{'fraction_non_alphanumeric': 0.032169811320754714, 'fraction_numerical': 0.011721698113207547, 'mean_word_length': 4.465455014178912, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnetic field oscillations of the critical current in long ballistic graphene Josephson junctions

We study the Josephson current in long ballistic superconductor-monolayer graphene-superconductor junctions. As a first step, we have developed an efficient computational approach to calculate the Josephson current in tight-binding systems. This approach can be particularly useful in the long junction limit, which has hitherto attracted less theoretical interest but has recently become experimentally relevant. We use this computational approach to study the dependence of the critical current on the junction geometry, doping level, and an applied perpendicular magnetic field. In zero magnetic field we find a good qualitative agreement with the recent experiment of Ben Shalom et al. (Reference[12]) for the length dependence of the critical current. For highly doped samples our numerical calculations show a broad agreement with the results of the quasiclassical formalism. In this case the critical current exhibits Fraunhofer-like oscillations as a function of the magnetic field. However, for lower doping levels, where the cyclotron orbit becomes comparable to the characteristic geometrical length scales of the system, deviations from the results of the quasiclassical formalism appear. We argue that due to the exceptional tunability and long mean free path of graphene systems a new regime can be explored where geometrical and dynamical effects are equally important to understand the magnetic field dependence of the critical current.

The theoretical work has mainly focused on short SGS junctions to date [32][33][34]41,43,44 , where the length of the normal region L is smaller than the superconducting coherence length ξ 0 = vF ∆0 . In addition, it was usually assumed that the width W of the junction is much larger than L. Although the long junction regime has been studied theoretically for superconductor-normal metalsuperconductor (SNS) systems [45][46][47][48][49][50][51] , the physics of long SGS junctions is less explored. An experimental study of long (L ξ 0 ) and wide (W ≫ L) diffusive SGS junctions was presented in Reference 25. However, in recent experiments different transport regimes have become accessible, where ballistic propagation was achieved in graphene samples where L ξ 0 10,12 and/or W/L ≈ 1 10 . Furthermore, the dependence of the superconducting critical current I c on a perpendicular magnetic field B has also been measured [10][11][12] in these SGS junctions. While References 11,12 have found that the oscillations of I c as a function of B can be described, at least in doped samples, by a Fraunhofer-like interference pattern, in Reference 10 deviations from the Fraunhofer pattern have been observed for samples that are in the long junction limit and have an aspect ratio W/L ≈ 1. Previously, deviations from the Fraunhofer-like I c (B) dependence were also observed in SNS junctions both in the diffusive 52 and in the quasi-ballistic limit 53 , and the subsequent theoretical work have elucidated the role of the junction geometry [54][55][56] using the quasiclassical Green's function approach. It is not immediately clear, however, if these theoretical results are directly applicable to SGS junctions, especially in the low doping regime.
Our aim in this work is twofold. Firstly, we want to present a newly developed computational approach to calculate the Josephson current in tight-binding (TB) systems. The method is general and can be implemented for many TB systems, not only for graphene. It takes into account on equal footing the contributions coming from both the Andreev bound states (ABS) and the scattering states (ScS), the latter being especially important in long Josephson junctions, where it is known that cancellation between different supercurrent contributions occur 47,57 . Since our method accounts for both contributions, it can be used for efficient simulations of recent experimental systems [10][11][12] . Secondly, using the above computational method, we study the length L and magnetic field dependence of the critical current in long SGS junctions. Although the length dependence of I c has been studied before using various theoretical approaches 28,30,58,59 we revisit this question because the recent observations in Reference 12 offer the possibility to directly compare theory and experiments. Encouragingly, we find a good qualitative agreement between our results and the experimental observations of Reference 12 , indicating that our approach can capture important aspects of the physics of long SGS junctions. Regarding the magnetic field effects in long SGS junctions, to our knowledge no detailed study is available at present. We study the magnetic field oscillations of I c as function of the doping of the normal region. For high doping we find that the semiclassical formalism [54][55][56] , developed for ballistic Josephson junction where the normal region is a two-dimensional electron gas, can also describe the oscillations I c in SGS junctions. However, for lower doping, where the cyclotron radius R c becomes comparable to W and/or L, orbital effects can no longer be neglected and deviations from the quasiclassical results appear.
The paper is organized as follows. In Section II we briefly introduce the model system that we used in our calculations. In order to make the paper accessible to a broad audience, this is followed by the presentation of our main results for the critical current I c . First, in Section III we discuss the length dependence of I c and also the current-phase relation. The effect of the magnetic field on I c is treated in Section IV. A general numerical approach to calculate the Josephson current for TB Hamiltonians is presented in Section V, while some of the relevant details of the TB model used in this work is given in Section VI. Finally, we conclude with Section VII.

II. THE MODEL
We first briefly describe the model we employed to calculate the Josephson current in SGS junctions, further details can be found in Section VI. In the normal conducting region of length L and width W we use the nearest-neighbour TB model of graphene 60,61 with the HamiltonianĤ Here ε i is the on-site energy on the atomic site i, γ ij is the hopping amplitude between the nearest neighbour atomic sites ij in the graphene lattice, and c † i (c i ) creates (annihilates) an electron at site i. The magnetic field can be incorporated by means of the Peierls substitution 62 : where φ 0 = h/e is the flux quantum, A(r) denotes the vector potential and the vector R i points to the ith atomic site in the lattice. The spatial dependence of A(r) is such that it yields a homogeneous perpendicular magnetic field B = (0, 0, B z ) T in the normal region and zero field in the superconducting regions, see Section VI C. The superconducting regions are modelled by a highly doped graphene region 43 of width W and open boundary condition in the transport direction. It is assumed that a finite on-site pair-potential ∆ 0 e iϕL,R is induced by proximity effect in the left (L) and right (R) electrodes. We note that our methodology would allow for other models of the superconducting regions as well 63 . For the superconducting pair-potential we assume a step-like change at the normal-superconductor (NS) interfaces: ( Here ϕ L (ϕ R ) is the phase of the pair-potential in the left (right) lead and we will denote by δϕ = ϕ R − ϕ L the phase difference. Our main interest in this work is to study SGS junctions where there is a significant difference between the doping levels of the S and N regions: is the Fermi wavelength in the superconductor (normal) region. Moreover, the junction is long L/ξ 0 1 with respect to the coherence length ξ 0 . In this case we expect that the detailed spatial dependence of ∆ 0 in the vicinity of the normal-superconductor interface is not very important and therefore the above approximation should give qualitatively correct results. Indeed, References 28,29,31,64,65 have shown that the selfconsistent calculation of ∆(x) in clean SGS junctions is most important for (i) short junctions, (ii) no Fermi-level mismatch between the S and N regions, (iii) temperatures close to T c . The simulation of realistic samples on the micrometer length scale is quite challenging in the TB framework due to the huge number of the atomic sites. Part of the problem can be circumvent by using an efficient numerical approach, see Section VI A for details. Moreover, we expect that experimentally relevant informations can be extracted from TB systems that follow certain scaling laws but imply significantly lower computational costs. Such an approach has proved to be very useful recently in the calculation of normal transport 66,67 for mesoscopic graphene structures. We expect that as long as the characteristic dimensions W and L of the system are much larger than the lattice constant of graphene, the same physical behavior should be observed in systems with the same L/ξ 0 , W/L, 2W λN , ξ0 λN , T Tc and L/l B control parameters, where λ N = EF vF is the Fermi wave number in the normal region l B = |eB| is the magnetic length and T c ≈ ∆ 0 /(1.76k B ) is the critical temperature. We are interested here in the bulk properties of the supercurrent, i.e., we need to ensure that edge effects do not play role. In most of our calculations we used zigzag nanoribbons, however, we have checked that we obtain very similar results for armchair nanoribbon as well. (see Section III). Therefore we expect that our results would not change for more general edges either.

A. Length dependence of Ic
An important property of long Josephson junctions is the dependence of the critical current I c on the junction length L, which was measured recently in Reference 12 . In general, at zero temperature I c is given by the relation 68 where N is the number of open channels in the normal region: N = EF W π vF = W kF π = 2W λF . The dimensionless coefficient α can depend on a number of factors, such as the junction transparency, the presence of a p − n junction or other disorder induced by the superconducting contacts, the doping level in the normal region 43 etc. However, for the case considered here (no disorder at the SG interface) and in the long junction limit α is expected to be a function of the ratio L/ξ 0 only 46,47 .
First, in Figure 1(a) we present the zero temperature calculations for α. The good agreement between the calculations for different widths and edge types indicates that these results are free of finite size and edge effects. One can see that the numerical results can be fitted with a function A/(L/ξ 0 + C), where A and C are fitting parameters. For L/ξ 0 ≫ C the critical current falls off as 1/L, whereas in the short junction limit L/ξ 0 ≪ 1 it goes to a finite value and reproduce the analytical prediction of Reference 43: for L/ξ 0 → 0 and λ S ≪ λ N the value of α is ≈ 1.22. Next, in Figures 1(b)-(d) we show the I c (L) curves for low, but finite temperatures T . As the temperature is increased from T = 0, both fit parameters A and C decrease and eventually there is a temperature range where the parameter C becomes very small so that I c falls off as ∼ 1/L for L/ξ 0 . Such a ∼ 1/L dependence of I c was recently observed in Reference 12 . It was suggested that this peculiar I c (L) dependence is a signature of the SGS junctions being truly ballistic and in the long junction limit. Our calculations support this conclusion, nevertheless, it would be interesting to measure I c (L) at lower temperatures in order to map out experimentally the T dependence of the parameters A(T ) and C(T ) and compare it to our results. It is important to note the following: (i) Taking the T c of bulk niobium, the experimental data of Reference 12 correspond to T /T c = 0.045, i.e., to lower temperature than shown in Figure 1(b). However, if one assumes that the induced pair potential in the graphene is smaller than the bulk value of ∆, i.e., it corresponds to a smaller T c , then the agreement with our calculations becomes better. (ii) The values of α in Figure 1(c) significantly overestimate the experimental ones reported in Reference 12. We think this is due to the fact that we have assumed perfectly transparent SG interface and no disorder in the normal part of the junction, while in Reference 12 the SG interface had a finite transparency and a p-n junction was probably formed due to the doping effect of the contacts. We leave the study of these effects for a future work. Finally, as shown in Figures 1(d), as the temperature is further increased such that the energy scale k B T becomes non-negligible with respect to the ballistic Thouless energy vF L 46,48 , I c becomes exponentially suppressed: I c ∼ e −cL/ξ0 , where c is a fitting parameter.

B. Current-phase relation
Regarding the current-phase relation (CPR), it has long been known that the CPR in long SNS Josephson junctions [45][46][47] is substantially different from the ∝ sin(δφ) relation valid in the short junction limit 69 . In the one-dimensional TB model studied in Reference 70 the deviation from the ∝ sin(δφ) dependence was explained by the contribution of the scattering states to the total current which becomes comparable to the contribution of the Andreev bound states as the length of the junction increases. A partial cancellation effect between the different contributions to the Josephson current in SNS junctions was also pointed out by References 46 and 47. Recently, anharmonic CPR was found in the calcu-lations of References 28,29,31 for SGS junctions. We have also calculated the CPR using our method for a long SGS junction. As one can see in Figure 2(a), at zero temperature the CPR deviates significantly from the ∝ sin(δφ) dependence and for δφ/π 0.3 it is a linear function of δφ/π. The contribution coming from the ScS (dashed- dotted line) is of the same magnitude as the contribution of the ABSs (dashed line). In Figure 2(a) the two contributions have the same sign, however, this is not always the case: calculations not shown here indicate that depending on the L/ξ 0 ratio the supercurrent due to the ScS can be either positive or negative. For finite temperatures [see Figure 2(b)], similarly to SNS junctions 46 , the CPR acquires a harmonic dependence on δφ. In these wide ribbons the major characteristics of the CPR do not depend on whether zigzag or armchair nanoribbons are used in the calculations.
To briefly summarize the results presented in this section, we find that our results show a good qualitative agreement with i) the measurements of Reference 12 for the length dependence of I c , ii) with previous theoretical works 29,31 for the CPR, even though we do not calculate the pair potential ∆ self-consistently. This suggests that in long SGS JJs, if there is a finite Fermi-level mismatch at the SG interface, the self-consistent calculation of ∆ is less important than in the short junction limit with no Fermi-level mismatch.

