Vortex Dynamics and Dissipation Under High-amplitude Microwave Drive

In this paper, we describe the vortex dynamics under high-amplitude microwave drive and its effect on the surface resistance of superconductors. The vortex surface resistance is calculated with a Montecarlo approach, where the vortex motion equation is solved for a collection of vortex flux lines each oscillating within a random pinning landscape. This approach is capable of providing a detailed description of the microscopic vortex dynamics and in turn important insights into the microwave field amplitude dependence of the vortex surface resistance. The numerical simulations are compared against experimental data of vortex surface resistance at high microwave amplitude measured by means of bulk niobium superconducting-radio frequency cavities operating at 1.3 GHz. The good qualitative agreement of simulations and experiments suggests that the non-linear dependence of the trapped flux surface resistance with the microwave field amplitude is generated by progressive microwave depinning and vortex jumps.

In this paper, we describe the vortex dynamics under high-amplitude microwave drive and its effect on the surface resistance of superconductors. The vortex surface resistance is calculated with a Montecarlo approach, where the vortex motion equation is solved for a collection of vortex flux lines each oscillating within a random pinning landscape. This approach is capable of providing a detailed description of the microscopic vortex dynamics and in turn important insights into the microwave field amplitude dependence of the vortex surface resistance. The numerical simulations are compared against experimental data of vortex surface resistance at high microwave amplitude measured by means of bulk niobium superconducting-radio frequency cavities operating at 1.3 GHz. The good qualitative agreement of simulations and experiments suggests that the non-linear dependence of the trapped flux surface resistance with the microwave field amplitude is generated by progressive microwave depinning and vortex jumps.

I. INTRODUCTION
Upon cooldown below critical temperature in presence of an external magnetic field, magnetic flux quanta-socalled vortexes-can exist in thermodynamic equilibrium in the mixed state of type-II superconductors. 1,2 Below the lower critical field vortexes are not stable in the superconductor, however, because of the occurrence of defects in real materials, vortexes get pinned and survive even in the Meissner state. Incomplete Meissner effect where 100% of the applied field is trapped by the superconductor can be achieved with superconducting radio-frequency (SRF) cavities by field-cooling slowly across transition, as experimentally observed by many studies. [3][4][5][6] Driven by radio-frequency (RF) and microwave fields, vortexes oscillate dissipating power and contributing to the total surface resistance. Other than vortex oscillation, others contributions to the surface resistance are: the temperature dependent surface resistance due to thermally-excited quasi-particles above the gap, 7 surface resistance due to sub-gap states, 8 proximitycoupled normal-conducting inclusions (or oxides), 9 nanohydrides, 10 and dielectric dissipation due to spurious two-level systems at low microwave amplitudes and low temperatures. 11,12 Vortex dissipation under microwave drive is a topic that has been under investigation for many decades. Initial studies on type-II superconducting films 13 brought about the development of models to describe the vortex surface resistance dependence on frequency. 14,15 Since then, several other studies investigated the microwave behavior of conventional superconductors, [16][17][18] cuprates, [19][20][21] and iron-based superconductors [22][23][24] in presence of vortexes. In parallel, many theoretical works allowed to describe the expected penetration depth change due to vortex oscillation, 25 the vortex surface impedance temperature-dependence in high-T c superconductors, 26 and the precise electrodynamics of vortexes lattices in the mixed state. 27 In the framework of modern particle accelerators em-ploying SRF cavities as accelerating devices, dissipation due to vortex motion under microwave drive is a topic of central importance, inasmuch trapped-fields in the order of ∼1% of the earth magnetic field can increase the surface resistance of about 1 nΩ, 6 a large value for stateof-the-art cavities capable of achieving surface resistance values as low as 4 nΩ at 2 K. 28,29 Early studies on S-band SRF bulk Nb cavities 30 recognized the effect of vortexes on the surface resistance, identifying the microwave amplitude dependence of the latter to follow a non-linear behavior. Later studies showed that up to moderate gradients the microwave field amplitude dependence of the vortex surface resistance in Nb-sputtered on copper SRF cavities is linear. 31 Studies of vortex dynamics in cavity-grade Nb and theoretical description of the phenomenon observed in the zero microwave field amplitude approximation were also studied in detail. 32 More recently, intensive studies of the trapped flux surface resistance at low and moderate accelerating gradients as a function of the electron mean-free-path (l) discovered a bell-like trend as a function of l, 6,33 which is well described by the interplay of pinning-and flux-flowlimited dissipation. [34][35][36] The studies also reported an almost linear dependence of the vortex dissipation with the microwave field amplitude in the resonator up to moderate values.
The microwave field amplitude dependence can be understood as the occurrence of non-quadratic terms in the pinning potential. 37 The linear dependence was shown to be well described by adding a cubic term to a parabolic pinning potential in the vortex motion equation. 37 However, this approach is limited since built around the single-vortex dynamics, while in the material there are several vortexes contributing to the overall dissipation, each one interacting with its surrounding.
A second approach, 38 that overcomes the latter issue, implements a mean-field method to solve the vortex motion equation for collective weak pinning. In this case the microwave field amplitude dependence has two regimes: i) hysteretic losses at low amplitudes where the depen- dence is linear, and ii) viscous losses at higher amplitudes where the surface resistance saturates to a constant value.
In the present work, the high-amplitude microwave field dependence of the vortex surface resistance is explored by comparing the experimental data acquired for a high-performing single cell Nb accelerating cavity to the surface resistance calculated by means of numerical simulations. The model developed is based on a Montecarlo approach that assumes a random distribution of pinning potentials in the vortex oscillation plane. The simulation scheme here presented allows to gather insights into the dynamics of vortexes oscillation in random pinning potentials such us vortex jumps an microwave depinning, which are not captured by previous models, but key to correctly interpret the experimental data.

