Soliton formation in the FFLO phase

M. D. Croitoru, Alexandre I. Buzdin

► To cite this version:

M. D. Croitoru, Alexandre I. Buzdin. Soliton formation in the FFLO phase. Physical Review B: Condensed Matter and Materials Physics (1998-2015), 2016, 94 (21), pp.214512 (1-8). 10.1103/Phys-RevB.94.214512. hal-01441825

HAL Id: hal-01441825

https://hal.science/hal-01441825

Submitted on 20 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Soliton formation in the FFLO phase

M. D. Croitoru¹,² and A. I. Buzdin¹,³

¹University of Bordeaux, LOMA, UMR-CNRS 5798, F-33405 Talence, France

²Institute of Theoretical Physics III, University of Bayreuth, D-9440 Bayreuth, Germany

³Institut Universitaire de France, Paris

There is increasing body of experimental evidences of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase in quasi-low-dimensional organic and heavy-fermion superconductors. The emergence of the FFLO phase has been demonstrated mainly based on a thermodynamic quantity or microscopically with spin polarization distribution that exhibit anomalies within the superconducting state in the presence of the in-plane magnetic field. However, the direct observation of superconducting order parameter modulation is so far (still) missing. Within the quasiclassical approach and Ginzburg-Landau formalism we study how the orbital effect of the in-plane field influences the FFLO instability in quasi-one-dimensional superconductors with a sufficiently weak interlayer coupling locking the magnetic flux to Josephson-type vortices. By making use of the continuum limit approximation of the Frenkel-Kontorova model for competing periodicities, we find and characterize the locking behavior of the modulation wave vector, when it remains equal to the magnetic length through some range of values of the external field.

# I. INTRODUCTION

If the superconducting state in high magnetic fields is para-magnetically limited, the phase transition from a homogeneous superconducting state to a normal state may occur either directly, at the Pauli limiting field $H_p = \Delta_0/\sqrt{2}$ with $\Delta_0$ superconducting gap, or via an interposed wedged spatially inhomogeneous state, characterized by the simultaneous coexistence of normal and superconducting phases and known as the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [1]. The latter is predicted for clean spin-singlet superconductors in the vicinity of the upper critical magnetic field $H_{c2}$, beyond the field set by the Pauli paramagnetic limit as a result of the competition between pairing correlations, favoring antiparallel spin alignment, and the Zeeman effect, favoring parallel spin alignment along the field, when the Pauli pair breaking dominates the orbital pair breaking [2], with the limiting field $H_{c2}^{orb} = \Phi_0/2\pi\xi^2$. In this state the order parameter varies periodically in space and Cooper pairs acquire a nonzero momentum. Another requirement for the FFLO state is that the superconductor should be in the clean limit, $l \gg \xi_0 = \hbar v_F / \pi \Delta_0$, where $l$ is the mean free path of the quasiparticles and $\xi_0$ is the superconducting coherence length at $T = 0$ and $H = 0$ [3-5]. Furthermore, the highly anisotropic Fermi surface [6,7] favors the FFLO phase formation. This is an inherent property of layered conductors that exhibit a highly anisotropic structure and hence have features of systems with reduced dimensionality [8].

There has been much attention focused on discovering the FFLO state in various superconductors under a high magnetic field since its prediction in 1964 [9,10]. Nevertheless, the conclusive evidence for a nonhomogeneous state has been provided only recently. Several signatures indicating the experimental realization of the FFLO state in organic superconductors were observed for the in-plane external magnetic field when the flux penetrates between the layers in the form of Josephson vortices, thus limiting orbital suppression: (i) the anomaly in the thermal conductivity [11], and finer details in the phase diagram obtained by the tunnel diode oscillator rf penetration depth measurements and pulsed field techniques

[12] for the clean organic sample $\lambda$-(BETS)$_2$GaCl$_4$; (ii) the calorimetric and magnetic torque evidence for the appearance of an additional first-order phase transition line within the superconducting phase in the in-plane high-field regime for the organic sample $\kappa$-(BEDT-TTF)$_2$Cu(NCS)$_2$ [13,14]; (iii) evidence for phase transition within the superconducting phase obtained from the local electronic spin polarization in a $^{13}$C NMR spectroscopy experiment [15], NMR detection of spin-polarized quasiparticles [16] forming the Andreev bound states spatially localized in the nodes of the order parameter [6,17,18], and phase transitions in the $(H-T)$ phase diagram that are consistent with the FFLO phase obtained from rf penetration depth measurements [19] for an organic sample $\kappa$-(BEDT-TTF)$_2$Cu(NCS)$_2$; (iv) NMR relaxation rate evidence for an additional phase transition line [20] and the clear upturn beyond the Pauli limit in the magnetic field and angular-dependent high-resolution specific-heat measurements for the organic material $\beta''$-(ET)$_2$SF$_5$CH$_2$CF$_2$SO$_3$ [21]; (v) magnetic torque evidence for the tricritical point between the FFLO, homogeneous superconducting, and paramagnetic metallic phases in the 2D magnetic-field-induced organic superconductor $\lambda$-(BETS)$_2$FeCl$_4$ [22]; and (vi) measurements of the temperature and angular dependencies of $H_{c2}$ in pnictide superconductor LiFeAs [23] and KFe$_2$As$_2$ [24] probably suggesting that these systems may have realized the FFLO phase as well.

Similarly, quasi-one-dimensional (1D) superconducting compounds have been studied extensively for exhibiting the inhomogeneous state. The field-amplitude and field-angle dependence of the superconducting transition temperature $T_c(H)$ of the organic superconductor (TMTSF)$_2$ClO$_4$ in a magnetic field applied along the conduction planes have been reported [25]. The authors observed an upturn of the curve of the upper critical field at low temperatures. (An enhancement of almost two times over $H_p \simeq 27$ kOe is observed, $H_{c2} \simeq 50$ kOe [25].) Moreover, an unusual in-plane anisotropy of $H_{c2}$ in the high-field regime was observed. Both observations were interpreted as an evidence of the FFLO state formation [26].

The above overview of the recent experimental developments in hunting inhomogeneous states shows that so far
mainly thermodynamic and microscopic due to spin relaxation (NMR) evidences for the existence of an FFLO state have been provided. However, the direct evidence of FFLO states has not been observed yet in experiments. In Ref. [26] we have demonstrated that the temperature evolution of the in-plane critical magnetic field anisotropy as well as of its fine structure may give important information about the FFLO state and unambiguously prove its existence. Our subsequent analysis [27] of the properties of the organic superconductors revealed that in some of them the resonance between the period of the FFLO modulation and the period of the interlayer coupling modulation induced by the external field is possible, which results in anomalous cusps in the temperature and angular dependencies of the in-plane critical field. Furthermore, at the resonances the interplay between the orbital and the paramagnetic effects may result in a commensurability effect between the FFLO and magnetic wave vectors [22,28]. These effects are directly based on the main feature of the inhomogeneous state, the spatial modulations of the order parameter.

