Phase Slips in Superconducting Weak Links

TL;DR

This paper investigates the nonlinear dynamics of phase slips in superconducting weak links, revealing that their onset occurs through an infinite period bifurcation, with analytical and numerical results aligning well.

Contribution

It introduces a weakly nonlinear approximation for phase slips in weak links and identifies the bifurcation mechanism behind their onset, which was previously not well understood.

Findings

Phase slips occur via an infinite period bifurcation.

The voltage-current dependence shows characteristic signatures of this bifurcation.

Analytical results agree with numerical simulations.

Abstract

Superconducting vortices and phase slips are primary mechanisms of dissipation in superconducting, superfluid, and cold atom systems. While the dynamics of vortices is fairly well described, phase slips occurring in quasi-one dimensional superconducting wires still elude understanding. The main reason is that phase slips are strongly non-linear time-dependent phenomena that cannot be cast in terms of small perturbations of the superconducting state. Here we study phase slips occurring in superconducting weak links. Thanks to partial suppression of superconductivity in weak links, we employ a weakly nonlinear approximation for dynamic phase slips. This approximation is not valid for homogeneous superconducting wires and slabs. Using the numerical solution of the time-dependent Ginzburg-Landau equation and bifurcation analysis of stationary solutions, we show that the onset of phase slips occurs via an infinite period bifurcation, which is manifested in a specific voltage-current dependence. Our analytical results are in good agreement with simulations.

Equations

$$\Gamma \left({\partial }_t+i\frac{2e}{\hslash }\mu \right)\Psi =a_0\nu (x)\Psi -b|\Psi {|}^2+\Psi +\frac{{\hslash }^2}{4m}{\partial }_x^2\Psi ,$$

$$u({\partial }_t+i\mu )\Psi ={\partial }_x^2\Psi +[\nu (x)-|\Psi {|}^2]\Psi$$

$$j_0=\Im ({\Psi }^{\ast }{\partial }_x\Psi )-{\partial }_x\mu .$$

$$\Psi$$

$$\mu$$

$$u[{\partial }_t+ix({\partial }_{t}K-A)+i\tilde{\mu }]\tilde{\Psi }=({\partial }_x+iK{)}^2\tilde{\Psi }+[\nu (x)-|\tilde{\Psi }{|}^2]\tilde{\Psi }.$$

$$j_0=\Im ({\tilde{\Psi }}^{\ast }{\partial }_x\tilde{\Psi })+|\tilde{\Psi }{|}^{2}K+{\partial }_{t}K-{\partial }_x\tilde{\mu }.$$

$${\partial }_{t}K+⟨|\tilde{\Psi }{|}^2⟩K=j_0-\Im ⟨{\tilde{\Psi }}^{\ast }{\partial }_x\tilde{\Psi }⟩\equiv j_n.$$

$$u({\partial }_t+i\mu )\Psi$$

$${\mu }_x$$

$$j_n$$

