text:
# Ripples and Ripples: from Sandy Deserts to Ion-Sputtered Surfaces

## I Introduction

When an ion hits a surface liberates locally a large amount of energy that melts a region of the solid immediately below. For geometrical reasons, the sputtering effect depends on the surface-curvature: the energy concentrates on regions of positive curvature and this favorites the excavation of valleys and the growth of hills. On the other hand, thermal diffusion and surface tension tend to smoothen the irregularities by flattening the surface. It has been observed that under the combined action of these mechanisms the surface tends to create spontaneously ripples Ekl91 ; Krim93 ; Chra94 ; Valbusa02 ; Buat00 ; Rusp98 ; Rusp98a ; Rusp97 . In nature, ripples are commonly observed in sandy-deserts as the result of dynamical instability of the sand surface under the action of a sufficiently strong wind Bagnold . In this case, the formation of ripples is commonly associated with the effect produced by some grains that are lifted from the sand-bed and accelerated by the wind. These grains, when re-impact with the bed, splash up a number of other grains. Most of these grains return to the bed leading to a local rearrangement, whereas some other are accelerated by the wind and impact again on the bed after a certain ‘saltation’ length. In the literature, many studies have been devoted to understanding the mechanism of ripple formation Andersen02 ; Andersen01 ; Kurtze00 ; Hoyle99 ; Prigozhin99 ; Werner93 . In particular, an hydrodynamical model for aeolian ripple formation, based on a continuum dynamical description with two species of grains (immobile and rolling grains), was proposed with success by Bouchaud et al. Bouch95 ; Terz98 ; Val99 ; Csa00 . The main ingredient of such a model is a bilinear differential equation, for the population of the two species of grains, which shows the instability of a flat bed against ripple formation.

In this paper we show that the same reasoning which has been used to describe the sand ripples formation in deserts, applied to the studies of dynamical surface roughening, leads to an accurate description of the morphogenesis and evolution of ripples on crystal and amorphous surfaces during ion sputtering. The present approach contains not only the Bradley-Harper approach Bra88 ; Par99  (based on a Kardar-Parisi-Zhang type equation KPZ ), but it is also able to describe some of the crucial experimental features recently observed in these systems Rusp98 . In particular, by means of this approach we can explain the two distinct dynamical behaviors experimentally observed in amorphous/semiconductors systems and in metals Valbusa02 . In the first case (amorphous/semiconductors) we find that the ripples growth exponentially fast at constant wavelength $`\widehat{\lambda }`$ up to a critical roughening $`W_c`$ at which the growing process interrupts. On the other hand, in metals (when the Erlich-Schwoebel barrier is active), we find a transition between an initial exponential to a slower power-law growth of the roughness (the root mean square of the height profile). In this regime the ripple-wavelength tends also to growth with time.

## II Particle mobility and Ripple dynamics

When the surface of a solid is taken under ion sputtering some atoms in the proximity of the surface receive energy from the sputtered ions and pass from a bounded - *‘immobile’* \- solid state to a *‘mobile’* melted unbounded state. The opposite mechanism is also allowed: some *mobile* atoms can gain in energy by becoming *immobile* and bounding in a given position in the solid. A certain fraction of atoms might also be dispersed into the atmosphere. Let us call $`h(𝐫,t)`$ the height of surface profile made of immobile -bounded- atoms and call $`R(𝐫,t)`$ the height of mobile -melted- atoms. In analogy with the theory developed to explain the dynamical evolutions of dunes in deserts Bouch95 ; Terz98 ; Val99 ; Csa00 , we describe the mechanisms of excavation, exchange between mobile and immobile atoms and surface displacement of mobile atoms in term of the following differential equation:

$`{\displaystyle \frac{h}{t}}`$ $`=`$ $`\mathrm{\Gamma }(R,h)_{ex}+\mathrm{\Gamma }(R,h)_{ad}`$
$`{\displaystyle \frac{R}{t}}`$ $`=`$ $`𝐉(R,h)+(1\varphi )\mathrm{\Gamma }(R,h)_{ex}\mathrm{\Gamma }(R,h)_{ad}.`$ (1)
Where $`\mathrm{\Gamma }(R,h)_{ex}`$ is the rate of atoms that are excavated under the action of the sputtering, and $`(1\varphi )`$ is the part of them that pass from immobile to mobile, whereas $`\varphi `$ is the fraction that is dispersed into the atmosphere.

Let us now write in details the various terms contained in Eq.II.

### II.1 Excavation

The excavation effect must clearly depend on the number and velocity of the sputtered ions (i.e. its flux), but also the local shape and orientation of the surface might play an important role. Indeed, the energy transmitted by the impacting ions concentrate more in regions of the surface with positive curvature. Moreover, part of the surface facing the flux are likely to experience a different erosion respect to others which are less exposed to the flux. Crystalline orientation and anisotropies might be also taken into account. We can write:

$$\mathrm{\Gamma }(R,h)_{ex}=\eta (1+𝐚h+b^2h);$$
(2)
here $`\eta `$ is the sputtering flux, whereas $`𝐚`$ and $`b`$ are respectively associated with the flux-direction-dependent and with the curvature-dependent sputtering erosion.

### II.2 Adsorption

The rate of adsorption of mobile atoms into immobile solid positions must be dependent on the quantity of mobile atoms in a given spatial position. Similarly to the excavation process, the adsorption is also dependent on the local curvature and orientation. We can write:

$$\mathrm{\Gamma }(R,h)_{ad}=R(\gamma +𝐜h+d^2h),$$
(3)
where the parameter $`\gamma `$ is the recombination rate and $`𝐜`$ and $`d`$ are associated to the different probabilities of recombination in relation with the local orientation and shape of the surface.

Note that Eqs. 2, 3 contains the same terms as the ones proposed in the literature for the formation of aeolian dunes in the so-called hydrodynamical model Bouch95 ; Terz98 ; Val99 ; Csa00 ; Mak00 . Indeed, in deserts, sand grains are lifted from the sand-bed and readsorbed into it with a probability which is dependent on the local shape and orientation of the dunes. Eqs. 2, 3 represent the simplest analytical expressions which formally take into account these shape and orientation dependences. In the search for simple explanations, such equations are therefore rather universal.

### II.3 Mobility

Mobile atoms will move on the surface, and the quantity $`𝐉(R,h)`$ in Eq.II is the ‘current’ of these atoms. In surface growth, there are two main mechanisms that are commonly indicated as the responsible for the surface mobility of atoms Bar95 . The first is a current, driven by the variations of the local chemical potential, which tends to smoothen the surface asperity moving atoms from hills to valleys. The second is the current induced by the Erlich-Schwoebel barrier which -on the contrary- moves atoms uphill. In addiction to these main mechanisms we might also have to take into account a drift velocity and a random thermal diffusion, obtaining:

$$𝐉(R,h)=𝒦R(^2h)+sR\frac{h}{1+(\alpha _dh)^2}+𝐯R𝒟R.$$
(4)
In this equation, the first term describes a deterministic diffusion driven by the variations of the chemical potential which depends on the local shape of the surface; the second term is associated with the uphill current due to the Erlich-Schwoebel barrier and $`\alpha _d`$ is a constant associated with the characteristic length. The quantity $`𝐯`$ is a drift velocity of the mobile atoms on the surface, whereas $`𝒟`$ is the dispersion constant associated with the random thermal motion.

Note that Eq.4 is substantially different from the one proposed in the literature to describe ripples in granular media Bouch95 ; Terz98 ; Val99 ; Csa92 ; Csa93 . Here the current is supposed to be dependent on the local shape and orientation of the surface (the $`h(𝐫,t)`$ profile). The equations describing sand deserts can be retrieved from Eq.4 by imposing $`𝒦=0`$ and $`s=0`$, but -on the contrary- in surface growth these two parameters are the leading terms of the equation and play the role of control parameters in the dynamics of ripple formation. Nonetheless, these terms describe a rather simple dependence of the dynamic of particles on a surface on the geometrical shape of the surface itself. Again, in our seek for universality, we expect that similar terms can be profitably introduced in the context of aeolian sand ripples in order to describe specific phenomena (associated, for instance, with packing properties Aste00  or granular flow Komatsu01 ) which relate the current of grains with the dune-shapes.

