text:
# Phase dynamics after connection of two separate Bose-Einstein condensates

Present address:\] Neworldcity, C/ Almagro 2, E-28010 Madrid, Spain

## Abstract

We study the dynamics of the relative phase following the connection of two independently formed Bose-Einstein condensates. Dissipation is assumed to be due to the creation of quasiparticles induced by a fluctuating condensate particle number. The coherence between different values of the phase, which is characteristic of the initial Fock state, is quickly lost after the net exchange of a few atoms has taken place. This process effectively measures the phase and marks the onset of a semiclassical regime in which the system undergoes Bloch oscillations around the initial particle number. These fast oscillations excite quasiparticles within each condensate and the system relaxes at a longer time scale until it displays low-energy, damped Josephson plasma oscillations, eventually coming to a halt when the equilibrium configuration is finally reached.

The way in which two Bose condensates, originally separately formed, behave when they are connected poses a fundamental problem in non-equilibrium macroscopic quantum mechanics. Two separate condensates 1 and 2 can be described by a relative number eigenstate $`|n`$ which in the phase representation may be written comment1

$$|n=_0^{2\pi }\frac{d\phi }{\sqrt{2\pi }}e^{in\phi }|\phi ,$$
(1)
where $`|\phi `$ is an eigenstate of the relative phase comment11 ; legg01  $`\phi =\phi _1\phi _2`$, with normalization $`\phi |\phi ^{}=\delta (\phi \phi ^{})`$, and $`n=(N_1N_2)/2`$ is the number of transferred atoms. The Fock state (1) may be viewed as resulting from a coherent superposition of all possible values of the relative phase. This poses no problem as long as the condensates are physically separated, e.g. by a barrier so high as to be effectively impenetrable, since then no physical observable depends on $`\phi `$. The situation changes when the two condensates are brought into connection, e.g. by lowering the barrier to the point where Josephson tunnelling becomes appreciable. Then a nonzero supercurrent develops which depends on the phase through the relation

$$\dot{n}=I(\phi )=(E_J/\mathrm{})\mathrm{sin}\phi ,$$
(2)
where $`E_J`$ is the Josephson coupling energy. If we interpret the Josephson relation (2) as an operator identity, we conclude that a system initially prepared in the state (1) will display a coherent superposition of states with macroscopically different values of the Josephson current (2). After some time, these different values of the current will give rise to macroscopically different values of the relative particle number

$$n(t)=n_0+_0^tI𝑑t.$$
(3)
The resulting macroscopic coherence is fragile and likely to be quickly lost by a tiny perturbation from the environment. After decoherence sets in, the relative phase is effectively measured and the system is expected to be describable by a reduced density matrix that is asymptotically diagonal in the phase representation. We will show that such a density matrix is a uniform mixture of pure state density matrices, each of them lying in the semiclassical regime where phase and number may be taken as simultaneously well-defined. In general the system finds itself very far from equilibrium and, due to the interaction with the dissipative environment, it eventually relaxes to the equilibrium configuration.

The purpose of this work is to analyze such a dynamical process. We assume that dissipation is caused by the spontaneous creation of quasiparticles within each well due to a fluctuating number of particles in the condensate that is caused in turn by the coherent exchange of atoms. Such an exchange can, of course, occur only in the presence of a physical coupling, e.g. of the Josephson type, between the two condensates. When the Josephson coupling is operating, the Hamiltonian of the system at low energies may be written zapa98 ; sols99

$$H=E_J(1\mathrm{cos}\phi )+E_cn^2/2,$$
(4)
where $`E_c`$ represents (up to a constant) the effects of the interatomic interaction for small values of $`n`$. When $`E_J`$ becomes nonzero, $`n`$ is no longer a constant of motion, as explicitly indicated by the equation of motion (2), and its dynamic fluctuations will induce excitation of quasiparticles within each well. The description of such a process will require going beyond the pendulum Hamiltonian (4) by including the dissipative coupling to the quasiparticle field. In practice, there may be other dissipative mechanisms, such as incoherent exchange of single atoms zapa98  or coupling to the walls of a cavity jaks01 , but quasiparticle creation is an intrinsic mechanism which lends itself to an approximate analytical treatment.

The problem which we address in this article is fundamental in several respects. First, it presents an important case of (phase) quantum measurement that can be treated in detail. Second, the connection between independent BEC’s has the essential features of a macroscopic interference experiment andr97 ; ande98 ; ande84 ; legg91 ; cast97 ; kohl01 , since it involves an initially random relative phase which becomes well defined, in this case due to the introduction of a Josephson coupling with the help of the quasiparticle field. Finally, the establishment of phase coherence between two recently connected condensates has bearing on the “Kibble problem” of formation of cosmological singularities in the early stages of the Universe kibb76 . Recent work on condensate coherence in optical lattices grei02  suggests that the study of phase dynamics after nonadiabatic switching of the Josephson coupling is within experimental reach.

The analytical treatment which we present in this article requires the introduction of some approximations. We have already said that dissipation is assumed to be caused by the creation of intrawell quasiparticles caused by a fluctuating condensate number. Following ref. meie00 , we take box normalization for the one-atom wave functions within each well. This is not an important limitation in the context of laser confined alkali condensates. The initial relative particle number is assumed to be sufficiently large for the semiclassical approximation to be applicable just after decoherence takes over between distant phase values. We also assume that the two independently prepared condensates are identical and neglect quasiparticle tunneling. For conceptual simplicity, we will focus here on the case in which macroscopic quantum self-trapping takes place smer97 ; milb97 . This requires $`E_cn_0^2>4E_J`$. In the future it will be desirable to extend the present analysis by relaxing some of these approximations.

The picture that emerges from our theoretical study is schematically depicted in (Fig. 1). There are at least three time scales in the connection process. First, decoherence between different phase values takes place after a few atoms have been exchanged. This process effectively measures the macroscopic phase. The second time scale is the Bloch oscillation period around the initial value of the relative particle number. These oscillations are a manifestation of the ac Josephson effect due to a chemical potential difference. Finally, the system relaxes in a longer time scale towards the equilibrium configuration. Damping is due to the irreversible energy transfer from the quickly oscillating number-phase degree of freedom to the quasiparticle field within each well.

First we perform a semiclassical study in which the relative phase and number are taken as $`c`$-numbers while the quasiparticle field is treated quantum mechanically. The semiclassical assumption is later justified from a full quantum treatment which explicitly shows how the quasiparticle field transforms the reduced density matrix of the phase-number coordinate from that of a pure Fock state to a mixture of Gaussian wave packets. This is the process that effectively “measures” the phase and makes it well-defined in a particular experimental realization.

At the semiclassical level (where the question of identity or not between the states $`|\phi `$ and $`|\phi +2\pi `$ can be neglected), it is clear from Eq. (4), that the problem maps on to that of the the “Bloch oscillations” of a particle in a tight-binding lattice which is subject to a harmonic potential, with $`n`$ playing the role of (discrete) coordinate and $`\phi `$ that a of a wave vector. A simple semiclassical analysis reveals that a particle placed with initial “wave vector” $`\phi _0`$ at a high point $`n_0`$ of the harmonic tight-binding lattice, undergoes periodic motion. Specifically, for $`E_cn_0^2E_J`$, and with an appropriate choice of time origin, we may write

$$n(t)=\overline{n}+p\mathrm{sin}(\omega _0t),$$
(5)
where $`p(E_J/\overline{n}E_c)`$, $`\omega _0E_c\overline{n}/\mathrm{}`$, and, initially, $`\overline{n}n_0`$. Equation (5) describes Bloch oscillations undergone by a quantum particle in a tilted lattice. The resulting self-trapping is not stable because the rapidly oscillating $`n(t)`$ induces excitations in the quasiparticle field, which acts as a dissipate environment for the macroscopic degree of freedom ruos98 ; legg77 . These oscillations are very fast and one is rather interested in the slower dynamics of the time average $`\overline{n}`$ as it experiences damping.

Since quasiparticles within each well are sensitive to the number of particles in the condensate to the extent that they are but linear deviations from the mean field description dalf99 , it is clear that (5) yields a time-dependent perturbation that will excite quasiparticles. The time-dependent Hamiltonian in e.g. the left well is fett71

$`H(t)={\displaystyle \underset{k}{}}\left[\epsilon _k+{\displaystyle \frac{g}{V}}\left({\displaystyle \frac{N}{2}}+n(t)\right)\right]c_k^{}c_k`$ (6)
$`+{\displaystyle \underset{k}{}}{\displaystyle \frac{g}{2V}}\left({\displaystyle \frac{N}{2}}+n(t)\right)(c_k^{}c_k^{}+c_kc_k),`$ (7)
where $`g`$ is the coupling constant, $`N`$ is the total particle number, and $`V`$ is the volume of one box. The Hamiltonian for the right well looks like (6) with $`n(t)`$ replaced by by $`n(t)`$. Eq. (6) may be rewritten as

$$H(t)=H_0+Ap\mathrm{sin}(\omega t),$$
(8)
where $`H_0`$ is the unperturbed “bath” Hamiltonian \[with $`n(t)`$ replaced by $`\overline{n}`$\] and

$$A=\frac{g}{2V}\underset{k}{}(2c_k^{}c_k+c_k^{}c_k^{}+c_kc_k).$$
(9)
By applying linear response theory forster  with $`\overline{n}`$ treated as constant, we obtain for the dissipated power within one well:

$$W=(p^2/2)\omega _0\chi ^{\prime \prime }(\omega _0),$$
(10)
where

$`\chi ^{\prime \prime }(\omega )`$ $``$ $`{\displaystyle \frac{1}{2\mathrm{}}}{\displaystyle 𝑑te^{i\omega t}[A(t),A(0)]}`$
$`=`$ $`{\displaystyle \frac{\pi g^2}{2\mathrm{}V^2}}{\displaystyle \underset{k}{}}|u_kv_k|^4\mathrm{coth}\left({\displaystyle \frac{\mathrm{}\omega _k}{2k_BT}}\right)\delta (\omega 2\omega _k)`$
$``$ $`{\displaystyle \frac{\mathrm{}V}{128\pi N^2c^3}}\omega ^4\mathrm{coth}\left({\displaystyle \frac{\mathrm{}\omega }{4k_BT}}\right).`$
In the third line, we have taken the low frequency approximation $`\mathrm{}\omega mc^2`$, $`c`$ being the sound velocity and $`m`$ the atom mass fett71 . Here $`u_k`$ and $`v_k`$ are the coherence factors of the Bogoliubov transformation, so that, in this limit fett71 , $`u_kv_k=(\mathrm{}Vck/Ng)^{1/2}`$, where $`N/2V`$ is the atom density in one well. We note that $`\chi ^{\prime \prime }(\omega )`$ goes like $`\omega ^4`$ at low temperatures and like $`\omega ^3`$ at high temperatures. The second line of Eq. (Phase dynamics after connection of two separate Bose-Einstein condensates) indicates that the quasiparticles are excited in pairs, which can be ascribed to the pairing structure of the depletion cloud fett72 ; legg01 .

The energy that is dissipated into the quasiparticle field has to be taken from the energy of the macroscopic pendulum (4) which, under the self-trapping assumption, is dominated by the interaction term. Thus, we write

$$2W=E_c\overline{n}\dot{\overline{n}},$$
(12)
where the factor of 2 accounts for the presence of two wells. We derive the differential equation

$$\dot{\overline{n}}=C\overline{n}^2\mathrm{coth}\left(\frac{E_c\overline{n}}{4k_BT}\right),$$
(13)
where $`C2\pi \eta E_J^2/\mathrm{}N^4E_c`$, with $`E_c=2g/V8\pi \mathrm{}^2a/mV`$ and $`\eta (Na^3/2\pi V)^{1/2}`$ the small parameter of the dilute gas. At low temperatures ($`k_BTE_c\overline{n}`$) we get

$$\overline{n}(t)=\frac{n_0}{1+\gamma _rt},$$
(14)
while, at high temperatures ($`k_BTE_c\overline{n}`$),