In this work we keep studying the effects of orbital contribution to the FFLO modulation, leading to unambiguous proof of the FFLO phase. We study the orbital effect in quasi-1D superconductors under the influence of the external magnetic field and demonstrate that in the FFLO phase there are two competing length scales: the period of the magnetic field potential and the characteristic length associated with the FFLO modulation [29]. The interplay between these two lengths results in locking phenomena, when the modulation wave vector remains equal to the period (or matches m periods) of the potential through some range of the magnetic field values. To study effects of competing periodicities in the system we use the Frenkel-Kontorova (FK) model [30] and show that the solution in the proximity of the locking transition can be described by $\cos[Qx + \varphi(x)]$, with $\varphi(x)$ a phase soliton lattice. The surprising ability of this model to describe many physically important phenomena, such as the dynamics of absorbed layers of atoms on crystal surfaces [31], charge-density waves [32], ferro- or antiferromagnetics [33], Josephson junctions chains [34], transport properties of vortices [35], and superionic conductors [36], has attracted much attention in condensed matter physics and nonlinear dynamics. With this work we further extend this series of described phenomena.

## II. THEORETICAL MODEL

We consider a quasi-1D conductor with the following electron spectrum:

$$E_{\mathbf{p}} = p_x^2/2m_x + 2t_y \cos(p_y d_y) + 2t_z \cos(p_z d_z), \quad (1)$$

where $d_y$ and $d_z$ are the interchain distances along the y and z axis, respectively. We assume that the couplings between chains are small, i.e., $t_z < T_{c0}$ and $t_y < T_{c0}$, but sufficiently large to suppress the charge- and spin-density-wave transitions, to stabilize the superconducting long-range order, and to make the mean-field treatment justified, $T_{c0}^2/E_F \ll t_z$, $T_{c0}^2/E_F \ll t_y$ [37]. Here $T_{c0}$ is the critical temperature of the system at $H = 0$.

In quasi-1D superconductors the orbital effect is extremely weak for the magnetic field applied along the chains [38].

Indeed, near $T_{c0}$ the parallel (along the x axis) upper critical field is $H_{c2}^x = \frac{\phi_0}{2\pi\xi_y\xi_z} \frac{T_{c0}-T}{T_{c0}}$, where $\xi_y = d_y \frac{t_y}{T_{c0}}$ and $\xi_z = d_z \frac{t_z}{T_{c0}}$, while the perpendicular upper critical field (i.e., along the y axis) is $H_{c2}^y = \frac{\phi_0}{2\pi\xi_x\xi_z} \frac{T_{c0}-T}{T_{c0}}$, with $\xi_x \sim \frac{t_y}{T_{c0}}$. We see that $H_{c2}^x \gg H_{c2}^y$ and the orbital effect for the parallel component of the magnetic field is weakened by the factor $\frac{t_y}{T_{c0}} \frac{d_y}{\xi_x} \ll 1$. Therefore, in our analysis of the orbital effect it is enough to take into account only the perpendicular component of the magnetic field because it provides the dominant contribution to the orbital effect. For better clarity we consider the magnetic field to be aligned in the xy plane, making the angle $\theta$ with the x axis. The perpendicular component of the magnetic field $H_y = H \sin\theta$ and the corresponding vector potential may be chosen as $A_z = -xH_y = -xH \sin\theta$, $A_y = 0$, $A_x = 0$. The vector potential varies slowly at the interchain distances, because the diamagnetic screening currents are very weak and the magnetic induction is practically uniform and equal to the external field, $H$. With the help of the Eilenberger equations we obtain the following linearized Ginzburg-Landau equation (for details, see the Appendix):

$$\begin{aligned} \alpha \Delta(x) = & \beta \partial_x^2 \Delta(x) - \delta \partial_x^4 \Delta(x) - 48 t^4 \nu \sin^4(Qx) \Delta(x) \\ & - \frac{1}{2} t^2 \mu [1 - 7 \cos(2Qx)] Q^2 \Delta(x) \\ & + 2t^2 \mu [3 \sin(2Qx) Q \partial_x \Delta(x) \\ & + 3 \sin^2(Qx) Q^2 \Delta(x)] \\ & - 4t^2 \gamma [1 - \cos(2Qx)] \Delta(x). \end{aligned} \quad (2)$$

Here $Q = -\pi d_z H_y / \phi_0 = -\pi d_z H \sin\theta / \phi_0$, with $\phi_0 = \pi \hbar c / e$, $t \equiv t_z$, $\alpha \equiv \ln \frac{T_c}{T_{c0}} - [K_1(T_c) - K_1^0(T_c)] = \frac{T_c - T_{c0}}{AT_c}$, and $\alpha < 0$ in the uniform superconducting state. Equation (2) contains the magnetic-field-induced potential, which is periodic in real space with the period $\lambda_H = 2\phi_0/d_z H_y$. The corresponding functional reads as

$$F_{sn} = \frac{1}{L_x} \int dx \{ \alpha |\Delta(x)|^2 + \beta |\partial_x \Delta(x)|^2 + \delta |\partial_x^2 \Delta(x)|^2 \\ + \frac{1}{2} t^2 \mu [1 - 7 \cos(2Qx)] Q^2 |\Delta(x)|^2 \\ + 6t^2 \mu \sin^2(Qx) |\partial_x \Delta(x)|^2 \\ + 4t^2 \gamma [1 - \cos(2Qx)] |\Delta(x)|^2 \}, \quad (3)$$

where the coefficients are given in the Appendix [39]. The functional (3) is obtained for $t < v_F Q$ and in this limit the term $48t^4 \nu \sin^4(Qx) \Delta(x)$ can be simply neglected.

Without the orbital contribution in 1D, the solution of the Ginzburg-Landau equation is an elliptic sine function (see Appendix) [40], which in the proximity of the line of the transition takes the simple form $\Delta_0 \cos(qx)$, where $q$ is the absolute value of the modulation wave vector. The linearized Ginzburg-Landau functional (3) and Eq. (2) in the paramagnetic limit (when neglecting the orbital contribution) reduce to

$$F_{sn}^{1D} = \frac{1}{L_x} \int dx \{ \alpha |\Delta(x)|^2 + \beta |\partial_x \Delta(x)|^2 + \delta |\partial_x^2 \Delta(x)|^2 \} \quad (4)$$
and

$$ \alpha\Delta(x) = \beta\partial_x^2\Delta(x) - \delta\partial_x^4\Delta(x). \quad (5) $$

The expression for the modulation wave vector is readily obtained from the solution of the linearized equation in the paramagnetic limit, $\Delta(x) = \Delta_0 \cos(qx)$, $q^2 = -\beta/2\delta = 2K_3(T_{cP})/K_5(T_{cP})v_F^2 > 0$. This wave vector increases from 0 at the tricritical point to a very large value at $T \to 0$.

