Nonlocal electrodynamics of normal and superconducting films

Electrically conducting films in a time-varying transverse applied magnetic field are considered. Their behavior is strongly influenced by the self-field of the induced currents, making the electrodynamics nonlocal, and consequently difficult to analyze both numerically and analytically. We present a formalism which allows many phenomena related to superconducting and Ohmic films to be modelled and analyzed. The formalism is based on the Maxwell equations, and a material current-voltage characteristics, linear for normal metals, and nonlinear for superconductors, plus a careful account of the boundary conditions. For Ohmic films, we consider the response to a delta function source-field turned on instantly. As one of few problems in nonlocal electrodynamics, this has an analytical solution, which we obtain, in both Fourier and real space. Next, the dynamical behaviour of a square superconductor film during ramping up of the field, and subsequently returning to zero, is treated numerically. Then, this remanent state is used as initial condition for triggering thermomagnetic avalanches. The avalanches tend to invade the central part where the density of trapped flux is largest, forming dendritic patterns in excellent agreement with magneto-optical images. Detailed profiles of current and flux density are presented and discussed. Finally, the formalism is extended to multiply connected samples, and numerical results for a patterned superconducting film, a ring with a square lattice of antidots, are presented and discussed.


Introduction
The flux dynamics in electrically conducting films experiencing a time-varying transverse applied magnetic field is governed by the Maxwell equations, with the material characteristics supplied as an additional relation between the electric field and current density. The systems to be addressed in this work range from superconductors to Ohmic materials. To solve these equations it is necessary to determine the currents induced in the film as the magnetic field varies. This is a nontrivial task since one must also account for the significant self-field of the induced currents, which makes the final relations nonlocal [1].
The electromagnetic behaviour of type-II superconductors is often well described by Bean's critical-state model [2]. For bulk samples initially zero-field-cooled below the transition temperature, T c , and then exposed to an increasing applied magnetic field, H a , the model tells that the material sets up lossfree shielding currents of critical density, j c . This current flows in the same macroscopic regions as where the magnetic flux is allowed to penetrate, while the inner unpenetrated part remains free of currents. In films, on the other hand, the electromagnetic non-locality implies that induced currents flow in the entire sample [3,4]. Thus, the film behaviour is qualitatively different from that of bulks, and magneto-optical imaging (MOI) of thin superconductors has revealed strong piling up of the magnetic field around the sample edges, where values far above H a are reached [5]. At internal boundaries, such as the inner edge of a planar ring, the field can, due to the nonlocal electrodynamics, be in opposite direction of the applied field [6,7]. Strongly modified behaviour is found also in films patterned with regular arrays of small holes (antidots), which tend to guide the flux into the superconductor [8,9,10,11].
The response of Ohmic films exposed to varying transverse magnetic fields is also described by nonlocal electrodynamics, but here the material responds linearly. Numerical solutions for strip and disk geometries have shown that the combination of non-locality and dissipation causes a rapid penetration of a suddenly applied magnetic field [12,13]. Different from superconductors, even regions deep inside an Ohmic film are quickly penetrated by the magnetic field.
A phenomenon that involves both the critical-state and Ohmic properties is the occurance of flux avalanches or flux jumps. These are commonly observed in type-II superconductors at low temperatures, and are caused by a thermomagnetic instability which drives the superconductor from the critical-state to a high resistivity state [14]. The instability is triggered, e.g., by a small temperature fluctuation which reduces the flux pinning locally, and some quantized flux lines, or vortices, will start moving. This creates local heat dissipation and the temperature will increase even further, thus forming a positive feedback loop. The result can be an exponential growth in the temperature and a large-scale runaway of magnetic flux. In superconducting films the thermomagnetic instability is seen by MOI to manifest as abrupt avalanches of magnetic flux which form complex branching filamentary structures, socalled dendritic avalanches [15,16,17,18,19,20,21].
These avalanches can be modeled using the equations describing nonlocal and nonlinear electrodynamics coupled with an equation for the production and propagation of heat [22]. Linearization of the equations has been highly succesful in parametrizing the conditions for onset of the instability, confirming that the nonlocal electrodynamics makes a significant difference between bulk [23] and film geometries [24,25,19,26,27]. Numerical simulations of the full time evolution of the avalanches Figure 1. An electrically conducting film in time-varying transverse applied magnetic field Ha(t). Due to the induced current, J, the sample partly expels the transverse field component Hz.
have produced dendritic flux patterns in excellent agreement with the experimental MOI results [26,28]. The propagation of these avalanches is extremely fast -velocities up tp 180 km/s have been measured [29], and the process is driven by adiabatic heating [30]. During an avalanche the local temperature is expected to rise above T c , thus bringing for a very short time interval the superconductor to the normal conducting state. In such cases, the process is governed by the interplay between supercurrents and Ohmic normal-state currents.
In this work, we consider the electrodynamics of normal and superconducting films in transverse applied magnetic field; see figure 1. The basic idea is that a wide range of physical problems in this geometry can be described by the same formalism based on the Maxwell equations, only by supplying a relation between electric field E and sheet current J to characterize the material. Thus, we will describe the formalism in detail, with particular focus on enforcement of the boundary conditions. Having described the formalism, we apply the solution method to various physical problems. As a first example, we calculate the remanent flux in a superconducting square, after the applied magnetic field has been increased to reach full penetration, and then decreased back to zero. Due to the nonlocal electrodynamics, the field and current distributions in the square are highly nontrivial. From the remanent state, we consider the evolution of dendritic flux avalanches, which means that we must model the propagation of heat in the system, in addition to the electrodynamics. Our numerical solution is compared with a magneto-optical imaging experiment which maps the magnetic flux distribution in a NbN superconductor in descending magnetic field. As a separate problem, we considers the response of an infinite Ohmic sheet to a delta function source field. This problem is analytically solvable and the solution sheds light on the dynamics of Ohmic films, as well as the dynamics of dendritic flux avalanches, which is driven by a normal domain invading a superconducting phase. Finally, we consider a superconducting ring patterned with a regular array of antidots. This system is interesting due to the conflicting symmetries of the sample and the antidot array, but rather difficult to handle numerically due to the complicated sample layout. In total, all these problems demonstrate that our formalism is powerful and flexible, as it can be applied to a wide range of physical problems in the thin-film transverse geometry. This paper is organized as follows. Section 2 describes the transverse geometry. Section 3 finds the remanent flux distribution in a square superconductor, and calculates also the analytical solution for the field and currents in a normal metal film subjected to delta function source field. Section 4 considers dendritic flux avalanches in the remanent state, both numerically and by a magneto-optical imaging experiment. Section 5 considers the dynamics of a superconducting ring patterned with antidots. Finally, section 6 provides the conclusion.

