# Transport Anisotropy as a Probe of the Interstitial Vortex State in Superconductors with Artificial Pinning Arrays

## Abstract

We show using simulations that when interstitial vortices are present in superconductors with periodic pinning arrays, the transport in two perpendicular directions can be anisotropic. The degree of the anisotropy varies as a function of field due to the fact that the interstitial vortex lattice has distinct orderings at different matching fields. The anisotropy is most pronounced at the matching fields but persists at incommensurate fields, and it is most prominent for triangular, honeycomb, and kagomé pinning arrays. Square pinning arrays can also show anisotropic transport at certain fields in spite of the fact that the perpendicular directions of the square pinning array are identical. We show that the anisotropy results from distinct vortex dynamical states and that although the critical depinning force may be lower in one direction, the vortex velocity above depinning may also be lower in the same direction for ranges of external drives where both directions are depinned. For honeycomb and kagomé pinning arrays, the anisotropy can show multiple reversals as a function of field. We argue that when the pinning sites can be multiply occupied such that no interstitial vortices are present, the anisotropy is strongly reduced or absent.

###### pacs:

74.25.Qt## I Introduction

Vortices in superconductors interacting with artificial arrays of periodic pinning exhibit a wide range of commensurability and dynamical effects that can be observed readily in critical current, transport, and other bulk measurements (1); (2); (3); (4); (5); (6); (7); (8); (9); (10). Advances in lithography techniques permit the creation of pinning arrays in which the size, shape, and composition of the individual pinning sites and the global geometry can be well controlled (2); (3); (4); (5); (6); (7); (8); (9); (10); (11); (12); (13); (14). Commensurability effects in these systems occur when the number of vortices equals an integer multiple of the number of pinning sites, resulting in peaks or anomalies in bulk measurements as a function of field. At the first matching field, there is one vortex per pinning site, and as the field is further increased, additional vortices can be located either at the pinning sites in the form of pinned multi-quanta vortices (3); (4); (15); (16); (17); (18); (19), or in the interstitial regions between the pinning sites. The interstitial vortices can be effectively pinned by the repulsive interactions from the vortices at the pinning sites, which create a caging potential (6); (5); (20); (21); (22); (19); (23); (24). It is also possible for mixed vortex pinning to occur in which the first few matching fields have only pinned multi-quanta vortices until the pinning sites are saturated, while for higher matching fields the additional vortices are located in the interstitial regions (3); (6); (13); (15); (16); (19). Conversely, it is also possible that interstitial vortices appear at the lower matching fields, but that as the vortex-vortex interactions increase at higher matching fields, multi-quanta vortices will begin to form at the pinning sites (19).

Interstitial vortex lattice crystals in square periodic pinning arrays have been observed directly with Lorentz microscopy, which revealed that there are several distinct types of interstitial vortex structures that have symmetries different from that of the triangular vortex lattice (6). The same types of vortex structures have been produced in simulations of square pinning arrays (20); (21), while simulations have also shown that similar vortex structures can form in triangular (20), rectangular (25); (24), honeycomb (26), and kagomé pinning arrays (26); (27). Other numerical works indicated that a rich variety of composite lattices with multiple and interstitial vortex configurations are possible (22); (19) and that new types of interstitial vortex configurations can occur for arrays of antipinning sites (28).

Vortex imaging experiments provide direct evidence for both multi-quanta vortex pinning and the formation of ordered interstitial vortex lattice structures (15); (16). Anomalies at matching fields found in bulk measurements occur for both multi-quanta vortex pinning and interstitial vortex pinning, so without direct imaging it can be difficult to determine whether multi-quanta or interstitial vortex pinning is occurring (18). In some cases, the presence of interstitial vortices can only be inferred from the shapes and characteristics of the current voltage curves or from phase locking experiments (5); (29). It would be highly desirable to identify additional clear signatures in transport measurements that can distinguish between interstitial vortex pinning and multi-quanta vortex pinning and that can also reveal the types of vortex lattice symmetries that are present.

An anisotropic response was recently measured for the critical current applied in two perpendicular directions to a triangular pinning array in recent experiments and simulations (30); (31). The experiments were performed on several different samples and the same anisotropic response appeared in each one, while the anisotropy observed in the simulations agreed with that seen in the experiments, suggesting that the behavior is due to distinct intrinsic features of the vortex dynamics. The simulations show that the vortex flow patterns are different for the two directions of applied current, which could account for the anisotropic response. The anisotropy is particularly pronounced at the second matching field but is absent at the first matching field, which suggests that depinning of the interstitial vortices is responsible for the anisotropy. Interestingly, for lower temperatures the experiments showed that the critical current anisotropy vanished at both the first and second matching fields but persisted at the third matching field, suggesting that multi-quanta vortex pinning occurs at the second matching field at low temperatures. These results indicate that the presence of anisotropy can be a useful way to probe the interstitial vortex state and the type of vortex ordering that occurs at different matching fields.

The experiments and simulations of Refs. (30); (31) only examined a triangular pinning array up to the third matching field. In this work we study anisotropic transport for a much wider range of fields and system parameters for triangular, square, honeycomb, and kagomé pinning arrays. For the triangular array we find a variety of distinct anisotropic behaviors due to the differing symmetries of the interstitial vortex lattice at different matching fields. For example, at certain matching fields the vortex lattice is disordered and the anisotropy vanishes. For honeycomb and kagomé pinning arrays, an even richer anisotropic behavior occurs due to the formation of vortex molecular crystal states (26) with additional rotational degrees of freedom, resulting in a series of reversals in the anisotropy as a function of field. Remarkably, we find that it is also possible for square pinning arrays to show anisotropic transport even though the two perpendicular directions of the pinning lattice are identical. This occurs at certain matching fields where a triangular vortex structure forms, such that one driving direction is oriented with the easy-shear direction of the vortex lattice. Anisotropic transport measurements for vortices in periodic pinning arrays have already been been shown to be experimentally feasible, and experiments have been performed in which the current is injected in two directions for samples with rectangular pinning arrays or asymmetric pinning shapes, revealing anisotropic depinning thresholds (32); (33); (34). Although our work is focused on superconducting vortices, our results should be general to the class of systems of particles interacting with periodic substrates where both pinned and interstitial particles are present. Examples of this type of system include colloidal particles in periodic pinning (35); (36) and charged metallic balls (37).

