Josephson effects in twisted cuprate bilayers

Twisted bilayers of high-$T_c$ cuprate superconductors have been argued to form topological phases with spontaneously broken time reversal symmetry $\mathcal{T}$ for certain twist angles. A key outstanding challenge is to identify unambiguous signatures of these topological phases in experiments. With this goal in mind we theoretically investigate a suite of Josephson phenomena between twisted cuprate layers. At intermediate twist angles we find an unusual nonmonotonic temperature dependence of the critical current, which we attribute to the unconventional sign structure of the $d$-wave order parameter. The onset of the $\mathcal{T}$-broken phase near the ${45}^{\circ }$ twist is marked by a crossover from the conventional $2\pi$-periodic Josephson relation $J\left(\phi \right)\simeq J_{c}sin\phi$ to a $\pi$-periodic function as the single-pair tunneling becomes dominated by a second-order cotunneling process. Despite this fundamental change, the critical current remains a smooth function of the twist angle $\theta$ and temperature $T$ implying that a measurement of $J_c$ alone will not be a litmus test for the $\mathcal{T}$-broken phase. We show that clear signatures of the $\mathcal{T}$-broken phase appear when $J_c$ is measured in the presence of an applied magnetic field or radiofrequency drive: The resulting Fraunhofer oscillations acquire additional nodes and fractional Shapiro steps become visible. We discuss these results in light of recent experiments on twisted bilayers of the high-$T_c$ cuprate superconductor ${\mathrm{Bi}}_2{\mathrm{Sr}}_2\mathrm{Ca}{\mathrm{Cu}}_2{\mathrm{O}}_{8+\delta }$.

Figure 1

(a) Schematic of the Josephson effect measurement on a twisted cuprate bilayer. Magnetic field $\mathit{B}$ parallel to the interface and electromagnetic radiation with frequency $\omega$ can be used to probe the current-phase relationship. (b) Typical phase diagram of the system calculated from the microscopic model discussed in Sec. 2. It is assumed that $d+is$ order is nucleated near the magic angle ${\theta }_M$ with the $s$ component driven by a weak on-site attractive interaction. Model parameters are ${\epsilon }_c=60\mathrm{meV}$, $g=10.5\mathrm{meV}$, $N_F\mathcal{V}=0.12$, and $N_F{\mathcal{V}}_s=0.043$.

Figure 2

[(a),(b)] Typical current-phase relations obtained from the continuum model. Data points represented as black dots are obtained by evaluating Eq. (20) at fixed twist angle $\theta$ for a range of temperatures between zero and just above $T_c$. The color-coded solid lines represent the best fit to the relation Eq. (6) amended by a third-order term $-J_{c3}sin\left(3\phi \right)$ to achieve a close fit. Model parameters are ${\epsilon }_c=60\mathrm{meV}$, $g=7\mathrm{meV}$, $N_F\mathcal{V}=0.12$. Panel (c) shows the temperature dependence of the best-fit coefficients $J_{c\alpha }(\alpha =1,2,3)$ for $\theta ={44}^{\circ }$.

Figure 3

(a) Critical current $J_c$ and (b) superfluid stiffness ${\rho }_s$ as a function of temperature for various twist angles. The traces for twist angles above ${40}^{\circ }$ are shown in $1^{\circ }$ increment. Red lines at $\theta ={10}^{\circ }$ correspond to the magic angle ${\theta }_M$. Blue lines show results for $\theta =42-{45}^{\circ }$ where time reversal is spontaneously broken at low $T$. Panel (c) shows the results for the superfluid stiffness calculated in the presence of the secondary $s$-wave order parameter that arises at low $T$ close to ${\theta }_M$ as indicated in the phase diagram Fig. 1. In these plots the same model parameters are used as in Fig. 2.

Figure 4

Behavior of the critical current for various accessible commensurate twist angles in the lattice model. For a detailed discussion of the model, see Ref. [15]. The parameters used are $\mu =-1.3t$, $g_0=10\mathrm{meV}$ in (a) and $\mu =-1.4t$, $g_0=30\mathrm{meV}$ in (b).

Figure 5

Critical current $I_c(\mathrm{\Phi },T)$ vs applied magnetic flux $\mathrm{\Phi }$ in the twist junction for various temperatures $T$, with twist angle $\theta ={42}^{\circ }$ (top) and $\theta ={44}^{\circ }$ (bottom). The curves are computed from a Lawrence-Doniach type model of the bilayer junction, taking into account the first and second harmonic terms of the current in Eqs. (6) and (31) with coefficients $J_{c1,c2}$ fitted from Fig. 2, closely following the discussion of Ref. [23].

Figure 6

(a) Current-Voltage curves in the ${45}^{\circ }$ twisted configuration in the presence of an external rf drive with amplitude $J_{\mathrm{rf}}=0.7$ at three different frequencies. The voltage is scaled in units of $\hslash \omega /2e$ to highlight the Shapiro steps and the lowest four fractional steps are indicated with the horizontal dashed lines. [(b),(c)] Voltage as a function of drive current for two twists that are in the topological and trivial phases respectively, with $\omega =0.6$. The color bar indicates step amplitudes. Since the numerical values of $J_{c1,2}$ at different twists vary widely, to aid comparison, $J\left(\phi \right)$ is normalized such that the critical current in the absence of rf drive ($J_{\mathrm{rf}}=0$) is unity (or equivalently the zeroth-order step has width 2). Clearly, fractional steps appear only for twist angles close to ${45}^{\circ }$, as in panel (b). At high rf currents, one enters a regime where the step amplitudes oscillate. The temperature is fixed at $T=0.5T_c$ and $R=0.7$ in all simulations.

Figure 7

(a) Current-phase relation computed for two flakes, each with $M=12$ monolayers, for twist angles near ${45}^{\circ }$. Various curves are for temperatures $T$ between zero and $T_c$ in increments $T_c/20$. (b) Critical current $J_c\left(T\right)$ and (c) phase stiffness ${\rho }_s\left(T\right)$ extracted from the calculated $J\left(\phi \right)$. We obtained results for $M$ up to 24 but we find that the curves become essentially independent of $M$ when $M\ge 4$.

Figure 8

Evolution of the spectral $E_{\mathrm{gap}}$ with the flake thickness $M$. Note that $E_{\mathrm{gap}}$ values below $\sim 0.1\mathrm{meV}$ are comparable to the finite-size gap in our system ($500\times 500$ grid of $\mathit{k}$ points) and are therefore consistent with gapless behavior. Inset shows scaling of the maximum gap $E_{\mathrm{gap}}\left({45}^{\circ }\right)$ with $M$, suggesting exponential decay, and the critical twist angle ${\theta }_c$ beyond which the $\mathcal{T}$-broken phase sets in.

Figure 9

Temperature dependence of the superconducting gap amplitude (blue) calculated by numerically iterating the gap equation (18), and the in-plane superfluid stiffness $\gamma \left(T\right)$ (red) evaluated as the coefficient of the $q^2$ term in Eq. (C8).