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

**Benjamin Spetzler, Elizaveta V. Golubeva, Ron-Marco Friedrich, Sebastian Zabel, Christine Kirchhof, Dirk Meyners, Jeffrey McCord and Franz Faupel \***

> Institute of Materials Science, Faculty of Engineering, Kiel University, 24143 Kiel, Germany<EMAIL_ADDRESS>(B.S<EMAIL_ADDRESS>(E.V.G<EMAIL_ADDRESS>(R.-M.F<EMAIL_ADDRESS>(S.Z<EMAIL_ADDRESS>(C.K<EMAIL_ADDRESS>(D.M<EMAIL_ADDRESS>(J.M.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +49-431-880-6229

**Abstract:** 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 *C*66 component has the promising potential of significantly increasing the magnetic sensitivity and the maximum normalized frequency change Δ*f*r. 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 Δ*f*r of TMs on the mode number, which differs fundamentally from BMs and is well explained by our model. Because the contribution of *C*11 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.

**Keywords:** delta-E effect; magnetoelectric; magnetoelastic; resonator; torsion mode; bending mode; magnetic modeling; MEMS; FEM

#### **1. 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–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–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–18]. It was used for magnetic

**Citation:** Spetzler, B.; Golubeva, E.V.; Friedrich, R.-M.; Zabel, S.; Kirchhof, C.; Meyners, D.; McCord, J.; Faupel, F. Magnetoelastic Coupling and Delta-E Effect in Magnetoelectric Torsion Mode Resonators. *Sensors* **2021**, *21*, 2022. https://doi.org/10.3390/ s21062022

Academic Editor: Nicolò Marconato

Received: 10 February 2021 Accepted: 9 March 2021 Published: 12 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

field sensing with magnetoelectric plate resonators [19–22] and beam structures [23–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–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.

#### **2. 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 (Fe90Co10)78Si12B10 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.

**Figure 1.** 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.

#### **3. Sensitivity**

#### *3.1. 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–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]:

$$\left.S\_{\mathbf{Y},\mathbf{r}} \frac{\partial |Y|}{\partial f} \right|\_{f=f\_{\mathbf{r}}, \ H=H\_0} \cdot f\_{\mathbf{r}}; \quad \left.S\_{\Phi,\mathbf{r}} \frac{\partial \Phi}{\partial f} \right|\_{f=f\_{\mathbf{r}}, H=H\_0} \cdot f\_{\mathbf{r}}; \quad \left.S\_{\mathbf{H},\mathbf{r}} \frac{1}{f\_{\mathbf{r}}} \frac{\partial f\_{\mathbf{r}}}{\partial \mu\_0 \mathbf{H}} \right|\_{H=H\_0},\tag{1}$$

with the magnetic vacuum permeability *μ*0 ≈ 4*π*·10−7 N/A2. 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*.

#### *3.2. 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 *Cij* is a function of the stiffness tensor components *Cij*. Depending on the respective resonance mode, different *Cij* 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/*∂Cij* that contains the resonance properties of the structure and a purely magnetoelastic part *∂Cij*/*∂μ*0*H*:

$$\left. \left. S\_{\mathrm{H,r}} \frac{1}{f\_{\mathrm{r}}} \frac{\partial f\_{\mathrm{r}}}{\partial \mu\_{0} H} \right|\_{\mathrm{H} = \mathrm{H}\_{0}} = \sum\_{i=1}^{3} \sum\_{j=1}^{3} \frac{1}{f\_{\mathrm{r}}} \frac{\partial f\_{\mathrm{r}}}{\partial \mathrm{C}\_{ij}} \left. \frac{\partial \mathrm{C}\_{ij}}{\partial \mu\_{0} H} \right|\_{\mathrm{H} = \mathrm{H}\_{0}} \sum\_{i=1}^{3} \sum\_{j=1}^{3} \partial\_{\mathrm{C}} f\_{\mathrm{r},ij} \, \partial\_{\mathrm{H}} \mathrm{C}\_{ij}. \tag{2}$$

If treated separately, the factors *∂ f*r/*∂Cij* and *∂Cij*/*∂μ*0*H* must be normalized to remove the dependency on the absolute value of *Cij* that cancels out in *S*H,r. We define:

$$\left. \partial\_{\mathbf{C}} f\_{\mathbf{r},ij} \frac{\mathbf{C}\_{ij}}{f\_{\mathbf{r}}} \frac{\partial f\_{\mathbf{r}}}{\partial \mathbf{C}\_{ij}} \right|\_{H=H\_0}; \qquad \left. \partial\_{\mathbf{H}} \mathbf{C}\_{ij} \frac{1}{\mathbf{C}\_{ij}} \frac{\partial \mathbf{C}\_{ij}}{\partial \mu\_0 H} \right|\_{H=H\_0}. \tag{3}$$

From Equation (3), the factor *∂*C *f*r,*ij* represents a normalized measure for the influence of the stiffness tensor component *Cij* 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*Cij*, includes the delta-E effect and describes the normalized influence of the applied flux density *μ*0*H* on *Cij*. 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.

#### **4. 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.

#### *4.1. 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])

$$
\rho \frac{\partial^2 \overline{u}}{\partial t^2} = \overline{\nabla} \cdot \sigma \, , \tag{4}
$$

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]:

$$\begin{aligned} \overline{E} &= -\overline{\nabla} \mathcal{U}, \\ \overline{\nabla} \cdot \overline{D} &= \rho\_{\mathbb{C}} \end{aligned} \tag{5}$$

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]:

$$\begin{array}{l} \sigma = \mathcal{C}^\* \varepsilon - \mathfrak{e}\_{\text{c}}^T \overline{E} \\ \overline{D} = \mathfrak{e}\_{\text{c}} \varepsilon - \mathfrak{e}\_{\text{el}} \overline{E} \end{array} \tag{6}$$

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 *U*app = *U*0· exp(*i*[*ωt* + *ϕ*v]) is applied, 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* = *u*ˆ· exp(*i*[*ωt* + *ϕ*u]) and a solution for the electrical potential of *U* = *U*ˆ · exp(*i*[*ωt* + *ϕ*v]). 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.

**Figure 2.** Coordinate system used for the electromechanical and the magnetic model. All three components *mi* of the reduced magnetization vector *m* are described by the polar angle *θ* and the azimuthal angle *ϕ*.

#### *4.2. 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:

$$\mu = \mathcal{K}\left(1 - \left(m\_i E A\_i\right)^2\right) - \mu\_0 M\_\text{s} m\_i H\_i - \frac{1}{2} \mu\_0 M\_\text{s} m\_i H\_{\text{d},i} - \sigma\_j \lambda\_j \quad \text{with } i = 1, 2, 3 \quad j = 1, \dots, 6 \text{ } . \tag{7}$$

In this equation, the components of the reduced magnetization vector are denoted by *mi*, 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 *EAi* and the effective first-order uniaxial anisotropy energy density constant *K*. The components of the external magnetic field vector are given by *Hi* and the components of the mean demagnetizing field by *H*d,*i* = −*DiimiM*s, with the main diagonal components *Dii* 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 *mi*. The exact definition of all vector and tensor components is given in the appendix. The linearized change of the elastic compliance components *Sij* with the magnetic field and stress is derived from the expression

$$S\_{i\bar{j}}(H,\sigma) = \frac{\partial \varepsilon\_i}{\partial \sigma\_{\bar{j}}} = \frac{\partial (\varepsilon\_i + \lambda\_{\bar{i}})}{\partial \sigma\_{\bar{j}}} S\_{\mathbf{m},i\bar{j}} + \Delta S\_{i\bar{j}} \,. \tag{8}$$

where the first summand *S*m,*ij* is the constant, fixed magnetization elastic compliance tensor component. The magnetization dependent part Δ*Sij* can be obtained from the equilibrium conditions that are given by the first-order derivatives of *u*:

$$
\mu\_{\varphi} \frac{\partial \mu}{\partial \varphi} = 0 \quad \text{and} \; \mu\_{\theta} \frac{\partial \mu}{\partial \theta} = 0 \; . \tag{9}
$$

From these equilibrium conditions a general expression for the linearized change Δ*Sij* of the compliance tensor can be derived (Appendix A). Denoting the second-order derivatives as *uϕϕ* and *uθθ* it is:

$$
\Delta S\_{i\bar{j}} \frac{\partial \lambda\_i}{\partial \sigma\_{\bar{j}}} = -\frac{\partial \lambda\_i}{\partial \varrho} \frac{\partial u\_{\varrho}}{\partial \sigma\_{\bar{j}}} \frac{1}{u\_{\varrho\varphi}} - \frac{\partial \lambda\_i}{\partial \theta} \frac{\partial u\_{\theta}}{\partial \sigma\_{\bar{j}}} \frac{1}{u\_{\theta\theta}} \,. \tag{10}
$$

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:

$$
\Delta S\_{11} = \Delta S\_{22} = -\Delta S\_{12} = \frac{9\lambda\_s^2 \cos[\varphi]^2 \sin[\varphi]^2}{\mu\_{\varphi\varphi}},\tag{11}
$$

$$
\Delta S\_{16} = -\Delta S\_{26} = -\frac{9\lambda\_s^2 \cos[\varrho] \cos[2\varrho] \sin[\varrho]}{u\_{\varrho\varphi}}\,,\tag{12}
$$

$$
\Delta S\_{44} = \frac{9\lambda\_s^2 \sin[\varphi]^2}{\mu\_{\theta\theta}} \,, \tag{13}
$$

$$
\Delta S\_{45} = \frac{9\lambda\_s^2 \cos[\varphi] \sin[\varphi]}{u\_{\theta\theta}} \,, \tag{14}
$$

$$
\Delta S\_{55} = \frac{9\lambda\_s^2 \cos[\varphi]^2}{\mu\_{\theta\theta}} \,, \tag{15}
$$

$$
\Delta S\_{66} = \frac{9\lambda\_s^2 \cos[2\ \phi]^2}{\mu\_{\phi\Phi}} \,. \tag{16}
$$

The final compliance tensor for in-plane magnetization as a function of magnetic field and stress is:

$$S(H,\sigma) = \begin{bmatrix} S\_{11} & S\_{12} & S\_{m,13} & 0 & 0 & \Delta S\_{16} \\ S\_{12} & S\_{22} & S\_{m,23} & 0 & 0 & \Delta S\_{26} \\ S\_{m,13} & S\_{m,23} & S\_{m,33} & 0 & 0 & 0 \\ 0 & 0 & 0 & S\_{44} & \Delta S\_{45} & 0 \\ 0 & 0 & 0 & \Delta S\_{45} & S\_{55} & 0 \\ \Delta S\_{16} & \Delta S\_{26} & 0 & 0 & 0 & S\_{66} \end{bmatrix} \text{with } S\_{ij}(H,\sigma) = S\_{m,ij} + \Delta S\_{ij}. \tag{17}$$

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.

#### **5. Implications of the Magnetic Model**

In the following, results for the *Cij* 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(*H* → *H*K) → ∞ (Equation (4)) it is *C*66(*H* = *H*K) = 0. At |*H*| > |*H*K|, it is *C*66 < *C*m,66 with *C*66 = *C*m,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 *C*66(*ϕ* = 45◦, 135◦) = *C*m,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 *Sij* 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.

**Figure 3.** (**a**) Magnetic field dependent components *Cij* of the effective stiffness tensor for an ideal hard axis magnetization process of a macrospin. The external magnetic field with magnitude *H* is applied along the *x* axis and normalized the anisotropy field *H*K; (**b**) component *C*66 for different angles *ϕ*EA of the magnetic easy axis to the x axis; (**c**) maximum value *∂*H*Cij*,max of *∂*H*Cij*(*H*) (Equation (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 *Cij* components on the relative magnetic sensitivity *S*H,r using *∂*H*Cij* as defined in Equation (3). Calculating *∂*H*Cij* requires forming the derivate *∂Cij*/*∂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*Cij*(*H*) numerically from *Cij*(*H*) and extract the maximum *∂*H*Cij*,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*Cij*,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.

#### **6. Results**

#### *6.1. 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/m3 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.

**Figure 4.** (**a**) Magneto-optical Kerr effect microscopy image of the analyzed structure, demagnetized along the x axis and composed of a series of different images. The region of the magnetic layer is marked by a white square and the approximated orientation of the magnetic easy axis is 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**).

#### *6.2. 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 μ0H = − 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).

**Figure 5.** (**a**) Example comparison of measurement and FEM simulation of the sensor admittance around the first torsion mode (TM1), close to magnetic saturation at −10 mT; (**b**) comparison of measured and modeled maximum phase shift of the first three torsion modes (TM, samples 1–3) and the first and second bending modes (BMs, samples 4–9) with various electrode configurations. The BM measurements are published in Ref. [28].

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 μ0H = −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 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 *E*1 and *E*2, phase shifted by 180◦. Additionally, alternative piezoelectric materials with larger piezoelectric coefficients, such as AlScN [63–65] could increase the electrical sensitivity significantly and result in Δ*φ* comparable with bending modes.

#### *6.3. 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.

**Figure 6.** Measurements of the resonance frequency *f*r as a function of the applied magnetic flux density *μ*0*H* along the long cantilever axis (x axis) starting at *μ*0*H* = −10 mT: (**a**) normalized resonance frequency *f*r/ *f*r,max of the first three torsion modes (TMs) (left) and the first three transversal bending modes (BMs) (right). The maximum resonance frequencies *f*r,max of the TMs are 26.256, 87.478, 175.150 kHz, and of the BMs: 7.649, 47.182, 121.400 kHz; (**b**) relative magnetic sensitivities *S*H,r = *S*H/ *f*r,max calculated from the data in (**a**) with Equation (1).

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*Cij*, 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.

#### *6.4. 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/m3, 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 *Cij* 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).

**Figure 7.** (**a**) Simulated resonance frequency *f*r normalized to its maximum, as a function of the magnetic bias flux density *μ*0*H* for the first three torsional modes (TMs); and (**b**) the first three bending modes (BMs); and (**c**) the magnitude of the displacement vector of the first three torsion modes plotted and calculated with COMSOL(*r*) Multiphysics v. 5.3a (COMSOL AB, Stockholm, Sweden).

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.

#### **7. 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

delta-E effect of different *Cij* can have opposite effects on Δ*f*r, depending on the resonance mode. This was demonstrated in the example of torsion modes. Because the contribution of *C*11 increases with the torsion mode number, the first-order torsion mode shows the highest magnetic sensitivity. In addition to fundamental insights on the delta-E effect in higher resonance modes, a model for the electrical and the magnetic sensitivity was presented. The results are not only relevant for the development of magnetoelastic magnetic field sensors, but also for frequency tunable devices based on the delta-E effect.

**Author Contributions:** Conceptualization, B.S.; methodology, B.S., E.V.G., R.-M.F.; validation, B.S.; formal analysis, B.S., E.V.G.; investigation, B.S., E.V.G., S.Z., C.K.; writing—original draft preparation, B.S.; writing—review and editing, B.S., E.V.G., R.-M.F., S.Z., C.K., D.M., J.M., F.F.; visualization, B.S., E.V.G.; supervision, D.M., J.M., F.F.; project administration, D.M., J.M., F.F.; funding acquisition, D.M., J.M., F.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Research Foundation (DFG) via the collaborative research center CRC 1261.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors thank Christoph Elis for performing the density measurements on FeCoSiB.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **Appendix A. Magnetoelastic Model**

#### *Appendix A.1. Definition of Vectors*

In the following, a detailed definition of all vectors is given in the spherical coordinate system. All polar angles are denoted by *θ* and all azimuthal angles by *ϕ*, with an index if it is not the angles of the reduced magnetization. The reduced magnetization vector is denoted as

$$
\overline{m} = \begin{bmatrix} \cos\varphi \sin\theta & \sin\varphi \sin\theta & \cos\theta \end{bmatrix}^T,\tag{A1}
$$

The easy axis unit vector *EA* is given by

$$\overline{EA} = \begin{bmatrix} \cos\varphi\_{\text{EA}}\sin\theta\_{\text{EA}} & \sin\varphi\_{\text{EA}}\sin\theta\_{\text{EA}} & \cos\theta\_{\text{EA}} \end{bmatrix}^T. \tag{A2}$$

The vector *H* of the external applied field and *H*d of the demagnetizing field are given by

$$\overrightarrow{H} = H \begin{bmatrix} \cos\varphi\_{\text{H}}\sin\theta\_{\text{H}} & \sin\varphi\_{\text{H}}\sin\theta\_{\text{H}} & \cos\theta\_{\text{H}} \end{bmatrix}^{T} \tag{A.3}$$

$$\overline{H}\_{\rm d} = -M\_{\rm S} \begin{bmatrix} D\_{11}m\_1 & D\_{22}m\_2 & D\_{33}m\_3 \end{bmatrix}^T. \tag{A4}$$

For all higher-order tensors, we use the Voigt notation. The magnetostriction tensor is then given by

$$
\lambda = \frac{3}{2} \lambda\_5 \begin{bmatrix} m\_1^2 - \frac{1}{3} & m\_2^2 - \frac{1}{3} & m\_3^2 - \frac{1}{3} & 2m\_2 m\_3 & 2m\_1 m\_3 & 2m\_1 m\_2 \end{bmatrix}^T. \tag{A5}
$$

*Appendix A.2. General Expression for* Δ*Sij*

From the equilibrium conditions, one can write

$$\begin{array}{l} \Delta S\_{i\bar{j}} = \frac{\partial \lambda\_{i}}{\partial \sigma\_{\bar{j}}} = \frac{\partial \lambda\_{i}}{\partial \varphi} \frac{\partial \varphi}{\partial \sigma\_{\bar{j}}} + \frac{\partial \lambda\_{i}}{\partial \theta} \frac{\partial \theta}{\partial \sigma\_{\bar{j}}} = \frac{\partial \lambda\_{i}}{\partial \varphi} \left( -\frac{\partial u\_{\varphi}}{\partial \sigma\_{\bar{j}}} / \frac{\partial u\_{\theta}}{\partial \varphi} \right) + \frac{\partial \lambda\_{i}}{\partial \theta} \left( -\frac{\partial u\_{\theta}}{\partial \sigma\_{\bar{j}}} / \frac{\partial u\_{\theta}}{\partial \theta} \right) \\\ -\frac{\partial \lambda\_{i}}{\partial \varphi} \frac{\partial u\_{\varphi}}{\partial \sigma\_{\bar{j}}} \frac{1}{u\_{\varphi\varphi}} - \frac{\partial \lambda\_{i}}{\partial \theta} \frac{\partial u\_{\theta}}{\partial \sigma\_{\bar{j}}} \frac{1}{u\_{\theta\theta}} . \end{array} \tag{A6}$$

*Appendix A.3. Derivatives of the Energy Density Functional*

For in-plane magnetization (*θ* = *π*/2), easy axis (*θ*EA = *π*/2), and magnetic field (*θ*H = *π*/2), the second-order derivatives of *u* (Equation (7)) are given by

$$\begin{split} \mu\_{\theta\theta} \frac{\partial^2 u}{\partial \rho^2} &= 3\lambda\_s [\sigma\_{11} \cos(2\rho) - \sigma\_{22} \cos(2\rho) + 2\sigma\_{12} \sin(2\rho)] + \mu\_0 M\_s H \cos(\rho - \varphi\_{\text{H}}) + 2K \cos(2[\varphi - \varphi\_{\text{E}\text{A}}]) \\ &+ \mu\_0 M\_s^2 (D\_{22} - D\_{11}) \cos(2\rho) \end{split} \tag{A7}$$

*uθθ ∂*2*u ∂θ*2 = 3*λ*s *σ*11 cos2 *ϕ* + *σ*22 sin2 *ϕ* − *σ*33 + *σ*12 sin(2*ϕ*) + *μ*0*M*s*H* cos(*ϕ* − *ϕ*H) + 2*K* cos2(*ϕ* − *ϕ*EA) −*μ*0*M*2 s *D*11 cos2 *ϕ* + *D*22 sin2 *ϕ* − *D*33 (A8)

#### **Appendix B. Resonance Frequencies and Sensitivities**

A summary of the maximum measured resonance frequencies *f*r,max at *μ*0*H* = −10 mT is given in Table A1, together with the corresponding quality factor *Q*max and the maximum magnetic sensitivities from Figure 6. The bending modes (BMs) and the torsion modes (TMs) were all measured using only electrode E1. Hence, the electrodes are not optimized for the bending modes. Additionally, the maximum quality factor of the BMs is a factor of three smaller than in TM1. Due to both factors, the electrical sensitivity of BMs is significantly smaller than that of TMs for our cantilever.

**Table A1.** Measured resonance frequencies *f*r,max of the six modes analyzed and the quality factor *Q*max, both measured in magnetic saturation at *μ*0*H* = −10 mT. The maximum magnetic sensitivity *S*H and the maximum relative magnetic sensitivity *S*H,r are obtained from Figure 6b. The maximum relative electrical sensitivities found are given by *Sφ*,r and *S*Y,r.

#### **Appendix C. Geometry and Material Parameters**

#### *Appendix C.1. Geometry*

The poly-Si cantilever was measured with an optical microscope to be *L* = 3.12 mm long and *W* = 2.15 mm wide. For the simulations, the length in the model was slightly adjusted within the measurement accuracy to 3.116 mm. The magnetostrictive layer was deposited directly at the clamping on the bottom side of the poly-Si cantilever and has a width of *W*mag = 2 mm and a length of *L*mag = 3.05 mm. The AlN layer is of the same geometry but deposited on top of the poly-Si. The electrodes *E*1 and *E*2 are positioned at the clamping on top of the AlN layer and the left and right edge, respectively. An additional parallel capacitance of *C*0 = 17.7 pF is used. It is consistent with the area of the bond pads, the conduction lines, and the relative electrical permittivity used for the AlN.

#### *Appendix C.2. Substrate (Poly-Si)*

For the poly-silicon substrate, we use isotropic material parameters, with Young's modulus *E*Si = 160 GPa [66,67] a Poisson's ratio of *v*Si = 0.22 [67] and a mass density of *ρ*Si = 2300 kg/m3 [68].

#### *Appendix C.3. Magnetic Material (FeCoSiB)*

The mass density of FeCoSiB was experimentally determined by estimating the volume with profilometer measurements and the mass with a microbalance. The measurements were performed on a 6-inch wafer with a mean FeCoSiB layer thickness of approximately 1.5 μm. A density of *ρ*FeCoSiB = (7870 ± 1350) kg/m3 was obtained. For the

simulation, a corresponding mass density of *ρ*FeCoSiB = 7700 kg/m3 was used. For the stiffness tensor of the mechanically isotropic magnetic film, we use:

$$
\mathbf{C}\_{\mathbf{m},ij} = \begin{bmatrix}
\mathbf{C}\_{\mathbf{m},11} & \mathbf{C}\_{\mathbf{m},12} & \mathbf{C}\_{\mathbf{m},12} & 0 & 0 & 0 \\
\mathbf{C}\_{\mathbf{m},12} & \mathbf{C}\_{\mathbf{m},11} & \mathbf{C}\_{\mathbf{m},12} & 0 & 0 & 0 \\
\mathbf{C}\_{\mathbf{m},12} & \mathbf{C}\_{\mathbf{m},12} & \mathbf{C}\_{\mathbf{m},11} & 0 & 0 & 0 \\
0 & 0 & 0 & \mathbf{C}\_{\mathbf{m},44} & 0 & 0 \\
0 & 0 & 0 & 0 & \mathbf{C}\_{\mathbf{m},55} & 0 \\
0 & 0 & 0 & 0 & 0 & \mathbf{C}\_{\mathbf{m},66}
\end{bmatrix} \tag{A9}
$$

Using a Young's modulus of *E* = 150 GPa and a Poisson's ratio of *v* = 0.3, both at fixed magnetization, the non-zero components of the fixed magnetization stiffness tensor are:

$$\mathcal{C}\_{\mathfrak{m},11} = \mathcal{C}\_{\mathfrak{m},22} = \mathcal{C}\_{\mathfrak{m},33} = \frac{E(1-v)}{(1+v)(1-2v)} = 201.92 \text{ GPa} \tag{A10}$$

$$\mathbb{C}\_{\mathfrak{m},12} = \mathbb{C}\_{\mathfrak{m},13} = \mathbb{C}\_{\mathfrak{m},23} = \frac{E(v)}{(1+v)(1-2v)} = 86.54 \,\text{GPa} \tag{A11}$$

$$\mathcal{C}\_{\text{m},66} = \mathcal{C}\_{\text{m},55} = \mathcal{C}\_{\text{m},44} = \frac{E}{2(1+v)} = 57.69 \text{ GPa} \tag{A12}$$

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:

$$\mathbf{C}\_{\text{AlN}} = \begin{bmatrix} 410.2 & 142.2 & 110.1 & 0 & 0 & 0\\ 142.2 & 410.2 & 110.1 & 0 & 0 & 0\\ 110.1 & 110.1 & 385.0 & 0 & 0 & 0\\ 0 & 0 & 0 & 122.9 & 0 & 0\\ 0 & 0 & 0 & 0 & 122.9 & 0\\ 0 & 0 & 0 & 0 & 0 & 134.0 \end{bmatrix} \text{GPa} \tag{A13}$$

The piezoelectric coupling tensor *c*e = *dC*AlN results as

$$\mathbf{c}\_{\mathbf{e}} = \begin{bmatrix} 0 & 0 & 0 & 0 & -0.27828 & 0 \\ 0 & 0 & 0 & -0.27828 & 0 & 0 \\ -0.4496 & -0.4496 & 1.41 & 0 & 0 & 0 \end{bmatrix} \frac{\mathbf{p}\mathbf{C}}{\mathbf{m}^2}.\tag{A14}$$

For the density, we use *ρ*AlN = 3300kg/m3 and for the electrical permittivity *ε*el = *ε*r*ε*0 with the electrical vacuum permittivity *ε*0 and the relative electrical permittivity *ε*r, given by

$$
\mathfrak{e}\_{\mathbf{t}} = \begin{bmatrix}
9.2081 & 0 & 0 \\
0 & 9.2081 & 0 \\
0 & 0 & 10.1192
\end{bmatrix} \tag{A15}
$$

#### **References**


**2. Properties of biocomposites**

The most important characteristic feature of selection materials for various applications is depending on its properties. The properties of materials are often dependent on the isotropic and anisotropic nature of the materials. The properties of materials that relate to different physical phenomena often behave linearly (or approximately so) in a given operating range. Modeling them as linear can significantly simplify the differential constitutive equations that the property describes. On the other hand, the relevant equations are also used to determine the material properties. If we know the original length of a material, then we can determine the gain or loss of its original length by calculating change of the length. Material properties are most reliably measured by standardized test methods. Many such test methods have been documented by their respective user communities and published through ASTM International [6]. We noted the most important properties of the renewable biocomposites from several researchers' investigations that are listed below.

(1) Mechanical properties, (2) thermal properties, (3) optical properties, (4) degradable properties, and (5) electrical properties.


**1. Introduction**

With a new climate system that replaces the Kyoto Protocol, which expired in 2020, the Paris Climate Agreement has been applied since January 2021, which includes strengthening global warming control targets and reducing greenhouse gas emissions. As follow, various policies have been enacted to encourage the development of new renewable energy sources to reduce greenhouse gas emissions.

In addition to international treaties, South Korea is promoting a "Renewable Energy 3020" policy to increase the ratio of renewable-energy power generation by 2030. The goal of the policy is to increase the renewable energy power generation ratio by 20% or more, and to supply over 95% of the capacity of the newly installed platform with clean energy such as solar and wind power [1]. Since 2012, the Renewable Portfolio Standard (RPS), which mandates a certain percentage of power generation as new renewable energy for industrialization such as solar power, wind power, and hydrogen, has also been enforced [2].

Recently, in the case of solar power generation, the demand for installation on water (e.g., lakes, coasts, and on the sea) is increasing rapidly due to the lack of installation area on land. It is more convenient to secure space on water and in the sea than on land, minimizes the effect of natural disasters, and eases the development large-scale power generation complexes [3].A floating solar power farm system is frequently used as a unified platform considering its installation convenience, expandability, and mobility.

On the other hand, since the coast and sea are exposed directly to natural disturbance (e.g., seawater, wind, and waves), technologies that can minimize the effects of waves and actively reduce fluctuations are needed [4–6]. Owing to the periodic motion of waves, fatigue loads are generated on structures (e.g., floating bodies, supporting frames, panels,

**Citation:** Lee, J.-H.; Paik, K.-J.; Lee, S.-H.; Hwangbo, J.; Ha, T.-H. Experimental and Numerical Study on the Characteristics of Motion and Load for a Floating Solar Power Farm under Regular Waves. *J. Mar. Sci. Eng.* **2022**, *10*, 565. https://doi.org/ 10.3390/jmse10050565

Academic Editor: Eugen Rusu

Received: 31 March 2022 Accepted: 19 April 2022 Published: 21 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and coupling devices), and a design is needed for supporting structures considering structural stability. To consider the stability of a floating solar power generation system, experiments and numerical studies on the motion performance and connections between adjacent platforms are still underway.

Oliveira-Pinto et al. [7] provided a literature review on the applicability of floating solar power farms and introduced the technologies currently available, as well as the technical challenges and risk factors when designing studies in the marine environment.

Sahu et al. [8] illustrated the concept and advantages of floating solar power systems and studied each component of solar power systems. The environmental loads can cause deformation and stress in the module, causing microcracks and reducing productivity and durability. The technology of thin-film research was proposed to withstand harsh environmental conditions. It was noted that remote sensing and GIS (geographic information system)-based technologies can be utilized to determine the potential of solar panels.

The floating structure uses individual HDPE (high-density polyethylene), and the floating bodies were connected to a pin. It was said that the connection part was weak due to environmental factors. When this connector is exposed to high wave heights, stress concentration could occur and cause system problems. It was shown that the behavior of floating systems is complex [9,10].

Shi et al. [11] introduced a network modeling method for the dynamic prediction of multimodule floating structures and conducted it experimentally in the wave. The connector was made on three-axis motion, and RAOs (response amplitude operator) were compared at various frequencies through the experiment. Recently, the SCOTRA company developed a connector system that allows floating bodies to have their own motility in all directions.

When waves interact with a floating solar power farm (FSP), the motions of the FSP generate inertial forces and dynamic loads on the structures. Non-linear effects make it difficult to predict what will happen on the platforms. It could increase the motion response complexity exponentially of floating bodies. Therefore, recent studies are trying to predict the motion of the system using the using the computational fluid dynamics (CFD) [12,13].

All previous studies have stated that it is important to predict the load at the design stage due to the non-linear environmental factors that affect the FSP depending on the wave conditions. Also, it is necessary to consider motion response with different FSP components (e.g., connecting method, the type of mooring system, the array of floating bodies, and the position and angle of the solar panel).

In this study, we focused on the motion response and tried to understand the relative motion that may occur in a unit platform. Owing to the establishment of numerical techniques, the simulations were performed for a floating solar power farm in the form of a unit under wave conditions. Uniaxial hinge and catenary techniques were applied to the connector and mooring system to determine the motion performance of the floating solar power farm for six degrees of freedom (6DOF). This was verified by an experiment conducted in the Inha University towing tank (IUTT). It investigated fluid forces such as load and pressure distribution, which are difficult to determine using the experiment.


**4. Discussion**

In the summer of 2017, southern Europe and the Euro-Mediterranean were hit by an exceptional heat wave [6,42]. After an outstandingly warm June in western Europe, the heat returned to southern Italy in July. It contributed to more than 400 wildfires, which destroyed approximately 800 km2 of forest and vegetated areas. The number of fires has been unprecedented in the last 20 years. Furthermore, early August saw a particularly intense heat wave described as the "worst heat wave since 2003", with the air temperature above 40 ◦C in many parts of Italy [42].

In the quest for possible satellite indices to assess and possibly mitigate the effect of long-lasting drought on vegetation, we have devised an index, *wdi*, which takes advantage of the IASI capability to retrieve surface data simultaneously with atmospheric parameters. Other methods that use satellite data exploit the visible region of the electromagnetic spectrum (e.g., NDVI, NDMI, and related indices) or the microwave band (e.g., *ssm* and LAI). The normalized difference moisture index (NDMI) (e.g., [43]) is mainly intended to detect humidity in vegetation using a combination of near-infrared (NIR) and shortwave infrared (SWIR) spectral bands. The index NDMI and the original greenness index, NDVI, with the same *ssm*, have also been used coupled to surface temperature and air temperature (e.g., see [8,11–14,18]). Other tools have tried to couple surface temperature and the humidity field, e.g., [6].

In contrast, our *wdi* exploits the thermal band of the Earth's emission spectrum and simultaneously uses the surface temperature, air temperature, and humidity. To our knowledge, this combination is unique. In effect, the water deficit index we have defined can monitor water deficit and assess vegetation stress, as the comparison with in situ measurements has demonstrated. It can be used complementary to *ssm*, LAI, and the set of NDVI-related indices to better understand the intensity and danger of heatwaves for the vegetation. Sequences of increasing *wdi* can help to identify the onset of water deficit for the vegetation, hence the increased risk of fire, especially in forests.

The *wdi* index is meant to identify regions where particular weather conditions can produce water deficits. The index is not intended as an estimate or an estimator of evapotranspiration. This process is also affected by vegetation/crop characteristics, environmental conditions, and cultivation types. Therefore, there is too much variability, which cannot be condensed into a single index. The *wdi* parameter is a bulk index, which can help to monitor forest and wood regions suffering from long-lasting droughts because of adverse weather conditions. It can be mapped on a regional and even global scale, allowing us to monitor drought processes at a glance. The *wdi* maps could be important to monitor and evaluate the risk of fires in the large forested area, which is otherwise inaccessible. In addition, we have shown that in regions where the vegetal ecosystem has a particular fragility to water deficit, the index can soon quantify the possible danger and require more accurate in situ observations.

In this respect, *wdi* is most effective in the case of a heatwave. In the wintertime, for example, large values of *wdi* could be linked to a dry atmosphere and low air temperature. In effect, this is the case in January 2017 for the more southern area on the map of Figure 5, which belongs to the high mountains of the Sila chain. Additionally, in summer, very humid and warm conditions could lead to *wdi* ∼ 0. For example, this is the case for the coastal regions in July–August 2017, as seen again in the map in Figure 5. For these cases, it is better to look separately at the maps of *Ts* and *Td*. In this respect, we observe that the dew point temperature has been individuated as a key parameter to compute sophisticated indicators of health stress for human beings during heatwaves [44].

Some words of caution should also be said about the temporal sampling of *wdi*. The occasional occurrence of a high *wdi* for one day should be of no concern. Drought is a process that takes several days or months. The severity of the process depends on its time continuity and persistence. Therefore, it is crucial to assess the persistence of the process, which can be done by looking at the time series. Averaging over several days can help to understand the persistence of the phenomenon. Another important point concerns the

capability of the retrieval system to solve the daily cycle, which cannot be done with the present polar satellite IASI instrument. During the night, the surface temperature could go below the dew point temperature and cause water vapour to condense at the surface. Therefore, it could be interesting to examine day and night separately. Hopefully, this could be the case when the MTG-IRS (https://www.eumetsat.int/mtg-infrared-sounder (accessed on 15 August 2022)) is put in orbit.