Finite momentum Cooper pairing in three-dimensional topological insulator Josephson junctions

Unconventional superconductivity arising from the interplay between strong spin–orbit coupling and magnetism is an intensive area of research. One form of unconventional superconductivity arises when Cooper pairs subjected to a magnetic exchange coupling acquire a finite momentum. Here, we report on a signature of finite momentum Cooper pairing in the three-dimensional topological insulator Bi2Se3. We apply in-plane and out-of-plane magnetic fields to proximity-coupled Bi2Se3 and find that the in-plane field creates a spatially oscillating superconducting order parameter in the junction as evidenced by the emergence of an anomalous Fraunhofer pattern. We describe how the anomalous Fraunhofer patterns evolve for different device parameters, and we use this to understand the microscopic origin of the oscillating order parameter. The agreement between the experimental data and simulations shows that the finite momentum pairing originates from the coexistence of the Zeeman effect and Aharonov–Bohm flux.

I n the conventional Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity, Cooper pairs form an isotropic condensate with a zero center-of-mass momentum 1 . However, introducing magnetism can change the stability of the BCS superconducting state, thereby destroying superconductivity or, as in unconventional superconductors, altering the pairing symmetry 2 . The potential for unconventional superconductivity at the confluence of magnetism and superconductivity has made it an area of great theoretical and experimental interest. One such example of an unconventional superconducting state is Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superconductivity, which was proposed as a way for maintaining superconductivity even beyond the critical Zeeman field 3,4 . Despite the intensive search for an FFLO superconductor in various types of materials such as heavy fermion compound CeCoIn 5 5,6 and BEDTTTFbased organic superconductors [7][8][9] , FFLO superconductivity still remains a controversial subject 2,10-12 .
To better hunt for unconventional superconductivity, there have been proposals for utilizing materials with strong spin-orbit interaction coupled to a conventional s-wave superconductor. This is predicted to stabilize an FFLO superconducting state: the spin-orbit coupling lifts the degeneracy in the Fermi surfaces of the material and introduces an anisotropy to the surfaces that makes it more amenable to a finite momentum phase [13][14][15] . In particular, time-reversal invariant topological insulators (TIs), whose surface states are massless Dirac fermions, are proposed to be an attractive candidate for unconventional superconductivity that carries finite momentum pairing 16 . To the best of our knowledge, experimental signatures of finite momentum Cooper pairs in TIs have mainly been sought after in the electron-doped regime of the two-dimensional (2D) TI HgTe quantum wells 17 , but the surface states of three-dimensional (3D) TIs also provide unique advantages to engineering finite momentum pairing. The Dirac cones on the surfaces are non-degenerate and have spinmomentum locking. As a consequence, the Fermi surface of the Dirac cone shifts uni-directionally under the application of an in-plane magnetic field to the surface, which can lead to an FFLO state 13,16 . Even though transport measurements in normal 3D TIs are often complicated by the presence of bulk carriers, there is experimental consensus that the metallic surface state dominates transport in a proximity-coupled TI even when the bulk is not depleted [18][19][20][21][22] .
To this end, we study the experimental signatures of Cooper pairs in a superconductor (S)-3D TI-S Josephson junction subjected to in-plane and out-of-plane magnetic fields. We probe the phase of Cooper pairs by generating Fraunhofer patterns with an out-of-plane field, and we find that adding an in-plane field distorts the Fraunhofer patterns by (1) transferring the intensity of superconductivity from the central Fraunhofer peak at B z = 0 out to finite magnetic field values and (2) introducing asymmetries between positive and negative values of B z in the Fraunhofer patterns. We show that the intensity transfer is suggestive of a spatially oscillating superconducting order parameter phase, which we call a finite momentum shift; we propose two potential origins for this finite momentum shift; and we demonstrate that asymmetries in the transport signal come from the sample geometry. Simulations show a close match between experimental data and finite momentum Cooper pair theory.