It should be noted that the factors $`𝐚`$, $`𝐜`$ and $`𝐯`$ in Eqs.2,3 and 4 are *vectors* (i.e. they have -in general- different components in the two horizontal directions). Indeed, crystal surfaces are in general anisotropic and therefore one must take into account the dependence of the parameters on the relative orientation of the crystal-surface and the sputtering direction.

## III Dispersion Relation

A trivial solution of Eq.II can be written for a completely flat surface: $`h(𝐫,t)=h_0(t)`$ and $`R(𝐫,t)=R_0`$. In this case, we obtain $`R_0=(1\varphi )\eta /\gamma `$ and $`h_0(t)=\varphi \eta t+const.`$. This describes a surface that rests flat and it is eroded with a speed equal to $`\varphi \eta `$. But this behavior is only hypothetical since -in general- the dynamics of the surface-profile presents instabilities against spontaneous roughening and therefore its evolution is more complex. For instance, a numerical solution of Eq.II, is shown in Fig.2 (for the 1-dimensional case). We observe that, in a certain range of the parameters, the surface is unstable and periodic ripples are formed spontaneously.

### III.1 Stability analysis

In order to infer indications about the amplification or the smoothing of small perturbations and to deduce an analytical expression for the ripples wave-length at their beginning, we performe a stability analysis on Eq.II. For this purpose we assume that the surface-profile is made by the combination of a flat term plus a rough part:

$`R(𝐫,t)`$ $`=`$ $`R_0+R_1(𝐫,t)`$
$`h(x,t)`$ $`=`$ $`h_0(t)+h_1(𝐫,t),`$ (5)
with $`R_1(𝐫,t)=\widehat{R}_1\mathrm{exp}(i\omega t+i\mathrm{𝐤𝐫})`$ and $`h_1(𝐫,t)=\widehat{h}_1\mathrm{exp}(i\omega t+i\mathrm{𝐤𝐫})`$. We substitute these quantities into Eq.II and linearize the equation by taking only the first order in $`R_1`$ and $`h_1`$. A Fourier analysis (see Appendix A) shows that such a linearized equation admits solutions when the frequencies $`\omega `$ and the wave vectors $`𝐤`$ satisfy:

$`[i\omega +\gamma +i\mathrm{𝐤𝐯}+k^2𝒟]`$
$`\left\{i\omega +i𝐤\left[𝐯_1(1\varphi )𝐯_2\right]k^2\left[𝒟_1(1\varphi )𝒟_2\right]\right\}`$
$`\gamma (1\varphi )[i𝐤(𝐯_1𝐯_2)k^2\left(𝒟_1𝒟_2s_1\right)k^4𝒦_1]=0;`$
where, to simplify the equations, we have introduced the following notation:

| $`𝐯_1=\eta 𝐚`$ | $`𝐯_2=\eta 𝐜/\gamma `$ |
| --- | --- |
| $`𝒟_1=\eta b`$ | $`𝒟_2=\eta d/\gamma `$ |
| $`s_1=\eta s/\gamma `$ | $`𝒦_1=\eta 𝒦/\gamma `$ |

Equation III.1 establishes a *dispersion relation* $`\omega (𝐤)`$ that is a complex function with two branches corresponding to the solutions of the quadratic Equation III.1.

## IV Surface Instabilities

The kinetic growth of the surface instability is related to the immaginary part of $`\omega (𝐤)`$. Indeed, $`Im(\omega (𝐤))`$ corresponds to modes with amplitudes that change exponentially fast in time, and negative values correspond to unstable modes that increase with the time. We can therefore study $`Im(\omega )`$ from the solution of Eq.III.1 and search for the range of $`k`$ in which $`Im(\omega )`$ is negative. The most unstable mode is the one that grows faster and it corresponds to the value of $`𝐤`$ at which $`Im(\omega )`$ reaches its most negative value (see Fig.3).

The solution of Eq.III.1 for $`Im(\omega )`$, is

$`2Im(\omega )_\pm `$ $`=`$ $`\gamma +\left[𝒟𝒟_1+(1\varphi )𝒟_2\right]k^2`$ (7)
$`\pm `$ $`\sqrt{{\displaystyle \frac{\mathrm{\Delta }_1+(\mathrm{\Delta }_1^2+4\mathrm{\Delta }_2^2)^{1/2}}{2}}}`$
where we have

$`\mathrm{\Delta }_1`$ $`=`$ $`\gamma ^2\left\{\left[𝐯𝐯_1+(1\varphi )𝐯_2\right]𝐤\right\}^2`$
$`+`$ $`2\gamma \left[𝒟(12\varphi )𝒟_1+(1\varphi )𝒟_2+2(1\varphi )s_1\right]k^2`$
$`+`$ $`\{[(𝒟+𝒟_1(1\varphi )𝒟_2]^24\gamma (1\varphi )𝒦_1\}k^4`$
and

$`\mathrm{\Delta }_2`$ $`=`$ $`\gamma \left[𝐯+(12\varphi )𝐯_1(1\varphi )𝐯_2\right]𝐤`$
$`+`$ $`\left[𝒟+𝒟_1(1\varphi )𝒟_2\right]\left[𝐯𝐯_1+(1\varphi )𝐯_2\right]𝐤^3.`$
Let us first observe that in absence of sputtering (i.e. when $`\eta =0`$ and therefore, $`𝐯_1=0`$, $`𝐯_2=0`$, $`𝒟_1=0`$, $`𝒟_2=0`$, $`s_1=0`$, $`𝒦_1=0`$) the solutions of Eq.III.1 are $`\omega (𝐤)=0`$ and $`\omega (𝐤)=\mathrm{𝐤𝐯}+i(\gamma +k^2𝒟)`$. In this case, the imaginary part of $`\omega (𝐤)`$ is non-negative, therefore we -correctly- expect no spontaneous corrugation of the surface. On the contrary, when the sputtering is active ($`\eta 0`$), the immaginary part of $`\omega (𝐤)`$ can assume negative values as shown in Fig.3 where a plot of $`Im(\omega )_{}`$ is reported (along one direction of the vector $`𝐤`$). As one can see, typically the branch $`Im(\omega )_{}`$ takes negative values for $`k`$ between 0 and a critical value $`k^{}`$ at which it passes the zero typical . The critical point $`k^{}`$, fixes the minimal unstable wavelength. We therefore expect to find unstable solutions associated with the formation and evolution of ripples with wavelengths $`\lambda \lambda ^{}=2\pi /k^{}`$.

## V Ripple wavelength

Several analytical solutions of Eq.III.1 can be found in special cases which are discussed in Appendix B. But the study of the surface instabilities can be highly simplified if we consider the first order effects when the sputtering flux $`\eta `$ is small.

### V.1 Approximate equation

In the case of small spattering fluxes, the branch of $`Im(\omega (𝐤))`$, with negative values can be approximated to:

$$Im(\omega )_{}\frac{P_1k^6+P_2k^4+P_3(𝐤)k^2+P_4k^2+P_5(𝐤)}{𝒟^2k^4+2\gamma 𝒟k^2+(\mathrm{𝐯𝐤})^2+\gamma ^2}$$
(10)
with

