text:
# Non-Observation of the Production of Superheavy Nuclei in the Reaction of 86Kr with 208Pb

## I Retraction

In 1999, our group reported  the synthesis of element 118 in the 208Pb(86Kr,n) reaction. Prompted by the absence of similar decay chains in subsequent experiments  the original data files from our experimentshave been re-analyzed. The reported chains are not in the original data files.

Therefore, the previously reported results on the production and decay of element 118 are not valid.


format: text

text:
# Information entropy in fragmenting systems

## Abstract

The possibility of facing critical phenomena in nuclear fragmentation is a topic of great interest. Different observables have been proposed to identify such a behavior , in particular, some related to the use of information entropy as a possible signal of critical behavior . In this work we critically examine some of the most widespread used ones comparing its performance in bond percolation and in the analysis of fragmenting Lennard Jones Drops.

When nuclei collide at energies of the order of hundreds of MeVs, the resulting highly excited system breaks up into several medium-size fragments. This phenomenon is usually referred to as multifragmentation and attracts a lot of attention in the nuclear community because it is conjectured that the fragmenting system might undergo a phase transition. In particular the possibility of facing critical behavior (originally triggered by the seminal work of the Purdue Group ) has induced the development of different signatures to characterize such a behavior from the analysis of the observed mass spectra. Starting with the work of Bauer  and Campi  the percolation model has emerged as a quite useful test-bed for the analysis of the tools developed to study nuclear multifragmentation.

In  the information entropy $`(S_1)`$ has been proposed as a signature of critical behavior and two systems has been analyzed: a lattice gas model and a molecular dynamics model.

In this work we test the information entropy together with other proposed signatures of critical behavior in a system (percolation on a simple cubic lattice) that it is known to undergo a second order phase transition. For this model, the behavior near the critical point (including the critical exponents) is well known. Besides the information entropy, we also analyze the normalized variance of the size of the maximum fragment ($`NVM`$) , the second moment of the distribution ($`M_2`$) , the best fitting power law  and minimum $`\tau `$ .

We have calculated the above mentioned signatures from the mass distributions resulting from a bond percolation process on two simple cubic lattices of sizes $`30\times 30\times 30`$ and $`6\times 6\times 6`$. It is known that for the larger one, the finite size effects are negligible, meanwhile for the smaller one (whose size is of the order of the one corresponding to nuclear system’s size) the finite size effects are important .

Once the various criticality signatures have been tested we use them to analyze excited Lennard-Jones drops undergoing fragmentation. We focus on whether the non-equilibrium process of fragmentation can be characterized as a critical phenomenon.

In percolation, the number of clusters of size $`A`$ per lattice site, $`n_A`$, at a given bond activation probability, $`p`$, was proposed in  to be given by :

$$n_A(p)=q_0A^\tau f(z),$$
(1)
where $`q_0`$ is a normalization factor, $`\tau `$ a critical exponent and $`f(z)`$ a scaling function. This equation shows that $`n_A(p)/q_0A^\tau `$ depends on $`p`$ and $`A`$ via the combination $`z=((pp_c)/p_c)A^\sigma =ϵA^\sigma `$. Thus it is clear that at the critical point $`n_A(p_c)=q_0A^\tau f(0)`$ with $`f(0)=1`$. All this is valid for not too small fragments and infinite systems for which the percolating cluster has been removed. In finite systems we won’t get a pure power law because finite size effects introduce distortions. Nevertheless, there exists a range of masses (depending of the size of the system) for which the mass spectra presents a power law like behavior, as can be seen in figure (1f).

There are two main methods to find the critical point from the mass spectra. The first one consists in fitting the mass spectra with a simple power law, with two independent parameters, $`q_0`$ and $`\tau `$, and then look for the value of $`p`$ (or the quantity taken as control parameter) for which the exponent $`\tau `$ is minimum .

The second one was developed in  and takes care of the constraint imposed by the normalization, i.e. $`q_0=q_0(\tau )`$. This dependence can be found by analyzing the first moment of the cluster distribution at the critical point :

$$M_1(ϵ=0)=\underset{A}{}n_A(ϵ=0)A=q_0\underset{A}{}A^{1\tau }=1,$$
(2)
with the sum running over all clusters. Equation (2) shows that the normalization constant, $`q_0`$, depends on $`\tau `$ via a Riemann $`\zeta `$ function:

$$q_0=\frac{1}{_AA^{1\tau }}.$$
(3)
Therefore, if we take into account the dependence of $`q_0`$ with $`\tau `$ via equation (3), we have to fit the cluster distribution with a single parameter power law. In this way, the critical point corresponds to the best fitted spectra, i.e. minimum $`\chi ^2`$.

The critical point can also be found by looking for the maximum (divergence in an infinite system) of the moments of the distribution  or of the normalized variance of the size of the maximum fragment (NVM).

The moments of the cluster distribution, $`M_k`$, are defined as:

$$M_k=\underset{A}{}A^kn_A,$$
(4)
where $`n_A`$ is the cluster distribution.

The normalized variance of the size of the maximum fragment (NVM)  is defined as:

$$NVM=<\frac{<A_{max}<A_{max}>>^2}{<A_{max}>}>,$$
(5)
where $`A_{max}`$ is the largest fragment of a given event and $`<\mathrm{}>`$ is an average over an ensemble of events at a given probability $`p`$.

Finally, the information entropy, as proposed in , is defined as:

$$S_1=\underset{A}{}p_A\mathrm{ln}(p_A),$$
(6)
where $`p_A`$ is the probability of detecting a fragment of mass $`A`$ defined as $`p_A=N_A/N_t,`$ with $`N_A`$ the number of fragments of size $`A`$ detected in a set of experiments and $`N_t=_AN_A`$ the total number of fragments detected in the same set of experiments. This magnitude is referred to as $`S_1`$ because it belongs to the family of generalized Renyi entropies defined by :

$$S_q=\{\begin{array}{cc}_ppln(p)\hfill & q=1\hfill \\ \frac{1}{q1}ln(_pp^q)\hfill & q1.\hfill \end{array}$$
(7)
In figure (1) we show the results of our calculations for percolation three dimensional lattices of two different sizes specified above. The $`NVM`$ in figure (1a) , $`M_2`$ in figure (1b), $`\chi ^2`$ in figure (1c) (the minimum of $`\chi ^2`$ denotes the best fitting power law with one-parameter fit), the $`\tau _{eff}`$ in figure (1d) (corresponding to a two-parameter fit of the cluster distribution. The minimum value of $`\tau `$ would correspond to the critical point) and $`S_1`$ in figure (1e).

We can see that for the mass spectra corresponding to a lattice of size $`30\times 30\times 30`$ (filled symbols of figure (1)) the signals $`NVM`$, $`M_2`$ and $`\chi ^2`$ give the same critical probability ($`p_c=0.25\pm 0.01`$). The extremes are very sharp and the critical probability is the one corresponding to the infinite system limit, indicating that the finite size effects are negligible. On the other hand, the maximum of the information entropy is not as sharp and it is also shifted (at $`p=0.225\pm 0.015`$) with respect to the extreme of the other signatures, showing that it is not a good signature of the critical point. Finally, we can see that when we fit the cluster distribution with a power law with two independent parameters the exponent $`\tau `$ is minimum at $`p=0.24\pm 0.01`$, slightly shifted from the other signals (but with overlapping error bars) and with a minimum less pronounced than when we fit the mass spectra with $`\tau `$ and $`q_0(\tau )`$ (minimum $`\chi ^2`$).

For the smaller ($`6\times 6\times 6`$) lattices (open symbols in figure (1)), we see that the peaks of the various signatures analyzed are slightly shifted from one another but with overlapping error bars, even for the case of the information entropy. The best fitted mass spectra according to the condition given by equation (3) is at $`p=0.31\pm 0.01`$ (figure (1c), the $`NVM`$ is maximum at $`p=0.30\pm 0.01`$ (figure (1a)), the $`M_2`$ at $`p=0.29\pm 0.01`$ (figure (1b)) and $`S_1(p_i)`$ at $`p_i=0.27\pm 0.05`$ (figure(1e)).

As we can see, when the system is small, finite size effects can introduce distortion and smear out the signature of the critical point. On the other hand, $`S_1`$ is shifted with respect the critical point when finite size effects are negligible.

In order to understand the behavior of the information entropy in this kind of systems, we have chosen to generate a family of clusters distributions according to equation (1) with the critical exponents corresponding to percolation in an infinite lattice in three dimensions. In this case, the normalization constant $`q_0`$ is taken as the inverse of the zeta Riemann function evaluated at $`\tau =2.18`$ and the scaling function $`f(z)`$ is the one fitted in . We have also chosen $`p_c=0.25`$ i.e. the infinite size value. These cluster distributions have the same functional form as the ones corresponding to a percolation problem, but without the percolating cluster and without finite size effects.

The information entropy can then be calculated by replacing $`n_A(p)`$ in the definition of $`S_1`$ (equation (6)) by its expression according to equation (1) . The results of this calculation are plotted in figure (2c) and we can see that the maximum in the entropy does not correspond to the pure power law distribution ($`p_c=0.25`$ in this case) as was claimed in  but to the distribution corresponding to $`p=0.22`$ showed in figure (2d). We can check this result by calculating the derivative of $`S_1`$ against $`ϵ=(pp_c)/p_c`$:

$$\frac{dS_1}{dϵ}=q_0^2\underset{s_1s_2}{}s_1^{(\tau )}s_2^{(\tau )}s_1^\sigma f^{}(s_1,ϵ)f(s_2,ϵ)log[(s_1^\tau f(s_1,\tau ))/(s_2^\tau f(s_2,\tau ))],$$
(8)
with $`f(s,ϵ)=f(z)`$ the scaling function and $`f^{}(z)=\frac{df}{dz}`$. If we evaluate equation (8) at the critical point, i.e. $`ϵ=0`$, we get:

$$\left(\frac{dS_1}{dϵ}\right)_{ϵ=0}=\frac{\tau }{\zeta ^2(\tau )}f^{}(0)\underset{s_1s_2}{}s_1^\sigma s_1^{(\tau )}s_2^{(\tau )}[log(s_1)log(s_2)]1.6.$$
(9)
It is clear from this expression that the information entropy does not show an extreme in the critical point for this family of distribution functions and, therefore, $`S_1`$ is not a signal of critical behavior.

In Figure (2a) we also plot the $`\chi ^2`$ coefficient for the cluster distribution fitted with one parameter ($`q_0=q_0(\tau )`$ according to equation (3)) and in figure (2b), the second moment of the distribution, $`M_2`$. It can be seen that these magnitudes signal the critical point very accurately.

Finally, we apply these tools to molecular dynamics simulations of fragmentation. These simulations have been described in previous works (see  for details). We follow the time evolution of excited drops composed of 147 particles interacting via a Lennard-Jones potential that undergo multifragmentation.

Using this model and taking the excitation energy as a control parameter we perform the same analysis as in the case of the percolation model. The corresponding results are displayed in figure (3). In this figure we can see that $`M_2`$ and $`NVM`$ peak at $`E=0.3\pm 0.2ϵ`$, meanwhile the best fitted mass spectra according to equation (3) (minimum $`\chi ^2`$) is, approximately, at $`E=0.15\pm 0.15ϵ`$. On the other hand the information entropy monotonically increases in the range of energies analyzed and therefore does not show a maximum when the system displays a power law mass spectra, i.e. in the energies around $`E=0.2ϵ`$.

In conclusion, we have analyzed different quantities proposed in the literature as signatures of critical behavior: i) the information entropy $`(S_1)`$, ii) the normalized variance of the size of the maximum fragment ($`NVM`$) , iii) the second moment of the distribution ($`M_2`$), iv) the best fitting power law and minimum $`\tau `$ . We found that $`S_1`$ is not a reliable signature of a second order phase transition as was proposed in . For this purpose, we have analyzed systems that display critical behavior: percolation on lattices of two sizes. For the larger one, where the finite size effects can be neglected, it is clear that the signature corresponding to information entropy is shifted with respect to the other quantities analyzed. But when the size is much smaller and the finite size can not be neglected, the signatures analyzed in this work display peaks that do not coincide in a given value of the probability. Instead, they are smeared in a small range (see figures (1)). In this case, the broad peak of the entropy can turn out to be rather close to this range of values (even though the error is large) and then the information entropy can be confused with the signatures.