Connecting magnetic field and current distributions
A key element in solving magnetic flux dynamics problems in films placed in a transversely applied field, H a , both Ohmic and superconducting ones, is the relation between the distributions of electrical current and transverse magnetic field H z (x, y) over the (x, y)-plane defined by the film. To establish the formalism used in this work, we assume the film thickness to be much smaller than any length characterizing the patterns of flux and currents. The current density in the film can then be expressed as where J is the sheet current. It is convenient to introduce the local magnetization g = g(x, y) as whereẑ is the unit vector transverse to the sample plane. The total magnetic moment of the film can then be expressed as Neglecting the displacement field, the Ampère law becomes and Fourier transforms along the Cartesian axes give ik y H [3] z − ik z H [3] y = ik y g [2] , −ik x H [3] z + ik z H [3] x = −ik x g [2] , ik x H [3] y − ik y H [3] x = 0. Here H [3] = H [3] (k x , k y , k z ) is the three-dimensional Fourier transform of H and g [2] = g [2] (k x , k y ). Conservation of magnetic flux, ∇ · H = 0, yields ik x H [3] x + ik y H [3] y + ik z H [3] z = 0, so that where k = k 2 x + k 2 y . Thus, H x is nonzero, and the same holds for H y , which is a general feature of films in the transverse geometry. Isolating H [3] z gives H [3] z = k 2 k 2 z + k 2 g [2] , and inverse Fourier transform in z direction results in the final expression where H [2] z = H [2] (k x , k y , z). For inversion, e.g., of magneto-optical images [31,32,33] one often uses a finite z to account for a small gap between the sample and the field sensing, i.e., Faraday rotating layer. However, for the flux dynamics calculations in this work we only consider the expressions at z = 0.
The H z − g relation will henceforth be denoted as the Biot-Savart law, and it can be written as [34] where F and F −1 is forward and inverse Fourier transform, respectively. The inverse relation is equally simple, The above equations are exact on an infinite sheet. For films having a finite area, they are good approximations for short wavelengths [1].

