Chip-based superconducting traps for levitation of micrometer-sized particles in the Meissner state

We present a detailed analysis of two chip-based superconducting trap architectures capable of levitating micrometer-sized superconducting particles in the Meissner state. These architectures are suitable for performing novel quantum experiments with more massive particles or for force and acceleration sensors of unprecedented sensitivity. We focus in our work on a chip-based anti-Helmholtz coil-type trap (AHC) and a planar double-loop (DL) trap. We demonstrate their fabrication from superconducting Nb films and the fabrication of superconducting particles from Nb or Pb. We apply finite element modeling (FEM) to analyze these two trap architectures in detail with respect to trap stability and frequency. Crucially, in FEM we account for the complete three-dimensional geometry of the traps, finite magnetic field penetration into the levitated superconducting particle, demagnetizing effects, and flux quantization. We can, thus, analyze trap properties beyond assumptions made in analytical models. We find that realistic AHC traps yield trap frequencies well above 10 kHz for levitation of micrometer-sized particles and can be fabricated with a three-layer process, while DL traps enable trap frequencies below 1 kHz and are simpler to fabricate in a single-layer process. Our numerical results guide future experiments aiming at levitating micrometer-sized particles in the Meissner state with chip-based superconducting traps. The modeling we use is also applicable in other scenarios using superconductors in the Meissner state, such as for designing superconducting magnetic shields or for calculating filling factors in superconducting resonators.


Introduction
Superconducting magnetic levitation [1,2] is a fascinating phenomenon. Its applications range from demonstration Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 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, micrometer-sized superconducting or magnetic particles levitated by magnetic fields are proposed to lead to a new generation of quantum experiments that enable spatial superposition states of levitated particles [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].
We consider levitation of superconducting particles in the Meissner state, inspired by references [5,6,8]. Their stable levitation requires traps that generate a local magnetic field minimum accompanied by a field gradient [17]. Superconducting chip-based trap structures have already been developed in the context of atom optics for trapping atomic clouds on top of superconducting chips [18][19][20][21]. 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 levitation of a perfect diamagnetic sphere in a quadrupole field [5] or in a field of four parallel wires [8], for a superconducting sphere in a quadrupole field [22], for a perfect diamagnetic ring in a quadrupole field [23] or can be derived for symmetric geometries and perfect diamagnetic objects using the image method [24]. However, in the general case when considering realistic three-dimensional trap geometries with 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).
In our work, we present the fabrication and modeling of two promising chip-based trap architectures suitable for levitation of micrometer-sized superconducting objects of spherical, cylindrical or ring shape. We focus on multi-layer anti-Helmholtz coil-like traps (AHC) and single-layer double-loop traps (DL). We first demonstrate fabrication of the traps using thin films of Nb [25] and of particles made from Nb or Pb of spherical, cylindrical or ring shape. We then use FEM-based simulations to numerically calculate crucial trap parameters, such as stability, frequency and levitation height, for realistic geometries incorporating the finite extent of the wires and the non-symmetry of the traps.
Our FEM simulations are based on implementing Maxwell-London equations in the static regime using the A-V formulation under the assumption that the levitated particles are in the Meissner state [26][27][28][29]. We specifically assume levitation of a particle in the Meissner state, which has been proposed to minimize mechanical loss [5,8], a limiting factor for performing quantum experiments. We compare the numerical FEM results to idealized situations of increased symmetry, where analytical results can be obtained [22,24,30]. While the analytical results are indicative of the underlying physics, numerical modeling yields predictions independent of most idealizing assumptions. Finally, we apply FEM modeling to estimate the signal induced by the motion of a levitated particle in a nearby pick-up loop. This signal would be used to manipulate the center-of-mass motion of the particle in subsequent quantum experiments [5].

Microfabrication of traps and particles
In the following, we describe the microfabrication of chip-based traps from superconducting Nb films and of  superconducting particles from Pb and Nb. Note that other superconducting materials, such as Al, can also be used. The choice of material determines the maximal allowed temperature of the cryogenic environment. While Pb and Nb, for example, allow levitation at liquid He temperatures, Al requires temperatures below 1.2 K. Further, the particles need to be in the Meissner state to avoid mechanical loss [5,8]. Hence, the magnetic field close to the particle surface must be smaller than the first critical field of the chosen material.

