Spin-polarized superconductivity: order parameter topology, current dissipation, and double-period Josephson effect

We discuss transport properties of fully spin-polarized triplet superconductors, where only electrons of one spin component (along a certain axis) are paired. Due to the structure of the order parameter space, wherein phase and spin rotations are intertwined, a configuration where the superconducting phase winds by $4\pi$ in space is topologically equivalent to a configuration with no phase winding. This opens the possibility of supercurrent relaxation by a smooth deformation of the order parameter, where the order parameter remains non-zero at any point in space throughout the entire process. During the process, a spin texture is formed. We discuss the conditions for such processes to occur and their physical consequences. In particular, we show that when a voltage is applied, they lead to an unusual alternating-current Josephson effect with double the usual period. These conclusions are substantiated in a simple time-dependent Ginzburg-Landau model for the dynamics of the order parameter. Our analysis is potentially applicable to twisted double-bilayer graphene, where evidence for triplet superconductivity in proximity to ferromagnetism was recently discovered.


I. INTRODUCTION
Two-dimensional moiré materials, such as twisted bilayer graphene (TBG) heterostructures based on other two-dimensional van der Waals materials, have recently emerged as a fertile ground for novel many-body phenomena . In particular, recent experiments in twisted double bilayer graphene (TDBG) with a perpendicular electric field revealed a ferromagnetic ground state at half filling of the moiré lattice [23,24]. Upon changing the density away from half filling, evidence for superconductivity was found [23,25]. Most strikingly, the superconducting onset temperature increases linearly as a function of an in-plane magnetic field for small fields, a clear indication for triplet superconductivity [23,[26][27][28][29][30]. Moreover, since the superconductor (SC) emerges from a ferromagnetic state, it is natural to expect that the spin polarization persists into the superconducting state [31].
These remarkable findings motivate us to reexamine the physics of triplet SCs in systems with negligible spinorbit coupling. We focus on a particular state which is natural for TDBG: a fully spin-polarized triplet superconductor, in which only electrons of one spin component (along a spontaneously chosen quantization axis) are paired. This phase is analogous to the A 1 and β phases discussed in the context of superfluid 3 He [32][33][34]. The order parameter of such a superconductor is specified by the spin direction of the condensate and its phase; topologically, the order parameter space is equivalent to the space of three-dimensional rotation matrices, SO(3) [35]; see Fig. 1(b). This property crucially determines the topology of the order parameter space, and hence the ability of the superconductor to carry stable supercurrents.
In a ring geometry, the superconducting phase of an ordinary (e.g., singlet) superconductor can wind any integer number of times around the hole of the ring. Math-  1. (a) A depiction of the proposed experimental setup; a fully spin-polarized triplet superconductor is connected to a DC voltage source and drain via the left and right leads. We overlay a current-carrying configuration of the order parameter. (b) A visualization of the order parameter triad where d (1) and d (2) are depicted by the red and blue arrows and the pairing polarization m is depicted by the gray arrow. (c) Schematics of the current response of the system. ematically, this is expressed by the homotopy group π 1 [U (1)] = Z. Changing the winding number requires creating a topological defect at which the magnitude of the order parameter is suppressed to zero, such as a phase slip or a vortex [36,37]. Since such defects are energetically suppressed, a configuration with a non-zero winding number (corresponding to a supercurrent) is metastable and may persist over a very long time.
In contrast, the homotopy group the order parameter of a fully spin-polarized triplet SC on a ring is π 1 [SO(3)] = Z 2 . Physically, this means that in such a SC, two vortices can always annihilate each other, regardless of their vorticity. A configuration in which the SC phase winds twice, either around a point or in a ring geometry, can be "untwisted" back to a uniform configuration continuously, without diminishing the magnitude of the order parameter at any location. The untwisting process involves a temporary change in the system's magnetization; thus, the timescale for this process to occur depends on the coupling of the SC condensate to either an intrinsic or extrinsic bath with which it can exchange spin angular momentum.
The purpose of this work is to explore the unusual transport properties of fully spin-polarized triplet SCs that result from the topology of their order parameter space. We find that in such SCs, a supercurrent-carrying state is much more fragile than in an ordinary SC. The critical current density, J c , that the system can sustain in a metastable state scales inversely with the system size along the current direction, and depends on both the spin and phase stiffnesses. A small voltage applied across the superconductor results in a time-dependent current with a direct-current (DC) component of magnitude close to this critical current, and an alternating-current (AC) component with a fundamental frequencyω J = eV / , i.e., half of the usual Josephson frequency across a SC weak link; see Fig. 1(c). This doubled periodicity is directly related to the Z 2 topological structure of the order parameter space.
We demonstrate these phenomena within a timedependent Ginzburg-Landau model, where we assume that the system is coupled to a bath with which it can exchange both energy and spin angular momentum. For currents larger than the critical current density mentioned above, J c , the model provides a prediction for the current-voltage relation: V ∝ (J − J c ) 2 , up to logarithmic corrections, where J is the DC component of the current. When an external Zeeman field is applied (e.g., an in-plane magnetic field in a two-dimensional system), it pins the direction of the spin magnetization, and the properties of the system rapidly cross over to those of an ordinary superconductor. We hope that these predictions will provide guidance to experiments in novel exotic superconductors, such as TDBG, where they can be used to confirm or invalidate the existence of a fully spin-polarized triplet SC. This paper is organized as follows. In Sec. II we review the properties of a fully-spin polarized SC, its distinction from other triplet SC phases, and the topology of its order parameter space. We then describe the physi-cal picture and summarize the main results of the paper. Sec. III describes the Ginzburg-Landau (GL) model. In Sec. IV we study the energy landscape of the model, either with or without an applied Zeeman field, deriving the maximum supercurrent the system can carry in a metastable state. Sec. V introduces a time-dependent extension of the GL model. This model allows us to study dynamic phenomena, such as the properties of the system in the presence of a finite applied voltage. The results are discussed in Sec. VI. The Appendices contain technical details of the solution of the time-dependent Ginzburg-Landau (TDGL) equations.

II. PHYSICAL PICTURE & MAIN RESULTS
In order to gain physical intuition of spin-polarized triplet superconductivity, we begin with a short review of the relevant order parameters. Those who are already familiar with this material may opt to skip to the results which are presented in Sec. II B.
A. Review of triplet superconductivity and application to TDBG

Order parameter
The order parameter of a triplet SC can be represented in terms of the so-called d-vector, which is a complex vector defined as , and the Pauli matrices σ = (σ x , σ y , σ z ) act on the spin degrees of freedom (↑, ↓). The Pauli principle forces d k = −d −k . [38] The existence of multiple valleys in systems such as TDBG enables the possibility of a valley-singlet state which does not require the order parameter to be momentum dependent within each valley. In the context of TDBG, we denote creation operators of electrons at the K + and K − valleys by the two-component spinors respectively. An inter-valley, valley-singlet superconducting order parameter corresponds to where d is independent of momentum within each valley. For simplicity, we will focus on this type of state henceforth. The d-vector may be conveniently described using two real vectors d (1,2) , We define the pairing polarization m and superfluid density ρ as FIG. 2. Continuous trajectory for supercurrent relaxation. We present the trajectory of Eq. (30) relaxing a supercurrent of J = 2J0 at t = 0 to a supercurrent of J = 0 at t = τ . The order parameter components, d (1) and d (2) , are depicted by the red and blue arrows, the gray arrows depict the pairing polarization, m; see Fig. 1 see Fig. 1(b). At temperatures far below the mean-field superconducting phase transition, we may treat the problem within the London limit where the magnitude of the order parameter is essentially fixed, Note that the pairing polarization m arises from the anomalous expectation value d in Eq. (1). At equilibrium, the pairing polarization is coupled to the physical magnetization (spin polarization) density, and proportional to it. This relation also holds for slowly varying magnetization textures. Thus, below we will discuss the behavior of the pairing polarization interchangeably with that of the magnetization density.

Symmetry classification of valley-singlet states
Triplet superconducting phases have been studied in great detail in the context of superfluid 3 He; see Refs. 32 and 34 for a review. These works naturally focused on pwave phases. However, as discussed above, in materials such as TDBG there is a possibility of a valley-singlet phase. Nevertheless, a similar symmetry classification of the possible phases can be carried out.
Neglecting spin-orbit coupling, the free energy of the system is invariant under transformations in the SO(3) S group of global spin rotations and in the group of U (1) φ gauge transformations. We thus distinguish three distinct triplet phases, according to their patterns of broken symmetries: (i) Non-spin-polarized, unitary spin nematic phase.-In this phase, the pairing polarization vanishes, m = 0, corresponding to d (1) || d (2) . The symmetry is reduced to the group of SO(2) S spin rotations around a preferred direction and the [Z 2 ] S+φ operation of a simultaneous π spin rotation around an axis perpendicular to d 1,2 and a π gauge transformation. The order parameter degeneracy space is given by This phase has the same pattern of spin and gauge symmetry breaking as the A and B 2 phases of 3 He [32,33]. We note, however, that the phases of 3 He differ by the breaking of orbital symmetry which is absent for the valley-singlet phase described here.
(ii) Fully spin-polarized non-unitary triplet.-In this phase d · d = 0, corresponding to d (1) ⊥ d (2) and |d (1) | = |d (2) |. The condition d·d = 0 can be interpreted as the vanishing of a charge 4e scalar order parameter ∆ 4e = d·d. The symmetry is reduced to that of U (1) S+φ rotations around the pairing polarization direction, m, matched with gauge transformations by the same phase. The order parameter degeneracy space is given by This situation is similar to the A 1 and β phases of 3 He [32,33], up to orbital symmetries (see above).
(iii) Partially spin-polarized non-unitary triplet.-In this phase neither m nor ∆ 4e vanish. The symmetry is completely broken and the degeneracy space is given by This situation is similar to the B and A 2 phases of 3 He [32,33], up to orbital symmetries (see above).
Recent theoretical works [26,27] have pointed out that the order parameter of TDBG might be in a fully spinpolarized state, at least in a part of the phase diagram. Motivated by this possibility, together with a fundamental interest in the novel properties that we uncover below, in the remainder of this work we will focus on the properties of fully spin-polarized triplet SCs [phase (ii) above].

B. Supercurrent decay mechanisms
The robustness of persistent supercurrents in ordinary SCs has a topological origin: in a superconducting ring, the winding number of the phase of the order parameter around the ring is an integer that cannot be changed by a small, local perturbation. For fully spin-polarized triplet SCs, the SO(3) topology of the order parameter configuration space, has important consequences for the stability of the supercurrent. Specifically, consider an initial current-carrying configuration with a constant pairing polarization, and a certain number n of windings of the superconducting phase across the system. A continuous relaxation event, such as the one illustrated in Fig. 2, can reduce the phase winding to n − 2. The unwinding event may involve a finite free energy barrier, which depends on the initial number of windings, n, on the ratio between the phase and pairing polarization stiffnesses, and on the applied Zeeman field. The barrier originates from the fact that the unwinding process involves twisting the pairing polarization in space and changing its average direction. We analyze the energy landscape of the continuous unwinding process in Sec. IV. The result is that the energy barrier vanishes above a certain current density, which we define as the critical current J c . The critical current has two contributions, where µ is the magnetic moment density in the material, L x is the length of the sample along the current direction, and B is the applied Zeeman field. Here, c is some non-universal constant which is nevertheless limited to the range 1 ≤ c ≤ 2. Intuitively, in the absence of a Zeeman field, we find only a mesoscopic contribution that scales inversely with system size, corresponding to a finite (fixed) number of windings. On the contrary, the presence of a Zeeman field acts to lock the pairing polarization, increasing the critical current. For a large enough Zeeman field, J c exceeds the intrinsic (microscopic) critical current density; continuous unwinding events then become irrelevant, and the system behaves as an ordinary SC. A direct physical manifestation of the order parameter topology can be observed when a voltage is applied across a fully spin-polarized SC [ Fig. 1(a)]. The voltage causes the superconducting phase at the right lead to wind relative to the phase at the left lead, at a constant rate ω J = 2eV / . Once the phase difference across the system has increased by 4π, there is no topological obstruction to continuously "unwind" the phase twist, relaxing the supercurrent. After the unwinding event, the phase difference (and hence the supercurrent) starts growing linearly in time again. Assuming that the dissipation rate is large compared to the Josephson frequency (a condition to be discussed further in Sec. VI C), we find that the current undergoes periodic oscillations with frequencyω J = eV / , i.e., half of the ordinary AC Josephson frequency [see Fig. 1(c)]. This is one of our main results. We substantiate it by solving a time-dependent Ginzburg-Landau model in Sec. V.

A. Complex vector description
The starting point of our analysis of a spin-polarized triplet SC is most generic gauge-invariant form of the of the free energy, expanded in long wavelengths up to second order in derivatives, as a function of the complex d-vector order parameter: Here, κ d is a generalized phase stiffness, κ m is the excess spin stiffness, B is the applied Zeeman field, µ is the magnetic moment density, and the potential U (d) can be any gauge-invariant function [39]. In order for the free energy to be bounded from below, the generalized phase stiffness, κ d , and the excess spin stiffness, κ m , must satisfy κ d ≥ 0 and κ m ≥ −κ d . We have opted for the Weylgauge, i.e., zero electric scalar potential, φ = 0, such that the system is coupled to the magnetic vector potential A; for generality we keep the charge, q, generic (in electronic SCs q = −2e). For the purpose studying dynamics, we postulate the following Poisson bracket between the components of d: These give rise to the following relations: consistent with m and ρ being the generators of spin rotations and gauge transformations, respectively. The supercurrent density is given by

B. Unitary matrix description
A very useful description of the SO(3) topology for the purpose of explicit computations is through its SU (2) double cover. This description is applicable only for the fully spin-polarized phase (see Sec. II A 2), where it captures the three-dimensional nature of the degeneracy space. For every matrix u ∈ SU (2) we can associate a complex d-vector The double cover is evident by the identical order parameters corresponding to both u and −u. The homotopy expressed in Eq. (8) follows directly from this property (see discussion in Sec. IV). This description holds in the London limit where |d| 2 = 2d 2 0 . Gauge transformations are implemented as where Λ(r) is a scalar function. On the other hand, a spin rotation of angle θ around an axisn acts by The most general form of the free energy invariant under gauge and spin rotation transformations is given by where κ d and κ m are the same coefficients as in Eq. (10). Within this description, the supercurrent density [Eq. (13)] and the pairing polarization [Eq. (3)] are given by

IV. ENERGY LANDSCAPE
We now consider the energetics of the continuous unwinding process for supercurrent relaxation. In general, a supercurrent carrying configuration would be energetically less favorable than a configuration with a smaller supercurrent. This less favorable configuration may be either a metastable or an unstable state. As long as the higher-energy configuration is metastable, either thermally activated or quantum tunneling events would act to relax the current; however, such processes are expected to be exponentially suppressed. Therefore, as the supercurrent gradually increases (e.g., in response to an applied DC voltage; see Fig. 1), the system should stabilize on a value of the current close to its lowest unstable state.

A. Linear stability
We seek the lowest possible supercurrent for which a supercurrent carrying state ceases to be metastable (i.e., where an instability develops).
Without loss of generality we set B = Bẑ. In the presence of this Zeeman field, the system always has as a fixed point, δF/δd = 0, where the d-vector exhibits a uniform configuration described by This configuration may be set to carry arbitrarily high currents by applying a vector potential A = − 1 2q 2 κ d J; see Eq. (13). In search of an instability, we study small deviations, d = d 0 + δd, from the configuration d 0 in Eq. (19), with where η = (η 1 , η 2 , η 3 ) is real, and |η| 1. Here we have taken the most general deviation respecting the fully spin-polarized SO(3) structure, d * ·d = 2d 2 0 , and d·d = 0, up to second order in η. Expanding the free energy to second order in the deviation η yields When considering a system on a cylinder of height L y and circumference L x with a current J = Jx, one immediately finds the least stable direction (defined with respect to the curvature of δF ) in configuration space to be For a deviation of this form, the corresponding change of the free energy is given by We thus find the critical current for the onset of instability, J c1 , as the value of J for which the curvature of the quadratic dependence of δF on η becomes negative: For J ≥ J c1 , the uniform pairing-polarization fixed-point configuration in Eq. (19) becomes a saddle point, and a spontaneous polarization texture should form. Nevertheless, this does not guaranty the existence of a trajectory in configuration space which connects the system to a new stable fixed-point without traversing an energy barrier. In the next section we thus estimate the possible energy barriers and seek out these trajectories.

B. Unwinding trajectories and energy barriers
We move on to estimate the free energy barrier for relaxing two windings of the superconducting phase continuously, as allowed by the topology of the order parameter [Eq. (8)].
This is easiest to describe using the unitary matrix presentation and by taking the gauge of A = 0. In this presentation, a uniform persistent supercurrent carrying state is given by As before, we keep the cylindrical geometry, whereby the order parameter, d, has periodic boundary conditions along x, which force either periodic or anti-periodic boundary conditions for u. These correspond to the two Z 2 classes of maps from the cylinder to the order parameter space [see Eq. (8) and discussion in Sec. III B]. As a consequence, the supercurrent is quantized to J = J 0 n, where n ∈ Z and is the fundamental current; this configuration is overlaid in Fig. 1(a). From the Z 2 order parameter topology we know there exist trajectories u(r, t) with n windings at t = 0 and n − 2 windings at some later time t = τ , in which the order parameter magnitude remains constant throughout the deformation. Among these, we seek the trajectory with the minimal energy barrier.
We keep B = Bẑ, and look for trajectories with uniform pairing polarization m = d 0ẑ at initial and final times t = 0 and t = τ , respectively: We study the following variational family of trajectories u(r, s 1 , s 2 ), which connect these initial and final states, and parametrically depend on time via s 1,2 (t = 0) = 0 and s 1,2 (t = τ ) = 1: Comparing Eqs. (28)- (30) with Eq. (26), it is evident that this trajectory connects a supercurrent-carrying state with n windings to a state with n − 2 windings. For example, the trajectory s 1 = t/τ , s 2 = 0 is shown in Fig. 2. By inserting the ansatz for u(r, t) in Eq. (30) into Eq. (18), one finds that the variation of s 1 (for s 2 held fixed at 0) acts to dissipate the current but flips the pairing polarization relative to its initial orientation parallel to the external field. The variation of s 2 acts to re-align the pairing polarization. An upper bound on the free energy barrier can be obtained by considering an example trajectory where s 1 rises from 0 to 1 strictly before s 2 does so. Such a trajectory has the appealing property that the free energy density is spatially homogeneous, and is given by We find a simple expression for the maximal energy barrier ∆F [u] along this trajectory: Note that J c1 [Eq. (25)] is in general a lower bound on J c2 , and indeed we find that J c1 ≤ J var c2 . Particularly, for B = 0 our variational bound is fully saturated and our estimate for J c2 is exact.
In terms of this critical current estimate, J var c2 , the free energy barrier in Eq. (32) gives For J < J c2 , current relaxation by continuous unwinding processes is possible, but requires traversing a barrier. In this case, relaxation is expected to be much slower than relaxation for J ≥ J c2 , especially at low temperatures.

C. Mean supercurrent
Upon driving the system, e.g., by applying a DC voltage, the system would oscillate back and forth between a current close to the critical current J ∼ J c2 and the relaxed state with J ∼ J c2 − 2J 0 ; see Fig. 1. Hence, the observable mean current,J c , should be well approximated byJ c J c2 − J 0 . Moreover, since J c1 ≤ J c2 ≤ J var c2 , we may use our results, Eq. (25) and Eq. (33), to get an estimate for the mean critical current: where c is limited to the range 1 ≤ c ≤ 2. Note, that this estimate only holds for J c2 > J 0 , since otherwise, J c2 < |J c2 − 2J 0 |, and there is no energy gained by dissipating two windings. As we show in Sec. V, using dynamical simulations, Eq. (30) provides a good approximation for the continuous unwinding trajectories that naturally emerge from the dynamics of the system, and the critical current is indeed given by (35).

D. Competing mechanisms
Throughout this analysis, we have assumed that the order parameter is uniform in the y direction. There also exist trajectories where the unwinding process occurs in a "domain wall" whose width in the y direction is smaller than L y . The domain wall then propagates along y and unwinds the phase twist in the entire system. However, the optimal width of the domain wall is proportional to L x , and thus its formation is only favorable when L y is sufficiently larger than L x .
Another competing mechanism for current relaxation is by vortex motion perpendicular to the current direction. As discussed in Sec. VI, we find the continuous unwinding mechanism is dominant at sufficiently low currents and low temperatures.

V. DYNAMIC TRANSPORT
So far, we have explored current relaxation from an energetic point of view. In this section, we support the predictions from the analysis of Sec. IV using an unbiased model for the order parameter dynamics that enables the system to explore its configuration space.

A. Time dependent model
To capture the dynamics of the current relaxation process, we use a Time-Dependent Ginzburg-Landau (TDGL) formulation supplemented by a stochastic noise term [40][41][42]. Crucially for the continuous unwinding process, we must include a mechanism for the system to exchange magnetization with its environment (relaxing the constraint of angular momentum conservation). To this end, we assume that the system is coupled to an external spin bath; for further discussion see Sec. VI B.
Within the stochastic TDGL approach, the order parameter evolves according to a Langevin-type equation of motion of the form Here, we have utilized the Einstein summation convention, with α, β ∈ {1, 2} and i, j, k ∈ {1, 2, 3}. The three terms on the right hand side of Eq. (36) are the kinetic (dissipationless) term, the dissipative term, and the source (noise) term, respectively. We take ζ j (r, t) to be a real Gaussian-distributed stochastic field with zero mean and correlation function We choose g α ik (d) = ε ijk d α j , which corresponds to a spin bath acting as a spatiotemporally-fluctuating Zeeman field ζ that induces precession of the order parameter. The form of Eq. (36) guarantees that the fluctuationdissipation theorem is satisfied at thermal equilibrium; see Appendix A.
In terms of the complex d-vector, The Langevin equation takes the form As our system is in the fully spin-polarized state, the potential in the free energy [Eq. (10)] must include a polarizing term, U (d) = λ |d · d| 2 + . . ., with large coupling λ. When taking the limit λ → ∞ we generate the correction to the kinetic term in the first line of Eq. (38).
Importantly, Eq. (38) is manifestly gauge invariant, and, as constructed, leads to the conservation laws By design, the pairing polarization m = 1 2i d * × d, which is proportional to the spin magnetization at equilibrium, is not conserved.

B. Solutions
To compute transport in the system, we consider a cylindrical geometry with circumference L x and height L y , where a current flows due to the presence of an applied electromotive force, equivalent to a DC voltage, V . Within our stochastic model, we thus seek the mean current, J(t) ζ (averaged over all realizations of the noise, ζ), that flows due to a uniform, constant electric field Here, E = Vx/L x . Since the TDGL equations [Eq. (38)] are nonlinear stochastic partial differential equations, their solutions are cumbersome to write. Therefore, we leave the details of the solution to Appendices B and C and present here the main results.
We study the dependence on the model parameters, i.e., the stiffnesses κ d and κ m , sample dimensions L x and L y , fluctuation coefficient Γ, temperature T , and applied voltage V . We focus on certain analytically tractable limits of physical interest. First, in order to avoid the domain-wall formation discussed in Sec. IV, we keep L y sufficiently small compared to L x . Current relaxation by transverse motion of vortices is neglected, assuming that the temperature is sufficiently low (see discussion in Sec. VI C). Under these conditions, the order parameter may be assumed to be independent of y. Second, for simplicity, we study the case of large fluctuation and dissipation, Γ d −2 0 , such that the kinetic term in Eq. (38) is negligible compared to the dissipation and noise terms. The applied Zeeman field is moreover set to zero. Finally, in order to work in the regime where relaxation occurs through well-separated-in-time individual unwinding events, we focus on low temperatures such that k B T L y κ d qV /Γ; see Eq. (B27) in Appendix B. Below, we present the time dependent solution of the TDGL equations, Eq. (38), as well as the the long-timeaveraged mean supercurrent, Note that, as the system is self averaging, the long-time averaged current takes the valueJ for each generic realization,J = lim τ →∞  x ; exact details of the simulation are found in Appendix C. Remarkably, there are no fitting parameters.

Energy landscape predictions
Our analysis of the continuous supercurrent dissipation mechanism and estimates of the critical supercurrent in Sec. IV enable us to give predictions for the DC transport setup discussed above. Suppose we start from a configuration with a uniform pairing polarization and apply a voltage. As long as this configuration is locally stable, the current rises linearly with time. This proceeds until the current exceeds J c . Then, the energy barrier for spontaneous dissipation vanishes and a dissipation event would initiate. This process should repeat itself, thus maintaining the average current close to the critical value, J c , estimated in Eq. (35). At a finite voltage, the current may overshoot the critical value, since the uniform configuration is a saddle point of the free energy (despite being unstable). The escape process is triggered by thermals fluctuations and takes a finite time to occur. The larger the voltage, the more the current overshoots J c , and hence the average current increases. Hence we anticipateJ with some exponent β.
In the following sections, we present both analytic and numeric results (see Figs. 3 and 4) showing these generic predictions are indeed realized in our specific TDGL model.

Large dissipation limit
The simplest case to analyze is the case of Γκ d L 2 x qV → ∞, where the dissipation rate is large compared with the Josephson frequency. In this case, the dependence of the current on time can be entirely found analytically (see Appendix B).
In Fig. 4 we show time-dependent supercurrent profiles J(t) for various values of κ m /κ d , together with the corresponding values of the time-averaged supercurrentJ. For κ m ≤ κ d the current profile is an asymmetric triangle wave pattern which averages out to 0, but for κ m > κ d we have a shifted saw-tooth wave pattern which averages to a finite value. These patterns always exhibit a period of 2t 0 = 2(h/2eV ), exactly matching our predictions for the doubling of the AC Josephson effect; see Sec. II B.
In this limit of large dissipation, the time-averaged supercurrent is independent of voltage and temperature, Due to the effectively infinite dissipation coefficient, the unwinding event begins immediately when the current reaches J c1 , exactly matching our predictions of Eq. (25): there is no "overshoot," even at a finite voltage. There-fore, the finite voltage correction discussed in Sec. V B 1 vanishes and we indeed getJ J c ; cf. Eq. (35). Note that the cusps in current vs. time are special to the infinite dissipation limit, and are smoothed out for finite dissipation; see Fig. 5 in Appendix B.

Finite dissipation
In the generic case of finite dissipation rate, we obtain a power-law relation between the time-averaged current J and the applied voltage. We find that

Transition point
In the limit of infinite dissipation rate, Eq. (44) displays a non-analytic dependence ofJ on the ratio κ m /κ d , with a critical (transition) point at κ m = κ d . For finite dissipation rates, the full equations forJ are somewhat simplified at this transition point. At this point we can obtain a complete solution; we find a characteristic voltage and temperature dependence of the average supercurrent,J This result is in excellent agreement with numerical simulations of Eq. (38), as seen in Fig. 3; the details of the simulations are given in Appendix C. Note that our result for the average current in Eq. (46) naïvely diverges at T → 0. Although the current exceeds a level where there is no longer an energy barrier for the continuous unwinding process, at T = 0 there are no fluctuations to trigger decay within the classical model of Eq. (38). However, at exponentially low temperatures whereJ J 0 , spontaneous quantum decay process uncaptured by our analysis are expected to dominate the relaxation.

VI. DISCUSSION AND CONCLUSIONS
We end with a discussion of various points regarding the implications of our results in the context of polar-ized triplet SCs in general, and in TDBG in particular. In the following, we discuss the thermal phase diagram of two-dimensional spin-polarized SCs, possible mechanisms for relaxation of the magnetization (necessary for our topological mechanism for continuous current relaxation), the role of conventional supercurrent relaxation by vortex-antivortex dissociation, and estimated parameters for triplet superconductivity in TDBG. Finally, we point out some topics worthy of future studies.

A. Thermodynamics of 2D spin-polarized superconductors
Throughout our discussion we have assumed that, at equilibrium, our system is essentially long-range ordered, i.e., d = 0. This is, of course, never strictly true in two spatial dimensions at any non-zero temperature. In the absence of a Zeeman field and any other breaking of spin rotational symmetry, the system is always in a disordered phase at all T > 0 [43], in accordance with the Mermin-Wagner theorem. However, the spin correlation length grows rapidly with decreasing temperature, ξ ∼ e 2π(κm+κ d )/(k B T ) ; hence, at sufficiently low temperatures, ξ max(L x , L y ) and our treatment should apply. If κ m is larger than κ d , there is a finite temperature crossover for finite systems that resembles a Berezinskii-Kosterlitz-Thouless (BKT) transition. In the presence of a Zeeman field that pins the direction of the pairing polarization, the system undergoes a genuine finite-temperature BKT transition.
We note in passing that, in contrast to a fully polarized triplet SC [phase (ii) in Sec. II A 2], a partially polarized phase [phase (iii) in Sec. II A 2] undergoes a true BKT transition at finite temperature even for B = 0, at which ∆ 4e becomes quasi-long range ordered [43]. This phase is interesting in its own right, as it supports half quantum vortices (carrying flux h/4e). We leave a more detailed discussion of this phase and its physical implications for future work.

B. Mechanism for magnetization relaxation
The continuous supercurrent relaxation process discussed in this work requires a mechanism for the pairing polarization to dynamically change in time (both locally and globally). Since the pairing polarization is pinned to the physical spin polarization, this requires a mechanism for the system to exchange spin angular momentum with its environment. In our time-dependent model (Sec. V), we assumed the existence of a spin bath on empirical grounds. In practice, the maximum rate at which the continuous relaxation process can proceed, 1/τ s , depends on the physical mechanism of spin relaxation. (Within our model, 1/τ s corresponds to the dissipation rate, Eq. B2.) The system's spin can relax either through spin-orbit coupling (which also introduces anisotropic terms in spin space into the free energy), a spin bath due to internal or substrate impurities or nuclear spins, or through the boundaries, in cases where the system is coupled to metallic leads.
In the context of TDBG, spin-orbit coupling is expected to be weak, and nuclear spins are rare (unless 13 C impurities are introduced intentionally). In the absence of other magnetic impurities, the most dominant source of spin relaxation is likely by spin transport to the metallic leads. This mechanism is absent in our simple dynamical model of Sec. V and we leave its detailed study to follow up work. An appropriate description of this mechanism might be the use of either solitons of the degeneracy space [44,45] or the Goldstone modes of the order parameter in order to carry the magnetic texture to the leads.
Here, we make a rough estimate of the spin relaxation time due to this mechanism. Given that the spin stiffness in our model is κ m + κ d [see Eq. (B4)], we expect on dimensional grounds that 1/τ s ∼ where ξ m is a microscopic "magnetic coherence length" (essentially, µ B /ξ 2 m is the magnetization density) and L x is the distance between the leads. This estimate for τ s can be understood as the time it takes for a quadratically dispersing magnon with wavevector |Q| ∼ 1/L x to travel across the system.
For a sufficiently small applied voltage such that eV 4π /τ s , the continuous unwinding process can take place within one (doubled) Josephson period, with spin carried in and out of the system through its connection to the leads. For larger voltages, the continuous unwinding process is limited by the rate of spin relaxation. The study of this case goes beyond our present analysis. A possible scenario, in this case, is that the phase accumulates more than two windings between the continuous unwinding events, decreasing the effective Josephson frequency belowω J = eV / .
To get a rough estimate of 1/τ s due to spin dissipation at the boundaries in TDBG, we assume that the spin stiffness κ m + κ d is of the order of 0.1 − 1 meV and ξ m is of the order of a few times the moiré lattice spacing a ≈ 10 nm. For example, suppose that κ m + κ d = 1 meV and ξ m = 5a. Then, for a system of size L x = 1 µm, the above considerations give that the crossover voltage is V ≈ 4π /(eτ s ) ≈ 30 µV.

C. Phase unwinding due to vortex motion
In our analysis of the supercurrent relaxation through continuous unwinding, we have neglected the ordinary mechanism of vortex-antivortex dissociation and motion of free vortices. In two spatial dimensions, this mechanism gives rise to the celebrated nonlinear I-V characteristics [36,37,[46][47][48] associated with the BKT transition: V ∝ I α with α(T ) = πκs(T ) k B T + 1, where the phase stiffness κ s satisfies κ s (T BKT ) = 2 π k B T BKT such that α(T BKT ) = 3 (in our model κ d = 1 2 κ s ). In order to determine which mechanism dominates the supercurrent relaxation, we need to compare the rate of phase unwinding (i.e., the voltage) generated by vortex motion to the rate of the continuous unwinding process. The voltage due to vortex motion is estimated as [36] (47) Here, ρ n is the normal state resistivity, J is the current density, ξ is the coherence length, L = min(L x , L y ), and J c,0 ∼ e k B TBKT ξ is microscopic critical current density. The continuous unwinding mechanism occurs when J J c [whereJ c is estimated in Eq. (35)], whereby it does not involve passing through an energy barrier. Hence, the continuous unwinding mechanism is dominant when V V vm (J c ). Therefore, as long asJ c < J c,0 , one sees that V vm (J c ) → 0 at low temperature, since the vortexantivortex unbinding process is thermally activated.

D. TDBG parameters
Specifically for the case of TDBG we may estimate the values of the various parameters in our model using the experimental data of Ref. 23. The sample dimensions in these experiments were of the order of a few µm, with a carrier density of n ≈ 2 × 10 12 cm −2 . As discussed in Sec. II A, we focus on dynamics far below the critical temperature T BKT ≈ 3.5 K. For an order of magnitude estimate, we set κ d ∼ 1 2 κ s (T BKT ) = 1 π k B T BKT ≈ 0.1 meV. We furthermore estimate the magnetic moment density as µ ≈ nµ B .
Using these estimated parameter values, we assess the two contributions to the critical current in Eq. (9). For small applied Zeeman field we have I micro

E. Conclusions and future directions
In this paper, we have studied the transport properties of fully spin-polarized triplet SCs. Interestingly, we found that the topology of the order parameter degeneracy space can have dramatic consequences, given that a mechanism for spin dissipation is available. Basically, two windings of the order parameter are topologically equivalent to zero windings, and hence a supercurrent can decay via a smooth deformation of the order parameter that involves a transient inhomogeneous spin texture but does not involve creating topological defects (such as vortices). We expect that this mechanism of supercurrent relaxation may become dominant at sufficiently low temperatures, where vortex configurations are suppressed.
A direct observable manifestation of the continuous supercurrent relaxation mechanism is that a sample subjected to a DC voltage exhibits an oscillatory current with a fundamental frequency which is half the Josephson frequency. Furthermore, an applied Zeeman field suppresses the continuous supercurrent decay mechanism, since it pins the direction of the system's magnetization.
A possible extension of our analysis is the study of the equilibrium properties of the system in the presence of a perpendicular magnetic field. At high applied field one would find the standard vortex lattice groundstate. However, at low magnetic fields, one would suspect a stabilization of a vortex-free spin texture; this would be reminiscent of Fig. 2 in space rather than in time. We leave the analysis to a followup work. In this Appendix, we show that the Langevin equation [Eq. (36)] satisfies the fluctuation-dissipation theorem (FDT) at equilibrium. Recall that a sufficient condition for a stationary stochastic process to satisfy the FDT is that the Gibbs distribution, W ∝ e −F/k B T , is time-independent under the stochastic process [49]. This means that the Gibbs distribution is a solution of the corresponding Fokker-Planck equation.
We consider a set of dynamical variables X I which satisfy general Langevin equations of the form (using the Stratonovich convention [42]) where repeated indices are summed over; h I , g Iµ are differentiable functions, and ζ µ (t) are independent Gaussian noise sources with zero mean, satisfying ζ µ (t)ζ ν (t) = k B T Γδ µν δ(t − t ). The corresponding Fokker-Planck equation that describes the evolution of the distribution function W ({X}, t) is [42]: Applying this rule to our Langevin equation (36) with i (r), we obtain: In this appendix we present the details of the analytic solution to Eq. (38).

Preliminaries
We focus on the case of large dissipation Γ d −2 0 . In this case there are three frequency scales in the system. The inverse of the Josephson period, the dissipation rate, and the fluctuation rate, Whenever possible we shall express our equations and results in terms of these quantities.

The equations
We explicitly take d · d = 0, characteristic of the fully spin-polarized phase, as discussed in Sec. II A 2; this condition is conserved by the equations of motion. We furthermore take the London limit of d · d * = 2d 2 0 . These conditions imply U (d) = const and ∇ · J = −ρ = 0. This allows one to rewrite the free energy as and to explicitly write the equations of motion [Eq. (38)], This nonlinear partial differential equation is best solved by the method of guessing the solution. We use an ansatz inspired by the unwinding trajectory of Sec. IV, Here, n and ∆ are integers, while φ n and φ ∆ are some arbitrary phases, and u 0 is an arbitrary SU (2) matrix. The current is given by This trajectory connects two families of constant configurations; f − = − π 2 + 2π with n − ∆ windings, and f + = π 2 + 2π with n + ∆ windings, where ∈ Z. Note that these configurations have uniform pairingpolarization and hence the winding number is unambiguously defined. Both configurations have supercurrent profiles growing linearly with time, The dynamics of the system consist of a series of unwinding events with different values of parameters n, ∆, φ n,∆ , u 0 . The pairing-polarization and supercurrent at the end of an unwinding event determine u 0 and n − ∆ of the next event. We assume the system is initialized with J(t = 0) = 0 such that n = ∆ for the first unwinding event.

Zero temperature
First, we analyze the case of T = 0, whereby our ansatz [Eq. (B6)] reduces the equations of motion to an ordinary differential equation (ODE), (B8) There are two families of fixed-point solutions to this equation, corresponding to the constant configurations discussed above, i.e., f ± = ± π 2 + 2π , where ∈ Z. However, when t → ∞ only f + are stable fixed-points under small perturbations while f − are unstable; when t → −∞ only f − are stable fixed-points while f + are unstable. Therefore, a solution to this equation describes a trajectory from f − to f + . This trajectory describes relaxation of 2∆ windings, and the fundamental 2-windings relaxation trajectory is attained for ∆ = 1.
At an intermediate time, t init , when J(t init ) = J c1 [see Eq. (25)], the stable fixed point f − at t → −∞ becomes a saddle-point. Any small deviation would thus send the system on a trajectory towards the stable fixed point f + at t → ∞, which itself ceases being a saddle-point only at some other intermediate time (t fin discussed below).
However, at T = 0, there is no source for small deviations and the system would thus remain frozen at f − . Moreover, even given some initial perturbation, it could only trigger the first unwinding event and the system would eventually remain asymptotically close to f + with nothing to trigger the next unwinding event. Therefore it is crucial to study the effects of finite temperature.

Effective model for finite temperatures
At finite temperatures, T > 0, the thermal fluctuations would nudge the system from its saddle-point and initiate an unwinding trajectory. The system may spontaneously choose any random fluctuation direction v(x) within the manifold of configurations escaping the saddlepoint. Without loss of generality, we pick such a direction v(x) compatible with our ansatz [Eq.
This is supported by numerical simulations of Eq. (B5) for finite temperatures.
This choice of direction enables us to drastically simplify the functional vector noise term in Eq. (36) and replace it with an effective uniform scalar noise term, g α ij (d)ζ j (r, t) → ε ijk d α j v k (x)ζ(t), with The corresponding equations of motion are modified accordingly (see Appendix A), As we show in Sec. B 4 a, consistency requires Γ ∆ = Γ/L x L y . Remarkably, this effective noise term propagates well into our ansatz [Eq. (B6)], thus reducing the equations of motion to a stochas-tic ODE. As seen in Fig. 6 in Appendix C, this effective model exquisitely captures the numerical detail of Eq. (38).
a. Evaluation of the diffusion rate Given our original 2+1 dimensional vector model, ζ i (r, t)ζ j (r , t ) = 2k B T Γδ ij δ(r − r )δ(t − t ), (B13) and an effective 0+1 dimensional scalar model fluctuating in the v(x) direction, we wish to find the value of Γ ∆ such that the effective scalar model best approximates the original vector model. Hence, we project ζ(r, t) onto v(x), i.e.
such that and thus find Γ ∆ = Γ/L x L y .

Solutions of the effective model
We are now finally in a position to derive the results of Sec. V B. We solve Eq. (B12) and use Eq. (B7) and Eq. (42) to obtain the average supercurrentJ for various ratios of the rates t −1 0 , γ, E and various ratios of the stiffnesses κ d , κ m .
In Fig. 6 we normalize Fig. 3 by the leading powerlaw behavior and find a very good agreement with our predicted logarithmic corrections.