Fractional Josephson effect in nonuniformly strained graphene

TL;DR

This paper predicts a fractional Josephson effect in strained graphene with induced superconductivity, where pseudo-Landau level inversion leads to gapless modes and 4π-periodic tunneling, revealing novel topological phenomena.

Contribution

It introduces a mechanism for realizing the fractional Josephson effect in strained graphene via pseudo-Landau level inversion and gapless mode formation.

Findings

Prediction of 4π-periodic fractional Josephson effect.

Identification of gapless modes along the junction.

Demonstration of pseudo-Landau level inversion effects.

Abstract

Nonuniform strain distributions in a graphene lattice can give rise to uniform pseudomagnetic fields and associated pseudo-Landau levels without breaking time-reversal symmetry. We demonstrate that by inducing superconductivity in a nonuniformly strained graphene sheet, the lowest pseudo-Landau levels split by a pairing gap can be inverted by changing the sign of the pairing potential. As a consequence of this inversion, we predict that a Josephson $\pi$ junction deposited on top of a strained graphene sheet exhibits one-dimensional gapless modes propagating along the junction. These gapless modes mediate single electron tunneling across the junction, giving rise to the $4\pi$-periodic fractional Josephson effect.

Equations

$$H_0$$

$$H_s$$

$$H_0=-t\sum _{\mathbf{k}\sigma }\sum _{m=1}^3{e}^{i\mathbf{k}\cdot {\mathbf{a}}_m}{c}_{B\sigma }^{†}(\mathbf{k})c_{A\sigma }(\mathbf{k})+H.c.$$

$${\mathcal{H}}_0=v_F({\tau }_z{\sigma }_x{p}_x+{\sigma }_y{p}_y),$$

$$H_s^{K_{\pm }}=\sum _{\sigma }\int d^{2}r\delta t(y)e^{i{K}_{\pm }\cdot {\mathbf{a}}_3}{c}_{B\sigma }^{K_{\pm }}(\mathbf{r}{)}^{†}{c}_{A\sigma }^{K_{\pm }}(\mathbf{r})+H.c.,$$

$${\mathcal{H}}_s=\delta t(y){\sigma }_x.$$

$$\mathcal{H}={\mathcal{H}}_0+{\mathcal{H}}_s=v_F\left[{\tau }_z{\sigma }_x\left(-i{\partial }_x+\delta t(y){\tau }_z/v_F\right)-i{\sigma }_y{\partial }_y\right].$$

$$\mathcal{H}=v_F[{\tau }_z{\sigma }_x(p_x+{\tau }_{z}eBy)-i{\sigma }_y{\partial }_y].$$

$$\mathcal{H}=\frac{\sqrt{2}{v}_F}{l_B}\left(\begin{array}{cccc}0& a^{†}& 0& 0\\ a& 0& 0& 0\\ 0& 0& 0& -{\overline{a}}^{†}\\ 0& 0& -\overline{a}& 0\end{array}\right).$$

$$H_{\Delta }=$$

$${\mathcal{H}}_{\text{BdG}}=\frac{v_F}{\sqrt{2}{l}_B}\left(\begin{array}{cccc}0& a^{†}& \tilde{\Delta }& 0\\ a& 0& 0& \tilde{\Delta }\\ {\tilde{\Delta }}^{\ast }& 0& 0& -a^{†}\\ 0& {\tilde{\Delta }}^{\ast }& -a& 0\end{array}\right),$$

$${\psi }_0^{\text{BdG}}(\theta ,E_0=+|\Delta |)$$

$${\psi }_0^{\text{BdG}}(\theta ,E_0=-|\Delta |)$$

$${\psi }_0^{\text{BdG}}(\theta =0,E_0=+|\Delta |)$$

$${\psi }_0^{\text{BdG}}(\theta =0,E_0=-|\Delta |)$$

$$H=$$

$$\displaystyle-\sum_{j=2}^{N}\bigl{(}t-\delta t(y)\bigr{)}c^{\dagger}_{B,j-1,\sigma}(k_{x})c_{A,j\sigma}(k_{x})+\textrm{H.c.}\Biggr{]},$$

$$H_{\Delta }=$$

$$+(\uparrow \leftrightarrow \downarrow ).$$

Full text

I INTRODUCTION