Fabrication of traps
The AHC-type trap is formed by two coils arranged in an anti-Helmholtz-like configuration. This trap yields a large magnetic field gradient in the trap center, resulting in trap frequencies above 10 kHz, see section 3.2. Figure 1(a) shows a schematic and figure 1(b) scanning electron microscope (SEM) images of a trap with 3 µm inner coil radius and 1 µm vertical coil separation fabricated in a three-layer process. The three layers used are Nb/Si/Nb, which are 300/1000/300 nm thick, respectively. The lower Nb layer is sputtered first and subsequently patterned by optical lithography and etched using inductively coupled plasma-reactive ion etching (ICP-RIE) [31]. Then, the Si layer is sputtered and subsequently etched via RIE to expose the contact pads of the lower Nb coil. The upper Nb layer is sputtered on top of this Si layer and structured. An electrical connection between the lower and upper Nb layer is facilitated by the Nb material sputtered on the sidewalls of the openings in the Si layer. Finally, a hole is etched through the three layers via ICP-RIE, which becomes the trapping region.
An alternative trap arrangement consists of two concentric and co-planar coils that carry counter-propagating currents. A schematic of such a DL trap is shown in figure 1(c), which can be regarded as a AHC-type trap in the plane. Figure 1(d) shows a microfabricated DL trap made from a 300 nm thick Nb film and patterned via electron beam lithography (EBL). This trap generates a local energy minimum above the plane of the coils, where a particle will be stably levitated with trap frequencies below 1 kHz, see section 3.3. The DL trap has the advantages of a simple single-layer microfabrication process and that the trap region is not restricted by a vertical separation between coils like in the AHC-trap.
We determined the properties of the 300 nm thick Nb film from R-T, I-V and Hall effect measurements 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, similar to previously reported values [32][33][34]. For the analysis of the traps, we will assume a current density in the coil wires of 1 · 10 11 A/m 2 (unless otherwise stated), which is close to the measured critical current density.

Fabrication of particles
The particles can be obtained from particle powders or can be microfabricated directly in the trap. Figure 2(a) shows a spherical Pb particle individually selected from Pb powder. Note, however, that most particles in the powder are nonspherical and one has to pick-and-place the desired particles into the trap region. A systematic approach towards fabricating particles can rely on etching of thin superconducting layers. To this end, we fabricated cylinder-and ring-shaped particles directly on the trap chip by sputtering a 300 nm thick Nb layer on top of a sacrificial layer of hard-baked resist, see figure 2(b). The particle shape is patterned via EBL followed by ICP-RIE etching. The sacrificial resist layer is removed using oxygen plasma, releasing the particles onto the chip.

Numerical analysis of superconducting trap architectures
In the following, we systematically analyze the presented trap architectures with respect to the stability of the trap and achievable trap frequencies for different trap sizes and geometries of the levitated particle. Before we proceed with this analysis, we recall the conditions for achieving stable levitation and present the different models we are going to use.

