A Combined Magnetoelectric Sensor and Human Head Model for Biomagnetic FEM Simulations

We investigate a coupled multiscale 3-D finite-element model consisting of an orthogonal array of thin-film cantilever magnetoelectric (ME) magnetic field sensors and a simplified human head model. Electric point dipole sources are placed inside the head to generate an electromagnetic field. This field propagates through the tissue layers and outside of the head, where it reaches the sensor array. The investigated sensors are based on a 300-<inline-formula> <tex-math notation="LaTeX">$\mu \text{m}$ </tex-math></inline-formula>-thick silicon substrate layer of 26.25 mm length and 2.45 mm width, with a 20-<inline-formula> <tex-math notation="LaTeX">$\mu \text{m}$ </tex-math></inline-formula>-thick aluminum nitride (AlN) layer and a 20-<inline-formula> <tex-math notation="LaTeX">$\mu \text{m}$ </tex-math></inline-formula>-thick FeCoSiB layer as piezoelectric (PE) and magnetostrictive materials, respectively, located on opposite sides of the substrate. We position three sensors orthogonally to obtain a vector field sensor. The head model is based on the three-shell approach and consists of concentric spheres representing white matter, skull, and skin. All three layers are assigned specific conductivity and relative permittivity values from the literature, aiming to approximate the propagation of electromagnetic fields through different tissue types. We observe the propagation of an electric field generated by the dipole source and subsequent induction of a magnetic field inside the head structure, propagating outwards to the sensor array, where the generated voltage in the PE layer is evaluated for different ME sensor array and source positions and orientations. We show different behaviors of the sensors for a dipole source inside the head tissue versus a simple air environment, highlighting the benefits in accuracy and specificity of a combined head and sensor model with realistic material parameters.


I. INTRODUCTION
M EASURING biomagnetic signals, such as in mag- netocardiography (MCG), magnetoencephalography (MEG), or magnetomyography (MMG) applications, has become increasingly important in medical diagnostics and patient care.Researchers have demonstrated several viable approaches to measure biomagnetic fields, such as fluxgate sensors [1], [2], optically pumped magnetometers (OPMs) [3], or superconducting quantum interference device (SQUID)based sensor systems [4].While each of these technologies has its own unique operation principle, each one also has specific drawbacks, such as the limitation of SQUID-based systems to require cooling with liquid helium or limitations of OPM sensors to be sensitive to ambient background magnetic fields.In this work, we consider strain-coupled composite magnetoelectric (ME) magnetic field sensors as a promising technology for conducting such measurements in clinical environments.As opposed to existing measurement systems such as OPMs, fluxgate sensors, or SQUID-based sensor systems, ME sensors offer highly sensitive measurements of biomagnetic fields in the pT/Hz 1/2 range, while enabling passive, uncooled, and unshielded operation [5], [6].Recent investigations characterized resonant ME sensor systems with a signal-to-noise ratio of −90 dB, an application-specific capacity of 9.8 × 10 −7 dB•Hz, and a limit-of-detection of 22 pT, while the inherent noise for the sensor itself (i.e., without readout electronics) is given as ≈4 pT/ √ Hz at 7.684 kHz at room temperature [6].Research into ME devices also includes studies with devices such as novel ME antennas that can be used for energy harvesting and magnetic field sensing purposes simultaneously, including energy harvesting experiments on mice brain tissue and simu-lation of the magnetic field and tissue interactions [7].Experts in the mathematical and medical fields have been studying the relationship between electric sources inside the human brain and the resulting electric potential on the scalp for many years.This article has led to numerous important findings and the development of highly detailed and precise so-called electroencephalography (EEG) forward models.Publications in this field cover many regions of interest, such as the efficient computation of EEG forward models [8], the effects of forward model errors on EEG source localization [9], or improved forward head models for specific applications such as ear-EEG [10].Although this work uses a simplified three-shell approach for the head geometry, which is already well-known in the literature regarding the EEG forward problem, the novel approach we present here is to include a fully coupled ME sensor model in combination with the simplified head model.With this strategy, we can focus on the behavior of ME sensors in the presence of a head model with reduced complexity that is less computationally taxing and offers readily available literature for validation, while the combined model retains the potential to be increased in complexity for further investigations, e.g., with realistic MRI data.While, in past studies, the performance characterization of ME sensors has often been conducted under simplified homogeneous magnetic field excitation, this work focuses on investigating sensor behavior under the influence of inhomogeneous magnetic fields in order to more accurately model sensor behavior in application-oriented scenarios.Our theoretical investigation is based on the approach of 3-D finite-element method (FEM) simulations and introduces a combined simplified three-shell human head and an array of three orthogonal ME sensors in a single FEM model [11], [12], [13].The ME sensor array is positioned inside an air volume surrounding the human head, which in turn contains one or multiple electric point dipole sources.Depending on the position and orientation of the point dipole inside the head, a different magnetic field propagation through the head tissue and the surrounding air volume is observed and propagates outwards and into the ME sensors.Previous studies show strong directionality of the sensor response depending on incident angle of the magnetic field relative to the cantilever's long axis [14], [15], [16].By using this approach, we can calculate the inhomogeneous magnetic flux reaching the sensors and the subsequent electric sensor response based on the interaction of one or multiple electric sources inside of the biological tissue, the propagation of fields through the tissue, and surrounding air environment.This novel approach of a combined model offers new insights into the behavior of ME sensors, but it also requires a careful modeling setup due to the difficulties of combining very small structures, such as the sensor's 20-µm-thin layers, with comparatively large geometries such as the head.We will describe this work in detail in the following sections.In Section II, we will give a detailed overview of the methodology used in our FEM simulation models.This section will first illustrate the ME sensor model and the head model separately and then describe the combined sensor and head model.Relevant material and geometry parameters, as well as equations for the FEM calculations, are also given.Section III covers the results of our investigation.These include the electric and magnetic fields created inside the head domain and the surrounding air environment, as well as the sensor response to the propagation of the magnetic field outwards from the head into the sensor array.We conclude our findings in Section IV, emphasizing the benefits of using a combined ME sensor and head model with realistic material parameters.

