Vortex jamming in superconductors and granular rheology

We demonstrate that a highly frustrated anisotropic Josephson junction array(JJA) on a square lattice exhibits a zero-temperature jamming transition, which shares much in common with those in granular systems. Anisotropy of the Josephson couplings along the horizontal and vertical directions plays roles similar to normal load or density in granular systems. We studied numerically static and dynamic response of the system against shear, i. e. injection of external electric current at zero temperature. Current-voltage curves at various strength of the anisotropy exhibit universal scaling features around the jamming point much as do the flow curves in granular rheology, shear-stress vs shear-rate. It turns out that at zero temperature the jamming transition occurs right at the isotropic coupling and anisotropic JJA behaves as an exotic fragile vortex matter : it behaves as superconductor (vortex glass) into one direction while normal conductor (vortex liquid) into the other direction even at zero temperature. Furthermore we find a variant of the theoretical model for the anisotropic JJA quantitatively reproduces universal master flow-curves of the granular systems. Our results suggest an unexpected common paradigm stretching over seemingly unrelated fields - the rheology of soft materials and superconductivity.

Physics continues to thrive on analogy [1]. Rheological properties of matters [2] and electric transport properties of superconductors [3] exhibit intriguing analogies. The flow curves in rheology, the shear-stress vs the shear-rate, correspond to the current-voltage curves in superconductors [4]. Exploring this analogy further we here demonstrate by computer simulations that there exist the zero-temperature jamming transitions and glassy non-linear rheology, originally found in granular and other materials [5,6,7,8,9,10,11,12,13,14,15,16], in a class of highly frustrated anisotropic Josephson junctions arrays (JJA) on a square lattice. Our key observation is that anisotropy of Josephson coupling plays the role of normal load or density in granular system such that a jamming transition takes place in the limit of isotropic Josephson coupling at zero temperature. Combined with accumulating evidences that the (vortex) liquid-glass transition occurs at zero temperature in isotropic JJA, our result provides a strong evidence that the isotropic coupling point at zero temperature is an ideal example of the so called J-point (Jamming point) [6] and that the anisotropic JJA is a promising system which allows explorations of both (athermal) unjamming-jamming transition and (thermal) liquid-glass transitions in a unified manner in a single system as originally proposed in the context of granular and glassy materials [6]. Furthermore, we show that a variant of the original JJA model emphasizing elastic nature in a particular direction can quantitatively reproduce scaling features of granular jamming transitions observed near the J-point.
The JJA [3,17,18] is a network of superconducting islands as depicted in Fig. 1 (a). The phase of superconducting order parameter of the islands at site i, θ i , is coupled to its nearest neighbors by the Josephson junctions. External transverse magnetic field B thread the cells in the forms of flux lines each of which carrying a flux quanta φ 0 . On average each unit cell of the square lattice carries f = Ba 2 /φ 0 flux lines. A flux line threading a unit cell induces a vortex of the phases around the cell much as a dislocation in a crystal. Here A ij is the vector potential. The static and dynamic properties of the JJA under transverse magnetic field are known to be described to a good accuracy by the energy term associated with the Josephson couplings [3,17], where K ij is the strength of the Josephson coupling. In mid 80's a tantalizing possibility of a superconducting glass in the JJA has been raised by Halsey [19]: if the transverse magnetic field is tuned in such a way that the number density f of the vorticies takes irrational values, a glassy state may be realized at low temperatures because the vorticies may not be able to form periodic structures called vortex lattices which are analogous to ordered structures of dislocations in Frank-Kasper phases. It has been argued that the frustration due to the gauge field like A ij in Eq. (1) mimics geometric frustration in structural glasses (see [20] for a review on the perspective on the frustration-based point of view on the glass transition). Indeed equilibrium relaxations were similar to the primary relaxation observed in typical fragile supercooled liquids [21,22]. Now recent studies appear to convincingly suggest that the putative (thermal) liquid-glass transition is actually taking place only at zero temperature exhibiting diverging length scale(s) [23,24,25]. We note that all these observations are made on systems in which the Josephson couplings K ij are isotropic. As stated at the beginning our key observation is that the anisotropy of the Josephson coupling K ij is a relevant variable and plays the role similar to the particle normal load or density in the jamming of granular materials [7] ( Fig. 1 (a)(b)). We parametrize the anisotropic coupling as, Such an anisotropic JJA can be realized experimentally by controlling the width of the junctions [26]. Let us also note that similar anisotropic couplings arise naturally also in cuprate high-Tc superconductors [27] and charge-density-wave systems [28].
Another key observation is that the model Eq. (1) can be slightly modified to build an effective Eulerian model for rheology of granular systems under horizontal shear as shown in Fig. 1 (b). We assume that particles move predominantly parallel to the direction of driving (x axis) confined within the horizontal layers without passing over others in the same layer. Then the variable θ i can be interpreted as a phase variable [4] by which the number density of particles within the i-th cell is described, for instance, as ρ i = (1/2)(1 + cos(θ i )). Here the size a of the cell corresponds to the typical scale of particles and their mutual distances [46]. The sinusoidal intra-layer couplings E ij (ψ) = − cos(ψ) in Eq. (1) are replaced by elastic couplings E ij (ψ) = ψ 2 /2 while the sinusoidal form is kept for the inter-layer couplings to allow phase slips between different layers. Let us call such a model as semi-elastic model. Here the vorticies represent dislocations. We assume that the geometrical frustration induced by the gauge field A ij mimics real frustrations in granular and other glassy systems [20]. In the absence of the frustration (A ij = 0) the semi-elastic model exhibits a non-linear rheology associated with a Kosterlitz-Thouless transition at finite temperatures [4].
The dynamics of the models can be described by the equation of motion, Here the frictional force F i is given by F i = −γ j (v i − v j ) with the summation taken over the 4 nearest neighbours of i. This equation of motion is nothing but the standard resistively and capacitively shunted junction (RCSJ) dynamics [3,24], which can also be viewed as a model for rheology of the layered systems [4]. For simplicity we choose mass (capacitance) m = 1, damping constant (resistance) γ = 1. Here we focus only on the zero temperature dynamics. Appropriate thermal noise can be added to Eq. (3) for finite temperature dynamics. For the semi-elastic model we assumed two different types of constitutive relations for the frictional force F i in Eq. (3): (1) Newtonian viscous (2) Bagnold's friction [38,39] where the sum is taken over the nearest neighbours on the two adjacent layers. The anisotropic coupling Eq. (2) can be motivated by recalling a well known problem in the science of friction. A class of friction models related to the Frenkel-Kontorova model [29] is known to exhibit the so called Aubry's transition [30] which is a kind of jamming transition at zero temperature. In Fig. 1 (d) we display a friction model proposed by Matsukawa and Fukuyama [31] which consists of two layers of atoms representing surfaces of two different solids. In general the ratio (winding number) f = l/a of the mean atomic spacings on the two different materials takes irrational values. The atoms in the same layers are connected to each other by springs while those on different layers interact with each other via short-ranged interactions of strength λ, which mimic the normal load. In the weak coupling regime λ < λ c , the two chains of atoms slide smoothly with respect to each other thanks to the incommensurability. On the other hand in the strong coupling regime λ > λ c the system is pinned into amorphous metastable states and a finite static frictional force or yield stress emerges. Note that an anisotropic JJ ladder [32] shown in Fig. 1 (c) which consist of two horizontal layers can be viewed as an Eulerian formulation of the friction model. The irrational winding number f can be identified with the irrational number density f of vorticies in a unit cell in the JJ ladder.
An important consequence of the anisotropic coupling Eq. (2) in 2 (and higher) dimensions is that the effective repulsive long-ranged interactions between the vorticies become anisotropic. For λ < 1 the vorticies will tend to align vertically since the repulsive force is stronger along the x axis, which make it much harder for the vorticies to move along the y axis, i.e. the direction with weaker coupling. Of course the situation becomes reversed for λ > 1.
The rigidity of the system can be probed by applying an external current just as external shear stress is applied on a solid. What corresponds to the shear stress Σ xy along the x axis ( Fig. 1 (b)) is the vertical external electric current J y ( Fig. 1 (a)) [4]. Then vorticies (dislocations) are driven along x axis by the Lorentz force. The resultant electric field E y which is proportional to the average velocity of the vorticies corresponds to the shear-rateΓ xy which measures the rate of plastic deformations in rheology. If the vorticies don't move significantly resisting against the Lorentz force, the energy dissipation is negligible and the system remains macroscopically superconducting. In practice, we apply shear to the system by forcing the top and bottom layers (walls) to move along the opposite directions at constant velocities. We measure the resultant electric field E y (shear-rateΓ xy ) defined as the slope of phase velocity v i developed in the system along the y axis. Electric current flowing through a junction from site i to j is defined as sin(θ i − θ j − A ij ). The currents running through the junctions parallel to x and y axes correspond to the shear stresses Σ xx and Σ xy respectively in rheology.
By construction of the system, static and dynamic resposes to J x at anisotropy λ at T = 0 is just the same as those to J y with λ ′ = 1/λ. Thus in the following we only display results of resposes to J y .
We numerically solved the equation of motion Eq. (3) by the 4th order Runge-Kutta method [47]. Periodic boundary condition is imposed along the x axis only. For a given irrational vortex density f we used its rational approximations p/q with integer p and q in systems of sizes L = nq (with n = 1 or 2). To explore larger length/time scales we use systematically better approximants to prevent commensurability (or matching) effects. (See APPENDIX A) Before starting measurements, we checked that the velocity profile becomes linear into the y axis without shear-bands and that observed quantities do not depend on the prior shear histories.
Shown in Fig. 2 are snapshots of the vorticies and local currents under shear of the anisotropic JJA. Chains of electric currents reminiscent of "force chains" [5,7] in granular materials can be noticed. The configurations at the isotropic point λ = 1 appear to manifest the diagonal stripe structures found in the ground states [33,34,35]. As expected the vorticies and the electric currents flowing along the trains of vorticies tend to align into the direction with weaker coupling. This observation strongly suggests that a jamming transition takes place at the isotropic point λ c = 1. The system behaves as a fluid (unjammed phase) for λ < λ c and amorphous solid (jammed phase) for λ > λ c with respect to J y as depicted in Fig. 1 (e). Furthermore it is interesting to note that the jammed state is inevitably fragile in somewhat similar sense as proposed in the context of granular matters [5]: the system with a given λ can resist against shear only into one direction (i. e. superconducting). Thus we may call such a state of matter as fragile vortex matter in the same spirit of [5]. The qualitative features are essentially the same in the semi-elastic model except that vorticies move only into x axis in the latter model.
The current-voltage curves obtained at different values of the coupling λ are displayed in Fig. 3 (a). At stronger coupling λ > 1 it appears that a non-zero critical current J c (λ) = lim E→0 J(E, λ) exists, which becomes larger with increasing λ. This means that the Lorentz force does not drive the vorticies significantly so that the system remains macroscopically superconducting along the y axis at strong enough coupling λ. The finite critical current corresponds to the yield stress Σ c in rheology. The disordered configurations of vorticies shown in Fig. 2 suggests that the system is an amorphous glassy state of vorticies. On the other hand, at smaller coupling λ < 1 and low enough E the Ohm's law J = σ(λ)E holds with finite linear conductivity σ(λ) which becomes larger with increasing λ. Thus the vorticies can flow easily producing significant energy dissipation at weak enough coupling λ. At the isotropic point λ = 1, we find a power law J ∝ E 1−α with 1 − α = 0.34 (3). This corresponds to the so called shear-thinning behaviour (α > 0) in rheology [2].
The above results strongly indicate that λ c = 1 is the critical point of a 2nd order phase transition at zero temperature. This is supported by a good scaling collapse of the data onto a master curve as shown in Fig. 3

