format: text

text:
# Coherent 𝜔 photoproduction from nuclei and 𝜔 properties in nuclear matter.

## 1 Introduction

The properties of $`\rho `$ and $`\omega `$ mesons in a nuclear medium remain a topic of extreme interest since the first experimental observations by the CERES and HELIOUS-3 Collaborations report of a strong enhancement of low mass dilepton production from heavy ion collisions Agakichiev1 ; Agakichiev2 ; Agakichiev3 ; Masera . This enhancement was interpreted as a modification of the mass and width of the $`\rho `$ and $`\omega `$ mesons in dense hadronic matter Brown ; Li ; Rapp ; Cassing . The $`\rho `$-meson is quite broad in the vacuum and an additional increase of its width in the hadronic environment leads to the production of an almost continuum dilepton spectrum, which does not indicate a resonance structure and can be identified as perpetual spectral strength of the isovector current.

While the $`\rho `$-meson seems to be absolutely melted in matter, one expects that the signal of the $`\omega `$-meson might survive as a quasi resonance structure, even if being strongly modified. More dedicated experiments aim to isolate the $`\omega `$-meson produced in the nuclear medium either by high resolution dilepton spectroscopy Friese1 ; Friese2 ; Stroth ; Morrison ; Zajc  or through the $`\omega \pi ^0\gamma `$ decay mode Messchendorp ; Sibirtsev4 . The goal of these experiments is to identify the structure of the $`\omega `$-meson and to measure its width and pole position.

The experimental studies of the in-medium modification of the $`\omega `$-meson mass and width are planned by utilizing pion, proton and heavy ion beams with HADES at GSI Friese1 ; Friese2 ; Stroth ; Metag1 ; Bratkovskaya3 , with photo-nuclear reactions Messchendorp ; Effenberger  at TAPS Gabler  and Crystal Barrel Aker  at ELSA, by heavy ion collisions with PHENIX at RHIC Morrison ; Zajc  and with proton-nucleus reactions by ANKE at COSY Sibirtsev4 ; Golubeva .

Very recently the E325 experiment at KEK-PS reported Ozawa  new results on the $`\rho /\omega `$ meson modification in nuclear matter. It was found that the dielectron spectra produced in $`p+C`$ and $`p+Cu`$ collisions at proton beam energy of 12 GeV indicate a significant difference below the $`\omega `$-meson mass. The enhancement of the invariant mass spectra of electron-positron pairs produced from a heavy $`Cu`$ target might suggest the modification of the $`\rho /\omega `$ properties at normal nuclear density and for average $`\omega `$-meson momenta $`k_\omega `$1 GeV.

Obviously, apart of the heavy ion reactions, in the above mentioned experiments the $`\omega `$-mesons are produced with a certain laboratory momenta $`k_\omega `$, that can be quite large. The scale of $`k_\omega `$ is defined by the threshold of the relevant elementary $`\omega `$-meson production reaction. On the other hand, currently available predictions Tsushima ; Friman1 ; Friman2 ; Klingl1 ; Klingl2 ; Saito  for the in-medium modification of the $`\omega `$-meson mass and width are given for the $`\omega `$-meson at rest, i.e. for momentum $`k_\omega =0`$.

Only some estimates are presently available Friman1 ; Friman2 ; Lykasov  for finite $`k_\omega `$ and they are model dependent. These predictions Friman1 ; Friman2 ; Lykasov  are obtained by calculations of a forward $`\omega N`$ scattering amplitude $`f_{\omega N}`$, which within the low density approximation is related to an additional in-medium collisional $`\omega `$-meson width $`\mathrm{\Delta }\mathrm{\Gamma }_\omega `$ given by Lenz ; Dover ; Friman1 ; Friman2

$$\mathrm{\Delta }\mathrm{\Gamma }_\omega =4\pi \rho _B\frac{m_n+m_\omega }{m_nm_\omega }\mathrm{}f_{\omega N}(0),$$
(1)
where $`m_N`$ and $`m_\omega `$ are the free nucleon and $`\omega `$-meson masses, respectively, while $`\rho _B`$ is the baryon density. At the same time the $`\omega `$-meson mass shift $`\mathrm{\Delta }m_\omega `$ in nuclear matter is given by Lenz ; Dover ; Klingl2

$$\mathrm{\Delta }m_\omega =2\pi \rho _B\frac{m_n+m_\omega }{m_nm_\omega }\mathrm{}f_{\omega N}(0).$$
(2)
Here we propose phenomenological method to evaluate the forward $`\omega N`$ scattering amplitude by coherent $`\omega `$-meson photoproduction off nuclei and to reconstruct the momentum dependence of the $`\omega `$-meson mass and width in nuclear matter. Our paper is organized as follows. In Sect. 2 we provide the formalism of the $`\omega `$-meson coherent photoproduction from nuclei. The analysis of available data as well as an estimation of the incoherent background are given in Sect. 3. In Sect. 4 we evaluate the forward scattering amplitude and compare our results with the predictions from available sources. The section ends with summary of our results.

## 2 Coherent $`\omega `$ photoproduction

The total cross section of $`\omega `$-meson photoproduction from nuclei, $`\gamma +A\omega +X`$ is given by the sum of the coherent and incoherent processes. The differential cross section $`d\sigma _{\gamma A}^{coh}/dt`$ for the coherent $`\omega `$-meson photoproduction from nuclei can be written in the eikonal form as Ross ; Drell ; Kolbig ; Bauer ; Trefil1 ; Trefil2

$`{\displaystyle \frac{d\sigma _{\gamma A}^{coh}}{dt}}=|2\pi f_{\gamma N}^{dif}(0){\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dbJ_0(k_tb)b`$
$`\times {\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}dz\rho (b,z)\mathrm{exp}[ik_lz+i\chi _\omega (b)]|^2,`$ (3)
where an integration is performed over the impact parameter $`b`$ and $`z`$ coordinate along the direction of the incident photon, $`\rho (r=\sqrt{b^2+z^2})`$ is the nuclear density function normalized to the total nucleon number $`A`$, while $`k_l`$, $`k_t`$ are the longitudinal and transverse component of the momentum transfered to the nucleus, respectively, given by

$$k_l=k_\omega \sqrt{k_\omega ^2m_\omega ^2},k_t=2k_\omega \mathrm{sin}(\theta /2),$$
(4)
with $`k_\omega `$ is the total momentum, $`m_\omega `$ is the pole mass of the $`\omega `$-meson and $`t`$ is the squared four momentum transfered from the photon to $`\omega `$-meson. In (3) $`J_0`$ is the zeroth order Bessel function and $`\chi _\omega `$ is the nuclear phase shift related to $`\omega `$-meson distortion in nucleus. Within the $`t\rho `$-approximation this phase shift can be approximated by

$$\chi _\omega (b)=\frac{2\pi f_{\omega N}(0)}{k_\omega }\underset{z}{\overset{\mathrm{}}{}}\rho (b,y)𝑑y,$$
(5)
where $`f_{\omega N}(0)`$ is complex amplitude for the forward $`\omega N`$ elastic scattering. The imaginary part of $`f_{\omega N}(0)`$ is related by an optical theorem to the total cross section $`\sigma _{\omega N}`$ of the $`\omega N`$ interaction as

$$\mathrm{}f_{\omega N}(0)=\frac{k_\omega }{4\pi }\sigma _{\omega N}.$$
(6)
Now we introduce the ratio of the real to imaginary part of the forward scattering amplitude $`\alpha _\omega `$ defined as

$$\alpha _\omega =\frac{\mathrm{}f_{\omega N}(0)}{\mathrm{}f_{\omega N}(0)}$$
(7)
and finally express the phase shift $`\chi _\omega `$ as

$$\chi _\omega =\frac{\sigma _{\omega N}(i+\alpha _\omega )}{2}\underset{z}{\overset{\mathrm{}}{}}\rho (b,y)𝑑y.$$
(8)
The quantity $`f_{\gamma N}^{dif}(0)`$ in (3) is the forward diffractive $`\omega `$-meson photoproduction amplitude for a single nucleon related to the diffractive cross section $`d\sigma _{\gamma N}^{dif}/dt`$ at $`t=0`$ as

$$\frac{d\sigma _{\gamma N}^{dif}}{dt}|_{t=0}=|f_{\gamma N}^{dif}(0)|^2,$$
(9)
ant can be expressed in a vector dominance model VDM in terms of the forward $`\omega N`$ scattering amplitude $`f_{\omega N}(0)`$.

The total invariant amplitude $`_{\gamma \omega }`$ for the reaction $`\gamma +N+\omega N`$is given by a VDM as

$$_{\gamma \omega }=[\underset{V=\rho ,\varphi ,J/\mathrm{\Psi },\mathrm{}}{}\frac{\sqrt{\pi \alpha }}{\gamma _V}_{V\omega }]+\frac{\sqrt{\pi \alpha }}{\gamma _\omega }_{\omega \omega },$$
(10)
where the summation is performed over the available vector meson states $`V`$. The quantity $`\alpha `$ represents the fine structure constant, $`\gamma _V`$ denotes the photon coupling to the vector meson state and $`_{V\omega }`$ is the amplitude for the transition $`V+N\omega +N`$ on a nucleon. Moreover, in (10) $`\gamma _\omega =8.24\pm 0.24`$ is the photon-$`\omega `$ coupling constant evaluated Sibirtsev1  from the leptonic decays of the $`\omega `$-meson, and $`_{\omega \omega }`$ is an invariant amplitude for the elastic scattering $`\omega +N\omega +N`$.

In principle, (10) might also contain the sum given by a non-diagonal transition of non-vector meson states to the $`\omega `$-meson as well as coupling of these possible intermediate continuum states to the photon. However, only the second term of (10) provides diffractive $`\omega `$-meson photoproduction on a nucleon and enters (3).

The differential cross section $`\gamma +N\omega +N`$ is given by an invariant amplitude $``$ as

$$\frac{d\sigma _{\gamma N}}{dt}=\frac{|_{\gamma \omega }|^2}{64\pi sq_\gamma ^2},$$
(11)
where $`s`$ is the squared invariant collision energy and $`q_\gamma `$ is photon momentum in the center of mass system. Now, the forward diffractive $`\omega `$-meson photoproduction cross section on a single nucleon, $`d\sigma _{\gamma N}/dt`$ at $`t=0`$ can be expressed in terms of the $`\omega N`$ forward scattering invariant amplitude $`_{\omega \omega }(0)`$ as

$$\frac{d\sigma _{\gamma N}^{dif}}{dt}|_{t=0}=\frac{1}{64sq_\gamma ^2}\frac{\alpha }{\gamma _\omega ^2}|_{\omega \omega }(0)|^2.$$
(12)
Furthermore, the scattering amplitude $`f_{\omega N}`$ of (6) is related to the invariant scattering amplitude $`_{\omega \omega }`$ as

$$f_{\omega N}=\frac{1}{8\pi \sqrt{s}}\frac{k_\omega }{q_\omega }_{\omega \omega },$$
(13)
where $`k_\omega `$ and $`q_\omega `$ are the $`\omega `$-meson momenta in the laboratory and the $`\omega N`$ center of mass system, respectively.

Finally, substituting (6) and (13) into (12) the differential cross section for the forward $`\omega `$ meson photoproduction on a nucleon is given by

$$\frac{d\sigma _{\gamma N}^{dif}}{dt}|_{t=0}=\frac{\alpha }{16\gamma _\omega ^2}\frac{q_\omega ^2}{q_\gamma ^2}(1+\alpha _\omega ^2)\sigma _{\omega N}^2.$$
(14)
The differential cross section $`d\sigma _{\gamma A}^{coh}/dt`$ for the coherent $`\omega `$-meson photoproduction from nuclei can then be written as

$`{\displaystyle \frac{d\sigma _{\gamma A}^{coh}}{dt}}={\displaystyle \frac{\pi \alpha \sigma _{\omega N}^2}{4\gamma _\omega ^2}}{\displaystyle \frac{q_\omega ^2}{q_\gamma ^2}}|(i+\alpha _\omega ){\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑bJ_0(k_tb)b{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}𝑑z\rho (b,z)`$
$`\times \mathrm{exp}[ik_lz]exp\left[{\displaystyle \frac{\sigma _{\omega N}(i\alpha _\omega 1)}{2}}{\displaystyle \underset{z}{\overset{\mathrm{}}{}}}\rho (b,y)dy\right]|^2.`$ (15)
The parameters $`\alpha _\omega `$ and $`\sigma _{\omega N}`$ can be extracted from experimental data on coherent $`\omega `$-meson photoproduction from nuclei and through (6) and (7) should be converted to the real and imaginary part of the forward $`\omega N`$ scattering amplitude. Furthermore, the coherent photoproduction can be measured at various photon beam momenta, thus providing the evaluation of the momentum dependence of the complex forward scattering amplitude $`f_{\omega N}(0)`$. In this way the momentum dependent potential of the $`\omega `$-meson in nuclear matter can be reconstructed