II. METHODOLOGY
We implement the joint model and calculate the sensor behavior with 3-D FEM simulations utilizing the commercial software COMSOL Multiphysics 6.1 including the physics packages solid mechanics, magnetic fields, and electrostatics.We perform the simulations with the frequency-domain study and an excitation frequency of 848 Hz, corresponding to the cantilever sensor's previously established resonance frequency for operation in direct detection mode [13].Sections II-A-II-C will briefly describe the modeling process from the separate models used for the ME sensors and the human head, leading to the setup for our combined sensor and head model.The complexity of this combined model gave rise to some difficulties that we encountered during the stages of our investigation.We will also use this section to describe the problems that these difficulties lead to, as well as the simplifications and solution strategies that were necessary to achieve reliable results.

A. ME Sensor Model
Each of the three modeled ME sensors consists of three materials: the substrate layer, magnetostrictive (MS) layer, and piezoelectric (PE) layer.These layers form a composite with slipless and ideal interfaces.The coupled system is described in detail in previous works based on the constitutive equations given in [17], [18], [19], [20], [21], and [22], while this work utilizes the built-in physics interfaces of COMSOL Multiphysics.The solid mechanics interface governs the mechanical behavior of our model and is used to solve the equations of motions together with a constitutive model for solid materials.The following equations define the mechanical behavior of our model: In the linear elastic material node, S and ε stand for the stress and strain, respectively, u is the displacement vector, ω is the angular frequency, and ρ is the density.C describes the coupling between stress and strain for the silicon substrate and is a function of its Young modulus and Poisson's ratio.
The MS material node contains the relation between stress, strain, magnetic field, and magnetic flux density in the stressmagnetization form Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Here, c H and e HS are the elasticity matrix and coupling matrix (Voigt notation), respectively, B and H are the magnetic field and magnetic flux density, respectively, and µ rS is the relative permeability.The PE material node of our model includes equations determining the coupling between the mechanical and electric properties of the sensors D is the electric displacement field and ρ v is the volume charge density according to Gauss's law.c E and e ES denote the elasticity matrix and coupling matrix (Voigt notation) in the stress-charge form, respectively.E and ε rS represent the electric field and the relative permittivity, respectively.In conjunction with solid mechanics described in (1)-( 9), our model utilizes the interfaces of magnetic fields and electrostatics.While the first interface uses Maxwell's equations to primarily solve for the vector potential with ( 10)-( 12), the latter utilizes Gauss' law and solves for the electric potential related to ( 13) and ( 14) Here, σ refers to the specific conductivity, A is the magnetic vector potential, V is the electric potential, and P is the polarization.Equations ( 10)-( 14) combined with the mechanical interrelations in (1)-( 9) allow us to simulate MS and PE coupling and enable us to model the physical effects that govern the functioning principle of ME sensors.Fig. 1(a) shows a schematic of a single ME sensor.The cantilever structure's substrate layer has a total length of 26.25 mm, out of which 3 mm consists of a fixed structure to enable resonant bending modes along the long axis in a fixed-free setup.The width of the substrate layer is 2.45 mm and its thickness is 0.3 mm.The red layer on the top of the substrate illustrates the PE layer with a length of 25.6 mm, width of 1.6 mm, and thickness of 20 µm.On the opposite side of the substrate, the MS layer is located with a length of 22.9 mm, width of 1.8 mm, and thickness of 20 µm.Fig. 1(b) illustrates the array consisting of three orthogonal ME sensors.Fig. 1(c) shows a cross section of the sensor with the layer thicknesses to scale.The geometry of this sensor is based on a fabricated ME sensor, which was computationally and experimentally characterized in [13].Lastly, we illustrate moving from previous studies based on homogeneous magnetic excitation fields in Fig. 1(d) toward ME sensors in inhomogeneous magnetic fields, such as those propagating outwards from a source inside the human head depicted in Fig. 1(e).Due to the strong directionality of the sensor and the definition of the materials parameters in our setup, the elasticity and coupling tensors for the PE and MS layers must be rotated according to the orientation of the respective cantilever in the coordinate system.We perform this tensor rotation within COMSOL Multiphysics by defining additional base vector systems and assigning the MS and PE properties of each sensor to the corresponding coordinate system, based on the sensor's orientation in the modeling space.
The material parameters for the nonrotated base vector system (which corresponds to the sensor oriented in x-direction) are given in the Appendix and used to implement and compute the sensor behavior.The dimensions and material types of the sensors are given in Table I.The material used for the substrate layer is silicon.The PE layer consists of aluminum nitride (AlN) [17].Lastly, the MS material of choice is the compound material FeCoSiB, commercially known as Metglas [17].The sensor model does not consider any additional layers such as electrode layers.

