Magnetoelastic Coupling and Delta-E Effect in Magnetoelectric Torsion Mode Resonators

Magnetoelectric resonators have been studied for the detection of small amplitude and low frequency magnetic fields via the delta-E effect, mainly in fundamental bending or bulk resonance modes. Here, we present an experimental and theoretical investigation of magnetoelectric thin-film cantilevers that can be operated in bending modes (BMs) and torsion modes (TMs) as a magnetic field sensor. A magnetoelastic macrospin model is combined with an electromechanical finite element model and a general description of the delta-E effect of all stiffness tensor components Cij is derived. Simulations confirm quantitatively that the delta-E effect of the C66 component has the promising potential of significantly increasing the magnetic sensitivity and the maximum normalized frequency change Δfr. However, the electrical excitation of TMs remains challenging and is found to significantly diminish the gain in sensitivity. Experiments reveal the dependency of the sensitivity and Δfr of TMs on the mode number, which differs fundamentally from BMs and is well explained by our model. Because the contribution of C11 to the TMs increases with the mode number, the first-order TM yields the highest magnetic sensitivity. Overall, general insights are gained for the design of high-sensitivity delta-E effect sensors, as well as for frequency tunable devices based on the delta-E effect.


Introduction
In recent years, thin-film magnetoelectric sensors have been studied, frequently envisioning biomedical applications in the future [1,2]. Such applications often require the measurement of small amplitude and low frequency magnetic fields [1][2][3]. With the direct magnetoelectric effect, such small detection limits are only obtained at high frequencies and in small-signal bandwidths of a few Hz [2,4]. One way to overcome these limitations is by using a modulation scheme based on the delta-E effect. The delta-E effect is the change of the effective elastic properties with magnetization due to magnetoelastic coupling [5][6][7][8]. It results from inverse magnetostriction that adds additional stress-induced magnetostrictive strain to the purely elastic Hookean strain. The delta-E effect can occur generally in various elastic moduli and several components of the elastic stiffness tensor C [9,10]. Hence, it is sometimes referred to as the delta-C effect [11]. Typically, delta-E effect sensors are based on magnetoelectric resonators that are electrically excited via the piezoelectric layer at or close to the resonance frequency f r . Upon the application of a magnetic field, the magnetization changes and the delta-E effect alters the mechanical stiffness tensor of the magnetostrictive layer. If the altered stiffness tensor components contribute to the resonance frequency of the excited mode, the resonance frequency changes, which can be read-out electrically. The delta-E effect of the Young's modulus has especially been studied thoroughly in soft magnetic amorphous materials [12][13][14][15][16][17][18]. It was used for magnetic field sensing with magnetoelectric plate resonators [19][20][21][22] and beam structures [23][24][25][26][27][28][29][30][31][32]. Such resonators are operated in bending or bulk modes and some have achieved limits of detection down to the sub-nT regime at low frequencies. Microelectromechanical systems (MEMS) cantilever sensors based on the delta-E effect were recently used for the mapping of magnetically labeled cells [33], and have shown promising properties for sensor array applications [34].
In contrast to the delta-E effect of the Young's modulus, the delta-E effect of the shear modulus has been studied less extensively [35] and mainly in amorphous wires [36,37]. It has been used for a different kind of delta-E effect sensors where shear waves, traveling through the magnetoelastic material, are influenced by the delta-E effect. This concept was realized with bulk acoustic shear waves in amorphous ribbons [38] and recently with surface acoustic shear waves in magnetic thin film devices [10,[39][40][41][42]. Only very few studies investigate torsion modes in beam structures [43,44], either with electrostatically actuated cantilevers [43] or double-clamped beams [44]. Both studies are limited to specific configurations of the magnetic system and consider neither the full tensor relations of the mechanics and the delta-E effect nor higher resonance modes. Until now, a comprehensive experimental and theoretical analysis has been missing as well as a discussion of implications for the design of delta-E effect-based devices.

MEMS Torsion Mode Sensors
In this study, all measurements and models are made for a microelectromechanical system (MEMS) technology-fabricated cantilever with an electrode design that permits the excitation of torsion modes. A sketch including dimensions and layer structure and a top-view photograph of the design are shown in Figure 1. The approximately 3.1 mm-long and 2.15 mm-wide cantilever consists of a ≈ 2 µm-thick piezoelectric layer of AlN [45] on a 50 µm-thick poly-Si substrate. A 2 µm-thick amorphous magnetostrictive multilayer is deposited on the rear side. A magnetic field is applied during the deposition to induce a magnetic easy axis along the short cantilever axis. For actuation and read-out, three top electrodes (E 1 , E 2 and E 3 ) of 100 nm-thick Au with lengths L 1 = L 2 ≈ 1 mm and L 3 ≈ 0.6 mm and widths W 1 = W 2 ≈ 0.5 mm and W 3 ≈ 1 mm contact the AlN layer on the top. The counter electrode (150 nm Pt) covers the whole beam area and is located between the AlN layer and the substrate. All measurements are performed with electrode E 1 . As a magnetostrictive material, we use a 2 µm multilayer of 20 × (100 nm (Fe 90 Co 10 ) 78 Si 12 B 10 and 6 nm Cr). It is covered by a top Cr-layer that serves as a protection against corrosion. More information about the layer structures and the fabrication process can be found elsewhere [27]. In contrast to the sensors in Ref. [27], the sensor presented here is significantly wider and the adapted electrode design additionally permits the excitation of torsion modes. Details on the geometry are given in the appendix.
In contrast to the delta-E effect of the Young's modulus, the delta-E effect of the shear modulus has been studied less extensively [35] and mainly in amorphous wires [36,37]. It has been used for a different kind of delta-E effect sensors where shear waves, traveling through the magnetoelastic material, are influenced by the delta-E effect. This concept was realized with bulk acoustic shear waves in amorphous ribbons [38] and recently with surface acoustic shear waves in magnetic thin film devices [10,[39][40][41][42]. Only very few studies investigate torsion modes in beam structures [43,44], either with electrostatically actuated cantilevers [43] or double-clamped beams [44]. Both studies are limited to specific configurations of the magnetic system and consider neither the full tensor relations of the mechanics and the delta-E effect nor higher resonance modes. Until now, a comprehensive experimental and theoretical analysis has been missing as well as a discussion of implications for the design of delta-E effect-based devices.

