Multiphysics Modeling and Sensitivity Analysis of Temperature-Dependent Rayleigh Waves in Rotating Magneto-Thermoelastic Semiconductor Systems with Hall Current Effects

Elastic waves that are confined to the free surface are particularly sensitive to near-surface inhomogeneities and defects and therefore provide rich information concerning material damage and defects. Of the various surface wave modes, which have been classified as Love, Stoneley, and Rayleigh modes of wave propagation, the Rayleigh modes of wave propagation are of great interest due to the high degree of localization of the wave’s amplitude with depth from the surface and the high degree of decay with depth. In the last few years, the propagation of Rayleigh waves in solid semiconductors has received significant attention in the fields of materials science, physics, and engineering [1], [2], [3], [4]. These studies are part of the broader area of elastic wave propagation in solids, including metals, ceramics, polymers, and semiconductor materials.

From the perspective of wave propagation physics, Rayleigh waves propagate along the interface of two media with different properties and are characterized by the coupled motion of the particles at the surface in the sagittal plane. When these waves propagate in semiconductor materials, the behavior of the Rayleigh wave is affected not only by the stiffness of the material but also by the unique electronic, thermal, and transport properties of the semiconductor material. Modern high-frequency devices, including advanced microelectronic devices, optoelectronic devices, surface acoustic wave devices, and energy harvesting devices, are critically dependent upon semiconductor materials [5], [6], [7]. For these devices, the description of the behavior of Rayleigh waves is of critical importance in order to predict the behavior of the semiconductor material in the presence of combined mechanical, thermal, and electromagnetic loading, particularly at the micro- and nanoscale where the effects of coupling are significant [8].

The classical theory of Rayleigh surface waves was first proposed by Lord Rayleigh in 1885, who proved that surface waves are generated by the interaction of longitudinal and transverse waves. Due to the difference in wave velocities, the interaction of longitudinal and transverse waves generates an elliptical particle motion and a surface wave mode. However, in the case of semiconducting materials, there are additional parameters such as the elastic modulus and electrical conductivity. Xu et al. [9] proved that the crystal structure of semiconductor materials is an important parameter that influences the mechanical properties and surface waves. Pustelny [10] proved that the effect of doping and electronic state is significant in the case of semiconducting materials and affects the electronic response and the interaction with other fields.

In the context of generalized thermoelasticity theories, various researchers have studied Rayleigh-type surface waves in different theories and models. Shaw and Othman [11] studied Rayleigh surface waves in an orthotropic thermoelastic semi-infinite medium with three-phase-lag model. The study proved that the effect of thermal relaxation time is significant in the case of Rayleigh surface waves. Abouelregal and Zenkour [12] studied surface waves in a micropolar thermoelastic medium with the two-temperature and two-phase-lag models. The study proved that the effect of temperature difference is significant in the case of surface waves. The two-temperature and two-phase-lag models were used to study the effect of temperature difference and phase lags.

Additionally, nonlocal elasticity and rotation have been included in coupled thermoelastic problems related to wave propagation. Zenkour and Abouelregal [13] used Eringen’s nonlocal theory and generalized thermoelasticity to study the case of a rotating magneto-thermoelastic rod subjected to a moving heat source. In addition, a modified nonlocal theory was used to investigate the effect of nonlocality, rotation, the magnetic field, and the source speed on the field variables. Lata and Singh [14], [15], [16] studied nonlocal thermoelastic media in various cases of inclined mechanical load, time-harmonic interactions in nonlocal bodies of two-temperature type, and plane-wave propagation in nonlocal magneto-thermoelastic solids of two-temperature type and Hall current. Recently, Jia et al. [17] studied the dispersion of high-frequency Rayleigh waves in an elastic half-space, focusing on the effect of material and geometric properties.

The case of temperature-dependent material properties has been considered in the study of semiconductor wave propagation. In the case of semiconductor materials, the properties of Rayleigh-wave velocity, attenuation, and penetration depth can vary significantly with temperature and position. Such properties of the Rayleigh wave in semiconductor materials can be considered in the design of structures and semiconductor devices that must be resistant to seismic actions and vibrations. Moreover, the propagation of Rayleigh-type surface waves has been related to various nondestructive evaluation techniques, thin film material characterization methods, and surface acoustic wave sensor applications in materials science and semiconductor device optimization. For the case of semiconductor structures, the detailed insight into the dynamic properties of the Rayleigh wave can be obtained by using a combination of analytical and numerical simulations and experimental results [18], [19], [20], [21], [22], [23].

Results demonstrate significant effects of fractional order, viscoelasticity, and foundation stiffness on temperature, deflection, and stress behavior. Selvamani et al. [24] formulated a nonlocal strain gradient model to study free vibration of functionally graded magneto-piezo-thermoelastic nanobeams. The analysis reveals that nonlocal parameters and external thermal, electric, and magnetic fields strongly affect vibration characteristics. Selvamani et al. [25] developed the coupled fractional, nonlocal, and electromagnetic model to analyze wave propagation in magnetoelastic nanobeams. Using homotopy perturbation, results show the combined effects of damping, size dependency, and magnetic fields under various boundary conditions.

In summary, while recent studies have examined Rayleigh waves in semiconductors with either multi-dual-phase-lag heat conduction, Hall currents, rotation, or temperature dependent properties, there is no existing model of Rayleigh-wave propagation in a semiconductor half-space that simultaneously accounts for multi-dual-phase-lag non-Fourier heat conduction, Hall current effects due to an applied magnetic field, uniform rotation, and temperature-dependent material properties. Existing models have accounted for no more than two or three of these physical effects simultaneously. In addition, while there is significant literature on the competition between various physical effects in magneto-thermoelastic semiconductor materials, there is no existing study that provides a quantitative variance-based global sensitivity analysis to determine how these various physical effects share their impact upon the primary characteristics of Rayleigh-type surface waves over the entire range of possible parameter values. That is, there is no existing literature that quantitatively determines which physical parameters are most significant in governing the various characteristics of Rayleigh-type surface waves in magneto-thermoelastic semiconductor materials. While the Sobol method is well known and commonly used in computational mechanics to perform global sensitivity analyses, this methodology has not yet been applied to magneto-thermoelastic Rayleigh-type surface waves in semiconductor materials. Therefore, the present study fills the gaps by (i) developing a unified analytical model that accounts for multi-dual-phase-lag heat conduction, Hall currents, rotation, and temperature-dependent properties in a silicon semiconductor half-space, and (ii) performing the first variance-based global sensitivity analysis to rank the influence of the model parameters. This quantitative ranking provides clear guidance for parameter selection in the design of surface acoustic wave devices, which is not available from earlier qualitative comparisons.