Eq.(15) illustrates that the absolute value of the $`\omega `$-meson coherent photoproduction differential cross section is very sensitive to the size of the total $`\omega N`$ cross section, $`\sigma _{\omega N}`$, and thus to the imaginary part of the forward scattering amplitude $`\mathrm{}f_{\omega N}(0)`$. Therefore $`\mathrm{}f_{\omega N}(0)`$ can be uniquely fixed by measuring the coherent $`\omega `$-meson photoproduction from nuclei.

To investigate the sensitivity of the coherent $`\omega `$-meson photoproduction on the real part of the forward scattering amplitude we calculate the coherent $`\omega `$-meson photoproduction cross section from $`Cu`$ nuclei at a photon beam energy of $`E_\gamma `$=6.8 GeV. In the calculations we use the nuclear density function $`\rho (r)`$ given by a Wood-Saxon distribution as

$$\rho (r)=\frac{\rho _0}{1+\mathrm{exp}[(rR)/d]},$$
(16)
with the nuclear radius being parameterized as

$$R=1.28A^{1/3}0.76+0.8A^{1/3}\text{f}m,$$
(17)
and the diffusion parameter given by $`d=\sqrt{3}/3`$ fm. Furthermore, for the moment we fix $`\sigma _{\omega N}`$ at 35 mb, as will be motivated later.

Fig.1 shows $`d\sigma _{\gamma A}^{coh}/dt`$ calculated by (15) as a function of the square of four momentum transfer $`t`$. The solid line indicates the result calculated with the ratio of the real to imaginary part of the $`f_{\omega N}(0)`$ taken as $`\alpha _\omega `$=–0.4, while the dashed line is the result for $`\alpha _\omega `$=0.4. It is clear that the ratio $`\alpha _\omega `$ can be fixed only at the vicinity of the diffractive minima. Moreover, at the diffractive minima the calculated coherent $`\omega `$-meson photoproduction cross section indicates a substantial dependence on both, the absolute value and the sign of the ratio $`\alpha _\omega `$.

Furthermore, Fig.2 shows the differential cross section for the coherent photoproduction of $`\omega `$-mesons from a $`Cu`$ target calculated around the diffractive minima at $`t`$= –0.03 GeV2 as a function of ratio $`\alpha _\omega `$. We found that around the difractive minima the $`d\sigma _{\gamma Cu}^{coh}/dt`$ is very sensitive to the $`\alpha _\omega `$ and therefore, the real part of the forward scattering amplitude $`\mathrm{}f_{\omega N}(0)`$ can be directly measured by coherent $`\omega `$-meson photoproduction. It is worthwhile to notice that this method allows a determination of the sign of $`\mathrm{}f_{\omega N}(0)`$.

The solid circles in Fig.1 show experimental results for $`\omega `$-meson photoproduction at an average photon energy $`E_\gamma `$=6.8 GeV obtained by Behrend et al. Behrend  at the Cornell electron synchrotron. The produced $`\omega `$-mesons were detected through the $`\omega \pi ^+\pi ^{}\pi ^0`$ decay mode, thus the measurements Behrend  provide the data on inclusive $`\omega `$-meson photoproduction. The contributions from coherent and incoherent $`\omega `$-mesons photoproduction were not separated experimentally Behrend .

The experimental results on $`t`$-distribution shown by Fig.1 indicate a forward peak typical for the coherent $`\omega `$-meson photoproduction. The photoproduction at large $`|t|`$ is dominated by incoherent production mechanism. Fig.1 illustrates that our calculations based on (15) reasonably reproduce the experimental results Behrend  at small $`|t|`$ and manifistate the dominance of the coherent photoproduction at $`t`$–0.01 GeV2. The data Behrend  can be used only for the evaluation of $`\sigma _{\omega N}`$, or imaginary part of the forward $`\omega N`$ scattering amplitude.

## 3 Data analysis

As it was shown in Fig. 1 the comparison between our calculations and the experimental results requires an estimate of the contribution from the incoherent $`\omega `$-meson photoproduction. Moreover, the evaluation of the incoherent background is important for an determination of the experimental efficiency necessary for the separation of the coherent and incoherent $`\omega `$-meson photoproduction, especially near the difractive minima, where the coherent cross section is sensitive to the real part of the $`\mathrm{}f_{\omega N}(0)`$.

Actually, the results shown in Fig. 1 indicate that for $`Cu`$ targets and a photon energy of $`E_\gamma `$=6.8 GeV the contribution from incoherent $`\omega `$-meson photoproduction around $`t`$–0.03 GeV2 may exceed the coherent one by two or three orders of magnitude. Whether this large difference holds for other nuclear targets and photon energies requires an additional calculation of the incoherent $`\omega `$-meson photoproduction.

Taking into account only the leading term in $`t`$ given by $`\omega `$-meson photoproduction and the absorption in the nucleus the differential cross section $`d\sigma _{\gamma A}^{inc}/dt`$ for incoherent $`\omega `$-meson photoproduction is given by Kolbig

$`{\displaystyle \frac{d\sigma _{\gamma A}^{inc}}{dt}}={\displaystyle \frac{d\sigma _{\gamma N}^{tot}}{dt}}|_{t=0}\mathrm{exp}(at){\displaystyle \frac{2\pi }{\sigma _{\omega N}}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}b𝑑b`$
$`\times (1\mathrm{exp}\left[\sigma _{\omega N}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\rho (b,y)𝑑y\right]),`$ (18)
where $`a`$ is the slope of $`t`$-distribution from the elementary reaction $`\gamma +N\omega +N`$. Furthermore, the elementary forward $`\omega `$-meson photoproduction cross section $`d\sigma _{\gamma N}^{tot}/dt`$ at $`t=0`$ now should account not only for the diffractive process due to the $`\omega +N\omega +N`$ transition as is given by Eq.14, but also for the total $`\gamma N\omega N`$ amplitude of (10).

However, as will be shown later, at high energies (14) accounts for the substantial part of a total forward $`\gamma +N\omega +N`$ cross section. On the other hand, $`d\sigma _{\gamma N}^{tot}/dt`$ at $`t=0`$ as well as slope $`a`$ can be evaluated from experimental data LB .

Moreover, higher order corrections to (18) due to the single elastic scattering after the $`\omega `$-photoproduction on a target nucleon are proportional to the powers of $`\mathrm{exp}(at)`$ and might be important only at large $`|t|`$. However, this is beyond the region of our interest. We also do not consider photon shadowing, which might play certain role in the incoherent $`\omega `$-meson photoproduction.

The differential cross section for $`\gamma +Cu\omega +X`$ reaction at $`E_\gamma `$=6.8 GeV. is shown in Fig. 3. Solid circles represent experimental results from Ref. Behrend . The lines show the calculations performed with $`\sigma _{\omega N}`$=46 mb, the ratio $`\alpha _\omega `$=0 and a slope $`a`$=9.5 GeV-2 LB . Here we also take $`d\sigma _{\gamma N}^{tot}/dt`$ at $`t=0`$ as given by (14). The dashed line in Fig. 3 shows our result for the differential cross section for coherent $`\omega `$-meson photoproduction calculated by (15). The dotted line indicates the result for the incoherent $`\omega `$-meson photoproduction as calculated by (18). The solid line is the sum of coherent and incoherent processes and reasonably describes the data Behrend .

Fig. 4 shows the differential cross section for $`\omega `$-meson photoproduction from $`Be`$ at a photon energy $`E_\gamma `$=6.8 GeV. Again, the experimental results are taken from Ref. Behrend  and the lines show the calculations with $`\sigma _{\omega N}`$=46 mb, $`\alpha _\omega `$=0 and $`a`$=9.5 GeV-2 LB . The results for coherent and incoherent $`\omega `$-meson photoproduction are shown by the dashed and dotted lines, respectively. The solid line indicates their sum and describes the data reasonably well. Note that the position of the diffractive minima for light nuclei lies at substantially larger $`|t|`$ Sibirtsev2 .

In the calculations for $`Be`$ we use the density function given by Dalkarov

$$\rho (r)=\frac{1}{(R\sqrt{\pi })^3}\left[4+\frac{2(A4)r^2}{3R^2}\right]\mathrm{exp}\left[\frac{r^2}{R^2}\right],$$
(19)
where $`R`$=1.58 fm.

The $`\omega `$-meson photoproduction from nuclei at average photon energies $`E_\gamma `$=8.2 GeV was studied experimentally by Aramson et al. Abramson  with results compatible to their previous measurements given at $`E_\gamma `$=6.8 GeV in Ref. Behrend .

In addition, the $`\omega `$-meson nuclear photoproduction at average photon energies $`E_\gamma `$=5.7 GeV was measured at DESY Braccini  by detecting the $`\omega \pi ^0\gamma `$ decay. The coherent and incoherent processes were not separated and the data were published without absolute normalization.

Fig.5 shows the differential cross section for $`\omega `$-meson photoproduction from $`C`$ and $`Al`$ targets at photon beam energy $`E_\gamma `$=5.7 GeV. The DESY experimental results are taken from Ref. Braccini  and normalized by our calculations. The solid lines show the sum of the coherent and incoherent photoproduction cross sections. The dotted lines show the calculations for incoherent $`\omega `$-meson photoproduction, which is dominant at large $`|t|`$. The contribution from coherent $`\omega `$-meson photoproduction is given by the dashed lines. The calculations were performed using (15) and (18) with parameters $`\sigma _{\omega N}`$=46 mb, $`\alpha _\omega `$=0 and $`a`$=9.5 GeV-2 LB .

The data Braccini  are quite reasonably reproduced by the calculations. However, they are insensitive to the ratio $`\alpha `$. Our choice of $`\alpha _\omega `$=0 for the comparison between the calculations and the data Behrend ; Braccini  is indeed arbitrary. We conclude that the ratio $`\alpha `$ as well as the real part of the forward $`\omega N`$ scattering amplitude $`\mathrm{}f_{\omega N}(0)`$ can not be constrained by the available data Behrend ; Abramson ; Braccini , since these measurements did not isolate the coherent photoproduction.

The total coherent and incoherent $`\omega `$-meson photoproduction cross sections from complex nuclei at $`E_\gamma `$=3.9 GeV were measured by Brodbeck et al. Brodbeck  and extrapolated to $`t`$=0. Experimental results for $`\omega `$-meson photoproduction from nuclei at $`E_\gamma `$=6.8 GeV and at $`t`$=0 were also reported by Behrend et al. Behrend1 .

Fig.6 shows results for $`\omega `$-meson photoproduction cross section at $`t`$=0 as a function of the target mass number. The solid circles indicate data Brodbeck  collected at $`E_\gamma `$=3.9 GeV, while the squares show the results Behrend1  measured at photon energy $`E_\gamma `$=6.8 GeV. The data do not show any dependence on the photon beam energy.

The solid line in Fig.6 indicates our calculations with $`\sigma _{\omega N}`$=37 mb and ratio $`\alpha _\omega `$=0. We note that our results do not depend on the ratio $`\alpha `$. The dashed line in Fig. 6 shows the $`A^2`$-dependence which is given by (15) assuming that $`\sigma _{\omega N}`$=0 and neglecting the $`k_l`$-dependence. The dashed line is arbitrary normalized at $`A=Be`$.

Moreover, the difference between the solid and dashed lines in the Fig. 6 is due to the finite value of $`\sigma _{\omega N}`$. Therefore, the data Brodbeck ; Behrend1 ; Behrend1  on forward $`\omega `$-meson photoproduction provide a reasonable evaluation of the imaginary part of the $`\omega N`$ forward scattering amplitude $`\mathrm{}f_{\omega N}(0)`$ .

## 4 Evaluation of the forward scattering amplitude

The comparison between our calculations and experimental results Behrend ; Abramson ; Braccini ; Brodbeck  allows us to extract the imaginary part of the forward $`\omega N`$ scattering. Finally the data on $`\omega `$-meson photoproduction from nuclei were fitted in order to evaluate $`\sigma _{\omega N}`$. Although we found that the available data are absolutely insensitive to the ratio $`\alpha _\omega `$ we allow its variation by a minimization procedure Minuit .

The results of the minimization are shown in the Table. 1. None of the data sets allows to extract a specific value for the ratio $`\alpha _\omega `$. Some of the experimental results allow an extraction of $`\sigma _{\omega N}`$ with sufficiently large uncertainty.

However, we should clarify, that experimental results are insensitive to the ratio $`\alpha _\omega `$ in case it is less than one. As is shown by (15) the large values of $`\alpha _\omega `$ should immediately effect our results. Roughly, a large ratio $`\alpha _\omega `$ should cause a decrease of the total $`\omega N`$ cross section $`\sigma _{\omega N}`$. The fitting procedure does not indicate such a tendency. Our analysis suggest that $`\alpha _\omega <1`$ or $`\mathrm{}f_{\omega N}(0)>`$$`\mathrm{}f_{\omega N}(0)`$.