(d).
Our scaling ansatz is similar in spirit to the ones used for the usual normal-to-superconducting phase transition at finite temperatures [36,37] which can also be reinterpreted in the context of rheology [4], The scaling function (master flow curve) is expected to behave asymptotically asJ(x) ∝ x for small enough x in the Ohmic phase (λ < λ c ) and lim x→0J (x) → const in the superconducting phase (λ > λ c ). The scaling ansatz Eq. (4) implies 1) the linear conductivity diverges as σ(λ) ∝ (λ − λ c ) −(∆−β) for λ → λ − c , 2) the critical current vanishes as J c (λ) ∝ (λ − λ c ) β for λ → λ + c and 3) the critical behaviourJ(x) ∝ x 1−α = x β/∆ sets-in for large x. Here we used the notations reflecting the analogy with the equilibrium critical behaviour of ferro-magnets under magnetic field as noticed by Wolf, Gubser and Imry [36]: the shear plays the role of symmetry breaking field like the magnetic field and the critical current emerges as an order parameter like the magnetization (see [11] for a similar argument in the context of rheology). As shown in Fig. 3, we find our scaling ansatz works well with λ c = 1. We have checked that the universality does not depend on the use of different irrational values of f as demonstrated in Fig. 3 (d). Let us note that the critial current discussed above J c (λ) = lim E→0 J(E, λ) corresponds to the dynamical yield stress in rheology. On the other hand one can also define the quasi-static critical current(s) needed to move out of a generic metastable state (or the ground state), corresponding to the static yield stress in rheology, as studied by Teitel and Jayaprakash [17] on the JJA. In practice one can consider athermal quasi-static processes similar to those used in some recent studies on amorphous solids [40,41,42]: starting from a metastable state reached from a random initial configuration by a qunech, i. e. deterministic energy descent process, the system is subjected to externally induced uniform strain δ (See Eq. (B1)) which is increased step by step. The system relaxes down to an energy minimum by the energy descent process after each small increment of δ. As the result one finds that the current J(δ) is a sawtooth-like function of δ: piecewise linear lines corresponding to elastic deformations broken by yield points at plastic events, as observed in amorphous solids [40,41,42]. A critical current can be defined as the value of the current just before reaching a yield point. In [19] Halsey has found that typical value of such quasi-static critical current J static c (1) is finite in the isotropic system λ = 1. Moreover we found that J static c (λ) varies smoothly with λ and remains finite even in the unjammed phase λ < 1 found above [43]. Our data of J(E, λ) in the flow curves becomes smaller than J static c (λ) meaning that the quasi-static current needed to move out of a generic metastable state and dymamic critical currents are distinct in the present system. Quite interestingly we observed that the flow curves become strongly dependent on strain histories below J static c (λ). In practice we had to use an annealing procedure to obtain stationary data: decrease E (shear rate) very slowly down to the target one. Slower shear rates are needed to investigate smaller current J (shear stress) regions. Note that no annealing is performed in the athermal quasi-static process discussed above. These observations suggest ruggedness of the energy landscape of the present system.
We also investigated static response to shear (See APPENDIX B for the details.). As shown in Fig. 4, the static helicity (shear) modulus is very sensitive to the anisotropy. The figure shows that the helicity modulus remains zero down to T → 0 for λ < 1 and becomes finite for λ > 1. Thus at T = 0 the isotropic point appears to be the critical point λ c = 1 being consistent with the dynamic response discussed so far. Here let us recall again that just by symmetry, helicity modulus Y x /λ at anisotropy λ is identical to Y y /λ ′ at λ ′ = 1/λ. Then the fact that λ c = 1 means that this system is quite exotic: fragile vortex matter which behaves as a solid with respect to shear along one direction but liquid for the other direction.
Another remarkable feature is the difference of the static helicity modulus between T = 0 and T → 0 limit which can be seen in Fig. 4 b). The helicity modulus at T = 0 only reflects local stability of energy minima. This difference suggests the existence of certain softmodes with vanishingly small energy gaps as discussed in APPENDIX B. This observation suggests that the vortex liquid behaviour at T = 0 is realized by some non-trivial softmodes due to frustrations.
Based on the above results we obtain the jamming phase diagram of the anisotropic JJA as shown in Fig. 1 (e) which is surprisingly similar to that of granular systems shown in Fig. 1 (f) [6,7,8]. Now let us turn to the semi-elastic model. In Fig. 3 (b,c), we display the flow curves of the semi-elastic model. The shear-stress due to the inter-layer coupling terms also obeys the Newtonian or Bagnold scaling for small λ at low enoughΓ. Flow curves at large λ suggests existence of non-zero yield stresses Σ c (λ) = limΓ →0 Σ(Γ, λ). We find again a power law behaviour Σ ∝ (Γ) 1−α with α > 0 (shear-thinning) at λ = 1. Indeed the scaling ansatz Eq. (4) (with E →Γ and J → Σ) works well again assuming λ c = 1 and the universality does not depend on the different irrational values of f (Fig. 3 (e,f)).
Recent numerical simulations of granular materials with/without strong dissipation at the particle level (coefficient of restitution smaller than/equal to 1) have found Bagnold/Newtonian scalings in the fluid phase and different critical exponents [10,12,13,14,44]. Quite interestingly the values of the shear-thinning exponent 1 − α found in our semi-elastic model with Bagnold/Newtonian frictions are 0.63 and 0.42 respectively in agreement with the exponents of the corresponding two-dimensional granular systems with repulsive linear spring forces between the particles. For a comparison master flow curves of the granular systems [44] are displayed in Fig. 3 (e,f). Quite remarkably the functional forms of the master flow-curves themselves agree very well.
In the present paper we focused on the responses of the systems to shear at zero temperature. Recent studies at the isotropic point λ = 1 at finite temperatures suggest critical behaviour with T → 0, i. e. T c (λ = 1) = 0 with diverging length scales [23,24,25]. If this would be confirmed, our system would provide a fascinating example where both the (thermal) liquid-glass transition and the (athermal) unjamming-jamming transition take place at the same thermodynamic point, T = 0 and λ c = 1, demonstrating deep connection between the two transitions. At least in this system, the jamming and glass transitions appears to be the two sides of a coin. It would be interesting to further explore this connection to make this statement more substantial. A related interesting question is whether or not the jamming (glass) phase survives at finite temperatures. Then the possibilities are 1) the jamming point at λ = 1 at T = 0 is an isolated critical point or 2) a critical line T c (λ) starts from the jamming point as shown in Fig. 1 (e). Our preliminary study points to the latter possibility.
We emphasize that the jamming transition here is purely due to geometrical frustration which is free from any quenched disorder in sharp contrast to the conventional vortex glasses [37] for which presence of random pinning centers are crucial. This is a much awaited, concrete example of a jamming-glass transition purely due to geometrical frustration [20]. It is tempting to speculate that similar phenomena may exist in frustrated magnets such as antiferromagnets on triangular, kagome and pyrochlore lattices [45]. We also note that it will be important and interesting to clarify how quenched disorders, which may not be completely avoided in experimental JJAs, come into play. by Teitel and Jayaprakash [17] (see also [33]). The latter behaviour suggests that a periodic vortex lattice (crystal) [17,18,33,34] associated with a given rational vortex density f n is formed and that the latter dictates physical properties of the system at length/time scales larger than its lattice spacing. Thus we expect the so called 'Bingham fluid' behaviour (fluid with finite yield stress) cannot be avoided for any small λ for fixed n and that the genuine fluid phase at T = 0 is realized only for truly irrational vortex densities. On the other hand the above results shown in Fig. 5 (b) suggest that physical properties at short enough length/time scales, which are n independent, reflect those properties of irrational f . Thus our strategy in the present work is to choose large enough n such that the n dependency do not become relvant within the range of E (shear rateΓ) we choose to work on. To study static response to shear along, say y axis, it is useful to consider a modified Hamiltonian, with periodic boundary conditions along both x and y axes. Here e x and e y are unit vectors along x and y axes. This is equivalent to consider a system with twisted boundary condition (shear) with total phase difference Lδ forced across the system along the y axis [17]. In equilbrium at temperature T the free-energy F (δ) of the system under shear strain δ can be defined. Then the static helicity modulus Y y is defined as, −λβ <i,j>||ey <k,l>||ey ( sin(φ ij ) sin(φ kl ) − sin(φ ij ) sin(φ kl ) ) (B2) where β ≡ 1/(k B T ) and . . . are equilibrium thermal averages at finite temperature T . The 1st term on the r. h. s reflects direct elastic response around energy minima with respect to a small externally induced shear strain. On the other hand, the 2nd term reflects relaxation of the system against the external strain at finite temperatures.
Here the distinction between T = 0 and T → 0 is important. At T = 0 only the 1st term exists. However, the contribution of the 2nd term can remain, in principle, in T → 0 limit if the strength of the thermal fluctuation of the current sin(φ ij ) is O(T ). Such a situation can arise if there are soft modes with vanishingly small energy gap such that they remain thermally active at arbitrarily low temperatures.
To obtain the equilibrium ensemble to evaluate the helicity modulus at δ = 0, we performed simulations of relaxational dynamics by numerically solving the Langevin equation with the Hamiltonian H given in Eq. (1) and ξ i (t) being Gaussian noise with zero mean and unit variance, by the 2nd order Runge-Kutta method. The system is cooled with cooling rates dT /dt = 10 −10 − 10 −9 starting from initial temperature at T = 0.3 − 1.0. We have checked that the cooling rate is slow enough by comparing with the results of 4 times faster cooling rate.