Tunneling Hamiltonian analysis of DC Josephson currents in a weakly-interacting Bose-Einstein condensate

Atomtronics experiments with ultracold atomic gases allow us to explore quantum transport phenomena of a weakly-interacting Bose-Einstein condensate (BEC). Here, we focus on two-terminal transport of such a BEC in the vicinity of zero temperature. By using the tunnel Hamiltonian and Bogoliubov theory, we obtain a DC Josephson current expression in the BEC and apply it to experimentally relevant situations such as quantum point contact and planar junction. Due to the absence of Andreev bound states but the presence of couplings associated with condensation elements, a current-phase relation in the BEC is found to be different from one in an s-wave superconductor. In addition, it turns out that the DC Josephson current in the BEC depends on the sign of tunneling elements, which allows to realize the so-called $\pi$ junction by using techniques of artificial gauge fields.


I. INTRODUCTION
Quantum transport between macroscopic reservoirs reveals colorful quantum phenomena. When macroscopic reservoirs are separated by a two-dimensional barrier, a variety of quantum tunneling phenomena that depends on quantum states of matter in reservoirs emerge [1]. A prototypical example is a superconducting tunneling junction where an insulating barrier is sandwiched by superconductors. In such a junction also known as the Josephson junction, electrons can flow even in the absence of a voltage bias, which is also an important building block of superconducting qubits [2]. When macroscopic reservoirs are connected by mesoscopic samples, in addition, quantum nature in transport is known to be enhanced [3]. If reservoirs are filled in superconductors, for instance, the mesoscopic supercurrent is essentially determined by properties of Andreev bound states [4]. In order to reveal such quantum transport phenomena, typical condensed matter systems composed of semiconductors and superconductors have conventionally been examined.
When it comes to a purely fundamental physics point of view, one may take advantage of controllable artificial systems to deepen understanding of quantum transport. Ultracold atomic gases are one of such systems [5], and now allow to realize typical quantum transport systems such as the tunneling junction [6][7][8] and mesoscopic systems [9,10]. An interesting perspective in transport research with ultracold atomic gases is to explore regimes hard to attain with solid-state materials.
A weakly-interacting Bose-Einstein condensate (BEC) realized in ultracold atomic gases is a distinctive example in that one can examine a variety of transport properties of a BEC in a controlled way. Since the realization of a weakly-interacting BEC in ultracold atomic gases, numerous theoretical studies on tunneling phenomena between BECs have been done in connection with the celebrated Josephson effect [11][12][13], which now becomes one chapter of the typical textbooks [14,15]. In addition to the theoretical works, the corresponding exper-imental works on the Josephson physics have also been done [6,7,9,[16][17][18][19][20][21][22][23][24][25]. Josephson effects emerging in junction systems can be formulated in terms of the so-called tunneling Hamiltonian [26]. Then, a standard assumption is that tunneling amplitudes are so small that an analysis up to the linear response theory is justified. However, it is not always trivial on whether or not experiments look at the linear response regime in reality, and it may also be important to estimate effects beyond the linear response theory [27][28][29]. In the typical mesoscopic setups such as a quantum point contact, moreover, it is easy to reach the regime where higher-order tunneling effects are nonnegligible [30].
The purpose of this paper is to obtain a DC Josephson current of a weakly-interacting BEC that includes higher-order tunneling processes in the tunneling Hamiltonian. To this end, we apply the Bogoliubov theory that explains a BEC near zero temperature to the tunneling Hamiltonian, and adopt the Keldysh formalism to take into account higher-order tunneling processes in an efficient manner. As discussed in Refs. [31,32], dominant contributions in low-energy transport between BECs are tunneling processes associated with condensation elements. Therefore, we focus on effects of such contributions, which is reasonable near zero temperature, since the occupancy effect of Bogoliubov modes is negligible. In consideration of cold-atom experiments, our formulation is specifically applied to quantum point contact and planar junction (see also Fig. 1). We find that as increasing the tunnel coupling, the DC Josephson current in a BEC deviates from the conventional sinusoidal arXiv:2106.12750v2 [cond-mat.quant-gas] 21 Oct 2021 curve in a different way from one in an s-wave superconductor. We note that our formulation is also applicable to other BEC junction systems such as spinor BECs in ultracold atomic gases [33,34] and magnon BECs in condensed matter [35] if the corresponding quantum point contact and planar junction systems are realized. This paper is organized as follows. Section II discusses the Bogoliubov theory and tunneling Hamiltonian, which give a basis of formulation of DC Josephson currents in a weakly-interacting BEC. In Sec. III, the Keldysh formalism that allows us to examine higher-order tunneling processes is applied to the tunneling Hamiltonian. Section IV obtains dominant contributions of DC Josephson currents in a BEC and analyzes quantum point contact and planar junction cases in a concrete manner. Section V discusses related subjects such as internal Josephson effect, thermal current, trapped systems, and possible experiments. For a self-contained description, Green's function formulae used in quantum point contact and planar junction systems and DC Josephson current formula in an s-wave superconductor are respectively given in Appendices A and B.

II. FORMULATION OF THE PROBLEM
We consider two-terminal transport systems where two macroscopic reservoirs are weakly coupled via transport channels. In our formulation, we assume that two reservoirs are in equilibrium in that the thermodynamic variables such as temperature, chemical potential, and phase degree of freedom are fixed in time. If the reservoirs are infinite, this assumption is reasonable. In ultracold atomic gases where a finite system is concerned, however, in the presence of a coupling between reservoirs, the particle number in each reservoir can change in time. Thus, the applicability of the above assumption may not be trivial.
At the same time, in the typical mesoscopic setups with ultracold atomic gases such as quantum point contact and tunneling junction, we encounter a situation that the thermalization time in each reservoir is much faster than the transport time that is an equilibrated time of the total system including both reservoirs and channel. In such a situation, one may use the so-called quasi-steady approximation where two reservoirs are in equilibrium in each time of the evolution. Since an analysis based on the quasi-steady approximation has a great success in cold atoms [10], in what follows, two-terminal transport properties in the presence of a fixed phase bias between reservoirs are discussed.

A. Bogoliubov theory
First, we focus on bulk properties of a weaklyinteracting BEC. The system of interest is that a number of non-relativistic identical bosons are trapped in each reservoir, where each atom experiences an s-wave collision. We can write down the corresponding reservoir Hamiltonian as follows: where α = L or R is the reservoir index, k = k 2 /2M , and µ and g are respectively the chemical potential and coupling constant, each of which takes a same value between reservoirs. For the stability reason, the coupling g is assumed to be positive, that is, repulsive interaction. In this paper, we set = k B = 1 and the volume in each reservoir is taken to be unity. We now consider a situation that a temperature of the system is near absolute zero and therefore bosons in each reservoir form a BEC. Since each boson is weakly interacting, we can apply the so-called Bogoliubov prescription to the above Hamiltonian [36]. Namely, we assume that the k = 0 component gives the dominant contribution and is replaced by c-number such that b 0,α = √ v α e −iφα . Here, v α and φ α are physically the condensate fraction and phase degree of freedom in a BEC, respectively. In the leading order analysis, the Hamiltonian is replaced by a c-number, which is a function of v α . Notice that up to the leading order, φ α dependence in the Hamiltonian is absent, since the original Hamiltonian consists of the combination b † b (U(1) symmetry). Then, v α itself is determined from the stationary condition. By assuming the uniform solution, we obtain We next consider a fluctuation effect from the cnumber analysis above. Namely, we take into account k = 0 components that are indeed generated in the presence of the interaction. By using the stationary condition above, one can show that the term linear in b k =0,α vanishes. Thus, the first nontrivial fluctuation effect comes from the term quadratic in b k =0,α . Up to this order, which is indeed the matter in the Bogoliubov theory, we obtain the following effective Hamiltonian: where we neglect the constant term. Since the above Hamiltonian contains b k b −k and b † k b † −k where the explicit phase dependence appears, one can introduce the following Bogoliubov transformation: with Then, the effective Hamiltonian is diagonalized as H α = k E k β † k,α β k,α , where we again neglect the constant term. We note that as a result of the common chemical potential, the Bogoliubov excitation energy E k does not depend on reservoirs. Since the Bogoliubov transformation is canonical, it follows that β k,α obeys the following commutation relation: Below, Eqs. (3)-(6) are used to obtain DC Josephson currents.

B. Tunneling Hamiltonian
We now introduce a coupling between reservoirs that induces transport. We are particularly interested in a situation that such a coupling is described by the following tunneling term: Here t pk is the tunneling amplitude that obeys t * pk = t −p−k in systems with time reversal symmetry. Thus, the Hamiltonian of the total system is described by sum of reservoirs' Hamiltonian and tunneling term H = H L + H R + H T , which is the so-called tunneling Hamiltonian. Josephson current analyses based on tunneling terms are successfully applied to discuss low-energy transport of tunnel junction [31,37], quantum point contact [32,38,39], and quantum dot systems [40].
In the presence of H T , the particle number in each reservoir is not conserved and a nonzero particle current is induced in the presence of a bias between reservoirs. By taking the convention that a current is positive when particles flow from left to right reservoirs, the particle current operator in the tunneling Hamiltonian is expressed as where N L = k b † k,L b k,L is the particle number operator in reservoir L, and we use the Heisenberg equation. By using the boson commutation relation, the particle current operator is reduced to We note that the above expression of I N is an operator relation and therefore is available regardless of quantum states in detail [41]. We now take into account the fact that each reservoir is filled by a BEC, where there is the global U(1) phase degree of freedom φ L(R) . Our interest is the current induced by a relative phase between BECs, ∆φ = φ L − φ R . To this end, it is convenient to perform the U(1) gauge dependence present in Eq. (3) vanishes, and instead the tunneling term is transformed as follows: where the relative phase dependence is present. It is then straightforward to confirm that due to this gauge transformation, I N also gains the corresponding phase factor e ±i∆φ .
In what follows, we calculate the currents based on this transformed Hamiltonian [42].

III. CURRENT EXPRESSION WITH REAL TIME GREEN'S FUNCTIONS
In order to examine quantum transport properties of a BEC, we use the following Nambu representation of the field operator:b By using this representation, the tunneling term is expressed as where we use the property of time reversal symmetry. In this work, we adopt the Keldysh formalism [43,44] to calculate the DC Josephson current. Under the Nambu representation, contour-ordered Green's function has the following 2 × 2 matrix structure: where T c is the contour ordering, and a, b = ± denotes upper (+) and lower (−) branches on the contour [43,44]. We note that the above average · · · is taken under the total Hamiltonian H including the tunneling term. In particular, the component of a = + and b = − is lesser Green's function given bŷ By means of lesser Green's function introduced above, the average particle current can be expressed as whereσ z = 1 0 0 −1 and the trace is taken over the Nambu space. For later use, we also introduce retarded and advanced Green's functionŝ In the absence of H T , each reservoir is decoupled and Green's functions can be determined by using the Bogoliubov theory. If we denote the corresponding uncoupled Green's functionsĝ, whose derivation is given in Appendix A, the Dyson equation allows us to relateĜ withĝ by treating H T as perturbation.
Since H T is the single-particle potential term, the Dyson equation for G R(A) has the following simple structure [43]:Ĝ Here −t pk e i∆φ ort † kp , and • represents summation over the internal momentum variable and integration over the internal time variable from minus infinity to plus infinity. On the other hand, by using the so-called Langreth rules [43], the Dyson equation forĜ < is obtained aŝ In order to obtain a convenient expression of the current, we now make the reservoir indices explicit. Then, the Dyson equations forĜ R,A are expressed aŝ (24) wheret is the short hand notation oft pk . In addition, the Dyson equation forĜ < is rewritten aŝ wherê In the above expressions, odd and even in the superscript respectively represent current components of odd-order and even-order terms in tunneling amplitudes, which can be confirmed by looking at the Dyson equation forĜ R(A) and Eq. (16) that contains an additional tunneling amplitude dependence. We note thatĜ <,odd andĜ <,even stem fromĝ < LR(RL) andĝ < LL(RR) contributions, respectively. What is peculiar in the BEC system is the presence ofĝ <

LR(RL)
that arises from the c-number component of the field operator, and thereforeĝ consists of the commutator of the field operators (see also Appendix A). The presence ofĜ <,odd is allowed by virtue of the tunnel coupling between condensation elements and one between condensation and non-condensation elements, and thereforeĜ <,odd is absent in fermionic superconductors (see also Appendix B). This difference between fermionic and bosonic systems can also be understood as the fact that while in fermionic superconductors a phase degree of freedom responsible for Josephson currents emerges from condensation of Cooper pairs whose transport is associated with even-order processes in the tunneling amplitude, in the bosonic condensate discussed in this paper, odd-order terms in the tunneling amplitude are allowed owing to the property b = 0 [45].
For convenience, we also introduce renormalized tunneling matrices defined aŝ By using above, it is straightforward to show the following relations: Thus, we find Tr σ zt •Ĝ <,even The equations above give a foundation of DC Josephson currents in a weakly-interacting BEC.

IV. RESULTS
By using Eqs. (33) and (34), we show typical behaviors of DC Josephson currents in the vicinity of zero temperature, where occupation effects of the Bogoliubov mode are negligible. As theoretically discussed in Refs. [31,32,46,47] and experimentally verified in Refs. [23,25], the dominant contribution of a mesoscopic BEC in such a regime originates from couplings associated with the condensation elements. We point out that in the Keldysh formalism presented in this paper such a dominant contribution is extracted from the replacement g < k,αβ →ĝ < 0,αβ . In addition, since there is no chemical potential bias between reservoirs, an AC Josephson current is absent and the current does not depend on time.
Thus, convenient expressions can be obtained by using the Fourier transformation. In frequency space, the DC Josephson current is expressed as where → represents the prescription that the dominant contribution is picked up. The odd-order and even-order terms in tunneling amplitudes are respectively expressed as where1 = 1 0 0 1 andt pk = −t pk e −i∆φ 0 0 −t pk e i∆φ in the frequency-space representation. Equations (35)- (37) are the main results of this paper, which are applicable to generic tunneling amplitudes within the perturbative meaningful regime. In what follows, DC Josephson currents in point-contact and planar junction systems, which are relevant to cold-atom experiments, are separately discussed.

A. Quantum point contact case
Here we discuss the case of a quantum point contact where a particle transfer between reservoirs occurs through a one-dimensional wire. The tunneling Hamiltonian description of a quantum point contact is known to be reasonable when short one-dimensional wire with channel transmittance whose energy dependence is negligible is concerned [29,32,38,39,[47][48][49][50][51]. In particular, the tunneling term in a quantum point contact is modeled in such a way that the tunneling occurs at a single point in each reservoir, e.g. at origin [48]. In the momentum space representation, the corresponding tunneling term is expressed as with the tunneling amplitude t. Thus, compared with the general expression (11), a momentum-independent tunneling amplitude is now concerned in transport of a quantum point contact.
Since the momentum dependence in the tunneling amplitude is absent, the renormalized tunneling matrices become momentum independent. By using the frequencyspace representation, we find that Eqs (28) and (29) are expressed aŝ αα (ω) is uncoupled local retarded (advanced) Green's function whose expression is given in Appendix A. Since the above are 2 × 2 matrix equations, they can be solved as , f R(A) are respectively 11, 22, 12 matrix elements ofĝ R(A) . In addition, by using lesser Green's function expression shown in Eq. (A2), we find that the dominant contribution of the DC Josephson current in a quantum point contact is obtained as wheret −1 = − 1 t e i∆φ 0 0 e −i∆φ , and Φ(0) and Ψ(0) are real parts of g R 1 (0) and f R (0), respectively (see also Appendix A).
We now comment on basic features of the DC Josephson current formula above. First, as is expected from generic properties of Josephson currents, Eq. (44) is an odd function in ∆φ, that is, I N (−∆φ) = −I N (∆φ). Second, in contrast to fermionic superconductor and superfluid systems, where Josephson currents consist of evenorder terms in t, I N above contains odd-order terms in t. The presence of odd-order terms in t shows that the current depends on the sign of t. However, we point out that the effect of the sign of t can be absorbed by the following shift of the relative phase: ∆φ → ∆φ + π. Thus, it turns out that the result of a negative t can be inferred from one of a positive t.
To see some intuition of Eq. (44), we next perform the series expansion in t. The leading term in Eq. (44) is linear in t and is given by Up to this order, the current obeys the conventional sinusoidal current-phase relation. We note that the result above can also be obtained with substitution of the macroscopic wave function for b † k=0,L(R) and b k=0,L(R) in Eq. (10). We next look at the sub-leading term in Eq. (44), which is quadratic in t. From the series expansion of Eq. (44), it is obtained as Thus, the argument of the sine function in this term is not ∆φ but 2∆φ. We point out that the essentially same result can also be obtained by using the linear response theory, where the tunnel current up to t 2 is concerned. It may also be clear from the analysis above that higherorder terms in t generate current components proportional to sin n∆φ with n ≥ 3. We now directly look at Eq. (44), which contains all order effects in t. Figure 2 (a) shows typical behaviors of the current-phase relation at different tunneling amplitudes. Here the current is plotted in units of I 0 = 2µt/g with a positive t, and the horizontal axis is restricted as 0 ≤ ∆φ < 2π due to the 2π periodicity. In addition, we set Φ(0) = −Ψ(0), since Φ(0) contains the cutoff Λ but we check that the value of Λ is irrelevant to the qualitative discussion of the current-phase relation. Owing to this setting, the current-phase relation in a BEC can be characterized by a single dimensionless parameter |tΨ(0)|. We note that the condition |tΨ(0)| < 1/2 is necessary to obtain a perturbatively meaning result. In terms of the physical model parameters, the current can also be expressed as with the speed of sound c = µ M . In a small tunneling amplitude (blue curve), the current-phase relation obeys the conventional Josephson relation Eq. (45). In larger tunneling amplitudes such as purple, orange, and red curves, on the other hand, we can see clear deviations from the conventional relation. There, the current takes large values in the vicinity of ∆φ = 0 and 2π, and is suppressed except for them. Such a tendency is somewhat akin to one obtained by the Gross-Pitaevskii analysis with a barrier potential [52][53][54][55]. At the same time, we note that due to symmetry of the sign in t mentioned before, the current-phase relation for a negative t is shifted by π and its amplitude takes large values in the vicinity of ∆φ = π. In ultracold atomic gases, the current-phase relation for a negative t that shows the π-junction behavior may also be achieved by using techniques of synthetic gauge fields [56]. For comparison, we show the DC Josephson current in an s-wave superconductor in Fig. 2

(b) (derivation with the tunneling Hamiltonian is given in Appendix B).
In the case of a superconductor, the DC Josephson current is carried by tunneling of Cooper pairs via multiple Andreev reflections. Due to the absence of a chemical potential bias, the loop of multiple Andreev reflections is closed and leads to the presence of Andreev bound states. As in the case of a BEC, the Josephson current in a superconductor deviates from the conventional sinusoidal curve as increasing the tunneling amplitude. This deviation from the conventional curve now originates from multiple Andreev reflections. In contrast to the BEC case, as increasing the channel transmission parameter, the Josephson current becomes maximum near ∆φ = π, which does not depend on the sign of the tunneling amplitude due to the absence ofĜ <.odd in a superconductor.

B. Planar junction case
We next consider the planar junction system where there is a translational invariance along the interface. As pointed out in Ref. [31], one can introduce three independent tunneling couplings in such a junction: coupling between condensation elements, one between condensation and non-condensation elements, and one between non-condensation elements. The corresponding tunneling amplitudes can respectively be expressed as [31] t 00 = t cc δ 00 , where t cc , t nc , and t nn are real constants, and p and k are momenta perpendicular to the interface [57]. In t 00 , δ 00 represents the tunneling between zero momentum elements. In addition, the factor δ p k in the coupling between non-condensation elements originates from the momentum conservation along the interface. Our interest is the dominant contribution where k = 0 component of lesser Green's function is concerned. Then, the momentum conservation along the interface leads to the fact that the three dimensional momentum summation appearing in the current expression is reduced to one dimensional one perpendicular to the interface. This actually simplifies the problem and allows us to obtain the DC Josephson current expression of the planar junction in an analytic manner.
At the same time, the generic expression in the planar junction is still complex due to the presence of the three parameters t cc , t nc , and t nn . From a practical point of view, furthermore, an order-by-order analysis in tunnel couplings may be sufficient to discuss transport properties of the junction system. Therefore, we turn to obtain the DC Josephson current of the planar junction, which contains first few series in tunneling amplitudes. To be specific, we consider the current expression up to fourth order in tunneling amplitudes. For this purpose, we can use the following approximations: p,RR (ω).
By substituting above into Eqs. (35)-(37), we obtain Here, Φ and Ψ are same as ones used in the point contact case (see also Appendix A). We note that the first two terms of the right hand side in Eq. (51) corresponds to the linear response theory and are consistent with [31] [58].
The terms proportional to t 2 nc t nn and t 2 nc t 2 nn are the new results obtained owing to the nonlinear response analysis in this work. Finally, we point out that it is also possible to obtain a current expression including higher order in t pk by incorporating higher tunneling processes inT LR(RL) andĜ LR .

V. DISCUSSIONS
By using the tunneling Hamiltonian and Bogoliubov theory, together with the Keldysh formalism, we have derived the DC Josephson current expression of a weaklyinteracting BEC. We applied the derived formula to quantum point contact and planar junction systems and discussed how the current-phase relations are modified from the conventional sinusoidal Josephson relation as increasing the tunnel coupling.
The formulation presented in this work is based on the previous work [32], where AC Josephson and direct currents in a BEC quantum point contact system have been examined under the presence of chemical potential and temperature biases. As technical differences, this work uses Green's functions without the phonon approximation and takes into account the momentum dependence in the tunnel coupling, which allow to obtain more general current expressions including the planar junction. In addition, a connection between AC and DC Josephson currents can be understood as follows. In the presence of a chemical potential bias, oscillating current as a function of time that contain both sine and cosine components is generated. As shown in [32], the amplitudes in the cosine components vanishes in the small bias limit while those in the sine components remain nonzero, which is also consistent with the AC Josephson current in a fermionic superconductor. Then, the zero bias limit of the AC Josephson current must be converged to the DC Josephson current by replacing the time oscillating argument with the relative phase.
Our formulation can also be applied to the so-called internal Josephson effect [16], where a two-component BEC with an internal coupling is concerned. In such a system, the coupling with a Rabi frequency Ω plays a role of the tunneling term. In the typical situation that the spatial dependence of Ω is negligible, however, the above coupling term only produce the momentum-conserving process between ↑ and ↓. Thus, the dominant DC Josephson current in this situation can be obtained by t nc = t nn = 0 in the planar junction case. Due to the absence of the momentumnonconserving processes, the renormalization of the tunnel coupling does not occur. Thus, it turns out that the internal Josephson current is described by the conventional sinusoidal curve, which is consistent with the Gross-Pitaevskii analysis [15]. At the same time, when Ω in Eq. (52) contains a nonnegligible spatial variation, higher-order tunneling processes are generated. Recently, such a coupling term is realized in a mixture system of 1 S 0 and 3 P 0 with Yb atoms [59]. If the similar system can also be realized in bosonic mixtures, one may harness Eqs. (35)-(37) obtained in this paper.
We also point out that the tunneling Hamiltonian approach can also be applied to the heat current. The heat current operator is defined as The commutator above consist of two parts: one between kinetic energy and tunneling terms and one between interaction energy and tunneling term. Especially, the latter part leads to a term quartic in the field operators whose analyses beyond the linear response regime is not so simple. To obtain a more convenient expression, we consider the Heisenberg equation The explicit evaluation of the above commutator shows that the heat current operator can be rewritten as [49] As in the case of the particle current, the average heat current is expressed in terms of lesser Green's function as follows: Since the form ofĜ < has already been obtained in the analyses of the particle current in this paper, it is straightforward to confirm that the heat current originating from the coupling associated with the condensation elements is obtained as where I odd N and I even N are given by Eqs. (36) and (37).
This physically means that the tunneling processes associated with condensation elements do not carry heat. Thus, if any, a nonzero heat current comes from the Bogoliubov-mode tunneling, which contains the Lee-Huang-Yang factor √ va 3 with the s-wave scattering length a [31,32] and therefore is expected to gives rise to a small effect in a weakly-interacting BEC [60].
So far, we have analyzed the two-terminal DC Josephson current where each reservoir is assumed to be uniform. By using box potential techniques available in ultracold atomic gases [61], the current in such a geometry can directly be measured and in fact has been measured in Ref. [24]. At the same time, by taking into account the fact that all experiments do not always implement such a scheme, we briefly comment on trapped systems. Usually, spatial inhomogeneities in ultracold atomic gases can be explained with the local density approximation [62]. If the spatial inhomogeneity in each reservoir is indeed enough smooth that the local density approximation is available, the analyses presented in this paper can be applied to trapped systems. In the case of the point contact system, the current is purely expressed in terms of local Green's functions at the center of the trap. Thus, Eq. (44) can directly be applied to the trapped case, provided that µ or v is now one at the center of the trap. Such an analysis has indeed been used in Ref. [48]. In the case of the planar junction system, on the other hand, one may locally introduce the hopping amplitude t pk (r). As discussed in [46], the total current should then be determined from the spatial integration of the corresponding local current.
Finally, we mention possible experiments to measure the DC Josephson current. In ultracold atomic gases, such a current can be measured by making current-biased junctions [37], which have been realized in a bosonic system [7] and in a molecular BEC system [23,25]. There, an external current in current-biased junctions is generated by shifting a trapping potential [7] or potential barrier [23,25]. It is then interesting to see the effects of higher order tunneling processes, which have already been done in a superconducting quantum point contact in condensed matter [63]  Here, we derive formulae of uncoupled Green's functions used in the main text.
Due to the absence of the tunnel coupling, uncoupled Green's functions are diagonal in the momentum index g pk =ĝ p δ pk . When it comes to retarded and advanced Green's functions, moreover, they are diagonal in the reservoir indexĝ αα δ αβ , since the commutator between the condensation elements k = 0 vanishes. Thus, the retarded Green's functions in the frequency space within the Bogoliubov theory is given bŷ where · · · 0 means an average without H T . We note that since the same chemical potential between reservoirs is assumed,ĝ R does not depend on α. Similarly, advanced Green's function can simply be obtained by using the relation g A k,αα (ω) = [ĝ R k,αα (ω)] † . In contrast, lesser Green's function contain an off-diagonal component in the reservoir index, which is associated with the presence of a BEC. By using the Bogoliubov theory, lesser Green's function is determined with the Bose distribution function n α (ω) = 1 e ω/Tα −1 . We note thatĝ < 0,αβ takes the same expression irrespective of reservoir indices α, β, and the factor δ(ω) appearing in k = 0 component represent the fact that the condensation elements do not evolve in time [36].

Quantum point contact case
In the case of quantum point contact transport, uncoupled Green's functions that take the sum over momentum index are concerned.
We note that they are essentially local Green's functions at origin since Since local Green's function may contain a UV divergence, we introduce a momentum cutoff Λ to regularize the problem. The introduction of the momentum cutoff is physically sound in that after all the tunneling Hamiltonian describes low-energy transport of the system. Then, local regarded Green function at the frequency space is given bŷ Here, Φ and Ψ are real parts of retarded Green's function that are given by withω = ω/µ. In addition, ρ and σ are the density of states associated with the imaginary parts of retarded Green's function. For ω > 0, they are determined as where Again, advanced Green function can be determined from g A (ω) = [ĝ R (ω)] † . Finally, uncoupled lesser Green function in frequency space is obtained aŝ

Planar junction case
For calculation of the dominant DC Josephson current in the planar junction, it is important to evaluate the following form of retarded or advanced Green's function: where k is the momentum perpendicular to the interface and k = 0 is assumed above. Sinceĝ

R(A)
−k,αα (ω) = g R(A) k,αα (ω), we see that the similar calculation to one in the point contact case is allowed. A straightforward calculation leadŝ where g R 1 , g R 2 , and f R are introduced in Eq. (A3). Thus, retarded and advanced Green's functions relevant to the dominant DC Josephson effect in the planar junction turns out to be obtained by multiplying ones in the point contact by 2π. In this appendix, we show a derivation of the DC Josephson current in conventional BCS superconductors by means of the tunneling Hamiltonian, which has already be done by using the fluctuation-dissipation relation in the total system [64]. For future reference, here we give another derivation without use of the fluctuationdissipation relation in the total system, which allows us to extend the system with a temperature imbalance between reservoirs.
In what follows, we assume that the tunneling amplitude is momentum-independent, which is usually done in fermionic systems [26]. Then, the analysis is reduced to one used in the point contact system. Namely, we can consider the following tunneling term where we introduce the Nambu spinor and the tunneling matrix with the relative phase ∆φ. We note that ψ ↑/↓ obeys the fermionic anti-commutation relation. The presence of e ±i∆φ/2 in the above tunneling matrix comes from the fact that Cooper pairs carry the phase φ L(R) and fermions carry half of it. By introducing local lesser Green's func-tion, G < LR (τ, τ ) = i ψ † ↑,R (0, τ )ψ ↑,L (0, τ ) ψ ↓,R (0, τ )ψ ↑,L (0, τ ) ψ † ↑,R (0, τ )ψ † ↓,L (0, τ ) ψ ↓,R (0, τ )ψ † ↓,L (0, τ ) , (B4) the particle current is expressed as I N = −Tr σ z tĜ < LR (τ, τ ) −t †Ĝ< RL (τ, τ ) . (B5) We now evaluate I N within the BCS theory. Then, uncoupled retarded (advanced) Green's function in the frequency-space representation is determined as [64] g R αα (ω)= whereĝ A αα (ω) = [ĝ R αα (ω)] † , W = 1/[πρ( F )] with the density of state at the Fermi energy ρ( F ), and ∆ α is the superconducting energy gap in reservoir α. We note that in order to obtain the expression above, we use the approximation that the density of states in free fermions is replaced by a constant value determined at the Fermi energy. Similarly, uncoupled lesser Green's function in the frequency space is obtained from the fluctuationdissipation relation in each reservoir, with the Fermi distribution function n F,α (ω) = 1 e ω/Tα +1 . We note that the above usage of the fluctuationdissipation relation is limited within each reservoir and we do not assume that the total system is in equilibrium.
Based on uncoupled Green's functions introduced above, we can obtain Dyson equations relevant to the Josephson current. Dyson equations for retarded and advanced components are given bŷ On the other hand, Dyson equations for the lesser component arê