text:
# Tunneling spectra for quasi-one-dimensional organic superconductors

## I Introduction

Almost $`20`$ years after the discovery of the organic superconductors , the problem of determination of their pairing state has not yet found a definite solution. The critical magnetic field $`H_{c2}`$ exceeds the Pauli paramagnetic limit and indicates that the pairing symmetry is triplet . The Knight shift does not change between the normal and superconducting states and it is a signature for triplet pairing state . However the absence of Hebel-Slichter peak and the power-low decay of $`1/T_1`$ below $`T_c`$  is an indication of the presence of nodes of the pair potential on the Fermi surface.

The scattering theory can be used to distinguish the symmetry of the pair potential . In $`d`$-wave superconductors the pair potential changes sign under a $`90^o`$-rotation. So under appropriate orientation of the $`a`$-axis within the $`ab`$ plane of $`d`$-wave superconductor the transmitted quasiparticles feel different sign of the pair potential. This results in the formation of bound states within the energy gap, which are detected as zero energy peaks (ZEP) in the conductance spectra .

The scattering theory has also been used for the determination of the pairing symmetry in (TMTSF)2X by Sengupta et al. . However their calculation is restricted to the presence or absence of ZEP at the surface, which can not distinguish the $`p`$ from the $`f`$-wave case. In a more realistic calculation Tanuma et al.  used the extended Hubbard model on a quasi-one-dimensional lattice at quarter-filling, to study the quasiparticle states near the surface of a quasi-one dimensional organic superconductor, where the pairing symmetry can actually be distinguished from the overall line shape of the surface density of states (SDOS) and the presence or absence of ZEP.

In this paper we extend the BTK formula to calculate the tunneling conductance in a normal-metal/insulator/quasi-one-dimensional organic superconductor, where the structure of the Fermi surface is taken into account. The pairing symmetry of the superconductor is assumed to be triplet $`p`$, singled $`d`$, $`d_{xy}`$, and triplet $`f`$-wave. In particular, we find that the ZEP appears in the electron tunneling along the $`a`$ or $`b`$ axis, when the transmitted quasiparticles experience different sign of the pair potential. Also the line shape of the spectra is sensitive to the presence or absence of nodes of the pair potential on the Fermi surface. The present calculation is useful for (TMTSF)2PF6, although this salt needs pressure to show superconductivity, but is not very suitable for (TMTSF)2ClO4, which is the only material exhibiting superconductivity in the Bechgaard salts at ambient pressure, but its unit cell is doubled due to anion ordering . Also tunneling spectroscopy on the quasi two-dimensional organic superconductor $`\kappa `$-(BEDT-TTF)2Cu(NCS)2 has been performed by using scanning tunneling microscopy . In this compound the conducting Cu(NCS)2 layers run parallel to the $`bc`$ plane, along the $`a`$-axis, and the in-plane tunneling data strongly suggest $`d`$-wave pairing state. These features can be used to distinguish the pairing symmetry in quasi-one-dimensional organic superconductor.

## II The model for the NS interface

The motion of quasiparticles in inhomogeneous superconductors is described by the four component Bogoliubov de Gennes (BdG) equations. The BdG equations are decoupled into two sets of (two component) equations, one for the spin up electron, spin down hole quasiparticle wave functions ($`u_{}(𝒓),v_{}(𝒓)`$) and the other for ($`u_{}(𝒓),v_{}(𝒓)`$). The BdG equation for spin index $`s(\overline{s})=()`$ or $`s(\overline{s})=()`$, read

$$\begin{array}{cc}(_e(𝒓)s\mu _BH)u_s(𝒓)+𝑑𝒓^{}\mathrm{\Delta }_{s\overline{s}}(𝒔,𝒙)v_{\overline{s}}(𝒓^{})\hfill & =E_su_s(𝒓)\hfill \\ 𝑑𝒓^{}\mathrm{\Delta }_{\overline{s}s}^{}(𝒔,𝒙)u_s(𝒓^{})(_e^{}(𝒓)s\mu _BH)v_{\overline{s}}(𝒓)\hfill & =E_sv_{\overline{s}}(𝒓)\hfill \end{array}.$$
(1)
The single-particle Hamiltonian is given by $`_e(𝒓)=\mathrm{}^2_𝒓^2/2m_e+V(𝒓)E_F`$, $`H`$ is the external magnetic field, $`E_s`$ is the energy measured from the Fermi energy $`E_F`$. For a given spin projection $`s`$ the magnetic field shifts the energy by $`s\mu _BH`$. $`\mathrm{\Delta }_{s\overline{s}}(𝒔,𝒙)`$ is the pair potential, after a transformation from the position coordinates $`𝒓,𝒓^{}`$ to the center of mass coordinate $`𝒙=(𝒓+𝒓^{})/2`$ and the relative vector $`𝒔=𝒓𝒓^{}`$. After Fourier transformation the pair potential depends on the relative wave vector $`𝒌`$ and $`𝒙`$. In the weak coupling limit $`𝒌`$ is fixed on the Fermi surface ($`|𝒌|=k_F`$), and only its direction is variable. After applying the quasi-classical approximation, i.e.,

