Modelling of Cavity Optomechanical Magnetometers

Cavity optomechanical magnetic field sensors, constructed by coupling a magnetostrictive material to a micro-toroidal optical cavity, act as ultra-sensitive room temperature magnetometers with tens of micrometre size and broad bandwidth, combined with a simple operating scheme. Here, we develop a general recipe for predicting the field sensitivity of these devices. Several geometries are analysed, with a highest predicted sensitivity of 180 pT/Hz at 28 μm resolution limited by thermal noise in good agreement with previous experimental observations. Furthermore, by adjusting the composition of the magnetostrictive material and its annealing process, a sensitivity as good as 20 pT/Hz may be possible at the same resolution. This method paves a way for future design of magnetostrictive material based optomechanical magnetometers, possibly allowing both scalar and vectorial magnetometers.


Introduction
Magnetometers with high spatial resolution are required for many applications such as magnetoencephalography [1], measurements of topological spin configurations [2] and nuclear magnetic resonance spectroscopy to identify chemical composition, molecular structure and dynamics [3]. Optical readout of magnetometers can offer high sensitivity for a given resolution, while being well decoupled from the magnetic signal. Among optical magnetometers, an ensemble of nitrogen-vacancy (NV) centres with a volume size of 8.5×10 5 µm 3 pushes the sensitivity down to 1 pT/ √ Hz [4]. However, NV magnetometry generally requires high optical power for excitation (e.g., 400 mW in Ref. [4]), as well as complicated microwave decoupling sequences in NMR spectroscopy, and is limited by the sample fabrication reproducibility [5]. A magnetometer based on micro-sized Bose-Einstein condensates has a volume of 90 µm 3 , but its quantum-enhanced sensitivity is limited to 1.86 nT/ √ Hz [6]. It is crucial yet challenging to reduce the size of magnetometers while maintaining competitive sensitivities.
Among various types of magnetometers, optomechanical magnetometers [7,8] reach sensitivities in the high pT/ √ Hz range at room temperature with sizes of tens of micrometres, comparable to the best cryogenic SQUID-magnetometer of the same size [9]. The principle of an optomechanical magnetometer is illustrated in Figure 1a. A magnetostrictive material converts the magnetic field to a force as a result of mechanical deformation.
The magnetostrictive response has a nonlinear component, a property that has been utilised in previous work to mix low frequency magnetic fields up to megahertz frequencies and therefore evade low frequency noise [8]. However, in general, it is far smaller than the linear component, so that the force may be well approximated by F field = c act B sig , where c act (N/T) is the actuation parameter and B sig (T) is the magnetic field to be measured. The amplitude of the mechanical response to this force is greatly enhanced when the magnetostrictive material is driven resonantly at its mechanical eigenfrequency by a modulated magnetic field. The mechanical response changes the path length of the optical cavity to which the magnetostrictive material is attached, allowing the magnetic field to be read out optically from the shift of the optical resonance [10]. While significant successes have been achieved in experimental demonstrations of optomechanical magnetometers [7,8], modelling and sensitivity-prediction for these devices have been somewhat ad hoc [11,12]. Better modelling techniques are needed to both enhance understanding of previous experimental results and for design of future magnetometers.
In this work, we present a model of magnetostrictive magnetometers that accounts for arbitrary mechanical mode shape and device geometry. We modify the elastic wave equation, which describes the small-amplitude motion of elastic materials, by including magnetostrictive stress. This modified elastic wave equation is then numerically solved by finite element analysis (using COMSOL Multiphysics). Magnetomechanical overlap, describing the overlap between the magnetostrictive deformation induced by the signal magnetic field and the excited mechanical eigenmode, is intrinsically included in the matrix form of the modified elastic wave equation, with each matrix element containing directional information. Mechanical properties are extracted from the solution to the modified elastic wave equation from COMSOL to be further combined with optomechanical analysis [10] to predict the sensitivity of a magnetometer for a given geometry. The coupling of magnetostrictive material to an optical cavity is quantified by the effective cooperativity C eff . The magnetostrictive material converts a magnetic field to a force F field = c act B sig with B sig being an oscillating magnetic field. Thermal force and optical shot noise act as noise terms.
κ (rad·s −1 ), Γ (rad·s −1 ), ω 0 (rad·s −1 ), and Ω M (rad·s −1 ) are optical and mechanical decay rate, optical and mechanical resonance frequency, respectively; (b) sketch of a magnetometer with micro-toroidal structure coupled to a tapered optical fibre; (c) homodyne detection scheme. The signal arm couples a coherent light source in and out from a magnetometer via a tapered optical fibre through an evanescent optical field, and is mixed with a strong reference beam (local oscillator field) by a 3 dB coupler.
The magnetometer is embedded in the signal magnetic field.
We apply this analysis to study the effect of the position of the magnetostrictive material on the sensitivity of devices similar to those reported in Ref. [7]. Using the piezomagnetic constant measured from a rod of the magnetostrictive material Terfenol-D [13], we model a magnetometer design, where the Terfenol-D is deposited directly on top of a standard silica toroid. From there, we employ single mode analysis (Appendix A) and discover that a bimetallic-stripe-like bending effect, similar to the bimetallic bending effect in a cantilever [14], greatly enhances the sensitivity when the magnetostrictive material is positioned off-centre. Optimisation of this effect may allow substantial improvements in sensitivity in future devices. Furthermore, we investigate the sensitivity achievable from a device comprised of a toroidal structure with a centre hole that is filled with the Terfenol-D, as studied experimentally in Ref. [8] and sketched in Figure 1b . We predict a peak sensitivity of 180 pT/ √ Hz over a broad spectrum by using multi-mode analysis under optimised operational conditions, in good agreement with current experimental observations. This numerical model allows specification of the orientation of a sample to maximally enhance the magnetomechanical overlap, thus amplifying the detected magnetic field signal, as well as characterization of the magnetomechanical overlap in response to the variation of the magnetic field direction. This is crucial to vectorial magnetometers that measure not only the intensity but also the direction of the magnetic field.