$$\nu (x)\equiv \{\begin{array}{ll}1,& x\notin I\\ -C,& x\in I\end{array}.$$

$$0$$

$$j_0$$

$$0={\partial }_x^{2}F+[\nu (x)-j_0^2{F}^{-4}-F^2]F.$$

$$j_c=\frac{1}{2\sqrt{C}}{e}^{-2r\sqrt{C}},$$

$$j_c=\frac{1}{2}{e}^{-2r}.$$

$$j_c\approx \frac{1}{2}{e}^{-2\underset{k}{\max }{r}_k}.$$

$$u{\partial }_t\eta$$

$$i(2{\partial }_x\varphi {\partial }_x\eta +{\partial }_x^2\varphi \eta -uF\mu )-F^2{\eta }^{\ast }$$

$$0$$

$$\begin{array}{rl}0& ={\partial }_x^{2}U+\left(\nu -({\partial }_x\varphi {)}^2-3F^2-\lambda \right)U-\\ & (2{\partial }_x\varphi {\partial }_{x}V+V{\partial }_x^2\varphi )\end{array}$$

$$\begin{array}{rl}0& ={\partial }_x^{2}V+\left(\nu -({\partial }_x\varphi {)}^2-F^2-\lambda \right)V+\\ & (2{\partial }_x\varphi {\partial }_{x}U+U{\partial }_x^2\varphi )-uF\mu \end{array}$$

$${\partial }_x\mu$$

$$⟨|\Psi {|}^2⟩$$

$$\approx \frac{1-j_0^2}{L}(L-r)+\frac{1}{L}{\int }_0^r{F}^2+2F(\eta +{\eta }^{\ast })dx$$

$$\approx 1-j_0^2+O\left(\frac{r}{L}\right).$$

$$K=\frac{ϵ}{1-j_0^2}+O\left(\frac{r}{L}\right).$$

$$\begin{array}{rl}u{A}_{\tau }& =B+c_1{A}^2\\ u{B}_{\tau }& =c_{2}AB+c_3{A}^3,\end{array}$$

$$x=\frac{A}{L_A},Y=\frac{B}{L_B},t^{\prime }=\frac{t}{u{L}_t}.$$

$$L_A=\frac{1}{c_2{L}_t},L_B=\frac{1}{c_2{L}_t^2},$$

$$\begin{array}{rl}\dot{X}& =Y+a{X}^2\\ \dot{Y}& =XY+b{X}^3,\end{array}$$

$$H(X,Y)=\frac{1}{2}{Y}^2-\frac{1}{2}{X}^{2}Y-\frac{1}{4}b{X}^4.$$

$$Y=\frac{U}{2a+1}-a{X}^2,$$

$$\dot{X}$$

$$\dot{U}$$

Full text

I Introduction

The motion of Abrikosov vortices is recognized as the main cause of dissipation in type-II superconductors Blatter et al. (1994). Conversely, in thin nanowires, the motion of vortices is impeded and phase-slip events are responsible for the dissipation. Phase slips, changing the phase difference of the superconducting order parameter by $2\pi$, may be caused by different physical mechanisms. Thermally activated phase slips at high temperatures and small applied currents are well understood Tinkham (1996). At very low temperatures, phase slips can be caused by quantum fluctuations (aptly called quantum phase slips) Mooij and Nazarov (2006); Lau et al. (2001); Glatz and Nattermann (2002). Phase slips are not unique to superconductors, they also occur in superfluid systems Anderson (1966); Langer and Fisher (1967); Schwarz (1990), and more recently, dissipation due to phase slips were studied in cold atom systems McKay et al. (2008); Scherpelz et al. (2014, 2015). In particular, phase slips can be triggered in a superfluid cold atom system by a rotating weak link Wright et al. (2013).

Even without thermal and quantum fluctuations, the phase slip phenomena and dissipative (or resistive) states can be induced by an applied current Skocpol et al. (1974); Meyer (1973). Magnetic field penetrates type-II superconductors in the form of Abrikosov vortices. If an external current is applied, the Lorentz force induces motion of the vortices. This motion is the main cause of dissipation in 2D and 3D superconductors. However, in quasi-one dimensional nanowires with the coherence length $\xi(T)$ and the penetration depth $\lambda(T)$ large compared to the wire diameter, vortex motion is suppressed. In this situation the transition to the normal state was made through successive voltage jumps which are attributed to the appearance of phase slip centersSkocpol et al. (1974); Meyer (1973). A study of this phenomenon was given first by Kramer and Baratoff who found that slightly below the depairing current, there is a dissipative state which consists of localized phase slips occurring in the superconducting filament Kramer and Baratoff (1977). In a narrow range of currents close to the depairing current, the material is superconducting except in narrow regions where phase slip centers (PSCs) occur. The period of these PSCs diverge as the external current approaches the lower bound in this narrow region. It was also shown that random thermal fluctuations allow for phase slips Little (1967), but these did not persist indefinitely. Further numerical study of the one-dimensional time-dependent Ginzburg-Landau equation revealed periodic phase slips existing in a narrow range of currents close to the depairing current Kramer and Rangel (1984); Rangel and Kramer (1989). Follow-up numerical studies of narrow two-dimensional superconducting strips discovered a transition from a phase-slip-line to vortex pairs Weber and Kramer (1991). Periodic lattices of the phase slip centers were studied in the context of vortex penetration in thin superconducting films near the third critical magnetic field Aranson and Vinokur (1998). Using a saddle-point approximation for the Ginzburg-Landau energy in narrow superconducting strips, the dependence of voltage drop vs temperature and bias current (neglecting thermal fluctuations) was studied in [Ovchinnikov and Varlamov, 2015].

The situation is different, however, for spatially inhomogeneous systems, such as superconductors with macroscopic defects or weak links Langer and Ambegaokar (1967). Perhaps the most famous examples are Dayem bridges and Josephson junctions Josephson (1962); Anderson and Dayem (1964). The mechanism for dissipation in these cases is the quantum tunneling of Cooper pairs between the two superconductors, which is caused by a phase difference between the weakly-linked superconductors. When the current is below some threshold $j_{c}$, the phase difference is fixed in time and a stationary superconducting state persists. Above this threshold, the solution exhibits oscillations, which lead to a finite voltage. In a review paper by Ivlev and Kopnin, inhomogeneities were analyzed, but in regards to the stability of the normal state Ivlev and Kopnin (1984). Thus, their analysis involved currents much closer to the GL critical threshold $j_{GL}=2/\sqrt{27}$. A lower bound $j_{1}$ at which the normal state was globally unstable (i.e. arbitrary small perturbation lead to instability of the normal state), and above which there was a critical-sized perturbation which separated the normal and superconducting states was estimated. Also, an upper critical current $j_{2}$ such that the normal state was absolutely stable for an external current $j_{0}>j_{2}$ was found. An inhomogeneity much smaller than the coherence length, $\xi(T)$, was used and was approximated by a $\delta$ function, simplifying the algebra. Here we consider a more realistic situation for the type-II high-temperature superconductors: an inclusion on the scale of $\xi(T)$. The transition we are interested in analyzing, occurs between the non-uniform superconducting state and the oscillatory state with phase slips. Therefore, the steady state and linearization in this paper are much more complex then in analyzing the normal state. The authors of [Van Dover et al., 1981] have shown experimental results of weak-links with non-hysteric behavior.

The phase slip state of homogenous systems have recently been analyzed in much greater detail Baranov et al. (2011). Using bifurcation analysis, Baranov et. al. extract the normal form of a saddle-node bifurcation when the current is near the critical current. They then correctly determine the characteristic scaling law and show its agreement with numerical simulations. The period diverges in an infinite-period saddle-node bifurcation as $j_{0}\to j_{c}$. These authors further expanded upon their analysis by showing the important role that the material parameter $u$ plays in the type of bifurcation that can occur Baranov et al. (2013) ($u$ is related to the electric field penetration depth). They observed that for finite lengths and values of $u$ above some critical threshold $u_{c2}$, numerical simulations showed hysteresis in the I-V curve. However, our work focuses on analytical methods for the inhomogeneous system, which as stated previously makes the steady state and linearization to much more difficult to handle. We show that a simplified system can be obtained through weakly nonlinear analysis and that this system contains the normal form obtained in [Baranov et al., 2011] as the size of the weak link shrinks to zero. We also demonstrate that in addition to the infinite period bifurcation for small $u$, a hysteresis exists in our system for large $u$ values, similar to that in Ref. [Baranov et al., 2011]. However, in contrast to previous studies, our reduced two-dimensional nonlinear system exhibits evolution of periodic orbits and a transition between superconducting and normal states that are not properly captured by the one-dimensional model in Ref [Baranov et al., 2011].

A work by Michotte et. al. in [Michotte et al., 2004] have found that the condition for PSCs to occur is based on the competition between two relaxation times: the relaxation time for the magnitude of the order parameter $t_{|\Psi|}$ and the relaxation time for the phase of the order parameter $t_{\phi}$. They observed that phase slips are possible only when $t_{\phi}<t_{|\Psi|}$. A linearized Eilenberger equation in the dirty limit was studied, resembling a generalized TDGL equation with additional parameters related to inelastic electron-phonon collisions, which was first derived in [Kramer and Watts-Tobin, 1978]. They derived an approximate critical current via this equation and their results implied that there was a finite maximal oscillation period for the order parameter. In contrast, for weak links all oscillation periods diverged. The generalized GL equation used contained an additional parameter $\gamma$ characterizing relative superconducting phase relaxation time (for us, $\gamma=0$). For large $\gamma$ values hysteresis was observed in the I-V curve. On a qualitative level, the effect of increasing parameter $\gamma$ is similar to an increase in parameter $u$ Baranov et al. (2011). Correspondingly, we observed hysteresis when $u\gg 1$. The authors of [Berdiyorov et al., 2012] have done numerical analysis of a periodic array of weak links using the generalized TDGL equation. They showed I-V curves for different magnetic fields, however no analysis of the divergence of the period of vortices was presented.

We focus on a 1D superconductor, separated by a normal or weakly superconducting inhomogeneity. The complete system is modeled by a spatially dependent critical temperature $T_{c}(x)$. The weak link is created by a lower transition temperature inside an interval $I=[-r,r]$, which leads to a suppression of the order parameter. Here $r$ is the inclusion radius. Below some critical current, this system relaxes to a stationary superconducting state, but above it, the superconductor exhibits a finite voltage with oscillatory behavior. Thermal fluctuations are initially not considered in this model and therefore does not cause a finite voltage in the superconducting state. The Josephson junction analysis is not applicable here. Indeed, since there is no dielectric contact between the two superconductor pieces, the phase should always be the same, implying zero voltage. We will show via simulations of the time-dependent Ginzburg-Landau equation, that the oscillations in the voltage is caused by phase slips in the center of the inclusion. The system approaches this state via a saddle-node bifurcation of two superconducting states, which occur at the critical current (at a saddle-node bifurcation stable and unstable stationary superconducting states annihilate and a periodic resistive state appears). The suppression of the order parameter in and near the weak link allows us to employ analytical methods in the vicinity of the critical current. We derive a reduced two-dimensional system governing the time evolution of the phase slip solution and describe a sequence of transitions between superconducting and dissipative states.

The paper is organized as follows: section II describes the model, section III deals with the stationary case and estimates the critical current which is obtained from the saddle-node bifurcation condition. Sections IV-VII deal with the time periodic solutions, extracting a time-dependent system via weakly nonlinear analysis and then studying the simplified model to show that it exhibits the same qualitative behavior. In section VIII, we interpret our analytical results, show the correspondence to the parameters of the superconductor and its effects on the phase slip state. Finally, section IX gives closing remarks and ideas for further study.

II Governing equations

The time-dependent Ginzburg-Landau equations (TDGLE) are obtained by minimization of the GL free energyAranson and Kramer (2002). In the absence of a magnetic field this results in

[TABLE]

where $\Gamma,a_{0},b$ are phenomenological parameters that can be found from the microscopic theory Bardeen et al. (1957), $e,m$ are the electron charge and mass, $\mu$ is the scalar potential, and $\nu(x)$ a spatially dependent linear coefficient modeling inhomogeneities in the system. Following Sadovskyy et al.Sadovskyy et al. (2015), we define the $+x$ direction be the direction of the external current and obtain the following dimensionless form:

[TABLE]

Here $\Psi$ is the complex order parameter, satisfying $|\Psi|=1$ in the purely superconducting state, and $|\Psi|=0$ in the normal state. The parameter $u=\Gamma/a_{0}\tau_{\text{GL}}$ with time $\tau_{\text{GL}}=4\pi\sigma\lambda_{0}^{2}/c^{2}$, $\lambda_{0}=\sqrt{\frac{mc^{2}}{8\pi e^{2}\psi_{0}^{2}}}$ is the magnetic penetration depth ($c$ the speed of light) and $\psi_{0}=\sqrt{a_{0}/b}$ is the equilibrium value of the order parameter when spatial variations are neglected, i.e., $\nu(x)=1$. The zero temperature coherence length $\xi_{0}=\sqrt{\frac{\hbar^{2}}{4ma_{0}}}$ is used for the unit of length. For more details see [Sadovskyy et al., 2015].

We apply periodic boundary conditions for $\Psi$. Since $\mu$ is on average an increasing function of $x$, there is necessarily a discontinuity at the boundary. This is resolved by making the following transformations:

[TABLE]

Here, $\tilde{\mu}$ is a periodic function in $x$. Essentially, we are moving the growth of $\mu$ to the phase of $\Psi$. The growth in $K$ now does not affect the magnitude. Indeed, this also allows us to rewind $K$ through $K\to K-(2\pi/\Delta x)\lfloor K\Delta_{x}/2\pi\rfloor$ which will remove any error from rapid phase oscillations Sadovskyy et al. (2015). Inserting this into (2a) gives

[TABLE]

Setting $\partial_{t}K=A$ eliminates the linear term. Now inserting this into (2b), we have

[TABLE]

Averaging this equation over space and noting that $\langle\tilde{\mu}_{x}\rangle=0$ results in an ordinary differential equation (ODE) for $K$

[TABLE]

For clearer notation, we now suppress the tildes, and we arrive at our modified TDGLE

[TABLE]

The integration domain is periodic with the period $L$. For the numerical integration, we generally took $L=20$ and $u=1$, however this was relaxed to see if the qualitative behavior changed. We verified that increasing $L$ does not affect the results, however changing $u$ can have a large effect (see section VIII.3). To make the analysis simpler, we placed the weak link of length $2r$ symmetrically at the origin in the interval $I$. The inclusion’s effect enters through the term $\nu(x)\Psi$ defined by

[TABLE]

Numerical analysis has shown that for $L\gg r$ there exists a critical current $j_{c}$, which is a function of $r$ that separates the dynamics of this system. For $j_{0}<j_{c}$, the system goes to a stationary superconducting state, while for $j_{0}>j_{c}$ the system exhibits a dissipative state represented by periodic phase slips occurring in the center of the inclusion via a stable limit cycle. In the following sections, we explain these results analytically. We first provide an analytical approximation of the critical current. Next, we extract a coupled two-dimensional nonlinear system of ODEs from (2a) which describes qualitatively, the correct behavior for suitable choices of the coefficients of the simplified system.

III The stationary case $j_{0}<j_{c}$

In the superconducting state with an applied current of $j_{0}<j_{c}$, it can be shown that $\mu=0$, (see appendix .1 for details). To proceed, we rewrite (2a) in terms of amplitude and phase of the order parameter, i.e., $\Psi=Fe^{i\phi}$. Inserting this into (2a) and (2b) gives for the stationary equation

[TABLE]

Plugging (6b) into (6a) gives the nonlinear ODE

[TABLE]

III.1 Large $C$ approximation

We now assume a large $C$ approximation, that is, the weak link strongly suppresses superconductivity in the inclusion (i.e. $C\gg j_{0}^{2}F^{-4}$). This allows us to neglect the nonlinear term and obtain a first order approximation of the solution of (6). From this we notice that (6a) has a first integral for both the inclusion domain and the superconducting domain. Asymptotic analysis of the size of these coefficients gives us a condition for $j_{c}$ given by

[TABLE]

for details see appendix .2. Setting $C=1$, we have that

[TABLE]

Comparing this approximation with numerical simulations, we see that the large $C$ approximation with $C=1$ is in good agreement with the numerical solution (see Fig. 1). Thus, we derived that a weak link results in a exponential suppression of the critical current as a function of the inclusion width $2r$ and strength $C$. A similar result was obtained through a different method in Ref. Rangel and Kramer, 1989. However, our method is appealing for the simple generalization to multiple inclusions.

III.2 Multiple inclusions

Let $r_{1},\dots r_{k}$ be the radii of $k$ inclusions in the domain. We have $k+1$ superconducting domains and $k$ normal domains, each with their own first integral constant. The analysis from appendix .2 carries over and we expect the inclusion domain’s first integral constant $E_{I_{k}}$ to be approximately 0 for each $k$. This holds at the center of each respective inclusion, which each give different critical currents. However, when one is no longer satisfied, the system will no longer be satisfied and the global $j_{c}$ is determined by the lowest local $j_{c}$, which appears at the longest inclusion:

[TABLE]

III.3 Linear stability analysis of the stationary state

Consider now a perturbation $\eta$ of the stable state in the form $\Psi=(F+\eta)e^{i\phi}$. Inserting this into (2a)-(2b) and linearizing in $\eta$, we obtain with (6a) and (6b)

[TABLE]

Separating $\eta(x,t)=(U+iV)e^{\lambda t}$ we obtain the following system (here $\lambda$ is the growth rate)

[TABLE]

This system along with (6a)-(6b) represents a 7 dimensional boundary-value eigenvalue problem which must be solved with appropriate boundary conditions. First, we note from (6a) that replacing $x\to-x$ leaves the differential equation unchanged. This with the reflection symmetry implies that $F$ is an even function in $x$. This symmetry implies from (2b) that $\partial_{x}\phi$ and $\partial_{x}\mu$ are even in $x$. Thus $x\to-x$ changes $\Psi\to\Psi^{*}$. The action of this must be retained in the linearization implying that $\eta(-x)$ and $\eta^{*}(x)$ are both solutions. Hence $U$ is even and $V$ is odd in $x$. Furthermore, by symmetry it suffices to solve the equations only on the half interval $(0,L/2)$ with the obtained natural boundary conditions from symmetry and the remaining conditions to be found by matching-shooting algorithm. To solve this we used a technique developed in Tsimring and Aranson (1997); Aranson and Vinokur (1998). In order to do so, we used a numerical matching-shooting solver for ODEs by beginning with a small domain (typically $L\sim 3$). We extracted the appropriate shooting boundary conditions and approximation for $\lambda$ and used these as guesses for a larger system size. Iterating this process, we continued to $L$ sufficiently large until the boundary conditions and $\lambda$ were not changing significantly. The results are plotted in Fig. 2. We note here that $j_{c}\approx 0.0637$ obtained by the solver is only 6% away from the value obtained through direct numerical solution of the Ginzburg-Landau model. The step size used in the dynamic simulations were much larger ($\Delta x=0.05$ compared to shooting solver with $\Delta x=0.001$) and each had an associated numerical error. Therefore, $j_{c}\approx 0.0637$ is more accurate. We checked if the error is independent of the solvers by analyzing the dynamic simulations $j_{c}$ as a function of $\Delta x$ in appendix .3. We found that as $\Delta x\to 0$, we approached a similar value to that found from shooting. Thus, from Fig. 2 one sees that at the critical current, when stable ($\lambda<0$) and unstable ($\lambda>0$) solutions merge and annihilate, the corresponding linear system becomes degenerate. At the critical point it possesses two zero eigenvalues $\lambda_{1,2}=0$. This degeneracy is taken into account through weakly nonlinear analysis.

IV Analysis of time-periodic solutions for $j_{0}>j_{c}$

When the current is above the critical threshold, the above analysis breaks down. Numerical simulations indicate that the superconductor exhibits oscillations in the order parameter, where phase slips are now present (i.e. $|\Psi(0,t)|=0$ for some $t$). In figure 3 we have estimated the period of oscillation $T$ as a function of $j_{0}-j_{c}\ll 1$. Numerical simulations indicate that the period $T\sim O(|j_{0}-j_{c}|^{-1/2})$, which is indicative of an infinite-period bifurcation (IPB) at the point $j_{0}=j_{c}$. In general for a bifurcation parameter $R$ (e.g. current $j$) the period of oscillations $T\sim O(|R-R_{c}|^{-1/2})$ for $|(R/R_{c})-1|\ll 1$ for an IPBStrogatz (2014). We can see from figure 3 that an IPB is occurring at the critical value. In section VIII.3, we show that for $u\gg 1$, we also observe hysteresis, behavior which is typical of a homoclinic bifurcation, a different mechanism through which a limit cycle can be destroyedStrogatz (2014).

Figure 4 shows time-voltage curves for $j_{0}>j_{c}$. One clearly sees the period diverging as we approach the critical value. To calculate the current-period relationship, we ramped the current from an initial amount (typically $j_{\text{init}}<j_{c}$). If the system was stationary for a certain number of iterations, we increased the current. Once the system started oscillating, we calculated peaks in voltage, while skipping the first few to account for system equilibration. Then we averaged over the remaining peaks to obtain the period. We then used linear extrapolation to find the new current. For example, at the $n^{\text{th}}$ step, we have the current $j_{n}$ and corresponding period $T_{n}$. Let $m_{n}=\Delta T_{n}/\Delta j_{n}$, then suppose we want to find the current corresponding to a new period $T_{n+1}=(1+\alpha)T_{n}$, with $\alpha>0$. This is given by $j_{n+1}=j_{n}+\frac{\alpha T_{n}}{m_{n}}$. Figure 5 shows a similar period divergence of the oscillations of $\Psi$ and the simplified model (see section VII).

V Weakly nonlinear analysis

We now extract a coupled ODE system, which exhibits two dynamical possibilities. In the case $j_{0}<j_{c}$, we show that the stationary (fixed) solution is stable, while in the opposite case, a stable limit cycle exists. It is of course possible that a bistability region can exist, which would lead to hysteric effects. Such effects have been observed in homogeneous superconductors Weber and Kramer (1991); Baranov et al. (2011, 2013). For large $u$, we have also observed hysteric I-V curves and we show that our extracted system contains both possibilities. The process is standard and is broken into these steps:

• •

Find stationary (basic) state $\Psi_{0}=Fe^{i\phi}$ (it is already shown in Fig. 2)

• •

Perturb solution and solve linearized system.

• •

Extract weakly nonlinear effects from orthogonality condition.

• •

Show that certain conditions allow for a stable limit cycle to exist.

Though standard, the difficulty in this problem is that the basic state and linearization cannot be solved in closed form. Though we can approximate it to a certain degree, its region of validity is dependent on the radius of the inclusion $r$, the current $j_{0}$ and to a smaller extent, the system size $L$. Indeed, it is impractical to obtain it numerically since the solutions are sensitive to these choices. However, our analysis will assume that these are all known a priori and proceed through the framework. The simplified system is then obtained generally, and we show that the system exhibits the appropriate behavior for certain values in parameter space.

We expand Eqs. (4a)-(4c) near the stationary solution and near the critical point $j_{0}=j_{c}+\epsilon$ with $\epsilon\ll 1$. The first order solution will be given by $\Psi_{0}=Fe^{i\phi}$ (since $K=0$, there is no electric potential in the super conducting state), in fact the initial transient would show exponential decay of $K\to e^{-\left\langle|\Psi|^{2}\right\rangle t}$ and so $\mu=0$ as expected. Let $\Psi=(F+\eta)e^{i\phi}$, where $\eta$, and time will now both slowly vary and be controlled by a small parameter $0<\delta\ll 1$, whose size will be related to $\epsilon$. The proper scaling will be determined from the ODE for $K$. Based on numerical simulations, we assume $K=O(\delta^{2})$. We claim that we may regard $K$ as constant in the relevant order of the perturbation method by the following argument. The perturbation $\eta$ at first order is highly localized inside the inclusion and from this we argue that

[TABLE]

For $L\gg r$, we can regard $\langle|\Psi|^{2}\rangle$ as a constant. In a similar way all averaged quantities in the voltage equation can be neglected in the large superconductor domain limit. This analysis shows that the time-dependence of the voltage is slaved to the behavior of the order parameter $\Psi$. Therefore, we set $K$ to a constant by

[TABLE]

From this, we extract the relation $\epsilon=\alpha\delta^{2}$ where $\alpha=\pm 1$. The linearized system at $\epsilon=0$ has a degenerate eigenvalue as was shown previously in Fig. 2. Therefore we expand $\eta(x,\tau)=A(\tau)\eta_{1}(x)+\sqrt{\delta}[B(\tau)\eta_{2}(x)+z_{1}(x,\tau)]+\delta z_{2}(x,\tau)$ where $\mathcal{L}\eta_{1}=0$, $\mathcal{L}\eta_{2}=\eta_{1}$ and $\mathcal{L}$ is the linear operator from (11). Using orthogonality conditions, we arrive at the coupled system

[TABLE]

where the coefficients $c_{k}$ can be found through evaluating the integrals (see appendix .4). We will show in section VI why we chose to not include the constant $K$ at this order. The general behavior is only captured correctly at $\epsilon=0$. When $\epsilon\neq 0$ (i.e. $K\neq 0$) we do not see a saddle-node bifurcation. To correct for this deficiency, higher order terms will be included. However, we can still gain some insight by analyzing this simplified system first.

VI Dynamical System Analysis

We begin with (13) by making a dimensionless system to analyze it more easily. We introduce the dimensionless variables

[TABLE]

Inserting this into the system and defining the characteristic variables

[TABLE]

we arrive at the dimensionless system

[TABLE]

where $a\equiv c_{1}/c_{2}$ and $b\equiv c_{3}/c_{2}^{2}$. The characteristic scale for time is arbitrary and is a consequence of the degeneracy in the system. The culprits are the $X^{2}$ term and $XY$ terms whose combination of characteristic scales simultaneously vanish.

VI.1 Fixed points and stability

There is only one fixed point located at the origin, provided that $a\neq b$. In this case there is a family of non-isolated fixed points along the parabola $Y=-aX^{2}$, however this case is not physical so we omit it. Next, we note the symmetry $t\to-t$ and $X\to-X$ of (14), which implies that the linearized center located at the origin is robust. We wish to see if this system exhibits closed orbits. The system is conservative if $a=-1/2$. In this case, a first integral can be obtained

[TABLE]

This has closed orbits provided that $b<-1/2$. So now that we have established the existence of closed orbits, we seek to gain insight if $a\neq-1/2$. We replace $Y$ via the transformation

[TABLE]

and rescale $X\to\frac{X}{2a+1}$ and obtain

[TABLE]

This leaves us with one independent parameter $\gamma$. We have already analyzed the case where $a=-1/2$ which, if $b<-1/2$ corresponds to $\gamma\to-\infty$ and has a family of closed orbits. If $b>-1/2$ then $\gamma\to\infty$ and we know this does not have closed orbits. Therefore, there must be some critical value of $\gamma$ where this behavior changes. We seek a solution of (15) of the form $X=\tilde{C}t^{-1}$ with $\tilde{C}$ to be determined. Plugging this into the equation gives the condition

[TABLE]

These two solutions form a saddle-type connection only when they are equal which occurs at $\gamma_{c}=-1/8$ or in the original coefficients

[TABLE]

This critical curve separates closed orbit solutions in the $(a,b)$ parameter space. We have shown that the simplest (first order) system obtained, demonstrates a saddle-type infinite period bifurcation, however this creates an infinite family of closed orbits and a unique stable limit cycle is not obtained. The bottleneck is created near the origin (see Fig 6). Additionally, it does not have a saddle-node bifurcation which we expect to occur at $j_{0}=j_{c}$. We note also that introducing $K$ at this order, which adds a nonzero constant term to the second ODE would still only have one fixed point and a constant at this order would destroy the degeneracy (and also any closed orbits) in a degenerate Hopf-type bifurcation when that constant crosses through zero. This should be corrected by including the next higher order cubic terms which will saturate and force the system to select one unique closed orbit.

The bottleneck created near the emergence of the saddle-node bifurcation is apparent in both the physical and simplified system (see figure 5). Note that the time scales need not be the same and careful treatment of the parameters in the simplified system (see section IV) would lead to the relation between the GL time and the time scale of the simplified system.

VII Full Dynamical System

We modify the system to include the next order cubic terms. In principle, we could obtain the next order terms by continuing the perturbation expansion, however, we chose to include the generic next higher order terms $X^{3},X^{2}Y,XY^{2}$, and so on. We then found that the removal of some cubic terms e.g. $XY^{2},Y^{3}$ slightly shifts the transitions boundaries but does not qualitatively change the bifurcation sequence. Therefore, we chose to keep the following system for our analysis:

[TABLE]

where we have introduced the new coefficients $c,w_{1},w_{2}$. We will enforce $w_{1},w_{2}<0$ to ensure the phase flow cannot escape to infinity, which would be a nonphysical state for this system.

VII.1 Analysis

The fixed points cannot be found analytically in general since the equation involves a quintic polynomial. Instead we look to find the two critical curves which correspond to our system. We wish to find a saddle-node bifurcation curve and an infinite-period bifurcation as the current is varied. The saddle-node bifurcation involves the merging and annihilation of the stable and unstable stationary solutions. An infinite-period bifurcation is a saddle-node bifurcation which occurs on the limit cycle in the phase plane Strogatz (2014).

We first find the fixed points of (16). Using (16a), we obtain $Y^{*}=-(X^{*})^{2}(a+w_{1}X^{*})$, which leads to the quintic equation

[TABLE]

A saddle-node bifurcation occurs provided that $f(X^{*})=f^{\prime}(X^{*})=0$. The curve exists only if $X^{*}$ is real which leads to the requirement that

[TABLE]

To motivate our choice of parameters, we write this in terms of $\gamma$

[TABLE]

If we set $w_{2}=2w_{1}$ we can eliminate $a$ from the dependence on $\gamma$. Thus, we have that the saddle-node bifurcation exists only if $\gamma\geq-\frac{2}{15}$.

Writing the quintic now with $a=-1$ allows us to cast the quintic function solely in terms of $w_{1},\gamma\text{ and }c$.

[TABLE]

The saddle node bifurcation then occurs along the curve

[TABLE]

where $X^{*}$ is given by

[TABLE]

The Jacobian of this system is

[TABLE]

A necessary condition for a Hopf bifurcation to occur is for a (un)stable spiral to change stability. This occurs when the trace of the Jacobian $\tau=X^{*}[2a+1+(3w_{1}+w_{2})X^{*}]=0$ and the determinant $\Delta>0$. For our analysis this implies that $X^{*}=0$ or $X^{*}=(5w_{1})^{-1}$. Of course our fixed point $X^{*}$ must also satisfy the quintic equation. Inserting this gives a necessary condition and curve in $(\gamma,c)$ space for a Hopf bifurcation

[TABLE]

The determinant is

[TABLE]

Thus, the first Hopf bifurcation curve exists only when $\gamma<-2/15$. The second Hopf bifurcation is more complicated since $\Delta=0$ and so nonlinear terms are important. The existence of that curve was found numerically.

VII.2 Phase Diagram

In general, this system has many different ways in which a limit cycle is destroyed. Numerical experiments indicate that this can occur via a Hopf, cycle bifurcation, infinite period or homoclinic bifurcation. Slightly changing the parameters can change which bifurcation we obtain. From the preceding section, we motivated the choices $w_{1}=-0.05,w_{2}=-0.1,a=-1$ to keep our parameter space $(\gamma,c)$. This leads to a generalized phase diagram of section VI.

The Hopf and saddle-node bifurcation curves of figure 7 were obtained analytically. The IPB curve $c_{\text{IP}}=c_{\text{IP}}(\gamma_{\text{IP}})$ was found numerically and for comparison is compared to the observed physical limit cycle in figures 3 and 5. Additionally, it was found numerically that the HB in region III, did not exhibit the birth of a stable limit cycle. Possible trajectories of the superconductor through this phase diagram is shown with purple lines.

A more generic phase diagram with $w_{2}\neq 2w_{1}$ is given in figure 8. Here, both an IPB and homoclinic bifurcation can destroy the limit cycle. The existence of the homoclinic bifurcation changes the morphology of the phase diagram to now include a bistability region in which the limit cycle (phase slips) and fixed point (superconducting state) coexist. This is particularly encouraging since we also found hysteresis for $u\gg 1$ (see section VIII.3). Possible trajectories of the superconductor through this phase diagram is shown with purple lines.

VIII Discussion

VIII.1 Sensitivity to temperature

To test the sensitivity of these phase slips to small thermal noise we modified (4) to include a small random noise term uniformly distributed between $[-T_{f},T_{f}]$ at each point in space. Numerical simulations indicate that the system is stable to small fluctuations. The qualitative change is the existence of finite voltage in the superconducting state, however the critical current at which phase slips begin is unchanged.

VIII.2 Effect of parameter $u$

The parameter $u$ characterizes the penetration of the electric field. In homogeneous wires, it has been found that hysteresis of the phase slip state exists for finite domains with $u\gg 1$ Baranov et al. (2013). We analyzed $u=0.01,1,10,100$ with $L=20$ and $r=1$ (see figure 9). Another important quantity not yet discussed is that of the retrapping current $j_{r}$. The authors of [Baranov et al., 2013] discuss the effect of $u$, numerically simulating the GL equation and finding a curve separating the hysteresis region of the I-V curve through some length dependent critical curve $u_{c2}(L)$. For our simulations of weak links, $u$ small (for $r=1$, $u<30$ is small enough), $j_{r}=j_{c}$. However for $u\gg 1$, $j_{r}<j_{c}$, this leads to hysteresis in the I-V curve (see figure 10).

VIII.3 Physical quantities in simplified system

The phase diagram is in $(\gamma,c)$ space. We can relate the important physical quantities $u,r,j_{0},L$ to $\gamma,c$ by using appendix .4. The coefficient $c$ is strongly affected by the parameter $u$ and the current $j_{0}$. Consider $j_{0}<j_{c}$ and $u\to 0$, then we know that there is no voltage (i.e. $K=0$), and $\alpha=-1$. This implies that $c\sim-u\zeta^{2}$ for some $\zeta(r,j_{0},L)$ for small $u$. At a significantly large enough $u$ we expect our initial trajectory to begin from a region in figure 8 where hysteresis is possible. Increasing the current $j_{0}>j_{c}$ switches $\alpha=1$ and $K\neq 0$, as $j_{0}$ continues to increase, $F$ decreases and we expect $c$ to change sign as we continue to increase it, which explains our motivation for the direction of trajectories. Increasing $r$ lowers $j_{c}$ and so we expect the trajectories to spend more time in the phase slip state, which leads us to expect that $c$ decreases. A similar argument, leads us to assume the same holds for $\gamma$ (see figure 7 and 8). The effects on $\gamma$ are more complicated for the current and probably non-monotonous in a general case. From physical arguments we know that the trajectories must begin in the superconducting state and move into the phase slip state via either an IPB or homoclinic bifurcation. Comparing this to the phase diagrams, we see that as $j_{0}$ increases, $\gamma$ must decrease. We also attempt to justify this from the terms in appendix .4. We consider the scaling from section VI, which implied $b=c_{3}/c_{2}^{2}$. We noted that $F$ is decreasing as $j_{0}$ increases (where $F^{\prime}$ is relatively unchanged). Again, employing appendix .4, we see that $c_{3}$ is decreasing with the current since the positive terms involve $F$ and the negative terms involve $F^{\prime}$. Finally we use the fact that $b=c_{3}/c_{2}^{2}$ to deduce that $b$ must be decreasing and since $\gamma=(b-a)/(2a+1)^{2}$, we see that $\gamma$ is also decreasing with the current.

IX Conclusion

We have considered a weak-link superconductor qualitatively similar to other weak-link systems, but fundamentally different in mechanism. We demonstrated the existence of a superconducting state and a PSC periodic state separated by a critical current $j_{c}$. This current was calculated asymptotically and agrees very well with numerical simulations. We then extracted a coupled ODE system from the TDGL equations using weakly nonlinear theory and showed under certain choices of parameters, an infinite period bifurcation and homoclinic bifurcations can occur. This demonstrates that the dynamics of phase-slip behavior in weak links described by the TDGL equations can be correctly captured by a simpler system of two coupled ordinary differential equations.

Further research is to extend this analysis to two dimensions. We anticipate additional transitions from phase slips occurring instantly inside the weak link to a more complicated dynamic regime involving phase slips and nucleation of vortex pairs, similar to that in [Weber and Kramer, 1991]. Another interesting generalization is to include disorder in the transverse direction inside the weak link. Possibly, some of the vortices will be pinned in the weak link. It may. in turn, lead to further suppression of the critical current.

The work was supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy for computations. The Office of Science, Advanced Scientific Computing Research and Basic Energy Science, Division of Materials Science and Engineering for analysis.

X Appendix

.1 No voltage in the superconducting state

We begin by multiplying (2a) by $\Psi^{*}$ and we differentiate (2b) with respect to $x$. This gives

[TABLE]

Taking the imaginary part of (17) and substituting this result into (18), we obtain

[TABLE]

Far from the inclusion, all the applied current is supercurrent and so if $L\gg r$, we expect $j_{0}=\Im{(\Psi^{*}\partial_{x}^{2}\Psi|)_{x=\pm L}}$, which implies that $\partial_{x}\mu(\pm L)=0$. Multiplying (19) by $\mu$ and integrating over the domain gives

[TABLE]

Noting the boundary conditions for $\mu$ and the fact that $\partial_{t}\Psi=0$ (stationary state), we see that $\mu\equiv 0$.

.2 Critical current calculation

We separate (7) by region (superconducting vs. normal metal) and then take the first integral to obtain the equations

[TABLE]

Now, far from the inclusion (near the boundary of the superconductor), $F\to F_{\infty}$ a constant. Assuming the relevant approximation that $j_{0}\ll 1$, we see that $F_{\infty}^{2}\approx 1-j_{0}^{2}$. Inserting this into (20), implies that $E_{S}\approx\frac{1}{2}+j_{0}^{2}$. We now use the large $C$ approximation that $C\gg j_{0}^{2}F^{-4}$. Proceeding, we obtain

[TABLE]

where we have introduced the radius $r$ of the inclusion. Solving the outer region at first order is given by

[TABLE]

The two constants are determined by the continuity conditions at the boundary of the inclusion. By symmetry, we may analyze just one side of the boundary, then our conditions are

[TABLE]

Solving for $K_{1}$ and $K_{2}$, we obtain

[TABLE]

Note the identity $E_{S}-E_{I}=(1+C)F^{2}(r)\geq 0$. This implies that

[TABLE]

Motivated by this, we assume that $E_{I}$ is a small parameter. At first order then $E_{I}=0$ and looking at $x=0$ we see that

[TABLE]

where the derivative has vanished by symmetry. Since $F$ is small in the inclusion, the last term can be neglected and we are left with $j_{0}\approx\sqrt{C}F^{2}(0)$. This leads to Eq. (8).

.3 Numerical analysis of $j_{c}$

To analyze the error associated with calculating $j_{c}$ numerically, we took $L=20$ and varied $\Delta x$. The results are shown in figure 11. Assuming the error is linear, we extrapolate the critical current to be $j_{c}\approx 0.06366$, which is in excellent agreement with the linear system solved using the shooting method with $\Delta x=0.001$.

For fixed $\Delta x=0.05$, we measured the sensitivity of $L$ on $j_{c}$ and found no significant change for $L\gg r$ (typically $L>5r$ was sufficient).

.4 Weakly nonlinear calculation

To obtain the weakly nonlinear system, we analyze near $j_{0}=j_{c}+\epsilon$ where $|\epsilon|\ll 1$. Linearizing about the base state near $\epsilon=0$ with $\Psi=(F+\eta)e^{i\phi}$. From before, we saw that $\epsilon=0$ leads to a degenerate zero eigenvalue implying that the linearized system has a generalized eigenvector solution where $\mathcal{L}\eta_{1}=0$ and $\mathcal{L}\eta_{2}=\eta_{1}$. We use Ansatz $\eta=A\delta\eta_{1}+\delta^{2}B(\eta_{2}+z)+\delta^{3}\zeta$ where $\eta_{k}=\left(\begin{matrix}U_{k}\\ V_{k}\end{matrix}\right)$ and $\epsilon=\alpha\delta^{2}$. Inserting this into (4a)–(4c) with the aid of mathematica and obtain at first order the ODE for $A$

[TABLE]

At next order, we obtain the ODE for $B$ (where we have already projected onto the eigenvector)

[TABLE]