$$\left(\begin{array}{cc}\overline{u}_s(𝒓)\hfill & \\ \overline{v}_{\overline{s}}(𝒓)\hfill & \end{array}\right)=e^{i𝒌𝒓}\left(\begin{array}{cc}u_s(𝒓)\hfill & \\ v_{\overline{s}}(𝒓)\hfill & \end{array}\right),$$
(2)
so that the fast oscillating part, of the wave function is divided out, the BdG equations are reduced to the Andreev equations

$$\begin{array}{cc}E_s\overline{u}_s(𝒓)\hfill & =iv_F𝒌\mathbf{}\overline{u}_s(𝒓)+\mathrm{\Delta }_{s\overline{s}}(𝒌,𝒓)\overline{v}_{\overline{s}}(𝒓)\hfill \\ E_s\overline{v}_{\overline{s}}(𝒓)\hfill & =iv_F𝒌\mathbf{}\overline{v}_{\overline{s}}(𝒓)+\mathrm{\Delta }_{\overline{s}s}^{}(𝒌,𝒓)\overline{u}_s(𝒓)\hfill \end{array},$$
(3)
where the quantities $`\overline{u}_s(𝒓)`$ and $`\overline{v}_{\overline{s}}(𝒓)`$ are electron-like and hole-like quasiparticles with spin index $`s`$, and $`\overline{s}`$ respectively, and $`v_F`$ is the Fermi velocity. We consider the normal-metal/insulator/superconductor junction shown in Fig. 1. The electron momentum parallel to the interface $`k_{}`$ is conserved. For the interface that is normal to the $`a`$-axis, the insulator is modelled by a delta function, located at $`x=0`$, of the form $`V\delta (x)`$. The temperature is fixed to $`0`$ K. We take the pair potential as a step function, i.e., $`\mathrm{\Delta }_{s\overline{s}}(k_{},𝒓)=\mathrm{\Theta }(x)\mathrm{\Delta }_{s\overline{s}}(k_{})`$. For the geometry shown in Fig. 1, Eqs. 3 take the form

$$\begin{array}{cc}E_s\overline{u}_s(x)\hfill & =iv_Fk_{Fx}\frac{d}{dx}\overline{u}_s(x)+\mathrm{\Delta }_{s\overline{s}}(k_{})\overline{v}_{\overline{s}}(x)\hfill \\ E_s\overline{v}_{\overline{s}}(x)\hfill & =iv_Fk_{Fx}\frac{d}{dx}\overline{v}_{\overline{s}}(x)+\mathrm{\Delta }_{\overline{s}s}^{}(k_{})\overline{u}_s(x)\hfill \end{array}.$$
(4)
When a beam of electrons is incident from the normal metal to the insulator, with momentum $`k`$, the general solution of Eqs. (4), is the two-component wave function $`\mathrm{\Psi }_I=(u_{[]},v_{[]})`$ which for $`x<0`$ is written as

$$\mathrm{\Psi }_I=\left(\begin{array}{cc}1\hfill & \\ 0\hfill & \end{array}\right)e^{ixk_{Fx}}+a_{[]}\left(\begin{array}{cc}0\hfill & \\ 1\hfill & \end{array}\right)e^{ixk_{Fx}}+b_{[]}\left(\begin{array}{cc}1\hfill & \\ 0\hfill & \end{array}\right)e^{ixk_{Fx}},$$
(5)
where $`a_{[]},b_{[]}`$, are the amplitudes for Andreev and normal reflection for spin-up(-down) quasiparticles, $`k_{Fx}`$ is the $`x`$-component of the Fermi wave vector $`k_F=\sqrt{\frac{2mE_F}{\mathrm{}^2}}`$. The wave vector of quasiparticles in the normal-metal, and the wave vector of the electron-like and hole-like quasiparticles in the superconductor are set equal. Since the translational symmetry holds in the $`y`$-axis direction, the momentum parallel to the interface is conserved. Using the matching conditions of the wave function at $`x=0`$, $`\mathrm{\Psi }_I(0)=\mathrm{\Psi }_{II}(0)`$ and $`\mathrm{\Psi }_{II}^{}(0)\mathrm{\Psi }_I^{}(0)=(2mV/\mathrm{}^2)\mathrm{\Psi }_I(0)`$, the Andreev and normal reflection amplitudes $`a_{[]},b_{[]}`$ for the spin-up(-down) quasiparticles are obtained as

