text:
# Navier-Stokes-alpha model: LES equations with nonlinear dispersion

## 1 Introduction

The effects of subgrid-scale (SGS) fluid motions occurring below the available resolution of numerical simulations must be modeled. One way of modeling these effects is simply to discard the energy that reaches such subgrid scales. This is clearly unacceptable, though, and many creative alternatives have been offered. A prominent example is the large eddy simulation (LES) approach, see, e.g.,  \- . The LES approach is based on applying a spatial filter to the Navier-Stokes equations. The reduction of flow complexity and information content achieved in the LES approach depends on the characteristics of the filter that one uses, its type and width. In particular, the LES approach introduces a length scale into the description of fluid dynamics, namely, the width of the filter used. The LES approach is conceptually different from the Reynolds Averaged Navier-Stokes (or, RANS) approach, which is based on statistical arguments and exact ensemble averages, rather than spatial and temporal filtering. After filtering, however, just as in the RANS approach, one faces the classic turbulence closure problem: How to model the effects of the filtered-out subgrid scales in terms of the remaining resolved fields? In practice, this problem is compounded by the requirement that the equations be solved numerically, thereby introducing further approximations.

A turbulence modeling scheme — called here the Navier-Stokes-alpha model, or NS$`\alpha `$ model (also called the viscous Camassa-Holm equations in  \- ) — was recently introduced. This model imposes an energy “penalty” inhibiting the creation of smaller and smaller excitations below a certain length scale (denoted alpha), but still allows for nonlinear sweeping of the smaller scales by the larger ones. This energy penalty implies a nonlinearly dispersive modification of the Navier-Stokes equations. The alpha-modification appears in the nonlinear convection term, it depends on length scale and we emphasize that it is dispersive, not dissipative. This modification causes the translational kinetic energy wavenumber spectrum of the NS$`\alpha `$ model to roll off rapidly below the length scale alpha as $`k^3`$ in three dimensions, instead of continuing to follow the slower Kolmogorov scaling law, $`k^{5/3}`$, thereby shortening its inertial range and making it more computable .

This roll-off in the energy spectrum means that the inertial range is curtailed and the length scale $`\mathrm{}_\alpha `$ at which viscous dissipation takes over in the NS$`\alpha `$ model is larger than for the Navier-Stokes equations with the same viscosity. Hence, for a given driving force and viscosity, the number of active degrees of freedom $`N_{dof}^\alpha `$ for the NS$`\alpha `$ model is smaller than for Navier-Stokes. A rigorous analytical estimate was derived in  for the fractal dimension $`D_{frac}`$ of the global attractor for the NS$`\alpha `$ model. Namely,

$$D_{frac}(N_{dof}^\alpha )^{3/2}\text{and}N_{dof}^\alpha (L/\mathrm{}_\alpha )^3\frac{L}{\alpha }Re^{3/2},$$
(1)
where $`L`$ is the integral scale (or domain size), $`\mathrm{}_\alpha `$ is the end of the NS$`\alpha `$ inertial range and $`Re=L^{4/3}\epsilon _\alpha ^{1/3}/\nu `$ is the Reynolds number (with NS$`\alpha `$ energy dissipation rate $`\epsilon _\alpha `$ and viscosity $`\nu `$). The Kolmogorov dissipation length scale is $`\mathrm{}_{Ko}=(\nu ^3/\epsilon _\alpha )^{1/4}`$. The three length scales $`\alpha >\mathrm{}_\alpha >\mathrm{}_{Ko}`$ are related by

$$\mathrm{}_\alpha ^{\mathrm{\hspace{0.17em}3}}\alpha \mathrm{}_{Ko}^{\mathrm{\hspace{0.17em}2}}\text{or}\frac{\alpha }{\mathrm{}_\alpha }\frac{\mathrm{}_\alpha ^{\mathrm{\hspace{0.17em}2}}}{\mathrm{}_{Ko}^{\mathrm{\hspace{0.17em}2}}}.$$
(2)
Thus, these length scales $`\alpha >\mathrm{}_\alpha >\mathrm{}_{Ko}`$ stand in relation to each other in the same way as the length scales $`(K^3/\epsilon ^2)^{1/2}>(K\nu /\epsilon )^{1/2}>(\nu ^3/\epsilon )^{1/4}`$ in a K-epsilon model of Navier-Stokes turbulence.

The number of degrees of freedom for a corresponding Navier-Stokes flow with the same viscosity and energy dissipation rate is

$$N_{dof}^{NS}(L/\mathrm{}_{Ko})^3Re^{9/4},$$
(3)
The implication of these estimates of degrees of freedom for direct numerical simulations that access a significant number of them using the NS$`\alpha `$ model is an increase in computational speed relative to Navier-Stokes of

$$\left(\frac{N_{dof}^{NS}}{N_{dof}^\alpha }\right)^{4/3}\left(\frac{\alpha }{L}\right)^{4/3}Re.$$
(4)
Thus, if $`\alpha `$ tends to a constant value, say $`L/100`$, when the Reynolds number increases – as found in  \-  by comparing steady NS$`\alpha `$ solutions with experimental data for turbulent flows in pipes and channels – then one could expect to obtain a substantial increase in computability by using the NS$`\alpha `$ model at high Reynolds numbers. An early indication of the reliability of using these estimates to gain a relative increase in computational speed in direct numerical simulations of homogeneous turbulence in a periodic domain is given in .

In this paper, we introduce the NS$`\alpha `$ model as one of a class of equations obtained by filtering the Kelvin circulation theorem for Navier-Stokes. These equations possess the vortex stretching and transport properties emphasized in . We then compare a spatially filtered version of the NS$`\alpha `$ equations with a generalization of the LES similarity models. We also calculate the kinematic tranformation properties of the NS$`\alpha `$ model under changes of reference frames and discuss its nonlinear spectral-mode dynamics.

## 2 Kelvin-filtered turbulence models

Following , the Navier-Stokes-alpha (NS$`\alpha `$) model may be introduced by starting from Kelvin’s circulation theorem for the original Navier-Stokes (NS) equations,

$$\frac{𝐯}{t}+𝐯𝐯+p=\nu ^2𝐯+𝐟,\text{with}𝐯=0,$$
(5)
for an appropriate forcing $`𝐟`$ and constant kinematic viscosity $`\nu `$. These equations satisfy Kelvin’s circulation theorem,

$$\frac{d}{dt}_{\gamma (𝐯)}𝐯𝑑𝐱=_{\gamma (𝐯)}(\nu ^2𝐯+𝐟)𝑑𝐱,$$
(6)
for a fluid loop $`\gamma (𝐯)`$ that moves with velocity $`𝐯(𝐱,t)`$, the Eulerian fluid velocity.

### Kelvin-filtering the Navier-Stokes equations.

The equations for the NS$`\alpha `$ model may be introduced by modifying the NS Kelvin circulation theorem (6) to integrate around a loop $`\gamma (\widehat{𝐯})`$ that moves with a spatially filtered Eulerian fluid velocity given by $`\widehat{𝐯}=g𝐯`$, where $``$ denotes the convolution,

$$\widehat{𝐯}=g𝐯=g(𝐱𝐱^{})𝐯(𝐱^{})d^{\mathrm{\hspace{0.17em}3}}x^{}L_g𝐯(𝐱).$$
(7)
The “inverse” operation is denoted

$$𝐯=𝒪\widehat{𝐯}L_g^1\widehat{𝐯},$$
(8)
thereby defining an operator $`𝒪`$ whose Green’s function is the filter $`g`$ and which we shall assume is positive, symmetric, isotropic, translation-invariant and time-independent. Under these assumptions the quantity (kinetic energy)

$$E=\frac{1}{2}𝐯\widehat{𝐯}d^3x=\frac{1}{2}𝐯g𝐯d^3x=\frac{1}{2}\widehat{𝐯}𝒪\widehat{𝐯}d^3x,$$
(9)
defines a norm.

We modify the Navier-Stokes equations (5) by replacing in their Kelvin’s circulation theorem (6) the loop $`\gamma (𝐯)`$ with another loop $`\gamma (\widehat{𝐯})`$ moving with the spatially filtered velocity, $`\widehat{𝐯}`$. Then we have,

$$\frac{d}{dt}_{\gamma (\widehat{𝐯})}𝐯𝑑𝐱=_{\gamma (\widehat{𝐯})}(\nu ^2𝐯+𝐟)𝑑𝐱.$$
(10)
After taking the time derivative inside the Kelvin loop integral moving with filtered velocity $`\widehat{𝐯}`$ and reconstructing the gradient of pressure, we find the Kelvin-filtered Navier-Stokes equation,

$$\frac{𝐯}{t}+\widehat{𝐯}𝐯+\widehat{𝐯}^\mathrm{T}𝐯+p=\nu ^2𝐯+𝐟,$$
(11)
with

$$\widehat{𝐯}=0,\text{and}𝐯=𝒪\widehat{𝐯}.$$
(12)
The transport velocity $`\widehat{𝐯}(𝐱,t)`$ is the spatially filtered Eulerian fluid velocity in Eq. (7). Note that the continuity equation is now imposed as $`\widehat{𝐯}=(g𝐯)=0`$. If filtering commutes with differentiation, then $`g(𝐯)`$ also vanishes and the transported velocity $`𝐯`$ is also divergenceless. The energy balance relation derived from the NS$`\alpha `$ equation (11) is

$$\frac{d}{dt}\frac{1}{2}\widehat{𝐯}𝐯d^3x=\widehat{𝐯}𝐟d^3x\nu \mathrm{tr}(\widehat{𝐯}^T𝐯)d^3x,$$
(13)
where we have used incompressibility and have dropped any boundary terms that appear upon integrating by parts.

### Vortex transport and stretching.

Let $`𝐪=\times 𝐯`$ be the vorticity of the unfiltered Eulerian fluid velocity. The curl of the Kelvin-filtered Navier-Stokes equation (11) gives the vortex transport and stretching equation,