## Ii Simulation and System

We consider a two-dimensional system with periodic boundary conditions in the and -directions of size . The magnetic field is applied out of the plane in the -direction, while vortices and pinning sites are placed within the system for a vortex density of and a pinning density of . The matching field is defined as the field at which the number of vortices equals the number of pinning sites, . In Fig. 1 we show representative examples of the triangular, square, honeycomb, and kagomé pinning geometries used in this work, and indicate the and directions along which current is applied. The initial vortex positions are prepared by simulated annealing with no applied drive, and then the vortex velocities are measured in the presence of a driving force that is applied in the direction. The driving force corresponds to the Lorentz force generated by an applied current, while the vortex velocities are proportional to the voltage response that would be measured experimentally. We repeat the simulation from the same initial vortex positions with the driving force applied in the direction, and compare the velocity response and and the critical depinning force and for the two driving directions.

The time evolution of the vortex dynamics is governed by integrating coupled overdamped equations of motion. The equation of motion for a single vortex at position is given by

(1) |

Here the damping constant , where is the sample thickness, is the coherence length, is the normal-state resistivity, and is the elementary flux quantum. The vortex-vortex interaction force is

(2) |

where is the modified Bessel function, is the London penetration depth, the unit of force is , , and . The vortex-vortex interaction force falls off sufficiently rapidly that a cutoff can be placed at . Use of a longer cutoff of produces identical results. An additional short range cutoff is placed at to avoid a divergence in the force. The pinning sites are modeled as attractive parabolic wells of radius and strength with . Here is the location of pinning site , , , and is the Heaviside step function. The pinning sites are arranged in a triangular, square, honeycomb, or kagomé array. The external drive or is a constant force that is uniformly applied to all of the vortices. The thermal force is used during the simulated annealing procedure and has the following properties: and . We decrement the temperature by every 1000 simulation time steps. After the initialization, the applied drive is imposed in increments of every simulation time steps. The velocity-force curves are obtained by averaging the velocity every simulation time steps: , where , . Here . The critical depinning forces in the and directions, and , are determined by the criterion .

## Iii Anisotropy in Triangular Pinning Arrays

We first consider the depinning forces for the and -directions in the triangular pinning lattice illustrated in Fig. 1(a) with , , and . In Fig. 2 we plot and vs along with a detail of the region from . In general at most of the commensurate fields. The anisotropy at incommensurate fields is most pronounced for , as shown in the inset of Fig. 2. In previous simulations a similar anisotropy was observed for this range of fields (31). For both and are small, so it is difficult to determine whether the depinning is anisotropic at incommensurate fields. At the matching fields , 9, and 12, where the depinning forces are high, and can be measured accurately. In Fig. 3 we plot representative velocity-force curves for driving in both the and directions at , 3, 4, and 5. There is a clear anisotropy in the depinning force with at , 3, and 4 in Figs. 3(a,b,c), while at in Fig. 3(d), the depinning is isotropic. This can be seen more clearly in the plot of the anisotropy at the matching fields, shown in Fig. 4. The anisotropy is largest for 4, and 9, weaker for and 12, and essentially absent at , 5, 6, 7, 8, 10, and 11. Corresponding to this, there are no peaks in the depinning force at , 6, 8, 10, and 11 in Fig. 2. At the fields with isotropic depinning , the overall vortex lattice is disordered, while at the other matching fields where anisotropic depinning occurs, ordered vortex lattices form. Previous numerical work for triangular pinning arrays at fields up to showed the same features in the critical depinning force as well as the existence of disordered lattices at the matching fields , 6, and 8 (20). From a geometric construction, it can be shown that a triangular vortex lattice can be placed on a triangular pinning lattice at the integer matching fields , where and are integers (20). This predicts the formation of triangular vortex lattices at fields with , 3, 4, 7, 9, and 12, in agreement with our observation of a peak in at each of these fields. The geometric construction does not predict the formation of a triangular lattice at ; however, a strong matching peak appears at this field both in Fig. 2 and in the previous work (20). The peak at occurs due to the formation of an ordered honeycomb vortex lattice structure, rather than a triangular vortex lattice. In general, we expect to find a peak in the critical current at fields where a triangular or other ordered vortex lattice structure forms. Figures 2 and 4 also show that although there is a peak in at , where a triangular vortex lattice forms, there is no anisotropy for this field and .

In order to explain the different degrees of anisotropy that appear at different fields, we analyze the vortex positions and trajectories for driving in the and directions. Figure 5(a) shows the ordered honeycomb vortex lattice structure that forms at for on a triangular pinning lattice. When driven in the -direction, the vortices can easily form slightly undulating channels of flow that pass between the filled pinning sites, as illustrated in Fig. 5(b). For driving in the -direction, the pinned vortices act as barriers that prevent the formation of simple flow channels, and the interstitial vortices can only move by making significant excursions into the -direction, as illustrated in Fig. 5(c) for a drive just above depinning. As a result, a larger external force is required to cause depinning in the -direction, and . The onset of motion in the -direction is very sharp, as seen by the jump in the velocity-force curve in Fig. 3(a), and after a brief initial period of disordered motion, the vortex flow quickly organizes into the pattern shown in Fig. 5(c). We observe similar motions just above depinning in the and -directions at the incommensurate fields for ; however, the presence of vacancies in the honeycomb vortex lattice causes certain rows of interstitial vortices to depin at a lower value of than at the commensurate fields. The velocity-force curves at shown in Fig. 3(a) indicate that , meaning that the vortex velocity has a lower slope as a function of increasing drive in the -direction than in the -direction, even though the same number of vortices are moving for either direction of drive. This difference is a result of the fact that the vortex motion for -direction driving has much larger excursions transverse to the driving direction since the interstitial vortices must go out of their way to pass around the pinned vortices.

