Designing high-performance superconductors with nanoparticle inclusions: comparisons to strong pinning theory

One of the most promising routes for achieving unprecedentedly high critical currents in superconductors is to incorporate dispersed, non-superconducting nanoparticles to control the dissipative motion of vortices. However, these inclusions reduce the overall superconducting volume and can strain the interlaying superconducting matrix, which can detrimentally reduce $T_c$. Consequently, an optimal balance must be achieved between the nanoparticle density $n_p$ and size $d$. Determining this balance requires garnering a better understanding of vortex-nanoparticle interactions, described by strong pinning theory. Here, we map the dependence of the critical current on nanoparticle size and density in (Y$_{0.77}$,Gd$_{0.23}$)Ba$_2$Cu$_3$O$_{7-\delta}$ films in magnetic fields up to 35 T, and compare the trends to recent results from time-dependent Ginzburg-Landau simulations. We identify consistencies between the field-dependent critical current $J_c(B)$ and expectations from strong pinning theory. Specifically, we find that that $J_c \propto B^{-\alpha}$, where $\alpha$ decreases from $0.66$ to $0.2$ with increasing density of nanoparticles and increases roughly linearly with nanoparticle size $d/\xi$ (normalized to the coherence length). At high fields, the critical current decays faster ($\sim B^{-1}$), suggestive that each nanoparticle has captured a vortex. When nanoparticles capture more than one vortex, a small, high-field peak is expected in $J_c(B)$. Due to a spread in defect sizes, this novel peak effect remains unresolved here. Lastly, we reveal that the dependence of the vortex creep rate $S$ on nanoparticle size and density roughly mirrors that of $\alpha$, and compare our results to low-$T$ nonlinearities in $S(T)$ that are predicted by strong pinning theory.


I. INTRODUCTION
Recently, the demand for high-performance superconductors -having large current-carrying capacities and resiliency to strong magnetic fields -has increased. These materials are necessary constituents in a multitude of applications, 1,2 including next-generation fusion reactors for clean energy production, 2-6 cables for power transmission, particle accelerators for high-energy physics, 7,8 magnetic resonance imaging for radiology, [9][10][11][12] and nuclear magnetic resonance spectroscopy for biological applications. 11,13 Meeting established metrics for operating conditions requires careful design of the conductor material, necessitating substantial changes from the intrinsic crystalline structure.
Interest in nanoparticle inclusions spans beyond their ability to boost J c , but also encompasses accompanying increases in the irreversibility field 22 H irr and decreases in the rate of thermally activated vortex motion (creep). 23,24,28,29 Most notably, YBa 2 Cu 3 O 7−δ (YBCO) films containing nanoparticles exhibit a non-monotonic temperature-dependent creep rate S(T ), dipping to very low values at T ∼ 0.3T c -0.5T c , which has not otherwise been achieved in YBCO single crystals 30 nor iron-based films containing nanoparticles. Understanding the mechanism behind this dip is key to further minimizing S while maximizing J c in YBCO and other superconductors. Remarkably, the beneficial effects of nanoparticles cannot be presumed. Their effects depend specifically on the intrinsic microstructure and measurement conditions (e.g., temperature T and applied magnetic field H). For example, BaZrO 3 inclusions unfavorably increased the creep rate S in BaFe 2 (As 0.67 P 0.33 ) 2 thin films exposed to low magnetic fields. 31 Despite the vast studies demonstrating remarkable advances in growth techniques and accompanying enhancements of J c , [15][16][17][18][19][20][21][22][23][24][25][26] many open questions remain largely unanswered: Why do nanoparticles induce non-monotonic temperature-dependent creep rates in YBCO films and what controls the local minimum in S(T )? Under what conditions do nanoparticles act collectively, versus as strong pinning centers? What is the maximum non-superconducting volume fraction that can be consumed by nanoparticle inclu-sions to maximize J c without compromising the critical temperature T c ? Furthermore, what is the optimal combination of nanoparticle size d and density n p to constitute this nonsuperconducting volume? Answering these questions will enable optimized design of the microstructure in superconductors.
Here, we present a comparative study of vortex pinning in (Y 0.77 ,Gd 0.23 )Ba 2 Cu 3 O 7−δ films containing different sizes and densities of nanoparticle inclusions and compare the results to the expectations of strong pinning theory. [32][33][34][35] A previous study 23 mapped the dependence of nanoparticle size and density on J c at high temperatures T /T c ∼ 0.7-0.9. We build upon this work by presenting a similarly systematic study at significantly lower temperatures T /T c ∼ 0.05-0.5 as well as at high fields up to 35 T. Mapping the trends at low temperatures enables direct comparisons to strong pinning theory formalisms that do not consider thermal activation, whereas high field measurements explore the high vortex density regime. We present four main results: First, at low fields, we find that J c increases roughly linearly with nanoparticle density n p , then J c (n p ) flattens out at higher densities. Second, we observe the two-stage power-law behavior of the critical current J c ∝ B −α that is expected, and previously observed, 24,28,30,[36][37][38][39][40][41][42] for the strong-pinning scenario. We subsequently find that α decreases with increasing nanoparticle density and increases with nanoparticle diameter d/ξ , normalized to the coherence length ξ ≡ ξ ab . Third, we study J c (B) at high vortex densities to look for the theoretically predicted peak effect that may occur when individual nanoparticles capture more than one vortex. Last, we find that samples containing denser nanoparticles produce slower creep rates, whereas larger nanoparticle diameters d/ξ result in increased creep rates at low temperatures T T c .