The orbital effect weakly modifies this solution except for the resonance region, $q \approx Q$, where it may strongly influence the FFLO structure [28]. To investigate the influence of the orbital effect on the solution we apply the theory developed by Dzyaloshinskii [41–43]. Hence we look for the solution of the type $\Delta_0 \cos[Qx + \varphi(x)]$, where $\varphi(x)$ is a slowly varying function $\varphi'(x) \ll Q$. This implies that the orbital effect weakly modifies the exact one-dimensional solution, which is justified for $(T_c - T)/T_{c0} \gg (t/T_{c0})^2$. We insert this solution into the functional (3) and perform the average over the period $Q^{-1}$ as

$$ \langle \dots \rangle = \frac{Q}{2\pi} \int_{-\pi/Q}^{\pi/Q} [\dots] dx, \quad (6) $$

and taking into account that $\varphi(x) \ll Qx$, we obtain

$$ \begin{aligned} \langle F_{sn} \rangle_Q = & \frac{\Delta_0^2}{L_x} \int dx \frac{\alpha}{2} + \frac{1}{2} \beta [Q + \varphi']^2 \\ & + \frac{1}{2} \delta[Q + \varphi']^4 + t^2 \gamma [2 - \cos 2\varphi(x)] \\ & + \frac{1}{8} t^2 \mu Q^2 [2 - 7 \cos 2\varphi(x)] \\ & + \frac{6}{8} t^2 \mu [2 - \cos 2\varphi(x)][Q + \varphi']^2. \end{aligned} \quad (7) $$

In the adopted approximation the functional $F_{sn}$ depends only on the space-varying phase, $\varphi(x)$. Since we are interested in the behavior of the system in the vicinity of the resonance $q \approx Q$, we can write the functional in a simpler form (normalized by $-\pi K_5(T_c)\Delta_0^2$, note that $K_5(T_c) < 0$)

$$ \langle F_{sn}^{\varphi} \rangle_Q = \frac{\Delta_0^2}{L_x} \int dx \frac{-\alpha}{2} + \frac{1}{8} v_F^4 Q^2 (\varphi' - \delta q)^2 + \frac{5}{8} t^2 v_F^2 Q^2 - \frac{11}{16} t^2 v_F^2 Q^2 \cos[2\varphi(x)], \quad (8) $$

where $\delta q \equiv Q - q$ is the detuning, or the relative misfit between two periodicities in the system. The phase $\varphi(x)$ is the shift of the FFLO modulation relative to the minima in the potential $\sim \cos[2\varphi(x)]$. The state $\varphi(x) = 0$ is the commensurate phase between $Q$ and $q$. Equation (8) is a continuum limit approximation of the Frenkel-Kontorova Hamiltonian introduced by Frank and Van der Merwe [44].

### III. FFLO SOLITON LATTICE

Taking into account the periodicity of the function $\varphi(x)$ we obtain (for the purpose of clarity the expression is normalized

by $-\pi K_5(T_c)v_F^2 Q^2 \Delta_0^2$

$$ \langle F_{sn}^{\varphi} \rangle_{Q, \varphi} = \frac{1}{K(\kappa^2)} \int_0^{K(\kappa^2)} dx \frac{v_F^2}{8} \varphi'^{,2}(x) - \frac{v_F^2}{4} \delta q \varphi'(x) + \frac{11t^2}{16} \cos[2\varphi(x)]. \quad (9) $$

Therefore, the ground state which minimizes $\langle F_{sn}^{\varphi} \rangle_{Q, \varphi}$, Eq. (9), is given by the solutions of the exactly integrable sine-Gordon equation for the phase $\varphi(x)$

$$ 2 \frac{\partial^2 \varphi}{\partial x^2} + v \sin[2\varphi(x)] = 0, \quad (10) $$

where $v \equiv 11t^2/v_F^2$ is the effective interlayer coupling parameter. In the absence of this coupling, $v = 0$, Eq. (10) has a solution, $\varphi(x) = Cx$, that describes unperturbed one-harmonic incommensurate phase, $\Delta_0 \cos(Q + C)x$. On the other hand, Eq. (10) has a trivial solution, $\varphi = 0$, that corresponds to a commensurate structure $\Delta_0 \cos Qx$. For a finite $v$ Eq. (10) should describe some inhomogeneous distribution of the phase of the order parameter. The first integral of this equation is given by

$$ \left(\frac{\partial \varphi}{\partial x}\right)^2 = \frac{v}{\kappa^2} \{1 - \kappa^2 \sin^2[\varphi(x)]\}, \quad (11) $$

where $\kappa^2 \equiv 2v/(\varepsilon + v)$ and $\varepsilon$ is the constant of integration. Then the exact solution of Eq. (11), expressed in terms of the Jacobi amplitude, is

$$ \varphi(x) = \operatorname{am} \left[ \frac{\sqrt{v}}{\kappa} x, \kappa^2 \right], \quad (12) $$

where $\kappa$ is the constant of integration and must be found from the energy minimum of the system. The solution is the regularly spaced solitons, a soliton lattice. The soliton lattice is a compromise between the term $\cos[2\varphi(x)]$, which favors $\varphi(x) = \text{const}$, and the derivative part that favors $\varphi(x) = Cx$. For $\kappa \to 1$ it reduces to

$$ \varphi(x) = 2 \tan^{-1} \left[ \tanh \frac{\sqrt{v}}{2\kappa} x \right], \quad (13) $$

that describes a domain wall, which separates two commensurate regions, one with phase $\varphi = -\pi/2$ and other with $\varphi = \pi/2$. Within this wall (soliton) a $\pi$ change of the phase occurs. The length of the soliton can be small, meaning the fast change of the phase.

Substituting the solution $\varphi(x)$ for the order parameter into Eq. (8), the normalized energy of the state with the spatial distribution of the order parameter given by Eq. (12) is expressed in terms of complete elliptic integrals $E$ and $K$ as

$$ \langle F_{sn}^{\varphi} \rangle = \frac{\delta q \pi}{4 \kappa K(\kappa^2)} - \frac{\sqrt{v}}{16} \left[ 1 - \frac{2}{\kappa^2} + 4 \frac{E(\kappa^2)}{\kappa^2 K(\kappa^2)} \right]. \quad (14) $$

Here it is necessary to take into account the fact that the function $K_5(T_c) < 0$ for the region of interest. Minimization of the energy with respect to $\kappa$ leads to the following equation
<figure><img src="image_2.png" /><figcaption>FIG. 1. Minimization parameter as a function of the anisotropy parameter.</figcaption></figure>

for $κ$:

$$ \frac{\kappa}{E(\kappa^2)} = \sqrt{\frac{v}{v_c}} \equiv \zeta, \qquad (15) $$

where $v_c = \pi^2 \delta q^2 / 4$. The numerical solution of this equation is given in Fig. 1. The parameter $κ$ varies in the range $0 \le κ \le 1$ as long as $v$ varies in the range $0 \le v \le v_c$. In this plot we use the notation $\zeta = \sqrt{v/v_c}$. The variation of the parameter $κ$ results in drastic changes in the behavior of the phase $\varphi(x)$. Figure 2 displays the phase for several values of the parameter $\zeta \in \{0.2, 0.6, 0.8, 0.9, 0.99, 0.999\}$. As $\zeta$ is increased a plateau section appears within the period of the soliton lattice $L = 2K(\kappa^2)\kappa/\sqrt{v}$ or $\delta qL = (2/\pi)^2 E(\kappa^2)K(\kappa^2)$. With increasing $\zeta$, not only the period of the soliton lattice increases but also the relative fraction of the plateau width. On this plateau the phase remains almost constant, but it strongly changes beyond the plateau at the end of the periods. Thus the system can be represented as a periodic structure of domains of the commensurate/locked phase separated by solitons. We can give to Eq. (14) another form taking into account the distance $L$ between the domain walls,

$$ \langle F_{sn}^{\varphi} \rangle = \frac{\pi}{4} \frac{\delta q}{L} - \frac{v}{16} \left[ 1 - \frac{2}{\kappa^2} + \frac{4 E(\kappa^2)}{\kappa^2 K(\kappa^2)} \right], \quad (16) $$

<figure><img src="image_3.png" /><figcaption>FIG. 2. The phase shift of the FFLO modulation for several values of the anisotropy parameter $\zeta$. The straight line corresponds to an unperturbed incommensurate structure.</figcaption></figure>

or since $q(\kappa^2) = \pi \Delta_0 / \kappa K(\kappa^2)$ and $h(\kappa^2) = 2E(\kappa^2)\Delta_0 / \pi \kappa$ is the parametric definition of the modulation wave vector and the critical field, respectively [45],

$$ \langle F_{sn}^{\varphi} \rangle = \frac{\delta q}{8} \sqrt{vq(\kappa^2)} - \frac{v}{16} \left[ 1 - \frac{2}{\kappa^2} + 2q(\kappa^2)h(\kappa^2) \right]. \quad (17) $$

The first term is proportional to the soliton density (inversely proportional to the distance between the solitons) and the second one describes the interaction with the field-induced potential.

One can see that as $\delta q \to 0$ the parameter $\zeta$ increases its value. In the limit $\zeta \to 1$ the period of the soliton lattice diverges,

$$ \lim_{\zeta \to 1} L = \frac{2\kappa}{\sqrt{v}} \ln \frac{4}{\sqrt{1-\kappa^2}}, \quad (18) $$

and the modulation wave vector is negligible. For $\zeta \ge 1$ there is no real valued solution of Eq. (15). Consequently, there exists a limiting value of $\delta q_c \equiv 2\sqrt{11t/\pi v_F}$ beyond which the phase is always locked, and the commensurate phase remains stable as $\delta q$ decreases further. That means that the modulation vector coincides with the modulation induced by the field. Therefore at $\delta q_c$ we have a lock-in of the FFLO modulation to the field-induced potential, i.e., an incommensurate/commensurate-type phase transition. That corresponds to $L \to \infty$ and the FFLO modulation is described by $\Delta_0 \cos Qx$. Figure 3 illustrates the behavior of the spatially varied order parameter in the vicinity of the IC-C phase transition when $\zeta \to 1$. One can see the domain walls (phase solitons) in the spatial distribution of the order parameter.

<figure><img src="image_4.png" /><figcaption>FIG. 3. The spatial distribution of the order parameter in the vicinity of the IC-C phase transition for $Q = 0.7$ (upper panel) and $Q = 1.7$ (lower panel).</figcaption></figure>
So far we have discussed the quasi-1D case. However, the theory developed is directly applicable to the quasi-2D materials as well. The orientation of the FFLO modulation vector is arbitrary in the pure Pauli-limited s-wave superconductor in the case of a Fermi surface with an elliptical cross section, which can be mapped by scaling transformation to the isotropic case [46]. However, the direction of the FFLO wave vector in real crystals is determined by the crystal field effects [47,48] and/or the pairing symmetry [49]. Therefore, when the in-plane magnetic field is perpendicular to the FFLO wave vector **q** we reduce the quasi-2D case to the above-described quasi-1D case, Eq. (3), but with the renormalized coefficients: $\beta \equiv \frac{\pi}{8}K_3(T)v_F^2$, $\delta \equiv -\frac{3\pi}{128}K_5(T)v_F^4$, $\mu \equiv -\frac{\pi}{4}K_5(T)v_F^2$, $\gamma \equiv \frac{\pi}{4}K_5(T)$. The coefficients $\gamma$, $v_F$ have the same values. The resonance conditions could be realized, for example, in $(TMTSF)_2ClO_4$ with in-plane perpendicular to **q** $H_{c2}^{b'}(0) = 3$ T, $d_z = 1.31$ nm, $\xi_0^a = 45$ nm that result in $\lambda_H = 2\pi/Q = 1050$ nm and $\lambda_{FFLO}^0 = 2\pi/q = \pi\hbar v_F/\Delta_0 = \pi^2\xi_0^a = 444$ nm, i.e., $\lambda_{FFLO}^0 < \lambda_H^0$ at $T=0$. Therefore it guarantees that while approaching a tricritical point the resonance $\lambda_{FFLO}(T) = \lambda_H(T)$ will occur. In the case of the quasi-2D magnetic-field-induced organic superconductor $\lambda$-(BETS)$_2FeCl_4$ the perpendicular $H_{c2}(0) = 4$ T, $d_z = 1.85$ nm, $\xi_0 = 6.5$ nm, and $\lambda_{FFLO}^0 = 65$ nm. The superconductivity appears at the in-plane field $H \gtrsim 18$ T and reaches the maximum critical temperature, $T_C = 4$ K, at $H \sim 32$ T. This means that $\lambda_H$ decreases from $\lambda_H = 556$ nm at $T=0$ to 312 nm at $T_C$. This again guarantees the fulfillment of the resonance condition in the FFLO phase.

In conclusion: Starting from the Eilenberger equations for a quasi-1D superconductor, which takes into account the orbital effects of an external field, in the proximity to the resonance between the period of the FFLO modulation and the period of the interlayer coupling modulation induced by the external field, we have reduced them to the continuum limit approximation of the Frenkel-Kontorova model equation. This equation is obtained by minimization of the Ginzburg-Landau functional (analogous to the extended Lawrence-Doniach theory for layered superconductors [27]) for an inhomogeneous distribution of the order parameter. As a main result we show that, besides the anomalous cusps in the temperature and angular dependencies of the in-plane critical field, when the FFLO modulation wave vector is in the resonance with the magnetic wave vector, predicted in our previous work, the following occurs: (i) The order parameter is defined by $\Delta_0 \cos[Qx + \varphi(x)]$, with $\varphi(x)$ describing regularly spaced phase solitons. (ii) The proximity to the resonance condition is characterized by the field-induced phase transition from the incommensurate to the commensurate structure (C-I) of the spatial modulation of the order parameter with respect to the magnetic wave vector structure. The periodic potential tends to “lock-in” the FFLO modulation into a commensurate configuration. (iii) In the incommensurate phase near the C-I transition the periodicity of the phase solitons or phase-domain walls depends continuously on the temperature and external field strength. This soliton lattice is a compromise between two periodicities of the system. The soliton represents a “phase slip” of the order parameter modulation. (iv) The distance between the phase slips logarithmically diverges at critical detuning or misfit, $\delta q_c$.

It would be promising to study anomalies of the superconducting properties accompanying the field-induced passage of the modulation wave vector of the order parameter through the commensurate value. The formation of the phase soliton lattice could be observed in the NMR spectroscopy measurements similar to those performed in Ref. [15]. The NMR signal, which is proportional to the local electronic spin susceptibility, may be very different in the phase domain wall as compared to the commensurate phase. The appearance of the domain walls in the FFLO phase may enhance its pinning and therefore could be revealed in transport measurements.

## ACKNOWLEDGMENTS

We acknowledge fruitful discussion with L. Bulaevskii, S. Uji, and J. Wosnitza. This work was supported by the French ANR SUPERTRONICS: Grant No. ANR-15-CE24-0008-02.

## APPENDIX: DERIVATION OF THE EXTENDED NONLINEAR LAWRENCE-DONIACH EQUATION

We start from the Eilenberger equations for a quasi-1D superconductor [50]:

$$ big[\Omega_n + \frac{1}{2}\hat{\Pi}(Q)]f_\omega(x, k_\perp, p_\perp) = \Delta(x, k_\perp)g_\omega(x, k_\perp, p_\perp), \quad (A1) $$

$$ [\Omega_n - \frac{1}{2}\hat{\Pi}(Q)]f_{\omega}^{\dagger}(x, k_{\perp}, p_{\perp}) = \Delta^{*}(x, k_{\perp})g_{\omega}(x, k_{\perp}, p_{\perp}), \quad (A2) $$

$$ \frac{1}{2}\hat{\Pi}(0)g_{\omega}(x, k_{\perp}, p_{\perp}) = \Delta^{*}(x, k_{\perp})f_{\omega}(x, k_{\perp}, p_{\perp}) - f_{\omega}^{\dagger}(x, k_{\perp}, p_{\perp})\Delta(x, k_{\perp}), \quad (A3) $$

where

$$ \hat{\Pi} \equiv \hbar v_F \frac{\partial}{\partial x} + 4it_y \sin(p_y d_y) \sin \left( Q_y x - \frac{k_y d_y}{2} \right) + 4it_z \sin(p_z d_z) \sin \left( Q_z x - \frac{k_z d_z}{2} \right), \quad (A4) $$

obtained for an arbitrary direction of the external magnetic field with the chosen gauge (see Sec. II) resulting in $A_z = -xH_y = -xH \sin\alpha \sin\theta$, $A_y = xH_z = xH \cos\alpha \sin\theta$, $A_x = 0$, where $\alpha$ accounts for the field direction, with amplitude $H$, in the $y-z$ plane from the positive $y$ axis, while $\theta$ is the angle accounting the field direction from the positive $x$ axis. This choice leads to $Q_y = \pi d_y H_z / \phi_0 = \pi d_y H \cos\alpha \sin\theta / \phi_0$, $Q_x = -\pi d_z H_y / \phi_0 = -\pi d_z H \sin\alpha \sin\theta / \phi_0$, with $\phi_0 = \pi \hbar c/e$. Here $k_\perp = (k_y, k_z)$, $p_\perp = (p_y, p_z)$, $v_F = v_{Fx}$ **i** is the Fermi velocity along the $x$ axis, $\Omega_n \equiv \omega_n + 1/2\tau - i\hbar\operatorname{sgn}(\omega_n)$ with $\hbar = \mu_B H$ is the Zeeman energy, and we have assumed that the vector potential varies slowly at the inter-chain distances. (This assumption means that we neglect the diamagnetic screening currents and take the magnetic field as uniform and given by the external field $H$.) These equations depend on $k_\perp$, which takes into account the dependence of the order parameter on the center-of-mass coordinate $(z+z')/2$ and $(y+y')/2$. The order parameter is defined self-consistently
as

$$
\frac{1}{\lambda} \Delta(x, k_{\perp}) = 2\pi T \text{Re} \sum_{\omega>0} \langle f_{\omega}(x, k_{\perp}, p_{\perp}) \rangle, \quad (A5)
$$

where $\lambda$ is the pairing constant and the brackets denote
averaging over $p_y$ and $p_z$ and $v_{F_x}$,

$$
\langle \dots \rangle \equiv N^{(1)} \int_{-\frac{\pi}{d_y}}^{\frac{\pi}{d_y}} \frac{d_y dp_y}{2\pi} \int_{-\frac{\pi}{d_z}}^{\frac{\pi}{d_z}} \frac{d_z dp_z}{2\pi} \frac{1}{2} \sum_{\pm v_F} (\dots), \quad (A6)
$$

with $N^{(1)} \equiv m/2\pi p_x d_y d_z$. The temperature unit is so that the Boltzmann constant $k_B = 1$. The maximum critical temperature corresponds to $k_y = 0$, and we have freedom in choosing $k_z$, for example, we can make $k_z = 0$. For the purpose of simplicity, choosing the external magnetic field in the x-y plane, finally we obtain

$$
\hat{\Pi} \equiv \hbar v_F \frac{\partial}{\partial x} + 4it_z \sin(p_z d_z) \sin(\chi_A), \quad (A7)
$$

where $\chi_A \equiv Q_x x - \frac{k_z}{2} d$.

In sequel we make two basic assumptions. First, it is
required that the system is near the phase transition that allows
expansion of Green functions in powers of Δ/ωn. However, the
assumption used is that not only the order parameter is small
but also the spatial gradients of the order parameter are small.
If we assume that the system resides near the phase transition,
then g<sub>ω</sub>(x,k<sub>⊥</sub>,p<sub>⊥</sub>) ≃ sgnωn, and this condition gives rise to

$$
(\Omega_n + \frac{1}{2} \hat{\Pi}) f_{\omega}(x, k_{\perp}, p_{\perp}) = \Delta, \quad (A8)
$$

where $\Omega_n = \omega_n - i\mu_B H$. It is convenient to write the procedure as the following iterative scheme:

$$
f_{\omega}^{(k+1)}(x, k_{\perp}, p_{\perp}) = \frac{\Delta}{\Omega_n} g_{\omega}^{(k)} - \frac{1}{2\Omega_n} \hat{\Pi} f_{\omega}^{(k)}, \quad (A9)
$$

$$
f_{\omega}^{\dagger(k+1)}(x, k_{\perp}, p_{\perp}) = \frac{\Delta^*}{\Omega_n} g_{\omega}^{(k)} + \frac{1}{2\Omega_n} \hat{\Pi} f_{\omega}^{\dagger(k)}. \quad (A10)
$$

The zero-order approximation (the absence of the su-
perconducting state) results in $g_{\omega}^{(0)}(x,k_{\perp},p_{\perp}) = 1$ and
$f_{\omega}^{(0)}(x,k_{\perp},p_{\perp}) = 0$. Then for the first-order approxima-
tion the iterative scheme produces $f_{\omega}^{(1)}(x,k_{\perp},p_{\perp}) = \Delta/\Omega_n$,
$f_{\omega}^{\dagger(1)}(x,k_{\perp},p_{\perp}) = \Delta^*/\Omega_n$. For the iterative scheme we need
the iterative expression for the $g_{\omega}^{(k)}$ function. If we write it, for
example, for $g_{\omega}^{(5)}$, as $g_{\omega}^{(5)} = \tilde{g}_{\omega}^{(0)} + \lambda \tilde{g}_{\omega}^{(1)} + \lambda_{\omega}^2 \tilde{g}_{\omega}^{(2)} + \lambda_{\omega}^3 \tilde{g}_{\omega}^{(3)} +$
$\lambda_{\omega}^4 \tilde{g}_{\omega}^{(4)} + \lambda_{\omega}^5 \tilde{g}_{\omega}^{(5)}$ and $f_{\omega}^{(5)} = \tilde{f}_{\omega}^{(0)} + \lambda \tilde{f}_{\omega}^{(1)} + \lambda^2 \tilde{f}_{\omega}^{(2)} + \lambda^3 \tilde{f}_{\omega}^{(3)} +$
$\lambda^4 \tilde{f}_{\omega}^{(4)} + \lambda^5 \tilde{f}_{\omega}^{(5)}$, then making use of the normalization condi-
tion,

$$
g_{\omega}^{2} + f_{\omega} f_{\omega}^{\dagger} = 1, \quad (A11)
$$

we obtain the following system of coupled expressions for the
kth iteration,

$$
\tilde{g}_{\omega}^{(0)}(x, k_{\perp}, p_{\perp}) = 1, \quad (A12)
$$

$$
\sum_{i=0}^{k} \tilde{g}_{\omega}^{(i)} \tilde{g}_{\omega}^{(k-i)} = - \sum_{i=0}^{k} \tilde{f}_{\omega}^{(k-i)} \tilde{f}_{\omega}^{(i)}. \quad (A13)
$$

The second iteration of the iteration scheme leads to

$$
f_{\omega}^{(2)}(x, k_{\perp}, p_{\perp}) = \frac{\Delta(x, k_{\perp})}{\Omega_n} - \frac{1}{2\Omega_n^2} \hat{\Pi} \Delta(x, k_{\perp}), \quad (A14)
$$

$$
f_{\omega}^{\dagger(2)}(x, k_{\perp}, p_{\perp}) = \frac{\Delta^{*}(x, k_{\perp})}{\Omega_{n}} + \frac{1}{2\Omega_{n}^{2}} \hat{\Pi} \Delta^{*}(x, k_{\perp}). \quad (A15)
$$

Therefore here $\tilde{f}_{\omega}^{(2)} = -\hat{\Pi}\Delta(x, k_{\perp})/2\Omega_{n}^{2}$, $\tilde{f}_{\omega}^{\dagger(2)} = \hat{\Pi}\Delta^{*}(x, k_{\perp})/2\Omega_{n}^{2}$. From the normalization condition we can obtain also

$$
\tilde{g}_{\omega}^{(1,2)} + 2\tilde{g}_{\omega}^{(0)}\tilde{g}_{\omega}^{(2)} = -\tilde{f}_{\omega}^{(2)}\tilde{f}_{\omega}^{\dagger(0)} - \tilde{f}_{\omega}^{(1)}\tilde{f}_{\omega}^{\dagger(1)} - \tilde{f}_{\omega}^{(0)}\tilde{f}_{\omega}^{\dagger(2)}, \quad (\text{A16})
$$

$$
2\tilde{g}_{\omega}^{(0)}\tilde{g}_{\omega}^{(2)} = -\tilde{f}_{\omega}^{(1)}\tilde{f}_{\omega}^{\dagger(1)}, \quad (\text{A17})
$$

which gives $\tilde{g}_{\omega}^{(2)} = -\tilde{f}_{\omega}^{(1)}\tilde{f}_{\omega}^{\dagger(1)}/2 = -f_{\omega}^{(1)}\tilde{f}_{\omega}^{\dagger(1)}/2$ and $\tilde{g}_{\omega}^{(1)} = 0$.
Therefore the normal Green function $g_{\omega}^{(2)}(x, k_{\perp}, p_{\perp})$ acquires
the following form:

$$
\begin{align*}
g_{\omega}^{(2)}(x, k_{\perp}, p_{\perp}) &= \tilde{g}_{\omega}^{(0)} + \tilde{g}_{\omega}^{(1)} + \tilde{g}_{\omega}^{(2)} = 1 - \frac{f_{\omega}^{(1)} f_{\omega}^{\dagger(1)}}{2} \\
&= 1 - \frac{|\Delta(x, k_{\perp})|^2}{2\Omega_n^2}.
\end{align*}
$$

Following similar steps, then performing averaging over $p_y$ and $p_z$ and Fermi surface $v_{F_x}$, we obtain on the fifth iteration for $f_{\omega}^{(5)} \equiv f_{\omega}^{(5)}(x, k_{\perp}, p_{\perp})$,

$$
\begin{align*}
f_{\omega}^{(5)} &= \frac{\Delta(x, k_{\perp})}{\Omega_n} - \frac{\Delta(x) |\Delta(x, k_{\perp})|^2}{2\Omega_n^3} \\
&\quad + \frac{1}{4\Omega_n^3} \hat{\Pi}^2 \Delta(x, k_{\perp}) + \frac{3\Delta(x, k_{\perp}) |\Delta(x, k_{\perp})|^4}{8\Omega_n^5} \\
&\quad + \frac{1}{16\Omega_n^5} \hat{\Pi}^4 \Delta(x, k_{\perp}) - \frac{4}{8\Omega_n^5} |\Delta(x, k_{\perp})|^2 \hat{\Pi}^2 \Delta(x, k_{\perp}) \\
&\quad - \frac{1}{8\Omega_n^5} \Delta^2(x, k_{\perp}) \hat{\Pi}^2 \Delta^*(x, k_{\perp}) \\
&\quad - \frac{3}{8\Omega_n^5} \Delta^*(x, k_{\perp}) \hat{\Pi} \Delta(x, k_{\perp}) \hat{\Pi} \Delta(x, k_{\perp}) \\
&\quad - \frac{2}{8\Omega_n^5} \Delta(x, k_{\perp}) \hat{\Pi} \Delta(x, k_{\perp}) \hat{\Pi} \Delta^*(x, k_{\perp}).
\end{align*}
$$

Substituting the found averages into the self-consistency re-
lation given by Δ<sup>(5)</sup>(x, k<sub>⊥</sub>)/λ = πT Σ<sub>n</sub> ⟨f<sub>ω</sub><sup>(5)</sup>(x, k<sub>⊥</sub>, p<sub>⊥</sub>)⟩ and
using the standard regularization rule

$$
\frac{1}{\lambda} = \ln \frac{T}{T_{c0}} + 2\pi T \sum_{\omega_n > 0} \frac{1}{\omega_n}, \quad (A18)
$$

we obtain the following extended version of the Ginzburg-
Landau equation (for clarity we dropped down k⊥ depen-
dence), if by $\hat{\Pi}$ one understands its three-dimensional version,

$$
\begin{align*}
0 = \Delta(x) \Biggl\{ & \pi [K_1(T) - K_1^0(T)] - \ln \frac{T}{T_{c0}} \Biggr\} \\
& + \frac{\pi K_3(T)}{4} \hat{\Pi}^2 \Delta(x) - \frac{\pi K_3(T)}{2} \Delta(x) |\Delta(x)|^2 \\
& + \frac{3\pi K_5(T)}{8} \Delta(x) |\Delta(x)|^4 + \frac{\pi K_5(T)}{16} \hat{\Pi}^4 \Delta(x)
\end{align*}
$$
$$ -\frac{\pi K_5(T)}{8} \{4|\Delta(x)|^2 \hat{\Pi}^2 \Delta(x) + \Delta^2(x) \hat{\Pi}^2 \Delta^*(x) \\ + 3\Delta^*(x) \hat{\Pi} \Delta(x) \hat{\Pi} \Delta(x) + 2\Delta(x) \hat{\Pi} \Delta(x) \hat{\Pi} \Delta^*(x)\}, \quad (A19) $$

where we have introduced the notations [39] $K_m(T) = 2T \sum_{\omega>0} (\omega_n - i\mu_B H)^{-m}$ and $K_m^0(T) = 2T \sum_{\omega>0} \omega_n^{-m}$. This result is similar to the one obtained by Houzet and Mineev [5]. Now we recall the original notations of Eq. (A7) and after the averaging procedure we can write (neglecting the $k_⊥$ dependence)

$$ \begin{aligned} \alpha \Delta(x) = & \beta \partial_x^2 \Delta(x) - \delta \partial_x^4 \Delta(x) - 2\gamma \Delta(x) |\Delta(x)|^2 \\ & - 3\nu \Delta(x) |\Delta(x)|^4 + \mu |\Delta(x)|^2 \partial_x^2 \Delta(x) \\ & + \frac{\mu}{4} \Delta^2(x) \partial_x^2 \Delta^*(x) + \frac{3\mu}{4} \Delta^*(x) [\partial_x \Delta(x)]^2 \\ & + \frac{\mu}{2} \Delta(x) \partial_x \Delta(x) \partial_x \Delta^*(x) - \frac{t^2 \mu}{2} \\ & \times [1 - 7 \cos(2\chi_A)][\partial_x \chi_A]^2 \Delta(x) - 80t^2 \nu \sin^2(\chi_A) \\ & \times \Delta(x) |\Delta(x)|^2 + 2t^2 \mu [3 \sin(2\chi_A) \partial_x \chi_A \partial_x \Delta(x) \\ & + \sin(2\chi_A) \partial_x^2 \chi_A \Delta(x) + 3 \sin^2(\chi_A) \partial_x^2 \Delta(x)] \\ & - 8t^2 \gamma \sin^2(\chi_A) \Delta(x) - 48t^4 \nu \sin^4(\chi_A) \Delta(x), \end{aligned} $$

where $\alpha \equiv \ln \frac{T}{T_{c0}} - \pi [K_1(T) - K_1^0(T)]$, $\beta \equiv \frac{\pi}{4} K_3(T) \langle v_{F_x}^2 \rangle$, $\delta \equiv -\frac{\pi}{16} K_5(T) \langle v_{F_x}^4 \rangle$, $\gamma \equiv \frac{\pi}{4} K_3(T)$, $\mu \equiv -\frac{\pi}{2} K_5(T) \langle v_{F_x}^2 \rangle$,

$\nu \equiv -\frac{\pi}{8} K_5(T)$. Here, the derivative signs mean the following: $\partial^2 \Delta(x) = \langle v_{F_x}^2 \rangle \partial_x^2 \Delta(x)$, $\partial^4 \Delta(x) = \langle v_{F_x}^4 \rangle \partial_x^4 \Delta(x)$, $\partial \chi_A \partial \Delta(x) = \langle v_{F_x}^2 \rangle \partial_x \chi_A \partial_x \Delta(x)$, and $[\partial \chi_A]^2 = \langle v_{F_x}^2 \rangle [\partial_x \chi_A]^2$. If one neglects the interchain hopping, then the previous equation reduces to

$$ \begin{aligned} 0 = & \alpha \Delta(x) - \beta \partial_x^2 \Delta(x) + \delta \partial_x^4 \Delta(x) \\ & + 2\gamma \Delta(x) |\Delta(x)|^2 + 3\nu \Delta(x) |\Delta(x)|^4 \\ & - \mu |\Delta(x)|^2 \partial_x^2 \Delta(x) - \frac{\mu}{4} \Delta^2(x) \partial_x^2 \Delta^*(x) \\ & - \frac{3\mu}{4} \Delta^*(x) [\partial_x \Delta(x)]^2 - \frac{\mu}{2} \Delta(x) \partial_x \Delta(x) \partial_x \Delta^*(x). \end{aligned} \quad (A20) $$

As shown in Ref. [39], a solution of this equation is the Jacobi elliptic sine function $\Delta(x) = \frac{v_F}{\xi(k)} \text{ksn}(\frac{x}{\xi(k)}, k)$, expressed in terms of the modulus $k$, which is determined upon minimizing the free energy at the fixed external field parameter $h$. The linearized version of this equation,

$$ \begin{aligned} \alpha \Delta(x) = & \beta \partial_x^2 \Delta(x) - \delta \partial_x^4 \Delta(x) - 48t^4 \nu \sin^4(\chi_A) \Delta(x) \\ & - \frac{t^2}{2} \mu [1 - 7 \cos(2\chi_A)][\partial_x \chi_A]^2 \Delta(x) \\ & + 2t^2 \mu [3 \sin(2\chi_A) \partial_x \chi_A \partial_x \Delta(x) \\ & + 3 \sin^2(\chi_A) \partial_x^2 \Delta(x)] - 8t^2 \gamma \sin^2(\chi_A) \Delta(x), \end{aligned} \quad (A21) $$

is used in the main part of the article.

[1] G. Zwicknagl and J. Wosnitza, in *BCS: 50 Years*, edited by L. N. Cooper and D. Feldman (World Scientific Publishing, Singapore, 2011).

[2] A. I. Buzdin, *Rev. Mod. Phys.* **77**, 935 (2005).

[3] L. G. Aslamazov, Zh. Eksp. Teor. Phys. **55**, 1477 (1968) [Sov. Phys. JETP **28**, 773 (1969)].

[4] S. Takada, *Prog. Theor. Phys.* **43**, 27 (1970).

[5] M. Houzet and V. P. Mineev, *Phys. Rev. B* **74**, 144522 (2006).

[6] Y. Matsuda and H. Shimahara, *J. Phys. Soc. Jpn.* **76**, 051005 (2007).

[7] H. Shimahara, *J. Phys. Soc. Jpn.* **67**, 736 (1998); **68**, 3069 (1998); H. Shimahara and K. Moriwake, *ibid.* **71**, 1234 (2002).

[8] A. I. Buzdin and L. N. Bulaevskii, *Sov. Phys. Usp.* **27**, 830 (1984).

[9] A. I. Larkin and Yu. N. Ovchinnikov, Zh. Eksp. Teor. Phys. **47**, 1136 (1964) [Sov. Phys. JETP **20**, 762 (1965)].

[10] P. Fulde and R. A. Ferrell, *Phys. Rev.* **135**, A550 (1964).

[11] M. A. Tanatar, T. Ishiguro, H. Tanaka, and H. Kobayashi, *Phys. Rev. B* **66**, 134503 (2002).

[12] W. A. Coniglio, L. E. Winter, K. Cho, C. C. Agosta, B. Fravel, and L. K. Montgomery, *Phys. Rev. B* **83**, 224507 (2011).

[13] R. Lortz, Y. Wang, A. Demuer, P. H. M. Böttger, B. Bergk, G. Zwicknagl, Y. Nakazawa, and J. Wosnitza, *Phys. Rev. Lett.* **99**, 187002 (2007).

[14] B. Bergk, A. Demuer, I. Sheikin, Y. Wang, J. Wosnitza, Y. Nakazawa, and R. Lortz, *Phys. Rev. B* **83**, 064506 (2011).

[15] J. A. Wright, E. Green, P. Kuhns, A. Reyes, J. Brooks, J. Schlueter, R. Kato, H. Yamamoto, M. Kobayashi, and S. E. Brown, *Phys. Rev. Lett.* **107**, 087002 (2011).

[16] H. Mayaffre, S. Krämer, M. Horváti, C. Berthier, K. Miyagawa, K. Kanoda, and V. F. Mitrović, *Nat. Phys.* **10**, 928 (2014).

[17] A. I. Buzdin and V. V. Tugushev, Zh. Eksp. Teor. Fiz. **85**, 735 (1983) [Sov. Phys. JETP **58**, 428 (1983)].

[18] A. B. Vorontsov, J. A. Sauls, and M. J. Graf, *Phys. Rev. B* **72**, 184501 (2005).

[19] C. C. Agosta, J. Jin, W. A. Coniglio, B. E. Smith, K. Cho, I. Stroe, C. Martin, S. W. Tozer, T. P. Murphy, E. C. Palm, J. A. Schlueter, and M. Kurmoo, *Phys. Rev. B* **85**, 214514 (2012).

[20] G. Koutroulakis, H. Kühne, J. A. Schlueter, J. Wosnitza, and S. E. Brown, *Phys. Rev. Lett.* **116**, 067003 (2016).

[21] R. Beyer, B. Bergk, S. Yasin, J. A. Schlueter, and J. Wosnitza, *Phys. Rev. Lett.* **109**, 027003 (2012).

[22] S. Uji, K. Kodama, K. Sugii, T. Terashima, Y. Takahide, N. Kurita, S. Tsuchiya, M. Kimata, A. Kobayashi, B. Zhou, and H. Kobayashi, *Phys. Rev. B* **85**, 174530 (2012).

[23] K. Cho, H. Kim, M. A. Tanatar, Y. J. Song, Y. S. Kwon, W. A. Coniglio, C. C. Agosta, A. Gurevich, and R. Prozorov, *Phys. Rev. B* **83**, 060502(R) (2011).

[24] D. A. Zocco, K. Grube, F. Eilers, T. Wolf, and H. v. Löhneysen, *Phys. Rev. Lett.* **111**, 057007 (2013).
[25] S. Yonezawa, S. Kusaba, Y. Maeno, P. Auban-Senzier, C. Pasquier, K. Bechgaard, and D. Jérome, Phys. Rev. Lett. 100, 117002 (2008).

[26] M. D. Croitoru, M. Houzet, and A. I. Buzdin, Phys. Rev. Lett. 108, 207005 (2012); M. D. Croitoru and A. I. Buzdin, Phys. Rev. B 86, 224508 (2012); J. Supercond. Nov. Magn. 26, 1657 (2013).

[27] M. D. Croitoru and A. I. Buzdin, Phys. Rev. B 86, 064507 (2012); J. Phys.: Condens. Matter 25, 125702 (2013).

[28] M. D. Croitoru and A. I. Buzdin, Phys. Rev. B 89, 224506 (2014); J. Supercond. Nov. Magn. 28, 1305 (2014).

[29] P. Bak, Rep. Prog. Phys. 45, 587 (1982).

[30] T. A. Kontorova and Ya. I. Frenkel, Zh. Eksp. Teor. Fiz. 8, 89 (1938); 8, 1340 (1938); Ya. I. Frenkel, Introduction Into the Theory of Metals (Nauka, Leningrad, 1972) (in Russian).

[31] L. M. Floria and J. J. Mazo, Adv. Phys. 45, 505 (1996).

[32] J. B. Sokoloff, Phys. Rev. B 16, 3367 (1977).

[33] O. M. Braun and Y. S. Kivshar, The Frenkel-Kontorova Model, Concept, Method and Application (Springer-Verlag, Berlin, 2004); Phys. Rep. 306, 1 (1998).

[34] M. Lucci, D. Badoni, V. Merlo, I. Ottaviani, G. Salina, and M. Cirillo, A. V. Ustinov, and D. Winkler, Phys. Rev. Lett. 115, 107002 (2015).

[35] R. Besseling, R. Niggebrugge, and P. H. Kes, Phys. Rev. Lett. 82, 3144 (1999).

[36] W. Dieterich, P. Fulde, and I. Peschel, Adv. Phys. 29, 527 (1980).

[37] T. Tsuzuki, J. Low Temp. Phys. 9, 525 (1972).

[38] L. A. Turkevich and R. A. Klemm, Phys. Rev. B 19, 2520 (1979).

[39] A. I. Buzdin and H. Kachkachi, Phys. Lett. A 225, 341 (1997).

[40] A. I. Buzdin and S. V. Polonskii, JETP 66, 422 (1987) [Zh. Eksp. Teor. Fiz. 93, 747 (1987)].

[41] E. Dzyaloshinskii, Zh. Eksp. Teor. Fiz. 47, 992 (1964) [Sov. Phys. JETP 20, 665 (1964)].

[42] Yu. A. Izyumov, Sov. Phys. Usp. 27, 845 (1984).

[43] Yu. A. Izyumov, in Neutron Diffraction in Long-Periodic Structures (Energoatomizdat, Moscow, 1987).

[44] F. C. Frank and J. H. Van der Merwe, Proc. Roy. Soc. (London) A 198, 205 (1949); 198, 216 (1949).

[45] K. Machida and H. Nakanishi, Phys. Rev. B 30, 122 (1984).

[46] J. P. Brison et al., Phys. C (Amsterdam, Neth.) 250, 128 (1995).

[47] D. Denisov, A. Buzdin, and H. Shimahara, Phys. Rev. B 79, 064506 (2009).

[48] A. Buzdin, Y. Matsuda, and T. Shibauchi, Europhys. Lett. 80, 67004 (2007).

[49] H. Maki and H. Won, Czech. J. Phys. B 46, 1035 (1996).

[50] N. B. Kopnin, Theory of Nonequilibrium Superconductivity (Clarendon Press, Oxford, 2001).