Velocity fluctuations of vortices in driven two-dimensional vortex matter

A molecular dynamic simulation is conducted for a dense two-dimensional vortex matter. At the critical moment when depinning events were concentrated and dynamic phase transition occurred, the simulation shows unexpectedly slow fluctuation in collective velocity of vortices along the direction of driving force (longitudinal direction). A detailed analysis exhibited heavy braiding of vortex flow-channels, and it was revealed that the velocity fluctuation is well synchronized with fluctuation of local density of flow-channels upon the depinning transition. In addition, vortex transport properties of a superconducting MgCNi3 single crystal seemed to show consistency between the experiment and the theory. At last, the density of flow-channels is proposed as a dynamic order parameter of current-driven vortex phase.


Introduction
Decades of studies on charge density wave dynamics [1], two-dimensional (2D) motion of electrons at high magnetic fields with disorder [2,3], and collective behavior of vortices in type-II superconductors [4] have expanded our understanding about many-particle nonlinear dynamics in a condensed matter. Among the above subjects, vortex dynamics is probably the most extensively studied subject because it is directly related to the application-oriented research about controlling critical current density. Also, the vortex matter (VM, vortices with quenched disorder) is a convenient model system as it can be easily generated by applying an external magnetic field to a superconductor. So far, there have been much effort to understand phenomena such as a peak effect (a sudden increase of depinning crictical current at the verge of superconductivity) [5][6][7][8][9][10] and a melting transition ofVM [11][12][13]. Valuable information regarding collective motion of vortices has been gathered by analyzing flux-flow noise [14][15][16][17][18] while real-time imaging of a driven VM has been a formidable task until very recently [19]. Therefore, numerical methods have been applied to understand experimentally observed fluxflow phenomena [20][21][22][23][24][25][26][27][28][29][30][31][32] On the other hand, heuristic approaches are needed to explain slow vortex dynamics in a driven VM [7,15,17,18]. For example, fast injection of vortices could introduce disordered phase spreading from the sample edge. Then, the average velocity is gradually reduced and the disordered phase is relaxed towards an ordered phase. In an ordered state, the velocity increases and injection of disorder is promoted again. The repetition of these procedures has long been proposed as the origin of low-frequency voltage (or collective vortex velocity) modulation in VM of Bi 2 Sr 2 CaCu 2 O y [15] and YBa 2 Cu 3 O 7 [33]. Also, the time-of-flight effect in Nb, which presumes ballistic propagation of vortices between sample edges, has been reported [17]. Lastly, observation of sharp sub-hertz oscillations of vortices in direct-current (dc) driven VM of NbSe 2 still remains enigmatic [18].
In this paper, we present an additional perspective of extremely slow fluctuations in collective vortex velocity without considering thermal stochastic process or injection of disorder from the edge. The model system is dense 2D VM at zero temperature. The result of molecular dynamic (MD) simulation well reproduces usually observed current-voltage (I-V ) characteristics, formation of flow-channels, and the maximum velocity fluctuation accompanied by the depinning transition. Because our particular interest lies in understanding lowfrequency phenomena, simulations are conducted over many more integration steps than previous MD simulations. To clarify meaningful signatures, white noise in velocity-time - ( ) v t series is removed by obtaining its autocorrelation function, Ψ(t), and Ψ(t) is Fourier transformed. It is notable that very slow fluctuations of Ψ(t) appears only in the depinning transition regime. Detailed analysis has indicated heavy braiding of flow-channels. In the critical region, local density fluctuations of flow-channels in transverse direction (which is perpendicular to the driving force) is well synchronized with longitudinal collective velocity fluctuations. An experimental result obtained from MgCNi 3 single crystals roughly followed a driving force dependent spectral distribution of velocity at low frequencies. Finally, it is argued that an order parameter of current-driven vortex phase can be the number density of flow-channels.