MEMS Torsion Mode Sensors
In this study, all measurements and models are made for a microelectromechanical system (MEMS) technology-fabricated cantilever with an electrode design that permits the excitation of torsion modes. A sketch including dimensions and layer structure and a top-view photograph of the design are shown in Figure 1. The approximately 3.1 mmlong and 2.15 mm-wide cantilever consists of a ≈ 2 µm-thick piezoelectric layer of AlN [45] on a 50 µm-thick poly-Si substrate. A 2 µm-thick amorphous magnetostrictive multilayer is deposited on the rear side. A magnetic field is applied during the deposition to induce a magnetic easy axis along the short cantilever axis. For actuation and read-out, three top electrodes ( , and ) of 100 nm-thick Au with lengths = ≈ 1 mm and ≈ 0.6 mm and widths = ≈ 0.5 mm and ≈ 1 mm contact the AlN layer on the top. The counter electrode (150 nm Pt) covers the whole beam area and is located between the AlN layer and the substrate. All measurements are performed with electrode . As a magnetostrictive material, we use a 2 µm multilayer of 20 (100 nm (Fe Co ) Si B and 6 nm Cr). It is covered by a top Cr-layer that serves as a protection against corrosion. More information about the layer structures and the fabrication process can be found elsewhere [27]. In contrast to the sensors in Ref. [27], the sensor presented here is significantly wider and the adapted electrode design additionally permits the excitation of torsion modes. Details on the geometry are given in the appendix.  Delta-E effect sensor analyzed in this study: (a) schematic top view of the cantilever, with three different electrodes E 1 , E 2 and E 3 of lengths L 1 = L 2 ≈ 1 mm, L 3 ≈ 0.6 mm and widths W 1 = W 2 ≈ 0.5 mm and W 3 ≈ 1 mm; (b) schematic side-view of the cantilever with the thickness of the functional layers and the poly-Si substrate; (c) top-view photograph of the fabricated structure.

Definition of the Sensitivity
An important parameter that characterizes a magnetic field sensor is its sensitivity. During sensor operation, an alternating voltage is applied to excite the cantilever at its mechanical resonance frequency f r . Applying a magnetic field, shifts f r via the delta-E effect and correspondingly the sensor's admittance characteristic on its frequency axis f. Hence, the magnitude |Y| = abs{Y} and phase angle φ = arg{Y} of the sensor admittance Y depend on the magnetic field. Consequently, the ac magnetic field to be measured causes an amplitude modulation (am) and phase modulation (pm) of the current through the sensor. Detailed information on the operation and read-out can be found elsewhere [46][47][48]. The linearized change of |Y| and φ with the magnetic field can be described by the amplitude sensitivity S am = S Y,r ·S H,r and the phase sensitivity S pm = S φ,r ·S H,r [49], respectively. Both sensitivities have a magnetic part S H,r that includes the delta-E effect and an electric part S Y,r or S φ,r , which can be determined from the admittance. We refer to the three sensitivities as relative sensitivities, because they are normalized to the excitation frequency f ex = f r . The normalization is required to eventually compare the electrical and magnetic sensitivities of sensors with different geometries operated at different f r or in different resonance modes. Usually a magnetic bias field H 0 is applied to operate the sensor at optimum conditions. The relative sensitivities are then defined in linear approximation as derivatives [49]: with the magnetic vacuum permeability µ 0 ≈ 4π·10 −7 N/A 2 . From Equation (1), the relative magnetic sensitivity S H,r is the linearized and normalized change of the resonance frequency f r with the applied magnetic flux density µ 0 H.

Magnetic Sensitivity of Arbitrary Resonance Modes
The delta-E effect is included in the relative magnetic sensitivity S H,r because the resonance frequency f r = f r C ij is a function of the stiffness tensor components C ij . Depending on the respective resonance mode, different C ij dominate f r and depending on the magnetoelastic properties they might result in non-zero S H,r . To describe S H,r for arbitrary resonance modes, it can be separated into a purely mechanical part f −1 r ∂ f r /∂C ij that contains the resonance properties of the structure and a purely magnetoelastic part ∂C ij /∂µ 0 H: If treated separately, the factors ∂ f r /∂C ij and ∂C ij /∂µ 0 H must be normalized to remove the dependency on the absolute value of C ij that cancels out in S H,r . We define: From Equation (3), the factor ∂ C f r,ij represents a normalized measure for the influence of the stiffness tensor component C ij on the resonance frequency f r of the considered resonance mode. It is a purely mechanical quantity and hence determined by the geometry, the resonance mode, and the effective mechanical properties of the resonator. The second factor ∂ H C ij , includes the delta-E effect and describes the normalized influence of the applied flux density µ 0 H on C ij . Hence, the two factors quantify the mechanical and the magnetoelastic parts of the relative magnetic sensitivity S H,r . They will be used later to analyze the sensitivity and the frequency detuning of higher bending and torsion modes of the cantilever.