At on a triangular pinning lattice, the overall vortex lattice ordering at is triangular, as shown in Fig. 5(d). The vortices move in slightly winding channels upon application of a drive in the -direction, as seen in Fig. 5(e); these structures are similar to the channels that form at in Fig. 5(b). For depinning in the -direction at , the pinned vortices again create a barrier to the formation of simple channels of interstitial vortex flow; however, unlike the ordered and strongly winding channels that form at , for the vortices move in a disordered fashion just above depinning, as illustrated in Fig. 5(f). An ordered flow state similar to that shown in Fig. 5(c) does not occur for until a much higher value of is applied. This result indicates that in addition to the anisotropy in , there are also strongly different vortex velocity fluctuation characteristics for driving along the and -directions at . Narrow band noise signatures should arise from the synchronized vortex motion that occurs for driving along the -direction, while for driving in the -direction the velocities are more random and a broad band noise signature should appear. In mode-locking experiments, where an external ac drive is imposed along with an applied dc drive, Shapiro type steps would appear for the ordered motion along the -direction, while Shapiro steps would be absent for driving along the -direction. Shapiro steps could be induced for both driving directions at since both vortex flow patterns show synchronized motion; however, some of the characteristics of the Shapiro steps might differ for the two directions since the meandering of the vortices is distinct in the and -directions.

A triangular vortex lattice with a single row of interstitial vortices between each row of pinning sites forms for at , as shown in Fig. 6(a). Just above the depinning transition for driving in the -direction, a portion of the interstitial vortices depin into the one-dimensional channels illustrated in Fig. 6(b). Here, one-third of the interstitial vortices remain immobile in the interstitial regions between pinning sites; these immobile interstitial vortices depin at a higher value of that is outside the range of driving forces shown in Fig. 3(c). For driving in the -direction, all of the interstitial vortices depin simultaneously and flow in winding channels around the occupied pinning sites, as shown in Fig. 6(c). The velocity-force curves shown in Fig. 3(c) indicate that although for , . As a result, the two velocity-force curves cross near . The slope is steeper for the -direction driving since all the interstitial vortices take part in the motion, whereas for the -direction driving, only of the interstitial vortices are moving.

At and the vortex lattice is disordered, as seen in Fig. 6(d), and the depinning is isotropic. Figure 6(d,e) shows that the same type of disordered vortex motion occurs for depinning in both the and -directions. In general, we observe disordered flow states in both driving directions at the other matching fields where a disordered vortex lattice forms, including , 8, 10, and 11.

At , the depinning thresholds and velocity-force curves are very similar to those found for since the triangular vortex lattices that form at these two fields have the same orientation. In Fig. 7(a), the vortex configuration at for contains two rows of interstitial vortices between adjacent pairs of pinning rows, whereas at , Fig. 6(a) shows that there is only one row of interstitial vortices between each pair of pinning rows. Just above the depinning transition for -direction driving at , Fig. 7(b) indicates that the two rows of interstitial vortices flow in one-dimensional channels between the pinning rows while two interstitial vortices remain trapped behind every pinning site so that of the interstitial vortices are moving. This is similar to the -direction depinning that occurs for , where a single row of interstitial vortices flows between each pair of pinning rows and a single interstitial vortex is trapped behind each pinning site. Just above the depinning transition for -direction driving at , shown in Fig. 7(c), all of the interstitial vortices are depinned and a combination of ordered and disordered vortex flow occurs. The velocity-force curves for driving in the and directions show a similar crossing at as that illustrated in Fig. 3(c) for .

For , the vortex configurations and depinning dynamics are similar to those seen for . In Fig. 7(d) the vortex configuration at consists of a triangular lattice that is aligned with the -direction. Since there is now an interstitial column of vortices that could depin and flow between adjacent columns of pinning sites, while there are no straight rows of interstitial vortices, it might be expected that for this field. Instead, Fig. 4 shows that , so the depinning is still easier in the -direction than in the -direction; however, the depinning anisotropy is much smaller than that which appears at . Vortex motion just above depinning at for driving in the -direction occurs in the form of three ordered winding rows passing between each pair of pinning rows, as shown in Fig. 7(e). This is similar to the motion at shown in Fig. 5(e) where one winding row of interstitial vortices moves between the pinning rows. For depinning in the -direction at , Fig. 7(f) shows that the vortex trajectories are disordered in a manner similar to that found for -direction depinning at , as seen in Fig. 5(f).

At , a peak in occurs as shown in Fig. 2, but there is almost no anisotropy in the depinning thresholds, as seen in Fig. 4. Figure 8(a) illustrates that at this field for , a triangular vortex lattice forms which is not aligned with either the or directions, unlike the configurations found at 4, 9, and , but is at an angle to the -axis. The vortices flow in ordered patterns just above depinning in both the and directions, as shown in Fig. 8(b,c). Since a portion of the interstitial vortices remain pinned just above depinning for driving in the -direction, whereas all of the interstitial vortices are flowing for driving in the -direction, at . The absence of the anisotropy at this field is likely due to the fact that the main symmetry axis of the triangular vortex lattice is not aligned with either the or directions, as is the case at the other matching fields which show anisotropy.