Concept of Optomechanical Magnetometry
Optomechanical magnetometry can be schematically explained via the example of a Fabry-Pérot optical resonator coupled to a spring-mass mechanical oscillator as depicted in Figure 1a. An applied magnetic field B sig causes a deformation to a magnetostrictive material attached to the mechanical oscillator (see Appendix B for details of how this field is generated in COMSOL). This induces a F field = c act B sig on one movable end mirror of the optical resonator, changing the optical path length and thus the optical resonance frequency. The shift in the optical resonance frequency is therefore proportional to the applied magnetic field. The transduction from magnetic field to mechanical motion is determined by the actuation parameter c act depending on magnetomechanical overlap and magnetostrictive coefficient. The magnetic field signal encoded on the motion of the mechanical element is read out by optically probing the the optical resonance frequency. This can be achieved with high precision by coupling a coherent optical field into the cavity, collecting the output field, and measuring the change in its amplitude or phase due to the modulation of the optical resonance frequency. For instance, directly detecting the output field, as in several reported experiments [7,8], measures changes to the amplitude of the output optical field and enables simple operation. Alternatively, a homodyne scheme can be used, allowing an arbitrary quadrature of the optical field to be accessed as shown in Figure 1c. Here, the output field is interfered with a bright local oscillator field prior to detection. The transduction from mechanical displacement to optical signal can be quantified by the effective cooperativity C eff [10].
The magnetic field sensitivity is limited by noise consisting primarily of thermal force and shot noise on the optical field. Thermal noise is explained by the equipartition theorem, which states that each mechanical degree of freedom of an object has a mean energy of k B T/2 (k B is the Boltzmann constant and T is the temperature). This energy excites incoherent mechanical vibration near mechanical eigenfrequencies. The bandwidth of the magnetometer depends on the visibility of the thermal noise over the optical shot noise. For the case of a single mechanical resonance, the sensitivity is flat over the frequency range where thermal noise dominates shot noise, and degrades outside of this region. Consequently, in this case, the bandwidth is given simply by the thermal-noise-dominant frequency band, which is typically on the order of a few megahertz [15]. The case of multiple mechanical modes is more complex due to variations in actuation constants, effective cooperativities and mechanical parameters, and due to interferences in the coherent response of the mechanical modes.
In this paper, as a test geometry for our model, we choose optomechanical magnetometers of the form reported in Refs. [7,8]. They utilise a silica microtoroid as the optical resonator. The magnetostrictive material is embedded in or deposited onto the microtoroid as sketched in Figures 1b and 2a, respectively. Combined, the silica microtoroid, the magnetostrictive material and the silicon pedestal serve as the mechanical oscillator. Using a tapered optical fibre placed next to the toroid, the optical field can be coupled in and out of the microtoroid through an evanescent optical field. This optomechanical magnetometry platform offers a simple operational scheme and low energy consumption with state-of-the-art field sensitivity for a micro-magnetometer.