Sensor Modelling
The model used to describe and analyze the sensor consists of two parts. With a semi-analytical magnetoelastic macrospin model, the delta-E effect is obtained, i.e., the effective mechanical stiffness tensor C(H) as a function of the applied field H. It is used as an input for an electromechanical finite element mechanics (FEM) model that describes the resonance frequency and the sensor's impedance response. In addition to a macrospin approximation, we assume a quasi-static magnetization behavior. Consequently, it is only valid for operation frequencies and magnetic field frequencies far below the ferromagnetic resonance frequency (FMR). The FMR generally depends on the geometry and the magnetic properties of the thin-film [50] and can cause a frequency dependency of the delta-E effect [51]. For the soft-magnetic material and thin-film geometry used here, it is in the GHz regime [51,52]. Because the operation frequencies are of the order of several kHz, magnetodynamic effects and the frequency dependency of the delta-E effect are neglected [51]. Due to the low frequencies, we assume that also electrodynamic effects can be omitted in the electromechanical model. In the following, both parts of the model are discussed in detail.

Electromechanical Model
In the electromechanical part of the model, we consider a simplified cantilever geometry, reduced to the poly-Si substrate, the piezoelectric AlN layer, and the magnetic FeCoSiB layer. Details on the geometry used are given in Appendix C. We assume all materials to be mechanically linear, which is a good approximation at sufficiently small excitation voltages. The material parameters used are given in the appendix. The cantilever is oriented in a cartesian coordinate system as illustrated in Figure 2, used throughout this paper. The mechanical equation of motion is given by (e.g., [53]) if no external forces are present. It includes the displacement vector u, the time t, the mass density ρ, and the divergence ∇·σ of the mechanical stress tensor σ. For sufficiently small excitation frequencies, eddy current effects can be neglected and the electrostatic equations [54]: are valid in good approximation. They include the electrical vector field E, the gradient ∇U of the electrical potential U and the divergence ∇·D of the electric flux density D with the free charge density ρ c . The electrostatic equations are coupled to the mechanical equation of motion via the constitutive piezoelectric equations, here in the stress-charge form [54,55]: with the linear strain tensor ε and the complex mechanical stiffness tensor C * = C(1 + iη). Its real part is the material's stiffness tensor C and its imaginary part Cη includes the isotropic loss factor η, which is used to consider damping in the materials [56]. The electromechanical coupling tensor is denoted as e c and the electrical permittivity tensor as ε el . For the calculation, we set fixed boundary conditions (u = 0) at the left face of the beam. For the piezoelectric material, we assumed at the boundaries nD = 0 (with surface normal vector n), and an initial value for the electric potential of U = 0, except for the area covered by the electrodes. The electrodes are modeled with a fixed potential boundary condition, where an alternating voltage with amplitude U 0 , the angular frequency ω and phase angle ϕ v . To calculate the electrical admittance Y = I/U the current I is obtained from integrating the surface charge density over the electrode areas. For the solution, a linear response of the system is assumed, with a displacement of the form u =û· exp(i[ωt + ϕ u ]) and a solution for the electrical potential . The equations are solved within a frequency domain study in COMSOL (r) Multiphysics v. 5.3a (COMSOL AB, Stockholm, Sweden) [56]. All material parameters used are given in the appendix.
. For the calculation, we set fixed boundary conditions ( = 0) at the left face of the beam. For the piezoelectric material, we assumed at the boundaries = 0 (with surface normal vector ), and an initial value for the electric potential of = 0, except for the area covered by the electrodes. The electrodes are modeled with a fixed potential boundary condition, where an alternating voltage = • exp ( [ + ]) is applied, with amplitude , the angular frequency and phase angle . To calculate the electrical admittance = ⁄ the current is obtained from integrating the surface charge density over the electrode areas. For the solution, a linear response of the system is assumed, with a displacement of the form = • exp ( [ + ]) and a solution for the electrical potential of = • exp ( [ + ]). The equations are solved within a frequency domain study in COMSOL ® Multiphysics v. 5.3a (COMSOL AB, Stockholm, Sweden) [56]. All material parameters used are given in the appendix.

Magnetoelastic Model
For the magnetic model, we consider the enthalpy density function of a macrospin with a uniaxial anisotropy energy density, an external magnetic field, a demagnetizing term, and magnetoelastic energy density. Using Einstein's summation convention, the enthalpy density term we use is: In this equation, the components of the reduced magnetization vector are denoted by , the magnitude of the magnetization vector by and the magnetic vacuum permeability by . The effective easy axis of magnetization is characterized by its orientation vector and the effective first-order uniaxial anisotropy energy density constant . The components of the external magnetic field vector are given by and the components of the mean demagnetizing field by , = − , with the main diagonal components of the demagnetizing tensor. For the magnetoelastic energy density, we use the coupling term − with the stress tensor components and the components of the isotropic magnetostrictive strain tensor. Both are given in Voigt's notation. The coupling term results from omitting magnetostrictive self-energy and incorporating the term constant with stress into K [57]. In the following, the polar angle and the azimuthal angle of in the spherical coordinate system ( Figure 2) are used to define its components . The exact definition of all vector and tensor components is given in the appendix. The linearized change of the elastic compliance components with the magnetic field and stress is derived from the expression where the first summand , is the constant, fixed magnetization elastic compliance tensor component. The magnetization dependent part ∆ can be obtained from the equilibrium conditions that are given by the first-order derivatives of :