II. EXPERIMENTAL DATA
Experimental data of vortex surface resistance was estimated from microwave measurements of a singlecell TESLA-type 39 bulk niobium cavity operating at 1.3 GHz.
The cavity surface was prepared with a series of treatments designed to increase the maximum accelerating gradient achievable by the resonator, so that to capture the power dissipation fingerprint due to the highamplitude microwave behavior of vortexes. The cavity was initially electropolished to remove 200 µm and subsequently baked at 800 • C for 3 hours to degas hydrogen and release internal stress. After this initial surface preparation the cavity was high-pressure water rinsed in a class 10 clean-room to remove field-emitting particles form the inner surface and assembled with antennae needed for microwave characterization. The cavity was then evacuated and baked at 120 • C for 48 hours in situ.
The resonator was equipped with three single-axis flux gates magnetometers from Bartington Instruments located symmetrically around the cavity equator, a set of Helmholtz coils aligned coaxially to the cavity beam pipe, and a temperature mapping system (Tmap) 40 to precisely monitor the cavity temperature rise during the measurements.
The cavity was initially zero-field-cooled (ZFC) to 1.5 K and data Q-factor versus accelerating gradient acquired. In such conditions, the resonator reached maximum peak magnetic field (B p ) of 204.5 mT (48 MV/m) with Q-factor of 9.6 · 10 9 . The same cavity was then field-cooled (FC) very slowly to 1.5 K four times trapping respectively 7.5 mG, 30 mG, 50 mG, and 100 mG upon the superconducting transition.
As described in detail in Ref. 6, the vortex surface resistance is calculated by subtracting ZFC data from FC data acquired at the same temperature. By performing the measurement at 1.5 K, the thermal contribution to the surface resistance is negligible, and the systematic error introduced by the accuracy of the temperature control in the dewar minimized.
In Fig. 1(a), the vortex surface resistance (R f l ) as a function of the peak magnetic field is reported. The normalized vortex surface resistance over the magnetic field trapped during the cooldown (B 0 )-so-called trappedflux sensitivity-is reported in Fig. 1(b).
As shown in Fig. 1(b), once the vortex surface resistance is normalized by B 0 , the field dependence appears the same within the measurement error (∼10% of R f l value) 6,41 . This means that, within the range of B 0 and B p investigated, the peak magnetic field amplitude dependence of R f l is dictated by vortex dynamics and not by thermal feedback.
In agreement with previous measurement of vortex surface resistance as a function of accelerating gradient performed in thin film and bulk Nb cavities, 6,31 R f l increases almost linearly up to moderate fields (∼80 mT). For higher field amplitudes the vortex surface resistance deviates from the almost linear growth and the slope increases with B p , leading to higher values of resistance. Experimental data was also acquired by means of the Tmap system, composed of 36 boards arranged symmetrically around the cavity, each incorporating 15 thermometers in contact with the outer surface of the cavity as sketched in Fig. 2(a).
The data was acquired at B p = 166 mT both in the ZFC experiment and in the FC experiment with 100 mG applied during the cooldown. The ZFC tmap data was then subtracted from the FC data in order to obtain the heating pattern due to trapped-vortexes only. Temperature mapping data is shown in Fig. 2(c). The board number is reported in the abscissa, while the thermometer number along the ordinate.
Interestingly, the vortex dissipation shows a specific pattern with heating concentrated symmetrically with respect to the equator of the cavity (the equator location corresponds to thermometer number 8), where, instead, the temperature differential ∆T out is much lower than neighboring areas and tends to zero.
In this configuration, the trapped magnetic field B 0 is parallel to the cavity beam axis (as shown in Fig. 2(a)) and at the equator the microwave currents are parallel to the trapped field. Because of the angular dependence of the Lorentz force that drives vortex oscillation, this geometrical effect leads to lower dissipation at the equator compared to areas of the cavity where the flux has an angle almost perpendicular to the microwave currents.
Important to point out that this behavior is intertwined with the modulation of the microwave currents that varies with the local amplitude of the surface magnetic field B at the cavity surface. B is almost constant from thermometer 4 to thermometer 12-surface peak field B p at thermometer 6 and 10, and it decays rapidly to zero from thermometer 3 to 1 and from 13 to 15 (as shown in Fig. 2(b)).
This modulation of the surface field B explains why the vortex dissipation is almost zero at the cavity iris (thermometers 1 and 15) where flux and currents are instead close to orthogonality.