Moreover, the $`\omega `$-meson photoproduction data on $`A`$-dependence at $`t`$=0 shown in the Fig.6 show substantial sensitivity to $`\sigma _{\omega N}`$. Obviously, the situation might improve substantially when separate data for coherent photoproduction will be available.

Final results for $`\sigma _{\omega N}`$ are converted by (6) to imaginary part of $`f_{\omega N}(0)`$ and are shown by the solid circles in Fig.7 as a function of the total laboratory $`\omega `$-meson energy $`E_\omega `$. In the case, when the results for $`\sigma _{\omega N}`$ are evaluated from different data sets but at the same photon energy, we average them in calculating the imaginary part of the forward scattering amplitude.

Now our results can be compared with estimations given by (6) and (14) and available data on forward $`\omega `$-meson photoproduction on a free nucleon. The squares in Fig.7 show the $`\mathrm{}f_{\omega N}(0)`$ obtained from the data LB  on $`\gamma +p\omega +p`$ reaction cross section at $`t`$=0 and neglecting the ratio $`\alpha _\omega `$ in (14).

An arrow at Fig. 7 indicates the lower limit for the imaginary part of the $`\omega N`$ scattering amplitude at threshold $`E_\omega =m_\omega `$ given at the pole mass of $`\omega `$-meson by

$$\mathrm{}f_{\omega N}(0)=0.164\text{f}m,$$
(20)
as recently evaluated Hanhart  from $`\pi ^{}+p\omega +n`$ measurements Karami .

Moreover, within an additive quark model the forward $`\omega N`$ scattering amplitude is given as the average of the forward $`\pi ^{}N`$ and $`\pi ^+N`$ amplitudes Lipkin ; Kajantie ; Donnachie1 . Using the Donnachie and Landshoff results Donnachie2  for a Regge theory fit to the $`\pi N`$ total cross section the complex forward $`\omega N`$ scattering amplitude can be written as

$`f_{\omega N}(0)={\displaystyle \frac{k_\omega }{4\pi \mathrm{}c}}[(\mathrm{\hspace{0.17em}0.173}\left[{\displaystyle \frac{s}{s_0}}\right]^ϵ2.726\left[{\displaystyle \frac{s}{s_0}}\right]^\eta )`$
$`+i(\mathrm{\hspace{0.17em}1.359}\left[{\displaystyle \frac{s}{s_0}}\right]^ϵ+3.164\left[{\displaystyle \frac{s}{s_0}}\right]^\eta )]\text{fm},`$ (21)
where the $`\omega `$-meson momentum $`k_\omega `$ and squared invariant collision energy $`s`$ are given in GeV/c and $`s_0`$=1 GeV2. The effective powers $`ϵ`$=0.08 and $`\eta `$=0.45 are given by the pomeron and $`\rho `$, $`\omega `$, $`f`$, $`a`$ exchanges, respectively. The solid line in Fig. 7 shows the prediction by additive quark model calculated by (21). We find that the prediction from the additive quark model and the Regge theory systematically underestimates other results.

Although we do not consider the results for ratio $`\alpha _\omega `$ listed in the Table 1 at some confidence level, it seems worthwhile to discuss a possible estimate for the ratio of the real to imaginary part of the forward $`\omega N`$ scattering amplitude. Moreover, as is illustrated by Fig.2 the magnitude of the coherent differential $`\omega `$-meson photoproduction cross section at diffractive minima substantially depends on the both sign and size of the $`\alpha _\omega `$ and thus the estimate for this ratio may be useful for valuation of an experimental accuracy necessary for measurements around the diffractive minima.

Very rough estimates of the ratio $`\alpha _\omega `$ may be given by an additive quark model through (21) and is shown in Fig.8 by the solid line as a function of total $`\omega `$-meson laboratory energy $`E_\omega `$. The solid circles in the Fig.8 show the experimental results Foley  for the ratio of the real to imaginary part for the $`\pi ^{}p`$ forward scattering amplitude while the squares indicate the data for $`\pi ^+p`$ scattering. The $`\pi N`$ data were taken at the same invariant collision energies $`\sqrt{s}`$.

Furthermore, the results from the coupled channel analysis Friman1  of the meson-nucleon scattering data predicts at the $`\omega N`$ threshold

$$f_{\omega N}=0.5+i\mathrm{\hspace{0.17em}0.2}\text{fm},$$
(22)
providing the negative ratio $`\alpha _\omega `$=–2.5. Moreover, most recently within the same model approach Friman2  it was deduced that $`\mathrm{}f_{\omega N}`$=0.013 fm, while the real part $`\mathrm{}f_{\omega N}`$ stands the same as in Eq.22 resulting finally in the ratio $`\alpha _\omega `$=–38.5.

At the same time, the calculations Klingl1  based on an effective Lagrangian approach combined with chiral SU(3) dynamics predicts the forward $`\omega N`$ scattering amplitude as

$$f_{\omega N}=3.34+i\mathrm{\hspace{0.17em}2.1}\text{fm},$$
(23)
with a positive defined ratio $`\alpha _\omega `$=0.62. More recent results of the calculation Klingl2  of the scattering amplitude is given as

$$f_{\omega N}=1.6+i\mathrm{\hspace{0.17em}0.3}\text{fm},$$
(24)
with positive and large ratio $`\alpha _\omega `$=5.3.

Theoretical predictions available for $`\alpha _\omega `$ differ substantially both in sigh and the absolute value of the ratio. Obviously, the coherent $`\omega `$-meson photoproduction from nuclei can not directly resolve the existing discrepancies between available predictions for the forward $`\omega N`$ scattering amplitude exactly at $`E_\omega =m_\omega `$. However the $`\gamma +A\omega +A`$ reactions might be considered as a model independent measurements of $`f_{\omega N}`$ at final $`\omega `$-meson momenta as well as provide an orientation for analytical continuation of the forward scattering amplitude close to the threshold Eletsky ; Kondratyuk ; Sibirtsev3 .

It is important to note that both, the coherent and incoherent cross sections for forward $`\omega `$-meson photoproduction from nuclei will be influenced when large values of $`\alpha >`$1 are used.

## 5 Conclusions

The coherent $`\omega `$-meson photoproduction from nuclei is studied as a possible phenomenological method to evaluate the complex forward $`\omega N`$ scattering amplitude at finite $`\omega `$-meson momenta.

We found that the real part of the forward scattering amplitude $`\mathrm{}f_{\omega N}(0)`$ can be only fixed by the data on coherent $`\omega `$-meson photoproduction around the diffractive minima. At the same time, the imaginary part of the forward scattering amplitude $`\mathrm{}f_{\omega N}(0)`$ can be well evaluated by an absolute value of the cross section at small $`|t|`$, where the contribution from incoherent processes is relatively small. Moreover, the sign and the magnitude of the ratio $`\alpha _\omega `$ of the real to imaginary part of the forward $`\omega N`$ scattering amplitude can be identified by the magnitude of the coherent photoproduction cross section at diffractive minima.

We analyze available experimental data Behrend ; Abramson ; Braccini ; Brodbeck ; Behrend1  on $`\omega `$-meson photoproduction from nuclei. Unfortunately the contributions from coherent and incoherent processes are not separated experimentally. We find reasonable agreement between the data Behrend ; Abramson ; Braccini ; Brodbeck ; Behrend1  and our calculations including both, coherent and incoherent $`\omega `$-meson photoproduction. By fitting the experimental results we evaluate the imaginary part of the forward $`\omega N`$ scattering amplitude. Within our approach we do not variate the $`\omega `$-photon coupling constant, but fix it by the dileptonic decay. We find that data, which are not separated into coherent and incoherent production, can not provide reliable results for the real part of the scattering amplitude $`\mathrm{}f_{\omega N}(0)`$.

Finally, our results for $`\mathrm{}f_{\omega N}(0)`$ are compared with predictions given by a vector dominance model relation between the forward $`\gamma +N\omega +N`$ reaction cross section and imaginary part of the forward $`\omega N`$ scattering amplitude. The experimental data LB  for forward $`\omega `$-meson photoproduction from the proton provide an $`\mathrm{}f_{\omega N}(0)`$ that agrees reasonably well with our estimates.

Furthermore, the predictions by an additive quark model together with a Regge model fit to the data for hadronic cross sections systematically underestimate our results for $`\mathrm{}f_{\omega N}(0)`$ as well as that extracted from the $`\gamma +p\omega +p`$ reaction. This discrepancy might indeed indicate that, while experimental results on exclusive $`\rho `$, $`\varphi `$ and $`J/\mathrm{\Psi }`$ photoproduction excellently confirm Donnachie1 ; Donnachie2 ; Sibirtsev3  the Regge theory, the $`\omega `$-meson photoproduction could be considered as an exception. It is important that $`\mathrm{}f_{\omega N}(0)`$ is evaluated precisely from available data on the $`A`$-dependence of $`\omega `$-meson photoproduction at $`t`$=0.

We also discuss the current status of the predictions for the real part of the $`\omega N`$ scattering amplitude and detect a strong disagreement between the results from different models.

###### Acknowledgements.

This work was performed in part under the auspices of the U. S. Department of Energy under contract No. DE-FG02-93ER40756 with the Ohio University.


text:
# Pairbreaking Without Magnetic Impurities in Disordered Superconductors

## Abstract

We study analytically the effects of inhomogeneous pairing interactions in short coherence length superconductors, using a spatially varying Bogoliubov-deGennes model. Within the Born approximation, it reproduces all of the standard Abrikosov-Gor’kov pairbreaking and gaplessness effects, even in the absence of actual magnetic impurities. For pairing disorder on a single site, the T-matrix gives rise to bound states within the BCS gap. Our results are compared with recent scanning tunneling microscopy measurements on Bi2Sr2CaCu2O8+δ with Zn or Ni impurities.

The effect of disorder in superconductors has long been a subject of considerable interest. A generally accepted physical picture is that magnetic impurities destroy superconductivity by locally breaking the pairs,AG ; deg64 ; amb65 ; pa69  whereas non-magnetic impurities are not pairbreaking, according to Anderson’s theorem.pa59  This is true for an isotropic s-wave BCS superconductor, in which the order parameter is uniform and momentum independent. Since in most high transition temperature $`T_c`$ cuprates, the suppression of $`T_c`$ with Zn or Ni doping is comparable, mah94  there were proposals to explain this in terms of $`d`$-wave superconductivity.sal96  However, nuclear magnetic resonance (NMR) experiments indicated that the nominally non-magnetic Zn++ ions polarize the spin background in the CuO2 planes upon substitution of the nominally $`S=1/2`$ Cu++ sites.mah94 ; julien00  Moreover, recent scanning tunneling microscopy (STM) measurements directly above the Zn or Ni impurity sites observed strong resonance peaks.pan00 ; hudson01  Very recently, several groups noticed from STM measurements that underdoped Bi2Sr2Cu2O8+δ (BSCCO) and Bi2-xPbxSr2Cu2O8+δ are extremely disordered on a scale of a few nm. pan01 ; lang02 ; HFK ; Cren  This disorder is characterized by two gaps, one corresponding to the superconducting gap, with characteristic superconducting peaks, and a non-superconducting gap.

There is now a large body of evidence that the pseudogap observed in cuprates above $`T_c`$ is not superconducting. The most convincing of these experiments determined that the magnetic field dependence of the resistivity, NMR spin-lattice relaxation rate, magnetization, and the gaps observed in intrinsic tunneling experiments for the superconducting and pseudogaps are qualitatively different.Shibauchi ; SW ; lang02 ; Krasnov ; Zheng ; Klemm  In particular, the pseudogap regime is field independent until one reaches the Zeeman field for breaking up chargeless spin zero pairs.Shibauchi  This is precisely consistent with the pseudogap being particle-hole pairs, which could be either of the charge-density wave (CDW) or spin-density wave (SDW) form. As in one dimension, such excitations are expected to arise from repulsive interactions between like charges, whereas superconductivity can arise when the interactions between like charges are attractive.

Thus, if indeed the disorder involves a percolation problem between superconducting and density-wave regions on a nanometer scale, then the phase coherence can only arise from Josephson coupling of the superconducting grains. In addition, one would expect the $`c`$-axis tunneling to be very incoherent, as inferred in BSCCO from $`c`$-axis twist Josephson junction experiments. Li

In this letter, we assume the superconductor is electronically disordered on the scale of the coherence length. We further assume the essential ingredient in the disorder is not one of impurities, but rather disorder in the pairing interaction itself. Thus, we expect that $`T_c`$ varies from site to site, as does the resulting order parameter amplitude.OWK  We treat this type of disorder using a Bogoliubov-deGennes procedure, assuming that the order parameter amplitude varies locally. HS ; Nandini