Magnetoelastic Model
For the magnetic model, we consider the enthalpy density function of a macrospin with a uniaxial anisotropy energy density, an external magnetic field, a demagnetizing term, and magnetoelastic energy density. Using Einstein's summation convention, the enthalpy density term we use is: In this equation, the components of the reduced magnetization vector are denoted by m i , the magnitude of the magnetization vector by M s and the magnetic vacuum permeability by µ 0 . The effective easy axis of magnetization is characterized by its orientation vector EA i and the effective first-order uniaxial anisotropy energy density constant K. The components of the external magnetic field vector are given by H i and the components of the mean demagnetizing field by H d,i = −D ii m i M s , with the main diagonal components D ii of the demagnetizing tensor. For the magnetoelastic energy density, we use the coupling term −σ i λ i with the stress tensor components σ i and the components λ i of the isotropic magnetostrictive strain tensor. Both are given in Voigt's notation. The coupling term results from omitting magnetostrictive self-energy and incorporating the term constant with stress into K [57]. In the following, the polar angle θ and the azimuthal angle ϕ of m in the spherical coordinate system (Figure 2) are used to define its components m i . The exact definition of all vector and tensor components is given in the appendix. The linearized change of the elastic compliance components S ij with the magnetic field and stress is derived from the expression where the first summand S m,ij is the constant, fixed magnetization elastic compliance tensor component. The magnetization dependent part ∆S ij can be obtained from the equilibrium conditions that are given by the first-order derivatives of u: From these equilibrium conditions a general expression for the linearized change ∆S ij of the compliance tensor can be derived (Appendix A). Denoting the second-order derivatives as u ϕϕ and u θθ it is: This expression permits a quick calculation of the compliance tensor for different magnetic systems described by an enthalpy density u. From Equation (10), the non-zero components of ∆S for in-plane magnetization (θ = π/2) are: The final compliance tensor for in-plane magnetization as a function of magnetic field and stress is: Because in our case both, S m and ∆S. are symmetric, and S is also symmetric. Note that S m,16 = S m,26 = S m,45 = 0 in our isotropic magnetic material and consequently S 16 = ∆S 16 , S 26 = ∆S 26 and S 45 = ∆S 45 . Finally, the stiffness tensor C is obtained by numerically calculating the inverse C(H, σ) = S(H, σ) −1 . It has the same non-zero components and symmetry. All equations (Equations (11)-(17) are obtained from Equation (10) assuming in-plane magnetization (θ = π/2) and are valid for the isotropic magnetoelastic coupling used in the enthalpy density function (Equation (7)). For all the following simulations, we additionally assume in-plane magnetic fields (θ H = π/2) and an in-plane easy axis (θ EA = π/2).
These two assumptions influence and simplify u and its derivatives, which are given in the appendix.

Implications of the Magnetic Model
In the following, results for the C ij of the magnetoelastic model are discussed at the example of a thin-film geometry. For the calculations, we assumed zero static stress (σ i = 0) and D 33 = 1. The large shape anisotropy results in C 44 ≈ C m,44 and C 55 ≈ C m, 55 . We limit the discussion to the C 11 , C 12 and C 66 components as they are most relevant for torsion and bending modes. In Figure 3a, the normalized C 11 , C 12 and C 66 components are plotted for a macrospin and ϕ EA = 90 • . Because u ϕϕ (H = H K ) = 0 and so ∆S 66 66 only for H → ∞ . Hence, for finite H even a small shear stress σ 6 can always tilt the magnetization vector out of the applied magnetic field direction. It occurs, because the magnetoelastic energy density contribution −σ 6 λ 6 of the shear stress σ 6 is asymmetric around ϕ = 0 • . Its minimum is shifted by 45 • compared to the minimum of the one-component at ϕ = 0 • . Consequently, at the two local maxima it is 66 . The C 11 component shows two distinct minima but unlike the delta-E effect in the Young's modulus (e.g., [6,14,49]) no discontinuities at |H| = |H K |. Although the discontinuities are present in S 11 (not shown), they vanish during the inversion due to contributions of other S ij components to C 11 . In contrast to C 11 , C 12 stiffens with applied magnetic bias field because ∆S 12 = −∆S 11 . The signs are a direct consequence of the positive isotropic magnetoelastic coupling. As the macrospin rotates towards the x axis, magnetostrictive expansion occurs along the x axis, but contraction occurs along the y axis. Compared with C 11 , the maximum relative change of C 12 is larger because S m,12 < S m,11 , which results in a different weighting in Equation (8). In Figure 3b, C 66 is shown for three different angles of the easy axis ϕ EA = 90 • , 85 • , and 75 • . It is apparent that a change of ϕ EA strongly influences C 66 . Relative to ϕ EA = 90 • , the two minima at H = ±H K shift to a larger |H| and the minimum value increases strongly by more than 85% at ϕ EA = 85 • and about 95% at ϕ EA = 75 • . The center minimum shifts due to the single domain hysteresis and decreases slightly with decreasing ϕ EA . A singularity occurs at ϕ EA = 85 • due to the magnetic discontinuity at the switching field of the single-spin model. Due to the strong impact of small deviations from ϕ EA = 90 • on C 66 (H), the magnetic sensitivity is expected to change notably with ϕ EA .
Sensors 2021, 21, x FOR PEER REVIEW 7 of 18 limit the discussion to the , and components as they are most relevant for torsion and bending modes.
Hence, for finite even a small shear stress can always tilt the magnetization vector out of the applied magnetic field direction. It occurs, because the magnetoelastic energy density contribution − of the shear stress is asymmetric around = 0°. Its minimum is shifted by 45° compared to the minimum of the one-component at = 0°. Consequently, at the two local maxima it is ( = 45°, 135°) = , . The component shows two distinct minima but unlike the delta-E effect in the Young's modulus (e.g., [6,14,49]) no discontinuities at | | = | |. Although the discontinuities are present in (not shown), they vanish during the inversion due to contributions of other components to . In contrast to , stiffens with applied magnetic bias field because ∆ = −∆ . The signs are a direct consequence of the positive isotropic magnetoelastic coupling. As the macrospin rotates towards the x axis, magnetostrictive expansion occurs along the x axis, but contraction occurs along the y axis. Compared with , the maximum relative change of is larger because , < , , which results in a different weighting in Equation (5). In Figure 3b, is shown for three different angles of the easy axis = 90°, 85°, and 75°. It is apparent that a change of strongly influences . Relative to = 90°, the two minima at = ± shift to a larger | | and the minimum value increases strongly by more than 85% at = 85° and about 95% at = 75°. The center minimum shifts due to the single domain hysteresis and decreases slightly with decreasing . A singularity occurs at = 85° due to the magnetic discontinuity at the switching field of the single-spin model. Due to the strong impact of small deviations from = 90° on ( ) , the magnetic sensitivity is expected to change notably with . ( ) (Equation (3)) for the and components as functions of the easy axis angle . For the calculations in (c), a distribution of effective anisotropy energy density is used as in Ref. [49] with a standard deviation of = 15 % for a more quantitative estimation and to prevent singularities.
In the following, we quantify the influence of the components on the relative magnetic sensitivity , using as defined in Equation (3). Calculating requires forming the derivate ⁄ , which results in singularities for the 66-component  (3)) for the C 11 and C 66 components as functions of the easy axis angle ϕ EA . For the calculations in (c), a distribution of effective anisotropy energy density is used as in Ref. [49] with a standard deviation of δ K = 15 % for a more quantitative estimation and to prevent singularities.
In the following, we quantify the influence of the C ij components on the relative magnetic sensitivity S H,r using ∂ H C ij as defined in Equation (3). Calculating ∂ H C ij requires forming the derivate ∂C ij /∂H, which results in singularities for the 66-component at ϕ EA = 90 • and |H| = |H K |. A finite derivative can be estimated by including the distribution of the effective anisotropy energy density K in a mean-field approach [15,49,58]. With such a distribution, inhomogeneities in the magnetization response are considered that can occur, e.g., from spatially varying stress or internal stray fields. We use a normal distribution of K with a standard deviation δ K = 15 % as a representative example value that has been used previously for a similar device [49]. We calculate ∂ H C ij (H) numerically from C ij (H) and extract the maximum ∂ H C ij,max (H) for H > 0, at various angles ϕ EA of the easy axis. They are plotted in Figure 3c. As a result of the distribution, both ∂ H C ij,max are finite at ϕ EA = 90 • with ∂ H C 66,max ≈ 10 × ∂ H C 11,max ≈ 4.5 × ∂ H C 12,max . This is reduced to ∂ H C 66,max ≈ 4 × ∂ H C 11,max ≈ 2 × ∂ H C 12,max at ϕ EA = 80 • . In conclusion, the C 66 component potentially offers a significantly larger magnetic sensitivity than the C 11 and C 12 components.