$$a_{[]}=\frac{4n_+}{4+z_{[]}^2z_{[]}^2n_+n_{}\varphi _{}\varphi _+^{}},$$
(6)
$$b_{[]}=\frac{(2iz_{[]}+z_{[]}^2)+(2iz_{[]}+z_{[]}^2)n_+n_{}\varphi _{}\varphi _+^{}}{4+z_{[]}^2z_{[]}^2n_+n_{}\varphi _{}\varphi _+^{}},$$
(7)
where $`z_0=\frac{mV}{\mathrm{}^2k_s}`$, $`z_{[]}=\frac{2z_0}{k_{Fx}}`$. The BCS coherence factors are given by

$$u_\pm ^2=[1+\sqrt{E^2|\mathrm{\Delta }_\pm |^2}/E]/2,$$
(8)
$$v_\pm ^2=[1\sqrt{E^2|\mathrm{\Delta }_\pm |^2}/E]/2,$$
(9)
and $`n_\pm =v_\pm /u_\pm `$. The internal phase coming from the energy gap is given by $`\varphi _\pm =[\mathrm{\Delta }_\pm /|\mathrm{\Delta }_\pm |]`$, where $`\mathrm{\Delta }_+`$ ($`\mathrm{\Delta }_{}`$), is the pair potential experienced by the transmitted electron-like (hole-like) quasiparticle.

According to the BTK formula the conductance of the junction, $`\overline{\sigma }_{[]}(E_s,k_{})`$, for up(down) spin quasiparticles, is expressed in terms of the probability amplitudes $`a_{[]},b_{[]}`$ as

$$\overline{\sigma }_{[]}(E_s,k_{})=1+|a_{[]}|^2|b_{[]}|^2.$$
(10)
The tunneling conductance, normalized by that in the normal state is given by

$$\sigma (E)=\sigma _{}(E_{})+\sigma _{}(E_{}),$$
(11)
$$\sigma _{[]}(E_s)=\frac{1}{R_N}_{k_{}^{max}}^{k_{}^{max}}𝑑k_{}\overline{\sigma }_{[]}(E_s,k_{}),$$
(12)
where

$$R_N=_{k_{}^{max}}^{k_{}^{max}}𝑑k_{}[\sigma _N_{}(k_{})+\sigma _N_{}(k_{})],$$
(13)
$$\sigma _{N_{[]}}(k_{})=\frac{4\lambda _1}{4+z_{[]}^2}.$$
(14)
## III pairing states and Fermi surface line shape

For the spin triplet pairing state the Cooper pairs have spin $`1`$ degree of freedom. The gap function is a $`2\times 2`$ symmetric matrix which in the spin space can be written as

$$\widehat{\mathrm{\Delta }}(𝒌)=i\sigma _y(𝒅\widehat{𝝈}),$$
(15)
where $`\widehat{𝝈}`$ denotes the Pauli matrices and $`𝒅`$ is a vector which defines the axis along which the Cooper pairs have zero spin projection. In the following we will take $`𝒅\widehat{𝒂}`$, i.e., parallel to the chains. In that case $`\mathrm{\Delta }_{}=\mathrm{\Delta }_{}=0`$, while $`\mathrm{\Delta }_{}=\mathrm{\Delta }_{}=\mathrm{\Delta }(k_{})`$.

The Fermi surface (FS) consists of two branches and is open in the $`k_y`$ direction, as shown in Fig. 2.

We consider the following pairing cases

a) In case of $`p_x`$-wave, $`p_y`$-wave superconductor

$$\mathrm{\Delta }_{p_x}=\mathrm{\Delta }_0\mathrm{sin}(2k_xa),\mathrm{\Delta }_{p_y}=\mathrm{\Delta }_0\mathrm{sin}(k_ya).$$
(16)
$`\mathrm{\Delta }_{p_x}`$ changes its sign along the FS in the $`k_x`$ direction while the $`\mathrm{\Delta }_{p_y}`$ changes its sign along the FS in the $`k_y`$ direction as seen in Figs. 2(a) and 2(d).

b) In case of $`d_x`$-wave, $`d_y`$-wave superconductor

$$\mathrm{\Delta }_{d_x}=\mathrm{\Delta }_0\mathrm{cos}(2k_xa),\mathrm{\Delta }_{d_y}=\mathrm{\Delta }_0\mathrm{cos}(k_ya)$$
(17)
$`\mathrm{\Delta }_{d_x}`$ changes its sign along the FS in the $`k_x`$ direction while the $`\mathrm{\Delta }_{d_y}`$ changes its sign along the FS in the $`k_y`$ direction as seen in Figs. 2(b) and 2(e).

c) In case of $`f_x`$-wave, $`f_y`$-wave superconductor