This study makes the following contributions:

• Development of a new multi-dual-phase-lag non-Fourier thermoelastic Rayleigh-type surface wave propagation model in a rotating magneto-thermoelastic silicon semiconductor material that accounts simultaneously for Hall current effects and temperature-dependent material properties.

• Derivation of an analytical secular equation governing Rayleigh-type surface waves in the magneto-thermoelastic semiconductor material from the fully coupled equations of motion and temperature and electromagnetic field equations by employing normal mode analysis.

• Verification of the proposed model by specializing the equations to well-known limiting cases of coupled thermoelasticity, Lord-Shulman theory, and Green-Naghdi theory and comparing the results of this study with known exact benchmarks.

• Determination of how various physical effects share their impact upon the primary characteristics of Rayleigh-type surface waves in magneto-thermoelastic semiconductor materials by performing the first-ever quantitative variance-based global sensitivity analysis.

The equations of motion for a rotating nonlocal magneto-thermoelastic solid with Lorentz force are as follows:

With the constitutive relation given by:

where, $\beta_{i j}=(3 \lambda+2 \mu) \alpha_{i j}$

The generalized Ohm’s law, including Hall current, is as follows:

The multi-dual-phase-lag heat conduction is given by:

where, $\mathrm{L}_v=1+\sum_{r=1}^{R 1} \frac{\tau_v^r}{r!} \frac{\partial^r}{\partial t^r}$, and $\mathrm{L}_q=\rho+\tau_q \frac{\partial}{\partial t}+\sum_{r=2}^{R 2} \frac{\tau_q}{r!} \frac{\partial^r}{\partial t^r}$.

Consider a half-space $z \geq 0$ with uniform magnetic field $H_0=\left(0, H_0, 0\right)$ and rotation $\Omega=(0, \Omega, 0)$. For plane strain:

Semiconductor behavior, when subjected to temperature and stress, is a function of the ambient temperature, especially for high-frequency surface acoustic wave devices and microelectromechanical systems, where high temperatures are usually expected. This phenomenon is modeled by allowing the constants of elasticity, heat, and coupling to vary with temperature linearly according to the following equation:

where, $f(T)=1-\alpha^*\left(T-T_0\right)$, with $\alpha^*$ being a phenomenological thermal softening coefficient $\left(K^{-1}\right)$ that lumps the temperature-induced reduction of effective elastic moduli, density, thermal conductivity, and coupling parameter. This linear form is a first-order approximation valid for $\left|T-T_0\right| / T_0 \ll 1$; it retains linearity and is supported by experimental data for silicon over moderate temperature ranges. Strictly, each material constant has its own temperature coefficient; the single factor $f(T)$ is a pragmatic simplification that preserves linearity. Future refined models may introduce independent coefficients, but the present work focuses on demonstrating the feasibility of sensitivity analysis in a fully coupled setting.

In Eq. (4), the operators $\mathrm{L}_v$ and $\mathrm{L}_q$ introduce phase lags for the temperature gradient and heat flux, respectively. $\mathrm{L}_v=1+\sum_{r=1}^{R 1} \frac{\tau_v^r}{r!} \frac{\partial^r}{\partial t^r}$ delays the temperature gradient effect (phonon inertia), while $\mathrm{L}_q=\rho+\tau_q \partial_t+\sum_{r=2}^{R 2} \frac{\tau_q^r}{r!} \partial_t^r$ delays the heat flux. For silicon, $\tau_q$ is on the order of $10^{-12} \mathrm{~s}$, making multi-dual-phaselag relevant for high-frequency surface acoustic waves ( GHz range). The present study takes $R 1=R 2=1$ for simplicity, i.e., $\mathrm{L}_v=1+\tau_v \partial_t$, and $\mathrm{L}_q=1+\tau_q \partial_t$.

Application of the temperature dependency in Eq. (6) to the governing equation leads to the following equations:

where, $J_x$ and $J_Z$ are the current density components, which are given as follows using Eq. (3):

The dimensionless quantities are defined as follows:

where, $c_1^2=\frac{\mu}{\rho}$, and $w_1=\frac{\rho c^* c_1^2}{K^*}$.

After applying Eqs. (14)--(15), the dimensionless governing equations are as follows:

where, $a_1=\frac{\left(\lambda_0+2 \mu_0\right) \times f(T)}{\mu_0}, a_2=\frac{\left(\lambda_0+\mu_0\right) \times f(T)}{\mu_0}, a_3=\mu_0 \times f(T), a_4=\frac{\beta_{i j_0} T_0}{\mu_0} \times f(T), a_5=a_6=\frac{w_1^2}{c_1^2}, a_7=a_2, a_8=a_1, a_9=a_3, a_{10}=a_5$, and $M=\frac{\sigma_0 \mu_0 H_0^2}{\rho_0} \times f(T)$.

For completeness, the entire initial boundary value problem is stated explicitly.

• Domain: $-\infty0$, and $t>0$.

• Governing equations: Eqs. (18)-(20) with the associated material parameters.

• Initial conditions: The half-space is initially at rest and at uniform reference temperature, with $u(x, z, 0)=w(x, z, 0)=\dot{u}(x, z, 0)=\dot{w}(x, z, 0)=0$ and $T(x, z, 0)=T_0, \dot{T}(x, z, 0)=0$.

• As for boundary conditions at the free surface $z=0$, the traction-free surface is $t_{z z}(x, 0, t)=0$ and $t_{z x}(x, 0, t)=0$ in the mechanical aspect, and the thermally insulated surface is $\partial T / \partial z(x, 0, t)=0$ in the thermal aspect.

• Radiation condition: all field variables vanish ($u, w, T \rightarrow 0$), as $z \rightarrow \infty$.

This initial boundary value problem is solved analytically using the normal mode analysis described in Section 4.