Magnetization Measurements
Magneto-optical Kerr effect (MOKE) microscopy [59] was used to analyze the magnetic multilayer. The picture in Figure 4a shows the rear side of the cantilever and is composed of a series of images. For each image, the magnetic multilayer was demagnetized along the x axis and the MOKE sensitivity axis was set along the y axis. The region of the magnetic multilayer is marked with a white frame and the estimated easy axis orientation is indicated with white arrows. In a large region around the left, top, and bottom edge, no magnetic response is visible. A comparison with light microscopy images reveals possibly corroded regions. They might have formed due to incomplete Cr-coverage at the edges. At the time, a particularly thin Cr-layer was deposited to ensure good magnetooptical contrast. Close to the clamping region (blue rectangle in Figure 4a), the layer is partially delaminated. Despite these nonidealities, the overall magnetic response in the magnetically active region is quite homogeneous. The average easy axis orientation is approximately ϕ EA = −75 • ± 5 • relative to the x axis. An effective uniaxial anisotropy energy density of K = (1.2 ± 0.1) kJ/m 3 is estimated with the magnetoelastic model. We used the ballistic demagnetizing tensor [60] in the center of the film and assumed σ j = 0. A representative magnetization curve of the center region, recorded along the x axis, is shown in Figure 4b, and compared with one recorded at the clamping region. The difference between these curves indicates a different alignment of the effective anisotropy. However, due to the magnetic multilayer structure and the partial delamination additional effects cannot be excluded. From previous investigations [49], we expect that the deteriorated magnetic properties at the clamping especially reduce the resonance frequency detuning and the magnetic sensitivity of the first bending mode.

at
= 90° and |H| = |HK|. A finite derivative can be estimated by including the distribution of the effective anisotropy energy density K in a mean-field approach [15,49,58]. With such a distribution, inhomogeneities in the magnetization response are considered that can occur, e.g., from spatially varying stress or internal stray fields. We use a normal distribution of with a standard deviation = 15 % as a representative example value that has been used previously for a similar device [49]. We calculate ( ) numerically from ( ) and extract the maximum , ( ) for > 0, at various angles of the easy axis. They are plotted in Figure 3c. As a result of the distribution, both

Magnetization Measurements
Magneto-optical Kerr effect (MOKE) microscopy [59] was used to analyze the magnetic multilayer. The picture in Figure 4a shows the rear side of the cantilever and is composed of a series of images. For each image, the magnetic multilayer was demagnetized along the x axis and the MOKE sensitivity axis was set along the y axis. The region of the magnetic multilayer is marked with a white frame and the estimated easy axis orientation is indicated with white arrows. In a large region around the left, top, and bottom edge, no magnetic response is visible. A comparison with light microscopy images reveals possibly corroded regions. They might have formed due to incomplete Cr-coverage at the edges. At the time, a particularly thin Cr-layer was deposited to ensure good magnetooptical contrast. Close to the clamping region (blue rectangle in Figure 4a), the layer is partially delaminated. Despite these nonidealities, the overall magnetic response in the magnetically active region is quite homogeneous. The average easy axis orientation is approximately = −75°± 5° relative to the x axis. An effective uniaxial anisotropy energy density of = (1.2 ± 0.1) kJ/m³ is estimated with the magnetoelastic model. We used the ballistic demagnetizing tensor [60] in the center of the film and assumed = 0. A representative magnetization curve of the center region, recorded along the x axis, is shown in Figure 4b, and compared with one recorded at the clamping region. The difference between these curves indicates a different alignment of the effective anisotropy. However, due to the magnetic multilayer structure and the partial delamination additional effects cannot be excluded. From previous investigations [49], we expect that the deteriorated magnetic properties at the clamping especially reduce the resonance frequency detuning and the magnetic sensitivity of the first bending mode.

