Activated vortex lattice transition in a superconductor with combined sixfold and twelvefold anisotropic interactions

Magnetic flux enters a type-II superconductor in the form of quantized vortices [1]. In the absence of strong pinning by impurities or thermally driven disordering, the mutual repulsion between vortices causes them to crystallize into an ordered vortex lattice (VL) [2]. In a perfectly isotropic superconductor the VL would be triangular and without a preferred orientation relative to the host lattice [3]; however, real materials always posses some degree of anisotropy that causes an alignment of the VL relative to the crystalline axes, often giving rise to field and/or temperature driven symmetry transitions [4]. Such transitions lead to the formation of competing VL domains with different orientations. The VL can thus serve as a simple model system for domain nucleation and growth, with relevance for a wide range of physical systems including martensitic phase transitions [5], domain switching in ferroelectrics [6], and the skyrmion lattice in chiral magnets [7–9].

For magnetic fields applied along the c axis of superconductors with a hexagonal crystal structure, the VL is empirically found always to have a triangular symmetry, but will often undergo continuous and sometimes even non-monotonic rotations as seen in MgB₂ [10] and UPt₃ [11, 12]. In addition, the VL in these materials exhibit significant metastability when cooled or heated across equilibrium phase transitions in a constant field [11, 13]. To reach the equilibrium configuration the VL must be excited by an external perturbation, which causes vortex motion and allows the system to find the global minimum in the free energy [4]. In MgB₂ the transition from the metastable to the equilibrium state has been studied extensively using small-angle neutron scattering (SANS), revealing an activation barrier that increases as the metastable state is gradually suppressed [14]. This 'work hardening' is speculated to arise from a proliferation of grain boundaries in the VL. Although spatially resolved SANS measurements provide some support for such a scenario [15], a real space study of the VL is desirable.

Spatial information about the VL can be obtained by molecular dynamics (MD) methods in which the vortices are treated as point particles with repulsive interactions. In addition to static VL configurations, MD simulations provide information about the dynamics of a large number of vortices over long times. This technique has been widely used to study how pinning affects the structural [16] and dynamical [17–19] properties of vortex matter. Due, however, to the added complexity introduced by a fully anisotropic interaction potential and the resulting nonradial forces, such studies were limited to systems with isotropic vortex–vortex interactions. Recently we implemented a phenomenological model that allows MD simulations of the VL in the presence of anisotropic interactions [20–22]. This was used to study systems with combined sixfold and twelvefold anisotropies relevant for MgB₂ [21]. By varying the ratio of the six- and twelvefold contributions to the interaction potential, it was possible to reproduce the three different triangular VL phases observed experimentally in this material [10, 13]. The VL phases, denoted F, L and I, differ only in their orientation relative to the MgB₂ crystalline axes. Specifically, F and I phases are oriented along the two high symmetry directions within the hexagonal basal plane ( a along the nearest neighbor direction and $\boldsymbol{a}^*$ along the next nearest neighbor direction), and are connected by an intermediate, continuously rotating L phase [10, 23]. The orientation of the VL phases relative to a and $\boldsymbol{a}^*$ is shown in [23], and the directions of the crystalline axes relative to the interaction potential are indicated in figure 1(a). The simulations also spontaneously produced VL dislocations as well as a range of vacancy and interstitial defects that were analyzed in terms of their energy cost [21].

Zoom In Zoom Out Reset image size

In this work we extend our MD simulations to study the gradual VL transition from a metastable L phase towards the equilibrium F phase. The metastable L phase is obtained by first annealing the system with a vortex–vortex interaction potential that has a twelvefold anisotropy and subsequently changing this to a sixfold anisotropy. Weak pinning is included to stabilize the initial twelvefold symmetric phase, and the temperature is kept low enough that thermal fluctuations alone are insufficient to overcome the activation barrier for a transition to the ground state. We do not observe metastable effects in simulations that does not include pinning, suggesting that at least some weak quenched disorder plays a role in the SANS experiments [13–15]. Here, the quenched disorder will decrease thermal hopping effects so that only the mechanical driving from changing the effective field, as discussed below, will play a role.