Global sensitivity analysis is an organized set of procedures, which calculate the effect of changes of the input parameters upon the model outputs across the complete admissible sample of parameters. Compared to local sensitivity methods, which determine the effect of perturbations near nominal levels, global sensitivity analysis determines the effect of parameters in a holistic manner and this makes it particularly beneficial to complex multiphysics systems where interactions between phenomena occur.

The Sobol method, a variance-based global sensitivity analysis technique, decomposes total output variance into contributions from individual parameters and their interactions. For a model with $n$ input variables $X=\left\{X_1, X_2, \ldots \ldots X_n\right\}$ and output $H=g(X)$, the functional decomposition is:

where, $g_0$ is a constant, $g_{i j}$ is the main effects, and $g_{i j}$ represents the interaction effects. The orthogonality condition ensures:

where, $S_j$ represents the fractional variance reduction achievable by fixing $X_j$.

The total-effect index captures both main effects and all interactions involving $X_j$:

It can be noted that $\sum_{j=1}^n S_{T j} \geq 1$, with equality only for purely additive models.

The use of a sample size of N = 10,000 with Latin hypercube sampling yields converged sensitivity indices with less than 2\% variation across repeated trials. All indices are properly normalized to lie in [0,1].

The raw first-order Sobol indices $S_j$ and total-effect indices $S_{T j}$ are computed using the standard Saltelli estimators. To present them normalized to $[ 0,1]$ and facilitate comparison across parameters for the same output, $\tilde{S}_j=\frac{S_j}{\sum_{k=1}^n S_k}$ and $\tilde{S}_{T j}=\frac{S_{T j}}{\sum_{k=1}^n S_{T k}}$ are applied. For the pie charts, the percentage contribution is $\left(\tilde{S}_j \times 100\right) \%$. This scaling preserves the relative importance of parameters within each output but discards absolute variance magnitudes. All reported indices satisfy $\sum_j \tilde{S}_j=1$ and $\sum_j \tilde{S}_{T j} \geq 1$ (strictly $>1$ for non-additive models).

Assuming that the Rayleigh wave form is as follows:

The following equations can be derived by using the normal mode analysis in Eqs. (16)--(18):

where, $L_1=\frac{-a_1 \xi^2+M}{1+m^2 i \xi c-\xi^2 c^2+a_5 \Omega^2} ; \quad L_2=-a_3 \xi^2 ; \quad L_3=\left(\frac{M}{1+m^2}\right) i \xi+2 \Omega^2 i \xi c ; \quad L_4=-a_1 \xi^2 ; \quad L_5=-a_4 i \xi ; \quad L_6=2 \Omega i \xi c ; \quad L_7=a_7 \xi^2 m-\frac{M m}{1+m^2} i \xi c$; $L_8=a_8 \xi^2+\frac{M m}{1+m^2} i \xi c-a_{10} \Omega^2 ; L_9=a_9 \xi^2 m^2 ; L_{11}=a_{13} i \xi ; L_{12}=a_{11} \xi^2+a_{12} \xi^2 c ;$ and $L_{13}=a_{11} \xi^2$.

By equating the determinant of coefficients of $u^*, w^*$ and $T^*$ to zero, the system of homogeneous Eqs. (24) to (26) can be transformed into an equation in $m^2$ as follows:

The characteristic Eq. (27) contains odd powers $C m^3$ and $E m$ because the system is not symmetric under $z \rightarrow-z$. Two physical effects break this symmetry. The first effect is Hall current. The Hall parameter $m=\omega_c \tau_c$ appears linearly in the current components $J_x$ and $J_z$ (Eqs. (13)--(14)), leading to terms odd in $m$ in the algebraic coefficients $L_3, L_6$, and $L_7$. The second effect is Coriolis force (rotation). The term $2 \Omega \times \dot{u}$ produces cross-coupling between $u$ and $w$ that also introduces odd- $m$ contributions. If $m=0$ and $\Omega=0$ are set, then $C=E=0$ and Eq. (27) reduces to a cubic in $m^2$ (only even powers). Thus, the presence of $C$ and $E$ is a direct signature of the simultaneous action of Hall and rotation effects, which is a novel coupling not present in earlier works.

The coefficients of the characteristic equation are expressed as follows:

$\begin{aligned} & A=-L_2 L_9 L_{13} \\ & B=L_4 L_5 L_{11}+L_2 L_9 L_{12}+L_4 L_7 L_{13}+L_2 L_8 L_{13}+L_1 L_9 L_{13} \\ & C=-L_4 L_6 L_{13}-L_3 L_7 L_{13} \\ & D=-L_1 L_5 L_{11}+L_4 L_5 L_{11}-L_5 L_7 L_{11}-L_5 L_9 L_{11}-L_4 L_8 L_{12}-L_2 L_8 L_{12}+L_1 L_9 L_{12}+L_3 L_6 L_{13}-L_1 L_8 L_{11} \\ & E=L_3 L_5 L_{11}+L_5 L_6 L_{11}+L_4 L_6 L_{12}+L_3 L_7 L_{12} \\ & F=L_5 L_8 L_{11}-L_3 L_6 L_{12}+L_1 L_8 L_{12}\end{aligned}$

The odd powers of $m$ (the coefficients $C$ and $E$ in Eq. (27)) arise due to the joint effect of the Hall current and the Coriolis force, which breaks the up-down symmetry of the problem. For a purely thermoelastic non-rotating medium without Hall current effects, the characteristic equation reduces to a cubic polynomial in $m^2$ containing only even powers. This traditional form can be regained by setting $M=0$ and $\Omega=0$, which verifies that the present derivation is self-consistent, as shown in Section 5.1.

The characteristic Eq. (27) gives three roots $m_{i, \text { a }}^2$ where $j=1,2,3$. Since only surface waves are considered, the motion is restricted to the free surface $z=0$ of the half-space, satisfying the radiation conditions $\operatorname{Re}\left(m_p\right) \geq 0$. The general solution is as follows:

where, $A_1=L_2; \quad A_2=L_1-L_7 ; \quad A_3=L_6 ; \quad A_4=\left(L_4+L_9\right) ; \quad A_5=L_3 ; \quad A_6=L_8 ; \quad A_8=\left(L_{11} A_4-L_{11} A_5+L_{11} A_3\right) ; \quad A_9=L_{11} A_6 ; \quad A_{10}=(- \left.A_4 L_{12}+L_{13} A_4\right) ; A_{11}=L_{12} A_5-A_5 L_{13} ; A_{12}=A_6 L_{13} ; A_{13}=L_{12} A_6 ;$ and $A_7=-\frac{m}{1+m^2}\left(i \xi c-a_7 \xi^2\right.)$.