An interesting feature that we observe at which does not occur at the other matching fields we have investigated is a switching dynamics that can be induced within the pinned phase. This is illustrated in Fig. 9 where in the initial state, shown in Fig. 9(a), the vortex lattice is tilted at to the -axis. The ground state is twofold degenerate since it could also have been aligned at to the -axis. When a small driving force is applied to the ground state shown in Fig. 9(a) along to the -axis, a structural rearrangement of the vortices occurs accompanied by a domain wall that passes through the system from right to left, as shown in Fig. 9(b,c). The vortices move by approximately one lattice constant in the -direction as the domain wall sweeps past, and in the final state the vortex lattice is tilted at to the -axis. The system can be switched back to if an external drive is applied along to the -axis. We expect that for higher matching fields beyond the fields that we consider here, a switching effect will be present for matching fields containing two or more degenerate ground states where the vortex lattice could be arranged in several possible ways. The application of an external drive will lower the energy of one of the orientations.

### iii.1 Effect of Pinning Strength

We next study the effect of varying the pinning strength on the anisotropy at the different fields. In Fig. 10(a) we plot and versus for a sample with a triangular pinning lattice at , and in Fig. 10(b) we show the corresponding anisotropy ratio . Since the pinning induces the honeycomb vortex lattice structure shown in Fig. 11(a) at this field, as decreases the vortices shift into configuration that is closer to a triangular lattice. To accommodate this shift, some of the pinning sites are vacated with decreasing . For , the honeycomb lattice (H) structure is stabilized, all of the pinning sites are occupied, the anisotropy is fixed near , and the depinning thresholds and do not change significantly with . For , the pinning is not strong enough to stabilize the honeycomb lattice and a partially pinned (PP) lattice forms in which only some of the pinning sites are occupied, as illustrated in Fig. 11(b) for . The depinning in this regime is still plastic with the interstitial vortices depinning first. There is a strong enhancement of the anisotropy in the PP phase, with reaching values as large as 5. For , the vortices form a distorted triangular (DT) lattice illustrated in Fig. 11(c) which becomes increasingly triangular with decreasing . In the DT phase, only a small fraction of the pinning sites are occupied. Here the depinning transition is elastic and both the interstitial and pinned vortices depin simultaneously. The anisotropy of the depinning is lost and the critical depinning forces are isotropic in the DT phase. In general, and decrease with decreasing ; however, near the transition between the PP and DT phases, increases with decreasing as the depinning changes from plastic to elastic. We believe that this effect is similar to the peak in the critical depinning force observed for periodic (38) and random arrays (39) above the first matching field when the vortex-vortex interaction strength is varied. For periodic pinning, Ref. (38) illustrated that the depinning force decreases with increasing vortex-vortex interaction strength since the system depins elastically and all of the vortices begin to move at the same time. For weak vortex-vortex interactions, the interstitial vortices depin easily and flow plastically past the vortices at the pinning sites, so the depinning force decreases with decreasing vortex-vortex interaction strength in this regime. Between these two extremes, a peak in the depinning force occurs. A similar effect is observed in vortex systems with very dilute random pinning arrays (39). In Fig. 10(a) the vortex-vortex interaction strength is fixed; however, there is still a transition from elastic to plastic depinning. The peak appears only for -direction driving, and this may be due to the existence of an easy shear mode in the -direction which is absent for -direction driving.

For , Fig. 10(c) shows that and saturate at large , while Fig. 10(d) indicates that the anisotropy ratio also saturates at . For , the depinning threshold decreases rapidly with decreasing as the system enters the elastic depinning regime. We note that since the vortex lattice at is triangular, there is no elastic energy cost for occupying the pinning sites and one-third of the vortices will always be located at the pinning sites for arbitrarily small , unlike the situation at . There is a very small peak in near the plastic-elastic depinning transition.

Figure 12(a,b) shows the depinning thresholds and as a function of for along with the anisotropy ratio . Both and saturate with increasing , and the anisotropy for is fixed at . For , the depinning threshold decreases rapidly with decreasing , and at the same time there is a drop in the anisotropy as the system passes from the plastic to the elastic depinning regime. In Fig. 12(c,d) we illustrate , , and for , where the vortex lattice lattice is disordered. Here the depinning is isotropic for all values of . A peak in the depinning thresholds occurs for both directions of driving near at the transition between elastic and plastic depinning.

We observe the same general trends for , including a saturation of the anisotropy with increasing and a crossover from plastic to elastic depinning at low . For higher values of and , multiple vortices are trapped at each pinning site. In this case, the vortices form vortex molecular crystals within the pinning sites (40), which changes both the depinning threshold and the anisotropy. In general, the anisotropy is reduced when multiple vortex pinning is present. This means that the anisotropy is only observable in a regime where the pinning is strong enough to allow for plastic depinning or channeling of vortices between the pinning sites, but not strong enough to permit multiple vortices to occupy each pinning site.

### iii.2 Effect of Changing

We next examine the depinning thresholds and anisotropy at different matching fields for the triangular pinning lattice with fixed but with increasing , achieved by increasing the density of pinning sites . The vortex-vortex interactions become more important for higher values of . For , illustrated in Fig. 13(a,b), both and show a peak feature, while the anisotropy remains near . The depinning thresholds increase with increasing for low due to the fact that the increasing strength of the vortex-vortex interactions raises the interstitial pinning barriers. The threshold does not continue to monotonically increase with increasing since the vortex configuration at is a honeycomb lattice. As a result, the elastic energy cost of maintaining the honeycomb structure increases with increasing , and for , a portion of the vortices shift out of the pinning sites to form a distorted triangular lattice similar to that seen for low in Fig. 11(c). This causes a drop in the depinning thresholds with increasing . We expect similar behavior for other matching fields at which an ordered but non-triangular vortex lattice forms.