In the SANS experiments, the transition to the ground state is driven by small-amplitude magnetic field oscillations, leading to a temporary variation in the vortex density akin to a breathing mode [14]. Such a scenario is not straightforward to implement in the MD simulations, due to the periodic boundary conditions, and we instead employ a periodic variation of the superconducting penetration depth. This changes the effective vortex–vortex interaction strength and thus mimics a field oscillation. As successive penetration depth cycles are applied, the metastable L phase gradually rotates towards the ground state F phase. Initially, the angle between the two degenerate L phase domain orientations decreases as a power law with an exponent that increases with the amplitude of the penetration depth cycles, similar to what was observed in the SANS experiments [14]. For the larger amplitudes there is a cross-over to a logarithmic behavior as the number of penetration depth cycles increases. Here, the VL relaxation occurs in bursts or avalanches that are associated with correlated dislocation annihilation events along the grain boundaries. We discuss how these behaviors compare to two-dimensional (2D) coarsening models [24–27], and show that the $\log(t)$ decay could arise due to coarsening or glassy motion with quenched disorder [28].

Our approach can be applied to other particle-based systems with anisotropic interactions, including skyrmion lattices [29–31] or colloidal particles with anisotropic interactions [32, 33]. Furthermore, it will allow modeling of the dynamics associated with the generation and recombination of dislocations, grain boundaries and domain formation both in the VL [14, 15] and in general [34, 35]. We discuss the relevance of our results to the decay of metastable states in other systems with periodic perturbations, such as granular matter [36, 37], colloidal systems [28], magnetic systems with return point memory [38, 39], and other types of periodically sheared systems that relax to a steady state under ac driving [40–42].

We follow the approach of Olszewski et al and use a phenomenological vortex-vortex interaction potential that includes six- and twelvefold anisotropies [21]:

$$\begin{equation} U\left(r,\theta\right) = A_v K_0\left(r\right) \left[ 1 + K_6 \cos^2 3 \theta + K_{12} \cos^2 6 \theta \right]. \end{equation} \tag{ 1 }$$

Here r and θ are the distance and angle, respectively, between two vortices at positions r _i and r _j , namely $r = |\boldsymbol{r}_i - \boldsymbol{r}_j|$ and $\theta = \tan^{-1}(r_y / r_x)$, with $\boldsymbol{r} = \boldsymbol{r}_i - \boldsymbol{r}_j$, $r_x = \boldsymbol{r} \cdot \boldsymbol{\hat{x}}$, and $r_y = \boldsymbol{r} \cdot \boldsymbol{\hat{y}}$. The radial part of the repulsive interaction potential is given by the amplitude A_v times the zeroth order Bessel function $K_0(r)$ [43]. The parameters K₆ and K₁₂ are the magnitudes of the six- and twelvefold contributions to the potential, with an anisotropy ratio defined by $\kappa = K_6/K_{12}$. When $|\kappa| \gt 4$, the potential is dominated by the six-fold term, with minima, for a fixed r, along either the horizontal ($\kappa \lt -4$) or vertical (κ > 4) direction and at every 60^∘. This is illustrated in figure 1(a). For $|\kappa| \unicode{x2A7D} 4$, there are 12 minima in the interaction potential as shown in figure 1(b), and consequently there are two degenerate triangular VL configurations. The locations of the minima change as κ is varied and are given by

$$\begin{equation} \theta_{\textrm{min}}\left(\kappa\right) = \pm \frac{1}{6} \arccos \left( - \frac{\kappa}{4} \right). \end{equation} \tag{ 2 }$$

Six of the minima rotate clockwise as κ is reduced while the other six rotate counterclockwise, as shown in figure 1(c). Changing the sign of κ is equivalent to a 30^∘ rotation, which exchanges the locations of the minima and maxima.

It should be noted that the interaction in equation (1) is not based on an ab-initio derivation, and, within Ginzburg–Landau (GL) theory, it is strictly speaking not possible to formulate an interaction potential for a single vortex [44]. Instead equation (1) is motivated by an expression for the entire VL free energy obtained from GL theory [45]. As shown by Olszewski et al, this potential reproduces the experimentally observed macroscopic VL phases in MgB₂ [10] and UPt₃ [11, 12] for a suitable choice of parameters K₆ and K₁₂ [21]. It can thus be considered a minimal and conceptually simple phenomenological model for MD simulations of the VL in superconductors with a hexagonal crystal structure.

In the MD simulations the dynamics of each vortex is governed by an overdamped equation of motion:

$$\begin{equation} \eta \frac{\mathrm{d}\boldsymbol{r}_i}{\mathrm{d}t} = \boldsymbol{F}_i^{vv} + \boldsymbol{F}_i^p +\boldsymbol{F}_i^T. \end{equation} \tag{ 3 }$$

Here η is the damping constant which is set equal to unity. The force field from the surrounding vortices is given by $\boldsymbol{F}^{vv} = -\boldsymbol{\nabla}(U) = (-\partial U/\partial x, -\partial U/\partial y)$, yielding components

$$\begin{align} \frac{F_x}{A_v} & = \cos \theta \, K_1\left(r\right) \left[ 1 + K_6 \cos^2 \theta + K_{12} \cos^2 6 \theta \right] - \frac{1}{r} \sin \theta \, K_0\left(r\right) \left[ 3 K_6 \sin 6 \theta + 6 K_{12} \sin 12 \theta \right] \end{align} \tag{ 4a }$$

$$\begin{align} \frac{F_y}{A_v} & = \sin \theta \, K_1\left(r\right) \left[ 1 + K_6 \cos^2 3 \theta + K_{12} \cos^2 6 \theta \right] + \frac{1}{r} \cos \theta \, K_0\left(r\right) \left[ 3 K_6 \sin 6 \theta + 6 K_{12}\sin 12 \theta \right]. \end{align} \tag{ 4b }$$

We consider a system of size L × L with $L = 108\lambda$, and apply periodic boundary conditions in both the x and y directions. Here λ is the London penetration depth. The sample contains $N_v = 5130$ vortices interacting with $N_p = 4000$ pinning sites placed randomly in non-overlapping positions and generating a pinning force $\boldsymbol{F}_i^p$. Each pinning site is modeled as a parabolic trap of range $r_p = 0.25$ that can exert a maximum pinning force of $F_p = 1.15$ on a vortex. The vortex-pin interaction is isotropic, directed towards the center of the pinning site, and given by $\boldsymbol{F}_i^p = F_p \sum_j^{N_p} (\boldsymbol{r}_j^p - \boldsymbol{r}_i) \, \Theta(r_p - |\boldsymbol{r}_i - \boldsymbol{r}_j^p|).$ For both $\boldsymbol{F}_i^{vv}$ and $\boldsymbol{F}_i^p$ the unit of length is λ. The thermal forces $\boldsymbol{F}_i^T$ are modeled by Langevin kicks with the properties $\langle \boldsymbol{F}^T \rangle = 0$ and $\langle \boldsymbol{F}_i^T(t) \boldsymbol{F}_j^T(t^{^{\prime}}) \rangle = 2 \eta k_B T \, \delta_{ij} \, \delta(t-t^{^{\prime}})$, where k_B is the Boltzmann constant. The Langevin kicks are implemented by applying random forces in the x and y directions drawn from a Gaussian distribution with zero mean and a standard deviation of F^T . In our dimensionless units F^T is the strength of the thermal forces and referred to as the temperature of the system. The values of r_p and F_p permit at most one vortex to occupy each pinning site, and were chosen by trial and error in order to stabilize a metastable VL configuration yet still allow the system to be driven towards the equilibrium configuration.

In the 2D MD approximation, the vortices are assumed to be stiff lines spaced far enough apart that their cores do not overlap. This is appropriate as MgB₂ is a three dimensional (i.e. not layered) materials with weak pinning and therefore little or no vortex bending is expected. It should also be noted that MD simulations have been used extensively to model the dynamics of stiff vortices in systems with pinning [46–48].

The results are organized as follows. First, the protocol used to prepare the metastable VL configuration is described, as is how the bulk behavior is characterized by the peak splitting in the structure factor. Second, it is demonstrated how the VL can be driven towards the equilibrium configuration by successive applications of the effective field oscillation, and the dynamics of this process is explored as a function of the oscillation amplitude. Finally, the evolution of VL grain boundaries and energetics is studied at the local scale as the system is driven towards the equilibrium state.

We fix the twelvefold anisotropy amplitude at $K_{12} = 0.0125$ and vary the anisotropy ratio κ by changing K₆, which was previously established to produce reliable results [21]. Varying κ changes the magnitude of the interaction potential experienced by each vortex, which is equivalent to introducing a change in the effective vortex density. To eliminate this variation, the effective magnetic field, defined by a 2D integral of the interaction potential

