text:
# 1 The gap-shape function 𝜒ᴳ_F(𝑋) calculated in the model of pairing with the Gogny force in the slab of nuclear matter within the hard wall potential for 𝜀_F=42 MeV (panel a) and 𝜀_F=46 MeV (panel b). The half-width 𝐿 of the slab is given (in fm) by the numbers close to the curves.

Surface behaviour of the pairing gap
in a slab of nuclear matter

M. Baldo1, M. Farine2, U. Lombardo3,4, E. E. Saperstein5, P. Schuck6,

and M. V. Zverev5

1INFN, Sezione di Catania, 57 Corso Italia, I-95129 Catania, Italy
2Ecole Navale, Lanvéoc-Poulmic, 29240 Brest-Naval, France
3INFN-LNS, 44 Via S.-Sofia, I-95123 Catania, Italy
4 Dipartimento di Fisica, 57 Corso Italia, I-95129 Catania, Italy
5 Kurchatov Institute, 123182, Moscow, Russia
6Insitut de Physique Nucléaire, IN2P3-CNRS,
Université Paris-Sud, F-91406 Orsay-Cédex, France

Abstract

The surface behavior of the pairing gap previously studied for semi-infinite nuclear matter is analyzed in the slab geometry. The gap-shape function is calculated in two cases: (a) pairing with the Gogny force in a hard-wall potential and (b) pairing with the separable Paris interaction in a Saxon-Woods mean-field potential. It is shown that the surface features are preserved in the case of slab geometry, being almost independent of the width of the slab. It is also demonstrated that the surface enhancement is strengthened as the absolute value of chemical potential $`|\mu |`$ decreases which simulates the approach to the nucleon drip line.

1. Introduction

Recently, the surface behavior of the pairing gap in the $`{}_{}{}^{1}S_{0}^{}`$-channel in semi-infinite nuclear matter was investigated independently from two quite different approaches , . A rather sophisticated approach was used in  which starts from the microscopic gap equation for semi-infinite nuclear matter with the separable representation  of the Paris potential . The effective pairing interaction $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ adopted in the gap equation was previously found within the Bethe-Goldstone formalism for semi-infinite nuclear matter without any form of local approximation . All the calculations were made for two values of the chemical potential: $`\mu =16`$MeV and $`\mu =8`$MeV. A surface enhancement in the gap $`\mathrm{\Delta }`$ was found, the effect being more pronounced for $`\mu =8`$ MeV.

In  a more simplified model was used in which nuclear matter was embedded in a semi-infinite hard wall potential and the pairing problem was considered in the BCS approximation with the Gogny force. Such a simple approach makes it possible to solve the problem to a great deal analytically and to examine the coordinate dependence of the pairing gap, pairing tensor and correlation energy density in a rather transparent way. A relatively mild surface enhancement of all the quantities under consideration was found. As to $`\mathrm{\Delta }`$, it is of the same size as in .

In this paper we carried out an analogous analysis for a slab of nuclear matter within both approaches with the hope that a direct comparison of results can help to clarify the general features of the phenomenon under consideration. The slab system is much closer to real atomic nuclei than the semi-infinite one and many results can be qualitatively related to them.

The structure of the article is as follows. In Section 2 we extend the model with the Gogny force  to the case of a hard-wall slab potential. Section 3 contains the extension of the model of  with the Paris force to the case of slab geometry. The results obtained in both models are discussed in Section 4.

2. Pairing with the Gogny force in the hard wall slab potential.

Let us consider a slab of nuclear matter embedded in a hard wall potential of thickness $`2L`$ along the $`x`$-direction. Let the origin be at the center of the slab. We start from expanding the gap operator in $`r`$-space $`\mathrm{\Delta }(𝐫,𝐫^{})`$ in terms of the wave functions $`\phi _𝐤(𝐫)`$ of the hard-wall slab potential,

$$\phi _𝐤(𝐫)=\frac{1}{L}\theta (L+x)\theta (Lx)\mathrm{sin}k_x(xL)e^{i𝐤_{}𝐫_{}},$$
(1)
where $`k=\{k_x,𝐤_{}\}`$, the quantum number $`k_x`$ running over the discrete set of eigenvalues $`k_n=\pi n/(2L)`$, $`n=1,2,\mathrm{}`$. Within the usual BCS approximation the expansion reads

$$\mathrm{\Delta }(𝐫_1,𝐫_2)=\underset{𝐤}{}\phi _𝐤(𝐫_1)\phi _𝐤(𝐫_2)\mathrm{\Delta }(𝐤),$$
(2)
where the state $`|𝐤>`$ is time-reversed with respect to the state $`|𝐤>`$. The gap $`\mathrm{\Delta }(𝐤)`$ obeys the BCS equation

$$\mathrm{\Delta }(𝐩)=\underset{𝐤}{}<𝐩,𝐩\left|V\right|𝐤,𝐤>\frac{\mathrm{\Delta }(𝐤)}{2\sqrt{\xi ^2(𝐤)+\mathrm{\Delta }^2(𝐤)}},$$
(3)
where $`\xi (𝐤)=\epsilon (𝐤)\mu `$ is the single-particle energy relative to the chemical potential and

$$<𝐩,𝐩\left|V\right|𝐤,𝐤>=d𝐫_1d𝐫_2\phi _𝐩^{}(𝐫_1)\phi _𝐩^{}(𝐫_2)V(𝐫_1,𝐫_2)\phi _𝐤(𝐫_2)\phi _𝐤(𝐫_1)$$
(4)
are the matrix elements of the pairing interaction $`V(𝐫_1,𝐫_2)`$. For the purpose of performing calculations analytically we used here, as in , the Gogny force D1

$$V(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐})=\underset{c=1}{\overset{2}{}}(W_cB_cH_c+M_c)e^{(𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐})^2/\alpha _c^2}$$
(5)
with the values of the parameters given in .

Using the explicit form of the eigenfunctions (1) in Eq. (4) one obtains

$$<𝐩,𝐩\left|V\right|𝐤,𝐤>=\frac{1}{2L^2}\underset{c=1}{\overset{2}{}}(W_cB_cH_c+M_c)\pi \alpha _c^2e^{\alpha _c^2(p_{}^2k_{}^2)/4}$$

$$\times dx_1dx_2\theta (x_1+L)\theta (Lx_1)\theta (x_2+L)\theta (Lx_2)$$

$$\times \mathrm{sin}\left(k_x(x_1L)\right)\mathrm{sin}\left(k_x(x_2L)\right)\mathrm{sin}\left(p_x(x_1L)\right)\mathrm{sin}\left(p_x(x_2L)\right)e^{(x_1x_2)^2/\alpha _c^2}.$$
(6)
Upon substituting this expression into Eq. (3) and, as the $`s`$-wave pairing is considered, averaging over the angle between vectors $`𝐩`$ and $`𝐤`$, the gap equation can be rewritten as follows

$$\mathrm{\Delta }(p_m,p_{})=\frac{1}{4L^2}\underset{c=1}{\overset{2}{}}(W_cB_cH_c+M_c)\alpha _c^2e^{\alpha _c^2p_{}^2/4}$$

$$\times \underset{n}{}\underset{0}{\overset{\mathrm{}}{}}k_{}dk_{}e^{\alpha _c^2k_{}^2/4}I_0\left(\frac{\alpha _cp_{}k_{}}{2}\right)\frac{\mathrm{\Delta }(k_n,k_{})}{\sqrt{\xi ^2(k_n,k_{})+\mathrm{\Delta }^2(k_n,k_{})}}$$

$$\times dx_1dx_2\theta (x_1+L)\theta (Lx_1)\theta (x_2+L)\theta (Lx_2)$$