This behavior is confirmed in the analysis of highly excited Lennard-Jones drops that undergo multifragmentation.

Finally, according to these results, we would choose as signatures of critical behavior the following signals: the best fit mass spectra with a one-parameter power law (minimum $`\chi ^2`$), the normalized variance of the size of the maximum fragment ($`NVM`$) or the second moment of the cluster distribution ($`M_2`$) and we would discard the information entropy ($`S_1`$).

Acknowledgments : We acknowledge partial financial support from the University of Buenos Aires via grant TW98 and Conicet via grant 4436/96. C.O.D is member of the Carrera del Investigador (CONICET). P.B is fellow of the CONICET .

Figure Captions

Fig.1: The normalize variance of the size of the maximum fragment, $`NVM`$ (a), the second moment of the cluster distribution $`M_2`$ (b), the $`\chi ^2`$ coefficient of a one parameter power law fitting (with $`q_0=q_0(\tau )`$) (c), the $`\tau `$ exponent from the fitting of two parameters power law (d) and the information entropy ($`S_1`$) (e) against bond probability for percolation lattices of sizes $`30\times 30\times 30`$ (full triangles) and $`6\times 6\times 6`$ (empty circles). In (f) we plot the cluster distribution for the critical probabilities at both sizes.

Fig.2: The $`\chi ^2`$ coefficient (a) and $`M_2`$ (b) against probability for the cluster distributions generated from equation (1) using the scaling function $`f(z)`$ obtained in . In (c) we plot the information entropy calculated by replacing the cluster distribution given by eq.(1) in the definition of $`S_1`$ (equation (6)). In (d) we can see both: the cluster distribution corresponding to the critical probability (chosen at $`p_i=0.25`$, a pure power law) and the one that maximize the information entropy ($`p_i=0.22`$). Here it is clear that the distribution that maximize the entropy is not a power law.

Fig.3: The $`\chi ^2`$ coefficient (a), the $`M_2`$ (b), the $`NVM`$ (c) and $`S_1`$ against the energy per particle for a molecular dynamics multifragmentation process for 147 particles interacting via a Lennard-Jones potential.


text:
# Selection of dune shapes and velocities. Part 2: A two-dimensional modelling.

## 1 Introduction

The beauty of the crescentic barchan dunes have recently attracted the interest of physicists for a better understanding and modelling of sand transport, as well as ripples and dunes formation and propagation. E. Guyon had this witty remark about them: ‘barchans are our drosophila’, which means that beyond the scientific and fundamental works on these dunes, we all keep in mind that such studies may lead to potential applications in the fight of saharan countries against sand invasion. One of the first reference work in the field is certainly the famous book of Bagnold which dates back from 1941 B41 . Since then, a great effort of measurement and modelling has been done which we have reviewed in details in the first part of these twin papers.

Our aim here is to discuss and model the selection of two-dimensional dune shape and velocity. For that purpose, we will simplify the model proposed by Sauermann et al. S01 ; SKH01 ; KSH01 . We will show that although rather severe approximations, we are able to recover their main results, in particular that dune profiles are not scale invariant, and that the inverse of the velocity grows almost linearly with the dune size. Besides, analytical expressions of dome and dune propagative profiles can be obtained, but whose coefficients have to be numerically computed. We furthermore get the existence of a critical mass below which no stable dune exists. The apparition of dunes can however be understood with the linear stability analysis of a flat sand sheet which is unconditionally unstable towards large wavelengths perturbations.

The paper is organized as follows. Section 2 is devoted to the equations of the model. The linear stability of a uniform sand sheet is treated in section 3. In section 4 we simplify further the equations and show what is the general shape of the purely propagative solutions of the model. The specific case of domes and ‘actual’ dunes are discussed in sections 5 and 6 respectively. At last, we conclude with a discussion of the relevance of these results, and the possible extension of the model to three-dimensional situations and dynamical studies.

## 2 Basic equations

We wish to give a description of the shape and evolution of two-dimensional dunes in terms of two fields: the profile $`h(x,t)`$ and the volumic sand flux $`q(x,t)`$ which is the volume of sand transported through an infinite vertical line per unit time. $`x`$ denotes the horizontal coordinate and $`t`$ is the time. We are going to write a set of three equations for these quantities in order to include into the model (i) the mass conservation, (ii) the progressive saturation of sand transport and (iii) the feedback of the topography on the sand erosion/deposition processes. Although we shall keep in this paper to two-dimensional situations which correspond to transverse dunes (invariant in the direction perpendicular to the wind), our ultimate goal is of course to be able to describe three dimensional dunes and barchans in particular.

Because barchans in dune fields – so-called ‘ergs’ – organize themselves like gooses or ducks during their migration flights, we named this class of models for purely graphical reasons the $`\text{C}_\text{c}^\text{c}`$ modellings. This denomination includes the approach of Sauermann et al. as well as the different variations and simplifications we derived in this paper from their work. At this stage of the modelling, we are however far from being able to take into account such ‘interactions’ between dunes which are necessary to explain the $`\text{C}_\text{c}^\text{c}`$ spatial organisation, and we shall focus on isolated objects only.

A simple balance calculation shows that the erosion rate $`_th`$ is directly related to the divergence of the flux $`q`$. This gives the common continuity equation:

$$_th+_xq=0.$$
(1)
The saturation effect of the sand transport has been already evoked in the first part of these twin papers: consider a patch of sand on which the wind is blowing – see figure 1. The flux of transported sand $`q`$ first increases and, because of the feedback of the grains on the wind velocity profile, get saturated after a typical length $`l_{sat}`$. In the first part of the paper, we showed that $`l_{sat}=\xi d\rho _{sand}/\rho _{air}`$ ($`d`$ is the grain diameter, $`\rho `$ the densities and $`\xi `$ a non dimensional prefactor), i.e. is almost independent of the wind shear velocity $`u_{}`$ – slow logarithmic dependencies hidden in $`\xi `$ only. This phenomenon has been reported and studied by several authors, e.g. B41 ; AH88 . The real shape of $`q(x)`$ is certainly more complicated than the one drawn on figure 1. In particular, oscillating or overshooting features were reported in B41  when $`q`$ reaches its asymptote. However, what is important for our purpose is only that a saturated value $`q_{sat}`$ is reached after a length $`l_{sat}`$. This space lag is satisfactorily described by the following equation:

$$_xq=\frac{q_{sat}q}{l_{sat}}.$$
(2)
This charge equation can be seen as a simplification of that proposed by Sauermann et al. in their continuum saltation model SKH01 . An important remark is that this equation is valid only if *some grains are available on the sand bed*. On a firm soil indeed, the flux cannot increase to become saturated. As suggested by Peer and Hakim PH01 , the right hand side of the relation (2) must be therefore multiplied by some matching function which quickly tends to zero when $`h`$ is decreased below, say, $`d`$ the grain diameter, and which is equal to unity above this value – the altitude of the firm soil is $`h=0`$. Then, the equation (2) becomes non linear but no boundary conditions have to be specified at the edges of the sand covered region. For example, if the matching function tends to zero like $`h`$, the dune will always keep a thin sand sheet at its back. But if it varies as $`\sqrt{h}`$, the dune will have a finite extension and will join the firm soil with an horizontal tangent. We shall ignore at present these subtleties but keep them in mind to invoke them later when necessary.

Another important remark is that the time scale on which the dune profile $`h`$ evolves is incomparably larger than that of the sand flux $`q`$. We then assume that $`q`$ adapts its profile instantaneously according to equation (2) and makes $`h`$ change slowly through equation (1). Therefore any term $`_tq`$ is irrelevant in this modelling.

The saturated flux $`q_{sat}`$ is uniform for a flat sand bed only. To the first order, the saturated flux $`q_{sat}`$ is a function of the local shear stress $`\tau =\rho _{air}u_{}^2`$ which itself depends – non locally – on the topography: basically, bumps and upwind slopes get more eroded than dips and downwind faces. A classical relationship between the saturated flux and the shear velocity that can be recovered with the scaling arguments of the part 1 of the paper is:

$$q_{sat}\frac{\rho _{air}}{\rho _{sand}}\frac{u_{}^3}{g}.$$
(3)
In principle, such a relationship is valid far from the velocity threshold $`u_{thr}`$ under which no sand can be eroded by the wind, i.e. $`q_{sat}=0`$ for $`u_{}u_{thr}`$. Refined formulas can be obtained which essentially smooth the step from $`0`$ to the previous asymptotic expression such as that obtained in part 1:

$$q_{sat}\frac{\rho _{air}}{\rho _{sand}}\frac{u_{}}{g}\left(u_{}^2u_{thr}^2\right).$$
(4)
In the whole range $`u_{thr}<u_{}<u_{thr}`$, $`q_{sat}`$ is null so that in practice, $`q_{sat}`$ cannot become negative on a dune. This condition will be used in section 5.

To close the equations, we have to explicit the spatial variations of the turbulent wind velocity due to the dune profile. The simplest model which verifies the basic requirements (see part 1) is certainly the perturbative calculation by Jackson and Hunt JH75 ; Weng91 . Neglecting logarithmic scale dependencies, Kroy et al. KSH01  have extracted the main features of their work by expressing the shear velocity as:

$$\frac{u_{}^2(x)}{U_{}^2}=1+A\frac{ds}{\pi s}_xh(xs)+B_xh(x),$$
(5)
where $`U_{}`$ is the shear velocity exerted on a flat bed. More precisely, as discussed in details in the part 1, $`A`$ and $`B`$ in principle depend on the size of the dune $`D`$ with $`\mathrm{ln}D/z_0`$ factors, where $`z_0`$ is the roughness of the sand surface. For $`D`$ varying between $`20`$ and $`200m`$ and a roughness of order of the grain size, such a logaritmic factor does not change by more than $`20\%`$ over the whole range. This justifies the fact that it is reasonable to take constant effective values for the coefficients $`A`$ and $`B`$. Several further important remarks must be made on this expression. First, it must be noted that the convolution integral acts on $`_xh`$ roughly like a derivative, leading to a term which encodes curvature effects. But this curvature is dimensionless and thus does not depend on the dune size – in other words, this term can be seen as a curvature rescaled by the dune size. It reflects the observation that the wind velocity increases on bumps (negative curvature) and decreases on hollows (positive curvature). Second, it is a non local term, meaning that the shear velocity depends on the whole shape of the dune. Of course sharp variations of the dune profile will also have a strong local effect. At last, the second term simply takes into account slope effects: positive slopes are more eroded than negative ones. Again, this term does not introduce any new lengthscale.

Expression (5) can be used to close up the set of equations as was done by Kroy et al. in KSH01  who used besides a more sophisticated – and non-linear – charge equation than (2). It is useful to simplify further the equation linking $`q_{sat}`$ to $`h`$ without loosing too much physics, in order to let more analytical developments. The linear expansion of the expression (5) rewritten in terms of $`q_{sat}u_{}^3`$ gives:

$$\frac{q_{sat}(x)}{Q_{sat}}=1+\frac{3}{2}A\frac{ds}{\pi s}_xh(xs)+\frac{3}{2}B_xh(x),$$
(6)
where $`Q_{sat}=q_{sat}(U_{})`$ is the saturated flux on a flat sand bed submitted to a shear velocity $`U_{}`$.

Furthermore, we shall use the saturation length $`l_{sat}`$ and flux $`Q_{sat}`$ to make our variables dimensionless. Thus, for a given wind shear velocity, all relevant scales of the problem are fixed. For instance, $`Q_{sat}/l_{sat}`$ is the velocity scale, $`l_{sat}^2/Q_{sat}`$ the time scale, $`Q_{sat}/l_{sat}^2`$ the frequency scale, etc. Note that the strength of the wind is completely encoded in these rescalings.

## 3 Stability of a flat sand bed

Before going further in the modelling of dune shape and propagation, we wish to investigate the problem of dune initiation. Indeed, there are two striking field observations. First, no persistent barchan dunes exist smaller than say, $`1m`$ high, $`20m`$ large and $`20m`$ long. Second, any small conical sandpile blown by the wind disappears even when a sand supply is provided. Then, how can barchan dunes appear?

To investigate this problem, we have integrated numerically equations (1,2,6) with $`A=6`$ and $`B=4`$. Two initial conditions were tested. First, a small triangular sandpile at the repose angle is prepared on the firm soil (figure 2). It can be seen that it is rapidly eroded and disappears, as observed in the field and in wind tunnel experiments – see part 1. Second, we look at the evolution of a thin sand sheet (figure 3) disturbed by a flat bump. This initial conditions mimics a sand beach on which sand is deposed by water. It can be seen that the bump propagates downwind and induces a strong erosion of the sand bed in front of it. A second bump nucleates from the initial perturbation which itself induces a strong erosion in front of it, and so on. After some time, a series of growing oscillations is generated. The amplification of this phenomenon stops when the oscillations eventually reach the firm soil from which no sand can be eroded, and/or when recirculation bubbles – see section 6 – appear and make these bumps interact. As a conclusion, depending on their spatial extension, small sand bumps either disappear or grow and initiate dunes.

It is then instructive to make the linear stability analysis of a flat sand bed. Let us consider an infinite uniform sand bed blown by a uniform wind. The sand flux is everywhere saturated: $`q=1`$ (in units of $`Q_{sat}`$). To investigate its stability, we can consider, without loss of generality a small perturbation of the form:

$`h(x,t)`$ $`=`$ $`He^{\sigma ti\omega t+ikx},`$ (7)
$`q(x,t)`$ $`=`$ $`1+Qe^{\sigma ti\omega t+ikx},`$ (8)
where $`Q`$ and $`H`$ are related one to the other by the conservation of matter,

$$(i\omega \sigma )H=ikQ.$$
(9)
¿From the relation (6), we get the following expression for the saturated flux: $`q_{sat}=1+3/2(|k|A+ikB)h`$. Once replaced in the saturation equation, it gives:

$$(1+ikl_{sat})Q=\frac{3}{2}(A|k|l_{sat}+Bik)H.$$
(10)
Combining equations (9) and (10) we finally obtain:

$`\sigma `$ $`=`$ $`{\displaystyle \frac{3k^2(BA|k|)}{2(1+k^2)}},`$ (11)
$`\omega `$ $`=`$ $`{\displaystyle \frac{3k|k|(A+B|k|)}{2(1+k^2)}}.`$ (12)
Note that a more complicated relation between $`q_{sat}`$ and $`u_{}`$ would only affect the prefactor, i.e. the time scale, in this calculation.

The growth rate $`\sigma `$, shown for $`A=6`$ and $`B=4`$ on figure 4 (top), is positive for small wavenumbers ($`k<B/A`$) and negative for large wavenumbers ($`k>B/A`$). A flat sand bed thus exhibits a large wavelength instability which can explain the initiation of dunes. But it is stable towards small wavelength disturbances, so that a small sandpile on a firm soil is quickly eroded.

The fastest growing mode – that which maximises $`\sigma `$ – is obtained for $`3k+k^3=2B/A`$. Neglecting the $`k^3`$ term, we obtain a still good approximation of the most unstable wavelength: $`\lambda =2\pi /k3\pi l_{sat}A/B`$. For $`A=6`$ and $`B=4`$ we get $`\lambda 14l_{sat}`$. For $`l_{sat}9m`$ we get a wavelength which is reasonable compared to what is observed on transverse dune fields in deserts. Under water $`l_{sat}`$ rather scales on the grain size and we have $`l_{sat}1cm`$ which gives again a good estimation of what has been measured by Betat et al. BKFR01  or Andersen et al. AGRL01  for example.

Finally, it can be seen from equation (12) that disturbances propagate downwind. We can compute the group velocity of these surface waves:

$$c=\frac{d\omega }{dk}=\frac{3[A|k|+Bk^2(3+k^2)]}{2(1+k^2)^2}.$$
(13)
$`c`$ is plotted on figure 4 (bottom) for $`A=6`$ and $`B=4`$. For small wavenumbers $`k`$, the group velocity $`c`$ increases linearly with $`k`$. This means that for asymptotically large wavelengths $`\lambda `$, the propagation speed scales as $`Q_{sat}/\lambda `$. We thus recover in the limit of large bumps, the scaling proposed by Bagnold. We also see that the velocity deviates significantly from this law, when the wavelength becomes comparable to the saturation length $`l_{sat}`$. This is also the case for actual dunes – see part 1.

## 4 A simple 2d modelling

### 4.1 Simplified equations

The expression (6) can be kept as it is, mixed with the mass conservation (1) and the charge (2) equations, and treated in Fourier space, but the idea is rather to replace the convolution expression by a simpler anzats. Since it is a ‘scaleless and non-local curvature’, we can replace it by $`D_{xx}h`$, where $`D`$ is the dune length. $`D`$ will be determined by the boundary conditions. We then end up with the following set of – linear – equations of our model:

$`_th+_xq`$ $`=`$ $`0,`$ (14)
$`_xq`$ $`=`$ $`q_{sat}q,`$ (15)
$`q_{sat}`$ $`=`$ $`1\alpha D_{xx}h+\beta _xh.`$ (16)
$`\alpha `$ and $`\beta `$ are two of the parameters of the model – a third one, $`\mu _b`$, will be introduced later on. They could be in principle measured independently. The $`\alpha `$ term has a stabilizing role on small perturbations, while the $`\beta `$ one make them growing. The rescalings $`h\alpha h`$ and $`t\alpha t`$ (with $`q`$ and $`x`$ unchanged) shows that only the ratio $`\beta /\alpha `$ is actually important. As stated in the figure captions, most of the curves of the paper have been plotted with $`\alpha =1.`$ and $`\beta =4`$. At last, as we already mentioned, the charge equation (15) is valid only if there is some sand to be eroded $`h>0`$. Note besides that, to be fully consistent, $`q_{sat}`$ should remain positive everywhere.

### 4.2 Boundary conditions

Before getting into the process of solving this system of equations, we need to specify what the boundary conditions at the edges of the dune are. This rises several important questions. The first thing to notice is that, when we transformed the relation (6) which links $`q_{sat}`$ to the profile $`h`$ with a convolution term into the simpler expression (16), we have incremented the order of the differential equation. Therefore, such a simplification requires an additional boundary condition than what is needed to integrate equations (6,14,15) for example. First order equations like the later need to know the profile $`h`$ at some position. However, when using the second order equations (14-16), the slope $`h^{}`$ should be also specified. As we already mentioned, one way to solve this problem is to regularise equation (15) by a non linear factor when $`h`$ tends to zero which gives conditions on $`h(x)`$ and $`h^{}(x)`$ as $`x\mathrm{}`$, see discussion after equation (2) in section 2.

An alternative is to consider the surface profile – the firm soil as well as the dune – as a whole. This is precisely what is done when the convolution (6) is used. If the slope $`h^{}`$ was not continuous at the dune edge, in some sense, $`h^{\prime \prime }`$ would be infinite. Then, from equation (16), the saturated flux would formally tend to $`\mathrm{}`$ over a region around this edge, i.e. would be forced to be null because it cannot be negative. Thus in this region, only sand deposition is permitted. This means that a discontinuity in the slope at the dune edge, immediately reacts to prevent the motion of this boundary point, which starts moving again only when the profile becomes flat, as shown on figure 5. Therefore, considering that the wind and thus the saturated flux are not sensitive to the transition from the firm soil to the dune back, gives a ‘natural’ upwind boundary condition for propagative solutions: the profile $`h`$ as well as its slope $`h^{}`$ should vanish at the upwind boundary of the dune. This argument actually also applies at the downwind edge, and we shall use it in the dome and dune next sections.

In the sequel, we are going to look for propagative solutions, i.e. functions of the type $`h(xct)`$ and $`q(xct)`$. We shall see that the explicit form of these functions can be found for given upwind conditions. However, such solutions are parametrized by coefficients such as this length $`D`$ or the propagative velocity $`c`$ which must be fixed by the right – downwind – conditions. Two kinds of right boundary conditions will be considered, leading to so-called ‘dome dunes’, i.e. without avalanche slip face, or ‘actual’ ones for which a ‘recirculation bubble’ will be introduced.

### 4.3 General form of the propagative solutions

For functions of $`xct`$, the continuity equation (14) can be easily integrated and gives

$$q=q_0+ch.$$
(17)
We now describe everything in the propagating frame referential, and rename $`x`$ as the new space coordinate. We look for isolated propagative objects. The point $`x=0`$ is the beginning of the dune where $`h=0`$. $`q_0`$ is thus the incoming sand flux – the sand supply. In the region $`x<0`$, no grains are available on the ground: $`h=0`$ and $`q=q_0`$ everywhere. In the region $`x>0`$, using the relation between $`q`$ and $`h`$ (17) in the charge equation (15) with $`q_{sat}`$ given by its expression (16), we get an ordinary differential equation for the dune profile $`h`$:

$$1\alpha Dh^{\prime \prime }+(\beta c)h^{}chq_0=0,$$
(18)
where $`h^{}`$ and $`h^{\prime \prime }`$ denote the first and second derivatives of the dune profile.

Under the two conditions $`h(0)=0`$ and $`h^{}(0)=0`$, the solution of equation (18) is

$$h(x)=\frac{1q_0}{c}\left[1+\left(\frac{s}{k}\mathrm{sin}(kx)\mathrm{cos}(kx)\right)e^{sx}\right],$$
(19)
where the coefficients $`s`$ and $`k`$ are given by

$`s`$ $`=`$ $`{\displaystyle \frac{\beta c}{2\alpha D}},`$ (20)
$`k`$ $`=`$ $`{\displaystyle \frac{1}{2\alpha D}}\sqrt{4c\alpha D(\beta c)^2}.`$ (21)
Figure 6 shows the comparison between the central longitudinal profile of an actual barchan and the theoretical form (19) which is parametrised by $`k`$, $`s`$ and $`c`$. It can be seen that the shape of the dune as well as the boundary conditions are well captured by the model. In fact, an oscillatory function as (19) does not look like an isolated dune yet. In particular, this $`h(x)`$ takes negative values. Besides, as already mentioned, the two coefficients $`D`$ and $`c`$ are not determined yet. In the next sections we show how it is possible to cut this solution at some point with an adequate right boundary condition in order to get a genuine propagative profile.

## 5 Domes