$$\mathrm{\Delta }_{f_x}=\mathrm{\Delta }_0\mathrm{sin}(4k_xa),\mathrm{\Delta }_{f_y}=\mathrm{\Delta }_0\mathrm{sin}(2k_ya)$$
(18)
$`\mathrm{\Delta }_{f_x}`$ changes its sign along the FS in the $`k_x`$ direction while the $`\mathrm{\Delta }_{f_y}`$ changes its sign along the FS in the $`k_y`$ direction as seen in Figs. 2(c) and 2(f).

d) In case of $`d_{xy}`$-wave, $`d_{x^2y^2}`$-wave superconductor

$`\mathrm{\Delta }_{d_{xy}}=`$ $`\mathrm{\Delta }_0\mathrm{sin}(2k_xa)\times \mathrm{sin}(k_ya),`$ (19)
$`\mathrm{\Delta }_{d_{x^2y^2}}=`$ $`\mathrm{\Delta }_0(\mathrm{cos}(2k_xa)\mathrm{cos}(k_ya)),`$ (20)
$`\mathrm{\Delta }_{d_{xy}}`$, $`\mathrm{\Delta }_{d_{x^2y^2}}`$ change sign along the FS as seen in Figs. 3(a) and 3(b).

## IV Tunneling conductance characteristics

In Figs. 4-7 we plot the tunneling conductance $`\sigma (E)`$ as a function of $`E/\mathrm{\Delta }_0`$ for various values of $`z_0`$, for different orientations normal to the $`a`$ or $`b`$ axis. The pairing symmetry of the superconductor is $`p_x`$-wave, $`p_y`$-wave, in Fig. 4, $`d_x`$-wave, $`d_y`$-wave, in Fig. 5, $`f_x`$-wave, $`f_y`$-wave, in Fig. 6, $`d_{xy}`$-wave, $`d_{x^2y^2}`$-wave, in Fig. 7.

The conductance peak is formed in the electron tunneling along the $`a`$ or $`b`$ axis, when the transmitted quasiparticles experience different sign of the pair potential. Also the line shape of the spectra is sensitive to the presence or absence of nodes of the pair potential on the Fermi surface.

For the $`p_x`$-wave case, when the $`a`$-axis of the crystal is at right angle to the interface, the horizontal line in Fig. 2 (a), representing the scattering process, connects points of the FS with opposite sign of the pair potential for $`\pi /a<k_y<\pi /a`$. As a result a peak exists in the conductance spectra, at $`E=0`$ as seen in Fig. 4 (a) for $`z_0=2.5`$. The height of the ZEP is proportional to the range of $`k_y`$ for which sign change occurs and it is expected to have its maximum value when the $`a`$-axis of the crystal is normal to the interface, since for this orientation the transmitted quasiparticles feel a different sign of the pair potential for all $`\pi /a<k_y<\pi /a`$. On the other hand, when the $`a`$ axis is along the interface, then the vertical line in Fig. 2 (a) connects points of the FS with the same sign of the pair potential for $`\pi /2a<k_x<\pi /2a`$, i.e., there is no $`k_x`$ for which, the transmitted quasiparticles feel the sign change of the pair potential, and no ZEP is formed as seen in Fig. 4 (b). For the $`p_x`$-wave case the nodes of the pair potential do not intersect the FS and the line shape of the spectra is U-like as in the case of the $`s`$-wave superconductor.

The situation is opposite in the $`p_y`$-wave case, where the scattering process for the surface orientation normal to the $`a`$ ($`b`$)-axis, described by the horizontal (vertical) line in Fig. 2 (d) connects points of the FS with the same (opposite) sign. As a consequence for the surface orientation normal to the $`a`$-axis, there is no $`k_y`$ for which the transmitted quasiparticles feel the sign change of the order parameter, and no ZEP is formed, as seen in Fig. 4 (c), while for the interface along the $`a`$-axis the transmitted quasiparticles feel the different sign of the pair potential, for all $`\pi /2a<k_x<\pi /2a`$, and a ZEP exists as seen in Fig. 4 (d). For the $`p_y`$-wave case the nodes of the pair potential intersect the FS and the spectra has a V-shaped form as in the case of $`d`$-wave superconductors.

For the $`d`$-wave case, and for tunneling along the $`a`$ or $`b`$-axis, the scattering process, described with the horizontal or vertical line in Figs. 2(b) and 2(e) connects points of the FS with the same sign in all cases. This means that the pair potential does not change sign and no ZEP is formed as seen in Fig. 5. Also due to the presence of nodes of the pair potential along the FS, the line shape of the spectra is V-like.