Models and assumptions
Two requirements have to be met to achieve stable levitation [2,17], see the more detailed discussion in appendix A. First, the magnetic and gravitational force have to balance each other, such that the particle is levitated in free space above the chip surface. Second, the levitation position needs to be stable, i.e. the particle needs to experience a restoring force along each spatial direction. If these two conditions are met, we can calculate a trap frequency, ω t , from the gradient of the force, F, at the levitation position, x lev , via where m is the mass of the particle and k t is the spring constant of the trap. A non-spherical particle also requires rotational stability and, thus, we also analyze torques, τ i , rotating the particle around an axis i by an angle θ j . If stable at θ lev , we calculate a corresponding angular frequency, ω τ i , from where I is the moment of inertia of the particle and k τ i is the angular spring constant. Equation (1) and equation (2) yield accurate trap frequencies as long as the force and torque depend linearly on displacement and angle, respectively, to which we restrict our analysis. Deviations can occur for larger particle amplitudes, see, e.g. references [14,35].
Knowing the magnetic field distribution of a particle in the trap allows calculating the necessary forces and torques, for details see appendix A. Table 1 summarizes the analytical [22,24,30] and FEM models we use for calculating magnetic field distributions of the traps. We consider different levels of FEM modeling, which allow us comparing to the analytical models that necessarily make assumptions about the trap geometry or neglect the finite magnetic field penetration into the particle.
The FEM modeling we use is based on the following assumptions. First, the particle is assumed to be in the Meissner state, which is motivated by the proposals of references [5,8] and implemented in FEM via the A-V formulation of the Maxwell-London equations [26][27][28][29], for details see appendix B, for validation examples see appendix C and for the FEM meshing see appendix D (discretized using quadratic mesh discretization). We, thus, only consider trap fields that remain below the first critical field on the particle surface (we are restricting us to B c = 0.08 T of Pb). Second, we account for flux quantization when considering levitation of a ring ad hoc by defining an area in the FEM model over which the flux should be constant. We neglect the flux in the interior of the material caused by the finite magnetic field penetration depth of the external field. This approximation is valid [36] for Λ/R ≪ 1 (we have Λ/R < 0.04), where Λ = λ 2 L /d is the two-dimensional effective penetration depth, λ L is the London penetration depth, R is the lateral size of the superconducting object and d its thickness. Third, for simplicity we model the wires as very low resistivity, diamagnetic normal Table 1. Different models we apply for calculating the trap architectures. The first three models are analytical models, while the other three are implemented in FEM. Parameters: current through wire I, magnetic field gradient at trap center b, wire radius r, wire thickness t, wire dimensions in 3D [r 3D ], sphere radius R, dimensions of rotational symmetric particle [R 2D ], particle dimensions in 3D [R 3D ], London penetration depth λ L .

Model
Trap Particle Parameters Comment Point particle [30] 1D closed current loops point particle I, r point particle Perfect diamagnet [24] 1D closed current loops superconducting sphere I, r, R, λ L = 0 image method Superconducting sphere [22]   conducting material carrying a uniform current across the wire geometry. The latter assumption is inspired by the situation of using a rectangular type-II superconducting film as wire material transporting a current under self-field that is close to its critical current density [37]. Future extensions could model the wires using the critical state model [38][39][40], which would also allow analysis of various loss mechanisms [41,42]. Note that hysteresis or AC losses are negligible for the cases we are going to consider in section 4 [5]. Finally, we need to consider that the magnetic field and the current density are gauge invariant. The gauge is fixed in the utilized FEM software COMSOL Multiphysics [43] by implementing the Coulomb gauge at the cost of adding an extra variable and by solving the model in the quasi-static regime, see appendix B for details.

Anti-Helmholtz coil-trap
We first analyze the magnetic field distribution of the AHC trap. Figure  for an AHC with a superconducting sphere. As expected, the field distributions depend on the modeling used and, thus, will affect the trap frequency and levitation point.

Trap stability for translational degrees of freedom.
The force acting on the spherical particle can now be calculated from the field distributions. Figure 4 shows the force acting on a superconducting sphere close to the center of the realistic AHC-trap. At the center of the trap, the force equals zero as the magnetic force balances the gravitational force. 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 position from which the spring constants k i , their uncertainties and trap frequencies ω i are calculated.   Table 3. Trap frequencies of a cylinder (1 µm diameter, 300 nm height) and a ring (300 nm thickness, inner and outer diameters of 0.5 µm and 1 µm, respectively) in the AHC trap from figure 3. Note, ωx and ωy for FEM-2D were simulated with FEM-3D and a symmetric trap. The uncertainty on ω for the (cylinder, ring) is below (0.13%, 0.13%) and (1.3%, 0.5%) for FEM-2D and FEM-3D, respectively.