Results
Experimental set-up. Our devices consist of Bi 2 Se 3 flakes that are mechanically exfoliated from crystals with a bulk carrier density n~5 × 10 17 cm −3 (see Supplementary Note 1). Angle-resolved photoemission (ARPES) measurements have been made in a previous work on crystals similar to the ones used in our experiments, and they show that the Fermi energy is close to the bulk gap 23 . The flakes are contacted with two superconducting electrodes, forming a Josephson junction device, and the measured devices vary in flake thicknesses and junction dimensions ( Table 1). A representative atomic force microscope (AFM) image is shown in Fig. 1a for device 1, which has flake thickness t~9 nm, average junction width W~860 nm, and electrode spacing of 140 nm. Because we utilize high in-plane fields to tune the behavior of the junction, we choose NbTi/NbTiN (T c~1 2.5 K and H c2,in-plane > 9 T) as the superconducting material.
Out-of-plane magnetic field B z is applied to the superconducting junction to generate a Fraunhofer pattern. In the absence of any in-plane fields, the devices exhibit a central peak with maximum critical current I c and decaying side peaks, as is expected for a conventional Fraunhofer pattern. The nodes of our Fraunhofer pattern are at B z ¼ nΦ 0 Wd ; where magnetic flux quantum Φ 0 ¼ h 2e , d is the effective electrode spacing that takes into account flux focusing (see Supplementary Note 2), and n is an integer 24 . The Fraunhofer pattern of device 1, shown in Fig. 1b, is representative of our devices.
When an in-plane field along the current direction B y is introduced to the devices, the conventional Fraunhofer pattern is modulated to an anomalous Fraunhofer pattern: the Fraunhofer pattern is shifted along the B z direction and the critical current of the side lobes increases as the critical current of the central lobe disappears. To measure the evolution of the Fraunhofer pattern as a function of B y , we apply a small AC excitation (with zero DC current) and measure the differential resistance dV/dI as a function of B z and B y , where lower resistance corresponds to higher critical current, similar to in ref. 17 . The evolution of the Fraunhofer pattern for device 1 is shown in Fig. 2a. As B y is applied, the Fraunhofer pattern is shifted along B z . This is evident as an overall tilt of the 2D differential resistance map in Fig. 2a. Although a tilt could be caused by misalignment of the sample with respect to the B y -B z plane, this type of misalignment would cause a similar shift in the Fraunhofer pattern when a field is applied in any in-plane direction. Because we do not see a corresponding shift when an in-plane field perpendicular to current direction B x is applied, we can exclude sample misalignment as a cause for the shift of the Fraunhofer pattern.
Besides the shift to the Fraunhofer pattern, we also observe additional side branch features in the evolution of the Fraunhofer pattern: at B y = 0, the intensity of superconductivity is maximum at the central lobe, but as B y is increased, the intensity is transferred outwards to higher values of B z . The emergence of this anomalous side branch becomes more evident if the tilt in the 2D differential resistance map is removed by rotating the graph until the lobe minima are vertical, as shown in Fig. 2b. The approximate locations of the minima of the lifted side lobes are marked as a guide to the eye, and a slope for the side branch can be approximated, as illustrated by the dashed black line in Fig. 2b (see Supplementary Note 3). This unique Fraunhofer evolution is evidence of finite momentum pairing, as discussed below. Modeling Josephson junction with finite momentum pairing.
To determine the origin of the evolution of the Fraunhofer pattern, we begin by considering the mesoscopic effects of magnetic fields B z and B y on the superconducting order parameter phase. For the following discussion we look primarily at the surface contribution because we find that the bulk order parameter decays much more rapidly than the surface contributions in the junction, which means the supercurrent is predominantly carried by the surface states (see Supplementary Notes 4 and 5). The phase modulation manifests itself as a spatially oscillating current distribution in thex-direction (I x ð Þ ¼ i c sinðΔϕ À 2eB z xd c Þ) 24 . By summing up all the oscillating components of the current, the conventional Fraunhofer diffraction pattern arises, which can be derived from the following equation: While the Fraunhofer pattern is generated by B z , the additional high in-plane magnetic field, B y , generates a Zeeman effect within the surface bands and adds magnetic flux inside the TI flake.
When the in-plane Zeeman effect is present, the low-energy Hamiltonian of the TI surface can be written as where v f is the Fermi velocity of the Dirac cone, g is the g-factor, and μ is the Bohr magneton 25 . By examining the Hamiltonian, we find that the location of the Dirac node is shifted from the Γ-point along thex-direction by gμB y hv f , resulting in a shift of the corresponding Fermi surface. Figure 3a shows a schematic of the shifted Fermi surface. As a result of this shift, when the electrons on the TI Fermi surface form spin singlet Cooper pairs, the Cooper pairs gain a finite center of mass momentum 2gμB y hv f . As a consequence, the superconducting order parameter at the end of each junction has a phase modulation in thex-direction. The order parameter is given as is the phase modulation due to a finite momentum shift. The finite momentum of the Cooper pair under these conditions is similar in nature to FFLO states 2,3 . Unlike for the Fermi surface, the bulk energy eigenvalue is degenerate. The Zeeman effect will therefore only lift the spin degeneracy rather than shifting the bulk bands, so the bulk does not acquire a phase from the Zeeman effect (see the Methods section).
Besides the Zeeman effect contribution to the order parameter, there is also a contribution from the finite flux that is inserted along the in-plane, orŷ, direction of the flake. This magnetic flux from B y results in a phase modulation given by an Aharonov-Bohm phase 2π 26 . For a bulk electron, the phase will be determined by its trajectory inside the flake, and we find that the bulk electrons will acquire a negligible net phase b Fraunhofer evolution for device 1 that has been rotated so that the lobe minima are vertical, making it easier to compare across samples. There is a side branch feature that develops as B y is applied to the junction, which can be quantified as a line with slope m (dashed line). Black dots mark approximate locations of minimum resistance at different side lobes as a guide to the eye modulation (see the Methods section and Supplementary Note 6). For a surface electron, on the other hand, the phase modulation is given as . As a result, we see that there are two finite momentum shift contributions to the surface order parameter of the 3D TI: the Zeeman effect and the Aharonov-Bohm effect, which we call the Zeeman modulation effect (ZME) and flux modulation effect (FME), respectively.
By summing up the three relevant contributions to the phasethe out-of-plane magnetic field B z , the ZME, and the FME-we get the total phase difference between the two junctions: Here, ϕ 1 (x 1 ) and ϕ 2 (x 2 ) are the phases of the order parameters of superconductor 1 and 2 at the coordinates x 1 and x 2 , respectively, along the width of the junction. Based on the phase difference between the two leads, we can model the total transport current along theŷ-direction in the Josephson junction using quasi-classical methods 17,27 . This analysis is equivalent to summing up all possible quasi-classical trajectories of electron transport, so the total transport current is where W 1(2) is the width of the superconducting lead 1(2), Δϕ is the overall phase difference between the superconductors, and d is the distance between the two superconductors. Using Eq. (4), we can calculate the critical current I c ðB y ; B z Þ ¼ max ϕ Iðϕ; B y ; B z Þ as a function of B z and B y to derive the evolution of the Fraunhofer pattern. Figure 3b shows simulations of the Fraunhofer pattern for various values of B y calculated using Eq. (3) and illustrates how the intensity of superconductivity is transferred from the central peak out to the side peaks as B y is increased. This transfer of intensity is proportional to the momentum shift of the Cooper pair. In terms of the differential resistance as a function of B z and B y , the transfer of the superconducting intensity can be seen as the formation of two side branches (formed by evolving side peaks) in the differential resistance map with slope m. These side branches are visible in the simulation for a symmetric Josephson junction in Fig. 3c and in the data for device 1 (Fig. 2b). The agreement between the simulation and the experimentally observed pattern indicates that the formation of the side branches is a result of the transfer of superconducting intensity due to the additional phase modulation generated by B y . This feature is known to be a key signature of finite momentum pairing and distinguishes the system from typical BCS superconductivity 16,28,29 .
Additionally, the slope m of the side branches reflects the relative contributions of finite momentum pairing due to ZME and FME as a function of B y . In the simulations, the slope is defined by a line between the origin and the nth side lobe as n becomes large (green dashed line in Fig. 3c). By calculating the We find that the intensity of superconductivity is transferred outwards to higher values of B z as B y is increased. c Evolution of the Fraunhofer pattern for a symmetric Josephson junction due to finite momentum pairing. The differential resistance is calculated and normalized to 1. The slope of the side branch is indicated by the dashed green line integral in Eq. (4), the slope of the side branches is estimated as In Eq. (5), the first and the second terms in the denominator are the contributions of the ZME and the FME to the slope, respectively. Because of the inverse relation, larger slopes reflect smaller ZME and FME contributions. Looking at the FME and ZME contributions separately, we can see that the FME contribution is proportional to the thickness of the flake t since the flux through B y will increase as the thickness increases. The ZME, on the other hand, is proportional to the intrinsic material parameters of the TI, g v f , rather than on an external parameter, such as the thickness.
We examine the slope dependence on TI thickness across multiple samples, where devices 2-4 are shown in Fig. 4a-c, respectively, with approximate minima marked and the 2D differential resistance maps rotated so that the lobe minima are vertical for ease of comparison. The slope m for each device is extracted from the minima (see Supplementary Note 3) and is illustrated for device 2 by the dashed line in Fig. 4a. To compare the experimental data with theory, we also calculate m for each device based on t and d using Eq. (5). Figure 4d shows the dependence of m (normalized by an effective d that takes into account flux focusing effects) on thickness, where m is extracted from the data (black) and calculated using theoretical predictions for finite momentum pairing due to ZME alone (red), FME alone (blue), and ZME and FME together (purple).
It is clear that the dominant contribution to the finite momentum shift comes from the FME. In the thickest sample, in particular, the FME closely predicts the slope in the experiment since more flux is enclosed in the flake and the orbital effect therefore has a more significant contribution. However, there are some deviations between the data and the FME curves. First, some of the normalized slopes do not follow the monotonically decreasing trend predicted by the FME theory. The deviations are caused by variations in the sample: in real samples, there are variations in chemical potentials, which would result in the g v f ratio changing from sample to sample, and fabrication imperfections, which would result in distortions in the transport signal. Nevertheless, the overall decreasing trend in the experimental slopes mirrors the decreasing trend predicted by the FME model.
We also find that the FME is not enough to predict the observed slope from the experiment on its own even when the error bars of slope calculation are taken into account. In fact, the ZME contribution needs to be considered for the theory to more closely match the data, which means that the finite momentum pairing due to the shifted Dirac cones is non-negligible. We calculate a curve that takes into account an additional ZME value (purple) that demonstrates that by introducing the ZME contribution to the theory, the theory curve is shifted closer to the experimental data. Therefore, our data are generally better explained by the coexistence of both FME and the unconventional ZME, which is more closely related to the FFLO state.
We also looked at the potential phase contribution from inserting a flux through the superconducting leads. In the supplementary section of a previous work 17 , it was found that a flux through the electrode seemed to contribute to the phase Effect of Josephson junction asymmetries on the evolution of the Fraunhofer pattern. In addition to simulating the Fraunhofer evolution as B y is applied to an ideal junction, we also consider the effect of asymmetries in the junction geometry and on the evolution of the Fraunhofer pattern. Due to the fact that typical sample fabrication can result in imperfect device configurations and flux focusing effects, it is important to understand what happens to the Fraunhofer pattern as the devices deviate from the ideal junction. For example, as reported in ref. 30 , asymmetric features between positive and negative values of B z often appear in Fraunhofer patterns and can be attributed to a combination of device-dependent factors such as disorder and the microscopic structure of the device. We consider some sources of asymmetries in the device configuration to model the effect of these asymmetries on the transport signal. One form of geometric asymmetry that arises in a Josephson junction is the asymmetry in the width of the two superconducting leads W 1 and W 2 . To model this effect, we introduce the width asymmetry factor α, which is the ratio of the two superconducting lead widths and satisfies W 1 = αW 2 in Eq. (4). Figure 5a, b shows how the Fraunhofer evolution changes as we increase the asymmetry between W 1 and W 2 by increasing α. Because finite α breaks the symmetry of Iðϕ; B y ; B z Þ upon reversing the sign of B z in Eq. (4), we find that for increasing α, the amplitude of the left and right side branches becomes more asymmetric. Another form of asymmetry, as quantified by flux asymmetry factor β, comes from the flux focusing effect 30 . Due to the screening of the magnetic field inside the superconductor, magnetic field B y may bend and cause contributions to B z . Since we apply large B y compared to B z , a very small bending of B y can cause a large tilt in the Fraunhofer pattern (see Supplementary  Information). We model this effect by replacing B z with (B z −βB y ) in Eq. (3), which generates the tilt seen in Fig. 5c, d. After understanding possible origins for anomalous features in the evolution of the Fraunhofer pattern, we can now compare our experimental data with simulations that take into account asymmetry factors. The results are shown in Fig. 6 (see Supplementary Note 8 for the detailed numerical methods). The width asymmetry factor α is extracted from scanning electron microscope images of the devices (summarized in Table 1). Due to the difficulty of quantifying the magnitude of the flux asymmetry, we add in an artificial β factor to simulations that generates a tilt that best matches with the data for better comparison. As discussed earlier, the overall structure and inplane field dependence of the Fraunhofer pattern is determined by the ZME, the FME, and the geometry of the junction. However, we find that by also incorporating the sample width asymmetries into the finite momentum pairing model and adding an artificial tilt in the data to take into account flux focusing effects, the theoretical prediction and the experiment agree very well.