III. VORTEX DYNAMICS
Numerical simulations of vortex dynamics were performed in order to shed light on the field dependence of the vortex surface resistance.
Let's assume a superconducting surface at z = 0 parallel to the plane xy, where the semi-space z > 0 is filled by superconducting niobium and the semi-space z < 0 is vacuum free space. The magnetic field trapped (B 0 ) is parallel to the z axis, and we assume that upon transition the magnetic flux is subdivided in individual vortexes not interacting with each other B 0 << B c2 .
The microwave excitation-with frequency ω/2π-at the superconducting surface has magnetic field oscillating along the x axis, inducing current parallel to the y axis, which decays exponentially in the material along z, with characteristic decay length equal to the penetration depth λ. 42 The amplitude of the microwave field coincides with the peak magnetic field B p .
Neglecting the inertial term, the vortex motion equation has form: where u (t, z) is the vortex displacement as a function of time and depth z,u the time partial derivative, u the second partial derivative along z, and γ = φ 0 /(µ 0 λ). is the single-vortex line tension, 43 defined as: with κ = λ/ξ. The vortex viscosity η is equal to σφ 0 B c2 , 44 and σ is the niobium normal-state conductivity at low temperature. The pinning force per unit of length is introduced through f p (t, z).

A. Pinning Landscape
Vortex pinning in superconductors is explained as the tendency of the system to minimize the total loss in condensation energy. Being a vortex itself a singularity in the condensation energy, it interacts with areas in the material where superconductivity is suppressed or absent in order to minimize the overall condensation energy loss.
This mechanism intrinsically implies the minimum dimension for a defect to be an efficient pinning center as the coherence length ξ, which represents the minimum length over which a given change in condensed electron density can be made. 45 The coherence length in niobium is of the order of 38 nm or lower when scattering centers are present. This implies that any defect in the material with dimension of the order of ξ or larger could pin vortexes. Due to their small transverse dimension (∼ 1 − 5 nm), grain boundaries and single dislocations cannot be consider efficient pinning centers, whereas dislocations tangles, and material non-uniformity can extend for even hundreds of nanometers, allowing for more efficient pinning. 46,47 Due to the defective nature of real materials, many defects may pin simultaneously a single vortex along its length, making virtually impossible to measure directly the pinning potential, especially for bulk materials. Theoretical calculations 48 suggests that isolated pinning centers can be approximated by Lorentzian functions of the type U (r) = −U 0 /(1 + (r/ξ) 2 ), with r the radial distance form the vortex on the plane orthogonal to φ 0n , while experimental measurements on thin films 49 showed that real pinning potentials can be approximated by a summation of several of these Lorentzian functions.
In the case under study, the pinning landscape is defined in the oscillation plane of vortexes-xz, and it is constructed through a random bottom up approach based on the summation of N Lorentzian pinning potentials described above: Here, U i , X i , and Z i represent respectively the i-th pinning potential energy and (x,z) coordinates in the vortex oscillation plane, the variable a i is a dimensionless parameter that accounts for the anisotropy of the pinning potential along the z direction, while b is a multiplicative factor to account for the pinning center dimension. The pinning force is defined by the summation of the negative derivative along x of the i-th pinning potential, over the total number of pinning centers N .
The coordinates (X i , Z i ) and the anisotropy values a i are generated randomly within a certain domain, usually −1 < X i < 1 µm, 0 < Z i < 3 µm, and 0.001 < a i < 1. The pinning energy values U i are generated randomly following a normal distribution centered around an average value (U ave ). U ave is defined as the pinning energy that translates to a maximum pinning force per unit of length (f p,ave ) comparable to values found in literature. Usually for niobium the pinning force per unit of length varies depending on the material impurity content and crystallographic structure in between 10 −6 − 10 −4 N/m. 50,51 The input value f p,ave , together with the pinning potential depth distribution normal deviation (σ U = 5%U ave ), the coordinates (X i , Z i ), and parameters a i and b, describe completely the random pinning landscape used in the calculation.

