text:
# On the Zero-Bias Conductance Peak in the Tunneling Spectroscopy

## I Introduction

Since the pioneering work of Giaever,giaever  the tunneling measurement has been a major experimental method for investigations into the electronic states of condensed matter systems.wolf  In the simplest setup, a metal with known spectral property is made in contact with a material X, forming an NX junction so that the electronic states of X can be probed. For many years, despite the fact that many insights into the spectral properties of many states have been gained from the differential conductance ($`\mathrm{d}I/\mathrm{d}V`$) curves obtained from tunneling measurements, nonetheless, unlike many other experiments, it is fair to say that there is no clear and solid statement as to exactly what bulk properties are being probed in tunneling measurements. For example, it is known that in neutron scattering experiments, the neutron intensity is a measure of the imaginary part of the bulk spin susceptibility, $`Im\chi (𝐤,\omega )`$; no similar statement has ever been firmly established for tunneling measurements.

The difficulty for establishing the relation between the tunneling conductance and the bulk quantities can be traced back to the very existence of the junction interface. It has been realized that the presence of the interface can change the conductance curve dramatically. A well-known example is the zero-bias conductance peak (ZBCP) observed in the tunneling spectra when X is a $`d`$-wave superconductor (ND junction) in (110) direction.hu  The appearance of the ZBCP is entirely tied up with the presence of the interface and its orientations, and therefore can not be obtained by simple calculations based on bulk density of states.

Recent theoretical analyses of the ZBCP have been mostly concentrated on the ND junctions. Furthermore, they are based largely on the standard BTK theory.btk  In the continuum limit, analytic expressions of the differential conductance for general orientations of the interface were obtained. Numerical calculations were later carried out for the BTK theory in the lattice version.Tanuma ; datta ; cl  While these works have supplied insights into the ZBCP, they are, however, specifically designed for studying the ND junction. Moreover, because the relation of the conductance curve to the bulk quantities was not clearly manifested, essentially the numerical computation had to be done individually for each interface orientation. Another technical inconvenience is that the BTK theory is a mean-field theory based on solving the mean-field quasi-particle wavefunctions, it is thus difficult in this formulation to take into account the effects of interaction systematically. To extend into the study of other systems, especially those with strong correlations where almost all relevant models are on discrete lattices, it is therefore an urgent need to have a formulation which can go beyond the mean-field BTK formulations. As an illustration of our approach, in this paper we will focus on mean-field analysis of several tunneling problems. The effects of fluctuations and interactions will be discussed elsewhere.

In this paper, we shall adopt an approach that is based on the non-equilibrium Keldysh-Green’s function formalism. In the lowest order approximation, we will be able to express the differential conductance entirely in terms of bulk Green’s functions and include the interface effects. Thus, the relation of the conductance curve to the bulk quantities is clearly manifested. The tunneling between N and X will be treated as a perturbation, so that in the zeroth order the Green’s function is the mean-field half-space Green’s function that resides only on the semi-infinite plane and satisfies the boundary conditions to be specified later. Based on the half-space Green’s function, $`g`$, higher order corrections can be systematically constructed.cl ; luke ; cuevas  In particular, a class of infinite series in $`g`$, which consists of all elastic tunneling processes in the perturbation theory, will be considered and summed to all orders for calculating the current across the junction.cl ; caroli ; yan  To fully take into account the tight-binding nature of the problem, we shall employ discrete models for both the materials N and X and the tunnel junctions. Thus the essential quantity to be calculated is the half-space lattice Green’s function for the X state. In resemblance to the conventional method of image, we express the half-space Green’s function in terms of the bulk Green’s functions propagating from the real source and a fictitious image source

$`g=G_{\text{real}}G_{\text{imag}}\times \alpha `$ (1)
with the factor $`\alpha `$ accounting for the boundary conditions. In this picture the half-space Green’s function is decomposed into two parts: the real-source part comes solely from the bulk and hence reveals purely the bulk properties, the image part contains all interface effects which are encoded in the factor $`\alpha `$. In this way, the interface effects are clearly identified in the course of the analysis and one can pinpoint any departure from the bulk property.

The factor $`\alpha `$ can be expressed in terms of bulk Green’s function. Right on the interface, it is found

$`\alpha _0=G^1(d)G(d).`$ (2)
Here $`d`$ is an effective lattice constant whose precise meaning will be explained in below. Clearly, the tunneling spectrum can be classified according to whether the reflection symmetry is broken or not. In the case when reflection symmetry is broken with respect to the interface, one has $`G(d)G(d)`$, hence $`\alpha _0`$ is not unity, possible zero modes may arise due to the presence of zeros in the denominator of the left hand side. The number of localized zero mode is thus determined by the order of zeros in the bulk Green’s function $`G(d)`$. In the lowest order approximation, the differential conductance is given by the local density of state at the interface

$`\mathrm{d}I/\mathrm{d}V`$ $``$ $`{\displaystyle \underset{𝐤\sigma }{}}Im\{g_0(𝐤,eV)\},`$ (3)
where $`g_0`$ is $`g`$ of Eq. (1) evaluated at the interface and $`e`$ is the charge of an electron. Since $`\alpha _0`$ can be expressed entirely in terms of bulk Green’s functions, this is then the relation between the bulk quantities and the differential conductance alluded to earlier.

This paper is organized as follows. In Sec. II, we outline the theoretical formulation and derive the generalized method of image for discrete lattices. In Sec. III this method is applied to the study of tunneling spectroscopies for various systems. We first study the ND junctions at various surface orientations and examine the doping dependence of the ZBCP using mean-field slave boson theory. We then study the effects of applied magnetic fields perpendicular to the $`ab`$ plane. A one dimension model, based on the structure of polyacetylene, is then studied in Sec. III.3. On the basis of this model, we further apply this method to investigate tunneling into $`d`$-density-wave states and graphite sheets. We conclude in Sec. IV with some comments on the significance and further applications of our formulation. The Appendix describes techniques for deriving the current expressions for the tunnel junctions studied in the text.

## II Theoretical Formulation and Generalized method of image

### II.1 Theoretical model

We start by modeling the planar junction. As illustrated in Fig. 1, the tunnel junction consists of two truncated two dimensional (2D) lattices connected through a tunnel barrier, with the left half the normal (N) electrode ($`\mathrm{}<xa_L`$, $`a_L`$ is the lattice constant) and the right the test (X) electrode ($`0x<\mathrm{}`$).3D  We take the interface the $`y`$ direction. The total Hamiltonian of the system thus comprises two parts: the Hamiltonian $`H_0=H_L+H_R`$ for the left and right electrodes, and the tunneling Hamiltonian which connects the surface points at $`x=a_L`$ and $`x=0`$

$`H_T={\displaystyle \underset{y_l,y_r,\sigma }{}}t(|y_ly_r|)c_{l\sigma }^{}c_{r\sigma }+\mathrm{h}.\mathrm{c}.`$ (4)
Here $`\sigma `$ are spin indices, and $`y_l`$, $`y_r`$ are the $`y`$-coordinates of the surface sites on the left and right electrodes; $`c_{l\sigma }`$, $`c_{r\sigma }`$ are the corresponding electron annihilation operators. $`t`$ is the tunneling amplitude whose magnitude models the barrier height in the tunnel junction. Since all points over the interface layers contribute to the tunneling process, one has to sum over all interface sites. Suppose the chemical potentials on the left and the right electrodes are $`\mu _L`$ and $`\mu _R`$, respectively, the total grand Hamiltonian is then given by

$`K`$ $`=`$ $`(H_L\mu _LN_L)+(H_R\mu _RN_R)+H_T`$ (5)
$``$ $`K_0+H_T.`$
The difference $`(\mu _L\mu _R)`$ is fixed to be the voltage drop $`eV`$ across the junction.

To calculate the tunneling current, we choose the unperturbed state to be the ground state of $`K_0`$ and adiabatically turn on $`H_T`$. In the Heisenberg picture, the tunneling current is obtained from the time rate of change of the particle number $`N_L`$ of the left electrodemahan  (we set $`\mathrm{}=1`$ throughout)

$`I(t)`$ $`=`$ $`+ie[N_L,H_T]`$ (6)
$`=`$ $`+ie{\displaystyle \underset{y_l,y_r,\sigma }{}}\left\{tc_{l\sigma }^{}c_{r\sigma }t^{}c_{r\sigma }^{}c_{l\sigma }\right\}.`$
The expectation values $`\mathrm{}`$ here represent the ensemble average $`Tr[Z^1\mathrm{exp}(\beta K_0)\mathrm{}]`$.

In actual experiments, the normal metal on the left electrode could be in any orientations, and the detail connection between the two lattices may also cause complications in the tunneling spectroscopy. To be definite, however, in our model we fix the lattice on the left side at (100) orientation and connect its boundary sites to those of the right at $`x=0`$ (Fig. 1$`(a)`$). As one can observe easily, the system is translational invariant along the interface direction with period $`a_L`$, the lattice constant of the left side. We exploit this symmetry by making a partial Fourier transformation along the interface direction in Eq. (6) and arriving at

$`I(t)=+ie{\displaystyle \underset{k_y,\sigma }{}}t(k_y)c_{l\sigma }^{}(k_y)c_{r\sigma }(k_y)+\mathrm{h}.\mathrm{c}.`$ (7)
The range of $`k_y`$ is determined by the periodicity of the interface sites along $`y`$ direction, hence

$`\pi /a_L<k_y\pi /a_L.`$ (8)
We emphasize that the problem is now effectively one dimensional: in Eq. (7) different $`k_y`$ modes are decoupled completely. Moreover, only surface quantities are involved. These are very appealing features especially for the feasibility of our method of image, as we will discuss in the following section.

In the Keldysh Green’s function formulation the time evolution of the density matrix can be formally solved as a closed time-ordered path integral,keldysh  the expectation value $`c_{l\sigma }^{}(k_y)c_{r\sigma }(k_y)`$ in Eq. (7) is then related to the components of Keldysh’s Green’s functions over the close time-path. One can then calculate perturbatively the average current $`I`$ in terms of the zeroth order Green’s function. Details of this calculation can be found in Ref. cl,  and an outline is presented in the Appendix. Here an essential difference from earlier work is that previously the Green’s functions were obtained through directly solving the equation of motion, while here we shall make use of the method of image elucidated in the following section. In this way the current approach is more general and versatile, and can be easily applied to various hybrid systems.

### II.2 Generalized method of image

In our scheme for the calculation of the tunneling current, the building blocks are the zeroth order half-space Green’s functions (see the Appendix). Because in the zeroth order, lattices on the left and right sides are disconnected, the Green’s functions are defined only for each semi-infinite plane. Therefore, lattice points on the interface will have “dangling bonds”. Effectively, as shown in Figs. 1$`(b)`$ and $`(c)`$, we are imposing hard-wall boundary conditions at the end points of these dangling bonds. One thus envisages a method of image similar to that in electrostatics.

In the usual practice, the method of image is done for the continuum differential equations. It is based on the principle of superposition and the uniqueness of the solutions.jackson  When applying it to the discrete lattice, one encounters the difficulty that the image point to any source point $`𝐫`$, may not locate right at the allowed lattice points. To overcome this difficulty, we note that for each semi-infinite lattice there is discrete translational invariance along the surface direction, which we choose as the $`y`$ direction. Furthermore, since in the analysis of tunneling problems each electrode is considered to be in steady states, time-translational invariance is preserved in the individual half-space. In the subsequent sections we shall fully exploit these symmetries and thus will be concerned with the half-space Green’s function in its Fourier space representation $`g(x,x^{};k_y;\omega )`$, which effectively propagates from $`x^{}`$ to $`x`$. For each $`k_y`$ and $`\omega `$ one is therefore dealing with an effective one dimensional (1D) system (Fig. 1).

As a demonstration of the method, let us consider a 2D semi-infinite square lattice with lattice constant $`a`$ extending over the region $`x0`$ at orientation $`(hk0)`$. The hard-wall boundary condition prescribes the half-space Green’s functions to vanish over the hard walls, which consist of all points where the boundary sites can reach away from the bulk lattice (Figs. 1$`(b)`$, $`(c)`$). For general surface orientations $`(hk0)`$ and with only nearest-neighbor (n.n.) hopping one can find that the number of hard walls is given by $`\mathrm{max}\{|h|,|k|\}`$. Let us consider first the single hard-wall configurations, which includes the $`(100)`$ and the $`(110)`$ orientations (when there is no next n.n hopping in the latter). As we shall discuss later, the multi-hard-wall problem are simple generalizations to the single hard-wall cases.