$$\begin{equation} B_{\mathrm{eff}} \propto \int_0^{\infty} r \, \mathrm{d}r \int_0^{2\pi} \mathrm{d}\theta \; U\left(r,\theta\right) \propto A_v \left( 1 + \frac{K_6}{2} + \frac{K_{12}}{2} \right), \end{equation} \tag{ 5 }$$

is kept at a constant value. Here, $K_6 = 0$ and $A_v = 2.0$ is used as a reference, and the magnitude of the isotropic contribution to the interaction potential is adjusted such that $A_v = 2 (2 + K_{12})/(2 + K_6 +K_{12})$.

3.1. Preparing a metastable VL

The protocol used to prepare the VL in a metastable configuration is shown schematically in figure 2. The system begins with κ = 0 and a temperature $F^T = 4.0$, which is sufficiently high to put it in a molten state. It is then cooled gradually to $F^T = 0.0$ by reducing the temperature by $\Delta F^T = -0.05$ every $2 \times 10^4$ simulation time steps, which is sufficiently slow to ensure that the system reaches an equilibrium state [21]. The anisotropy ratio is then increased to κ = 5 in increments of $\Delta \kappa = 0.25$ every $2 \times 10^4$ simulation time steps, while keeping $F^T = 0$. Finally, the temperature is increased to $F^T = 0.35$ using the same rate as for the initial cooling. The final temperature increase provides the system with the necessary thermal energy to overcome activation barriers and be driven towards the equilibrium configuration by the effective field oscillations described in section 3.2.

Zoom In Zoom Out Reset image size

The global VL configuration is characterized by computing the structure factor, $S(\boldsymbol{k}) \propto |\sum^{N_v}_i \exp(-i \boldsymbol{k} \cdot \boldsymbol{r}_i)|^2$. Figure 3(a) shows $S(\mathbf{k})$ at κ = 0 following the initial annealing to $F^T = 0$. Here the intensity is concentrated in 12 (Bragg) peaks lying on a circle, indicating the presence of VL domains that are oriented along one of two different directions corresponding to the minima in the interaction potential in equation (1). The radius agrees with the calculated value for a triangular lattice, $k_0 = (2/\sqrt{3})^{1/2} \; 2\pi (\sqrt{5130}/108\lambda) = 0.71 (2\pi/\lambda)$. Figure 3(b) shows the azimuthal dependence of the structure factor folded into one 60^∘ segment, fitted by a two-peak Voigt function. Voigt fits (a convolution of a Gaussian function and a Lorentzian) are used throughout this work because they were found to provide better fits to the azimuthal intensity distribution than simple Gaussian or Lorentzian. The peak splitting, $\Delta \theta$, is defined as the difference between the two fitted peak positions.

Zoom In Zoom Out Reset image size

A metastable VL is established when κ is increased at zero temperature, since without thermal fluctuations the system is no longer able to evolve enough to remain in the changing equilibrium state. This is illustrated in the plot of $\Delta \theta$ in figure 4(a). Although $\Delta \theta$ decreases it always remains higher than the equilibrium value $\Delta\theta_{\textrm{eq}} = 60^{\circ} - 2 \theta_{\textrm{min}}(\kappa)$, indicating that the system becomes trapped in a metastable state. The measured splitting angle has a kink near κ = 4 where $\Delta\theta_{\textrm{eq}}$ reaches zero. Once the final value of κ = 5 has been reached, the system is allowed to evolve for an additional $2 \times 10^4$ simulation time steps. This only gives rise to a modest additional reduction of the splitting angle, which stabilizes at a value of $\Delta \theta \sim 20^{\circ}$. During the subsequent temperature increase at κ = 5, $\Delta \theta$ remains largely unchanged as shown in figure 4(b), indicating that the thermal fluctuations alone are insufficient to permit the system to escape the metastable state.

Zoom In Zoom Out Reset image size

Real space information about the VL is obtained from Voronoi constructions of the vortex positions. Within each Voronoi cell, all points are closer to a given vortex than to any other. Furthermore, the number of edges yields the coordination of the vortex within the cell. The Voronoi constructions provide local information about domain formation, VL defects and grain boundaries. Figure 5(a) shows the Voronoi construction for an equilibrium VL at κ = 0 and $F^T = 0$, corresponding to figure 3. Here the system breaks up into domains with two different orientations. Within each triangular domain, the vortices are sixfold coordinated, but due to the twelvefold anisotropy of the vortex–vortex interaction potential, there are two degenerate orientations. The grain boundaries between domains are decorated by dislocations consisting of fivefold and sevenfold coordinated vortices.