$`P_1=𝒟[(1\varphi )\gamma 𝒦_1𝒟(𝒟_1(1\varphi )𝒟_2)]`$
$`P_2=(1\varphi )\gamma [𝒟(𝒟_2s_1)+\gamma 𝒦_1](1+\varphi )\gamma 𝒟𝒟_1`$
$`P_3(𝐤)=[𝒟_1+(1\varphi )𝒟_2](\mathrm{𝐯𝐤})^2`$
$`P_4=\gamma ^2[(1\varphi )s_1\varphi 𝒟_1]`$
$`P_5(𝐤)=(1\varphi )\gamma (\mathrm{𝐯𝐤})[(𝐯_1𝐯_2)𝐤]`$ (11)
When $`k`$ is sufficiently small ( $`k\gamma /\eta `$ ), we can develop Eq.10 at the 4th order obtaining:

$$Im(\omega )_{}Ak^4Bk^2,$$
(12)
with

$`A`$ $`=`$ $`(1\varphi )[𝒦_1+(s_1+𝒟_2𝒟_1){\displaystyle \frac{\gamma 𝒟+v^2}{\gamma ^2}}+`$
$`v(v_1v_2){\displaystyle \frac{2\gamma 𝒟+v^2}{\gamma ^3}}]`$
$`B`$ $`=`$ $`\varphi 𝒟_1+(1\varphi )\left[s_1+{\displaystyle \frac{v(v_1v_2)}{\gamma }}\right].`$ (13)
Here $`v`$, $`v_1`$ and $`v_2`$ are respectivelly the components of $`𝐯`$, $`𝐯_1`$ and $`𝐯_2`$ in the direction parallel to $`𝐤`$). (In Fig.3 a comparison between this approximate solution and the exact one is given.)

### V.2 Solutions

The expected wavelength of the ripples is associated with the fastest growing mode, which corresponds to the value of $`𝐤`$ at which $`Im(\omega )_{}`$ reaches its most negative point. Here the minimum of $`Im(\omega )`$ is at

$$\widehat{k}=\sqrt{\frac{B}{2A}}.$$
(14)
Therefore, at the beginning, the roughness will grow exponentially fast as $`W\mathrm{exp}(B^2t/(4A))`$ with associated ripple-wavelength at:

$$\widehat{\lambda }2\pi \sqrt{\frac{2A}{B}}.$$
(15)
### V.3 Some special cases

Let us first observe that, when $`𝒦_1`$, $`s_1`$ and $`\varphi =0`$, the ripple wavelength, given by Eq.15, coincides with the one found for sand dunes in deserts (see for instance Val99 ). In our notation the ‘reptation length’ is $`l_0=v/\gamma `$, the ‘cut-off length’ is $`l_c=(𝒟_2𝒟_1)/v`$, whereas $`v_1v_2`$ is the collective drift velocity of the dunes. The approximations usually applied in this context Terz98 ; Val99 , imply: $`l_c\sqrt{𝒟/\gamma }`$, and $`\gamma l_cv_1v_2`$. Giving, from Eq.15

$$\widehat{\lambda }2\pi \sqrt{\frac{2vl_0l_c}{v_1v_2}}.$$
(16)
Let us now consider the dynamical evolution of a surface under ion sputtering and in particular the case when the effect of the Erlich-Schwoebel barrier is not present (as for semiconductors and glasses). In this case, $`s=0,s_1=0`$ and we also expect that the drift velocity $`v`$ and the dispersion constant $`𝒟`$ are equal to zero or infinitesimally small. Indeed, here the current of mobile atoms on the surface is mainly induced by the differences in the chemical potential. Under these assumptions, from Eq.15, the wavelength of the most unstable ripple is:

$$\widehat{\lambda }2\pi \sqrt{\frac{2𝒦}{\nu }}.$$
(17)
where we called $`\nu =\gamma b\varphi /(1\varphi )`$, a quantity which plays the role of an effective surface tension. Note that Equation 17 is the same result as from the Bradley and Harper theory Bra88 ; Par99 ; Bar95 ; Cuer95 ; Gill01 .

When the Erlich-Schwoebel barriers are active ($`s,s_10`$), effects can be observed on the ripple-wavelength at their beginning, which becomes:

$$\widehat{\lambda }2\pi \sqrt{\frac{2𝒦}{\nu +s}}.$$
(18)
## VI Exponential/power-law growth and critical roughness

In metals, when the Erlich-Schwoebel barrier is active, there is an important non-linear contribution in the current of mobile atoms which becomes sizable when the roughness becomes sufficiently large and therefore $`(\alpha _dh)^21`$ (see Eq.4) where the average is over all the surface positions. We observe numerically that this changes the law of growth of ripples: from an exponential to a power-law behavior. This effect is shown in Fig.4. Numerically, all the computed exponents follow in the range between $`0.65`$ and $`0.85`$. The theoretical evaluation of this exponent is under current investigation. In this regime the ripple-wavelength tends also to grow with time. Experimentally, power law growth of the roughness and growth of characteristic wavelengths were observed in erosive sputtering on Ag(001) Valbusa02 .

In semiconductors or glasses, when no Erlich-Schwoebel barrier is present, it is physically intuitive that the exponential growth of the surface roughness (which is a characteristic of the beginning of the surface instability) cannot continue indefinitely. Indeed, from the expression $`R(𝐫,t)=R_0+\widehat{R}_1\mathrm{exp}(i\omega t+i\mathrm{𝐤𝐫})`$, which we used to derive Eq.15, we can immediately observe that when $`\widehat{R}_1>R_0=(1\varphi )\eta /\gamma `$, the amount of mobile atoms might become negative. Since a negative amount of atoms is physically impossible, the process of exponential roughness-growth described above must necessarily finish around a critical roughness given by:

$$W_c(1\varphi )\frac{\eta }{\gamma }.$$
(19)
This behavior is confirmed by numerical solutions of Eq.II and it is expected to be observable in semiconductors and glasses after sufficiently long times.

## VII Conclusions

We have shown that the same theoretical approach introduced to describe the formation of aeolian sand ripples can be conveniently applied to the study of the formation of periodic structures on surfaces under ion sputterning. Although the two phenomena are rather different, they can be described within the same conceptual framework by using rather general ideas that relate mobility, excavation and adsorption rates with the surface shape and orientation.

We have obtained general expressions for the ripples wavelength in term of the system parameters. It has been shown that in some particular cases such a solution coincide with the ones already known in the literature for sand dunes and surface instability Bra88 ; KPZ ; Par99 ; Terz98 ; Val99 . We have discussed the effect of the Erlich-Schwoebel barriers and compared the result with numerical solutions. We pointed out the Erlich-Schwoebel barrier can be responsible for a dramatic change in the system dynamics: from the exponential growth to a power law. Finally, we have demonstrated that the occurrence of a critical roughness is predictable within the present theoretical framework.

It should be noted that the main purpose of this paper is to point out a relevant example of universality: two processes which have completely different scales present a dynamical evolution which obeys to the same geometrical constraints and thus can be described by using the same phenomenological model. On the other hand, we must observe that the class of solutions of Eq.II is rich and complex - even in the linear approximation. Exhaustive, systematic studies of the classes of solutions of this equation and their dependence on the set of parameters will be the subject of future studies and publications.

## Appendix A Fourier transform of the linearized equation

By substituting Eqs.III.1, 4, 3 and 2 into Eq.II and by neglecting the second order terms (in $`R_1`$ and $`h_1`$), we obtain the following linearized equation:

$`{\displaystyle \frac{h_1}{t}}`$ $`=`$ $`\gamma R_1\left[𝐯_1(1\varphi )𝐯_2\right]h_1`$
$`\left[𝒟_1(1\varphi )𝒟_2\right]^2h_1`$
$`{\displaystyle \frac{R_1}{t}}`$ $`=`$ $`\gamma R_1𝐯R_1+𝒟^2R_1+`$ (20)
$`(1\varphi )[(𝐯_1𝐯_2)h_1+`$
$`(𝒟_1𝒟_2s_1)^2h_1𝒦_1^4h_1].`$
A Fourier analysis of Eq.A leads to