B. Human Head Model
The implemented simplified head for this study is a basic three-shell head model [11], [12].It consists of three concentric spheres of different radii representing the tissue types of the head: skin, skull, and white matter (from outermost to innermost sphere).Fig. 2 illustrates the layer structure of the head in a cross section of the xy plane.The three tissue types with their radii 88, 92, and 100 mm are assigned material parameters for their relative permittivity and specific conductivity according to the literature [23], [24], [25], with our white matter layer encompassing the geometry for separate cortex and fluid layers from other related investigations [23].The relevant geometry parameters and electric properties of the head are given in Table II.

C. Combined Model
Following the individual design of the sensor model and head model, we combine both geometries in a joint model.The center of the head is placed in the origin of the coordinate Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.system, and the sensor array can be positioned at arbitrary locations around the head.This setup allows us to investigate the behavior of each of the three orthogonal sensors for different positions in space, relative to the source inside the head.Fig. 3(a) shows a cross section of the head model with its three tissue layers skin, skull, and white matter, as well as the ME sensor and an electric point dipole source located in the white matter domain.Fig. 3(a) illustrates the method of positioning the ME sensor in proximity of the head and calculating its response for various positions inside the air surrounding, e.g., performing a position sweep in the 3-D space around the head.The electric point dipole can be positioned at arbitrary locations inside the head geometry to study the effects of its position in the material frame on the propagation of electric and magnetic fields inside and outside of the head.It is also possible to have multiple active dipole sources, and the model is compatible with different source configurations altogether, e.g., deep brain stimulation (DBS) electrodes.
One of the challenges in combining both the sensor and head in a joint 3-D FEM model is that different geometries can have vastly different aspect ratios.While the simplified human head typically has a circumference in the cm-range, depending on the tissue compartment, the ME sensor's material layers have thicknesses in the µm-range, i.e., up to four orders of magnitude smaller.In order to calculate both relatively large and substantially smaller geometries in the same model, the 3-D mesh must be designed carefully.Therefore, the mesh elements constituting the head geometry are chosen to be small enough to achieve sufficiently accurate results, while avoiding sizes so small that they would lead to excessive computational load.Fig. 4(a)-(c) shows the mesh setup for the head, airbox, and ME sensor with their differently sized and structured tetrahedral mesh elements, respectively.The mesh is divided Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.into three separate areas, with one area representing the head, one for the sensor array, and one for the air environment.The mesh for the sensor array has an extremely fine element size, where the head and airbox geometries have comparatively smaller amounts of mesh elements which are also larger in size and include a higher maximum element growth rate and slightly smaller resolution of narrow regions.With this setup, the volume between the head and sensor array can be very finely meshed in order to enable a progression from the head geometry to the much finer sensor geometry, while the sections of less interest further away from the head and sensor array can be meshed with much larger elements.The mesh always extends between the head and ME sensor array continuously through the surrounding air environment.Over the course of our mesh optimization, however, we achieved the best results using a physics-controlled mesh that was automatically generated by COMSOL Multiphysics and slightly adjusting mesh parameters maximum element size, minimum element size (i.e., longest edge of tetrahedral mesh element), and maximum element growth rate for specific regions.This way, we were able to optimize the mesh settings while keeping the computational load as low as possible.The number of degrees of freedom solved for is ∼19.4 × 10 6 .Table III gives an overview of the mesh parameters used for different regions of the model.The surrounding air environment has the dimensions of (W × D × H ) 65.63 × 65.63 × 65.63 cm (corresponding to 25-times substrate length) and is shown in Fig. 4(b).This size was determined through a parameter sweep and found to be a good compromise between computational efficiency and the effect of the boundary condition (magnetic isolation) of the air domain on the magnetic field inside of it.An essential numerical limitation of this combined model is the large range of values for certain material parameters, such as the specific conductivity value.While this value can be in the range of 1-0.05 S/m for the tissues of the head as stated in Table II, it is close to 0 S/m for the surrounding air.In FEM modeling, and specifically 3-D models, these high contrasts of material parameters can lead to numerical stability problems, potentially resulting in random high-amplitude maxima in the magnetic field propagating through the air environment or nonconverging solutions (this was confirmed by COMSOL) [26].These instabilities also occurred in this work, and to combat these effects, we set the electrical conductivity of air to a much greater value, using 1 mS/m for the simulations in this work.We speak of this artificial medium as "highconductivity air."