$$\frac{𝐪}{t}+\widehat{𝐯}𝐪𝐪\widehat{𝐯}=\nu ^2𝐪+\times 𝐟.$$
(14)
We note that the coefficient $`\widehat{𝐯}`$ in the vortex stretching term is the gradient of the spatially filtered Eulerian velocity $`\widehat{𝐯}`$. Thus, Kelvin-filtering tempers the vortex stretching in the modified Navier-Stokes equation (11), while preserving the original form of the vortex dynamics. This tempered vorticity stretching is reminicent of Leray’s approach in regularizing the Navier-Stokes equations. Except for the term $`(\widehat{𝐯})^\mathrm{T}𝐯`$, the Kelvin-filtered Navier-Stokes equation (11) is otherwise quite similar to Leray’s regularization of the Navier-Stokes equations proposed in 1934 . Extension of the Leray regularization to satisfy the Kelvin circulation theorem was cited as an outstanding problem in Gallavotti’s review . Indeed, for the case that $`𝒪`$ is the Helmholtz operator, full advantage of this regularized vorticity stretching effect was taken in  to prove the global existence and uniqueness of strong solutions for the three dimensional NS$`\alpha `$ model. Paper  also proves the estimate (1) for the fractal dimension of the global attractor for the NS$`\alpha `$ model.

The Navier-Stokes-alpha model emerges from the Kelvin-filtered Navier-Stokes equation (11) when one chooses the operator $`𝒪`$ to be the Helmholtz operator, thereby introducing a constant $`\alpha `$ that has dimensions of length,

$$𝒪=1\alpha ^2^2,\text{with}\alpha =const.$$
(15)
In this case, the filtered and unfiltered fluid velocities in Eq. (8) are related by

$$𝐯=(1\alpha ^2^2)\widehat{𝐯},$$
(16)
and the reconstructed pressure $`p`$ is related to the hydrodynamic pressure $`P`$ as

$$p=P\frac{1}{2}|\widehat{𝐯}|^2\frac{\alpha ^2}{2}|\widehat{𝐯}|^2.$$
(17)
The original derivation of the ideal Euler-alpha model (the $`\nu =0`$ case of the NS$`\alpha `$ model) was obtained by using the Euler-Poincaré approach in , . The physical interpretations of $`\widehat{𝐯}`$ and $`𝐯`$ as the Eulerian and Lagrangian mean fluid velocities were given in .

The kinetic energy norm corresponding to (9) for the NS$`\alpha `$ model is given by

$$E_\alpha =\left[\frac{1}{2}|\widehat{𝐯}|^2+\frac{\alpha ^2}{2}|\widehat{𝐯}|^2\right]d^3x.$$
(18)
This kinetic energy is the sum of a translational kinetic energy based on the spatially filtered Eulerian velocity $`\widehat{𝐯}`$, and a gradient-velocity kinetic energy, multiplied by $`\alpha ^2`$. By showing the global boundedness, in time, of the kinetic energy (18) one concludes that the coefficient $`\widehat{𝐯}`$ in the vortex stretching relation (14) is bounded in $`L^2`$ for the NS$`\alpha `$ model. Thus, the second term in the kinetic energy (18) imposes an energy penalty for creating small scales. The spatial integral of $`|\widehat{𝐯}|^2`$ in the second term has the same dimensions as filtered enstrophy (the spatial integral of $`|\times \widehat{𝐯}|^2`$, the squared filtered vorticity). For a domain with boundary, or when considering transformations to rotating frames, the spatial integral of $`\frac{1}{2}|\widehat{𝐯}|^2`$ in the second term in the kinetic energy norm (18) should be replaced by the integral of $`\mathrm{trace}(\widehat{𝐞}\widehat{𝐞})`$, where $`\widehat{𝐞}`$ is the strain rate tensor, $`\widehat{𝐞}=(1/2)(\widehat{𝐯}+\widehat{𝐯}^\mathrm{T})`$ of the filtered velocity in Euclidean coordinates . For the case we treat here, in Euclidean coordinates and in the absence of boundaries, all three of these norms are equivalent for an incompressible flow.

## 3 Comparison of the NS$`\alpha `$ model with LES equations

The large eddy simulation (LES) equations for an incompressible flow are conventionally written as

$$\frac{}{t}\overline{u}_i+\frac{}{x_j}\overline{u}_i\overline{u}_j=\frac{}{x_i}\overline{P}+\nu \frac{^2}{x_jx_j}\overline{u}_i\frac{}{x_j}\tau _{ij}$$
(19)
$$\text{with}\frac{\overline{u}_i}{x_i}=0.$$
(20)
The overbar denotes spatial filtering which, for a quantity $`f(𝐱)`$, is defined by the integral

$$\overline{f}(𝐱)=G(𝐱𝐱^{})f(𝐱^{})d^3x^{}L_Gf(𝐱),$$
(21)
where $`G`$ is a given kernel function, not necessarily the same as $`g`$ in Eq. (7). The quantities in Eq. (19), $`u_i,P,`$ and $`\nu `$ are the velocity components, pressure, and the kinematic viscosity, respectively, and the effects of subgrid-scale quantities on the resolved velocity are described by a subgrid-scale (SGS) stress tensor

$$\tau _{ij}=\overline{u_iu}_j\overline{u}_i\overline{u}_j,$$
(22)
which must be modeled in terms of the resolved quantities to close Eq. (19). Note that the unfiltered velocity $`u_i`$ is assumed to be a solution of the incompressible Navier-Stokes equation.

In contrast, the NS$`\alpha `$ equation (11) with definition (16) and (17) is written in index notation as

$$\frac{}{t}v_i+\widehat{v}_j\frac{v_i}{x_j}+v_k\frac{\widehat{v}_k}{x_i}=\frac{}{x_i}\left(P\frac{1}{2}\widehat{v}_k\widehat{v}_k\frac{\alpha ^2}{2}\frac{\widehat{v}_m}{x_l}\frac{\widehat{v}_m}{x_l}\right)+\nu \frac{^2}{x_jx_j}v_i,$$
(23)
where $`\widehat{v}_i`$ is the filtered quantity related to the NS$`\alpha `$ velocity $`v_i`$ through the equation

$$v_i=\widehat{v}_i\alpha ^2\frac{^2}{x_jx_j}\widehat{v}_i=\widehat{v}_i\alpha ^2^2\widehat{v}_i$$
(24)
and $`\alpha `$ is the filter length scale. Both fields $`v_i`$ and $`\widehat{v}_i`$ are incompressible since the ‘hat’ filter commutes with differentiation for a constant alpha.

The NS$`\alpha `$ equation (23) can be rearranged into a form similar to the LES equation (19)

$$\frac{}{t}v_i+\frac{}{x_j}v_iv_j=\frac{P}{x_i}+\nu \frac{^2}{x_jx_j}v_i\frac{}{x_j}\left((\widehat{v}_jv_j)v_i\alpha ^2\frac{\widehat{v}_k}{x_j}\frac{\widehat{v}_k}{x_i}\right).$$
(25)
### Two interpretations of the NS$`\alpha `$ velocity.

Both the NS$`\alpha `$ velocity $`v_i`$ and the filtered NS$`\alpha `$ velocity $`\widehat{v}_i`$ can be represented numerically using a lower resolution than required to represent the full NS velocity for the same problem. This is the property that they share with the filtered velocity in traditional LES. Nevertheless, it must be recognized from the derivation in the previous section that there is no LES filtering procedure that would produce the NS$`\alpha `$ equations from the full NS equations. Nevertheless because of the similar behavior both NS$`\alpha `$ velocities can be formally interpreted as the traditional LES velocity. If the NS$`\alpha `$ velocity $`v_i`$ in equation (25) is interpreted as the filtered Navier-Stokes velocity $`\overline{u}_i`$ in the LES equation (19) then the subgrid-stress tensor associated with the alpha model is

$$\tau _{ij}^\alpha =(\widehat{v}_jv_j)v_i\alpha ^2\frac{\widehat{v}_k}{x_j}\frac{\widehat{v}_k}{x_i}=\alpha ^2\left(\frac{\widehat{v}_k}{x_j}\frac{\widehat{v}_k}{x_i}\widehat{v}_i^2\widehat{v}_j\right)\alpha ^4(^2\widehat{v}_i)(^2\widehat{v}_j)$$
(26)
where the last equality is obtained using Eq. (24). This is a straightforward identification of terms.

Alternatively, one may propose to interpret the filtered NS$`\alpha `$ velocity $`\widehat{v}_i`$ as the filtered Navier-Stokes velocity $`\overline{u}_i`$. To obtain the corresponding SGS stress in this case we first filter the NS$`\alpha `$ equation (25) with the ‘hat’ filter. Since this filtering commutes with differentiation one finds

$$\frac{}{t}\widehat{v}_i+\frac{}{x_j}\widehat{v}_i\widehat{v}_j=\frac{\widehat{P}}{x_i}+\nu \frac{^2}{x_jx_j}\widehat{v}_i\frac{}{x_j}\left(\left(\widehat{\widehat{v}_i\widehat{v}_j}\widehat{v}_i\widehat{v}_j\right)\alpha ^2\left(\widehat{\frac{\widehat{v}_k}{x_j}\frac{\widehat{v}_k}{x_i}}+\widehat{\widehat{v}_j^2\widehat{v}_i}\right)\right),$$
(27)
which involves doubly filtered quantities. The corresponding SGS stress tensor is

$`\tau _{ij}^{f\alpha }`$ $`=`$ $`(\widehat{\widehat{v}_i\widehat{v}_j}\widehat{v}_i\widehat{v}_j)\alpha ^2\left(\widehat{{\displaystyle \frac{\widehat{v}_k}{x_j}}{\displaystyle \frac{\widehat{v}_k}{x_i}}}+\widehat{\widehat{v}_j^2\widehat{v}_i}\right)`$ (28)
$`=`$ $`(\widehat{\widehat{v}_i\widehat{v}_j}\widehat{\widehat{v}}_i\widehat{\widehat{v}}_j)\alpha ^2\widehat{{\displaystyle \frac{\widehat{v}_k}{x_j}}{\displaystyle \frac{\widehat{v}_k}{x_i}}}\alpha ^2\left(\widehat{\widehat{v}_j^2\widehat{v}_i}\widehat{\widehat{v}}_j^2\widehat{\widehat{v}}_i\right)`$
$`+\alpha ^2\widehat{\widehat{v}}_i^2\widehat{\widehat{v}}_j\alpha ^4(^2\widehat{\widehat{v}}_i)(^2\widehat{\widehat{v}}_j),`$
where the last equality is obtained using Eq. (24). We shall next consider the transformation properties of the NS$`\alpha `$ model. In section 5 we shall discuss how the SGS stress tensor $`\tau _{ij}^{f\alpha }`$ for the filtered NS$`\alpha `$ model is related to LES similarity models for NS.

## 4 Transformation properties

### Galilean invariance.