Trap stability for angular degrees of freedom.
When a non-spherical particle, such as a cylinder or ring, is placed in the field of the realistic AHC-trap, torques also act on the particle, see figure 5. Equilibrium orientations are found when the torque is zero and its slope negative, whereby the orientation with the largest slope is the stable and all others are metastable orientations. For a cylinder, a stable and metastable orientation are found at a tilt angle of 0 and π/2 with respect to the y axis, respectively. For the orientation with respect to the The geometric parameters of the trap and particle are taken from table 2 and scaled by a factor while the current density in the coils and λ L are kept constant. The vertical lines indicate the values for the initial geometry. The black points in the insets indicate the location of the 1D-current loops. The grey area represents geometries in which the particle is subject to magnetic fields above 80 mT (Bc of lead) with a maximal field of up to 230 mT. In appendix E we also consider the case when the 1D-current loops are centered in the wire.
x axis, the stable orientation is close to 0, with a slight shift in angle due to the coil openings. For a ring with no trapped flux, a stable and metastable orientation are found at a tilt angle of 0 and π/2 with respect to the y-axis, respectively. However, for the other orientation, there is only one stable orientation close to π/6. This asymmetry is caused by the coil openings and flux quantization that generates an additional current in the ring. A torque acts to minimize this current, orienting the ring towards the coil openings, where the field is weaker. If the AHC-trap had no openings, a stable and metastable orientation would appear at an angle of 0 and π/2 with respect to the y-axis, respectively.

Trap frequency.
The previous analysis confirms that particles of different shapes can be stably levitated in a realistic AHC. We now systematically study the trap frequency and consider first particles of different shape in the same AHC trap, see table 2 and table 3. We observe in table 2 that the trap frequency for a spherical particle along z is by a factor of two larger than along x or y for the analytical models, which is expected due to the ideal anti Helmholtz coil arrangement in the trap. In FEM, however, this factor is reduced, due to the deviation from a quadrupole field caused by the finite extent of the coil wires. We observe further that when accounting for the volume of the particle and treating it as a superconductor in the Meissner state, the magnetic field gradient around the particle is decreased and, thus, also the trap frequency. When also accounting for the opening of the coils via FEM-3D, the magnetic field distribution becomes asymmetric and leads to different trap frequencies along x and y. Table 3 shows that particles of non-spherical shape result in higher trap frequencies along the z axis. This difference can be attributed to the lower mass m of the non-spherical particles as ω t = k t /m (the diameter of all particles is the same). Additionally, the spring constant k t is also different due to the varying demagnetizing effect of each particle shape, for details see appendix C.
We now analyze the dependence of the trap frequency on the size of a spherical particle in a trap with unaltered dimensions. In figure 6(a) we observe that for large particles the perfect diamagnetic sphere model yields similar results as FEM-2D-1D, since the normal conducting volume fraction of the particle is negligible compared to its superconducting volume fraction. Deviations occur when the particle radius is decreased to a size where magnetic field penetration becomes relevant, i.e. for λ L /R sphere ⪆ 0.1. When comparing FEM-2D-1D to a superconducting particle in a quadrupole field [22], we observe that for small particle sizes FEM gives similar results. However, for larger particle sizes (λ L /R sphere ⪅ 0.15), the two methods give different results, which we attribute to the difference between a quadrupole field and the field generated by the wires, becoming more pronounced for larger particles (see also figure E9). When accounting for coils of finite extent via FEM-2D, the gradient of the field decreases compared to FEM-2D-1D and, thus, the trap frequency also decreases. Also in this case, assuming a superconducting sphere in a quadrupole field gives similar results for small particle sizes, but deviates for larger ones. When accounting for the opening of the trap wires via FEM-3D, the trap frequency further decreases, as expected.
In figure 6(b) we analyze a scaled AHC-trap architecture, whereby the dimensions of the particle and trap are simultaneously scaled, while keeping the current density in the coils and λ L constant. For large geometries, i.e. when the penetration depth is small compared to the particle size, the perfect diamagnetic particle method is in agreement with FEM-2D-1D. The decrease of the trap frequency for FEM-2D-1D when scaling down the system (for scaling factors ⪅ 3, i.e. 1/scaling factor ⪆ 0.3) is due to the fact that for particles with a radius approaching λ L a portion of the sphere's volume becomes a normal conductor, and, thus, the magnetic force on the particle weakens. As before, when modeling the finite extent of the wires via FEM-2D the trap frequency decreases compared to FEM-2D-1D. For a superconducting sphere in a quadrupole field, we get similar results for small geometries, but deviations for large geometries. We attribute this behaviour as in figure 6(a) to the deviation of the field of the trap from a quadrupole field.
Levitation of a ring in the AHC trap is particularly interesting. Figure 7 shows that the trap frequency and levitation height depend on the amount of trapped flux, Φ t , in the ring. The trap frequency decreases with increasing number of trapped flux, regardless of its orientation. The levitation height, however, increases monotonously with flux. 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.
To summarize, we find that FEM gives useful predictions for the stability, orientation and trap frequencies of different particle shapes levitated in realistic AHC traps. In contrast, analytical models tend to overestimation of trap frequencies and deviating predictions when scaling the trap geometry, which can be traced back to the assumptions made by these models.