Electromechanical Properties
To analyze the electromechanical properties of the sensor, the sensor admittance Y( f ex ) is measured over a large range of excitation frequencies f ex . Six resonance modes are characterized in detail by fitting a modified Butterworth van Dyke (mBvD) model (e.g., [61]) with the equivalent circuit configuration from [47] to the measurements. The resonance frequencies f r and quality factors are calculated from the mBvD parameters of each admittance curve and compared with the eigenfrequencies obtained from the finite element method (FEM) model. With this comparison, the eigenmodes are identified to be the first three bending modes (BM1-3) and the first three torsion modes (TM1-3). The FEM model was fitted to admittance measurements of the first torsion mode (TM1) close to magnetic saturation at µ 0 H = − 10 mT. It matches the measurements very well as shown in Figure 5a. The material parameters match excellently with literature values. Details on the material parameters and on the geometry are given in Appendix C. A comparison of the measured resonance frequencies in magnetic saturation with the FEM simulations results in extremely small deviations <2% for all six modes (Appendix B). indicated with arrows at approximately −75°± 5° relative to the x axis; (b) magnetization curve close to the clamping and in the center of the magnetic film. The two evaluated regions are indicated with squares in (a).

Electromechanical Properties
To analyze the electromechanical properties of the sensor, the sensor admittance ( ) is measured over a large range of excitation frequencies . Six resonance modes are characterized in detail by fitting a modified Butterworth van Dyke (mBvD) model (e.g., [61]) with the equivalent circuit configuration from [47] to the measurements. The resonance frequencies and quality factors are calculated from the mBvD parameters of each admittance curve and compared with the eigenfrequencies obtained from the finite element method (FEM) model. With this comparison, the eigenmodes are identified to be the first three bending modes (BM1-3) and the first three torsion modes (TM1-3). The FEM model was fitted to admittance measurements of the first torsion mode (TM1) close to magnetic saturation at μ 0 H = − 10 mT. It matches the measurements very well as shown in Figure 5a. The material parameters match excellently with literature values. Details on the material parameters and on the geometry are given in Appendix C. A comparison of the measured resonance frequencies in magnetic saturation with the FEM simulations results in extremely small deviations < 2% for all six modes (Appendix B). The set of material parameters found is used to predict, and compare the impedance characteristic of other cantilever delta-E effect sensors published previously [28]. The sensors differ in their geometry from our torsion mode sensor. They were designed to excite the first and second bending mode with various electrode geometries. For the simulations, we used the same material parameters found for the torsion mode sensors but adjusted the geometry.
As a figure of merit for the electromechanical model, we compared the absolute difference ∆ = − of the phase angle of the electrical admittance . The simulation results are plotted in Figure 5b and compared with values of the torsion modes (Appendix B) measured here, and the bending mode from Ref. [28]. The TMs were measured close to magnetic saturation at μ 0 H =−10 mT to reduce the influence of the delta-E effect. Slight deviations between the measurement and simulation might result from effectively different magnetoelectric coupling factors, e.g., due to the slightly different material parameters, geometric inaccuracies, or stress [62]. In conclusion, the model can estimate the electromechanical properties of the device and the effect of different electrode configurations well. For the application of magnetoelastic resonators as delta-E effect sensors, a high ∆ and hence a high electrical sensitivity is desirable. In comparison to the The set of material parameters found is used to predict, and compare the impedance characteristic of other cantilever delta-E effect sensors published previously [28]. The sensors differ in their geometry from our torsion mode sensor. They were designed to excite the first and second bending mode with various electrode geometries. For the simulations, we used the same material parameters found for the torsion mode sensors but adjusted the geometry.
As a figure of merit for the electromechanical model, we compared the absolute difference ∆φ = φ max − φ min of the phase angle φ of the electrical admittance Y. The simulation results are plotted in Figure 5b and compared with values of the torsion modes (Appendix B) measured here, and the bending mode from Ref. [28]. The TMs were measured close to magnetic saturation at µ 0 H = −10 mT to reduce the influence of the delta-E effect. Slight deviations between the measurement and simulation might result from effectively different magnetoelectric coupling factors, e.g., due to the slightly different material parameters, geometric inaccuracies, or stress [62]. In conclusion, the model can estimate the electromechanical properties of the device and the effect of different electrode configurations well. For the application of magnetoelastic resonators as delta-E effect sensors, a high ∆φ and hence a high electrical sensitivity is desirable. In comparison to the bending modes, the ∆φ of the torsion modes is systematically smaller, which is also reflected in the electrical sensitivities. With S Y,r ≈ 0.85 mS of TM1, the maximum relative electrical amplitude sensitivity S Y,r ≈ 5.8 mS of sample No. 7 (BM2) [28] is almost a factor of seven larger, despite a similar quality factor. Hence, the large factor potentially gained in the magnetic sensitivity from utilizing the C 66 component can be diminished by a reduced electrical sensitivity.
Additional simulations show that further optimization of the electrode design and reduction in the parasitic capacity from bond pads and wires could improve ∆φ of TM1 to ∆φ = 10 • . Alternatively, the parasitic effect of the sensor capacitance could be neutralized with additional electronics to utilize the phase-modulated signal for magnetic field detection [48]. A further improvement by a factor of two could be obtained by exciting both electrodes E 1 and E 2 , phase shifted by 180 • . Additionally, alternative piezoelectric materials with larger piezoelectric coefficients, such as AlScN [63][64][65] could increase the electrical sensitivity significantly and result in ∆φ comparable with bending modes.