Coupling quantum Hall states with superconductivity van Ostaay et al. (2011); Amet et al. (2016); Allen et al. (2016, 2015); Lee et al. (2016) has been intensively studied recently as it provides a platform for exotic excitations such as Majorana fermions Lindner et al. (2012); San-Jose et al. (2015) and parafermions Clarke et al. (2013); Mong et al. (2014); Alicea and Stern (2015); Alicea and Fendley (2016). In contrast to the supercurrent in conventional Josephson junctions that is mediated by Cooper pairs, these elusive modes give rise to single electron or fractional charge tunneling, which results in a fractional Josephson effect Kitaev (2001); Alicea et al. (2011) with flux periodicity larger than $2\pi$.

Magnetic field strengths normally required for the quantum Hall effect typically suppress proximity-induced superconductivity, which leads to challenges in attempting to couple these two phenomena. Instead of coupling superconductivity and quantum Hall states, an alternative is to couple superconductivity with quantum spin Hall states Maciejko et al. (2011), i.e., 2D time-reversal invariant topological insulators. The topological superconductor that emerges from this coupling gives rise to Majorana modes and the fractional Josephson effect Fu and Kane (2009). Strong spin-orbit coupling, an essential ingredient of the quantum spin Hall effect, is crucial for the occurrence of the fractional Josephson effect in this proposed realization.

In this paper we propose an alternative way of realizing the fractional Josephson effect that requires neither spin-orbit coupling nor a magnetic field. Our proposal relies on the recently demonstrated ability to engineer large uniform pseudomagnetic fields in graphene by applying nonuniform distributions of strain Guinea et al. (2010a); Levy et al. (2010); Ghaemi et al. (2012); Guinea et al. (2010b); Masir et al. (2013); Low and Guinea (2010); Roy and Juričić (2014). Just like a physical magnetic field in a 2D electron gas gives rise to the Hall effect and, upon quantization, Landau levels, in strained graphene the pseudomagnetic field leads to the valley Hall effect and pseudo-Landau levels Vaezi et al. (2013). Experimentally, strain-induced pseudomagnetic fields in excess of 300 Tesla have been reported in graphene nanobubbles Levy et al. (2010); Gomes et al. (2012). Strain-induced pseudo-Landau levels have also been observed in graphene grown by chemical vapor deposition Yeh et al. (2011).

Coupling superconductivity to strained graphene modifies the pseudo-Landau levels Uchoa and Barlas (2013) in a manner similar to how the Landau levels of a massless Dirac Hamiltonian are affected by the addition of a Dirac mass. In the Dirac-Landau problem a mass domain wall (i.e., an interface at which the Dirac mass changes sign) gives rise to a gapless chiral 1D mode that propagates along the domain wall Callan Jr. and Harvey (1985); Qi et al. (2008) and disperses in the otherwise empty intra-Landau-level bulk gap. Likewise here, we find that an interface at which the pairing term changes sign, i.e., a Josephson $\pi$ junction, is accompanied by 1D gapless modes dispersing within the bulk pairing gap between the two lowest (zeroth) pseudo-Landau levels. These gapless modes in turn lead to a $4\pi$ energy-phase relation and the fractional Josephson effect.

II MODEL OF STRAINED GRAPHENE

The device we consider is a Josephson junction between two $s$-wave superconducting electrodes deposited on top of a graphene sheet with nonuniform (linearly increasing) uniaxial strain (Fig. 1), for instance using the methods described in Ref. Heersche et al., 2007. There are three vectors, ${\mathbf{a}}_{1}=a(\sqrt{3}/2,1/2)$, ${\mathbf{a}}_{2}=a(-\sqrt{3}/2,1/2)$ and ${\mathbf{a}}_{3}=a(0,-1)$ that connect any atom on an $A$ sublattice of the graphene lattice to the nearest neighboring atom on a $B$ sublattice, where $a$ is the closest distance between carbon atoms. In strained graphene, we assume the hopping amplitudes along ${\mathbf{a}}_{1}$ and ${\mathbf{a}}_{2}$ remain unmodified with value $t$, but the hopping amplitude along the ${\mathbf{a}}_{3}$ direction is reduced as $t-\delta t(y)$ Ghaemi et al. (2012); He and He (2013) where $\delta t(y)\propto y$ represents the effect of linearly increasing strain in the $y$ direction. The Hamiltonian with nearest neighbor hopping is composed of two terms, $H=H_{0}+H_{s}$. The first is the usual nearest neighbor hopping term on the honeycomb lattice,

