Numerical Implementation

In this section we describe our numerical implementation. Simulations were carried out on a square region with 7250 randomly placed pins. Having constructed a random pinscape, we placed N = 65 flux lines at random and solved numerically for stable equilibrium configurations of Eq. (5) using a three stage procedure: (i) Starting from a random initial configuration, the vortex system is allowed to relax dynamically according to

$$\frac{d{\bf {r}}_i}{dt} = {\bf F}_i({\bf r}_1,{\bf r}_2, \ldots, {\bf r}_N),$$ (9)

until the magnitudes of all forces on the RHS are less than a prescribed target value. (ii) Once this has been achieved, the resulting configuration is used as an initial guess in a Newton-Raphson algorithm to find a solution to Eq. (5). Solutions found this way are not guaranteed to be stable with respect to small perturbations of the particle locations. Therefore, (iii), a linear stability analysis is carried out to check for stability. If the configuration obtained this way does not turn out to be stable, stage (i) is resumed with a reduced target value. The whole procedure is repeated until a stable configuration has been found. The combination of a molecular dynamics algorithm with a Newton-Raphson scheme, as described above, turns out to be far more computationally efficient than using molecular dynamics alone [9].

Lengths were chosen so that the square corresponds to a $5 \times 5~\mu$m² field of view in a magnetic field of roughly 50 Gauss. With this choice, all lengths will be reported in $\mu$m. The London penetration depth was taken to be $\lambda = 0.04~{\mu}$m. These choices resemble the experimental conditions in Ref. [3]. The parameters varied in the simulations were the range of the pinning forces r_p, with values $r_p \in [0.015,0.05]~\mu$m, and the distribution of pinning strengths $V_{\alpha}$, cf. Eq. (4). We considered the following distributions:

$${\cal P}(V_\alpha) = \delta(V_\alpha - V_o),$$ (10)

for identical pinning strengths, and

$${\cal P}_\nu (V_\alpha) = \frac{c_1}{V_o} e^{-(c_2 \frac{V_\alpha}{V_o})^\nu}, V_\alpha > 0,$$ (11)

for stretched-exponentially distributed pinning strengths, with the constants c₁ and c₂ given by $c_1 = \nu c_2/\Gamma(1/\nu)$, and $c_2 = \frac{\Gamma(2/\nu)}{\Gamma(1/\nu)}$, where $\Gamma(x)$ is the Gamma function. In particular, the distributions used were $\nu = 2$ (half-gaussian), $\nu = 1$ (exponential), and $\nu = 1/2$ (stretched exponential). The mean of the above distributions is given by V_o. The unit of energy was chosen such that

$$\frac{{\phi_{o}}^2 r_p}{8 \pi^2 \lambda^2 V_o} = 200.$$ (12)

This choice is consistent with experimental observations [3], and ensures that the maximum force exerted by a pin is independent of the range of the pinning potential.

The numerical integration in stage (i) was carried out using a 4^th order Runge-Kutta method [10] with adaptive step size and an initial force resolution of 10^-6 in the above units with subsequent decrements by factors of 10, as needed. Algorithms of the software libraries LINPACK and EISPACK [11] were used in putting together stages (ii) and (iii). The maximum residual force on any particle in a stable configuration was found to be smaller than $\Delta F = {\cal O}(10^{-12})$.

Free boundary conditions were employed, and flux lines leaving the square (leaks) were put back in random locations. Typically, only a few leaks occurred in a given run.

In order to speed up computations further, several additional measures were taken: If a stable configuration was not found after a fixed number of integration steps (5000), a new initial configuration was generated. Initial configurations were such that flux lines were randomly positioned, but no two were allowed to be closer initially than $0.2~\mu$m. These measures were found to reduce the computation time without changing the results significantly. The particle-pin interaction, Eq. (4), was cut off at a distance of 4r_p. Prior to each run, the area was subdivided into a mesh of squares and a look-up table was generated providing each square with the location and strength of pins within the effective interaction distance 4r_p. This eliminates the need to search all pinning sites for each computation of the pinning forces.

Figure 2: Pair correlation function g(r) obtained from averaging over 212 stable flux line configurations on a pinscape, for $\lambda = 0.04$ and r_p = 0.015. The pinning strengths are exponentially distributed. The inset shows a comparison with the liquid theory prediction, Eq. (8).