Delta-E Effect and Sensitivities
The f r (H) plots extracted from the modified Butterworth van Dyke (mBvD) fits of the first three bending modes (BM1-3) are shown in Figure 6a (right). They are normalized to f r,max ∆ f r (−10 mT) and have a respective minimum resonance frequency f r,min . As a measure for the maximum resonance frequency detuning, we defined the normalized resonance frequency change f r ( f r,max − f r,min )/ f r,max . All three curves are w-shaped and ∆ f r increases with increasing mode number. This effect was reported previously and explained with a strong weighting of the magnetic properties at the clamping in BM1 [49]. Here, the difference between the BM1 and BM2 is significantly larger, which is consistent with the deteriorated material around the clamping region, visible in the magneto-optical Kerr effect microscopy (MOKE) images ( Figure 4a). Correspondingly, the relative magnetic sensitivity S H,r ≈ 3.5 T −1 is smallest in BM1 and increases up to S H,r ≈ 9 T −1 in BM2.
bending modes, the ∆ of the torsion modes is systematically smaller, which is also reflected in the electrical sensitivities. With , ≈ 0.85 mS of TM1, the maximum relative electrical amplitude sensitivity , ≈ 5.8 mS of sample No. 7 (BM2) [28] is almost a factor of seven larger, despite a similar quality factor. Hence, the large factor potentially gained in the magnetic sensitivity from utilizing the C66 component can be diminished by a reduced electrical sensitivity.
Additional simulations show that further optimization of the electrode design and reduction in the parasitic capacity from bond pads and wires could improve ∆ of TM1 to ∆ = 10°. Alternatively, the parasitic effect of the sensor capacitance could be neutralized with additional electronics to utilize the phase-modulated signal for magnetic field detection [48]. A further improvement by a factor of two could be obtained by exciting both electrodes and , phase shifted by 180°. Additionally, alternative piezoelectric materials with larger piezoelectric coefficients, such as AlScN [63][64][65] could increase the electrical sensitivity significantly and result in ∆ comparable with bending modes.

Delta-E Effect and Sensitivities
The ( ) plots extracted from the modified Butterworth van Dyke (mBvD) fits of the first three bending modes (BM1-3) are shown in Figure 6a (right). They are normalized to , ≔ (−10 mT) and have a respective minimum resonance frequency , . As a measure for the maximum resonance frequency detuning, we defined the normalized resonance frequency change ∆ ≔ ( , − , )/ , . All three curves are w-shaped and ∆ increases with increasing mode number. This effect was reported previously and explained with a strong weighting of the magnetic properties at the clamping in BM1 [49]. Here, the difference between the BM1 and BM2 is significantly larger, which is consistent with the deteriorated material around the clamping region, visible in the magneto-optical Kerr effect microscopy (MOKE) images ( Figure 4a). Correspondingly, the relative magnetic sensitivity , ≈ 3.5 T is smallest in BM1 and increases up to , ≈ 9 T in BM2. The normalized ( ) plots of the torsion modes (TMs) and their corresponding , are shown in Figure 6 (left). Although the sample is close to magnetic saturation at = −10 mT, all three ( ) curves still exhibit a non-zero slope as expected from the presented theory. The three ( ) curves have a global minimum around = 0, two local The normalized f r (H) plots of the torsion modes (TMs) and their corresponding S H,r are shown in Figure 6 (left). Although the sample is close to magnetic saturation at µ 0 H = −10 mT, all three f r (H) curves still exhibit a non-zero slope as expected from the presented theory. The three f r (H) curves have a global minimum around H = 0, two local minima at around ±2 mT, and two local maxima at about ±1 mT. With an increasing mode number, the local maxima are almost unaffected, whereas ∆ f r significantly decreases. Consequently, the maximum S H,r also decrease with the increasing mode number, here from S H,r = 12.6 T −1 in TM1 to S H,r = 3.0 T −1 in TM3. The trend is notably opposed to the corresponding behavior of the bending modes and will be analyzed and explained in detail in the next section using the magnetoelastic and electromechanical models. Overall, the magnetic sensitivities are in the range of ≈ 10 T −1 also measured with other magnetoelastic resonators in bending and bulk resonance modes [22,49]. At first glance, the similarity of BM and TM in S H,r ∝ ∂ H C ij , might contradict the magnetoelastic model results in Figure 3c.
To resolve this and explain the dependency of the torsion modes on the mode number, the second factor ∂ C f r,ij that contributes to S H,r must be considered.