$$\overline{n}(t)=n_0\mathrm{exp}(\gamma _rt).$$
(15)
In (14) and (15), $`\gamma _r`$ is, respectively, the low and high temperature limit of the general expression

$$\gamma _r=\frac{\mathrm{}C}{E_c}\omega _0\mathrm{coth}\left(\frac{\mathrm{}\omega _0}{4k_BT}\right),$$
(16)
where $`\omega _0`$ refers to $`n_0`$.

As the dissipation process continues, the energy of the macroscopic pendulum decreases until it eventually undergoes low energy harmonic oscillations. In this limit, the analysis of the dissipation process is similar to that described above, with $`\omega _0`$ replaced by the Josephson plasma frequency $`\omega _{\mathrm{JP}}=\sqrt{E_JE_c}/\mathrm{}`$ and with the pendulum energy written as $`\stackrel{~}{n}^2E_c/2`$, where $`\stackrel{~}{n}`$ is the amplitude of the number oscillations. Since $`\omega _{\mathrm{JP}}`$ does not depend on $`\stackrel{~}{n}`$ for $`\stackrel{~}{n}N`$, we obtain for both low and high temperatures

$$\stackrel{~}{n}(t)=\stackrel{~}{n}(0)\mathrm{exp}(\gamma _rt),$$
(17)
$`\gamma _r`$ being given by (16) with $`\omega _0`$ replaced by $`\omega _{\mathrm{JP}}`$.

It remains to justify the semiclassical approximation. For that we have performed a full quantum calculation of the phase dynamics right after the connection has been established. The direct reversal process $`nn`$ may be safely neglected. We generalize quantum dissipation techniques weis99  to the case of quadratic coupling in the bath coordinates \[see Eq. (9)\]. After some lengthy algebra, we find that, in the phase representation, the reduced density matrix simply evolves as

$$\left|\rho (\phi ,\phi ^{};t)\right|=\frac{1}{2\pi }\mathrm{exp}\left[4\gamma _dt\mathrm{sin}^2\left(\frac{\phi \phi ^{}}{2}\right)\right],$$
(18)
where

$$\gamma _d=\gamma _rn_0\mathrm{coth}(\frac{E_cn_0}{2k_BT},).$$
(19)
For $`\gamma _dt1`$, Eq. (18) becomes effectively equivalent to

$$|\rho (\phi ,\phi ^{};t)|=(2\pi )^1\mathrm{exp}[(\phi \phi ^{})^2\gamma _dt],$$
(20)
and hence $`\gamma _d`$ is identified as the decoherence rate. It is possible to prove that the reduced density matrix (20) is a mixture of pure state Gaussians $`\mathrm{exp}[2(\phi \theta )^2\gamma _dt]`$ with $`\theta `$ uniformly distributed between 0 and $`2\pi `$. Thus we assert that the density matrix (20) describes a quantum system in which the variable $`\phi `$ is measured with precision $`(\gamma _dt)^{1/2}`$ with the average result chosen among a menu of uniformly distributed possible values. Here the role of the “measurement apparatus” is played by the quasiparticle field with its many degrees of freedom.

At low temperatures, we obtain $`\gamma _d/\gamma _r=n_0`$. Now we note that $`2\gamma _r`$ is the relaxation of the total energy $`E_cn_0^2/2`$ and that $`E_cn_0`$ is the variation in energy due to the exchange of a single atom. We conclude that the ratio between the decoherence and the energy relaxation rates is the ratio between the total energy and the energy variation associated to the transfer of one atom. Thus decoherence is achieved as soon as $`|\overline{n}(t)|`$ decreases by a few units comment12 . Once this has happened, Eq. (20) applies and the system may be treated semiclassically with an initial value of the phase that is chosen at random and which is defined with increasing precision as time evolves kohl01 ; comment3 . In other words, the decoherence rate can be identified with the rate of net atom exchange. On the contrary, an appreciable energy relaxation requires the exchange of many atoms. This picture is consistent with the idea that the net exchange of a single atom plays the role of the excitation of one quantum of the dissipative bath cald85 .

Once we can describe our system as a mixture of semiclassical wave packets, linear evolution applies to these even in the presence of dissipation. Since the translational motion of a Gaussian wave packet is semiclassical, we apply the method initially described with the confidence that it has been justified by a full quantum treatment.

The final picture of the connection process is one which involves several time scales. Decoherence in the phase representation sets in after the exchange of a few atoms has been completed. For a given experimental run this can be interpreted as the measurement of the phase in which a given value is chosen with an accuracy that improves with time. After the phase is defined, a semiclassical treatment suffices to conclude that the macroscopic degree of freedom will undergo very fast Bloch oscillations that will dissipate energy into the many degrees of freedom of the quasiparticle field. Ultimately this Bloch oscillations become low energy Josephson plasma oscillations that are also damped by the quasiparticle field. At long times we expect that the system density matrix will be the thermal equilibrium one which corresponds to a total energy of $`E_cn_0^2/2`$. For large $`V`$ the corresponding temperature is much smaller than the critical temperature.

The main message is that the depletion cloud, described by the quasiparticle field, can be an intrinsic agent for the decoherence and measurement of the macroscopic phase. We have studied the most fundamental intrinsic mechanism. Incoherent exchange of thermal atoms zapa98  and coupling to the walls of a cavity jaks01  will yield additional contributions to both the decoherence and the energy relaxation rates.

This work has been supported by NSF Grant no. DMR-99-86199, by the Dirección General de Investigación Científica y Técnica under Grants No. PB96-0080-C02 and BFM2001-0172, and by the Ramón Areces Foundation. One of us (I.Z.) acknowledges support from the Fulbright Commission and Spain’s Ministerio de Educación y Ciencia.


format: text

text:
# Using the Continuous Spectrum to “Feel” Integrability: the Effect of Boundary Conditions

## Abstract

The scope of this work is to propose a method for testing the integrability of a model partial differential (PDE) and/or differential difference equation (DDE). For monoparametric families of PDE/DDE’s, that are known to possess isolated integrable points, we find that very special features occur in the continuous (“phonon”) spectrum at these “singular” points. We identify these features in the case example of a PDE and a DDE (that sustain front and pulse-like solutions respectively) for different types of boundary conditions. The key finding of the work is that such spectral features are generic near the singular, integrable points and hence we propose to explore a given PDE/DDE for such traits, as a means of assessing its potential integrability.

Integrable models of partial differential (PDE) and differential difference (DDE) equations have been a topic of intense investigation over the past few decades . The main reason for that, except for the wide variety of physical applications that can be described by integrable or near-integrable systems, is that the special case of integrable models can be analyzed completely by means of the inverse scattering transform . This can then serve as a starting point for perturbative treatment of near-integrable systems.

In the process of these developments, a number of techniques have been developed for assessing integrability in continuous  or discrete  settings (or applicable to both ). An interesting feature of these “tests” is that they are necessary (but not sufficient) conditions for integrability. Hence, if a model equation fails such a criterion, it is non-integrable, but if it passes, it may or may not be integrable. In a sense, this suggests that we still do not understand the essential ingredients that render a system completely integrable. Of course, should a Lax pair be identified and the inverse scattering mechanism be applied, we know that the system is integrable, but it would certainly be desirable (as is clear from all the above effort to create “integrability tests”) to have a mechanistic (“black box”) type of criterion to assess that.

We, of course, do not claim to be providing a full answer to this question in the present work. However, we will try to give a number of useful hints that may lead to partial answers to the above questions and may provide some new intuition in the effort to construct such mechanistic criteria.

Our tool of choice will be the use of different sets of boundary conditions (BC) to “feel” the continuous spectrum (CS) of the linearization around the nonlinear coherent structure that the PDE/DDE of interest supports. Notice that the effect of boundary conditions in related contexts has been studied in a number of references; see e.g.,  and references therein. However, in all of these works the effects of the BC to the point spectrum (rather than the CS) were assessed and moreover this was not done in direct connection with issues of integrability.

In the present work, we focus on two model problems, to establish our findings and demonstrate their generality. The models are selected as one parameter families of equations such that one member of the family is an integrable system. Moreover, in illustrating the generality of the conclusions, they are selected in a form such that the one model corresponds to a PDE, while the other to a DDE and so that the one is kink bearing, while the other is pulse bearing. The models of interest will be the parametrically modified sine-Gordon equation (often also called the Peyrard-Remoissenet (PR) model)  and a modified version of the discrete nonlinear Schrödinger (DNLS) model (occasionally called the Salerno model) . The former PDE reads:

$`\varphi _{tt}\varphi _{xx}`$ $`=`$ $`{\displaystyle \frac{dU}{d\varphi }},U(\varphi ,r)={\displaystyle \frac{(1r)^2[1\mathrm{cos}(\varphi )]}{1+r^2+2r\mathrm{cos}(\varphi )}},|r|<1,|x|<\mathrm{},`$ (1)
while the latter DDE is of the form:

$`i\dot{u}_n=\mathrm{\Delta }_2u_n|u_n|^2\left[2ϵu_n+(1ϵ)(u_{n+1}+u_{n1})\right].`$ (2)
The most well-known among these monoparametric families of models are the sine-Gordon equation (Eq. (1), for $`r=0`$) which is relevant to superconductivity and charge density waves among other applications  and the experimentally realizable discrete nonlinear Schrödinger equation  of $`ϵ=1`$, as well as its integrable, so-called Ablowitz-Ladik  counterpart for $`ϵ=0`$ in the case of Eq. (2).

Notice that for the PDE, the subscripts denote partial derivatives of the field, while for the DDE, the overdot denotes temporal derivative, $`\mathrm{\Delta }_2u_nC(u_{n+1}2u_n+u_{n1})`$ where $`C=1/(\mathrm{\Delta }x)^2`$ is a constant determined by the lattice spacing $`\mathrm{\Delta }x`$; the subscript $`n`$ denotes the lattice site index. In the former case, there exist kink-like solutions which have been detailed in , while in the latter, the field is complex and there exist pulse-like solutions of the form $`u_n=\mathrm{exp}(i\mathrm{\Lambda }t)v_n`$, where $`\mathrm{\Lambda }`$ is the frequency of the solutions and $`v_n`$ its (real) exponentially localized spatial profile .

In the PDE, linearization around a state $`\varphi _0(x)`$, using the ansatz $`\varphi _0(x)+\delta \mathrm{exp}(i\omega t)f(x)`$ into Eq. (1), yields to $`O(\delta )`$ the linearization equation

$`f_{xx}+\left(\omega ^2{\displaystyle \frac{(r^21)^2[(r^2+1)\mathrm{cos}[\varphi _0(x)]+r(3\mathrm{cos}[2\varphi _0(x)])]}{[1+r^2+2r\mathrm{cos}(\varphi _0(x))]^3}}\right)f`$ $`=`$ $`0,|x|<\mathrm{}.`$ (3)
Notice that when $`r=0`$, $`\varphi _0(x)=4\mathrm{arctan}[\mathrm{exp}(x)]`$ (static kink solution of the sG equation) and for this function, the Sturm-Liouville problem (3) can be exactly solved  yielding one discrete mode (Goldstone mode) at $`\omega =0`$ and the continuous spectrum represented by the phonons. For $`r0`$, neither the static solution nor the linearization spectrum are explicitly available. Hence, to investigate the linear spectrum of Eq. (3), we should find its numerical solution, by discretizing the equation in a numerical mesh for a finite domain. The mesh consists of the $`N+1`$ points $`x_j=\{L/2+j\mathrm{\Delta }x;j=0,1,2,\mathrm{},N\}`$ defined in the finite length $`L`$ of the system ($`\mathrm{\Delta }x=L/N`$). Notice that since, in this case, we wish to emulate the behavior of the PDE, $`\mathrm{\Delta }x`$ is very fine (typically $`0.05`$), and the robustness of the findings upon variation of the (small) $`\mathrm{\Delta }x`$ has been verified.