IV. OSCILLATIONS OF THE CRITICAL CURRENT IN PERPENDICULAR MAGNETIC FIELD
We now turn to the properties of the critical current in the presence of a perpendicular magnetic field B. Oscillations of I c have been measured in several recent experiments 10-12 but have not yet received much theoretical attention. It is well known that in tunnel junctions I c exhibits Fraunhofer-like oscillations as a function of the piercing magnetic flux Φ with a period φ 0 = h/2e and oscillation minima at integer multiples of the flux quantum φ 0 : , where I c (0) is the zero-field critical current. The long junction limit in SNS systems was studied in References 54-56,71 using the quasiclassical Green's function formalism. It was pointed out that the magnetic field oscillation of I c also depend on the geometry of the junctions. For wide and long ballistic junctions, where W ≫ L ≫ ξ 0 , the critical current oscillates as 55,71 {x} denoting the fractional part of x. The oscillation pattern given by Eq. (5) is very similar to the Fraunhofer-like oscillations in tunnel junctions, except for Φ/φ 0 ≪ 1. However, deviations from Eq. (5) were observed in the measurements of Reference 53 in quasi-ballistic SNS junctions. Subsequent theoretical work 54,55 showed that geometrical effects become important when W ∼ L ≫ ξ 0 . In particular, Reference 54 found that the periodicity of I c as a function of magnetic field changes from φ 0 [see Eq. 5] to 2φ 0 as the flux through the junction increases and at low temperatures the crossover to the 2φ 0 periodicity appears at a flux ∼ φ 0 W/L. Regarding the recent experiments in SGS junctions, Fraunhofer-like pattern for I c (B) was found in Reference 12 for wide (W ≫ L) junctions, whereas Reference 10 measured periodicity that was longer than φ 0 for junctions with aspect ratios W/L ≈ 1. An important assumption behind the quasiclassical formalism [54][55][56]71 is that the effect of magnetic field can be taken into account through a phase factor that the wave function of the quasiparticles acquires along classical trajectories that are straight lines. This is a good approximation when the cyclotron radius of the particles is much larger than other characteristic length scales in the system.
However, compared to traditional SNS systems, the exceptional tunability of the state-of-the-art graphene systems, combined with the use of Nb 12 or MoRe 10 superconductors allows, in principle, to reach a regime where the size of the cyclotron radius R c is comparable to W and/or L already for relatively low magnetic fields such that the effect of magnetic field on the superconducting electrodes can still be neglected. Considering the semiclassical cyclotron radius R c in graphene R c = EF l 2 B vF one can see that for a charge density n e ≈ 10 12 1 cm 2 the cyclotron radius is ≈ 1µm for a B ≈ 0.1T . Thus R c becomes an important length scale when the chemical potential approaches the charge neutrality point. This regime is expected to harbour rich physics, because geometrical effects, discussed above, and dynamical effects related to the curved semiclassical trajectories of the quasiparticles are equally important. It is then not clear if the quasiclassical formalism 54,55,71 , which neglects the dynamical effects, can still give a good description. The possible importance of this regime has also been discussed in Reference 12.
In the following we take the chemical potential µ N in the normal region as a tuning parameter discuss separately the high doping limit, when R c ≫ L, W holds for the considered magnetic field range, and the low doping range, where R c L, W can be reached for relatively small magnetic fields. In all our calculations we assume that the effect of magnetic field on the pair potential in the superconducting contacts can be neglected.