In this section, we look for so-called ‘dome’ solutions, i.e. dunes which do not show any avalanche plane. This is possible if the local slope of the dome is everywhere not steeper than some threshold $`\mu _b`$. In the next section, it will become clearer that this threshold corresponds to the slope at which the wind stream lines detach the dune profile, creating a backward wind flow – or a ‘recirculation bubble’ – which leads to a slip face. As will be explained in the conclusion, these dome solutions may play an important role when extending the present model to three dimensional situations. One of the major results of this section is that the shape of the dome, when it exists, is selected by the value of the incident flux $`q_0`$. Because we look for purely propagative solutions, $`q_0`$ of course is also the flux of the sand which leaves the dune.

Let us first take a very simple example. For $`c=\beta `$ the coefficient $`s`$ defined in equation (20) vanishes, and the solution (19) reduces to $`h(x)=\left[1\mathrm{cos}(kx)\right](1q_0)/\beta `$, where $`k=[\beta /(\alpha D)]^{1/2}`$. As the general solution (19) was constructed for, $`h`$ and $`h^{}`$ vanish at $`x=0`$, but also at $`x=2\pi /k`$. This value sets the length of the dune: $`D=4\pi ^2\alpha /\beta `$. Interestingly, this is precisely the cut-off $`\lambda _c`$ of the linear stability analysis. As discussed in the previous section, this simple solution satisfies the ‘natural’ downwind boundary conditions $`h(D)=0`$ and $`h^{}(D)=0`$. Another important restriction is that on the saturated flux $`q_{sat}`$ which must remain positive everywhere. As a matter of fact, there is a minimum value of $`q_0`$ below which $`q_{sat}`$ takes negative values and makes this solution inconsistent. It is easy to show that this lower bound for $`q_0`$ reads

$$q_0^{}=1\frac{1}{\sqrt{1+\left(\frac{\beta }{2\pi \alpha }\right)^2}}.$$
(22)
For the values $`\alpha =1.`$ and $`\beta =4.`$, we get $`q_0^{}0.16`$. Two such dome profiles for $`q_0=0.5`$ and $`q_0=q_0^{}`$ are displayed on figure 7. At last, note that the steepest slope of this solution is $`h_{min}^{}=\frac{1q_0}{2\pi \alpha }`$. For $`q_0=q_0^{}`$ and the same numerical values as above, it gives $`|h_{min}^{}|0.13`$. For any slope threshold $`\mu _b`$ larger than $`|h_{min}^{}|`$, the whole set of these solutions such that $`q_0^{}q_01`$ is acceptable. For $`\mu _b`$ below this (actually rather small) value however, $`q_0`$ must be larger than $`q_0^c=1.2\pi \alpha \mu _b`$, and no fully consistent dome solution can be constructed for a smaller incident sand flux. In the sequel, the typical value of $`\mu _b`$ that we shall use will be $`\mu _b=0.25`$.

How can we construct a dome solution for $`q_0<q_0^{}`$? Suppose the values of the velocity $`c`$ and the length $`D`$ are given, from equation (19) we can compute $`h`$ and its derivatives as well as $`q_{sat}`$. Because $`q_0`$ is smaller than $`q_0^{}`$, $`q_{sat}`$ will reach zero at some position $`x=L<D`$. Negative values of $`q_{sat}`$ are not permitted, and we therefore set it to zero for $`xL`$. Then, equation (15) is very easy to integrate and gives and exponential branch for the flux $`q(x)`$, which, using the linear relation (17) between $`q`$ and $`h`$ leads to:

$$h_r(x)=\frac{1}{c}\left[(cH_r+q_0)e^{(xL)}q_0\right].$$
(23)
The subscript $`r`$ is used to avoid any confusion with $`h(x)`$ given by equation (19), but both are part of the same dome solution. In particular, the dome must be smooth and in this expression $`H_r=h(L)`$ – we shall come back in the next paragraph on the continuity conditions at $`x=L`$. We cannot end the dome solution with $`h_r`$ for two reasons. First, we argued in the previous section that the dome profile must ends with a horizontal slope, which is not possible with an exponential function at finite distance $`D`$. Second, if one compute what would be the saturated flux $`q_{sat}^r`$ calculated from the profile $`h_r`$ with a relation like (16), one sees that it is negative as it should close to $`x=L`$ (remember that $`q_{sat}`$ is set to zero), but crosses zero at some other position $`x=L+R<D`$. From that point and for larger $`x`$ we thus need to come back to the original profile $`h`$. Then, a natural way to end up the dome profile is to use

$$h_y(x)=h(Dx),$$
(24)
which is consistent if $`q_{sat}(Y)=0`$, with $`L+R+Y=D`$. This choice ensures that both $`h(D)`$ and $`h^{}(D)`$ vanish as required. These three regions of the dome solution are illustrated on figure 8.

Let us now be more explicit about the continuity conditions at the two matching points. At $`x=L`$, by construction of the relation (23), the dome profile, the sand flux, and the saturated flux are continuous. Because equation (15) holds everywhere, it implies that $`_xq`$, and therefore the slope $`h^{}`$ are also continuous. By contrast, the curvatures $`h^{\prime \prime }(L)`$ and $`h_r^{\prime \prime }(L)`$ are different, and therefore $`q_{sat}^r(L)q_{sat}(L)=0`$. The position $`x=L+R`$ is defined by $`q_{sat}^y(L+R)=q_{sat}(Y)=0`$. At this point of course we do not want any step in the dome profile, such that $`H_yh_y(L+R)`$ and $`h_r(L+R)`$ must be equal. Again, because of equation (15) and the fact that $`q_{sat}`$ has been built to be continuous at $`x=L+R`$, the continuity of the profile makes the slope continuous too. However this position is also the point where the pseudo saturated flux $`q_{sat}^r`$ crosses zero, such that the curvature of the dome profile is also continuous. All these continuity conditions can be shown on figure 9.

To sum up, we have in practice four coefficients to determine: $`c`$, $`L`$, $`R`$ and $`Y`$ ($`D=L+R+Y`$) with the four non-linear following equations:

$`q_{sat}(L)`$ $`=`$ $`0,`$ (25)
$`q_{sat}(Y)`$ $`=`$ $`0,`$ (26)
$`h_r(L+R)`$ $`=`$ $`h(Y),`$ (27)
$`q_{sat}^r(L+R)`$ $`=`$ $`0,`$ (28)
which can be done numerically. One can find such a solution for any $`q_0<q_0^{}`$. However, to be consistent with the next section, and as already mentioned in the $`c=\beta `$ case, such a dome is acceptable only if its steepest slope is larger than $`\mu _b`$. This fixes a lower bound $`q_0^c`$ under which there are no solution. With $`\alpha =1.`$, $`\beta =4.`$ and $`\mu _b=0.25`$, we find that $`q_0^c0.09`$. The solution for this particular value of the incident sand flux is plotted on figure 9, and the dome profile for $`q_0=0.1`$ is compared, on figure 7, with two cases of the $`c=\beta `$ solutions.

When we have all coefficients of a particular solution, we can compute the mass of the corresponding dome:

$$M=_0^D𝑑xh(x).$$
(29)
For the $`c=\beta `$ solutions, this mass can be simply expressed as $`M=4\pi ^2\alpha (1q_0)/\beta ^2`$. The interesting point is that the function $`M(q_0)`$ decreases with its argument, see figure 10. Then, take a propagating dome under a flux $`q_0`$, and suppose that, for some reasons, the incident flux is suddenly increased by some small amount $`\delta q_0`$. The shape of the profile is such that the flux of the sand which leaves the dome is exactly $`q_0`$. Then, its mass gets larger. But a larger mass means that the flux of sand which can leave the dome is even smaller than before, which makes the mass of the dome increase again, and so on. In the same way, if the initial flux $`q_0`$ is suddenly decreased, the dome will loose more and more sand because less mass means a higher leaving sand flux. This negative feedback mechanism suggests that domes are unstable objects that may either quickly disappear or reach the point where their slope is steep enough to generate a recirculation bubble and create an avalanche slip face to become an actual dune. One could however observe them, for example under water with periodic boundary conditions BKFR01 , which can stabilize such an instability. Then an interesting prediction on the velocity of these domes which could be compared to experiments is shown in the inset of figure 10: it is constant for small domes and $`1/c`$ grows almost linearly with their size $`M^{1/2}`$ for larger ones. The dome length $`D`$ behaves also very much like $`1/c`$ as a function of $`M^{1/2}`$.

## 6 Dunes

It is possible to get another type of propagative profiles – transverse dunes – with a quite more sophisticated downwind boundary condition, namely a recirculation bubble. It is a common field observation (see part 1) that the wind streamlines on a dune follow exactly the shape of its back profile but separate at the point where the avalanche slip face begins – the brink – and reattach further downwind, as shown on figure 11. This phenomenon creates an eddy in the ‘shadow’ of the dune, where the wind is much less strong than anywhere else. As a consequence, all the sand eroded on the back is deposited around the top of the slip face which avalanches when the slope becomes too steep. As explained in the first part of these twin papers, it is fortunate that an accurate description of these avalanches is not necessary due to the fact that they do not have any feedback on the back profile of the dune: they simply relax their equilibrium slope $`\mathrm{tan}\phi `$. Because no grains can escape the dune from the slip face through the recirculation bubble, the net out flux is zero, which fixes $`q_0=0`$. Note that this is particular to two-dimensional situations that we are focused on: three dimensional barchan dunes loose sand from their horns.

### 6.1 The recirculation bubble

The simplest way to model the feedback of the recirculation bubble on the whole dune has been proposed by Zeman and Jensen JZ85  and used more recently by Kroy et al. KSH01 . The idea is to build an envelope of the dune which prolongs the dune profile by the separation streamline – see figure 11. To the first order, the wind on the back of the dune is the same as that would have been obtained if the envelope was solid. For example, the convolution integral in equation (5) used to calculate the shear stress on the soil should be applied not to the dune profile but to the dune + bubble envelope.

In our approach, the effect of the recirculation bubble on the wind is simply to modify the total length of the dune $`D`$ which becomes the dune length $`L`$ plus the bubble length $`R`$: $`D=L+R`$. The importance of the recirculation bubble becomes very clear: because it makes the apparent length of the dune be larger, it increases the erosion at the top of the dune, due to the curvature effect.

We are going to construct a very simple description of the recirculation streamline. Because it extends smoothly the dune shape, it will be possible to link its parameters to the dune ones, and to impose this way what could be called ‘non-local right boundary conditions’ to the solution expressed by equation (19).

Let us call $`h_b(x)`$ the separation streamline profile. As shown on figure 11, $`x=L`$ is the point at which this line starts and $`x=L+R`$ is its reattachment point with the soil. At both matching points, the extension of the dune or the soil by the separation streamline should be smooth:

$`h_b(L)`$ $`=`$ $`h(L)B,`$ (30)
$`h_b^{}(L)`$ $`=`$ $`h^{}(L)p,`$ (31)
$`h_b(L+R)`$ $`=`$ $`0,`$ (32)
$`h_b^{}(L+R)`$ $`=`$ $`0.`$ (33)
The recirculation bubble will be essentially governed by one important parameter: a slope that we already named $`\mu _b`$ in the previous section. It corresponds to the maximum angle that the wind streamlines can follow behind an obstacle. In other words, there is a flow separation if the slope is locally steeper than $`\mu _b`$. It is natural to impose that the steepest slope of the separation streamline should then be $`\mu _b`$:

$$h_b^{}(x_b)=\mu _b\text{with}h_b^{\prime \prime }(x_b)=0.$$
(34)
As a matter of fact, this steepest slope should be in principle less than or equal to $`\mu _b`$ otherwise this would mean that the streamlines actually detach at a higher threshold value. Besides, the bubble initially appears exactly at this slope on critical domes. The relation (34) can be then taken as the simplest self consistent slope condition.