For single hard-wall case the only hard wall is located at $`x=d`$, where $`d=a/\sqrt{h^2+k^2}`$ is the spacing between two consecutive $`(hk0)`$ planes. Since the Green’s function must vanish on the hard wall regardless the position of the source point $`x^{}`$, one imposes the boundary condition

$`g(d,x^{};k_y;\omega )=0.`$ (9)
To implement the method of image, we construct the half-space Green’s function $`g(x,x^{};k_y;\omega )`$ from the full-space Green’s function $`G(x,x^{};k_y;\omega )`$ as

$`g(x,x^{};k_y;\omega )=`$ $`G(x,x^{};k_y;\omega )`$ (10)
$``$ $`G(x,x_1^{};k_y;\omega )\alpha (x^{};k_y;\omega ),`$
where $`x_1^{}=2dx^{}`$ is the image point of the point source $`x^{}`$ with respect to the hard wall $`x=d`$. The Green’s function $`G(x,x^{})`$ describes direct propagation from the point source to the point $`x`$, while $`G(x,x_1^{})`$ propagates from the image point to $`x`$. The factor is determined by fitting the boundary condition (9), we find

$`\alpha (x^{})`$ $`=`$ $`G^1(d,2dx^{})G(d,x^{})`$ (11)
$`=`$ $`G^1(d+x^{})G(dx^{}).`$
Here and in the following we suppress the $`k_y`$ and $`\omega `$ dependence whenever no confusion would arise. In going from the first to the second expressions above, we have used $`G(x,x^{})=G(xx^{})`$, namely that the full-space Green’s functions are translational invariant along the $`x`$ direction. However, this is not essential for establishing the method of image. It is used here only for brevity. For systems without translational symmetry along $`x`$ direction (such as the $`d`$-density-wave state to be discussed in Sec. III.4), the following discussion still proceeds with only minor modification.

From the half-space Green’s function $`g(x,x^{})`$, one obtains the surface Green’s function $`g_0`$ by setting $`x=x^{}=0`$. In the Fourier space, $`g_0`$ can be expressed by

$`g_0={\displaystyle \underset{\pi /dk_x<\pi /d}{}}G(k_x)\times [1\mathrm{exp}(2ik_xd)\alpha _0].`$ (12)
Here $`\alpha _0=\alpha (0)`$ does not depends on $`k_x`$ and the sum over $`k_x`$ extends over the first Brillouin zone of the effective 1D lattice. The advantage of this formulation is clearly seen from (12): the surface Green’s function is obtained from combinations of full-space Green’s functions. For different surface orientations, one simply rotates the full-space Green’s function to the appropriate angle. Furthermore, it is also clear that here we have a scheme for studying the effects of interactions and fluctuations in tunneling problems. Essentially one can take these effects into account through the bulk Green’s function. Here, however, we shall concentrate on mean-field treatments and defer correlation effects to a separate publication.

It is when dealing with lattices with an anisotropic order parameter that one could most easily appreciate the power of the present formulation. For instance in dealing with $`d`$-wave superconductors, apart from fitting the boundary conditions (9), $`\alpha `$ also takes care of the different gap structures for propagation along the reflected path and the fictitious path (such as $`AO`$ and $`A^{}O`$ depicted in Fig. 2). In the presence of reflection symmetry (such as an $`s`$-wave superconductor, or a $`d`$-wave superconductor at $`(100)`$ orientation), since the gap structure as seen by these two paths are identical, the full-space bulk Green’s function possesses the symmetry $`G(d)=G(d)`$. Therefore $`\alpha `$ becomes (independent of $`k_y`$ and $`\omega `$) universally equal to the identity matrix and Eq. (12) reduces to the familiar formcl ; yan

$`g_0={\displaystyle \underset{\pi /dk_x<\pi /d}{}}G(k_x)\times 2\mathrm{sin}^2(k_xd).`$ (13)
For general orientations or when taking into account next nearest neighbor hopping, as noted earlier, there could be more than one hard walls. In these circumstances the surface Green’s function must satisfy the boundary condition that it vanishes on all these hard walls simultaneously. This is a simple generalization of the single hard-wall problem. For instance, let us consider the $`(210)`$ case with n.n hopping: as depicted in Fig. 1$`(c)`$ there are two hard walls located at $`x=d`$ and $`2d`$, where $`d=a/\sqrt{5}`$. Analogous to the single hard-wall problem, we write the half-space Green’s function

$`g(x,x^{})=G(x,x^{})G(x,x_1^{})\alpha _1(x^{})G(x,x_2^{})\alpha _2(x^{})`$ (14)
with $`x_1^{}=2dx^{}`$, $`x_2^{}=4dx^{}`$ being the location of the image sources, and $`\alpha _1`$, $`\alpha _2`$ determined by the boundary conditions

$`g(d,x^{})=0=g(2d,x^{}).`$ (15)
In other words, for the point source at $`x^{}`$, each hard wall “generates” an image source on the other side of the surface and introduces an $`\alpha `$ factor which accounts for the additional boundary conditions. The half-space Green’s function is a superposition of contributions from the real and all image sources. To obtain the surface Green’s function, one again substitutes $`x=x^{}=0`$ into Eq. (14).

Before proceeding to the applications in the following sections, we comment that the present method is not restricted to square lattices. In Sec. III.5 we will apply this method to systems involving honeycomb lattices. Indeed our generalized method of image relies only on the possibility of reducing 2D lattices into 1D structures through a Fourier transformation in the transverse direction.

## III Tunneling spectroscopy in hybrid systems

### III.1 Normal metal–$`d`$-wave Superconductors

We first study the $`ab`$-plane tunneling between a normal metal and a $`d`$-wave superconductor. The superconductor occupies the half-space $`x>0`$ and is described by the mean-field Hamiltonian

$`H_R`$ $`=`$ $`{\displaystyle \underset{ij,\sigma }{}}t_Rc_{i\sigma }^{}c_{j\sigma }{\displaystyle \underset{ij^{},\sigma }{}}t_R^{}c_{i\sigma }^{}c_{j\sigma }`$ (16)
$`+{\displaystyle \underset{ij}{}}\mathrm{\Delta }_{ij}(c_ic_jc_ic_j)+\mathrm{h}.\mathrm{c}.,`$
where $`ij`$ denotes the nearest-neighbor (n.n.) bond, $`ij^{}`$ the next nearest neighbor (n.n.n.) bond, $`t_R`$ and $`t_R^{}`$ are hopping amplitudes between n.n and n.n.n sites, respectively; $`\mathrm{\Delta }_{ij}`$ is the mean-field pairing amplitude which possesses the $`d`$-wave symmetry

$`\mathrm{\Delta }_{ij}=\{\begin{array}{cc}\mathrm{\Delta }_0& \text{for }𝐫_j=𝐫_i+𝐚\text{,}\\ \mathrm{\Delta }_0& \text{for }𝐫_j=𝐫_i+𝐛\text{.}\end{array}`$ (19)
The normal metal on the left is modeled by a Hamiltonian similar to $`H_R`$ but with only n.n hopping terms.

To obtain the corresponding 1D structure, we Fourier transform the Hamiltonian along the $`y`$ direction. For example, at (110) orientation if including only n.n. hopping $`H_R`$ becomes

$`H_R={\displaystyle \underset{x_i,k_y,\sigma }{}}2t_R\mathrm{cos}\left({\displaystyle \frac{k_ya}{\sqrt{2}}}\right)c_{i\sigma }^{}(k_y)c_{i+1\sigma }(k_y)`$
$`+{\displaystyle \underset{x_i,k_y}{}}2i\mathrm{\Delta }_0\mathrm{sin}\left({\displaystyle \frac{k_ya}{\sqrt{2}}}\right)[c_i(k_y)c_{i+1}(k_y)`$
$`+c_i(k_y)c_{i+1}(k_y)]+\mathrm{h}.\mathrm{c}.`$ (20)
Here $`c_i`$ are the electron annihilation operators for the 1D lattice at the $`i`$-th site (see Fig. 1$`(b)`$) and $`a`$ is the lattice constant of the original 2D lattice. The lattice constant of the 1D lattice is identical to the distance $`d`$ between two consecutive $`(hk0)`$ planes (cf. Fig. 1). This 1D structure of the problem is very helpful for us since each of the lattice planes $`(hk0)`$ now becomes a point on the $`x`$ axis. This enables us to define for each lattice site the corresponding image sites with respect to the hard walls.wm

It is convenient to use the Nambu notation which distinguishes particle and hole components. We define the spinor field-operatornote\_spinor

$`\mathrm{\Psi }_i(k_y,t)=\left(\begin{array}{c}c_i(k_y,t)\\ c_i^{}(k_y,t)\end{array}\right).`$ (23)
The upper and the lower components of $`\mathrm{\Psi }_i`$ correspond to the particle and hole components, respectively. For (110) orientation one can then write $`H_R`$ of (20) in the form

$`H_R`$ $`=`$ $`{\displaystyle \underset{x_i,k_y}{}}\left(\mathrm{\Psi }_i^{}H_{i,i+1}\mathrm{\Psi }_{i+1}+\mathrm{\Psi }_i^{}H_{i,i1}\mathrm{\Psi }_{i1}\right)`$ (24)
with

$`H_{i,i\pm 1}=\left(\begin{array}{cc}2t_R\mathrm{cos}(\frac{k_ya}{\sqrt{2}})& \pm 2i\mathrm{\Delta }_0\mathrm{sin}(\frac{k_ya}{\sqrt{2}})\\ \pm 2i\mathrm{\Delta }_0\mathrm{sin}(\frac{k_ya}{\sqrt{2}})& 2t_R\mathrm{cos}(\frac{k_ya}{\sqrt{2}})\end{array}\right).`$ (27)
In order to apply the Keldysh formulation to calculating the tunneling current, as detailed in the Appendix, the basic quantity one shall need is the (bare) half-space retarded Green’s function. According to our method of image this can be obtained from superposition of the full-space Green’s functions. Therefore our remaining task is to find the full-space Green’s function. For general interface orientations $`(hk0)`$ all that we need is to rotate the full-space Green’s function to the appropriate angle and then build up the half-space Green’s function based on the recipe outlined in Sec. II.

To find the full-space Green’s function, we go over to the momentum space and express the Hamiltonian (16) in Nambu’s representation

$`H_R`$ $`=`$ $`{\displaystyle \underset{𝐤,\sigma }{}}ϵ_𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }+{\displaystyle \underset{𝐤}{}}(\mathrm{\Delta }_𝐤c_𝐤^{}c_𝐤^{}+\mathrm{\Delta }_𝐤^{}c_𝐤c_𝐤)`$ (33)
$`=`$ $`{\displaystyle \underset{𝐤}{}}(\begin{array}{cc}c_𝐤^{}& c_𝐤\end{array})(\begin{array}{cc}ϵ_𝐤& \mathrm{\Delta }_𝐤\\ \mathrm{\Delta }_𝐤^{}& ϵ_𝐤\end{array})\left(\begin{array}{c}c_𝐤\\ c_𝐤^{}\end{array}\right).`$
Here the quasiparticle dispersion $`ϵ_𝐤`$ and the gap function $`\mathrm{\Delta }_𝐤`$ are given by

$`ϵ_𝐤`$ $`=`$ $`2t_R[\mathrm{cos}(𝐤𝐚)+\mathrm{cos}(𝐤𝐛)]`$
$`4t_R^{}\mathrm{cos}(𝐤𝐚)\mathrm{cos}(𝐤𝐛),`$
$`\mathrm{\Delta }_𝐤`$ $`=`$ $`2\mathrm{\Delta }_0[\mathrm{cos}(𝐤𝐚)\mathrm{cos}(𝐤𝐛)].`$ (34)
For $`(hk0)`$ orientation the lattice vectors $`𝐚=a(\mathrm{cos}\theta ,\mathrm{sin}\theta )`$, $`𝐛=a(\mathrm{sin}\theta ,\mathrm{cos}\theta )`$, where $`a`$ is the lattice constant and $`\theta `$ is the angle between the $`a`$-axis and the $`x`$-direction (thus $`\mathrm{tan}\theta =k/h`$). The Hamiltonian $`H_R`$ in the form (33) is readily diagonalized, the quasiparticle excitation energy is found to be $`\pm E_𝐤=\pm \sqrt{ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2}`$.

The full-space retarded Green’s function can be obtained from