$`\gamma \widehat{R}_1\{i\omega +i𝐤[𝐯_1(1\varphi )𝐯_2]`$
$`k^2[𝒟_1(1\varphi )𝒟_2]\}\widehat{h}_1=0`$
$`[i\omega +\gamma +i\mathrm{𝐤𝐯}+k^2𝒟]\widehat{R}_1(1\varphi )[i𝐤(𝐯_1𝐯_2)`$
$`k^2(𝒟_1𝒟_2s_1)k^4𝒦_1]\widehat{h}_1=0,`$ (21)
with $`\widehat{R}_1`$ and $`\widehat{h}_1`$ the Fourier components of $`R_1`$ and $`h_1`$ respectively. This equation is a simple linear equation in two variables. It admits a non-trivial solution when the determinant of the coefficients is equal to zero. This leads to Eq.III.1.

## Appendix B Exact solutions

Analytical expressions for the value of $`𝐤`$ at which $`Im(\omega )=0`$ ($`𝐤^{}`$) can be calculated from Eq.7 in some special cases.

In particular, when $`\varphi =0`$, $`s_1=0`$, $`𝒦_1=0`$ and $`𝒟=0`$, we obtain

$$k^{}=\sqrt{\frac{\gamma (v_1v_2)}{(vv_1+v_2)(𝒟_2𝒟_1)}},$$
(22)
where $`v`$, $`v_1`$ and $`v_2`$ are the components of $`𝐯`$, $`𝐯_1`$ and $`𝐯_2`$ in the direction of $`𝐤^{}`$.

On the other hand when, $`𝒦_1`$, $`s_1`$, $`𝒟_1`$ and $`𝒟_2=0`$, we find

$$k^{}=\sqrt{\frac{\gamma (v\varphi v_1)}{𝒟(v_1v(1\varphi )v_2)}}.$$
(23)
The effect of the deterministic diffusion induced by the chemical potential can be studied from the solution

$$k^{}=\sqrt{\frac{\varphi \gamma 𝒟_1}{(1\varphi )\gamma 𝒦_1𝒟(𝒟_1(1\varphi )𝒟_2)}},$$
(24)
which holds when $`𝐯=0`$, $`𝐯_1=0`$, $`𝐯_2=0`$, $`s_1=0`$ and $`𝒟𝒟_1+(1\varphi )𝒟_2>0`$.

Whereas, when $`𝐯_1`$, $`𝐯_2`$, $`𝒟_1`$ and $`𝒟_2=0`$, we find

$$k^{}=\sqrt{\frac{s_1}{𝒦_1}},$$
(25)
which implies that the uphill current due to the Erlich-Schwoebel barrier can generate instability even when the shape-dependent erosion and recombination terms are inactive.

## Appendix C Numerical Solutions

The numerical solutions of Eq.II presented in this paper and in particular the ones shown in Figs.2 and 4 have been performed as follows. We considered a one-dimensional flat substrate ($`h(x,0)=h0`$) of length $`L`$, with periodic boundary conditions. An infinitesimal quantity of mobile atoms were added randomly to the substrate (with $`0<r(x,0)<L10^{10}`$). We then computed the profile-evolution using Eq.II with the derivative substituted with finite differences. To this purpose, substrate has been divided into $`N`$ discrete points. The -adimensional- time indicated in Fig.4 is the number of numerical steps. The height is in unit of $`L/N`$ and the roughness is defined as $`w(t,L)=[h(x^{},t)h(x,t)_x]^2_x^{}^{1/2}`$ (see, for instance, Csa92 ).

Several computations with a number of points equal to $`N=100`$, 200 and 300 (the one published have $`N=300`$) have been performed to verify the effect of boundary and discretization. Moreover, simulations with no periodic boundary conditions and with the sputtering term (Eq.2) applied only to a central mask, have also been performed obtaining very similar results. The robustness of the present approach has been verified varying the parameters, the time steps, the initial roughness of the substrate, etc. Comparable results have been always found but, we must stress that, under some conditions, numerical instabilities (in particular small surface-deformations with $`\lambda L/N`$) can be trigged on depending on the protocol utilized.


format: text

text:
# Transition Region Emission From Very Low Mass Stars

## 1 Introduction

Surface magnetic fields on low mass stars give rise to a variety of interesting astrophysical phenomena, including a persistent hot outer atmosphere (which may be contained within magnetic loop structures) and transient events such as starspots and flares. Non-radiative heating processes transfer energy from the magnetic field into the outer stellar atmosphere to power these phenomena, but a physical description of the heating processes remains elusive even on the Sun. Observationally, the heating produces chromospheric (T$`10^4`$K), transition region (T$`10^5`$K) and coronal (T$`10^6`$K) emission, primarily in emission lines of the hydrogen Lyman and Balmer series, and in ionized species of abundant elements such as helium, carbon, oxygen, calcium, magnesium, and iron.

Recently, there has been renewed interest in the atmospheric properties of the lowest mass stars and brown dwarfs. Chromospheric observations of H$`\alpha `$ emission (Hawley et al., 1996; Gizis et al., 2000; Burgasser et al., 2002) indicate that both the frequency of active stars at a given spectral type, and the activity level (measured by L$`{}_{H\alpha }{}^{}/L_{bol}`$) decrease at spectral types later than approximately M7. Recent investigations of magnetic activity in these very late type dwarfs have led to a new, yet already widely accepted, paradigm: that hot gas in the atmospheres of very low mass main sequence stars and brown dwarfs is present only during flares (Fleming et al., 2000; Rutledge et al., 2000; Mohanty et al., 2002). Here, hot gas refers to gas at transition region and coronal temperatures. The M9 brown dwarf, LP944-20, observed with Chandra by Rutledge et al. (2000) does show persistent chromospheric (H$`\alpha `$) emission during quiescence (i.e. outside of flares), as do many of its counterparts.

This new paradigm is based on the following evidence. Linsky et al. (1995) observed the M8 dwarf VB 10 with HST/GHRS and found no quiescent transition region emission, but did observe C IV emission during a flare. Fleming et al. (2000) re-analyzed ROSAT data on VB10 and found that the observed X-ray emission came only from a flare event, and no X-rays were detected during six hours of quiescent observation. The X-ray flare luminosity was substantially larger than expected based on the upper limit of the quiescent X-ray luminosity, when compared to flares on earlier type (M0 - M6) stars. An H$`\alpha `$ flare on an M9.5 2MASS dwarf showed a similarly substantial flare luminosity compared to quiescence (Liebert et al., 1999). Finally, the afore-mentioned Chandra investigation of the M9 brown dwarf LP944-20 detected no X-rays in a quiescent observation of nine hour duration (Rutledge et al., 2000). By contrast, Berger et al. (2001) found both persistent and flaring radio emission from LP944-20, while Berger (2002) recently reported radio observations of both flaring and quiescent emission for 3 additional low mass dwarfs of spectral types M8.5, M9.5 and L3.5. These authors speculate that perhaps the radio emission mechanism differs in the low mass dwarfs compared to earlier type active stars. On the theoretical side, both Fleming et al. (2000) and Mohanty et al. (2002) have proposed mechanisms related to the low ionization state of very cool atmospheres which would prevent the occurrence of magnetic activity in very cool (hence late spectral type) dwarfs, in order to explain the various observations described above.

Recent X-ray studies have probed into the brown dwarf regime in young clusters (Comerón, Neuhäuser, & Kaas, 2000; Feigelson et al., 2002) and find both quiescent and flaring emission. However, the extreme youth of these objects means that most have early-mid M spectral types. Further, fossil magnetic fields may play a significant role in their activity (Feigelson et al., 2002), making them unsuitable as a testbed for theories of the evolution of dynamo-generated magnetic fields in older, main sequence dwarfs.