A. Electrical Excitation and Propagation of Electromagnetic Fields
To test functionality and stability, we compute several studies on electric and magnetic field propagation through the entire modeling domain and repeatedly validate results against a homogeneous high-conductivity air environment.A fine tetrahedral mesh is chosen for the sensor and its thin, high aspect ratio layer structure.The head and the surrounding air volume with their comparatively large structures and lower aspect ratios utilize a successively coarser tetrahedral mesh to reduce computational load.The dimensions of the sensor and head are listed in Tables I and II, respectively.
An electric point dipole is placed inside the head geometry with an arbitrary dipole direction and magnitude.Based on the electromagnetic coupling in the model, this dipole introduces an electric field and a current density inside the head geometry, leading to the creation and propagation of a magnetic field.Fig. 5 shows a comparison of the electric current density norm, electric field density norm, and magnetic flux density norm based on electromagnetic field propagation from a dipole source inside the head.As a validation step, Fig. 5(a)-(c) demonstrates the results for a case in which all three layers of the head have their material parameters set to those of highconductivity air, while Fig. 5(d)-(f) shows the results with the material parameters for the head structure described in Table II.Fig. 5(a) and (d) illustrates the computed electric current density norm inside the head with an electric point dipole as the source.In this case, the electric point dipole is oriented in the y-direction and has a magnitude of 1 mA•m.For reference, a current amplitude of 1 mA at a standard tissue amplitude range is described as falling within the therapeutic range of DBS, and similar magnitudes within head models are described in [24] and [27].Furthermore, magnetic flux density measurements with fluxgate magnetometers obtained magnitudes in the nT and pT ranges for monopolar and bipolar DBS stimulations, respectively, validating a range of excitation amplitudes for our calculations [28].Due to the linear behavior, the simulated excitation source can be scaled in amplitude to investigate different excitation sources, such as various DBS electrode configurations or brain activity during epileptic seizures.The dipole is placed at distances of x = y = −4 cm and z = +2 cm from the center of the head geometry.To increase the numerical stability of the solution in setups with highly symmetrical source positions, an additional dipole source with the same direction and a magnitude of 1 µA•m was placed at coordinates x = z = 0 and y = 1 mm in the center of the white matter region.Furthermore, the logarithmically scaled arrows represent the current density norm vector inside different tissue types of the head.This illustrates how the current is guided by the shape of the head geometry and the tissue layers' conductivity and permittivity parameters.Because of the geometry-and material-based field propagation inside the head, the difference in propagation of the electric field density norm is apparent in Fig. 5(b) and (e), while the propagating magnetic flux density for the two cases can be compared in Fig. 5(c) and (f).The magnetic flux density is computed based on Ampere's law in the entire model domain and couples the electric field propagation with the material parameters and geometries to calculate the propagation of a magnetic field in the modeling space.The magnetic field that propagates through the high-conductivity air environment and ultimately reaches the sensor array dictates its behavior.Fig. 5 illustrate how the magnetic flux density is dependent on the head geometry through the electric material parameters, even though constant vacuum permeability for white matter, skull, and skin is assumed [6], [29].