Iteration scheme
Consider a planar conducting film surrounded by vacuum, and with H z known inside the sample area defined by its boundary. Given the task to determine the local magnetization, g, the intuitive approach is to use (7). However, this fails to give correct result unless H z is known over the entire plane. An approach allowing g to be found correctly was invented by Brandt [1,35], and is based on a matrix inversion scheme. The approach proved to work very well for simple geometries which can be represented by a fairly small number of discrete grid points. Later, the numerical performance of the matrix inversion was improved by using a congruent gradient method [36,37]. An alternative approach is to try and extrapolate H z (x, y) to the outside area, and then apply (7). For an infinitely long strip, this can be done by symmetry considerations, as shown by Aranson et al. [26]. In this work we consider far more general geometries, and will calculate H z making use of the fact that outside the sample one has g = 0. Our scheme is iterative, and as will be demonstrated, computationally efficient [28].
To describe our approach, it is convenient to define a function representing the projection on the sample, S(x, y) = 1, inside the sample boundary, 0, outside the sample boundary.
The corresponding projection on the outside region is 1 − S(x, y). For brevity, the argument (x, y) is omitted in the next expressions. The iterations start by defining a trial function H   The local magnetization, the correct and anticipated one, is then expressed respectively as is an improved approximation to H z , and we repeat the whole procedure s times untill (1 − S)g (s) becomes vanishingly small. In this case g (s) gives the correct magnetization distribution, and we have successfully inverted the Biot-Savart law.

Test case: Array of superconducting strips
In order to illustrate the iteration scheme, the algorithm will first be applied to a reference case with known analytical solution. We consider a periodic arrangement of infinitely long superconducting strips in the Bean critical-state, where an exact solution was obtained by Mawatari [38]. The configuration is seen in figure 2, where three strips in an infinite array are shown. Each strip has width 2w, thickness d and center-to-center distance 2L. Due to the magnetic field applied in the z-direction, the magnetic flux penetrates from both sides of the strips. For the strip centered at x = 0, the flux front position is at |x| = a, and the magnetic flux distribution is given by where H c = J c /π. The corresponding sheet current is where the function ϕ(x) is The width of the fluxfree area, 2a, shrinks with the increasing applied field according to Let us now assume that the magnetic field distribution, equation (13), is known over the area of the strip, |x| < w, and based only on that, set out to determine both the sheet current, J, and local magnetization, g. We use (10) and iterate over an area 2L × 2L, discretized on a 256 × 256 equidistant grid. The calculations were performed using L = 1.5 and H a = 0.5 in units where J c = w = 1. As initial guess we set H The results obtained after 0, 1 and 5 iterations, are presented in figure 3. In spite of a poor initial guess for the outside field, already after one iteration, the result is very much improved. The largest deviation is that significant currents flow in the region between the strips. However, after 5 iterations this unphysical feature is negligible, and the numerical and exact solutions are practically the same. Thus, we conclude that our iterative inversion scheme is very rapidly converging towards the exact solution in this non-trivial test case.

Superconducting films
Consider now the more general situation where a superconducting film of finite size is experiencing a time-varying transverse homogeneous applied magnetic field, H a (t). We want to calculate numerically the electrodynamics as the field is gradually changing. In such cases, electrical currents will be induced in the sample, setting up their own magnetic self-field. The total transverse field, H z , has therefore two contributions, where the last term represents the induced field (6). Taking the time derivative and inverting this equation, one getṡ Outside the sample,Ḣ z is found by a boundary condition, as described in section 2.2. Inside the sample area,Ḣ z is found using the Faraday law, µ 0Ḣz = −(∇ × E) z , which combined with a material law E = ρJ /d giveṡ where the resistivity ρ represents the material characteristics. The conventional material characteristic used to describe a superconductor in the slow dynamics, or flux creep regime, is a power law where ρ 0 is a resistivity constant, H c2 is the upper critical field, and J c the critical sheet current. The exponent m is typically small, while the creep exponent n 1. For high-T c superconductors, e.g., YBa 2 Cu 3 O x , one commonly finds n = 10 − 70 [39,40], while for MgB 2 exponents as high as n = 78 were found at T = 25 K [41]. In conventional superconductors flux creep is not observed unless very close to T c , so in simulations one may then set n sufficiently large to make creep negligible.
In this work we present simulation results for both stable flux creep dynamics, and for the far more dramatic flux avalanche dynamics. We illustrate first the numerical scheme by applying it to the smooth dynamics when the applied field is ramped from zero and up to a value giving essentially full flux penetration, and then back again. This produces a remanent state which contains trapped flux, and is the state used in Section IV as starting point for simulations of avalanches.
To solve the dynamical equations numerically we convert them to dimensionless form, assuming that J c and |Ḣ a | are constants. Based on the sample half-width w, and the parameter we choose dimensionless quantities asQ In these units the ramp rate satisfies |dH a /dt| = 1, and the only free parameters are the exponents m and n. We will henceforth omit the tildes when writing the dimensionless quantities. The sample of size 2×2 is embedded in a 2.6×2.6 square which is discretized on a 512 × 512 equidistant grid. The results are obtained Shown in figure 4 (left) is the flux distribution and current stream line patterns after increasing the applied field to H a = 1. The field along the sample edge is much higher, reaching values nearly twice as large close to the mid-point of the sides. The flux front reaches almost to the center of the sample forming a flux density pattern often observed in MOI experiments [5]. The current streamlines are in the fluxpenetrated regions essentially equidistant, and display sharp turns at the diagonals, as typical for a square sample in the critical state [1]. A slight overall convexity of the streamline loops is due to the finite creep exponent.
The two panels on the right show the state after the applied field was ramped down to zero. In this remanent state the edge field is reversed. As seen in the current map below the regions of maximum current are now tongue-shaped extending from each side of the square. Here, the nearly equidistant stream lines represent current flow in opposite direction as compared to those in the left panel. Only in the central part of the sample is the J circulating the square in the same direction as at maximum applied field. However, the magnitude J is in a large central area far below J c , and the current flows in a different pattern. The result is in good agreement with scanning Hall probe measurements [42], MOI and previous numerical simulations [43,44].
We return to this remanent state, when reporting simulations of avalanche dynamics. Since the transient electromagnetic behavior during flux avalanches in superconductors involves rapid localized variations in the field taking place in normally conducting regions, we present next, as reference, a useful exact solution to a generic dynamical problem for an Ohmic film.