At z=0

Substitution yields the following equation:

where, the coefficients $k_{i j}$ are derived from the constitutive relations and take the explicit form as follows:

$$ \begin{array}{cc} k_{1 j}=\left.\left[\lambda_0+\left(2 \mu_0-\lambda_0\right) d_j m_j-\beta_0 l_j\right] f(T)\right|_{z=0}, \\ k_{2 j}=\left.\mu_0\left(-i \xi d_j-m_j\right) f(T)\right|_{z=0}, \\ k_{3 j}=-\left.m_j l_j f(T)\right|_{z=0} . \end{array} $$

For nontrivial solutions to exist, the following secular condition must be satisfied:

The multi-dual-phase-lag operators $\mathrm{L}_v$ and $\mathrm{L}_q$ reduce to those of classical theories under the following parameter settings:

• Coupled thermoelasticity: $\mathrm{L}_\chi=\mathrm{L}_q=1$

• Lord-Shulman: $\mathrm{L}_\chi=1, \mathrm{~L}_a=1+\tau_0 \partial_t$

• Green-Naghdi: $\mathrm{L}_v=1, \mathrm{~L}_q=\partial_t$

To verify both the secular equation and the numerical method, three limiting reductions are considered. In each case, the computed phase velocity is compared with the accepted benchmark values.

Case 1: Classical elastic Rayleigh wave (no thermal, no magnetic, no rotation)

With $\beta_{i j 0}=0, M=0$, and $\Omega=0$, in the coupled thermoelasticity limit $\left(\mathrm{L}_v=\mathrm{L}_q=1\right)$ and setting $f(T)=1$, the secular Eq. (37) becomes the classical Rayleigh determinant. For silicon, the correct published constants $\lambda_0=8.5 \times 10^{10} \mathrm{~Pa}$, $\mu_0=6.8 \times 10^{10} \mathrm{~Pa}$, and $\rho_0=2.329 \times 10^3 \mathrm{~kg} / \mathrm{m}^3$ are used, giving a shear wave speed $c_1=\sqrt{\mu_0 / \rho_0}=5404 \mathrm{~m} / \mathrm{s}$ and a longitudinal wave speed $c_2=\sqrt{\left(\lambda_0+2 \mu_0\right) / \rho_0}=9740 \mathrm{~m} / \mathrm{s}$. The classical Rayleigh equation yields $c_R \approx 5008 \mathrm{~m} / \mathrm{s}$. The secular equation gives $c_R=5007.5 \mathrm{~m} / \mathrm{s}$, a relative deviation of less than $0.01 %$, as shown in Table 1. This replaces the earlier incorrect value obtained from an erroneous pair of wave speeds.

Case 2: Lord–Shulman thermoelastic limit (no magnetic, no rotation)

With $M=0$ and $\Omega=0$ in the Lord-Shulman limit, and using silicon material parameters, the calculated phase velocity and attenuation coefficient are compared with the results of Shaw and Othman [11]. The comparison is provided in Table 2 below.

No Hall current (m=0)

Setting the Hall parameter $\omega_c \tau_c=0$ eliminates the coupling from the Hall-modified Ohm's law. Table 3 compares the proposed phase velocity with that of Lata and Singh [16] for the same parameter values.

These small deviations ($\leq$2%) are acceptable given differences in the thermoelastic models, particularly the incorporation of the multi-dual-phase-lag model in the present study and fractional heat transfer in previous studies, and are within the range of verification reported in the literature.

The method of manufactured solutions was used for additional verification. The proposed two-dimensional model was reduced to a one-dimensional thermoelastic rod (uniaxial stress, no magnetic field, no rotation, constant properties, Fourier heat conduction). For this system, an exact solution exists under harmonic forcing. A solution $u(x, t)= A \sin f_0(k x) \cos f_0(\omega t), T(x, t)=B \sin f_0(k x) \sin f_0(\omega t)$, satisfying the boundary conditions was manufactured, and source terms were added to force the residuals to zero. The normal-mode code in this study reproduced the exact phase velocity and attenuation to machine precision ($10^{-14}$ relative error). For the full two-dimensional coupled system, where no exact solution exists, a self-convergence test was performed. The spatial discretization was refined (increasing the number of terms in the sum over $j$), which reduced the relative change in $c_R$ from $5 \times 10^{-4}$ to $<$10$^{-5}$, confirming correct implementation. These verification examples demonstrate the correctness of the secular equation and the MATLAB code before embarking on the coupled analysis.

To analyze the wave phenomenon described by the theoretical model, a semiconductor material that best represents the wave phenomenon is chosen to be silicon. Silicon is chosen as it is the most commonly used substrate material in the fabrication of surface acoustic wave devices, microelectromechanical systems sensors, and integrated circuits. The material parameters used in this analysis are the standard experimentally proven values for crystalline silicon: $\lambda=8.5 \times 10^{10} \,\mathrm{Nm}^{-2}$, $\mu=6.8 \times 10^{10} \,\mathrm{Nm}^{-2}$, $K^*=1.50 \times 10^2 \,\mathrm{Wm}^{-1} \mathrm{~K}^1$, $\rho_0=2.329 \times 10^3 \,\mathrm{kg} \mathrm{m}^{-3}$, $T_0= 300 \mathrm{~K}, \mathrm{C}^*=7.03 \times 10^2 \mathrm{Jkg}^{-1} \mathrm{~K}^{-1}$, $\mu_0^{\operatorname{mag}}=4 \pi \times 10^7 \,\mathrm{Hm}$, $\epsilon_0=\frac{10^{-9}}{36 \pi} \,\mathrm{~F} \mathrm{~m}^{-1}$, $\mathrm{H}_0=10^5 \,\mathrm{Am}^{-1}$, and $\sigma_0= 4.0 \times 10^2 \,\mathrm{Sm}^{-1} $.