B. Sensor Response and Array Behavior
After computing and analyzing the propagation of the magnetic field originating from the point dipole source through the head and high-conductivity air environment, we calculate the magnetic field arriving at the sensor array.The sensor's Metglas layer exhibits a deformation based on the MS effect, and via the mechanical coupling of the composite layers, a PE potential is generated between the AlN layer's top and bottom surfaces, which defines the sensor's electric output signal.The voltage between these surfaces is the transduced electric response for a given magnetic excitation field.The ME exhibits frequency-and angle-dependent responses regarding the incident external magnetic field.To visualize these dependencies within our combined sensor and head model, we performed a frequency sweep with varying excitation frequencies for the dipole inside the head and the full combined model.Fig. 6(a) shows the frequency characteristic for the sensor in the x-direction and excitation frequencies between 825 and 875 Hz.The blue line corresponds to the mean absolute potential generated over the PE layer of the sensor and shows a distinct peak at the first resonant bending mode at 848 Hz.The black line corresponds to the mean magnetic flux density norm inside the volume of the MS layer when this layer's material parameters are set to air.This configuration confirms that the magnetic flux density propagating through the air volume and into the sensor geometry is constant over the excitation frequencies and no frequency-dependent influence from the head's tissue layers is affecting the sensor's frequency characteristics.Important to note is that the sensor model in this work does not include damping effects due to sensor materials or the surrounding air volume.Previous investigations on damping effects reveal relevant contributors such as viscoelastic, molecular, and thermoelastic damping and compare sensor performance for air damping under atmospheric pressure and vacuum [30].It is also possible to implement experimentally obtained quality factors in postprocessing of the simulation results and to prescribe arbitrary damping coefficients based Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.on the Fano resonance profile of the sensor frequency characteristic.This enables dynamic investigations of different sensor damping configurations and calculating ME coefficients after running the FEM simulation [31], [32].This work omits damping effects and calculating the ME coefficient, as our focus is to highlight the relative behavior of the ME sensors in the orthogonal array for different dipole orientations and positions and investigate the position adaptability of the magnetic field detection in highly inhomogeneous fields.In Fig. 6(b), we show the angular dependence of the normalized PE voltage of a single sensor in the x-direction obtained in a homogeneous magnetic field with varying angle of incidence relative to the x-(azimuthal angle) and z-axes (polar angle).We observe a strong sinusoidal dependence for both azimuthal and polar angles of incidence and a maximum sensor response for incident angles parallel to the long axis of the sensor in the x-direction.Computing the total magnetic flux density in the MS layer and the potential across the PE layer geometry can be performed for arbitrary combinations of dipole positions and orientations, as well as ME sensor positions around the head.Fig. 7 illustrates the previously discussed dipole setup and the electric behavior of the three orthogonal ME sensors for two different positions of the array relative to the source inside the head.Fig. 7(a) shows the tangential component of the magnetic flux density in the xy plane at the location of the MS layer for the sensor oriented in the y-direction at approximately 10 cm distance to the surface of the head.Red arrows (logarithmically scaled) illustrate the vector field of the magnetic flux density.The ME sensor array is positioned at y = −4 cm (position marked as "A").Fig. 7(b) shows the same setup as in Fig. 7(a), but with the sensor array positioned at y = +5 cm (position marked as "B").
The plot of the magnetic flux density in the xy plane visualizes its inhomogeneity.Depending on their position along the y-axis, the sensors are exposed to different compositions of parallel and orthogonal magnetic flux vector components, relative to the long axis of the cantilever.As an example, we take a closer look at the sensor in y-direction.In cases of strongly orthogonal components, such as in position A, the sensor's MS layer exhibits low magnetic flux densities.At position B, the vector field has increasing parallel components leading to an enhanced magnetic flux density inside the MS layer, which ultimately leads to a stronger sensor response.This effect can be extended and analyzed for different sensor  locations and all sensors in the array, as shown in Fig. 8. Here, we demonstrate the strong location-dependent behavior of the three sensors, including the sensor in y-orientation (red line) with its minimum at y = −4 cm and maximum at y = +5 cm, corresponding to the respective locations "A" and "B" highlighted in Fig. 7.With increasing position in positive y-direction, this effect is counteracted by the growing distance to the dipole source.Note that the sensors in xand z-orientations show similarly constant behavior as opposed to distinct extrema, while the latter exhibits the strongest magnetic flux density, due to the dominant z-component of the magnetic field in this dipole configuration (the arrows in Fig. 7 only depict the tangential component of the magnetic flux density to exemplary highlight the behavior of the sensor in the y-direction).The increased magnetic flux density leads to a higher PE potential in the PE layer through the ME effect and, ultimately, a stronger sensor response.
To verify that this effect is explained by the interaction of the sensors' behavior and the head model, as opposed to solely the propagation of the magnetic field through the head and air environment, we perform a validation step.We replace the entire sensor array geometry with the same volume of air and compute the magnetic flux density norm in the air volumes previously corresponding to the three MS layers.The results for these air volumes are displayed as the dotted lines in Fig. 8, whereas the magnetic flux density norm inside of the actual MS layers corresponds to the solid lines.These results confirm that all three air volumes, without any sensor material present, have near identical magnetic flux density norms, with the difference being slight distance variations to the source due to the orthogonal sensor positioning.In contrast to that, the simulations containing the sensor model with MS layers show distinct behavior based on the orientation of each sensor.We conclude that each of the three sensors is exposed to a very similar external magnetic field, but exhibits vastly different results based on its previously described directionally dependent behavior [14], [15], [16].
In order to perform another validation step to confirm the directional dependence of the ME sensors, as well as the dependence of the propagated magnetic field on the electric parameters of the head, we execute the same simulation with a variable y-position sweep of the sensor array (as shown in Fig. 8) in combination with setting the material parameters of the head structure to the properties of high-conductivity air [analogous to the setup described in Fig. 5(a)-(c)].The resulting plots are shown in Fig. 9 where each of the three sensors is compared for both cases of the head material parameters set to tissue values according to Table II (solid lines) or high-conductivity air (dashed lines) for the setup with an electric point dipole source oriented in the y-direction and the integrated absolute potential over the PE layer is given for each array position.The behavior for the sensor in x-orientation is shown in Fig. 9(a) and is close in both cases, with a slightly higher potential generated for the case of the head geometry set to air material parameters.Fig. 9(b) shows the plot for the sensor oriented in y-direction.The influence of the material parameters of the tissue on the propagation of the magnetic field leads to the curves showing significant differences.In the case of the head geometry consisting of air, the sensor shows a steadily declining potential for the location sweep in a positive y-direction.In the case of the head geometry with assigned tissue parameters, the previously discussed behavior for the y-sensor becomes apparent again based on the effects of the electrical properties of the tissue on the magnetic field propagation.Due to the dominant z-component of the magnetic field, the integrated absolute potential over the PE layer is highest in Fig. 9(c) in both cases compared to the other sensors, with the dipole in the air environment generating a potential which is roughly one order of magnitude higher.
After the validation, a synoptic overview with all three sensors in the array combined with three different excitation dipole orientations is given in Fig. 10.Here, the effects on each orthogonal sensor for the dipole excitation with x-, y-, and z-orientations are given in Fig. 10(a)-(c), respectively.The different reactions of the sensors based on the dipole orientation are apparent, e.g., the minima for the sensors in x-, and y-orientations for different dipole moment directions, or the dominating potential of the sensor in z-direction in Fig. 10(a) and (b) compared to its near constant and minimal values in Fig. 10(c).We can pinpoint the differences in behavior seen in Fig. 10 on the strong anisotropy in the coupling parameters of the PE and MS layers.The materials' coupling is at its maximum when the incident magnetic field is parallel to the long axis of the cantilever, while it is minimal for angles of incidence that are perpendicular to it.Fig. 10 emphasizes the importance of analyzing the sources and sensors in a combined model, as slightly different dipole configurations can already lead to strong variation in the resulting sensor signals, based on the position and orientation of the sensor array relative to the source and the electric parameters of the head's tissues.Although only a limited number of examples can be presented in this work, theoretical investigations of different positions and orientations of the sensors, as well as a wide array of source configurations, can be modeled.