B. Numerical Simulations
Numerical simulations of the vortex motion are defined within the xz plane domain where the pinning force landscape is defined by means of Eq. 3, assuming distributions of pinning energy, location, and anisotropy as described above.
The initial condition and the boundary conditions used in the calculation are the following: where z max represent the maximum depth of the solution domain defined from 0-the microwave surface, to z max .
In several works, 32,38 the boundary condition at the far edge of the solution domain is set to be u(t, z max ) = 0, implying a vortex strongly pinned in the bulk. On the opposite, in this work vortexes are allowed to find their own equilibrium position for any randomly-generated pinning landscape and we set u (t, z max ) = 0. As we will see later, this boundary conditions choice is mandatory to allow the vortex line displacement in the bulk to evolve with time until an equilibrium position is found, condition reached after several oscillations when the simulation converges to a steady-state solution.
The vortex motion equation is rewritten to be solved numerically, by defining lengths in units of λ and times in unit of 1/f : u =ũλ, z =zλ, X i =X i λ, Z i =X i λ, and t =t/f . The motion equation then becomes: The motion equation in Eq. 5 was solved numerically by means of the method of lines. The spatial derivatives along z are discretized into a grid of M points with a finite elements method, while the time variable is kept continuous. Equation 5 was thus approximated by a system of M differential equations each describing the time evolution ofũ evaluated on one of the M grid points. To solve Eq. 5, M was chosen in between 1000 and 10000 depending on the complexity of the pinning landscape.
The calculation is let evolve for several periods of oscillation, until the solution converges to a steady-state condition, identified by a close trajectory in the vortex motion phase portrait,u versusũ.
In Fig. 3(a), an example of the vortex motion phase portrait at the surface (z = 0) calculated for surface peak field B p = 50 mT, f p,ave = 50 µN/m, and l = 300 nm is shown. The motion starts with the vortex in the initial condition (ũ(0,z) = 0) and evolves with time until the steady-state condition is reached.
In Fig. 3(b), ment within the pinning landscape is shown. The vortex initiates its motion from the initial condition represented by the straight dashed line at x = 0. Then, it starts oscillating at the surface while moving towards negative x values attracted by a nearby pinning center (darker region). Because of this, the initial motion of the vortex is chaotic-as shown by the phase space in Fig. 3(a) and by the dotted lines in Fig. 3(b)-and eventually converges to steady-state oscillations when the vortex find its equilibrium position attracted by a pinning center, as shown by the solid lines in Fig. 3(b) on the left side of the solution domain. As briefly described above, the self-stabilization of the vortex to equilibrium according to the pinning landscape is allowed by the boundary condition chosen at the far edge of the solution domain. The conditionũ (t,z max ) = 0 allows the vortex to move and reach an equilibrium position without being constrained to a specific location.
Simulations were also performed preparing the initial condition by letting the vortex relaxing to equilibrium before applying the driving force. The employment of such simulation strategy resulted in an equivalent outcome as ifũ(0,z) = 0 was used, but required longer time to converge to a steady-state solution. Therefore, all the simulations here presented were performed using the initial condition reported in Eq. 4.