Discussion
In conclusion, we observe an anomalous Fraunhofer pattern that is indicative of the presence of finite momentum pairing in 3D TI Josephson junctions that are subjected to in-plane magnetic fields. We identify the two microscopic origins of the finite momentum pairing to be the ZME and the FME. By comparing the slope of the side branches in the anomalous Fraunhofer pattern with the theoretical predictions, we conclude that the measured slope can only be explained by the coexistence of both the ZME and the FME. In particular, the ZME-the contribution associated with the FFLO phase-becomes significant when there is less phase accumulation due to enclosed flux, which occurs for thinner samples. We believe that this is a promising start for the hunt for Simulation of effect of junction asymmetry and field inhomogeneity. The width and field asymmetry dependence of the Fraunhofer pattern. a, b When a width asymmetry factor α is added to the model, we find that the signal becomes asymmetric between positive and negative B z . Here, the amplitude of the left side lobes increases relative to the amplitude of the right side lobes. We used the values α = 0.3, 0.6 respectively. c, d When the asymmetry factor β is introduced to the model, the Fraunhofer patterns are shifted along B z , which can be seen as a tilt introduced to the 2D differential resistance map. We used β = 0.01, 0.02 respectively unconventional superconductivity in a proximity-coupled 3D TI in the presence of in-plane fields, but further work can be done to mitigate the effect of the FME. For example, besides finding thinner flakes, the ZME can be enhanced by increasing g/v f , which can be done by tuning the TI flake closer to the Dirac point. Furthermore, as demonstrated by our work and others, a continued understanding of the effect of non-ideal junctions is conducive to identifying anomalous asymmetric signatures, like the Fraunhofer asymmetry across B z , in transport signals.

