Superconducting microsphere magnetically levitated in an anharmonic potential with integrated magnetic readout

Magnetically levitated superconducting microparticles offer a promising path to quantum experiments with picogram to microgram objects. In this work, we levitate a 700ng $\sim 10^{17}$amu superconducting microsphere in a magnetic chip trap in which detection is integrated. We measure the particle's center-of-mass motion using a DC-SQUID magnetometer. The trap frequencies are continuously tunable between 30 and 160 Hz and the particle remains stably trapped over days in a dilution refrigerator environment. We characterize motional-amplitude-dependent frequency shifts, which arise from trap anharmonicities, namely Duffing nonlinearities and mode couplings. We explain this nonlinear behavior using finite element modelling of the chip-based trap potential. This work constitutes a first step towards quantum experiments and ultrasensitive inertial sensors with magnetically levitated superconducting microparticles.

We levitate a superconducting microsphere in a fully chip-based system. The chip-based approach [44,58,59] enables higher magnetic field gradients and trapping frequencies, as well as the potential to scale up the system to levitate multiple particles on the same chip. We measure the particle's motion using magnetic pickup loops which are coupled to a SQUID magnetometer. The pickup loops are integrated in the chip; this allows for precise positioning and enhanced measurement sensitivity. In the future we will replace the SQUID by a flux-tunable superconducting microwave cavity [60][61][62][63][64][65] to achieve quantum control over the COM motion of the levitated microparticle [10,[47][48][49].
In this work we demonstrate stable levitation of a 48 µm-diameter (700 ng) superconducting microsphere over days. We smoothly tune the particle's COM frequencies between 30 and 160 Hz by varying the trap current. We observe that the COM frequencies depend on the motional amplitudes. This arises from trap anharmonicities [32,[66][67][68][69][70]. The observed behavior is consistent with estimations of the trap anharmonicities extracted from finite element modelling (FEM) of our system. In the future we will employ cryogenic vibration isolation [71,72] and feedback cooling [26,27,73] to reduce the motional amplitudes, then the effects of trap anharmonicities will be mitigated.

II. EXPERIMENTAL SETUP
The magnetic trap is produced by a current flowing in two superconducting coils, which are patterned on two silicon chips. The chips are stacked on top of each other to form an anti-Helmholtz-like configuration (see Fig. 1), for details see Ref. [44]. The superconducting particle stably levitates near the minimum of the trap's magnetic field. The particle is confined within a closed container given by the side walls of the hole in the top chip, the top surface of the bottom chip, and a glass slide on top of the hole. The force due to the magnetic trapping field is restoring within the container's bounds, with a trap depth larger than 1 × 10 10 K. To detect the particle motion, we use two pickup loops which are integrated on the same chips, see Fig. 1.
As the particle moves it changes the magnetic flux threading the pickup loops and thus the current induced in the loops. The pickup loops are connected to a commercial DC-SQUID magnetometer, which transduces the flux into a measurable voltage. Typical pickup efficiencies are {η x , η y , η z } = {1.58, 3.3, 19.4} mϕ 0 µm −1 for the three COM modes, in terms of the flux coupled into the SQUID, and where ϕ 0 is the magnetic flux quantum. η z > η x , η y due to the system's geometry (see Appendix S2). The measurement noise floor (0.32 mϕ 0 Hz −0.5 ) is limited by magnetic field fluctuations caused by trap current fluctuations and corresponds to a noise floor of {200, 97, 17}nm Hz −0.5 for displacements along the x, y, z directions, respectively. We connect the pickup loops in series to reduce the sensitivity to these field fluctuations. In the future, we will mitigate this noise by driving the trap using a persistent current [56], which should enable reaching the intrinsic noise floor of the SQUID (∼ 1 µϕ 0 Hz −0.5 ).
The magnetic trap and the SQUID magnetometer are thermally connected to the mixing stage of a dilution refrigerator. This allows the experiment to operate between temperatures as low as 50 mK and as high as the critical temperature of the superconducting particle (6.2 K for lead). Before levitating, the particle thermalizes on the bottom chip surface to the temperature of the chip substrate [see Fig. 1(b)]. To lift the particle off the chip surface, we ramp the trap current up to 0.8 A to overcome the adhesive force between the particle and the chip surface. With a current of 0.8 A, the lift force is ≈ 300 nN, which is usually sufficient to lift the particle. After that, the current is ramped down to the operating trap current.