II. THEORETICAL BACKGROUND
The dynamics of vortices depends on the interplay between current-induced forces, thermal energy, vortex elasticity, and pinning forces exerted by defects. Broadly, the complexity of the problem is dictated by three constituents: First, the superconductor's intrinsic properties defining the 'pure' vortex state, i.e. penetration depth λ , coherence length, and electronic mass anisotropy γ 2 = m c /m a . Second, artificial or naturally occurring inhomogeneities in the superconducting matrix provide a defect structure (pinning landscape) that impedes vortex motion for sufficiently small driving forces. The third contribution originates from the sample geometry that is known to define a panoply of geometric and surface effects. [43][44][45][46][47][48][49] Segments of a vortex line or vortex bundles can be pinned by the independent action of large defects (strong pinning) or by the collective action of many weak, small pins (weak pinning). [50][51][52][53] In the absence of thermal effects, the optimal pinning configuration results from balancing the energy cost to deform a vortex state away from its lattice arrangement by an energy gain when vortices accommodate to neighboring pinning sites. Along a vortex line, segments of character-istic length L c are independently pinned, while in the transverse direction, pinning of vortices separated by more than R c is uncorrelated. The correlation volume thus becomes V c ∝ L c R 2 c ; its parametric dependence has been calculated for the weak and strong pinning cases. [50][51][52][53] Though pinning disrupts long-range order of the vortex lattice, short-range order exists within V c , which can contain a single vortex or a bundle of vortices, depending on the vortex density, where the vortex spacing is a 0 = (4/3) 1/4 (Φ 0 /µ 0 H) 1/2 . At finite temperatures, these single vortices or bundles hop out of their pinning sites due to thermal activation, such that the vortex arrangement relaxes (creeps) from its out-ofequilibrium arrangement. This is often captured by measuring the relaxation rate of the persistent currents, 54 creep-induced curvature in IV characteristics 17,35 near I I c for V ∼ I n , or imaging studies (e.g. scanning tunneling microscopy). 55 Moreover, at sufficiently high temperatures, the vortex state undergoes a melting transition in which pinning becomes ineffective. [56][57][58][59][60] A. Strong pinning theory: predictions for the critical current The limit of dilute, strong defects has been studied analytically in recent years. Contrary to the qualitative scaling arguments of weak collective pinning theory, the strong pinning framework provides quantitative results for various observables (J c 32,50,52,61,62 , Campbell penetration [63][64][65][66] , zero-and finite-temperature current-voltage characteristics. 34,35,67,68 ) Specifically, the critical current has been shown 52 to follow J c = J dp (n p a 0 ξ 2 )( f p /ε 0 ) 2 , over a wide field range, with J dp being the depairing current, n p the defect density, f p the elementary force provided by a single defect, and ε 0 the vortex line energy. This result implies a power-law field dependence J c ∝ B −α with α = 1/2, in agreement with an early work by Larkin and Ovchinnikov. 50 A recent study 32 demonstrated that the effect of defects trapping vortices transverse to the direction of the force increases the power α by 1/8 to 5/8. Furthermore, numerical simulations of the strong-pinning regime have revealed that the power-law behavior J c ∝ B −α persists beyond the limit of diluted defects, albeit with a reduced power α < 5/8. Specifically, it was found that the power α depends primarily on the volume fraction V occupied by the defect, only weakly depends on the defect size, and follows an empiric law α = 0.026-0.083 ln(V ) (for V > 10 −3 ), see also Ref. [69].
As the magnetic field increases, so does the vortex density, such that vortex-vortex interactions become increasingly relevant and vortices compete for pinning sites. When the vortex density surpasses the number of available pinning sites, some vortices may be effectively immobilized by caging effects (pinned vortices repel interstitial vortices) [70][71][72] or pinning centers may host more than one vortex. 32 This change in pinning behavior can evoke a change in the magnetic and transport properties. For example, Ref. [33] predicts a novel type of peak-effect in J c (B) that appears at high fields as the system transitions from nanoparticles capturing one to capturing multiple vortices. Despite rapid improvements in controlling the size of incorporated nanoparticles, experimental evidence for this effect remains elusive. Finally, at large fields, where the inter-vortex distance falls below the typical defect size, the critical current is expected to decay faster J c ∝ B −(1+β ) , where β > 0 accounts for the competition between vortices for the same defect.