Here we show that in the Born approximation, this problem maps exactly onto that of pairbreaking in a superconductor, with all of the features of that model, including gapless superconductivity. AG ; deg64 ; pa69 ; amb65  At a particular defect site, the T-matrix gives rise to bound states within the gap, even without magnetic impurities.

We use the Nambu representation,

$`\mathrm{\Psi }(𝐫)`$ $`=`$ $`\left(\begin{array}{c}c_{}(𝐫)\\ c_{}(𝐫)\\ c_{}^{}(𝐫)\\ c_{}^{}(𝐫)\end{array}\right),`$
$`\mathrm{\Psi }^{}(𝐫)`$ $`=`$ $`\left(\begin{array}{cccc}c_{}^{}(𝐫),& c_{}^{}(𝐫),& c_{}(𝐫),& c_{}(𝐫)\end{array}\right),`$ (1)
where $`c_{}(𝐫)`$ \[$`c_{}^{}(𝐫)`$\] annihilates \[creates\] a quasiparticle with spin eigenstate $``$ \[$``$\] at the position $`𝐫`$. We set $`\mathrm{}=c=k_B=1`$. The Hamiltonian under study is $`H=H_0+H_1+H_2+H_3`$, where

$`H_0`$ $`=`$ $`{\displaystyle d^3𝐫\mathrm{\Psi }^{}(𝐫)[\widehat{\xi }(𝐫)\rho _3\sigma _0+\mathrm{\Delta }_0\rho _2\sigma _2]\mathrm{\Psi }(𝐫)},`$ (2)
$`H_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle d^3𝐫\mathrm{\Psi }^{}(𝐫)\widehat{V}_i(𝐫)\mathrm{\Psi }(𝐫)},`$ (3)
$`\widehat{V}_1(𝐫)`$ $`=`$ $`U_1(𝐫)\rho _3\sigma _0,`$ (4)
$`\widehat{V}_2(𝐫)`$ $`=`$ $`U_2(𝐫)𝐒\stackrel{}{\alpha }/[S(S+1)]^{1/2}`$
$`\stackrel{}{\alpha }`$ $`=`$ $`\widehat{𝐱}\rho _3\sigma _1+\widehat{𝐲}\rho _0\sigma _2+\widehat{𝐳}\rho _3\sigma _3,`$ (5)
$`\widehat{V}_3(𝐫)`$ $`=`$ $`U_3(𝐫)\rho _2\sigma _2,`$ (6)
where $`\rho _j\sigma _j^{}\rho _j\sigma _j^{}`$ is a rank 4 tensor composed of two Pauli matrices for $`j,j^{}=1,2,3`$ and $`\rho _0`$, $`\sigma _0`$ are rank 2 identity matrices, respectively. $`H_0`$ is the Bogoliubov-deGennes version of the BCS Hamiltonian, with momentum space quasiparticle energy dispersion $`\xi _𝐤`$ relative to the Fermi energy $`\mu `$, $`\mathrm{\Delta }_0(T)`$ is the real bare uniform BCS order parameter, and $`H_1`$ and $`H_2`$ are the interactions due to scattering off random non-magnetic and magnetic impurities with effective potentials $`U_1(𝐫)`$ and $`U_2(𝐫)/[S(S+1)]^{1/2}`$ , respectively. $`H_3`$ with effective potential $`U_3(𝐫)`$ is the new interaction arising from random variations in the pairing interaction.OWK ; Nandini  In $`H_2`$, $`𝐒`$ and $`S`$ are the spin vector and quantum number of the magnetic impurities, respectively, and $`\stackrel{}{\alpha }`$ represents the quasiparticle spin eigenvector. We assume the spatial averages of all of the random potentials vanish, $`U_i(𝐫)=0`$ for $`i=1,2,3`$. In the absence of all defects, the order parameter $`\mathrm{\Delta }_0(T)=Vc_{}(𝐫)c_{}(𝐫)`$ satisfies the standard BCS gap equation,

$`\mathrm{\Delta }_0`$ $`=`$ $`VT{\displaystyle \underset{|\omega _n|\omega _0}{}}{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\mathrm{Tr}[\rho _2\sigma _2\widehat{G}_0(𝐤,\omega _n)]},`$
$`\widehat{G}_0^1(𝐤,\omega _n)`$ $`=`$ $`i\omega _n\rho _0\sigma _0\xi _𝐤\rho _3\sigma _0\mathrm{\Delta }_0\rho _2\sigma _2,`$ (7)
where $`\widehat{G}_0`$ is the bare Green’s function matrix, $`V<0`$ is the uniform (BCS) part of the pairing interaction, $`N(0)`$ is the single-spin quasiparticle density of states, $`\omega _0`$ is a BCS-like cutoff, and $`\omega _n`$ are the Matsubara frequencies.

We assumed a real bare uniform order parameter $`\mathrm{\Delta }_0`$, and restricted our consideration in $`H_3`$ to spatial fluctuations of the amplitude of $`\mathrm{\Delta }_0`$. The model can also treat spatial fluctuations of the phase of $`\mathrm{\Delta }_0`$ by letting $`\mathrm{\Delta }_0\rho _1\sigma _2`$ and $`U_3(𝐫)\rho _2\sigma _2`$ be generalized to $`\mathrm{\Delta }_{01}\rho _1\sigma _2+\mathrm{\Delta }_{02}\rho _2\sigma _2`$ and $`U_{31}(𝐫)\rho _1\sigma _2+U_{32}(𝐫)\rho _2\sigma _2`$, respectively.

Our main interest lies in studying $`H_3`$. Using quantum Monte Carlo techniques to study a two-dimensional square lattice with an on-site attractive Hubbard pairing interaction in $`H_0`$,Nandini  Ghosal et al. obtained interesting results in excellent qualitative agreement with those obtained from STM measurements.lang02  We also consider $`H_1`$ and $`H_2`$ for comparison, because the combination of one or both of them with $`H_3`$ can lead to interesting novel behavior. In the Born approximation, these interactions add or subtract in a simple fashion. However, in the T-matrix approximation for a single defect site, these interactions interact in a highly non-trivial manner.

In the self-consistent Born approximation, the quasiparticle self-energy $`\widehat{\mathrm{\Sigma }}=\widehat{\mathrm{\Sigma }}_1+\widehat{\mathrm{\Sigma }}_2+\widehat{\mathrm{\Sigma }}_3`$, where

$`\widehat{\mathrm{\Sigma }}_i(𝐤,\omega _n)`$ $`=`$ $`n_i{\displaystyle \underset{𝐤^{}}{}}\widehat{V}_i(𝐤𝐤^{})\widehat{G}(𝐤^{},\omega _n)\widehat{V}_i(𝐤^{}𝐤),`$ (8)
$`\widehat{G}^1(𝐤,\omega _n)`$ $`=`$ $`\widehat{G}_0^1(𝐤,\omega _n)\widehat{\mathrm{\Sigma }}(𝐤,\omega _n),`$ (9)
$`\widehat{V}_i(𝐤),U_i(𝐤)`$ are the spatial Fourier transforms of $`\widehat{V}_i(𝐫),U_i(𝐫)`$, respectively. Neglecting any possible anisotropy arising from Fermi surface integrations, the effective rates of the three processes are

$`1/\tau _i`$ $`=`$ $`2\pi n_iN(0)|U_i(𝐤_F)|^2,`$ (10)
where $`n_i`$ is the density of defects of type $`i`$.

As in the usual pairbreaking theory, AG ; deg64 ; amb65 ; pa69 , $`\widehat{G}`$ has the same form as does $`\widehat{G}_0`$, except that $`\omega _n`$ and $`\mathrm{\Delta }_0`$ are replaced by their renormalized equivalents $`\stackrel{~}{\omega }_n`$ and $`\stackrel{~}{\mathrm{\Delta }}`$, respectively. We then obtain the standard equations for the renormalized gap and Matsubara frequency,

$`\stackrel{~}{\omega }_n`$ $`=`$ $`\omega _n+(1/\tau _1+1/\tau _{\mathrm{pb}}){\displaystyle \frac{\stackrel{~}{\omega }_n}{2[\stackrel{~}{\omega }_n^2+\stackrel{~}{\mathrm{\Delta }}^2]^{1/2}}},`$ (11)
$`\stackrel{~}{\mathrm{\Delta }}`$ $`=`$ $`\mathrm{\Delta }_0+(1/\tau _11/\tau _{\mathrm{pb}}){\displaystyle \frac{\stackrel{~}{\mathrm{\Delta }}}{2[\stackrel{~}{\omega }_n^2+\stackrel{~}{\mathrm{\Delta }}^2]^{1/2}}},`$ (12)
$`1/\tau _{\mathrm{pb}}`$ $`=`$ $`1/\tau _2+1/\tau _3`$ (13)
is the total pairbreaking rate. The new physics arises from $`H_3`$. Evidently, within the self-consistent Born approximation, the effects of the random interactions are exactly equivalent to those of magnetic impurities.

Using standard pairbreaking theory, AG ; deg64 ; amb65 ; pa69  one finds

$`{\displaystyle \frac{\omega _n}{\mathrm{\Delta }_0}}`$ $`=`$ $`u(1{\displaystyle \frac{\zeta }{\sqrt{1+u^2}}}),`$ (14)
$`u`$ $`=`$ $`\stackrel{~}{\omega }_n/\stackrel{~}{\mathrm{\Delta }},`$ (15)
$`\zeta `$ $`=`$ $`1/(\tau _{pb}\mathrm{\Delta }_0).`$ (16)
The spatial average gap $`\mathrm{\Delta }(T)`$ is then

$`\mathrm{\Delta }`$ $`=`$ $`|V|T{\displaystyle \underset{|\omega _n|\omega _0}{}}{\displaystyle \frac{d^2𝐤}{(2\pi )^3}\mathrm{Tr}[\rho _2\sigma _2\widehat{G}(𝐤,\omega _n)]},`$ (17)
$`=`$ $`\pi |V|N(0)T{\displaystyle \underset{|\omega _n|\omega _0}{}}{\displaystyle \frac{1}{\sqrt{1+u^2}}},`$
leading to the standard equation for $`T_c/T_{c0}=t`$,

$`0`$ $`=`$ $`\mathrm{ln}(t)+\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\alpha _{\mathrm{pb}}}{2\pi t}}\right)\psi \left({\displaystyle \frac{1}{2}}\right),`$ (18)
$`\alpha _{\mathrm{pb}}`$ $`=`$ $`1/(\tau _{\mathrm{pb}}T_{c0}),`$ (19)
where $`\psi (x)`$ is the digamma function. For small $`\alpha _{\mathrm{pb}}`$, $`T_cT_{c0}\pi /4\tau _{\mathrm{pb}}`$. We note that $`T_c/T_{c0}`$ can be suppressed to zero even without any magnetic impurities, for $`1/\tau _31/\tau _{3c}=\pi T_{c0}/2\gamma `$, where $`\gamma =1.781`$ is the exponential of Euler’s constant. In addition, the superconductivity becomes gapless for $`1>\tau _{3c}/\tau _312e^{\pi /4}0.912`$. Thus, even an isotropic, $`s`$-wave superconductor can become gapless, as observed in the cuprates with STM.pan00 ; hudson01 ; HFK

This can only occur in short coherence length superconductors with strong local inhomogeneities in the pairing interaction, as is a likely explanation for the vanishing of the $`T_c`$ in the highly underdoped region of the cuprate phase diagram, although that region is also complicated by the simultaneous appearance of local SDW/CDW order at the non-superconducting regions not included in this calculation.lang02

In order to make direct comparison with STM experiments, we solve the T-matrix for a single defect site. We assume the site has all three types of defects associated with it. We approximate the effects of the magnetic impurity by assuming that its spin behaves classically.shi68  Then, the T-matrix equation can be solved exactly,

$`\widehat{T}(\omega _n)`$ $`=`$ $`{\displaystyle \frac{\widehat{V}(0)}{\widehat{1}\widehat{g}_0(\omega _n)\widehat{V}(0)}},`$ (20)
$`\widehat{V}(0)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}\widehat{V}_i(0),`$ (21)
$`\widehat{g}_0(\omega _n)`$ $`=`$ $`{\displaystyle \frac{\pi N(0)}{[\omega _n^2+\mathrm{\Delta }_0^2]^{1/2}}}[i\omega _n\rho _0\sigma _0+\mathrm{\Delta }_0\rho _2\sigma _2],`$ (22)
is the bare Green’s function at the origin, $`\widehat{1}=\rho _0\sigma _0`$ is the rank 4 identity matrix, and the effects of a finite quasiparticle energy bandwidth have been neglected. Bound states within the gap at $`T=0`$ at the frequency $`\omega `$ are obtained from