[TABLE]

where ${\mathbf{r}}$ denotes sites of the underlying triangular Bravais lattice and $\sigma=\uparrow,\downarrow$ is spin. The second represents a modification due to strain, and is given by

[TABLE]

Eq. (1) can be expressed in momentum space as

[TABLE]

Expanding the momentum ${\bf{k}}$ near the Dirac points $K_{\pm}\equiv\pm 4\pi/(3\sqrt{3}a)(1,0)$ as ${\bf{k}}=K_{\pm}+{\bf{p}}$ where $|{\mathbf{p}}|$ is small compared to the dimensions of the Brillouin zone allows us to rewrite $H_{0}$ as a Dirac Hamiltonian Castro Neto et al. (2009); Ghaemi et al. (2012); Hou et al. (2007) in the basis of $\Psi_{\sigma}({\bf{p}})=(c^{K_{+}}_{A\sigma}({\bf{p}}),c^{K_{+}}_{B\sigma}({\bf{p}}),c^{K_{-}}_{A\sigma}({\bf{p}}),c^{K_{-}}_{B\sigma}({\bf{p}}))^{T}$ where $c_{A,B,\sigma}^{K_{\pm}}({\mathbf{p}})\equiv c_{A,B,\sigma}(K_{\pm}+{\mathbf{p}})$. More specifically, the effective low-energy Hamiltonian near the Dirac points is expressed as $H_{0}\approx\sum_{\sigma}\int[d^{2}p/(2\pi)^{2}]\Psi_{\sigma}^{\dagger}\mathcal{H}_{0}\Psi_{\sigma}$, where

[TABLE]

describes linearly dispersing Dirac fermions with velocity $v_{F}\equiv(3/2)at$. We use $\sigma_{i}$ and $\tau_{i}$, $i=x,y,z$ to denote the Pauli matrices in sublattice ($A$ and $B$) and valley ($K_{+}$ and $K_{-}$) space, respectively.

We now turn to the term $H_{s}$ in Eq. (2) that is due to strain. Assuming $\delta t(y)$ varies slowly on the scale of the lattice constant $a$, the two valleys remain approximately decoupled and one can write $H_{s}\approx H^{K_{+}}_{s}+H^{K_{-}}_{s}$ with Ghaemi et al. (2012); He and He (2013); Masir et al. (2013)

[TABLE]

where ${\mathbf{r}}=(x,y)$ now denotes the continuum position and $c_{A,B,\sigma}^{K_{\pm}}({\mathbf{r}})$ is the continuum Fourier transform of $c_{A,B,\sigma}^{K_{\pm}}({\mathbf{p}})$. In our choice of basis $K_{\pm}\cdot{\mathbf{a}}_{3}=0$, thus the strain Hamiltonian becomes $H_{s}\approx\sum_{\sigma}\int d^{2}r\,\Psi_{\sigma}^{\dagger}\mathcal{H}_{s}\Psi_{\sigma}$ with

[TABLE]

The full Hamiltonian for nonuniformly strained graphene in the continuum limit is given as

[TABLE]

Comparing Eq. (4) and Eq. (7), the variation in hopping strength along the ${\mathbf{a}}_{3}$ direction is equivalent to the substitution $-i\partial_{x}\rightarrow-i\partial_{x}+\delta t(y)\tau_{z}/v_{F}$. In other words, the strain in graphene generates a pseudo-vector potential Ghaemi et al. (2012); He and He (2013) along the $x$ direction given by ${{\mathbf{A}}_{s}}({\mathbf{r}})=\tau_{z}\delta t(y)\hat{{\mathbf{x}}}/(ev_{F})$. The Pauli matrix $\tau_{z}$ in ${{\mathbf{A}}_{s}}$ signifies that electrons in different valleys experience a pseudomagnetic field of the same magnitude but opposite sign. This comes from the fact that the pseudomagnetic field preserves time-reversal symmetry. One can indeed check that Eq. (7) is invariant under a time-reversal transformation $\mathcal{T}$. Momentum changes sign under time reversal, which implies $\mathcal{T}\tau_{z}\mathcal{T}^{-1}=-\tau_{z}$ because we exchange the valley $K_{+}$ to $-K_{+}=K_{-}$ (modulo a reciprocal lattice vector) under time reversal. The sublattices $A$ and $B$ remain the same under time reversal, giving $\mathcal{T}\sigma_{x}\mathcal{T}^{-1}=\sigma_{x}$ and $\mathcal{T}\sigma_{y}\mathcal{T}^{-1}=-\sigma_{y}$ because $\mathcal{T}i\mathcal{T}^{-1}=-i$. Combining these results gives $\mathcal{T}\mathcal{H}\mathcal{T}^{-1}=\mathcal{H}$.