For the $`f`$-wave case, although the pair potential has a different structure, the lines in Figs. 2(c) and 2(f) connect points of the FS with the same sign-change of the pair potential as in the $`p`$-wave case. As a result the tunneling shows ZEP as in the $`p`$-wave case, due to the sign change of the transmitted quasiparticles. However, in all cases in the $`f`$-wave case the pair potential intersects the FS and nodes are formed. As a consequence the line shape of the spectra is V-like, as seen in Fig. 6 unlike to $`p_x`$-wave case, where the pair potential is nodeless and the tunneling spectra has a U-shaped form. The conclusion is that the line shape of the spectra can be used to distinguish the $`p_x`$ from the $`f`$-wave pairing state. Our results are comparable to that of Tanuma et al.  although their calculation was based on a different model, i.e., the extended Hubbard model on a quasi-one dimensional lattice at quarter-filling.

For the $`d_{xy}`$-wave ($`d_{x^2y^2}`$-wave) case, and for tunneling along the $`a`$ or $`b`$-axis, the scattering process, described with the horizontal or vertical line in Figs. 3(a) and 3(b) connects points of the FS with different (the same) sign. This means that the pair potential felt by the transmitted quasiparticles changes (conserves) sign during the scattering process and a ZEP (no ZEP) is formed at $`E=0`$ for the $`d_{xy}`$-wave ($`d_{x^2y^2}`$-wave) case as seen in Fig. 7. Also due to the presence of nodes of the pair potential along the FS, the line shape of the spectra is V-like. In addition the tunneling spectra for the $`d_{xy}`$-wave, for tunneling direction along the $`a`$ or $`b`$ axis, is equivalent to the tunneling spectra for $`d_{x^2y^2}`$-wave with surface orientation tilted by $`\pi /4`$, i.e., along the nodal direction.

In the metallic limit ($`z_0=0`$), $`\sigma (E)`$ has the same form for each pairing potential independently from the orientation.

In table I we summarize the results concerning the presence or absence of ZEP, and also the overall line shape of the tunneling spectra for various pairing potentials.

## V Bound state energies

These features are explained if we calculate the energy of the midgap state, which is given for large $`z_0`$ by the value in which the denominator of Eqs. (6) and (7) vanishes. The equation giving the energy peak level is written as