A. High doping limit
To see if the quasiclassical theory 54,55 also applies to SGS junctions when the normal region is strongly doped, we calculate I c (B) for different W/L ratios. We have found that the oscillations depend only weakly on the exact value of L/ξ 0 and on the temperature (in the low temperature limit relevant for recent experiments [10][11][12] ). Thus, we present our results only for particular temperature and L/ξ 0 values. Figure 3 shows the magnetic oscillations calculated for different aspect ratios W/L 72 . These results were obtained for zigzag nanoribbons but we have checked that the results do not change qualitatively for armchair nanoribbons. For wide junctions [ Figure 3(a)], we recover the Fraunhofer-like pattern of the oscillations with minimums at integer multiples of the flux quantum φ 0 , see Eq. (5). Deviations from the ideal curve Eq. (5) start to appear only for magnetic fluxes Φ 4φ 0 . As the aspect ratio W/L decreases [ Figures 3(b)-(c)], the periodicity of I c becomes longer than φ 0 , even for smaller magnetic fields. Interestingly, for W/L ∼ 1 [ Figure 3(c)] the oscillation period is roughly constant in the considered magnetic field range. Finally, in Figure 3(d) one can see that in narrow samples the first minimum in I c is at 2φ 0 and the current does not go to zero. These results are in broad agreement with the quasiclassical calculations of Reference 54 indicating that I c (B) in highly doped graphene samples can essentially be described by the theory used previously for SNS junctions in References 54, 55, and 71. As mentioned above, deviations from the Fraunhofer-like pattern, similar to the ones shown in Figures 3(c) and (d) have been measured recently in Reference 10 for samples with W/L ≈ 1. However, a more quantitative comparison of our calculations with Reference 10 is difficult because i) due to the inevitable disorder at the edges, the effective width of the samples may be smaller than the geometrical width, and ii) the current oscillations are likely to depend on the properties of the p-n junction formed at the SG interface which is not taken into account in our calculations.