$`G(x_i,x_j)={\displaystyle \underset{\pi /dk_x<\pi /d}{}}G(k_x,k_y,\omega )\times \mathrm{e}^{ik_x(x_ix_j)},`$ (35)
where $`G(k_x,k_y,\omega )=[\omega +i\eta \widehat{H}_R(k_x,k_y,\omega )]^1`$, with $`\widehat{H}_R`$ the matrix in the second line of (33) and $`\eta `$ an infinitesimal positive number. The half-space bare Green’s function $`g_0^r`$ is then obtained from the method of image. In the tunneling problem, the tunneling Hamiltonian brings in tunneling events between the two sides of the tunnel junctions which “renormalize” the half-space Green’s functions (see the Appendix). In the Keldysh formulation this is expressed as a perturbation series which can be re-summed to all orders in the tunneling amplitude $`t`$, yielding the renormalized half-space Green’s functions.caroli  With the assumption that the renormalized advanced and retarded half-space Green’s functions satisfy

$`[g_{\alpha \beta }^r]^{}=g_{\beta \alpha }^a`$ (36)
($`\alpha ,\beta =\{L,R\}`$ labels the electrodes), we can express the tunneling current as

$`I=I_1+I_2+I_3+I_A,`$ (37)
where

$`I_1`$ $`=`$ $`{\displaystyle \underset{k_y}{}}4\pi e{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\omega t^2[f(\omega eV)f(\omega )]A_{L,11}(\omega eV)A_{R,11}(\omega )|1+tg_{RL,11}^r(\omega )|^2,`$ (38)
$`I_2`$ $`=`$ $`{\displaystyle \underset{k_y}{}}8\pi e{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\omega t^2[f(\omega eV)f(\omega )]A_{L,11}(\omega eV)Re\{A_{R,12}(\omega )[tg_{LR,21}^a(\omega )(1+tg_{RL,11}^r(\omega ))]\},`$ (39)
$`I_3`$ $`=`$ $`{\displaystyle \underset{k_y}{}}4\pi e{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\omega t^4[f(\omega eV)f(\omega )]A_{L,11}(\omega eV)A_{R,22}(\omega )|g_{RL,12}^r(\omega )|^2,`$ (40)
$`I_A`$ $`=`$ $`{\displaystyle \underset{k_y}{}}4\pi e{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\omega t^4[f(\omega eV)f(\omega +eV)]A_{L,11}(\omega eV)A_{L,22}(\omega +eV)|g_{RR,12}^r(\omega )|^2.`$ (41)
Here the range of the $`k_y`$ is given by (8),

$`A_\alpha =i/(2\pi )(g_{0,\alpha \alpha }^rg_{0,\alpha \alpha }^{r}{}_{}{}^{})`$ (42)
are the spectral weight matrices for the electrode $`\alpha =\{L,R\}`$, and $`f(\omega )`$ is the Fermi function (at zero temperature it is simply the step function $`\mathrm{\Theta }(\omega )`$). The indices 1, 2 in the Green’s functions and the spectral weight matrices refer respectively to the particle and the hole components in the Nambu representation. $`t=t(\omega ,k_y)`$ is the tunneling amplitude between the two electrodes. It is remarkable that the expression for $`I_2`$ here generalizes that found in Ref. cl,  and is applicable to any interface orientation. For the special cases considered in Ref. cl, , where the surface Green’s functions are symmetric (for $`(100)`$ orientation) or antisymmetric (for $`(110)`$ orientation), Eq. (39) reproduces previous results. From these formulas one can clearly identify the contributions from each channel in the tunneling process. In particular, $`I_1`$ is the contribution from single particle tunneling and $`I_A`$ the Andreev reflection (thus $`I_A`$ depends on the particle and hole components of the spectral weight matrix $`A_L`$).

We now present some of our results. Fig. 3 shows the tunneling spectra for (110) and (210) orientations at the doping levels $`\delta =0.08`$, 0.14, and 0.20. Here we study the doping dependence by resorting to the mean-field slave boson theory for the $`t`$-$`t^{}`$-$`J`$ model. The electron operators $`c`$ and $`c^{}`$ are then essentially the spinon operators and the Green’s function for spinons as well. The holons condense so that $`b=\sqrt{\delta }`$. The mean-field parameters $`t_R`$, $`t_R^{}`$, $`\mathrm{\Delta }_0`$, and the chemical potential $`\mu _R`$ for each doping are calculated self-consistently.cl  It is obvious from Fig. 3 that the ZBCP is significantly reduced in the (210) orientation. Interestingly, for (110) orientation the ZBCP decreases upon increasing doping while for (210) case it grows and then falls with doping. Another interesting feature in the tunneling spectra is the subgap structures near $`\pm 2\mathrm{\Delta }_0`$ in the (210) case. These may have originated from resonances due to broken surface pairs, resulting from the dangling bonds in (210) orientations.barash

The ZBCP originates from zero-energy surface states (or the midgap states) due to Andreev reflections. In our formulation these states arise from singularities in the image contributions which manifest as poles in the $`\alpha `$ factors. In the presence of a single hard wall, the poles are determined by the zeros of the following factor when $`\eta =0`$

$`\beta (k_y)=det[G(d;k_y,\omega =0)].`$ (43)
This produces singular behavior in the Green’s functions and results in the ZBCP. In the (100) case, since $`\alpha _0`$ is simply the identity matrix the surface Green’s function (13) is regular at $`\omega =0`$ thus there is no ZBCP.

### III.2 Tunneling into current-carrying superconductors

We now consider the $`ab`$-plane tunneling from normal metals to current-carrying superconductors. In experiments one applies magnetic field along the $`c`$-axis of the superconductor, so that a screening current is generated over the $`ab`$ plane. When a quasiparticle tunnels across the surface layer, it acquires additional energy from the supercurrent. Thus the zero-energy surface state evolves in this case into two surface states with non-zero energy. In the tunneling spectra this appears as “splitting” of the ZBCP (Fig. 4). Fogelström et al. have analyzed splittings of the zero-energy peak in the surface density of states under applied fields in the continuum limit.fogels  Here we examine the tunneling spectra base on our discrete model.

To marry with formulations in the previous section, we note that in the presence of supercurrent the gap function is modified asdeGennes

$`\mathrm{\Delta }_{ij}\mathrm{\Delta }_{ij}\mathrm{exp}(i𝐪(𝐫_i+𝐫_j)),`$ (44)
where $`𝐪=(0,q_y)`$ is the superfluid momentum and is proportional to the magnetic field. We shall assume that tunneling events take place only within a shallow layer of order about the penetration depth from the surface, so that $`q_y`$ is approximately uniform in the region of our concern. This additional phase can be absorbed into the electron operator by the transformation $`c_{i\sigma }c_{i\sigma }\mathrm{exp}(i𝐪𝐫_i)`$. In Fourier space the Hamiltonian becomes

$`H_R`$ $`=`$ $`{\displaystyle \underset{𝐤,\sigma }{}}ϵ_{𝐤+𝐪}c_{𝐤,\sigma }^{}c_{𝐤,\sigma }+{\displaystyle \underset{𝐤}{}}(\mathrm{\Delta }_𝐤c_𝐤^{}c_𝐤^{}+\mathrm{\Delta }_𝐤^{}c_𝐤c_𝐤)`$ (50)
$`=`$ $`{\displaystyle \underset{𝐤}{}}(\begin{array}{cc}c_𝐤^{}& c_𝐤\end{array})(\begin{array}{cc}ϵ_{𝐤+𝐪}& \mathrm{\Delta }_𝐤\\ \mathrm{\Delta }_𝐤^{}& ϵ_{𝐤𝐪}\end{array})\left(\begin{array}{c}c_𝐤\\ c_𝐤^{}\end{array}\right).`$
Here the quasiparticle dispersion $`ϵ_𝐤`$ and the gap function $`\mathrm{\Delta }_𝐤`$ are those of Eq. (34). After diagonalizing $`H_R`$ above one finds the quasiparticle excitation energy becomes

$`E_𝐤^{(\pm )}=({\displaystyle \frac{ϵ_{𝐤+𝐪}ϵ_{𝐤𝐪}}{2}})\pm \sqrt{\left({\displaystyle \frac{ϵ_{𝐤+𝐪}+ϵ_{𝐤𝐪}}{2}}\right)^2+\mathrm{\Delta }_𝐤^2}.`$ (51)
The momentum space Green’s function that is fed into Eq. (12) is obtained in the same way: $`G(k_x,k_y,\omega )=[\omega +i\eta \widehat{H}_R(k_x,k_y,\omega )]^1`$, with $`\widehat{H}_R`$ the matrix in the second line of (50).

Fig. 4 shows typical tunneling spectra for the splitting of the ZBCP when increasing $`q_y`$. Note that the slightly asymmetric splitting originates from the particle-hole asymmetry in $`ϵ_𝐤`$. Fig. 5 plots the magnitude of the splitting versus the applied magnetic field for underdoped case. For small $`𝐪`$ , the expansion of Eq (51) leads to linear splitting in the lowest order terms: $`E_𝐤^{(\pm )}=\pm E_𝐤+𝐪\frac{ϵ_𝐤}{𝐤}`$, where $`E_𝐤=\sqrt{ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2}`$, as observed in small applied fields. For higher fields, one has to retain the full $`𝐪`$ dependence, resulting in the bending of the splitting. This is purely due to the lattice effect. Also shown in Fig. 5 are the results taking into account suppression of the superconducting gap under magnetic fields self-consistently. The curve is seen to be “pushed” inwards while maintaining similar features. Note that quantitative agreement with experimental observationsaprili ; deutsch  can be obtained by fitting scales of our results to the experimental data. Nevertheless, we did not observe any zero-field splitting at overdoping. This is in contrast with the experiment of Ref. deutsch, . The mechanism inducing this splitting might have eluded from our simple model.

The doping dependence of splitting is also shown in Fig. 6 for $`q_y=0.15`$ and $`0.80`$, which are respectively in the linear and the saturated regimes in Fig. 5. Note that the splitting increases with doping, in agreement with Ref. deutsch, .

In passing we point out that the splitting depends sensitively on the Fermi surface topology. Indeed for $`\mu _R=0`$ we find no splitting of the ZBCP whatever the value of $`q_y`$ is. One can confirm this analytically by making an asymptotic expansion of the Green’s function around $`\omega =0`$. At $`\mu _R=0`$ one finds the conductance peak invariantly stays at the zero bias.

### III.3 Polyacetylene

Up to this point, we have considered tunnel junctions with superconducting test electrodes, where particles of opposite spins form pairs. In this and the following sections we will consider systems which exhibits particle-hole pairing over bipartite lattices. To start with we shall consider first a simple 1D model based on the structure of polyacetylene.su  This will turn out to be very helpful for understanding results in the following sections. Most importantly, it provides the criteria for the formation of midgap states in semi-infinite bipartite systems.

The model we shall examine here is a 1D chain with alternating hopping amplitudes $`t_1`$, $`t_2`$ as shown in Fig. 7. The separation between the lattice points is taken to be a constant $`a`$.note\_poly  It is convenient to categorize the lattice points into $`A`$ and $`B`$ sublattices and express the Hamiltonian for this “$`t_1`$-$`t_2`$ model” as

$`H_R={\displaystyle \underset{i_B,\sigma }{}}t_1c_{i1}^Ac_i^Bt_2c_i^Bc_{i+1}^A+\mathrm{h}.\mathrm{c}.`$ (52)
Here $`c_i^\alpha `$ annihilates electrons over site $`i`$ on the $`\alpha =\{A,B\}`$ sublattice (spin indices $`\sigma `$ will be omitted throughout), and the sum run over sites $`i`$ in the $`B`$ sublattice only. Going over to momentum space, one finds

$`H_R`$ $`=`$ $`{\displaystyle \underset{k,\sigma }{}}\mathrm{\Lambda }_kc_{k}^{A}{}_{}{}^{}c_k^B+\mathrm{h}.\mathrm{c}.`$ (58)
$`=`$ $`{\displaystyle \underset{k,\sigma }{}}(\begin{array}{cc}c_{k}^{A}{}_{}{}^{}& c_{k}^{B}{}_{}{}^{}\end{array})(\begin{array}{cc}0& \mathrm{\Lambda }_k\\ \mathrm{\Lambda }_k^{}& 0\end{array})\left(\begin{array}{c}c_k^A\\ c_k^B\end{array}\right),`$
where

$`\mathrm{\Lambda }_k`$ $`=`$ $`t_1\mathrm{exp}(ika)t_2\mathrm{exp}(ika)`$ (59)
$`=`$ $`(t_1+t_2)\mathrm{cos}(ka)i(t_1t_2)\mathrm{sin}(ka).`$
The Hamiltonian (58) can be diagonalized easily and the quasiparticle excitation energy is found to be $`\pm |\mathrm{\Lambda }_k|`$. Note that the real part of $`\mathrm{\Lambda }_k`$ is similar to the usual hopping energy in one dimensional chain. Therefore, when $`Im\{\mathrm{\Lambda }_k\}(t_1t_2)0`$ a single particle excitation gap opens at the chemical potential.