B. Strong pinning theory: predictions for vortex creep
Material disorder creates a landscape of potential energy wells in which vortices localize to reduce their core energies by a pinning energy U c , defined at low currents J → 0. This energy landscape is then tilted by a current, reducing the energy barrier that a pinned vortex must surmount to escape from the well to U(J). The approximate time for thermal activation over the barrier is where t 0 is a microscopic attempt time for vortex creep (typically 1 µs). 54 Generally, the creep energy barrier is thought to depend on current as U(J) ≈ U c (1 − J/J c ) n , where n is determined by the pinning landscape. The Anderson-Kim model 51,73 predicts n = 1; as it disregards vortex elasticity and vortex-vortex interactions, it is often a useful approximation at low temperatures (T T c ), low magnetic fields, and in the early stages of the relaxation process (J J c0 ), where J c0 is the creep-free critical current. For currents J J c0 , the weak collective creep paradigm, which considers vortex elasticity, models the energy barrier as U(J) = U c (J c0 /J) µ , where µ is the glassy exponent that depends on the size of the vortex bundle that hops during the creep process and its dimensionality. 51,74 Interpolating between the two regimes, it is common to model the energy barrier as U(J) = (U c /µ)[(J c0 /J) µ − 1], which consequently covers a broad range of currents. Combining Eq. 1 with this interpolation formula, we expect J c to decay over time as The vortex creep rate then parameterizes this decay as such that, in the collective pinning scenario, the creep rate depends on the temperature and U c as Recent results from strong pinning theory 35 , however, predict that n ≈ 3/2, resulting in for J < J c . This leads, instead, to a thermal creep rate of 35 .