B. Low doping limit
We now consider the case when the doping of the normal region approaches the charge neutrality point, but it is still large enough so that the formation electron-hole puddles can be neglected. The most interesting results in this parameter range are shown in Figure 4. As in the high doping case, the current oscillations depend on the aspect ratio W/L. For wide junctions, as in Figure 4(a), the oscillations are similar to the strong doping case shown in Figure 3(a) but the oscillation period is somewhat longer and the amplitude, relative to I c (0), is larger. For L > R c , however, the current is strongly suppressed. In Figures 4(b)-(d) we calculate the current for a narrower junction as the doping is decreased. Comparing Figure 4(b) and Figure 3(c), where the same geometrical parameters were used, one can see that for smaller doping dynamical effects influence the period of the current oscillations. (As indicated in Figure 4(b), the ratio L/R c ≈ 1 is obtained for Φ/φ 0 ≈ 7.8 in this case. ) While in the strongly doped case [ Figure 3(c)] the period of oscillations is ≈ 1.1φ 0 in the shown magnetic field range, for smaller doping the distance between consecutive minima grows with the magnetic field. As the doping level is further reduced the oscillations of I c become rather complex for magnetic fields where L > R c is fulfilled, see illustrates that geometrical effects (due to W ∼ L) and dynamical effects (R c ∼ L) can both strongly affect the oscillation pattern. Note, that for the parameters used in Figure 4(b)-(d) the diameter of the cyclotron motion is still smaller than the geometrical parameters, i.e., 2R c > W, L and therefore no quantum Hall edge states are formed.
The regime R c L has recently been considered in Reference 12, where the sample dimensions were W ≫  , the situation is somewhat different from the wide junction case because here the electron-hole trajectories cannot drift away to large distances and are more likely to form (nearly) closed orbits and hence bound states. This may explain why the supercurrent is not suppressed.

