Modeling and fabrication of chip-based superconducting traps for levitation of micrometer-sized superconducting particles

We describe the finite-element modeling and fabrication of chip-based superconducting traps for levitating micrometer-sized superconducting particles. Such experiments promise to lead to a new generation of macroscopic quantum experiments and of force and acceleration sensors. An accurate modeling of the utilized trap architectures is crucial for predicting parameters of the traps, such as trap stability, frequency and levitation height, in realistic situations accounting for the finite extent of the involved superconducting objects. To this end, we apply a modeling method that is applicable to arbitrary superconducting structures in the Meissner state. It is based on Maxwell-London equations in the static regime using the A-V formulation. The modeling allows us to simulate superconducting objects with arbitrary geometry subject to arbitrary magnetic field distributions and captures finite volume effects like magnetic field expulsion. We use this modeling to simulate two chip-based trap architectures: an anti-Helmholtz coil-type trap and a planar double-loop trap. We calculate important parameters of the superconducting traps for the cases of levitating micrometer-sized particles of either spherical, cylindrical or ring shape. We compare our modeling results to analytical test cases for idealized geometries. We also model detection of the motion of the levitated particle by measurement of flux-changes induced in a nearby pick-up loop. We demonstrate the fabrication of these chip-based traps and particles using thin Nb films. Our modeling is generic and has applications beyond the one considered, such as for designing superconducting magnetic shields or for calculating filling factors in superconducting resonators.


I. INTRODUCTION
Superconducting magnetic levitation 1,2 is a fascinating phenomenon and its applications range from demonstration experiments 3 to precise measurements of gravity using the superconducting gravimeter 4 . Recently, theoretical proposals suggest the use of superconducting magnetic levitation as a means to enable new experiments in the field of quantum optics 5,6 . Specifically, micrometersized superconducting or magnetic objects levitated by magnetic fields are proposed to lead to a new generation of quantum experiments that enable spatial superposition states of levitated objects [5][6][7][8] or ultra-high sensitivities for measurement of forces or accelerations 7,9,10 , with recent experiments along these lines [11][12][13][14][15][16] .
Superconducting magnetic levitation experiments require superconducting traps that generate magnetic field distributions that are capable of levitating superconducting particles. As a superconducting particle is a perfect diamagnet, the traps require a local magnetic field minimum accompanied by a field gradient for achieving stable levitation 17 . For the cases we consider, the levitated particles and the associated traps have dimensions on the order of micrometers. Such micrometer-sized traps can be fabricated using standard techniques developed for superconducting thin films, as we exemplify using thin films of Nb 18 . Superconducting chip-based trap structures have already been developed in the context of atom * witlef.wieczorek@chalmers.se optics for trapping atomic clouds on top of superconducting chips [19][20][21][22] . However, in contrast to trapped atomic clouds, a levitated particle has a finite extent and, thus, requires accounting for its volume and the finite magnetic field penetration in the levitated object such that trap properties can be accurately predicted. Analytical formulas exist for idealized geometries, such as for a diamagnetic sphere in a quadrupole field 5 or in a field of four parallel wires 8 , for a superconducting sphere in a quadrupole field 23 , for a diamagnetic ring in the field of an anti-Helmholtz coil 24 or can be derived for symmetric geometries using the image method 25 . However, in the general case when considering realistic trap geometries of reduced symmetry, trap wires of finite extent or arbitrary shapes of the levitated particle, analytical formulas do not exist and one has to resort to modeling using finite-element methods (FEM).
structures, micrometer-sized multilayer anti-Helmholtz coil-like (AHC) and single layer double-loop traps. We compare the modeling results to more ideal situations of increased symmetry, where analytical results exist 23,30 or results can be obtained using the image method 25 . We note that the presented FEM modeling can be also used for simulating larger trap structures and particles, which are interesting in the context of acceleration sensing 9,13 . We consider different particle shapes and find that higher trap frequencies can be achieved for non-spherical particles that locally increase the magnetic field gradient. Our models are generic and can be applied to other settings involving superconductors in the Meissner state.
Our work is structured as follows. In Sec. II we summarize the underlying simulation model and in Sec. III state the principles for realizing stable magnetic levitation. We describe the fabrication of the two considered trap architectures and of different particle shapes in Sec. IV. We turn to the simulation of the traps in Sec. V and discuss the results in the light of achieving stable levitation. In Sec. VI we analyze measurement of the centerof-mass motion of a particle levitated in a superconducting trap and calculate the expected signal detected in a pick-up loop. We conclude our work in Sec. VII. The appendix contains details on the modeling method (Appendix A), the validation of our modeling against test cases (Appendix B), dependence of the results on the FEM mesh (Appendix C) and the relation between the Maxwell stress tensor, forces and fields (Appendix D).