III. CHARACTERIZATION OF CENTER-OF-MASS MOTION
Near the trap center, the particle experiences a harmonic trapping potential. The COM frequencies depend on the trap geometry, the particle's density, and the trap current. The penetration depth (∼ 40 nm for lead) is much smaller than the particle radius (24 µm) and so the particle can be modelled as an ideal diamagnet with magnetic susceptibility of −1 and its trapping frequencies ω i are given by [74] where ∇ i B is the magnetic field gradient along the i direction at the trap center, µ 0 is the vacuum permeability, I is the trap current, N is the number of trap coil windings, R is the trap coil inner radius, ρ is the particle's density and ζ is a geometric factor. At the center of an ideal anti-Helmholtz configuration 2ζ x = 2ζ y = ζ z = 0.86. In our system ζ x = 0.04, ζ y = 0.06 and ζ z = 0.12. (our coils are separated by 280 µm and have inner radii ≈ 125 µm as shown in Appendix S2). The trap axes are indicated in Fig. 1(a). Fig. 2(a) shows the power spectrum of a levitated microsphere. Throughout this work, unless otherwise stated, we levitated a 48 µm-diameter lead sphere using 0.5 A trap current and the cryostat temperature was 4 K. The peaks corresponding to the COM motion are colored. We identified these modes by comparing the peak frequencies with predictions from FEM simulations of our system [44,45]. We find good agreement between the measured and simulated COM frequencies. Peaks at the second harmonic of the COM frequencies are pronounced, particularly the second harmonic of the ∼40 Hz peak at ∼80 Hz. These peaks arise from the nonlinear pickup efficiency, rather than due to actual particle motion at these frequencies. We describe these peaks further in Appendix S3. In the future we will feedback cool the particle, then effects of the nonlinear pickup will be negligible.
The linear relation between the COM frequencies and the trap current given by Eq. (1) is shown in Fig. 2 We observed no significant difference between measuring this effect with the cryostat temperature of 50 mK or 4 K. We find good agreement between measurement results and FEM simulation results for the x and y modes. We attribute the 4% discrepancy for the z-mode frequency to simulation uncertainties in the FEM.
We confirm the inverse relation between trapping frequency and particle density of Eq. (1) by comparing the trapping frequencies of a lead particle and a tin-lead particle as we vary the trap current, see Fig. 2(c). The ratio of the frequencies is given by the square root of the ratio of the material densities. We can stably levitate the superconducting sphere for days in the chip-based magnetic trap. Fig. 3(a) shows the fluctuations of the COM frequencies of a levitated sphere over a 35 h measurement. We have yet to observe an upper limit to the levitation time, provided that the particle is not illuminated. When the particle is illuminated, as in Ref. [44], it heats up, loses superconductivity and falls on the bottom chip. When the particle is kept in the dark, as in this work, we have not yet observed any upper limit to the levitation time; we have measured up to 48 h.
Around once per day, we observe a sudden jump of all the COM frequencies of around 1 Hz (see Appendix S5). We attribute these jumps to changes in the magnitude or orientation of trapped flux in the particle. Fig. 3(a) shows 35 h of data in which we do not discern any such frequency jumps. The trapped flux will be the topic of a dedicated future investigation.
The particle's COM motion does not thermalize at the cryostat temperature, since it is strongly driven by the cryostat's vibrations. We estimate the mean amplitudes of the three modes to be 24, 10, and 7 µm for the x, y, and z modes, respectively, corresponding to an effective temperature of about 1 × 10 9 K. The energy in each mode (mω 2 i ⟨r 2 i ⟩) fluctuates on a time scale ∼ 80 s, as shown by the autocorrelation functions in Fig. 3(b). The mode energy fluctuations are due to the mode amplitude fluctuations, which occur over this time scale. By fitting the autocorrelation functions to exponential decay functions, we extract quality factors 3400, 4500, and 9300 for the x, y, and z modes, respectively [39]. We expect the quality factors to be limited by eddy current damping caused by normal-conducting metals in the vicinity of the levitated particle. In the future, we will mitigate this damping mechanism by surrounding the particle with a superconducting shield, with no normal-conducting metal within the shield.