Ohmic films
When a uniform magnetic field is suddenly applied transverse to a normally conducting film, electrical currents will be induced everywhere in the specimen. This global character of the response has similarities to that of superconducting films. The case considered here, is a delta function source-field applied instantaneously to an infinite sheet of normal conductor, with resistivity ρ 0 . Let the applied field be described by where H 0 is the field strength, δ 2 is the two dimensional delta function, and Θ is the Heaviside step function. The dynamical response is described by (18), combined with Ohm's law, E = ρ 0 J /d. One then getṡ whereQ −1 is the inverse Biot-Savart operator (7) and v 0 = ρ 0 /(dµ 0 ) is a constant of dimension velocity. Fourier transforms yield where g [2+1] is the Fourier transform of g in two spatial dimensions plus time. Isolation of g and transforming back to time domain for t > 0 gives This means that the eddy currents and magnetic fields decay with characteristic time where l = 2π/k is the wavelength. Thus, the longest decay times are found for the largest wavelengths. Note that the characteristic time for films is shorter by a factor d/l compared to the slowest decaying modes in bulk Ohmic samples [45]. Interestingly, (26) gives results in fairly close agreement with the numerical evaluation of relaxation times after a uniform field is abruptly applied to conducting strips and disks [12,13]. This shows that the decay time is only weakly sensitive to the spatial profile of the applied field excitation as well as the shape of the Ohmic film. Inverse Fourier transform in space of (25) yields the result At t → 0 + , the self-field is proportional to a delta function, which means that it shields exactly the applied field. At all times, (27) conserves flux, since d 2 rH z = 0. . An Ohmic film exposed to a delta function source-field applied at t = 0. At t = 0 the self field Hz − Ha shields the applied field completely, while at times 0.01, 0.1 and 1 the shielding is gradually reduced due to the decay of the shielding currents Jϕ.
The corresponding decaying sheet current is given by The sheet current has a maximum at i.e. the peak moves with a constant velocity √ 2v 0 , similar to the eddy current front in disks after a uniform magnetic field is instantly applied [13].
Shown in figure 5 are the spatial profiles of the self field, H z −H a , and the shielding current, J ϕ , plotted at times t = 0.01, 0.1, and 1, in units where H 0 = v 0 = 1. In the beginning the self field is focused near x = 0, and with time it decays and becomes increasingly uniform. The result shows that (23) produces a solution that is very different from that of a diffusion process. In particular, there is no well defined diffusion front, since both magnetic field and currents decay algebraically as J ∼ 1/r 2 and H ∼ 1/r 3 at large r. This contrasts the heat kernel solution of the ordinary diffusion equation, which decays exponentially.
Although this calculation is exact only for an infinite film, the behavior of the shorter wavelengths should be reasonable approximations for finite normal domains. One such example is the normal parts in the center of dendritic flux avalanches. The time τ and velocity √ 2v 0 are thus also expected to be characteristic for the propagation of dendritic flux avalanches.