In Fig. 13(c) we plot and versus for the same system with , and in Fig. 13(d) we show the corresponding . At this matching field, the vortex lattice is triangular as seen in Fig. 5(d). Therefore, the vortex positions do not shift as increases, unlike the case for , and the depinning thresholds increase monotonically with increasing . The anisotropy does not vary strongly with . In Fig. 14(a,b) we show , , and versus for the same system at , where the vortex lattice is triangular as indicated in Fig. 6(a). The depinning thresholds increase monotonically with increasing while the anisotropy remains constant at . At , where Fig. 6(d) shows that the vortex lattice is disordered, the depinning thresholds are isotropic and increase monotonically with increasing , as illustrated in Fig. 14(c,d). For any matching field , if is increased above the range of values considered here, multiple vortex pinning by the pinning sites eventually occurs when the vortex lattice constant becomes of the order of the pinning radius and very little distortion of the vortex lattice is required to shift the vortices into the pinning sites. The occurrence of multiple vortex pinning would alter both the depinning thresholds and the anisotropy.

### iii.3 Anisotropy for

For , interstitial vortices are present and the depinning threshold is determined by the strength of the vortex-vortex interactions, provided that the pinning is strong enough to produce a plastic depinning transition. For , the depinning threshold is controlled by a combination of the vortex-vortex interaction strength and the strength of the pinning sites. For a triangular pinning array at the matching and submatching fields of , 1/3, and , the vortex lattice is triangular and every vortex is trapped by a pinning site. In this case, the depinning thresholds are determined by the maximum force exerted by the pinning sites and are independent of the direction of the external drive, so the depinning is isotropic. This is in agreement with previous work, where an isotropic depinning threshold was observed for vortices in triangular pinning arrays at (30); (31). In Fig. 15 we plot and versus at and for a system with . The depinning thresholds are isotropic for both fields, and at , as expected. At , the vortex lattice is disordered since a triangular vortex lattice cannot match the pinning array at this filling, as shown in Fig. 11(c) of Ref. (41). Although all of the vortices are pinned, due to the disorder of the vortex structure, some vortices experience stronger vortex-vortex repulsion from neighboring vortices than other vortices, and as a result the vortex-vortex interactions do not cancel at as they do at . This lowers the depinning thresholds and causes both and to decrease with increasing due to the dependence of the depinning force on the vortex-vortex interaction strength. In general, we observe little or no anisotropy in the depinning thresholds for . If multiple vortex pinning occurs at , so that no interstitial vortices are present, we expect that the same type of depinning phenomena seen for will appear and there will be little anisotropy. This suggests that the appearance of anisotropic depinning forces and the existence of voltage-current curves with different values of in different directions are indicators of the presence of interstitial vortices in the system.

### iii.4 Temperature Effects

We next consider how robust the anisotropy is to thermal fluctuations. We anneal the system down to a finite temperature and then measure the depinning forces in the and -directions. The temperature is given in terms of the melting temperature at which the vortices begin to diffuse significantly and the system is in a molten state. In Fig. 16(a) we plot and versus for a triangular pinning lattice system with at , and in Fig. 16(b) we show the corresponding anisotropy . Both depinning thresholds monotonically decrease with increasing temperature and reach zero at . The anisotropy decreases with increasing temperature for low temperatures; however, just below the melting temperature, passes through a peak. A similar trend is seen for , as shown in Fig. 16(c,d). These results indicate that the anisotropy should be robust against temperature. This agrees with the experiments in Ref. (31), which found that the anisotropy became more pronounced at higher temperatures. The enhancement of the anisotropy at higher temperatures could also result from the presence of weak random intrinsic pinning in the sample, which would reduce the magnitude of the anisotropy but which can be washed out by thermal fluctuations.

## Iv Anisotropy in Square Pinning Arrays

For square pinning arrays, we find that the depinning thresholds at most fields are isotropic. This is a result of the fact that the same type of vortex motion occurs in both directions for most fields, as illustrated for in Fig. 17. In Fig. 18(a) we show the isotropic velocity-force curves at for and -direction driving. Since the perpendicular directions of the square pinning array are identical, unlike the perpendicular directions of the triangular pinning array, it is not surprising that most matching fields have the same depinning thresholds in both driving directions for the square array. Nevertheless, strongly anisotropic depinning occurs at and , as shown in Fig. 18(b,c).

In Fig. 19(a) we plot the vortex positions and pinning site locations for , where anisotropic depinning occurs. Here a triangular vortex lattice that is aligned with the -axis forms. Under -direction driving, the interstitial vortices flow in one-dimensional channels between adjacent columns of pinning sites, as seen in Fig. 19(b) where of the interstitial vortices have depinned. A simple channeling motion cannot occur for direction driving, so is much higher than . Just above the depinning transition for -direction driving, as illustrated in Fig. 19(c), every other column of vortices shifts in the direction in order to permit every other vortex to join a one-dimensional flowing channel passing between adjacent rows of pinning sites. The remaining interstitial vortices remain trapped between the pinning sites. After this rearrangement, of the interstitial vortices flow in the steady-state one-dimensional channels shown in Fig. 19(d), which are similar to the channels that form for -direction driving. The ground state shown in Fig. 19(a) is degenerate, and a state with the vortex lattice aligned along the -direction has equal energy. Thus, the realignment process seen in Fig. 19(c) is simply a shift of the vortices into the other ground state prior to the onset of flow. Figure 18(b) shows that is approximately 1.75 times higher than . Due to the symmetry of the square lattice, either the or the -direction can show a higher depinning threshold depending on the initial configuration of the vortex lattice. This is distinct from the triangular pinning lattice, where the higher depinning force always occurs in the same direction at a given field. If the annealing process for the square pinning lattice is repeated with different initial conditions, the vortex lattice has a 50% chance of aligning with the -direction, in which case the anisotropy will be reversed from that shown in Figs. 18, 19. If the driving force is cycled, the velocity-force curve in the hard driving direction is hysteretic during the first cycle due to the realignment effect, while there is no hysteresis for the easy driving direction or for subsequent cycles in the initially hard driving direction.