Numerical Methods
The primary objective of this work is to develop a versatile technique to numerically obtain a meaningful estimation of the magnetic field sensitivity for a wide range of sensor geometries. We consider the case of phase quadrature detection in a homodyne scheme and on-resonance optical probing of the cavity resonance, which maximises the signal-to-noise. We note, however, that simpler direct detection with off-resonance probing and an optimal detuning of κ , where κ is the optical cavity linewidth, only degrades the sensitivity by a factor of 8 3 3/2 ∼ 1.5. The sensitivity as a function of magnetic field frequency Ω can be determined from the finite-time sensor power spectrum S(Ω), which can be separated into a stochastic noise term S noise (Ω) and coherent signal term S signal (Ω) as where i is the photocurrent, normalised so that the optical shot noise contribution to S noise (Ω) is equal to 1/2 [10], and τ is the measurement time. At frequencies Ω 2π/τ and considering j mechanical modes, S noise (Ω) is given by [10] S noise (Ω) = 1 2 where the first term is the optical shot noise and the second term constitutes the combination of mechanical thermal noise and quantum back-action noise. The detection efficiency η, consisting of the loss in the fibre-device coupling and detection process, is ideally taken to be 1 in the model. However, in the non-back-action dominated regime relevant here, reductions in efficiency can be exactly modelled by a proportionate decrease in the optomechanical cooperativity. Γ j is the mechanical decay rate of mode j and C eff,j is its effective cooperativity, which depends on the input laser power used, the decay rate of the optical field and mechanical excitation, and the radiation pressure coupling rate between them. The mechanical susceptibility of mode j is defined as χ j (Ω) ≡ Ω M,j /(−Ω 2 − iΩΓ j + Ω 2 M,j ), with Ω M,j its mechanical resonance frequency. k B T/hΩ M,j is the number of phonons thermally excited at room temperature, withh being the reduced Planck constant. The mechanical motion induced by an alternating-current (AC) magnetic field is quantified by the finite-time power spectrum S signal (Ω). This is calculated by replacing the thermal environment forcing F th in the input momentum fluctuation P in = x zpf F th /h √ Γ [10], which leads to Equation (2), with a coherent sinusoidal driving force F field (t) = c act B(t) sig at frequency Ω, and neglecting the incoherent noise terms (laser shot noise in amplitude and phase quadrature). This results in the expression where B signal,rms is the root-mean-square amplitude of B sig (t), m eff,j is the effective mass of mode j, and c act,j is the actuation constant associated with that mode. This finite-time power spectrum takes into account mechanical interference, as experimentally observed, for example in optoelectromechanical systems coherently driven by an electric field [16].
The frequency dependent signal-to-noise ratio (SNR) of the magnetic field measurement is given simply by The minimum detectable field in the measurement time τ is defined as the field that produces a signal-to-noise ratio SNR of one, i.e., B min,τ = B sig (SNR = 1). It should be noted that the stochastic noise power spectral density S(Ω) noise of Equation (2) is independent of integration time, whereas the integral of a coherent band-limited signal power spectrum, as described by S signal (Ω) in Equation (3), increases linearly with time. Consequently, B min,τ improves with measurement time as τ −1/2 . To obtain a minimum detectable field in the conventional units of Tesla per root Hertz, independent of time, we multiply through by τ 1/2 with the result To determine the minimum detectable field via finite element simulations, we use COMSOL Multiphysics. Simulations detailed in the appendices allow us to extract each of the parameters in Equations (2) and (3) and therefore predict the sensitivity. These simulations involve both mechanical eigenmode solving to determine the resonance frequency, effective mass and effective cooperativity of each mechanical mode of a given device geometry; and magnetic field driving to determine the coherent response of the mechanical modes to a magnetic field and the interferences between them. The approach is briefly sketched in what follows.
The spatio-temporal mechanical modeshape is described by a separable function u(r, t) = Ψ(r)x(t). The effective mass m eff,j for one mechanical resonance at an eigenfrequency Ω M,j is calculated from the maximum physical displacement max r [|Ψ(r)|] as m eff,j = V ρ n |Ψ(r)| 2 dV [17], with normalization max r,t [|u(r, t)|] = max t [x(t)] and therefore max r [|Ψ(r)|] 2 = 1. ρ n is the density of the material and the subscript n denotes different parts of the device (for instance, silica for the optical resonator and Terfenol-D for the transducing medium). Note that, while this definition of effective mass is the convention for microelectromechanical systems, an alternative definition-where the effective mass is defined with respect to the optical path length-is commonly used in the optomechanical community [18]. This choice of convention has no effect on the ultimate predictions of our model.
The magnetic field response S signal (Ω) of the sensor is determined by the eigenmode-dependent actuation parameter c act . For a single mechanical eigenmode, the equation of motion is At the resonance frequency of each mechanical eigenmode, c act can be extracted as a fitting parameter in the mechanical signal frequency response spectrum obtained from COMSOL. Taking the Fourier transform of Equation (6), we see that This allows c act to be determined for each mechanical mode. Due to the magnetostrictive energy stored within compressed magnetostrictive materials, the extraction of m eff and displacement from COMSOL for such materials requires modification of the elastic wave equation. To treat the magnetostrictive material in COMSOL, we built upon a previously used method [19][20][21], including the magnetic field in a driving stress σ driv and adding a damping stress σ Γ to the elastic stress σ ela , which describes the mechanical properties without driving force in the elastic wave equation [22], resulting in The modulated driving stress is linked to the magnetic field via the piezomagnetic constant [20], and a low value for the damping stress σ Γ is chosen manually to avoid an artefactual infinity in the mechanical displacement at resonance (see Appendix C for technical details). Simulations reveal that the influence of a particular value chosen for σ Γ on numerical results is negligible (Appendix D).
To obtain the value of effective cooperativity, we quantify the effectiveness of transduction of mechanical motion to measurable optical path length change as the geometrical factor, as where δL is the change of the optical path length due to the mechanical displacement. The extraction of the value of ξ from COMSOL is detailed in Appendix E. Within one mechanical mode, ξ is directly linked to the effective cooperativity (Appendix F) by where κ is the optical decay rate, ω 0 is the optical resonance frequency, N in (photons·s −1 ) is the input optical photon number flux, L is the optical path length, and η esc is the escape efficiency counting fibre-device coupling. The front part of the right hand side of Equation (10) is arranged to be mechanical mode dependent. The calculation of the magnetic field sensitivity from Equations (2), (3) and (5) can then be obtained based on the value of the geometrical factor ξ.