When we compute the solution either of the PDE or of the linearization equation, we consider three different types of BC. Free BC are expressed through $`f(0)=f(1)`$ and $`f(N)=f(N1)`$ (notice that $`f(j)=f(x_j)`$). Fixed BC are given by $`f(0)=0`$ and $`f(N)=0`$, while (anti-)periodic BC are given by $`f(0)=f(N1)`$ and $`f(N)=f(1)`$.

Analogously to the PDE, for the linear stalibity analysis of DDE (2) we insert $`\mathrm{exp}(i\mathrm{\Lambda }t)[v_n+\delta (U_n\mathrm{e}^{i\omega t}+W_n\mathrm{e}^{i\omega ^{}t})]`$ into Eq. (2). We thus obtain to order $`O(\delta )`$ the following eigenvalue problem for $`\{\omega ,\{U_n,W_n^{}\}\}`$

$`\omega \left(\begin{array}{c}U_n\hfill \\ W_n^{}\hfill \end{array}\right)`$ $`=`$ $`\left(\begin{array}{c}U_n\hfill \\ W_n^{}\hfill \end{array}\right),=\left(\begin{array}{cc}A& B\\ & \\ B& A\end{array}\right),`$ (10)
$`A_{mn}`$ $`=`$ $`[\mathrm{\Lambda }+2C(4ϵv_n^2+(1ϵ)v_n[v_{n+1}+v_{n1}])]\delta _{m,n}+[(1ϵ)v_n^2C](\delta _{m,n+1}+\delta _{m,n1}),`$ (11)
$`B_{mn}`$ $`=`$ $`v_n[2ϵv_n+(1ϵ)(v_{n+1}+v_{n1})]\delta _{m,n},`$ (12)
where the stars denote complex conjugation. We compute the solutions of DDE (2) and Eq. (10) using $`200`$ points, $`\mathrm{\Delta }x=0.75`$. The BC are defined analogously through $`U_0=U_1`$, $`U_N=U_{N1}`$, $`W_0=W_1`$, and $`W_N=W_{N1}`$ for free BC. For fixed BC: $`U_0=0`$, $`U_N=0`$, $`W_0=0`$, and $`W_N=0`$, while for periodic BC: $`U_0=U_{N1}`$, $`U_N=U_1`$, $`W_0=W_{N1}`$, and $`W_N=W_1`$.

We have also tried no-flux boundary conditions that yield essentially the same results as in the free BC case.

Our results when the parameter of the PR potential or $`ϵ`$ in the DDE are varied can be summarized in Figs. 1-3.

From the above results the following conclusions can be drawn:
1. For fixed BC, the CS band edge frequency is prohibited. Hence, we compare $`\omega _k^{free}`$ with $`\omega _{k1}^{fixed}`$. We find that for small wavenumbers, fixed and free BC eigenfrequencies practically coincide only in the integrable case, whereas for the non-integrable case we observe an oscillatory behaviour of this function.
2. For periodic BC, the CS comprises of modes coming alternately from the free and fixed BC. This seems natural as the free boundary conditions select eigenmodes symmetric at the boundary, the fixed ones select modes antisymmetric at the boundary, while the periodic BC allow for both.
3. An additional feature, equally important as 1. (especially in view of its potential predictive power) is the fact that for the integrable case of $`r=0`$, periodic BC essentially imply the presence of double eigenvalues. The difference between the two eigenvalues is of the $`O(10^9)`$ for all pairs (except for the cutoff, discretization induced phenomena at the upper end of the spectrum). This is in sharp contrast (in particular for small wavenumbers), to even mild breakings of integrability, as can be inferred from Fig. 2.
4. Statements 1. and 3. above can be used in predictive form and constitute the criterion (algorithm) set forth in this work: for a given PDE/DDE model, find the steady state coherent structure (e.g., by finding the exact solution of an ODE, or numerically performing a Newton type algorithm). Linearize around the exact solution and study, in particular, the small wavenumbers, close to the lower edge of the CS (we assume that the problem is monoparametric in what follows, but it is clear that the application of the criterion does not require that). If for a critical/singular value of the parameter the fixed BC and free BC (small $`k`$) eigenvalue spectra essentially coincide and the multiplicity of periodic BC eigenvalues becomes double, then the model for this unique value of the parameter can be “strongly suspected” to be integrable. We use the above expression, as we provide no rigorous proof, but only supporting (but rather universal in distinct models with distinct features/solutions) numerical evidence for this statement.
5. We now attempt to give a qualitative explanation for the criterion set forth. One can solve explicitly (by means of Chebyshev polynomials) the discretized linearization problem of the PR model, to find: $`\omega _{k+1}^{free}=\sqrt{\omega _{min}^2+(4/(\mathrm{\Delta }x)^2)\mathrm{sin}^2(k\pi /(2N2))},k=0,1,\mathrm{},N2`$; for fixed BC $`\omega _k^{fixed}=\sqrt{\omega _{min}^2+(4/(\mathrm{\Delta }x)^2)\mathrm{sin}^2(k\pi /(2N))},k=1,2,\mathrm{},N1`$; for periodic BC $`\omega _0=\omega _{min},\omega _k=\sqrt{\omega _{min}^2+(4/(\mathrm{\Delta }x)^2)\mathrm{sin}^2(k\pi /(N1))},\omega _{N1k}=\omega _k,k=1,\mathrm{},(N2)/2,`$ respectively. From the above expressions, $`\omega _k^{free}\omega _{k1}^{fixed}`$ for small $`k`$. The difference between these two sets of frequencies as a function of $`k`$ practically does not depend on $`r`$ and coincides with the spectrum for the nonlinear integrable case $`r=0`$, except for the oscillatory behavior at the edge. For periodic BC, except for the first phonon frequency, $`\omega _{k+1}=\omega _{N2k}`$, all the eigenfrequencies are of algebraic multiplicity $`2`$. In both situations (for fixed-free and for periodic BC), the linearization around the nonlinear wave (kink) shows the exact same features as the linearization around the uniform steady state only in the integrable case. That is because, it is only in that case, that the potential is reflectionless  and, hence, there is no reflection in the scattering from the nonlinear wave. Such a reflection would lead to an eigenfrequency adjustment in the non-integrable case (in order to satisfy the boundary value eigenvalue problem). However, it is only for the integrable case that the CS (used here to “perceive” the presence of the nonlinear wave through the scattering) will not “feel” the wave (as only in this singular case there is no reflection from it).
7. We have also tested the validity of these results also in equation (2), in the vicinity of the integrable limit $`ϵ=0`$, with similar conclusions. For brevity, we only show the case of periodic BC. In Fig. 3, it can be clearly seen that it is only for the integrable case that the double eigenvalue multiplicity is obtained.

In conclusion, we have proposed and used a new test for revealing the potential integrable nature of a given model problem. By varying the boundary conditions and “feeling” the effects of such variations in the CS, we have revealed that the small wavenumbers (particularly, but also the CS more generally) have singular ways of responding to the unique parameter values for which the model is integrable. These singular features (such as an approximate identification of fixed with free BC, small $`k`$ eigenvalues and the double multiplicity of eigenvalues for periodic BC) can be used to identify and single out the integrable behavior. We have provided two model examples respectively for kinks and pulses and for a PDE and a DDE. Independently of the detailed structure of the model these properties have been identified as universal. It is naturally interesting to explore the potential usefulness of such a criterion in various more complex settings. It would also be of value to rigorously elucidate the mathematical structure that underlies the findings presented herein (we have only qualitatively attempted to justify this structure here).


text:
# Perturbative QCD analysis of the nucleon’s Pauli form factor 𝐹₂(𝑄²)

## Abstract

We perform a perturbative QCD analysis of the nucleon’s Pauli form factor $`F_2(Q^2)`$ in the asymptotically large $`Q^2`$ limit. We find that the leading contribution to $`F_2(Q^2)`$ has a $`1/Q^6`$ power behavior, consistent with the well-known result in the literature. Its coefficient depends on the leading- and subleading-twist light-cone wave functions of the nucleon, the latter describing the quarks with one unit of orbital angular momentum. We also derive at the logarithmic accurary the asymptotic scaling $`F_2(Q^2)/F_1(Q^2)(\mathrm{log}^2Q^2/\mathrm{\Lambda }^2)/Q^2`$ which describes recent Jefferson Lab data well.

preprint: DOE/ER/40762-271 UM-PP#03-031

The electromagnetic form factors are fundamental observables of the nucleon containing important information about its internal nonperturbative structure. Since the first measurement in the mid 1950s, experimental studies of these observables have become an active frontier in nuclear and particle physics. More recently, with the development of novel experimental techniques, the magnetic form factor of the proton, and electric and magnetic form factors of the neutron have been measured with unprecedented precision at Jefferson Lab and other facilities around the world Gao03 . In particular, the Pauli form factor of the proton $`F_2`$ extracted from the recoil polarizations at Jefferson Lab reveals a significant difference from previous data and theoretical expectations Jon00 ; Gay02 . The data shows that the ratio of Dirac to Pauli form factor $`Q^2F_2(Q^2)/F_1(Q^2)`$ continues climbing at the largest $`Q^2`$’s measured and seems to scale as $`\sqrt{Q^2}`$. The result has sparked myriad speculations on its implication about the underlying microscopic structure of the proton BraLenMahSte01 ; Mil02 ; RalJai02 ; MaQinSch02 ; Bro02 .

A simple dimensional counting rule was devised in Refs. BroFar73 ; MatMurTav73  to determine the dominant power contribution to hadronic form factors at large momentum transfer. For the hadron helicity-conserving form factor $`F_1(Q^2)`$, it predicts a dominant scaling behavior $`1/Q^4`$. The power counting can be justified by QCD factorization theorems which separate short-distance quark-gluon interactions from soft hadron wave functions LepBro80 ; EfrRad80 ; CheZhi77 ; DunMue79 ; AznEsaAtsTer80 ; AvdKorChe81 ; CarGro87 ; Ste99 , see Fig. 1. Recently, a more sophisticated formalism has been developed to include higher-order perturbative QCD (PQCD) resummation, or Sudakov form factor, which treats contributions from the small-$`x`$ partons more appropriately DunMue79 ; Rad84 ; IsgLle89 ; LiSte92 ; Li93 . The procedure extends the applicability of leading-order PQCD predictions to moderate $`Q^2`$, but does not change the dominant power-law behavior significantly BolJacKroBerSte95 ; Li93 ; Ste99 .

The literature on PQCD studies of $`F_2(Q^2)`$, however, is much meager. Since $`F_2(Q^2)`$ is related to the hadron helicity-flip amplitude, its power behavior is suppressed compared to $`F_1(Q^2)1/Q^4`$, and a natural expectation is then $`1/Q^6`$. This is confirmed by calculations where the quark masses are the mechanism for the helicity flip BroFar73 ; JiSil86 . Since the up and down current quark masses are negligible, the dominant mechanism for the spin flip in QCD come from the quark orbital angular momentum and the polarization of an extra gluon which must be included to maintain gauge invariance BroDre80 . A generalized power counting including the parton orbital angular momentum validates the expected scaling law for $`F_2(Q^2)`$ JiMaYua03 . An actual PQCD calculation for $`F_2(Q^2)`$ requires the extension of the usual formalism to include the quark orbital angular momentum, a technology has not yet been systematically developed in the literature. The recent development in classification of the light-cone wave functions BurJiYua02 ; JiMaYua02  provides the necessary ingredient to perform these calculations.

In this paper we report on a first PQCD calculation of $`F_2(Q^2)`$ at asymptotically large momentum transfer $`Q^2`$. Modulo logarithms, we confirm that the leading contribution to $`F_2(Q^2)`$ goes like $`1/Q^6`$. The coefficient of the leading-power term depends on leading order (twist-three) and next-to-leading order (twist-four) light-cone wave functions. The latter is the probability amplitude for one of the quarks to carry one unit of orbital angular momentum. We compute the leading-order perturbative kernel, and estimate the coefficient using models for the light-cone wave functions. We comment on the role of the Sudakov form factor in regulating possible end-point singularities in the phase space integrals. We also derive with the logarithmic accuracy the asymptotic scaling $`F_2(Q^2)/F_1(Q^2)(\mathrm{log}^2Q^2/\mathrm{\Lambda }^2)/Q^2`$ which we find to describe the data unexpectedly well.