Because the wind velocity profile is scale invariant, a further property is that the length $`R`$ of the recirculation bubble should scale on the dune size and no new length scale should be introduced by the bubble. Therefore, the curvature at $`x=L`$ must also be continuous:

$$h_b^{\prime \prime }(L)=h^{\prime \prime }(L)C.$$
(35)
Such a condition is also naturally observed when numerically integrating equations (6,14,15) – i.e. with the convolution term of equation (6) rather than the simplified relation (16). Furthermore, we noticed that the stability of the numerical scheme is very sensitive to the way the slope is computed at $`x=L`$: the only choice which does not create any instability is that which uses a formula which mixes the dune and bubble profiles $`h^{}(L)=[h_b(L+\mathrm{\Delta }x)h(L)]/\mathrm{\Delta }x`$, where $`\mathrm{\Delta }x`$ is the discretisation step.

In the absence of available data or models on flow separation behind a dune, the explicit form of the bubble profile we chose is the simplest that satisfies the requirements written above, namely a polynomial of $`3^{\text{rd}}`$ degree:

$$h_b=B+p(xL)+\frac{1}{2}C(xL)^2+\frac{1}{6}G(xL)^3.$$
(36)
Such a choice was also that of Kroy et al. in KSH01 . Using the conditions (30-35), $`G`$, $`C`$ and $`R`$ can be computed as functions of the slope $`p`$ and the height $`B`$ at the brink:

$`G`$ $`=`$ $`{\displaystyle \frac{9\mu _b}{R^2}}\left(1{\displaystyle \frac{2B}{3R\mu _b}}\sqrt{1{\displaystyle \frac{4B}{3R\mu _b}}}\right)`$ (37)
$`C`$ $`=`$ $`RG+\sqrt{2G\mu _b}`$ (38)
$`p`$ $`=`$ $`\mu _b+{\displaystyle \frac{C^2}{2G}}`$ (39)
We have checked that other similar choices for the parametrization of the separation streamline do not change the qualitative conclusions that are going to be presented in the next subsections. In a general way, the recirculation bubble conditions can be expressed as $`R/B=f_1(p)`$ and $`CB=f_2(p)`$, where the functions $`f_1`$ and $`f_2`$ encode the particular choice of the streamline separation profile $`h_b`$.

These two relations put together with the explicit expression of the dune profile (19) and the differential equation that it verifies (18), where $`D`$ has been set to the whole dune + bubble length $`D=L+R`$, let one to get three relations which link together the four parameters $`c`$, $`L`$, $`R`$ and $`B`$. Instead of plotting, say, the three first with respect to the last one, we rather chose to express, as we did for the domes in the previous section, all of them as functions of the total mass $`M`$ of the dune:

$$M=_0^L𝑑xh(x)+\frac{B^2}{2\mathrm{tan}\phi }.$$
(40)
We then get a continuous set of dune solution, from very large values of the mass, down to some cut-off value $`M_c`$ below which no stable recirculation bubble can be constructed – see below.

### 6.2 Profiles, dimensions and velocities

Let us fix a value of the mass of sand $`M`$ available to construct a dune. As for the domes, even for the simple choice of the recirculation bubble profile (36), a numerical resolution is required at this point to get the corresponding values of the parameters $`c`$, $`L`$, $`R`$ and $`B`$. Three examples of such solutions are plotted on figure 12. As evidenced on figure 13, when rescaled by a typical dimension of the dune – here the square root of the mass –, the different profiles do not collapse on a single curve. In particular the slope at the brink $`p=h^{}(L)`$ varies with respect to $`M`$, and even changes its sign for this choice of the parameters $`\alpha `$, $`\beta `$ and $`\mu _b`$. We get $`p<0`$ for the smallest dune, while $`p>0`$ for the largest one.

The inset of figure 13, as well as the curves of figures 14 and 15, which all show the dune and bubble features (lengths, slope, aspect ratios and velocity) as a function of the mass of the dune $`M`$, are all cut off at some value $`M_c`$ (big dot). This point corresponds to the smallest solution that can be constructed. Technically, one can see that the equations (37-39) which link bubble and dune parameters loose their sense when the quantity $`3R\mu _b2B`$ becomes negative. This point precisely gives the cut-off mass under which no bubble, and therefore no dune solution can be constructed. It is important to notice that this smallest dune has an avalanche slip face of finite size. Another choice for the function $`h_b(x)`$ would of course have given different values for $`M_c`$ as a function of the parameters of the model. But physically, this cut-off scale corresponds to the fact that such a recirculation bubble must have a minimum spatial extension to accommodate all continuity constrains. Not surprisingly, this extension is of the order of unity, that is to say of the order of $`l_{sat}`$ in non rescaled length units – which is the only length scale of our system. At last, it must be mentioned that the parameters of this smallest dune are pretty close to that of the largest dome, i.e. the dome for which the steepest slope is precisely equal to $`\mu _b`$, see table 1.

All results of figures 14 and 15 are remarkably simple and resemble very much field observations: the lengths $`L`$, $`R`$, $`H`$ and $`B`$ are almost straight lines as a function of the dune size $`M^{1/2}`$. Similarly, $`1/caM^{1/2}+b`$ with a very good precision. At last, it must be noted that the dune and bubble aspect ratios are not constant: large dunes are more compact with a proportionally larger bubble than small ones.

### 6.3 Parametric study

To complete the results of the previous subsection, we present here a parametric study of few quantities, namely the minimal dune height $`H_c`$ and the dune slope at the brink of asymptotically large dunes ($`M\mathrm{}`$). The parameter $`\alpha `$ is kept to unity, but $`\beta `$ has been varied between $`2`$ and $`6`$, and $`\mu _b`$ from $`0.2`$ to $`0.4`$.

Figure 16 shows that the variation of $`H_c`$ as a function of the parameters of the model is not very strong. Similarly, the aspect ratio of very large dunes is always of the same order of magnitude. The most interesting fact is perhaps that for small values of $`\beta `$, or for large values of $`\mu _b`$, the slope at the brink of very large dunes can remain negative – see figure 17. In other words, the brink and the crest can be always distinct, even for asymptotically large dunes.

One can understand intuitively the variations of $`H_c`$ and $`p`$ with $`\beta `$ and $`\mu _b`$. Increasing $`\beta `$ gives more strength to the destabilizing process, which lets small dunes appear at lower critical scale (smaller $`H_c`$) and makes large dunes more bumpy (larger $`p`$). If $`\mu _b`$ gets larger, the situation is reversed.

## 7 Conclusion

We have shown in this paper how, inspired from the work of Sauermann et al., one can build a simpler two dimensional model for the formation and the propagation of dunes. This modelling is based on two main variables: the dune profile $`h`$ and the volumic sand flux $`q`$. It includes three effects: (i) the mass conservation, (ii) the space lag over a length $`l_{sat}`$ for the sand flux to become saturated at some value $`q_{sat}`$, and (iii) the feedback of the profile on the saturated flux. In this third phenomenological equation, erosion and deposition processes are the result of the competition between two antagonist mechanisms: a stabilizing non local curvature term ($`\alpha `$) and a destabilizing slope one ($`\beta `$) which breaks the upwind-downwind symmetry.

Two kinds of solutions have been found: so-called ‘domes’ which do not show any avalanche slip face, and ‘real’ dunes for which a downwind recirculation bubble has been introduced. We were able to predict an analytical form for their propagative profiles, but whose coefficients have to be computed numerically.

The results of the model resemble very much field observations. We found for example that, due to this length scale $`l_{sat}`$, the dune profiles are not scale invariant: small dunes are flatter than large ones. Another point is that the inverse of the propagative velocity $`c`$ is, to a very good precision, almost linear with the size of the dune. This is consistent with Bagnold’s argument that $`cq_{sat}/H`$ for a dune of height $`H`$. In fact, this relation overestimates the velocity of small dunes for which the sand flux may not be already saturated at the crest and also because, compared to large dunes, the value of $`q_{sat}`$ is reduced due to a smaller curvature.

An important point discussed all through the paper was the issue of the boundary conditions. In particular, an important physical input was the ‘recirculation bubble’ behind the dune. This bubble makes the dune effectively look larger to the wind, and, due to the non local term in the relation between the saturated flux and the dune profile, has a stabilizing role. In fact, very little is known and well established about this recirculation bubble, but most of the dune features (position of the slip face, cut-off size, etc.) precisely depend on fine interactions and feedbacks between the dune and the bubble. More studies on this point are then urgently needed.

Another central result of the paper is the existence of this cut-off scale just mentioned, below which no dune can form. It corresponds to the fact that the bubble must have a minimum spatial extension to accommodate all continuity constrains. Not surprisingly, this scale is of order of $`l_{sat}`$. This result then rises the question of dune initiation and formation. The two scenari usually proposed by geophysicists for the formation of dunes are the following: first possibility, a small bump (of the size of ripples) grows continuously and forms a dune; second one, the sand accumulates on a solid obstacle like a rock or a bush and, when the size of the accumulation becomes larger than the obstacle, a dune forms and starts propagating downwind. However, observations show that ripples are stable and no structures between dunes and ripples can be seen. Similarly, rocks and bushes creates lee dunes of the size of the obstacle which remain anchored to the obstacle. An alternative explanation can be proposed, following the results of the stability analysis of the equations of the model, as well as that of the dome solutions. We found that large wavelengths perturbations get amplified, and that the dome profiles, selected by their incident flux $`q_0`$, are unstable to changes of that flux. Then, a possibility is that first domes form with a small height but directly with large length and width, and second that these domes progressively become more and more compact, and eventually reach the point where their slope is steep enough to generate a bubble and create an avalanche slip face to become an actual dune.

Several extensions to the present work can be thought of. First, we would like to go beyond the calculation of purely propagative solutions, and study the full dynamics of a given dune profile. In particular, as just said, an important point is the evolution of the dome solutions when submitted to incident sand flux variations. A second point is to go from a 2d description – transverse dunes – to real three dimensional situations. The idea is to ‘cut a barchan into longitudinal 2d slices’. As a matter of fact, a barchan slice close to the center of the dune looks like our dune solution, while a slice made at the edges where the horns are present rather have a dome shape. Suppose these slices are completely decoupled. Because the small ones go faster than the large ones, an initial conical sandpile will soon get a crescentic shape. However, when equilibrium is reached, all the slices should move at the same velocity. There should thus be a coupling between them, namely a lateral sand flux from the centre towards the horns. When the flux is saturated at the crest the velocity at the crest is $`c=(q_{sat}q_0)/H`$. This suggests than an equilibrium can indeed be achieved if $`q_0`$ increases in the small slices. Eventually the 3D dune slip face will be the sum of the contributions of all 2D slices whose brinks depend on their heights. Note that this scenario is consistent with the field observation that barchan horns are more elongated at strong winds which make the lateral sand flux less important, and consequently slices less coupled.

Finally, quantitative comparisons between experimental dune profiles and our theoretical predictions will be performed. This idea is to use barchan longitudinal slices as that shown on figure 6, but also sand structures under water, such as those obtained by Betat et al. BKFR01  or Andersen et al. AGRL01  for which, in principle, this model should be also valid.

###### Acknowledgements.

Acknowledgments
We are grateful to Gerd Sauermann and Klaus Kroy for precise explanations of all the details of their modelling. The measurements of the dune profiles (figure 6) were performed by B. Andreotti, S. Douady, P. Hersen and L. Quartier. We wish to thank Stéphane Douady, Christof Krülle and Brad Murray for useful discussions.