Bending Effect
To verify the numerical model, we apply it to the first experimentally realized optomechanical magnetometer [7]. For simplicity, we begin the analysis considering only a single mechanical eigenmode (Appendix A). The magnetometer as sketched in Figure 2a consists of a silica micro-toroidal cavity with major radius of 33 µm. The Terfenol-D is glued on top of the silica and is modelled as a semi-sphere with a transverse radius of 18.5 µm and a height of 15 µm. The optical quality factor Q o = ω o /κ is taken to be 2×10 7 from the experiment. The mechanical quality factor Q M = Ω M /Γ is assumed to be 200 for all modes which is a simplification, but is roughly in line with the experimentally observed quality factors. A continuous input laser is locked to the optical cavity resonance in the homodyne detection scheme, and the input laser power ensures that on mechanical resonances thermal noise dominates over optical shot noise.
From available optical microscopic images, it is not clear whether the Terfenol-D is centred on the toroid or not. Therefore, we sweep the position of the Terfenol-D from the centre. Without loss of generality, we analyse the magnetic response for a second order crown mode because this mode has been commonly observed in experiments [15,16,23]. For the magnetometer with centred Terfenol-D, the effective motional mass is m eff = 3.9 pg with eigenfrequency at 10.1 MHz. As the Terfenol-D is moved away from the centre as illustrated in Figure 2a, the mechanical eigenmode changes (Figure 2b). Generally, the top of the Terfenol-D stretches more than the bottom part attached to a silica disk during a mechanical oscillation. This is also the case for silica where the top layer experiences the force from the Terfenol-D and the bottom layer is clamped to the silicon pedestal. Therefore, a bimetallic-like strain gradient is formed vertically. In the second order crown mode, as the major motion takes place at the silica layer instead of the Terfenol-D, the strain gradient can be viewed inside the silica disk at the edge of bottom Terfenol-D and top facet of silicon. With the centred Terfenol-D, the strain at the top layer of the silica is nearly two orders of magnitude smaller than that with 4 µm Terfenol-D position offset as shown in the red areas in Figure 2c. This local maximum strain leads to the maximum displacement of the device (pointed by the arrows on tori in Figure 2b) in the radial direction. Figure 2d shows the best sensitivity of 78 nT/ √ Hz, when driven by an in-plane magnetic field, takes place when the Terfenol-D offset is at 4 µm, nearly two orders of magnitude better than that of Terfenol-D centred (the same order of magnitude difference as that of the strain). We therefore see that the position of the Terfenol-D on the silica layer has strong influence on the bimetallic-like strain effect, and consequently the sensitivity. Moreover, the effective cooperativity |C eff | of the crown mode experiences four orders of magnitude enhancement with only a few micrometres Terfenol-D offset as plotted in Figure 2d. |C eff | is chosen for evaluating mechanical mode shape induced characteristics. Terfenol-D with offset breaks the axial symmetry of the crown mode, creating a first order circumference difference of the toroid as the mechanical mode oscillates, and thus improves the value of |C eff |.
The numerical results show that asymmetry and the bimetallic-like bending effect helps to enhance the sensitivity. With an optimal offset of Terfenol-D, low nT/ √ Hz sensitivity is predicted, which is five times better than the experimental result [7]. It is likely that the experimental results were degraded not only due to a lack of Terfenol-D offset, but also by the epoxy used to fix the Terfenol-D on top of the toroid, reducing the expansion of the silica disk.

Effect of the Size of the Terfenol-D
The single mechanical mode analysis is then applied to a proposed [7] thin disk structure: 1 µm sputter coated Terfenol-D film on top of a 400 nm-thick silica disk. Magnetometers with sputter coated Terfenol-D have the advantage of a reproducible fabrication process. The silica disk has a radius of 30 µm and the pedestal has a top facet of 15 µm (sketched in Figure 3a inset top). The optical quality factor is taken to be 1 × 10 6 [24], a coherent laser source is again used to probe the system with zero detuning and measured via homodyne detection. The effective mass extracted from numerical simulation varies from 1 pg to 3.8 pg depending on the Terfenol-D size, for the radial breathing modes of the device.  Figure 3a shows the relation of sensitivity to the size of the Terfenol-D for the first order radial breathing-type mode. The signal magnetic field drives the radial breathing mode in the axial direction to create a magnetic field induced deformation profile as shown in Figure 3b. Unlike an isotropic magnetostrictive material breathing radially under axial magnetic field driving, the spatial profile from the non-isotropic Terfenol-D stretches only in one direction. When the size of the Terfenol-D is larger than the top facet of the silicon pedestal, the silica disk is also significantly affected by the motion from the Terfenol-D. The part of Terfenol-D inside the top pedestal facet (Terfenol-D is highlighted with the white dashed line in the mechanical eigenmode simulation in Figure 3a inset) is motionless because it is obstructed by the silicon pedestal. When the rim of the Terfenol-D reaches outside the top pedestal facet, the device mechanical motion is hybridized with mechanical modes of the Terfenol-D. This leads to a bi-metallic-strip-like effect close to the edge of the top facet of the silicon pedestal across the silica layer, increasing the silica displacement and thus allowing for better sensitivity than in the cases where the Terfenol-D is confined inside the silicon pedestal. Generally, the sensitivity scales with the size of the motional part of the Terfenol-D. A sensitivity of 2.9 nT/ √ Hz is predicted when the diametre of the Terfenol-D disk covers more than 2/3 of the silica disk in Figure 3a. A power-law fit (y(x) = a · x b with fitting results of a = 7.9 × 10 −7 and b = −1.7) is applied to the data with the Terfenol-D radius larger than that of the pedestal, predicting a 300 µm radius of Terfenol-D may lead to 50 pT/ √ Hz sensitivity. To achieve better sensitivity, the size of the Terfenol-D must be larger than the pedestal so as to have large portion of motional Terfenol-D and large bi-metallic-strip-like bending effect, which could be realised by decreasing the size of the silicon pedestal and by increasing the size of the Terfenol-D.