$$\times \mathrm{sin}((k_n(x_1L))\mathrm{sin}\left(k_n(x_2L)\right)\mathrm{sin}\left(p_m(x_1L)\right)\mathrm{sin}\left(p_m(x_2L)\right)e^{(x_1x_2)^2/\alpha _c^2},$$
(7)
where $`I_0(z)`$ is the modified Bessel function.

Integrating then over $`x_1`$ and $`x_2`$ and introducing the function

$$g_c(p,k)=\frac{i\sqrt{\pi }\alpha _c}{2L^2(p+k)}\{e^{\alpha _c^2k^2/4}[(e^{2i(p+k)L}+1)\mathrm{erf}\left(\frac{ik\alpha _c}{2}\right)\mathrm{erf}(\frac{ik\alpha _c}{2}\frac{2L}{\alpha _c})$$

$$e^{2i(p+k)L}\mathrm{erf}(\frac{ik\alpha _c}{2}+\frac{2L}{\alpha _c})]+e^{\alpha _c^2p^2/4}[(e^{2i(p+k)L}+1)\mathrm{erf}\left(\frac{ip\alpha _c}{2}\right)$$

$$\mathrm{erf}(\frac{ip\alpha _c}{2}\frac{2L}{\alpha _c})e^{2i(p+k)L}\mathrm{erf}(\frac{ip\alpha _c}{2}+\frac{2L}{\alpha _c})\left]\right\},$$
(8)
we arrive at the gap equation in the following form

$$\mathrm{\Delta }(p_m,p_{})=\frac{1}{2}\underset{c=1}{\overset{2}{}}\alpha _c^2(W_cB_cH_c+M_c)e^{\alpha _c^2p_{}^2/4}\underset{0}{\overset{\mathrm{}}{}}k_{}𝑑k_{}e^{\alpha _c^2k_{}^2/4}I_0\left(\frac{\alpha _cp_{}k_{}}{2}\right)$$

$$\times \underset{n}{}\frac{\mathrm{\Delta }(k_n,k_{})}{\sqrt{\xi ^2(k_n,k_{})+\mathrm{\Delta }^2(k_n,k_{})}}G_c(p_m,k_n),$$
(9)
where

$$G_c(p,k)=\frac{1}{8}\mathrm{Re}\{g_c(p+k,p+k)+g_c(pk,pk)+g_c(p+k,pk)$$

$$+g_c(pk,p+k)2g_c(p+k,p+k)2g_c(p+k,pk)\}.$$
(10)
To reveal spatial behavior of the nonlocal pairing gap operator the Wigner transform of the gap is very useful. It reads

$$\mathrm{\Delta }(𝐑,𝐤)=𝑑𝐬\mathrm{\Delta }(𝐑,𝐬)e^{i\mathrm{𝐤𝐬}},$$
(11)
where $`𝐑=(𝐫_\mathrm{𝟏}+𝐫_\mathrm{𝟐})/2`$ and $`𝐬=𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}`$. In the case of slab geometry, the Wigner transform (11) of the gap operator depends only on $`X`$ which is the component of $`𝐑`$ perpendicular to the surface. Considering $`\mathrm{\Delta }(X,k_x,k_{})`$ only for $`X>0`$, since it is an even function of $`X`$, one can easily obtain the following series:

$$\mathrm{\Delta }(X,k_x,k_{})=\frac{1}{\pi L}\theta (LX)\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{\Delta }(k_n,k_{})[2\mathrm{cos}\left(k_n(XL)\right)\frac{\mathrm{sin}\left(2k_x(XL)\right)}{k_x}$$

$$\frac{\mathrm{sin}\left(2(k_xk_n)(XL)\right)}{k_xk_n}\frac{\mathrm{sin}\left(2(k_x+k_n)(XL)\right)}{k_x+k_n}].$$
(12)
In the bulk, the gap $`\mathrm{\Delta }`$ depends mainly on the total momentum $`k=\sqrt{k_x^2+k_{}^2}`$. Approximately, this is true also for the surface region. Within this approximation, we can treat the gap at the Fermi surface $`\mathrm{\Delta }(X,k_\mathrm{F}(X))`$ as the Wigner transform $`\mathrm{\Delta }(X,k_x,k_{})`$ taken at $`k=k_\mathrm{F}(X)`$, where $`k_\mathrm{F}=\left(3\pi ^2\rho (X)/2\right)^{1/3}`$ is the local Fermi momentum. We define the gap-shape function as the gap at the Fermi surface normalized to one at the center of the slab:

$$\chi _\mathrm{F}^\mathrm{G}(X)=\mathrm{\Delta }(X,k_\mathrm{F}(X))/\mathrm{\Delta }(0,k_\mathrm{F}(0)).$$
(13)
The function $`\chi _\mathrm{F}^\mathrm{G}(X)`$ is drawn in Fig. 1 (a and b) for four different values of the half-width of the slab: $`L=4,6,8`$ and 10 fm. In Fig. 1 a the gap-shape function is calculated for the Fermi momentum $`k_\mathrm{F}=1.4232`$ fm-1 corresponding to the Fermi energy $`\epsilon _\mathrm{F}=42`$ MeV, while in calculation of Fig. 1 b the value of Fermi momentum, $`k_\mathrm{F}=1.4894`$ fm-1, corresponds to $`\epsilon _\mathrm{F}=46`$ MeV.

Let us discuss the salient features of the shape function $`\chi _\mathrm{F}^\mathrm{G}(X)`$. The main observation is that in the slab, almost independent of its thickness, the surface enhancement is not very different from the semi-infinite matter case previously studied . The enhancement of the surface peak is rather moderate, not exceeding 20% - 30%. In fact it is not clear how much of the surface peak comes from Friedel oscillations and how much is a genuine enhancement of $`\mathrm{\Delta }`$. From this result and our previous study  we are therefore inclined to conclude that at least with the Gogny force pairing is only very moderately surface enhanced. As to the dependence of the gap-shape function on the size of the slab, in eq. (13) there is a competition between two effects. First, with diminishing of the size parameter $`L`$, the peak value of the numerator $`\mathrm{\Delta }(XL)`$ becomes higher. Second, the denominator (the gap value in the center) also gets enhanced. As a result, the $`L`$-dependence of the gap-shape function turns out to be rather smooth.

Of course, finite nuclei are rather different from a slab but qualitatively things should be similar. Some results for finite nuclei with the Gogny force do exist  in the $`Sn`$-region but no definite conclusions about the surface features of $`\mathrm{\Delta }`$ can be drawn from this limited number of values. In fact, in finite nuclei like the tin isotopes, we should look at a great number of nuclei in the isotopic chain because the behavior of $`\mathrm{\Delta }(r)`$ can fluctuate a great deal passing through the open neutron shell in question.

A detailed study of several long isotopes chains  with density dependent effective pairing forces confirmed that, in the framework of a phenomenological approach, it is rather difficult to distinguish between two opposite possibilities, the volume pairing and the pairing with pronounced surface enhancement. For this purpose, the analysis of some binding energy characteristics, such as separation energies, should be accompanied by the study of the variations of nuclear radii along the chain. The odd-even staggering phenomenon is especially sensitive to the coordinate dependence of $`\mathrm{\Delta }`$. The microscopic calculation of $`\mathrm{\Delta }`$ should help to solve this problem.

Of course, the conclusion may strongly depend on the employed pairing force and below we will investigate the surface behavior of pairing with a separable version of the Paris force.

3. Pairing with the Paris force in a slab of nuclear matter

within a Saxon-Woods potential.

Now we shall adopt the more realistic one-dimensional Saxon-Woods potential $`U(X)`$ for a slab with the width of $`2L`$ symmetrical around the origin $`X=0`$:

$$U(X)=\frac{U_0}{1+\mathrm{exp}((XL)/d)+\mathrm{exp}((X+L)/d)}.$$
(14)
Here $`U_0`$ is the potential depth in the central region and $`d`$ is the diffuseness parameter (to be more exact, the maximum potential depth is $`U(0)=U_0/\left(1+2\mathrm{exp}(L/d)\right)`$). The two parameters ($`U_0=50`$MeV and $`d=0.65`$fm ) are taken to be close to those of real atomic nuclei. The half-width parameter $`L`$ will be changed to examine the size dependence of the effect under consideration.

To avoid a rather cumbersome resolution of the Bogolyubov equations for the nonlocal gap , as in , we use a powerful method  (we refer to it as KKC) of solving the gap equation for the case of a nonlocal interaction. This method was originally suggested for infinite matter where the gap $`\mathrm{\Delta }`$ can be represented as a product $`\mathrm{\Delta }(p)=\mathrm{\Delta }_\mathrm{F}\chi (p)`$ of the constant $`\mathrm{\Delta }_\mathrm{F}=\mathrm{\Delta }(p_\mathrm{F})`$ and the ”gap-shape” function $`\chi (p)`$ normalized to $`\chi (p_\mathrm{F})=1`$. Basically the KKC method is a transformation of the gap equation to a set of two coupled equations: an integral equation for $`\chi (p)`$, which is almost independent of the value of $`\mathrm{\Delta }_\mathrm{F}`$, and an algebraic equation for the value $`\mathrm{\Delta }_\mathrm{F}`$. This significantly simplifies the solution of the gap equation in infinite matter. In  the extension of the KKC method to nonzero temperatures  was used in the case of semi-infinite nuclear matter where the spatial dependence of the gap-shape function was also taken into account. Within the KKC method for infinite matter, the temperature and momentum parts of the gap function can be factorized as follows: $`\mathrm{\Delta }(𝐩,T)=\mathrm{\Delta }_\mathrm{F}(T)\chi (𝐩)`$. In  it was supposed that a similar separation of the temperature factor can be made for the semi-infinite system:

$$\mathrm{\Delta }(x_1,x_2,k_{}^2;T)=\mathrm{\Delta }_\mathrm{F}(T)\chi (x_1,x_2;k_{}^2).$$
(15)
An additional advantage from using the KKC method in this case comes from the possibility of finding the normalization factor $`\mathrm{\Delta }_\mathrm{F}(T)`$ by solving the gap equation in infinite matter.

In this paper we use the same ansatz (15) for the slab geometry. Unfortunately, in this case no direct relation to infinite nuclear matter exists and there is no evident way to find the normalization factor without solving the Bogolyubov equations. However, since the gap-shape function is of main importance for an analysis of the surface enhancement of the gap, we do not calculate here the normalization factor postponing the solution of the Bogolyubov equation to a forthcoming publication.

As far as in the case of the slab geometry all the equations are very similar to those for semi-infinite matter, we write down explicitly only those which are necessary for explaining our calculations and refer the reader to  for details. In symbolic notation, the gap equation has the form :

$$\mathrm{\Delta }(T)=𝒱_{\mathrm{eff}}^\mathrm{p}A_0^s(T)\mathrm{\Delta }(T)$$
(16)
where $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ is the effective pairing interaction acting in the model space $`S_0`$ in which the superfluid two-particle propagator $`A_0^s`$ is defined.

The separable $`3\times 3`$ form  of the Paris potential  is used

$$V(𝐤,𝐤^{})=\underset{ij}{}\lambda _{ij}g_i(k^2)g_j(k^2),$$
(17)
where $`𝐤`$ and $`𝐤^{}`$ are the relative momenta before and after scattering. The effective interaction has a similar separable form which, in notation of , is as follows:

$$𝒱_{\mathrm{eff}}^\mathrm{p}(x_1,x_2,x_3,x_4,k_{}^2,k_{}^2;E)=\underset{ij}{}\mathrm{\Lambda }_{ij}(X,X^{};E)g_i(k_{}^2,x)g_j(k_{}^2,x^{}),$$
(18)
where $`E=2\mu `$, $`X=(x_1+x_2)/2`$, $`x=x_1x_2`$, etc., and $`g_i(k_{}^2,x)`$ stands for the inverse Fourier transform in the $`x`$-direction of the form factor $`g_i(k_{}^2+k_x^2)`$. The gap-shape factor can be also written as

$$\chi (x_1,x_2;k_{}^2)=\underset{i}{}\chi _i(X)g_i(k_{}^2,x).$$
(19)
After substituting Eqs. (15), (17) – (19) into Eq. (16) at $`T=T_c`$ we obtain the following equation for the components $`\chi _i`$:

$$\chi _i(X)=\underset{lm}{}𝑑X_1𝑑X_2\mathrm{\Lambda }_{il}(X,X_1;E)B_{lm}(X_1,X_2,E;T_c)\chi _m(X_2),$$
(20)
where

$`B_{lm}(X,X^{},E;T)={\displaystyle \underset{n_1n_2}{}}{\displaystyle }{\displaystyle \frac{d𝐤_{}}{(2\pi )^2}}{\displaystyle \frac{1N_{\lambda _1}(T)N_{\lambda _2}(T)}{E\epsilon _{\lambda _1}\epsilon _{\lambda _2}}}\times `$ (21)
$`\times G_{n_1n_2}^l(k_{}^2,X)G_{n_1n_2}^m(k_{}^2,X^{}),`$
$$G_{n_1n_2}^l(k_{}^2,X)=y_{n_1}(X+\frac{x}{2})g_l(k_{}^2,x)y_{n_2}(X\frac{x}{2})𝑑x,$$
(22)
$$N_\lambda (T)=\left(1+\mathrm{exp}\frac{\epsilon _\lambda \mu }{T}\right)^1,$$
(23)
In Eqs. (21) – (23), $`\lambda =(n,𝐤_{})`$, $`\epsilon _\lambda =\epsilon _n+k_{}^2/2m`$, $`\epsilon _n`$ and $`y_n(x)`$ stand for the energies and wave functions, respectively, of the 1-dimensional Schrödinger equation with the potential (14).

Following the recipe of Ref.  for the effective interaction $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ the propagator $`A_0^s`$ embodies all the single-particle states with negative energies only. Thus, the summation over $`n_1,n_2`$ and the integration over $`𝐤_{}`$ are limited by the conditions: $`\epsilon _\lambda <0,\epsilon _\lambda ^{}<0`$.

The coefficients $`\mathrm{\Lambda }_{ij}`$ obey the set of integral equations

$`\mathrm{\Lambda }_{ij}(X_{12},X_{34};E)=\lambda _{ij}\delta (X_{12}X_{34})+`$ (24)
$`+{\displaystyle \underset{lm}{}}\lambda _{il}{\displaystyle 𝑑X_{56}B_{lm}(X_{12},X_{56};E)\mathrm{\Lambda }_{mj}(X_{56},X_{34};E)},`$
where $`B_{lm}`$ are defined by an expression similar to Eq. (21), but without the temperature factor, including the states $`\lambda _1,\lambda _2`$ from the complementary subspace.

To simplify the calculations we used the local potential approximation (LPA) which has turned out to be accurate for semi-infinite nuclear matter  and for nuclear slabs . The LPA prescription consists in using for the 2-particle propagator $`B`$ of the complementary space the local momentum approximation, very similar to the Local Density Approximation (LDA) , for each particle separately: $`\epsilon _np_x^2/2m+V(X)`$, $`\epsilon _n^{}(p_x^{})^2/2m+V(X)`$. This type of approach, where the individual particles are treated in semiclassical approximation was used in  for examining the response function. This approximation has been shown to be very accurate, if one is not interested in fine details. It should be stressed that the LPA is only applied to the equation for the effective interaction $`𝒱_{\mathrm{eff}}^\mathrm{p}`$, while no local approximation is used in the renormalized gap equation (16). That is, the local approximation is used only for two-particle states belonging to the complementary space, for which the corresponding energy denominators in eq.(21) are large. Therefore the individual contribution of each state is negligible, and only the sum of a number of such contributions is important. For such a sum the semiclassical and local approximations are expected to be accurate. In this respect LPA is different from the standard LDA, since in the latter the local approximation is used for all two-particle states.

Within the LPA, the exact values of $`B_{lm}(X_1,X_2;E)`$ are replaced by the set of $`B_{lm}^{\mathrm{inf}}(t,E;U[X])`$ calculated for infinite nuclear matter put into the homogeneous potential well of the depth $`U[X]`$. Here $`t=X_1X_2`$, $`X=(X_1+X_2)/2`$ are the relative and average values of the CM coordinates. In the first step of the LPA procedure we calculate the set of vectors $`B_{lm}^{\mathrm{inf}}(t,E=2\mu ;U_i)(V_i=\delta V(i1))`$ for a fixed value of the chemical potential $`\mu `$. In the second step, for each value of $`(X,t)`$, we find $`B_{lm}^{\mathrm{LPA}}(X_1,X_2)`$ interpolating the values of $`B_{lm}^{\mathrm{inf}}(t;V_i)`$ by the values of $`V_i`$ nearest to $`V(X)`$. After substitution of $`B_{lm}^{\mathrm{LPA}}(X_1,X_2)`$ into Eq. (24) we find the LPA prescription $`\mathrm{\Lambda }_{lm}^{\mathrm{LPA}}(X_1,X_2)`$ for the effective interaction which should be substituted into the homogeneous Eq.(20) for the gap-shape function. If the critical temperature $`T_c`$ was known, this equation could be solved directly. To find $`T_c`$ a more general integral equation must be considered

$$\chi _i(X)=\lambda (T)\underset{lm}{}𝑑X_1K_{il}(X,X_1)\chi _m(X_1),$$
(25)
which involves the eigenvalue $`\lambda (T)`$. Here the abbreviation $`K=\mathrm{\Lambda }B^0`$ is introduced for the kernel. The critical temperature can be found from the evident condition $`\lambda (T_c)=1`$.

The entire calculation scheme is similar to that for semi-infinite matter except for some details. First, in the slab case we are dealing with the discrete spectrum $`\epsilon _n`$ in Eq. (21). Second, due to the obvious reflection symmetry of the slab system in the $`x`$-direction, all the integral equations under consideration can be readily reduced to a form including positive $`X`$ only. Just as in , instead of the direct solution of Eq. (25) in the coordinate space, we use the Fourier expansion within the symmetrical interval $`(L_0,L_0)`$,

$$\chi _i(X)=\underset{n}{}\chi _i^nf_n(X),$$
(26)
where only the even functions must be retained, $`f_n(x)=\mathrm{cos}(\pi n(XX_c)/L_0)`$. The kernels $`K_{ij}(X,X^{})`$ of Eq. (25) are also expanded in a double Fourier series. Finally we arrive at a set of homogeneous linear equations for the coefficients $`\chi _i^n`$ :

$$\chi _i^n=\underset{j=1}{\overset{3}{}}\underset{n^{}=1}{\overset{N}{}}K_{ij}^{nn^{}}\chi _j^n^{},$$
(27)
which can be solved by standard numerical methods. Then the components of the gap-shape function are found from Eq. (26).

Instead of analyzing the separate components $`\chi _i(X)`$ for the separable form (17) of Paris force it is more useful to display the local form of the gap-shape function which enters the matrix elements of the gap for states nearby the Fermi surface:

$$\chi _\mathrm{F}^\mathrm{P}(X)=\underset{i}{}\chi _i(X)g_i\left(k^2=k_\mathrm{F}^2(X)\right),$$
(28)
where $`k_\mathrm{F}(X)=\sqrt{2m(\mu U(X))}`$ is the local Fermi momentum ($`k_\mathrm{F}(X)=0`$ for $`\mu U(X)<0`$). We calculated the gap-shape function $`\chi _\mathrm{F}^\mathrm{P}(X)`$ for the same values of the half-width of the slab, $`L=4,6,8`$ and 10 fm, as in the case of the hard wall potential and for two values of the chemical potential $`\mu =8`$ MeV and $`4`$ MeV. As the depth of the Saxon-Woods potential was taken $`U_0=50`$ MeV, these two values of $`\mu `$ correspond just to the same values of the Fermi energy as in the hard-wall case. The results are shown in Fig. 2 a and b.

One observes a surface bump which is much more pronounced than in the previous model with the Gogny force and the hard wall potential. The enhancement is now around 80% - 100% and, as for the Gogny force, it is quite similar to the semi-infinite matter case . Two values of the chemical potential $`\mu =8`$ MeV and $`\mu =4`$ MeV have been chosen which in account of the depth of the Saxon-Woods potential of $`V_0=50`$ MeV correspond precisely to the Fermi energies $`\epsilon _\mathrm{F}=42`$ MeV and $`\epsilon _\mathrm{F}=46`$ MeV, respectively, of the previous model. Going from $`\mu =8`$ MeV to $`\mu =4`$ MeV a rather important increase of the enhancement of the order of 30% is observed which is much larger than in the case of the hard-wall potential with the Gogny force. On the other hand, there is little variation with the thickness of the slab, the maximum of $`\mathrm{\Delta }`$ being a couple of percent larger for small slab size. This situation is analogous to the previous model.

The surface effect in $`\mathrm{\Delta }`$ for the Paris force can be qualitatively explained from the properties of the effective interaction $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ which were analyzed in . There it was found that $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ undergoes a sharp variation in the surface region, from almost zero in the bulk to very strong attraction in vacuum. In the asymptotic region, the latter coincides with the off-shell $`T`$-matrix of free $`NN`$-scattering $`T(E=2\mu )`$ which exhibits a resonant behavior at small $`E`$. The strong surface attraction and the sharp variation in the surface region are mostly responsible for the surface effect of the gap. The $`\mu `$-dependence of the surface effect can be explained by the increase of the jump $`\delta 𝒱_{\mathrm{eff}}^\mathrm{p}`$ from inside to outside as $`|\mu |`$ is decreasing. There are two reasons for such an increase. The first one is the $`k^2`$-dependence of the form factors in Eqs.(17), (18) leading to a reduction of $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ in the inner region with increasing values of $`k_F`$. The second one is a pole-like behavior of $`T(E)`$ at small $`E`$ which results in an increase of $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ with decreasing $`|\mu |`$ in the exterior, due to the approach to the virtual pole. One sees that both effects work in the same direction resulting in strengthening the surface effect at small values of $`|\mu |`$.

4. Discussion, and conclusions.

In this work we continued our effort to understand the surface behavior of the nuclear gap. Previous investigations considered semi-infinite nuclear matter embedded in i) a hard wall potential with the Gogny force  and ii) a potential of Saxon-Woods shape with the separable version of the Paris force . In both cases a surface enhancement was found but which is relatively modest in view of what one could expect from LDA. In the present study we addressed the question whether finite size effects can strongly alter this situation and repeated the former calculation  in a slab configuration.