$$\varphi _{}\varphi _+^{}n_+n_{}|_{E=E_p}=1.0.$$
(21)
In the $`p_x`$-wave case, for surface orientation normal to $`a`$-axis, this equation has the solution $`E=0`$, for $`\pi /a<k_y<\pi /a`$, since $`n_+n_{}|_{E=0}=1`$, and also the transmitted quasiparticles feel a different sign of the pair potential, i.e., $`\varphi _{}\varphi _+^{}|_{E=0}=1`$. When a midgap state exists the tunneling conductance $`\overline{\sigma }_s(E,k_y)`$ is equal to $`2`$ and the peek in $`\sigma (E)`$ seen in Fig. 4(a), is due to the normal state conductance $`R_N`$ in Eq. 12, which depends inversely on the $`z_0^2`$ for large $`z_0`$. For surface orientation normal to $`b`$-axis the range of $`k_x`$ for which Eq. 21 has solutions collapses to zero, and no bound states are formed. Then $`\sigma (E)`$ goes to zero as $`1/z_0^2`$ and there is no conductance peak. In the $`p_y`$-wave case for surface orientation normal to $`a`$-axis Eq. 21 has no solutions because $`\varphi _{}\varphi _+^{}|_{E=0}=1`$, while for surface orientation normal to $`b`$-axis a bound state is formed for all $`k_x`$ in the interval $`\pi /2a<k_x<\pi /2a`$. For the $`d`$-wave case and for tunneling along the $`a`$ or $`b`$ axis of the crystal, the condition $`\varphi _{}\varphi _+^{}|_{E=0}=1`$ is not satisfied and no bound states are formed. For the $`f`$-wave case that condition is the same as in the $`p`$-wave case, and thus the ZEP is formed as in the $`p`$-wave pairing state. For the $`d_{xy}`$-wave ($`d_{x^2y^2}`$-wave) case that condition $`\varphi _{}\varphi _+^{}|_{E=0}=1`$ is (is not satisfyied for tunneling along the $`a`$ and $`b`$ axis and hence the ZEP (no ZEP) is formed for both tunneling directions.

## VI magnetic field effect

In this section we describe the effect of the external magnetic field $`H`$ in the direction parallel to the chains, in the spectra. The effect of the magnetic field depends on the spin, of the quasiparticles. The tunneling conductance is given by

$$\sigma (E)=\sigma _{}(E\mu _BH)+\sigma _{}(E+\mu _BH).$$
(22)
In Fig. 8 the tunneling conductance $`\sigma (E)`$ is plotted for a fixed magnetic field $`\mu _BH/\mathrm{\Delta }_0=0.2`$, and barrier strength $`z_0=2.5`$. The pairing symmetry of the superconductor is $`p`$, $`d`$, $`f`$, $`d_{xy}`$ respectively.

The magnetic field splits symmetrically the tunneling spectrum which is a linear superposition of the spectra for spin up(down) quasiparticles. The amplitude of the splitting depends linearly on the magnetic field $`H`$. For the case of $`\mu _BH/\mathrm{\Delta }_0=0.2`$, seen in Fig. 8 the spin up(down) part of the spectra partially overlaps while for larger values of the magnetic field the spin up and down branches are well separated. In the latter case the right(left) branch of the spectra originates from spin up(down) quasiparticle spectra $`\sigma _q_{}(E\mu _BH)(\sigma _q_{}(E+\mu _BH))`$.

The condition for the formation of bound states is slightly modified under the presence of a magnetic field to $`E\mu _BH=0`$, for the spin-up region, and $`E+\mu _BH=0`$, for the spin-down, from the corresponding $`E=0`$ in the absence of any field. So the multiplication of the bound states and the presence of magnetic field results into the appearance of double peak in the conductance spectra.

## VII Conclusions

We calculated the tunneling conductance in normal-metal/insulator/quasi-one-dimensional superconductor, using the BTK formalism. We showed that the ZEP appears in the electron tunneling along the $`a`$ or $`b`$ axis, when the transmitted quasiparticles experience a different sign of the pair potential. Also the line shape of the spectra is sensitive to the presence or absence of nodes of the pair potential on the Fermi surface, and results to a U-shaped structure for the $`p_x`$-wave case and to a V-shaped one for the $`d`$ and $`f`$-wave cases.

The ZEP are due to the formation of bound states within the gap. The calculation of the conductance $`\overline{\sigma }_s(E,k_y)`$, for which bound state occurs, shows an enhancement at the bound state energy. The effect of the magnetic field in the direction parallel to the chains is to split linearly the ZEP.

Throughout this paper the spatial variation of the order parameter near the surface, which depends on the boundary orientation, is ignored for simplicity. However, since the features presented here are intrinsic which are generated by the existence of surface bound states, the essential results do not change qualitatively. Also we assumed perfectly flat interfaces in the clean limit, so any impurity scattering and the effect of the surface roughness are ignored.

2


format: text

text:
# Inverse cascade in Charney-Hasegawa-Mima turbulence

## Abstract

The inverse energy cascade in Charney-Hasegawa-Mima turbulence is investigated. Kolmogorov law for the third order velocity structure function is shown to be independent on the Rossby number, at variance with the energy spectrum, as shown by high resolution direct numerical simulations. In the asymptotic limit of strong rotation, coherent vortices are observed to form at a dynamical scale which slowly grows with time. These vortices form an almost quenched pattern and induce strong deviation form Gaussianity in the velocity field.

The existence of an inverse cascade is the most remarkable property of two dimensional turbulence. It was predicted by Kraichnan  for Navier–Stokes equation: as a consequence of inviscid enstrophy conservation, energy is forced to flow to large scales. The inverse cascade can be sustained only in presence of an external forcing injecting energy at a characteristic scale into the system. At scales larger than the forcing the turbulent flow is essentially random with Gaussian velocity difference statistics following Kolmogorov scaling . Thus, the presence of the inverse energy cascade prevents the formation of the large scale coherent structures observed in the case of decaying turbulence .

Large scale coherent structures in presence of the inverse cascade have been observed only in presence of a characteristic scale breaking scale invariance. A well known example is the so-called Bose-Einstein condensation, when the energy accumulates at the largest available scale  forming vortices at the system size.

Another example of vortex formation is the quasi-crystallization phenomenon observed in Charney-Hasegawa-Mima (CHM) turbulence, a paradigm for both geostrophic motion in planetary atmospheres  and drift-wave turbulence in magnetically confined plasma . For the stream function $`\psi (𝐱,t)`$ CHM equation is written as:

$$\frac{}{t}(^2\psi \lambda ^2\psi )+J(^2\psi ,\psi )=D+F$$
(1)
where $`J`$ denotes the Jacobian and $`D`$ and $`F`$ are damping and forcing respectively. In this case, vortices have been observed to form, in a quasi-crystal structure, at the intrinsic scale $`1/\lambda `$, corresponding to the Rossby deformation radius in the atmosphere  or to the effective ion Larmor radius in plasma .

In this Letter we focus on dynamics on scales much larger than $`\lambda ^1`$. In this regime there is no intrinsic scale involved in the evolution of the system, which nevertheless exhibits formation of strong vortices. In this case we observe that the scale of vortices is a dynamical one which increases in time as a consequence of vortex merging. The characteristic time of evolution slows down leading, in the limit of large Reynolds numbers, to a disordered pattern of quenched vortices. Despite the presence of strong vortices, we find that the two-dimensional $`3/2`$ Kolmogorov law for the third-order velocity structure function holds, independently on the value of $`\lambda `$. As a consequence, the kinetic energy spectrum follows Kolmogorov scaling but with a different constant with respect to the Navier–Stokes turbulence.

The CHM equation (1) has two quadratic inviscid invariants corresponding to total energy

$$E=E_k+E_\lambda =\frac{1}{2}(\psi )^2+\lambda ^2\psi ^2$$
(2)
where $`\mathrm{}`$ denotes spatial average, and total enstrophy

$$Z=Z_k+Z_\lambda =\frac{1}{2}(^2\psi )^2+\lambda ^2(\psi )^2$$
(3)
Both the inviscid invariants consist of two terms, the first corresponding to kinetic contribution and the second to potential one. The kinetic terms are, by definition, the only ones which survive in the Navier–Stokes limit $`\lambda 0`$.

The range of scales are separated by the characteristic wavenumber $`\lambda `$. For $`k\lambda `$ the kinetic contributions dominate in (2-3), at very large scales $`k\lambda `$ the leading terms are the potential ones. In the following we will assume that the forcing $`F`$ is limited to a narrow band of wavenumber around $`k_f`$. This will be the other relevant wavenumber in our problem.

Kolmogorov-like dimensional analysis can be easily extended to the present problem . If $`k_f\lambda `$ we recover the well-known energy spectra for two-dimensional Navier–Stokes turbulence with energy spectrum $`E(k)k^{5/3}`$ and $`E(k)k^3`$ for $`\lambda <k<k_f`$ and $`k>k_f`$ respectively. When $`k_f\lambda `$ one obtains the prediction $`E(k)k^{11/3}`$ for $`k<k_f`$ and $`E(k)k^5`$ for $`k_f>k>\lambda `$. Dimensionally predicted spectra have been confirmed by direct numerical simulations  but little is known about structure functions and probability density functions.

The starting point for a statistical approach to turbulent cascade is the energy flux, written in the physical space as

$$\epsilon (\mathrm{})=\frac{}{t}\frac{1}{2}\psi (𝐱)\psi (𝐱+\mathrm{})+\lambda ^2\psi (𝐱)\psi (𝐱+\mathrm{})|_{NL}$$
(4)
where the subscript $`NL`$ stands for the nonlinear contribution of (1) to the time derivative. Making use of (1) and of integration by parts we easily obtain

$$\epsilon (\mathrm{})=\frac{1}{2}\psi (𝐱)J(\psi (𝐱+\mathrm{}),^2\psi (𝐱+\mathrm{}))+(𝐱𝐱+\mathrm{})$$
(5)
in which the parameter $`\lambda `$ has formally disappeared. As a consequence we can make use of the well known result for Navier–Stokes (i.e. $`\lambda =0`$) and write

$$(\delta u_L)^3=\frac{3}{2}\epsilon \mathrm{}$$
(6)
where $`\delta u_L`$ represents the longitudinal increment of the velocity $`𝐯=(_y\psi ,_x\psi )`$ and $`\epsilon `$ is the energy input due to the forcing. We thus have a new degree of universality for the 3/2 Kolmogorov law in two-dimensional turbulence, with respect to the class of equations (1) parameterized by $`\lambda `$. From a dimensional point of view, (6) implies a scaling exponent $`h=1/3`$ for velocity increments and thus a scaling exponent $`4/3`$ for $`\psi `$, as in Navier–Stokes turbulence. From (2) one obtains the different predictions for the spectrum discussed above. In particular, the kinetic spectrum has the form

$$E_k(k)=C_\lambda \epsilon ^{2/3}k^{5/3}$$
(7)
for any value of $`\lambda `$, but with a Kolmogorov constant $`C_\lambda `$ which, in principle, can depend on $`\lambda `$. A simple physical argument for this dependency is as follows. The scaling of the eddy turnover time depends on the scaling exponent. In the kinetic limit $`k\lambda `$ one has the standard Kolmogorov scaling $`\tau (k)\epsilon ^{1/3}k^{2/3}`$ . On the other hand, in the potential limit $`k\lambda `$, (1) gives $`\tau (k)(\lambda /k)^2\epsilon ^{1/3}k^{2/3}`$. Thus, for $`k\lambda `$, the efficence of energy transfer is reduced and one expects a larger value of the constant in (7).

We have numerically investigated the inverse cascade in the potential energy regime by direct numerical simulations. In order to avoid complications induced by the crossover from the kinetic domain to the potential domain, we study the system in the limit $`\lambda \mathrm{}`$. This is to be seen only as a formal procedure, equivalent to considering wavenumbers much smaller than $`\lambda `$, which physically might be the case for magnetized plasma in presence of a strong magnetic field. Indeed this limit provides us with a model suitable for any $`\lambda k_f`$: because the energy is transferred to large scale, the dominance of the potential term will be assured in all the inertial range. In the limit $`\lambda k_f`$, rescaling the time $`tt/\lambda ^2`$, one obtains the so-called asymptotic model

$$\frac{\psi }{t}=J(^2\psi ,\psi )+D+F$$
(8)
for which the conserved quantities become

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2}}\psi ^2`$ (9)
$`Z`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\psi )^2`$ (11)
We have integrated (8) with a standard pseudo-spectral code in a double periodic domain of size $`L=2\pi `$ at resolution $`N=512`$. The forcing is white in time in a narrow band of wavenumbers around $`k_f=160`$. The dissipative term in (8) has the role of removing potential enstrophy at small scales and, as customary, it is numerically substituted by a hyperviscous term (of order $`8`$ in our simulations). Time evolution is obtained by a standard second-order Runge-Kutta scheme starting from a zero initial condition. The run is stopped at a given time $`T`$ at which the energy containing scales are still much smaller than the computational box in order to avoid condensation effects  (see Figure 1). All the results discussed in the following are taken after averaging over $`14`$ independent realizations.