Zoom In Zoom Out Reset image size

Figure 5(b) shows the Voronoi construction plotted after the anisotropy ratio has been increased to κ = 5 at zero temperature. Although the number of dislocations has decreased, the system still shows clear domain formation with easily discernible grain boundaries. Some grain boundaries have shifted relative those in figure 5(a) and indicated by the black lines. Furthermore, the VL domains have rotated towards each other as indicated by a reduced peak splitting in the structure factor. In other words, the relaxation towards the equilibrium configuration occurs primarily through a continuous rotation of the metastable VL domains rather than the nucleation and growth of separate ground state domains. The ground state for this κ would be a single domain containing sixfold coordinated vortices. Heating the system to $F^T = 0.35$ at κ = 5 causes some additional annihilation of dislocations and grain boundary motion, but the VL remains in a metastable configuration as seen from the Voronoi construction and structure factor in figure 5(c). This confirms that thermal effects at this temperature are not strong enough to permit vortices to hop over the activation barrier between the metastable equilibrium states.

3.2. Effective field oscillations

In the SANS experiments, the transition from the metastable to the equilibrium state is induced by the application of magnetic field oscillations, which change the number of vortices in the sample [14]. Changing the number of vortices is not possible in the MD simulations with periodic boundary conditions, and instead an effective change in the vortex density is modeled by varying the penetration depth $\lambda = \lambda_0 + \Delta \lambda$. This corresponds to a change in both the system size $L = 108 \lambda$ and the range of the vortex–vortex interaction in equations (4a ) and (4b ) as r is measured in units of λ. As such, it is equivalent to a change of the magnetic field

$$\begin{equation} B = \frac{N_v \Phi_0}{\left(L \lambda\right)^2} = \frac{B_0}{\left(1 + \Delta \lambda/\lambda_0\right)^2}, \end{equation} \tag{ 6 }$$

where Φ₀ is the flux quantum and B₀ is the field for $\Delta \lambda = 0$. The corresponding amplitude of the field variation, $\Delta B = B - B_0$, is given by

$$\begin{equation} \frac{\Delta B}{B_0} = \frac{1}{\left(1 + \Delta \lambda/\lambda_0\right)^2} - 1. \end{equation} \tag{ 7 }$$

Changing the range of the vortex–vortex interactions by varying λ modifies the force balance between the unevenly spaced vortices near the grain boundaries. This creates mechanical instabilities that allow for vortex jumps or rearrangements to occur. It is important to note that whereas varying λ in an experiment is associated with a change in temperature, all the effective field oscillations are performed with a constant F^T . An alternative to the effective field oscillations would have been to drive the system would be with an applied current. Such an approach was previously used to study the dynamical annealing of disordered vortex states [49].

Each effective field oscillation is performed by changing λ in steps of 0.05 every 2000 simulation time steps until the desired value of $\Delta \lambda$ is reached. The cycle is completed by changing the penetration depth back to the original value $\lambda = \lambda_0$ at the same rate. It was verified that this rate is slow enough that the results do not change if the frequency is further reduced. Values of $\Delta \lambda$ between $-0.1 \lambda_0$ and $-0.5 \lambda_0$ were studied, corresponding to $\Delta B/B_0$ between 0.23 and 3.0. Figure 5(d) shows the Voronoi construction and structure factor after 1000 penetration depth oscillations were applied to the mestastable state in figure 5(c). The number of dislocations are reduced and the VL domains have rotated closer to the equilibrium orientation as indicated by the smaller structure factor peak splitting. Furthermore, more grain boundary motion is evident and also some fracturing of the VL into smaller domains. This shows that the penetration depth oscillations are able to move the grain boundaries and thus mimic the behavior observed in the SANS experiments [14]. Notably, while all simulations were performed with a different configuration of the pin locations the system evolve toward the same equilibrium state.

