{'fraction_non_alphanumeric': 0.0273774282071782, 'fraction_numerical': 0.006828150997212559, 'mean_word_length': 4.747640337804272, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Simple Model for Plastic Dynamics of a Disordered Flux Line Lattice

We use a coarse-grained model of superconducting vortices driven through a random pinning potential to study the nonlinear current-voltage (IV) characteristics of flux flow in type II superconductors with pinning. In experiments, the IV relation measures flux flow down a flux density gradient. The work presented here treats this key feature explicitly. As the vortex repulsion weakens, the vortex pile maintains a globally steeper slope, corresponding to a larger critical current, for the same pinning potential. In addition, the magnitude of the peak in the differential resistance falls as the resistance peak shifts to higher currents. The model also exhibits so-called"IV fingerprints", and crossover to Ohmic (linear) behavior at high currents. Thus, many of the varieties of plastic behavior observed experimentally for soft flux line systems in the ``peak regime'' are reproduced in numerical simulations of the zero temperature model. The nonlinear transport behaviors are related to the self-organized, large scale morphologies of the vortex river flow down the slope of the vortex pile. These morphologies include isolated filamentary channels, braided rivers, and flooded rivers. We propose that these self-organized morphologies of flux flow down a flux gradient govern the various plastic flow behaviors, including nonlinear IV characteristics, observed in type II superconductors.


I. INTRODUCTION
Collective transport in disordered media is a widespread and poorly understood phenomena. A great deal of experimental and theoretical effort in this area has been devoted to studying the nonlinear dynamics of the disordered flux line lattice (FLL) in type II superconductors. The FLL exhibits a threshold behavior due to the competition between pinning and flux line repulsion [1,2]. In response to a force, such as that associated with a transport current, the three dimensional FLL can move smoothly via elastic deformations, maintaining its integrity and order. However, in another regime the FLL deforms plasticly. In that regime disorder becomes more important, and the moving FLL manifold tears as some flux lines (in two-dimensions, vortices) move while others do not. As a result, the flow pattern breaks up in a nonuniform way. It is generally believed that the microscopic structure of the FLL, or defects in it, is fundamental to transport behavior both in the plastic and elastic regimes.
Part of the attention to the dynamics of a moving FLL has been motivated by interests in possible, exotic phase transitions and glassy phases, melting and other complicated scenarios associated with structural order in the FLL. Most experiments, however, are essentially transport studies, and are, quoting Higgins and Bhattacharya [3], "notoriously ill-suited for the study of thermodynamic phase transitions. These experiments yield direct information only about the mobility of flux lines, i.e. on dynamics and pinning, which would then have to be connected, through highly model-dependent ways, to the structure of the state they pertain to." Here, we show that many of the empirical results found in transport studies of plastic flux flow in type II superconductors may be obtained with an extremely simple model [4]. It has been recognized for many years that, in the presence of pinning, magnetic flux in type II superconductors forms a pile with an overall, global slope, akin to a sandpile. Penetration of magnetic flux into superconductors driven solely by the flux density gradient has been described using molecular dynamics (MD) simulations [5][6][7][8][9][10], and by the model used here [4]. However, the flux gradient has not been taken into account in any previous numerical simulation studies of the current-voltage characteristic. One possible reason is that previous numerical studies of flux motion at the scale of the vortex cores have not been able to reach a sufficiently large sys-tem size. The Bassler-Paczuski (BP) model [4], on the other hand, is a coarse-grained model and describes the magnetic flux dynamics at the much larger scale of the London length, making the large system size limit much more accessible.
The numerical simulations presented here of the BP model [4] show that nonlinear behaviors, characteristic of experimental transport measurements of plastic flow in superconductors, arise as a result of vortex flow down a vortex density gradient. The effect of the transport current is modelled by a shift in boundary conditions, which leads to a generalized "tilt" of the vortex pile. Eventually the "tilt" is sufficiently great so that some vortices can flow down the pile in the steady state, leading to the onset of a finite voltage. The vortex flow forms a variety of river morphologies depending on the interaction strength between the vortices, compared to pinning, and the overall rate of flow. The flow patterns of magnetic flux are self-organized together with the magnetic flux profile, which is the substrate on which the flow takes place. Since the BP model does not contain any detailed information on the positions of the vortex cores at the micro-scale of the FLL, it cannot exhibit any structural ordering or disordering behavior.
Although the BP model can also be studied at finite temperatures, here we use the zero temperature limit where thermal fluctuations of the flux motion may be ignored. The zero temperature approximation seems reasonable to describe the plastic transport dynamics of the low temperature superconductors and, perhaps, some aspects of the high temperature superconductors as well. Also, we consider the limit where the depairing current density, j 0 , is extremely large compared to the critical current density, j c . This corresponds to the so called "weak pinning" regime.
We argue that many of the varieties of collective transport dynamics observed in superconductors may be generic to repulsive particle systems driven through a disordered media. These behaviors are directly related in our coarse-grained model to large scale morphologies of flow down a vortex density gradient and changes from filamentary strings, to a braided river, to uniform flow at high applied currents. Since the BP model arguably contains the essential physics of the disordered flux line system, at a coarse-grained level, and reproduces a wide variety of experimental results on transport properties, we propose that these self-organized, large scale flow morphologies are also governing the nonlinear, plastic dynamics in the actual physical system: flux lines driven through a superconductor with a disordered pinning landscape by an applied transport current.

A. Summary
In the next section, we focus attention on an extensive study of 2H − N bSe 2 as summarized by Higgins and Bhattacharya [3]. The reader may already compare Fig. 1 and Fig. 2, which are schematic reproductions of experimental results, with Figs. 6 and 7, which are results from numerical simulations presented in this work. Section III presents a short review the standard theoretical picture, associated with the work of Shi and Berlinsky [11], on tearing of the FLL and the resultant nonlinear IV as found via molecular dynamics (MD) simulations [12][13][14][15][16][17][18]. These MD simulations do not, however, take into account a vortex density gradient, and do not agree with some the varieties of plastic dynamics observed in experiments. Section IV describes the original coarse-grained model, first used to describe vortex avalanches [4] and vortex rivers [19] as seen in experiments [20,21], and explains the exact model used here to describe the transport experiments summarized in Section II. To eliminate all potentially spurious sources of noise, a completely deterministic variant of the original BP model is used. To describe the IV experiments, we apply a shift in boundary conditions or generalized "tilt" of the vortex pile to represent the effect of a transport current. The resultant vortex flow represents the measured voltage. Section V contains the main numerical results and comparisons with both experiments and MD simulations. The last section summarizes our conclusions.

II. BRIEF SUMMARY OF TRANSPORT MEASUREMENTS
In a series of papers, Bhattacharya and Higgins [3,[22][23][24] describe experiments on the nonlinear transport properties of the FLL in the anisotropic superconductor 2H − N bSe 2 . In this system, the pinning is extremely weak (j c /j 0 ∼ 10 −6 ), the lattice is well-formed and a robust "peak effect" occurs slightly below H c2 . Since the London length, λ, is much less than the film thickness, the system operates in the three dimensional regime, where the flux line interactions decay exponentially for lengths larger than λ. They use the strong magnetic field dependence of the critical current (which they interpret in terms of a changing rigidity of the FLL) to explore the crossover between different type of dynamics including "elastic" flow, "plastic" flow, and "fluid" flow. In this regard, the material they study is an ideal experimental system, allowing the exploration of very different regimes in a well-controlled manner.
Near the upper critical magnetic field, they observe a pronounced peak in the pinning force. This is referred to as a "peak effect". (It should be distinguished from the peak in the differential resistance to be described later.) Equivalently, the critical current, I c , where some flux lines start to move, increases as the external magnetic field, H, increases. For I < I c the entire flux system is pinned, while for I > I c , some magnetic flux lines flow across the sample leading to the onset of a finite voltage, V . We will focus mostly on results associated with the current-voltage (IV) relation as the external magnetic field, H, is varied. As shown in Fig. 1b of Ref. [22], in the "peak regime" the IV characteristics of the superconductor vary enormously. For clarity, these experimental observations are illustrated schematically in Fig. 1. First, on increasing H, the critical current, I c , increases. This is the "peak effect". Second, below some threshold magnetic field, the IV curves always rise concave upward from I c . This is the generic form of the IV for an FLL that is usually reported in the literature, and it is associated with an "elastic" regime. Third, when the external magnetic field enters the "peak regime", the IV curves change drastically, starting as concave upwards but then bending over after a pronounced inflection point associated with a change of curvature. In the plastic regime, one obtains a characteristic S-shape IV curve. Close to onset, I c , it is concave upward, but then bends over as I increases further, saturating to a finite slope at large currents. Between the S-shaped and elastic IV curves may be a special curve which is always convex for I > I c . This appears at approximately 5.8T in the experiments where the inflection point has moved to onset. Finally, above another magnetic field value, the inflection point is at currents larger than those used in the experiments, and there is no saturation in the slope of the IV curves that can be observed. The numerical simulations presented here reproduce the entire progression of these curves in the plastic regime and their changes as a parameter in the model, representing vortex interactions, is varied. It does not reproduce the behavior in the "elastic" regime, as explained later. The elastic regime is also not observed in twodimensional MD simulations of driven vortices near the onset of flow. Experimental measurements of the differential resistance, R = dV /dI, reveal a peak in R, corresponding to the inflection point for the S-shaped IV curves. The magnitude of this peak diminishes as the external magnetic field approaches H c2 , and the position of the peak shifts to higher currents. Except for the largest magnetic fields used, the differential resistance eventually saturates at constant values at large enough currents, indicating Ohmic or fluid like behavior for sufficiently high driving. Also, in the plastic regime corresponding to the S-shaped IV curves, the usual scaling ansatz associated with dynamic critical phenomena V ∼ (I − I c ) β does not appear to hold. These experimental observations are illustrated schematically in Fig. 2. The numerical simulations presented here reproduce this precise pattern of behavior for the differential resistance in the "peak" regime, including the changing position (to higher currents) and decreasing magnitude of the peak resistance as vortex interactions weaken, and the saturation to Ohmic or fluid-like behavior at high currents.
Over a narrow range of parameters in the plastic regime, Bhattacharya and Higgins observed jaggedness in the differential resistance. This corresponds to secondary peaks in addition to the main peak in the differential resistance. As the external magnetic field is varied, the peaks can be made to appear and disappear, but for a given value of magnetic field and for a given sample, the peaks are reproducible, and thus act as "fingerprints" of the underlying pinning disorder. Our model also exhibits such IV fingerprints.

III. SOME COMMENTS ON CONVENTIONAL THEORY
The usual interpretation of the peak effect follows the original proposition of Pippard [25], that softer systems are pinned more strongly than more rigid ones. This picture has been put on a formal basis by Larkin and Ovchinnikov's theory of collective pinning [26]. It was found that critical current increases as the external magnetic field increases near H c2 because the FLL elastic moduli soften. However, collective pinning theory doesn't explain the shape of the IV curve outside of the elastic regime. Some advances have recently been made by Le Doussal and Giamarchi to probe the transition to plastic flow and other features [31].
MD simulations have indicated that, at least for the two dimensional case, the pinning forces and IV characteristics are determined by plastic deformations (tearing of the FLL) which fall outside the region of validity of collective pinning theory [12,13]. It is believed that the two dimensional elastic system is always unstable at onset, when a finite voltage first appears, and flow must take place in the plastic regime. These simulations can also be interpreted as describing two dimensional cross-sections of the three dimensional FLL in the plastic regime. The scale of the MD simulations is that of the vortex cores, which is much smaller than the London length, λ.
Shi and Berlinsky [11] presented a limited analytic treatment of the way in which "lattice defects alter the flow characteristics of the lattice under the influence of an external drive". However both the density of lattice defects and the IV relation were determined by numerical MD simulations. These simulations showed that at high currents where the IV relation is linear, the defect density drops. It increases sharply at lower currents where the IV relation develops an S-shaped curve.
Further simulations have shown that the two dimensional plastic flow takes place in terms of rivers of moving vortices separated by islands where the vortices do not move [27,28]. This channel flow behavior has been observed experimentally using Lorentz microscopy [21]. Previous simulations of the coarse-grained BP model revealed results consistent with these and also showed that the channels can form a braided river which exhibit selfaffine (multifractal) behavior similar to fluvial braided rivers (see for example Ref. [29]).
Nori and collaborators first studied, using MD simulations, flux driven into superconductors with random pinning, with the driving force solely due to the flux-density gradient. They elucidated many properties of the Bean state including the magnetic field profile, magnetization hysteresis loops, critical currents, vortex avalanches, and vortex rivers [5][6][7][8][9][10]. None of these studies using an open system with an overall density gradient reported the IV characteristics, though.
Recently, Nori et al [16] have simulated the IV curve as the FLL softens by varying the vortex interaction parameter (see also Ref. [30]). As in all previous MD studies of IV behavior of the FLL [12][13][14][15][16][17][18], the vortices are contained in a periodic system, where they can neither enter nor leave. The initial condition is an ordered vortex lattice. Motion takes place via MD updates, with a uniform force applied to all vortices.
The most important difference, besides the scale of the model, with our description of the IV experiments, is that, in all of these MD studies of the current-voltage relation, no overall vortex density gradient can develop owing to the periodic boundary conditions. This makes the vortices travel perpetually around the system and forces the paths to circle, which is unphysical, and does not occur in the real system. The actual physical situation is an open system with flux pushed in and out, rather than a periodic one. The physical system adjusts its overall magnetic flux profile in response to applied forces. This is not possible in a periodic, closed system. Nevertheless, the basic result that the critical current increases as the vortex interactions weakens is obtained. However, all the IV curves measured via MD simulations of periodic systems, fall on top of each other, or overlap, at high currents, for different values of the vortex interaction strength. This is inconsistent with the experiments of Bhattacharya and Higgins (see Fig. 1), and reflects the fact that the artificially periodic system is not able to adjust its profile in response to applied forces. An even more significant, but related, difference is that the peak in the differential resistance grows monotonically and gets sharper as the resistance peak shifts to higher currents, the exact opposite of what happens in experiments.
In short, the systematically changing morphology of the IV curves and differential resistance for a soft FLL has not been theoretically described.

A. Motivation
We take a different approach to modeling flux flow in superconductors. In order to describe collective transport in disordered media, it is reasonable to consider coarse-grained models where the details of the precise interactions between flux lines (or vortices) are lost but the general effects of repulsive interaction between granular or discrete objects, pinning, and over-damped motion leading to stick-slip dynamics (tearing) are preserved.
The BP model is an interacting sand-pile model of vortices in a type-II superconductor, where the "sand" grains, representing magnetic vortices, repel each other. It was originally motivated by the observation of (possibly) self-organized critical [32] avalanches in field ramping experiments [20], where the distribution of flux packages falling into the interior coil of a hollow cylindrical superconductor were measured. Actually, the similarity between the Bean state [33] (or vortex pile) and sandpiles was first pointed out by de Gennes [34]. Later, Vinokur, Feigel'man, and Geshkenbein [35] suggested that thermally induced flux creep would lead to a self-organized critical state in a type II superconductor, as did Tang [36]. The key observation is that flux line flow always takes place on a flux pile which has an overall density gradient. This pile may be in a self-organized critical, or some other nontrivial state.
In addition to describing vortex avalanches in field ramping experiments [20], the BP vortex model has also been used to describe flux noise [37], vortex avalanches in the presence of a periodic, dense array of pinning centers [38], thermally activated flux creep [39,40], and magnetization loops [41].

Coarse-graining
Consider a transverse two-dimensional slice of a superconducting slab at T = 0. The BP model [4] results from a coarse-grained description to the scale of the London length, λ, of the microscopic dynamics, and incorporates the features that are essential to produce the observed complex behavior: Repulsive interactions between vortices, variations in the pinning potential, and variations in the vortex density -all at the scale of λ. One can imagine imposing a grid of cells on the system; vortices in the model correspond to a vortex number in an extended region of the actual physical system. Pinning in the model corresponds to a number of point pins in an extended cell. Each lattice site in the model can hold many vortices, and can have a different, albeit quenched, pinning potential, due to the underlying randomness in the positions and strengths of the microscopic pinning centers. The model allows many vortices to interact with each other while maintain locality (at the scale of λ) in interactions. It is the only model of this sort that has been proposed to describe flux dynamics in superconductors. It enables numerically studies, using ordinary workstations, of the steady state and transient properties of systems larger than (500λ) 2 containing tens of millions of vortices.

B. Definition
The BP model is defined as follows (see Fig. 3). Consider a two-dimensional honeycomb lattice [42] where each cell, x, has three nearest neighbors, and is occupied by an integer number of vortices, m(x). The total energy of the vortex system includes the repulsive pairwise interaction between the vortices and the attractive interaction of vortices with the pinning potential,V . For a given configuration of vortex number {m(x)}, the total energy of the system is: Since the model describes a system coarse-grained to the scale of the London length, the repulsive interactions J ij are short-ranged. This usually includes an on-site interaction, and a weaker nearest-neighbor interaction.
As in the microscopic case described by MD simulations, the change in the total energy of the model when moving a unit vortex from one site to a nearest neighbor site is determined. This yields the force to move a unit vortex from x to y, which is As indicated in Fig. 3, the nearest neighbor cells of x are y, x1, and x2, and the nearest neighbors cells of y are x, y1, and y2, and 0 ≤ r < 1. A slightly different implementation of the disorder is used than before. The normalized pinning potential V pin (x) is a random number taken from a uniform distribution in the interval between zero and V max . In each time step, all cells of the lattice are updated in parallel. A single vortex moves from a cell to a neighboring cell if the force in that direction is positive, or equivalently if the total energy of the system is lowered. In Eq. 2, the units of force on a vortex have been normalized so that the on-site term is unity. Thus there are two dimensionless parameters remaining, V max and r.
x y x1 x2 y1 y2 FIG. 3. The two-dimensional honeycomb lattice. Each cell, x, has three nearest neighbors, and is occupied by an integer number of vortices, m(x). The force pushing a vortex from cell x to cell y is calculated by taking the discrete gradient of the sum of two potentials, one representing the repulsive interaction between vortices occupying the same and nearest neighbor cells, and the other representing the attractive interaction between vortices and pinning centers.
Many alternatives exist to handle the situation when more than one unstable direction appears for a vortex to move. In the previous implementation of the model, one unstable direction was chosen at random. In order to simplify the model and eliminate all potentially spurious sources of noise, the most unstable direction that has the largest force is chosen and the vortex moves to that site. This represents an extremal process. In fact the entire model is now completely deterministic, corresponding to a T = 0 limit of the dynamics.

C. The External Magnetic Field: Building a Vortex Pile
Flux lines enter the superconductor from the edges, pushed in by the external magnetic field. This is represented by putting all sites on the left edge of the model in contact with a reservoir of vortices at some potential, corresponding to the external magnetic field on the left side of the sample. The same is done for the right side of the sample. For simplicity, periodic boundary conditions are used for the top and bottom. If the two reservoirs are set equal, representing embedding the sample in an external magnetic field, vortices enter the system generating the classic V-shaped flux density curve as the external magnetic field is increased (see below).

Details about Boundary Algorithm
In order to calculate the force on a vortex to move to or from a boundary site, a special algorithm must be used because one of the nearest neighbor sites of each boundary site is not on the lattice. The rule used here simply assumes that the "virtual" off lattice site neighboring each boundary site is occupied with an equal number of vortices as that boundary site. More specifically, at the beginning of each lattice update, all of the sites on the boundary are set to be occupied with the same number of vortices. The lattice update then proceeds, during which vortices can move off the boundary sites into the system, or from the system onto the boundary sites, thereby changing the number of vortices occupying a boundary site. However, at the beginning of the next lattice update all of the sites on the boundary are reset to their original value. Through this process, vortices can be removed or added to the system. In general, the left and right boundaries are held at different values. There is no pinning at the boundary sites.

Details about Parallel Update
An artifact of parallel updating is the existence of local instabilities in which two, or more, vortices oscillate back and forth between neighboring sites. These local instabilities disappear if the model is coarse-grained, because then the neighboring sites are incorporated into a single one, and therefore are not important to the large scale behavior of the system. The instabilities can be eliminated by keeping track of the direction from which the last vortex moved onto each site and always forbidding a return movement backwards in that direction. Otherwise, backward jumps are rare and therefore disallowing them does not change the large scale behavior of the system. A similar rule applies to the boundary sites. The advantage of the parallel update is that it is numerically more efficient, than other update schemes and does not introduce any uncontrollable spurious effects.

D. The Transport Current: Shifting the Boundary Conditions
We consider an infinite slab of finite thickness in a parallel applied magnetic field, carrying a current perpendicular to the field. Depending on the direction of the applied current, the magnetic field on one side of the superconductor (e.g. the right hand side) will be decreased and the magnetic field on the other side will be increased (e.g. the left hand side). Assuming the applied magnetic field is sufficiently strong, this corresponds just to changing the heights of the magnetic flux pile on the left and right edge, or a general kind of "tilt" of the pile. See, for example, Ref. [43]. Of course, the internal currents and forces inside the superconductor will readjust to accommodate this new boundary condition. The flux lines are considered to be perfectly stiff, and described by a two-dimensional slice of the three dimensional slab. An applied magnetic field in the geometry described above will produce the well known "V profile" of the magnetic flux density discussed, for example, in Orlando [44] and observed in many experiments, such as those of Behnia, et al [45]. This is shown in Fig. 4. The actual magnetic profile depends on the history of the sample and how magnetic field has been applied in the past. Applying a finite current shifts the boundary conditions resulting in a hysteretic profile, also shown in Fig. 4.

Onset
Eventually, as the applied current is increased further, a critical current is reached, where the shift of the boundary conditions is so large that steady vortex flow occurs down the gradient spannning the entire sample. A profile just above the critical current is also shown in Fig. 4. Note that this profile also retains some hysteretic properties. For example, it has a bump corresponding to where the increasing and decreasing portions of the magnetic field profile merged as the external current was increased above the threshold. This bump disappears if the boundary conditions on the sample are shifted further and then lowered back to the previous value.

E. Making IV Measurements
The IV characteristic is determined by the relation between applied current and the vortex flow, which induces a voltage. To our knowledge, this type of numerical measurement, made by shifting the boundary conditions on the Bean state, has not been investigated before. Here the IV characteristic is simply the relation between the magnitude of the shift (representing an applied current) and the average flow rate of vortices (representing a voltage) when the critical current (tilt) is exceeded. In order not to confuse the reader we use the term current to refer to the applied electrical transport current and the term flow to refer to the motion of magnetic vortices.
In general, the boundary sites in the model can be set to any real value, including non-integer values. This is because vortex number on the boundary sites describes the external magnetic field density at the boundary of the sample. However, only integer units of magnetic flux can enter the interior of the system, and thus only integer numbers of vortices occupy interior lattice sites. Obviously, the magnitude of the difference between the heights of the left and right boundaries, can also be set to non-integer values. However, if the boundary heights are shifted by less than whole integers, steps can appear in the IV data. These steps are caused by the fact that only whole numbers of vortices can occupy interior lattice sites. Since one vortex unit in the model represents many actual physical vortices, this is to some extent an unphysical artifact of the model, and the discreteness effect should be less apparent in experiments. To eliminate this effect, all of the IV data presented in this paper were calculated by shifting the boundary heights by only integer numbers of vortices.
In simulations of the model, there is almost no backward movement of vortices. Thus, the vortex flow can be determined by measuring the average number of vortices moving per lattice update (the average activity). Furthermore, the velocity of each moving vortex is one lattice site per update. Thus, the experimentally measured voltage, which is equal to the amount of moving flux times the velocity of that flux, is therefore proportional to the average vortex activity. IV data from simulations of different size systems collapse nicely in a scaling plot, Fig. 5, if the voltage is measured as vortices moving per lattice update per lattice site (the average activity per site), and the current is measured as the average slope across the system (the magnitude of the height difference between the left and right boundaries divided by the length of the system). In the following sections, IV data is presented in these scaling units. In Fig. 5, the IV data was produced by repeatedly increasing the height of the left boundary by 1 vortex and lowering the height right boundary by 1 vortex. The vortex interaction strength was r = 0.1.
The IV data presented in the following sections was calculated in a similar fashion. All of the data is for systems of size 200x400. Each IV data point was calculated by first shifting the left boundary height up 1 vortex and shifting the right boundary height down 1 vortex. As for Fig. 5, the lattice was updated 20000 times to eliminate transient behavior, and finally the lattice was updated another 20000 times during which the average vortex activity was measured.

V. COLLECTIVE TRANSPORT BEHAVIOR OF THE MODEL
The IV relation was measured for different values of the parameters and for different system sizes. Our main result is shown in Fig. 6, where the parameter V max = 5 and the parameter r is varied. The parameter r represents the strength of repulsion between vortices at nearest neighbor cells (of size λ). The first result is that as r decreases, the critical current, which is the slope of the pile where vortices first start to flow, increases. Clearly, applying an increased r to the pile in the steady state lowers the slope of the pile since formerly stable local slopes will now become unstable due to the increased repulsion between vortices at neighboring sites. Thus, fixing all other parameters, we can identify the parameter r in the model as a way of controlling the critical current, or slope of the pile. In the superconductor, the critical current can be controlled by the applied magnetic field. In the peak effect regime, it turns out that increasing the applied magnetic field leads to a softer FLL and thus a higher critical current. Therefore, the regime where the critical current is an increasing function of the applied magnetic field is represented in our model by a critical current which is a decreasing function of the parameter r. This is made evident by comparing Figs. 1,2 and Figs. 6,7.
The second result is that, except for the largest r, all of the IV curves have a characteristic S shape. They start out at some I c increasing concave upward and then bend over saturating to a finite slope at large currents. The IV relations for a given realization of disorder, however, do not overlap at high currents, unlike the results obtained with previous MD simulations [16]. Again, all of this agrees with the experimental results. In Fig. 7, the differential resistance dV /dI is calculated by taking numerical derivatives of the curves shown in Fig. 6. At the largest r value, the resistance peak is very large. As the parameter r decreases, the peak moves to higher currents and decreases in magnitude. In fact the same behavior is observed in experiments as the external magnetic field approaches H c2 . As mentioned before, the MD simulations give the opposite result of a peak that increases in magnitude as it shifts to higher currents [16].
In the actual experiments, at fields below the peak effect regime, the vortex flow is believed to be elastic and there is no observed peak in the differential resistance. As the magnetic field increases a peak in the differential resistance starts to develop, which then reaches a maximum, decreasing again for larger fields. Our simulations appear to describe the experiments once the peak in differential resistance has reached its maximum.
In order to obtain an elastic regime, we could consider a three dimensional coarse-grained model of repelling flux lines rather than point vortices representing a two-dimensional cross section of that system. Work is in progress along those lines.
The model discussed here does not describe the behavior of the superconductor in fields greater than that which gives the largest critical current, where the critical current decreases as the magnetic field increases. This may be due to the fact that we take the depairing current to be strictly infinite (i.e. j c /j 0 = 0) and there is no transition to a non-superconducting state in the model presented here.
Note that the differential resistance curves for small r contain secondary peaks, in addition to the main peak. This is similar to the jaggedness or "fingerprint" found in experiments. This jaggedness in our results occurs in the filamentary channel regime, discussed below, and is due to filaments opening and closing as the applied current increases.

A. Relation to River Morphology
So far we have only characterized the model using measurements analogous to those that experimentalists typically have available. However, numerical simulations can also provide a great deal of easily accessible information about the morphologies of flow patterns associated with different configurations. In fact, different morphologies of vortex flow patterns have been observed using MD simulations [8][9][10]12,15,30,46]. To understand the nature of the differences in the shape of the IV curves, we have examined how the flow morphologies change as one increases the applied current for a system that has an Sshaped IV curve, and also for the singular IV curve that occurs at large r.
The paths that the vortices take as they cross the sample can be determined by measuring the average activity at each site. For example, Fig. 8 shows a series of grey scale images of the vortex flow patterns. These images represent "time-lapsed photographs" of the vortex activity. The different images in Fig. 8 show the vortex flow patterns for the case r = 0.1 as a function of increasing external current, I. In these images, sites with no vortex activity are blank, sites with an activity greater than or equal to 0.5 (one vortex moves every other lattice update) are black, and sites with activity between 0 and 0.5 are indicated by grey dots with a darkness proportional to their activity. These patterns are fixed and do not change in time with fixed external driving conditions.
The nature of the flow can be quantified by the distribution of activity at the different lattice sites. Figure  9 shows histograms of the activity corresponding to the images of Fig. 8. These histograms are constructed with 1000 bins and normalized so that the area under the curve is equal to one. Note that there are peaks corresponding to zero activity not included in the figure. Those peaks at zero activity decrease in size as the current is increased.
As can be seen in Figs. 6 and 7, the IV curve for r = 0.1 is S-shaped. Figure 8a shows the vortex flow pattern just above threshold for that case. The vortex flow takes place only on a single filamentary string with some small side branches. This behavior is also evident in the histogram Fig. 9a, which shows a small number of isolated peaks. As the tilt of the pile increases the number of filaments of the vortex flow increases, and they begin to merge. This process can be seen in Figs. 8b and 8c. The corresponding histograms of activity shown in Figs. 9b and 9c indicate an increasing number of peaks, and the development of a continuous distribution. During this merging process, while the vortex flow remains filamentary, the IV curve remains concave upward, and the differential resistance rises.
Eventually the filaments of vortex flow merge to form a braided river, as shown in Fig. 8d. This occurs around the peak in differential resistance, corresponding to a change in curvature of the IV relation. In this case, the corresponding histogram of activity (Fig. 9d) shows a continuous distribution of activity peaked at zero. Thus the peak in differential resistance signals a change in the underlying vortex flow morphology from filamentary strings to a braided river. This appears to be consistent with the observation of Kolton et al that at the peak in differential resistance "all the vortices are moving in a seemingly isotropic channel network with maximum interconnectivity" [18].
As the tilt of the pile is increased even further, the vortex flow becomes spatially more and more uniform, as the braided river floods. This uniform flow region corresponds to linear, or Ohmic behavior. These results can be seen in the river flow pictures in Figs. 8e and 8f. The corresponding histograms shows a peaked continuous distribution, which has separated from zero activity. As the tilt increases further the continuous distribution of activity obtains an increasing mean and narrowing width. The transition to Ohmic behavior comes about as the braided river floods at higher flow rates than can be supported by such a structure. A similar progression of flow morphologies, from isolated filamentary channels, to a braided river, to a flooded river with uniform flow, occurs for other values of r with S-shaped IV curves. However, a different scenerio occurs at large values of r where the IV curve is not S-shaped, but is instead always concave downward from the onset of vortex motion. The IV curve for r = 0.4 in Fig. 6 is an example. We will refer to this as the critical curve; its functional form is discussed in the next subsection. The differential resistance for the critical curve is discontinuous, as can be seen in Fig. 7.
For the critical curve, there is no filamentary region of vortex flow. Instead, right at the onset of vortex motion, the flow has the structure of a braided river. This can be seen in Fig. 10a, and in the corresponding activity histogram in Fig. 11a which shows a continuous distribution of activity. (Figures 10 and 11 were produced in the same manner as Figs. 8 and 9, respectively.) As the tilt of the pile is increased even further, the flow again increases and becomes more spatially uniform, similar to behavior of the S-shaped IV case after the peak in the differential resistance, where the braided river floods. This behavior can be seen in the river flow pictures in Figs. 10b, 10c and 10d, and the corresponding histograms in Figs. 11b, 11c and 11d. Since the critical curve has a diverging differential resistance at onset, this further supports our observation that the peak in the differential resistance is associated with a braided river structure.

B. Scaling of the Critical IV Curve
The critical IV curve for r = 0.4 is well described by with β = 0.6±0.1, as shown in Fig. 12. For that measurement, the critical current, I c , was measured to be 0.095, which was the largest current found with zero voltage. Since the current was sampled at values spaced by 0.005, the uncertainty in I c is approximately ±0.005. Varying the value of I c over that range changes the value of β that best fits the data near the onset of vortex motion, and allows an estimate of the error on β.
The exponent β can be related via scaling arguments to the exponents characterizing the distribution of avalanches in the self-organized critical state. In the self-organized critical state, the average flow rate, V , of vortices is controlled rather than the overall slope of the system. The scaling argument is similar to that used in Corral and Paczuski [47] to describe the transition from avalanche to continuous flow in the one dimensional Oslo model [47]. The excess slope above the critical slope at onset is ∆m = I − I c = (∆N )L 3 , where ∆N is the excess number of vortices in the system, and L is the system size. If vortices are added very slowly then there will be distinct avalanches separated by intervals of no activity. In that regime, superposition applies and ∆N ∼ V L z , where L z is the cutoff in the duration of the avalanches. In the rapidly driven regime, the avalanches overlap, and the excess slope becomes independent of system size, thus ∆m ∼ V 1/β . These two limits can be combined into a single cross-over scaling function: The exponent x in this expression measures the average duration of avalanches t ∼ L x . Since avalanches in this case come about from adding an entire row of L vortices to the system, the average size of avalanches, measured in terms of the number of topplings is s ∼ L 2 , rather than s ∼ L, as in the Oslo model. Obviously, on average each site in the system topples once when a row of vortices is added in the self-organized critical state. Using this result together with conservation of probability gives a scaling relation between x and z; x = z + 2 − D .
Combining all these results gives Using the exponent values z = 1.5 and D = 2.7 obtained in Ref. [4] gives β = 0.6, in good agreement with the numerical result presented in Fig. 12.

VI. SUMMARY
We have studied the nonlinear current-voltage characteristics of flux flow in type II supeconductors with random pinning using a simple, cellular model. As in the physical system, vortices flow down a flux density gradient. Because the coarse-grained model does not include any degrees of freedom at scales much smaller than the London length, λ, it can not describe structural ordering or disordering of the flux line lattice. Despite this fact, simulations of the model reproduce many of the empirically observed transport characteristics of plastic flux flow in type II superconductors.
In particular, our results reproduce many of the features attributable to plastic flow of the FLL in the "peak regime". By weakening the vortex interaction strength, we find an increase of the critical current and a falling of the magnitude of the peak in the differential resistance. The model also exhibits IV fingerprints, and crossover to Ohmic or linear behavior at high currents. Also, the IV curves for different vortex interaction strengths do not merge at large external currents. All of these features are completely consistent with experimental results. However, presumably due to the two-dimensionality of the model, the elastic behavior of the FLL is not reproduced.
The success of these efforts to describe the plastic transport behavior of magnetic flux in superconductors with a coarse-grained cellular model suggests a possibly generic explanation for and ubiquity of plastic flow phenomena observed in superconductors. It may not depend on the degree of disorder or defects in an underlying microscopic flux line lattice. Instead the behavior may be common to driven repulsive particle systems in a disordered media. The varieties of plastic flow behaviors in the model studied here result from the changing morphologies of the vortex flow pattern down a vortex density gradient, self-organizing into distinct large scale patterns. These include isolated filaments, which merge at higher flow rates to give a braided river, and lead to uniform flow, or Ohmic behavior, at the highest flow rates. The filamentary structure is associated with a concave IV characteristic; the braided river structure is associated with the peak in the differential resistance; the change to Ohmic behavior comes about as the braided river floods as it can not support a higher level of flow.