IV. VORTEX SURFACE RESISTANCE
During a FC experiment, we assume that all the magnetic flux applied during cooldown gets trapped in the material and divided into N v vortexes. Once an electromagnetic field is applied to the material, inducing currents parallel to the surface, vortexes oscillates absorbing part of the electromagnetic energy dissipating power.
The power per unit of surface dissipated by vortex motion is defined as: Considering a situation in which the magnetic field B p and R f l are constant over the surface, R f l can be defined as: where P represents the total power dissipated by the collection of the N v oscillating vortexes. Because of the random nature of the pinning landscape in real materials, we assume that each vortex oscillates within its local pinning landscape, which differs from vortex to vortex except for its average pinning force. Under this assumption, the power P represents the sum of the single-vortex dissipated powers P total number of N v vortexes: withu i (t, z) the i-th vortex velocity and T the microwave field oscillation period. The number of vortexes trapped in the material is estimated by assuming that each vortex carries a single flux quantum φ 0 and that vortexes are equally distributed over the area Σ. In turn, the applied magnetic flux during the cooldown B 0 Σ is divided into N v = B 0 Σ/φ 0 vortexes. Thus Eq. 7 can be rewritten as: where P (1) represents the average of the N v vortexes dissipated power. By applying the same definitions used in Eq. 5, we can rewrite the single-vortex average dissipated power as: A Montecarlo approach was implemented to estimate P (1) by considering a statistical sample of 100 vortexes. In the range of field probed-tens to hundreds of mG-and areas considered-order of the square meter, vortexes exceeds the hundred million unities, making Eq. 10 impossible to compute by means of a desktop computer if the total number of vortexes involved was considered.
The vortex motion equation in Eq. 5 was computed for each of the 100 vortexes generating each time an unique random pinning landscape. P  The colored solid lines represents the single-vortex power P (1) i as a function of B p calculated each time for a random pinning landscape, for a total of 100 vortexes. The solid line with black dots represents instead the average single-vortex dissipated power P (1) .
In Fig. 4(b) and (c) the average phase portrait of vortex oscillation calculated over 100 vortexes, respectively at B p = 0.001 mT and 300 mT, is shown by the thick black solid line. The thin coloured lines represent the individual single-vortex phase portraits. Fig. 4(d) and its inset show instead the calculated normalized vortex surface resistance with abscissa in semi-logarithmic and linear scales, respectively.
At very low microwave fields, the vortex performs oscillations with small amplitude, since well constrained inside a local pinning potential. In this regime, the pinning potential can be approximated by a parabolic function, and the vortex motion equation (Eq. 5) becomes linear. Thus, the calculated vortex response is linear, implying an elliptical phase portrait and constant surface resistance ( Fig. 4(b) and (d)).
At high microwave fields, the vortex surface resistance tends to saturate to a constant value equal to the surface resistance in absence of pinning. In this case, even if each vortex shows a non-linear response, the average phase portrait tends to an ellipse (Fig. 4(c)), meaning that the global response becomes linear for higher microwave fields. Pinning is overcame when vortexes are driven at higher microwave fields and the surface resistance saturates to a constant value equal to the vortex surface resistance in absence of pinning, as shown in the inset of Fig. 4(d).
Such saturation behavior of the vortex surface resistance above a certain microwave field level was also predicted as generated by viscous losses in Ref. 38.
In summary, for low microwave fields pinning can be treated as a parabolic potential-the pinning force is linear and proportional to the vortex displacement, while for high microwave fields pinning becomes negligible. In both cases, the vortex dynamics is approximated to be linear leading to a constant surface resistance.