Molecular dynamic simulations
Collective motion of rigid flux tubes can be traced by observing 2D slices of the flux tubes. For MD simulations, the following viscous equation of motion is solved.
where η,  r i , r ij , r kl , U vv , U vp , and  F d are viscosity, position of i-th vortex, vortex-vortex distance, vortex-pin distance, vortex-vortex interaction potential, vortex-pin interaction potential, and driving force, respectively. If the vortex core size (or ξ, coherence length) can be ignored compared with the vortex-vortex distance, U vv would be expressed with the modified Bessel function. However, in a dense VM, ξ is comparable to the Abrikosov lattice constant (  a ) and ´ B around a vortex core, which is proportional to the repulsive force between vortices, becomes small. Therefore, although derivation of a general analytic expression for U vv is a formidable task, we can conjecture that vortex-vortex interaction would be short-ranged in a dense VM. All the length parameters are normalized by ξ, the energy scale is normalized by vortex line energy per unit length of A vv , the driving force is normalized by A vv /ξ, and the time scale is normalized by t hx = A 2 vv . Simulations are conducted by substituting a lattice parameter  x = a 4 and a Gaussian form of interaction potentials, and R α is a normalized interaction range. To check the size effect, the ratio between numbers of pinning centers ( ) N p and vortices ( ) N v is kept constant for two different sizes of simulations. The initial state is prepared by waiting until the system finds itself an equilibrium configuration in which vortices present zero mean velocity over time with To mention values of parameters which were used to obtain representative results, = ( ) R R vv vp was 1, A vp / A vv was −1, and N p / N v was approximately 1.4. Most importantly, in order to capture low-frequency behaviors, the number of iterations as well as the iteration interval (dt) are largely increased compared to those found from previous works [20-25, 27, 28]. In this work, 4th order Runge-Kutta method has been applied and the total period of a simulation is 10 5 τ with dt=0.1 τ.
In figure 1(a), representative velocity-force curves (black symbols) and corresponding differential curves are shown (red symbols). The left scale of the graph is time and number average of the longitudinal velocity, á ñ The size effect was checked by running simulations for two different VM configurations, one contains 72 vortices with 100 pinning centers and the other contains 242 vortices and 340 pinning centers. The result from the small system is marked by dotted squares and the result from the large system is marked by dotted upper triangles.
shows sharp peak in the depinning transition regime, and the critical driving force F cr is defined at the peak position. In figure 1(a), it is clear that F cr =0.87. For both systems, extensive analysis in time domain has been followed, and the result from , is recorded for a given F d . Then, to find meaningful signals other than the white noise superposed in the data, the following autocorrelation procedure has been conducted.
yt and • t is the total length of the data. When t+t′ exceeds • t during the ¢ t -average process, v y is set to zero. This is why the magnitude of Ψ(t) reduces as  • t t . Finally, Ψ(t) is Fourier transformed to reveal frequencies unique to collective motion of vortices. Notations for the autocorrelation function and its fast Fourier transformed function comply with those introduced in [34].
The leftmost column in figure 1 displays flow patterns of vortices as the driving force increases. Each figure is constructed by overlapping snapshots of the system at every 50 iterations for the last 10 4 iterations. Panels in the middle column show evolution of Ψ(t) with F d from t=2000 τ to t=5000 τ. The rightmost column exhibits fast Fourier transform (FFT), w( f ), of Ψ(t)s. Far below F cr , not even a single vortex flows ( figure 1(b)) but the subtle spectral weight in figure 1(d) is attributed to the diffusive motion of vortices which are strongly bound to pinning centers. As the system is in the vicinity of the depinning transition regime, vortices start to flow through channels: see figure 1(e). The channeling is energetically favorable because it is easy to tear a VM rather than to compress it. Say, the shear modulus C 66 is smaller than the bulk modulus C 11 [4]. When a certain vortex is driven The similar set of graphs as in along a channel, it encounters a random pinning center and experiences velocity damping. In the same channel, following vortices are almost equally spaced with an average distanceā and the velocity damping occurs regularly. This is the well known origin of a wahsboard frequency but v y (t) itself does not show clear oscillations because the pinning center is randomly distributed. The random feature due to scattered pinning centers is wiped out by autocorrelation, and figure 1(f) shows oscillations in Ψ(t) caused by the washboard effect. In addition, there exist beating because oscillation frequencies are slightly different between the two major channels shown in figure 1(e). In figure 1(g), sharp peaks are seen and the arrow points the estimated value of the washboard frequency, = á ñf v a y t N wsh , v . When the system is driven with F d =F cr =0.87, channels are heavily braided ( figure 1(h)). The autocorrelation function no longer shows sharp oscillations but exhibits very slow and large amplitude fluctuations ( figure 1(i)). The washboard frequency which must be proportional to á ñ v y t N , v cannot be found figure 1(j)). For F d >F cr , independent flow-channels begin to appear as shown in figures 1(k) and (n). Also, major spectral weight in w( f ) is strongly focused around f wsh . In figure 1(p), washboard frequencies for each flow channels are distributed since the channels become more and more independent at large F d andā for each flow channels are also distributed. What is singular in figures 1(h)-(j) is that the fluctuation time-scale of Ψ(t) in the critical regime is approximately three orders of magnitude longer than the period inherent to the washboard effect. To figure out the origin of this slow longitudinal velocity fluctuation, v y (t) and flow patterns are scrutinized. Figure 2(a) shows a part of v y (t) (black line, left vertical scale) when the driving force is F cr . Note that the graph shows fluctuation over very long periods. Technically, 1 τ corresponds to 10 iterations, and Fig 2(a) shows results during 15000 iterations. In figure 2(b), 20 snapshots of the driven VM are superposed for Δ t A =100 τ to see a flow pattern when vortcies are moving fast. Similarly, in figure 2(c), the same number of snapshots are superposed for Δ t B =100 τ to see a flow pattern when vortcies are moving slow. This is nothing but disassembling the heavily braided flow pattern ( figure 1(h)) at the dynamic critical point. From the comparison between figures 2(b) and (c), it is easy to discover that channels are closer when vortices are moving fast and vice versa. To explain, the effective number of pinning centers seen by a vortex in a narrow bundle of flow channels (say, a narrow stream) is smaller than the number of pinning centers seen by a vortex in a wide stream. Thus, the narrow stream ( figure 2(b)) experiences less resistance and flows faster than the wide stream ( figure 2(c)). The above feature is most clearly seen for ranges of parameters, 0.8<R vv <1.2 (R vp =1) and 5 . The narrowness of a flow stream is related to the local density of moving vortices or flow-channels as shown in figures 2(b) and (c). Focusing on the region where channels are strongly fluctuating in the x-direction (transverse direction), the local area of interest is chosen between x=10 and x=20 (L loc =10). The local density, n loc , is calculated as the time averaged number of flowing vortices for 100 τ inside the area L y ×L loc (L y is the full length in the y-direction). In figure 2(a), n loc for every 100 τ is overlapped (red dots), and n loc (t) is well synchronized with v y (t). In the end, the delicate competition between channel-channel repulsion due to vortexvortex repulsion and channel-channel attraction mediated by vortex-pin attraction makes the width of the local flow stream and v y to fluctuate over a very long period.
To provide plausible physical interpretations, the concept of second-order phase transition is reminded. Given that the number of flow-channels is drastically increasing as F d exceeds F cr , it seems reasonable to take the number density of flow-channels as a dynamic order parameter. In this context, figure 1(h) might be regarded as a snapshot of critical opalescence while critical exponents should be found to discuss critical phenomena in a strict sense. In the meantime, features observed in figure 2 have not been anticipated from a general theory of critical phenomena. Slow density fluctuations of flow-channels upon dynamic phase transition is analyzed long after the system is relaxed from the initial static phase. Thus, it is physically different phenomenon compared to the critical slowdown of recovery rate of perturbed VM [29][30][31][32].