Because of the first term in Eq. (26), or equivalently the expression $`\widehat{v}_i^2\widehat{v}_j`$ in the second part of (26), the SGS tensor $`\tau _{ij}^\alpha `$ is neither symmetric nor Galilean invariant, in contrast to the ordinary SGS stress tensor (22). However, the SGS force which is the divergence of the SGS stress tensor is Galilean invariant for both representations of $`\tau _{ij}^\alpha `$. Therefore, the dynamics of the alpha model is Galilean invariant. Similarly, the SGS stress tensor $`\tau _{ij}^{f\alpha }`$ is not symmetric and lacks the Galilean invariance because of the term $`\widehat{\widehat{v}}_i^2\widehat{\widehat{v}}_j`$ but the Galilean invariance is restored in its divergence since $`\widehat{\widehat{v}}_j`$ is incompressible.

### Rotating frame.

It is of interest to consider a transformation of the SGS stress tensor to a frame of reference rotating with a uniform angular velocity $`\mathrm{\Omega }_i`$. The full SGS stress tensor (22) transforms in such a way that the SGS force, i.e. $`\tau _{ij}/x_j`$, in an inertial and in a rotating frame are the same. Horiuti  has recently analyzed several SGS models and showed that some of them do not satify this condition.

The velocity in a rotating frame, denoted by the asterisk, and the velocity in an inertial frame are related as follows

$$u_i=u_i^{}+ϵ_{imn}\mathrm{\Omega }_mx_n^{}.$$
(30)
Transforming the substantial derivative to a rotating frame and using (30) gives

$$\frac{}{t}u_i+\frac{}{x_j}u_iu_j=\frac{}{t}u_i^{}+\frac{}{x_j^{}}u_i^{}u_j^{}+2ϵ_{imn}\mathrm{\Omega }_mu_n^{}+ϵ_{imn}ϵ_{nkl}\mathrm{\Omega }_m\mathrm{\Omega }_kx_l^{},$$
(31)
where the last two terms are the Coriolis force and the centrifugal force, respectively. Transforming the NS$`\alpha `$ equation (25) to a rotating frame gives

$$\frac{}{t}v_i^{}+\frac{}{x_j^{}}v_i^{}v_j^{}=2ϵ_{imn}\mathrm{\Omega }_mv_n^{}ϵ_{imn}ϵ_{nkl}\mathrm{\Omega }_m\mathrm{\Omega }_kx_l^{}\frac{P}{x_i^{}}+\nu \frac{^2}{x_j^{}x_j^{}}v_i^{}\frac{}{x_j^{}}(\tau _{ij}^\alpha +\tau _{ij}^c),$$
(32)
where $`\tau _{ij}^\alpha `$ is the expression (26) for ‘star’ quantities and the additional correction term appears

$$\tau _{ij}^c=\alpha ^2\left(ϵ_{imn}\mathrm{\Omega }_mx_n^{}^2\widehat{v}_j^{}ϵ_{kmj}\mathrm{\Omega }_m\frac{\widehat{v}_k^{}}{x_i^{}}ϵ_{kli}\mathrm{\Omega }_l\frac{\widehat{v}_k^{}}{x_j^{}}\right).$$
(33)
The divergence of this correction stress tensor does not vanish; in fact,

$$\frac{}{x_j^{}}\tau _{ij}^c=\alpha ^2\left(2ϵ_{imn}\mathrm{\Omega }_m^2\widehat{v}_n^{}ϵ_{kmj}\mathrm{\Omega }_m\frac{^2\widehat{v}_k^{}}{x_i^{}x_j^{}}\right).$$
(34)
However, the first term in (34) combines with the Coriolis force in Eq. (32) and the second term combines with the centrifugal force to yield the following form of the NS$`\alpha `$ equation in the rotating frame:

$$\frac{}{t}v_i^{}+\frac{}{x_j^{}}v_i^{}v_j^{}=2ϵ_{imn}\mathrm{\Omega }_m\widehat{v}_n^{}\frac{\mathrm{\Pi }}{x_i^{}}+\nu \frac{^2}{x_j^{}x_j^{}}v_i^{}+\frac{}{x_j^{}}\tau _{ij}^\alpha ,$$
(35)
where the total pressure $`\mathrm{\Pi }`$ is

$$\mathrm{\Pi }=P+\frac{1}{2}(\mathrm{\Omega }_mx_m^{})(\mathrm{\Omega }_lx_l^{})\frac{1}{2}(\mathrm{\Omega }_m\mathrm{\Omega }_m)(x_l^{}x_l^{})\alpha ^2ϵ_{mjk}\mathrm{\Omega }_m\frac{\widehat{v}_k^{}}{x_j^{}}.$$
(36)
### Remarks.

This slight pressure shift is of no importance for incompressible fluid dynamics. Had we replaced $`\frac{1}{2}|\widehat{𝐯}|^2`$ in the second term of the kinetic energy norm (18) by $`\mathrm{trace}(\widehat{𝐞}\widehat{𝐞})`$, as mentioned earlier, we would have removed this pressure shift at the cost of making the motion equation and stress tensor slightly more complicated. We shall decline this option here; however, it shall be required in the compressible case.

Note that in the transformed motion equation (35) the Coriolis force is computed using filtered velocity $`\widehat{𝐯}^{}`$. This is the physical velocity of the fluid parcels in the rotating frame.

One may also obtain the nondissipative part of the NS$`\alpha `$ motion equation written in a rotating frame as in (35) by first substituting $`\widehat{v}_i=\widehat{v}_i^{}+ϵ_{imn}\mathrm{\Omega }_mx_n^{}`$ into Hamilton’s principle in which the Lagrangian is given by the kinetic energy (18) and then taking variations in the Euler-Poincaré framework, as discussed in .

## 5 Relation to generalized similarity models

We shall discuss the relation of the NS$`\alpha `$ model to other SGS models in terms of the inversion of the filtering operation (21). In generalized similarity models for the NS equations the full NS velocity $`u_i`$ in the expression for the SGS stress tensor (22) is obtained by an inversion of the filtering operation (21), sometimes called deconvolution. The deconvolution is frequently encountered in image processing and its properties are discussed in many textbooks on this subject, e.g. . In general, the inversion of the continuous filter is a ill-conditioned and often a singular problem precluding a unique solution. This difficulty leads to a concept of an approximate or regularized deconvolution such as a pseudo-inverse through a singular value decomposition approach . In the context of subgrid-scale modeling the concept of deconvolution was discussed previously by several authors. Germano  introduced an exponential filter for which the filtering operation can be inverted exactly, i.e. the function $`f(x)`$ in (21) can be represented as a closed form expression in terms of the filtered function $`\overline{f}(x)`$. Using this representation one obtains a formal expression for the subgrid-scale stress tensor as a function of the filtered field. Germano  shows that for this filter the dominant term in the expression for the subgrid scale stress has the form of a nonlinear velocity gradient model

$$\tau _{ij}\frac{\overline{u}_i}{x_k}\frac{\overline{u}_j}{x_k}.$$
(37)
This form, considered first by Leonard , has been more recently investigated by Liu et al. , Borue and Orszag , Leonard , and Winckelmans et al. . In addition to the exact deconvolution procedure of Germano  and Leonard  other approaches have been proposed and utilized recently. Geurts  constructs an approximate inverse operator by requiring that polynomials up to a certain order are recovered exactly from their filtered counterparts. Stolz and Adams  use a formal power series expansion of the inverse operator and obtain the approximate inverse by truncating the expansion. Similar methods are used by Horiuti  leading to a multi-level filtered model. The deconvolution is also an integral part of the SGS estimation model of Domaradzki et al. . In the spectral space version of the estimation model an exact deconvolution is possible . In the physical space representation the approximate deconvolution is accomplished by solving a sequence of tridiagonal equations .

Here, following Stolz and Adams  and Horiuti  we shall consider a formal inverse of the linear operator $`L_G`$ in (21) as a power series expansion

$$L_G^1=(I(IL_G))^1=I+(IL_G)+(IL_G)^2+\mathrm{}$$
(38)
where $`I`$ is the identity operator. To first order in the expansion (38) the velocity $`u_i`$ is approximated as

$$u_i=L_G^1\overline{u}_i\overline{u}_i+(\overline{u}_i\overline{\overline{u}}_i).$$
(39)
The second term in the approximation (39) is obtained by filtering of the difference $`(u_i\overline{u}_i)`$ which can be written using the Taylor expansion as

$`u_i\overline{u}_i`$ $`=`$ $`{\displaystyle \left(u_i(𝐱)u_i(𝐱+𝐱^{})\right)G(𝐱^{})d^3x^{}}`$ (40)
$`=`$ $`{\displaystyle \left(\frac{u_i}{x_j}x_j^{}\frac{1}{2}\frac{^2u_i}{x_kx_m}x_k^{}x_m^{}+\mathrm{}\right)G(𝐱^{})d^3x^{}}{\displaystyle \frac{\mathrm{\Delta }^2}{24}}^2u_i(𝐱)`$ (41)
where the derivatives are evaluated at $`𝐱`$ and the above relation is valid for both the symmetric top hat filter or the Gaussian filter with the width $`\mathrm{\Delta }`$ in each Cartesian direction. Using this relation the approximate inversion (39) becomes

$$u_i\overline{u}_i\frac{\mathrm{\Delta }^2}{24}^2\overline{u}_i,$$
(42)
upon again invoking commutation of filtering and differentiation. By comparing (42) for the approximate inversion and (24) for the definition of $`\widehat{v}_i`$ in the NS$`\alpha `$ model, we recognize that $`v_i`$ and $`\widehat{v}_i`$ are related by the approximate inverse (39) with $`\alpha ^2=\mathrm{\Delta }^2/24`$. We shall pursue this relation farther by comparing the NS$`\alpha `$ model with LES generalized similarity models.

Employing the approximate inverse (42) in the expression for the SGS stress tensor (22) one finds

$`\tau _{ij}^{sim}`$ $`=`$ $`\overline{u_iu}_j\overline{u}_i\overline{u}_j(\overline{\overline{u}_i\overline{u}_j}\overline{\overline{u}}_i\overline{\overline{u}}_j)`$
$`\alpha ^2\left[(\overline{\overline{u}_j^2\overline{u}_i}\overline{\overline{u}}_j^2\overline{\overline{u}}_i)+(\overline{\overline{u}_i^2\overline{u}_j}\overline{\overline{u}}_i^2\overline{\overline{u}}_j)\right]`$
$`+\alpha ^4(\overline{^2\overline{u}_i^2\overline{u}_j}^2\overline{\overline{u}}_i^2\overline{\overline{u}}_j).`$
The generalized similarity tensor $`\tau _{ij}^{sim}`$ in expression (5) is symmetric and Galilean invariant. On comparing the last equation with the expression for $`\tau _{ij}^{f\alpha }`$ in Eq. (3), we find as the most prominent difference the product derivative term in the filtered alpha-model equation