V. TIGHT-BINDING APPROACH TO CALCULATE THE DC JOSEPHSON CURRENT
In this section we describe the TB approach we used to calculate the DC Josephson current. It is a generalization of the approach developed to calculate two-terminal normal transport 73,74 to the case when both terminals are superconducting and the main interest is not the scattering matrix (and the differential conductance) but the current flowing due to the superconducting phase difference between the terminals. As already mentioned, the method is general and can be implemented for many TB system. It is a generalization of the one-dimensional TB work by References 70 and 75 to nanoribbon geometries and would also allow for extension to multi-terminal systems studied, e.g., in References 36, 37, and 76. The actual calculations presented in Sections III and IV were performed with the EQuUs software 77 .
A. The general setup The studied system, including the central region and the electrodes, is schematically depicted in Fig. 5. The Let us label the eigenvalues and eigenstates of the Hamiltonian (6) by E n and Ψ n , respectively. The index n labels both the bound and scattering states that are formed in the system. In the latter case n = (m, k) stands for a pair made of the discrete transverse quantum number m and the wave vector k describing a propagating incoming state in one of the leads. Generally, the Green's function of the studied system can be written as: The normalization of the eigenstates Ψ n in Eq. (7) is straightforward for the bound states, namely Ψ n1 |Ψ n2 = δ n1n2 . The scattering states, on the other hand, are normalized to unit incoming current 73,78 .