text:
# New Multicritical Random Matrix Ensembles

## 1 Introduction

Random matrix ensembles arise in an overwhelming number of diverse applications. Their utility stems from the fact that even though one can construct a multitude of different random matrix models, when one studies properties on the scale of eigenvalue spacing the results become universal and independent of the detailed structure of the random matrix measure. Hence one can extract predictions even without knowing the details of the appropriate microscopic model of the phenomenon being studied.

The type of universality class depends in general on the symmetry properties of the random matrix model and, of course, on the scaling properties of eigenvalue spacing. Thus one has different bulk and edge universality regimes, as well as a discrete series of multicritical universality regimes. Consequently only a discrete series of (rational) scaling exponents could be obtained.

In this paper we would like to show that there exists a real continuous range of universality regimes with generic real scaling exponents. We will perform an explicit construction for ensembles with eigenvalue spacing behaving like $`1/N^{1/(1+\alpha )}`$ with $`\alpha (0,1)`$.

A similar continuous series of universality classes appeared naturally in , for random matrix models with noncompact eigenvalue support and power-law tails in the eigenvalue density. There the universal behaviour appeared for large eigenvalues. In general, however, the appropriate measures were difficult to obtain in an explicit way. This was the main motivation for this investigation.

In this paper we would like to perform a construction of a more conventional class of models with compact support and where the universal regime is located close to $`\lambda =0`$.

We strongly suspect that a mapping of this behaviour to infinity would correspond to the universal scaling for the Lévy ensembles of . Another possible application of similar random matrix models (in their chiral variant) might be to model the behaviour of low lying eigenvalues of the QCD Dirac operator at the chiral phase transition. In this context the exponent $`\alpha `$ in $`\rho (\lambda )|\lambda |^\alpha `$ is the inverse of the $`\delta `$ critical exponent of the chiral phase transition .

The plan of this paper is as follows. In section 2 we briefly recall the main features of the old and new universality regimes in random matrix theory. In section 3 we construct the random matrix measures of the new multicritical ensembles and compute the eigenvalue density in the bulk of the spectrum. In section 4 we move to a discussion of the microscopic (presumably universal) scaling behaviour and propose a simple approximate method based on the distribution of zeroes of orthogonal polynomials. In section 5 we compare these approximate formulas with exact expressions calculated numerically for $`\alpha =1/3`$. We close the paper with a discussion.

## 2 Universality regimes

The universality of a correlation function in a random matrix model means that it does not change (up to a trivial rescaling) when the probability measure is modified. Typically the universal quantities involve properties of the eigenvalues and their correlations considered on the scale of eigenvalue spacing in the large $`N`$ limit. For a multicritical point (see below), the modifications of the probability measure are restricted to those that do not destroy multicriticality. However, in general, this still leaves an infinite dimensional space of possible deformations.

The type of universal behaviour in random matrix models depends on the part of the spectrum that one is studying. In general the scaling of eigenvalue spacing with $`N`$, in the vicinity of a fixed eigenvalue $`\lambda _0`$, can be easily obtained by looking at the local behaviour of the (bulk) eigenvalue density (normalized to $`N`$):

$$N\rho (\lambda )N(\lambda \lambda _0)^\alpha $$
(1)
Consequently the number of eigenvalues between $`\lambda _0`$ and $`\mathrm{\Lambda }`$ is approximately $`nN(\mathrm{\Lambda }\lambda _0)^{\alpha +1}`$. Reexpressing $`\mathrm{\Lambda }`$ in terms of $`n`$ shows that the eigenvalue spacing in the vicinity of $`\lambda _0`$ scales like $`1/N^{1/(1+\alpha )}`$.

In the bulk the eigenvalue spacing is of the order $`1/N`$ — this is the classical Wigner-Dyson regime . At the edges of the spectrum the spacing is like $`1/N^{2/3}`$ and one observes universal Airy-like oscillations . When the random matrix posseses an additional chiral structure as in applications to QCD i.e.

$$=\left(\begin{array}{cc}0& 𝒜\\ 𝒜^{}& 0\end{array}\right)$$
(2)
there appears a different universal regime close to $`\lambda =0`$ (with eigenvalue spacing $`1/N`$) .

Apart from these universal regimes, when one finetunes the potential so that $`\rho (\lambda )\lambda ^{2m}`$, a discrete series of multicritical points appear (with scaling $`1/N^{1/(2m+1)}`$). The behaviour on the scale of eigenvalue spacing in this regime is very difficult to extract . A different (but also discrete) class of multicritical models appears when adding an appropriate fixed matrix to the random matrix. This behaviour and novel scaling has been analyzed both for ordinary  and chiral  random matrix models.

From the above discussion we see that in order to construct a class of random matrix models with arbitrary real scaling exponents we have to find a random matrix ensemble whose bulk eigenvalue spacing behaves like

$$\rho (\lambda )|\lambda |^\alpha $$
(3)
for real $`\alpha `$. We will perform the construction for $`\alpha (0,1)`$ but an extension to other $`\alpha `$ should not be difficult.

## 3 Eigenvalue density in the bulk

For a hermitian random matrix model with measure

$$e^{N\text{tr}V()}$$
(4)
we consider the following class of potentials

$$V()=\frac{t}{p}||^p+\frac{g}{2p}||^{2p}$$
(5)
For hermitian matrices we may define $`||^p`$ in the following way. By a unitary transformation $``$ can be rewritten in the form $`=U\mathrm{\Lambda }U^{}`$ with $`\mathrm{\Lambda }=diag(\lambda _1,\mathrm{},\lambda _N)`$ a diagonal matrix. Then we may define $`||^pU|\mathrm{\Lambda }|^pU^{}`$ where $`|\mathrm{\Lambda }|^p=diag(|\lambda _1|^p,\mathrm{},|\lambda _N|^p)`$. In fact we only need its trace which is given by

$$\text{tr}|M|^p=\underset{i=1}{\overset{N}{}}|\lambda _i|^p$$
(6)
We note that although random matrix models with potentials of the type (5) are perfectly well defined, they seem to lack an evident diagrammatic perturbative expansion. Potentials with a single power-like term were first considered in the mathematical literature  and in . The parameter $`p`$ will be related to $`\alpha `$ by

$$p=1+\alpha $$
(7)
In the bulk we solve for the eigenvalue density by standard saddle point method. We are looking for a solution with a single cut $`(a,a)`$. Later we will choose $`a=1`$. The expression for the Green’s function is then

$$G(z)=\frac{1}{2\pi }\sqrt{z^2a^2}_a^a\frac{d\lambda }{z\lambda }\frac{V^{}(\lambda )}{\sqrt{a^2\lambda ^2}}$$
(8)
We note that for $`\alpha >0`$, the derivative exists. The cut endpoint is fixed through the constraint

$$\frac{1}{\pi }_0^a𝑑\lambda \frac{\lambda V^{}(\lambda )}{\sqrt{a^2\lambda ^2}}=1$$
(9)
The eigenvalue density can be extracted from the imaginary part through

$$\rho (z)=\frac{1}{2\pi ^2}\sqrt{a^2z^2}PV_a^a\frac{d\lambda }{z\lambda }\frac{V^{}(\lambda )}{\sqrt{a^2\lambda ^2}}$$
(10)
Because of the nonanalytic power like structure of the potential (5), it is convenient to get rid of the principal value and rewrite the formula as an ordinary integral. First using the fact that the potential is even we may use

$$\frac{1}{z\lambda }\frac{1}{z+\lambda }=\frac{2\lambda }{z^2\lambda ^2}$$
(11)
to obtain

$$\rho (z)=\frac{1}{\pi ^2}\sqrt{a^2z^2}PV_0^a\frac{d\lambda }{z^2\lambda ^2}\frac{\lambda V^{}(\lambda )}{\sqrt{a^2\lambda ^2}}$$
(12)
Now we use the fact that

$$PV_0^a\frac{d\lambda }{z^2\lambda ^2}\frac{1}{\sqrt{a^2\lambda ^2}}=0$$
(13)
Therefore we may rewrite (12) as

$$\rho (z)=\frac{1}{\pi ^2}\sqrt{a^2z^2}_0^a\frac{zV^{}(z)\lambda V^{}(\lambda )}{z^2\lambda ^2}\frac{1}{\sqrt{a^2\lambda ^2}}$$
(14)
where we could erase the principal value as the integrand is now nonsingular.

We will now fix the endpoint of the cut to $`a=1`$. Substitution of (5) into (9) yields

$$t\frac{\mathrm{\Gamma }\left(\frac{1+p}{2}\right)}{\sqrt{\pi }p\mathrm{\Gamma }\left(\frac{p}{2}\right)}+g\frac{\mathrm{\Gamma }\left(\frac{1+2p}{2}\right)}{\sqrt{\pi }2p\mathrm{\Gamma }\left(p\right)}=1$$
(15)
Since we want to locate a multicritical point we have to require that $`\rho (0)=0`$. This gives the second equation for $`t`$ and $`g`$:

$$_0^1\frac{V^{}(\lambda )d\lambda }{\lambda \sqrt{1\lambda ^2}}=0$$
(16)
i.e.

$$t\frac{\mathrm{\Gamma }\left(\frac{p1}{2}\right)}{\mathrm{\Gamma }\left(\frac{p}{2}\right)}+g\frac{\mathrm{\Gamma }\left(\frac{2p1}{2}\right)}{\mathrm{\Gamma }\left(p\right)}=0$$
(17)
The solution is

$`t`$ $`=`$ $`{\displaystyle \frac{2\sqrt{\pi }p\mathrm{\Gamma }\left(p\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{p}{2}\right)}{\mathrm{\Gamma }\left(p+\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{p1}{2}\right)2\mathrm{\Gamma }\left(p\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{1+p}{2}\right)}}`$ (18)
$`g`$ $`=`$ $`{\displaystyle \frac{2\sqrt{\pi }p\mathrm{\Gamma }\left(\frac{p1}{2}\right)\mathrm{\Gamma }\left(p\right)}{\mathrm{\Gamma }\left(p+\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{p1}{2}\right)2\mathrm{\Gamma }\left(p\frac{1}{2}\right)\mathrm{\Gamma }\left(\frac{1+p}{2}\right)}}`$ (19)
#### Bulk eigenvalue density

Once the potential is fixed let us compute explicitly the bulk eigenvalue density

$$\rho (z)=\frac{1}{\pi ^2}\sqrt{1z^2}\left[t\varphi _p(z)+g\varphi _{2p}(z)\right]$$
(20)
where

$$\varphi _p(z)_0^1\frac{z^p\lambda ^p}{z^2\lambda ^2}\frac{1}{\sqrt{1\lambda ^2}}$$
(21)
This integral can be explicitly expressed in terms of hypergeometric functions

$$\varphi _p(z)=\frac{\sqrt{\pi }\mathrm{\Gamma }\left(\frac{p1}{2}\right)}{2\mathrm{\Gamma }\left(\frac{p}{2}\right)}{}_{2}{}^{}F_{1}^{}(1,1\frac{p}{2},\frac{3p}{2};z^2)+\frac{\pi }{2}\mathrm{tan}\frac{\pi p}{2}\frac{z^{p1}}{\sqrt{1z^2}}$$
(22)
We see that it is the last term which gives the multicritical behaviour that we wanted to obtain $`\rho (z)z^{p1}=z^\alpha `$. Owing to equation (16) the constant terms cancel. We could have chosen of course a different power for the second term in (5), however with this choice the numerical construction of the relevant orthgonal polynomials is made easier.