IV. FREQUENCY PULLING
We can explain the COM frequency fluctuations of Fig. 3(a) as resulting from the fluctuating mode amplitudes [ Fig. 3(b)] together with frequency pulling. Frequency pulling describes the dependence of the COM frequencies on the mode amplitudes. It arises from quartic terms in the trapping potential of the form which cause the motional frequencies to be shifted depending on the motional amplitudes according to [75] ∆ω The γ ii terms in the potential are called Duffing nonlinearities, while the γ ij terms describe couplings between the modes. Experimental data showing frequency pulling is shown in Fig. 4. The spectral peak corresponding to the xmode (y-mode) shifts depending on the amplitude of the x-mode motion in Fig. 4 Fig. 4(c) and (d), in which the spectral peak frequencies ω x and ω y are plotted against the spectral peak area corresponding to the x-mode (the remaining graphs of the same form are included in Appendix S4). The slopes of these lines depend on the values of γ xx and γ yx as well as the pickup efficiency η x . We extract estimates of γ ij from FEM simulations of our system (see Appendix S1); this allows us to use the nine gradients (i.e. the nine linear relations between ω i and ⟨x 2 j ⟩, for i ∈ {x, y, z} and j ∈ {x, y, z}) to estimate the three efficiencies η i (the values are quoted earlier). The estimated efficiencies are used to convert the lower x-axes of Fig. 4(c)-(d) into the upper x-axes, and yield the secondary y-axis of Fig. 2(a).
The estimation of the three efficiencies is an overconstrained problem, it yields fairly consistent results for the nine graphs. For instance, the estimation of η x together with the FEM estimations of γ xx and γ yx describes well both the slope in Fig. 4(c) and Fig. 4(d).
To observe the frequency pulling effect presented in Fig. 4 we did not control the particle's motional amplitudes, instead, we filtered a long 35 h dataset in which the motional amplitudes randomly fluctuated: We separated the data into 10 s chunks and extracted the mode frequencies and mode areas from each chunk. To investigate the dependence of the mode frequencies on, e.g., the x-mode amplitude, we filtered out the data in which the y-and z-mode amplitudes were high. Each point in Fig. 4(c)-(d) corresponds to one 10 s chunk. To produce the seven spectra in Fig. 4(a) and (b) we binned the power spectrum of each 10 s-long dataset into seven bins, based on the x-peak areas, then we averaged the power spectra for each of the seven bins.
The mode amplitude fluctuations together with frequency pulling describe the frequency fluctuations of Fig. 3(a). We represent the frequency distributions by histograms in Fig. 4(e)-(g); we note they are asymmetric. We also plot histograms based on the frequency pulling model; the model histograms were constructed by predicting the mode frequencies for each 10 s interval of data based on the measured mode amplitudes, the nonlinear coefficients γ ij and the pickup efficiencies η i . The model histograms describe the observed frequency fluctuations well.
Because the modes are coupled via the trap anharmonicity [Eq. (3)] we expect the average energies of the modes to be similar. We can estimate the mean energies of the modes E i ≈ mω 2 i ⟨r 2 i ⟩ using the estimated pickup efficiencies η i , and we find E x = 2.1E y = 1.6E z . Because the estimated average energies rely on the estimated pickup efficiencies, the similarity of the average energies lends credence to our estimates of η i .