$`\mathrm{det}|\widehat{1}\widehat{g}_0(i\omega )\widehat{V}(0)|`$ $`=`$ $`0.`$ (23)
Solving Eqs. (21) - (23) exactly, we generally obtain four bound states within the gap at $`\omega =\overline{\omega }_b\mathrm{\Delta }_0`$, where

$`\overline{\omega }_b`$ $`=`$ $`\pm [A+sR]^{1/2}/B,`$ (24)
$`A`$ $`=`$ $`16v_2^2v_3^2+a_+^2(a_{}^2+4v_1^2),`$ (25)
$`B`$ $`=`$ $`a_+^2+4v_2^2,`$ (26)
$`C`$ $`=`$ $`a_+^24v_3^2,`$ (27)
$`R`$ $`=`$ $`[A^2B^2C^2]^{1/2},`$ (28)
$`a_\pm `$ $`=`$ $`v_1^2+v_3^2v_2^2\pm 1,`$ (29)
$`v_i`$ $`=`$ $`\pi N(0)U_i(0),`$ (30)
and $`s=\pm 1`$.

It is useful to rewrite Eq. (24) for the three special cases of two defects only. For non-magnetic defects, $`v_2=0`$, there are only two bound states symmetric about zero bias,

$`\overline{\omega }_b`$ $``$ $`\pm {\displaystyle \frac{\left([(1+v_3)^2+v_1^2][(1v_3)^2+v_1^2]\right)^{1/2}}{(1+v_1^2+v_3^2)}}.`$ (31)
For the trivial case of a site with just a non-magnetic impurity, Eq. (31) shows that there are no poles inside the gap. However, for a pairing interaction defect alone, there are poles at $`\overline{\omega }_b=\pm |1v_3^2|/(1+v_3^2)`$.

For impurities alone, $`v_3=0`$, there are also only two bound states symmetric about zero bias, at

$`\overline{\omega }_b`$ $``$ $`\pm {\displaystyle \frac{|1+u_+u_{}|}{[(1+u_+^2)(1+u_{}^2)]^{1/2}}},`$ (32)
where $`u_\pm =v_1\pm v_2`$. For $`v_1=0`$, this reduces to the result of Shiba for a classical magnetic impurity, $`\overline{\omega }_b=\pm |1v_2^2|/(1+v_2^2)`$, even though we did not average over the classical spin direction before summing the T-matrix.shi68

The most interesting case arises when both $`v_2,v_30`$ and $`v_3`$ and $`v_2v_3`$. In the limit $`v_1=0`$, there are four bound states symmetric about zero bias at

$`\overline{\omega }_b`$ $``$ $`\pm {\displaystyle \frac{(v_+^2v_{}^2+s|1v_+^2v_{}^2|)}{(1+v_+^2)(1+v_{}^2)}},`$ (33)
where $`s=\pm 1`$, $`v_\pm =v_2\pm v_3`$. If $`v_2,v_30`$ and $`v_2v_3`$, there are four bound states, exhibiting reflection symmetry about zero bias. Otherwise, if either $`v_2`$ or $`v_3=0`$, or if $`v_2=v_3`$, there are only two bound states that are symmetric about zero bias. We note that Eq. (33) for either $`v_2=0`$ or $`v_3=0`$ reduces to Eq. (31) and (32) with $`v_1=0`$, respectively. In addition, for $`v_2=v_3`$, it reduces to Eq. (32) with $`v_1=0`$ and $`v_22v_2`$, etc. Moreover, setting $`v_2=v_3`$ in Eq. (24) with $`v_1`$ arbitrary leads to only two bound states symmetric about zero bias, at $`\overline{\omega }_b=\pm [|b_{}|/b_+]^{1/2}`$, where $`b_\pm =(1+v_1^2)^2\pm 4v_2^2`$. Thus, we conclude that in order to obtain four bound states, two on each side of zero bias, one requires $`0v_2v_30`$.

When the defect is a quantum spin with a single component normal to the surface, the spin operator $`S_z`$ commutes with the Hamiltonian, and the spin states are easily described by $`|SM`$, with $`S_z|SM=M|SM`$. Then the magnetic impurity in the presence of the non-magnetic potential and the pairing disorder can all be solved exactly. There are bound states for each of the $`2S+1`$ eigenstates.

In Figs. 1 and 2, we have illustrated how this solution can aid in understanding the STM results obtained from the Bi sites directly above Zn and Ni impurity sites in the presumed top underlying CuO2 plane of BSCCO, respectively.pan00 ; hudson01  In Fig. 1, we fit the Zn STM data with a single pole obtained from Eq. (24) with $`𝐯(v_1,v_2,v_3)=(0,0,\pm 1)`$ \[or equivalently with $`𝐯=(0,\pm 1,0)`$\] and a width $`\delta \omega `$ chosen to fit the data. The experimental peak center appears at a slight offset from zero, which can be understood quantitatively by adjusting $`v_1<<1`$ and a width that merges the two resulting peaks into a single peak broader than the offest.

In Fig. 2, we used Eq. (33) with $`𝐯=(0,\pm 0.51,\pm 1.01)`$ \[or equivalently with $`𝐯=(0,\pm 1.01,\pm 0.51)`$\] to fit the two peak positions, and adjusted the widths to those of the data obtained above a Ni site.hudson01  We did not show the other two peaks that one expects from Eq. (33). That is because the data not pictured here show that the two peaks present in these data appear at equivalent negative biases on adjacent Bi sites,hudson01  for reasons that are not particularly clear. In any event, our theory suggests that the STM data for Zn are consistent with it behaving either as a strong pairing fluctuation defect or as a strong magnetic impurity, and the Ni data suggests it behaves as both a strong magnetic impurity and a strong pairing defect in BSCCO. The similar strengths of the two pairbreaking defects might explain why the $`T_c`$ suppressions obtained in doping BSCCO with these elements are nearly identical.

Thus, we solved the T-matrix in this modified BCS model of a local, on-site attractive pairing interaction with three types of defects on a site.Nandini  For a superconductor with local, near-neighbor pairing of $`d_{x^2y^2}`$-wave symmetry, the local gap $`\mathrm{\Delta }_{ij}`$ at the site $`(i,j)`$ on a tetragonal lattice is coupled heirarchically to the $`\mathrm{\Delta }_{i^{}j^{}}`$ at every site $`(i^{},j^{})`$. Hence, a generalization of our results to $`d`$-wave superconductors would not be straightforward.

In summary, we have shown that disordered short coherence length superconductors can exhibit pairbreaking from spatial fluctuations in the pairing interaction, in a manner very similar to that found with magnetic impurities. We studied the effects of a single site with up to three types of defects using the T-matrix approximation, and found bound states within the superconducting gap arising from either pairing fluctuations or magnetic impurities. Our best fits to the scanning tunneling microscopy data above the sites of Zn and Ni impurities suggest that Zn behaves either as a strong pairing fluctuation defect or as a strong magnetic impurity, whereas Ni behaves as both a strong magnetic impurity and a strong pairing defect.


text:
# Coarsening of Sand Ripples in Mass Transfer Models with Extinction

## I Introduction and Motivation

When a surface of sand is exposed to wind or water flow, patterns like ripples or dunes are commonly formed. The physics of this process is extremely complex because it involves the interaction of a granular medium with a possibly turbulent hydrodynamic flow Bagnold . It is therefore desirable to develop simplified models that capture some of the key features of the pattern formation.

In this paper we are concerned with a class of models which focus on the role of the mass transfer in the evolution of the pattern. Along a one-dimensional cut perpendicular to the ripples, the pattern is described by a set $`\{\lambda _i\}`$ of ripple lengths, where the index $`i`$ labels the ripples in the array. The $`\lambda _i`$ are used here as a general measure of ripple size, without reference to the detailed geometry of individual ripples (see Fig.1). In particular, we do not distinguish between the linear size of a ripple and the mass it contains (for further discussion of this point see Andersen:PRE63 ).

During the evolution of the patterns, the flow transfers mass between neighboring ripples. The central assumption of the model is that the mass transferred to ripple $`i`$ from ripple $`i+1`$ or $`i1`$ (per unit time) is a function $`\mathrm{\Gamma }(\lambda _i)`$ of the size of the ripple which gains the mass. Further motivation for this assumption will be given below. We refer to $`\mathrm{\Gamma }(\lambda )`$ as the *robber function* 111This term was suggested to us by Ko van der Weele..

Depending on the characteristics of the flow, the mass transfer between ripples can be symmetric or asymmetric. In the symmetric case the balance between loss and gain processes for a given ripple leads to the evolution equation Andersen:PRE63

$$\frac{\mathrm{d}\lambda _i}{\mathrm{d}t}=\frac{1}{2}[\mathrm{\Gamma }(\lambda _{i1})+2\mathrm{\Gamma }(\lambda _i)\mathrm{\Gamma }(\lambda _{i+1})],$$
(1)
while in the asymmetric case (assuming, say, mass transfer only to the left) one has

$$\frac{\mathrm{d}\lambda _i}{\mathrm{d}t}=\mathrm{\Gamma }(\lambda _{i1})+\mathrm{\Gamma }(\lambda _i).$$
(2)
The factor $`1/2`$ in Eq. (1) makes the time scales for both the dynamics equal.

A homogeneous state of equally sized ripples, $`\lambda _i\overline{\lambda }`$, is stationary under (1) and (2), but its stability depends on the derivative of the robber function: The pattern is stable for $`\mathrm{\Gamma }^{}(\overline{\lambda })<0`$ and unstable for $`\mathrm{\Gamma }^{}(\overline{\lambda })>0`$ Andersen:PRE63 . In the unstable case the dominant mode is a modulation of period 2, in which every second ripple grows and every second one shrinks. As the size of the shrinking ripples reaches zero in a finite time, the evolution equations (1) and (2) have to be supplemented by an *extinction rule*: When the size of a ripple vanishes, it is removed from the system and the remaining ripples are relabeled such that the previous neighbors of the removed ripple become neighbors of each other. Extinction events contribute to the *coarsening* of the pattern, i.e., to an increase of the mean wavelength. In this work the reverse process of ripple creation is not considered, hence coarsening is irreversible.

The symmetric mass transfer model (1) was first proposed as a description of vortex ripples in coastal waters, which are created under the oscillatory flow of surface waves Andersen:PRE63 . In that context the dependence of the robber function on the size of the gaining ripple is motivated by the observation that the mass transfer is effectuated mostly by a separation vortex which appears in the wake of that ripple. Numerical simulations Andersen:PRE63  and experiments Andersen:preprint  show that $`\mathrm{\Gamma }(\lambda )`$ is nonmonotonic, with a maximum near $`\lambda =a`$, where $`a`$ is the amplitude of the fluid oscillations. Thus patterns of wavelength $`\overline{\lambda }<a`$ ($`\overline{\lambda }>a`$) are unstable (stable), and the main interest is in the wavelength selection process starting from a short wavelength, unstable state Andersen:PRE63 ; Andersen:preprint ; Krug:ACS .

A related, asymmetric mass transfer model for wind driven sand ripples was introduced in Werner:PRL71 . The basic hypothesis of the model is that wind ripples wander with a speed that is inversely proportional to their size. This implies that a leading ripple (ripple $`i+1`$) is eroded by the trailing ripple (ripple $`i`$) at a rate which is proportional to $`1/\lambda _i`$, so the resulting evolution equation is of the type (2) with $`\mathrm{\Gamma }(\lambda )1/\lambda `$. Since $`\mathrm{\Gamma }^{}(\lambda )<0`$, the homogeneous pattern is stable. However, when fluctuations are included by discretizing the ripple sizes and implementing a stochastic mass transfer rule, a fluctuation-driven coarsening mechanism becomes effective and leads to an increase of the mean wavelength with time $`t`$ as $`\mathrm{ln}t`$.

In this article we consider a class of stochastic models whose noiseless counterparts are described by (1) or (2). We concentrate on monotonic, algebraic robber functions $`\mathrm{\Gamma }(\lambda )\lambda ^\gamma `$ and study the coarsening process regarding $`\gamma `$ as a variable parameter. For $`\gamma <0`$ this extends the results of Werner:PRL71  on fluctuation-driven coarsening. The case $`\gamma >0`$ is a simple realization of linearly unstable ripple evolution, and it is studied here as a first step towards a better understanding of models with nonmonotonic robber functions Andersen:PRE63 ; Andersen:preprint ; Krug:ACS . Although the models are defined using the terminology of sand ripples, they are connected to other problems in nonequilibrium statistical physics. For example, for $`\gamma =0`$ the system maps to coalescing random walks and is therefore exactly solvable. Other equivalences include exclusion processes, zero range processes, urn models, and cluster-cluster aggregation.