In the first case of pairing with the Gogny force (D1) within the slab of nuclear matter, the hard wall potential allows to perform most part of calculations analytically. The second rather sophisticated case demanded a lot of numerical effort due to the use of the Paris interaction and the Saxon-Woods shape of the mean field potential. In both cases a noticeable surface effect for the pairing gap was obtained of the same order of magnitude as it was previously found in semi-infinite nuclear matter. The shape of the gap in coordinate space turned out to be qualitatively similar in both cases, with a significant surface enhancement. For the value of the chemical potential $`\mu =8`$MeV which simulates stable atomic nuclei, the enhancement is of the order of 30% for the first model and is almost 100% for the second one. A general feature of both models is the rather smooth dependence of the enhancement on the slab thickness which is approximately 10% in the first case and only 5% in the second one. In both cases, a $`\mu `$-dependence of the surface effect is found: the enhancement coefficient increases as the absolute value of the chemical potential $`|\mu |`$ decreases. The latter effect is more pronounced for the Paris interaction and Saxon-Woods potential reaching 30% with diminishing $`|\mu |`$ from 8 MeV to 4 MeV. For the Gogny force and box potential, the corresponding $`\mu `$-effect is approximately 10%.

To understand in detail the possible common physical origin of this surface enhancement as well as the mentioned differences, it is instructive to consider the gap equation in the form (16) in which $`\mathrm{\Delta }`$ is expressed in terms of the effective interaction $`𝒱_{\mathrm{eff}}^\mathrm{p}`$ given by Eq. (18). Properties of the effective interaction generated by the Paris force were analyzed in , where it was found to undergo a sharp variation in the surface region. It is this sharp variation that is mostly responsible for the surface effect in the gap. It is natural that in the case of the hard wall box potential the influence of the surface interaction is smaller which makes the surface enhancement weaker.

The $`\mu `$-dependence of the surface effect is explained by the increase of the jump $`\delta 𝒱_{\mathrm{eff}}^\mathrm{p}=𝒱_{\mathrm{eff}}^{\mathrm{in}}𝒱_{\mathrm{eff}}^{\mathrm{ex}}`$ with decreasing $`|\mu |`$. This increase is caused by two reasons. The first one is the strong $`k^2`$-dependence of the form factors in Eq.(17) leading to a reduction of $`𝒱_{\mathrm{eff}}^{\mathrm{in}}`$ for the larger values of $`k_\mathrm{F}`$ in the bulk. The second one is an increase of $`𝒱_{\mathrm{eff}}^{\mathrm{ex}}`$ with decreasing $`|\mu |`$ caused by the pole-like behavior of the $`T`$-matrix at small $`E`$. Thus, both reasons jointly work towards making the surface effect stronger at small values of $`|\mu |`$. Qualitatively, the two reasons should work also for Gogny force, but the hard wall potential strongly suppresses the second reason, diminishing the surface effect itself and its $`\mu `$-dependence as well. Finally, it should be mentioned that a large value of the coherence length of pairing which is comparable with the size of the slab makes all the effects under considerations rather smooth. In particular, it results in the weak dependence of the surface enhancement of the pairing gap on the slab thickness.

It is worth to point out that the mechanism of the surface enhancement in $`\mathrm{\Delta }`$ considered in this paper is different from the one of ref. . The latter is related to a contribution to the effective pairing interaction of the virtual exchange by collective surface vibrations. In a consistent description of the surface behavior of the pairing gap in nuclei these two effects should be considered on equal footing.

This research was partially supported by Grants No. 00-15-96590 and No. 00-02-17319 from the Russian Foundation for Basic Research. Two of the authors (E.E.S. and M.V.Z.) thank INFN (Sezione di Catania and LNS) and Catania University for hospitality during their stay in Catania.


format: text

text:
# Shot noise in tunneling transport through molecules and quantum dots

## I Introduction

Single-molecule electronics provides an exciting possibility for further miniaturization of electronic devices. Nonlinear electrical transport through single molecules deposited between mechanically controlled break junctions has been measured in several experiments. reed-etal ; kergueris-etal ; reichert-etal ; platin-complex  Furthermore, recently observed diode behaviormetzger ; vilan  memory effects and negative differential conductancechen-reed-cp  in devices which make use of molecular films have attracted much interest.

