High-Q magnetic levitation and control of superconducting microspheres at millikelvin temperatures

We report the levitation of a superconducting lead-tin sphere with 100 micrometer diameter (corresponding to a mass of 5.6 micrograms) in a static magnetic trap formed by two coils in an anti-Helmholtz configuration, with adjustable resonance frequencies up to 240 hertz. The center-of-mass motion of the sphere is monitored magnetically using a dc superconducting quantum interference device as well as optically and exhibits quality factors of up to 2.6e7. We also demonstrate 3D magnetic feedback control of the sphere's motion. The setup is housed in a dilution refrigerator operating at 15 millikelvin. By implementing a cryogenic vibration isolation system we can attenuate environmental vibrations at 200 hertz by approximately seven orders of magnitude. The combination of low temperature, large mass and high quality factor as well as adjustable resonance frequencies provides a promising platform for testing quantum physics in previously unexplored regimes with high mass and long coherence times.

Diamagnets and superconductors partially expel magnetic fields, allowing stable levitation in field minima [1]. A prominent application of magnetically levitated systems is the use as ultraprecise acceleration sensors, most notably in the superconducting gravimeter [2], which relies on the levitation of centimeter-sized hollow superconducting spheres in a stable magnetic field generated by superconducting coils in persistent current mode. More recently, several proposals [3,4] have highlighted the potential of magnetic levitation for tests of quantum physics with macroscopic, micrometer-sized objects.
Successfully preparing a magnetically levitated particle in a quantum state requires a magnetic trap with low damping, low heating rates and the ability to control the mechanical motion. Some of these features are already present in recent demonstrations of magnetically levitated systems, such as the levitation of permanent magnets above a superconducting surface [5,6] and the levitation of diamagnets in the field of permanent magnets at cryogenic temperatures [7] or room temperature [8]. However, as of yet no system has combined these features. Furthermore, as typical levitation frequencies for these systems range from subhertz to kilohertz, the heating rate is dominated by environmental vibrations.
One of the most promising avenues toward the quantum regime is the levitation of a type-I (or zero-fieldcooled type-II) superconductor in a magnetic field produced by persistent currents, as this avoids dissipation due to hysteresis and eddy currents, which is inherent to levitation schemes involving ferromagnets or flux-pinned type-II superconductors [9,10]. Using persistent currents should result in an extremely stable trap, as both the drift and noise can be kept extremely low [2,11,12]. Working at millikelvin temperatures also naturally enables coupling to nonlinear quantum systems such as superconducting qubits [3], and miniaturizing the trap architecture [13][14][15] might lead to fully on-chip coupled quantum systems [16]. In comparison with optical levitation, which has already been established as a means for studying massive objects in the quantum regime [17][18][19], magnetostatic levitation has the potential to access a new parameter regime of even larger mass [20,21] and longer coherence times [3,16]: While the size of optically levitated particles is, in practice, limited by the available laser power to the micrometer scale [22], magnetic levitation can support train-scale objects [23] and due to the static nature of the fields there is no heating from photon absorption.
A first step toward a fully superconducting platform was presented in [24], where a lead particle was levitated in the magnetic field of a superconducting coil. This approach, however, was limited in performance by the presence of cryogenic exchange gas at 4 K and environmental vibrations. Here, we report an experiment that provides orders of magnitude improvements in both damping and vibration isolation. We demonstrate the levitation of a superconducting sphere in the magnetic field generated by a superconducting coil. We detect the center of mass (c.m.) motion of the particle magnetically by inductive coupling to a superconducting quantum interference device (SQUID) via a pickup coil (cf. Fig. 1) and demonstrate feedback control of all three c.m. translational degrees of freedom by real-time processing of the SQUID signal. We find that the motion of the particle has a low damping rate, corresponding to quality factors of up to 2.6 × 10 7 . We also implement a custom vibration isolation system to reduce heating from environmental vibrations. Finally, we discuss the limitations of our current setup and the improvements that are necessary for ground state cooling.
Trapping.-The energy of an object of volume V with permeability µ in an applied magnetic field B 0 in free space can be approximated by −1/2V (µ − µ 0 )B 2 0 [1], where µ 0 is the permeability of free space. Because ∇ 2 B 2 0 ≥ 0, stable levitation is possible in a field minimum for µ < µ 0 , i.e. for diamagnets. A superconducting sphere can act as a perfect diamagnet (µ = 0), because it reacts to a change in the applied magnetic field by forming screening currents on its surface, counteracting the applied field and keeping the interior of the sphere field free.
For a sphere trapped in the field of an anti-Helmholtz coil, the resulting potential close to the field minimum is harmonic and two of the c.m. modes are degenerate [3,20]; to lift the degeneracy we use elliptical coils [25]. The z direction is along the coil axis, which coincides with the vertical direction, while x and y are perpendicular to z, with x along the major axis of the coils. We use numerical simulations [25] to model the magnetic field created by such an arrangement and confirm that the magnetic field near the center between the coils is well described by (b x x, b y y, b z z), where the field gradients b i are pro-portional to the trap current with |b x | < |b y | < |b z | and b z = −(b x + b y ). This field shape constitutes a threedimensional harmonic trap for a levitated superconducting sphere with the c.m. frequencies [25] where ρ is the density of the sphere. Our setup is housed in a dilution refrigerator (Bluefors BF-LD400) with a loaded base temperature of approximately 15 mK. The microspheres are commercial (Easy-Spheres) solder balls with diameter 100(6) µm, made of a 90-10 lead-tin alloy with density ρ = 10.9 × 10 3 kg/m 3 . Lead-tin is a type-II superconductor with a critical temperature of approximately 7 K [26]. A single sphere is initially placed in a 3D-printed polylactide bowl glued to the lower trap coil, such that it rests around 1 mm below the trap center. The bowl helps to prevent accidental loss of the sphere during assembly and between measurement runs. After cooling down the system to millikelvin temperatures, lift-off is achieved by increasing the coil current until the upward magnetic force on the sphere is stronger than the downward gravitational and adhesive forces. The trap coils are surrounded by an aluminum shield with small openings for optical access and wires, which screens the sphere from magnetic field fluctuations when the aluminum is superconducting. To suppress field fluctuations from the trapping coils and feedback coil, we use a low pass filter and attenuation stages, respectively [25].
Readout.-A planar Niobium thin film gradiometric pickup coil consisting of two square loops positioned approximately 400 µm above the trap center is used to magnetically probe the particle motion. Each loop is 145 µm × 145 µm and the separation between the loops is 3 µm. As the particle moves, it induces a current in the pickup loop [25], which is connected with Nb wires to the input coil of a SQUID current sensor [27]. The SQUID has an intrinsic flux noise level of S Wb is the flux quantum. When the SQUID is incorporated into the setup and the pickup loop is connected, the noise rises to S 1/2 Φ ≈ 10 µΦ 0 / √ Hz. We believe this is due to external field fluctuations inductively coupling into the wires connecting the pickup loop and the SQUID. With a commercial SQUID (Supracon CSblue), that was used for some of the measurements, the noise floor in the setup further increased by an order of magnitude.
To calibrate the SQUID signal, we optically record the motion of the sphere by illuminating it and imaging its shadow on a camera, while simultaneously recording the motion with the SQUID [25]. The coupling strengths (flux induced in the SQUID per sphere displacement) for the different motional modes depend on the position of the pickup loop in relation to the sphere as well FIG. 2. The PSD of the SQUID signal displays three peaks which correspond to the c.m. modes of the particle. The motion of the particle is excited to approximately one micrometer rms displacement in this figure -measurement noise prevents us from resolving the equilibrium motion of the particle. The inset shows the linear dependence of the c.m. frequencies on the trap current. The solid lines are zero-intercept linear fits.
as the trap gradients, i.e. mechanical frequencies. The highest measured coupling strength was approximately 13 × 10 3 Φ 0 /m for the z mode and 3 × 10 3 Φ 0 /m for the x and y modes. Together with a flux noise floor of 10 µΦ 0 / √ Hz these coupling strengths correspond to a displacement sensitivity of approximately 1 nm/ √ Hz for the z mode. During these measurements the separation between pickup loop and sphere was approximately 400 µm. We estimate that, with a smaller separation and an optimized pickup loop geometry, an improvement of 4 orders of magnitude is possible [25].
-The c.m. frequencies of the particle are visible in the power spectral density (PSD) of the SQUID signal, as shown in Fig. 2. The c.m. frequencies depend linearly on the magnetic field gradient of the trap [see Eq. (1)] and, thus, the current in the trap coils, as shown in the inset. The highest frequency reached in the present setup is 240 Hz (corresponding to b z ≈ 150 T/m), limited by the maximum current output of the low-noise current source. We can trap at lower frequencies, down to f z ≈ 20 Hz, where the lower limit stems from the gravity-induced shift of the trapping position and our trap geometry. However, at low trapping frequencies, the sphere is strongly excited by environmental vibrations.
We can apply a magnetic feedback force on the levitated sphere by processing the SQUID signal on an FPGA (STEMlab 125-14) and applying a feedback current to a small coil with approximately 20 windings positioned approximately 1 mm below the trap center [25]. The upper right inset shows a linear plot for the latter part of the data, starting at t=2.5 h. The red line is a fit to an exponential decay. The spikes correspond to a brief loss of the SQUID lock, they were removed for the amplitude plot. The lower left inset shows a ringdown while direct feedback is applied (note that the x axis scale is in seconds).
As measurement noise prevents us from resolving the equilibrium c.m. motion of the particle, we use feedback here only as a means to quickly adjust the amplitudes and prepare other measurements (cf. Fig. 3, inset).
Quality factor.-We measure the quality factor by performing ringdown measurements, where the initial starting amplitude is set by applying an appropriate feedback signal before the start of the measurement. We do see occasional jumps in the damping rate (cf. Fig. 3), which we attribute to the detaching and attaching of flux lines to a pinning center. Unpinned flux lines can move within the superconductor (flux creep), a mechanism that is known to cause damping in the levitation of flux-pinned superconductors or magnets [10,28]. Although we do not apply a magnetic field during the cooldown, there is a nonzero background field, and hence, we expect some frozen-in flux to be present in the levitated sphere. As a consequence, the measured quality factors vary with time and between measurement runs, but they are generally between 1 × 10 7 and 2.5 × 10 7 for the z mode and about half as high for the x and y modes. The highest quality factor we have measured is 2.6 × 10 7 at 212 Hz, corresponding to a dissipation rate of 5 × 10 −5 s −1 , and likely still limited by flux creep inside the particle. We also consider other possible contributions to the damping [25], but find that the expected contributions cannot explain the measured values. In future experiments we plan to use a particle made of a type-I superconductor such as monocrystalline lead, which will expel all magnetic flux via the Meissner effect and prevent flux creep.
Sensing.-A massive system acting as a harmonic oscillator can be used for force and acceleration sensing, with the sensitivities depending on the mass m and the dissipation rate γ of the oscillator. In equilibrium with a thermal bath at temperature T 0 , the thermal force noise is given by S th , where k B is Boltzmann's constant. The dissipation rate γ is related to the quality factor Q and the angular mechanical resonance frequency ω 0 by γ = ω 0 /Q. The corresponding acceleration sensitivity is S th F F /m. For a particle of mass m ≈ 5.6 µg, assuming it is in thermal equilibrium with its surroundings in the dilution refrigerator at 15 mK, this sets the limits for force and acceleration sensing at 5 × 10 −19 N/ √ Hz and 9 × 10 −12 g/ √ Hz, using our highest measured quality factor Q = 2.6 × 10 7 at 212 Hz. Our current experimental sensitivities are approximately 1 order of magnitude above the thermal noise, limited by drifts in the trap current and measurement noise [25].
Cryogenic vibration isolation.-Vibrations of the cryostat accelerate the trapping coils and, thus, the trap center, effecting a force onto the levitated sphere. Denoting the vibrational PSD by S ǫǫ , the displacement PSD of our sphere becomes |χ(ω)| 2 m 2 ω 4 0 S ǫǫ (ω), where χ denotes the mechanical susceptibility. From independent accelerometer measurements on the cryostat, we estimate √ S ǫǫ ≈ 1 × 10 −10 m/ √ Hz at 200 Hz, which would result in a peak height of 2.6 × 10 −3 m/ √ Hz and a root-meansquare displacement of 9 µm, corresponding to an effective temperature of 6 × 10 10 K. To mitigate the effects of vibrations, we implement a passive cryogenic vibration isolation system: the aluminum shield containing the trap is hung from the 4 K stage of the dilution refrigerator via 38 µm-thick stainless steel wires and two intermediate stages. The system acts like a triple pendulum, offering isolation from horizontal (x and y) vibrations, while the elasticity of the wires provides isolation from vertical (z) vibrations. To prevent coupling of external vibrations into the experimental setup via electrical wires and the copper braids used for thermalization, we connect these wires and braids to an additional vibration isolation platform, as represented in Fig. 4. Vibrations of the experimental system could also be induced by fluctuating magnetic fields acting on the aluminum shield. To prevent this, we surround the setup by a second aluminum shield which is rigidly mounted to the cryostat. We typically operate at vertical (z) trap frequencies above 200 Hz and initially found that the vibration isolation system attenuates vibrations at these frequencies by approximately 5 orders of magnitude. We have further improved the system by optimizing the mass and wire length for each stage, leading to an attenuation of almost 7 orders of magnitude at 200 Hz. Theoretically, vibrations along the horizontal axis as well as librational fluctuations around the axes are suppressed even more, but in practice, we expect vibrations from the vertical axis to couple into all Prospects for ground state cooling.-Now, we discuss the potential of this system to reach the ground state, which would be a first step toward accessing the quantum regime. SQUID noise can be separated into flux noise S φφ and noise in the circulating current S JJ , which may be partially correlated [31] and fulfill S φφ S JJ − S 2 φJ ≥h [32, 33]. The flux noise results in the measurement noise S φφ η 2 , while the current noise corresponds to a back action η 2 S JJ , where η is the coupling strength. Applying direct feedback with optimum gain and assuming vanishing correlations between flux and current noise, one can estimate the final phonon occupation [25] as SJJ depends on the working point of the SQUID and is typically on the order of the SQUID inductance L S . The first term describes heating from coupling to an effective thermal reservoir with temperature T 0 , while the second term is the backaction limit under continuous feedback. We can reach a coupling strength of η ≈ 5.5 × 10 7 Φ 0 /m [25], so assuming T 0 ≈ 15 mK and a typical SQUID inductance of L S ≈ 15 pH, we need γ ≈ 1 × 10 −6 s −1 to keep the former negligible. Ground state cooling then requires S φφ S JJ < 3h. In addition to feedback cooling with a SQUID, we will investigate other potential avenues for ground state cooling such as coupling the levitated sphere to a flux qubit or a microwave cavity [3,[34][35][36].
Discussion and outlook.-We built a platform based on the diamagnetic levitation of a superconductor in a magnetostatic trap in a dilution refrigerator at 15 mK and demonstrate quality factors exceeding 1 × 10 7 for the three c.m. modes of a 5.6 µg oscillator. The high quality factors result in excellent sensing capabilities, with achievable force (acceleration) sensitivities that are, otherwise, only reached with much smaller (larger) masses [2,37]. We also demonstrated simultaneous feedback cooling of all c.m. modes using either direct or parametric feedback. While we have focused, here, on a 100 µm-sized sphere, the setup allows for particles from submicron size to millimeter size to be levitated [20]. We have identified improvements to mitigate technical issues: First, we plan to use a type-I superconducting sphere to avoid trapped flux and dissipation caused by flux creep, which we expect is limiting our quality factor. Second, we can improve the magnetomechanical coupling by 4 orders of magnitude by better positioning of the pickup coil and by using a pickup coil with multiple windings matched to the input coil of the SQUID. Together with an optimized SQUID shielded from environmental noise (such that our readout is limited by the intrinsic SQUID noise), this would result in an improvement of the measurement noise floor by 6 orders of magnitude [25]. Third, we plan to implement a persistent current switch to improve the trap current stability. With these improvements, our system offers a promising approach for bringing microgram objects to the quantum regime, which opens a potential avenue for quantum-limited acceleration sensing or even probing quantum effects of gravity [38][39][40][41]. The data of this study are available at the Zenodo repository [42].
We thank G. Kirchmair, O. Romero-Isart, and A. Sanchez for discussions. We thank Matthias Rudolph for assisting with the initial setup of a SQUID measurement system. We thank F. Wulschner for setting up a helium flow cryostat and assisting with measurements during the initial phase of the project.

Trapping field
We performed numerical simulations in comsol multiphysics to predict the field configuration generated by the arrangement of our trapping coils. The design of the coils and the simulated magnetic field density are shown in Fig. 1a,b. The field gradients for a trap current of 2.5 A are obtained as (b x , b y , b z ) = (57, 90, 147)T/m, which would correspond to frequencies of (92, 144, 236)Hz. Comparing this to the measured values of (109, 127, 236)Hz we can see that the agreement between the numerical and measured z-mode frequency is excellent, but that the simulations predict a much stronger split for the x-and y-mode than what we see in the experiment. This is likely due to the fact that the spooling process resulted in coils that were more circular than the design (cf. Fig. 1c). For a perfectly circular coil the x-and y-mode frequencies are degenerate [1], so a smaller frequency split is to be expected. Close to the center the field is very well approximated by a quadrupole field: Within a cubic volume of 200 µm×200 µm×200 µm the maximum relative root-mean-square deviation, defined as |B sim −(b x x, b y y, b z z)|/|B sim |, is less than 0.1 %.

Mechanical frequencies
In this section we build on the calculations in [1] to analytically calculate the magnetic field distribution for a su-perconducting sphere in an axially-asymmetric magnetic quadrupole field. We consider a superconducting sphere at the origin of the coordinate system and displaced from the center the quadrupole field by x 0 = (x 0 , y 0 , z 0 ). The applied field then reads As our sphere's radius R is much larger than the expected penetration depth λ [2], we can use the approximation λ/R → 0, i.e. we assume B = 0 inside the sphere and B r = 0 on the sphere's surface [1], where B r denotes the radial component of the field. As there are no currents outside the sphere, there exists a scalar potential Φ s.t.
where r = x 2 + y 2 + z 2 and the Y m n are spherical harmonics. The coefficients a n,m are determined by the boundary condition B 0,r = (∇Φ) r on the sphere's surface as all other coefficients are zero. From this field distribution we obtain the force on the sphere as where V is the volume of the sphere. The magnetic field thus creates a harmonic trapping potential for a superconducting sphere, with trapping frequencies described by

Coupling strength
The linear coupling strength with respect to a pickup loop that is described by a closed path γ is defined as ν i = ∂ i Φ γ , where Φ γ is the flux through the area enclosed by the pickup loop and i ∈ {x 0 , y 0 , z 0 }. We can express this as where A denotes the magnetic vector potential, defined by ∇ × A = B and ∇A = 0. While the resulting expressions are generally bulky, they can be easily evaluated using a computer algebra system, we use Wolfram Mathematica. We measure the current that is induced in the pickup loop using a SQUID current sensor, i.e. the pickup loop is connected to an input coil on the same wafer as the SQUID and inductively coupled to the SQUID. The coupling strength with respect to the SQUID is thus obtained as η i = ν i M/L, where L is the inductance of the circuit formed by pickup loop, input coil of the current sensor, and the connecting wires, while M is the mutual inductance between input coil and SQUID. In our case M/L ≈ 0.02. In Fig. 2 we plot the dependence of the coupling strength on the vertical separation between trap center and pickup loop, when the x-and y-position of the pickup loop relative to the trap center (as shown in the inset) is approximately the same as what we estimate for our setup (from pictures taken at room temperature, cf. Fig. 1c). The horizontal dotted lines correspond to the measured coupling strengths, while the vertical dotted line corresponds to the vertical separation extracted from camera images (approx. 400 µm) of the levitated sphere.
Trap displacement due to gravity The rest position of the sphere is slightly displaced from the field minimum due to gravity. In our case the setup is aligned such that the direction of gravity coincides with the z-axis. The expected displacement is thus where g ≈ 9.81 m/s 2 . In general this displacement needs to be taken into account when calculating the coupling strength, but for larger trapping frequencies it becomes negligible. We usually operate at f z > 200 Hz, for which z g < 6 µm.

CURRENT FLUCTUATIONS
Since the sphere's motional frequencies depend linearly on the trap current, current fluctuations result in frequency fluctuations, and thus can cause heating as well as broadening of the mechanical spectral peaks. We stabilize the trap current both passively and actively, as described in the following. We implement a first-order low-pass RL-filter by adding a short copper wire (resistance R C ) in parallel to the trap coils (inductance L C ). We characterize the filter's cut-off frequency by measuring its step response in an empty trap (Fig. 3). The cutoff frequency is R C /L C ≈ 0.036 s −1 , which results in an attenuation of approximately −180 dB at 200 Hz, while the added thermal current noise from the copper wire is only on the order of 1 × 10 −11 A/ √ Hz. We actively stabilize the filtered trap current by monitoring the mean current source output using a digital ammeter and updating the current source output based on an estimation of the trap current. This lessens current drifts and the associated resonance broadening, as seen by comparing Fig. 4(a) with (b). Even with active current stabilization we observe drifts of the particle's resonance frequency shown in Fig. 4(d), which we attribute to low-frequency noise in the ammeter reading.

SENSING
Taking into account measurement noise, we can calculate the force sensitivity as Hz is the measurement noise and χ(ω) = 1/[m(ω 2 0 − ω 2 − iγω)] denotes the mechanical susceptibility. In our case, using Q = 2.6 × 10 7 at 212 Hz and assuming T 0 = 15 mK, this would result in a force sensitivity of √ S F F = 6.3 × 10 −19 N/ √ Hz on resonance, close to the thermal limit of 4.9 × 10 −19 N/ √ Hz, and a noise equivalent temperature of T n = T 0 + |χ(ω 0 )| −2 S nn /(4k B mγ) = 24 mK (cf. Fig. 5). However, due to the stochastic frequency fluctuations described in the last section the susceptibility is stochastic as well. We briefly sketch how this reduces our sensitivity: Note that a measurement with frequency drifts can be approximated as a series of shorter measurements at different, but constant, resonance frequencies. For a finite measurement duration T the measured force power spectral density at resonance is approximately where the integration runs over the bin width ∆f = 1 T of the bin containing the resonance. Our force sensitivity thus becomes √ S F F = 4k B T 0 mγ + |χ| −2 avg S nn with |χ| 2 avg = 1 ∆f ∆f df |χ(f )| 2 . For ∆f > γ 2π , we get |χ(f0)| 2 |χ| 2 avg ∝ ∆f γ . In our case this boosts the measurement noise above the thermal noise and, for a typical measurement, our force sensitivity worsens by at least an order of magnitude, corresponding to an noise equivalent temperature of at least 2.5 K. Both improving the coupling, which corresponds to an effective decrease of the measurement noise, and improving the stability of the resonance frequency will allow us to perform measurements limited by thermal noise. We show below that we can increase the coupling by almost four orders of magnitude. Regarding the latter, we note that persistent current coils have reached a relative current stability better than 2 × 10 −11 /h [3].

FEEDBACK
We can apply a magnetic feedback force on the levitated sphere by processing the SQUID signal and applying a feedback current to a small coil with approximately 20 windings positioned approximately 1 mm below the trap center. The feedback current generates an additional magnetic field, thereby shifting the field minimum of the trapping field and enabling us to apply direct feedback. The gradient of the trapping field is changed as well, such that we also have the possibility to apply parametric feedback at twice the particle's resonance frequencies. Since the COM modes are separated in frequency, we can apply feedback to all modes simultaneously. When we apply direct feedback, we pass the SQUID signal to a FPGA (STEMlab 125-14), which, for each of the mechanical frequencies, applies a bandpassfilter, gain and phase shift [4]. The resulting signals are then recombined on the FPGA and fed back to the feedback coil. When we apply parametric feedback at twice the mechanical frequencies we additionally use a phaselocked loop and a clock divider on a second FPGA. As described in the last section, our current noise floor corresponds to a noise equivalent temperature that is higher than the base temperature of our dilution refrigerator. This noise equivalent temperature also sets the limit for feedback cooling. Due to the large ambient noise reduction provided by our vibration isolation, magnetic shielding and trap current filter, the equilibrium temperature of the particle is below the noise equivalent temperature, meaning that even without applied feedback the particle's motion undergoes exponential decay until it is not detectable anymore. We thus cannot cool the particle's modes below their equilibrium temperature and use feedback currently only as a way to quickly adjust the amplitudes and prepare other measurements.

CHARACTERIZATION OF THE VIBRATION ISOLATION SYSTEM
Our vibration isolation system consists of several cylindrical plates connected by straight wires, with the top plate mounted to the 4K-platform of the dilution refrigerator and the bottom plate mounted to the aluminum can containing the magnetic trap. Electrical and thermal connections going to the setup are loosely coiled up and guided to the setup via a separate vibration isolation stage hung from the Mixing Chamber platform (cf. Fig. 6). We first characterized a single stage at room temperature, by exciting it mechanically and recording its motion with a camera. We find that the horizontal and vertical resonance frequencies are f h = 1.1 Hz and f v = 18.7 Hz, respectively, while the librational resonance frequencies are f l,1 = 0.6 Hz (around vertical axis), f l,2 = 8.7 Hz and f l,3 = 9.0 Hz (around horizontal axes). We are thus limited by vertical vibrations, which in general couple to the other degrees of freedom as well [5]. We then assembled the setup as shown in Fig. 6, with the bottom plate mounted to an accelerometer (instead of the setup) used to measure acceleration along the vertical axis. The mass of the bottom plate (including the accelerometer) was 0.29 kg in this assembly. The intermediate stages are each supported by three wires of equal length, while the bottom stage is supported by a single wire. The wires are made from type 304 stainless steel with a diameter of D = 38 µm and a yield load of 0.37 kg, as specified by the manufacturer (Fort Wayne Metals). The wires are connected to the stages with clamping connections. The vertical spring constant for the i-th stage is given by k v,i = (N i Y D 2 π)/(4L i ), s.t. f v,i = 1/(2π) k v,i /m i is the resonance frequency of a single stage. Here Y is the Young's modulus, N i the number of supporting wires, L i the wire length and m i the mass of the stage. The normal mode frequencies f n,i of the assembled system are obtained by solving the characteristic equation of the system [5], resulting in the transfer function ( In Fig. 7 we compare the vibrations measured without vibration isolation (with the accelerometer mounted to the Mixing Chamber stage of the dilution fridge) and with vibration isolation. The highest normal mode frequency is approximately 30 Hz, above which the measured accelerations quickly drop below the sensitivity of the accelerometer, when it is mounted to the vibration isolation system. At lower frequencies the measured transfer function corresponds very well to the analytical values and we use the analytic expression to extrapolate the attenuation at higher frequencies, resulting in an attenuation of 1 × 10 −5 at 100 Hz and 1.5 × 10 −7 at 200 Hz. The peaks background in the lab, which is verified by an independent measurement in which the accelerometer is replaced by an 1 MΩ resistor inside the fridge. This measurement reproduces the peaks at precisely the same frequencies as in the accelerometer measurement, although with a lower magnitude -this we attribute to the additional (unshielded) wiring going from the SMA port to the accelerometer, which can act as an antenna. In an initial approach we used different masses and wire lengths, such that the attenuation was approximately 5 orders of magnitude at 200 Hz.

QUALITY FACTOR
In this section we estimate contributions to the dissipation due to gas damping and eddy currents. We also consider dissipation in the SQUID and flux creep in the pickup loop, which we rule out experimentally.
Collisions with background gas result in an additional damping of [6] γ P = β P ρRv th , where β ≈ 1.8 andv th is the mean thermal velocity of the molecules. We do not have a direct measure of the pressure at the trap location, but we can use the measured pressure at the room temperature side of the vacuum can, P ≈ 1 × 10 −6 mbar, to set an upper limit for the expected damping. Assuming the background gas consists mostly of Helium, we get γ P ≈ 3 × 10 −7 s −1 , two orders of magnitude below the measured damping rates. We avoid eddy current damping by not placing any normal conductors within the inner shield. Eddy currents can still be induced in conductors outside the shield, either via openings in the shield or via the pickup circuit, part of which is placed outside the shield. In the first case, we can estimate possible eddy current losses by imagining a conductive loop with resistance R o and inductance L o placed around the openings in the shield. The dissipation in mode i ∈ {x 0 , y 0 , z 0 } then takes a maximum value for R o = 2πf i L o , corresponding to a damping Here φ denotes the flux in the loop induced by the oscillation of the sphere, which can be calculated from Eq. 1. Using the dimensions of our windows and estimating L o > 10 nH for a similarly sized conductive loop, we get γ i < 1 × 10 −9 s −1 for all modes.
We now consider eddy current damping mediated by the pickup circuit or damping due to flux creep in the pickup circuit, as well as dissipation in the SQUID due to a real part in the SQUID input impedance [7]. This implies γ i ∝ ν 2 i , while in the experiment the damping is approximately equal for all three COM modes. In addition, during initial measurements the pickup loop was positioned farther from the trap center and coupling strengths were approximately one order of magnitude lower than the ones reported here, but we did not see any change in the quality factor. of light, which is further collimated by small windows in the cryostat and the magnetic shields. After passing through the cryostat, the light is focused by a teleobjective into the camera (DFK 33UX287). Any object in the light path thus appears as a shadow on a bright background, corresponding to the projection of the object onto the image sensor. We use a frame rate of 596fps and 0.5ms exposure time to record the videos and then process them as follows. Each frame is first converted to grayscale, smoothed with Gaussian blurring, thresholded and then analyzed using the particle tracking software Trackpy [8]. This results in datasets for the horizontal and vertical (as defined by the camera's orientation) position of the sphere. The camera is aligned such that the vertical axis is along the trap coil axis and the direction of the particle's z-motion, while the particle's motion along x and y is projected onto the camera's horizontal axis. The trap coils are aligned such that the elliptic axes are at around 45 • to the optical axis, so both the x and y motion can be imaged. To get separate trajectories for each mode we implement a digital bandpass filter around the respective resonance frequency. Calibration of the sphere's displacement relies on knowledge of the optical system's magnification, which we determine using knowledge of the sphere's size (100(6) µm, as specified by the manufacturer). Including the uncertainty from the alignment of the trap coils with respect to the optical axis in addition to the ±6 % uncertainty from the size of the sphere, we estimate a total uncertainty of ±13 % for the calibrated COM displacement.
The light source is turned on only for calibration measurements, as the illumination heats the particle and limits the levitation time to approximately one hour. Further, optical band-pass filters (approximately 400 nm -650 nm) on the millikelvin stage reduce ambient light reaching the particle and heating it. Without the filters (and the light source turned off) levitation times are limited to around four hours, while we have not yet found a limit (>30 h) using the filters.

OPTIMIZATIONS AND GROUND STATE COOLING
In this section we provide more detail on how to implement the improvements that are necessary to perform feedback cooling to the ground state. We will consider only cooling of the z-mode, cooling of the other modes and 3D-cooling can be characterized and optimized in an analogous manner. We first show that the coupling strength can be increased by almost four orders of magnitude, s.t. the measurement noise (in units of m/ √ Hz) with a state-of-the-art SQUID decreases to 1 × 10 −15 m/ √ Hz. We then take into account the back action from the SQUID due to the circulating current and determine the standard quantum limit (SQL) of our system. We proceed to evaluate the requirements for ground state cooling as well as contributions to heating from vibrations and frequency fluctuations.

Coupling strength
We can increase the coupling strength by better positioning of the pickup loop relative to the superconducting sphere and by increasing the number of turns of the pickup loop, thus matching its inductance to that of the input coil. For the calculation, we will assume a planar pickup loop with circular windings, positioned coaxial with the z-axis. Then Eq. 1 leads to for the coupling strength w.r.t. a single turn of the pickup loop with radius R P and z-position Z P . The coupling strength to a multi-turn pickup loop can be written as a sum of terms of this form, ν = i ν 1 (R P,i , Z P,i ), and the coupling strength to the SQUID (cf. Fig. 9) is obtained as Here L P , L I and L S are the inductances of pickup coil, SQUID input coil and SQUID, respectively. We have used that the the mutual inductance M between input coil and SQUID can be written as M = k √ L I L S with |k| < 1. L W denotes the stray inductance of the pickup circuit and is in our case dominated by the wires connecting the pickup coil to the SQUID, L W ≈ 100nH.
A relevant figure of merit for the measurement noise of the SQUID is the so-called energy resolution S EE = S φφ 2LS , with state-of-the-art SQUIDs reaching values on the order ofh [9,10]. In terms of position resolution this can be written as In order to get a realistic estimate for this we will use the SQUID parameters presented in [10], i.e. L S = 15 pH, L I = 0.53 µH, M = 2.3 nH and S φφ 2k 2 LS = 5.5h. For the pickup coil we assume a wire width of 0.3 µm and a distance between the wires of 0.45 µm, which is also the resolution of the input coil of our existing SQUID. We use the modified Wheeler formula [11] to determine L P , s.t. we get a fully analytical expression for S nn . We also restrict us to Z P ≥ R. Given these parameters, as well as b z = 147 T/m, we minimize S nn with respect to R P , Z P and N , where R P now denotes the inner radius of the pickup coil and N is the number of turns. This results in η ≈ 5.5 × 10 7 φ 0 /m and S nn ≈ 1 × 10 −15 m for Z P = R, R P ≪ R and N ≈ 87. We note that the scaling of the coupling (Eq. 3) means that the measurement noise can be further decreased by working with higher gradients or using larger particles (a scale-independent expression for the coupling is easily obtained by measuring distance in units of R, magnetic fields in units of b z R and thus the coupling in units of b z R 2 ).
The standard quantum limit Equation 4 for the readout signal measured with the SQUID is valid only when the back-action stemming from the circulating SQUID current onto the levitating particle is negligible. This is a reasonable assumption given the weak coupling in our current setup, but for increased coupling strength the back-action should be considered. All relevant parameters are stated for the SQUID coupled to the input circuit and can be different than the parameters of the uncoupled SQUID, depending on the capacitive coupling at the Josephson frequency [7]. The readout noise is set by the flux noise, while the backaction noise is determined by fluctuations of the circulating current J. For a shunted SQUID the origin of both flux and current fluctuations is current noise in the shunt resistors and the quantum limit that arises from zero-point-fluctuations is S φφ S JJ − S 2 φJ ≥h [12,13]. We will assume in the following that correlations between flux and current noise are negligible, S φJ ≈ 0 [13,14].
A current J around the SQUID loop will effect a force −ηJ on the levitated sphere and thus the measured displacement PSD of the sphere becomes S m zz = S φφ η 2 + |χ| 2 (S th F F + η 2 S JJ ).
The product of measurement noise and back-action force noise thus fulfills the Heisenberg-like inequality S φφ η 2 η 2 S JJ ≥h 2 and S m zz is minimized for η 2 = 1 |χ| S φφ SJJ , such that S m zz = 2|χ| S φφ S JJ + |χ| 2 S th F F . For S φφ S JJ =h this corresponds to the standard quantum limit.

Feedback cooling
We now assume we provide direct feedback s.t. the effective damping of the oscillator becomes γ + Γ, where Γ is the cold damping added by the feedback [15] and χ(ω) = 1/[m(ω 2 0 − ω 2 − i(γ + Γ)ω)] is the effective susceptibility. The measurement noise will also enter the feedback system, resulting in an additional force, and the displacement power spectral density becomes S zz = |χ| 2 S th F F + η 2 S JJ + (mωΓ) 2 S φφ η 2 .
For Γ = η 2 mω0LS ≫ γ this corresponds a mean phonon number where we have introduced the shorthandL S = S φφ SJJ . Close to the optimum working pointL S ≈ L S [14]. Using η = 5.5 × 10 7 Φ 0 /m as well as assuming T 0 = 15 mK, m = 5.6 µg and L S = 15 pH, we need γ ≈ 1 × 10 −6 s −1 to keep the left term negligible. Ground state cooling then requires S φφ S JJ < 3h. While current noise in SQUIDs is challenging to measure [16], a DC-SQUID with S φφ S JJ ≈ 10h has been reported [17]. Further, forL S ≈ L S we would expect S φφ S JJ ≈ 2S EE . The lowest energy resolution reported up to date is S EE ≈ 1.7h [9], so ground state cooling by continuous feedback will require a state-of-the-art SQUID.

Heating
We now consider fluctuations of the trap center with spectral density S ǫǫ as well as fractional fluctuations of the spring constant with spectral density S δδ . We note that the heating rate Γ δ E due to frequency fluctuations is proportional to the energy E of the oscillator. The average energy of the oscillator is then described bẏ whereQ ǫ is the heating rate due to fluctuations of the trap center. For Γ δ < γ this results in an effective temperature of The heating rates are given by [18] Q ǫ = 1 4 mω 4 0 S ǫǫ (ω 0 ), and Γ δ = 1 4 ω 2 0 S δδ (2ω 0 ).
Using the values from the preceding paragraph, ground state cooling requires approximately S δδ (2ω 0 ) < 1 × 10 −7 / √ Hz and S ǫǫ (ω 0 ) < 1 × 10 −19 m/ √ Hz (corresponding to T ef f < 20 mK). Regarding the former we note that relative fluctuations on the order of one part per billion can be reached in superconducting coils [19], while the latter requires further improvements to our vibration isolation system. Suspending the system from the top plate of the cryostat would allow us to add two more stages and achieve the required vibration suppression.