To start with, we choose a coordinate system for the process $`P+\gamma ^{}P^{}`$ such that $`P`$ and $`P^{}`$ form conjugate light-cone vectors. We let the light-cone components $`P^+`$ and $`P_{}^{}{}_{}{}^{}`$ to be large in the asymptotic limit, $`Q^2=2P^+P_{}^{}{}_{}{}^{}`$. The $`F_2(Q^2)`$ form factor can be extracted from the helicity-flip matrix element

$$P_{}^{}|J^\mu |P_{}=F_2(Q^2)\overline{u}_{}(P^{})\frac{i\sigma ^{\mu \alpha }q_\alpha }{2M}u_{}(P),$$
(1)
where $`J^\mu `$ is the quark electromagnetic current and $`u(P)`$ is an on-shell spinor of the nucleon with momentum $`P`$. In the above coordinate system, the spin of the nucleon actually remains the same in the initial and final states.

A perturbative analysis of the above matrix element is done by computing leading-order Feynman diagrams with two gluon exchanges, for a typical graph see Fig. 1, in which the initial (final) nucleon contains three quarks with longitudinal momenta $`x_iP`$ ($`y_iP^{}`$) and transverse momenta $`𝒌_i`$ ($`𝒌_i^{}`$) of order $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. Two hard-gluon exchanges ensure that the three quarks in the final state propagate collinearly after the injection of a large momentum transfer $`q^\mu `$. Since up and down quarks are light, quark helicity is approximately conserved during the hard scattering. Therefore, to produce a nucleon-spin flip, the quark orbital angular momentum in the initial and final states must differ by one unit. The leading contribution to $`F_2(Q^2)`$ comes from the configurations in which the quarks in the initial or the final state carry zero unit, and those in the other state carry one unit, of orbital angular momentum. For definiteness, let us assume that the final state quarks are symmetric in azimuth.

To isolate the leading contribution, we expand the hard part of the diagram in $`𝒌_i^2/Q^2`$ in the limit of large $`Q^2`$ and $`x_i0`$ (a procedure dubbed the collinear expansion). We will comment on the $`x_i0`$ case later. Since the final state nucleon has no orbital angular momentum, we throw away all the subleading terms in $`𝒌_i^{}`$. On the other hand, the quarks in the initial state nucleon have one unit of orbital angular momentum, and hence the hard part must have linear terms in the quark transverse momenta. Therefore, the leading hard part, we are interested in, has a structure $`𝒌_iT(x_1,x_2,x_3,y_1,y_2,y_3;Q^2)`$.

Let us consider the simplification of the nucleon wave functions after the collinear expansion. For definiteness, we consider the proton form factor. The light-cone wave function for the final state is

$`|P_{}_{1/2}={\displaystyle \frac{1}{12}}{\displaystyle \frac{[dx][d^2𝒌]}{\sqrt{x_1x_2x_3}}\psi _1(\kappa _1,\kappa _2,\kappa _3)}`$ (2)
$`\times \epsilon ^{abc}u_a^{}(\kappa _1)\left\{u_b^{}(\kappa _2)d_c^{}(\kappa _3)d_b^{}(\kappa _2)u_c^{}(\kappa _3)\right\}|0,`$
where the argument $`\kappa _i`$ is a shorthand notation for $`(x_i,𝒌_i)`$, and the integration measures for the quark momenta are

$`[dx]dx_1dx_2dx_3\delta \left(x_1+x_2+x_31\right),`$
$`[d^2𝒌]d^2𝒌_1d^2𝒌_2d^2𝒌_3\delta ^{(2)}\left(𝒌_1+𝒌_2+𝒌_3\right).`$
If the hard part has no dependence on $`𝒌_i^{}`$, we can ignore the transverse momentum dependence in the quark creation operators and define the twist-three amplitude

$`\mathrm{\Phi }_3(x_1,x_2,x_3)=2{\displaystyle [d^2𝒌]\psi _1(\kappa _1,\kappa _2,\kappa _3)}.`$
Then the final nucleon state can be simplified to,

$`|P_{}_{1/2}={\displaystyle \frac{1}{24}}{\displaystyle \frac{[dx]}{\sqrt{x_1x_2x_3}}\mathrm{\Phi }_3(x_1,x_2,x_3)}`$ (3)
$`\times \epsilon ^{abc}u_a^{}(x_1)\left\{u_b^{}(x_2)d_c^{}(x_3)d_b^{}(x_2)u_c^{}(x_3)\right\}|0.`$
In the first quantization formalism, the second line can be replaced by the standard $`SU(6)`$ wave function.

Since the quarks in the initial nucleon state must have a total helicity $`1/2`$, only the following wave function component is relevant JiMaYua02 ; BurJiYua02

$`|P_{}_{1/2}={\displaystyle \frac{1}{12}}{\displaystyle \frac{[dx][d^2𝒌]}{\sqrt{x_1x_2x_3}}\left\{\overline{k}_1\psi _3+\overline{k}_2\psi _4\right\}(\kappa _1,\kappa _2,\kappa _3)}`$
$`\times \epsilon ^{abc}u_a^{}(\kappa _2)\left\{u_b^{}(\kappa _1)d_c^{}(\kappa _3)d_b^{}(\kappa _1)u_c^{}(\kappa _3)\right\}|0,`$ (4)
where we use the notation $`\overline{k}_{}k^xik^y`$.

The spin-isospin structure of this wave function is exactly the same as that in Eq. (Perturbative QCD analysis of the nucleon’s Pauli form factor $`F_2(Q^2)`$). After the collinear expansion, the above wave function will be convoluted with transverse momenta of quarks. We define the twist-four amplitudes via

$`\mathrm{\Phi }_4(x_2,x_1,x_3)`$ $`=`$ $`2{\displaystyle \frac{[d^2𝒌]}{Mx_3}𝒌_3\left\{𝒌_1\psi _3+𝒌_2\psi _4\right\}(\kappa _1,\kappa _2,\kappa _3)},`$
$`\mathrm{\Psi }_4(x_1,x_2,x_3)`$ $`=`$ $`2{\displaystyle \frac{[d^2𝒌]}{Mx_2}𝒌_2\left\{𝒌_1\psi _3+𝒌_2\psi _4\right\}(\kappa _1,\kappa _2,\kappa _3)}.`$
Apart from a gluon-potential dependent term, the above light-cone amplitudes are the same as those defined in Ref. BraFriMahSte00 . In the following, we include the former so that the result is gauge invariant. To the order that we are working, the gluon potential terms arise from perturbative diagrams with explicit gluons attached to the nucleon blobs. These contributions are generally considered as dynamically suppressed Pol80 .

We introduce effective nucleon wave functions for the initial state after integrating over the transverse momenta weighted with a momentum factor from the hard part,

$`|P_{}_{1/2}[𝒌_1]={\displaystyle \frac{M}{24}}{\displaystyle \frac{[dx]}{\sqrt{x_1x_2x_3}}x_1\mathrm{\Psi }_4(x_2,x_1,x_3)}`$ (5)
$`\times ϵ^{abc}u_a^{}(x_1)\left\{u_b^{}(x_2)d_c^{}(x_3)d_b^{}(x_2)u_c^{}(x_3)\right\}|0,`$
$`|P_{}_{1/2}[𝒌_3]={\displaystyle \frac{M}{24}}{\displaystyle \frac{[dx]}{\sqrt{x_1x_2x_3}}x_3\mathrm{\Phi }_4(x_1,x_2,x_3)}`$ (6)
$`\times ϵ^{abc}u_a^{}(x_1)\left\{u_b^{}(x_2)d_c^{}(x_3)d_b^{}(x_2)u_c^{}(x_3)\right\}|0,`$
where the argument on the left-hand side indicates the momentum being averaged. Averaging over $`𝒌_2`$ does not yield independent information.

Using the twist-3 and -4 amplitudes, we can write the $`F_2(Q^2)`$ form factor in the following factorized form,

$`F_2(Q^2)`$ $`=`$ $`{\displaystyle }[dx][dy]\{x_3\mathrm{\Phi }_4(x_1,x_2,x_3)T_\mathrm{\Phi }(\{x\},\{y\})`$
$`+`$ $`x_1\mathrm{\Psi }_4(x_2,x_1,x_3)T_\mathrm{\Psi }(\{x\},\{y\})\}\mathrm{\Phi }_3(y_1,y_2,y_3),`$
where $`\{x\}=(x_1,x_2,x_3)`$. Let us see how to extract the hard part, $`T`$, in Fig. 1.

We use the nomenclature and strategy of Ref. LepBro80  by calling the top quark line 1, the middle 2, etc., without committing them to a specific flavor. Given the spin-isospin structure of the initial and final state nucleon wave functions, we assume that the first and third quarks have spin up and the second one spin down. Call the amplitudes $`T_i`$ when the electromagnetic current acts on particle $`i`$. Clearly because of symmetry, $`T_3`$ can be obtained from $`T_1`$ by a suitable exchange of variables (see below). Using the $`SU(6)`$ wave function and the quark charge weighting, we find the flavor structure of the hard part for the proton,

$`T^p(\{x\},\{y\})`$ $`=`$ $`{\displaystyle \frac{2e_u}{3}}T_1+{\displaystyle \frac{e_u+e_d}{3}}(T_2+T_3)`$ (8)
$`+`$ $`{\displaystyle \frac{e_u}{3}}(T_1^{}+T_3^{})+{\displaystyle \frac{e_d}{3}}T_2^{}.`$
where we have omitted the argument of $`T_i`$ on the right-hand side, and $`T_i^{}`$ has $`y_1`$ and $`y_3`$ interchanged. For the neutron, $`ud`$.

According to the above convention, we find the contribution to $`T_1`$ from Fig. 1,

$`T_{1\mathrm{\Psi }}^{\mathrm{fig}.1}(\{x\},\{y\})`$ $`=`$ $`{\displaystyle \frac{M^2C_B^2}{Q^6}}(4\pi \alpha _s)^2{\displaystyle \frac{1}{\overline{x}_1^2x_3\overline{y}_1y_3^2}},`$
$`T_{1\mathrm{\Phi }}^{\mathrm{fig}.1}(\{x\},\{y\})`$ $`=`$ $`{\displaystyle \frac{M^2C_B^2}{Q^6}}(4\pi \alpha _s)^2{\displaystyle \frac{1}{\overline{x}_1x_3^2\overline{y}_1y_3^2}},`$ (9)
where here and everywhere $`\overline{x}1x`$. $`C_B=2/3`$ is a color factor and the contribution from the twist-4 wave functions of the final-state nucleon has also been included. The corresponding contribution to $`T_3`$ can be obtained by interchanging the labels, $`T_{3\mathrm{\Psi }}^{\mathrm{fig}.1}=T_{1\mathrm{\Phi }}^{\mathrm{fig}.1}(13)`$, $`T_{3\mathrm{\Phi }}^{\mathrm{fig}.1}=T_{1\mathrm{\Psi }}^{\mathrm{fig}.1}(13)`$. This in fact is a general feature.

To find a complete expression for the hard part, we must calculate perturbative diagrams displayed in Fig. 2. A straightforward evaluation yields,

$`T_{1\mathrm{\Psi }}(\{x\},\{y\})`$ $`=`$ $`{\displaystyle \frac{M^2C_B^2}{Q^6}}(4\pi \alpha _s)^2x_3(G_{11}G_{12}),`$ (10)
$`T_{1\mathrm{\Phi }}(\{x\},\{y\})`$ $`=`$ $`{\displaystyle \frac{M^2C_B^2}{Q^6}}(4\pi \alpha _s)^2[x_1G_{11}\overline{x}_1G_{12}],`$
$`T_{2\mathrm{\Psi }}(\{x\},\{y\})`$ $`=`$ $`{\displaystyle \frac{M^2C_B^2}{Q^6}}(4\pi \alpha _s)^2`$
$`\times `$ $`[x_3(G_{22}\stackrel{~}{G}_{21}\stackrel{~}{G}_{22})\overline{x}_3G_{21}],`$
and $`T_{2\mathrm{\Phi }}=T_{2\mathrm{\Psi }}(13)`$, $`T_{3\mathrm{\Psi }}=T_{1\mathrm{\Phi }}(13)`$, and $`T_{3\mathrm{\Phi }}=T_{1\mathrm{\Psi }}(13)`$. The functions $`G_{ij}`$ are defined as