The limitation in the resolution ($`N=512`$) is due to the discussed scaling of the characteristic time. Even with this moderate resolution, the ratio of the large scale characteristic time with the forcing scale time is about $`2000`$ and thus time evolution is very expensive ($`10^6`$ time steps for each realization). In the case of Navier–Stokes turbulence ($`\lambda =0`$) this would correspond to an integration covering about $`5`$ decades of inertial range.

In Figure 1 we plot the potential energy spectrum $`E_\lambda (k)`$ at two different times. The scaling exponent $`11/3`$ is clearly visible even if some accumulation at the largest mode is evident. This accumulation is not due to condensation as it is still well below the largest mode and it moves in time. We think that the existence of this “bump” is a genuine effect, probably due to the rapid growth of characteristic times and to the presence of intense vortices, as discussed below. The energy flux $`\epsilon 1.24`$ is estimated by the plateau of the energy flux shown in the inset.

The “$`3/2`$” law for the longitudinal velocity structure function is shown in Figure 2, also plotted at two different times. The compensation with the theoretical prediction (6) is remarkable, taking into account the limited resolution of our runs. As expected, the extension of the inertial range increases with time without changing small scale statistics. The oscillations observed at small scales are due to the contamination of the forcing. A similar effect was observed also in NS simulations.