Extensive observations of earlier M dwarfs (type M0-M6) have led to scaling relations between the persistent H$`\alpha `$ (chromospheric) and soft X-ray (thermal coronal) emission outside of flares (Hawley et al., 1996). We have used H$`\alpha `$ observations (Hawley et al., 1996) and C IV (transition region) measurements (Johns–Krull et al., 2000) for the dMe stars AD Leo (M3), EV Lac (M3.5), and YZ CMi (M4) to derive empirical C IV to H$`\alpha `$ line flux ratios of 0.06, 0.07, and 0.13 respectively. There may be a tendency for this ratio to increase toward later spectral types. A conservative estimate suggests that the ratio is typically $``$ 0.1. In terms of the observed fluxes, the relations between chromospheric, transition region, and coronal diagnostics for active early-mid M dwarfs are therefore:

$$F_{CIV}0.1F_{H\alpha }$$
$`(1)`$
$$F_X3F_{H\alpha }$$
$`(2)`$
If these relations are applied to the observed H$`\alpha `$ fluxes for VB 10 and LP944-20, the upper limits for the reported non-detections in C IV and X-rays are near or even above the predicted values. For example, Linsky et al. (1995) found an upper limit for VB 10 of $`F_{CIV}<6\times 10^{16}`$ ergs cm-2 s-1 while the scaling relation gives a predicted C IV line flux of $`9\times 10^{16}`$ ergs cm-2 s-1 based on the observed H$`\alpha `$ line flux from Hawley et al. (1996). Recently, Martín & Bouy (2002) reobserved LP944-20 with XMM-Newton, again failing to detect quiescent emission; however, this observation is only a factor of 3 more sensitive than the previous Chandra result of Rutledge et al. (2000). Thus, only slight deviations from the above relations would explain the current non-detections, without requiring a major paradigm shift with regard to magnetic activity on late M dwarfs compared to earlier M dwarfs.

The case for the absence of a quiescent transition region and corona in the lowest mass stars and brown dwarfs thus rests on tenuous observational evidence. The goal of this investigation is to place stringent limits on the presence of transition region emission using the sensitive HST/STIS spectrograph with the FUV-MAMA detector. These observations provide the first rigorous test of the presence of persistent quiescent emission on very late type main sequence dwarfs at temperatures well above chromospheric levels, and by inference the ubiquitous existence of a hot outer atmosphere.

## 2 Observations

Our target list comprises the brightest active stars known at spectral types M7 (VB 8), M8 (VB 10) and M9 (LHS 2065). These targets were chosen to give us the best possible detection limits for persistent quiescent emission. The total exposure time on each object corresponds to a C IV detection limit ten times fainter than the expected C IV line strength scaled from earlier type active M dwarfs (equation 1).

The Hubble Space Telescope was used with the STIS spectrograph and the FUV-MAMA detector to obtain 15 consecutive exposures of each target. All observations were made using the 52′′ long by 0.2′′ wide slit and the G140L grating to give a spectral resolution $`R\lambda /\delta \lambda 3600`$ at 1500Å. For VB 8, spectral observations started at 12:38 UT on 2001 October 6. Each G140L spectrum was integrated for approximately 5 minutes, and the total exposure sequence spanned two orbits. For LHS 2065, spectral observations started at 5:23 UT on 2002 April 27. Each G140L spectrum was integrated for approximately 10 minutes, and the total exposure sequence spanned four orbits. The observing sequence for VB 10 was identical to that of LHS 2065 with G140L spectral observations starting at 5:03 UT on 2002 September 7. Table 1 gives a detailed log of the STIS observations.

The data were downloaded from the HST archive, and data reduction was performed with the CALSTIS package (Lindler, 1999) utilizing the STIS pipeline. However, due to the relative faintness of our targets, the pipeline package was unable to properly locate the spectrum of the star on the detector and extract it. Therefore, we performed the spectral extraction from the flatfielded science image manually. To accomplish this, we coadded all 15 exposures for each star and examined the resulting image by eye, where it was easy to identify bright emission lines such as C IV and Si IV. We then used the order trace from the STIS pipeline, applying an offset such that the trace went through the emission lines from our object. Examining each of the 15 exposures individually, we verified that the order location did not change substantially (more than a pixel) during any one exposure. We then extracted the spectrum and background around this offset order location using the same extraction width (11 pixels) used by the pipeline to determine a count rate spectrum. The count rate spectrum was converted to a flux spectrum via the transformation provided in the CALSTIS pipeline reduction. We verified our procedure by extracting spectra around the nominal order locations provided by the CALSTIS pipeline for the 15 exposures of each source, recovering the same (noise) spectrum produced by the STIS pipeline to within the (somewhat large) photon statistics.

## 3 Results

The spectra extracted for each target from the sum of all 15 exposures are shown in Figure 1. The emission lines we have identified are typical of those found in active late-type stars (Ayres et al., 2001), and are very similar in relative strength to earlier type M dwarf transition region spectra, as shown by comparison to the recent STIS/FUV-MAMA spectrum of the dM1e star AU Mic (Pagano et al., 2000). The C IV doublet at 1550Å is the strongest feature typically observed from low mass stars in the wavelength range of study here and is also the strongest line observed in these 3 stars. For presentation purposes, we have applied Gaussian smoothing to the spectra with a FWHM of 3Å. This smoothing is not used in the flux measurements given in Table 1. The 1-$`\sigma `$ uncertainty limits are also indicated for the smoothed spectra in Figure 1.

Table 1 gives the observed C IV doublet line flux for each exposure in all three stars. The fluxes were determined by integrating the line over 6.4Å centered on the doublet; the integration limits were set by visual examination of the individual spectra. The significance of these detections can be estimated from the off-source background in the images, as discussed by Kraft et al. (1991); see also Babu & Feigelson (1996). The image of each STIS spectrum shows a clear variation in the background along the dispersion direction, likely due to scattered light from the strong Ly-$`\alpha `$ airglow present. The background shows no variation in the spatial axis at the position of the C IV line. To estimate the significance of the C IV detection in each exposure, we step above and below the position of the C IV line on the spectral image and extract the number of background counts in a box the same size as that used to extract the actual C IV line. This produces approximately 100 estimates of the background, to which we fit a Poisson distribution. We then compare the net number of counts detected in the C IV line to this distribution in order to determine the significance of the detection. In all but one exposure on LHS 2065, the C IV line is detected above the 3$`\sigma `$ level. For VB 8, the weakest detection is at the 3.1$`\sigma `$ level, the strongest at the 15.0$`\sigma `$ level, with a mean significance of 7.7$`\sigma `$. For VB 10, these three values are 3.3$`\sigma `$, 314$`\sigma `$, and 39.1$`\sigma `$ respectively. For LHS 2065, the corresponding detections are 2.9$`\sigma `$, 17.6$`\sigma `$, and 10.1$`\sigma `$. The continuum adjacent to the C IV line is not detected in the same sized integration box at a significant level in any exposure.

Although these estimates show that the C IV line is detected in all exposures, they do not represent the uncertainties in the measured flux. Using the Bayesian method described by Kraft et al. (1991), we have estimated the 68.3% confidence limits111The 68.3% confidence limits correspond to $`\pm 1\sigma `$ for a Gaussian distribution. We have not assumed a Gaussian distribution here, but these confidence limits represent a typical benchmark.. The flux confidence limits for each spectrum are reported in Table 1.

Time sequences of the C IV measurements, together with the uncertainties, for each target are shown in Figure 2. There is no definite sign of long-lived flaring (C IV rise and decay) in the time sequences observed. The second orbit on VB 10 shows a strong flare at the end of the orbit (exposure 7), so the duration of the event cannot be measured. The C IV flux does appear to be variable between exposures on all three stars. Quiescent emission does not typically occur at a constant flux level on active M dwarfs; similar variability is commonly seen in the quiescent H$`\alpha `$ emission fluxes (Gizis, Reid, & Hawley, 2002). This may be a sign of small scale flaring as a cause of the quiescent emission. Our C IV data suggest a common origin for the chromospheric and transition region emission, in conflict with the hypothesis that quiescent chromospheric emission is present but that transition region emission is absent except during strong flares.