III. SAMPLES AND EXPERIMENTAL PROCEDURE
In this work, we studied one reference film of (Y, Gd)Ba 2 Cu 3 O 7−x and seven films containing BaMO 3 inclusions (for M = Zr, Hf, or Sn). The reference film was 900 nm thick and the thicknesses of the films containing BaMO 3 nanoparticles varied from 290 nm to 910 nm. Table I summarizes the parameters of each sample and assigns a sample ID that is referenced in the data. All films were grown epitaxially on buffered Hastelloy substrates by means of the trifluoroacetate metal-organic deposition process from Y-, Gd-, and Ba-trifluoroacetates and Cu-naphthenate solutions with the cation ratio of 0.77:0.23:1.5:3. For the films containing BaMO 3 nanoparticles, Zr-, Hf-, or Sn-naphthenate was introduced into the precursor solution. The interposing buffer was a stack of NiCrFeO, Gd 2 Zr 2 O 7 , Y 2 O 3 , MgO (deposited using ion beam assisted deposition), LaMnO 3 , and CeO 2 layers. Note that all samples contain a sparse distribution (0.3 × 10 21 /m 3 ) of Y 2 Cu 2 O 5 (225) nanoparticle precipitates that naturally form during the growth process. Further details of the growth process are elaborated elsewhere. 21 To map the average size and density of nanoparticle inclusions, the films were characterized via transmission electron microscopy (TEM) and energy dispersive spectroscopy. 23 For example, Fig. 1(a,b) shows TEM images of a (Y,Gd)BCO film containing, on average, 8 nm-sized BaHfO 3 nanoparticles that was processed under the same conditions as films in this study. As shown in the histogram in Fig. 1(c), all films contain a finite distribution of nanoparticle sizes and Table I indicates the mean diameter. Note that this study benefits from films grown using a modified process, described in Ref. [21], that results in a relatively narrow distribution range.
Incorporation of BaMO 3 nanoparticles in these films is accompanied by an increased density of two types of planar defects: c-axis oriented twin boundaries and stacking faults along the ab-plane. TEM studies showed that the stacking faults are short (50-100 nm) and do not segment the twin boundaries, which maintain their integrity throughout the film thickness (contrary to the films studied in Ref. [75]). The twin boundary spacing ranges from~25 nm in the film containing  the highest density of nanoparticles to~45 nm in the reference film. For more details, see Ref. [22].
To determine the temperature and field dependencies of the critical currents and creep rates, magnetization studies were performed using Quantum Design SQUID magnetometers (fields up to 7 T) at Los Alamos National Lab and the Colorado School of Mines, and a Vibrating-sample Magnetometer (VSM) at the National High Magnetic Field Lab (fields up to 35 T). In all cases, the magnetic field was applied parallel to the c-axis (perpendicular to the film plane). For the highfield VSM measurements, to observe measurably large signals over broad field ranges, we measured a stack of two identical cuts of each sample, effectively doubling the thickness and magnetization.
Vortex creep rates were determined by capturing the timedependent magnetization M(t) using standard methods, 54 which is expected to decay quasi-logarithmically over time t as per Eq. 2. To measure creep, first, the critical state is established by sweeping the field ∆H > 4H , where H is the minimum field at which magnetic flux fully penetrates the sample. Second, the field is set to the value of interest for the measurement and M is then recorded every~15 s for an hour. We then subtract the background produced by the substrate and adjust the time to account for the difference between the initial application of the field and first measurement. This difference is found by maximizing the correlation coefficient of a linear fit to ln M − lnt. Finally, the creep rate S is extracted from the slope of the linear fit, per Eq. 3.
The critical current J c (T, H) was extracted from the magnetization data using the extended Bean critical-state model for rectangular samples. 76,77 Here, ∆M is the difference between the upper and lower branches of the magnetization loops M(H), h is the film thickness, w = 3-5 mm and l = 3-5 mm specify sample width and length.