Experiment
Now that anomalously slow longitudinal velocity fluctuation of driven dense VM is anticipated as a consequence of the local density fluctuation of flow-channels during the dynamic phase transition, an experimental verification will be required. One of the most important requirements would be preparing 2D superconducting specimen. But, we can also test 3D isotropic superconductor. Flux lines in MgCNi 3 single crystals are examined because these are considered to be rigid and isotropic [9]. In this case, the flux motion can be treated as practically two-dimensional. The purpose of the experiment is to investigate possible appearance of slow fluctuations in voltage developed by moving vortices.
The experiment was performed on 200 × 50 × 20 μm 3 rectangular shaped single crystal [35]. The superconducting transition in resistivity appears around 7K [35]. The upper critical field ( ) H c2 is about 10 T, and all of the superconducting properties are isotropic. The standard 4-wire method was used to probe vortex transport properties. By using photo-lithography, current terminals were deposited along short edges of the sample, and voltage terminals in strip shape were deposited between the current terminals [35]. The distance between the voltage terminals was 50μm and all four terminals were parallel to each other. The Lorentz-like force exerted on a vortex is parallel to the short side of the sample and the voltage difference is developed along the long side. This geometry of metallic terminals rules out possible appearance of flicker noise (pink noise or 1/f noise) [14] which may not be relevant to low frequency spectra due to the critical slowdown. A customized cryostat equipped with a 9 T superconducting magnet was used. In this apparatus, the sample is placed in the middle of uniform flow of low-pressure 4 He gas. Because the stream of cold 4 He gas immediately takes the excessive heat away from the specimen, an artifact related to self-heating was negligible.
The recording time for V(t) (not shown) was 300 s and the sampling interval was 50 ms. Therefore, frequency window of the Fourier space is 10 Hz wide and a washboard frequency is missing in this frequency range. For instance, with á ñ = V 1 t μV and  x = a 4 , kHz. Roughly 1 μs of time resolution is required to measure a washboard frequency and the required time resolution should increase with á ñ V t since the vortices are moving faster. However, the voltage resolution is inversely proportional to the time resolution of a dc voltmeter. With 1 μs of time resolution, voltage resolution becomes about 1 mV, which largely exceeds a level of the signal we need to measure. Admittedly, it is difficult to measure f wsh with a dc circuit, but it is also beyond the scope of this article to provide precise measurement of f wsh . nm (F • , flux-quantum) and  x  a 4 [35,36]. Therefore, vortices are as densely packed as the simulated system in section 2. Existence of three distinct dynamic regimes are deduced from figure 3(a). At currents well below I c , which value is defined at the peak position of a dV/dI, almost all the vortices are presumed to be pinned. As I is increased close to the I c , vortices begin to move and nonlinearity of the I-V curve characterizes the depinning transition. A linear I-V curve denotes that almost all the vortices are moving at high enough driving current (free-flow state). Usually, V(t) exhibits randomly looking fluctuation centered at an average value á ñ V t . A root-mean-square (rms) fluctuation was not so sensitive to the length of data acquisition time, but the current dependence of the rms value was similar to that of dV/dI showing the maximum in the critical regime. In the mean time, when figures 3(b) and (h) are compared with figures 1(c) and (o), the driving force dependence of Ψ(t) seems to be reversed. The discrepancy can be understood as follows. First, at I=3.6 mA (pinned state), the voltage level is fluctuating around zero and the signal to noise ratio is the minimum. So, contribution of the instrumental noise will be maximized and small signals characterizing diffusive motions will not be observable in w( f ). Second, remind again that the experiment is not designed for capturing washboard oscillations which would make Ψ(t) in figure 3(h) to show large amplitude oscillations as in figure 1(o). With 50 ms of sampling rate, it is not possible to observe a featured Ψ(t) with I=7.6 mA (figure 3(h)) because f wsh ∼4 MHz.
In figures 3(b), (d), (f), and (h), we present Ψ(t)s with I=3.6, 6.0, 6.4, and 7.6 mA, respectively at 2.3 K (see black symbols). Except I=3.6 mA, all other values of current are marked by blue arrows in figure 3(a). In the critical regime ( figure 3(d), I=6.0 mA), substantial amplitude in Ψ(t) is seen and large low frequency spectral weight is observed ( figure 3(e)). Similar results were obtained with slightly increased current (I=6.4 mA, figure 3(g)). Although the length of the data may not be enough to discuss possible long-term oscillations, major fluctuation time intervals (50 s in figure 3(d), 10 s in figure 3(f)) are incompatible with approximate value of 1/f wsh ∼10 −6 s. Thus, an intrinsic fluctuation time scale of the VM practically diverges inside the critical regime. At elevated temperatures, a low frequency spectral weight was not so pronounced although it is still recognizable.
At last, for field values between 4 T and 5 T at 2.3 K, experimental Ψ(t) shows similar driving force dependence as in the simulation. For a dilute VM with ã 10 ξ (B=1.3 T), slow fluctuations of Ψ(t) in a critical region were not so distinct (not shown). In the meantime, characteristic frequencies in a critical region are distributed over low-frequencies. Therefore, it is difficult to specify I and B dependent characteristic fluctuation time which might contain information about critical exponents. For now, a qualitative argument that the characteristic fluctuation time of Ψ(t) increases drastically in certain area of vortex phase diagram seems plausible. To clarify such critical region, extensive experiments are now underway with various materials encompassing heavy-fermion, iron-based, and high-temperature superconductors.

Conclusion
To summarize, molecular dynamic (MD) simulations has been conducted for a dense vortex matter (VM). Drastically increasing number of flow-channels upon the depinning transition suggests that the density of flowchannels can be a dynamic order parameter. The spectral analysis of vortex velocity in different dynamic regimes indicates that extremely slow longitudinal velocity fluctuation appears only in the depinning regime. Detailed real space and time-domain analysis has revealed that transverse fluctuations of heavily braided flow-channels are well synchronized with collective longitudinal velocity fluctuations of vortices. Hence, at least for dense VM, the braiding of flow-channels is a prerequisite for the unforeseen sluggish vortex dynamics. The driving current dependency of low frequency voltage spectrum observed from VM of MgCNi 3 single crystals roughly follows theoretical prediction. The above result might be a guidance for additional findings related to plastic flow and dynamic criticality of various kinds of type-II superconductors.