Resonance Frequency Simulations
In the following, we use the stiffness tensor components from the magnetic model as input in the finite element method (FEM) model to describe and analyze the frequency detuning and the magnetic sensitivity of the bending and torsion modes measured before. The demagnetizing tensor is approximated with the ballistic demagnetizing tensor in the center of the magnetic layer [60]. Consistently with the measurements the easy axis angle is set to ϕ EA = −75 • and the effective anisotropy energy density constant to K = 1.2 kJ/m 3 , assuming σ j = 0. Results for the normalized resonance frequencies f r (H) of the torsion modes are shown in Figure 7a and of the bending modes in Figure 7b. Despite the simplifying assumptions, a striking similarity with the measurements is apparent. All simulated torsion mode (TM) curves in Figure 7a exhibit two local maxima around one global minimum. Due to the single-domain hysteresis, the local minimum is shifted slightly leftwards away from µ 0 H = 0. The frequency difference between the local maxima and the global minimum decreases significantly with increasing mode number, as also observed in the measurements.
Within the model, this phenomenon can be explained with the mode shapes of the higher torsion modes (Figure 7c). Due to the multiple twisting of the cantilever in higher modes, the resonance nodes are closer together. This results in an increasing contribution of the stiffness tensor components C 11 and C 22 to f r relative to the C 66 component. Quantitatively, we can explain the contribution of C ij to f r with the 11-and the 66-components of the normalized frequency factors ∂ C f r,ij (Equation (3)). They are estimated with the FEM model and summarized in Table 1. Whereas ∂ C f r,11 increases by almost a factor of three, ∂ C f r, 66 shows the opposite trend and decreases by a factor of approximately two, from TM1 to TM3. Because the minima and maxima of C 11 and C 66 occur at similar magnetic bias fields they increasingly compensate each other in higher torsion modes. This causes similar magnetic sensitivities of TMs and BMs in our sensor, although ∂ H C 66,max > ∂ H C 11,max in Figure 3c. If the delta-E effect of C 66 is to be utilized, consequently, the first torsion mode is preferable to higher modes. Table 1. Normalized frequency factors ∂ C f r, 11 and ∂ C f r,66 (Equation (3)) of the first three torsion modes (TMs) and bending modes (BMs), calculated with the electromechanical finite element model. In contrast to the measured bending mode curves (Figure 6), the corresponding modeled curves (Figure 7b) are almost independent of the mode number. Consistently, the ∂ C f r,11 of the BMs are approximately constant with the mode number. The other frequency factor ∂ C f r,12 is very small and ∂ C f r,66 ≈ 0. A different effect dominates the mode dependency observed in the measured bending modes. This corroborates the hypothesis stated earlier in Section 6.3 that the reduced maximum normalized resonance frequency change ∆ f r (as defined in Section 6.3) of BM1 is likely caused by the deteriorated magnetic layer present around the clamping (Figure 4a). the simplifying assumptions, a striking similarity with the measurements is apparent. All simulated torsion mode (TM) curves in Figure 7a exhibit two local maxima around one global minimum. Due to the single-domain hysteresis, the local minimum is shifted slightly leftwards away from = 0. The frequency difference between the local maxima and the global minimum decreases significantly with increasing mode number, as also observed in the measurements. Within the model, this phenomenon can be explained with the mode shapes of the higher torsion modes (Figure 7c). Due to the multiple twisting of the cantilever in higher modes, the resonance nodes are closer together. This results in an increasing contribution of the stiffness tensor components C11 and C22 to relative to the C66 component. Quantitatively, we can explain the contribution of Cij to with the 11-and the 66-components As shown earlier in Figure 3a, the minima of C 11 (H) occur at the same magnetic bias fields as the maxima of C 12 (H). Whereas C 11 (H) softens upon the application of a magnetic bias field, C 12 (H) increases. However, upon application of a magnetic field, they both reduce the resonance frequency of bending modes. Consequently, their corresponding frequency factors have opposite signs and ∂ C f r,12 < 0.

Summary and Conclusions
In summary, we provide an experimental and theoretical study on the delta-E effect, the normalized resonance frequency change ∆ f r (defined in Section 6.3) and the sensitivity of first and higher-order bending modes (BMs) and torsion modes (TMs). The study was conducted on a magnetoelectric thin-film cantilever with a soft magnetic FeCoSiB-Cr multilayer and an electrode design that enables the excitation of various resonance modes. A general expression was developed that permits the detailed analysis of the magnetic sensitivity of arbitrary resonance modes. An electromechanical finite element method model was set up to describe the resonator and the electrical sensitivity. It was combined with a magnetoelastic macrospin model to include the tensor of the linearized delta-E effect for isotropic magnetostriction in the approximation of negligible magnetostrictive self-energy. The models are valid for moderately high-operation frequencies, where electrodynamic and magnetodynamic effects can be omitted.
The delta-E effect model is discussed in detail for here the most relevant components C 66 , C 11, and C 12 of the magnetic field-dependent stiffness tensor. Simulation results imply that the C 66 component potentially offers a ten-fold higher contribution to the magnetic sensitivity than the C 11 component. With an increasing tilt of the magnetic easy axis, this factor reduces to approximately four at an easy axis angle aligned at 80 • relative to the long axis of the cantilever. However, the measurements and simulations of the current design confirm that the TMs exhibit a systematically smaller electromechanical response compared to BMs, which can significantly diminish the potential gain in sensitivity. Possible ways of improvement are sketched out. From simulated and measured resonance frequency curves f r (H) we found that the maximum normalized resonance frequency change ∆ f r and the magnetic sensitivity of TMs reduce with the increasing mode number due to the increasing contribution of C 11 to the resonance frequency. Hence, the dependency of TMs on the mode number is opposite to the one observed for BMs and caused by a different mechanism.
In conclusion, the delta-E effect of the C 66 component shows the promising potential of significantly increasing the magnetic sensitivity and the maximum normalized resonance frequency change ∆ f r . However, the efficient electrical excitation of TMs remains challenging for achieving high electrical sensitivity. Generally, the results imply that the simulation, a corresponding mass density of ρ FeCoSiB = 7700 kg/m 3 was used. For the stiffness tensor of the mechanically isotropic magnetic film, we use: For the magnetoelastic simulations, we use a saturation magnetic flux density of µ 0 M s = 1.5 T [17] and saturation magnetostriction of λ s = 35 ppm [17].

Appendix C.4. Piezoelectric Material (AlN)
For the stiffness matrix C AlN and the piezoelectric stress-charge coefficient tensor d we use values based on ab initio calculations [69]. Those tend to overestimate d and are here slightly adjusted. The stiffness tensor is: For the density, we use ρ AlN = 3300kg/m 3 and for the electrical permittivity ε el = ε r ε 0 with the electrical vacuum permittivity ε 0 and the relative electrical permittivity ε r , given by