$$\tau _{ij}^{f\alpha }=\tau _{ij}^{\alpha sim}\alpha ^2(\widehat{\frac{\widehat{v}_k}{x_j}\frac{\widehat{v}_k}{x_i}}\widehat{\widehat{v}_i^2\widehat{v}_j})\alpha ^4\widehat{(^2\widehat{v}_i)(^2\widehat{v}_j})=\tau ^{\alpha sim}_{ij}+\widehat{\tau _{ij}^\alpha },$$
(44)
where $`\tau _{ij}^{\alpha sim}`$ is expression (5) written in terms of $`v_i`$ and the ‘hat’ filter, instead of $`u_i`$ and the ‘bar’ filter. Also $`\widehat{\tau }_{ij}^\alpha `$ is the ‘hat’-filtered version of $`\tau _{ij}^\alpha `$ in Eq. (26). According to Eq. (44) the SGS stress tensor for the filtered NS$`\alpha `$ equation (27) is the sum of the similarity stress tensor for NS and the filtered SGS stress tensor (26) for the original NS$`\alpha `$ equation (25). The last equality is the Germano identity  and follows from the the definitions of $`\tau _{ij}^{\alpha sim}`$, $`\tau _{ij}^\alpha `$, and $`\tau _{ij}^{f\alpha }`$.

We note that the important new term appearing in all SGS stress tensor expressions for the alpha-model, namely the product of derivatives, is different from the similar product of derivatives in the nonlinear models (37).

The NS-$`\alpha `$ model is reminiscent, but not equivalent to the generalized similarity models. This is encouraging because it may avoid difficulties encountered by the pure similarity models. One of the most serious problems encountered by the pure similarity models is insufficient SGS dissipation caused by neglecting in the modeling scales which cannot be resolved on a given LES grid. This fact has been long recognized, e.g. by Zhou et al. , and more recently by Domaradzki and Loh . Indeed, consider the most favorable situation when the filtering operation can be inverted exactly, i.e. the function $`f(x)`$ in (21) can be represented as a closed form expression in terms of the filtered function $`\overline{f}(x)`$. In such a case one immediately obtains an expression for the subgrid-scale stress tensor as a function of the filtered field from the definition (22). However, the apparent simplicity of such an exact approach to subgrid-scale modeling is deceptive. This is because in the procedure $`f(x)`$ and $`\overline{f}(x)`$ have the same spectral support since the same LES mesh is used to represent both functions. Therefore, no modeling of the subgrid scales is actually performed in the deconvolution procedure. As a result the similarity SGS stress tensor is unable to account for the dissipative effects of the nonlinear interactions between the resolved scales and the unknown, subgrid scales which cannot be represented on a given LES mesh. This implies that a good SGS model should either explicitly model subgrid scales or possess extra terms in addition to the similarity part. The former approach is used in the SGS estimation modeling  and the additional terms are encountered in the mixed models  as well as in the NS-$`\alpha `$ model (44).

### Section 5 summary.

We have identified the ‘hat’-filtered NS$`\alpha `$ equations with a class of generalized LES similarity models in which the SGS stress tensor possesses additional terms involving products of derivatives, as in $`\widehat{\tau }_{ij}^\alpha `$. Applying the approximate-inverse methods of deconvolution to this class of models restores the Kelvin circulation theorem to the deconvolved solution. The corresponding dynamics of the deconvolved solution is the NS$`\alpha `$ model, whose analytical and numerical properties are summarized in .

## 6 Spectral space interpretation

The nonlinear term in the NS$`\alpha `$ equation (25) transforms to spectral (Fourier) space as

$$F.T.\left[\frac{}{x_j}v_l(𝐱)\widehat{v}_j(𝐱)\right]=ik_jv_l(𝐩)\widehat{v}_j(𝐤𝐩)d^3p=ik_jv_l(𝐤𝐩)\widehat{v}_j(𝐩)d^3p.$$
(45)
Here the spectral representation is indicated by the dependence on wavenumbers. The incompressibility condition requires projection of (45) onto a plane normal to the vector $`𝐤`$ which is accomplished by applying the tensor

$$P_{lm}(𝐤)=(\delta _{lm}k_lk_m/k^2).$$
(46)
In spectral space one has

$$v_j(𝐤)=\widehat{v}_j(𝐤)+\alpha ^2k^2\widehat{v}_j(𝐤),$$
(47)
so the projection of the nonlinear term can be written as

$$\frac{1}{2}i\left(k_jP_{lm}d^3p\frac{v_j(𝐩)}{1+\alpha ^2p^2}v_m(𝐤𝐩)+k_mP_{lj}d^3pv_j(𝐩)\frac{v_m(𝐤𝐩)}{1+\alpha ^2|𝐤𝐩|^2}\right).$$
(48)
This expression allows us to interpret the alpha-model in terms of triad interactions. Qualitatively, the factor $`1/(1+\alpha ^2k^2)`$ attenuates amplitudes of modes with wavenumbers $`k>1/\alpha `$ and leaves modes with $`k<1/\alpha `$ largely unaffected.

Consider first a mode with $`k<1/\alpha `$. Its triadic interactions with all other modes $`𝐩`$ and $`𝐪=𝐤𝐩`$ for which $`p,q<1/\alpha `$ are only slightly affected by the attenuation factors and so are essentially the same as for the Navier-Stokes dynamics. If one of the wavenumbers $`p,q`$ is greater than $`1/\alpha `$ its effect on the mode $`𝐤`$ is weakened by one of the attenuation factors in (48). In particular, the highly nonlocal interactions of the eddy viscosity type, $`k<1/\alpha `$ and $`p,q1/\alpha `$, are removed from the dynamics.

Consider now a mode with $`k>1/\alpha `$. In this case all local triads with $`p,qO(k)`$ are attenuated and the nonlocal triads with $`p,q1/\alpha `$ are again removed from the dynamics. However, not all nonlocal triads affecting the mode $`𝐤`$ are neglected. Indeed, for $`k>1/\alpha `$ it is possible to have triads with $`p<1/\alpha <k`$ and $`qO(k)`$ and then at least one of the terms in (48) is not directly affected by the attenuation factors. The NS$`\alpha `$ model thus retains nonlocal sweeping effects of large scales (either $`p`$, or $`q`$ $`<1/\alpha <k`$) acting on small scales $`k>1/\alpha `$.

### Section 6 summary.

The NS$`\alpha `$ model is equivalent to the Navier-Stokes dynamics for all triads with $`k,p,q<1/\alpha `$, it consistently neglects the nonlocal interactions of the eddy viscosity type for modes $`p,q>1/\alpha `$, and retains nonlocal sweeping effects of large eddies on small eddies ($`p`$ or $`q`$ $`<1/\alpha <k`$ ).

## 7 Acknowledgements

It is our pleasure to acknowledge instructive and stimulating discussions on many of these matters with P. Constantin, G. Eyink, C. Foias, B. J. Geurts, N. Sandham, K. Sreenivasan, E. Titi, C. Vassilicos and V. Yakhot. We are grateful for hospitality at the UC Santa Barbara Institute for Theoretical Physics where this work was initiated during their Hydrodynamic Turbulence program in spring 2000. ITP is supported by NSF. One of the authors (JAD) was also partially supported by the NSF Grant No. CTS-9704728.


format: text

text:
# References

Do neutrino oscillations allow an extra phenomenological parameter?

I.S. Tsukerman

State Research Center

“Institute for Theoretical and Experimental Physics”

Moscow, 117218, Russia

E-mail:zuckerma@heron.itep.ru

Submitted to JETP Letters

## Abstract

The quantity $`\xi `$ introduced recently in the phenomenological description of neutrino oscillations is in fact not a free parameter, but a fixed number.

The literature on phenomenology of neutrino oscillations is vast (see, e.g.,  \-  and references therein). In a recent paper  Giunti and Kim in the case of two-flavour mixing have introduced a new phenomenological parameter $`\xi `$. According to , $`\xi =0`$ corresponds to the so-called equal momentum assumption , while $`\xi =1`$ corresponds to equal energy assumption . Authors of  emphasize that $`\xi `$ disappears from final expressions for the neutrino oscillation probability.

The aim of this note is to indicate that parameter $`\xi `$ is fixed by energy-momentum conservation in the process which is responsible for neutrino emission, as explicitly assumed in ref. .

Following ref.  we will consider the decay $`\pi \mu \nu `$ in the framework of two-flavour toy model. The parameter $`\xi `$ is defined in  for the pion rest-frame by considering the auxiliary case of absolutely massless neutrinos and denoting the energy of such neutrinos as $`E`$,

$$\xi ={}_{}{}^{1}/_{2}^{}(1+m_\mu ^2/m_\pi ^2),$$
(1)
where $`m_\mu `$ and $`m_\pi `$ are the masses of the muon and the pion. Then for massive (but light!) neutrinos they get:

$$E_{1,2}=E+(1\xi )m_{1,2}^2/2E,$$
(2)
$$p_{1,2}=E\xi m_{1,2}^2/2E.$$
(3)
Here $`E_{1,2}`$, $`p_{1,2}`$ and $`m_{1,2}`$ are the energies, momenta and masses of the neutrinos, respectively. Wherefrom the above statement about $`\xi =0,1`$ follows:

$$E_1=E_2\mathrm{for}\xi =1\mathrm{and}p_1=p_2\mathrm{for}\xi =0.$$
(4)
Thus, the equal energy and equal momentum assumptions in the form $`\mathrm{\Delta }EE_1E_2=0`$ and $`\mathrm{\Delta }pp_1p_2=0`$, respectively, are treated by authors of ref.  as particular cases of the general kinematical relations (2) and (3):

$$\mathrm{\Delta }E=(1\xi )\mathrm{\Delta }m^2/2E=0\mathrm{for}\xi =1,$$
(5)
$$\mathrm{\Delta }p=\xi \mathrm{\Delta }m^2/2E=0\mathrm{for}\xi =0.$$
(6)
Unfortunately, both the treatment and the relations (4) - (6) are erroneous.