IV. CONCLUSION
In conclusion, we have demonstrated a combined multiscale 3-D FEM model consisting of an array of three orthogonal ME sensors and a simplified human head.Using electric point dipole sources located inside the head, we have simulated the excitation and propagation of electric and magnetic fields inside white matter, skull, and skin tissue domains, through the surrounding air environment, and into the ME sensor array.Inside the sensor layers, we have evaluated the influx of magnetic field, as well as the resulting electric potential in the PE layer, based on the ME effect, for different ME sensor and dipole source positions and orientations.This comprehensive investigation allows us to study how the head geometry with its three tissue types and corresponding material parameters influences the propagation of the electromagnetic Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
field through the tissues and outside into the air environment, where it is picked up by the ME sensor array.In addition to that, we can compare sensor behavior for different sensor positions and orientations relative to the source inside the head and observe different responses.Furthermore, validation studies with the electric parameters of the head or sensor geometries set to high-conductivity air demonstrate the importance of modeling the ME sensor array in conjunction with the appropriate head geometry and material properties to achieve distinct and realistic sensor responses, as opposed to using simplified homogeneous excitation fields or relying on simple air environments.The ability to analyze both sensor behavior and the influence of the source configuration in a single comprehensive model should be a key aspect of future theoretical studies.Comparing the responses of differently oriented sensors for arbitrary source configurations holds the potential to optimize localization applications with arrays of ME sensors, such as DBS electrode localization [33] or ultrasound head localization [34].Our setup also allows future investigations to take a closer look at the coupling between sensors in the array in order to optimize array geometries by minimizing crosstalk between sensors.A crucial question for upcoming investigations with the combined model is whether there is an interaction between ME sensors and the electromagnetic field inside the head.Combining sensor physics with the head model, follow-up investigations with a more advanced head geometry could, for example, answer the question whether the presence of the ME sensor introduces secondary effects, such as volumetric currents inside certain areas of the brain or modulation of the surrounding magnetic flux based on phenomena like sensor crosstalk or the magnetoelastic effect.Using different sensor geometries and the approach of frequency modulation, sensor systems that operate at relevant frequencies for their respective medical applications (e.g., MCG) could be modeled and further investigated with our approach.Finally, it is possible to utilize specialized head phantoms in measurement setups with corresponding FEM models to experimentally validate theoretical results [28], [35], which is a priority goal for further investigations.