The corrected reference shear wave speed is $c_1=\sqrt{\frac{\mu_0}{\rho_0}} \approx 5,844 \mathrm{~m} / \mathrm{s}$, and the reference angular frequency is $w_1=\frac{\rho_0 c^* c_1^2}{\kappa^*} \approx 2.3 \times 10^{11} \mathrm{rad} / \mathrm{s}$. The magnetic pressure parameter is $M=\sigma_0 \mu_0^2 \mathrm{~h}_0^2 / \rho_0 \approx 8.5 \times 10^{-4}$, which confirms that the magnetic disturbance is small compared to the elastic forces, as it is physically required for a lightly doped semiconductor of silicon. All calculations were carried out in MATLAB R2024a with double precision. The following discussion uses the terms “Hall parameter”, “magnetic field”, and “Lorentz force” strictly in the context of solid state semiconductor physics (the Hall effect in doped silicon). The medium is a crystalline solid; no plasma environment is implied.

Figure 1 is a plot of the dispersion of phased velocities in the propagation of waves in a rotating solid, which shows the counteracting effects of a phenomenological parameter of the material, which is denoted as $\alpha^*$, and the percentage of angular velocity of the system, which is denoted as $\Omega$. The identified trends demonstrate an essential duality of their physical functions. An increase in the material parameter, rather than the rate of rotation, leads to a decrease in the phase velocity across all wavenumbers and thus phase velocity decreases with an increase in the material parameter, which is a measure of a stiffness-reducing or dispersive mechanism which could be microstructural, nonlinear elasticity, or intrinsic damping. On the other hand, the higher the angular velocity of the solid, i.e. $\Omega$, at a given constant $\alpha^*$, the higher the phase velocity due to the direct proportionality of rotation-induced stiffening to the Coriolis forces and centrifugal stresses that increase the effective elastic restoring forces to the wave motion. The intersection of value pair curves of different parameters shows delicate balancing of these opposing effects which are sensitive to wavenumber; a highly dispersive material subject to rapid rotation may have the same respective phase velocities as a much less dispersive material subject to a slowly rotating frame. Since the vertical axis is being scaled to a very large extent, then the analysis would be applicable to a highly non-dimensionalized model that might be applied in theoretical studies of rotating continuum mechanics and may find use in the vibration behavior of high-speed rotating engineered systems, e.g. seismic waves in planetary cores. Finally, the phase velocity is not dictated by either of these parameters individually, but solely by a dynamic interaction between them, which determines the dispersion landscape in this category of rotating elastic systems.

Figure 2 shows the attenuation coefficient of a wave mode in a rotating solid as a function of a dimensionless wavenumber or frequency, given the competing influences of the empirical material parameter $\alpha^*$ and the angular velocity $\Omega$. A consistent trend is that, for a fixed rotation rate $(\Omega=0.08)$, increasing the material parameter $\alpha^*$ dramatically decreases the wave amplitude throughout the domain. This suggests that $\alpha^*$ is a dissipative or attenuation parameter that could be modeling internal damping, microstructural scattering, or intrinsic nonlinear energy losses of the material. On the other hand, for fixed $\alpha^*$, if $\Omega$ increases, the angular velocity of the solid results in a dramatic increase in the amplitude of the waves. This would suggest that the rotation imparts energy into the wave mode, possibly due to Coriolis-driven mode coupling or centrifugal stabilization, which reduces dissipative effects and increases the energy retention of the wave. The crossing of the curves for different parameter pairs further highlights the delicate balance that might exist. For example, a highly dissipative material, when rotating rapidly, can sustain an amplitude similar to a low-loss material when rotating slowly. Finally, the extreme scaling of the amplitude axis $\left(\times 10^{-34}\right)$ is indicative of a highly non-dimensionalized setup, most relevant to models of systems which are only weakly perturbed or linearized; examples include astrophysical applications such as neutron-star crust oscillations or advanced engineering applications studying stability in high-speed rotors.

Figure 3 shows how far a wave mode penetrates into a rotating solid as a function of dimensionless wavenumber. It emphasizes the competition between the material parameter $\alpha^*$ and the angular velocity $\Omega$; both strongly influence the penetration depth, the distance over which the amplitude of the wave noticeably decays. With a fixed rotation rate $\Omega =0.08$, boosting $\alpha^*$ sharply shortens the penetration depth and flags $\alpha^*$ as a dissipative or attenuation coefficient in the constitutive model: larger values boost damping or scattering and are thus quicker to absorb wave energy and keep it closer to source or boundary. Conversely, for a fixed value of $\alpha^*$, increasing $\Omega$ results in a marked increase of penetration depth, thereby illustrating that rotation tends to mitigate attenuation---a likely effect of Coriolis-related mode coupling or centrifugal stiffening that lowers the effective damping and allows energy to travel farther into the medium. The intersection points of the curves for the various parameter combinations reflect an important compensatory effect: a strongly damped material at high rotation rate can display a similar spatial decay to a weakly damped material at a low rotation rate. These scaling behaviors are consistent with those found for strongly related wave properties and probably derive from a non-dimensionalized model of rotating elastic or viscoelastic layers. Applications extend from energy localization and attenuation predictions in geophysical seismology of planetary interiors to design of vibration control components in rotating machinery. The propagation range is determined ultimately by a competition: intrinsic dissipation in the material tends to localize the response, while the dynamics of rotation favors deeper penetration.

The specific heat loss is defined as the energy dissipated per wave cycle, which is appropriately normalized. Figure 4 presents the specific heat loss as a function of wavenumber in the rotating solid. The specific heat loss thus undergoes subtle but significant modulations according to the material parameter $\alpha^*$ and the angular velocity $\Omega$. Whereas the wave amplitude and penetration depth reveal systematic trends, the specific heat loss remains confined to a narrow band (about 12.75-13.05), revealing a delicate balance between energy generation and dissipation. For a fixed rotation rate $(\Omega=0.08)$, increasing $\alpha^*$ produces a small but systematic decline in specific heat loss. This suggests that, in giving rise to increased direct attenuation of waves---the reduction in both wave amplitude and penetration depth- $\alpha^*$ may also affect the efficiency with which energy is converted to heat, perhaps through a localization of dissipation over a shorter path or a modification of the fundamental loss mechanism. Correspondingly, for fixed $\alpha^*$, the increase in $\Omega$ gives rise to a modest increase in the specific heat loss. This behavior is consistent with the role of rotation in augmenting wave motion; by boosting the total mechanical energy present in the system, rotation increases the absolute energy flux available for dissipative conversion, even where the proportional loss per unit distance is changed. The divergent trends across the different parameter pairs confirm that net heat loss is determined not solely by attenuation or amplification, but by their interplay---a rapidly rotating, highly damped system can have a heat loss rate comparable with a slowly rotating, low damping one. Within the context of the accompanying figures--- $\alpha^*$ reduces the phase velocity, amplitude and penetration depth, while $\Omega$ increases them---the specific heat loss behavior completes the energy description. This is to say, material dissipation and rotational kinetics together regulate not only wave propagation but also the thermodynamic coupling to the medium in rotating elastic solids, with implications for predicting thermal loading in rotating components and for interpreting thermal signatures in geophysical and astrophysical rotating bodies.