IV. RESULTS
A. Impact of film thickness: Elastic versus rigid vortex dynamics Beyond the bulk properties, it is important to recognize that the critical current J c and creep rate S acquire an explicit dependence on the sample thickness h, once the latter drops below the pinning length L c , i.e. h < L c . Phenomenologically, this represents a transition between elastic dynamics in thick (3D, bulk) samples and faster rigid vortex dynamics in thinner (quasi-2D) samples. That is, when h > L c , J c and S do not depend on thickness, whereas when h < L c , both S and J c become anticorrelated with thickness. Moreover, because the pinning length is temperature-dependent, so is the regime. When the superconductor contains strong pinning centers in addition to point defects, this pinning length deviates from the predictions of weak collective pinning theory, e.g. L c (T ) ≈ (ξ /γ)(J d /J c ) 1/2 for single vortex dynamics, to a mixed pinning length L mp c > L c that can be extracted from creep measurements. 78 From Fig. 2, we expect bulk behavior in films of all tested thicknesses at T = 5 K.
Our study contains three samples in the bulk-regime at all temperatures and five that host quasi-2D vortex dynamics at higher temperatures T 20 K. The thickness metric for categorizing the films was based off of a previous study 78 in which we extracted an effective pinning length L c from creep measurements collected on (Y,Gd)BCO films of a wide range of thicknesses as well as the data shown in Fig. 2 for (Y,Gd)BCO films containing BaHfO 3 nanoparticles. For clarity, throughout this work, we will analyze results from each category separately.
B. Low-field critical current dependence on nanoparticle diameter and density Figure 3(a) shows the dependence of critical current J c on applied magnetic fields up to 7 T, normalized by the upper critical field H c2 , for all 8 samples. In Fig. 3(b), a snapshot of this data (at two fields and temperatures) is re-plotted versus the nanoparticle density, showing that increasing the density is indeed effective in boosting J c , most dramatically at low temperature T = 5 K. Focusing on the field dependence of J c , we see that below the self-field µ 0 H s f = µ 0 γJ s f h/π, J c is field-independent, whereas above it, J c exhibits the expected power-law behavior J c ∼ B −α . Here, the exponent α can provide information regarding the pinning regime. Consequently, the dependence of α on defect sizes and densities can reveal whether the microstructure induces weak collective or strong pinning and provide information on vortex-defect interactions, as discussed in Sec. II. Hence, to understand the low-field H H c2 pinning regimes, we extract α from the slopes of linear fits to the data in Fig. 3(a) between the calculated H s f and H/H c2 = 0.01. The extracted α parameters are plotted versus the nanoparticle density n p and normalized diameter d/ξ in Figs. 3(c) and 3(d).
Samples 6, 7 and 9 are all approximately 900 nm thick, therefore host three-dimensional vortex dynamics. 78 For these samples, α ranges from 0.6 to 0.7 at 5 K, which is within the expectations of strong pinning theory. 32 Additionally, within the measured magnetic field range, we find that J c is more sensitive to nanoparticle density than pinning volume (see Table  I) and many small defects are preferable to fewer large ones.
In the thinner samples (#1-#5), the power-law exponent is more sensitive to the defect density and size than in the bulk samples. Figure 3(c) reveals a striking, systematic correlation between α and n p , increasing from 0.2 at a high density of 100 × 10 21 /m 3 to 0.53 for a lower density of 5 × 10 21 /m 3 , then only slightly increasing to~0.6 for the reference sample containing only 0.03 × 10 21 /m 3 225 precipitates. This is remarkably consistent with the expectations from the TDGL simulations for low versus high nanoparticle densities. Empirically, we also observe a strong dependence of α on the defect size at nearly fixed volume fraction (e.g. comparing all films with a volume fraction of approximately 8.7-8.6%). This result is in contrast with the weak dependence expected for three-dimensional strong pinning. 32

C. Critical current at high magnetic fields
To study strong pinning at high fields and look for the predicted peak effect, we measured samples 6 and 7 at the NHMFL in static fields up to 35 T (H/H c2 ≈ 0.3) using a VSM. Accordingly, magnetization loops were collected for each sample at temperatures up to 40 K, from which we calculated J c (H), shown in Figs. 4(a) and (b) for samples 7 and 6, respectively. Figure 4(c) shows how the enhancement in J c that is achieved by incorporating a high density of BaHfO 3 nanoparticles (sample 6) is maintained at high fields.
We extracted the high-field α (denoted α HF ) from linear fits to log J c versus log H/H c2 , restricted to H/H c2 = 0.05-0.15, for direct comparison with the strong pinning result in Ref. [32]. For data collected at T = 30-40 K, the upper bound of the fit was decreased to H/H c2 = 0.125 due to noise in the system as the signal approached its sensitivity limit. Figure 4(c) compares results for the two samples and the inset displays the dependence of α on temperature. From these comparisons, α appears independent of the nanoparticle diameter or density at these high magnetic fields. Furthermore, α ∼ 0.8 at 5 K and α ∼ 1.2-1.4 at 40 K. This is consistent with the expectations of increased α 1 at high fields as each pinning site is occupied by a vortex.
According to strong pinning theory, J c (H) should exhibit non-monotonic behavior at high fields as nanoparticles capture multiple vortices. For sample 6, we would expect the lowest field peak to occur around B∼Φ 0 /d 2 ≈ 10 T and additional peaks at higher fields as the nanoparticle pinning site occupancy increases. However, we observe no non-monotonic behavior over the entire field range. A few reasons may account for this discrepancy. First, the distribution of nanoparticles in our samples may not be as homogeneous as represented in the simulations and the size variations noted in Fig. 1(c) may smear the peak. In fact, variations in the diameter enter quadratically into the expected peak position. Second, all samples contain a sparse distribution of large 225 nanoparticle precipitates that likely capture multiple vortices at lower fields than the smaller BaMO 3 inclusions. However, the den-sity of 225 precipitates may be too small for multiple vortex occupancy of these defects to induce a measurable change in the J c (H) dependencies. Furthermore, the presence of twin boundaries may effectively increase the field at which the multiple vortices occupy the BaMO 3 inclusions to beyond our measurement range or further smear the expected peak. There is not yet a formalism that can consider mixed pinning landscapes -microstructures that consists of multiple different types of defects. This is partially because it is notably complicated to concurrently consider defects of varying dimensionality: point-like (e.g. vacancies), 2D planar (e.g. twin boundaries, stacking faults), and 3D (e.g. nanoparticles). Third, note that at higher temperatures, J c curves decreased to below the sensitivity of the VSM at fields around 25 T such that we were unable to resolve J c beyond H/H c2 ≈ 0.2.

D. Vortex creep dependence on nanoparticle diameter and density
Bulk samples. Figure 5(a,b) shows temperaturedependent creep rates for the three bulk films. At 1 T, we observe the typical non-monotonic dependence reported in multiple studies of YBCO films containing nanoparticle precipitates. 24,28,79,80 In agreement with these earlier studies, see e.g. Ref. [30], we find that the non-monotonicity tends to disappear upon increasing the magnetic field. The position of the peak in S is the same in all samples, therefore, may depend on a common property, such as the density of 225 precipitates. Additionally, from Figs. 5(a) and 5(b), we see that the BaMO 3 nanoparticles have far less effect on the creep rate at low temperatures than at high temperatures, consistent with Ref. [28]'s finding that defects generally affect creep more significantly at high-T .
At all fields, the film containing a high density of BaHfO 3 nanoparticles (sample 6) produces the lowest creep rates. The film containing BaSnO 3 (sample 9) and the reference film (sample 7) produce similar creep rates at 1 T, however, at 3 T, creep in sample 9 is notably slower than in sample 7. This suggests that the vortex density now is high enough that there is an insufficient number of 225 precipitates to pin all vortices and sample 9 benefits from the additional BaSnO 3 nanoparti- cles.
All samples. Figure 5(c) presents a comparison of creep in all samples. As expected based on results in Ref. [78], the thinner samples (1-5) exhibit little-to-no peak and faster creep than the bulk samples. Despite this difference, S(T ) at low temperatures is strikingly similar for all samples. Notably, creep increases approximately linearly with T , qualitatively adhering to the expectations of the Anderson-Kim model describing creep of rigid vortices S ∼ k B T /U c . 73 Moreover, all curves deviate from this linear behavior around the same temperature T = 20 K. Consequently, we can again see that that the effects of the inclusions are more significant at high temperatures than at low temperatures.
To view the dependence of the creep rate on nanoparticle properties and determine whether nanoparticles have an effect at low T , we plotted creep at 5 K versus the nanoparticle size and density in Fig. 5(d,e). A direct comparison of Figs. 5(d,e) and 3(c,d) reveals that, qualitatively, the dependence of S on the defect density and size is very similar to that observed for the exponent α. Specifically, creep shows a somewhat linear dependence on d/ξ for all samples containing added nanoparticle inclusions at both 0.3 T and 1.0 T. Data for the reference sample deviates from this trend, likely due to the large diameter of the dominate pinning sites-the 225 precipitates. We also observe a decrease in creep values as nanoparticle density decreases. Similar to the results for α(n p , d/ξ ) presented in Fig. 3(c), this decrease in the vortex creep rate suggests that a higher density of smaller defects may be preferable to a lower density of larger ones.

Low-temperature creep rates
It is common to compare creep data piecemeal to Eq. 4, such that at low temperatures T U c /[k B ln(t/t 0 )], the creep rate is expected to increase approximately linearly with temperature as S ≈ k B T /U c . Consequently, signatures of quantum creep are typically identified as a linear fit of S that extrapolates out to finite values at zero temperature S(T = 0) > 0 or a creep rate that deviates from linear behavior at low temperatures. However, Eq. 6 shows that in the strong pinning scenario, thermal creep should cause a nonlinear, convex increase in S at low temperatures, which becomes concave at the inflection point T = U c /[5 3/2 ln(t/t 0 )]. 35 Furthermore, a linear extrapolation at this inflection point indeed results in finite S 0 = S(T = 0) solely from thermal considerations. 35 Hence, the extrapolation alone is insufficient evidence of quantum creep, without measurements in the quantum-regime.
In Fig. 6, we plot a fit of our creep data for sample 6 to Eq. 6, where U c is the only free parameter and t 0 ∼ 10 −4 /ν (for H c2 /H ∼ 100 and ν ∼ 10 2 is related to the vortex velocity). From this fit, we extract U c /k B ≈ 655 K and T ≈ 2.7 K using t = 3600 s and demonstrate in the figure that a linear extrapolation around the inflection point leads to finite creep at zero temperature. Note that our data does not only deviate from linearity at low T , but also from the functional form predicted for thermal creep from strong pinning centers, likely due to a significant, unknown contribution from quantum creep within our measurement range.  6. (Yellow triangles) Low-temperature creep rate in a field of 1 T for sample 6. Blue solid curve is a fit of the data to the expectations for thermal creep from strong pinning theory (Eq. 6) using t 0 = 10 −6 . Deviations between experimental data and model possibly due to contributions from quantum creep. Dashed blue curve is a linear fit k B T /U c around the inflection point T . Notice that the dashed curve extrapolates to a finite S > 0 as T → 0, such that the conventional fit to a thermal region alone is insufficient evidence of quantum creep within the strong pinning paradigm.

Effect of Twin Boundaries
Though we have focused on the effect of nanoparticles, it is important to mention that twin boundaries also likely contribute to slowing creep. Understanding the role of twin boundaries, planar defects prevalent in many crystalline materials, is especially challenging --twin boundaries can have diametric roles in different regimes, either channeling or pinning vortices, 75,81,82 and the density of twin boundaries in films is difficult to control. Twin boundaries act as strong pinning centers (pinning a higher density of vortices than the bulk), perturb the vortex lattice, and form a barrier for vortex crossing such that vortices will accumulate on one side of the twin boundary. 83,84 A better understanding of the effects of twin boundaries on vortex motion is required to further improve the performance of superconductors.
As vortices are not rigid objects, vortex motion will not simply be either arrested or funneled by a twin boundary. The elasticity of a vortex means that thermal fluctuations can cause segments of the vortex to depin and form curved, half-loops that either slowly grow or rapidly expand if partially pinned by a neighboring defect. 51,85,86 These mechanisms for vortex motion will each have different effects on S.
To better understand this, let us consider the temperature dependence of the twin boundary pinning potential. Vortex attraction to twin boundaries is caused by suppression of the superconducting order parameter within the twin boundary, Ψ T B . This suppression extends a distance √ 2ξ (T ) = √ 2ξ 0 (1 − T /T c ) −1/2 around the twin boundary, where ξ 0 ≡ ξ (0), such that the breadth increases greatly as T → T c . The magnitude of this suppression is δ Ψ T B = 1 − |Ψ T B /Ψ bulk | 2 , where Ψ bulk is the order parameter in the bulk infinitely far from the twin boundary. While the exact temperature dependence of δ Ψ T B is unknown, suppression of the order parameter is weak (δ Ψ T B 1) for T T c 1 − [a/ξ (0)] 2 ) ≡ AT c and strong (δ Ψ T B ∼ 1 − (ξ 0 /a) 2 [1 − (T /T c ) 2 ] 1) in the opposite case, 51 where a = 3.8 Å is the lattice parameter. In YBCO samples, 51 A ≈ 0.97 so the suppression is mostly nearly temperature independent, but should become stronger and temperature dependent as we approach T c . This increased effectiveness of twin boundaries in pinning at high temperatures is consistent with the slowest high-T creep rates in the sample with the highest density of twin boundaries, and suggest that twin boundaries may contribute to lowering S at high temperatures.