Fig. 1 .
Fig. 1.(a) Schematic of a single ME sensor with given length, width, and thickness of the substrate layer (gray) and view of PE layer on the top surface of the substrate (red).(b) Sensor array with three orthogonal ME sensors oriented in x-, y -, and z-directions to form a vector sensor array.(c) Cross section of ME sensor with relative layer z-thicknesses to scale (but not the y -widths).Shown are the PE layer (red) on the top surface and the MS layer (blue) on the bottom surface of the substrate layer (gray).(d) Magnetic flux density norm and logarithmically scaled vector field for an ME sensor inside a homogeneous magnetic field parallel to the long axis of the sensor as an example for past investigations.(e) The aim of this work is to study the behavior of ME sensors inside inhomogeneous magnetic fields, such as those propagating outwards from a source inside of a human head, in a combined ME sensor and human head model.

Fig. 2 .
Fig.2.Cross section of the simplified human head based on the three-shell model.Three tissue types are represented with concentric spheres: skin, skull, and white matter (from outermost to innermost sphere).The radii are 88 mm for the white matter layer, 92 mm for the skull layer, and 100 mm for the skin layer.The relative permittivity and specific conductivity are assigned to the tissue types according to the literature.

Fig. 3 .
Fig. 3. (a) Overview of model geometry for the entire combined head and sensor model.A rectangular air surrounding (not depicted) contains the head model with its three tissue layers, as well as the sensor array with three orthogonally positioned sensors.(b) Cross section of the head model with adjacent ME sensor array.The concentric spheres represent human tissue layers of skin, skull, and white matter (innermost to outermost layer).Inside the head geometry, an electric dipole source can be placed at arbitrary positions and with variable excitation parameters, such as frequency and amplitude.

Fig. 4 .
Fig. 4. Model geometry with overlaid mesh elements: (a) zoomed-in view of head geometry with the sensor array positioned in front of it and a cross section of the skin and skull tissue layers; (b) full scale of the head and sensor with the surrounding rectangular air volume.This air volume has the dimensions of (W × D × H) 65.63 × 65.63 × 65.63 cm; and (c) magnified view of the base of the ME sensor in z-orientation.A fine tetrahedral mesh is chosen for the sensor and its thin, high aspect ratio layer structure.The head and the surrounding air volume with their comparatively large structures and lower aspect ratios utilize a successively coarser tetrahedral mesh to reduce computational load.The dimensions of the sensor and head are listed in Tables I and II, respectively.

