Nonlinear dynamics of diamagnetically levitating resonators

The ultimate isolation offered by levitation provides new opportunities for studying fundamental science and realizing ultra-sensitive floating sensors. Among different levitation schemes, diamagnetic levitation is attractive because it allows stable levitation at room temperature without a continuous power supply. While the dynamics of diamagnetically levitating objects in the linear regime are well studied, their nonlinear dynamics have received little attention. Here, we experimentally and theoretically study the nonlinear dynamic response of graphite resonators that levitate in permanent magnetic traps. By large amplitude actuation, we drive the resonators into nonlinear regime and measure their motion using laser Doppler interferometry. Unlike other magnetic levitation systems, here we observe a resonance frequency reduction with amplitude in a diamagnetic levitation system that we attribute to the softening effect of the magnetic force. We then analyze the asymmetric magnetic potential and construct a model that captures the experimental nonlinear dynamic behavior over a wide range of excitation forces. We also investigate the linearity of the damping forces on the levitating resonator, and show that although eddy current damping remains linear over a large range, gas damping opens a route for tuning nonlinear damping forces via the squeeze-film effect. Supplementary Information The online version contains supplementary material available at 10.1007/s11071-024-10018-x.

The charge model is a useful method to calculate the magnetic field distribution of permanent magnets.In this model, a magnet is taken as a distribution of equivalent 'magnetic charge', which is the source of magnetic field.For a current-free region and magnetostatic field, ∇ × H = 0 and ∇ • B = 0, where H, B are the magnetic field strength and magnetic flux density due to a magnetic source, respectively.Then, the irrotational magnetic vector field H can be written as the gradient of the magnetic scalar potential ϕ m : Because B = µ 0 (H + M), where µ 0 is the vacuum permeability and M is the magnetization, combining Eq. ( 1) and ∇ • B = 0 results in: For a free space without boundary conditions, we can solve Eq. ( 2) using Green's function G(x, x 0 ) for ∇ 2 and obtain a particular solution: where x is the observation point, x 0 is the source point, ∇ ′ operates on the primed coordinates, and the integration is over the volume for which the magnetization exists.Assuming M is confined in a volume V with boundary surface S and falls abruptly to zero outside of this volume, Eq. ( 4) can be written as: where n is the outward unit normal to S. Therefore, the volume charge densities ρ m and surface charge densities σ m can be defined as: For free space, B = µ 0 H, and subtituting Eq. ( 5) into Eq.( 1) obtains: For a permanent magnet as shown in Fig. S1, assuming that its magnetization is M = M r ẑ along z direction and its dimension can be denoted by two points P1(x 1 , y 1 , z 2 ) and P2(x 2 , y 2 , z 1 ), its magnetic field B can be calculated using the charge model.For permanent magnets, the volume charge density inside the magnets is zero ρ m = −∇ • M = 0, while the surface charge densities on top and bottom surfaces are: Based on Eq. ( 8), the magnetic field of the rectangular magnet is: Therefore, by integrating Eq. (10) with respect to x 0 and y 0 , we can obtain the x-component of the magnetic field: and the y-component of the magnetic field: Similarly, the z-component of the magnetic field can be obtained: the magnetic property of a permanent magnet is described by its remanent magnetic flux density B r , and B r = µ 0 M r .Therefore, with a known B r and dimensions, the magnetic field outside a permanent magnet can be calculated using Eq.(11), Eq. ( 12) and Eq.(13).For example, for a magnet with B r = 1.4 T and dimensions of 10 × 10 × 10 mm 3 (x 1 = −5, x 2 = 5, y 1 = −5, y 2 = 5, z 1 = −5, z 2 = 5), its magnetic field B x , B y , B z along a line (x = 2 mm, y = 3 mm, 5 mm < z < 10 mm) are calculated using Eq.(11-13) and shown in Fig. S2.For comparison, we also use a finite element method (COMSOL Multiphysics) to calculate the magnetic field of the same magnet and the results are also shown in Fig. S2, from which a good agreement between the two methods is observed.

C. Magnetic field of an array of magnets
For an array of multiple permanent magnets, the magnetic field in free space can also be determined using Eq.(11-13).For example, for a magnet array as shown in Fig. S3, it magnetic field can be calculated by: where B i (x, y, z) is the magnetic filed of the i−th magnet.