Dendritic flux avalanches in superconductors
Dendritic flux avalanches appear in descending as well as ascending magnetic field. Since avalanches in descending field appear on a highly nontrivial background, it is interesting to study their properties. However, a numerical simulation of the full process, from initially zero-field-cooled state, to the reentrant stability at full penetration, and then back down, with the thermal feedback turned on, is computationally demanding. Therefore, we will here consider a simpler scheme, where the remanent state is prepared with the thermal feedback turned off, i.e., we assume that there are no avalanches, the temperature is everywhere T = T 0 , and the flux and current distributions are as described in figure 4. From this background, we will explicitly trigger the dendritic flux avalanche by a heat pulse near the edge, and consider its development.
We will in this section describe the equations governing the flow of heat, the solution method, the units, and how to rescale the previous result for the remanent state to these units. We will consider the time evolution of dendritic flux avalanches nucleated at two different locations. For comparison, we show the flux distribution of a superconductor in descending field, mapped by the magneto-optical imaging method.

Preparation of the pre-avalanche state
Due to motion of vortices there is heat dissipation in type-II superconductors experiencing a varying external magnetic field. Since many of the material parameters, most notably J c and n, depend on temperature, the dissipation will interfere with the electrodynamics. Thus, in order to get a complete description of the dynamics it is necessary to model, in addition to the electrodynamics, also the propagation of the produced heat.
Consider a superconducting film in thermal contact with a substrate of constant temperature T 0 . The propagation of heat can then be described by the equation where c is the specific heat, κ is the lateral thermal conductivity, h is the coefficient of heat transfer to the substrate, and the last term represents the Joule heating.
To transform (30) into a dimensionless form, one needs to decide on convenient scales for normalization. The most difficult scale to decide is the time scale, since our problem is composite, with physical processes at many different time scales. Under the assumption that a dendritic flux avalanches mainly propagates due to a domain in the normal state invading a superconducting domain, it is natural to chose a time scale appropriate for the decay of normal currents, as discussed in section 3.2. Thus we let t = tρ 0 /dwµ 0 , where ρ 0 is the normal resistivity of the superconductor at T c . This means that the timet ∼ 1 is characteristic for the decay of modes with size 4πw. For the other quantities, we letT = T /T c ,J = J/dj c0 , andẼ = E/ρ 0 j c0 , where j c0 is the critical current density at T = 0. Since this set of units is appropriate for describing the fast decay of normal currents during the propagation of the dendritic flux avalanches, the rate of change of the applied field is typically very small in comparison, i.e., dH a /dt 1.

The heat propagation equation then becomes
Here α is dimensionless heat conductance, β is dimensionless constant for heat transfer to the substrate, and γ is a Joule heating parameter. These constant parameters are defined as where the material parameters at the right-hand-sides are evaluated at T c . The temperature-dependence of γ is taken asγ(T ) = c(T c )/c(T ). In this work only the phonon contribution to c, givingγ =T −3 at low temperatures, is taken into account.
We have also assumed that the fractions κ/c and h/c are temperature-independent. Henceforth we will skip the tildes when reporting the results in dimensionless units.
In order to simulate the thermomagnetic instabilities one must specify temperature dependencies of J c = dj c and n in addition to the thermal parameters α, β, and γ. We let The resistivity is Equation (34) describes a flux creep regime at J < J c and T < 1, normal resistivity at T > 1, and a high-resistivity flux flow regime at J > J c . The latter implies that we have taken into account the flux flow instability [46,47]. This instability, which must not be confused with the thermomagnetic instability, makes the flux flow non-linear at high electric fields. When the vortex velocity is higher than a critical value, v * , the resistivity jumps from the usual flux resistive ρ 0 H/H c2 to the much higher value ρ 0 .
In (34), we have assumed that inside the avalanche, the flux motion satisfies v > v * , wherever J > J c . The simulation of the evolution of dendritic flux avalanches will be based on the remanent state of figure 4. It must then be transformed from units whereḢ a = 1 to units withḢ a 1. The relevant conversion factors are The physical quantities will then transform as g → ug, J → uJ, H → uH, t → vt, and E → uE/v. In the runs, the dimensionless parameters that characterize the thermal properties of the sample, equation (32), are selected as α = 2 · 10 −5 , β = 0.05, γ = 10. For the ramp rate and substrate temperature, we choseḢ a = 10 −10 and T 0 = 0.2. This gives J c = 0.8, and n = 100, when n 0 = 20. The conversions factors become u = 0.64 and v = 8 · 10 9 .
The spatial disorder usually present in superconducting films manifests itself in a non-uniform J c0 . Hence, a random disorder is added to model by assigning each grid point with J c0 → 1 + ∆(r − 1/2), where ∆ = 0.05 and r ∈ (0, 1) are random numbers.