II. MODELING
The FEM simulations we use are based on a simplification of the Ginzburg-Landau equations for the superconducting state in the limiting case of small external magnetic field B. Ginzburg-Landau theory predicts for a superconductor under a magnetic field that its supercurrent J s is given by 31 : where e and m are the charge and mass of the electron, respectively, is the reduced Planck constant, A is the magnetic vector potential, Ψ is the wave-function that describes the cooper pair density in the superconductor, and |Ψ| 2 is the cooper pair density. The cooper pair density gets smaller once an external magnetic field B is applied and equals zero when B is larger than the critical field B c . Reversely, in the absence of a magnetic field, |Ψ| 2 is maximum and constant. For small magnetic fields |Ψ| 2 decreases only slightly, making ∇Ψ and ∇Ψ * negligible. Thus, the equation for the supercurrent reduces to the London model 31 : where λ L = m µ 0 |Ψ | 2 e 2 is the London penetration depth and µ 0 the magnetic permeability of vacuum. By implementing this equation in FEM software as an external contribution to the current density in the superconductor domains, one can model domains as type I superconductors in the Meissner state. The model solves the Maxwell-London equations using A-V formulation [26][27][28][29] . That is, the field equations are solved using the magnetic vector potential A and the voltage V as the dependent variables, for details see Appendix A. In our case, the field equations are solved in the quasi-static regime, so time derivatives of the equations describing the system are not involved. We would like to point out that describing dynamic systems is possible as shown in Ref. 32 . We note that, if B is larger than the first critical field, B c1 , magnetic flux vortices will start nucleating in the superconductor and the gradient terms in Eq. (1) cannot be neglected. Thus B c1 puts a bound on the maximal trap strength that is achievable in our modeling, along with the critical current density J c in the coil material. Future work could extend our modeling by, e.g., using the critical state model 33 applicable to hard type-II superconductors.
Another feature of superconductivity is fluxoid quantization, which should be accounted for to accurately describe superconducting objects with holes. In our case, this concerns the levitation of ring-like particles. Fluxoid quantization can be derived by integrating Eq. (1) over a closed loop in the superconductor, which contains a hole with magnetic flux Φ hole . This results in 31 : where θ is the phase of the wave function, n an integer, Φ 0 = h/2e is the magnetic flux quantum. Eq. (3) tells us that the supercurrent will preserve the magnetic flux threading the hole of the superconductor as the multiple of Φ 0 closest to Φ hole .
Since our model does not account for the contribution of the wave function's gradient, fluxoid quantization cannot emerge from the implementation of Eq. (2). We simplify our modeling by considering only flux quantization in the following and, thus, neglect the flux in the ring's interior material caused by the finite penetration depth of the external magnetic field. This approximation is reasonable for Λ/R 1, where Λ = λ 2 L /d is the two dimensional effective penetration depth, R is the size of the superconducting object and d the superconducting film thickness 34 . In our case we have Λ/R < 0.05, thus, fulfilling the approximation. We implement flux quantization ad hoc by defining the area of the hole in the superconductor over which Eq. (3) is integrated and impose an additional contribution to the current density of the superconductor such that the constraint is fulfilled within the defined area. In this way, a superconductor with trapped flux in a hole can be modeled.
To perform the modeling we use COMSOL Multiphysics 35 .
Equations (1) and (4) are implemented in the magnetic formulation module in the quasi-static regime. We verified our modeling against test cases, where analytical results exist and find perfect agreement, see Appendix B.