On one hand, the quantity $`\xi `$ is not a free parameter for a certain decay process. Indeed, it follows from (1) that $`\xi `$ has a fixed value ($`0.8`$) for the decay under consideration. On the other hand, it is evident from definitions of $`E`$ and $`\xi `$ that

$$E=m_\pi (1\xi ).$$
(7)
The parameter $`\xi `$ determines sharing of the decay energy. As seen from (7), the values $`\xi =0`$ and $`\xi =1`$ are senseless ones because they refer accordingly to the limiting cases of $`E_\mu =0`$ and $`E=0`$. Therefore one cannot assume that $`\xi `$ can be equal to 1 or 0. Instead of that, the solution of the equalities (5) and (6) is the vanishing $`\mathrm{\Delta }m^2`$, that is absence of the oscillations.

Acknowledgements

The author is grateful to L.B. Okun for his friendly support. This work was supported by RFBR grant 00-15-96562.


text:
# Direct measurement of single fluxoid dynamics in superconducting rings

## Abstract

We have measured the dynamics of individual magnetic fluxoids entering and leaving photolithographically patterned thin film rings of underdoped high-temperature superconductor Bi2Sr2CaCu2O8+δ, using a variable sample temperature scanning SQUID microscope. These measurements can be understood within a phenomenological model in which the fluxoid number changes by thermal activation of a Pearl vortex in the ring wall. We place upper limits on the “vison” binding energy in these samples from these measurements.

Although there is a vast literature on vortex dynamics in superconductors , with a few notable exceptions  this work has involved indirect measurements of collective motions of vortices through, to cite a few recent examples, transport , voltage noise , magnetization loops , persistent currents in rings , or microwave impedance . In this Letter, we present direct measurements of the dynamics of individual magnetic fluxoids leaving and entering superconducting rings with a well defined geometry. These measurements represent a new tool for studying vortex dynamics.

Our measurements were made on 300 nm thick films of the high-temperature superconductor Bi2Sr2CaCu2O8+δ(BSCCO), epitaxially grown on (100) SrTiO3 substrates using magnetron sputtering. The oxygen concentration in these films was varied by annealing in oxygen or argon at 400-450 oC. The films were photolithographically patterned into circular rings using ion etching. The rings had outside diameters of 40, 60, and 80 $`\mu `$m, with inside diameters half the outside diameters. The film for the current measurements had a broad resistive transition (90% of the extrapolated normal state resistance at T=79 K, 10% at T=46 K) with a zero-resistance $`T_c`$ of 36 K before patterning. After patterning the rings had $`T_c`$’s from 25 K to 33 K, as indicated by inductive measurements. The rings were magnetically imaged using a variable sample temperature scanning Superconducting Quantum Interference Device (SQUID) microscope, which scans a sample relative to a SQUID with a small, well shielded, integrated pickup loop (a square loop 17.8 $`\mu `$m on a side for these measurements), the sample temperature being varied while the low-$`T_c`$ SQUID remains superconducting. Figure 1(a) shows a scanning SQUID microscope image of some of these rings, cooled in a magnetic field sufficient to trap one vortex in each ring. The dots in Figure 1(b) are a cross-section through the center of one ring. The curve in Figure 1(b) is modeled as follows:

Consider a thin film ring of thickness $`d\lambda _L`$ (the London penetration depth) with radii $`a<b`$ in the plane $`z=0`$. The London equations for the film interior read

$$𝐣=\frac{c\mathrm{\Phi }_0}{8\pi ^2\lambda _L^2}\left(\theta +\frac{2\pi }{\mathrm{\Phi }_0}𝐀\right),$$
(1)
where $`𝐣`$ is the supercurrent density, $`\mathrm{\Phi }_0=hc/2e`$ is the superconducting flux quantum, $`\theta `$ is the order parameter phase, and $`𝐀`$ is the vector potential. Since the current in the ring must be single valued, $`\theta =N\phi `$, where $`\phi `$ is the azimuth and the integer $`N`$ is the winding number (vorticity) of the state. Integrating $`𝐣`$ over the film thickness $`d`$, we obtain:

$$g_\phi g(r)=\frac{c\mathrm{\Phi }_0}{4\pi ^2\mathrm{\Lambda }}\left(\frac{N}{r}\frac{2\pi }{\mathrm{\Phi }_0}A_\phi \right),$$
(2)
where $`g(r)`$ is the sheet current density directed along the azimuth $`\phi `$, and $`\mathrm{\Lambda }=2\lambda _L^2/d`$ is Pearl’s film penetration depth . The vector potential $`A_\phi `$ can be written as

$$A_\phi (r)=_a^b𝑑\rho g(\rho )a_\phi (\rho ;r,0)+\frac{r}{2}H,$$
(3)
where the last term represents a uniform applied field $`H`$ in the $`z`$ direction and $`a_\phi (\rho ;r,z)`$ is the vector potential of the field created by a circular unit current of a radius $`\rho `$:

$`a_\phi (\rho ;r,z)`$ $`=`$ $`{\displaystyle \frac{4}{ck}}\sqrt{{\displaystyle \frac{\rho }{r}}}\left[\left(1{\displaystyle \frac{k^2}{2}}\right)𝐊(k)𝐄(k)\right],`$ (4)
$`k^2`$ $`=`$ $`{\displaystyle \frac{4\rho r}{(\rho +r)^2+z^2}}.`$ (5)
Here, $`𝐊(k)`$ and $`𝐄(k)`$ are the complete elliptic integrals in the notation of Ref. .

Substituting Eq. (3) and (5) into (2), we obtain an integral equation for $`g(r)`$:

$`{\displaystyle \frac{4\pi ^2\mathrm{\Lambda }}{c}}rg(r)+\pi r^2H\mathrm{\Phi }_0N`$ (6)
$`={\displaystyle \frac{4\pi }{c}}{\displaystyle _a^b}𝑑\rho g(\rho )\left[{\displaystyle \frac{\rho ^2+r^2}{\rho +r}}𝐊(k_0)(\rho +r)𝐄(k_0)\right],`$ (7)
where $`k_0^2=4\rho r/(\rho +r)^2`$. This equation is solved by iteration for a given integer $`N`$ and field $`H`$ to produce current distributions which we label as $`g_N(H,r)`$.

After $`g_N(H,r)`$ is found, the field outside the ring can be calculated using Eq. (5):

$`h_z(N;r,z)={\displaystyle \frac{2}{c}}{\displaystyle _a^b}{\displaystyle \frac{d\rho g_N(H,\rho )}{\sqrt{(\rho +r)^2+z^2}}}[𝐊(k)`$ (8)
$`+{\displaystyle \frac{\rho ^2r^2z^2}{(\rho r)^2+z^2}}𝐄(k)]+H.`$ (9)
The flux through the SQUID is obtained numerically by integrating Eq. (8) over the pickup loop area. The line in Fig. 1b is a two parameter fit of this integration of Eq. (8) to the data, resulting in $`z`$=3.5 $`\mu `$m, and $`\mathrm{\Lambda }`$=9 $`\mu `$m (corresponding to $`\lambda _L`$=1.1 $`\mu `$m).

The fluxoid number $`N`$ of a ring could be changed by varying an externally applied flux $`\mathrm{\Phi }_a=A_{eff}H`$, $`A_{eff}\pi (a^2+b^2)/2`$ , and monitored by positioning the SQUID pickup loop directly over it . We always observed single fluxoid switching events, as determined by the agreement (to within 10%) of the measured spacing in applied flux between vortex switching events with our calculations for $`\mathrm{\Delta }N=1`$, in the experiments reported here. Switching distributions $`P(\mathrm{\Phi }_{a,i})`$ were obtained by repeatedly sweeping the applied field, in analogy with experiments on Josephson junctions . The transition rates $`\nu `$ of the fluxoid states were determined from this data using

$$\nu (\mathrm{\Phi }_{a,m})=\frac{d\mathrm{\Phi }_a/dt}{\mathrm{\Delta }\mathrm{\Phi }_a}\mathrm{ln}\left\{\underset{j=1}{\overset{m}{}}P(\mathrm{\Phi }_{a,j})/\underset{i=1}{\overset{m1}{}}P(\mathrm{\Phi }_{a,i})\right\},$$
(10)
where $`m=1`$ labels the largest $`\mathrm{\Phi }_a`$ in a given switching histogram peak , and $`\mathrm{\Delta }\mathrm{\Phi }_a`$ is the flux interval between data points.

The dots in Figure 2 show the results from one such experiment. The lines are from a phenomenological model (a full theory will be published elsewhere): We take the energy of the ring in its initial or final state to be

$$E_r(N)=\frac{\mathrm{\Phi }_0^2}{2L}(N\varphi _a)^2,$$
(11)
where $`\varphi _a=\mathrm{\Phi }_a/\mathrm{\Phi }_0`$, and $`L`$ is the inductance of the ring. We take the ring maximum energy during the transition $`NN^{^{}}`$ to be

$$E_t(N,N^{^{}})=E_V\alpha (N+N^{^{}})^2/4+(E_r(N)+E_r(N^{^{}}))/2,$$
(12)
where the first two terms on the right hand side represent the energy required to nucleate a vortex in the ring wall. The maximum vortex energy in a straight thin film superconducting strip of width $`W\mathrm{\Lambda }`$ (carrying no transport supercurrent) is