For semi-infinite chain, there are two possible configurations with the terminating site being an $`A`$ or a $`B`$ sublattice point. In either case we choose the boundary point the origin $`x=0`$ and construct the surface Green’s function utilizing the method of image

$`g_0=G(0,0)G(0,2a)G^1(a,2a)G(a,0).`$ (60)
In this formula, the appropriate Green’s functions should be used depending the the type of the end point; for instance, in the case of an $`A`$-type boundary even/odd sites are attributed to the $`A`$/$`B`$ sublattices. The retarded Green’s function for infinite system is

$`G(x_i^\alpha ,t;x_j^\beta ,0)=i\mathrm{\Theta }(t)\{c_i^\alpha (t),c_{j}^{\beta }{}_{}{}^{}(0)\}`$ (61)
for $`x_i^\alpha `$ over the $`\alpha `$ sublattices; the braces denote the anti-commutators. After Fourier transformation in time, the Green’s function is obtained from

$`G(x_i^\alpha ,x_j^\beta ;\omega )={\displaystyle \underset{\pi /2ak_x<\pi /2a}{}}G_{\alpha \beta }(k,\omega )\times \mathrm{e}^{ik(x_i^\alpha x_j^\beta )}.`$ (62)
Utilizing $`\widehat{H}_R`$ the matrix in the second line of (58), we write for brevity the momentum space Green’s functions in matrix form

$`G_{\alpha \beta }(k,\omega )`$ $`=`$ $`[\omega +i\eta \widehat{H}_R(k,\omega )]_{\alpha \beta }^1`$ (65)
$`=`$ $`{\displaystyle \frac{1}{(\omega +i\eta )^2E_k^2}}\left(\begin{array}{cc}\omega +i\eta & \mathrm{\Lambda }_k\\ \mathrm{\Lambda }_k^{}& \omega +i\eta \end{array}\right)_{\alpha \beta }.`$
The matrix elements are assigned to the Green’s functions according to the convention in (58), namely $`\alpha =A,B`$ corresponds to $`\alpha =1,2`$ respectively.

Since in Eq. (60) the surface Green’s function $`g_0`$ is expressed as a combination of bulk Green’s functions, the only possible source of singular behavior in $`g_0`$ resides in the inverse part $`G^1(a,2a)`$. In other words, the existence of the zero-energy mode depends on the behavior of $`G(a,2a)`$ at $`\omega =0`$. This is analogous to the ND junctions where the ZBCP results from the zeros of the determinant $`\beta (k_y)`$, Eq. (43).

For example, in the case of $`A`$-type boundaries we find explicitly (setting $`\eta =0`$)

$`G_{BA}(a,2a;\omega =0)`$ $`=`$ $`{\displaystyle \frac{a}{2\pi }}{\displaystyle _{\pi /2a}^{\pi /2a}}dk{\displaystyle \frac{1}{t_1+t_2\mathrm{exp}(2ika)}}`$ (68)
$`=`$ $`\{\begin{array}{cc}0& \text{if }t_1<t_2\text{,}\\ 1/2t_1& \text{if }t_2<t_1\text{.}\end{array}`$
Here the index “$`BA`$” denotes the same meaning as in (65) and is used for emphasizing the correct Green’s function to be used. From (68), when $`t_1<t_2`$ a sharp singularity in the surface Green’s function $`g_0`$ arises at $`\omega =0`$ due to the divergent factor $`G_{BA}^1`$ in (60). On the other hand, for $`t_1>t_2`$ since $`G_{BA}^1`$ is finite at $`\omega =0`$, no singular behavior in $`g_0`$ could occur there. Thus one expects ZBCP in the former case while none in the latter. In the following sections we shall see that this provides for 2D bipartite systems a criterion for the range of transverse momenta where zero-energy states exist. For $`B`$-type boundary the analysis is identical, except an exchange in the roles of $`t_1`$ and $`t_2`$. Therefore when the ZBCP shows up in an $`A`$-type chain, it must be absent in a $`B`$-type chain, and vice versa. This is shown in Fig. 8 for the case of polyacetylene. The current expression here is identical to Eq. (79) given in the following section, except the extra sum over $`k_y`$ there.

### III.4 Normal metal–$`d`$-density wave states

In underdoped cuprate superconductors, it is observed in experiments that there are signatures of a “partial” gap well above the superconducting temperature $`T_c`$. This anomalous regime in the phase digram of the cuprate superconductors is thus termed the pseudogap phase.timusk  Experiments also find that the pseudogap is consistent with a $`d`$-wave structure. Recently Chakravarty et al. proposed that the pseudogap phase of the underdoped cuprate is possibly the $`d`$-density-wave (DDW) state.chakravarty  It is therefore of interest to examine the tunneling spectra of normal-metal–$`d`$-density-wave (N-DDW) junctions.

The DDW state is characterized by the staggered flux in the elementary plaquettes of the lattice. The bond currents circulating the unit cell of the underlying square lattice break, among other symmetries, the invariance of translation by one lattice spacing and lead to a bipartite structure (Fig. 9). Obviously, if the interface cuts at (110) direction, the reflection symmetry is broken – in contrast to the (100) case. Therefore, we shall examine the (110) direction with the following mean-field Hamiltonianaffleck

$`H_R={\displaystyle \underset{i_B,\sigma }{}}\{\chi (c_{i+𝐚}^Ac_i^B+c_{i𝐚}^Ac_i^B)`$
$`+\chi ^{}(c_{i+𝐛}^A^{}c_i^B+c_{i𝐛}^Ac_i^B)+\mathrm{h}.\mathrm{c}.\},`$ (69)
where $`c_i^\alpha `$ annihilates an electron at site $`i`$ over the $`\alpha `$ sublattice, and $`\chi `$ is the hopping amplitude on the bond (Fig. 9$`(b)`$). Making Fourier transformation along the $`y`$ direction in $`H_R`$, one finds

$`H_R={\displaystyle \underset{x_i^B,k_y,\sigma }{}}\{`$ $`\mathrm{\Lambda }_{i,i1}c_{i1}^A(k_y)c_i^B(k_y)`$ (70)
$`+`$ $`\mathrm{\Lambda }_{i,i+1}c_i^B(k_y)c_{i+1}^A(k_y)+\mathrm{h}.\mathrm{c}.\},`$
where $`\mathrm{\Lambda }_{i,i\pm 1}=2Re\{\chi \mathrm{e}^{\pm ik_ya/\sqrt{2}}\}`$ with $`a`$ the lattice constant of the square lattice; $`c_i^\alpha (k_y)`$ is the electron annihilation operator for fixed $`k_y`$ at site $`x_i^\alpha `$ over the 1D lattice. Going over to momentum space, one finds, similar to (58)

$`H_R={\displaystyle \underset{𝐤,\sigma }{}}(\begin{array}{cc}c_{𝐤}^{A}{}_{}{}^{}& c_{𝐤}^{B}{}_{}{}^{}\end{array})(\begin{array}{cc}0& \mathrm{\Lambda }_𝐤\\ \mathrm{\Lambda }_𝐤^{}& 0\end{array})\left(\begin{array}{c}c_𝐤^A\\ c_𝐤^B\end{array}\right)`$ (76)
with $`\mathrm{\Lambda }_𝐤=ϵ_𝐤+i\mathrm{\Delta }_𝐤`$. $`ϵ_𝐤`$ and $`\mathrm{\Delta }_𝐤`$ are given by Eq. (34) with $`t_R=Re\{\chi \}`$, $`\mathrm{\Delta }_0=Im\{\chi \}`$, and $`t_R^{}=0`$. The quasiparticle excitation energies is then obviously $`\pm E_𝐤=\pm |\mathrm{\Lambda }_𝐤|=\pm \sqrt{ϵ_𝐤^2+\mathrm{\Delta }_𝐤^2}`$.

To find the tunneling current we apply again the Keldysh formulation outlined in the Appendix. For fixed $`k_y`$ the full-space retarded Green’s function between the sites $`x_i^\alpha `$ and $`x_j^\beta `$ pertaining to the $`\alpha `$ and $`\beta `$ sublattices is given by

$`G(x_i^\alpha ,t;x_j^\beta ,0;k_y)=i\mathrm{\Theta }(t)\{c_i^\alpha (k_y,t),c_{j}^{\beta }{}_{}{}^{}(k_y,0)\},`$ (77)
which in Fourier space for fixed frequency $`\omega `$ becomes

$`G(x_i^\alpha ,x_j^\beta )={\displaystyle \underset{\pi /2dk_x<\pi /2d}{}}G_{\alpha \beta }(k_x,k_y,\omega )\times \mathrm{e}^{ik_x(x_i^\alpha x_j^\beta )}.`$ (78)
Here $`d=a/\sqrt{2}`$ is the lattice spacing of the 1D lattice and we have suppressed the $`\omega `$ and $`k_y`$ dependences on the left hand side. $`G_{\alpha \beta }`$ has the same form given in Eq. (65) except that now $`\alpha =A`$ or $`B`$.

From $`G_{\alpha \beta }`$, the half-space surface Green’s function is obtained again using the method of image. The current expressions here, however, are distinct from those of (38)–(41). Indeed since we are dealing with a single component Green’s function the calculation is much simpler than previously. As shown in the Appendix, the current expression is here

$`I={\displaystyle \underset{k_y,\sigma }{}}2\pi e{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\omega t^2[f(\omega eV)f(\omega )]`$
$`\times A_L(\omega eV)A_R(\omega )|1+tg_{RL}^r(\omega )|^2.`$ (79)
This is exactly the single-particle current $`I_1`$ of Eq. (38) for ND tunneling. There is no contribution from “Andreev reflections” in N-DDW tunneling. This is due to the fact that in the DDW state the pairing takes place between particles and holes of momenta $`𝐤`$ and $`𝐤+𝐐`$, with $`𝐐`$ the nesting vector of 2D square lattices. Thus the Andreev reflected particles are still electrons whose response to the bias voltage are the same as the incident particles; as a result their contributions to the tunneling current cancel exactly. In the ND junction, however, a particle is Andreev reflected as a hole, which behaves oppositely under applied bias. Fig. 10$`(a)`$ shows a typical plot for differential conductance versus voltage for N-DDW junctions. The conspicuous ZBCP agrees with recent calculations done by Honerkamp and Sigrist.honerkamp

The reason for the ZBCP here can be understood on the basis of the results in the previous section. Just like polyacetylene, the midgap states arises when $`g_0`$ is singular due to the zeros in the Green’s function such as in Eq. (68). For each $`k_y`$ Eq. (70) resembles the $`t_1`$-$`t_2`$ model with $`t_1=\mathrm{\Lambda }_{i,i1}`$ and $`t_2=\mathrm{\Lambda }_{i,i+1}`$. Therefore, for example, for $`A`$-type boundary one expects midgap states for the range of $`k_y`$ where

$`\mathrm{\Lambda }_{i,i1}>\mathrm{\Lambda }_{i,i+1}\mathrm{or}Im\{\chi \}\mathrm{sin}\left({\displaystyle \frac{k_ya}{\sqrt{2}}}\right)>0.`$ (80)
Since here $`Im\{\chi \}=\mathrm{\Delta }_0<0`$, the above equation leads to $`\sqrt{2}\pi /a<k_y<0`$.

We have so far considered only the case of vanishing chemical potential $`\mu _R`$ in the DDW state. At finite chemical potential the grand Hamiltonian for the DDW is $`K_R=H_R\mu _RN_R`$. Since the number operator

$`N_R={\displaystyle \underset{𝐤\sigma }{}}(\begin{array}{cc}c_{𝐤}^{A}{}_{}{}^{}& c_{𝐤}^{B}{}_{}{}^{}\end{array})(\begin{array}{cc}1& 0\\ 0& 1\end{array})\left(\begin{array}{c}c_𝐤^A\\ c_𝐤^B\end{array}\right).`$ (86)
Hence $`\mu _RN_R`$ is diagonal and it simply shifts the excitation energy $`E_𝐤`$ to $`E_𝐤\mu _R`$. It is easy to check that this change in the excitation energy induces a shift $`\omega \omega +\mu _R`$ in the Green’s function. This is in sharp contrast with the ND case; there the chemical potential shifts the quasiparticle energy $`ϵ_𝐤ϵ_𝐤\mu _R`$ in the Green’s function but not the frequency. This results in the distinct behavior of the ZBCP for ND and N-DDW junctions at finite $`\mu _R`$ .