N S N S
FIG. S3.An array of four permanent magnets with alternating magnetic poles, where 'N' represents the north pole and 'S' the south pole, respectively.
After the magnetic field is determined, the magnetic force applied on the graphite plate can be calculated with: where V is the volume of the plate, H x,y,z are the components of the magnetic field vector, M is the magnetization vector and B the magnetic flux density vector.In this analysis it is assumed that the plate does not significantly affect the magnetic field, since its relative magnetic permeability is close to 1.All the parameters used in our calculations for pyrolytic graphite plates are listed in Table I.

S2. MAGNETIC FIELD CALCULATION BY FEM MODELING
In order to compare with the analytical modeling and take into account the fillets on the edges of the magnets, we also use FEM method to calculate the magnetic force.The geometry model of the magnets and graphite plate is shown in Fig. S4.The magnets are with dimensions of 12 × 12 × 12 mm 3 and with rounded edges with a fillet radius of 1 mm.The remanent magnetic flux density of the magnets is B r = 1.4 T. The material properties of the 10 × 10 × 0.28 mm 3 graphite plate are listed in Table I.Using COMSOL Multiphysics 5.6, we can simulate the magnetic field and then calculate the magnetic force using Eq. ( 16).

S3. CHARACTERIZATION OF SHAKER
To determine the relations between the output motion and driving voltage of the shaker, we shine the laser directly on the magnets and measure their displacement amplitude d with different driving voltages V ac .The measurement results for 4 different driving frequencies are shown in Fig. S5.It can be seen that the displacement of the shaker is approximately linear with the driving voltage and d = C V V ac , where C V = 0.0158mm/V.We note that the slope of the conversion factor C V increases from 0.0152 to 0.0178 when the frequency increases from 14Hz to 20Hz, which is attributed to the shaker having a slightly frequency-dependent base amplitude.

S4. SOLVING PROCEDURES FOR NONLINEAR EQUATIONS OF MOTION
The equation of motion was solved using a pseudo arc-length continuation technique implemented in the numerical software AUTO, with detailed documentation available in the manual [3] and on GitHub [4].Here, we only summarize the key solving procedures: (1) Transform the equation of motion (Eq.(4) in the main text) into two first-order differential equations.Input these equations along with all normalized physical parameters into the parameter setting file of AUTO.
(2) Set the control file of AUTO.The important settings for our system are: • NDIM=1: Indicates only one degree of freedom.
• IPS=2: Specifies the computation of periodic solutions.• DSMIN=0.000001 and DSMAX=0.005:Set the minimal and maximal pseudo-arclength step sizes for the first attempted step along any branch.This is important since our system has a bifurcation point.
• Initialize the sweeping frequency and define the sweeping frequency range.
(3) Retrieve the amplitude of the periodic solution at each step and plot the frequency response curve.
The stability analysis of the periodic solutions is done also via AUTO and by checking the Floquet multipliers associated with the periodic solution.More details on these can be found in our earlier work [5], as well as [6] S4

. MODELING NONLINEAR FREQUENCY RESPONSE IN AIR
Taking the nonlinear damping coefficient as a fit parameter, we can model the nonlinear frequency response curves using Eq. ( 5) in the main text, which match the experiment results well, as shown in Fig. S7.Note that to obtain a good fit, we need to adjust the liner resonance frequency for each fitting.The linear resonance frequencies used for the frequency response curve fitting are listed in Table II.
FIG.S2.Schematic of one magnet with magnetization M and its corresponding coordinate system.
FIG. S4.Geometry model of one graphite plate levitating over four permanent magnets.
FIG. S5.Displacement amplitude of the shaker as a function of driving voltages for four driving frequencies.
FIG. S6.Frequency response curves with different excitation voltages measured in air and their corresponding modeled curves using Eq.(5) in the main text.

TABLE I .
Material properties used for the simulations of the levitating pyrolytic graphite.

TABLE II .
Linear resonance frequencies used for the frequency response curve fitting Hz) 16.78 16.80 16.82 16.86 16.89 16.92