The C IV line strengths we measure are: $`1.4\pm 0.1\times 10^{15}`$ ergs s-1 cm-2 for VB 8; $`2.8(1.4)\pm 0.1(0.1)\times 10^{15}`$ ergs s-1 cm-2 for VB 10; and $`7.1\pm 0.6\times 10^{16}`$ ergs s-1 cm-2 for LHS 2065. The values in parentheses for VB 10 represent the quiescent flux and uncertainty measured from the total of all spectra excluding the strong flare observed in exposure 7 (see Figure 2). Note that when integrated over the entire observation for each target, the C IV line is detected with sufficient counts that uncertainties can be estimated assuming a Gaussian distribution for the source and the background (Babu & Feigelson, 1996).

The scaling relation based on the H$`\alpha `$ emission in these stars, which assumes that the C IV to H$`\alpha `$ ratio is the same as for earlier type M dwarfs (equation 1), predicted values of 2 $`\times 10^{15}`$ (VB 8), 9 $`\times 10^{16}`$ (VB 10), and 1 $`\times 10^{15}`$ (LHS 2065) ergs s-1 cm-2. These are close to the observed values, with VB 8 and LHS 2065 having slightly lower observed fluxes, while the observed VB 10 (quiescent) flux is somewhat higher than predicted. The VB 10 (quiescent) flux is also a factor of two higher than the upper limit quoted in Linsky et al. (1995).

## 4 Summary

Our conclusion is that active, very low-mass, main-sequence stars (spectral types M7-M9) can maintain a quiescent chromosphere and transition region that is similar to those observed in active, earlier type M dwarfs. Flaring may still be ultimately responsible for the outer atmospheric heating in these stars, even during “quiescence”; the role of flares (particularly micro-flares) in coronal and transition region heating is still not understood for the Sun or for active stars (Kashyap et al., 2002; Güdel et al., 2003, for example). However, theories involving changes in atmospheric structure, dynamo heating and other mechanisms proposed to explain the lack of quiescent activity in very low mass stars and brown dwarfs (Fleming et al., 2000; Mohanty et al., 2002) must accommodate the persistent quiescent activity that we have observed at least through spectral type M9.

This research was made possible by grant HST-GO-9090 from the Space Telescope Science Institute, operated by AURA for NASA.


text:
# Critical field in a model with local pairs

## I Introduction

Almost twenty years after discovery of high–temperature superconductors (HTSCs) the mechanism responsible for their unusual properties remains unclear. The complex phase diagram of HTSCs suggests that there may be no single mechanism that dominates over the entire doping range. In particular, the normal–state properties in underdoped and overdoped regimes are different. Highly overdoped compounds in the normal state exhibit Fermi liquid behavior, whereas the superconducting state may be described within a weak–coupling BCS theoryFermi . On the other hand, in the underdoped regime the HTSCs exhibit unconventional features. The most remarkable of them are the extremely short coherence length and a pseudogap that opens in the normal state. The presence of the pseudogap has been confirmed with the help of various experimental techniques like: angle–resolved photoemission arpes1 ; arpes2 ; arpes3 , intrinsic tunneling spectroscopy krasnov1 ; krasnov2 , NMR nmr1 ; nmr2 , infrared infra  and transport transport  measurements. Although, there is no complete theoretical description of the pseudogap, one usually considers this phase as a precursor of the superconductivity. According to this hypothesis formation of Cooper pairs starts at temperature $`T^{}`$, higher than the superconducting transition temperature $`T_c`$. Then, at $`T_c`$, these preformed pairs undergo Bose-Einstein condensation.

This hypothesis seems to be supported by recent observations of the vortex–like Nernst signal above $`T_c`$ xu  that evolves smoothly into the analogous signal below the superconducting phase transition ochida . The Meissner effect does not occur in the pseudogap phase due to strong phase fluctuations rather than the vanishing of the superfluid density. Therefore, theoretical description of the suppression of the Meissner effect requires an approach beyond the mean–field level. Despite the absence of the Meissner effect above $`T_c`$, one can observe inhomogeneous magnetic domains that are interpreted as precursors to the Meissner state iguchi .

The short coherence length indicates that the pairing takes place in the real space, leading to boson–like objects. A few models are commonly used to describe systems with the local pairs. Namely, the attractive Hubbard (AH) model micnas , fermion–bosonFB  and purely bosonic modelsboson , as well as the Penson–Kolb (PK) modelPK , i.e., the tight–binging model with local pair hopping. These models should be considered as effective approaches which do not explain the microscopic origin of the pairing interaction.

Another unusual property of HTSCs is related to their behavior in the external magnetic field. In particular, temperature dependence of the upper critical field $`H_{c2}`$ has a positive curvatureosofsky ; mackenzie  in contradistinction to classical superconductors, where a negative curvature is observed. Moreover, $`H_{c2}`$ does not saturate even at genuinely low temperature. Recent experimentshpg  show that also the pseudogap is destroyed by sufficiently high magnetic field, $`H_{pg}`$. Although, the temperature dependence of $`H_{pg}`$ has a negative curvature, it significantly differs from the predictions of the standard Helfand–Werthammer theoryHW . Namely, $`H_{pg}(T)`$ has a large slope at temperatures close to $`T^{}`$ and saturates already at $`T0.7T^{}`$. These features may assist in verification of the preformed Cooper pairs hypothesis and, more generally, in choosing the most appropriate model of HTSC.

## II Model

In the present paper we show that opening of the pseudogap and its dependence on the magnetic field can be described within a model with local pair hopping. Our starting point is the two–dimensional (2D) Penson–Kolb model with the Hamiltonian given by:

$`H`$ $`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}t_{ij}\mathrm{e}^{i\mathrm{\Phi }_{ij}}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{i\sigma }{}}(g\mu _BH_z\sigma \mu )c_{i\sigma }^{}c_{i\sigma }`$ (1)
$`{\displaystyle \frac{1}{2}}J{\displaystyle \underset{i,j}{}}\mathrm{e}^{2i\mathrm{\Phi }_{ij}}c_i^{}c_i^{}c_jc_j.`$
Here, $`c_{i\sigma }^{()}`$ creates (annihilates) an electron with spin $`\sigma `$ at site $`i`$, $`t_{ij}`$ is the single electron hopping integral between sites $`i`$ and $`j`$, $`\mu `$ is the chemical potential and $`J`$ is the nearest neighbor pair hopping interactions. The external magnetic field perpendicular to the lattice $`H_z`$ shifts the energy levels by $`g\mu _BH_z\sigma `$ ($`g`$ is the gyromagnetic ratio and $`\mu _B`$ is the Bohr magneton) and modifies the hopping terms. The single electron hopping integral acquires the Peierls factor

$$\mathrm{\Phi }_{ij}=\frac{e}{\mathrm{}c}_{\text{R}_j}^{\text{R}_i}\text{A}𝑑\text{l},$$
(2)
whereas the phase factor in the pair hopping term is twice larger.