Symmetric nucleation
Based on the remanent state, an avalanche is nucleated centrally at one side of the sample by assigning T = 1.5 in a small area close to the edge. Figure 6 shows H z , J, and T at times t = 1, 5 and 40 after the nucleation.
At t = 1 (left column of figure 6) only the critical state region is affected, and the avalanche is mainly visible in H z and T as a long, thin filament with some tendency of branching. Some of the flux is negative, which means that the avalanche partly consists of positive flux leaving, partly of negative flux entering the sample. The heating is significant, with most of the avalanche already heated above the critical temperature. Yet, the tip is still superconducting, in a flux-flow state with high resistivity. The effect on the sheet current J is less visible, although the value drops locally inside the avalanche. As typical for the remanent state, the direction of the J along the edge is such that it favors positive flux leaving and negative flux entering the sample. At t = 5 (middle column) the avalanche spreads out into the inner parts of the sample. At this stage the avalanche prefers to invade the regions with highest flux density. The explanation of this behaviour is in the sheet current pattern, where the branch tips are seen to propagate transverse to the current stream lines, i.e., in the direction of the Lorentz forces density F L = µ 0 H zẑ × J . At the same time, due to the nonlocality of the equations, the propagating avalanche distorts the current density in a large portion of the sample.
At t = 40 (right column) the avalanche has essentially reached its largest extent, and due to the efficient heat removal to substrate, the branches are now colder. Because of the symmetric nucleation the avalanche is almost symmetric, but not entirely, since the state prior to the avalanche was seeded with randomly distributed disorder.
The avalanche is large and destructive as it affects the distribution of flux and currents in the entire sample. Another most dramatic effect is the strong change of the critical state region around the edge. Before the avalanche took place, the state was just as described by the critical state model, with constant current density and stream lines with almost equal spacing starting from the edge. After the avalanche the critical state has vanished completely, leaving a current density which is less than the half of the original value and stream lines that are no longer parallel. This means that the consequences are a lot more severe for the avalanches in the remanent state than in ascending field, where the critical state is destroyed only in the vicinity of the avalanche [30].
Worth noticing is also that at t = 40, there are small, embryonic avalanches appearing close to the edge at both sides of the large avalanche. However, due to the above mentioned destruction of the critical state, these are unable to develop into full avalanches, and therefore remain small.

Off-center nucleation
Here we investigate how the evolution of dendritic avalanches depends on the location where it is initiated. We explore this by nucleating an avalanche away from the center of the side of the square. We use the same remanent state and disorder configuration as for the symmetrically nucleated avalanche in figure 6.
The results of such an asymmetrically nucleated avalanche is shown in Figure 7, with H z , J and T obtained at t = 1, 5, and 40. The avalanche is nucleated close to the upper left corner, and spreads out and fills nearly the whole inner part of the sample. The size, shape and time evolution of the avalanche shows much resemblance with the avalanche in figure 6, but the symmetry of the final state is entirely different. The final state looks like a loop, also in this case, but it closes on the bottom right corner.
For the off-center triggered avalanche, all the main features discussed for the centrally triggered avalanche are present: the enormous size, the extensive spreading into the regions with highest flux density, the negative flux inside the avalanche, the destruction of the critical state, and finally the appearance of embryonic avalanches at the edge.
Some profiles of H z , J and T along the y-axis (vertical line through the center of the square) at times 1, 5 and 40 are shown in figure 8. At t = 1 all profiles are as expected for the remanent state in the critical state model. In H z (y) at t = 5 one sees the finger-like structures penetrating the places where the flux density was highest. The fingers consist of positive flux, while at t = 40 there is also significant amounts of negative flux in the avalanche. The overall |H z (y)| after the avalanche, both inside and outside, is much closer the zero than the state prior to the avalanche. The J x (y) profiles are complex as the currents of the fingering structures go in opposite directions on each side of the fingers. More than anything, the J x (y) shows that after the avalanche event the critical state has vanished completely. Moreover, there is no clearly preferred direction of the current. E.g., one sees that close to the edge there is a thin layer with reversed current direction. The T profiles at t = 5 shows individual hot branches with temperatures just below T c . At t = 40 it is no longer possible to distinguish the different branches as the thermal diffusion has smeared the temperature profiles.