In Fig. 20 we illustrate the effect of changing and on the depinning anisotropy for the square pinning lattice sample with . We show a case where the initial vortex lattice orientation is along the -axis, as in Fig. 19(a), which gives . Figure 20(b) indicates that the anisotropy saturates at for , while Fig. 20(a) shows that the depinning thresholds also saturate above this value of . For , the depinning is elastic and the anisotropy is lost. For fixed the depinning thresholds increase monotonically with increasing , as seen in Fig. 20(c), while in Fig. 20(d) the anisotropy passes through a shallow extremum of .

We find a similar anisotropic depinning behavior for . The vortex configurations for this field at are illustrated in Fig. 21(a). The vortex lattice is aligned in the direction, but just as at , there are two degenerate ground states, and a vortex lattice that is aligned in the -direction would have the same energy. There are two columns of interstitial vortices between adjacent columns of pinning sites. The overall vortex lattice structure is not triangular and there are dislocations present in the lattice. Along certain columns, neighboring interstitial vortices in neighboring columns lie along a line tilted by with respect to the -axis, while along other columns, neighboring interstitial vortices lie along a line tilted by with respect to the -axis. The dislocations in the vortex lattice are all aligned in the same direction, resulting in a smectic structure. This smectic state for the square pinning array has not been observed in previous work. There are two possible low-energy orientations for the vortex lattice, just as in the case: the -axis orientation shown in Fig. 21(a), or the same state rotated by and aligned with the -axis. For depinning in the -direction, Fig. 21(b) shows that the two columns of interstitial vortices depin into flowing one-dimensional channels, while two interstitial vortices remain trapped between adjacent pairs of pinning sites. For driving along the -direction, the same type of lattice reorientation found for occurs at , as illustrated in Fig. 21(c,d). The columns of interstitial vortices shift in such a way that two one-dimensional rows of interstitial vortices form, while two interstitial vortices remain trapped between adjacent pairs of pinning sites. The final vortex configuration for -direction driving, in Fig. 21(d), is a rotated version of the configuration for -direction driving seen in Fig. 21(b). The same hysteretic voltage-current response should occur for driving in the hard direction at as at . In Fig. 18(c) the velocity-force curves at show an anisotropy The more pronounced anisotropy at compared to is due to the fact that the overall structure at is smectic and the dislocations are aligned in the -direction, further decreasing the depinning force along this direction.

We expect that the depinning anisotropy in the square pinning arrays will be more difficult to observe experimentally than the anisotropy in the triangular pinning arrays due to the existence of twofold-degenerate ground states for the square system. If other forms of quenched disorder, such as intrinsic pinning, are present and the system is large, domains of the two different orientations could form which would render the depinning thresholds isotropic. On the other hand, we found that an applied drive can readily align the vortex lattice structures in the driving direction at and . Therefore, it may be possible to prepare the system in an aligned state using an external drive, and then measure the anisotropy of the depinning forces starting from this aligned state. This procedure should permit the anisotropy in the square pinning lattice system to be observed experimentally.

## V Anisotropy in Honeycomb and Kagomé Pinning Arrays

We next examine anisotropy in honeycomb and kagomé pinning arrays. Since the honeycomb and kagomé arrays can be constructed by removing selected pinning sites from a triangular array, it might be expected that the anisotropy would follow the same trend found for the triangular pinning arrays. In particular, one could expect that the depinning threshold would always be higher in the -direction. Instead, we find that the honeycomb and kagomé pinning arrangements show an anisotropy that undergoes reversals as a function of the applied magnetic field. In Fig. 22(a) we plot and as a function of for a honeycomb pinning system with , , and . Here, is the matching field for a honeycomb lattice constructed from a triangular lattice with matching field , and we have (26). Figure 22(b) shows versus for the honeycomb pinning array. For , and the depinning threshold is higher for -direction driving. This anisotropy is reversed from that seen in the triangular pinning arrays. The reversal can be understood by examining the vortex positions at in Fig. 23(a). Each interstitial vortex is located at the position where the pinning site was removed from the triangular lattice in order to create the honeycomb lattice. There is a pin-free channel of motion which the interstitial vortices can follow for driving in the -direction, as seen in Fig. 23(c). For -direction driving, the path of the interstitial vortex is blocked by pinned vortices, creating a much stronger barrier for depinning, and once the vortices begin to move, they follow the winding paths illustrated in Fig. 23(b).

In Fig. 24 we plot and versus for driving in the and -directions, respectively. The depinning in the honeycomb pinning array is anisotropic at both , shown in Fig. 24(a), and at , shown in Fig. 24(b). For the anisotropy is reversed compared to the lower fields. This is due to the formation of -mer states in the large interstitial spaces of the honeycomb lattice which permit a portion of the interstitial vortices to be aligned in the -direction between the pinning rows. In Fig. 23(d) we illustrate the vortex configuration at , where the interstitial vortices form triangular shapes within the large interstitial sites. The vortex at the top of each interstitial triangle forms a nearly one-dimensional channel in the -direction with the vortices at the base of the adjacent interstitial triangles, permitting easy depinning in the -direction into the flow pattern shown in Fig. 23(e). For depinning in the -direction, where the threshold is higher, the complex but ordered flow pattern in Fig. 23(f) appears. For , the anisotropy reverses again, as seen in Fig. 22.

Figure 22(c,d) shows that similar anisotropy reversals occur for depinning in kagomé pinning arrays, where is the kagomé matching field for a pinning lattice constructed from a triangular lattice with matching field and . For , Fig. 22(d) indicates that , but that this pattern reverses several times for increasing . Just as in the honeycomb pinning array, in the kagomé pinning array the reversals originate from the formation and alignment of interstitial vortex -mer states in the large interstitial sites. In general, the anisotropy is smaller for the honeycomb and kagomé pinning arrays than for the triangular or square pinning arrays.