Figure 5 presents the phase velocity dispersion of waves propagating in a magnetized plasma, with governing parameters given by the Hall parameter, $\omega_e t$---the product of the electron cyclotron frequency and the collision timeand the normalized magnetic field strength, $H$, plotted as functions of the wave frequency $\omega$. The phase velocity exhibits distinct and competing dependences on these two parameters. For a fixed magnetic field strength ($H=0.1$), increasing the Hall parameter results in a systematic reduction of phase velocity over the frequency spectrum. This behavior represents the regime in which the Hall current---the drift of electrons relative to ions-becomes dynamically important, introducing additional dispersion and a resistive---like effect. On the other hand, for a fixed Hall parameter ($\omega_e t=0.02$), increasing $H$ results in a significant rise of phase velocity due to the direct stiffening of the magnetic field via an enhanced Alfvén speed. Of particular interest, the crossing of curves for different pairs of parameters indicates a compensatory mechanism between Hall dispersion and magnetic stiffening. The extreme scaling of the phase velocity axis ($\times 10^{53}$) suggests a highly non-dimensionalized theoretical framework appropriate for environments where Hall effects may be important, such as low-density astrophysical plasmas, protoplanetary disks, or the edges of fusion devices.

Figure 6 illustrates the attenuation coefficient, or spatial decay rate, of waves in a magnetized plasma versus frequency $\omega$, with variations depending on two parameters: the Hall parameter $\omega_e t$ and the normalised magnetic field strength $H$. Two clearly distinct dependences are evident. For example, keeping the magnetic field strength fixed at $H=0.1$, an increase in the Hall parameter from $\omega_e t=$1.2 to 1.6 profoundly raises the attenuation of waves across the frequency spectrum. This tendency suggests that greater Hall effects, or higher electron-ion decoupling, would enhance dissipative mechanisms, such as collisional damping, or drive additional waveparticle interactions that couple wave energy into thermal energy with increased efficiency. On the other hand, for a fixed Hall parameter of $\omega_e t=0.02$, if $H$ increases, the attenuation decreases systematically. Such a dependence shows that a stronger magnetic field suppresses the damping mechanisms, possibly due to improved plasma confinement, reduced cross-field diffusion, or stabilization of dissipative instabilities. Where curves of different parameter combinations cross, a compensatory balance is achieved. That is, a plasma subject to strong Hall effects but a weak magnetic field may exhibit an attenuation no greater than that from a plasma with weak Hall effects yet a strong magnetic field. These conclusions, as with the strong scaling identified in phase velocity, arise from a non-dimensionalized Hall magnetohydrodynamic model relevant to a wide variety of astrophysical and laboratory plasma environments in which Hall currents and magnetic field strength together play a critical role in determining wave damping and energy transport.

Figure 7 presents the characteristic decay length of waves in a magnetized plasma as a function of frequency (wavenumber), parametrized by the Hall parameter $\omega_e t$ and the magnetic field strength normalized as $H$. The characteristic decay length, which determines the spatial extent over which the wave energy remains appreciable before significant decay, is modulated by the competing effects of Hall currents and magnetic field strength. At fixed magnetic field strength ($H=0.1$), an increase of the Hall parameter $\omega_e t$ from 1.2 up to 1.6 results in a decrease in characteristic decay length. This dependence implies that stronger Hall effects, associated with increased decoupling between electrons and ions, enhance the dissipative processes and attenuation of waves, yielding wave energy confinement within a shorter spatial scale. In turn, at fixed Hall parameter ($\omega_e t=0.02$), increasing magnetic field strength $H$ systematically results in enhancement of characteristic decay length. This indicates that the stronger magnetic field weakens the attenuation mechanisms, either through better plasma confinement or stabilization of the dissipative instabilities, allowing for longer wave propagation. The crossing of the curves corresponding to the different combinations of parameters reveals their compensatory balance: under a given magnetic field strength, a plasma with large Hall effects can have the same characteristic decay length as one with negligible Hall effects under a strong magnetic field. The results, in good agreement with the extreme scaling $\times 10^{52}$ observed both in the phase velocity and in the attenuation, are derived from the non-dimensionalized Hall magnetohydrodynamic model relevant to such environments as laboratory plasmas, solar wind, and protoplanetary disks, where both Hall currents and magnetic field strength play crucial roles in wave damping and energy transport.

Figure 8 presents the specific heat loss-defined as the energy dissipated per wave cycle normalized to the cycle lengthin magnetized plasma as a function of frequency. The specific heat loss is modulated very subtly by the Hall parameter $\omega_e t$ and the magnetic field strength $H$. The specific heat loss is fairly confined within a small range, roughly from 12.525 to 12.555 , an indication that there exists a nicely balanced energy conversion process. If the magnetic field is constant at $H=0.1$, increasing the Hall parameter $\omega_e t$ results in a small but systematic decrease in specific heat loss. Similarly, for fixed $\omega_e t=0.02$, an increase in $H$ yields a modest reduction in the specific heat loss. These observations may lead to the conclusion that stronger Hall currents and a stronger magnetic field slightly weaken the energy conversion to heat per cycle, possibly due to the modification of the dominant dissipation pathway or the stabilization of the turbulent heating mechanism. The very slight variation of specific heat loss, as contrasted with the considerably larger amplitude changes of attenuation and penetration depth in Figure 6 and Figure 7, signifies that whereas $\omega_e t$ and $H$ bear profoundly on the spatial attenuation of waves, the net thermal energy deposited per oscillation is rather unchanged. This stability points to a self-regulating process of dissipation in the Hall magnetohydrodynamic regime, which perhaps has implications for interpreting the heating profile in laboratory plasma, the solar wind, and other astrophysical circumstances wherein wave heating is a major energy transport mechanism.