Because a linearly increasing vector potential corresponds to a uniform magnetic field, we expect pseudo-Landau levels to appear in Eq. (7). Writing the strain-induced variation in hopping strength as $\delta t(y)=ev_{F}By$ where $B$ is the pseudomagnetic field Abanin and Pesin (2012), Eq. (7) becomes

[TABLE]

Translation invariance in the $x$ direction allows us to replace the momentum operator in the $x$ direction by its eigenvalue $p_{x}$, while translation invariance is broken by the strain pattern in the $y$ direction. Equation (8) is almost the same as the Dirac Hamiltonian in a uniform magnetic field in the Landau gauge, except for the presence of $\tau_{z}$ in the vector potential. To demonstrate the presence of pseudo-Landau levels, we define the dimensionless coordinate $\xi\equiv l_{B}p_{x}+y/l_{B}$ for the $K_{+}$ valley and $\overline{\xi}\equiv l_{B}p_{x}-y/l_{B}$ for the $K_{-}$ valley, where $l_{B}\equiv 1/\sqrt{eB}$ is a pseudomagnetic length. This allows us to define the raising operator as $a^{\dagger}\equiv\frac{1}{\sqrt{2}}(\xi-\partial_{\xi})$ for the $K_{+}$ valley and $\overline{a}^{\dagger}\equiv\frac{1}{\sqrt{2}}(\overline{\xi}-\partial_{\overline{\xi}})$ for the $K_{-}$ valley. The corresponding lowering operators are $a=\frac{1}{\sqrt{2}}(\xi+\partial_{\xi})$ and $\overline{a}=\frac{1}{\sqrt{2}}(\overline{\xi}+\partial_{\overline{\xi}})$. These operators satisfy the commutation relations $[a,a^{\dagger}]=1$, $[\overline{a},\overline{a}^{\dagger}]=1$, and allow us to write the first quantized Hamiltonian (8) in a simple form,

[TABLE]

By solving the eigenvalue problem $\mathcal{H}\psi_{n}=E_{n}\psi_{n}$, we obtain the pseudo-Landau levels as $E_{n}=\pm\sqrt{2}v_{F}\sqrt{n}/l_{B}$, where $n=0,1,2,\ldots$ The corresponding eigenstates $\psi_{n}=(\psi^{K_{+}}_{n},\psi^{K_{-}}_{n})^{T}$ are given in terms of the simple harmonic oscillator eigenstates $|n\rangle$ as $\psi^{K_{+}}_{n}=(|n\rangle,\pm|n-1\rangle)^{T}$ and $\psi^{K_{-}}_{n}=(|\overline{n}\rangle,\mp|\overline{n-1}\rangle)^{T}$ for the $K_{+}$ and $K_{-}$ valleys, respectively. Here the simple harmonic oscillator eigenstates $|n\rangle$ are given as $|n\rangle\propto\exp(-\xi^{2}/2)H_{n}(\xi)$ for the $K_{+}$ valley and $|\overline{n}\rangle\propto\exp(-\overline{\xi}^{2}/2)H_{n}(\overline{\xi})$ for the $K_{-}$ valley, where $H_{n}$ is the $n$th Hermite polynomial. The wavefunctions are plane waves in the $x$ direction, with an infinite degeneracy parametrized by the momentum eigenvalue $p_{x}$.

The wavefunction of the zeroth pseudo-Landau level is $\psi^{K_{+}}_{0}=(|0\rangle,0)^{T}$ for the $K_{+}$ valley and $\psi^{K_{-}}_{0}=(|\overline{0}\rangle,0)^{T}$ for the $K_{-}$ valley, thus the wavefunctions for both valleys have support only on the A sublattice. This is in contrast to the case of a real magnetic field where the zeroth Landau level wavefunctions for the $K_{+}$ and $K_{-}$ valleys have support on the A and B sublattice, respectively. Support on the same sublattice in a pseudomagnetic field is one of the key reasons why an on-site pairing potential opens a superconducting gap in the zeroth pseudo-Landau level, as will be seen shortly.