Let us note what happens when we reach the point $`\alpha =1`$. Then $`p=2`$ and the potential obtained here is the standard multicritical quartic one $`V(x)=4x^2+4x^4`$. However due to the vanishing of the coefficient of the last term in (22) the eigenvalue density behaves like $`\rho (z)z^2`$ and not $`z`$. When we increase $`\alpha `$ the expression for the eigenvalue density ceases to be positive and thus the construction fails. Presumably a modification of the power in the second term of the potential might cure the problem but we will not consider that here. To sum up, the expression (20) is nonnegative when $`\alpha (0,1)`$ as considered in this paper.

## 4 Behaviour near the origin — orthogonal polynomials

The ultimate interest in constructing the random matrix models with the new scaling properties is to extract universal properties which typically occur in the microscopic regime i.e. on the scale of eigenvalue spacing. It is well known that all the relevant properties are encoded in the orthogonal polynomials

$$_{\mathrm{}}^{\mathrm{}}𝑑xP_n(x)P_m(x)e^{NV(x)}=\delta _{nm}$$
(23)
It is convenient to introduce the wavefunctions

$$\psi _n(x)=P_n(x)e^{N\frac{V(x)}{2}}$$
(24)
Then the kernel which allows for the determination of all correlation functions is given by the expression

$$K(z,w)=\sqrt{R_N}\frac{\psi _N(z)\psi _{N1}(w)\psi _{N1}(z)\psi _N(w)}{zw}$$
(25)
where $`R_N`$ is the recursion coefficient entering $`z\psi _n(z)=\sqrt{R_{n+1}}\psi _{n+1}(z)+\sqrt{R_n}\psi _{n1}(z)`$.

The microscopic limit is defined through a rescaling

$$x=\frac{y}{N^{\frac{1}{1+\alpha }}}$$
(26)
and a limit $`N\mathrm{}`$ with $`y`$ kept fixed. The key quantity that determines the universal properties is the rescaled wavefunction

$$\underset{N\mathrm{}}{lim}C_N\psi _N\left(\frac{y}{N^{1/(1+\alpha )}}\right)=F(y)$$
(27)
(as well as a similar expression with $`\psi _{N1}`$), $`C_N`$ is a normalization constant chosen in such a way that the limit exists. In fact the existence of such a limit is a nontrivial and nongeneric property.

In general the determination of the scaling function $`F(y)`$ is a very difficult problem c.f. . One can formulate a differential equation satisfied by the wavefunction $`\psi _n(x)`$, but taking the scaling limit is extremely difficult and e.g. for the case of quartic multicritical ensemble involves data from an auxiliary Painlevé equation . Moreover the starting point of such considerations requires a polynomial potential, which is not the case for our class of models. It would be interesting in the future to explore the possibility of applying the Riemann-Hilbert methods of .

Here we would like to adopt a different approach and give a very simple but approximate method of constructing the scaling function $`F(y)`$. In section 5 we will compare the approximate solution with the numerically obtained exact result for $`\alpha =1/3`$, and with the mesoscopic approximation of .

### Zeroes of the orthogonal polynomials

The starting point of our construction is the elementary fact that we may reconstruct the orthogonal polynomial from the knowledge of its zeroes. Since in any case we normalize the (even) wave functions through $`\psi _{2n}(0)=1`$ the appropiate formula is222We always consider even potentials.

$$\psi _{2n}(x)=e^{N\frac{V(x)}{2}}\underset{i=1}{\overset{n}{}}\left(1\frac{x^2}{\lambda _i^2}\right)$$
(28)
Interestingly enough there is a theorem due to Ismail  which states that for an (almost) arbitrary potential $`V(x)`$ (see ) the zeroes of the $`n^{th}`$ orthogonal polynomial $`\lambda _i`$ solve the equations

$$\underset{1kn,ik}{}\frac{1}{\lambda _i\lambda _k}=\frac{NV^{}(\lambda _i)}{2}+\frac{1}{2}\left(\mathrm{log}\frac{A_n(\lambda _i)}{\sqrt{R_n}}\right)^{}$$
(29)
i.e. these are electrostatic equilibrium positions in an external field. The additional assumptions present in  serve only to prove the uniqueness of a solution to (29).

The correction term involves the function $`A_n(x)`$, given by the formula

$$\frac{A_n(x)}{\sqrt{R_n}}=N_{\mathrm{}}^{\mathrm{}}\frac{V^{}(x)V^{}(y)}{xy}P_n^2(y)e^{NV(y)}𝑑y$$
(30)
For polynomial potentials of degree $`m`$, $`A_n(x)`$ is a polynomial of at most degree $`m2`$. Note that the above equations are exact and valid for any $`N`$. We see that the last term in (29) is suppressed w.r.t. the ordinary potential term.

Now in order to develop our approximation scheme let us neglect the correction factor $`(\mathrm{log}A_n(x))^{}`$ (this is an approximation, the $`A_n`$ term may indeed give some contributions in the scaling limit). From (29) we see that the zeroes are distributed with a continuum density identical to the eigenvalue density $`\rho (\lambda )`$ of the corresponding random matrix model. The basis of our approximation scheme is the assumption that locally the zeroes are distributed uniformly w.r.t. $`\rho (\lambda )`$ (see the examples below).

#### Classical Wigner-Dyson scaling functions

Before we consider the case of immediate interest to us, in order to illustrate the method let us first rederive the formulas for the standard Wigner-Dyson scaling functions. To this end we assume that the eigenvalue density at $`\lambda =0`$ is nonvanishing and $`N`$ is very large. For simplicity we take the potential to be even. It is convenient to normalize the continuum density of zeroes by

$$_{\mathrm{}}^{\mathrm{}}\rho _{cont}(\lambda )𝑑\lambda =N$$
(31)
In our approximation $`\rho _{cont}(\lambda )`$ is simply given by $`N\rho (\lambda )`$ where $`\rho (\lambda )`$ is the bulk eigenvalue density of the random matrix model (normalized to 1). If the number of zeroes is even then they are distributed symmetrically around 0 and there is no zero at $`\lambda =0`$. The assumption of ‘uniform distribution’ means that the $`ith`$ zero is located at $`\lambda _i`$ where

$$_0^{\lambda _i}\rho _{cont}(\lambda )𝑑\lambda =i+\frac{1}{2}$$
(32)
Since we are interested only in zeroes which lie close to $`\lambda =0`$, in the above expression we may substitute $`\rho _{cont}(\lambda )N\rho (0)`$. Consequently $`\lambda _i=(i+1/2)/(N\rho (0))`$. Therefore the wavefunction is

$$\psi (x)=e^{N\frac{V(x)}{2}}\underset{i=0}{\overset{\mathrm{}}{}}\left(1\left(\frac{\rho (0)xN}{i+1/2}\right)^2\right)$$
(33)
Here we extended the upper limit of the product to infinity. After introducing the scaling variable $`y=xN`$ we obtain immediately

$$F(y)=\mathrm{cos}\left(\pi \rho (0)y\right)$$
(34)
the exact scaling function in the bulk. The odd case is similar but then there is a zero at $`\lambda =0`$ and consequently the relevant equation is

$$_0^{\lambda _i}\rho (\lambda )𝑑\lambda =i$$
(35)
From the infinite product representation one obtains then the sine function. The kernel (25) in the scaling limit follows immediately

$`K(y_1,y_2)`$ $``$ $`{\displaystyle \frac{\mathrm{sin}\left(\pi \rho (0)y_1\right)\mathrm{cos}\left(\pi \rho (0)y_2\right)\mathrm{cos}\left(\pi \rho (0)y_1\right)\mathrm{sin}\left(\pi \rho (0)y_2\right)}{y_1y_2}}`$ (36)
$`=`$ $`{\displaystyle \frac{\mathrm{sin}\left(\pi \rho (0)(y_1y_2)\right)}{y_1y_2}}`$
The fact that we obtained here the exact scaling functions does not mean, however, that we should expect to get exact results in general.

#### Approximate scaling function for $`\alpha (0,1)`$

The analysis of the scaling function for the new multicritical ensembles is very similar. For definiteness we will just consider the even case. We start from the behaviour of the eigenvalue density near zero:

$$\rho (\lambda )=c_\alpha |\lambda |^\alpha $$
(37)
where $`c_\alpha `$ can be easily extracted from (20). The positions of the zeroes according to the ‘uniform distribution’ approximation are determined by

$$_0^{\lambda _i}Nc_\alpha \lambda ^\alpha 𝑑\lambda =\frac{1}{\alpha +1}c_\alpha N\lambda _i^{\alpha +1}\frac{1}{\alpha +1}c_\alpha y_i^{\alpha +1}=i+\frac{1}{2}$$
(38)
where we introduced the rescaled variable $`y=\lambda N^{1/(1+\alpha )}`$. Hence

$$y_i=\left[\frac{1+\alpha }{c_\alpha }\left(i+\frac{1}{2}\right)\right]^{\frac{1}{1+\alpha }}$$
(39)
Consequently the scaling function normalized by $`F(0)=1`$ is

$$F(y)=\mathrm{exp}\left\{\frac{t}{1+\alpha }\frac{y^{1+\alpha }}{2}\right\}\underset{i=0}{\overset{\mathrm{}}{}}\left(1\frac{y^2}{y_i^2}\right)$$
(40)
In the next section we will compare this approximate result with the exact orthogonal polynomials for $`\alpha =1/3`$. However because the convergence properties of the infinite product are not very good we will truncate it at some $`n_{max}`$ (for numerical computation we use $`n_{max}=200`$) and approximate the rest of the terms through

$$\underset{i=n_{max}+1}{\overset{N}{}}\left(1\frac{y^2}{y_i^2}\right)\mathrm{exp}\left\{\mathrm{log}\left(1\frac{y^2}{y_i^2}\right)\right\}\mathrm{exp}\left\{y^2_{y_{n_{max}+1}}\frac{\rho (\lambda )}{\lambda ^2}\right\}$$
(41)
The result is

$$\mathrm{exp}\left(\frac{1}{1\alpha }c_\alpha \frac{y^2}{y_{n_{max}+1}^{1\alpha }}\right)$$
(42)
As we see when taking the $`n_{max}\mathrm{}`$ limit this term becomes equal to unity (because $`y_{n_{max}+1}\mathrm{}`$). However we include it just because of the slow numerical convergence of the infinite product.

## 5 An example — $`\alpha =\frac{1}{3}`$

The orthogonal polynomials for any potential can be constructed using the determinant formula:

$$P_n(x)det\left(\begin{array}{ccccc}\mathrm{\Delta }_0& \mathrm{\Delta }_1& \mathrm{\Delta }_2& \text{}& \mathrm{\Delta }_n\\ \mathrm{\Delta }_1& \mathrm{\Delta }_2& \mathrm{\Delta }_3& \text{}& \mathrm{\Delta }_{n+1}\\ \mathrm{\Delta }_2& \mathrm{\Delta }_3& \mathrm{\Delta }_4& \text{}& \mathrm{\Delta }_{n+2}\\ \text{}& \text{}& \text{}& \text{}& \text{}\\ x^0& x^1& x^2& \text{}& x^n\end{array}\right)$$
(43)
where $`\mathrm{\Delta }_i`$ are the moments

$$\mathrm{\Delta }_i=_{\mathrm{}}^{\mathrm{}}x^ie^{NV(x)}𝑑x$$
(44)
The advantage of the specific form of (5) is that the moments $`\mathrm{\Delta }_i`$ can be expressed analytically in terms of the confluent hypergeometric functions $`{}_{1}{}^{}F_{1}^{}`$.

We constructed orthogonal polynomials for $`\alpha =1/3`$ ($`p=4/3`$) for $`N=40,120`$ and $`N=240`$. In Fig. 1 we show the wavefunctions expressed in terms of the scaling variable $`y=xN^{3/4}`$