Our main results are the following. In general, one can identify two time scales in the dynamics: The one of ripple extinctions and the other at which the system would equilibrate to a steady state in the absence of extinctions. For $`\gamma <0`$ the loss of a ripple is a rare fluctuation when the mean ripple size is large. Therefore the two time scales are well separated, and the system has time to relax to a quasi-steady state between ripple extinctions. We show that this state is characterized by a product measure. This justifies the mean field assumption made in Werner:PRL71 , and allows us to calculate the stationary ripple size distribution. The product measure becomes exact only at the limit $`t\mathrm{}`$ as the correlations in the system decay as a power law. The average ripple size grows logarithmically at late times, with a prefactor $`\gamma ^1`$.

For $`\gamma >0`$ extinctions are frequent events which occur on the same time scale as the evolution of the surviving ripples. We find in this case that the noise is irrelevant, so that the dynamics can be described by (1) and (2). For $`0<\gamma <1`$ the mean ripple size grows algebraically with the exponent $`1/(1\gamma )`$, while the growth is exponential for $`\gamma =1`$. In the latter case the evolution equations become linear, and the problem can be solved exactly on the mean field level. The mean field theory reproduces the exponential growth for the mean ripple size, but incorrectly predicts a dependence of the ripple size distribution and the coarsening law on the initial conditions.

In the next section the model is introduced and its relations to other models are discussed. Algebraically decaying robber functions ($`\gamma <0`$) are considered in Sec. III. The product form of the mass distribution is derived in Sec. III.1, the coarsening law is calculated in Sec. III.2, and the approach to the product measure is analyzed in Sec. III.3. Section IV is devoted to algebraically growing robber functions ($`\gamma >0`$). The mean-field theory is first developed for $`\gamma =1`$ and then compared to simulations (Sec. IV.1). Section IV.2 examines the case $`0<\gamma <1`$. Conclusions and open questions are formulated in Section V.

## II The stochastic ripple model

### II.1 Definition and simulation algorithm

In the stochastic model a sand ripple is characterized by its mass $`m`$. The mass variables are integers such that each ripple consists of $`m_i`$ elementary mass units and occupies a site $`i`$ on a one-dimensional lattice. The mass is conserved, i.e., M := $`_im_i=\mathrm{const}.`$. The $`m_i`$ correspond to the length variables $`\lambda _i`$ used in equations (1) and (2). As mentioned in the Introduction, the mass and the length of ripples are here considered to be indistinguishable. We use different symbols for two reasons. We want to make a clear distinction between (i) the real and integer valued ripple sizes and (ii) between the deterministic and noisy dynamics.

Ripples interact only by exchanging mass with their nearest neighbors with an algebraic mass transfer rates $`\mathrm{\Gamma }(m)=\mathrm{\Gamma }_0m^\gamma `$. Since the constant $`\mathrm{\Gamma }_0`$ affects only the time scale it will be set equal to unity from now on. If ripples obtain mass only from one of their neighbors, say, from the right one, the mass transfer is called (totally) asymmetric. If the mass comes from both neighbors we call the dynamics symmetric. As was discussed in Section I, the asymmetric mass transfer naturally arises in the case of wind ripple formation Werner:PRL71  whereas the symmetric dynamics takes place for ripple patterns forming under an oscillatory flow Andersen:PRE63 .

In addition the model includes the removal of ripples when their mass becomes zero. This is done such that lattice sites containing no mass are eliminated from the system. In this way each ripple always has a neighbor from which it can gain mass. If we denote the number of lattice sites at time $`t`$ by $`N(t)`$, the average ripple mass is $`m(t)=M/N(t)`$.

In the simulations three different initial conditions are used. As random initial conditions we denote the case in which the probability to have a ripple of size $`m`$ is given by the geometric distribution $`(1q)q^{m1}`$, with $`0<q<1`$. The probability $`q`$ is related to the mean ripple size as $`m=(1q)^1`$. A distribution $`m_i=mi`$ is referred to as monodisperse. The third possibility is a Poisson distribution.

The dynamics is implemented as follows. First a ripple is selected randomly and time is incremented by $`N(t)^1\mathrm{\Gamma }_{\mathrm{max}}^1`$, where $`\mathrm{\Gamma }_{\mathrm{max}}`$ is the maximum of all the rates of the ripples in the system at time $`t`$. Denote the mass of the selected ripple by $`m`$. If $`x<\mathrm{\Gamma }(m)/\mathrm{\Gamma }_{\mathrm{max}}`$, where $`x`$ is an uniformly distributed random number in the interval $`[0,1]`$, the ripple gets a unit mass from its nearest neighbor. Otherwise a new ripple is selected and the process is repeated. For symmetric dynamics the neighbor is selected randomly whereas in the asymmetric case it is always the right one.

### II.2 Relation to other models

The model defined above is inspired by the worm model originally introduced to describe the coarsening of wind ripples Werner:PRL71 . Here it is generalized in two respects. First, the mass transfer rate in the worm model is inversely proportional to the ripple mass whereas in the generalized model it is given by $`\mathrm{\Gamma }(m)m^\gamma `$. Second, we consider also the symmetric mass transfer between neighbors. In the case of asymmetric dynamics and for $`\gamma =1`$ our model reduces to the worm model.

Apart from the extinction step, the sand ripple model is similar to a zero range process Spitzer:AIM5 ; Spohninkirja ; Evans:BJP . Both models are defined in terms of conserved, integer mass variables $`m_i`$ which interact through the (symmetric or asymmetric) exchange of unit masses between nearest neighbor sites of a lattice. The key difference is that in a zero range process the mass transfer rate is a function of the mass at the site of *departure*, while in the sand ripple model it depends on the mass at the *target* site. This reverses the sign of the right hand sides of (1) and (2), and hence the stability properties of the model: In a zero range process the homogeneous state is stable if $`\mathrm{\Gamma }^{}(\overline{m})>0`$ and unstable if $`\mathrm{\Gamma }^{}(\overline{m})<0`$. The coarsening behavior in zero range processes with nonmonotonic robber function is relevant to clustering in granular gases vanderMeer:PRE63 .

The occurrence of irreversible extinction events in our model is reminiscent of certain urn models that have been proposed in the context of glassy dynamics Godreche: . For example, consider the backgammon model Ritort:PRL75  which is defined by $`M`$ particles distributed among $`N`$ boxes with $`m_i`$ particles in the $`i`$th box. The Hamiltonian is $`H=_i\delta _{m_i,0}`$ so that the energy corresponds to minus the number of empty boxes. The $`N`$-fold degenerate ground state therefore consists of a condensate, where all the particles belong to one box. Associating the masses of ripples with the particle numbers and lattice sites with boxes, the ripple evolution becomes similar to the backgammon model at zero temperature, where the empty boxes are not refilled once they have become empty (ripple extinction).

In contrast to the ripple model, the urn models have no spatial structure, i.e., mass transfer is possible between any pair of boxes. In the standard dynamical scheme, originally due to Ehrenfest, in each time step one of the balls is chosen at random and a move to another box is attempted Godreche: . The probability for a box to be chosen is then proportional to its occupation number. In our setting this corresponds to a mass transfer rate $`\mathrm{\Gamma }(m_i)m_i`$, where $`i`$ is the site of *departure*; in this respect the urn models are related to zero range processes. Since $`\mathrm{\Gamma }^{}>0`$, the homogeneous state is linearly stable and coarsening (i.e., evolution towards the ground state) is very slow.

The sand ripple model can also be mapped to an exclusion process Spohninkirja ; Liggetinkirja ; DerridaEvans\_in\_Privmans\_book . The mapping can be done in two ways which differ in how the disappearance of ripples is taken into account. The mappings proceed along the lines of Majumdar:JSP99  and the first one is schematically presented in Fig. 2. One constructs a new lattice with $`L(t)=M+N(t)`$ sites. The mass variables $`m_i`$ of the ripples turn to $`m_i`$ consecutive holes separated by particles on the new lattice. More precisely, there exist particles on sites $`i+_{k=1}^im_k`$ ($`i=1,\mathrm{},N(t)`$) while the rest are empty. Moving one mass unit from one ripple to its neighbor corresponds to a hop of a particle in the exclusion process. Naturally the exclusion process is either symmetric or asymmetric as the original dynamics. For the exclusion process the loss of a ripple becomes a coalescence of particles at contact, which changes the length $`L(t)`$ of the system.

As the masses map to holes, the hopping rates of particles depend on the distances between them. For negative (positive) $`\gamma `$ there is a repulsive (attractive) interaction between the particles. For the marginal case of mass independent transfer rates ($`\gamma =0`$) the interaction vanishes and the exclusion process reduces to coalescing random walkers. This is a well-known problem which can be solved exactly e.g. with the method of interparticle distribution functions ben-Avraham\_in\_Privmans\_book . The most relevant results for our case are: (i) the average ripple mass grows asymptotically as $`m(t)\sqrt{t}`$ and (ii) the ripple size distribution (the probability of finding a ripple of mass $`m`$ at time $`t`$) can be written in a scaling form as $`p(m;t)=m^1G(m/m)`$, where the scaling function $`G(x)=\frac{1}{2}\pi x^2e^{\pi x^2/4}`$.

A length conserving mapping is obtained by considering the model with ripple extinction but without removal of empty lattice sites. Particles hop over these sites in order not to allow for creation of new ripples. Again the empty sites become particles in the exclusion picture but now the mass transfer corresponds to a hop of a particle over all the particles in the same cluster. This is illustrated in Fig. 3.

A long one particle hop can be considered as moving a cluster of particles as a whole. The loss of a ripple turns into the aggregation of clusters. In this way our model further maps to a cluster–cluster aggregation process where each cluster moves with a rate that depends on the distance to its neighbor(s). As the ripples map to holes, the main interest lies in the size distribution of distances between clusters and not in the cluster size distribution itself, as is usually the case. One-dimensional cluster-cluster aggregation models generally obey universal dynamical scaling (see Hellen:PRE62  and references therein), which will be seen to be the case also for the sand ripple model.

## III Noise-induced coarsening ($`\gamma <0`$)

It is known from the mean-field analysis of Werner:PRL71  that for $`\gamma =1`$ the average ripple size grows as $`m(t)\mathrm{ln}(t+e^{m(0)})`$. Intuitively, the slow growth follows from the fact that, for $`\gamma <0`$, the ripples near extinction are those with the highest incoming mass rates. Therefore the disappearance of a ripple involves a rare fluctuation; within the approach of Werner:PRL71 , the mass of a ripple evolving in a background of mean mass $`m(t)`$ performs a random walk that is biased away from zero.

We base our theoretical analysis on this observation. In what follows, we will assume that at long times, i.e., for large $`m`$, the extinction of a ripple is such a rare event, that it does not affect the ripple size distribution. Neglecting extinctions, we show that the steady state distribution is given by a product form. After validating the quasi-static approximation by simulations, we use it to show that to leading order in $`t`$ the average ripple mass grows as $`m\gamma ^1\mathrm{ln}(t)`$. Finally we consider the approach to the product measure by studying nearest neighbor correlations.

### III.1 Mass distribution

Without extinctions the ripple size distribution can be solved exactly. This is due to the short range of interactions between ripples: The mass transfer rate depends only on the mass at the target site. As was noted above in Sec. II.2, this is similar to the zero-range process, where the rate depends only on the site of departure. The most important characteristics of a zero range process is that its steady state is described by a product measure Evans:BJP . This was shown to generalize to processes where the transition rate is a product of functions of the occupation numbers at the site of departure and the target site Karimipour: . As our model is a special case of this class of models, the results of Karimipour:  apply here as well. For completeness we give a brief derivation.

The product measure property implies that the stationary probability distribution, $`P(\{m_i\})`$, of finding the system in configuration $`\{m_1,m_2,\mathrm{},m_N\}`$ factorizes as

$$P(\{m_i\})=\underset{i}{}p(m_i),$$
(3)
where $`p(m_i)`$ is the probability of finding mass $`m_i`$ at site $`i`$. In the steady state there are no correlations between the ripple sizes. Starting from the master equation for $`P(\{m_i\})`$ and using the product form (3), one obtains for the asymmetric case (the calculation is not presented here since up to index changes it is identical to that presented in Evans:BJP )

$$p(m_i)p(m_{i+1})\mathrm{\Gamma }(m_i)=p(m_i+1)p(m_{i+1}1)\mathrm{\Gamma }(m_{i+1}1).$$
(4)
The condition (4), known as pairwise balance Schutz:JPHYSA , generalizes the detailed balance condition familiar from equilibrium statistical mechanics. It has a simple interpretation. The left hand size of Eq. (4) represents the mass transfer to the site $`i`$ which has to balanced by a transfer out of this site (the right hand side) in order to be in the steady state. The first two terms give the probability to find a mass $`m_i`$ at site $`i`$ with a right neighbor with mass $`m_{i+1}`$ and the last term describes the rate at which the site $`i`$ gains mass from its neighbor. We emphasize that, provided a solution to (4) can be found, this proves that the product measure (3) is an *exact* stationary solution of the master equation; on the basis of general arguments, this solution is then also expected to be unique.