III PROXIMITY-INDUCED SUPERCONDUCTIVITY

To couple the pseudo-Landau levels with superconductivity, we add a proximity-induced pairing term $\Delta$ to the Hamiltonian,

[TABLE]

using $\uparrow,\downarrow$ to denote the physical spin. The full Hamiltonian with the pairing term included is $SU(2)$ spin rotationally invariant and all eigenstates have an exact two-fold degeneracy; in the following we factor out this degeneracy and focus on solving the Hamiltonian in one of the two nondegenerate subspaces. This Hamiltonian can be expressed in the Nambu/Bogoliubov-de Gennes (BdG) basis $\Psi^{\text{BdG}}(p_{x})=[c^{K_{+}}_{A\uparrow}(p_{x}),c^{K_{+}}_{B\uparrow}(p_{x}),c^{K_{-}\dagger}_{A\downarrow}(-p_{x}),c^{K_{-}\dagger}_{B\downarrow}(-p_{x})]^{T}$ as

[TABLE]

where $\tilde{\Delta}=\frac{\sqrt{2}l_{B}}{v_{F}}\Delta$ is a dimensionless measure of the pairing gap. Diagonalizing Eq. (11), we obtain the energy spectrum as $E_{n}=\pm\frac{v_{F}}{\sqrt{2}l_{B}}\tilde{E}_{n}$ where $\tilde{E}_{n}\equiv\sqrt{|\tilde{\Delta}|^{2}+n}$ is dimensionless. The corresponding (unnormalized) eigenstate is given as $\psi^{\text{BdG}}_{n}=(\tilde{\Delta}|n\rangle,\tilde{\Delta}|n-1\rangle,(\tilde{E}_{n}-\sqrt{n})|n\rangle,(\tilde{E}_{n}-\sqrt{n})|n-1\rangle)^{T}$. The energy spectrum $E_{n}$ is the same as the Landau level spectrum of a massive Dirac Hamiltonian with mass $\tilde{\Delta}$. In a 2D topological insulator, a change of sign of the mass in the Dirac Hamiltonian is accompanied by band inversion. We thus expect the zeroth pseudo-Landau level to become inverted as the sign of the pairing term in Eq. (11) is reversed.

We verify this intuition by explicit calculation. The $n=0$ pseudo-Landau level splits into dispersionless BdG bands $E_{0}=\pm|\Delta|$ in the presence of the pairing gap $|\Delta|$. Writing the pairing potential as $\Delta=|\Delta|e^{i\theta}$, the corresponding wavefunction $\psi(\theta,E_{0})$ of the zeroth pseudo-Landau level with energy $E_{0}$ and phase $\theta$ is

[TABLE]

The zeroth pseudo-Landau level wavefunction with energy $E_{0}=|\Delta|$ and phase $\theta=0$ is the same as the wavefunction with energy $E_{0}=-|\Delta|$ and phase $\theta=\pi$. In other words,

[TABLE]

Equation (13) is the main result of this work, and means that the zeroth pseudo-Landau level undergoes band inversion as the sign of the pairing potential is reversed. In experiment, reversing the sign of the pairing potential can be achieved by building a Josephson junction with a phase difference of $\pi$ across the junction. For this reason, we consider a Josephson junction that is built on top of strained graphene as shown in Fig. 1.

We now analyze the energy-momentum relation of the zeroth pseudo-Landau level in a Josephson $\pi$ junction. The zeroth pseudo-Landau level wavefunction is $|0\rangle\propto\exp(-\xi^{2}/2)$, which peaks at $\xi=0$. By setting $\xi=0$, we locate the guiding center of the wavefunction at $y_{0}=-l^{2}_{B}p_{x}$. Assuming the pairing potential changes sign at $y=0$, the pseudo-Landau levels at energy $\pm|\Delta|$ cross as the momentum $p_{x}$ goes from negative to positive (Fig. 2).

Figure 2 implies that when the zeroth pseudo-Landau level is inverted across the $\pi$ junction, gapless modes appear within the pairing gap. We confirm this result achieved with a continuum model with that obtained from a lattice model. The energy-momentum relation of strained graphene with a $\pi$ junction in a lattice model is shown in Fig. 3(a). We consider a graphene lattice with translation invariance in the $x$ direction but with a finite width with zigzag edges in the $y$ direction as in Fig. 1. We Fourier transform Eq. (1) and Eq. (2) along the $x$ direction to wavevector space $k_{x}$, which gives the tight-binding Hamiltonian of a strained graphene strip,