Fig. 5 .
Fig. 5. Computed electric current density norm, electric field density norm, and magnetic flux density norm for an electric point dipole source at coordinates x = y = −4 cm and z = +2 cm in the white matter domain.The dipole moment is oriented in y -direction and has a magnitude of 1 mA•m.(a) and (d) Electric current density norm inside the head geometry for the head material parameters set to air with a conductivity of 1 mS/m and realistic tissue values, respectively.The red arrows represent logarithmically scaled current density vectors (tangential components) in the xy plane and illustrate that the current density is strongly dependent on head geometry and electric material properties.(b) and (e) Resulting electric field norm and (c) and (f) magnetic flux density norm inside and outside of the head.The figures highlight the dependence of the magnetic flux density on the head geometry through the electric material parameters.

Fig. 6 .
Fig. 6.(a) Frequency characteristics of the sensor oriented in x-direction.The blue line corresponds to the mean absolute potential over the PE layer generated at different frequencies with a maximum at the first resonant bending mode at 848 Hz.The black line illustrates the resulting mean magnetic flux density norm inside the volume of the MS layer of the sensor with material parameters set to air, demonstrating that no variation in magnetic flux occurs based on the magnetic field propagation through tissue and air at different frequencies.(b) Angular dependence of the normalized PE voltage for a sensor in the x-direction in a homogeneous magnetic excitation field of incident azimuthal angle θ and polar angle φ.A strong sinusoidal dependence is observed with maximum sensor response at incident field angles parallel to the sensor's long axis.

Fig. 7 .
Fig. 7. Magnetic flux density norm and logarithmically scaled vector field (tangential component) for the excitation with an electric point dipole in y -orientation.Shown is the xy plane with a cross section through the head with a dipole source and the MS layer of the sensor in y -orientation.(a) Sensor is located at y = −4 cm on the y -axis, annotated by the position marker "A."The inset in (a) shows the magnetic vector field propagating orthogonally to the long axis of the sensor, resulting in a minimal magnetic flux density and therefore electric response of the sensor.(b) Sensor is located at y = +5 cm on the y -axis, annotated by the position marker "B."The inset in (b) shows the magnetic vector field propagating with a smaller angle relative to the long axis and therefore a larger parallel component.This leads to a higher magnetic flux density and therefore electric response of the sensor, compared to the case in (a).

Fig. 8 .
Fig. 8. Simulation results for integrated magnetic flux density in the MS layers for the ME sensor array and a dipole in y -orientation.The dipole is located at coordinates x = y = −4 cm and z = +2 cm in the white matter domain and oriented in the y -direction.The solid lines show the integrated magnetic flux density norm inside each individual sensor for the variable position of the sensor array on the y -axis.The dashed lines show the integrated magnetic flux density norm inside the air volume previously occupied by the sensor geometry, without any sensor material present.The graphs represent the directional behavior of the three orthogonal sensors based on their position in the magnetic vector field surrounding the head in the interval between y = −5 and +5 cm and the anisotropy of the PE and MS couplings.

Fig. 9 .
Fig. 9. Simulation results for the integrated absolute value of electric potential in the PE layers of each individual sensor and different material parameters for the head.The dipole is located at coordinates x = y = −4 cm and z = +2 cm and is oriented in y -direction.Plots with solid lines represent results for the head geometry with material parameters given in Table II, while the dashed lines show results for the head geometry set to material parameters of air.The graphs represent the directional behavior of the sensors oriented in: (a) x-direction; (b) y -direction; and (c) z-direction based on their position relative to the source inside of the head.

Fig. 10 .
Fig. 10.Simulation results for the integrated absolute value of electric potential in the PE layers for the ME sensor array and three dipole configurations.The dipole is located at coordinates x = y = −4 cm and z = +2 cm in the white matter domain.The orientation for this dipole is: (a) in x-direction; (b) in y -direction; and (c) in z-direction.Distinct sensor responses for different cases are visible, including distinct minima at specific locations based on the angle of incidence between the magnetic vector field and the ME sensors.

TABLE III MESH
PARAMETERS USED FOR THE SENSOR, HEAD, AND SURROUNDING AIR ENVIRONMENT