A crucial problem for understanding and designing molecular electronic devices is the lack of control of the metal-molecule interfaces. In molecular-film experiments, deposition of the covering electrode is very difficult and can easily destroy the film or induce short cuts. Once the molecules are deposited, no direct study and manipulation of the contact is possible. Even in a break junction experiment one has very limited influence on how a molecule bonds to the two electrodes. Information on the coupling strengths $`\mathrm{\Gamma }_\mathrm{L},\mathrm{\Gamma }_\mathrm{R}`$ between the molecule and electrodes has to rely on transport measurements. From the measurement of the current $`I`$, however, only a combination of $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$ can be deduced, whereas the asymmetry ratio $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}`$ remains undetermined.

In addition, it is impossible to judge from the $`IV`$ curve alone, whether charging effects such as Coulomb blockade are relevant to the transport. This is because any observed set of current steps expected for an $`IV`$ curve displaying Coulomb blockade can be mimicked by transport through a set of noninteracting levels. As so far there is no quantitative theory for describing hybrid metal-molecule functions, it remains a matter of opinion, whether such a fit is reasonable for a given experiment.

Further information can be gained from the shot noise $`S`$, allowing, in principle, to determine the ratio of the couplings. Furthermore, the combination of current and shot-noise measurement provides information about the presence and size of Coulomb charging energy. For transport through a noninteracting levelKhlus ; Lesovik ; blanter  the so-called Fano factor $`F=S/2eI`$ is known to be $`(\mathrm{\Gamma }_\mathrm{L}^2+\mathrm{\Gamma }_\mathrm{R}^2)/(\mathrm{\Gamma }_\mathrm{L}+\mathrm{\Gamma }_\mathrm{R})^2`$. However, transport through weakly coupled molecules can not be modeled by noninteracting molecular levels, since Coulomb blockade effects become important.hsw ; hettler\_prl ; cocomplex

Some results have been obtained before for shot noise through interacting systems.hershfield ; nazarov ; loss ; bagrets ; haug  However, the focus in these works was on transport through impurity states or quantum dots, where some questions pertinent to molecules are irrelevant.

In the present paper, we reformulate the “orthodox” theory of shot noiseblanter  in a way that allows a systematic perturbation expansion in the coupling strength. We study in detail the shot noise for a molecule accommodating one level, described by the single-level Anderson model, in the lowest-order perturbation theory. We perform a complete classification of all transport regimes, and evaluate the corresponding Fano factors, recovering known expressions for the cases previously studied in the literature. The Fano factor turns out to depend in a nonmonotonic way on the asymmetry ratio $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}`$. The $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}`$ dependence of the Fano factor can be used to identify the different transport regimes. Our results hold for transport through single-level quantum dots as well. In contrast to molecular devices, the contacts to leads in the semiconductor quantum dots are usually well controllable.

## II The model

The simplest model for electron transport through a molecule with Coulomb interaction is the Anderson-impurity model described by the Hamiltonian $`\widehat{H}=\widehat{H}_\mathrm{L}+\widehat{H}_\mathrm{R}+\widehat{H}_{\mathrm{Mol}}+\widehat{H}_{\mathrm{T},\mathrm{L}}+\widehat{H}_{\mathrm{T},\mathrm{R}}`$ with

$`\widehat{H}_r={\displaystyle \underset{k\sigma }{}}ϵ_{k\sigma r}a_{k\sigma r}^{}a_{k\sigma r},`$ (1)
$`\widehat{H}_{\mathrm{Mol}}={\displaystyle \underset{\sigma }{}}ϵ_\sigma c_\sigma ^{}c_\sigma +Un_{}n_{},`$ (2)
$`\widehat{H}_{\mathrm{T},r}={\displaystyle \underset{k\sigma }{}}(t_ra_{k\sigma r}^{}c_\sigma +h.c.)`$ (3)
and $`r=\mathrm{L},\mathrm{R}`$. Here, $`\widehat{H}_\mathrm{L}`$ and $`\widehat{H}_\mathrm{R}`$ model the left and right electrodes with noninteracting electrons, $`\widehat{H}_{\mathrm{Mol}}`$ describes the molecule with one relevant (spin-dependent) molecular level of energy $`ϵ_\sigma `$, and the Coulomb interaction $`U`$ on the molecule ($`n_{}`$ and $`n_{}`$ are the number operators for electrons with corresponding spin). Tunneling between leads and molecule level is modeled by $`\widehat{H}_{\mathrm{T},\mathrm{L}}`$ and $`\widehat{H}_{\mathrm{T},\mathrm{R}}`$. The coupling strength is characterized by the intrinsic linewidth $`\mathrm{\Gamma }_r=2\pi \rho _e|t_r|^2`$, where $`\rho _e`$ is the (constant) density of states of the leads. Furthermore, we define the total linewidth $`\mathrm{\Gamma }=\mathrm{\Gamma }_\mathrm{L}+\mathrm{\Gamma }_\mathrm{R}`$. The Fermi operators $`a_{k\sigma r}^{}(a_{k\sigma r})`$ and $`c_\sigma ^{}(c_\sigma )`$ create (annihilate) electrons in the electrodes and the molecule.

We are interested in transport through the molecule, in particular in the current $`I`$ and the (zero-frequency) current noise $`S`$. They are related to the current operator $`\widehat{I}=(\widehat{I}_\mathrm{R}\widehat{I}_\mathrm{L})/2`$, with $`\widehat{I}_r=i(e/\mathrm{})_{k\sigma }(t_ra_{k\sigma r}^{}c_\sigma h.c.)`$ being the current operator for electrons tunneling into lead $`r`$, by $`I=\widehat{I}`$, and

$$S=_{\mathrm{}}^{\mathrm{}}𝑑t\delta \widehat{I}(t)\delta \widehat{I}(0)+\delta \widehat{I}(0)\delta \widehat{I}(t),$$
(4)
where $`\delta \widehat{I}(t)=\widehat{I}(t)\widehat{I}`$.

## III Diagrammatic Technique

A diagrammatic technique to perform a systematic perturbation expansion of the current through localized levels has been developed in Ref. diagrams, . Here, we expand this technique to address the perturbation expansion of the current noise as well. After this general formulation, we specify the results for lowest order in $`\mathrm{\Gamma }`$ and discuss implications for transport through molecules in this regime.

All transport properties are governed by the nonequilibrium time evolution of the density matrix. The leads are described as reservoirs of noninteracting electrons, which remain in thermal equilibrium. A finite transport voltage enters the difference of the electrochemical potentials $`\mu _r`$ for the left and right leads, $`r=\mathrm{L},\mathrm{R}`$. These degrees of freedom can be integrated out making use of Wick’s theorem, such that the Fermion operators are contracted in pairs. As a result, we end up with a reduced density matrix for the molecule’s degrees of freedom only, labeled by $`\chi `$. In our case, four molecule states are possible: the level is empty, $`\chi =0`$, singly occupied with either spin up or down, $`\chi =`$ or $``$, or doubly occupied $`\chi =d`$. The time evolution of the reduced density matrix, described by the propagator $`\mathrm{\Pi }_{\chi ^{}\chi }(t^{},t)`$ for the propagation from state $`\chi `$ at time $`t`$ to state $`\chi ^{}`$ at time $`t^{}`$, can be visualized by diagrams.

An example is shown in Fig. 1. Vertices representing tunneling, $`\widehat{H}_\mathrm{T}`$, are connected in pairs by tunneling lines. The full propagation is expressed as a sequence of irreducible blocks $`W_{\chi ^{}\chi }(t^{},t)`$ containing one or more tunneling lines.

They are associated with transitions from state $`\chi `$ at time $`t`$ to state $`\chi ^{}`$ at time $`t^{}`$ (see Fig. 2). This leads to the Dyson equation for the propagator

$$𝚷(0,t)=\mathrm{𝟏}+_t^0𝑑t_2_t^{t_2}𝑑t_1𝐖(t_2,t_1)𝚷(t_1,t),$$
(5)
where the boldface indicates matrix notation related to the molecular state labels. To calculate transport properties in the stationary limit, we start at time $`t\mathrm{}`$, at which the system has a given diagonal but otherwise arbitrary distribution of probabilities $`p_\chi ^{\mathrm{init}}`$ to be in state $`\chi `$, comprised in the vector $`𝐩^{\mathrm{init}}`$. The stationary probability distribution $`𝐩^{\mathrm{st}}=𝚷(0,\mathrm{})𝐩^{\mathrm{init}}`$, however, does not depend on $`𝐩^{\mathrm{init}}`$, i.e., $`\mathrm{\Pi }_{\chi ^{}\chi }(0,\mathrm{})=p_\chi ^{}^{\mathrm{st}}`$, independent of $`\chi `$. To determine $`𝐩^{\mathrm{st}}`$ we use the Dyson equation Eq. (5) in the limit $`t\mathrm{}`$, multiply $`𝐩^{\mathrm{init}}`$ from the right, and define $`𝐖=_{\mathrm{}}^0𝑑t\mathrm{}𝐖(0,t)`$ to get $`𝐩^{\mathrm{st}}𝐩^{\mathrm{init}}=\frac{1}{\eta }\mathrm{𝐖𝐩}^{\mathrm{st}}`$, where $`\eta 0^+`$ comes from a convergence factor $`\mathrm{exp}(\eta |t_2|/\mathrm{})`$ in the time integral. This leads to $`\mathrm{𝐖𝐩}^{\mathrm{st}}=\mathrm{𝟎}`$ together with the normalization $`𝐞^T𝐩^{\mathrm{st}}=1`$, where the vector $`𝐞`$ is given by $`e_\chi =1`$ for all $`\chi `$. The matrix $`𝐖`$ has zero determinant and hence cannot be inverted (this can be seen from the sum rule $`𝐞^T𝐖=\mathrm{𝟎}`$). To write down a single matrix equation which determines $`𝐩^{\mathrm{st}}`$ we make use of the normalization of probabilities and introduce the matrix $`\stackrel{~}{𝐖}`$, which is identical to $`𝐖`$ but with one (arbitrarily chosen) row $`\chi _0`$ being replaced with $`(\mathrm{\Gamma },\mathrm{},\mathrm{\Gamma })`$. The stationary probabilities are the solution of the linear equation

$$(\stackrel{~}{𝐖}𝐩^{\mathrm{st}})_\chi =\mathrm{\Gamma }\delta _{\chi ,\chi _0}.$$
(6)
This can be solved for $`𝐩^{\mathrm{st}}`$ by inverting $`\stackrel{~}{𝐖}`$.

For a well-defined perturbation expansion in powers $`k`$ of the coupling strength $`\mathrm{\Gamma }`$ we write $`𝐖=_{k=1}^{\mathrm{}}𝐖^{(k)}`$, $`\stackrel{~}{𝐖}=_{k=1}^{\mathrm{}}\stackrel{~}{𝐖}^{(k)}`$, and $`𝐩^{\mathrm{st}}=_{k=0}^{\mathrm{}}𝐩^{\mathrm{st}(k)}`$. The order in $`\mathrm{\Gamma }`$ corresponds to the number of tunnel lines contained in the irreducible blocks. As a consequence, $`𝐖`$ and $`\stackrel{~}{𝐖}`$ start with first order in $`\mathrm{\Gamma }`$, whereas $`𝐩^{\mathrm{st}}`$ starts in zeroth order. The zeroth-order stationary probabilities are

$$p_\chi ^{\mathrm{st}(0)}=(\stackrel{~}{W}^{(1)})_{\chi \chi _0}^1\mathrm{\Gamma }$$
(7)
and higher-order corrections are obtained iteratively by

$$𝐩^{\mathrm{st}(k)}=\left(\stackrel{~}{𝐖}^{(1)}\right)^1\underset{m=0}{\overset{k1}{}}\stackrel{~}{𝐖}^{(km+1)}𝐩^{\mathrm{st}(m)}$$
(8)
for $`k=1,2,\mathrm{}`$. The diagrammatic rules for the calculation of the blocks $`𝐖^{(k)}`$ will be summarized at the end of this section.

The formulation of the perturbation expansion of the stationary probabilities is the first important step for a perturbation expansion of the current. The diagrammatic representation of contributions to the current involve a block $`𝐖^I`$, in which one (internal) vertex is replaced by an external vertex for the current operator $`\widehat{I}`$ divided by $`e/\mathrm{}`$. The current can also be expanded order by order in the coupling strength $`\mathrm{\Gamma }`$,

$$I^{(k)}=\frac{e}{2\mathrm{}}𝐞^T\underset{m=0}{\overset{k1}{}}𝐖^{I(km)}𝐩^{\mathrm{st}(m)}$$
(9)
for $`k=1,2,\mathrm{}`$, where the factor $`1/2`$ corrects for double counting of the external vertex being on the upper and lower branch of the Keldysh contour.

Now, we turn to the evaluation of the current noise $`S=2_{\mathrm{}}^0𝑑t[\widehat{I}(t)\widehat{I}(0)+\widehat{I}(0)\widehat{I}(t)2\widehat{I}^2]`$, which involves expectation values of two current operators at different times. These two current operators can either appear in the same irreducible block, denoted by $`𝐖^{II}`$, or in two different blocks. In total, the current noise is expressed as

$$S=\frac{e^2}{\mathrm{}}𝐞^T\left(𝐖^{II}𝐩^{\mathrm{st}}+𝐖^I\mathrm{𝐏𝐖}^I𝐩^{\mathrm{st}}\right),$$
(10)
(where the factor $`2`$ in the above definition has canceled against the correction factor $`1/2`$ for double counting) with

$$𝐏=_{\mathrm{}}^0𝑑t\frac{1}{\mathrm{}}\left[𝚷(0,t)𝚷(0,\mathrm{})\right].$$
(11)
The constant part $`𝚷(0,\mathrm{})`$ in the integrand is associated with the $`\widehat{I}^2`$ contribution to the noise. Furthermore, it guarantees that the integral is convergent and that $`𝐏`$ exists.

We also can formulate a perturbation expansion of $`𝐏`$ in orders of $`\mathrm{\Gamma }`$. To do so, we first relate $`𝐏`$ to the blocks $`𝐖`$ and the stationary probabilities $`𝐩^{\mathrm{st}}`$ by

$$\mathrm{𝐖𝐏}=𝐩^{\mathrm{st}}𝐞^T\mathrm{𝟏},$$
(12)
which follows from the Dyson equation and $`\mathrm{𝐖𝐩}^{\mathrm{st}}=\mathrm{𝟎}`$. As discussed above, $`𝐖`$ is not invertible, i.e., in addition to Eq. (12) we need some extra condition to determine $`𝐏`$. This extra condition is $`𝐞^T𝐏=\mathrm{𝟎}`$, which follows from the definition, Eq. (11), together with the Dyson equation, Eq. (5), and $`𝐞^T𝐖=\mathrm{𝟎}`$. Introducing the matrix $`𝐐`$ with elements given by

$$Q_{\chi ^{}\chi }=(p_\chi ^{}^{\mathrm{st}}\delta _{\chi ^{},\chi })(1\delta _{\chi ^{},\chi _0}),$$
(13)
Eq. (12) can be recast as $`\stackrel{~}{𝐖}𝐏=𝐐`$, which allows the determination of $`𝐏`$ by inversion of the matrix $`\stackrel{~}{𝐖}`$ introduced above. The solution for $`𝐏`$ can be expanded in orders of $`\mathrm{\Gamma }`$. We observe that, since the expansion of $`𝐐`$ starts in zeroth order in $`\mathrm{\Gamma }`$, $`𝐏`$ starts in order $`\mathrm{\Gamma }^1`$,

$$𝐏^{(1)}=(\stackrel{~}{𝐖}^{(1)})^1𝐐^{(0)}.$$
(14)
This reflects the fact that between the external current vertices the system evolves in time, as described by $`𝚷(0,t)`$. The time scale at which the dot relaxes and correlations decay is set by the inverse of the tunnel rates, i.e., by $`\mathrm{}/\mathrm{\Gamma }`$. As a consequence, the time integral in Eq. (11) is effectively cut at this time scale, i.e., the dominant contribution to $`𝐏`$ is of order $`\mathrm{\Gamma }^1`$. Furthermore, it is crucial that $`𝐏`$ starts in order $`\mathrm{\Gamma }^1`$, as otherwise the second term in Eq. (10) would not contribute to noise to first order in $`\mathrm{\Gamma }`$ (cf. Ref. korotkov, ), leading to a wrong result even for non-interacting systems ($`U=0`$), where exact expressions are available for comparison.hershfield

Corrections of higher order, $`k=0,1,\mathrm{}`$, can be computed from

$$𝐏^{(k)}=\left(\stackrel{~}{𝐖}^{(1)}\right)^1\left[𝐐^{(k+1)}\underset{m=1}{\overset{k1}{}}\stackrel{~}{𝐖}^{(km+1)}𝐏^{(m)}\right].$$
(15)
Eventually, we have all ingredients at hand to expand the noise in orders of $`\mathrm{\Gamma }`$, and we find

$$S^{(k)}=\frac{e^2}{\mathrm{}}𝐞^T\underset{m=0}{\overset{k1}{}}\left(𝐖^{II(km)}𝐩^{\mathrm{st}(m)}+\underset{m^{}=1}{\overset{km}{}}\underset{m^{\prime \prime }=1}{\overset{kmm^{}1}{}}𝐖^{I(kmm^{}m^{\prime \prime })}𝐏^{(m^{\prime \prime })}𝐖^{I(m^{})}𝐩^{\mathrm{st}(m)}\right)$$
(16)
for $`k=1,2,\mathrm{}`$. In lowest order this result simplifies to $`S^{(1)}=(e^2/\mathrm{})𝐞^T\left(𝐖^{II(1)}+𝐖^{I(1)}𝐏^{(1)}𝐖^{I(1)}\right)𝐩^{\mathrm{st}(0)}`$.

In summary of the technical part, we have derived a scheme to perform a systematic perturbation expansion of the current and the current noise. First, one calculates the irreducible blocks $`𝐖`$, $`𝐖^I`$, and $`𝐖^{II}`$. Second, the stationary probabilities $`𝐩^{\mathrm{st}}`$ are obtained from Eqs. (7) and (8) and $`𝐏`$ follows from Eqs. (14) and (15). Eventually, current and current noise are given by Eqs. (9) and (16).

The irreducible blocks $`𝐖`$ are calculated according to the rules specified in Ref. diagrams,  \[where the irreducible blocks were denoted by $`𝚺`$ and included a factor $`i`$ in the definition, $`𝚺=i𝐖`$\]. We complete this section with a summary of these diagrammatic rules, given as follows:

(1) For a given order $`k`$ draw all topologically different diagrams with $`2k`$ vertices connected by $`k`$ tunneling lines (an example of a first- and second-order diagram is shown in Fig. 2). Assign the energies $`ϵ_\chi `$ to the propagators and energies $`\omega _l`$ ($`l=1,\mathrm{},k`$) to the tunneling lines.

(2) For each of the ($`2k1`$) segments enclosed by two adjacent vertices there is a resolvent $`1/(\mathrm{\Delta }E_j+i0^+)`$, with $`j=1,\mathrm{},2k1`$, where $`\mathrm{\Delta }E_j`$ is the difference of the left-going minus the right-going energies.

(3) The contribution of a tunneling line of reservoir $`r`$ is $`(1/2\pi )\mathrm{\Gamma }_rf(\omega _l\mu _r)`$ if the line is going backward with respect to the closed time path and $`(1/2\pi )\mathrm{\Gamma }_r[1f(\omega _l\mu _r)]`$ if it is going forward. Here, $`f(x)`$ is the Fermi function.

(4) There is an overall prefactor $`(i)(1)^c`$, where $`c`$ is the total number of vertices on the backward propagator plus the number of crossings of tunneling lines plus the number of vertices connecting the state $`d`$ with $``$.

(5) Integrate over the energies $`\omega _l`$ of the tunneling lines and sum over all reservoir and spin indices.

The blocks $`𝐖^I`$ and $`𝐖^{II}`$ are determined in a similar way. The only difference to $`𝐖`$ is that in $`𝐖^I`$ ($`𝐖^{II}`$) one (two) internal vertices are replaced by external ones representing $`\widehat{I}\mathrm{}/e`$. This amounts to multiplying an overall prefactor, which arises due to the definition of the current operator and since the number of internal vertices on the backward propagator may have changed. We get a factor $`+1/2`$ for each external vertex on the upper (lower) branch of the Keldysh contour which describes tunneling of an electron into the right (left) or out of the left (right) lead, and $`1/2`$ in the other four cases. Finally, we have to sum up all the factors for each possibility to replace one (two) internal vertices by external ones.

## IV Results

In the following we discuss shot noise in the lowest- (first-) order perturbation theory in $`\mathrm{\Gamma }`$. We consider transport through a single level and allow for a finite spin splitting. Such a spin splitting could be realized due to the Zeeman effect of an external magnetic field (for semiconductor quantum dots) or an intrinsic exchange field due to atoms with magnetic moments (for molecules). The molecule can acquire four possible states: the level being unoccupied, occupied with either spin $``$ or $``$, or doubly occupied. The molecular level is characterized by level energies $`ϵ_{}`$ and $`ϵ_{}`$, the Coulomb repulsion $`U`$, and the coupling strengths $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$ to the electrodes. In addition, the electron distributions in the metal electrodes are governed by electrochemical potentials $`\mu _\mathrm{L}`$ and $`\mu _\mathrm{R}`$ and temperature $`T`$. We choose a symmetric bias voltage, such that $`\mu _\mathrm{L}=eV/2`$ and $`\mu _\mathrm{R}=eV/2`$. In this simple model without real spatial extent, the voltage is dropped entirely at the electrode-molecule tunnel junctions, meaning that the energies of the molecular states are independent of the applied voltage even if the couplings are not symmetric. For asymmetric coupling this might not be entirely realistic. However, we do not wish to include a heuristic parameter to describe the possibly unequal drop of the bias. In a comparison to an actual experiment this could be easily amended. We also do not include an explicit gate voltage, as its effects are straightforward to anticipate.

In Fig. 3 we plot the current and the current noise for a special choice of system parameters as a function of transport voltage. For this figure we consider a set of energy parameters ($`ϵ_{}=0.1\mathrm{eV}`$,$`ϵ_{}=0.2\mathrm{eV}`$, $`U=0.4\mathrm{eV}`$) that is comparable to energies encountered in small molecules. A rescaling of energies by a factor 1/100 would be necessary to obtain energies of the order achievable in semiconductor quantum dots and by external magnetic fields. If the same scaling is applied to the couplings and the temperature, the figure would remain unchanged. This is a special feature of first-order transport.

Electron transport becomes possible when charge excitations on the molecule become energetically allowed. Generally, at low bias, transport is exponentially suppressed, unless a degeneracy of states with different net charge is present (this could be tailored by application of a gate voltage, which we do not consider here). Each time when one of the four excitation energies $`ϵ_{}`$, $`ϵ_{}`$, $`ϵ_{}+U`$, or $`ϵ_{}+U`$ enters the energy window defined by the electrochemical potentials of the electrodes, a transport channel opens. This gives rise to plateaus, separated by thermally broadened steps. Due to the symmetric application of the bias, the steps occur at voltages of twice the corresponding excitation energy. The plateau heights depend on the coupling parameters $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$ only, i.e., they are independent of $`U`$ and $`T`$. The analytic expressions for current, noise, and Fano factor of these plateaus, which are labeled by $`i=0,\mathrm{},4`$, are given in Table 1.

The curves in Fig. 3 are normalized to $`I_{\mathrm{max}}=(e/\mathrm{})\mathrm{\Gamma }/2`$ and $`S_{\mathrm{max}}=(e^2/\mathrm{})\mathrm{\Gamma }/2`$, respectively, which is reached in the large-bias limit for symmetric coupling $`\mathrm{\Gamma }_\mathrm{L}=\mathrm{\Gamma }_\mathrm{R}`$. For asymmetric coupling, the plateaus are reduced in height. In Fig. 3, we show the results for $`3\mathrm{\Gamma }_\mathrm{L}=\mathrm{\Gamma }_\mathrm{R}`$ (dashed lines) and $`\mathrm{\Gamma }_\mathrm{L}=3\mathrm{\Gamma }_\mathrm{R}`$ (dotted lines) together with the case of symmetric coupling. The symmetry of our setup implies that all plateau heights are invariant under simultaneous exchange of $`\mathrm{\Gamma }_\mathrm{L}`$ with $`\mathrm{\Gamma }_\mathrm{R}`$ and $`\mu _\mathrm{L}`$ with $`\mu _\mathrm{R}`$. However, the plateau height can change if only $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$ are exchanged, or if only the bias voltage is reversed, as shown in our example for the two plateaus labeled by $`2`$ and $`3`$. This opens the possibility to access the asymmetry $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}`$ experimentally by reversing the bias voltage and comparing the plateau heights.