V. CONCLUSION
In conclusion, we conducted a thorough study of the effects of nanoparticle size and density on the critical current and vortex creep rates in (Y,Gd)BCO films. By conducting measurements in a wide range of temperatures 5-70 K and fields up to 35 T, we compared trends to the expectations of recent results from strong pinning theory. Though the J c ∝ B −α dependence has been noted in numerous previous studies, we systematically unveiled a clear and similar dependence of α and S on both the nanoparticle density n p and normalized size d/ξ , which is our main result. Our results are consistent recent advances in strong pinning theory, which predicts a decrease in alpha from 0.66 to 0.2 with increasing n p and we observed a decrease from 0.65 to 0.2. Lowering α signifies higher retention of J c upon exposure to magnetic fields, therefore, is favorable for applications. At higher fields, 10-35 T, we observed larger α parameters (0.8-1.4), consistent with pinning sites becoming fully occupied. Lastly, the study showed a clear advantage to incorporating a high density of small nanoparticles versus a lower density of large ones.

AUTHORS' CONTRIBUTIONS
S.E. conceived and designed the experiment. S.E. and S.J. collected and analyzed the magnetization data. M.M. grew the samples, analyzed the microstructure, and assisted with data interpretation. R.Y. and T.K. performed transmission electron microscopy and other microstructural analysis. L.C. contributed to data interpretation. R.W. provided theory support and assisted with data interpretation. S.E., R.W., and S.J. wrote the manuscript. All author's commented on the manuscript.

ACKNOWLEDGMENTS
We would like to thank Dr. Eun Sang Choi at NHMFL for assistance with high magnetic field measurements and Dr. Barry Zink and Mike Roos for use of their Quantum Design magnetometer to conduct preliminary measurements on stacked films in preparation for the NHMFL measurements. This material is based upon work supported by the National Science Foundation under Grant No. 1905909 (S.E., S.J.) and the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences (L.C.). A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by National Science Foundation Cooperative Agreement No. DMR-1644779* and the State of Florida. R.W. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) TRR 288-422213477 Elasto-Q-Mat (project B01), and by the Heidelberger Akademie der Wissenschaften (WIN08). The work at Seikei University was supported by JST FOREST (Grant Number JPMJFR202G, Japan). A part of this work is supported JSPS KAKENHI (18KK0414 and 20H02184), PMAC for Science Research Promotion Fund and NEDO.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.