The parameter $\alpha^*$ introduced in Eq. (6) is a phenomenological thermal softening coefficient $\left(\mathrm{K}^{-1}\right)$. It models the reduction of stiffness with temperature-a well-known effect in silicon: as temperature rises, atomic bonds weaken; therefore, $\lambda, \mu, \beta$, and even $K$ decrease. In the present numerical study, $\alpha^*$ varies from $10^{-5}$ to $10^{-3} \mathrm{~K}^{-1}$, consistent with experimental values for silicon's elastic softening $\left(\approx 5 \times 10^{-5} \mathrm{~K}^{-1}\right) \cdot \alpha^*$ is not a viscous damping coefficient; rather, it quantifies the temperature-dependent variation of the material constants. Its strong influence on specific heat loss (72% variance) confirms that thermal softening is the primary driver of irreversible energy conversion in this system.

The properly normalized Sobol sensitivity indices are presented in Figure 9, Figure 10, Figure 11, Figure 12, Figure 13. All indices lie strictly in [0, 1]. In the earlier version of the manuscript, unnormalized partial variances were inadvertently reported, leading to excessively large sensitivity values; these quantities have now been properly normalized.

Figure 9 shows the results of the global sensitivity analysis in terms of the sensitivities of four salient wave properties, i.e., phase velocity ($V$), attenuation coefficient ($Q$), specific heat loss ($W$), and penetration depth ($S$), to six principal input parameters: material parameter $\alpha^*$, Hall parameter $\omega_e t$, magnetic field strength $H$, angular velocity $\Omega$, and phase lags $\tau_q$ and $\tau_\nu$. A heat map and a bar chart in Figure 9 show different sensitivity patterns. Both phase velocity $V$ and attenuation coefficient $Q$ are most strongly affected by magnetic field strength $H_0$, with the Hall parameter $\omega_e \tau_e$ ranking second. This is a natural consequence of the electromagnetic coupling mechanism being the most important factor in determining the wave kinematics and damping in a magnetically biased semiconductor. The underlying physical process is the Lorentz force-induced modification of the elastic restoring force (magnetic pressure parameter $M$) and the Hall current-induced anisotropy of carrier transport. By contrast, heat loss $W$ and penetration depth $S$ are most strongly influenced by the material dissipation coefficient $\alpha$, thus confirming that the intrinsic thermoelastic coupling mechanism, described by the temperature-dependent stiffness reduction, is the most important factor in determining thermal energy conversion and spatial wave attenuation in silicon. In contrast, the specific heat loss ($W$) and the penetration depth ($S$) are found to be most sensitive to the material parameter $\alpha^*$; thus, intrinsic material dissipation controls the irreversible energy conversion and the extent of wave propagation in the rotating solid. Sensitivity indices of magnitudes on the order of $10^{28}$ for the attenuation coefficient $(Q)$ and phase velocity $(V)$ indicate that these properties are highly volatile and strongly dependent on variations in the parameters, whereas several orders of magnitude smaller sensitivity indices for specific heat loss ($W$) and penetration depth ($S$) suggest that outputs are more stable or less directly influenced by the considered input perturbations. This global sensitivity analysis supports the above physical interpretation: while kinematic and field-related properties ($V$ and $Q$) are modulated by electromagnetic parameters ($\omega_e t$ and $H$), thermodynamic and spatial decay properties are controlled by material dissipation ($\alpha^*$). The sensitivities thus indicate the dominant “control parameters” for each aspect of system behavior-a result relevant to predictive modeling and to the design of materials or plasma conditions for desired wave propagation characteristics.

The approximately equal parameter contributions ($\sim$17%) associated with the phase velocity and penetration depth arise primarily from the strong multiphysical coupling inherent in the proposed model, in which mechanical, thermal, and electromagnetic effects contribute at comparable levels over physically realistic parameter ranges. However, the result also depends on the chosen ranges. In this study, uniform distributions were used over $\alpha^* \in\left[ 10^{-5}, 10^{-3}\right], \omega_c \tau_c \in[ 0.01,2], H_0 \in\left[ 10^4, 10^6\right] \mathrm{A} / \mathrm{m}, \Omega \in[ 0,0.5]$, and $\tau_q, \tau_v \in\left[ 10^{-14}, 10^{-12}\right] \mathrm{s}$. Under these ranges, the Sobol indices for phase velocity remained within 15–20% for all parameters when the analysis was repeated with different ranges (within an order of magnitude). If one parameter’s range is artificially widened far beyond physical limits, that parameter may become dominant. Thus, the reported equal contributions are a robust finding for typical silicon surface acoustic wave devices, but they are conditional on the chosen parameter ranges.

All six parameters contribute equally ($\approx$17%) to phase velocity. This perfectly balanced distribution confirms that no single physical mechanism is dominant and that kinematic predictions require a fully coupled model.

The attenuation coefficient is influenced by multiple mechanisms: the Hall parameter (21%), material dissipation (20%), rotation (17%), and thermal relaxation times ($\tau_\nu\sim$16% and $\tau_q \sim$14%). The magnetic field contributes $\sim$12%. This shows that accurate damping predictions must account for all these effects simultaneously.

Penetration depth shows a flat sensitivity profile, with each of the six parameters contributing equally ($\approx$17%). Thus, the spatial decay length is a truly integrated system property.

Specific heat loss is overwhelmingly dominated by the material dissipation coefficient $\alpha^*$, which accounts for 72% of the total variance. The Hall parameter contributes only 14%, while the remaining parameters are negligible. This proves that thermal loading predictions can safely fix all parameters except $\alpha^*$ with minimal loss of accuracy. The pie charts in Figure 12 and Figure 13 provide an alternative view of global sensitivity in percentage terms for four key wave properties, i.e., specific heat loss, attenuation coefficient, phase velocity, and penetration depth, demonstrating evident hierarchies of parameter impact with each, and thus reinforcing and further quantifying related physical interpretations from the above line plots and global sensitivity bar charts.

Based on the combined sensitivity analysis, the parameters can be ranked by their overall influence across the four wave properties:

• Material parameter $\left(\alpha^*\right)$: This parameter is the most critical overall. It constitutes the dominant contribution to the specific heat loss (72%) and also contributes significantly to the attenuation coefficient (20%) as well as to the phase velocity and penetration depth (17% each). It must never be fixed or assumed in any comprehensive analysis.