For N-DDW junctions since $`\omega \omega +\mu _R`$ at finite chemical potential, the conductance peak is shifted from zero bias to the opposite value of the chemical potential $`\mu _R`$. For ND junctions, however, the midgap state stays at $`\omega =0`$ even at finite chemical potential, thus the conductance peak always position at zero bias (see Fig. 3). This shift has an obvious implication: the peak will split due to the Zeeman splitting (see Fig. 10). The orbital effects of magnetic fields can be included by changing $`\chi `$ into $`\chi \mathrm{e}^{i𝐪(𝐫_i𝐫_j)}`$ for any nearest neighbor sites $`i`$, $`j`$. This takes into account the current induced near the interface. Since under this change both $`ϵ_𝐤`$ and $`\mathrm{\Delta }_𝐤`$ undergo shifting of $`𝐤`$ by $`𝐪`$ which can be absorbed into the summation of $`𝐤`$, the peak does not split. Therefore, the splitting of ZBCP turn out the same for both in-plane and perpendicular magnetic fields. This is in contrast to the ND junction where orbital effects dominate for perpendicular fields.

In closing this section we note that since the next n.n term $`4t_R^{}\mathrm{cos}(𝐤𝐚)\mathrm{cos}(𝐤𝐛)`$ couples only lattice sites within each sublattice, it introduces only diagonal terms to (76). Therefore, similar to the chemical potential, the next n.n terms cause the ZBCP and the spectrum to migrate when $`t_R^{}0`$. This is displayed in Fig. 11.

### III.5 Graphite sheets

So far we have considered systems involving only square lattices. As commented in the end of Sec. II, our formulation is indeed quite general and can be applied to any systems which can be projected into 1D structures. As an example, we study in this section the in-plane tunneling from a normal metal into semi-infinite graphite sheets (NG junctions).

In the 2D graphite sheets, each carbon atom is bonded to its three n.n via three electrons in the $`sp^2`$ orbital, forming a honeycomb lattice. The fourth electron (the $`\pi `$ electron), which resides in the $`p_z`$ orbital perpendicular to the 2D layer, is not bonded and is free to hop from site to site. In the tight-binding limit the Hamiltonian for the bulk graphite sheet is thus

$`H_R={\displaystyle \underset{i_B,\sigma }{}}\gamma _0c_{i+𝐚}^Ac_i^B\gamma _1c_{i+𝐛}^Ac_i^B\gamma _2c_{i𝐚𝐛}^Ac_i^B+\mathrm{h}.\mathrm{c}.`$ (87)
Here the lattice is divided into $`A`$ and $`B`$ sublattices, and $`𝐚`$, $`𝐛`$ are the lattice vectors illustrated in Fig. 12. $`c_i^\alpha `$ are electron annihilation operators at site $`i`$ over the $`\alpha `$ sublattice, and $`\gamma _i`$ are the hopping integrals. For simplicity we shall take $`\gamma _0=\gamma _1=\gamma _2`$ in the following. We will be interested in two orientations of the lattice: one with zigzag and the other with armchair boundaries.

We first consider the zigzag case and choose the frame of coordinates as shown in Fig. 12$`(a)`$. Fourier transformation in the transverse direction leads to 1D Hamiltonian which resembles (52)

$`H_R={\displaystyle \underset{x_i^B,k_y,\sigma }{}}t_1c_{i1}^A(k_y)c_i^B(k_y)t_2c_i^B(k_y)c_{i+1}^A(k_y)+\mathrm{h}.\mathrm{c}.`$
with here

$`t_1=2\gamma _0\mathrm{cos}({\displaystyle \frac{\sqrt{3}}{2}}k_ya)\mathrm{and}t_2=\gamma _0.`$ (89)
Further $`k_x`$ transformation brings $`H_R`$ into the same form as (76) with

$`\mathrm{\Lambda }_𝐤=\gamma _0[\mathrm{e}^{ik_xa}+2\mathrm{cos}({\displaystyle \frac{\sqrt{3}}{2}}k_ya)\mathrm{e}^{ik_xa/2}].`$ (90)
In applying the method of image, we note that the projected 1D lattice for the zigzag case has alternating bond length, which breaks the reflection symmetry and hence implies the possible existence of the ZBCP. The alternating bond length, however, seems to cause difficulty in locating the image point of an arbitrary source site. For instance, the usual choice – the mirror image – does not always put the image point right on the lattice. Nevertheless, since in 1D the hard wall becomes a point, as long as the Green’s function propagating from the real source to the hard wall can be canceled by that from a fictitious source so that the boundary condition is satisfied, uniqueness of the half-space Green’s function implies that the location of the fictitious source can be chosen at will. Indeed, this can be explicitly checked numerically. To be definite, we shall place the fictitious source at $`x=(3/2)a`$ and apply the method of image. The boundary condition $`g(a,x^{})=0`$ for all $`x^{}`$ immediately leads to

$`g_0`$ $`=`$ $`G_{AA}(0,0)`$ (91)
$``$ $`G_{AA}(0,3/2a)G_{BA}^1(a,3/2a)G_{BA}(a,0).`$
Here we have labelled the attributes of the lattice points explicitly in the subscripts of the Green’s functions. Just like polyacetylene, the midgap states arises when $`g_0`$ is singular, namely at the zeros of $`G_{BA}(a,3/2a)`$ when $`\eta =0`$. From (89) the correspondence to the $`t_1`$-$`t_2`$ model indicates that midgap states exist for $`k_y`$ which satisfy $`\mathrm{cos}(k_y\sqrt{3}a/2)<1/2`$, ie. when setting $`\sqrt{3}a=1`$

$`\pi k_y<2\pi /3\mathrm{and}2\pi /3<k_y\pi .`$ (92)
This is exactly what is found in band structure calculations.fujita

For the zigzag orientation, apart from the zigzag boundary, there could also be the “bearded” boundary where the surface layer consists of $`B`$ sites. This is reminiscent of the case of $`B`$-type end point of the $`t_1`$-$`t_2`$ model. Similar analysis as above can also carry over here with minor change. We find in this case the zero energy state arises from the range ($`\sqrt{3}a=1`$)

$`2\pi /3<k_y<2\pi /3.`$ (93)
The current expression are the same as Eq. (79) for N-DDW junctions. The corresponding tunneling spectra are solid and dashed lines shown in Fig. 13.

Let us now consider the armchair case (Fig. 12$`(b)`$). After the Fourier transformation along the interface, one finds

$`H_R={\displaystyle \underset{x_i^B,k_y,\sigma }{}}\gamma _0`$ $`(`$ $`c_{i}^{B}{}_{}{}^{}(k_y)c_i^A(k_y)\mathrm{e}^{ik_ya}`$
$`+`$ $`c_{i}^{B}{}_{}{}^{}(k_y)c_{i+1}^A(k_y)\mathrm{e}^{ik_ya/2}`$
$`+`$ $`c_{i}^{B}{}_{}{}^{}(k_y)c_{i1}^A(k_y)\mathrm{e}^{ik_ya/2}+\mathrm{h}.\mathrm{c}.)`$
Note that here for each site $`x_i`$ there are both $`A`$ and $`B`$ components as shown in Fig. 12$`(b)`$. The propagation between $`x_i`$ and $`x_j`$ thus compose of four components and the full-space Green’s function form a $`2\times 2`$ matrix

$`G(x_i,x_j)=\left(\begin{array}{cc}G(x_i^A,x_j^A)& G(x_i^A,x_j^B)\\ G(x_i^B,x_j^A)& G(x_i^B,x_j^B)\end{array}\right).`$ (97)
Further $`k_x`$ transformation yields the Hamiltonian (76) with here

$`\mathrm{\Lambda }_𝐤=\gamma _0[\mathrm{e}^{ik_ya}+2\mathrm{cos}({\displaystyle \frac{\sqrt{3}}{2}}k_xa)\mathrm{e}^{ik_y\frac{a}{2}}].`$ (98)
Note that $`\mathrm{\Lambda }(k_x,k_y)=\mathrm{\Lambda }(k_x,k_y)`$ implies that the reflection symmetry is not broken. Similar to the case of polyacetylene, in momentum space, $`G(𝐤,\omega )`$ has exactly the same form as Eq. (65). However, now the hard-wall boundary condition becomes a matrix equation

$`g(x_i^\alpha ,x_j^\beta )|_{x_i=d}=0\text{for all }\alpha ,\beta =\{A,B\},`$ (99)
where $`d=(\sqrt{3}/2)a`$ is the lattice constant of the projected 1D lattice. The surface Green’s function then takes the form

$`g_0=G(0,0)G(0,2d)G^1(d,2d)G(d,0).`$ (100)
Since translational symmetry is preserved in this 1D lattice, we have $`G(d,2d)=G(d)`$ and $`G(d,0)=G(d)`$. Reflection symmetry implies $`G(d)=G(d)`$ and consequently

$`g_0=G(0)G(2d).`$ (101)
In this case, one thus expects no midgap states.

Without loss of generality, we connect the $`A`$ sublattice to the left side (Fig. 12$`(b)`$). The current expression is then the same as Eq. (79), where $`g_0`$ is replaced by the 11 (or $`AA`$) component of the right hand side of Eq. (101). The dotted line of Fig. 13 shows the conductance curve for the armchair case.

## IV Summary

In summary, the generalized method of image that we developed has allowed us to deal with various tunneling problems in a unified manner, with full tight-binding nature being taken into account. In particular, we have applied it in this paper to examine in-plane tunneling spectra of normal metal–$`d`$-wave superconductor junctions, with and without external magnetic fields, at arbitrary crystalline orientations. The doping dependence of the ZBCPs is also studied within the mean-field slave boson approach. Our results for the splitting of ZBCP under applied magnetic field agrees well with recent experiments. We further showed that tunneling into $`d`$-density-wave state at (110) orientation should display a sharp conductance peak at the chemical potential in the tunneling spectra. This peak will shift away from the chemical potential if the next nearest neighbor hopping $`t_R^{}`$ exists, which also offers a way to measure $`t_R^{}`$. Under in-plane magnetic fields, it also splits due to Zeeman splitting. These provide signatures to be looked for in experiments, especially in normal-metal–pseudogap-cuprate junctions for testing the proposal of Ref. chakravarty, . As a demonstration of the general applicability of our formulation, we further consider tunneling into graphite sheets at the zigzag and armchair orientations. ZBCP is found in the zigzag case while no ZBCP should be displayed in the armchair case, consistent with findings in the study of graphite ribbons.fujita  We analyze these results on the basis of the 1D $`t_1`$-$`t_2`$ model and obtain the criteria for the emergence of the ZBCP.

The merit of our formulation lies in two aspects. Firstly, it offers a unified method for theoretical study of the tunneling spectroscopy. Recently, applying this method, we discover a remarkable even-odd effect in semi-infinite 1D chains: for hopping amplitudes with even cycles there can be zero-energy localized edge-states, while for those with odd cycles there is none.wm2  Secondly, as already pointed out at the beginning of the paper, for the first time, our method allows us to express what is being measured in tunneling experiments in terms of bulk Green’s function. For instance, in a single hard-wall configuration, suppose the normal metal on the left is a wideband material and the junction barrier is very high ($`t1`$ in Eqs. (38)–(41) and (79)), tunneling experiments essentially measure the quantity

$`\mathrm{d}I/\mathrm{d}V`$ $``$ $`{\displaystyle \underset{𝐤\sigma }{}}Im\{g_0(𝐤,eV)\}`$
$`=`$ $`{\displaystyle \underset{𝐤\sigma }{}}Im\{G(𝐤,eV)\times [1\mathrm{exp}(2ik_xd)\alpha _0]\}.`$
Namely it probes the surface density of state, which is the bulk Green’s function modulated by the factor in the square bracket. In our formulation the surface Green’s function, and hence the surface density of states, decomposes into two parts: one from the real source and the other from the image source. The real-source part results purely from the bulk and hence reveals faithfully the bulk property; the image source part contains all surface effects and hence is responsible for any complications. In the presence of the reflection symmetry, we find that $`\alpha _0=1`$ and the image part contribute another bulk term. Thus in this case the conductance curve exhibits only the bulk property (with, however, the van Hove singularity “rounded” by a sine factor as in Eq. (13)). When the reflection symmetry is broken, singular behavior may arise from the image part. This singularity is contained in $`\alpha _0`$ and originates from the zeros of the Green’s function (Eq. (68)) or its determinant when considering superconductors (Eq. (43)).

###### Acknowledgements.

The authors would like to thank Profs. Sungkit Yip, Hsiu-Hau Lin, T. K. Lee, Hu Xiao and C. C. Chi for useful discussions. This research was supported by NSC of Taiwan.

## Appendix A Current expressions

In this appendix we outline techniques for calculating the tunneling currents for the ND and the N-DDW junctions. We start from the expression