[TABLE]

where $j=1,\ldots,N$ is the lattice index along the $y$ direction and $y\equiv(3/2)a[j-(N+1)/2]$ in $\delta t(y)$. The proximity-induced pairing term becomes

[TABLE]

To model a Josephson $\pi$ junction, the pairing potential is expressed as $\Delta=\Delta_{0}e^{i\theta}$ where $\theta=0$ for $y<0$ and $\theta=\pi$ when $y>0$. Diagonalizing the total Hamiltonian $H+H_{\Delta}$ numerically, we find a crossing of the zeroth pseudo-Landau levels at the both the $K_{+}$ and $K_{-}$ Dirac points [Fig. 3(a)], consistent with our analysis of the continuum model (Eq. (13) and Fig. 2). This crossing results in gapless 1D propagating modes localized near the $\pi$ junction [Fig. 3(b)-(d)] and appearing within the pairing gap. For energies near the Dirac point $|\tilde{E}|\ll|\tilde{\Delta}|$, the gapless modes disperse linearly with momentum $\tilde{E}\propto\pm p_{x}$. Since the guiding center of the zeroth pseudo-Landau wavefunctions is $-l_{B}^{2}p_{x}$, in this low-energy regime the spatial separation between the two counter-propagating modes is proportional to $\tilde{E}$.

IV DISCUSSION

The 1D gapless modes mediate single electron tunneling, which leads to a fractional Josephson effect with a $4\pi$ energy-phase relation (Fig. 4). The pseudo-Landau level inversion in the $\pi$ junction implies a zero-energy mode at the $K_{+}$ and $K_{-}$ points. Denoting by $\delta\theta$ the superconducting phase difference across the junction, the zero-energy mode $E(\delta\theta=\pi)=0$ at the $K_{+}$ or $K_{-}$ Dirac point implies $4\pi$ periodicity in the energy-phase relation of the junction. To demonstrate this effect, we fix the momentum $k_{x}=\pm 4\pi/(3\sqrt{3}a)$ at one of the Dirac points $K_{\pm}$ and diagonalize the total Hamiltonian $H+H_{\Delta}$ while varying the phase difference $\delta\theta$. This gives the spectrum shown in Fig. 4. Most of the pseudo-Landau levels do not disperse with $\delta\theta$, except the zeroth level which has a $4\pi$-periodic energy-phase relation. The supercurrent $I_{s}\propto\partial E/\partial(\delta\theta)$ is proportional to the derivative of the total energy $E$ with respect to the phase difference $\delta\theta$, which means only the zeroth pseudo-Landau level contributes, giving rise to a $4\pi$-periodic current-phase relation. In the absence of a $\pi$ junction, the non-dispersing pseudo-Landau levels correspond to a flat band superconductor Kauppila et al. (2016).

The $4\pi$-periodic Josephson effect is traditionally associated with unpaired Majorana fermions in 1D topological superconductors Kitaev (2001). While in the topological superconductor two Majorana modes at opposite ends of a wire form a single zero-energy electronic state which mediates single-electron tunneling, in our case this is achieved by two zero-energy electronic states in each valley which are degenerate by $SU(2)$ spin rotation symmetry. The zero-energy modes at the $K_{+}$ and $K_{-}$ Dirac points [see Fig. 3(a)] can be expressed in terms of eight Majorana operators. Unlike in topological superconductors, however, such Majorana modes are not topologically protected.

In conclusion, we have shown that the zeroth pseudo-Landau level induced by straining graphene in zero magnetic field can be inverted by a Josephson $\pi$ junction, which reverses the sign of the pairing potential. As a consequence of pseudo-Landau level inversion, a pair of gapless counter-propagating 1D modes appears near the center of the junction. These gapless modes in turn lead to a $4\pi$ energy-phase relation, giving rise to the possibility of observing the fractional Josephson effect in strained graphene.

Acknowledgements.

We thank Chien-Hung Lin for illuminating discussions. This research was supported by NSERC grants #RGPIN-203396 (F.M.) and #RGPIN-2014-4608 (J.M.), the Canada Research Chair (Program CRC), the Canadian Institute for Advanced Research (CIFAR), and the University of Alberta.