Double-loop trap
We now turn to analyze the properties of the DL trap and show in figure 8 its magnetic field distribution. In figure 8(d),(g) the trap region is visible as the region surrounded by high field intensity. As can be seen in figure 8(e),(f),(h),(i) a particle with a diameter similar to the trap size fills 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 DL trap is magneto-gravitational [11]. Hence, the simple layout of the DL trap comes at the expense of sacrificing magnetic field gradient and intensity.   The breaking of symmetry due to the openings of the coil wires has a significant effect in the DL trap. As shown in figure 8(i), the field 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. This effect can lead to the particle not being trapped. Thus, a careful design of the DL 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. In appendix E we also consider the case when the 1D-current loops are centered in the wire. Table 4 shows trap frequencies for a 10 µm spherical particle in a DL trap. The frequencies are below 1 kHz and, thus, lower compared to the AHC trap due to the DL trap being magneto-gravitational for this particle size. Note, the trap frequency will not change considerably for particles of a different shape, since any increase of the field gradient around the particle will push it higher up into regions of smaller magnetic field and, thus, smaller trap frequency. Figure 9(a) shows the trap frequency in the DL trap when changing particle size. For large particles, FEM-2D-1D agrees with the perfect diamagnetic particle method, while it deviates for smaller particles due to the finite field penetration. Modeling via FEM-2D and FEM-3D results in gradually smaller trap frequencies due to a reduced gradient of the trap. Interestingly, the trap frequency reaches a local maximum around λ L /R sphere ∼ 0.05. For larger particles, the trap frequency decreases due to the trap becoming more magnetogravitational, whereas for smaller particle sizes the magnetic field penetration into the particle leads to a reduction of the trap frequency.

Trap frequency.
In figure 9(b) we consider a scaled system, where both the trap and the particle change size while keeping the current density of the trap and λ L constant. Again, we find agreement between the perfect diamagnetic particle method and the FEM-2D-1D for large geometries and an increasing discrepancy for smaller geometries due to magnetic field penetration. The trap frequency decreases in FEM-2D compared to FEM-2D-1D due to reducing the field gradient and further decreases when modeling via FEM-3D due to accounting for the wire opening.
To summarize the analysis of the DL trap, we find that analytical models overestimate the trap frequency and may even fail to predict stability in case the wire coils have openings.

Numerical analysis of flux-based read-out of particle motion
Magnetic levitation of superconducting micrometer-sized objects promises to reach an exceptional decoupling of the levitated object from its environment [5,6]. To verify this decoupling, 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]. 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 pickup 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 analytically in previous work for the case of idealized situations [5,9].
Using FEM we can now calculate the expected signal for realistic geometries by accounting for extended volumes, field Table 5. Signal strength η i and noise power spectral density S ϕi on mechanical resonance detected by a pick-up coil with 2 µm radius located between the two coils of the AHC-trap. The dimensions of the trap and area of the pick-up coil are shown in figure 10(a). The trap and particle parameters are the same as in table 3. S ϕx,y,z (S ϕ0x,0y,0z ) denotes the signal assuming Q = 10 7 and T = 4 K (quantum ground state). The uncertainties are below 25% for the z direction and around 50% for the x and y directions. penetration and flux quantization. In the following, we first consider a 1 µm diameter spherical particle trapped in an AHC-trap (cf figure 3). We are interested in calculating the magnetic flux threading a pick-up coil for small particle displacements with respect to the trap center, see figure 10(a). In figure 10(b) we compare the analytical prediction for a perfect diamagnetic sphere in a quadrupole field from reference [5] with our numerical FEM-3D results and find similar behaviour.
The slope of the curve in figure 10(b) yields the signal strength per displacement along direction i (normalized by 10 −3 ϕ 0 = 1 mϕ 0 ) as Commonly, one measures the flux noise power spectral density S ϕi (ω), which is given as [44]: with S xi (ω) is the noise power spectral density of mechanical motion, x rms,i = k B T/mω 2 i (k B is Boltzmann's constant, T is temperature) the root mean square amplitude of the oscillation in direction i and γ i = ω i /Q i is the mechanical damping with Q i being the mechanical quality factor. On mechanical resonance, one obtains S ϕi (ω i ) = η i x rms,i / √ γ i . Table 5 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. x 0,i = ℏ/mω i , via measurement of flux, S ϕ0i (ω i ) = η i x 0,i / √ γ i (ℏ is the reduced Planck's constant). The values are on the order of mϕ 0 / √ Hz for thermally driven motion and some µϕ 0 / √ Hz for ground state motion. The former signals are well above the noise floor of stateof-the-art SQUID sensors, which are below 1 µϕ 0 / √ Hz for detection frequencies above 1 kHz [45][46][47]. While detection of ground state motion seems feasible, a further decrease in mechanical damping would be beneficial, as is predicted by theory [5,9]. Figure 10(c) shows the signal strength and noise power spectral density when varying the pick-up coil radius. For small radii, the FEM results correspond within their uncertainty to the values predicted by reference [5], but 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.