Multi-Mode Analysis
Single mode analysis is limited, in that it only correctly predicts the performance of devices over frequency ranges where only one mechanical mode contributes significantly to the dynamics. In reality, this is rarely the case, and often there is a dense spectrum of mechanical modes (see e.g., Ref. [7,8]). To extend our analysis to such situations, we use multi-mode analysis from Section 3. We first examine the limitations of the single mode analysis and then predict an optimal driving direction of the magnetic field leading to a best predicted sensitivity of an ensemble of mechanical eigenmodes.
We examine the limitations of single mode analysis by considering the magnetometer design reported in Ref. [8]. This type of magnetometer has a hole of 14 µm radius in the middle of a silica toroid, which has a 45 µm major radius. A cross-sectional view is shown in Figure 4a, where the outer silicon undercut is 15 µm. The Terfenol-D is modelled as an ellipsoid having the same transverse radius as the silica hole and an axial radius of 16 µm. Mechanical modes with resonant frequencies up to 45 MHz are selectively driven with the in-plane B sig in accordance with the experimental conditions of Ref. [8]. Three windows (∼7 MHz, ∼26 MHz and ∼43 MHz) of interest are selected. Mechanical modes in between are not taken into consideration due to their small optomechanical coupling resulting from their symmetrical mode shapes. The power spectral density S noise (Ω) and magnetic field sensitivity spectrum in Figure 4b are obtained, again choosing Q M = 200 for all modes, and setting Q o = 2 × 10 6 and a coherent laser with power of 1 µW at 1550 nm in an on-resonance homodyne detection scheme. With these parameters, the sensor noise floor is dominated by mechanical thermal noise close to the mechanical resonance frequencies, and optical shot noise at other frequencies (Figure 4b top). A single mechanical mode at Ω M /2π=23 MHz has the largest actuation parameter (see Appendix D for c act spectrum) due to a relatively large spatial mode overlap between the mechanical eigenmode (Figure 4b inset) and the magnetic field induced deformation profile (Figure 4c) compared with other modes. However, this particular mode has a very weak optomechanical coupling when the device is modelled uniformly and axial-symmetrically. This prevents the mode from being optically resolved from the thermal noise of others, causing a large difference of the magnetic field sensitivity between the single mode and multi-mode analysis, as shown in triangles and lines in Figure 4b bottom, respectively.   [8]. (a) cross-sectional view of the optomechanical magnetometer; (b) top: the power spectral density S noise (Ω) (blue) is the sum of individual thermal Brownian motion peaks (grey) and coherent laser shot noise on the optical phase quadrature (red); bottom: minimum detectable magnetic field from multi-mode (blue) and single mode (black triangles) analysis driven by in-plane magnetic field. The inset is the mechanical mode with the highest c act at Ω M /2π = 23 MHz; (c) deformation profile induced by in-plane magnetic field far away from mechanical resonance frequencies; (d) top: the power spectral density S noise (Ω) of the radial-breathing-like mechanical modes. The insets show the mechanical eigenmodes corresponding to each resolved thermal Brownian motion peaks; middle: the magnetic field response S signal (Ω)/τ to the axial magnetic field driving; bottom: the sensitivity spectrum from multi-mode (blue) and single mode (black triangles) analysis driven by the axial magnetic field; (e) deformation profile induced by axial magnetic field far away from mechanical resonance frequencies.
To achieve better sensitivity, the direction of the driving magnetic field needs to be optimised. The mechanical mode under magnetic field driving should have both relative large optomechanical coupling and relative good magnetomechanical overlap compared to other modes. As might be expected, and is shown in Figure 4b, top modes with radial-breathing-like motion (Figure 4d top insets show the eigenmodes) offer the largest optomechanical coupling. These mechanical modes are at 4.8 MHz, 26 MHz, 27 MHz, 43.2 MHz and 43.4 MHz. When driven axially, the deformation profile due to magnetostriction is also radial, as shown in Figure 4e. This suggests the magnetometer will perform well when axially driven near radial breathing modes. Choosing axial field magnetic field driving, we find the power spectral density, network response and sensitivity shown in Figure 4d. The radial breathing mode at Ω M /2π = 27 MHz, third from left in Figure 4d top inset, reaches a sensitivity of 180 pT/ √ Hz. We confirm that the result from multi-mode analysis (see Figure 4d bottom blue line) is consistent with single mode analysis (see Figure 4d bottom triangles) for this mode. The actuation parameter is 3200 times larger than if the same mechanical mode is driven by an in-plane magnetic field (see Appendix D for c act values), verifying a strong dependence of the magnetomechanical overlap on the magnetic field direction and the potential for vectorial magnetometry.
With in-plane magnetic field driving, the sensitivity observed in the experiment 200 pT/ √ Hz [8] surpasses the modelled sensitivity by around two orders of magnitude. This is likely due to the fact that the simulated mode at 23 MHz (Figure 4b bottom inset) is thermally resolved in the experiment, which is not the case in the model. This difference can be understood in terms of symmetry. In the model, the symmetry results in a very poor predicted optical transduction sensitivity. However, in the experiment, it can be expected that the symmetry is broken due to fabrication defects resulting in improved sensitivity [25].