$`I(t)=+ie{\displaystyle \underset{k_y,\sigma }{}}t(k_y)c_{l\sigma }^{}(k_y)c_{r\sigma }(k_y)+\mathrm{h}.\mathrm{c}.`$ (A$`\mathrm{\hspace{0.33em}103}`$)
In Keldysh’s formulation of non-equilibrium Green’s functions, the expectation values $`c_{l\sigma }^{}c_{r\sigma }`$ in the above equation can be expressed as a perturbation series. Under certain approximations one can resum this series to all orders in the tunneling amplitude $`t`$.

### A.1 ND junctions

We consider first ND tunneling. In dealing with superconducting phenomena it is convenient to use the Nambu representation which explicitly distinguishes particles and holes by assigning them to different components. One thus defines the spinor field-operators

$`\mathrm{\Psi }_\alpha (x_i,k_y,t)=\left(\begin{array}{c}\mathrm{\Psi }_{\alpha ,1}\\ \mathrm{\Psi }_{\alpha ,2}\end{array}\right)=\left(\begin{array}{c}c_\alpha (x_i,k_y,t)\\ c_\alpha ^{}(x_i,k_y,t)\end{array}\right),`$ (A$`\mathrm{\hspace{0.33em}108}`$)
where $`\alpha =\{L,R\}`$ labels the electrodes and the upper and lower elements are associated with, respectively, electrons and holes. The Keldysh non-equilibrium Green’s functions are then defined asyeyati ; cuevas

$`g_{\alpha \beta ,\mu \nu }^+(x_i,t;x_j,t^{})=+i\mathrm{\Psi }_{\beta ,\nu }^{}(x_j,t^{})\mathrm{\Psi }_{\alpha ,\mu }(x_i,t),`$ (A$`\mathrm{\hspace{0.33em}109}`$)
$`g_{\alpha \beta ,\mu \nu }^+(x_i,t;x_j,t^{})=i\mathrm{\Psi }_{\alpha ,\mu }(x_i,t)\mathrm{\Psi }_{\beta ,\nu }^{}(x_j,t^{}).`$ (A$`\mathrm{\hspace{0.33em}110}`$)
For brevity we have suppressed the $`k_y`$ dependence. The Green’s functions here carry the left right indices $`\alpha ,\beta =\{L,R\}`$, the Nambu (spinor) indices $`\mu ,\nu =\{1,2\}`$, and the Keldysh indices $`\{,+\}`$. For notational clarity we shall in the following frequently omit irrelevant indices and keep track of only those related to our discussion.

In this representation we define the tunneling matrix

$`\widehat{t}t\tau _3\sigma _3`$ (A$`\mathrm{\hspace{0.33em}111}`$)
where $`\tau _3`$ and $`\sigma _3`$ are the third Pauli matrices pertaining to the Nambu space and the Keldysh space, respectively. In particular, $`\sigma _3`$ is chosen so that in the Keldysh space

$`\sigma _3^{}=1=\sigma _3^{++},\mathrm{and}\sigma _3^+=0=\sigma _3^+,`$
since we have assigned the forward time-path the “$``$” time axis, and the return time-path the “+” time axis. In the following we will consider only real valued $`t`$ and hence $`t^{}=t`$.

The current expression (A$`\mathrm{\hspace{0.33em}103}`$) can now be written as

$`I(t)=+e{\displaystyle \underset{k_y}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\omega }{2\pi }}t(k_y)`$
$`\times \left\{Tr[g_{RL}^+(x_0,k_y,\omega )]Tr[g_{LR}^+(x_0,k_y,\omega )]\right\}.`$ (A$`\mathrm{\hspace{0.33em}112}`$)
where the trace is taken over the Nambu space and we have Fourier transformed over the time variables. In the presence of particle-hole symmetry the Nambu components 11 and 22 in the trace above contribute equally. Therefore, the trace yields twice the contribution from the 11 component. Applying the Dyson equationskeldysh

$`g=g_0+g_0\widehat{t}g=g_0+g\widehat{t}g_0`$ (A$`\mathrm{\hspace{0.33em}113}`$)
and noting that $`g_{0,RL}=0=g_{0,LR}`$, as there is no propagation between the two electrodes at the bare level, one finds

$`g_{RL}^+`$ $`=`$ $`\left[g_0+g\widehat{t}g_0\right]_{RL}^+`$ (A$`\mathrm{\hspace{0.33em}114}`$)
$`=`$ $`t(g_{RR}^{}g_{0,LL}^+g_{RR}^+g_{0,LL}^{++})`$
$`g_{LR}^+`$ $`=`$ $`\left[g_0+g_0\widehat{t}g\right]_{LR}^+`$ (A$`\mathrm{\hspace{0.33em}115}`$)
$`=`$ $`t(g_{0,LL}^+g_{RR}^{++}+g_{0,LL}^{}g_{RR}^+).`$
Since the components of Keldysh Green’s functions has the property $`g^{++}+g^{}=g^++g^+`$. Taking the difference between $`g_{RL}^+`$ and $`g_{LR}^+`$ using Eqs. (A$`\mathrm{\hspace{0.33em}114}`$) and (A$`\mathrm{\hspace{0.33em}115}`$) and inserting the result into (A$`\mathrm{\hspace{0.33em}112}`$) we obtain

$`I=2e{\displaystyle \underset{k_y}{}}t^2{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\omega }{2\pi }}\{g_{0,LL,11}^+(\omega eV)g_{RR,11}^+(\omega )`$
$`g_{0,LL,11}^+(\omega eV)g_{RR,11}^+(\omega )\}.`$ (A$`\mathrm{\hspace{0.33em}116}`$)
Note that the frequency arguments of the bare Green’s functions for the left electrode $`g_{0,LL}^{+/+}`$ has been shifted due to the applied bias $`eV`$ between the two sides ($`\mu _L\mu _R=eV`$). We emphasize that at this stage Eq. (A$`\mathrm{\hspace{0.33em}116}`$) is exact and our remaining task is to express the exact Green’s functions $`g_{RR,11}^{+/+}`$ in terms of the uncoupled Green’s functions $`g_0`$ via some approximations.

The Green’s functions $`g^{+/+}`$ can be expressed in terms of the bare Green’s functions $`g_0^{+/+}`$ and the exact Green’s functions $`g^{r,a}`$ by means of the following equationskeldysh

$`g^{+/+}(\omega )=[1+g^r(\omega )\widehat{t}]g_0^{+/+}[\widehat{t}g^a(\omega )+1].`$ (A$`\mathrm{\hspace{0.33em}117}`$)
To fully reduce to bare quantities, one can express the exact Green’s functions $`g^{r,a}`$ in terms of the bare Green’s functions $`g_0^{r,a}`$ by virtue of the Dyson equations

$`g^{r,a}(\omega )=g_0^{r,a}(\omega )+g_0^{r,a}(\omega )\widehat{t}g^{r,a}(\omega ).`$ (A$`\mathrm{\hspace{0.33em}118}`$)
Ignoring many-body effects, we are able to resum the series and obtain

$`g_{RR}^{r,a}(\omega )=𝒯_{RL}^{r,a}(\omega )g_{0,RR}^{r,a}(\omega ),`$ (A$`\mathrm{\hspace{0.33em}119}`$)
$`g_{LR}^{r,a}(\omega )=𝒯_{LR}^{r,a}(\omega )\left[tg_{0,LL}^{r,a}(\omega eV)\tau _3g_{0,RR}^{r,a}(\omega )\right],`$ (A$`\mathrm{\hspace{0.33em}120}`$)
$`g_{RL}^{r,a}(\omega )=𝒯_{RL}^{r,a}(\omega )\left[tg_{0,RR}^{r,a}(\omega )\tau _3g_{0,LL}^{r,a}(\omega eV)\right],`$ (A$`\mathrm{\hspace{0.33em}121}`$)
$`g_{LL}^{r,a}(\omega )=𝒯_{LR}^{r,a}(\omega )g_{0,LL}^{r,a}(\omega eV),`$ (A$`\mathrm{\hspace{0.33em}122}`$)
where sum over tunneling processes of all orders is signified by the factors

$`𝒯_{RL}^{r,a}(\omega )`$ $`=`$ $`\left[1t^2g_{0,RR}^{r,a}(\omega )\tau _3g_{0,LL}^{r,a}(\omega eV)\tau _3\right]^1,`$
$`𝒯_{LR}^{r,a}(\omega )`$ $`=`$ $`\left[1t^2g_{0,LL}^{r,a}(\omega eV)\tau _3g_{0,RR}^{r,a}(\omega )\tau _3\right]^1.`$
Since $`\mu _L\mu _R=eV`$ all frequency arguments of $`g_{0,LL}`$ are to be taken at $`\omega eV`$. In arriving at the above equations we have used for the normal electrode

$`g_{0,LL}(\omega eV)=(\begin{array}{cc}g_{0,LL,11}(\omega eV)& 0\\ 0& g_{0,LL,22}(\omega +eV)\end{array}).`$
Note that the 11 and the 22 elements of $`g_{0,LL}`$ represent particles and holes, respectively, and hence their response to applied bias is opposite. This is essential in giving rise to the Andreev contributions in the tunneling current.

Incorporating Eqs. (A$`\mathrm{\hspace{0.33em}119}`$)–(A$`\mathrm{\hspace{0.33em}122}`$) with (A$`\mathrm{\hspace{0.33em}117}`$), one can thus obtain $`g_{RR}^{+/+}`$ and substitute back into (A$`\mathrm{\hspace{0.33em}116}`$). Finally we make use of the following relations

$`g_0^+(\omega )`$ $`=`$ $`2\pi if(\omega )\widehat{A}(\omega ),`$ (A$`\mathrm{\hspace{0.33em}124}`$)
$`g_0^+(\omega )`$ $`=`$ $`2\pi i[1f(\omega )]\widehat{A}(\omega ),`$ (A$`\mathrm{\hspace{0.33em}125}`$)
where $`f(\omega )`$ is the Fermi function and $`\widehat{A}(\omega )`$ is the spectral weight matrix given by (42). This leads us to the current expressions (38)–(41).

### A.2 N-DDW junctions

We turn now to deriving the current expressions for N-DDW junctions which is also applicable to NG junctions. We shall also show that in this case Andreev-like processes do not contribute to the tunneling current. In the absence of external fields, spin degree of freedom merely introduces a factor of two. Thus the spin indices $`\sigma `$ will be omitted in the following.

We first define the Keldysh Green’s functions similarly to (A$`\mathrm{\hspace{0.33em}110}`$)

$`g_{\alpha \beta }^+(x_i,t;x_j,t^{})=+ic_\beta ^{}(x_j,t^{})c_\alpha (x_i,t),`$ (A$`\mathrm{\hspace{0.33em}126}`$)
$`g_{\alpha \beta }^+(x_i,t;x_j,t^{})=ic_\alpha (x_i,t)c_\beta ^{}(x_j,t^{}).`$ (A$`\mathrm{\hspace{0.33em}127}`$)
Here the subscripts $`\alpha `$, $`\beta =\{R,L\}`$ are labels for the electrodes (not to be confused with the labels for sublattices in the text). In terms of the Keldysh Green’s functions the tunneling current (A$`\mathrm{\hspace{0.33em}103}`$) can be written

$`I(t)=+e{\displaystyle \underset{l,r,\sigma }{}}[tg_{RL}^+(r,l)t^{}g_{LR}^+(l,r)].`$ (A$`\mathrm{\hspace{0.33em}128}`$)
Similar to the previous section, the renormalized Green’s functions $`g_{RL}`$ and $`g_{LR}`$ can be expressed as combinations of the bare Green’s functions $`g_{0,LL}`$ and the renormalized Green’s function $`g_{RR}`$. This results in the exact formula

$`I=e{\displaystyle \underset{k_y,\sigma }{}}t^2{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{d}\omega }{2\pi }}\{g_{0,LL}^+(\omega eV)g_{RR}^+(\omega )`$
$`g_{0,LL}^+(\omega eV)g_{RR}^+(\omega )\}.`$ (A$`\mathrm{\hspace{0.33em}129}`$)
Note that here, unlike the previous section, the tunneling matrix is

$`\widehat{t}=t\sigma _3,`$ (A$`\mathrm{\hspace{0.33em}130}`$)
where $`\sigma _3`$ is the third Pauli matrix in the Keldysh space.

