8 For v$_{F}$ k ≪ µ the this equation coincides with Eq. (A1) taken for µ ≪ Δ. Expanding it in small ˜ E$_{1}$ and ˜ E$_{2}$ we find The similar equations are obtained for the other state α = 3 which belongs to the same energy E$_{0}$ . We have j = − 2 e ∑ p [ ∂ ∂ k ( ˜ E 2 1 + ˜ E 2 ) d dE$_{0}$ tanh E$_{0}$ 2 T + v 2 $_{F}$p 2 E 3 0 v 2 $_{F}$[ p × k ] × p p 2 tanh E$_{0}$ 2 T ] H$_{αβ}$ = v$_{F}$ 2 ( ˆ u +(0) α ˆ σ · k ˆ u (0) β + ˆ v +(0) α ˆ σ · k ˆ v (0) β ) − µ ( ˆ u +(0) α ˆ u (0) β − ˆ v +(0) α ˆ v (0) β ) . Therefore H$_{11}$ = − H$_{33}$ = v$_{F}$ p · k 2 p − µv$_{F}$p E$_{0}$ , H$_{13}$ = − H$_{31}$ = iv$_{F}$ ([ z$_{0}$ × p ] k ) 2 p | Δ | E$_{0}$ , and which yields after integration over the momentum angle Subtraction of the normal-state current returns us to Eq. (42). j = − ev 2 $_{F}$k ∑ p [( 1 + | Δ | 2 E 2 0 ) d dE$_{0}$ tanh E$_{0}$ 2 T + v 2 $_{F}$p 2 E 3 0 tanh E$_{0}$ 2 T ] . (B14) 2. Undoped graphene The Bogoliubov–de Gennes equations can be solved ex- actly for the undoped case µ = 0 where the dispersion equation (20) becomes bi-quadratic with the energy spec- trum Eq. (21). We shall write down the solution of the Bogoliubov–de Gennes equations for the vector k parallel to the axis x ( k = k$_{x}$ ). Then Secular equations (B5), (B6) determine δE$_{1}$$_{,}$$_{3}$ together with the coefficients C (0) 11 , C (0) 13 . The first-order correc- tions to C$_{11}$ and C$_{13}$ are found from Eqs. (B1) and (B2) written up to the second-order terms. Using Eqs. (B9), (B10), and (B11) we obtain C (1) 11 = C (1) 13 = 0. The coef- ficients C (1) 12 and C (1) 14 are determined by Eqs. (B7) and (B8). In the same way we find C (1) 31 = C (1) 33 = 0 together with the coefficients C (1) 32 and C (1) 34 . Calculations for the states α = 2 , 4 can be done in exactly the same way. Using the obtained coefficients we can rewrite Eq. (5) for the current in the form H$_{21}$ = H$_{21}$ = H$_{43}$ = H$_{34}$ = − µ | Δ | E$_{0}$ , (B9) H$_{23}$ = − H$_{32}$ = − i v$_{F}$p E$_{0}$ v$_{F}$ ([ p × k ] z ) 2 p , (B10) H$_{41}$ = − H$_{14}$ = − i v$_{F}$p E$_{0}$ v$_{F}$ ([ p × k ] z ) 2 p . (B11) j = − ev$_{F}$ ∑ p [ A (0) 1 ( tanh E$_{0}$ − ˜ E$_{1}$ 2 T − tanh E$_{0}$ + ˜ E$_{1}$ 2 T ) − A (0) 2 ( tanh E$_{0}$ − ˜ E$_{2}$ 2 T − tanh E$_{0}$ + ˜ E$_{2}$ 2 T ) + 2 ( A (1) 1 − A (1) 2 ) tanh E$_{0}$ 2 T ] , (B12) where and A (0) 1 = ( | C (0) 11 | 2 − | C (0) 13 | 2 ) ˆ a † $_{↑}$ˆ σ ˆ a$_{↑}$ +2 uv ( C (0) ∗ 11 C (0) 13 − C (0) ∗ 13 C (0) 11 ) ˆ a † $_{↑}$ˆ σ ˆ a$_{↓}$ , A (0) 2 = ( | C (0) 22 | 2 − | C (0) 24 | 2 ) ˆ a † $_{↑}$ˆ σ ˆ a$_{↑}$ +2 uv ( C (0) ∗ 22 C (0) 24 − C (0) ∗ 24 C (0) 22 ) ˆ a † $_{↑}$ˆ σ ˆ a$_{↓}$ , where the phase factors are given by Using Eq. (B16) one finds the terms which determine the x component of the supercurrent: The current in Eq. (B12) takes the form A (1) 1 = ( u 2 − v $^{2}$) E − 1 0 H$_{23}$ ˆ a † $_{↑}$ˆ σ ˆ a$_{↓}$ , A (1) 2 = ( u 2 − v $^{2}$) E − 1 0 H$_{14}$ ˆ a † $_{↑}$ˆ σ ˆ a$_{↓}$ . e iφ ± = √ | Δ | $^{2}$/v 2 F + p 2 x ± k/ 2 ± ip$_{y}$ ( | Δ | $^{2}$/v 2 F + p 2 x ± k/ 2) 2 + p 2 y . j = 2 e ∑ p [ ∂ ˜ E$_{1}$ ∂ k ( tanh E$_{0}$ − ˜ E$_{1}$ 2 T − tanh E$_{0}$ + ˜ E$_{1}$ 2 T ) + ∂ ˜ E$_{2}$ ∂ k ( tanh E$_{0}$ − ˜ E$_{2}$ 2 T − tanh E$_{0}$ + ˜ E$_{2}$ 2 T ) v 2 $_{F}$p 2 v 2 $_{F}$[ p k ] p E$_{0}$ ˆ u = 1 2 √ 1 ± p$_{x}$ √ | Δ | $^{2}$/v 2 F + p 2 x ( ± E ± | E ± $_{|}$e − iφ ± / 2 e iφ ± / 2 ) , ˆ v = 1 2 √ 1 ∓ p$_{x}$ √ | Δ | $^{2}$/v 2 F + p 2 x ( ± e − iφ ± / 2 E | E ± $_{|}$e iφ ± / 2 ) , (B16) E 2 ± = ( √ | Δ | 2 + v 2 $_{F}$p 2 x ± v$_{F}$k/ 2) 2 + v 2 $_{F}$p 2 $_{y}$, (B15) and one may check by substitution that the four ortho- normalized solutions of the Bogoliubov–de Gennes equa- tions are given by the spinors