Conclusions
We have analyzed in detail using analytical [22,24,30] and FEM modeling two promising trap architectures for levitating micrometer-sized superconducting particles in the Meissner state. The FEM modeling that we used is based on the A-V formulation [26][27][28][29] and is generically applicable for superconductors in the Meissner state, such as for designing superconducting magnetic shields [48] or filling factors in superconducting resonators [49].
Crucially, we have shown that trap properties, like trap stability and frequency, can significantly differ from idealized, analytical models due to breaking of symmetry by coil openings, demagnetizing effects and flux quantization. We found that a chip-based AHC trap is capable of levitating micrometer-sized particles of spherical, cylindrical and ring shape with trap frequencies well above 10 kHz for a current density of 10 11 A/m 2 in the trap wires. However, the fabrication of such a trap on a single chip is complex and requires a three-layer process. A promising alternative would be to use a flip-chip architecture [50]. In contrast, the DL trap is straight forward to fabricate in a single layer process. However, it comes at the expense of considerably lower trap frequencies of below 1 kHz. Further, we confirmed numerically that read-out of the motion of the levitated particle using a pick-up loop in its vicinity [5,9] should lead to clearly detectable signals using presently available SQUID technology [45][46][47]. We, thus, conclude that the analyzed chip-based superconducting traps are a viable approach for future quantum experiments that aim at levitating superconducting particles in the Meissner state [5,6,8].
Extending our modeling by including flux pinning [51][52][53] via, for example, the critical state model [38,39] would allow studying alternative trap opportunities, which may offer chipbased traps with even higher trap frequencies. tion Laboratory at Chalmers. Simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at C3SE, Chalmers, partially funded by the Swedish Research Council through grant agreement no. 2016-07213.