To proceed, as in the ND case, we apply (A$`\mathrm{\hspace{0.33em}117}`$) and expand the exact Green’s functions $`g_{RR}^{+/+}`$ in terms of the bare Green’s functions $`g_0^{+/+}`$ and the renormalized retarded and advanced Green’s functions $`g^{r,a}`$. Again, ignoring many-body effects we can express the renormalized Green’s functions $`g^{r,a}`$ in terms of the bare Green’s functions $`g_0^{r,a}`$. Utilizing (A$`\mathrm{\hspace{0.33em}124}`$) and (A$`\mathrm{\hspace{0.33em}125}`$), one has

$`g_{RR}^+=2\pi i[f(\omega )M_R(\omega )+f(\omega eV)M_L(\omega )]`$
$`g_{RR}^+=2\pi i[(1f(\omega ))M_R(\omega )+(1f(\omega eV))M_L(\omega )]`$
with

$`M_R(\omega )`$ $`=`$ $`A_R(\omega )|1+tg_{RL}^r(\omega )|^2,`$ (A$`\mathrm{\hspace{0.33em}131}`$)
$`M_L(\omega )`$ $`=`$ $`t^2A_L(\omega eV)|g_{RR}^r(\omega )|^2.`$ (A$`\mathrm{\hspace{0.33em}132}`$)
In the last expressions we have used $`g_{RL}^r=(g_{LR}^a)^{}`$. Note that $`M_L`$ contains the spectral weight $`A_L`$ of the normal electrode. It is associated with tunneling processes where a particle is reflected back into the left side and at the same time a particle-hole pair is transmitted into the right; this is reminiscent of the Andreev channel in ND tunneling (cf. the integrand in (41)). Substituting the above results into the current expression (A$`\mathrm{\hspace{0.33em}129}`$), we obtain for the terms in the braces in the integrand

$`4\pi ^2[f(\omega eV)f(\omega )]A_L(\omega eV)M_R(\omega ),`$ (A$`\mathrm{\hspace{0.33em}133}`$)
which leads to the current expression (79).

It is remarkable that $`M_L`$ is canceled completely in the final current expression, and hence no contribution from the Andreev process remains in the tunneling current. The reason for this difference between the ND and the N-DDW (and likewise NG) tunneling is that for the former the bias voltage shifts the energy of particles and holes in opposite directions, while for the latter particles from different bands acquire the same shift under applied bias.


format: text

text:
# Josephson Current in S-FIF-S Junctions: Nonmonotonic Dependence on Misorientation Angle

## Abstract

Spectra and spin structures of Andreev interface states in S-FIF-S junctions are investigated with emphasis on finite transparency and misorientation angle $`\phi `$ between in-plane magnetizations of ferromagnetic layers in a three-layer interface. It is demonstrated that the Josephson current in S-FIF-S quantum point contacts can exhibit a nonmonotonic dependence on the misorientation angle. The characteristic behavior takes place, if the $`\pi `$-state is the equilibrium state of the junction in the particular case of parallel magnetizations.

The dc Josephson effect in junctions with ferromagnetic interfaces exposes remarkable features which have been intensively studied in recent years theoretically bul77 ; buz82 ; mrs88 ; buz92 ; demler97 ; yip00 ; heik00 ; fogel00 ; fogel01 ; koshina011 ; koshina012 ; volkov01 ; nazarov01 ; bb02 ; waintal02 ; gol021 ; gol022  and experimentally ryazanov01 ; kontos02 . Apart from interfaces with a fixed magnetization, considerable attention has been drawn also to more complicated cases, when the magnetization is spatially dependent inside the interface. An important particular model for this kind of interfaces is a three-layer FIF-interface, where two metallic ferromagnetic layers with in-plane magnetizations, making angle $`\phi `$ with each other, are separated by an insulating magnetically inactive interlayer koshina011 ; koshina012 ; volkov01 ; waintal02 ; gol021 . In the present article we identify theoretically spectra and spin polarization of Andreev states bound to the three-layer FIF-constrictions with finite transmission, separating clean $`s`$-wave superconductors. Then we determine the Josephson current in the S-FIF-S quantum point contacts. This problem was not studied previously in the literature. In the dirty limit, considered in koshina011 ; koshina012 ; gol021 , Andreev bound states are fully destroyed. Considerations of Ref. volkov01,  imply the absence of Andreev bound states in clean S-FIF-S junctions. This can be justified only for short superconductors, whose lengths are less than their coherence lengths. Spectra of Andreev states and the Josephson current in S-FIF-S junctions with $`\phi 0`$ have been found earlier only in the particular limit of fully transparent constriction which was characterized by the authors as a toy model waintal02 .

The misorientation angle $`\phi `$ can be considered, in general, as a variable quantity. Let, for instance, the magnetization axis be pinned in one ferromagnetic layer, while in the other one there is an easy in-plane magnetization layer. Then one can vary the misorientation angle (keeping other parameters fixed) by applying an external weak magnetic field to the interface. We find that the Josephson current as a function of the misorientation angle $`\phi `$ manifests a characteristic nonmonotonic behavior, if, at $`\phi =0`$, the $`\pi `$-state is the equilibrium state of the junction.

For our analysis, we examine a smooth plane interface between two superconductors which consists of two layers of the same ferromagnetic metal separated by an insulating nonmagnetic barrier. Two identical ferromagnetic layers are characterized by their thickness $`l`$ and internal exchange fields $`|𝒉_{1,2}|=h`$, which, being parallel to the layers, make an angle $`\phi `$ with each other.

The normal-state scattering $`𝒮`$ matrix is represented as $`𝒮=S(1+\widehat{\tau }_z)/2+\stackrel{~}{S}(1\widehat{\tau }_z)/2`$, where the Pauli-matrices $`\widehat{\tau }_j`$ are taken in particle-hole space and $`\stackrel{~}{S}(p_{})=S^T(p_{})`$. Each component $`\widehat{S}_{ij}`$ in matrix $`S=\widehat{S}_{ij}`$ ($`i(j)=1,2`$) is in its turn a matrix in spin space. Matrix $`\widehat{S}_{ii}=\left(\begin{array}{cc}r_i^{}& r_i^{}\\ r_i^{}& r_i^{}\end{array}\right)`$ contains, in general, spin-dependent interface reflection amplitudes for normal-state quasiparticles in $`i`$-th half-space, while $`\widehat{S}_{ij}=\left(\begin{array}{cc}d_{ij}^{}& d_{ij}^{}\\ d_{ij}^{}& d_{ij}^{}\end{array}\right)`$ with $`ij`$ incorporates spin-dependent transmission amplitudes for normal-state quasiparticles from side $`j`$. For the interface potentials conserving particle current, the scattering matrix has to satisfy the unitarity condition: $`𝒮𝒮^{}=1`$. If the interface Hamiltonian possesses time-reversal symmetry, one obtains an additional constraint on the scattering matrix: $`S(𝒑_f,𝒉_{1,2})=\widehat{\sigma }_yS^T(\underset{¯}{𝒑}_f,𝒉_{1,2})\widehat{\sigma }_y`$ mrs88 . Assume, for simplicity, the exchange fields to be much smaller compared to the Fermi energies. For the diagonalization of the $`S_{11}`$-matrix we choose the $`z`$-axis along the magnetization in the first (left) ferromagnetic layer. Then the other $`S_{ij}`$-matrices are nondiagonal unless $`\phi =0,\pi `$:

$`\widehat{S}_{21}=d\mathrm{exp}\left({\displaystyle \frac{i\mathrm{\Theta }}{4}}\left(\widehat{\sigma }_y\mathrm{sin}\phi +\widehat{\sigma }_z\mathrm{cos}\phi \right)\right)\mathrm{exp}\left({\displaystyle \frac{i\mathrm{\Theta }}{4}}\widehat{\sigma }_z\right),`$
$`\widehat{S}_{12}=d\mathrm{exp}\left({\displaystyle \frac{i\mathrm{\Theta }}{4}}\widehat{\sigma }_z\right)\mathrm{exp}\left({\displaystyle \frac{i\mathrm{\Theta }}{4}}\left(\widehat{\sigma }_y\mathrm{sin}\phi +\widehat{\sigma }_z\mathrm{cos}\phi \right)\right),`$
$`\widehat{S}_{11}=r\mathrm{exp}\left(i{\displaystyle \frac{\mathrm{\Theta }}{2}}\widehat{\sigma }_z\right),`$
$`\widehat{S}_{22}=\stackrel{~}{r}\mathrm{exp}\left(i{\displaystyle \frac{\mathrm{\Theta }}{2}}\left(\widehat{\sigma }_y\mathrm{sin}\phi +\widehat{\sigma }_z\mathrm{cos}\phi \right)\right).`$ (1)
Here $`\mathrm{\Theta }=\frac{4lh}{\mathrm{}\mathrm{v}_{f,x}}`$. Quantities $`r`$, $`\stackrel{~}{r}`$ and $`d`$ are the respective reflection and transmission amplitudes of the potential barrier $`V`$, satisfying the condition $`rd^{}=d\stackrel{~}{r}^{}`$.

We carry out calculations within the quasiclassical theory of superconductivity, based on the equations for the so-called Riccati amplitudes or coherence functions schopohl98 ; eschrig00 ; fogel00 . Considering a quantum point contact with FIF-constriction, the order parameter is taken spatially constant. We include interface exchange fields in the quasiclassical boundary conditions. Since they imply, as usually, that all interface potentials are much larger than the superconducting order parameter zai84 , one should assume $`|h_{1,2}|\mathrm{\Delta }`$.

With the above normal-state $`𝒮`$ matrix we get four branches of interface Andreev bound states:

$$\epsilon _{1,2}=|\mathrm{\Delta }|\mathrm{cos}\frac{\mathrm{\Phi }_{1,2}}{2},\epsilon _{3,4}=|\mathrm{\Delta }|\mathrm{cos}\frac{\mathrm{\Phi }_{1,2}}{2},$$
(2)
where the quantities $`\mathrm{\Phi }_{1,2}(\chi ,\mathrm{\Theta },\phi )`$ are defined as

$`\mathrm{cos}\mathrm{\Phi }_{1,2}(\chi ,\mathrm{\Theta },\phi )=\mathrm{cos}\mathrm{\Theta }2D\mathrm{cos}\mathrm{\Theta }\mathrm{sin}^2{\displaystyle \frac{\chi }{2}}+`$
$`+2D\mathrm{cos}\chi \mathrm{sin}^2{\displaystyle \frac{\mathrm{\Theta }}{2}}\mathrm{sin}^2{\displaystyle \frac{\phi }{2}}\pm 2\sqrt{D}\mathrm{sin}{\displaystyle \frac{\chi }{2}}\mathrm{sin}\mathrm{\Theta }\mathrm{cos}{\displaystyle \frac{\phi }{2}}\times `$
$`\times \sqrt{1D\mathrm{sin}^2{\displaystyle \frac{\chi }{2}}+D\mathrm{cos}^2{\displaystyle \frac{\chi }{2}}\mathrm{tan}^2{\displaystyle \frac{\mathrm{\Theta }}{2}}\mathrm{sin}^2{\displaystyle \frac{\phi }{2}}}.`$ (3)
Here, $`\chi `$ is the phase difference of the two superconductors. The energies $`\epsilon _i`$ ($`i=1,2,3,4`$) implicitly depend on quasiparticle momentum directions via the parameter $`\mathrm{\Theta }`$ and the transmission coefficient $`D`$.

For $`\phi =0`$ the spectra Eq. (3) reduce to those found in Ref. fogel01, . In the particular case of antiparallel orientation of the left and the right magnetization $`\phi =\pi `$, the spectra of Andreev interface states (2), (3) take the form

$$\epsilon _1=\epsilon _2=\epsilon _3=\epsilon _4=|\mathrm{\Delta }|\sqrt{D\mathrm{cos}^2\frac{\chi }{2}+R\mathrm{cos}^2\frac{\mathrm{\Theta }}{2}}.$$
(4)
Being symmetric with respect to the transformation $`\mathrm{\Theta }\mathrm{\Theta }`$, the spectrum (4) is doubly degenerated. In the limit of a nonmagnetic interface ($`\mathrm{\Theta }=0`$), our result, Eqs. (2) and (3), leads to a well known spectrum of spin-degenerate interface Andreev bound states fur90 ; fur191 ; been191 ; been291  $`\epsilon _B=\pm |\mathrm{\Delta }|\sqrt{1D\mathrm{sin}^2(\chi /2)}`$.

Quasiparticle spin is a qood quantum number in the BCS theory of superconductivity, if one can disregard spin-flip effects stimulated, for instance, by magnetic impurities, spin-orbit interactions or magnetically active interfaces. In the presence of a paramagnetic spin interaction with an externally applied magnetic field or an internal exchange field, spin degeneracy of quasiparticle states is lifted and only states with parallel or antiparallel spin-to-field orientations are still eigenstates of the problem. This can lead to effects having physics common with the Larkin-Ovchinnikov-Fulde-Ferrell state ff65 ; lo65  and, in particular, associated with opposite signs of the Zeeman terms for electrons forming a Cooper pair in singlet superconductors.