One difference between the penetration depth oscillations, used in the MD simulations, and the experimental magnetic field oscillations is that the latter produces a small gradient in the vortex density perpendicular to the sample boundary. This gradient gives rise to a current, which, in turn, causes vortex motion. In comparison, no current is applied in the MD simulations since the simulations use periodic boundary conditions and there is no sample boundary. It is therefore not possible to directly compare the values of the effective field oscillations used in the MD simulations and the much smaller values used in the SANS experiments ($\Delta B/B_0 = 0.001 - 0.003$) [14].

3.3. Metastable state relaxation dynamics

To study the relaxation dynamics of the metastable VL state, MD simulations were performed in which the system was subjected to successive penetration depth oscillations. Simulations were performed for $\Delta B/B_0 = 0, 0.2346, 0.5625, 1.0408, 1.7778,$ and 3.0. Ten simulation runs were performed for each value of $\Delta B/B_0$. Figure 6(a) shows $\Delta \theta$ determined from the structure factor peak splitting, averaged over the ten simulation runs, versus the cycle count n + 1, where n is the number of penetration depth oscillations that have been applied. Using n + 1 as the abscissa allows the inclusion of the initial metastable state (n = 0).

Zoom In Zoom Out Reset image size

For non-zero $\Delta B/B_0$ an initial rapid drop of $\Delta \theta$ is observed, with a magnitude that increases with the amplitude of the effective field oscillation. The data for $n \unicode{x2A7D} 25$ is well fitted by the functional form $\Delta \theta = a(n+1)^\gamma + b$, as shown by the full lines in figure 6(a). A similar form was also used to fit the experimental SANS data [14]. Fitted values of γ as a function of $\Delta B/B_0$ are shown in figure 6(b). For $\Delta B/B_0 = 0$ (no oscillations), $\Delta \theta$ decreases very slowly due only to thermal effects at $F^T = 0.35$ and $\gamma\sim 0$ although the uncertainty is large. For the lowest non-zero $\Delta B/B_0$, γ is still very small but the uncertainty is significantly reduced. As the amplitude of the field oscillations is increased, the magnitude of γ also grows until, for the highest value of $\Delta B/B_0$, $\gamma = -0.5$. The overall behavior of γ with increasing $\Delta B/B_0$ is similar to that found in the SANS experiments [14], where exponents ranging from $\gamma = -0.051$ to $\gamma = -0.22$ were obtained as shown in figure 6(b). It is likely that the largest values of $\Delta B/B_0$ used in the MD simulations exceed the range explored in the SANS in the experiments, but as discussed a direct comparison is not possible. As evident from figure 6(a), the power law fits of $\Delta \theta$ describe the data well for $n \lesssim 300$; however, at large n they do not perform as well. Rather, the system exhibits a logarithmic behavior at large n for all values of $\Delta B/B_0$.

To gain a better understanding of the longer time behaviors, figure 7(a) shows $\Delta \theta$ for $\Delta B/B_0 = 3.0$ versus n + 1 on a linear scale for a single simulation run. After the initial smooth decrease, a series of sharp downward jumps or avalanches appears for n > 100. Figure 7(b) shows the corresponding fraction of sixfold coordinated vortices P₆ versus n + 1. From this it is evident that the downward jumps in $\Delta \theta$ are accompanied by sudden increases of P₆, meaning that the avalanches are due to the annihilation of groups of dislocations. This is different from the behavior for $n+1 \unicode{x2A7D} 15$ where P₆ undergoes a smooth continuous increase.

Zoom In Zoom Out Reset image size

A more direct demonstration that the jumps in the splitting angle occur by means of avalanches, in which multiple dislocation pairs are annihilated, is shown in figure 8. This shows a portion of the sample immediately before and after the jump highlighted in the insets of figure 7. In the Voronoi construction prior to the jump in figure 8(a), eight 5–7 dislocations are located near a VL grain boundary within the circled region. After the jump, in figure 8(c), about four of the dislocations have annihilated. This is also reflected in the local VL energy shown in figures 8(b) and (d). After the jump, the overall energy of the system has dropped and become more uniform.

Zoom In Zoom Out Reset image size

In summary, the dynamics can be broken into two parts with a fast initial power law dependence on n followed by a slower universal logarithmic dependence dominated by bursts or avalanches.