Appendix A. Magnetic levitation, forces and torques
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, B the magnetic field, m the mass of the particle, g 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 perfect diamagnetic point particle with magnetic moment m = V M = −V B/µ 0 . Then, assuming B(r) depends linearly on r, the force acting on the particle is [17]: wherek is the unit vector in the z direction, and we see that levitation is achieved when F(r lev ) = 0, that is, when B · ∇B = −µ 0 gρk at r lev . In reality, we cannot make the above approximation and we need to evaluate equation (5) for an extended volume.
To this end, in our FEM model, the electromagnetic force and the torque on an object are calculated via the Maxwell stress tensor T, whose components T ij are given as: where ε 0 and µ 0 are the electrical permittivity and magnetic permeability, respectively, E i and B i are the vector components of the electric and the magnetic field and δ ij is the Kronecker delta. The knowledge of the field distributions E(r) and B(r) is sufficient to calculate electromagnetic forces and torques via surface integrals as [54] F =˛Ω nTdS, and where τ is the torque, n is the unit vector normal to the particle surface, Ω is the surface of the particle and r and r 0 are the application point of the torque and the center of mass of the particle, respectively. While balance of the gravitational and magnetic force 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 [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, so that the particle experiences a restoring force in the trap.

Appendix B. FEM Modeling
The FEM simulations we use are based on the London model [55] where, for small applied fields, the equation for the supercurrent in a superconductor can be written as [56] where λ L = m µ0|Ψ| 2 e 2 is the London penetration depth, |Ψ| 2 = n c is the squared amplitude of the order parameter's wave function Ψ(r, θ) = |Ψ(r)|e iθr with phase θ, n c is the Cooper pair density, and e is the electron charge. By implementing this equation in FEM software as an external contribution to the current density in the superconductor domains, one can model domains as superconductors in the Meissner state. Note that equation (10) is in general not gauge invariant under the transformation A ′ = A + ∇Φ s , where Φ s is here an arbitrary scalar potential. However, in the specific case we consider, charge is conserved and the potentials A and Φ s change slowly in time (i.e. in the quasi-static regime), such that we can use equation (10) in the Coulomb gauge ∇ · A = 0.
The FEM implementation 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. 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, however, possible as shown in reference [49]. We note that, if B is larger than the first critical field, B c1 , magnetic flux vortices will start nucleating in the superconductor. Thus B c1 puts a bound on the maximal trap strength that can be studied in our modeling.
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 the Ginzburg-Landau equation for the supercurrent [55] (m e is the mass of the electron, ℏ is Planck's constant) over a closed loop in the superconductor, which contains a hole with magnetic flux Φ hole . This results in [55] m e where n is an integer, and Φ 0 = h/2e is the magnetic flux quantum. Equation (12) 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 contributions of the wave function's gradient in equation (11), fluxoid quantization cannot emerge from the implementation of equation (10). We simplify our modeling by considering only flux quantization 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 (we have Λ/R < 0.04), where Λ = λ 2 L /d is the two dimensional effective penetration depth, R is the lateral size of the superconducting object and d its thickness [36]. We implement flux quantization ad hoc by defining the area of the hole in the superconductor over which equation (12) 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.

Appendix C. Validation of FEM modeling
In order to validate our specific FEM implementation, we compare its results to test case, where analytical results exist. To this end, we select the magnetic field expulsion of a superconductor and demagnetizing effects of superconducting objects with different geometries. We also look at flux quantization in a ring and calculate the torque acting on a ring in a homogeneous magnetic field. a. 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 thin superconducting film with infinite extension in the z axis, under a homogeneous magnetic field B 0 = B 0 ·k, see figure  C1. 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 ) [55] , 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  We simulate the structures for case (i) with a semi-infinite 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 ·k with B 0 = 100 mT applied parallel to the z axis. The results are shown in figure C1(c,d) and show excellent agreement between FEM modeling and analytical equations.
b. Demagnetizing Effects. Field expulsion concentrates field lines around the surfaces of the superconducting object parallel to the field. In these regions, an increase of 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 figure C2 we show the magnetic field distribution around a micrometer-sized sphere, cylinder and ring, under a homogeneous magnetic field B 0 = B 0 ·k with B 0 = 30 mT. The demagnetizing factors for a perfect diamagnet with such geometries are 1.5, 1.8 and 1.8, respectively [57]. Our modeling as shown in figure C2 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 figure C2(f), the thin section of the curve represents negative values of the magnetic field, which are generated by the supercurrent in the ring to keep Φ hole = 0.
c. Flux quantization: a ring in a homogeneous magnetic field. 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 figure C3, 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 . Figure D4. Trap frequency of a spherical particle in an AHC trap as considered in figure 3 as a function of the mesh element size l mesh on the particle's surface using FEM-3D. The insets show the mesh on the surface of the particle at the given element sizes. 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 ). d. Flux quantization: a ring and levitation. Reference [23] provides an analytical formula for the trap frequency along the vertical direction for levitating a ring in a quadrupole field, including flux quantization. We compared FEM-2D simulations to this formula for a ring with inner and outer radii of 0.4 µm and 0.5 µm, respectively, thickness of 50 nm and λ L = 50 nm in an AHC trap with coil radius and separation of 10 µm and a current of 3 A. Using FEM-2D and assuming zero flux trapped in the ring, we obtain (212 ± 0.6) kHz, which is in good agreement with the 209kHz predicted by reference [23]. We also calculated the inductance of such a superconducting ring with flux quantization with FEM and obtained Figure D6. Trap frequency of a spherical particle with 125 nm radius and λ L = 50 nm in an AHC trap as considered in figure 6 with a scaling factor of 10, as a function of the maximal mesh element size l mesh . Comparison of (a) FEM-2D to analytical results obtained for the configuration of a superconducting sphere in a quadrupole field [22] for two different meshing strategies: (i) triangular mesh only (red data) and (ii) triangular mesh combined with a shell mesh that meshes the outermost volume of 75 nm thickness of the sphere with onion-type layers of 1 nm thickness (blue data). (b) Degrees of freedom (dots) and average element area (crosses) in the sphere for each of the meshing strategies.

Appendix D. FEM meshing
Given that the model is based on FEM, the results are mesh dependent. Constructing a mesh fine enough at the surface of the superconducting domains is critical to get reliable results. This dependence is illustrated in figure D4, where the trap frequency along z for a 1 µm diameter sphere in an AHC trap (cf figure 3) is calculated via FEM-3D. For these simulations we changed the maximal allowed mesh element size, l mesh , on the surface of the particle resulting in gradually finer meshed particles, see the insets in figure D4 and figure G5. When reducing l mesh , the FEM meshing algorithm gradually increases the number of mesh elements in the sphere and, thus, reduces the average element area that one mesh element covers. This is reflected in the number of degrees of freedom (DOF) in the sphere, that is, the number of unknowns to solve for in the model, which in general equals the number of dependent variables (A x , A y , A z and the gauge fixing potential inside the sphere) times the number of nodes in the geometry. In all Figure E7. Trap frequency of a superconducting sphere in a realistic AHC predicted by the image method [24], by assuming a superconducting sphere in a quadrupole field [22], FEM-2D with 1D-wires, FEM-2D and FEM-3D. (a) The radius of the particle is scaled and the geometrical parameters of the trap and λ L are kept constant with parameters as given in table 2. (b) The geometrical parameters of the trap and particle are taken from table 2 and scaled by a scaling factor while the current density in the coils and λ L are kept constant. The vertical lines indicate the initial geometry. The black points in the schematics indicate the location of the 1D-current loops. The grey area represents geometries in which the particle is subject to magnetic fields above 80 mT (Bc of lead) with a maximal field of up to 100 mT. our simulations we use quadratic mesh discretization, which means the lines connecting the mesh nodes are not straight lines but polynomials of second order. For l mesh ⪅ 5 · λ L = 250 nm, we observe no clear trend of the trap frequency within its uncertainty. 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.
For FEM-2D we can decrease l mesh further as the computational cost is not as large as for FEM-3D simulations.  an AHC trap in dependence of l mesh . For fine enough meshing, i.e. l mesh ⪅ 10 nm corresponding to > 10 4 DOF, the FEM simulations converge to the analytical results obtained for a superconducting sphere in a quadrupole field. The small discrepancy is attributed to the difference between the field distribution of a quadrupole field and the field of the modeled trap.
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 for similar FEM models. To test this, we simulated the trap configuration as used for figure D4 for slightly different l mesh of (49.9, 49.95, 50.00, 50.05, 50.1) nm and get trap frequencies of (23.7, 23.6, 23.8, 23.6, 23.6) kHz, resulting in a mean value of (23.66 ± 0.09) kHz. Thus, the scatter of trap frequency of about ± 0.4% from using nearly similar meshes is smaller than the fit uncertainty of the trap frequency. Note that the computation time for obtaining a typical magnetic field distribution of a particle is 30 − 120 min and requires 50 − 600 GB of RAM on computing nodes with 20 core Intel E5-2650v3 CPUs with 2.30 GHz base frequency available via a computing cluster. Figure E7 and figure E8 show the dependence of the trap frequency with the scaling of the geometry of the respective trap. Here, we place the 1D current loops in FEM2D-1D and the perfect diamagnetic particle method at a position corresponding to the center of the wires. This data can be compared to the corresponding data shown in figure 6 and figure 9 when the 1D current loops are placed at the innermost corner of the coils. Figure E9 shows the magnetic field distribution in an AHC trap for the case of a superconducting sphere in a quadruple field [22] and in the field generated by quasi-1D wires (obatined via FEM-2D-1D). These field distributions are similar close to the particle surface, but deviate much more when approaching the coil wires.