V. MICROWAVE DEPINNING
Simulations of single-vortex dynamics allow to shed light into the mechanism behind the dependence of R f l against the amplitude of the microwave drive.
In Fig. 5, an example of single-vortex dynamics is reported. The single-vortex power dissipation (P (1) i ) curve shown in Fig. 5(a) was selected from Fig. 4(a). The curve presents a step like behavior (see highlighted points 1 and 2) as a function of B p . Such effect is consequence of a "vortex jump", where the vortex gets depinned from a location and captured by a stronger pinning center nearby, as a consequence of microwave depinning.
Simulations in Fig. 5(b), (c), and (d) show the vortex jump behavior just discussed and the microwave depinning phenomenon. The darker blue line in the phase portraits represents steady-state oscillations, while the  Fig. 5(a).
At B p = 63 mT (point 1), the vortex tip oscillates almost depinned from the local pinning center. As soon as the B p is increased to 75 mT (point 2), the oscillation widens and the vortex starts to interact with a nearby pinning center, which eventually captures it.
Comparing the two simulations (point 1 and 2), appears clear that because of the jump from one pinning center to a stronger one, the amplitude of the steadystate oscillations are reduced and so it is the oscillation velocity. As a result, the power dissipation is reduced accordingly, generating this abrupt step characteristic in the B p dependence of P (1) i . The progressive microwave depinning phenomenon is appreciable by comparing the simulations of point 2 and 3. At B p = 75 mT (point 2) the vortex oscillation is well constrained by the pinning center, while for higher peak magnetic field values (e.g 170 mT, point 3) the oscillation widens substantially reaching areas outside the local pinning potential. If the conditions are right, vortex jump induced by microwave depinning can repeat at higher fields, as shown in Fig. 5(a) around 180 mT.
Each vortex might experience jumps and depinning at different B p values depending on the local pinning landscape, producing an unique P In Fig. 6(a), (b), and (c), the trapped flux surface resistance measured experimentally is compared to several simulations of trapped flux surface resistance calculated for increasing values of average pinning force f p,ave , mean-free-path l, and pinning center size bξ, respectively.
As shown in Fig. 6(a), for larger average pinning force, the B p value at which saturation is reached increases and the overall R f l curve is shifted to lower values. The effect of different mean-free-path is instead shown in Fig. 6(b). The flux-flow viscosity decreases for low l values, the vortex moves more freely and the saturation value increases. The slope of the R f l microwave field amplitude dependence is instead dictated by the dimension of the pinning centers, as shown in Fig. 6(c). The smaller the pinning center, the steeper the trend and the earlier the saturation is reached.
Since the pinning landscape is unknown a-priori, it is virtually impossible to generate quantitative simulations of vortex surface resistance as a function of the microwave field amplitude. Nevertheless, the simulations can qualitatively describe the experimental data. Based on the simulations results, progressive microwave depinning is the most likely culprit for the microwave field dependence of the vortex surface resistance, the higher the driving force, the less constrained the vortexes and wider the oscillations, leading to higher resistance. orHowever, the experimental data does not present saturation at high B p values, contradicting the simulations. This might occur because: i) saturation is reached at higher fields, and ii) others non-linearity play a role in the vortex dynamics. The first possibility is related to the fact that pinning landscape and mean-free-path used in the simulations may not be accurate enough-pinning landscape is unknown a-priori and cavities treated as the one under study are characterized by l varying with depth, with much shorter values at the surface. 52 The second possibility is instead substantiated by the fact that in specific conditions the vortex flux-flow viscosity decreases with the vortex velocity, introducing another non-linearity in the motion equation. Such effect, first described by Larkin and Ovchinnikov, 53 will be discussed in Section VII.

B. Geometrical Considerations
Due to the not planar geometry of the SRF cavity, the vortex surface resistance is not constant across the surface, it instead varies as described previously in Section II.
The temperature variation ∆T out -calculated as the FC and ZFC Tmap data subtraction-is directly related to the surface resistance due to vortexes. The surface resistance is directly proportional to the power dissipated, as shown in Eq. 7, which, in turn, is directly proportional to the temperature variation on the outside surface of the cavity through the thermal transport law, i.e. R f l ∝ ∆T out /B 2 p . Figure 7 reports the normalized data shown in Fig. 2(c) averaged over the board number as a function of the angle θ between B 0ẑ and the microwave currents flowing parallel to the cavity surface. θ was estimated as the angle between B 0ẑ and the cavity surface tangent at each thermometer location.
The error bars are defined to be the standard deviation of the ∆T out /B 2 p average value in ordinate, while for the values in abscissa the error was estimated a-posteriori to be 10 • .
The experimental data is compared against the model here presented, where the Lorentz force sin(θ) dependence is taken into consideration, as well as the local surface current distribution at the cavity surface reported in Fig. 2(b).
The calculation was performed with the same Montecarlo approach described above, where every point is calculated averaging the dissipated power over 100 vortexes. The simulation was carried out assuming l = 100 nm, f p,ave = 200 µN/m, and b = 4, but any other set of parameters returned comparable results. In order to com- pare experimental data and simulation, all the quantities in ordinate are normalized between 0 and 1.
The angle dependence calculated is in good agreement to the experimental data. At θ = 0 • , the surface resistance is zero since the trapped field is parallel to the microwave currents, while at θ 80 • , the surface resistance decreases because the microwave current drops rapidly to zero approaching the cavity irises.