$`G_{11}`$ $`=`$ $`{\displaystyle \frac{1}{x_1x_2x_3^2y_2y_3^2\overline{y}_3}},`$
$`G_{12}`$ $`=`$ $`{\displaystyle \frac{1}{x_3^2\overline{x}_1^2y_3\overline{y}_1^2}}+{\displaystyle \frac{1}{x_2x_3\overline{x}_1^2y_2\overline{y}_1^2}}`$
$`+{\displaystyle \frac{1}{\overline{x}_1^2x_3^2\overline{y}_1y_3^2}}{\displaystyle \frac{1}{x_2x_3^2\overline{x}_1y_2y_3\overline{y}_3}},`$
$`G_{21}`$ $`=`$ $`{\displaystyle \frac{1}{x_1^2x_3\overline{x}_3y_1^2y_3\overline{y}_1}},`$
$`G_{22}`$ $`=`$ $`{\displaystyle \frac{1}{x_1x_3^2\overline{x}_2y_3^2\overline{y}_2}}+{\displaystyle \frac{1}{x_1^2x_3^2y_1y_3\overline{y}_1}},`$ (11)
and $`\stackrel{~}{G}_{21}=G_{21}(13)`$, $`\stackrel{~}{G}_{22}=G_{22}(13)`$.

To determine the normalization for $`F_2(Q^2)`$, we need to know the light-cone distribution amplitudes $`\mathrm{\Phi }_3`$, $`\mathrm{\Phi }_4`$ and $`\mathrm{\Psi }_4`$, which can only be obtained by solving QCD nonperturbatively. However, the scale evolution of these amplitudes selects at the asymptotically large $`Q^2`$ the leading component with a fixed small-$`x_i`$ behavior. For example, the asymptotic form of $`\mathrm{\Phi }_3`$ is $`x_1x_2x_3`$ LepBro80 , whereas that of $`\mathrm{\Phi }_4`$ is $`x_1x_2`$. In Ref. BraFriMahSte00  a set of phenomenological amplitudes satisfying these asymptotic constraints have been proposed on a basis of conformal expansion.

With the above wave functions, the integrals over momentum fractions $`x_i`$ and $`y_i`$ have logarithmic singularities, indicating that the factorization breaks down when one of the quarks in the wave function becomes soft Rad84 ; IsgLle89 . It has been suggested that the higher-order PQCD resummation, or the Sudakov form factor, suppresses the contribution at small-$`x`$ and provides an effective cut-off for the integrals at $`x\mathrm{\Lambda }^2/Q^2`$, where $`\mathrm{\Lambda }`$ is a soft scale related to the size of the nucleon Li93 ; LiSte92 ; BolJacKroBerSte95 ; MusRad97 ; DesSac01 . The outcome is that the $`x_i`$ integrations contribute an extra $`Q^2`$-dependent factor $`\mathrm{ln}^2Q^2/\mathrm{\Lambda }^2`$, compared to $`F_1(Q^2)`$. Physically, the end-point divergencies indicate that quarks with different rapidity contribute equally to the hard scattering. Since the contribution from quarks with very large rapidity (small-$`x`$) is suppressed by the Sudakov form factor, this $`Q^2`$-dependence reflects simply the kinematic broadening of the quark (and gluon) rapidity range with increasing nucleon momentum.

For an estimate, we use asymptotic wave functions BraFriMahSte00 ,

$`\mathrm{\Phi }_3=120x_1x_2x_3f_N,\mathrm{\Phi }_4=12x_1x_2(f_N+\lambda _1),`$
$`\mathrm{\Psi }_4=12x_1x_3(f_N\lambda _1),`$ (12)
with $`f_N=5.310^3\mathrm{GeV}^2`$ and $`\lambda _1=2.710^2\mathrm{GeV}^2`$. With a choice of $`\mathrm{\Lambda }=(0.3\mathrm{GeV})^2`$, $`Q^6F_2^p(Q^2)`$ is roughly $`0.6\mathrm{GeV}^6`$ for $`Q^2=(520)\mathrm{GeV}^2`$, about $`1/3`$ of the Jefferson Lab data at $`Q^2=5\mathrm{GeV}^2`$ Gay02 . Of course, to get a more realistic PQCD prediction in this regime, one must have the quark distribution amplitudes appropriate at this scale. However, from the comparison between the data and PQCD predictions for $`F_1(Q^2)`$ Bro02 , we believe that asymptotic PQCD is unlikely to be the dominant contribution to $`F_2(Q^2)`$ at $`Q^2=35\mathrm{GeV}^2`$: one must take into account higher-order corrections and higher-twist effects.

Coming back to the scaling behavior of the ratio $`F_2(Q^2)/F_1(Q^2)`$ for which the Jefferson Lab data has stimulated much discussion in the literature. PQCD predicts the power-law scaling $`1/Q^2`$. With the new result for $`F_2(Q^2)`$, we can determine its scaling up to logarithmic accuracy. The strong coupling constant in the ratio simply cancels. The wave function evolution yields a factor of $`\alpha _s^{32/(9\beta )}(Q^2)`$ for $`F_1(Q^2)`$ and $`\alpha _s^{8/(3\beta )}(Q^2)`$ for $`F_2(Q^2)`$ from the leading non-vanishing contribution, where $`\beta =112n_f/3`$. Thus PQCD predicts that $`(Q^2/\mathrm{ln}^{2+8/(9\beta )}Q^2/\mathrm{\Lambda }^2)(F_2(Q^2)/F_1(Q^2))`$ scales as a constant at large $`Q^2`$, $`8/(9\beta )1`$. Surprisingly, the Jefferson Lab data plotted this way (ignoring the small $`8/(9\beta )`$) exhibits little $`Q^2`$ variation for a range of choices of $`\mathrm{\Lambda }`$ as shown in Fig. 3. Since we do not expect the asymptotic predictions for $`F_{1,2}(Q^2)`$ to work at these $`Q^2`$, the observed consistency might be a sign of precocious scaling as a consequence of delicate cancellations in the ratio. A more detailed discussion on this issue along with more thorough phenomenological analyses will be given in a separate publication.

We thank V.M. Braun, H. Gao, G.P. Korchemsky, H.-N. Li, and A.V. Radyushkin for useful discussions. This work was supported by the U. S. Department of Energy via grant DE-FG02-93ER-40762.


text:
# Theory of charge transport in diffusive normal metal / conventional superconductor point contacts

## I Introduction

The electron coherence in mesoscopic superconducting systems is one of the important topics of solid state physics. The low energy transport in these systems is essentially influenced by the Andreev reflection Andreev , a unique process specific for normal metal/superconductor interfaces. The phase coherence between incoming electrons and Andreev reflected holes persists in the diffusive normal metal at a mesoscopic length scale and results in strong interference effects on the probability of Andreev reflection Hekking . These effects become prominent at sufficiently low temperatures where the thermal broadening is negligible. One of the remarkable experimental manifestations is the zero bias conductance peak (ZBCP) Giazotto ; Klapwijk ; Kastalsky ; Nguyen ; Wees ; Nitta ; Bakker ; Xiong ; Magnee ; Kutch ; Poirier . A calculation of tunneling conductance in a normal metal (N) / superconductor (S) junction is an interesting theoretical problem since quantum interference effects due to Andreev reflection are expected.

For a clean NS contact in the presence of the interface potential barrier the conductance was calculated by Blonder, Tinkham and Klapwijk BTK  (BTK) in terms of the corresponding transmission coefficients on the basis of the solution of the Bogoljubov - de Gennes equations. From the general set of boundary conditions connecting the quasiclassical Green’s functions on both sides of the interface for arbitrary transmission probabilities, Zaitsev Zaitsev  derived the expression for the conductance similar to that by BTK. The BTK method BTK  is confined to ballistic systems. The generalization of this method to systems with impurities has been performed by several authors (see the review Beenakker1 ). In a number of papers the transmission coefficients were directly calculated by numerical methods Lambert ; Takane . However, it is difficult to apply such methods to most of relevant experimental situations. Another approach, the so-called random matrix theory, was employed by Beenakker $`et`$ $`al.`$, where the total transmission coefficients are expressed in terms of those through the normal part of the system and the normal/superconductor interface separately Beenakker1 ; Beenakker2 . Within this theory, the ZBCP observed in experiments is understood as a resonance phenomenon related to reflectionless tunneling reflec . The scattering matrix approach was later generalized to finite voltage and temperature.Lesovik

On the other hand, a quasiclassical Green’s function calculation based on nonequilibrium superconductivity theories Larkin  is much more powerful and convenient for the actual calculations Volkov . In this approach, the impurity scattering is included in the self-consistent Born approximation and the weak localization effects are neglected. In the theory of tunneling conductance developed by Volkov, Zaitsev and Klapwijk (VZK) by solving the Usadel equations Usadel , the origin of the ZBCP observed in several experiments was clarified to be due to the enhancement of the pair amplitude in the diffusive normal metal by the proximity effect Volkov . VZK applied the Kupriyanov and Lukichev (KL) boundary condition for the Keldysh-Nambu Green’s function KL . The KL boundary condition is valid for the atomically sharp interface barrier dividing two diffusive metals. As shown by Lambert et al. Lambert1 , this condition is exact in two limits of either high or low barrier transparency, with small corrections in the intermediate transparency range. By applying the VZK theory, several authors studied the charge transport in various junctions Nazarov1 ; Yip ; Stoof ; Reentrance ; Golubov ; Takayanagi ; Bezuglyi ; Seviour  by solving the Usadel equations. (see review by Belzig et al. Belzig ).

The generalization of the KL boundary conditions for an arbitrary connector between diffusive metals was provided by Nazarov within the so-called ”generalized circuit theory” Nazarov2 . In this theory, the mesoscopic system is presented as a network of nodes and connectors. A connector is characterized by a set of transmission coefficients and can present anything from a ballistic point contact to a tunnel junction. A conservation law of matrix current holds in each node. The method to derive the relation between matrix current and Green’s functions puts the results of Ref.Zaitsev  to the framework of Landauer-Büttiker scattering formalism. The boundary condition for Keldysh-Nambu Green’s function was derived in Nazarov2  for an arbitrary connector including various situations from ballistic point contact to diffusive contact. Actually, this boundary condition is very general since the BTK theory is reproduced in the ballistic limit while in the diffusive limit with a low transmissivity of the interface, the KL boundary condition is reproduced.

Although a number of papers were published on charge transport in mesoscopic NS junctions, as far as we know, almost all of them are either based on the KL boundary conditions or on the BTK model. However in the actual junctions, transparency of the junction is not necessarily small and impurity scattering in the DN is important. Therefore, an interesting and important theoretical problem is the calculation of the tunneling conductance in normal metal / conventional superconductor junctions using the boundary condition from Ref.Nazarov2  since both the ballistic (BTK theory) and diffusive (VZK theory) cases can be covered simultaneously. In the present paper, we study the tunneling conductance in diffusive normal metal / insulator / conventional superconductor (DN/I/S) junctions for various parameters such as the height of the insulating barrier at the interface, resistance $`R_d`$ in DN and the Thouless energy $`E_{Th}`$ in DN. We concentrate on the normalized tunneling conductance of the junctions $`\sigma _T(eV)`$ as a function of the bias voltage $`V`$. The conductance $`\sigma _T(eV)`$ is given by $`\sigma _T(eV)=\sigma _S(eV)/\sigma _N(eV)`$ where $`\sigma _{S(N)}(eV)`$ is the tunneling conductance in the superconducting (normal ) state at a bias voltage $`V`$.

In the present paper the following points are clarified:

1. When the transparency of the junction is sufficiently low, $`\sigma _T(eV)`$ for $`eV<\mathrm{\Delta }_0`$ is enhanced with the increase of $`R_d`$ due to the enhancement of the proximity effect. The ZBCP becomes prominent for $`E_{Th}<<\mathrm{\Delta }_0`$ and $`R_d/R_b<1`$. In such a case, with a further increase of $`R_d/R_b`$ the ZBCP changes into a zero bias conductance dip (ZBCD). In the low transparent limit, the line shapes of $`\sigma _T(eV)`$ are qualitatively the same as those obtained by VZK theory Volkov ; Yip .

2. When the transparency of the junction is almost unity, $`\sigma _T(eV)`$ always exhibits a ZBCD except for the special case of $`R_d=0`$, i.e. the BTK limit.

3. The measure of the proximity effect, $`\theta `$, is mainly determined by $`R_d/R_b`$ and $`E_{Th}`$, where $`R_b`$ is the resistance from the insulating barrier. The proximity effect enhances (reduces) the magnitude of $`\sigma _T(eV)`$ for junctions with low (high) transparency.

4. Even for junctions between conventional $`s`$-wave superconductors, we can expect a wide variety of line shapes of the tunneling conductance, a ZBCP, ZBCD, $`U`$-shaped structure, and a rounded bottom structure.

5. We have clarified the parameter space where a ZBCP should be expected. Typically, small Thouless energy $`E_{Th}`$ is required for a ZBCP. If the magnitude of $`E_{Th}`$ is increasing up to $`\mathrm{\Delta }_0`$, a ZBCP is only expected for junctions with low transmissivity, $`R_d/R_b<<1`$.

The organization of this paper is as follows. In section 2, we will provide the detailed derivation of the expression for the normalized tunneling conductance. In section 3, the results of calculations of $`\sigma _T(eV)`$ and $`\theta `$ are presented for various types of junctions. In section 4, the summary of the obtained results is given.

## II Formulation

In this section we introduce the model and the formalism. We consider a junction consisting of normal and superconducting reservoirs connected by a quasi-one-dimensional diffusive conductor (DN) with a length $`L`$ much larger than the mean free path. The interface between the DN conductor and the S electrode has a resistance $`R_d`$ while the DN/N interface has zero resistance. The positions of the DN/N interface of the DN/S interface are denoted as $`x=0`$ and $`x=L`$, respectively. According to the circuit theory, the interface between DN and S is subdivided into two isotropization zones in DN and S, two ballistic zones and a scattering zone. The sizes of the ballistic and scattering zones in the current flow direction are much shorter than the coherence length. Although the generalized boundary condition of Ref. Nazarov2  is valid for arbitrary interfaces, here scattering zone is modelled as an infinitely narrow insulating barrier described by the delta function $`U(x)=H\delta (xL)`$. The resulting transparency of the junctions $`T_m`$ is given by $`T_m=4\mathrm{cos}^2\varphi /(4\mathrm{cos}^2\varphi +Z^2)`$, where $`Z=2H/(\mathrm{}v_F)`$ is a dimensionless constant, $`\varphi `$ is the injection angle measured from the interface normal to the junction and $`v_F`$ is Fermi velocity. Variation of the barrier shape will not change our results in the considered case of isotropic superconductivity in the S electrode.

We apply the quasiclassical Keldysh formalism in the following calculation of the tunneling conductance. The 4 $`\times `$ 4 Green’s functions in DN and S are denoted by $`\stackrel{ˇ}{G}_1(x)`$ and $`\stackrel{ˇ}{G}_2(x)`$ which are expressed in matrix form as

$`\stackrel{ˇ}{G}_1(x)`$ $`=`$ $`\left(\begin{array}{cc}\widehat{R}_1(x)& \widehat{K}_1(x)\\ 0& \widehat{A}_1(x)\end{array}\right),`$ (3)
$`\stackrel{ˇ}{G}_2(x)`$ $`=`$ $`\left(\begin{array}{cc}\widehat{R}_2(x)& \widehat{K}_2(x)\\ 0& \widehat{A}_2(x)\end{array}\right),`$ (6)
where the Keldysh component $`\widehat{K}_{1,2}(x)`$ is given by $`\widehat{K}_{1(2)}(x)=\widehat{R}_{1(2)}(x)\widehat{f}_{1(2)}(x)\widehat{f}_{1(2)}(x)\widehat{A}_{1(2)}(x)`$ with retarded component $`\widehat{R}_{1,2}(x)`$, advanced component $`\widehat{A}_{1,2}(x)`$ using distribution function $`\widehat{f}_{1(2)}(x)`$. In the above, $`\widehat{R}_2(x)`$ is expressed by

$$\widehat{R}_2(x)=(g\widehat{\tau }_3+f\widehat{\tau }_2)$$

with $`g=ϵ/\sqrt{ϵ^2\mathrm{\Delta }_0^2}`$ and $`f=\mathrm{\Delta }_0/\sqrt{\mathrm{\Delta }_0^2ϵ^2}`$, where $`ϵ`$ denotes the quasiparticle energy measured from the Fermi energy, $`\widehat{A}_2(x)=\widehat{R}_2^{}(x)`$ and $`\widehat{f}_2(x)=\mathrm{tanh}[ϵ/(2\mathrm{k}_\mathrm{B}\mathrm{T})]`$ in thermal equilibrium with temperature $`T`$. We put the electrical potential zero in the S-electrode. In this case the spatial dependence of $`\stackrel{ˇ}{G}_1(x)`$ in DN is determined by the static Usadel equation Usadel ,

$$D\frac{}{x}[\stackrel{ˇ}{G}_1(x)\frac{\stackrel{ˇ}{G}_1(x)}{x}]+i[\stackrel{ˇ}{H},\stackrel{ˇ}{G}_1(x)]=0,$$
(7)
with the diffusion constant $`D`$ in DN, where $`\stackrel{ˇ}{H}`$ is given by

$$\stackrel{ˇ}{H}=\left(\begin{array}{cc}\widehat{H}_0& 0\\ 0& \widehat{H}_0\end{array}\right),$$

with $`\widehat{H}_0=ϵ\tau _3.`$

The boundary condition for $`\stackrel{ˇ}{G}_1(x)`$ at the DN/S interface is given by Nazarov’s generalized boundary condition,

$$\frac{L}{R_d}(\stackrel{ˇ}{G}_1\frac{\stackrel{ˇ}{G}_1}{x})_{x=L_{}}=R_b^1<B>,$$
(8)
$$B=\frac{2T_m[\stackrel{ˇ}{G}_1(L_{}),\stackrel{ˇ}{G}_2(L_+)]}{4+T_m([\stackrel{ˇ}{G}_1(L_{}),\stackrel{ˇ}{G}_2(L_+)]_+2)}.$$

The average over the various angles of injected particles at the interface is defined as

$$<B(\varphi )>=_{\pi /2}^{\pi /2}𝑑\varphi \mathrm{cos}\varphi B(\varphi )/_{\pi /2}^{\pi /2}𝑑\varphi T(\varphi )\mathrm{cos}\varphi $$

with $`B(\varphi )=B`$ and $`T(\varphi )=T_m`$. The resistance of the interface $`R_b`$ is given by

$$R_b=R_0\frac{2}{_{\pi /2}^{\pi /2}𝑑\varphi T(\varphi )\mathrm{cos}\varphi }.$$

Here $`R_0`$ is Sharvin resistance, which in three-dimensional case is given by $`R_0^1=e^2k_F^2S_c/(4\pi ^2\mathrm{})`$, where $`k_F`$ is the Fermi wave-vector and $`S_c`$ is the constriction area. Note that the area $`S_c`$ is in general not equal to the crossection area $`S_d`$ of the normal conductor, therefore $`S_c/S_d`$ is independent parameter of our theory.

For $`T_m0`$ in Eq. (4), the quantity $`B`$ can be expressed as

$$B=\frac{T_m}{2}[\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2]$$

and we can reproduce the KL boundary condition. On the other hand, at $`x=0`$ $`\stackrel{ˇ}{G}_1(0)`$ coincides with that in the normal state.

The electric current is expressed using $`\stackrel{ˇ}{G}_1(x)`$ as

$$I_{el}=\frac{L}{4eR_d}_0^{\mathrm{}}𝑑ϵ\mathrm{Tr}[\tau _3(\stackrel{ˇ}{G}_1(x)\frac{\stackrel{ˇ}{G}_1(x)}{x})^K],$$
(9)
where $`(\stackrel{ˇ}{G_1}(x)\frac{\stackrel{ˇ}{G_1}(x)}{x})^K`$ denotes the Keldysh component of $`(\stackrel{ˇ}{G_1}(x)\frac{\stackrel{ˇ}{G_1}(x)}{x})`$. In the actual calculation it is convenient to use the standard $`\theta `$-parameterization when function $`\widehat{R}_1(x)`$ is expressed as

$$\widehat{R}_1(x)=\widehat{\tau }_3\mathrm{cos}\theta (x)+\widehat{\tau }_2\mathrm{sin}\theta (x).$$
(10)
The parameter $`\theta (x)`$ is a measure of the proximity effect in DN.

Functions $`\widehat{A}_1(x)`$ and $`\widehat{K}_1(x)`$ are expressed as $`\widehat{A}_1(x)=\widehat{R}_1^{}(x)`$ and $`\widehat{K}_1(x)=\widehat{R}_1(x)\widehat{f}_1(x)\widehat{f}_1(x)\widehat{A}_1(x)`$ with the distribution function $`\widehat{f}_1(x)`$ which is given by $`\widehat{f}_1(x)=f_l(x)+\widehat{\tau }_3f_t(x)`$. In the above, $`f_t(x)`$ is the relevant distribution function which determines the conductance of the junction we are now concentrating on. From the retarded or advanced component of the Usadel equation, the spatial dependence of $`\theta (x)`$ is determined by the following equation

$$D\frac{^2}{x^2}\theta (x)+2iϵ\mathrm{sin}[\theta (x)]=0,$$
(11)
while for the Keldysh component we obtain

$$D\frac{}{x}[\frac{f_t(x)}{x}\mathrm{cosh}^2\theta _{imag}(x)]=0.$$
(12)
At $`x=0`$, since DN is attached to the normal electrode, $`\theta (0)`$=0 and $`f_t(0)=f_{t0}`$ is satisfied with

$$f_{t0}=\frac{1}{2}\{\mathrm{tanh}[(ϵ+eV)/(2k_BT)]\mathrm{tanh}[(ϵeV)/(2k_BT)]\}.$$

Next we focus on the boundary condition at the DN/S interface. Taking the retarded part of Eq. (8), we obtain

$$\frac{L}{R_d}\frac{\theta (x)}{x}_{x=L_{}}=\frac{<F>}{R_b},$$
(13)
$$F=\frac{2(f\mathrm{cos}\theta _Lg\mathrm{sin}\theta _L)T_m}{(2T_m)+T_m[g\mathrm{cos}\theta _L+f\mathrm{sin}\theta _L]},$$

with $`\theta _L=\theta (L_{})`$.

On the other hand, from the Keldysh part of Eq. (8), we obtain

$$\frac{L}{R_d}(\frac{f_t}{x})\mathrm{cosh}^2\theta _{imag}(x)_{x=L_{}}=\frac{<I_{b0}>f_t(L_{})}{R_b},$$
(14)
with

$$I_{b0}=\frac{T_m^2\mathrm{\Lambda }_1+2T_m(2T_m)\mathrm{\Lambda }_2}{2(2T_m)+T_m[g\mathrm{cos}\theta _L+f\mathrm{sin}\theta _L]^2},$$

$$\mathrm{\Lambda }_1=(1+\mathrm{cos}\theta _L^2+\mathrm{sin}\theta _L^2)(g^2+f^2+1)$$