Proceeding similarly for the symmetric dynamics gives (transitions $`\{\mathrm{},m_{i1},m_i,m_{i+1},\mathrm{}\}\{\mathrm{},m_{i1},m_i+1,m_{i+1}1,\mathrm{}\}`$ and $`\{\mathrm{},m_{i1}1,m_i+1,m_{i+1},\mathrm{}\}\{\mathrm{},m_{i1},m_i,m_{i+1},\mathrm{}\}`$)

$`p(m_{i1})p(m_i)p(m_{i+1})\mathrm{\Gamma }(m_i)=`$
$`p(m_{i1}1)p(m_i+1)p(m_{i+1})\mathrm{\Gamma }(m_{i1}1).`$ (5)
Since $`p(m_{i+1})`$ cancels out we end up with Eq. (4). Therefore the steady state distribution is independent of the asymmetry of the dynamics.

Equation (4) can be recast as

$$\alpha ^1:=\frac{p(m_i)\mathrm{\Gamma }(m_i)}{p(m_i+1)}=\frac{p(m_{i+1}1)\mathrm{\Gamma }(m_{i+1}1)}{p(m_{i+1})},$$
(6)
where $`\alpha `$ must be a constant. Denoting $`p(0)=p_0`$ and recursively iterating equation (6), we obtain

$$p(m)=p_0\alpha ^m\underset{i=1}{\overset{m1}{}}\mathrm{\Gamma }(i)=p_0\alpha ^m[(m1)!]^\gamma ,$$
(7)
where the product for $`m=1`$ is defined to give unity and the last form follows from the definition $`\mathrm{\Gamma }(m)=m^\gamma `$.

The unknown constants $`p_0`$ and $`\alpha `$ can be determined by the normalization $`_{m=0}^{\mathrm{}}p(m)=1`$ and the expectation value $`m:=_{m=0}^{\mathrm{}}mp(m)`$. Explicit results for $`\gamma =1`$ and $`2`$ can be found in Appendix A; for $`\gamma =1`$, (7) is a (shifted) Poisson distribution. In general, the distribution for $`m1`$ can be written as

$$p(m)=C_2(\gamma )e^{\gamma m}m^{\gamma m(1\gamma )/2}[(m1)!]^\gamma ,$$
(8)
where the explicit form of $`C_2(\gamma )`$ is not important for our purposes. Using the form given in Eq. (8), it is easy to show that the width $`\sigma :=\sqrt{m^2m^2}`$ of the distribution behaves as $`\sigma \sqrt{m}`$ independent of $`\gamma `$.

The calculated distributions are compared to numerics in Figs. 4 and 5. The average ripple mass is not a constant as the simulations include also ripple extinction. The excellent agreement at long times shows that indeed these become so rare, that between subsequent extinctions the system has time to equilibrate to the steady state. Note that all initial distributions converge to the universal distribution $`p(m;t)`$ given by Eq. (8), where the time-dependence enters only through the mean ripple mass $`m(t)`$.

### III.2 Coarsening law

Next we proceed to calculate the mean ripple size $`m(t)`$ using an approach similar to the analysis of the backgammon model Godreche: . We assume that, at long times, the probability for a given ripple to vanish is equal to the probability $`p(0)`$ obtained by extrapolating the steady state probability distribution (8) to $`m=0`$. The number $`N`$ of ripples then decays according to $`dN/dtp(0)N`$. Since $`m(t)=M/N`$ we obtain

$$\frac{\mathrm{d}m(t)}{\mathrm{d}t}p(0)m(t)e^{\gamma m}m^{(1\gamma )/2},$$
(9)
which to leading order in $`t`$ gives

$$m(t)\gamma ^1\mathrm{ln}(t).$$
(10)
Simulations with different initial conditions are in accord with Eq. (10) (Fig. 6).

### III.3 Decay of correlations

The product measure for the ripple size distribution implies that there are no correlations between neighboring ripples. This is true only asymptotically. To study the approach to the product measure distribution we consider the normalized nearest neighbor time correlation function

$$g(t):=\frac{m_im_{i+1}m^2}{m^2}.$$
(11)
As is clear from Fig. 7, the early time behavior is sensitive to the details of the initial distribution. In this regime it is possible to have positive correlation between neighboring ripples but at long times there will always be anticorrelations, i.e., $`g(t)<0`$. The numerically observed correlations seem to be independent of the initial conditions and vanish in a universal manner as

$$g(t)t^{1/2}$$
(12)
for both symmetric and asymmetric dynamics.

At first sight one may be tempted to relate the decay of correlations to the extinction events, which perturb the product measure. However, as was shown in Sec. III.2, the probability of extinction events decays as $`p(0)e^{\gamma m(t)}t^1,`$ which is much faster than the numerically observed decay law (12). This implies that the power law (12) is associated with the dynamics between extinction events, which can be described using standard hydrodynamic fluctuation theory for a one-dimensional system with a single conserved density.

Let $`\varphi (x,t)`$ denote the coarse grained mass fluctuations in the (quasi-) steady state of mean mass $`m`$. The long wavelength behavior of $`\varphi `$ is governed by a Langevin equation of the generic form vanBeijeren:PRL54

$$\frac{\varphi }{t}=\nu \frac{^2\varphi }{x^2}\mu \varphi \frac{\varphi }{x}\frac{\eta }{x},$$
(13)
where $`\eta (x,t)`$ is Gaussian white noise with covariance $`\eta (x,t)\eta (x^{},t^{})=D\delta (xx^{})\delta (tt^{})`$. For $`\mu 0`$, Eq. (13) is the noisy Burgers equation Forster:PRA , which has been widely studied in the context of driven diffusive systems vanBeijeren:PRL54  and interface growth Kardar:PRL56 ; KrugSpohn\_in\_Godreches\_book ; Krug:ADVIP .

The coefficients $`\nu `$, $`\mu `$, and $`D`$ appearing in the long wavelength description can be related to the microscopic dynamics of the sand ripple model as follows. The nonlinear term on the right hand side of Eq. (13) is generated by the asymmetry, and its coefficient is given by $`\mu =j^{\prime \prime }(m)`$, where $`j`$ is the steady state mass current. Since in our case $`j=\mathrm{\Gamma }`$, we conclude that $`\mu m^{\gamma 2}`$. In the symmetric case $`\mu =0`$ and the diffusion coefficient $`\nu `$ is proportional to $`\mathrm{\Gamma }^{}m^{\gamma 1}`$ (this can be seen by expanding Eq. (1) around the homogeneous state). Finally, owing to a fluctuation-dissipation theorem Forster:PRA , the equal time correlations of (13) are Gaussian with covariance $`\varphi (x)\varphi (x^{})(D/\nu )\delta (xx^{})`$ independent of $`\mu `$. As we have shown above in Sec. III.1, in the ripple model the variance of the mass fluctuations is always of order $`m`$, hence $`D/\nu m`$.

We want to use equation (13) to describe the approach to the steady state, starting from some initial condition (e.g., the monodisperse state $`\varphi 0`$) specified at $`t=0`$. The analysis of Eq. (13) shows that at long times, and for $`xx^{}`$, the pair correlation function takes the scaling form KrugSpohn\_in\_Godreches\_book ; Krug:ADVIP

$$\varphi (x,t)\varphi (x^{},t)=\frac{D}{\nu }\frac{1}{\xi (t)}𝒢(\xi (t)^1|xx^{}|).$$
(14)
Here $`𝒢`$ is a scaling function, and $`\xi (t)`$ denotes the dynamic correlation length. The prefactor of the scaling function on the right hand side of Eq. (14) is fixed by the requirements that (i) the steady state density fluctuations are proportional to $`D/\nu `$, and (ii) the integral over the pair correlation function is constant due to mass conservation. The correlation length grows diffusively as $`\xi (t)(\nu t)^{1/2}`$ for $`\mu =0`$ and superdiffusively as $`\xi (t)[(D/\nu )^{1/2}\mu t]^{2/3}`$ for $`\mu 0`$.

Keeping $`|xx^{}|`$ fixed and taking $`t\mathrm{}`$, we see that the pair correlations (14) decay as $`(D/\nu )𝒢(0)\xi (t)^1`$. Expressing $`\nu `$ and $`D/\nu `$ in terms of the mean ripple mass, we conclude that in the symmetric case ($`\mu =0`$) the normalized correlation function (12) should decay as

$$g(t)\frac{m^{(1+\gamma )/2}}{t^{1/2}}\frac{(\mathrm{ln}t)^{(1+\gamma )/2}}{t^{1/2}}.$$
(15)
For $`\gamma =1`$ the logarithmic factor disappears and (15) becomes a pure power law with exponent $`1/2`$, in accordance with the simulation results shown in Fig. 7. Moreover, the explicit calculation for the diffusive case shows that the scaling function $`𝒢`$ in Eq. (14) is negative, hence $`g(t)<0`$ as is observed numerically.

In the asymmetric case the fluctuation theory predicts an asymptotic decay as $`g(t)1/\xi (t)t^{2/3}`$, with logarithmic corrections due to the growth of $`m(t)`$. Then why do we find $`g(t)t^{1/2}`$ also for asymmetric dynamics? To answer this question, we recall that the asymptotic, superdiffusive behavior predicted by the noisy Burgers equation sets in only beyond a crossover time scale $`t_\times `$, which increases rapidly with decreasing strength $`\mu `$ of the nonlinear term Krug:ADVIP . The crossover time is of the order of $`t_\times \nu ^5/(D^2\mu ^4)`$. Inserting the estimates for $`\mu `$, $`\nu `$, and $`D/\nu `$ derived above, we see that for the ripple model

$$t_\times m^{3\gamma }(\mathrm{ln}t)^{3\gamma }.$$
(16)
For the case $`\gamma =1`$ considered in Fig. 7, this implies that superdiffusive behavior can be expected only for times such that $`t/(\mathrm{ln}t)^41`$. The left hand side of this inequality becomes equal to unity for $`t5500`$ and reaches the value 10 only for $`t235000`$. Thus the asymptotic regime has not been reached in our simulations. The slight deviation of the simulation data from the $`t^{1/2}`$-behavior seen after $`t=10^5`$ may indicate the beginning of the crossover.

## IV Unstable coarsening ($`\gamma >0`$)

For $`\gamma >0`$ the homogeneous state is linearly unstable because the largest ripples are those with the highest growth rate. Ripple extinction is then no longer a rare event, and the product measure solution derived in Sec. III.1 becomes invalid. On the other hand, it is plausible (and will be confirmed by simulations, see below) that the linear instability supersedes the noise in the time evolution, so that the deterministic equations (1) and (2) and the stochastic ripple model show the same behavior.

In what follows, we first develop a mean-field theory for the deterministic model in the simplest case of a linear robber function ($`\gamma =1`$). Simulations show that the mean field theory is not quantitatively correct, presumably due to the neglect of spatial fluctuations. In the nonlinear regime $`0<\gamma <1`$ we use scaling analysis to derive the coarsening law.

### IV.1 Mean-field analysis for $`\gamma =1`$

We start our analysis from the deterministic equations (1) and (2). For $`\gamma =1`$ these become linear but the system is still non-trivial due to the ripple extinction. As the system is deterministic, the only randomness lies in the initial condition. We denote the initial ripple size distribution by $`P_0(\lambda _0)`$ and its average by $`\overline{\lambda }_0`$.

The mean field approximation consists of replacing the ripples surrounding an arbitrary ripple of size $`\lambda `$ by ripples of the average size $`\lambda `$, such that the evolution equation becomes

$$\frac{\mathrm{d}\lambda }{\mathrm{d}t}=\mathrm{\Gamma }(\lambda )\mathrm{\Gamma }(\lambda )=\lambda \lambda .$$
(17)
On this level there is no difference between symmetric and asymmetric mass transfer. The solution of Eq. (17) reads $`\lambda (\lambda _0,t)=e^t[\lambda _0F(t)]`$, where the function

$$F(t):=_0^td\tau e^\tau \lambda (\tau )$$
(18)
has to be calculated self-consistently. Note that at this point we do not explicitly restrict $`\lambda (t)`$ to be nonnegative (this constraint will enter later). Once $`F(t)`$ is known, the ripple size distribution at time $`t`$ can be obtained by inverting the solution for $`\lambda (\lambda _0)`$ and inserting this into the initial distribution, with the result

$$p(\lambda ;t)=e^tP_0(e^t\lambda +F(t)).$$
(19)
Thus in the mean-field approximation the ripple size distribution preserves its initial shape but gets scaled and shifted.