Magneto-optical imaging of avalanches
In order to validate the correctness of the numerical solution of the dendritic flux avalanches in the remanent state, magneto-optical imaging experiments were performed. The sample was a 180 nm thick NbN superconducting film shaped as a square of sides 5.35 mm. Placed on top of the sample was an in-plane magnetization ferrite garnet film used as Faraday-rotation sensor [48]. Since the Faraday rotation increases monotonously with the perpendicular component of magnetic field, one can by polarized light microscopy create a map of the magnetic field distribution above the film [31]. The sample was initially zero-field-cooled to 4 K and magnetic field was applied perpendicular to the film. During the ascending field ramp, there were many avalanches, but due to the reentrant stability in high fields, the full penetration state at 17 mT was critical state like [20]. In descending field, the flux dynamics was for a long time smooth and at 10.5 mT the flux distribution was as shown in the small image in figure 9. Then, suddenly, a large avalanche stroke, and in a short time, it entered a large portion of the sample. This large avalanche, seen in the main image of figure 9, is typical for avalanches in descending field, near the avalanche threshold temperature [16].
Although the avalanche did not strike exactly in the remanent state it is close enough to be used in a qualitative comparison with the simulation. First we note that the avalanche has a clear similarity to the simulated flux avalanche in figure 7, as it avoids the critical state region close to the edges and instead it invades the region with highest magnetic flux density. The size and extent of the avalanche is also similar. The majority of branches are dark meaning that the flux density is low. Some of the branches are white. In this case it is not clear if this means negative flux, as was reported in the simulations, or positive flux, since the image only shows the absolute value of H z . One more detail worth noticing, is the appearance of embryonic avalanches at the edge of the sample, just as predicted by the numerical simulations. The main discrepancy between the flux distribution of the simulation and the magneto-optical experiment is the width of the branches. In the experiment they are much more narrow than in the simulation. This is an indication that the NbN film has lower value of the effective heat diffusion parameter α, given in (32), than what was used in the simulation. However, due to the limited spatial resolution one cannot run the simulations with smaller values of α without at the same time increasing the number of grid points.

Film with antidots
The formalism described in section 2 for modelling the dynamics of thin-film superconductors in transverse field, is valid only for simply connected samples. We will now extend the formalism to multiply connected samples. This gives us the opportunity to study also the flux dynamics of superconducting films with antidots (non-conducting holes). This is of interest since, due to the nonlocal electrodynamics, the presence of antidots may strongly influence the distribution of flux and current in the films. For example, it has been reported that patterning with regular arrays of antidots makes the magnetic flux penetration anisotropic [8,10]. Currently, there are only few numerical simulation works that has considered the critical state flux penetration in samples patterned with antidots [49,50,9,51,11,52]. In order to improve the theoretical knowledge on the field, we will here consider the numerically challenging sample configuration of a superconducting ring patterned with a square array of disk-shaped antidots.
Let us first consider the boundary conditions. When the film contains holes of any shape and number, their presence can be implemented by an iterative scheme similar to that described in section 2.2. For each hole, labeled α = 1 . . . N , we then define the hole projection Equation (10) now becomes where which allows H z in each hole to be reconstructed. The constants C (i) α are fixed by the flux conservation condition The operatorQ α can be any implementation of the forward Biot-Savart law. In general, it is beneficial to use different implementations for large and small holes. For large holes, the best is to letQ m =Q, i.e., the full Biot-Savart law, equation (6). The drawback of this approach is that it runs over all grid points. The advantage is that the linear operatorQ can be moved outside the sum in (38) when there are more than one large hole. For small holes one can use an implementation ofQ which for each hole only loops over the grid points in the hole. For convergence of the procedure the input to the operatorQ n should first be shifted to minimize contributions from the edge of the hole. This will reduce the damaging effect of the sharp cut made by h m .
Let us consider the flux penetration in a superconducting ring with antidots patterned in a square grid. This layout allows us to illustrate the consequences of electromagnetic non-locality and non-trivial dynamics given the conflicting symmetries of sample and the array of antidots. In units where the outer radius is R = 1 anḋ H a = 1 (same as the remanent state of section 3.1), the inner radius is 0.8, and the antidots, 385 in total, have radii a = 0.013. The center-to-center distances of the antidots are 4a. In order to apply the boundary conditions, the ring is embedded in a square of size L x = L y = 1.3, which is discretized on a 1024 × 1024 equidistant grid.
The left panel of Figure 10 shows the flux distribution at H a = 0.2 when flux has fully penetrated the ring, starting to fill the central hole with flux. The outer edge is white indicating high flux density and the inner edge is dark indicating negative flux, as typical for the ring geometry [6]. This means that the currents flow in a clockwise direction everywhere in the ring, contrary to a strip where currents flow in both directions. The local flux distribution inside the sample is much distorted due to the presence of the antidots. The current stream lines has to bend around the antidots and this induces large amount of flux in the antidots. The flux is negative towards the outer edge and positive on the other side of the antidot. For the holes closest to the inner edge of the ring the situation is opposite. This means that the inner edge to large extent behaves like an outer edge subjected to a negative applied field.
The four right panels show close-up views of H z and J around 45 • and 90 • direction. The J-maps show that there are connected critical state region with J ≈ 1 extending from the outside to the inner edge of the ring. The critical state connected regions follow the symmetry of the antidot lattice and these act like channels for easy flux penetration [50,11]. Hence, the antidot lattice makes the flux penetration anisotropic, in good agreement with previous magneto-optical experiments on superconducting disks patterned with antidots [8,10].
Between the antidots there are places where J < 1. This is a feature that cannot be predicted by the Bean model and it shows that it is necessary to solve the time-dependent equations to get a correct description of the state. These sub-critical pockets may be of technological relevance because they imply that sample patterned with periodic antidot arrays has better shielding properties for local magnetic fields than unpatterned samples.