The Penson–Kolb model can be derived from a general microscopic tight–binding Hamiltonianhubbard , where the Coulomb repulsion may lead to the pair hopping interaction. In such a case $`J`$ is negative (repulsive Penson–Kolb model). However, treating the Penson–Kolb model as of a phenomenological nature, we assume $`J`$ to be an effective parameter, that can be negative as well as positive. It can be understood as a result of renormalization originating, e.g., from electron–phonon couplingrob-bul . For a nonzero single electron hopping integral $`JJ`$ is not a symmetry of the PK model japaridze . However, superconducting correlations occur in the Penson–Kolb model for attractive pair hopping interaction ($`J>0`$) as well as for the repulsive one ($`J<0)`$, provided that the pair hopping is large enough. The latter case is usually referred to as $`\eta `$–type pairing. Then, the total momentum of the paired electrons is $`\text{Q}=(\pi ,\pi )`$ and the phase of superconducting order parameter alters from one site to the neighboring one. It has been shown that there is a flux quantization and Meissner effect in this state.yang  Superconductivity survives also in the presence of on–site Coulomb repulsion (Penson–Kolb–Hubbard model), provided that this interaction is not too largeFab .

### II.1 Density of states

At the mean–field level, for $`J>0`$ one obtains an isotropic superconducting gap, identical to that obtained for AH model. On the other hand, in the case of $`\eta `$–type pairing ($`J<0`$), the density of states is finite for arbitrary energy. However, the density of states at the Fermi level may significantly be suppressed for some dopings. In the simplest case of the nearest neighbor hopping the density of states in the $`\eta `$–phase is of the form:

$$\rho (\omega )=\frac{1}{2}\left(1\frac{\mu }{\stackrel{~}{\mu }}\right)\rho _0(\omega \stackrel{~}{\mu })+\frac{1}{2}\left(1+\frac{\mu }{\stackrel{~}{\mu }}\right)\rho _0(\omega +\stackrel{~}{\mu }).$$
(3)
Here, $`\stackrel{~}{\mu }=\sqrt{\mu ^2+4|J\mathrm{\Delta }|^2}`$, $`\mathrm{\Delta }(1)^ic_ic_i`$ is the $`\eta `$–phase order parameter and $`\rho _0`$ is the density of states for $`\mu =\mathrm{\Delta }=0`$. One can see from Eq. (3) that the quasiparticle poles split when $`\mathrm{\Delta }`$ becomes finite. Therefore, a local minimum in the density of states may occur at the Fermi surface. Despite the presence of this minimum the density of states at the Fermi level remains finite provided that $`\stackrel{~}{\mu }`$ is small when compared to the band width. Inclusion of the next nearest neighbor hopping $`t^{}`$ leads to a more complicated expression for the density of states. However, the structure of $`\rho `$ in the $`\eta `$–phase remains unchanged. Fig. 1. shows the density of states calculated for $`t^{}0`$ and different values of the $`\eta `$–phase order parameter.

Gradual decreasing of $`\rho `$ at the Fermi level resembles opening of the pseudogap in HTSCs.

Another feature that could speak in favor of this interpretation is anisotropy of the gap arpes2 . More precisely, for $`t^{}0`$ the magnitude of splitting of the quasiparticle peaks depends on the direction in the Brillouin zone. The splitting of the spectral functions is presented in Fig. 2. As we consider isotropic order parameter the splitting is finite everywhere at the Fermi level, in contradistinction to a purely $`d`$–wave gap. However, this drawback may be removed when considering a nonlocal pairing.

### II.2 Response to magnetic field

In contradistinction to the AH model, the external magnetic field explicitly enters the term responsible for superconductivity, i.e., the pair hopping interaction. Therefore, the differences between AH and PK models may show up in the electromagnetic properties czart . Here, we investigate the temperature dependence of critical field $`H_{\mathrm{crit}}`$. It is defined as the highest magnitude of the magnetic field, for that there exists a non–zero solution for the superconducting order parameter:

$$\mathrm{\Delta }_i=c_ic_i.$$
(4)
As we carry out calculations at the mean–field level it is impossible to determine the phase coherence of the Cooper pairs. Therefore, the physical interpretation of the critical field is not unique. In the overdoped regime of cuprates it can directly correspond to the upper critical field. On the other hand, for underdoped systems, and within the precursor scenario of pseudogap, it can be interpreted as $`H_{pg}`$, i.e., the field at which the incoherent pairs appear. One can argue that the PK model can also be used as an effective model of other superconductors that are characterized by a short coherence length: barium bismuthates, fullerides, Chevrel phases, organic superconductors and heavy fermion systems. czart  In these cases $`H_{\mathrm{crit}}`$ should be considered as $`H_{c2}`$.

For the sake of simplicity we define:

$$\stackrel{~}{\mathrm{\Delta }}_i=\frac{1}{2}\underset{j}{}^{}\mathrm{e}^{2i\mathrm{\Phi }_{ij}}\mathrm{\Delta }_j,$$
(5)
where the prime means that the summation is carried out over the nearest neighbors of site $`i`$. Then, the mean–field Hamiltonian takes on the following form:

$`H`$ $`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}t_{ij}\mathrm{e}^{i\mathrm{\Phi }_{ij}}c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{i\sigma }{}}(g\mu _BH_z\sigma \mu )c_{i\sigma }^{}c_{i\sigma }`$ (6)
$`J{\displaystyle \underset{i}{}}(c_i^{}c_i^{}\stackrel{~}{\mathrm{\Delta }}_i+\mathrm{H}.\mathrm{c}.).`$
At the mean–field level the only difference between PK and AH models is the presence of $`\stackrel{~}{\mathrm{\Delta }}_i`$ in Eq. (6) instead of $`\mathrm{\Delta }_i`$. Therefore, in order to calculate the critical field one can follow an approach, that has previously been developed for the lattice gas with on–site attraction my . Then, one ends up with the lattice version of the Gor’kov equations:

$$\mathrm{\Delta }_i=\frac{J}{\beta }\underset{j,\omega _n}{}\stackrel{~}{\mathrm{\Delta }}_jG(i,j,\omega _n)G(i,j,\omega _n).$$
(7)
Here, $`G(i,j,\omega _n)`$ is the one–electron Green’s function in the presence of a uniform and static magnetic field. $`\omega _n`$ denotes the fermionic Matsubara frequency. With the help of Eq.(5) one can eliminate $`\stackrel{~}{\mathrm{\Delta }}_i`$ from the Gor’kov equations. Then, $`H_{\mathrm{crit}}(T)`$ can be calculated from

$$\mathrm{\Delta }_i=\frac{J}{2\beta }\underset{j,l,\omega _n}{}\mathrm{e}^{2i\mathrm{\Phi }_{lj}}\mathrm{\Delta }_jG(i,l,\omega _n)G(i,l,\omega _n)$$
(8)
or

$$\mathrm{\Delta }_i=\frac{J}{2\beta }\underset{j,l,\omega _n}{}(1)^{j+l}\mathrm{e}^{2i\mathrm{\Phi }_{lj}}\mathrm{\Delta }_jG(i,l,\omega _n)G(i,l,\omega _n).$$
(9)
Eqs. (8) and (9) are equivalent since $`(1)^{j+l}=1`$, for the neighboring sites $`j`$ and $`l`$. However, it is more convenient to use the first/second of them for attractive/repulsive pair hopping interaction.

In the following we consider only the nearest–neighbor one–particle hopping integral $`t`$ and use the Landau gauge $`\text{A}=H_z(0,x,0)`$. Then, the Harper equation

$`g(\overline{p}_x,p_y,x+1)+2\mathrm{cos}\left(hxp_ya\right)g(\overline{p}_x,p_y,x)`$
$`+g(\overline{p}_x,p_y,x1)=t^1E(\overline{p}_x,p_y)g(\overline{p}_x,p_y,x),`$ (10)
determines eigenvalues $`E(\overline{p}_x,p_y)`$ of the one–particle hopping term. The corresponding eigenstates are enumerated by $`\overline{p}_x,p_y`$ and are of the form:

$$U_{x,y}(\overline{p}_x,p_y)=e^{ip_yya}g(\overline{p}_x,p_y,x).$$
(11)
Here, $`x,y`$ are integers which enumerate the lattice sites in $`\widehat{x}`$ and $`\widehat{y}`$ directions, whereas $`h/(2\pi )`$ is a ratio of the flux through a lattice cell to one flux quantum. We refer to Ref. my  for the details.

The one–electron Green’s function can be expressed with the help of eigenvalues and eigenstates of the normal–state Hamiltonian. Then, the summation over Matsubara frequencies in Eqs. (8) and (9) can explicitly be carried out. In the Landau gauge the presence of magnetic field does not change the plane–wave behavior in $`\widehat{y}`$–direction \[see Eq. (11)\]. Therefore, the superconducting order parameter depends only on $`x`$ and pairing of electrons takes place for the same $`\widehat{y}`$–components of their momenta, as in the absence of magnetic field, i.e., $`(p_y,p_y)`$ for $`J>0`$ and $`(p_y,\pi p_y)`$ for $`J<0`$. Taking these features into account one can rewrite the Gor’kov equations for the attractive:

$`\mathrm{\Delta }_x^{}`$ $`=`$ $`{\displaystyle \frac{J}{2\sqrt{N}}}{\displaystyle \underset{x}{}}\mathrm{\Delta }_x{\displaystyle \underset{p_y,\overline{p}_x,\overline{k}_x}{}}\chi (\overline{p}_x,p_y,\overline{k}_x,p_y)`$ (12)
$`\times [2\mathrm{cos}(2hx)g(\overline{p}_x,p_y,x)g(\overline{k}_x,p_y,x)`$
$`+g(\overline{p}_x,p_y,x+1)g(\overline{k}_x,p_y,x+1)`$
$`+g(\overline{p}_x,p_y,x1)g(\overline{k}_x,p_y,x1)]`$
$`\times g(\overline{p}_x,p_y,x^{})g(\overline{k}_x,p_y,x^{}),`$
as well as for the repulsive pair hopping interaction:

$`\mathrm{\Pi }_x^{}`$ $`=`$ $`{\displaystyle \frac{J}{2\sqrt{N}}}{\displaystyle \underset{x}{}}\mathrm{\Pi }_x{\displaystyle \underset{p_y,\overline{p}_x,\overline{k}_x}{}}\chi (\overline{p}_x,p_y,\overline{k}_x,\pi p_y)`$ (13)
$`\times [2\mathrm{cos}(2hx)g(\overline{p}_x,p_y,x)g(\overline{k}_x,\pi p_y,x)`$
$`g(\overline{p}_x,p_y,x+1)g(\overline{k}_x,\pi p_y,x+1)`$
$`g(\overline{p}_x,p_y,x1)g(\overline{k}_x,\pi p_y,x1)]`$
$`\times g(\overline{p}_x,p_y,x^{})g(\overline{k}_x,\pi p_y,x^{}).`$
Here $`\mathrm{\Pi }_x\mathrm{\Delta }_{x,y}(1)^y`$ and the Cooper pair susceptibility $`\chi (\overline{p}_x,p_y,\overline{k}_x,k_y)`$ has the following form:

$`\chi (\overline{p}_x,p_y;\overline{k}_x,k_y)`$ $`=`$ $`[\mathrm{tanh}{\displaystyle \frac{E(\overline{p}_x,p_y)\mu g\mu _BH_z}{2k_BT}}`$
$`+`$ $`\mathrm{tanh}{\displaystyle \frac{E(\overline{k}_x,k_y)\mu +g\mu _BH_z}{2k_BT}}]`$
$`\times `$ $`\left[2\left(E(\overline{p}_x,p_y)+E(\overline{k}_x,k_y)2\mu \right)\right]^1.`$
The above equations determine the strength of magnetic field at which the local pairing disappears. We have carried out calculations for $`150\times 150`$ cluster with periodic boundary conditions (bc) along the $`\widehat{y}`$ axis. As the Landau gauge breaks the translation invariance along $`\widehat{x}`$ axis we have used fixed bc in this direction. Our previous calculations indicate that such a size of cluster is sufficient to obtain convergent results.my

Fig. 3 shows the temperature dependence of $`H_{\mathrm{crit}}`$ obtained for attractive and for the repulsive pair hopping interaction. These results are compared with $`H_{\mathrm{crit}}(T)`$ calculated from the 2D AH modelmy  with $`U=t`$. We have adjusted the strength of the pair hopping interaction to obtain the same critical temperature in the absence of magnetic field. For $`J>0`$ $`H_{\mathrm{crit}}(T)`$ in PK model is very close to that of AH model. It means that in the case of $`s`$–wave pairing the Peierls factor in the pair hopping term leads only to a small decrease of superconducting correlations. However, the temperature dependence of $`H_{\mathrm{crit}}`$ in the $`\eta `$–state differs qualitatively from the $`s`$–wave case. Namely, $`H_{\mathrm{crit}}(T)`$ has a very large slope for a weak magnetic field and saturates already at relatively high temperature. Such a behavior of the critical field resembles $`H_{pg}(T)`$, that has recently been observed in $`\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_20_{8+\mathrm{y}}`$.hpg  In Fig. 4 we compare our results and the experimental data.

Within the Helfand–Werthammer theory, the temperature dependence of critical field is predominantly determined by the diamagnetic pair breaking mechanism. The Zeeman coupling becomes important only for sufficiently strong magnetic field. This feature holds also in the case of the lattice gas.maciek  In order to investigate the role of the Zeeman and orbital contributions in PK model, we have repeated our calculations in the absence of the Zeeman term. The resulting $`H_{\mathrm{crit}}(T)`$ is shown in Fig. 5.

In contradistinction to the $`s`$–wave superconductivity, the diamagnetic pair breaking is of minor importance in the case of $`\eta `$–pairing. This feature is responsible for extremely high values of $`H_{\mathrm{crit}}`$ in the absence of Zeeman term. Experimental investigationshpg  show that the pseudogap closing field scales linearly with $`T^{}`$. In Ref. hpg  the value of the scaling factor has been interpreted in favor of the Zeeman coupling as a mechanism that closes the pseudogap. This pair breaking mechanism dominates also in the case of $`\eta `$–pairing.

## III Summary

To summarize, we have investigated the PK model with attractive as well as repulsive pair hopping interaction. We have shown that the repulsive pair hopping term may lead to the occurrence of local minimum in the density of states, that is characteristic for pseudogap phase of underdoped cuprates. It originates from the spitting of the quasiparticle peaks. Despite the on–site pairing the magnitude of the splitting is a direction–dependent quantity, provided that $`t^{}0`$. Anisotropy of the pseudogap is observed in ARPES experiments.arpes2  We have also calculated the temperature dependence of $`H_{\mathrm{crit}}`$, defined as the highest magnetic field for which there exists a non–zero solution for the order parameter. We have found that in the case of $`\eta `$–type pairing $`H_{\mathrm{crit}}(T)`$ reproduces the experimental data for the pseudogap closing field. These features do not occur for attractive pair hopping interaction. In this case the gap structure as well as $`H_{\mathrm{crit}}(T)`$ are similar to those obtained for AH model.

Our approach to the critical field accounts both for Zeeman and diamagnetic pair breaking mechanisms. In the case of $`s`$–wave pairing inclusion of the Zeeman coupling does not lead to any essential changes in $`H_{\mathrm{crit}}(T)`$. On the other hand, Zeeman term is of crucial importance for $`\eta `$–pairing, whereas the diamagnetic pair breaking is ineffective. According to the experimental data the pseudogap is closed by the Zeeman splitting.

As we have previouslymy  shown, other models appropriate for short coherence superconductors have ground states (with $`s`$–wave or $`d`$–wave symmetry) that are almost insensitive to the Zeeman interaction. Therefore, the PK model with $`J<0`$ is unique in that the gap is closed predominantly due to the Zeeman interaction.

Collecting the features: the presence of the pseudogap, its anisotropy, Zeeman origin of $`H_{\mathrm{crit}}`$ (in agreement with the experimental data), the presence of flux quantization and the Meissner effect, (consistent with the preformed Cooper pairs scenario), may lead to a tempting hypothesis that the pair hopping can be responsible for the pseudogap. However, in order to avoid the problem of interpretation of the critical field, it should be verified beyond the mean–field level, discussed in this paper.

###### Acknowledgements.

We acknowledge a stimulating discussion with S. Robaszkiewicz. This work was supported in part by the Polish State Committee for Scientific Research, Grant No. 2 P03B 050 23.