V. CONCLUSIONS
In conclusion, we have demonstrated an integrated superconducting chip device for magnetic levitation and detection of superconducting microparticles at mK temperatures. We have shown a good understanding of the trap potential, the COM motion, and inductive coupling by comparing trap models with measurements. Further, we have demonstrated that the COM frequencies can be tuned via the trap current and that the system is stable over days. For the large motional amplitudes in our experiments, we observe nonlinear behavior and mode coupling of the COM particle motion, in the form of amplitude-dependent frequency shifts.
Future experiments will employ a cryogenic passive vibration isolation system [71,72], which will decouple the particle motion from external mechanical vibrations. Such vibration isolation systems can attenuate our measured mechanical vibrations of about 1 × 10 −8 m Hz −0.5 by more than six orders of magnitude to achieve thermally driven COM motion at 50 mK. Further, the detection noise floor can be greatly reduced by the use of a persistent current trap [56] and improved magnetic shielding, which should allow us to reach the intrinsic noise floor of 1 µϕ 0 Hz −0.5 of our commercial SQUID. Alternatively, the particle motion can also be measured using flux-sensitive superconducting microwave resonators [60-62, 64, 65]. Furthermore, the pick-up efficiency can be increased by placing the pick-up coil closer to the particle, by using multi-winding coils, and by improving the flux-transfer efficiency to the SQUID. We estimate that these measures will yield an improvement in the pick-up efficiency of four orders of magnitude. This improvement is enabled by our understanding of the chip-based particle motion detection and the flexibility in the microfabrication of pick-up loops of our chip-based approach. Overall, the reduction of technical noise and increase in pick-up efficiency should enable feedback cooling to the quantum ground state [26,27,73] of the COM motion of microparticles, which would serve as a gateway to gen-erating quantum states of nano-to microgram masses [10,47,48].

S1. EXTRACTING THE TRAP ANHARMONICITY FROM FEM
We consider anharmonicities in the trapping potential up to quartic terms. The potential can then be described as The force acting on the particle is The FEM simulations give us the force acting on the particle due to the trapping potential at different displacements, as described in Ref. [45]. We fit the FEM results by Eq. (S2) with the coefficients of the trapping potential of Eq. (S1) as free parameters. In this way we extract the trapping potential coefficients from the FEM simulations. Note, in comparing γ ′ with γ from the main text: The potential energy which we calculate based on fitting the FEM simulations is represented in Fig. S1. We calculate the potential energy in 3D; in the figure 2D slices are shown. Fig. S1 also shows the relative contribution of the cubic and quartic terms to the potential.

S2. MEASUREMENT OF PARTICLE MOTION
As the particle displacement oscillates the magnetic flux threading the pickup loops oscillates in response. This oscillating flux of the pickup loops is transferred to the SQUID loop. The flux at the SQUID loop is converted into a voltage signal. In this section each of these steps are described in turn.