B. The expectation value of the current operator
Due to current conservation, the charge current is equal between any pair of slabs of the studied system. Therefore, as we will show, it is sufficient to calculate the Green's function only on a couple of slabs, leading to a numerical efficient method. To this end let us now introduce the operator P i projecting on the subspace of the ith slab of the system. The matrix form of the projector P i reads where I Ni is an N i × N i unity matrix. One can notice that the projectors P i = P i † are hermitian operators. In addition, the operator P {i} projecting on a set of slabs (i 1 , i 2 , . . . , i p ) can be given as a sum of the individual projectors Thus, the projected Green's function on the slabs (i 1 , i 2 , . . . , i p ) can be given as G P (z) = P G(z)P † . Using the Green's function given in Eq. (7), one would obtain the following projected Green's function are the projected wave functions normalized to unity by the factor Ω n = Ψ n |P |Ψ n . The expectation value of the current operatorÎ with respect to a projected state Φ n can be calculated as follows: Here E min and E max are arbitrary energies such that E max − E min ≤ BW , where BW is the bandwidth of the system. To obtain the total current I c,tot , we will need to sum over all projected states Φ n . Using the projected Green's function G P introduced in Eq. (10), one finds that Green's function. We remind, that the retarded Green's function may have poles corresponding to bound states which would lead to complications in the numerical evaluation of Eq. (13). Therefore, making use of the fact that retarded Green's function is analytical in the upper halfplane (Im(z) > 0) of the complex energy, the integral can be performed along a path Γ in the complex plane as shown in Fig. 6. Thus, the sum of the expectation values of the current operator between energies E min and E max can be calculated by the formula The total current I c,tot would also depend on the thermal occupation of the electronic states. If we have normal conducting terminals then I c,tot can be expressed as (15) where f (z) is the Fermi distribution function. In the case of superconducting terminal the situation is somewhat more complicated but we leave the discussion of superconductor-normal-superconductor systems to Section V C. Here we only note that for the evaluation of Eq. (14) one does not need to know explicitly the specific energies E n of the bound states. Choosing E min and E max appropriately the contributions of both the scattering and bound states are automatically taken into account. However, one has to choose such a contour Γ which does not enclose "unwanted poles", e.g., in the case of Eq. (15 it avoids the poles of f (z) located at the energies Z l = i(2l + 1)πk B T .
As mentioned, due to current conservation it is sufficient to calculate I c,tot between two arbitrary slabs, and for practical reasons we choose the surface slabs of the central region (see the 1st and the N th slabs in Fig.  5). To calculate the necessary projected Green's function G P (z), we use the Green's function technique of Refs. 73,78 . Namely, we account for the effect of the leads attached to the normal region by means of the Dyson's equation: 79 Here Σ L and Σ R are the self-energies of the left and right leads, respectively (see Section VI D for further details). H eff (z) is the effective energy dependent Hamiltonian describing the surface of the normal region (see Fig.5). (The energy dependence of H eff (z) can be interpreted as the effect of the inner sites located between the surface slabs 1 and N .) The effective Hamiltonian H eff (z) can be obtained via several methods. For example, one can eliminate all the sites inside the central region by the decimation method and keep only the sites of the surface slabs 74 . However, for long ballistic structures there is a more efficient method which we will briefly describe in Section VI A.

C. SNS system
The discussion in Sections V A and V B was general and would apply regardless of whether one assumes normal or superconducting leads. In this section we discuss those aspects of the problem which are specific to normal-superconductor-normal (SNS) systems, i.e., systems where there are two superconducting terminals and a central normal scattering region.
We describe this inhomogeneous superconducting system by the Bogoliubov de Gennes (BdG) model. Consequently, the Hilbert space of the superconducting system is constructed as the product of the Hilbert space of the normal system and the Nambu space describing the electron (u) end hole like (v) degrees of freedom. The matrix elements ofĤ in Eq. (6) can be written as and where, for simplicity, we omitted the indexes labeling the slabs of the system. Matrices H u 0 and H u ±1 describe the electron likes components. The Hamiltonians of the hole like components are given by Finally, the diagonal matrix ∆ contains the superconducting pair-potentials on the atomic sites of the slabs, and µ is the chemical potential. Since in the central region the superconducting pair potential is zero, the electron and hole like components of the BdG equations become uncoupled and H eff becomes diagonal in the Nambu space. Similarly, if one calculates the current in the central (normal region), the current operator can be written in a block-diagonal form whereÎ u andÎ v are the charge current operators of the electron and hole like states, respectively (their explicit form for our TB model is given in Section VI B ). The projected Green's function (see Eq. 16), on the other hand, has a matrix form where G P uv and G P vu are non-zero due to the fact that electron and hole components are coupled in the selfenergies of the superconducting leads.
Taking into account the thermal occupation of the electron-like and hole-like states, we obtain the following expression for the charge current: In general, the spectrum of the BdG Hamiltonian is symmetrical around E = 0 and therefore it is enough to consider either E > 0 or E < 0. Considering only the negative energies, the spectrum of the SNS junction can be divided into three spectral regions 70 . The first region corresponds to the states of energy −∆ 0 < E n < 0. Due to the energy gap in the superconducting leads, these states are bounded to the normal region since they are decaying exponentially in the superconducting leads. We refer to these states as the Andreev bound states (ABS). The next energy regime is given by −|BW L | < E n < −∆ 0 , where BW L is the bandwidth of the superconducting leads, i.e., the maximal energy of the propagating states in the leads. These scattering states form a continuous energy range in the spectrum. Finally, in the case when the bandwidth of the normal region BW N is larger than BW L , we can define a third energy regime. Namely, for energies −|BW N | < E n < −|BW L | the states formed in the system are decaying exponentially in the leads, but are still propagating in the normal region. According to Ref. 70 , we refer to these states as the normal bound states (NBS). Figure 7 shows the integration contours to be used to calculate the DC Josephson current due to the ScS's, the ABS's and the NBS's, respectively. In addition, one can also calculate the contribution from all these states at once by integration over the contour Γ ALL also shown in Figure 7. Note, that using this approach it is not necessary to find the zeros of a polynomial, as in Reference 70. Finally, due to the fact that in superconducting systems the spectral density distribution is symmetrical with respect to E = 0, the contour integration (21) needs to be evaluated only in the Re(z) < 0 half-plane. The contribution of the states in the Re(z) > 0 half-plane can be accounted for by a factor of 2 in the final result.
This completes the general discussion of the TB method that we used to calculate the Josephson current. In the next section we discuss certain model specific aspects of our calculations.

VI. DETAILS OF THE TB CALCULATIONS
In this section we give the details of our TB calculations for SGS junctions which are relevant for obtaining the results presented in Sections III and IV.

A. Calculation of the effective Hamiltonian
In this section we provide an efficient numerical method to calculate the effective Hamiltonian H eff (z) needed to evaluate the expectation value of the current operator in Eqs. (15) and (21). While the procedure described here is optimized for a ballistic scattering region containing identical slabs, it has been shown in Ref. 66 that this approach can be generalized also to more complex geometries as long as the system is ballistic and is numerically more efficient than the standard recursive Green's function techniques.
We assume that the central region is made of identical slabs described by the Hamiltonian H i 0 ≡ H 0 and coupled to each other by H i∓1,i ±1 ≡ H ±1 (1 ≤ i ≤ N ). Following the procedure described in Ref. 66 , we can obtain the H eff (z) using the Green's function g i,j of an infinite ribbon made of these slabs (here i, j are slab indices). The Green's function g i,j can be efficiently calculated using a semi-analytical formula introduced in Ref. 73 . In order obtain the elements of H eff (z) we need to calculate the propagators g i,j on slabs i ∈ {0, 1, N, N + 1} and between these slabs. They can be arranged into a matrix that reads Since the structure of the ribbon contains only nearest neighbor couplings between the slabs without long range interaction, the effective Hamiltonian defined as H(z) = zI − G(z) −1 has the following structure: Note that there is no coupling between the slabs 0 and N since these slabs are coupled via the slab 1, and therefore the matrix element H 0N vanishes.
The effect of the sites between the surface slabs are incorporated within the energy dependence of H eff . We note that using the decimation method 74 one would obtain an identical effective Hamiltonian. However, the described procedure involves only sites that are located on the surface of the central region and therefore the computational cost of calculating H eff is scaling only with the width of the ribbon. This is especially important in the case of long scattering regions.

B. The effective current-operator
In order to evaluate the contour integral given in (21) one needs to construct the matrix representation of the current operator between the surface slabs of the central region. Since the electron and hole like components are uncoupled in the central region, one can obtain individual charge current operators for the electron (Î u ) and hole (Î v ) like components. Each of them can be derived, similarly to the normal systems, from the corresponding effective Hamiltonian H eff by means of a discretized continuity equation. Thus, we obtain the current operator between the slabs 1 and N from the elements of the effective Hamiltonian (24). Making use of the block diagonal structure of H eff , the current operator readŝ Similarly to the effective Hamiltonian, the effective current operator is also energy dependent.
C. Implementing the magnetic field in the SNS junction We now discuss how the magnetic field is taken into account in our calculations. This is needed in order to calculate the oscillations of the critical current, see Section IV. We limit our considerations to low magnetic fields where the effect of screening currents on ∆ can be neglected. Consequently, the only effect of the induced supercurrents on the surface of the superconducting regions is that the magnetic field is expelled from the superconducting leads, but is considered to be homogeneous in the normal region. The corresponding vector potential can be given in a Landau gauge where As we have seen in section VI A, one can calculate the effective Hamiltonian H eff efficiently for long ballistic systems if the normal region consists of identical slabs. This translational invariance would be broken by the vector potential given in Eq. (27). To avoid this problem one may try to incorporate the magnetic field in the normal region by a vector potential A x that is translational invariant in the normal region [see Fig. 8.(a)]: Note however that A x cannot be fitted continuously to the constant vector potential in the superconducting leads. However, A y (x, y) is related to A x (x, y) by a gauge transformation A y (x, y) = A x (x, y) + gradχ(x, y) with a gauge field given as Since the magnetic field enters the calculations through the Peierls-substitution (see Eq. 2), one can show that the effect of this gauge transformation on the effective Hamiltonian H eff can be expressed as where is a matrix describing a unitary transformation. Here χ(x, y) is defined by Eq. (29) and R i is a lattice vector.

D. Calculation of the self energies
Finally, we briefly discuss the calculation of the selfenergies Σ L,R that enter Eq. (16). We obtained them using the model of References 43 and 68, i.e., assuming that the superconducting leads consists of highly doped semiinfinite graphene ribbons where a finite superconducting pair potential ∆ was induced by proximity effect. The surface Green's functions of the leads and the corresponding self-energies Σ L,R can be calculated as described in Reference 78. We note that other approaches to calculate Σ L,R , such as the "bulk-BCS" model discussed in Reference 63, could equally be used in the computational framework we introduced.

VII. SUMMARY AND OUTLOOK
In summary, we have studied theoretically the DC Josephson current in long SGS junctions. We developed a theoretical framework that can be applied to an arbitrary superconducting-normal-superconducting junction defined on a tight-binding lattice. By treating the bound and scattering states on an equal footing it presents an accurate and efficient numerical method to calculate the equilibrium Josephson current.
We used this theoretical approach to investigate the dependence of the critical current on the geometrical properties of the junctions and on the magnetic field in the ballistic transport regime. In the zero field and low temperature limit we have found that the critical current decays as ∼ ξ 0 /L in agreement with recent measurements. For temperatures comparable to T c , on the other hand, the critical current becomes exponentially suppressed. Furthermore, we have found that in the long junction limit the contribution of the ScS to the Josephson current is as important as are the contribution coming from the ABSs.
We have also studied the magnetic oscillations of the critical current. Generally, for a given magnetic field one can distinguish the high-doping and the low-doping limits which are defined in terms of the cyclotron radius R c as R c W, L and R c W, L, respectively. In the high-doping regime, when R c is much larger than other length scales of the junction, the the period of the oscillations depends on the geometry of the junction in a similar way as predicted by previous quasiclassical calculations for S-2DEG-S Josephson junctions. For wide junctions, i.e. W/L > 1, we recover the Fraunhofer-like oscillation pattern, however, deviations from this oscillation pattern start to appear at a magnetic flux Φ ∼ W/L φ 0 resulting in an increased oscillation period. For narrow junctions, in turn, one can observe a more complex interference pattern in the magnetic field dependence of the critical current without reaching the zero value. Close to the Dirac point R c becomes comparable to the dimensions of the junctions, thus the dynamical effects, due to the curved semiclassical trajectories, cannot be neglected any more. According to our results, the interplay of these dynamical effects and the magnetic field induced quantum mechanical interference has twofold effect on the critical current. Firstly, by decreasing R c the oscillation period starts to increase and the minimums of the critical current are lifted from zero. Secondly, for R c smaller than the length of the junction, the magnetic dependence of the critical current does not show any regular oscillations. We note, however, that in this work we did not address the case 2R c < L, W , i.e., when the formation of quantum Hall states is expected.
The methodology that we have introduced here is quite flexible and it would allow to address a number of fur-ther problems. As a brief outlook, we would mention the following ones. Firstly, the interplay of Landau quantization and Josephson current flow in a weak link, i.e., the regime 2R c < L, W . This question is timely, as the first report on the observation of supercurrent in the quantum Hall regime in a SGS junction has recently appeared 80 . Secondly, although we have focused on ballistic junctions in this work, disorder effects in the normal region can also be incorporated using, e.g., the recursive Green's function technique. We mention two problems in this regard: in Reference 11 an anisotropic supercurrent distribution was found at low dopings, where the role of disorder is expected to be especially important. This was explained by calculating the normal density of states and assuming that it also determines the supercurrent flow. One could verify this assumption by solving the BdG equations for a disordered normal region and calculating the supercurrent distribution. Finally, we note that Josephson vortices were predicted 50,51 and later experimentally verified 81 to exist in diffusive SNS junctions. One may expect that they also appear in diffusive SGS junctions.