$$E_V=\frac{\mathrm{\Phi }_0^2}{8\pi ^2\mathrm{\Lambda }}\mathrm{ln}\left\{\frac{2W}{\pi \xi }\right\}.$$
(13)
Since $`\lambda _L1/\sqrt{1t^4}`$ ($`t=T/T_c`$) , we take $`E_V=E_{V0}(1t^4)`$. $`\alpha `$ is a fitting parameter which describes the reduction of the vortex nucleation energy with increasing $`N`$. For thermal activation, the $`NN^{^{}}`$ transition rate is $`\nu =\nu _0e^{E_a(N,N^{^{}})/k_BT}`$, where $`\nu _0`$ is an attempt frequency and the activation energy $`E_a(N,N^{^{}})=E_t(N,N^{^{}})E_r(N)`$. For this ring (from fits to the data of Fig.’s 2, 3 and 4) we take $`E_{V0}`$=6$`\times 10^{13}`$ erg, $`\alpha `$=1.7$`\times 10^{16}`$ erg, and $`\mathrm{\Phi }_0^2/2L`$=4.6$`\times 10^{14}`$ erg. Taking $`\mathrm{\Lambda }=200\mu m`$ at T=31.5 K from fits of Eq. (8) to SSM data, $`\xi =3.2/\sqrt{1t}`$ nm, and $`W=20\mu m`$, we find $`E_{V0}=1.55\times 10^{12}`$ erg, a factor of 2.6 larger than the 6 $`\times 10^{13}`$ erg from our fits. It has been proposed that surface defects could reduce the barrier to entry of vortices in type-II superconductors . If we calculate an effective inductance for the ring as $`L^{}=\mathrm{\Phi }_0/I_s`$, where $`I_s=_a^b𝑑rg_N(H;r)`$, and use the solution of Eq. (6) for $`N=1`$, $`H=0`$, and $`\mathrm{\Lambda }=200\mu m`$, we find $`\mathrm{\Phi }_0^2/2L^{}`$=1.8$`\times 10^{14}`$ erg, a factor of 2.6 smaller than the 4.6$`\times 10^{14}`$ erg from our fit. For ring temperatures close to $`T_c`$, telegraph noise (Figure 3a) from thermally activated switching $`NN^{^{}}`$ occurs. This noise peaks when $`\mathrm{\Phi }_a=(N+1/2)\mathrm{\Phi }_0`$, $`N`$ an integer. Figure 3b shows the telegraph noise frequency (number of steps up per sec) as a function of applied field for the ring of Figure 2 at a temperature T=31.6K. The line is the prediction of the model described above, using the same fitting parameters, and writing the telegraph noise frequency as

$$\nu =2\nu _0\frac{e^{E_{in}/k_BT}e^{E_{out}/k_BT}}{e^{E_{in}/k_BT}+e^{E_{out}/k_BT}},$$
(14)
where $`E_{in}=E_a(N,N+1)`$ and $`E_{out}=E_a(N+1,N)`$. Finally, the symbols in Figure 4 show the telegraph noise frequency for the same ring as in Figures 2 and 3, as a function of temperature with $`\mathrm{\Phi }_a=\mathrm{\Phi }_0/2`$. The line is the prediction of our simple model, with the same fitting parameters as above.

The data shown in Figures 2, 3, and 4 are all consistent with an attempt frequency $`\nu _0`$ 3$`\times 10^8`$ sec-1. Glazman and Fogel , in a treatment of quantum tunneling of vortices, write $`\nu _0=\sqrt{\mathrm{\Phi }_0B_{c2}/4\pi \mathrm{\Lambda }m}`$, where the vortex mass $`m=\mathrm{}\eta /\mathrm{\Delta }`$, and the damping parameter $`\eta =\mathrm{\Phi }_0B_{c2}d/\rho _nc^2`$. Taking the second critical field $`B_{c2}`$ equal to 1T for $`T_c`$-T=1K, the normal state sheet resistivity $`\rho _n`$=1200 $`\mu `$Ohm-cm, and $`\mathrm{\Delta }`$ = 5k$`{}_{B}{}^{}T_{c}^{}`$ = 2.24$`\times 10^{14}`$ erg, we find $`\nu _0`$ = 9.7$`\times 10^8`$ sec-1, within a factor of three of our measurements. A second estimate, in a treatment of thermal activation of vortices , is $`\nu _0=6.96(D/a^2)\sqrt{E_V/k_BT}`$, where the diffusion constant $`D=k_BT/\eta `$. Using $`a`$=40$`\mu `$m, and the other parameters the same as before, this expression gives $`\nu _0=1.5\times 10^5`$sec-1. The disparity between these two estimates provides a measure of the uncertainties involved in calculating attempt frequencies.

Senthil and Fisher  and Sachdev  have proposed tests of the idea that the electron is fractionalized in the high-$`T_c`$ cuprate superconductors. The Senthil-Fisher proposal is to look for persistence of vorticity in underdoped cuprate cylinders as they are cycled through the superconducting transition temperature. This experiment has been performed by Bonn et al.  on single crystals of YBa2Cu3O6.3, and by us on the present ring samples, with no evidence to date for this persistent vorticity. The predicted effect depends on the existence of a gapped topological excitation, dubbed a “vison”, that is associated with a conventional vortex below the superconducting transition temperature, but which persists up to the pseudogap temperature. The present experiments provide another test of these ideas. It has been proposed that the vison should have a binding energy above the superconducting transition temperature of order $`E_{vison}k_BT^{}`$ . The pseudogap temperature $`T^{}`$ is estimated to be approximately 300 K  for BSCCO with a $`T_c`$ of 30 K. If we assume that the vortex thermal activation energy at $`T_c`$ is equal to the vison binding energy $`E_{vison}`$ (since the vortex core and supercurrent contributions to the vortex energy have gone to zero), then the telegraph noise frequency extrapolated to $`T_c`$ should be given by $`\nu (T=T_c)=\nu _0e^{E_{vison}/k_BT_c}`$, or $`E_{vison}=k_BT_c\mathrm{ln}(\nu _0/\nu (T=T_c))`$. The larger of our two estimates for the attempt frequency $`\nu _0`$ would lead to $`E_{vison}/k_B60K`$. The smaller estimate would indicate a smaller binding energy. A conservative estimate for an upper limit of the attempt frequency would be $`\nu _0c/b`$ = 4$`\times 10^{12}`$ sec-1, where $`c`$ is the speed of light, and $`b`$ is the outer radius of the ring, which would lead to an upper limit for $`E_{vison}/k_B`$ 300 K.

In conclusion, we have demonstrated a technique for measuring the dynamics of single vortices on a relatively short time scale, limited by the modulation frequency (100KHz) of our SQUID electronics (modulation frequencies 103 times faster have been demonstrated). These measurements were made possible by using cuprate superconductors that were highly underdoped, so that the penetration depths were longer than the ring wall lengths, and the vortex activation energies were comparable to the temperature, over an appreciable temperature range below $`T_c`$. Although these measurements were apparently made in a regime where the fluxoid transitions were mediated by thermally activated Pearl vortices, it may be possible to use similar techniques to study single vortex tunneling .

We would like to thank T. Senthil and M.P.A. Fisher for suggesting these experiments to us, and R. Koch, S. Woods, D. Bonn, W.A. Hardy, D.J. Scalapino, and F. Tafuri for useful conversations. Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This research was supported by the Director for Energy Research, Office of Basic Energy Sciences.


text:
# The energy gap of intermediate-valent SmB6 studied by point-contact spectroscopy

## I Introduction

Samarium hexaboride SmB6 is a homogeneous intermediate \- valence compound in which the electronic structure at low temperatures shows a narrow energy gap as well as spin gap, both originating from the hybridization between the narrow states formed by electrons of samarium 4$`f`$ \- shell and the wide conduction band formed by both boron $`p`$ \- states and Sm 6$`s`$ \- states . A review of this and similar materials can be found in Refs. .

Recent experiments on SmB6  have shown that in the low-energy excitation spectrum of this material several energy scales exist, and also several regimes of low-temperature electron kinetics. At least three different activation energies determine the behaviour of the conduction electrons. In the temperature range $`70\mathrm{K}>T>15`$K the properties of SmB6 are governed by the hybridization gap $`E_g1020`$meV. Between $`15`$K and $`5`$K, a narrow in - gap band separated from the bottom of the conduction band by a direct activation energy of $`E_d35`$meV has been observed. The properties of this narrow band seem to be influenced by the content of impurities and imperfections of the specific sample. Below about $`5`$K the electrical conductivity saturates, indicating a small conductivity channel within the $`E_d`$ in - gap states, where the Fermi level is pinned.

Various models have been proposed to explain the formation of the $`E_d`$ \- band and the origin of the residual conductivity , but so far no final conclusion could be obtained. While Ref.  favours hopping processes of electrons, Refs.  prefer a coherent metal-like state at low temperatures. In either case electrons are strongly localized ($`0.6`$nm localisation radius) at random impurities with a small concentration of $`N10^{23}\mathrm{m}^3`$ at low temperatures .

To search for anomalies in the quasi-particle density of states a number of experiments on small junctions with SmB6 has been performed. Frankowski and Wachter  brought a sharply etched molybdanum tip into direct contact with the surface of a cleaved SmB6 single crystal. The spectra of these low-resistance contacts showed a 4.6 meV wide anomaly (full width at half maximum) which was attributed to a gap in the density of states. The overall size of this anomaly was only about $`10\%`$ of the mean contact resistance. Güntherodt et al.  investigated Schottky-type tunnel contacts between SmB6 and GaAs. They found a huge but only 2.7 meV wide zero-bias anomaly. Such small gaps were obtained only when the SmB6 surface was sputter cleaned in situ. Without such a treatment the anomalies broadened to around 10 meV. Planar tunnel junctions with lead counter electrodes were investigated by Batlogg et al.  and by Amsler et al. . The latter experiments showed 14 meV wide spectral anomalies, roughly coinciding with the $`E_g`$ band gap, and a $`70\%`$ residual quasiparticle density of states at the Fermi level. No consistent picture could evolve from all these experiments, because the width as well as the size of the main zero-bias anomaly varied a lot. To obtain more and, hopefully, reliable information on the different energy scales involved we investigated direct junctions between two bulk pieces of SmB6 as function of contact size.

## II Experimental

Our SmB6 samples were cut from one batch which was grown by the zone-floating method as described in Ref. . About 800 ppm magnetic impurities, mostly magnetic lanthanide elements close to Sm in the periodic table, were found by magnetic susceptibility measurements . A similar amount of non-magnetic impurities, mainly lanthanum, was detected by induction-coupled plasma spectroscopy. sufficient size.

Two different type of point contacts were prepared. First, two mechanically polished pieces of SmB6 were brought into contact in a spear-anvil type setup and measured at $`4.2`$K in liquid helium. Second, bulk pieces of SmB6 were broken at a predefined notch in the ultrahigh vacuum region of a cold 3He-4He dilution refrigerator and measured mainly at 0.1 K. Reference  describes our break-junction apparatus in detail. With mechanically-controllable break junctions we avoid the oxidation of the interfaces. It also offers exellent mechanical stability of the contacts so that also very small junctions can be investigated. In both cases the $`dI/dU(U)`$ spectra were recorded in the standard four - terminal mode with current biasing.

## III Electrical conductivity and activation energies