The decay of a metastable state studied here and in experiments could have implications for many other non-thermal systems in metastable states. One example is granular matter confined in a vertical container. After the grains are initially poured in to the container, they have a fixed volume fraction; however, if the container is subjected to periodic vertical shaking or tapping, the grains become more compressed and the volume fraction V varies as $V \propto 1/(a + b \log(n))$ where n is the number of taps applied [36]. In 2D systems, there have been several studies examining the coarsening dynamics of topological defects in the case where the equilibrium system would be defect free. Over time, the density of defects decays according to a power law with an exponent in the range $1/4$ to $1/3$ [24–27]. The values of the exponents are consistent with some of the values we obtain at shorter times. Studies of systems that form 2D crystals in the absence of quenched disorder show that when quenched disorder is added, numerous defects appear in the lattice that form grain boundaries at zero temperature; however, if an ac driving is applied that shakes the particles back and forth over the substrate, the defects annihilate over time at a rate that becomes proportional to $\ln(t)$ at longer times [28, 40, 50]. If quenched disorder is present, there are generally some defects that remain permanently trapped.

One of the closest systems outside of superconducting vortices to which we can compare our results is the mechanical annealing of grain boundaries in 2D colloidal crystals, where experiments found that the density of grain boundaries decays as a power law as a function of time [50]. It may be possible that the coarsening dynamics has a faster time scale than the glassy dynamics which are typically associated with a $\ln(t)$ dependence, so that for the first 1–100 cycles the system behaves as if it is coarsening, but at longer times it reverts to $\ln(t)$ behavior when the dynamics becomes dominated by hopping or avalanches.

Although our results are not accurate enough to distinguish between the different types of coarsening, at early times they do show that the basic microscopic ingredients for the model we have presented capture the main features found in the experiments. Additionally, our results agree with the vortex–vortex interactions observed in the experiments, where a metastable twelvefold symmetric state is quenched into a sixfold symmetric state and the ac driving progressively reduces the splitting angle.

It could also be possible that in some cases an ordered system in the presence of quenched disorder would show increased defect formation under an applied drive because some particles become trapped while other particles shear past them. This also suggests that there could be an optimal ac amplitude or frequency for creating the maximum amount of dynamic annealing. In addition, introduction of quenched disorder can cause the system to remember the metastable states. If a cyclic driving current were applied, the vortices would move in a periodic manner, and it would be possible to examine whether they return to the same positions at the end of each cycle in a form of return point memory or random organization, which has been studied in numerous periodically driven systems [38, 39, 41, 42].

The relaxation of a metastable vortex state was studied numerically for a system quenched from a twelvefold to a sixfold anisotropy of the vortex–vortex interaction. The system is initially in an equilibrium configuration for the twelvefold anisotropy, consisting of triangular VL domains with two different orientations and separated by grain boundaries decorated by dislocations. In comparison, the ground state for the sixfold anisotropy is a single triangular VL domain. The bulk VL configuration is characterized by the structure factor. Following the quench from the twelvefold to the sixfold anisotropy, the initial twelve equally spaced structure factor peaks form six pairs with a decreasing splitting. However, the presence of grain boundaries and weak pinning traps the system in a metastable state with a non-zero splitting angle.

The system is driven towards the ground state by a periodic variation of the penetration depth, which mimics the magnetic field oscillations used in SANS experiments [14]. This causes a gradual decrease of the structure factor peak splitting with successive applications of the penetration depth cycles. Initially, the decrease follows a power law with an exponent that increases with the amplitude of the penetration depth oscillations, in agreement with the SANS experiments performed in a similar regime [14]. As the number of penetration oscillations grows, the decay of the peak splitting becomes logarithmic with a slope that is independent of the amplitude. This indicates a transition from annihilation of individual dislocation pairs at short times to coarsening dynamics at longer times, where the splitting angle decreases in bursts or avalanches that are correlated with the annihilation of larger groups of dislocations.

We are grateful to E Roe, M W Olszewski and D Spulber for assistance with the molecular dynamics simulations and analysis of the results. This research was supported in part by the Notre Dame Center for Research Computing. Work at the University of Notre Dame (DM, MRE: MD simulations, data analysis) was supported by the US Department of Energy, Office of Basic Energy Sciences, under Award No. DE-SC0005051. Part of this work (CR, CJOR: code development) was carried out under support by the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U. S. Department of Energy (Contract No. 892333218NCA000001).

All data that support the findings of this study are included within the article (and any supplementary files).