Figure 3 shows the result for one special choice of energy parameters and corresponding excitation energies. Nevertheless, Table 1 is complete in the sense that it contains all possible plateau values for any configuration of the excitation energies relative to the electrochemical potentials of the electrodes. The classification of the configurations and the algorithm to find the corresponding analytic expressions in Table 1 is given in Table 2. Without loss of generality we restrict ourselves to $`U0`$ and $`ϵ_{}ϵ_{}`$. The different configurations are classified by specifying which excitation energies lie within the energy window defined by the chemical potentials $`\mu _\mathrm{L},\mu _\mathrm{R}`$. We find 13 different possibilities, as listed in Table 2. For each case, the index $`i`$ indicates the column where the corresponding analytic expressions for current, noise and Fano factor can be found in Table 1. The indices $`2^{}`$ and $`3^{}`$ refer to columns $`2`$ and $`3`$ but with $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$ being exchanged.

In order to illustrate the use of the table we sketch the situation $`2a`$ in Fig. 4, realized in Fig. 3 in the region between $`0.4\mathrm{V}`$ and $`1\mathrm{V}`$. In this situation transport through both spin states is present as the molecule gets charged/uncharged in the sequential tunneling events. Double occupancy is still out of reach, since the excitation energies $`ϵ_{}+U`$ and $`ϵ_{}+U`$ are outside the energy window opened by the applied voltage.

We now turn to the discussion of the Fano factor $`F=S/2eI`$. In Fig. 5 we show the Fano factor as a function of bias voltage for the same parameters as in Fig. 3. Again we show three curves with $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}=1/1`$ (solid line), $`1/3`$ (dashed line) and $`3/1`$ (dotted line). At small bias, $`eVk_BT`$, the noise is dominated by thermal noise, described by the well-known hyperbolic cotangent behavior which leads to a divergence of the Fano factor.blanter ; loss  The plateau for bias voltages below $`0.2\mathrm{V}`$ corresponds to the Coulomb-blockade regime, where transport is exponentially suppressed (case $`0`$ in Table 2). In the region between $`0.2\mathrm{V}`$ and $`0.4\mathrm{V}`$ (case $`1a`$) transport through only one spin state ($``$) is possible. The Fano factor $`F_1`$ for this case has been derived earlier in Ref. nazarov, . For very large bias (region $`F_4`$), all states of the molecular level are involved in transport, and the Fano factor $`F_4`$ is again identical to the well-known blanter  formula for transport through a resonant level in the absence of Coulomb charging energy.

For the regions in between, corresponding to the cases $`2a`$ and $`3a`$, the Fano factor is different. The expression for $`F_2`$ has been very recently derived in Ref. bagrets, , while $`F_3`$, corresponding to region $`3a`$, is, to the best of knowledge, a new result.

On similar grounds as for the current and the noise, we find that the expressions for the Fano factors $`F_2`$ and $`F_3`$ are not invariant under exchange of $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$ alone, and also not invariant under reverse of the bias voltage alone. This is clearly seen in Fig. 5 in the very different plateau heights of the dotted and dashed curves, corresponding to exchange of $`\mathrm{\Gamma }_\mathrm{L}`$ and $`\mathrm{\Gamma }_\mathrm{R}`$. Furthermore, we see that all plateau heights lie between $`1/2`$ and $`1`$, and that the Fano factor, in general, is a nonmonotonic function of the bias voltage. We find that $`F_1=F_4`$ and $`F_3F_2`$ always hold, whereas $`F_1F_2`$ for $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}1/\sqrt{2}`$ and $`F_1F_3`$ for $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}1/2`$ only.

As a consequence, the pattern of the plateau sequence, in particular the relative height of the plateaus $`F_2`$ and $`F_3`$ compared to $`F_1=F_4`$, indicates not only the presence of an interaction or charging energy, but also gives overcomplete information on the ratio $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}`$ of the coupling strengths. This could be used in experiments to determine these parameters in a consistent way. In an experiment observing more than one plateau, the overcompleteness would give narrow constraints (due to experimental uncertainty) on whether a single interacting level can explain the measured values. Of course, it is always possible to fit $`n`$ plateaus with $`n`$ noninteracting levels and different couplings per level, so an absolute decision on the presence of interactions is not possible without the additional application of a gate voltage. However, if a fit with an interacting level is feasible, the principle of parsimony should favor the model with fewer parameters.