Figure 1 shows the temperature dependence of the bulk conductivity $`\sigma `$($`T`$) of the single crystalline SmB6 sample. $`\sigma `$$`(T`$) decreases in a rather complicated way with decreasing $`T`$, saturating at the lowest temperatures. Describing the conductivity by thermally activated mechanism, $`\sigma (T)\mathrm{exp}(W/k_\mathrm{B}T)`$, an activation energy

$$W=d(\mathrm{ln}\sigma (T)/d(1/k_\mathrm{B}T)$$
(1)
can be defined (Fig. 2). It shows a peak of about 5.6 meV in the temperature range between 70 K and 15 K. This value can be attributed to $`E_g^{}/2`$ when the Fermi level sits just in the center of the hybridization band gap of width $`E_g^{}11.2`$meV. Between 15 K and 5 K the second pronounced peak apppears. Its value of $`E_d^{}3.7`$meV corresponds to the formation of the narrow in - gap band within the hybridization gap located $`E_d^{}`$ below the conduction band edge. These two activation energies do not clearly level off. Therefore the index ’\*’ is used to distinguish the activation-derived values from the true energy gaps still to be determined.

With further decreasing temperature, the activation energy $`W`$ falls to very low values. The anomaly at about $`12`$K with $`E_a0.2`$meV originates probably from an additional narrow band formed inside the energy gap (close to the $`E_d`$ \- band) due to the relatively large content of impurities of the sample. Below about 1 K, however, the activation energy is lower than the available thermal energy $`k_\mathrm{B}T`$, indicating a metallic - like conduction mechanism at the lowest temperatures. This would agree with recent $`ac`$ \- conductivity  and specific heat measurements of SmB6 , which are interpreted as showing a transition into a coherent heavy-fermion like state below about 5 K.

## IV Point-contact spectroscopy

We have investigated more than 100 SmB6 \- SmB6 point contacts with zero-bias conductance ranging from about 0.01 $`\mu `$S to about 1 mS, most of them using the spear-anvil technique. These spear-anvil type junctions had a high zero - bias conductance and usually symmetric $`dI/dU(U)`$ spectra, similar to those obtained earlier on SmB6 \- Mo point contacts . But few of those junctions were asymmetric like that in Figure 3. We found anomalies, typically small kinks in the $`dI/dU`$ spectra indicating a change of slope, at $`U4`$mV and $`U12`$mV, respectively.

Figures 4, 5, and 6 show typical $`dI/dU`$ spectra of SmB6 break junctions at $`T=0.1`$K. Like the spear-anvil type contacts they had anomalies in $`dI/dU`$ and in $`d^2I/dU^2`$ at about 4 mV and 12 meV. Very often additional anomalies were observed close to 0.2 mV, 1.8 mV, and 22 meV. Several junctions also had anomalies at about 80 mV.

Spectra with high zero-bias conductance $`G_0=dI/dU(T=0.1K,U=0)`$ shown in Fig. 4 were symmetric with respect to the applied voltage. Figs. 5 and 6 show that when the zero-bias conductance $`G_0`$ falls below about $`30\mu `$S, the differential conductance becomes more and more asymmetric. This tendency towards asymmetry could result from an inhomogeneous distribution of impurities and imperfections in the sample which can lead to different local carrier concentration, and which become the more pronounced the smaller the contact interface is. The different possible crystallographic orientations of the SmB6 electrodes forming the contact could also contribute to the asymmetry.

The spectra of junctions with low conductance (Fig. 6) appear to be broader than those with higher zero-bias conductance (Fig. 5). However, according to Fig. 7 the position of the anomalies is not affected by any variation of $`G_0`$. This points then to different weighing factors of the anomalies at small and at large zero-bias conductance, respectively, possibly reflecting the two different regimes of charge transport discussed below. In Fig. 7 we have labeled the positions according to the activation energies $`E_a`$, $`E_d^{}`$, and $`E_g^{}`$. Anomalies around 1.8 mV as well as around 22 meV could be related with $`E_g^{}/2E_d^{}`$ and $`2E_g^{}`$, respectively. And one could speculate whether the 80 mV anomalies (not shown in the figure) correspond to the excitation of electrons into higher energy levels of the Sm - ions, considered for example to interprete magnetic properties of SmB6 .

The properties of the junctions depend only weakly on temperature (Fig. 8) and on magnetic field (Fig. 9). In both cases, the shape of the spectra does not change much, only the conductance increases slightly when either the temperature is raised from 0.1 K to 1 K or the magnetic field from 0 Tesla to 8 Tesla (the field was always perpendicular to the direction of current flow). Thus also the magneto resistance is small and negative like in the bulk sample. Junctions with large zero-bias conductance $`G_01`$mS had typically $`G_0(B=0)/G_0(B)11.3\times 10^3B^2`$, where $`B`$ is the applied field in Tesla (Fig. 10).

## V Regimes of charge transport

To evaluate the observed $`dI/dU(U)`$ dependencies not only qualitatively, one should know the regime of charge transport across the contacts. This is absolutely necessary if one wants to attribute a spectral anomaly at a bias voltage $`U`$ to an anomaly in the density of states or to some scattering mechanism at an energy $`eU`$. We face here a quite general – but very often ignored – problem of direct contacts (that means contacts without well-defined tunneling barrier) between conductors with a short electronic mean free path. Naidyuk and Yanson  have recently reviewed and dicussed this topic.

As one typical example to illustrate the situation we refer to point-contact experiments with the Kondo semiconductor CeNiSn. It is ususally believed that direct junctions with this short mean free path compound are in the tunnelling regime, implying that the spectra measure the density of states . However, it was demonstrated recently that those point contacts are more likely metallic . And this, in turn, makes it impossible to extract the density of states.

We emphasize that the shape of the spectra itself does not tell anything about the regime of charge transport, whether it is tunnelling or diffusive transport, for example, unless one works with a well-known substance for which the spectra can be predicted reliably. Experiments with a superconducting lead counter electrode like that in Ref.  allow the quality of a tunnel junction to be verified. For such an experiment quantum tunneling can be confirmed or disproved. Obviously, this does not work for SmB6 in contact with an other piece of SmB6, since the exact properties of this compound have yet to be determined. Thus for our experiments the situation seems hopeless. But there exists what we believe are strong arguments to identify the relevant transport mechanisms. We use a procedure that has been developed for contacts with heavy-fermion compounds  as well as for junctions with doped germanium , both materials with a typically short electronic mean free path (or hopping length in the latter case).

## VI Large junctions – Maxwell regime

SmB6 has – due to hybridization effects – an extremely short electronic mean free path at low temperatures, which amounts to not more than several lattice constants. Therefore it is difficult to evaluate our experiments in terms of the classical metallic point contact spectroscopy. When the diameter $`d`$, the characteristic dimension of the point-contact, is large, it can certainly be estimated using Maxwell’s formula $`G\sigma d`$, where $`G=dI/dU`$ is the point-contact conductance at zero bias (a circular shape of the contact is a useful simplifying assumption). Two approaches are possible: One could set $`G`$ of Maxwell’s formula equal to the measured $`dI/dU(U=0)`$ and use the measured $`\sigma `$ at the same temperature. Thus, with $`\sigma (T=0.1\mathrm{K})0.05\mu `$S/cm from Fig. 1, a $`G_0=1\mu `$S junction has a $`d200`$nm diameter. However, it is better – and this has been clearly demonstrated for junctions with the heavy-fermion superconductors  – to use the change of conductance $`\delta G`$ and conductivity $`\delta \sigma `$, respectively, when the temperature is increased. Thus the contact diameter in the Maxwell regime actually becomes

$$d\delta \sigma /\delta G.$$
(2)
We have noted that the conductance increases by roughly $`5\%`$ when the temperature is increased from 0.1 to 1 K. At the same time, the bulk conductivity of our sample increases by a factor of 4 to about 0.2 S/cm. If we attribute this increase of conductance to the increase of bulk conductivity in Fig. 1 then, according to Eq. 2, the diameter amounts to about 3 nm only. (Applying this method on SmB6 is by far less reliable than for junctions with the heavy-fermion compounds. The above value should be regarded as an estimate for the order of magnitude of $`d`$.) Such a small contact is still formed by $`100`$ atoms which contribute to its mechanical stability. In fact, we have found a lack of stability, indicating the transition to ’one atom’ junctions, only when the conductance is further reduced by a factor of at least 100. This qualitative agreement supports our estimate of the contact size. Nevertheless, Maxwell’s formula is valid only for large junctions, and we have yet to determine its lower bound.

The large temperature-independent background needs to be explained: It is not due to impurities of the bulk sample. More likely the contact distorts the crystal lattice locally, and these additional defects enhance the local conductivity. A similar phenomenom – but of opposite direction – has been observed at junctions with heavy-fermion materials . Lattice defects there enhance the local resistivity and, consequently, the contact resistance. Using the absolute value of the contact resistance itself instead of its temperature-dependent part would then considerably underestimate the contact diameter, contrary to the SmB6 junctions.

Joule heating of our SmB6 \- SmB6 contacts could be expected at high bias voltages. A simple heating model  predicts the effective temperature of the contact

$$T_{eff}=\sqrt{T_b^2+U^2/4L}$$
(3)
at a bath temperature $`T_b`$. Here the Lorenz number $`L=\kappa /\sigma T`$, and $`\kappa `$ denotes the thermal conductivity of the material in the contact region. In the Drude - Sommerfeld theory of metals $`L`$ equals $`L_0=2.45\times 10^8`$V$`{}_{}{}^{2}/`$K2, implying a significant overheating of the contact. However, SmB6 has a very low electrical conductivity below about 5 K (Fig. 1) but a rather high thermal conductivity  due to heat transport by phonons. This yields a more than five orders of magnitude higher $`L`$ and provides a very effective cooling of the contact. Electrons crossing such a large junction take part in many scattering events, but basically they keep their original kinetic energy as they diffuse through the contact region. At a bath temperature of $`T_b=0.1`$K a voltage of $`U=30`$mV would heat up the contact to not more than $`T_{eff}0.4`$K. Therefore the applied voltage $`U`$ is too small to account for any of the observed anomalies of the $`dI/dU`$ spectra in terms of overheating, and the observed anomalies should mark the correct energies at which certain energy bands can be occupied or scattering mechanisms start to work.

The small magnetoresistance, less than $`10\%`$ at $`B=8`$Tesla for the junction in Fig. 9, also indicates no overheating. Otherwise a much larger effect had to be expected from the strong magnetoresistivity of the bulk sample above 1 K.

## VII Small junctions – Tunnelling

On further reducing the lateral size or diameter of the junctions one should expect a transition to ballistic transport if SmB6 was an ordinary metal. Obviously, this is not the case. Because of the low carrier concentration of $`N10^{23}\mathrm{m}^3`$ , the junctions will instead undergo a transition to tunnelling. Again there are two different approaches to estimate the contact diameter at which this transition will take place, depending on whether conduction is due to electrons in a coherent metal-like state or due to hopping.

First, if the electron system at low temperatures is described by a coherent state, the Fermi wave number amounts to about $`k_\mathrm{F}=(3\pi N)^{1/3}=10^8\mathrm{m}^1`$. In the contact region, discrete energy levels can exist which carry electrical current as long as the contact diameter $`d`$ is larger than $`4/k_\mathrm{F}40`$nm. If the diameter is smaller, electrons have to cross the junction as evanescent waves, that means they tunnel. Such a situation has been observed at direct junctions with antimon, also a compound with a low carrier density .

Second, if electron transport is due to thermally activated hopping between localized impurity states as proposed in Ref., for example, then junctions with SmB6 should be compared with those of a real semiconductor like germanium. Junctions with doped Ge seem to be in a tunnelling state although there is a direct contact with the bulk material without dielectric barrier . Transport across such contacts can be called ’tunnelling’ because the hopping length is so large that the charge carriers cross the junction interface by one single hop process. These hop processes across the junction are the same as that in the bulk material. We believe that our SmB6 junctions with small zero-bias conductance can be in such a tunneling-like state. A reasonable estimate for the impurity density is the carrier density, and the average distance betwen the impurities equals the hopping length. The transition to tunneling has then to be expected when the contacts become smaller than $`N^{1/3}22`$nm.

For both approaches, using Eq. 2 to estimate the contact diameter, the transition towards tunnelling should take place when the zero-bias conductance gets smaller than about $`10\mu `$S. And this is just what we are observing: a transition in the size of the zero-bias anomalies when $`G_0`$ is reduced. We fit a parabola to the spectra a high voltages as indicated in Fig. 4 to define as reference the conductance at zero bias $`G_H`$ if there was no gap in the density o states. (According to Ref.  one would have to increase the temperature of a specific junction to above 40 K to get an experimental value for $`G_H`$. This was not possible with our apparatus). The relative size of the zero-bias anomaly at low temperatures is then $`G_0/G_H`$. For contacts between SmB6 and a normal metal with constant density of states around the Fermi level, like that in Ref. , the differential conductance $`dI/dU(U)`$ is proportional to the energy dependence of the local electronic density of states multiplied by the transmission probability of the contact . Thus $`g_0/g_HG_0/G_H`$ when $`g_0`$, the low-temperature density of states at the Fermi level, is normalized to the density of states without gap $`g_H`$. Since for SmB6 \- SmB6 homocontacts the same quasiparticle density of states on both sides of the contact contribute to the tunnel conductance, $`g_0`$ has to be described by

$$g_0/g_H\sqrt{G_0/G_H}$$
(4)
Fig. 11 shows that large junctions (large $`G_0`$) also have large $`\sqrt{G_0/G_H}`$ because of their large background signal. This is inconsist with the quasiparticle density of states of around $`70\%`$ measured using planar tunnel junctions at $`T=10`$K . Moreover there seems to be a systematic decrease of the relative size of the anomaly as $`G_01`$mS, possibly indicating a transition from diffusive to thermal transport when progressively more scattering processes are required for an electron to pass the contact region. Its initial kinetic energy is then not conserved any more. The signal size of smaller junctions, however, agrees much better with the expected density of states if we take into account our much lower temperatures of $`T=0.1`$K, letting us attribute the zero-bias conductance to the finite density of quasi-particles in the gap at the Fermi level. Thus bulk transport turns to tunneling below about $`G_02\mu `$S, in excellent agreement with our above estimate.

## VIII The energy gaps and the quasiparticle density of states

Most of the anomalies, for the break junctions as well as for the spear-anvil type contacts, were found at positions that coincide well with the characteristic activation energies of SmB6 (Fig. 7). The important point here is that, whatever the interpretation of these energies: the positions do not depend on the zero-bias conductance, that means on the lateral contact size. It implies that there is no additional voltage drop in the contact area, which could pose a problem at very large junctions with $`G_01`$mS. Thus Fig. 7 demonstrates that at least up to $`G_01`$mS we can still derive the kinetic energy of the electrons from the applied bias voltage because of the suppressed local heating.

Anomalies at $`E_g^{}`$ as well as at $`E_d^{}`$ could be observed at all junctions. These dominant anomalies could indicate when the top of the valence band and the narrow impurity band of one electrode, respectively, face the bottom of the conduction band of the other electrode. At some junctions we could also resolve a small anomaly which corresponds to the lowest activation energy $`E_a`$ of our sample. It may represent an additional impurity band close to the $`E_d`$ band.

Not all junctions showed anomalies at around 0.2 meV, 2 meV, and 22 meV. An obvious reason for the sometime missing low-energy anomalies could be the degradation of the interface region due to the contact, increasing the local conductivity and thereby suppressing the low-energy processes. Two of the anomalies, the one at 1.8 meV and the other at 22 meV, have no counterpart as an activation energy from Fig. 2. And it is unclear whether they just accidentally coincide with $`E_g^{}/2E_d^{}`$ and $`2E_g^{}`$, respectively. To speculate whether the anomaly at twice $`E_g^{}`$ results from a double junction, that is a contact with two junctions in series, can be safely discarded because there are no corresponding anomalies at lower energies at $`2E_d^{}`$ or at $`2E_a`$. The two above anomalies could be artefacts, created by the stress in the contact region, and which is not present in the bulk sample. But there is no reason why in such a case the position of the other anomalies should not be affected in a similar way. This open question leads us to reconsider our way of attributing certain energies derived form the activation energy to the observed anomalies.

The necessary additional information can be obtained from the tunnelling spectra. As discussed above, a residual quasiparticle density of states at the Fermi level of about $`45\%`$, derived from the size of the zero-bias anomlies of small junctions with $`G_01\mu `$S and shown in Fig. 11, is considerably smaller than the data obtained by conventional tunnel junctions . However, this difference is not unreasonable because the latter experiments, compared to ours, were carried out at much higher temperatures when the $`E_d`$ anomaly has not fully developed. The rather large scattering of our data points in Fig. 11 probably results, in a minor part, from our method of extrapolating $`G_H`$, the zero-bias conductance without gap. The major part is due to the fact that the break junctions are a local probe. Each junction has slightly different local properties, for example due to stress induced by the contact. At planar tunnel junctions, on the other hand, the density of states is sampled and averaged over a much wider contact area. Of course, uniform stress and a preferred direction could affect the spectra of those planar junctions as well.

To estimate the average density of states at our SmB6 tunnel junctions, that means junctions with $`G_01\mu `$S, we normalized their spectra with respect to the voltage-dependent background as described in Fig. 4, averaged the normalized spectra, and then fitted the average spectrum by various functional dependencies of the density of states, assuming a constant transmission probability. A good fit was obtained by assuming a constant background of $`45\%`$ and two different energy-dependent parts $`g_{1,2}(E)`$ of the density of states as shown in Fig. 12.

The one part, $`g_1(E)`$, has a rather wide gap of about 21 meV (full width at half maximum), while the other part, $`g_2(E)`$, has a much narrower gap of about 4.5 meV. These two gaps may be attributed to the hybridization gap and the $`E_d`$ anomaly, respectively. But in this case the width of the larger anomaly differs a lot, almost a factor of two, from the activation energy data. On the other hand, if the gap in $`g_2(E)`$ was suppressed, leaving only the gap in $`g_1(E)`$, the anomaly of the density of states found in Ref.  for planar SmB6 \- Pb tunnel junctions at 10 K would be reproduced quite well (the anomaly described by Ref.  had a somewhat smaller width of about 14 meV, possibly because it was also asymmetric). This supports our interpretation of tunnelling due to the lateral confinement at our junctions. It seems as if the gap in $`g_2(E)`$ develops only at low temperatures $`T10`$K. And this has to be expected in view of the bulk conductivity and the activation energy.

The second derivative $`d^2I/dU^2`$ of the average spectrum as well as of the fit in Fig. 13 clearly show anomalies at around 1.6 mV, 5.5 mV, 11 mV, and 21 mV, respectively. Thus except for $`E_a=0.2`$eV all anomalies found in the spectra as displayed in Fig. 7 can also be recovered from the fit curve. We believe that the activation-derived gaps differ from the real gap values $`E_g=21`$meV and $`E_d=4.5`$meV because even at low temperatures the density of states remains finite inside the gaps. Therefore electrons can be thermally excited into states with energy $`EE_d`$ or $`EE_g/2`$, respectively, smearing out the otherwise discrete activation-energy levels.

## IX Conclusion

The differential conductance of SmB6 \- SmB6 junctions has been investigated at low temperatures as function of contact size. A very wide range of data on the conductance scale down to very small junctions was made possible because of the excellent mechanical stability of the break junctions when compared to the spear-anvil type contacts. Two different regimes of charge transport were distinguished. Large junctions are in the diffusive Maxwell regime, in which the conductance is dominated by the bulk conductivity. Local heating is negligible because of the large phonon heat conductivity. Therefore, the applied bias voltage can still be attributed to the kinetic energy of the charge carriers, enabling spectroscopy. Small junctions are in the tunnelling regime, although there is no dielectric barrier. Depending on the model used, electrons tunnel through the contact either by a single hopping event like in the bulk material or as evanescent waves because of a large Fermi wavenumber. Both models predict the transition to tunnelling when the zero-bias conductance becomes smaller than about $`10\mu `$S, in good agreement with our experiments. The spectra of these tunnel junctions indicate a finite $`45\%`$ residual density of states at the Fermi level. In the two different regimes of transport, at large as well as at small junctions, anomalies can be resolved at the same energies. However, only the spectra in the tunnelling regime reveal the correct gaps in the density of states of $`E_g=21`$meV and $`E_d=4.5`$meV that are responsible for the observed anomalies. The absolute values of these two energies fit well those derived from other experiments.

## X Acknowledgments

This work was supported by the Slovak Scientific Grant Agency VEGA, contract 1148 - 01 and 7022 - 20, and by the SFB 252 Darmstadt / Frankfurt / Mainz. US Steel Košice sponsored part of the liquid nitrogen.