Conclusions
The macroscopic electrodynamics of thin films, either superconducting or Ohmic, in transverse applied field can be modeled by the Maxwell equations. The formalism is capable of handling a wide range of physical systems, where the material-specific properties are introduced as an E − J relation, which is linear for Ohmic conductors, nonlinear for superconductors. A challanging point in the formalism is to calculate the currents for a known distribution of the magnetic field. We solve this problem by a hybrid real space -Fourier space iterative scheme, which is both computationally efficient and is able to handle also samples with non-symmetric boundary.
When magnetic field is increased to reach full flux penetration and then decreased to zero, superconductors with strong flux pinning experience that a large amount of remanent flux is trapped inside the specimen. Both the distributions of current and magnetic field in this remanent state are highly nontrivial, as we showed by a numerical simulation on a film with square shape. In order to consider how dendritic flux avalanches evolve on the background of the remanent state, we developed the formalism for rescaling solutions and for calculating the flow of heat. The dendritic flux avalanche in the remanent state was found to develop as an irregular branching structure that enters the inner parts of the sample. The avalanche consisted partly of positive flux leaving the sample, partly of negative flux entering. It was found to be more destructive than avalanches in ascending field since, after the avalanche, the critical state had vanished completely from the entire film. The spatial extent of the avalanches was sensitive to the nucleation position, but the size and overall consequences were not. A magneto-optical imaging experiment on dendritic flux avalanche in descending field in a NbN film showed similar looking gigantic flux avalanches and supported the findings of the simulations.
Very few problems related to the nonlocal electrodynamics of thin films are analytically solvable. An exception is the response of an infinite Ohmic film to a delta function source field, turned on instantly. We calculated this solution, both in Fourier and real space, and the solution gave much insight into the behaviour of Ohmic films or Ohmic domains in transverse field. We found that, in Fourier space, the mode of wavelength l decays with characteristic time µ 0 dl/(4πρ 0 ). The solution in real space showed that there was no well-defined front of propagation, since both current and magnetic self-field decreased algebraically far from the source. Yet, the current had a maximum moving away from the source with constant velocity √ 2ρ 0 /(dµ 0 ). Finally, we generalized the numerical simulation formalism from simply connected to multiply connected geometry, i.e., we allow the samples to contain non-conducting holes. Thereupon, we consider the magnetic flux penetration in a superconducting ring with antidots (small holes) distributed in a square array. The magnetic flux distribution was locally much perturbed by the antidots, and also the large scale flux distribution was modified, as it became anisotropic when the magnetic flux was guided along the directions of the antidot array. Between the antidots there were localized regions with low flux traffic and J < J c . This is contrary to the situation in simply connected samples, where the regions with J < J c usually form large connected domains. The current distribution inside the ring patterned with antidots was thus highly nontrivial, even in the critical state.
In summary, we have shown that a wide range of apparently different phenomena related to the electrodynamics of superconducting and Ohmic films in transverse field can be described by one formalism based on the Maxwell equations and material-ofQ by a multiplication by a Gaussian in Fourier space, In real space, this implies a convolution with a Gaussian Φ σ (r) = 1 i.e., the result is an interpolation with a neighborhood of size σ. It is reasonable to let σ be a small number of order the grid size σ ∼ 2a/N x . Note that lim σ→0 Φ σ = δ(x)δ(y).
The same Gaussian smoothing should also be applied toQ −1 .