In the absence of a spin splitting, $`ϵ_{}=ϵ_{}`$, the number of plateaus is reduced. An example for the Fano factor as a function of bias voltage is displayed in Fig. 6. The solid curve is for symmetric, and the dotted and dashed curves are for asymmetric coupling. Again, a nonmonotonic sequence of plateau heights, $`F_4>F_2`$, indicates the presence of Coulomb charging and $`\mathrm{\Gamma }_\mathrm{L}/\mathrm{\Gamma }_\mathrm{R}<1/\sqrt{2}`$. For comparison we remark that for both negligible spin splitting and zero charging energy ($`U=0`$), only the plateaus $`F_0`$ and $`F_4`$ appear, the plateau heights are invariant under reverse of bias, and no nonmonotonic behavior is seen, as also shown in Fig. 6.

Finally, we briefly comment on the steps connecting the plateaus of the Fano factor. In the present first-order approximation the steps are solely broadened by temperature. In addition, we observe that sometimes a peak in the Fano factor shows up. This happens, for example, at the step between $`F_1`$ and $`F_2`$ for $`\mathrm{\Gamma }_\mathrm{L}=\mathrm{\Gamma }_\mathrm{R}`$, as seen in Fig. 5. The peak height, $`0.5625`$, exceeds that of the adjacent plateaus, $`0.5`$ and $`5/9`$. These features can even appear in the regime of negligible Coulomb charging energy, as previously shown in Ref. haug, . However, we point out that the behavior at the steps is going to be strongly affected by second-order tunneling events, so-called cotunneling. For example, the current steps will show additional broadening due to the intrinsic linewidth $`\mathrm{\Gamma }`$. On the other hand, the plateau values in regimes 1 - 4 will not be affected by second-order effects.

In summary, we presented a theory of transport through a molecule or quantum dot that allows for a systematic perturbation expansion of current and shot noise in the coupling strength between molecular orbital (dot level) and electrodes. We analyzed first-order transport in detail. In a complete overview, we derived analytic expressions for current, noise, and the Fano factor for all possible situations involving a single interacting level. In particular, we discussed how the sequence of Fano factor plateaus can provide information about the asymmetry ratio of the coupling strengths as well as the presence of Coulomb charging energy.

Acknowledgments. We enjoyed interesting and helpful discussions with A. Braggio, G. Johansson, H. Schoeller, H. Weber, and A. Zaikin, as well as financial support by the DFG via the Center for Functional Nanostructures and the Emmy-Noether program.


text:
# PHYSICS Version of Random Noise and Pole-Dynamics in Unstable Front Propagation(new version)

## I Introduction

In our last paper about flame front propagation we left two open problems. The first one is how to explain the existence of small dependence parameters of problem on the noise. The second problem is how to calculate numerically such values as excess number of poles in system , number of poles that appear in the system in unit of time, life time of pole.To solve these problem we write this paper.

## II Equations of Motion and Pole-decomposition in the Channel Geometry

It is known that planar flames freely propagating through initially motionless homogeneous combustible mixtures are intrinsically unstable. It was reported that such flames develop characteristic structures which include cusps, and that under usual experimental conditions the flame front accelerates as time goes on. A model in $`1+1`$ dimensions that pertains to the propagation of flame fronts in channels of width $`\stackrel{~}{L}`$ was proposed in . It is written in terms of position $`h(x,t)`$ of the flame front above the $`x`$-axis. After appropriate rescalings it takes the form:

$$\frac{h(x,t)}{t}=\frac{1}{2}\left[\frac{h(x,t)}{x}\right]^2+\nu \frac{^2h(x,t)}{x^2}+I\{h(x,t)\}+1.$$
(1)
The domain is $`0<x<\stackrel{~}{L}`$, $`\nu `$ is a parameter and we use periodic boundary conditions. The functional $`I[h(x,t)]`$ is the Hilbert transform of derivative which is conveniently defined in terms of the spatial Fourier transform

$`h(x,t)={\displaystyle _{\mathrm{}}^{\mathrm{}}}e^{ikx}\widehat{h}(k,t)𝑑k`$ (2)
$`I[h(k,t)]=|k|\widehat{h}(k,t)`$ (3)
For the purpose of introducing the pole-decomposition it is convenient to rescale the domain to $`0<\theta <2\pi `$. Performing this rescaling and denoting the resulting quantities with the same notation we have

$`{\displaystyle \frac{h(\theta ,t)}{t}}={\displaystyle \frac{1}{2L^2}}\left[{\displaystyle \frac{h(\theta ,t)}{\theta }}\right]^2+{\displaystyle \frac{\nu }{L^2}}{\displaystyle \frac{^2h(\theta ,t)}{\theta ^2}}`$ (4)
$`+{\displaystyle \frac{1}{L}}I\{h(\theta ,t)\}+1.`$ (5)
In this equation $`L=\stackrel{~}{L}/2\pi `$. Next we change variables to $`u(\theta ,t)h(\theta ,t)/\theta `$. We find

$$\frac{u(\theta ,t)}{t}=\frac{u(\theta ,t)}{L^2}\frac{u(\theta ,t)}{\theta }+\frac{\nu }{L^2}\frac{^2u(\theta ,t)}{\theta ^2}+\frac{1}{L}I\{u(\theta ,t)\}.$$
(6)
It is well known that the flat front solution of this equation is linearly unstable. The linear spectrum in $`k`$-representation is

$$\omega _k=|k|/L\nu k^2/L^2.$$
(7)
There exists a typical scale $`k_{max}`$ which is the last unstable mode

$$k_{max}=\frac{L}{\nu }.$$
(8)
Nonlinear effects stabilize a new steady-state which is discussed next.

The outstanding feature of the solutions of this equation is the appearance of cusp-like structures in the developing fronts. Therefore a representation in terms of Fourier modes is very inefficient. Rather, it appears very worthwhile to represent such solutions in terms of sums of functions of poles in the complex plane. It will be shown below that the position of the cusp along the front is determined by the real coordinate of the pole, whereas the height of the cusp is in correspondence with the imaginary coordinate. Moreover, it will be seen that the dynamics of the developing front can be usefully described in terms of the dynamics of the poles. Following  we expand the solutions $`u(\theta ,t)`$ in functions that depend on $`N`$ poles whose position $`z_j(t)x_j(t)+iy_j(t)`$ in the complex plane is time dependent:

$`u(\theta ,t)=\nu {\displaystyle \underset{j=1}{\overset{N}{}}}\mathrm{cot}\left[{\displaystyle \frac{\theta z_j(t)}{2}}\right]+c.c.`$ (9)
$`=\nu {\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{2\mathrm{sin}[\theta x_j(t)]}{\mathrm{cosh}[y_j(t)]\mathrm{cos}[\theta x_j(t)]}},`$ (10)
$$h(\theta ,t)=2\nu \underset{j=1}{\overset{N}{}}\mathrm{ln}\left[\mathrm{cosh}(y_j(t))\mathrm{cos}(\theta x_j(t))\right]+C(t).$$
(11)
In (11) $`C(t)`$ is a function of time. The function (11) is a superposition of quasi-cusps (i.e. cusps that are rounded at the tip). The real part of the pole position (i.e. $`x_j`$) is the coordinate (in the domain $`[0,2\pi ]`$) of the maximum of the quasi-cusp, and the imaginary part of the pole position (i.e $`y_j`$) is related to the depth of the quasi-cusp. As $`y_j`$ decreases the depth of the cusp increases. As $`y_j0`$ the depth diverges to infinity. Conversely, when $`y_j\mathrm{}`$ the depth decreases to zero.

The main advantage of this representation is that the propagation and wrinkling of the front can be described via the dynamics of the poles. Substituting (10) in (6) we derive the following ordinary differential equations for the positions of the poles:

$$L^2\frac{dz_j}{dt}=\left[\nu \underset{k=1,kj}{\overset{2N}{}}\mathrm{cot}\left(\frac{z_jz_k}{2}\right)+i\frac{L}{2}sign[Im(z_j)]\right].$$
(12)
We note that in (10), due to the complex conjugation, we have $`2N`$ poles which are arranged in pairs such that for $`j<N`$ $`z_{j+N}=\overline{z}_j`$. In the second sum in (10) each pair of poles contributed one term. In Eq.(12) we again employ $`2N`$ poles since all of them interact. We can write the pole dynamics in terms of the real and imaginary parts $`x_j`$ and $`y_j`$. Because of the arrangement in pairs it is sufficient to write the equation for either $`y_j>0`$ or for $`y_j<0`$. We opt for the first. The equations for the positions of the poles read

$`L^2{\displaystyle \frac{dx_j}{dt}}=\nu {\displaystyle \underset{k=1,kj}{\overset{N}{}}}\mathrm{sin}(x_jx_k)[[\mathrm{cosh}(y_jy_k)`$ (13)
$`\mathrm{cos}(x_jx_k)]^1+[\mathrm{cosh}(y_j+y_k)\mathrm{cos}(x_jx_k)]^1]`$ (14)
$`L^2{\displaystyle \frac{dy_j}{dt}}=\nu {\displaystyle \underset{k=1,kj}{\overset{N}{}}}({\displaystyle \frac{\mathrm{sinh}(y_jy_k)}{\mathrm{cosh}(y_jy_k)\mathrm{cos}(x_jx_k)}}`$ (15)
$`+{\displaystyle \frac{\mathrm{sinh}(y_j+y_k)}{\mathrm{cosh}(y_j+y_k)\mathrm{cos}(x_jx_k)}})+\nu \mathrm{coth}(y_j)L.`$ (16)
We note that if the initial conditions of the differential equation (6) are expandable in a finite number of poles, these equations of motion preserve this number as a function of time. On the other hand, this may be an unstable situation for the partial differential equation, and noise can change the number of poles. This issue will be examined at length in Section III. We turn now to a discussion of the steady state solution of the equations of the pole-dynamics.

### A Qualitative properties of the stationary solution

The steady-state solution of the flame front propagating in channels of width $`2\pi `$ was presented in Ref.. Using these results we can immediately translate the discussion to a channel of width $`L`$. The main results are summarized as follows:

1. There is only one stable stationary solution which is geometrically represented by a giant cusp (or equivalently one finger) and analytically by $`N(L)`$ poles which are aligned on one line parallel to the imaginary axis. The existence of this solution is made clearer with the following remarks.
2. There exists an attraction between the poles along the real line. This is obvious from Eq.(13) in which the sign of $`dx_j/dt`$ is always determined by $`\mathrm{sin}(x_jx_k)`$. The resulting dynamics merges all the $`x`$ positions of poles whose $`y`$-position remains finite.
3. The $`y`$ positions are distinct, and the poles are aligned above each others in positions $`y_{j1}<y_j<y_{j+1}`$ with the maximal being $`y_{N(L)}`$. This can be understood from Eq.(16) in which the interaction is seen to be repulsive at short ranges, but changes sign at longer ranges.
4. If one adds an additional pole to such a solution, this pole (or another) will be pushed to infinity along the imaginary axis. If the system has less than $`N(L)`$ poles it is unstable to the addition of poles, and any noise will drive the system towards this unique state. The number $`N(L)`$ is

$$N(L)=\left[\frac{1}{2}\left(\frac{L}{\nu }+1\right)\right],$$
(17)
where $`\left[\mathrm{}\right]`$ is the integer part. To see this consider a system with $`N`$ poles and such that all the values of $`y_j`$ satisfy the condition $`0<y_j<y_{max}`$. Add now one additional pole whose coordinates are $`z_a(x_a,y_a)`$ with $`y_ay_{max}`$. From the equation of motion for $`y_a`$, (16) we see that the terms in the sum are all of the order of unity as is also the $`\mathrm{cot}(y_a)`$ term. Thus the equation of motion of $`y_a`$ is approximately