III. PRINCIPLES OF MAGNETIC LEVITATION
The goal of the chip-based traps is to stably levitate a superconducting particle in a point r lev in free space above the surface of the chip. To this end, a local energy minimum in the potential energy landscape U(r) of the superconducting particle is required, with U(r) given by 17 : where M is the magnetization and m the mass of the particle, g is the gravitational acceleration and z is the height above the chip surface. The integration goes over the volume of the levitated particle. For illustration, let us assume the superconducting particle to be a point particle with magnetic moment m = V M = −V B/µ 0 . Then, the force acting on the particle is 17 : and we see that levitation is achieved when F(r lev ) = 0, that is, when B · ∇B = −µ 0 gρe z at r lev . In reality, we cannot make the above approximation and we need to evaluate Eq. (5) for realistic geometries, see Appendix D for details how this is implemented in FEM.
While balance of forces is a necessary condition, it is not sufficient. Additionally, the local energy minimum at r = (x, y, z) T = r lev = (x lev , y lev , z lev ) T must fulfill that 17 ∂ 2 U(r)/∂ x 2 > 0, ∂ 2 U(r)/∂ y 2 > 0 and ∂ 2 U(r)/∂z 2 > 0 in order to achieve stable levitation. Then, the particle experiences a restoring force in the trap.
For the traps under consideration, our simulations show that for small particle displacements (on the order of 100's of nanometers) the force acting on the particle can be approximated linearly, with F(r) = (F x , F y , F z ) T = (k x x, k y y, k z z) T and k i being the spring constant of the trap along direction i. Then, the frequency with which the particle oscillates is where (x 1 , x 2 , x 3 ) T ≡ (x, y, z) T . Assuming that thermal energy is driving the particle's motion, the root mean square oscillation amplitude x rms,i in a given direction i for the particle being in thermal equilibrium with the environment at temperature T is: where k B is Boltzmann's constant. Typically, this amplitude is much smaller than the particle or the trap region, e.g., for a realistic situation that we consider we assume a 1 µm particle at a temperature of T = 4 K in a trap with ω z = 2π · 1 kHz and obtain x rms,z = 25 nm.

IV. FABRICATION OF CHIP-BASED SUPERCONDUCTING TRAPS AND PARTICLES
We consider two superconducting trap architectures that can be used to stably levitate superconducting objects of spherical, cylinder or ring shape. We put particular focus on the capability to microfabricate these traps and particles in order to guarantee reproducible parameters and scalable device architectures.

A. Fabrication of anti-Helmholtz coil-type trap
The AHC-type trap is formed by two coils arranged in an anti-Helmholtz-like configuration. Fig. 1(a) shows a schematic and Fig. 1(b,c) scanning electron microscope (SEM) images of a microfabricated trap with 3 µm inner coil radius and 1 µm vertical coil separation. Note, that we deviated from the ideal coil separation of 3 µm due to limitations in growing thick layers of the Si spacing layer. Alternatively, one can pattern a hole of 2 µm diameter to approach an ideal AHC-trap, as modeled in Sec. V A. The layers Nb/Si/Nb, which are 300/1000/300 nm thick, respectively, were sequentially deposited on a silicon substrate by sputtering. Each of the layers was separately structured via electron beam lithography (EBL) or optical lithography and subsequently etched using inductively coupled plasma-reactive ion etching (ICP-RIE) 36 . The lowermost coil from Nb is fabricated first and covered by Si. The Si layer is etched via RIE to expose the contact pads of the lowermost coil. Then, the upper Nb layer is deposited and structured. The contact pads made via etching the Si layer ensure an electrical connection between the lower and upper Nb layers through the material sputtered on the sidewalls of the opening. Finally, a hole is etched from the top of the Si layer down into the Si substrate. The volume contained in the hole is the trap region. The yellowish region marks a 300 nm thick and 1 µm wide Nb wire on a Si substrate.

B. Fabrication of double-loop trap
An alternative trap arrangement consists of two concentric and co-planar coils with counter-propagating currents that generate a local energy minimum above the plane of the coils such that a particle will be stably levitated in this minimum as we show in Sec. V B. A schematic of such a planar double-loop trap is shown in Fig. 1(d), which can be seen as a planar version of an AHC-trap. Fig. 1(e) shows a microfabricated double-loop trap. We use a 300 nm thick Nb film to structure the wire. The double-loop trap has the advantage of a simple and straight-forward microfabrication process, since it can be fabricated with a single layer layout. Furthermore, the trap region is not restricted by the vertical separation between coils, unlike as is the case for the AHC-trap.

C. Superconducting properties of Nb films
We determined the superconducting properties of the 300 nm thick Nb film to have a T c ≈ 9 K, a critical current density up to j c = 5 · 10 11 A/m 2 and a critical field B c2 ≈ 0.4 T. For the simulations in our work, we assume a current in the trap wires that corresponds to a current density of 1 · 10 11 A/m 2 , which is below the critical current measured and, thus, leaves a safe margin for realistically running such a current through the trap. We model these wires as superconductors of finite extent, unless otherwise stated, with a thickness of 300 nm.

D. Fabrication of particles
The particles can be obtained from metallic powders or can be microfabricated directly in the traps. Fig. 2(a) shows a spherical Pb particle individually selected from Pb powder. Note, however, that most particles in the powder are non-spherical and one has to pick-and-place the desired particles into the trap region. A systematic approach towards fabricating particles of well defined shapes can rely on etching of thin superconducting layers. Fig. 2(b) shows such an approach, where we fabricated cylinder-and ring-shaped particles from a Nb layer. The particles are fabricated by sputtering a 300 nm thick Nb layer on top of a sacrificial layer of hard-baked resist. The particle shape is patterned via EBL followed by ICP-RIE etching in the same way as the trap architectures. The sacrificial resist layer is removed using oxygen plasma, releasing the particles onto the substrate.
Note, that all the following simulations do not exceed the first critical field of Nb (B c1 = 0.17 T) or of Pb (B c = 0.08 T) on the particle surface in order to conform with our modeling assumption of levitating superconductors in the Meissner state. Levitation of superconducting particles in the Meissner state is key for macroscopic quantum experiments, see Refs. 5,8 , to avoid undesired dissipation via flux penetration.

V. MODELING OF CHIP-BASED SUPERCONDUCTING TRAPS
In the following, we model the properties of the superconducting traps using different methods: Analytic: The simplest model assumes the coils to be one dimensional current loops and the superconducting object to be a point particle with an induced dipole moment, hence there is no dependence on the particle geometry. Analytical results can be derived for this case 30 .
Image: We consider the image method 25 that assumes the particle to be a perfect diamagnetic sphere in a magnetic field distribution that is rotationally symmetric. The image method allows obtaining results without resorting to FEM simulations.
Hofer: We use the results derived by Hofer et al. Ref. 23 that are applicable for trapping a superconducting sphere in a quadrupole field. This situation accounts for finite field penetration, but cannot account for the deviation of the magnetic field generated by an AHC from the idealized quadrupole field.
We use different levels of realistic FEM modeling: FEM-2D with 1D-wires: We account for the volume of the superconducting particle and field penetration. However, the magnetic field is generated by quasi one dimensional (1D) coils (modeled with a cross section of 1 × 1 nm 2 ). The coils are closed loops, resulting in a rotationally symmetric geometry. The corresponding magnetic field distributions are simulated with FEM using the modeling described in Sec. II.

FEM-2D:
A more realistic implementation accounts for the volumes of the particle and the coils while assuming that the coils are closed loops, resulting in a rotationally symmetric geometry. The corresponding magnetic field distributions are simulated with FEM, see Sec. II.
FEM-3D: The most realistic approach accounts for the symmetry breaking of the coil openings, which must be there to feed current into the coils. We use again the modeling as described in Sec. II.
We note that the most realistic modeling is FEM-3D, which is also the most computationally expensive in terms of runtime on a computer. The analytic results obtained via Ref. 30 , the image method 25 and Ref. 23 allow comparing the results obtained with FEM simulations to the trend observed when approximating the trap geometry or the finite penetration depth of the magnetic field in the superconducting regions. Note also that we keep for comparison the vertical distance between the coils as well as the inner diameter of the coil for all methods the same.

A. Anti-Helmholtz coil-trap
We first analyze the magnetic field distribution of an AHC trap obtained via the analytic formula for calculating the field of an empty AHC [ Fig . From the field distributions one observes that while there exist general similarities between the predictions of each model, the differences between them are also significant. Most importantly, the magnitude of the magnetic field and its distribution change, which affect the trap frequency and levitation point. In the following, we discuss the stability of the levitated particle in the trap and the resulting trap frequency.
a. Trap stability for translational degrees of freedom A stably levitated particle experiences a restoring force for all translational degrees of freedom. Fig. 4 shows the total force acting on the superconducting sphere close to the center of the AHC-trap. The force is calculated as explained in Appendix D from the magnetic field distributions shown in Fig. 3(f,i). At the center of the trap, the force equals zero as the magnetic levitation force and the gravitational force balance each other. The negative gradient of the force corresponds to a restoring force pushing the particle back to the trap center for small displacements. Thus, this parameter set results in a stably levitated particle. The thick solid lines are linear fits within ±100 nm of the trap region from which the restoring constants k i and, thus, trap frequencies ω i are calculated, see Eq. (7). In the following, we use the error of the linear fit on k for calculating the uncertainty on the trap frequency.
b. Trap stability for angular degrees of freedom For any non-spherical particle, the magnetic field distribution can result in a torque, which will align the particle in the field until its potential energy is minimized. Hence, just as there is a trap center for the translational degrees of freedom, there is a stable particle orientation for the angular degrees of freedom. The torque τ acting on the particle is calculated as detailed in Appendix D. The torque restores the particle to the stable orientation with a spring constant k τ i given by:  Figure 4. Force acting on a 1 µm diameter spherical particle in the AHC-trap considered in Fig. 3(c,f,i) as a function of the particle's displacement relative to the center of the AHCtrap in the x, y, z directions. Open circles show the FEM-3D simulation results, the thin solid line is a guide to the eye and the thick solid line represents a linear fit, from which the restoring constant k i is calculated.
where τ i is the torque component w.r.t. rotation around the x/y axis and θ j is the particle tilt with re-spect to the y/x axis.
In the following we consider the torque acting on a cylinder or a ring levitated in the realistic AHC-trap. Fig. 5 shows the x and y components of the torque versus the inclination of the particle with respect to the other axis, i.e., the x component of the torque tilts the particle with respect to the y axis and vice-versa. Equilibrium orientations are found for orientations in which the torque is zero and its slope is negative. From all orientations fulfilling these conditions, the one with the largest (absolute) slope is the stable orientation, while the others will be metastable. The stable and a metastable orientation for a cylinder are 0 and π/2 tilt with respect to the y axis, respectively. For the orientation with respect to the x axis, the stable orientation is close to 0, with a slight shift in angle due to the coil openings. For a ring with zero trapped flux the stable and metastable orientations are at 0 and π/2 tilt with respect to the y-axis, respectively. However, for the orientation with respect to the x-axis, there is only one stable orientation close to π/6. This asymmetry is caused by the coil openings and flux quantization in the ring hole, which generates an extra current that increases the energy of the system.   Sphere ω x,y,z /2π (kHz) Cylinder ω x,y,z /2π (kHz) Ring ω x,y,z /2π (kHz)  Table I. Trap frequencies in an AHC trap for levitating particles of different shape: a sphere of 1 µm diameter, a cylinder of 1 µm diameter and 300 nm height, and a ring with 300 nm thickness and inner and outer diameters of 0.5 µm and 1 µm, respectively. Other parameters of the trap as in Fig. 3. Note that trap frequencies along the x and y directions for the case FEM-2D were actually simulated using FEM-3D and a symmetric trap. The uncertainties for the frequencies are below 0.7% for the sphere, 1.3% for the cylinder and 0.5% for the ring.
Hence, a torque that works to minimize the current in the ring will appear, orienting it towards the coil openings, where the field is weaker. If the AHC-trap had no openings, such as the one shown in Fig. 3(b), the stable and a metastable orientation would appear at 0 and π/2 with respect to the y-axis, respectively. Note that the torque acting on the ring is approximately one order of magnitude stronger than that on the cylinder.
c. Trap frequency Tab. I shows the obtained trap frequencies of the three particle shapes, sphere, cylinder and ring, all with the same diameter, for a specific set of parameters, cf. with Ref. 9 . There is some general behavior observable. First, when accounting for the volume of the particle and treating it as a superconductor in the Meissner state instead of a perfect diamagnet (i.e. when comparing the analytic or image method to the methods of Hofer, FEM-2D and FEM-3D), the magnetic field gradient around the particle is smoothed, thus, decreasing the trap frequency. Second, accounting for the volume of the coils and their openings spreads out the magnetic field, decreasing its gradient at the trap position, thus, leading to different oscillation frequencies in the x and y directions (FEM-2D vs. FEM-3D).
Particles of non-spherical shape result in higher trap frequencies in the z axis, which can be traced back to the increased magnetic field gradient at the particle surface caused by the larger demagnetizing effect when compared to a spherical particle, see Appendix B 2.
Different trap and particle sizes have a systematic effect on the trap frequency. We first consider a scaling of the AHC-trap architecture by scaling the trap geometry with a scaling factor such that the dimensions of sphere and trap are simultaneously changing size, while keeping the current density in the coils and λ L constant. Fig. 6(a) shows the trap frequency in the z direction obtained using different methods. For large geometries, i.e., negligibly small penetration depths, the image method is in agreement with FEM-2D with 1D-wires, since the normal conducting volume fraction of the particle is negligible compared to its superconducting volume fraction and the trap wires are quasi one dimensional. The decrease of the trap frequency for FEM-2D with 1D-wires when scaling down the system (i.e. for scaling factors smaller than 2, i.e., 1/scaling factor larger than 0.5) is due to the fact that for particles with a radius approaching λ L , a significant portion of the sphere's volume becomes a normal conductor. This volume does not contribute strongly to the magnetic energy due to the much lower magnetic susceptibility. Thus, the magnetic force on the particle weakens, which decreases ω i leading to the difference between the image method and FEM-2D with 1D-wires. When modeling the coils of finite extent via FEM-2D, the gradient of the field decreases, thus decreasing the trap frequency. The FEM-3D modeling also incorporates the coil opening, thus further reducing the trap frequency. Assuming a superconducting particle in a quadrupole field from Hofer et al. 23 results in the same trend of decreasing trap frequency when scaling down the trap geometry. Interestingly, the prediction of trap frequencies from Hofer comes close to FEM-3D modeling.
We now consider the dependence of the trap frequency on the particle radius in a trap with unaltered size. In Fig. 6(b) we observe that for large particles the image method gives similar results as FEM-2D with 1D-wires. Deviations occur when the particle radius is decreased to a size that magnetic field penetration becomes relevant, i.e., for λ L /R sphere 0.1. As before, when modeling the finite extent of the wires via FEM-2D the trap frequency decreases compared to 1D modeling of the wires. When accounting for the opening of the trap wires via FEM-3D, the trap frequency further decreases. Again, the modeling from Hofer results in a similar trend as FEM-3D. Levitation of a ring is particularly interesting. Fig. 7 shows the dependence of trap frequency for the amount of trapped flux, Φ t , in the ring for an AHC trap with 100 mA of current (corresponding to increasing j c to 3.3 · 10 11 A/m 2 ). Strikingly, the trap frequency decreases with increasing number of trapped Φ t , regardless of its orientation.
The levitation height, however, increases monotonously with Φ t . This is because the ring seeks the region in the trap with a magnetic field strength that will generate the same flux as Φ t . As a result, the ring gets closer to one coil or the other depending on the orientation of Φ t , and, thus, further away from the trap center, where the field gradient is highest, reducing the trap frequency. Note, that if the 10μm 10μm 10μm trap does not generate a high enough magnetic field, a ring with non-zero number of trapped flux quanta will be ejected from the trap.

B. Double-loop trap
We now turn to modeling the properties of the doubleloop trap using the different modeling possibilities.
a. Trap stability The simple layout of the doubleloop trap comes at the expense of sacrificing magnetic field gradient and intensity. To see this, we show the magnetic field distributions for such a trap in Fig. 8. In Fig. 8(d,g) the trap region is visible as the region of the surface plot surrounded by high field intensity. As can be seen in Fig. 8(e,f,h,i) a particle with a diameter similar to the trap size will fill up the trap region and is stable in the z direction due to gravity, since there is no magnetic field from above pushing it down. For these particle sizes, the double-loop trap is magneto-gravitational.
Just as for the AHC-trap, the breaking of symmetry due to the openings of the coils has a significant effect in the double-loop trap. As shown in Fig. 8(i), the field at the opening of the coil on the side of the current feed lines interferes constructively with the field generated by the inner coil, creating a higher field intensity at the left side of the particle that pushes it towards the direction of positive x. At the same time, the field opening at the opposite side weakens the field, creating a lower field intensity at the right side of the particle, which weakens the push in the direction of negative x towards the coil center. These two factors combined weaken the restoring force along the x-axis, which can potentially lead to the particle not being trapped. Thus, a careful design of the double-loop trap is required in order to achieve stable levitation. As a rule of thumb, the opening left between the wires should be smaller than the wire width of the coil. The trap center position is displaced from the coil center towards the opening, and the trap frequency along that direction is smaller than for any other in-plane direction.
b. Trap frequency Tab. II shows trap frequencies for a 10 µm spherical particle in a double-loop trap. The frequencies are below 1 kHz and, compared to the AHC trap rather low, due to this double-loop trap being magnetogravitational. Note, the trap frequency will not change considerably for particles of different shapes, since any increase of the field gradient around the particle will push it higher up into regions of smaller magnetic field and, In the following, we look at systematic effects on the trap frequency. We first consider a scaled system in Fig. 9(a), where both the trap and the particle change size with a scaling factor while keeping the current density of the trap and λ L constant. The geometry corresponding to a scaling factor of one is that of Fig. 8. There is agreement between the image method and the FEM-2D with 1D-wires for all geometries shown. A discrepancy between the methods would be expected for smaller geometries than the ones shown, as was the case for the AHC trap. However, smaller geometries than the ones considered are not capable of levitation because the magnetic lift force does not overcome the weight of the particle. As for the AHC trap, when modeling via FEM-2D the trap frequency decreases due to reducing the field gradient and further decreases when modeling via FEM-3D due to accounting for the wire opening.
In Fig. 9(b) the trap frequencies of a particle with different size in a trap with unaltered dimensions are shown. Once again, the FEM-2D simulations agree with the image method for big particles, but deviate for smaller particles due to the finite field penetration accounted for in FEM-2D with 1D-wires. Modeling via FEM-2D and FEM-3D results in gradually smaller trap frequencies as discussed before. Interestingly, the trap frequency reaches a maximum around λ L /R sphere ∼ 0.05. For larger particles, the trap frequency decreases due to the trap be-coming more magneto-gravitational, whereas for smaller particle sizes the magnetic field penetration into the particle reduces the trap frequency.

VI. FLUX-BASED READ-OUT OF LEVITATED PARTICLE MOTION
A major motivation for undertaking magnetic levitation of superconducting micrometer-sized objects is the exceptional decoupling of the levitated object from its environment 5,6 . To verify this decoupling experimentally, one needs to detect the motion of the levitated particle. Motion detection can rely on flux-based read-out via a pick-up coil placed in the vicinity of the trap 9,14 . Then, particle oscillations around the trap center generate perturbations in the magnetic field distribution, which translate into a change of the magnetic flux threading through a pick-up coil. The pick-up coil could, in turn, be connected to a DC-SQUID, which converts the flux signal into a measurable voltage signal. The expected signal in a pick-up loop has been calculated in previous work for the case of idealized situations 5,9 .
We calculate the expected signal for realistic geometries by accounting for extended volumes and flux quantization based on our FEM modeling. In the following, we consider a 1 µm diameter spherical particle trapped in an AHC-trap (cf. Fig. 3). We are interested in the magnetic flux threading a pick-up coil (shaded area in Fig. 10(a)) for different particle displacements with respect to the trap center. We show in Fig. 10(b) that our model predicts magnetic flux changes similar to the ones of the analytical expressions for a perfect diamagnetic sphere in a quadrupole field 5 .
The slope of the curve in Fig. 10 Table III. Signal strength η i and noise power spectral density S φ i on resonance detected by a pick-up coil with radius of 2 µm located between the two coils of an AHC-trap. The dimensions of the trap and area of the pick-up coil are shown in Fig. 10(a). The trap and particle parameters are the same as in Tab. I. S φ x, y, z (S φ 0x,0y,0z ) denotes the signal assuming a Q of 10 7 and a temperature of 4 K (quantum ground state). The uncertainties are below 25% for the z direction and around 50% for x and y directions.
Commonly, one measures the flux noise power spectral density S φ i (ω), which is given as 37 with S x i (ω) is the noise power spectral density of mechanical motion, x rms,i the root mean square amplitude of the oscillation in direction i (cf. Eq. (8)) and γ i = ω i /Q i is the mechanical damping with Q i the mechanical quality factor. On mechanical resonance, one obtains Tab. III shows η i and S φ i (ω i ) for a sphere, cylinder and ring in an AHC-trap at a temperature of 4 K and for a conservative 5,6 Q = 10 7 . We also consider the case of detecting the ground state motion, i.e., The simulated signals are on the order of 10 3 µφ 0 / √ Hz for thermally driven motion and some µφ 0 / √ Hz for ground state motion. The former signals are well above the noise floor of state-of-the-art SQUID sensors (< 1 µφ 0 / √ Hz for detection frequencies above 1 kHz [38][39][40] ). While detection of ground state motion seems feasible, a further decrease in mechanical damping γ i would be beneficial, as is predicted by theory 5,9 .  Figure 11. Comparison of the signal strength η i and power spectral density S φ i detected by a pick-up loop of varying radius placed between the two coils of an AHC trap. We consider the same trap as the one in Fig. 10. The dashed line indicates the size of the pick-up loop as used in Tab. III. Fig. 11 shows the signal strength and noise power spectral density when varying the pick-up coil radius. For small radii, the FEM results match within their uncertainty the value predicted by the image method, but slightly deviate for larger radii. This is because as the radius of the pick-up loop grows, the FEM model integrates over more coarsely meshed regions of the model and numerical errors accumulate.

VII. CONCLUSIONS
In conclusion, we have used FEM modeling for simulating realistic type I superconducting structures in the Meissner state including flux quantization. We have based our modeling on the A-V formulation for implementing the Maxwell-London equations, which uses only a single parameter, the London penetration depth λ L .
We applied this modeling technique in the context of chip-based superconducting levitation of micrometersized particles to analyze in detail two trap architectures. Crucially, we have shown that trap properties, like trap stability and frequency, differ from idealized models due to breaking of symmetry by coil openings, demagnetizing effects and flux quantization. Our FEM-based simulations incorporate these effects and, thus, allow for a careful design of realistic trap architectures. We predict large trap frequencies in the 40 kHz range for ring-like particles trapped in a micrometer-sized AHC-trap compared to a couple of hundred Hz for a similar sized planar double-loop trap. We have demonstrated how to fabricate such traps and particles from thin Nb films.
Large trap frequencies are beneficial for reducing the impact of environmental vibrational noise on particle motion, both from the seismic background 41 as well as from the vibrational noise coming from the utilized cryostats necessary for cooling superconductors 42,43 . A lower vibrational background noise results in a smaller vibrationinduced heating rate of particle motion, which is beneficial for follow-up experiments targeting motional quantum ground state cooling of levitated particles 5,6 . Further, we have shown that for realistic trap geometries the motion of the particle in the trap can be measured by state-of-the-art SQUIDs, whose measurement signal would be used as feedback signal for performing feedback cooling of particle motion.
Including flux pinning in our modeling [44][45][46] via, for example, the critical state model 33 would allow studying alternative trap opportunities, which may offer chip-based traps with even higher trap frequencies. Our modeling is generic and applicable not only in the present context of simulating chip-based superconducting traps for quantum or sensing experiments, but also for designing superconducting magnetic shields 47 or filling factors in superconducting resonators 32 , when considering superconductors in the Meissner state. 10μm 150nm 50nm Figure 12. Surface mesh of a superconducting sphere and a double-loop trap. The maximum element size has been increased three-fold (to 150 nm in the figure) for the sake of visibility. Parameters as in Fig. 8. The insets show a magnified image of the surface mesh of the sphere with the same exaggerated surface mesh element size as the full image (150 nm) and the one used in the simulations (50 nm).
With that in mind, there are a couple of important things to consider when using FEM modeling for magnetic field distributions around superconducting structures. The mesh needs to be constructed very carefully in the regions close to the superconducting surfaces and regions with high magnetic field density and magnetic field gradient, see as an example Fig. 12. This has to do with the characteristic length scales of superconductors and field gradients in such systems, rather than the method we use to compute them. The magnetic field, current density and thus the magnetic forces on the superconductor are mainly determined by the regions of high magnetic field gradient, which in turn are closely related to demagnetizing effects and the field expulsion of the superconductor. Thus, to obtain valid results the meshing around the surfaces of the superconducting objects needs to be very dense, on the order of λ L , which can become computationally demanding. Furthermore, we need to consider that the magnetic field is gauge invariant. If the initial values for the magnetic vector potential are numerically too big or too small, they can become too hard to compute numerically and the system will not converge. Thus, especially when high values of B are involved, a guess within a few orders of magnitude for A is required for the software to converge. Alternatively, gauge fixing for A can be used if the software allows such possibility.

Appendix B: Validation of FEM modeling
In order to validate our model with test cases, we compare the results it yields to the predictions made by analytical equations. The cases used for validation are (i) the magnetic field expulsion of a superconductor and (ii) demagnetizing effects of superconducting objects with dif-ferent geometries. We also look at flux quantization in a ring and calculate the torque acting on a ring in a homogeneous magnetic field.

Magnetic Field Expulsion
To examine magnetic field expulsion we consider (i) a flat superconducting object with infinite extension in the z and positive x axes and (ii) a superconducting thin film with infinite extension in the z axis, under a homogeneous magnetic field B 0 = B 0 · e z , see Fig. 13. For the first case, the Maxwell-London equations predict that B 0 is expected to decay exponentially within the superconductor with the characteristic length scale λ L (for superconductors with sizes λ L ) 31 : where x is the distance from the superconductor's surface. For the second case, the magnetic field inside a superconducting thin film of thickness t is expected to also decay exponentially from both sides, but the tails of each exponential will overlap in the middle of the thin film, thus limiting the field expulsion 31 : (B2) We simulate the structures for case (i) with a semiinfinite superconductor that occupies the positive half space x > 0 and all z, and for case (ii) with a superconducting thin film with t = 1 µm in x direction centered at zero while y = z = ∞. In both cases, we use λ L = 100 nm and a homogeneous magnetic field B 0 = B 0 · e z with B 0 = 100 mT applied parallel to the z axis. The results are shown in Fig. 13(c,d) and show excellent agreement between FEM modeling and analytical equations.

Demagnetizing Effects
Field expulsion concentrates field lines around the surfaces of the superconducting object parallel to the field. In these regions, an increase of the magnetic field intensity appears. This increase can be calculated analytically as a multiplying factor called demagnetizing factor.
Demagnetizing effects arise naturally in our modeling. In Fig. 14 we show the magnetic field distribution around a micrometer-sized sphere, cylinder and ring, under a homogeneous magnetic field B 0 = B 0 · e z with B 0 = 30 mT. The demagnetizing factors for a perfect diamagnet with such geometries are 1.5, 1.8 and 1.8, respectively 48 . Our modeling as shown in Fig. 14 perfectly matches the analytically calculated values when λ L is close to zero, i.e., for an ideal diamagnet. In the case of the ring, flux quantization is partly responsible for the magnetic field distribution within the ring. As indicated in Fig. 14(f), the light section of the curve represents negative values of B, which are generated by the supercurrent in the ring to keep Φ hole = 0.
In case of chip-based superconducting traps, deviation from the spherical particle shape and flux quantization result in larger magnetic field gradients, which lead to higher trap frequencies, as we described in Sec. V.

Flux quantization
In general, generating a supercurrent has an energy cost. Then, it follows that the energy of the superconductor is minimized when the amount of supercurrent in it is smallest.
Such an effect is shown in Fig. 15, where we calculate the x component of the torque acting on a superconducting ring in a homogeneous magnetic field B 0 as a function of the ring's inclination with respect to the y axis. We consider the cases for a superconducting ring with (i) no flux quantization, (ii) flux quantization with zero flux trapped and (iii) one flux quantum trapped with the same orientation as B 0 .
The ring with no flux quantization experiences a torque because the field is less perturbed when B 0 is parallel to the area of the hole than when it is perpendicular. Hence, it takes less supercurrent to expel the field when θ = π/2 or 3π/2. When the area of the hole is perpendicular to B 0 the torque on the ring vanishes due to symmetry, since it is as likely to tilt clockwise or counter-clockwise, in other words, it is in an unstable equilibrium.
The stable configuration for the ring including flux quantization and no trapped flux, i.e., Φ t = 0, is to be oriented so that no flux is threading the hole, i.e., π/2 or 3π/2. The difference is that the torque is stronger due to additional current from flux quantization that keeps Φ t = 0 when θ π/2 or 3π/2.
For the case of a ring with one trapped flux quantum parallel to B 0 , the configuration in which the least supercurrent is generated is that where B 0 is parallel to the trapped flux quantum, since B 0 is chosen so that the flux through the hole equals Φ 0 when the ring is perpendicular to the field. Thus, the ring will experience a torque that will force it to θ = 0. For θ = π the ring will be unstable because the flux through the hole at this configuration is maximum (Φ = 2Φ 0 ).
Additionally, the inductance of the ring with flux quantization was calculated with the FEM model and found to be 2 pH, which is in reasonable agreement with the 1.6 pH predicted by Ref. 34 .

Appendix C: The effect of meshing
Given that the model is based on FEM, the results are bound to be dependent on the meshing of the system. As mentioned in Appendix A, constructing a mesh fine enough at the surface of the superconducting domains is critical to get reliable results. This dependency is best illustrated in Fig. 16, where the trap frequency along z for a 1 µm diameter sphere in an AHC trap (cf. Fig. 3) is calculated. For these simulations we changed the maximally allowed mesh element size, l mesh , on the surface of the particle. For maximal mesh element sizes l mesh 5 · λ L = 250 nm, we observe no clear trend of the trap frequency within its uncertainty. In this regime, the meshing is fine enough to resolve the spatially varying field penetration in the particle, on the order of λ L . However, for l mesh 0.5 · R sphere = 250 nm, the particle itself is not properly resolved and the magnetic field penetrates parts or the entire volume of the particle, which effectively increases the effect of field penetration and, thus decreases the trap frequency.
Furthermore, the trap frequency dependence on the mesh might not only be related to the mesh element size itself, but also on differences in the mesh being differently built by the FEM model for the very same calculation. To test this dependency the same trap as before was simulated for slightly different maximal surface mesh element sizes of (49.9, 49.95, 50.00, 50.05, 50.1) nm. We get trap frequencies of (23.7, 23.6, 23.8, 23.6, 23.6) kHz, respec-  Figure 16. Trap frequency in the z axis of a spherical particle in an AHC trap as considered in Fig. 3 as a function of the maximal mesh element size l mesh . The insets show the mesh on the surface of the particle at the given element sizes.
tively, resulting in a mean value of (23.66 ± 0.09) kHz. Thus, the scatter of trap frequency from using nearly similar meshes of about ±0.38% is smaller than the fit uncertainty for calculating trap frequencies ω or signal strengths η.

Appendix D: Magnetic force calculation
In our model, the electromagnetic force and the torque on an object are calculated via the Maxwell stress tensor T, whose components T i j are given as: where 0 and µ 0 are the electrical permittivity and magnetic permeability, respectively, and E i and B i are the vector components of the electric and the magnetic field. The knowledge of the field distributions E(r) and B(r) is sufficient to calculate electromagnetic forces and torques via surface integrals as 49 and where F is the force, τ is the torque, Ω is the surface of the particle and r and r 0 are the application point of the force and the center of mass of the particle, respectively.