As discussed above, the fact that the “$`3/2`$” law is independent on $`\lambda `$ (and thus the velocity scaling exponent has always the Kolmogorov value $`h=1/3`$) do not imply that the statistics, and in particular the form of the pdf of velocity differences, is the same as for Navier–Stokes equation. For example, in Figure 1 we also plot the kinetic energy spectrum $`E_k(k)`$ at final time $`t=10^2`$. The scaling exponent is compatible with the Kolmogorov value $`5/3`$ as predicted by (7), but the Kolmogorov constant $`C_\lambda 11`$ is about two times that of Navier–Stokes . A larger constant means a suppression of the energy flux which is a direct consequence of the dilatation of the characteristic times.

A more significant difference with respect to Navier–Stokes inverse cascade is the presence of strong vortices, as shown in Figure 3. Vortices in CHM turbulence are injected at the forcing scale and they organize themselves to form a random pattern on the characteristic scale $`\mathrm{}_E`$. This dynamical scale is associated to the peak of the spectrum of Figure 1. Vortex dynamics slows down as $`\mathrm{}_E`$ increases as $`\tau (\mathrm{}_E)\mathrm{}_E^{8/3}`$. Thus, in the limit of large Reynolds number the system will end in a disordered pattern of quenched vortices forming a kind of “turbulent glass”.

An important consequence of the presence of strong vortices is that the statistics of the velocity field strongly deviates from Gaussianity. In Figure 4 we plot the pdf of longitudinal and transverse velocity differences at three different scales within the inertial range. The effect of vortices is evident by the presence of large “wings” in the tails, in particular on the transverse velocity differences which are more sensible to a rotating structure.

In conclusion, we have shown that Kolmogorov “3/2” law for two-dimensional energy cascade in Charney-Hasegawa-Mima turbulence is independent on the value of the intrinsic scale $`\lambda `$. Velocity statistics satisfies Kolmogorov scaling with non-universal coefficients. In the asymptotic limit $`\lambda \mathrm{}`$ the Kolmogorov constant is found to be about $`2`$ times the Navier–Stokes case . Strong coherent vortices are found to emerge at the forcing scale and aggregate to form a pattern of quenched vortices at large scale. As a consequence of the presence of vortices, strong deviations from Gaussianity are observed in the velocity field.

###### Acknowledgements.

This work was partially supported by MURST Cofin 2001 n.2001023848. We acknowledge the allocation of computer resources from INFM Progetto Calcolo Parallelo.