Conclusions
We have developed a new versatile approach to model the sensitivity of optomechanical magnetometers, introducing magnetostriction into the elastic wave equation used to solve for mechanical eigenmodes. By numerically solving a modified elastic wave equation for a range of geometries, we model the sensitivity for magnetometers both experimentally demonstrated and not-yet fabricated. The modelling predicts that at least one order of magnitude improvement from previous experimental results [8] is possible. The sensitivity of optomechanical magnetometers can be significantly improved by optimising the size and the shape of the Terfenol-D, by utilizing the bending effect, which arises from a magnetic equivalent of the bi-metallic strip effect, and by optimizations of the composition and the annealing process of Terfenol-D, which may lead to sensitivity below 20 pT/ √ Hz using the piezomagnetic constant in Ref. [26] with micrometre-level resolution. The numerical method developed here is applicable to optomechanical magnetometers with a wide range of geometries and any magnetostrictive materials. A full characterization of the response of the magnetomechanical overlap to the variation of the signal magnetic fields direction may allow vectorial optomechanical magnetometry, complementary to vectorial optomechanical force sensors [27,28]. Micro-optomechanical magnetometers with pT/ √ Hz sensitivity can potentially be applied to detect signals from neurons, similar to recent results with nitrogen-vacancy centre based magnetometers [29] and atomic magnetometers [30], but with benefits of a simpler, silicon-chip fabricateable approach, as well as high bandwidth.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Derivation of the Sensitivity for a Single Mechanical Mode
For a single mechanical mode, the minimum detectable magnetic field can be obtained from the actuation parameter c act and the noise force spectral density S FF . Calibrated in the medium of air, the sensitivity at individual mechanical eigenfrequencies can be written as where a factor of 1/ √ 2π ensures B min having the unit of T/ √ Hz. Noise sources considered are thermal noise and noise from optical measurement including imprecision and back-action. Measurement imprecision comes from the laser shot noise in the optical phase quadrature. Back-action noise is due to the laser shot noise in the optical amplitude quadrature driving the mechanical oscillator.
The power spectral density in the unit of force (specifically rad·s·N 2 ) for individual mechanical modes driven by noise can be found in Ref. [10]. Here, we extend the calculation to include the back-action noise, as: in which a factor of 4 in front of the classical thermal force spectrum in Equation (A2) is due to the definition of Γ being the full-width-half-maximum of the mechanical oscillator. Inspection of Equations (A3) and (A4) shows that, despite optomechanical coupling, the effective cooperativity C eff also quantifies the trade-off between better measurement precision and large back-action noise due to the Heisenberg uncertainty relation. |C eff (Ω)| is given by where N is the intra-cavity photon number, N in (photons ·s −1 ) is the input photon number flux, and g 0 (rad·s −1 ) is the vacuum optomechanical coupling rate quantifying the optical resonance frequency shift by the mechanical displacement at zero energy excitation. Fibre-device coupling here is idealized to be lossless where the intra-cavity and end mirror loss due to the scattering and/or absorption of the light is neglected, leaving the optical decay only counted at the front mirror to be κ as shown in Figure 1a and thus making the cavity escape efficiency η esc = κ/(κ + 0) = 1.
To provide an idea of the relative magnitude of S therm FF , S imp FF and S ba FF , we choose the geometry and parameters (in Section 5 in the main text) of the magnetometer reported in Ref. [8]. Not surprisingly, on mechanical resonances at room temperature, we find that back-action noise is always smaller than the thermal noise by a large margin. For each mechanical mode, this is quantified, roughly, by the ratio of thermal phonon occupancy to effective cooperativity, with the former being in the range of 10 5 -10 6 for our mechanical frequencies at room temperature, and the later not exceeding 100 for typical parameters. Far from mechanical resonance, the backaction noise will eventually exceed the thermal noise [10]. However, in this regime, the optical shot noise dominates. As a consequence, backaction noise can be safely neglected at all frequencies.

Appendix B. COMSOL Implementation of the Magnetic Field
The magnetic field B sig is generated by a pair of Helmholtz coils whose axis can be freely rotated in a 4π solid angle as shown in Figure A1a for COMSOL layout. To enable the simulation of the magnetic field, the outermost sphere is filled with air. The amplitude of the magnetic field is controlled by inputting a known current in the pair of Helmholtz coils. The coils diametre is set to be more than 40 times larger than the lateral size of the Terfenol-D to ensure a uniform driving magnetic field. Therefore, the direction of B sig is along the axial axis of the pair of coils. The magnetic field can be viewed by the intersected orthonormal slices on which the magnetic field amplitude is projected, with the colour refers to the amplitude of the magnetic field as shown in Figure A1b. At high frequencies, the eddy current induced magnetic field opposes the external magnetic field and thus reinforces most of the magnetic field between the surface and skin depth, leaving most inner part of the magnetostrictive material unused. This undesired effect is evaluated by varying the frequency of the signal magnetic field B sig . The magnetometer [7] as illustrated in Figure 2a left is used to perform the eddy current simulation. Figure A1c shows that the skin depth effect starts to take place at frequencies above 1 GHz. The in-plane magnetic field B sig characterized at the location of Terfenol-D is 7.7 µT, and the colourmap shows the magnetic field inside the Terfenol-D is ∼25 µT at low frequencies as shown in Figure A1c left. The value from the colourmap is consistent with the calculation from B sig and the value of the relative permeability tensor in Table A1. The simulation agrees with the simple relation for skin depth δ skin = 1/ π f µ 0 µ r σ c , where the conductivity for Terfenol-D is σ c = 1.67×10 6 S/m. Mechanical eigenfrequencies of interest of modelled devices are below 50 MHz, far detuned from the range in which eddy currents pose a problem.