A Bogoliubov quasiparticle in the superconductor has well defined spin, although its electron and hole components are described with Zeeman terms of opposite signs. Also, an electron and its Andreev reflected partner (hole) at an interface, separating singlet superconductors and leading to no spin-flip processes, have identical spin orientations and opposite signs of Zeeman terms. With opposite velocity directions and electric charges, while in identical spin states, they carry jointly the electric supercurrent across the interface, but no equilibrium spin current. Hence, definite spin polarization of interface Andreev bound states is fully compatible with the fact that Cooper pairs in singlet supercoductors carry no spin current.

Andreev states bound to nonmagnetic interfaces are spin degenerated. For a ferromagnetic interface with uniformly oriented magnetization Andreev interface states are spin polarized, being parrallel or antiparallel to the magnetization. The ferromagnetic interface lifts spin degeneracy of the Andreev states, but still does not mix the spin-polarized channels, carrying the Josephson current bb02 . This is not the case, however, if various orientations of magnetization are present in the interface, as this takes place in the FIF-interface with $`\phi 0`$. Quasiparticle Andreev interface states with the spectra of Eqs. (2), (3) are characterized by a nontrivial spin structure, which substantially depends (together with the spectra themselves) on $`\phi `$, $`\mathrm{\Theta }`$ and $`D`$. In general, each of the two incoming and two outgoing parts of quasiparticle trajectories, forming an Andreev interface state, has its own specific spin polarization. This should be compatible with no spin current they carry, on the whole, across the interface.

Figs. 1 and 2 demonstrate the evolution of spectra and spin polarizations of four branches of Andreev interface states as functions of $`\mathrm{\Theta }`$ in tunnel junctions (with transparency $`D=0.1`$) and highly transparent junctions ($`D=0.95`$) respectively. Two particular relative orientations of magnetization $`\phi =0.1\pi `$ (left column) and $`\phi =2\pi /3`$ (right column) are chosen. For definiteness, we consider spin polarizations of Andreev states on the incoming part of the quasiparticle trajectory in the right superconductor. The spin polarization gradually changes with the parameter $`\mathrm{\Theta }`$ in all cases considered. A characteristic scale of $`\mathrm{\Theta }`$ for the spin reconstruction diminishes with decreasing the misorientation angle $`\phi `$. For vanishing $`\phi `$ the scale vanishes and there are abrupt jumps from parallel to antiparallel (or vice versa) spin orientations with respect to the magnetization, taking place at those values of $`\mathrm{\Theta }`$, where $`\epsilon _i(\mathrm{\Theta })=\pm \mathrm{\Delta }`$bb02 .

Only in the particular case $`\phi =0`$, when a single magnetization direction is present inside a symmetric magnetic interface, each of the Andreev interface states possesses, as a whole, a definite spin-up or spin-down polarization, identical for all incoming and outgoing quasiparticle trajectories forming the state. Then the spectra of the spin-up polarized Andreev bound states are bb02 : $`\epsilon _{1,2}^{}=|\mathrm{\Delta }|\mathrm{sgn}(\mathrm{sin}\frac{\mathrm{\Phi }_{1,2}\left(\chi ,\mathrm{\Theta },\phi =0\right)}{2})\mathrm{cos}\frac{\mathrm{\Phi }_{1,2}\left(\chi ,\mathrm{\Theta },\phi =0\right)}{2}`$. The energies $`\epsilon _{1,2}^{}`$ for spin-down Andreev bound states are obtained from $`\epsilon _{1,2}^{}`$ by substituting $`\mathrm{\Theta }\mathrm{\Theta }`$.

The spectra of Andreev states and their spin polarizations as functions of the misorientation angle $`\phi `$ are shown in Fig. 3. The spin polarization at $`\phi 0`$ makes a finite angle with both magnetization directions and differs on different incoming and outgoing trajectories related by the bound state. As already mentioned above, for antiparallel magnetizations ($`\phi =\pi `$) the spectra are doubly degenerated. Spin polarizations, shown in Fig. 3 for $`\phi =\pi `$, can be considered as correct eigenfunctions in zeroth order approximation in small deviations $`\phi \pi `$.

The spin structure of Andreev interface states at nonzero $`\phi `$ should be taken into account in producing nonequilibrium occupation of the states. For $`\phi =0`$ only the interlevel transitions accompanied with spin-flip processes are possible under certain conditions bb02 . On account of the complicated spin structure of the Andreev states at nonzero $`\phi `$, there are actually no strict restrictions to a change of quasiparticle spin in the transition.

The Josephson current is carried by the bound states (2), analogously to the situation in nonmagnetic symmetric junctions fur90 ; fur191 ; been191 ; been291 . Hence, in a quantum point contact with a FIF constriction the total Josephson current carried by four Andreev states (2) can be found as $`J(\chi ,T)=2e\underset{m}{}\frac{d\epsilon _m}{d\chi }n(\epsilon _m)`$ $`=2e\underset{\epsilon _m>0}{}\frac{d\epsilon _m}{d\chi }\mathrm{tanh}\frac{\epsilon _m}{2T}`$ . It is not difficult to calculate the current in this scheme numerically. The Josephson critical current as a function of the misorientation angle $`\phi `$, normalized to its value at $`\phi =0`$, is shown in Fig. 4 for various $`\mathrm{\Theta }`$ and for two values of the transparency $`D=0.01,0.8`$ (the upper and the lower panels respectively) and the temperature $`T=0.1T_c,0.8T_c`$ (the right and the left columns).

We define the critical current as a positive quantity, as it is usually determined experimentally. There are two qualitatively different regimes for the behavior of the Josephson critical current as a function of $`\phi `$ in all particular cases displayed in Fig. 4. The two regimes are separated by a characteristic value $`\mathrm{\Theta }^{}(T,D)`$, which depends on the temperature and the transparency. For $`\mathrm{\Theta }<\mathrm{\Theta }^{}`$ the current is a monotonous function of the misorientation angle, reaching the maximum for the antiparallel orientation of the magnetizations. For $`\mathrm{\Theta }>\mathrm{\Theta }^{}`$ the current manifests, however, a nonmonotonic dependence on the misorientation angle with well pronounced minimum at some intermediate value of $`\phi `$ and the maximum at $`\phi =\pi `$. In the case $`\mathrm{\Theta }=\pi `$ the currents at $`\phi =0`$ and $`\phi =\pi `$ are equal to each other. The parameter $`\mathrm{\Theta }^{}`$ admits a simple physical interpretation, associated with the junction at $`\phi =0`$. At zero misorientation angle the junction acquires a uniformly oriented ferromagnetic interface. Then the Josephson current is equivalent to that studied in fogel00 ; fogel01 ; nazarov01 ; bb02 . It can be shown, that for $`\phi =0`$ and $`\mathrm{\Theta }=\mathrm{\Theta }^{}(T,D)`$ the $`0\pi `$ transition takes place in the junction just at the given temperature $`T`$. Hence, for $`\mathrm{\Theta }>\mathrm{\Theta }^{}(T,D)`$ the equilibrium state of the junction with $`\phi =0`$ is the $`\pi `$-state, while for $`\mathrm{\Theta }<\mathrm{\Theta }^{}(T,D)`$ it is the $`0`$-state. We omit an analytical analysis of the total Josephson current in the case $`\phi =0`$, since the results exactly coincide with those obtained in Refs. nazarov01, ; bb02, .

Furthermore, there is no $`0\pi `$ transition in the junction with antiparallel magnetization, $`\phi =\pi `$, in the three-layer interface. Indeed, it is straightforward to get from Eq. (2) the Josephson current in the particular case $`\phi =\pi `$:

$`J(\chi ,T)={\displaystyle \frac{eD|\mathrm{\Delta }|\mathrm{sin}\chi }{\sqrt{D\mathrm{cos}^2{\displaystyle \frac{\chi }{2}}+R\mathrm{cos}^2{\displaystyle \frac{\mathrm{\Theta }}{2}}}}}\times `$
$`\times \mathrm{tanh}{\displaystyle \frac{|\mathrm{\Delta }|\sqrt{D\mathrm{cos}^2{\displaystyle \frac{\chi }{2}}+R\mathrm{cos}^2{\displaystyle \frac{\mathrm{\Theta }}{2}}}}{2T}}.`$ (5)
In contrast to the case $`\phi =0`$, the current (5) does not change its sign at any temperature. The same assertion is valid also for junctions with dirty superconductors, where Andreev states are fully destroyed koshina012 ; gol021 . We conclude, that the nonmonotonic dependence of the Josephson current on $`\phi `$ arises due to the $`0\pi `$ transition taking place with varying the misorientation angle at an intermediate value of $`\phi `$. This always occurs under the condition that there is a $`\pi `$-junction at $`\phi =0`$. If one defined the critical current in the $`\pi `$-junction as the negative quantity, then the nonmonotonic behavior would transform into the monotonic one, accompanied with a change of sign and discontinuity whenever $`\mathrm{min}|J_c|0`$.

The dependence of the Josephson current on the misorientation angle $`\phi `$ becomes especially simple and clear in the tunneling limit. In tunnel quantum point contacts the Josephson current takes the form $`J(T,\phi ,\chi )=J(T,\phi )\mathrm{sin}\chi `$ , where

$`J(T,\phi )=J^{(p)}(T)\mathrm{cos}^2{\displaystyle \frac{\phi }{2}}+J^{(a)}(T)\mathrm{sin}^2{\displaystyle \frac{\phi }{2}},`$ (6)
$`J^{(p)}(T)=eD|\mathrm{\Delta }|[\mathrm{cos}{\displaystyle \frac{\mathrm{\Theta }}{2}}\mathrm{tanh}\left({\displaystyle \frac{|\mathrm{\Delta }|\mathrm{cos}{\displaystyle \frac{\mathrm{\Theta }}{2}}}{2T}}\right)`$
$`{\displaystyle \frac{|\mathrm{\Delta }|}{2T}}{\displaystyle \frac{\mathrm{sin}^2{\displaystyle \frac{\mathrm{\Theta }}{2}}}{\mathrm{cosh}^2\left({\displaystyle \frac{|\mathrm{\Delta }|\mathrm{cos}\frac{\mathrm{\Theta }}{2}}{2T}}\right)}}],`$ (7)
$`J^{(a)}(T)={\displaystyle \frac{eD|\mathrm{\Delta }|}{\mathrm{cos}{\displaystyle \frac{\mathrm{\Theta }}{2}}}}\mathrm{tanh}\left({\displaystyle \frac{|\mathrm{\Delta }|\mathrm{cos}{\displaystyle \frac{\mathrm{\Theta }}{2}}}{2T}}\right).`$ (8)
The quantities $`J^{(p)}(T)`$, $`J^{(a)}(T)`$ are defined as $`J^{(p)}(T)J(T,\phi =0)`$$`J^{(a)}(T)J(T,\phi =\pi )`$, so that $`|J^{(p)}(T)|`$, $`|J^{(a)}(T)|`$ are critical currents in tunnel junctions with parallel and antiparallel orientations of the exchange fields in the three-layer interface. The dependence (6) on the misorientation angle has been derived in the tunneling limit earlier in Ref. volkov01, , disregarding the contribution from Andreev states and, hence, the $`0\pi `$ transition. As one can conclude from (6), the $`0\pi `$ transition can take place with varying $`\phi `$, if $`J_c^{(p)}(T)`$ and $`J_c^{(a)}(T)`$ have opposite signs. This is exactly the reason for a nonmonotonic dependence of the critical current on $`\phi `$ to show up. Eq. (6) is quite general within the tunneling limit and not applicable to highly transparent junctions. Spin polarizations of the eigenstates on any side of the impenetrable interface are aligned parallel or antiparallel to the respective magnetization direction. Making the projections of the spin polarized states from one side on the eigenstates on another side (with the spin polarization rotated by the angle $`\phi `$ with respect to the initial one), one confirms in the tunneling limit that the current is of the form (6).

In conclusion, we have investigated theoretically spectra and spin structures of interface Andreev states in S-FIF-S junctions. Both finite transparency and the misorientation angle between in-plane magnetizations of ferromagnetic layers were taken into account. We demonstrated that the Josephson critical current as a function of the misorientation angle always manifests a nonmonotonic behavior, if at $`\phi =0`$ the equilibrium state of the quantum point contact is the $`\pi `$-state.

Acknowledgments This work was supported by the Russian Foundation for Basic Research under Grant No. 02-02-16643 and by BMBF 13N6918/1.