Methods
Experimental techniques. 3D TI flakes were mechanically exfoliated from bulk Bi 2 Se 3 crystals onto Si/SiO 2 substrates, and thicknesses were identified using atomic force microscopy. After identifying suitable flakes, superconducting electrodes were defined using standard ebeam lithography techniques. Forty to 65 nm of NbTi/ NbTiN was then sputtered onto the device using an rf source following a brief ion mill to clean off the surface. Devices were then wire-bonded and measured in a dry dilution unit that reaches a base temperature of T = 25 mK and has a three-axis vector magnet.
Numerical methods. After calculating the Josephson current using Eq. (4), we numerically generate Fraunhofer patterns as a function of B y for each device. We find that a non-zero B y transfers the intensity of superconductivity from B z = 0 out to higher values of B z . To illustrate this feature, we present various slices of the Fraunhofer pattern as a function of B y in Supplementary Figure 6. Once we derive a set of Fraunhofer patterns as a function of B y , we can draw the surface contour plot of the critical current as a function of both B y and B z and map the critical current map to a differential resistance map with an effective thermal noise, as detailed in Supplementary Note 7.
Calculation of bulk and surface contributions to the ZME. In this section, we determine how the Zeeman effect affects surface and bulk states by calculating the normal propagator of the surface and the bulk of the TI. We first solve for the surface propagator. The Hamiltonian of the TI surface with a finite in-plane Zeeman term is given as and E f is the Fermi energy. The energy of the above Hamiltonian is given as After deriving the eigenstate and the energies and taking just the positive eigenstate (since we assume the chemical potential is near the conduction band), we can use these results to solve for the surface propagator, g surf (E,r): 0 dθ e ik f rcos θ k Àθ r ð Þ 1 e Àiθ k e iθ k 1 is the Fermi wave vector. The above integration is the Bessel function of the first kind, J 0 : In the canonical quasi-classical approximation, we assume that k f r ) 1, so the Bessel functions can be asymptotically approximated as This is the final expression for the surface propagator, where the e iax term represents oscillations in the surface propagator due to the Zeeman effect.
We now solve for the bulk propagator. To do so, we follow a similar procedure to the one used to solve for the surface propagator. Using Supplementary Eq. (2), the bulk Hamiltonian of the TI with the Zeeman effect can be written as The experimental data (a-d) agrees well with the normalized differential resistance of the simulations (e-h) that are based on a finite momentum shift model and take into account width asymmetries delineated in Table 1. An additional tilt is added to the simulations for better comparison with the data. (A possible origin for the tilt is a flux asymmetry, which is difficult to quantify.) In all the devices, side branches with slope m appear as an in-plane field B y is added. The appearance of these side branches is indicative of interference in the phase modulation due to B z and B y where β = gμB. The bulk energy eigenvalue is given as It is important to note that Eq. (6) differs from the surface energy eigenvalue because, in the bulk, the Zeeman effect does not shift the location of the Fermi surface but rather only lifts the degeneracy. As a consequence, we will see that the ZME is absent from the bulk band contribution.
We again solve for the normal propagator The eigenstate of the above Hamiltonian can be explicitly calculated by assuming k j j ( m. In this case, we can focus on just the conduction band to simplify the expression: Because the bulk band is initially degenerate and does not shift in one direction in momentum space, the bulk propagator does not have oscillations from the ZME.
Calculation of surface and bulk quasiparticles of the Aharonov-Bohm phase. In addition to the Zeeman effect, the electrons gain additional phases due to the FME. In this section, we calculate the phase accumulated by electrons in a magnetic field as they acquire an Aharonov-Bohm phase. The vector potential of the inplane magnetic field can be written as A ¼ ðB y ðz À t 2 Þ; 0; 0Þ. Under this gauge, the surface electrons, which travel along the outmost trajectory of the TI, gain the conventional Aharonov-Bohm phase. More specifically, the electrons at the top surface, z = t, gain momenta opposite from the bottom surface, z = 0, which is consistent with a current circulating around the surface loop. When the surface electron classically travels from x 1 to x 2 , the phase is given as ϕ x 2 ð Þ À ϕ x 1 ð Þ ¼ 2π total flux ϕ 0 travel distance circumference ¼ 2π tWB y ϕ 0 Unlike the surface electrons, the set of classical trajectories of the bulk electrons can be characterized by the number of reflections off the bottom surface (Supplementary Figure 4). Similarly the Green function can be decomposed as follows: F iω; r ð Þ¼F 0 iω; r ð ÞþF 1 iω; r ð Þþ ; where F n represent the normal propagators that consist of the trajectories having n reflections off the bottom surface. When the travel distance is longer than the magnetic length l B ¼ hc eB À Á 1=2 , the Green function gains the Aharonov-Bohm phase, 2π Á dl, where the integral is taken over the chord connecting r 1 and r 2 .
There are two distinct Aharonov-Bohm phases that come from the bulk quasiparticle trajectories: (1) a phase due to a trajectory with no reflections so that the quasiparticle is transferred directly from one superconducting contact into the other and (2) a phase due to a trajectory that reflects off of the bottom surface more than once before being transferred into the other contact. If there are no reflections, the bulk quasiparticle will follow the same trajectory as the surface quasiparticles. In this case, the bulk electron should gain the same phase as the surface electron: If the bulk quasiparticle is reflected from the bottom surface at least one time, it will follow a continuous path from z = 0 to z = t and the bottom half will have the opposite modulation from the top half of the flake. In this case, the net phase cancels each other out: As a result, bulk electrons with reflection-less trajectories gain the same FME effect as surface electrons while the bulk electrons with reflected trajectories do not experience the FME and only contribute to the conventional Fraunhofer pattern.
Data availability. The data that support the findings of this study are available from the corresponding author on reasonable request.