Appendix C. COMSOL Implementation of the Modified Elastic Wave Equation
COMSOL's Solid Mechanics module allows users to modify the elastic wave equation and solve it numerically. In the modified elastic wave equation given by Equation (8) displacement u = ∑ i u i (r, Ω)e i , and stress tensor σ = ∑ ij σ ij e ij = ∑ n σ n e n in three dimensions can be fitted into a 3 × 3 matrix as Up to the first term in the bracket on the right-hand side (RHS) of Equation (8) is the elastic stress σ ela . Elastic stress is connected to strain via a tensor coefficient λ ijkl as σ ij ela = λ ijkl kl . The tensor with λ ijkl being its elements is termed elasticity matrix in COMSOL. For isotropic materials, elements of the elasticity matrix are determined by an isotropic Young's modulus and an isotropic Poisson's ratio, while for anisotropic materials the number of independent elements can go up to 21 in the 6 × 6 matrix [22]. Terms σ ela + σ driv in the bracket on the RHS of Equation (8) incorporate both the elastic stress and the stress under the magnetic field driving. Assuming (1) the variation of the AC magnetic field is slow enough for the material to reach deformed equilibrium, (2) the magnetostrictive material exhibits reversibility, and (3) the operational point is far below the magnetostrictive saturated strain defined as the ratio of the maximum material elongation to its original length, the stress-magnetic field relation can be linearly approximated [19]. σ ela + σ driv with small modulation can be expressed via first order Taylor expansion, when projected onto one dimension, as where the modulated Maxwell stress tensor ∆σ ij Maxw [31] describes the stress caused by the interaction of a magnetized ferromagnetic material and the external magnetic field, and its value depends on the relative permeability of the magnetostrictive material. H k is the magnetic field strength inside the Terfenol-D when the external magnetic field strength H 0 = B sig /µ 0 (µ 0 is the vacuum permeability) magnetizes the magnetostrictive material due to the effect of the demagnetization field. The internal magnetic field H k decreases as the length of the magnetic rod is reduced from infinite length where H k = H 0 to zero-thickness where H k = H 0 /µ r , in which µ r is the relative permeability of the magnetostrictive material.
The last term in the bracket on the RHS of Equation (8) is the input damping σ Γ . It is chosen to be proportional to the time derivative of strain as σ Γ = η damp˙ (the unit of η damp is Pa·s) so that the integral form of the elastic wave equation along one dimension followed by a Fourier transformation to the frequency domain can have the damping coefficient Γ in front of the term linear with Ω in the mechanical susceptibility. Only the part of (σ ij ela + σ ij driv + σ ij Γ ) perpendicular to dS contributes to the integral. Roughly speaking (simplified to one dimension), m eff Ω 2 M is related to elasticity matrix and spatial profile of the mechanical eigenmode Ψ(r) and surface area S. The damping factor Γ includes the input damping η damp (Pa·s), m eff , Ψ(r) and surface area S. The actuation parameter is determined by c act = e S ⊥ /µ 0 µ r (S ⊥ being the area perpendicular to the magnetic field induced stress). Note that Ψ(r) comes from an ansatz of Equation (8) as the displacement u(r, Ω) = Ψ(r)x(Ω) can be decoupled into a spatial and a frequency dependent term.
Due to the transverse (axial) symmetry of the Terfenol-D, the independent elements of the elasticity matrix are reduced to only six [19] where small modulation ∆ kl and ∆H k in Equation (A8) is replaced with kl and H k in the frequency domain. The input elasticity matrix elements λ H ijkl and piezomagnetic constant elements e ijk are taken from an experimental measurement biased at 60 kA/m and prestressed at 20 MPa [13] as summarised in Table A1. In addition, the density of Terfenol-D and silica are taken as 9250 kg·m −3 and 2203 kg·m −3 , respectively. The Maxwell stress tensor in Equation (A8) is neglected in Equation (A10). This is because the contribution of the Maxwell stress tensor and magnetostrictive stress (taking the value from Table A1) is comparable only when the driving magnetic field is no less than the order of 10 Tesla under linear stress-magnetic-field approximation. Due to the large piezomagnetic constant of magnetostrictive materials, the influence from the Maxwell stress tensor in the Terfenol-D can be safely neglected in experimental condition where magnetic field is well below microtesla (magnetostrictive stress is ∼ 10 6 times larger than the Maxwell stress). The input external driving stress from Equation (A10) is fitted into a 3×3 matrix in the form of Equation (A6) in COMSOL as external stress after simple matrix product calculation. Table A1. Coefficients in the magnetomechanical coupling biased at 60 kA/m and prestressed at 20 MPa [13]. λ H is the elasticity matrix element, e is the piezomagnetic constant and µ σ is the relative magnetic permeability.