$$\frac{dy_a}{dt}\nu \frac{2N+1}{L^2}\frac{1}{L}.$$
(18)
The fate of this pole depends on the number of other poles. If $`N`$ is too large the pole will run to infinity, whereas if $`N`$ is small the pole will be attracted towards the real axis. The condition for moving away to infinity is that $`N>N(L)`$ where $`N(L)`$ is given by (17). On the other hand the $`y`$ coordinate of the poles cannot hit zero. Zero is a repulsive line, and poles are pushed away from zero with infinite velocity. To see this consider a pole whose $`y_j`$ approaches zero. For any finite $`L`$ the term $`\mathrm{coth}(y_j)`$ grows unboundedly whereas all the other terms in Eq.(16) remain bounded.
5. The height of the cusp is proportional to $`L`$. The distribution of positions of the poles along the line of constant $`x`$ was worked out in .

We will refer to the solution with all these properties as the Thual-Frisch-Henon (TFH)-cusp solution.

## III Acceleration of the Flame Front, Pole Dynamics and Noise

A major motivation of this Section is the observation that in radial geometry the same equation of motion shows an acceleration of the flame front. The aim of this section is to argue that this phenomenon is caused by the noisy generation of new poles. Moreover, it is our contention that a great deal can be learned about the acceleration in radial geometry by considering the effect of noise in channel growth. In Ref.  it was shown that any initial condition which is represented in poles goes to a unique stationary state which is the giant cusp which propagates with a constant velocity $`v=1/2`$ up to small $`1/L`$ corrections. In light of our discussion of the last section we expect that any smooth enough initial condition will go to the same stationary state. Thus if there is no noise in the dynamics of a finite channel, no acceleration of the flame front is possible. What happens if we add noise to the system?

For concreteness we introduce an additive white-noise term $`\eta (\theta ,t)`$to the equation of motion (6) where

$$\eta (\theta ,t)=\underset{k}{}\eta _k(t)\mathrm{exp}(ik\theta ),$$
(19)
and the Fourier amplitudes $`\eta _k`$ are correlated according to

$$<\eta _k(t)\eta _k^{}^{}(t^{})>=\frac{f}{L}\delta _{k,k^{}}\delta (tt^{}).$$
(20)
We will first examine the result of numerical simulations of noise-driven dynamics, and later return to the theoretical analysis.

### A Noisy Simulations