We illustrate the state for in Fig. 25(a), where there is one vortex per large interstitial site. In Fig. 25(b) we show the vortex trajectories just above depinning for driving in the -direction at . Each interstitial vortex diverts to the right or left around a pinned vortex, forming large asymmetric patterns of vortex flow around trios of occupied pins. The flow is slightly disordered, since in some cases two moving interstitial vortices approach the same pinned vortex simultaneously, causing one of the vortices to move outside of the flow pattern temporarily. For -direction driving at , shown in Fig. 25(c), the interstitial vortices flow around individual pinning sites in the sparse pinning rows, creating an asymmetric sinusoidal channel pattern. Figure 24(c) indicates that at this field, and depinning is easiest along the -direction. A similar anisotropy occurs for , as shown in Fig. 24(d). In Fig. 25(d) we illustrate the state at , where the large interstitial sites capture five vortices which form a pentagon structure. For depinning in the -direction at this field, only a portion of the interstitial vortices move in winding channels between the pinning rows while the large interstitial sites all capture three vortices, as seen in Fig. 25(e). The interstitial vortices which are not part of the channeling flow still undergo a small circular motion as the flowing interstitial vortices move past. At higher drives, all of the interstitial vortices depin and a step appears in the velocity-force curve. For depinning in the -direction at , shown in Fig. 25(f), the trajectories of the moving interstitial vortices are much more tortuous. Although a larger number of interstitial vortices spend at least part of the time flowing for depinning in the -direction than for depinning in the -direction, Fig. 24(d) shows that . This is because the vortex motion in Fig. 25(f) is not continuous. Instead, the interstitial vortices hop by one lattice constant and then repin so that a pulse of motion passes though the system. For depinning in the -direction, shown in Fig. 25(e), the vortices are continuously moving through the system.

## Vi Summary

We have shown that when interstitial vortices are present for fields beyond the first matching field in regular artificial pinning arrays, the transport response is anisotropic. For triangular pinning arrays, we find that the depinning thresholds are always higher in the -direction, defined to be perpendicular to a symmetry axis of the pinning lattice. This is in general agreement with previous numerical and experimental studies up to the third matching field for triangular pinning arrays. We show that the anisotropy also occurs for higher matching fields and that at certain matching fields the depinning is isotropic. The degree of anisotropy can be controlled by changing the strength of the pinning. For weak pinning, when the depinning transition is elastic and all of the vortices depin simultaneously, the anisotropy is reduced or destroyed. For strong pinning, the depinning is isotropic for matching fields at which the vortex lattice does not order. We find that the velocity-force curves can have different slopes for driving in the different directions depending on the number of interstitial vortices that initially depin. In some cases, although the depinning threshold is lower in one direction, the slope of the velocity-force curve is also lower in that direction, producing a crossing in the velocity-force curves for the two directions of driving. For fields at or below the first matching field, the anisotropy is strongly reduced, implying that the anisotropy should disappear if multiple vortex pinning rather than interstitial vortex pinning occurs above the first matching field. The vortex dynamics can be distinctly different for the two driving directions, with ordered flow states occurring for driving in one direction and disordered flow states appearing for driving in the other direction. The anisotropy is robust against temperature fluctuations and can be enhanced near the vortex lattice melting transition. For square pinning arrays, the two perpendicular driving directions are symmetric; nevertheless, at certain matching fields where triangular or smectic vortex structures form, a strong depinning anisotropy can occur. The easy-flow direction for square lattices is not fixed but depends on which of the two degenerate ground states is formed by the initial vortex lattice. A sufficiently large applied drive realigns the vortex lattice and sets the easy-flow direction in the direction of the drive. If the initial vortex lattice formed with domains of the two degenerate ground states, the domains can be eliminated and a pure ground state formed simply by sweeping the driving force. The anisotropic depinning behavior of this pure state can then be probed with small drives. Honeycomb and kagomé pinning arrays have a smaller anisotropy than that seen for the triangular and square arrays. The anisotropy for the honeycomb and kagomé arrays shows a series of reversals as a function of field due to the formation of vortex molecular crystals, which have orientational degrees of freedom that can lock to different angles.

###### Acknowledgements.

This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396.### References