$$+4\mathrm{I}\mathrm{m}\mathrm{a}\mathrm{g}[fg^{}]\mathrm{Imag}[\mathrm{cos}\theta _L\mathrm{sin}\theta _L^{}],$$
(15)
$$\mathrm{\Lambda }_2=\mathrm{Real}\{\mathrm{g}(\mathrm{cos}\theta _\mathrm{L}+\mathrm{cos}\theta _\mathrm{L}^{})+\mathrm{f}(\mathrm{sin}\theta _\mathrm{L}+\mathrm{sin}\theta _\mathrm{L}^{})\},$$
(16)
where $`\theta _{imag}(x)`$ denotes the imaginary part of $`\theta (x)`$. For $`T_m<<1`$, $`I_{b0}`$ is reduced to

$$I_{b0}=\frac{\mathrm{Real}[\mathrm{g}(\mathrm{cos}\theta _\mathrm{L}+\mathrm{cos}\theta _\mathrm{L}^{})+\mathrm{f}(\mathrm{sin}\theta _\mathrm{L}+\mathrm{sin}\theta _\mathrm{L}^{})]}{2}T_m,$$
(17)
which is the expression used in the VZK theory. After a simple manipulation, we can obtain $`f_t(L_{})`$

$$f_t(L_{})=\frac{R_bf_{t0}}{R_b+\frac{R_d<I_{b0}>}{L}_0^L\frac{dx}{\mathrm{cosh}^2\theta _{imag}(x)}}$$

Since the electric current $`I_{el}`$ can be expressed via $`\theta _L`$ in the following form

$$I_{el}=\frac{L}{eR_d}_0^{\mathrm{}}(\frac{f_t}{x})_{x=L_{}}\mathrm{cosh}^2[\mathrm{Imag}(\theta _L)]dϵ,$$

we obtain the following final result for the current

$$I_{el}=\frac{1}{e}_0^{\mathrm{}}𝑑ϵ\frac{f_{t0}}{\frac{R_b}{<I_{b0}>}+\frac{R_d}{L}_0^L\frac{dx}{\mathrm{cosh}^2\theta _{imag}(x)}}.$$
(18)
Then the total resistance $`R`$ at zero temperature is given by

$$R=\frac{R_b}{<I_{b0}>}+\frac{R_d}{L}_0^L\frac{dx}{\mathrm{cosh}^2\theta _{imag}(x)}$$
(19)
and the tunneling conductance in the superconducting state $`\sigma _S(eV)`$ is given by $`\sigma _S(eV)=1/R`$.

It should be mentioned that for $`R_d=0`$, $`\theta _L`$ becomes zero due to the absence of the proximity effect. Then $`I_{b0}`$ is given as follows

$$I_{b0}=\frac{(1+g^2+f^2)T_m^2+2T_m(2T_m)\mathrm{Real}(g)}{(2T_m)+T_mg^2}$$

$$=\frac{T_m[1+\mathrm{\Gamma }^2+(T_m1)\mathrm{\Gamma }^4]}{1(1T_m)\mathrm{\Gamma }^2^2}$$
(20)
with $`\mathrm{\Gamma }=\frac{\sqrt{ϵ\sqrt{ϵ^2\mathrm{\Delta }_0^2}}}{\sqrt{ϵ+\sqrt{ϵ^2\mathrm{\Delta }_0^2}}}`$ and the resulting $`\sigma _S`$ is given by

$$\sigma _S(eV)=\frac{1}{R_0}_{\pi /2}^{\pi /2}\frac{I_{b0}}{2}\mathrm{cos}\varphi d\varphi ,$$

and reproduces that by BTK theory.

It should be remarked that in the present circuit theory, $`R_d/R_b`$ can be varied independently of $`T_m`$, $`i.e.`$ independently of $`Z`$, since one can change the magnitude of the constriction area $`S_c`$ independently. In other words, $`R_d/R_b`$ is no more proportional to $`T_{av}(L/l)`$, where $`T_{av}`$ is the averaged transmissivity of the barrier and $`l`$ is the mean free path in the diffusive region, respectively. Based on this fact, we can choose $`R_d/R_b`$ and $`Z`$ as independent parameters.

In the following section, we will discuss the normalized tunneling conductance $`\sigma _T(eV)=\sigma _S(eV)/\sigma _N(eV)`$ where $`\sigma _N(eV)`$ is the tunneling conductance in the normal state given by $`\sigma _N(eV)=\sigma _N=1/(R_d+R_b)`$, respectively.

## III Results

### III.1 Tunneling conductance vs voltage: zero-bias anomalies

In this section, we focus on the line shape of the tunneling conductance. Let us first choose the relatively strong barrier $`Z=10`$ (Figs. 1 and 2) for various $`R_d/R_b`$. For $`E_{Th}=\mathrm{\Delta }_0`$, the magnitude of $`\sigma _T(eV)`$ for $`eV<\mathrm{\Delta }_0`$ increases with the increase of $`R_d/R_b`$. First, the line shape of the tunneling conductance remains to be $`U`$ shaped and only the height of the bottom value is enhanced (curve $`b`$ ). Then, with a further increase of $`R_d/R_b`$, a rounded bottom structure (curve $`c`$ and $`d`$) appears and finally it changes into a nearly flat line shape (curve $`e`$).

For $`E_{Th}=0.01\mathrm{\Delta }_0`$ (Fig. 2), the magnitude of $`\sigma _T(eV)`$ has a ZBCP once the magnitude of $`R_d/R_b`$ deviates slightly from 0. The order of magnitude of the ZBCP width is given by $`E_{Th}`$. When the magnitude of $`R_d/R_b`$ exceeds unity, the ZBCP splits into two (curve $`d`$) and finally $`\sigma _T(eV)`$ acquires a zero bias conductance dip (ZBCD) (curve $`e`$). In the limit of an extremely strong barrier with $`Z\mathrm{}`$, the results of the VZK theory are reproduced.

On the other hand, in the fully transparent case with $`Z=0`$ the line shape of the tunneling conductance becomes quite different. For $`E_{Th}=\mathrm{\Delta }_0`$ (Fig. 3), the magnitude of $`\sigma _T(eV)`$ decreases with the increase of the magnitude of $`R_d/R_b`$. The bottom parts of all curves are rounded and $`\sigma _T(eV)`$ always exceeds unity. On the other hand, for $`E_{Th}=0.01\mathrm{\Delta }_0`$, $`\sigma _T(eV)`$ has a ZBCD even for a small magnitude of $`R_d/R_b`$ except for the special case of $`R_d/R_b=0`$ where the BTK theory is valid. This feature is quite different from that shown in Fig. 2.

In the case of an intermediate barrier strength, $`Z=1`$, the shape of $`\sigma _T(eV)`$ becomes rather complex. For $`E_{Th}=\mathrm{\Delta }_0`$, $`\sigma _T(eV)`$ has a shallow gap structure similar to the case of the BTK theory (curve $`a`$ ). With the increase of $`R_d/R_b`$, the coherent peak structure at $`eV=\pm \mathrm{\Delta }_0`$ is smeared out and the voltage dependence becomes very weak as shown by curve $`e`$. For $`E_{Th}=0.01\mathrm{\Delta }_0`$, the ZBCP appears for a small magnitude of $`R_d/R_b`$ (see curve $`b`$). With increasing magnitude of $`R_d/R_b`$, the ZBCP changes into a ZBCD (curves $`c`$ $`d`$ $`e`$). As compared to the $`Z=10`$ case (see Fig. 2), the ZBCP is much more easier to change into a ZBCD by increasing the magnitude of $`R_d/R_b`$.

It is interesting to study how various parameters influence the proximity effect. The measure of the proximity effect at the S/N interface $`\theta _L`$ is plotted for $`Z=0`$ and $`Z=10`$ with corresponding parameters in Figs. 1 to 4. For $`R_d/R_b=0`$, $`\theta _L=0`$ is satisfied for any $`E_{Th}`$ and $`Z`$. Besides this fact, at $`ϵ=0`$, $`\theta _L`$ always becomes a real number. First, we study the case of $`E_{Th}/\mathrm{\Delta }_0=1`$ (Fig. 7) where the same values of $`R_d/R_b`$ are chosen as in Figs. 1 and 3. The real part of $`\theta _L`$ is enhanced with an increase in $`R_d/R_b`$ and is almost constant as function of $`ϵ`$. At the same time, the imaginary part of $`\theta _L`$ is an increasing function of $`ϵ`$ for all cases. There is no clear qualitative difference between the energy dependencies of $`\mathrm{Real}(\theta _L)`$ and $`\mathrm{Imag}(\theta _L)`$ at $`Z=0`$ and $`Z=10`$.

Next, we discuss the line shapes of $`\theta _L`$ for $`E_{Th}/\mathrm{\Delta }_0=0.01`$. Real($`\theta _L`$) has a peak at zero voltage and decreases with the increase of $`ϵ`$. Imag($`\theta _L`$) increases sharply from 0 and has a peak at about $`ϵE_{Th}`$, except for a sufficiently large value of $`R_d`$. Also in this case, there is no qualitative difference between the line shapes of $`\mathrm{Real}(\mathrm{Imag})(\theta _L)`$ for the $`Z=0`$ case and that for $`Z=10`$.

Although the magnitude of $`\theta _L`$, $`i.e.`$ the measure of proximity effect, is enhanced with increasing $`R_d/R_b`$, its influence on $`\sigma _T(eV)`$ is different for low and high transparent junctions. In the low transparent junctions, the increase in the magnitude of $`\theta _L`$ by $`R_d/R_b`$ can enhance the conductance $`\sigma _T(eV)`$ for $`eV0`$ and produce a ZBCP, whereas in high transparent junctions the enhancement of $`\theta _L`$ suppresses the magnitude of $`\sigma _T(eV)`$.

Finally, we focus on the condition where the ZBCP appears. We change $`Z`$ and $`R_d/R_b`$ for fixed $`E_{Th}`$. An upper critical value of $`R_d=R_{bu}`$ exists where the ZBCP vanishes for $`R_d>R_{bu}`$. For $`E_{Th}=0.01\mathrm{\Delta }_0`$, $`R_{bu}`$ increases with $`Z`$ and converges at nearly 1.4. For $`E_{Th}=0.8\mathrm{\Delta }_0`$, the lower critical value of $`R_d=R_{bl}`$ also appears where the ZBCP vanishes for $`R_d<R_{bl}`$. The ZBCP is expected for $`R_{bl}<R_d<R_{bu}`$. The magnitude of $`R_{bu}`$ is suppressed drastically as compared to that for $`E_{Th}=0.01\mathrm{\Delta }_0`$. For $`E_{Th}>\mathrm{\Delta }_0`$, a ZBCP region vanishes for $`0<Z<20`$.

In order to understand the crossover from the ZBCP to the ZBCD in much more detail, it is interesting to calculate the second derivative of the total resistance $`R`$ as a function of $`ϵ=eV`$ at $`ϵ=0`$. For simplicity, here we focus on the case of a sufficiently large $`Z`$. For $`ϵ<\mathrm{\Delta }_0`$, for simplicity total resistance is written as

$$R=R_1+R_2,$$
(21)
$$R_1=\frac{R_b}{<I_{b0}>},R_2=\frac{R_d}{L}_0^L\frac{dx}{\mathrm{cosh}^2\theta _i(x)},$$

$$<I_{b0}>=f\mathrm{sin}\theta _{L,r}\mathrm{cosh}\theta _{L,i},$$

with $`\theta _L=\theta _{L,r}+i\theta _{L,i}`$, where $`\theta _{L,r}`$ and $`\theta _{L,i}`$ are the real and imaginary parts of $`\theta _L`$, respectively. The second derivative of $`I_{b0}`$ at $`ϵ=0`$ is given by

$$\frac{^2I_{b0}}{ϵ^2}=\frac{\mathrm{sin}\theta _{L,r}}{\mathrm{\Delta }_0^2}+\mathrm{cos}\theta _{L,r}\frac{^2\theta _{L,r}}{ϵ^2}+\mathrm{sin}\theta _{L,r}(\frac{\theta _{L,i}}{ϵ})^2,$$

since $`\frac{\theta _{L,r}}{ϵ}=0`$ and $`\theta _{L,i}=0`$ is satisfied at $`ϵ=0`$.