Appendix D. Calculation of c act by Lorentzian Fit
Actuation parameter c act is extracted by fitting the equation of motion to the Lorentzian distribution. In COMSOL, the phase of the mechanical displacement spectrum is in the range 0 to −π/2 at frequencies below the resonance frequency and π/2 to 0 above it, this generates an artefact that we remove by taking the absolute value of the displacement. The modified mechanical frequency response given by the Fourier transform of Equation (6) for each single mechanical mode is given by All the fitting parameters are on the RHS of Equation (A11): c act , Γ and Ω M , while all the three parameters on the LHS of Equation (A11), maximum displacement, effective mass and signal magnetic field, can be drawn from COMSOL. The COMSOL syntax for extracting the maximum displacement is sqrt(abs(u^2)+abs(v^2)+abs(w^2)) under volume maximum analysis and the m eff is the quotient of solid.rho*(abs(w^2)+abs(v^2)+abs(u^2)) under volume integration analysis and maximum displacement where u,v,w are displacements in x, y, z directions.
Since damping is input manually, it is important to check whether the damping affects c act or not. By changing the input damping for a wide parameter range of 12500 times variation, it has been verified that the effect of the damping variation on c act is negligible. Figure A2a shows the fit to the Lorentzian distribution for an input damping factor 12.5 times smaller than that of in Figure A2b, and the fit in Figure A2c has the damping factor 12500 times larger than that of in Figure A2a. The fitting is performed on a mechanical mode simulated in Figure A2d. A smaller input damping allows for a clearer Lorentzian fit as shown in Figure A2a. Therefore, in the implementation, input damping is chosen to be as small as possible limited by COMSOL's convergence error.  Figure A2. (a) a fit to the Lorentzian distribution, and ξ spectrum for a small input material damping; (b) the fit and ξ spectrum with input damping factor 12.5 times larger than that of a); (c) the fit with damping factor 12500 times larger than that of a); (d) all the fits and ξ spectrum are performed on the mechanical mode with Terfenol-D position offset from the centre. Figure A3 shows the actuation parameter c act of the magnetometer [8] using the parameters in Section 5 in three frequency sections of mechanical modes. Blue dots present the c act for mechanical modes under in-plane magnetic field driving, while orange dots are for mechanical modes under axial driving. As can be seen, the actuation parameter varies over many orders of magnitude both for different mechanical modes with the same driving magnetic field direction (comparison of dots among the same colour) and for the same mode driving by the magnetic field in two orthogonal directions (in comparison between red and blue dots at the same mechanical eigenfrequency, connected via a black vertical line in Figure A3).

Appendix E. Calculation of Geometrical Factor ξ
This numerical modelling offers an accurate way to compute the effective cooperativity C eff from the geometrical factor ξ defined in Equation (9). The optical path length L corresponds to the outermost circumference in micro-toroidal resonators. The change of circumference δL(Ω) in a harmonic oscillation at a mechanical resonance can be obtained by synchronizing the output displacement to the maximum and minimum amplitudes. The implementation of the synchronization in COMSOL is to multiply a phase shift in all three Cartesian coordinates in the deformable mesh setting as dx=u*exp(-i*atan(imag(u)/(real(u)+1e-16))), dy=v*exp(-i*atan(imag(v)/(real(v)+1e-16))), dz=w*exp(-i*atan(imag(w)/(real(u)+1e-16))), where 1e-16 in the denominator is an example of adding a small value to eliminate the error of dividing by 0 when real(u)/real(v)/real(w) is 0, and dx,dy,dz are displacements after the synchronization. A line integral with integrand 1 along the outermost circumference of the devices is followed after the phase synchronization, which results in the circumference at a particular phase. Linear optomechanical coupling is evaluated through δL by taking the difference of circumference synchronized at phase 0 and phase π where Equation (A12) is multiplied by a phase factor exp(i*pi). Figure A2a,b shows the constant ξ across a mechanical resonance with 12.5 times variation of input damping, showing that ξ is insensitive to the mechanical quality. ξ is missing in Figure A2c due to the numerical error at large manually input damping where ξ spectrum is far away from constant. Plotting ξ spectrum offers a way of sanity check of possible numerical errors.

Appendix F. From ξ to Optomechanical Coupling
With the calculated geometrical factor ξ, the value of the parameters describing optomechanical coupling C eff is easy to achieve. The optomechanical coupling rate is defined as g 0 ≡ G · x zpf , where zero-point motion x zpf = h/(2m eff Ω M ), and G (rad·m −1 s −1 ) is the optomechanical coupling strength quantifying the shift of optical resonance frequency δω 0 by the mechanical displacement as For a Fabry-Pérot type and micro-toroidal structure cavity with length L, the shift of the optical resonance frequency is linked with the change of the cavity length by δL/L = δω 0 /ω 0 . Inserting the expression δL/L = δω 0 /ω 0 into Equation (A13) leads to Therefore, g 0 and thereby C eff can then be written as a function of the geometrical factor ξ. For an individual mechanical eigenmode at Ω M and based on Equation (A5), the expression for the effective cooperativity as a function of the geometrical factor ξ can be calculated to where the front part on the RHS is mode dependent, ξ(Ω M ) and m eff (Ω M ) can be accurately extracted from the numerical simulation. The optical resonance frequency ω 0 and length of the cavity L are input parameters, and the only empirical parameters left are the optical decay κ and the mechanical damping factor Γ with the assumption of an idealized lossless cavity escape efficiency η esc = 1.