VII. LARKIN-OVCHINNIKOV INSTABILITY
Larkin-Ovchinnikov (LO) instability 53 is a phenomenon of non-equilibrium related to vortex motion. Inside the vortex core superconductivity is suppressed and quasi-particles can occupy virtually any energy level. During vortex motion the energy of quasi-particles increases and when reaches ∼ ∆, quasi-particles diffuse out of the vortex into the superconductor, the vortex core gets depleted, shrinks in dimension, and the fluxflow viscosity decreases.
In the LO theory framework, the critical velocity is defined as [53][54][55] where τ −1 ε is the rate at which quasi-particles relax back into the vortex core. The shorter τ ε , the faster the vortex must travel to meet the instability onset.
Experimental studies on Nb thin films 54 showed that the LO instability critical velocity v * is function of the magnetic field trapped in the sample. The low B 0 region (B 0 ≤ 100 mG) which is of our interest is still poorly explored. However, it was shown that, in Nb, v * increases for very low fields, 54,56 reaching values of the order of couple of km/h for B 0 ∼ 1 G.
The theory developed by Larkin and Ovchinnikov is strictly applicable to T T c . However, the underlying decreasing of the flux-flow viscosity as a function of the vortex velocity is a feature observed also at low temperatures, 55 and therefore relevant in Nb SRF cavities operating at 2 K, especially at high gradients. In order to expand the applicability of such formulation to low temperatures, the quasi-particle relaxation rate τ −1 ε shall be defined as the sum τ −1 s + τ −1 r , respectively scattering and recombination rates. 55 Generally speaking, at low temperatures (T < T c /2), τ −1 s and τ −1 r decrease due to the lower number of phonons and quasi-particles, respectively. 57 However, if τ −1 s is governed by scattering with impurities or lattice defects, it becomes temperature-independent. The occurrence of LO instability should then be more relevant in clean niobium, where the quasi-particle relaxation time is not limited by scattering with impurities and results longer.
In the present study, we discuss trapped flux surface resistance in cavities with low mean-free-path, prepared with a mild baking treatment (Section II), that was proven to increase the concentration of interstitial oxygen and vacancies in the material. [58][59][60] We can then assume that the relaxation rate is most likely dominated by scattering with impurities, and therefore that v * increases with T as expected by Eq. 11.
Hence, if occurring, LO instability would be appreciable at higher microwave fields, introducing deviations from the expected plateau in the surface resistance predicted by this work. However, since no saturation in the experimental data was observed at high microwave field amplitudes, LO instability cannot be totally ruled out, inasmuch saturation might be masked by such phenomenon.

VIII. CONCLUSIONS
In this work, the vortex surface resistance as a function of the microwave field amplitude is simulated by means of a Montecarlo approach aimed to solve the motion equation for a collection of vortexes each oscillating within a randomly defined pinning landscape.
Simulations suggest that at low microwave amplitudes vortexes are well constrained within a local pinning potential and their response is linear, leading to a constant surface resistance. For increasing field amplitudes vortexes may oscillate out of their local pinning potential and possibly jump from their local to a neighbouring pinning center introducing non-linearity. The effect of such phenomena leads to an increase of the vortex surface resistance until the average response approaches a linear behavior for higher field amplitudes, where the average vortex motion becomes insensitive to the pinning landscape. When this happens the surface resistance saturates to a constant value.
Slope, saturation value, and saturation onset are closely related to the pinning landscape-not known apriori-within which vortexes oscillate. Quantitative description of R f l as a function of B p is then virtually impossible. Nevertheless, the experimental data are in qualitative agreement with the simulations especially the data as a function of the angle between microwave current and vortex magnetic flux direction.
The experimental data as a function of the microwave field amplitude do not show saturation as instead expected from simulations. As described above, this might be due to the not enough accurate choice of pinning parameters and mean-free-path value in the simulations, or to the onset of LO instability at high fields.
Concluding, this work allowed to gain insights into vortex dynamics at high amplitude microwave drive and on the transition from low to high field amplitudes. The microscopic nature the microwave field dependence of the vortex surface resistance was presented and attributed to microwave depinning of vortexes. Appreciable contributions from thermal feedback effects for the magnetic field values and microwave amplitudes explored experimentally-B 0 < 100 mG and B p < 200 mT-were instead ruled out.