A. Inductive coupling to particle motion
We estimate the response of the change of flux threading the pickup loops in response to displacement of the particle along x, y and z by treating the trapping field as a magnetic quadrupole field (which is a good approximation near the center), using Eq. (5) from [74], and using knowledge of the geometry of our system (represented in Fig. S3). The results are shown in Fig. S4. The response to displacement along all directions is nonlinear. That is to say, the response has strong quadratic components ({−7.5, −6.0, −3.7} mϕ 0 µm −2 ) which we describe by the parameter u i in Section S3 B. We note that the response of the pickup to displacement along the z direction is largely linear, whereas the response to x and y displacements has a clear quadratic component. The response can be made more linear along the x and y directions by changing the geometry such that the particle equilibrium position is displaced further from the center of the pickup loops. Because the pickup loops are connected in series, the overall pickup-loop flux is the sum of the fluxes in each loop.
where L input ≈ 24 nH is the input coil inductance, L parasitic ≈ 33 nH is the parasitic inductance which is dominated by the inductance of the twisted wire pairs, L pickup ≈ 0.72 nH is the pickup inductance, and M input ≈ 0.87 nH is the mutual inductance between the SQUID loop and the input coil inductance. Using these values, we estimate the flux transfer efficiency η flux = 3.1×10 −2 . By improving the inductance matching, η flux can be improved to 0.5.
Based on the dependence of the pickup-loop flux on the particle displacement (Fig. S4) and the flux transfer efficiency from the pickup loops to the SQUID loop Eq. (S4) we make independent estimates of the pickup efficiency in the pick-up loop η pickup i,quad−field = {40, 15.3, 600} mϕ 0 µm −1 . These values are close to the ones found by using the measurement data and the model for the potential energy, where we obtain η pickup i,model = {50, 106, 625} mϕ 0 µm −1 . We note that we obtain the latter values by dividing the pickup efficiencies in terms of flux in the SQUID η i = {1.58, 3.3, 19.4} mϕ 0 µm −1 by the flux transfer efficiency η flux , i.e., η pickup i,model = η i /η flux . The small discrepancy may arise because our estimates of the response of the pickup loops to displacements along x and y using the data in Fig. S4 depend critically on our knowledge of the particle's equilibrium position, due to the nonlinearity of the pickup efficiency. Additionally, we made the simplifying assumption in Section S2 A of a magnetic quadrupole field at the trap center.

C. Flux in SQUID loop converted to SQUID voltage
We operate the SQUID in flux-locked loop mode. Then the SQUID output voltage V out is related to the SQUID loop flux ϕ SQUID by where M Finput =38 pH is the mutual inductance between the feedback electronics and the SQUID, and R F =10 kΩ describes the resistor used in the feedback electronics.