Previous numerical investigations  did not introduce noise in a controlled fashion. We will argue later that some of the phenomena encountered in these simulations can be ascribed to the (uncontrolled) numerical noise. We performed numerical simulations of Eq.(6 using a pseudo-spectral method. The time-stepping scheme was chosen as Adams-Bashforth with 2nd order presicion in time. The additive white noise was generated in Fourier-space by choosing $`\eta _k`$ for every $`k`$ from a flat distribution in the interval $`[\sqrt{2\frac{f}{L}},\sqrt{2\frac{f}{L}}]`$. We examined the average steady state velocity of the front as a function of $`L`$ for fixed $`f`$ and as a function of $`f`$ for fixed $`L`$. We found the interesting phenomena that are summarized here:

1. In Fig.2 we can see two different regimes of behavior the average velocity $`v`$ as function of noise $`f^0.5`$ for fixed system size L. For the noise $`f`$ smaller then same fixed value $`f_{cr}`$

$$vf^\xi .$$
(21)
For these values of $`f`$ this dependence is very weak, and $`\xi 0.02`$. For large values of $`f`$ the dependence is much stronger
2. In Fig.1 we can see growth of the average velocity $`v`$ as function of the system size L. After some values of L we can see saturation of the velocity. For regime $`f<f_{cr}`$ the growth of the velocity can be written as

$$vL^\mu ,\mu 0.40\pm 0.05.$$
(22)

### B Calculation of Poles Number in the System

The interesting problem that we would like to solve here it is to find number of poles that exist in our system outside the giant cusp. We can make it by next way: to calculate number of cusps (points of minimum or inflexional points) and their position on the interval $`\theta :[0,2\pi ]`$ in every moment of time and to draw positions of cusp like function of time, see Fig.3.

We assume that our system is almost all time in ”quasi-stable” state, i.e. every new cusp that appears in the system includes only one pole. By help pictures obtained by such way we can find

1. By calculation number of cusp in some moment of time and by investigation of history of every cusp (except the giant cusp) , i.e. how many initial cusps take part in formation this cusp, after averaging with respect to different moments of time we can find mean number of poles that exist in our system outside the giant cusp. Let us denote this number $`\delta N`$. We can see four regimes that can be define with respect to dependence of this number on noise $`f`$:

(i) Regime I: Such small noise that no poles exist in our system outside the giant cusp.

(ii) Regime II Strong dependence of poles number $`\delta N`$ on noise $`f`$.

(iii) Regime III Saturation poles number $`\delta N`$ on noise $`f`$, so we see very small dependence of this number on noise

$$\delta Nf^{0.03}$$
(23)
Saturated value of $`\delta N`$ is defined by next formula

$$\delta NN(L)/2\frac{1}{4}\frac{L}{\nu }$$
(24)
where $`N(L)\frac{1}{2}\frac{L}{\nu }`$ is number of poles in giant cusp.

(iv) Regime IV We again see strong dependence of poles number $`\delta N`$ on noise $`f`$.

$$\delta Nf^{0.1}$$
(25)
Because of numerical noise we can see in most of simulations only regime III and IV. In future if we don’t note something different we discuss regime III.
2. By calculation of new cusp number we can find number of poles that appear in the system in unit of time $`\frac{dN}{dt}`$. In regime III

$$\frac{dN}{dt}f^{0.03}$$
(26)
Dependence on $`L`$ and $`\nu `$ define by

$$\frac{dN}{dt}L^{0.8}$$
(27)
$$\frac{dN}{dt}\frac{1}{\nu ^2}$$
(28)
And in regime IV

$$\frac{dN}{dt}f^{0.1}$$
(29)

### C Theoretical Discussion of the Effect of Noise

#### 1 The Threshold of Instability to Added Noise. Transition from regime I to regime II

First we present the theoretical arguments that explain the sensitivity of the giant cusp solution to the effect of added noise. This sensitivity increases dramatically with increasing the system size $`L`$. To see this we use again the relationship between the linear stability analysis and the pole dynamics.

Our additive noise introduces perturbations with all $`k`$-vectors. We showed previously that the most unstable mode is the $`k=1`$ component $`A_1\mathrm{sin}(\theta )`$. Thus the most effective noisy perturbation is $`\eta _1\mathrm{sin}(\theta )`$ which can potentially lead to a growth of the most unstable mode. Whether or not this mode will grow depends on the amplitude of the noise. To see this clearly we return to the pole description. For small values of the amplitude $`A_1`$ we represent $`A_1\mathrm{sin}(\theta )`$ as a single pole solution of the functional form $`\nu e^y\mathrm{sin}\theta `$. The $`y`$ position is determined from $`y=\mathrm{log}|A_1|/\nu `$, and the $`\theta `$-position is $`\theta =\pi `$ for positive $`A_1`$ and $`\theta =0`$ for negative $`A_1`$. From the analysis of Section III we know that for very small $`A_1`$ the fate of the pole is to be pushed to infinity, independently of its $`\theta `$ position; the dynamics is symmetric in $`A_1A_1`$ when $`y`$ is large enough. On the other hand when the value of $`A_1`$ increases the symmetry is broken and the $`\theta `$ position and the sign of $`A_1`$ become very important. If $`A_1>0`$ there is a threshold value of $`y`$ below which the pole is attracted down. On the other hand if $`A_1<0`$, and $`\theta =0`$ the repulsion from the poles of the giant cusp grows with decreasing $`y`$. We thus understand that qualitatively speaking the dynamics of $`A_1`$ is characterized by an asymmetric “potential” according to

$`\dot{A}_1`$ $`=`$ $`{\displaystyle \frac{V(A_1)}{A_1}},`$ (30)
$`V(A_1)`$ $`=`$ $`\lambda A_1^2aA_1^3+\mathrm{}.`$ (31)
¿From the linear stability analysis we know that $`\lambda \nu /L^2`$, cf. Eq.(18). We know further that the threshold for nonlinear instability is at $`A_1\nu ^3/L^2`$, cf. Eq(LABEL:nu3L2). This determines that value of the coefficient $`a2/3\nu ^2`$. The magnitude of the “potential” at the maximum is

$$V(A_{max})\nu ^7/L^6.$$
(32)
The effect of the noise on the development of the mode $`A_1\mathrm{sin}\theta `$ can be understood from the following stochastic equation

$$\dot{A}_1=\frac{V(A_1)}{A_1}+\eta _1(t).$$
(33)
It is well known  that for such dynamics the rate of escape $`R`$ over the “potential” barrier for small noise is proportional to

$$R\frac{\nu }{L^2}\mathrm{exp}^{\nu ^7/\frac{f}{L}L^6}.$$
(34)
The conclusion is that any arbitrarily tiny noise becomes effective when the system size increase and when $`\nu `$ decreases. If we drive the system with noise of amplitude $`\frac{f}{L}`$ the system can always be sensitive to this noise when its size exceeds a critical value $`L_c`$ that is determined by $`\frac{f}{L_c}\nu ^7/L_c^6`$. This formula defines transition from regime I (no poles) to regime II. For $`L>L_c`$ the noise will introduce new poles into the system. Even numerical noise in simulations involving large size systems may have a macroscopic influence.

The appearance of new poles must increase the velocity of the front. The velocity is proportional to the mean of $`(u/L)^2`$. New poles distort the giant cusp by additional smaller cusps on the wings of the giant cusp, increasing $`u^2`$. Upon increasing the noise amplitude more and more smaller cusps appear in the front, and inevitably the velocity increases. This phenomenon is discussed quantitatively in Section III.

#### 2 Verifing of asymmetric ”potential” form

If we know distribution of poles in the giant cusp we can find form of ”potential” and verify expressions for $`\lambda `$, $`A_{max}`$ and $`\frac{V(A_1)}{A_1}`$. For measurement of $`A_1`$ we use formula $`A_1=4\nu e^y`$ and $`\frac{V(A_1)}{A_1}=4\nu \frac{dy}{dt}e^y`$ where $`\frac{dy}{dt}`$ can be find from equation motion of poles. Numerical measurements were made for $`L=2n\nu ,ninteger,n>2`$. For graphs that we use for fixed $`\nu `$ and variable $`L`$;$`\nu =0.005`$ $`L:[1,150]`$. For fixed $`L`$ and variable $`\nu `$; $`L=1`$ and $`\nu :[0.005,0.05]`$. We can find $`A_{max}`$ as zero-point of $`\frac{V(A_1)}{A_1}`$. Obtain results are next. For $`\frac{A_{max}L^2}{\nu ^3}`$ as function of $`L`$ (Fig 10)

and $`\nu `$ (Fig 11)

we can see that $`\frac{A_{max}L^2}{\nu ^3}`$ is almost constant. Dependence $`\frac{A_{max}}{A_{N(L)}}`$($`A_{N(L)}`$ is defined by position of the upper pole.)as function $`L`$ (Fig 12)

and $`\nu `$ (Fig 13)

is almost $`Const0.5`$. From graphs Fig 14 (for $`L`$)

and Fig 15

(for $`\nu `$) demonstrate us that $`\frac{\lambda L^2}{\nu }`$ is almost constant.

Fig 16 give us full dependence of $`\frac{f}{L}`$ on $`L`$. These results are in good agreement with the theory.

#### 3 The Noisy Steady State and its Collapse with Large Noise and System Size

In this subsection we discuss the response of the giant cusp solution to noise levels that are able to introduce a large number of excess poles in addition to those existing in the giant cusp. We will denote the excess number of poles by $`\delta N`$. The first question that we address is how difficult is it to insert yet an additional pole when there is already a given excess $`\delta N`$. To this aim we estimate the effective potential $`V_{\delta N}(A_1)`$ which is similar to (31) but is taking into account the existence of an excess number of poles. A basic approximation that we employ is that the fundamental form of the giant cusp solution is not seriously modified by the existence of an excess number of poles. Of course this approximation breaks down quantitatively already with one excess pole. Qualitatively however it holds well until the excess number of poles is of the order of the original number $`N(L)`$ of the giant cusp solution. Another approximation is that the rest of the linear modes play no role in this case. At this point we limit the discussion therefore to the situation $`\delta NN(L)`$ (regime II).

To estimate the parameter $`\lambda `$ in the effective potential we consider the dynamics of one pole whose $`y`$ position $`y_a`$ is far above $`y_{max}`$. According to Eq.(18) the dynamics reads

$$\frac{dy_a}{dt}\frac{2\nu (N(L)+\delta N)}{L^2}\frac{1}{L}$$
(35)
Since the $`N(L)`$ term cancels against the $`L^1`$ term (cf. Sec. II A), we remain with a repulsive term that in the effective potential translates to

$$\lambda =\frac{\nu \delta N}{L^2}.$$
(36)
Next we estimate the value of the potential at the break-even point between attraction and repulsion. In the last subsection we saw that a foreign pole has to be inserted below $`y_{max}`$ in order to be attracted towards the real axis. Now we need to push the new pole below the position of the existing pole whose index is $`N(L)\delta N`$. This position is estimated as in Sec III C by employing the TFH distribution function (LABEL:dist). We find

$$y_{\delta N}2\mathrm{ln}\left[\frac{4L}{\pi ^2\nu \delta N}\right].$$
(37)
As before, this implies a threshold value of the amplitude of single pole solution $`A_{max}\mathrm{sin}\theta `$ which is obtained from equating $`A_{max}=\nu e^{y_{\delta N}}`$. We thus find in the present case $`A_{max}\nu ^3(\delta N)^2/L^2`$. Using again a cubic representation for the effective potential we find $`a=2/(3\nu ^2\delta N)`$ and

$$V(A_{max})=\frac{1}{3}\frac{\nu ^7(\delta N)^5}{L^6}.$$
(38)
Repeating the calculation of the escape rate over the potential barrier we find in the present case

$$R\frac{\nu \delta N}{L^2}\mathrm{exp}^{\nu ^7(\delta N)^5/\frac{f}{L}L^6}.$$
(39)
For a given noise amplitude $`\frac{f}{L}`$ there is always a value of $`L`$ and $`\nu `$ for which the escape rate is of $`O(1)`$ as long as $`\delta N`$ is not too large. When $`\delta N`$ increases the escape rate decreases, and eventually no additional poles can creep into the system. The typical number $`\delta N`$ for fixed values of the parameters is estimated from equating the argument in the exponent to unity:

$$\delta N\left(\frac{f}{L}L^6/\nu ^7\right)^{1/5}.$$
(40)
We can see that $`\delta N`$ depend on noise $`f`$ very seriously. It is not the case in regime III. Let us find conditions of transition from regime II to regime III, where we see saturation of $`\delta N`$ with respect to noise $`f`$.

(i) We use for the amplitude of pole solution that really equal to $`\frac{2\nu \mathrm{sin}\theta }{\mathrm{cosh}(y_{\delta N})\mathrm{cos}\theta }`$ expression $`A_{max}=4\nu e^{y_{\delta N}}`$ but it is right only for big number $`y_{\delta N}`$. For $`y_{\delta N}<1`$ better approximation is $`A_{max}=\frac{4\nu }{y_{\delta N}^2}`$. From equation (37) we can find that the boundary value $`y_{\delta N}=1`$ correspond to $`\delta NN(L)/2`$

(ii) We use expression $`y_{\delta N}2\mathrm{ln}\left[\frac{4L}{\pi ^2\nu \delta N}\right]`$ but for big value of $`\delta N`$ better approximation that can be find the same way is $`y_{\delta N}\frac{\pi ^2\nu }{2L}(N(L)\delta N)\mathrm{ln}\left[\frac{8eL}{\pi ^2\nu (N(L)\delta N)}\right]`$. These expressions give us nearly equal result for $`\delta NN(L)/2`$.

From (i) and (ii) we can make next conclusions

(a) Transition from regime II to regime III happens for nearly $`\delta NN(L)/2`$

(b) Using new expression in (i) and (ii) for amplitude $`A_{max}`$ and $`y_{\delta N}`$ we can find for noise $`\frac{f}{L}`$ in regime III:

$$\frac{f}{L}V(A_{max})\lambda A_{max}^2\frac{\nu \delta N}{L^2}(\frac{4\nu }{y_{\delta N}^2})^2\frac{L^2}{\nu }\frac{\delta N}{(N(L)\delta N)^4}$$
(41)
This expression define very slow dependence of $`\delta N`$ on noise $`\frac{f}{L}`$ for $`\delta N>N(L)/2`$ that explain noise saturation of $`\delta N`$ for regime III.

(c) Form of the giant cusp solution is defined by poles that are closely to zero with respect to $`y`$. For regime III $`N(L)/2`$ poles that have position $`y<y_{\delta N=N(L)/2}=1`$ stay on these place. This result explain why giant cusp solution is not seriously modified for regime III.

From Eq.(40) by help of boundary condition

$$\delta NN(L)/2$$
(42)
boundary noise $`f_b`$ between regime II and III can be found

$$f_b\nu ^2$$
(43)
The basic equation describing pole dynamics is next

$$\frac{dN}{dt}=\frac{\delta N}{T}$$
(44)
where $`\frac{dN}{dt}`$ is number of poles that appear in unit of time in our system, $`\delta N`$ is excess number of poles, T is mean life time of pole (between appearing and merging with giant cusp). Using result of numerical simulations for $`\frac{dN}{dt}`$ and (42) we can find for $`T`$

$$T=\frac{\delta N}{\frac{dN}{dt}}\nu L^{0.2}$$
(45)
So life time proportional to $`\nu `$ and depend on system size $`L`$ very weakly.

### D The acceleration of the flame front due to noise

In this section we estimate the scaling exponents that characterize the velocity of the flame front as a function of the system size. Our arguments in this section are even less solid than the previous ones, but nevertheless we believe that we succeed to capture some of the essential qualitative physics that underlies the interaction between noise and instability and which results in the acceleration of the flame front.

To estimate the velocity of the flame front we need to write down an equation for the mean of $`<dh/dt>`$ given an arbitrary number $`N`$ of poles in the system. This equation follows directly from (5):

$$\frac{dh}{dt}=\frac{1}{L^2}\frac{1}{2\pi }_0^{2\pi }u^2𝑑\theta .$$
(46)
After substitution of (10) in (46) we get, using (13) and (16)

$$\frac{dh}{dt}=2\nu \underset{k=1}{\overset{N}{}}\frac{dy_k}{dt}+2\left(\frac{\nu N}{L}\frac{\nu ^2N^2}{L^2}\right).$$
(47)
The estimates of the second and third terms in this equation are straightforward. Writing $`N=N(L)+\delta N(L)`$ and remembering that $`N(L)L/\nu `$ and $`\delta N(L)N(L)/2`$ we find that these terms contribute $`O(1)`$. The first term will contribute only when the current of poles is asymmetric. Since noise introduces poles at a finite value of $`y_{min}`$, whereas the rejected poles stream towards infinity and disappear at boundary of nonlinearity defined by position of highest pole

$$y_{max}2\mathrm{ln}\left[\frac{4L}{\pi ^2\nu }\right].$$
(48)
, we have an asymmetry that contributes to the velocity of the front. To estimate the first term let us define

$$d(\frac{dy_k}{dt})=\underset{l}{\overset{l+dl}{}}\frac{dy_k}{dt}$$
(49)
where $`_l^{l+dl}\frac{dy_k}{dt}`$ is sum over poles that are on the interval $`y:[l,l+dl]`$. We can write

$$d(\frac{dy_k}{dt})=d(\frac{dy_k}{dt})_{up}d(\frac{dy_k}{dt})_{down}$$
(50)
Where $`d(\frac{dy_k}{dt})_{up}`$ flux of poles moving up and $`d(\frac{dy_k}{dt})_{down}`$ flux of poles moving down.

For these flux we can write

$$d(\frac{dy_k}{dt})_{up},d(\frac{dy_k}{dt})_{down}\frac{dN}{dt}dl$$
(51)
So for the first term

$`0{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{dy_k}{dt}}=`$ (52)
$`{\displaystyle _{y_{min}}^{y_{max}}}{\displaystyle \frac{d(\frac{dy_k}{dt})}{dl}}𝑑l`$ (53)
$`={\displaystyle _{y_{min}}^{y_{max}}}{\displaystyle \frac{d(\frac{dy_k}{dt})_{up}d(\frac{dy_k}{dt})_{down}}{dl}}𝑑l`$ (54)
$`{\displaystyle \frac{dN}{dt}}(y_{max}y_{min})`$ (55)
$`{\displaystyle \frac{dN}{dt}}y_{max}`$ (56)
Because of slow($`\mathrm{ln}`$) dependence of $`y_{max}`$ on $`L`$ and $`\nu `$ $`\frac{dN}{dt}`$ term define oder of nonlinearity for first term. This term zero for symmetric current of poles and achieves maximum for maximal asymmetric current of poles. Comparison $`vL^{0.42}f^{0.02}`$ and $`\frac{dN}{dt}L^{0.8}f^{0.03}`$ confirm this calculation.

## IV summary and conclusions

The main two messages of this paper are: (i) There is an important interaction between the instability of developing fronts and random noise; (ii) This interaction and its implications can be understood qualitatively and sometimes quantitatively using the description in terms of complex poles.

The pole description is natural in this context firstly because it provides an exact (and effective) representation of the steady state without noise. Once one succeeds to describe also the perturbations about this steady state in terms of poles, one achieves a particularly transparent language for the study of the interplay between noise and instability. This language also allows us to describe in qualitative and semi-quantitative terms the inverse cascade process of increasing typical lengths when the system relaxes to the steady state from small, random initial conditions.

The main conceptual steps in this paper are as follows: firstly one realizes that the steady state solution, which is characterized by $`N(L)`$ poles aligned along the imaginary axis is marginally stable against noise in a periodic array of $`L`$ values. For all values of $`L`$ the steady state is nonlinearly unstable against noise. The main and foremost effect of noise of a given amplitude $`f`$ is to introduce an excess number of poles $`\delta N(L,f)`$ into the system. The existence of this excess number of poles is responsible for the additional wrinkling of the flame front on top of the giant cusp, and for the observed acceleration of the flame front. By considering the noisy appearance of new poles we rationalize the observed scaling laws as a function of the noise amplitude and the system size.

Theoretically we therefore concentrate on estimating $`\delta N(L,f)`$. The measurements do not test our theoretical consideration directly, but rather test the dependence of the velocity on $`L`$ and $`f`$. The only direct test for our theory is the critical line shown in Fig.7. The measured exponent is in accord with our analytic estimates. Nevertheless we note that some of our consideration are only qualitative. For example, we estimated $`\delta N(L,f)`$ by assuming that the giant cusp solution is not seriously perturbed. On the other hand we find a flux of poles going to infinity due to the introduction of poles at finite values of $`y`$ by the noise. The existence of poles spread between $`y_{max}`$ and infinity is a significant perturbation of the giant cusp solution. Thus also the comparison between the various scaling exponents measured and predicted must be done with caution; we cannot guarantee that those cases in which our prediction hit close to the measurement mean that the theory is quantitative. However we believe that our consideration extract the essential ingredients of a correct theory.

The “phase diagram” as a function of $`L`$ and $`f`$ in this system consists of three regimes. In the first one, discussed in Section III C 1 , the noise is too small to have any effect on the giant cusp solution. In the second the noise introduces excess poles that serve to decorate the giant cusp with side cusps. In this regime we find scaling laws for the velocity as a function of $`L`$ and $`f`$ and we are reasonably successful in understanding the scaling exponents. In the third regime the noise is large enough to create small scale structures that are not neatly understood in terms of individual poles. It appears from our numerics that in this regime the roughening of the flame front gains a contribution from the the small scale structure in a way that is reminiscent of stable, noise driven growth models like the Kardar-Parisi-Zhang model.

One of our main motivations in this research was to understand the phenomena observed in radial geometry with expanding flame fronts. A full analysis of this problem cannot be presented here. We note however that many of the insights offered above translate immediately to that problem. Indeed, in radial geometry the flame front accelerates and cusps multiply and form a hierarchic structure as time progresses. Since the radius (and the typical scale) increase in this system all the time, new poles will be added to the system even by a vanishingly small noise. The marginal stability found above holds also in this case, and the system will allow the introduction of excess poles as a result of noise. The results discussed in Ref. can be combined with the present insights to provide a theory of radial growth. This theory will be offered in a forthcoming publication.

Finally, the success of this approach in the case of flame propagation raises hope that Laplacian growth patterns may be dealt with using similar ideas. A problem of immediate interest is Laplacian growth in channels, in which a finger steady-state solution is known to exist. It is documented that the stability of such a finger solution to noise decreases rapidly with increasing the channel width. In addition, it is understood that noise brings about additional geometric features on top of the finger. There are enough similarities here to indicate that a careful analysis of the analytic theory may shed as much light on that problem as on the present one.

Acknowledgments This work was supported in part by the German Israeli Foundation, the US-Israel Binational Science Foundation, the Minerva Center for Nonlinear Physics, and the Naftali and Anna Backenroth-Bronicki Fund for Research in Chaos and Complexity.