- A.T. Fiory, A.F. Hebard, and S. Somekh, Appl. Phys. Lett. 32, 73 (1978).
- V.V. Metlushko, M. Baert, R. Jonckheere, V.V. Moshchalkov, and Y. Bruynseraede, Solid State Commun. 91, 331 (1994).
- M. Baert, V.V. Metlushko, R. Jonckheere, V.V. Moshchalkov, and Y. Bruynseraede, Phys. Rev. Lett. 74, 3269 (1995); M. Baert, V.V. Metlushko, R. Jonckheere, V.V. Moshchalkov, and Y. Bruynseraede, Europhys. Lett. 29, 157 (1995).
- A. Bezryadin and B. Pannetier, J. Low Temp. Phys. 102, 73 (1996).
- E. Rosseel, M. Van Bael, M. Baert, R. Jonckheere, V.V. Moshchalkov, and Y. Bruynseraede, Phys. Rev. B 53, R2983 (1996); L. Van Look, E. Rosseel, M.J. Van Bael, K. Temst, V.V. Moshchalkov, and Y. Bruynseraede, Phys. Rev. B 60, R6998 (1999).
- K. Harada, O. Kamimura, H. Kasai, T. Matsuda, A. Tonomura, and V.V. Moshchalkov, Science 274, 1167 (1996).
- V. Metlushko, U. Welp, G.W. Crabtree, Z. Zhang, S.R.J. Brueck, B. Watkins, L.E. DeLong, B. Ilic, K. Chung, and P.J. Hesketh, Phys. Rev. B 59, 603 (1999).
- U. Welp, Z.L. Xiao, J.S. Jiang, V.K. Vlasko-Vlasov, S.D. Bader, G.W. Crabtree, J. Liang, H. Chik, and J.M. Xu, Phys. Rev. B 66, 212507 (2002).
- J.I. Martín, M. Vélez, J. Nogués, and I.K. Schuller, Phys. Rev. Lett. 79 , 1929 (1997); A. Hoffmann, L. Fumagalli, N. Jahedi, J.C. Sautner, J.E. Pearson, G. Mihajlovic, and V. Metlushko, Phys. Rev. B 77, 060506(R) (2008).
- D.J. Morgan and J.B. Ketterson, Phys. Rev. Lett. 80, 3614 (1998).
- A.A. Zhukov, P.A.J. de Groot, V.V. Metlushko, and B. Ilic, Appl. Phys. Lett. 83, 4217 (2003).
- A.V. Silhanek, S. Raedts, M.J. Van Bael, and V.V. Moshchalkov, Phys. Rev. B 70, 054515 (2004).
- G. Karapetrov, J. Fedor, M. Iavarone, D. Rosenmann, and W.K. Kwok, Phys. Rev. Lett. 95, 167002 (2005).
- Y. Fasano and M. Menghini, Supercond. Sci. Technol. 31, 023001 (2008).
- S.B. Field, S.S. James, J. Barentine, V.V. Metlushko, G. Crabtree, H. Shtrikman, B. Ilic, and S.R.J. Brueck, Phys. Rev. Lett. 88, 067003 (2002).
- A.N. Grigorenko, G.D. Howells, S.J. Bending, J. Bekaert, M.J. Van Bael, L. Van Look, V.V. Moshchalkov, Y. Bruynseraede, G. Borghs, I.I. Kaya, and R.A. Stradling, Phys. Rev. B 63, 052504 (2001); A.N. Grigorenko, S.J. Bending, M.J. Van Bael, M. Lange, V.V. Moshchalkov, H. Fangohr, and P.A.J. de Groot, Phys. Rev. Lett. 90, 237001 (2003).
- C. Reichhardt and N. Grønbech-Jensen, Phys. Rev. Lett. 85, 2372 (2000).
- C. Reichhardt, G.T. Zimányi, R.T. Scalettar, A. Hoffmann, and I.K. Schuller, Phys. Rev. B 64, 052503 (2001).
- G.R. Berdiyorov, M.V. Milosevic, and F.M. Peeters, Phys. Rev. Lett. 96, 207001 (2006); Phys. Rev. B 74, 174512 (2006).
- C. Reichhardt, C.J. Olson, and F. Nori, Phys. Rev. B 57, 7937 (1998).
- Q.H. Chen, G. Teniers, B.B. Jin and V.V. Moshchalkov, Phys. Rev. B 73, 014506 (2006).
- G.R. Berdiyorov, M.V. Milosevic, and F.M. Peeters, Europhys. Lett. 74, 493 (2006).
- C. Reichhardt, C.J. Olson, R.T. Scalettar, and G.T. Zimányi, Phys. Rev. B 64, 144509 (2001).
- C.J. Olson Reichhardt, A. Libál, and C. Reichhardt, Phys. Rev. B 73, 184519 (2006).
- C. Reichhardt, G.T. Zimányi, and N. Grønbech-Jensen, Phys. Rev. B 64, 014501 (2001).
- C. Reichhardt and C.J. Olson Reichahrdt, Phys. Rev. B 76, 064523 (2007).
- M.F. Laguna, C.A. Balseiro, D. Domínguez, and F. Nori, Phys. Rev. B 64, 104505 (2001).
- G.R. Berdiyorov, V.R. Misko, M.V. Milosevic, W. Escoffier, I.V. Grigorieva, and F.M. Peeters, Phys. Rev. B 77, 024526 (2008).
- L. Van Look, E. Rosseel, M.J. Van Bael, K. Temst, V.V. Moshchalkov, and Y. Bruynseraede, Phys. Rev. B 60, R6998 (1999).
- T.C. Wu, P.C. Kang, L. Horng, J.C. Wu, and T.J. Yang, J. Appl. Phys. 95, 6696 (2004).
- R. Cao, T.C. Wu, P.C. Kang, J.C. Wu, T.J. Yang, and L. Horng, Sol. St. Commun. 143, 171 (2007).
- L. Van Look, B.Y. Zhu, R. Jonckheere, B.R. Zhao, Z.X. Zhao, and V.V. Moshchalkov, Phys. Rev. B 66, 214511 (2002).
- M. Velez, D. Jaque, J.I. Martín, M.I. Montero, I.K. Schuller, and J.L. Vicent, Phys. Rev. B 65, 104511 (2002).
- J.E. Villegas, E.M. Gonzalez, M.I. Montero, I.K. Schuller and J.L. Vicent, Phys. Rev. B 72, 064507 (2005).
- K. Mangold, P. Leiderer, and C. Bechinger, Phys. Rev. Lett. 90, 158302 (2003).
- P. Tierno, T.H. Johansen, and T.M. Fischer, Phys. Rev. Lett. 99, 038303 (2007).
- G. Coupier, M. Saint Jean and C. Guthmann, Phys. Rev. B 75, 224103 (2007).
- C. Reichhardt, K. Moon, R. Scalettar, and G. Zimányi, Phys. Rev. Lett. 83, 2282 (1999).
- M.-C. Cha and H.A. Fertig, Phys. Rev. Lett. 80, 3851 (1998).
- C. Reichhardt and N. Grønbech-Jensen, Phys. Rev. Lett. 85, 2372 (2000).
- C. Reichhardt and N. Grønbech-Jensen, Phys. Rev. B 63, 054510 (2001).