It is possible to derive a differential equation for the unknown function $`F(t)`$. The fraction $`\rho (t)`$ of surviving ripples is equal to the probability that $`\lambda (t)>0`$,

$$\rho (t)=_0^{\mathrm{}}\mathrm{d}\lambda p(\lambda ;t)=_{F(t)}^{\mathrm{}}\mathrm{d}xP_0(x)=:P_0^c(F(t)),$$
(20)
where the last equation defines the cumulative distribution $`P_0^c`$. The average ripple size is given by $`\lambda (t)=\overline{\lambda }_0/\rho (t)`$. Inserting this into the definition of $`F(t)`$ and differentiating once gives

$$\frac{\mathrm{d}F(t)}{\mathrm{d}t}=\frac{e^t\overline{\lambda }_0}{P_0^c(F(t))}.$$
(21)
Hence the problem reduces to solving the differential equation (21) for a given initial distribution $`P_0(\lambda )`$.

For example, for an exponential initial distribution $`P_0(\lambda _0)=\overline{\lambda }_0^1e^{\lambda _0/\overline{\lambda }_0}`$, we find $`F(t)=\overline{\lambda }_0t`$ and $`\lambda (t)=\overline{\lambda }_0e^t`$, whereas for a flat distribution

$$P_0(\lambda _0)=\{\begin{array}{cc}(2\overline{\lambda }_0)^1\hfill & \lambda _02\overline{\lambda }_0\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(22)
the solution is given by $`F(t)=2\overline{\lambda }_0(1e^{t/2})`$ and $`\lambda (t)=\overline{\lambda }_0e^{t/2}`$. As the rate of exponential growth is different in these two cases, we conclude that the coarsening behavior of the mean field model (17) is nonuniversal.

In general, the exponential growth rate of the mean ripple size is governed by the extremal statistics of the initial distribution $`P_0`$. If the initial ripple sizes are bounded by a maximal size $`\lambda _{\mathrm{m}ax}`$, and $`P_0(\lambda _0)(\lambda _{\mathrm{max}}\lambda _0)^a`$ for $`\lambda _0\lambda _{\mathrm{max}}`$, then the analysis of Eq. (21) shows that $`t^1\mathrm{ln}\lambda (t)(a+1)/(a+2)`$, while for fat initial distributions with a power law tail, $`P_0(\lambda _0)\lambda _0^{(b+1)}`$, we find $`t^1\mathrm{ln}\lambda (t)b/(b1)`$.

To compare the predictions of the mean-field theory to simulations we prefer to show the complement of the cumulative distribution

$$I(\lambda ;t):=_\lambda ^{\mathrm{}}\mathrm{d}xp(x;t)=:f\left(\frac{\lambda }{\lambda (t)}\right),$$
(23)
where the last equation defines the scaling function $`f(x)`$. A similar definition applies to $`p(m;t)`$ with the integral replaced by a sum. In the case of an exponential distribution $`p(\lambda ;t)`$ also the function $`f(x)`$ is exponential whereas for a flat $`p(\lambda ;t)`$ it is linear.

We solved the deterministic equations (1) and (17) using the fourth-order Runge-Kutta method NumericalRecipes . As a check of the algorithm, we reproduced the solution (19) of the mean-field equations. For the full noiseless system (1) the exponential initial distribution remains unchanged \[Fig. 8 a); dashed lines\] but also a flat initial distribution presumably approaches the exponential one \[Fig. 8 b); dashed lines\]. This is in conflict with the mean-field prediction. In both cases the average ripple size grows as $`m(t)e^t`$.

Similarly, in the discrete, noisy ripple model the random initial distribution quickly converges towards an exponential scaling function \[Fig. 8 a); solid lines\]. The monodisperse initial condition spreads out and approaches the same form \[Fig. 8 b); solid lines\]. Again, the mean ripple size grows as $`m(t)e^t`$ for both initial distributions.

Since the deterministic model behaves in a similar manner as the noisy one, we conclude that, in contrast to the case $`\gamma <0`$, the noise is irrelevant. The discrepancy between the mean-field theory and the full deterministic system suggests that the spatial fluctuations are important, as is often the case for low dimensional systems. In particular, the numerical results indicate that, in contrast to the mean field prediction, the behavior of the full system is universal with respect to the initial ripple size distribution.

### IV.2 Coarsening law for $`0<\gamma <1`$

As the mean-field equation is not readily solvable for $`\gamma 1`$ and probably would not describe the problem correctly anyway, here we present a simple scaling argument for the growth of the mean ripple size in the regime $`0<\gamma <1`$. We start from the observation that in the linearly unstable case ($`\mathrm{\Gamma }^{}(\lambda )>0`$) predominantly every second ripple grows and every second one shrinks Krug:ACS . Therefore we may consider a simplified system consisting of two ripples of initial sizes $`\lambda _1^0>\lambda _2^0`$. We calculate the time $`t^{}`$ at which the average size has doubled. It is given by the conditions $`\lambda _1(t^{})=\lambda _1^0+\lambda _2^0`$ and $`\lambda _2(t^{})=0`$. Since the mass is conserved we have $`\overline{\lambda }:=\lambda _1(t)+\lambda _2(t)=\mathrm{const}.`$.

Applying Eq. (1) gives

$$\{\begin{array}{cc}\dot{\lambda }_1\hfill & =\lambda _1^\gamma \lambda _2^\gamma \hfill \\ \dot{\lambda }_2\hfill & =\lambda _2^\gamma \lambda _1^\gamma ,\hfill \end{array}$$
(24)
where the dot denotes derivative with respect to time. The solution is implicitly given by $`\lambda _2(t)=\overline{\lambda }\lambda _1(t)`$ and

$$t=_{\lambda _1^0}^{\lambda _1(t)}\frac{dx}{x^\gamma (\overline{\lambda }x)^\gamma },$$
(25)
which together with the definition of $`t^{}`$ implies the homogeneity relation

$$t^{}(a\lambda _1^0,a\overline{\lambda })=a^{1\gamma }t^{}(\lambda _1^0,\overline{\lambda }).$$
(26)
Assuming that the evolving ripple size distribution is governed by a single size scale, it follows that the doubling time depends on the mean ripple size as $`t^{}\lambda ^{1\gamma }`$. The inverse of the doubling time is the growth rate of $`\lambda `$. Hence we may write

$$\frac{\mathrm{d}\lambda (t)}{\mathrm{d}t}\frac{1}{t^{}}\lambda (t)$$
(27)
which yields $`\lambda (t)t^z`$ with $`z=1/(1\gamma )`$. This is confirmed by simulations, which give $`z=1.32\pm 0.02`$ and $`1.98\pm 0.03`$ for $`\gamma =0.25`$ and $`0.50`$, respectively. We also numerically checked the universal scaling behavior of the ripple size distribution for $`0<\gamma <1`$, but in this region the scaling function is more complicated than a simple exponential.

## V Conclusions

In this paper we have studied a one-dimensional model for sand ripple evolution, where mass is transferred between neighboring sites with algebraic rates and sites containing no mass are removed from the system (ripple extinction). As the rates depend only on the site to which the mass is transferred, the system is similar to a zero range process. Thus the steady state in the absence of ripple extinction is characterized by a product measure. Asymptotically this continues to hold for algebraically decaying mass transfer rates ($`\gamma <0`$) since the extinctions are exponentially rare at late times. As a consequence the average ripple size grows to leading order logarithmically slowly with the prefactor $`1/\gamma `$.

For $`\gamma <0`$ the approach to the steady state product measure is algebraic. The correlations between masses of neighboring ripples decay universally, i.e., independent of the initial distribution, as $`t^{1/2}`$ and $`t^{2/3}`$ for symmetric and asymmetric mass transfer, respectively. In the asymmetric case the asymptotic regime is preceded by a long crossover where $`t^{1/2}`$-decay is observed.

For algebraically growing robber functions ($`\gamma >0`$) the coarsening is driven by the linear instability of the homogeneous state. Ripple extinctions become frequent, and the product measure is no longer relevant. The average ripple size grows algebraically as $`t^{1/(1\gamma )}`$ for $`0<\gamma <1`$. The behavior at $`\gamma =0`$ is discontinuous since $`m(t)t^{1/2}`$ for $`\gamma =0`$, which follows from the mapping to coalescing random walkers. For $`\gamma =1`$ $`m(t)e^t`$ and the scaling function of the ripple size distribution appears to be a simple exponential. The dynamical noise, which is necessary to have coarsening for $`\gamma 0`$, is irrelevant for $`\gamma >0`$. The mean-field theory developed for $`\gamma =1`$ reproduces the exponential growth of the mean ripple size, but it is insufficient to describe the universality of the growth law and the ripple size distribution which is observed numerically.

It is interesting to compare the results to the behavior in one-dimensional cluster-cluster aggregation. Recall that the model treated here can be mapped to cluster-cluster aggregation with hopping rates of clusters depending algebraically on the distance between them (Sec. II.2). When the hopping rates depend as $`\mathrm{\Gamma }(m)m^\gamma `$ on the masses of clusters, the growth of the average cluster size is algebraic with $`m(t)t^{1/(2\gamma )}`$ for all $`\gamma <2`$ Hellen:PRE62 . Thus the behavior for nonnegative values of $`\gamma `$ is rather similar in the two models, but for $`\gamma <0`$ one finds a drastic difference due to the repulsive interaction between clusters in the ripple model.

We conclude by adducing some open problems for future studies. One of the most interesting issues is to understand the coarsening and the final ripple size selection in the case of a nonmonotonic robber function. This has direct applications in the coarsening of vortex ripples, where the robber function has recently been measured Andersen:preprint . Initially these systems are in the unstable regime, where the transfer function is monotonically increasing. As we have seen in the present article, even this is a harder problem than the case where the dominant contribution to coarsening comes from the dynamical fluctuations. For nonmonotonic robber functions one needs an understanding of both coarsening mechanisms. Therefore the starting point into this direction would be to better understand the $`\gamma >0`$ case.

## Acknowledgements

E.K.O.H. thanks M.J. Alava for valuable discussions and a critical reading of a preliminary version of the manuscript. He is also indebted to the kind hospitality of the University of Essen where the main part of this research was performed. J.K. is grateful to K.H. Andersen for introducing him to the problem, and to M.R. Evans for mentioning the backgammon model. The financial support by Vilho, Yrjö, and Kalle Väisälä Foundation, DFG within SFB 237, and the Academy of Finland’s Center of Excellence program are gratefully acknowledged.

*

## Appendix A Calculation of mass distributions for $`\gamma =1`$ and $`2`$

Here we calculate the explicit form of the mass distribution

$$p(m)=p_0\alpha ^m\underset{i=1}{\overset{m1}{}}\mathrm{\Gamma }(i)=p_0\alpha ^m[(m1)!]^\gamma ,$$
(28)
for $`\mathrm{\Gamma }(m)=m^\gamma `$ in the cases $`\gamma =1`$ and $`2`$. For $`m=1`$ the product in Eq. (28) for $`m=1`$ is defined to be unity, and we set $`(1)!=1`$. The normalization condition $`_{m=0}^{\mathrm{}}p(m)=1`$ gives

$$p(m)=\{\begin{array}{cc}(1+\alpha e^\alpha )^1\alpha ^m/(m1)!\hfill & \text{for}\gamma =1\hfill \\ I_0^1(2\sqrt{\alpha })\alpha ^{m1}/[(m1)!]^2\hfill & \text{for}\gamma =2,\hfill \end{array}$$
(29)
where $`I_n(x)`$ is the modified Bessel function of the first kind. The parameter $`\alpha `$ is related to the expectation value $`m=_{m=0}^{\mathrm{}}mp(m)`$ by

$$m=\{\begin{array}{cc}\alpha (\alpha +1)e^\alpha /\left(1+\alpha e^\alpha \right)\hfill & \text{for}\gamma =1\hfill \\ 1+\sqrt{\alpha }I_1(2\sqrt{\alpha })/I_0(2\sqrt{\alpha })\hfill & \text{for}\gamma =2.\hfill \end{array}$$
(30)
Using the expansions $`I_n(x)=e^x/\sqrt{2\pi x}+𝒪(1/x)`$ for $`\alpha \mathrm{}`$, these formulae simplify to $`m\alpha `$ and $`m\sqrt{\alpha }`$, for $`\gamma =1`$ and $`2`$, respectively. Hence for $`m\mathrm{}`$ the distributions become

$$p(m)\{\begin{array}{cc}e^mm^{m1}/(m1)!\hfill & \text{for}\gamma =1\hfill \\ 2\sqrt{\pi }e^{2m}m^{2m3/2}/[(m1)!]^2\hfill & \text{for}\gamma =2,\hfill \end{array}$$
(31)
which are of the general form indicated in Eq. (8).