The resulting second derivative of the total resistance at zero energy is given as

$$\frac{^2R}{ϵ^2}=\frac{^2R_1}{ϵ^2}+\frac{^2R_2}{ϵ^2},$$

$$\frac{^2R_1}{ϵ^2}=$$

$$\frac{R_b}{\mathrm{sin}^2\theta _{L,r}}[\frac{\mathrm{sin}\theta _{L,r}}{\mathrm{\Delta }_0^2}+\mathrm{cos}\theta _{L,r}\frac{^2\theta _{L,r}}{ϵ^2}+\mathrm{sin}\theta _{L,r}(\frac{\theta _{L,i}}{ϵ})^2],$$
(22)
$$\frac{^2R_2}{ϵ^2}=\frac{2R_d}{L}_0^L[\frac{\theta _i(x)}{ϵ}]^2𝑑x,$$
(23)
where $`\theta _i(x)`$ is the imaginary part of $`\theta (x)`$. The sign of $`\frac{^2\theta _{L,r}}{ϵ^2}`$ becomes negative and it can induce $`\frac{^2R_1}{ϵ^2}>0`$ in some cases. The order of $`\frac{^2\theta _{L,r}}{ϵ^2}`$ and $`(\frac{\theta _{L,i}}{ϵ})^2`$ is proportional to the inverse of $`E_{Th}^2`$. The sign of $`\frac{^2R_1}{ϵ^2}`$ is crucially influenced by the relative magnitude of the second term at the right-hand side of Eq. (22). On the other hand, the sign of $`\frac{^2R_2}{ϵ^2}`$ is always negative.

For $`E_{Th}<<\mathrm{\Delta }_0`$ and small magnitude of $`R_d/R_b`$, the resulting $`\theta _{L,r}`$ is sufficiently small and the magnitude of $`\frac{^2R_1}{ϵ^2}`$ becomes positive. When this positive contribution overcomes the negative contribution from $`\frac{^2R_2}{ϵ^2}`$, we can expect a resistance minimum at zero energy, $`i.e.`$, a ZBCP. However, with increasing $`R_d/R_b`$, the magnitude of $`\frac{^2R_1}{ϵ^2}`$ decreases due to the enhancement of the third term in the right-hand side of Eq. (22), while $`\frac{^2R_2}{ϵ^2}`$ increases. Then, a critical value $`R_d=R_{bu}`$ appears, above which the ZBCP changes into a ZBCD with an increase of $`R_d/R_b`$. This is the mechanism of the crossover from a ZBCP to a ZBCD.

When the magnitude of $`E_{Th}`$ is enhanced, the magnitudes of the second and third terms in Eq. (22) are reduced and the first term can not be neglected. The resulting magnitude of $`\frac{^2R_1}{ϵ^2}`$ is suppressed and the value of $`R_{bu}`$ is reduced. This is the origin of the difference in $`R_{bu}`$ in Figs. 9(a) and (b).

It is also interesting to present similar plots as Fig. 9 using an averaged transparency of the junction $`T_{av}`$

$$T_{av}=\frac{_{\pi /2}^{\pi /2}\mathrm{cos}\varphi T(\varphi )𝑑\varphi }{2}.$$
(24)
The results are shown in Fig. 10. For $`E_{Th}=0.01\mathrm{\Delta }_0`$, the magnitude of $`R_{bu}`$ decreases monotonically with increasing $`T_{av}`$ and it vanishes about $`T_{av}0.8`$. For $`E_{Th}=0.8\mathrm{\Delta }_0`$, the magnitude of $`R_{bl}`$ increases for an increasing $`T_{av}`$ while that for $`R_{bu}`$ decreases. The ZBCP region is restricted to small values of $`R_d/R_b`$.

In any way, the preferred condition for the formation of a ZBCP is the combination of the low transparency of the junction and the smallness of the $`E_{Th}/\mathrm{\Delta }_0`$ ratio. This situation is understood as follows. It is well known from the BTK theory that the magnitude of the zero bias conductance is almost proportional to $`T_{av}^2`$ for $`R_d=0`$ for $`T_{av}<<1`$. With the increase in the magnitude of $`R_d`$, the measure of the proximity effect $`\theta _L`$ is enhanced for $`ϵ<E_{Th}`$ and the zero bias conductance is proportional to $`T_{av}`$. However, the magnitude of $`\theta _L`$ at finite energy in the range $`E_{Th}<ϵ<\mathrm{\Delta }_0`$ is drastically suppressed as shown in Fig. 8, due to the existence of the proximity induced minigap minigap  in the normal diffusive part DN of the order of $`E_{Th}`$. As a result, in this regime the conductance channel that provides a contribution proportional to $`T_{av}`$ is not available. Thus only the low voltage conductance is enhanced. On the other hand, for large $`E_{Th}`$ with $`E_{Th}\mathrm{\Delta }_0`$, the measure of proximity effect $`\theta _L`$ is insensitive to energy for $`ϵ<\mathrm{\Delta }_0`$ (see Fig. 7), then the resulting $`\sigma _T(eV)`$ is always enhanced for $`eV<\mathrm{\Delta }_0`$ and the degree of the prominent enhancement of $`\sigma _T(0)`$ is weakened.

### III.2 Tunneling conductance vs temperature: reentrance effect

Finally, we look at temperature dependence of conductance for various $`Z`$ and $`R_d/R_b`$ and focus on the relevance to the corresponding results in VZK theory based on the KL boundary condition. We calculate tunneling conductance at non zero temperature following

$$\sigma _S(eV,T)=dI_{el}/dV.$$
(25)
Then we define deviation of tunneling conductance from that at zero temperature given by

$$\delta \sigma _S=\sigma _S(eV,T)\sigma _S(eV,0).$$
(26)
In the following, we will plot $`\delta \sigma _S/\sigma _N`$ as a function of temperature $`T`$.

For $`R_d/R_b=0.1`$ (curve $`a`$) and $`R_d/R_b=1`$ (curve $`b`$), due to the existence of ZBCP as shown in Fig. 2, $`\delta \sigma _S/\sigma _N`$ takes negative value and decreases with $`T`$. While for $`R_d/R_b=2`$, $`\delta \sigma _S/\sigma _N`$ first increases and decreases again (curve $`c`$). With the further increase of the magnitude of $`R_d/R_b`$, $`\delta \sigma _S/\sigma _N`$ becomes positive and it has maximum at a certain temperature. This effect is known as ”reentrance effect”. It was predicted theoretically within the VZK theory in Stoof ; Reentrance ; Golubov  and observed experimentally Reentrance\_exp . According to the theory Stoof ; Golubov ; Reentrance , for $`R_b0`$, the maximum value of $`\delta \sigma _S/\sigma _N`$ is about 0.09 at temperature of the order of Thouless energy. Note that the reentrance of the metallic conductance occurs as a function of bias voltage as well, but here we concentrate on the temperature dependence since it was studied in most detail within the VZK theory.

In order to study the consequences of our theory for the reentrance effect, we have calculated the temperature dependence of $`\delta \sigma _S/\sigma _N`$ for various $`Z`$ for fixed $`R_d/R_b`$. For the comparison with the standard VZK theory, we also plot $`\delta \sigma _S/\sigma _N`$ using the KL boundary condition. For $`R_d/R_b=10`$, we always see the standard reentrant behavior. For $`Z=10`$, we can not see clear deviation from the VZK theory (compare curves $`a`$ and $`b`$). With the decrease of the magnitude of $`Z`$, the magnitude of $`\delta \sigma _S/\sigma _N`$ is enhanced much stronger (curves $`c`$ and $`d`$). Although there are quantitative difference between four curves, the qualitative line shapes are similar to those predicted by the VZK theory.

However, situation is different for decreasing $`R_d/R_b=1`$. For $`R_d/R_b=1`$ the resulting $`\delta \sigma _S/\sigma _N`$ takes negative value and similar feature is also obtained for $`Z=10`$ (see curves $`a`$ and $`b`$). At the same time, for small magnitude of $`Z`$, $`\delta \sigma _S/\sigma _N`$ takes positive value. As seen from these results, the deviation from VZK becomes significant for small magnitude of $`Z`$ with $`R_d/R_b<1`$.

## IV Conclusions

In the present paper, a detailed theoretical investigation of the tunneling conductance of diffusive normal metal / conventional superconductor junctions is presented. Even for conventional s-wave junctions, the interplay between diffusive and interface scattering produces a wide variety of line shapes of the tunneling conductance: ZBCP, ZBCD, U-shaped, and rounded bottom structures. There are several points which have been clarified in this paper.

1. When the transparency of a junction is sufficiently low and $`E_{Th}`$ is small, the ZBCP appears for $`R_d<R_b`$. With an increase in $`R_d`$, the ZBCP changes into a zero bias conductance dip (ZBCD). For large $`E_{Th}`$ with the same order of $`\mathrm{\Delta }_0`$, the ZBCP is only expected for small values of $`R_d`$. With increasing $`E_{Th}`$, the ZBCP vanishes when $`E_{Th}>\mathrm{\Delta }_0`$ is satisfied. For the low transparency limit, the results obtained by us are reduced to those by the VZK theory where the KL boundary condition is used.

2. When the transparency of the junction is almost unity, $`\sigma _T(eV)`$ always have a ZBCD except for the case of vanishing $`R_d`$.

3. The proximity effect can enhance (reduce) the tunneling conductance of junctions with low (high) transparency.

The above mentioned criteria for the existence of a ZBCP agree with available experimental data. However, we are not aware of an experimental observation of ZBCD in highly transparent junctions.

In the present paper, the superconductor is restricted to be a conventional s-wave superconductor. However, it is well known that a ZBCP also appears in unconventional superconductor junctions TK95 ; Kashi00 , the origin of which is the formation of midgap Andreev bound states (MABS) Buch . Indeed, a ZBCP has been reported in various superconductors that have an anisotropic pairing symmetry. Kashi00 ; Experiments  It should be remarked that the line shape of the ZBCP obtained in the present paper is quite different from that by MABS. In the present case, the height of $`\sigma _T(eV)`$ never exceeds unity and its width is determined by $`E_{Th}`$, while in the MABS case, the peak height is proportional to the inverse of the magnitude of $`T_m`$ and the width is proportional to $`T_m`$.

The proper theory of transport of unconventional junctions in the presence of MABS has been formulated TK95 ; Kashi00  only under the conditions of ballistic transport. Recently, this theory has been revisited to account for diffusive transport in the normal metal in Ref.Nazarov3 , where an extension of the circuit theory was provided for unconventional superconductor junctions. A general relation was derived for matrix current, $`B`$ in Eq. (4), which is available for unconventional superconductor junctions with MABS. An elaborated example demonstrated the interplay of MABS and proximity effect in a $`d`$-wave junction. It is actually quite interesting to apply this novel circuit theory to the actual calculation of $`\sigma _T(eV)`$ as in the present paper. Such a direction of study is important in order to analyze recent tunneling experiments where mesoscopic interference effects were observed in high-T<sub>C</sub> junction systems Hiromi . We will present the obtained results in the near future Golubov2 .

In the present study, we have focused on N/S junctions. The extension of Nazarov’s circuit theory to long diffusive S/N/S junctions has been performed by Bezuglyi et al Bezuglyi . In S/N/S junctions, the mechanism of multiple Andreev reflections produces the subharmonic gap structures on I-V curves M1 ; M2 ; M3 ; M4 ; M5 ; M6 ; M7 ; M8  and the situation becomes much more complex as compared to N/S junctions. Moreover, in S/N/S junctions with unconventional superconductors, MABS lead to the anomalous current-phase relation and temperature dependence of the Josephson current TKJ . An interesting problem is an extension of the circuit theory to S/N/S junctions with unconventional superconductors.

The authors appreciate useful and fruitful discussions with Yu. Nazarov, J. Inoue, and H. Itoh. This work was supported by the Core Research for Evolutional Science and Technology (CREST) of the Japan Science and Technology Corporation (JST). The computational aspect of this work has been performed at the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo and the Computer Center.