S3. HARMONICS IN SPECTRA
As well as observing the peaks in the spectra at frequencies ω x , ω y and ω z corresponding to COM motion along the x, y and z directions, we observe strong peaks at frequencies 2ω x , 2ω y and 2ω z (see the power spectrum in Fig. S6).
We observe that the peak areas of the second-harmonic peaks grow quadratically with the peak areas of the respective fundamental peaks, as shown in Fig. S2. The fundamental and harmonic peak areas A fi and A hi are related by The factor R obs i describes the observed quadratic relation between the peak areas, it is different for each mode. From the data in Fig. S2 we extract R obs = {1000, 870, 7.2} V −2 for x, y and z modes, respectively (in this section we work with the spectral peak areas in units of V 2 , that is why R has units V −2 ).
Constructing Fig. S2 involved breaking the 35 h-long dataset which was used in Fig. 3 and Fig. 4 of the main text into 10 s chunks, we then calculated the power spectrum for each 10 s chunk and extracted the peak areas for the three fundamental and three second harmonic peaks and compared them.
Below, we estimate that the actual motion of the particle at frequencies 2ω x , 2ω y and 2ω z arising from the trap anharmonicity is negligible. Instead, we believe the second harmonics arise primarily due to the pickup nonlinearity.
A. Particle motion at the second harmonic frequency arising from the trap anharmonicity In a potential of the form the particle will move with a component at frequency ω x and a component at frequency 2ω x . The amplitudes of these motions, at the fundamental frequency and at the harmonic frequency, are related by provided that the cubic term in the potential is a perturbation and provided that the motional amplitudes are small. We derive this equation by solving the equation of motion in this potential looking for a solution of the form in the regime when X 2 ≪ X 1 .
Because the peak areas are related to the mean-square amplitudes by ⟨x 2 i ⟩ = η i A i we expect the nonlinearity of the trapping potential to give And thus, using the values β iii extracted from FEM (see Section S1) we expect R trap = {6 × 10 −1 , 1.9 × 10 −4 , 2.1 × 10 −4 } V −2 for the x, y and z modes. These values are much smaller than the observed values R obs , and so the second harmonic peaks do not arise from motion of the particle at the second harmonic frequencies.
B. Second harmonic peaks caused by the pickup nonlinearity A particle in a 1D harmonic potential moves according to As described in Section S2 A and Fig. S4, the flux threading the pickup loops depends nonlinearly on the particle position x(t). This can be captured by the quadratic function (1 − cos 2ωt) + vX sin ωt + w (S13) Here v and u are conversion factors from particle displacement to pickup-loop flux; v describes the linear response of the pickup-loop, while u describes the quadratic response. w is an offset. The SQUID voltage signal is related to the flux threading the pickup loops by where we define the conversion factor k. And so, the voltage signal is related to the particle position by (1 − cos 2ωt) + qX sin ωt + r (S16) where q and p are conversion factors from particle displacement to SQUID voltage; q describes the linear response of the SQUID, while p describes the quadratic response, r is an offset. We see that  FIG. S6. The power spectrum displays peaks at the harmonics and at mixing frequencies, as well as peaks at the COM motional frequencies of the particle (the fundamental peaks are colored). The harmonics arise primarily from the nonlinearity of the pickup.
In the power spectrum of the voltage signal, the area of the fundamental peak is A f and the area of the harmonic peak is A h . And so Earlier we defined Overall, the nonlinearity of the pickup should give We extract u i and v i for the different directions x, y and z using the data in Fig. S4, and we extract η from the comparison of the observed frequency pulling (see Section S4) with the trap anharmonicities given by FEM (see Section S1). And so, we estimate R PU = {2.5 · 10 4 , 2.3 · 10 4 , 1.6 · 10 −3 } V −2 . These are order-ofmagnitude estimations, since the estimations of v i depend critically on our estimation of the position of the trap center. The estimates R PU are similar in magnitude to the observed values of R obs for the x and y modes, and so we believe the spectral peaks at the second harmonics arise from the nonlinear pickup.
S4. FREQUENCY PULLING Figures Fig. S7 and Fig. S8 employ the same dataset and analysis method as was used to construct Fig. 4(a-d) of the main text. While Fig. 4(a-d) shows the dependence of ω x and ω y on ⟨x 2 ⟩, Fig. S8 and Fig. S7 show the dependence of ω x , ω y and ω z on ⟨x 2 ⟩, ⟨y 2 ⟩ and ⟨z 2 ⟩.
In Fig. S7 the blue-colored straight lines represent linear fits to the data, they are meant as guides to the eye. The model curves are given by Eq. (3) in the main text. We extract estimates of the γ ij values from FEM (details in Section S1). In Fig. S8(a-i) the x-axes are the mean-square displacements in units of mϕ 2 0 , and thus the gradients of the nine model lines depend on the three pickup efficiencies η i [In (a-c) the gradients depend on η x , (d-f) depend on η y and (g-i) depend on η z ]. We treat the three pickup efficiencies as free parameters, and we obtain three estimates for each η i , from the data in each of the panels. We choose the values which provide the best agreement between the observations and the model. In this way we use the frequency pulling data to estimate the pickup efficiencies η i .

S5. FREQUENCY JUMPS
When we levitate a particle we occasionally observe sudden changes of all the COM frequencies at once. These changes appear to happen at random times. We expect these frequency jumps are due to reorientation of trapped flux in the particle, or due to a change of the amount of the trapped flux in the particle. In Fig. S9 three such frequency jumps are indicated by arrows. Here the particle was levitated for two days. The 35 h dataset used in Fig. 3 and Fig. 4 used the data between time 4 h and 39 h of this dataset. In each panel, lighter colored spectra have higher mode amplitudes. In (a-c) ⟨x 2 ⟩ is different for different spectra (both ⟨y 2 ⟩ and ⟨z 2 ⟩ were relatively low). Similarly, in (d-f) ⟨y 2 ⟩ is different for different spectra (both ⟨x 2 ⟩ and ⟨z 2 ⟩ were relatively low). Similarly, in (g-i) ⟨z 2 ⟩ is different for different spectra (both ⟨x 2 ⟩ and ⟨y 2 ⟩ were relatively low).