$$\psi _N^{sc}(y)\psi _N\left(y/N^{\frac{3}{4}}\right)=P_N\left(y/N^{\frac{3}{4}}\right)e^{\frac{N}{2}V\left(y/N^{\frac{3}{4}}\right)}$$
(45)
for the above values of $`N`$, normalized by $`\psi _N(0)=1`$. We see convergence towards a well defined scaling function.

In Fig. 2 we compare the exact wavefunction for $`N=240`$ with the approximate scaling function $`F(y)`$ obtained in the previous section. There is no free parameter in $`F(y)`$. The agreement is indeed surprisingly good. Even the small ‘bump’ close to $`y=0`$ is correctly reproduced. In fact it is difficult to judge from the numerical comparision whether the approximation is exact or not in this case. The small deviations might be caused either by true corrections which go beyond our approximation or by finite size effects. The reason why such a simple approximation scheme works so well certainly deserves further study.

Finally let us say a few words about the comparison with the ‘mesoscopic approximation’ of . An analysis of the differential equations for $`\psi _N(x)`$ for the potentials considered here is still lacking, so we will just take a suitable analytical continuation of the mesoscopic approximation from the discrete multicritical points considered in  to our case. The result is just

$$F_{meso}(y)=\mathrm{cos}\left(\frac{\pi }{\alpha +1}c_\alpha y^{\alpha +1}\right)$$
(46)
We see that the zeroes of this function on the real axis coincide with the zeroes of our approximation, but the analytical structure in the complex plane is certainly different. There are spurious cuts and zeroes coming from $`y^{4/3}`$. Nevertheless this approximation is also quite good but the substructure close to $`y=0`$ is not captured by the mesoscopic approximation (similarly as in the quartic multicritical case considered in ). The comparison is presented in Fig. 3.

## 6 Discussion

In this paper we constructed a continuous class of random matrix ensembles which exhibit a new type of universal behaviour close to the zero eigenvalue. The eigenvalue spacing in this region behaves like $`1/N^{1/(\alpha +1)}`$ with a real exponent $`\alpha (0,1)`$. We derived the bulk eigenvalue density and studied approximately the behaviour of the scaling functions in the microscopic limit. Our method for obtaining the scaling function was approximate and based just on the local behaviour of the bulk eigenvalue distribution. Nevertheless the approximation seems to be very good and captures even the fine structure of the scaling functions near the origin.

It would be very interesting to see if these new universality classes could appear in some physical systems. We conjecture that the scaling behaviour can be directly related to the large eigenvalue behaviour of Lévy random matrix models .

There are numerous open questions and directions for further study. Firstly it would be very interesting to try to obtain, even implicitly through e.g. a differential equation, the exact scaling function, especially as the vast majority of techniques dealing with multicritical random matrix ensembles requires from the outset a polynomial potential. Perhaps the most promising approach would be the Riemann-Hilbert method. On the technical side a more systematic analysis of the properties of the orthogonal polynomial zeroes starting from (29), in particular the nature and scaling properties of corrections would be very welcome. This is especially interesting in order to better understand the effectiveness and limitiations of the approximate method which seems to work so well here. In general, however, this seems to be quite a formidable problem. Secondly other quantities of interest, such as eigenvalue spacing distributions, would be important in view of possible applications. Other possible directions are an extension of the above considerations to chiral random matrix models (in view of possible applications to QCD), and a rigorous treatment of universality properties.

Note added: An independent study of a very similar class of multicritical random matrix model with real exponents was undertaken by Gernot Akemann and Graziano Vernizzi .

Acknowledgments. I would like to thank Graziano Vernizzi for interesting discussion on the bulk properties, Poul Damgaard for discussion on  and Maciej A. Nowak for discussions and comments. This work was supported in part by KBN grant 2P03B01917.


text:
# Current-induced dendritic magnetic instability in superconducting MgB2 films

## Abstract

Magneto-optical imaging reveals that in superconducting films of MgB2 a pulse of transport current creates avalanche-like flux dynamics where highly branching dendritic patterns are formed. The instability is triggered when the current exceeds a threshold value, and the superconductor, shaped as a long strip, is initially in the critical state. The instability exists up to 19 K, which is a much wider temperature range than in previous experiments, where dendrites were formed by a slowly varying magnetic field. The instability is believed to be of thermo-magnetic origin indicating that thermal stabilization may become crucial in applications of MgB2.

Resubmitted to Appl. Phys. Lett. on 15.04.2002

Following the discovery of superconductivity below 39 K in polycrystalline MgB2 a tremendous effort is now being invested to produce high-quality films for basic studies and for potential industrial use. Of primary importance for technological applications is the magnitude and stability of the critical current density $`J_c`$, or more generally, the static and dynamical behavior of the magnetic vortices. Whereas $`J_c`$ as high as 107 A/cm2 have already been reported for MgB2 thin films it has also been shown that this promising value and the technological potential may be challenged by a surprising complexity in the flux dynamics. Magneto-optical (MO) imaging recently revealed that below 10 K the flux penetration in MgB2 films is dominated by large dendritic structures abruptly entering the film, thus dramatically contrasting the smooth and regular flux dynamics normally found in type-II superconductors. This observation is consistent with the large fluctuations reported in the magnetization for various MgB2 samples. The abrupt drops appearing in the magnetization curve, much like at flux jumps only in very large numbers, give an effective 50% suppression of the apparent $`J_c`$.

Abrupt flux dynamics with dendritic penetration patterns has been observed earlier also in films of Nb and YBCO. In all the previous studies the dendrites or magnetization noise have been the result of a gradual increase or decrease of the applied magnetic field (or a short laser pulse perturbing the remanent state). It is therefore still an open question – one of vital importance for many practical applications – whether flux dendrites will be formed also in response to a transport current. In this work it is shown that in MgB2 thin films, a pulse of transport current indeed can trigger the dendritic instability, and we use MO imaging to reveal details of this behavior.

A film of MgB2 was fabricated on ($`1\overline{1}02`$) Al2O3 substrate using pulsed laser deposition. An amorphous B film was first deposited, and then sintered at high temperature in a Mg atmosphere. The film had transition at 39 K with width 0.7 K, and a high degree of $`c`$-axis alignment perpendicular to the film plane. The film thickness was 300 nm, and its lateral dimensions were 3$`\times `$10 mm2. The strip was at both ends equipped with a contact pad of size 3$`\times `$1.5 mm2 fabricated by sputtering of Au through a mask. Copper wires were attached to the pads using silver paste. The MO imaging was performed with the sample mounted on the cold finger in an Oxford Microstat–He optical cryostat. The MO sensing element, a plate of in-plane magnetization ferrite garnet film, was placed on the MgB2 film covering the area between the contacts. The pulses of transport current with rise time of 0.1 $`\mu `$s were sent through the MgB2 film using a power transistor controlled by a pulse generator. The MO images were recorded immediately after the current pulse, the camera exposure time being 250 ms. The time development of the $`IV`$ characteristics during the pulse measured by a two-point scheme was monitored via a TDS-420 oscilloscope.

First, MO images were taken of the MgB2 film after zero-field-cooling (ZFC) to 3.5 K and applying a perpendicular magnetic field. Results of this reference experiment (with no transport current) are shown in Fig. 1. In the images both for increasing (a) and decreasing field (b) one clearly sees that the flux penetration picture is dominated by dendritic patterns. The dendrites are formed instantaneously in response to a slowly varying applied field. Like in the thicker MgB2 film studied in Ref. , also this sample has a threshold temperature of 10 K, above which the flux penetration proceeds in a conventional way with a gradual advancement of flux front.

Shown in Fig. 2(a,b) are flux distributions recorded at 12.5 K, as the applied field cycled from 0 to 70 mT, and back to the remanent state. Evidently, no dendrites were formed, and the final state image, (b), shows all the characteristics of a critical remanent state, where the field is maximum at the so-called discontinuity lines. They are seen as the most bright lines in the image, one in the center of the strip and two pointing towards the corners. At these lines the current, which flows in rectangular loops, as indicated by the white arrows, changes its direction discontinuously. Even without dendrites the flux distribution is not smooth, especially near the edges, see also Fig. 1(a), indicating non-uniformity of the film.

With the sample prepared in the remanent state, (b), a pulse of transport current was applied. As a result, dendritic flux structures of the type seen in (c) and (d) were formed. In (c) and (d) the current pulse was passed through the strip in opposite directions, as indicated by black arrows. Interestingly, one observes that the dendrites always develop from the side of the film where the transport current adds up with the remanent state persistent current. Furthermore, the dendrites tend to expand to the middle of the film, in agreement with numerical simulations. This contrasts the virgin penetration seen in Fig. 1(a), where the dendrites develop more freely into the Meissner state region.

Similar experiments performed at different temperatures showed that a current pulse will trigger flux dendrites up to 19 K, but never above. This threshold temperature is approximately twice that of dendrite formation induced by a slowly changing ($``$ 0.01 T/s) applied field. The difference in threshold temperature can be qualitatively explained by the sensitivity of thermo-magnetic instability to the field (or current) ramp rate that determines magnitude of the electric field in the sample. The dendrite pattern was insensitive to a change in the pulse duration in the range 1.5 – 15 $`\mu `$s if the rise time was kept constant. This also suggests that all dendrites are formed during the rise time.

In another series of experiments the height of the current pulse, $`I`$, was varied. By applying the pulse to the same remanent state, Fig. 2(b), we found that dendrites occur only when $`I`$ exceeds a certain magnitude. In the present case, the threshold current equals 7.0 A at $`T=12.5`$ K, and depends only weakly on temperature. We also observe that the more $`I`$ exceeds the threshold, the larger the dendritic structure becomes in size. Eventually, for sufficiently large $`I`$ two dendrites are formed during the pulse. The full scenario is illustrated in Fig. 3,

which shows by solid circles the area covered by the dendrite(s). The vertical arrows indicate at which $`I`$ the number of dendrites is incremented from zero (below 7 A) to three (at 17 A). For very large currents heating near the contact pads affects the results. It was detected for $`I>15`$ A by the appearance of a non-stationary voltage across the sample which increased during the pulse, see Fig. 3. In addition, flux-free areas around the contacts were seen in the MO images.

Current-induced dendrites were formed not only by starting from the critical remanent state, but also from critical states like the one in Fig. 2(a). The dendrites develop then from the opposite side of the strip since the flow of screening currents is here reversed compared to the remanent state. Hence, all these results show consistently that the dendritic instability is triggered when the transport and the critical-state shielding current flow in the same direction. All attempts to create dendrites by passing a current pulse through a virgin ZFC film failed for $`I`$ up to 30 A, even at the lowest temperature of 3.5 K.

It is believed that the dendritic avalanche behavior results from a thermo-magnetic instability in the superconductor. Vortex dynamics simulations and experiments on MgB2 films under different thermal conditions suggest that the instability stems from the local heating produced by flux motion. The heating will facilitate flux motion nearby, which in turn can lead to a large-scale avalanche invasion of depinned flux lines. The same mechanism is responsible for flux jumps found in the Tesla-range of the magnetization loop of bulk materials. The present investigation is the first one to demonstrate that the dendritic type of the thermo-magnetic instability can be triggered by a transport current. Moreover, films of MgB2 are found more susceptible

for the current-induced instability than for instability induced by varying the applied field. This shows that careful design of thermal stabilization may become essential in order to make viable devices out of MgB2 films.

The work was financially supported by The Norwegian Research Council, and by the Ministry of Science and Technology of Korea through the Creative Research Initiative Program.