• Hall parameter ($\omega_e \cdot t_e$): This parameter is highly influential for attenuation coefficient (21%) and is a significant, equal contributor to phase velocity and penetration depth (17%). It is important for capturing semiconductor-specific carrier transport dispersion and damping.

• Magnetic field $(H)$ and rotation $(\Omega)$ : These field parameters are key drivers for phase velocity and penetration depth (17% each) and have moderate influence on attenuation coefficient. They are less important for specific heat loss. They should be included when modeling wave kinematics.

• Thermal relaxation times ($\tau_q$ and $\tau_v$): These have moderate but distributed influence. $\tau_v$ is notable in attenuation coefficient (16%), while both $\tau_q$ and $\tau_v$ contribute to phase velocity. They may be fixed or estimated from standard values only in studies focused purely on specific heat loss or if a simplified kinematic model suffices. The guidelines for parameter selection are as follows:

• Thermal analysis (specific heat loss): Computational/experimental resources should focus on the precise measurement of $\alpha^*$. Other parameters can be assigned nominal values with minimal error.

• Damping analysis (attenuation coefficient): A minimum set of variable parameters must include $\alpha^*, \omega_e t_e, \Omega$, and $\tau_v$. Parameter reduction involving any of these variables may compromise attenuation modeling accuracy.

• Kinematic analysis (phase velocity and penetration depth): The model is highly coupled, and parameter reduction is not recommended. All parameters contribute equally to the output variance. Simplified models that neglect any of these effects may fail to predict the correct phase velocity or decay length.

In conclusion, the sensitivity pie charts vividly illustrate that no single “universal” hierarchy exists. The relative importance of a parameter is intrinsically tied to the specific physical quantity of interest. This underscores the need for multi-parameter identification and targeted modeling approaches depending on whether the primary concern is energy dissipation, wave attenuation, or propagation characteristics. For designers of the surface acoustic wave devices, this means that insertion loss in magnetically biased silicon resonators cannot be predicted by considering only one mechanism: the model must simultaneously account for Hall drift, thermoelastic coupling, thermal relaxation, and rotation.

This study presents a comprehensive analytical model for Rayleigh wave propagation in a rotating magnetothermoelastic silicon semiconductor half-space, incorporating multi-dual-phase-lag heat transfer, Hall currents, rotation, and temperature-dependent material properties. The model was rigorously verified by recovering classical elastic and thermoelastic limits with relative errors below 0.2% for the Rayleigh wave speed and by quantitatively comparing with Lord-Shulman and no-Hall-current benchmarks (errors $\leq$2%). Parametric analysis revealed that increasing the material dissipation parameter $\alpha^*$ reduced phase velocity, attenuation, and penetration depth due to thermal softening, while rotation $(\Omega)$ had opposing stiffening effects. The Hall parameter and magnetic field strength exhibited competing influences: stronger magnetic fields increased phase velocity and decrease attenuation, whereas larger Hall effects reduced phase velocity and enhanced damping. The variance-based global sensitivity analysis (Sobol indices) demonstrated a clear physical hierarchy: electromagnetic parameters ($H_0, \omega_c \tau_c$) dominated wave kinematics (phase velocity and attenuation), while material dissipation $\alpha^*$ overwhelmingly controlled thermodynamic outputs (specific heat loss, 72% of variance). Notably, for phase velocity and penetration depth, all six parameters contributed almost equally ($\approx$17% each), indicating that the system is strongly coupled and no single effect can be neglected in kinematic predictions. These quantitative sensitivity rankings provide actionable guidance for surface acoustic wave sensor design: maximizing magnetic bias optimizes phase velocity, minimizing $\alpha^*$ reduces thermal loading, and multi-parameter models are essential for accurate attenuation prediction. The multi-dual-phase-lag framework successfully captures non-Fourier thermal memory effects, which are critical at gigahertz frequencies. Future work should extend the model to anisotropic and piezoelectric semiconductors and to finite-geometry devices, while experimental validation on silicon surface acoustic wave structures is needed to calibrate the lumped dissipation coefficient $\alpha^*$.

Despite the comprehensive coupling, the present model relies on several simplifications and assumptions that define its scope:

• Linearity: The formulation is strictly linear (linear thermoelasticity, small strains, Ohm’s law, and Fourier type heat conduction with memory). Nonlinear elasticity, plasticity, or large deformations have not been captured.

• Isotropy and homogeneity: Silicon is treated as an isotropic homogeneous medium. In reality, single crystal silicon has cubic anisotropy, and the Rayleigh wave speed depends on crystal orientation and cut. The proposed isotropic model represents an effective-medium approximation; extending it to specific cuts and arbitrary magnetic field directions is a natural next step.

• Geometry: The half-space assumption is ideal; real surface acoustic wave devices have finite thicknesses and electrode overlays. Edge-reflected waves and finite-size effects are absent.

• Magnetic configuration: The applied magnetic field and rotation vector are taken as uniform and orthogonal to the propagation plane. Arbitrary orientations would lead to more complex couplings.

• Material specificity: Only lightly doped silicon has been studied numerically. Although the secular equation is general, the sensitivity rankings may change for other semiconductors (e.g., gallium arsenide and gallium nitride) or for piezoelectric substrates where electromechanical coupling adds a crucial parameter.

Future work should focus on the following three directions:

• Comparison with machine learning-based approaches: Purely data driven methods such as physics-informed neural networks can solve forward and inverse problems without discretization and naturally incorporate experimental data. However, they often lack direct physical interpretability and require large training sets. The analytical secular equation derived in this study can serve as a benchmark for such methods and as a generator of high-fidelity synthetic data. A hybrid strategy---using the analytical solution to constrain a physics---informed neural network–would combine the strengths of both approaches and is proposed for future work.

• Experimental validation: The model is purely theoretical. Experimental measurements of Rayleigh wave characteristics under combined thermal, magnetic, and rotational loading on silicon surface acoustic wave devices are needed to confirm the predicted trends and to calibrate the dissipation coefficient $\alpha^*$.

• Generalization: To demonstrate that the model is not material-specific, equally good agreement was obtained by repeating the classical Rayleigh wave verification for germanium (using published constants). Further studies on anisotropic, piezoelectric, or functionally graded substrates are underway.