Supported data and manuscript "Understanding change in the sound wave frequency in a ferromagnet under magnetic field influence (Simon effect) in the low-field regime"

Contents lists available at ScienceDirect Results in Physics journal homepage: www.elsevier.com/locate/rinp Understanding change in the sound wave frequency in a ferromagnet under magnetic field influence (Simon effect) in the low-field regime I. Korniienkoa, P. Nievesb, D. Legutc,a,∗ aIT4Innovations, VŠB - Technical University of Ostrava, 17. listopadu 2172/15, Ostrava, 70800, Czech Republic bDepartamento de Física, Universidad de Oviedo, C. Leopoldo Calvo Sotelo, 18, Oviedo, 33007, Spain cDepartment of Condensed Matter Physics, Faculty of Mathematics and Physics, Charles University, Ke Karlovu 3, Prague 2, 121 16, Czech Republic A R T I C L E I N F O Keywords: Sound wave velocity Simon effect Magnetoelastic effect Ferromagnet Dispersion A B S T R A C T Studies of coupled magnetic and elastic effects have a long history, however, these effects still hold the scientific interest of researchers. In particular, in recent years, there has been an increased interest in research on the interaction of surface acoustic waves with spin waves in ferromagnetic films. In turn, this calls for better understanding of low field regime of Simon effect (the effect is related to the influence of the applied magnetic field on the sound wave velocity in a ferromagnetic crystal). In our work, based on the example of bcc Fe, we propose a refined formula to describe the Simon effect, which contains terms related to dispersion effects associated with exchange stiffness. We compare our analytical solutions with other alternative computational approaches and show that dispersion effects can be significant for the Simon effect and cannot be neglected in the low field regime. As a result, we propose a more accurate analytical formula, which, due to its relative simplicity, can become a convenient tool to estimate the magnitude of the expected magnetic field effect on the sound wave speed propagation in a cubic ferromagnetic crystal, as well as it explains observed deviations from analytically expected results in Simon effect at low magnetic field. Introduction Magnetoelastic coupling is a fundamental feature of magnetic materials. It is responsible for numerous interesting magnetostrictive and stress effects, among which there are magnetovolume effect, Joule magnetostriction, Wiedemann effect, Villari effect, Matteucci effect, anomalous thermal expansion and many others [1–8]. Due to magnetoelastic coupling, the elastic properties of magnetically ordered materials depend on the magnetic ordering of the sample and the magnitude and direction of the external magnetic field. In particular, this leads to the dependence of the ultrasonic wave velocity on the applied magnetic field beyond technical saturation as 1∕𝐻, which was explained by Simon in 1958 using linear theory of magnetoelasticity (Simon effect) [9]. The low magnetic field regime, below saturation, is difficult to understand because the ultrasonic wave velocity is influenced by applied magnetic field through the induced domain wall displacement (𝛥s-effect) [10]. In this regime, interesting anomalies on the fractional change of ultrasound frequency has been experimentally reported by Rouchy et al. in cubic Co-Pt alloy [10]. Recently, similar ∗Corresponding author at: Department of Condensed Matter Physics, Faculty of Mathematics and Physics, Charles University, Ke Karlovu 3, Prague 2, 121 16, Czech Republic. E-mail address: [email protected] (D. Legut). patterns to these anomalies have been obtained by means of atomistic spin-lattice simulations [11] for bcc Fe assuming a saturated magnetic state in the low magnetic field regime, suggesting that such anomalies could not arise only from domain wall propagation. This situation leads to an interesting question about the possible origin of these discrepancies between analytical and experimental results. In addition to Simon effect (linear theory of magnetoelasticity), there are higher order effects that one might consider, like morphic effects arising from second order magnetoelastic coupling combined with anharmonic third order elastic constants, or rotational-magnetostrictive effects due to finite strain theory and magnetoelastic free energy rotational invariance (not predicted by infinitesimal strain theory) [1,8]. However, typically these high order effects do not depend on external magnetic field, so that one might disregard them as possible mechanisms causing these anomalies. In this work, aiming to find a plausible explanation of the anomalies observed in ultrasound propagation at low applied magnetic field, we performed a theoretical analysis of this problem using linear theory of magnetoelasticity. To verify and extend our study, we also used three https://doi.org/10.1016/j.rinp.2025.108264 Received 15 January 2025; Received in revised form 1 April 2025; Accepted 14 April 2025 Results in Physics 73 (2025) 108264 Available online 28 April 2025 2211-3797/© 2025 Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ). I. Korniienko et al. computational approaches based on Finite Element Method (FEM), atomistic spin-lattice simulations, and effective corrections to elastic constants. Model Theory In this section, we present some details about the theory of magnetoelastic waves. Due to its complexity, it is convenient to present the linear theory where Simon effect 𝑆(𝐻) comes from. Higher order effects are described in Appendix ‘‘Higher order magnetoelastic effects and magnetic domain propagation’’. Linear magnetoelastic effects We consider a ferromagnetic material with cubic symmetry in an external magnetic field. The energy density of an inhomogeneously magnetized ferromagnet consists of the following components 𝐸=𝐸𝑒𝑥 +𝐸𝑎+𝐸𝑚𝑒 +𝐸𝑒𝑙 −𝐌𝐇,(1) where 𝐌 is the magnetization of the crystal. The first term of the right-hand side of the sum describes the contribution due to exchange interactions [12,13] 𝐸𝑒𝑥 =𝐴((∇𝛼𝑥)2+ (∇𝛼𝑦)2+ (∇𝛼𝑧)2),(2) written in terms of magnetization directional cosines 𝛼𝑖=𝑀𝑖∕𝑀𝑠 where 𝑀𝑠 is the saturation magnetization and 𝐴 is the exchange stiffness constant. Anisotropy energy density can be written as [14,15] 𝐸𝑎=𝐾1(𝛼2 𝑥𝛼2 𝑦+𝛼2 𝑦𝛼2 𝑧+𝛼2 𝑧𝛼2 𝑥) + 𝐾2(𝛼2 𝑥𝛼2 𝑦𝛼2 𝑧),(3) where 𝐾1, 𝐾2 are the magnetic anisotropy constants for the cubic ferromagnet. The magnetoelastic energy density in Eq. (1) is given by [11] 𝐸𝑚𝑒 =𝑏0(𝜖𝑥𝑥 +𝜖𝑦𝑦 +𝜖𝑧𝑧) + 𝑏1(𝛼2 𝑥𝜖𝑥𝑥 +𝛼2 𝑦𝜖𝑦𝑦 +𝛼2 𝑧𝜖𝑧𝑧)+ + 2𝑏2(𝛼𝑥𝛼𝑦𝜖𝑥𝑦 +𝛼𝑦𝛼𝑧𝜖𝑦𝑧 +𝛼𝑥𝛼𝑧𝜖𝑥𝑧)(4) where the constants 𝑏0, 𝑏1 and 𝑏2 are the magnetoelastic coupling constants, 𝜖𝑖𝑗 are the components of the strain tensor 𝜖. The fourth term in Eq. (1) is the energy density of elastic deformations [16]: 𝐸𝑒𝑙 =𝐶11 2(𝜖2 𝑥𝑥 +𝜖2 𝑦𝑦 +𝜖2 𝑧𝑧) + 𝐶12(𝜖𝑥𝑥𝜖𝑦𝑦 +𝜖𝑦𝑦𝜖𝑧𝑧 +𝜖𝑧𝑧𝜖𝑥𝑥)+ + 2𝐶44(𝜖2 𝑥𝑦 +𝜖2 𝑦𝑧 +𝜖2 𝑧𝑥), (5) where 𝐶11, 𝐶12, 𝐶44 are the second order elastic moduli. The last term in Eq. (1) is the energy density due to the Zeeman effect (applied field). The first step in the analysis of these systems is to find the equilibrium strain configuration from the energy local minima 𝜕𝐸 𝜕𝜖𝑖𝑗 = 0.(6) The solutions to these six equations give the components of the magnetostrictive strain tensor with respect to the equilibrium magnetization direction 𝜶0: 𝜖𝑥𝑥0= − 𝑏1 𝐶11 −𝐶12 (𝛼2 𝑥0−1 3)−1 3 3𝑏0+𝑏1 𝐶11 + 2𝐶12 , 𝜖𝑦𝑦0= − 𝑏1 𝐶11 −𝐶12 (𝛼2 𝑦0−1 3)−1 3 3𝑏0+𝑏1 𝐶11 + 2𝐶12 , 𝜖𝑧𝑧0= − 𝑏1 𝐶11 −𝐶12 (𝛼2 𝑧0−1 3)−1 3 3𝑏0+𝑏1 𝐶11 + 2𝐶12 , 𝜖𝑥𝑦0= − 𝑏2 2𝐶44 𝛼𝑥0𝛼𝑦0, 𝜖𝑦𝑧0= − 𝑏2 2𝐶44 𝛼𝑦0𝛼𝑧0, 𝜖𝑧𝑥0= − 𝑏2 2𝐶44 𝛼𝑧0𝛼𝑥0. (7) In order to investigate the effect associated with magneto-acoustic interaction, we use the system of equations for the magnetization and elastic displacement dynamics. Thus, the magnetic dynamics can be described by the Landau–Lifshitz–Gilbert equation 𝜕𝐌 𝜕𝑡 = −𝛾(𝐌×𝐇eff) − 𝛾𝛼𝐺 𝑀𝑠 𝐌×[𝐌×𝐇eff],(8) where 𝛾 is the gyromagnetic ratio, 𝐇eff = −𝜕𝐸∕𝜕𝐌 is the effective magnetic field, 𝛼𝐺 is the Gilbert damping parameter. For simplicity, in the following theoretical analysis we neglect the part related to damping, but note that the contribution of this term is analyzed in detail in the Refs. [17,18]. The elastic dynamics can be described by equation of motion for the displacement vector 𝐮 [16,19]: 𝜌𝜕2𝐮 𝜕𝑡2= ∇𝜎+𝐟𝑚𝑒,(9) where 𝜌 is the density of the crystal material, 𝜎 is the stress tensor with components 𝜎𝑖𝑗 related to the fourth-order elastic stiffness tensor with components 𝑐𝑖𝑗𝑘𝑙 and the second-order strain tensor with 𝜖𝑖𝑗 through the generalized Hooke’s law 𝜎𝑖𝑗 =𝑐𝑖𝑗𝑘𝑙𝜖𝑘𝑙.(10) To facilitate the manipulation of this equation it is convenient to replace 𝑐𝑖𝑗𝑘𝑙 by 𝐶𝑛𝑚 contracting a pair of cartesian indices into a single integer using Voigt notation: 𝑥𝑥 →1, 𝑦𝑦 →2, 𝑧𝑧 →3, 𝑦𝑧 →4, 𝑧𝑥 →5 and 𝑥𝑦 →6. The contribution from magnetic order for cubic crystal comes into the elastic wave equation in the form of effective force as 𝐟𝑚𝑒 = ∇ 𝛿𝐸𝑚𝑒 𝛿𝜖 =⎛⎜⎜⎜⎜⎝ 𝑏1 𝜕𝛼2 𝑥 𝜕𝑥 +𝑏2(𝜕𝛼𝑥𝛼𝑦 𝜕𝑦 +𝜕𝛼𝑥𝛼𝑧 𝜕𝑧 ) 𝑏1 𝜕𝛼2 𝑦 𝜕𝑦 +𝑏2(𝜕𝛼𝑥𝛼𝑦 𝜕𝑥 +𝜕𝛼𝑦𝛼𝑧 𝜕𝑧 ) 𝑏1 𝜕𝛼2 𝑧 𝜕𝑧 +𝑏2(𝜕𝛼𝑧𝛼𝑦 𝜕𝑦 +𝜕𝛼𝑥𝛼𝑧 𝜕𝑥 ) ⎞⎟⎟⎟⎟⎠ .(11) We consider small deviations of atoms and magnetization from their equilibrium positions 𝑢𝑗=𝑢𝑗0+𝑢𝑗𝑒𝑖(𝐤𝐫−𝜔𝑡), 𝛼𝑗=𝛼𝑗0+𝛼𝑗𝑒𝑖(𝐤𝐫−𝜔𝑡),(12) with the angular frequency 𝜔 and the wave vector of the collective wave 𝐤. The components of the strain tensor can also be written as the sum of a homogeneous part with small deviations 𝜖𝑙𝑚 =1 2(𝜕𝑢𝑙 𝜕𝑟𝑚 +𝜕𝑢𝑚 𝜕𝑟𝑙)= =𝜖𝑙𝑚0+𝑖 2(𝑘𝑚𝑢𝑙𝑒𝑖(𝐤𝐫−𝜔𝑡)+𝑘𝑙𝑢𝑚𝑒𝑖(𝐤𝐫−𝜔𝑡)). (13) It is easy to derive the dispersion relations for free acoustic waves from the Eq. (9) by considering only the elastic contribution (i.e. neglecting 𝐟𝑚𝑒). In the case where a magnetoacoustic interaction exists, these elastic waves can no longer propagate separately; under certain conditions, they will interact with the oscillations in the magnetic moment of the crystal [20]. To analytically study the Simon effect 𝑆(𝐻), we will analyze dispersion relations 𝜔(𝑘) of coupled magnetoacoustic waves for different values of the external magnetic field 𝐇 applied along [100] direction starting from the Eqs. (8) and (9). The bcc Fe, as an example of cubic material, for spin-lattice simulations and quantitative estimates is consider. Since the [100] is easy directions for single crystals of Fe [21], we may assume that the equilibrium value of 𝐌 will be also directed along this direction, i.e 𝜶0= (1,0,0) and as a result 𝜖𝑥𝑦0=𝜖𝑦𝑧0=𝜖𝑥𝑧0= 0, 𝜖𝑧𝑧0−𝜖𝑥𝑥0=𝜖𝑦𝑦0−𝜖𝑥𝑥0=𝑏1 𝐶11−𝐶12 . For given magnetic field and equilibrium magnetization orientations the system of dynamical equations (8), (9) after neglecting terms Results in Physics 73 (2025) 108264 2 I. Korniienko et al. quadratic in small amplitudes of deviations 𝛼𝑗 and 𝑢𝑗, leads to the following linearized system of equations: 𝑖𝜔𝛼𝑦=𝛾 𝑀𝑠 (𝛼𝑧(2𝑏2 1 𝐶11 −𝐶12 + 2𝐴𝑘2+ 2𝐾1+𝐻𝑀𝑠)+ +𝑏2(𝑖𝑘𝑧𝑢𝑥+𝑖𝑘𝑥𝑢𝑧)), −𝑖𝜔𝛼𝑧=𝛾 𝑀𝑠 (𝛼𝑦(2𝑏2 1 𝐶11 −𝐶12 + 2𝐴𝑘2+ 2𝐾1+𝐻𝑀𝑠)+ +𝑏2(𝑖𝑘𝑦𝑢𝑥+𝑖𝑘𝑥𝑢𝑦)), 𝜌𝜔2𝑢𝑥=𝐶11𝑘2 𝑥𝑢𝑥+𝐶12𝑘𝑥(𝑘𝑧𝑢𝑧+𝑘𝑦𝑢𝑦)+ +𝐶44(𝑘2 𝑧𝑢𝑥+𝑘𝑧𝑘𝑥𝑢𝑧+𝑘𝑥𝑘𝑦𝑢𝑦+𝑘2 𝑦𝑢𝑥)− −𝑏2(𝑖𝑘𝑧𝛼𝑧+𝑖𝑘𝑦𝛼𝑦), 𝜌𝜔2𝑢𝑦=𝐶11𝑘2 𝑦𝑢𝑦+𝐶12𝑘𝑦(𝑘𝑧𝑢𝑧+𝑘𝑥𝑢𝑥)+ +𝐶44(𝑘2 𝑧𝑢𝑦+𝑘𝑧𝑘𝑦𝑢𝑧+𝑘𝑥𝑘𝑦𝑢𝑥+𝑘2 𝑥𝑢𝑦)− −𝑏2𝑖𝑘𝑥𝛼𝑦, 𝜌𝜔2𝑢𝑧=𝐶11𝑘2 𝑧𝑢𝑧+𝐶12𝑘𝑧(𝑘𝑦𝑢𝑦+𝑘𝑥𝑢𝑥)+ +𝐶44(𝑘2 𝑦𝑢𝑧+𝑘𝑧𝑘𝑦𝑢𝑦+𝑘𝑥𝑘𝑧𝑢𝑥+𝑘2 𝑥𝑢𝑧)− −𝑏2𝑖𝑘𝑥𝛼𝑧. (14) As will be shown in Section ‘‘Results’’, from this theory it is possible derive the Simon effect 𝑆(𝐻) [9]. Simulations The general theory of coupled magnetoelastic waves has been known for a long time [12,19] and it has been constantly developed in an attempt to achieve greater accuracy [15,22]. In many cases, the analytical solution of the Eqs. (8) and (9) requires a number of simplifications but still can give relevant insights about this phenomenon. Such theoretical analysis can be supported and extended with the help of simulations. In this work, we compare the analytical solution of Eqs. (8) and (9) with three computational approaches: (i) Finite Element Method (FEM), (ii) atomistic spin-lattice simulations, and (iii) effective corrections to elastic constants. Here, we use the FEM to study this problem for the first time, while the other two approaches were already used in work [11]. Finite element method The FEM is based on continuum theory, and allows us to solve the partial differential equations of coupled elastic and magnetic dynamics Eqs. (8) and (9) numerically by subdividing a large system into smaller, simpler parts called finite elements. By solving Eqs. (8) and (9) means that in this approach we restrict our selves to the linear theory of magnetoelastic effects (Section ‘‘Linear magnetoelastic effects’’), so that we do not include high order magnetoelastic effects described in Appendix ‘‘Higher order magnetoelastic effects and magnetic domain propagation’’. These calculations are performed using COMSOL Multiphysics software (version 5.6) [23,24] combining the Solid Mechanic module with the Micromagnetic module developed by Weichao Yu et al. [25–27]. In Appendix ‘‘Magnetoelastic coupling in FEM’’ we provide details of the magnetoelastic coupling implementation in the FEM simulations. The frequency of the sound wave is extracted from the oscillation of the kinetic energy of a plane wave [11] 𝐮(𝐫, 𝑡) = 𝐮0cos (𝐤𝐫) cos (𝜔𝑡)(15) where 𝐮 is the displacement vector, 𝐮0 is the displacement amplitude, 𝜔 is the angular frequency of the sound wave, 𝐫 is the position and 𝑡 is the time. For such a plane wave, it turns out that the kinetic energy 𝐸𝑘𝑖𝑛 evolves with time as [11] 𝐸𝑘𝑖𝑛(𝑡) = 𝛬[1 − cos(4𝜋𝑓𝑡)],(16) where 𝛬 is a time-independent constant since we do not include energy dissipation, and 𝑓=𝜔∕2𝜋 is the frequency of the plane wave. Thus, Fig. 1. (Left) FEM mesh used in COMSOL to solve numerically coupled elastic and magnetic dynamics, which consists of 88 domain elements, 100 boundary elements, and 48 edge elements. (Right) The X-component of the displacement field (𝑢𝑥) for the initial state (𝑡= 0) of a transverse wave propagated along 𝐤∥ [001] with polarization 𝐮∥ [100]. 𝑓 can be easily extracted by fitting the simulated data of 𝐸𝑘𝑖𝑛 versus time to Eq. (16). In this method, we need to use periodic boundary conditions and the value of the wave number 𝑘= 2𝜋𝑛∕𝐿 where 𝑛 is an integer, and 𝐿 is the length of the simulated system along the 𝐤, see Fig. 1. For example, in the study of magnetic field effect (Section ‘‘Influence of magnetic field’’), we apply the parameter values 𝑢0= 0.01𝑎0, 𝑛= 1, 𝐿= 120𝑎0 and the equilibrium lattice parameter for magnetized iron 𝑎0= 2.83023 Å, leading to 𝑘𝑎0∕2𝜋= 0.0083, as in the atomistic spin-lattice model in Ref. [11]. For this geometry size, the FEM mesh consists of 88 domain elements, 100 boundary elements, and 48 edge elements, see Fig. 1, where we verified a good convergence of the presented results with the number of elements in the mesh. In the study of the influence of wave vector 𝑘 (Section ‘‘Influence of wave vector’’), the geometry length 𝐿 is modified accordingly. To calculate the kinetic energy as a function of time, we performed a FEM time-dependent study using the linear system solver SPOOLES with preordering algorithm Nested dissection, and pivot threshold equal to 0.01. We used the generalized-𝛼 method for time stepping [28] with initial time step 𝑑𝑡 = 10 fs, including linear predictor and Backward Euler consistent initialization, as implemented in COMSOL Multiphysics software [23,24]. Concerning material properties, we set the exchange stiffness constant 𝐴= 21 pJ/m according to Ref. [29], while the other necessary constants and parameters for the magnetized (i.e. collinear) state of bcc Fe are shown in Table 1. It is important to note that the applicability of FEM approach to study magnetoelastic waves is limited to long-wavelength spin waves (low 𝑘 values), where the discretization length (𝛥) of the micromagnetic mesh should be small enough to resolve spatial spin wave profile but include enough atoms to be valid as the continuous approximation. For computational studies of short-wavelength spin waves, one can use atomistic approaches, as described in Section ‘‘Atomistic spin-lattice simulations’’. Atomistic spin-lattice simulations The computational study of magnetic field effect on sound wave can be also done at the atomic scale by means of spin-lattice simulations, which allows to resolve coupled atom motion and magnetic moments dynamics. In this work, we use the same atomistic spin-lattice model of bcc Fe as described in our previous publication [11] using the SPIN Results in Physics 73 (2025) 108264 3 I. Korniienko et al. Table 1 Elastic constants (𝐶𝑖𝑗 ) for collinear oriented atomic magnetic moments, magnetoelastic constants (𝑏𝑖), magnetostrictive coefficients (𝜆), the magnetic anisotropy constant (𝐾1), saturation magnetization (𝑀𝑠), and density (𝜌) for bcc Fe at zero temperature [11]. 𝐶𝑖𝑗 (GPa) 𝑏(MPa) 𝜆(×10−6)𝐾1(kJ/m3)𝜇0𝑀𝑠(T) 𝜌(kg/m3) 𝐶11 252.56 𝑏1−4.41 𝜆001 26.08 55 2.26 8180 𝐶12 139.87 𝑏29.73 𝜆111 −30.33 𝐶44 106.95 Note that in table we use SI units as required for simulation in COMSOL although our theoretical formulas use CGS units. package of LAMMPS [30,31]. The model details are described in the Appendix ‘‘Atomistic spin-lattice model of bcc Fe’’ and the resulting properties of this model are shown in Table 1. The frequency of the sound wave is also extracted from the oscillation of the kinetic energy, as in the FEM approach. Higher order corrections arising from morphic coefficients 𝐺(𝑚) and finite strain theory (rotational-magnetostrictive effect) 𝑅(𝜆), described in Appendix ‘‘Higher order magnetoelastic effects and magnetic domain propagation’’ are naturally included in this type of simulations due to its atomic nature, but note that such effects are not included in the theoretical analysis presented in Section ‘‘Theory’’ and FEM simulations, based on infinitesimal strain theory. Hence, this fact should be taken into account for a fair comparison between these approaches, see Section ‘‘Results’’. The applicability of atomistic study based on oscillation of the kinetic energy is restricted to spin waves with non-extremely long wavelength, since simulating long wavelength requires very large atomic supercell which makes this task very computationally demanding. In such cases, the FEM approach could be more convenient and efficient. Effective corrections to elastic constants In this approach, we combine analytical formulas of group velocity for a general anisotropic media derived by Wang et al. [32] from elastic constants (𝐶𝑖𝑗 ) with magnetic corrections to elastic constants (𝛥𝐶𝑖𝑗 ) obtained by Rinaldi and Turilli [33]. Namely, the calculation of group velocity 𝑣 for any crystal symmetry is performed by calculating Christoffel matrix elements (𝛤), their partial derivative with respect to the phase velocity direction, and inserting them in analytical expression for group velocity derived by Wang et al. [32] 𝑣𝑙(𝐶𝑖𝑗 ) = 𝜕𝑐 𝜕𝑛𝑙 =𝜕𝜔 𝜕𝑘𝑙 , 𝑙 =𝑥, 𝑦, 𝑧 (17) where 𝑐=𝜔∕𝑘 is the phase velocity, and 𝒏=𝒌∕𝑘 is the unit vector of phase velocity direction (ray direction). On the other hand, the group velocity of a sound wave due to an external magnetic field can be obtained from effective magnetic corrections to elastic tensor (𝛥𝐶𝑖𝑗 ) [1,8,11,33]. Such corrections were derived by Rinaldi and Turilli and have the form [33] 𝛥𝐶𝑖𝑗 =1 𝑀2 𝑠 𝑏𝑝𝑞𝑖𝑏𝑟𝑠𝑗 𝛼𝑝𝛼𝑟𝜒𝑞𝑠,(18) where 𝑏𝑖𝑗𝑘 are magnetoelastic constants, 𝑀𝑠 is saturation magnetization, 𝛼𝑖 is direction cosine of magnetization at equilibrium, and 𝜒 is the susceptibility tensor, see Ref. [33] for details about subscript indexes. The corresponding fractional change in group velocity due to magnetic effects (𝑣−𝑣0)∕𝑣0 is obtained by computing Wang et al. formulas [32] using the elastic tensor including effective magnetic corrections (𝐶𝑖𝑗 + 𝛥𝐶𝑖𝑗 ), that leads to 𝑣, and elastic tensor without effective magnetic corrections (𝐶𝑖𝑗 ) which gives 𝑣0, that is 𝑣−𝑣0 𝑣0 =𝑣(𝐶𝑖𝑗 +𝛥𝐶𝑖𝑗 ) − 𝑣(𝐶𝑖𝑗 ) 𝑣(𝐶𝑖𝑗 ).(19) Note also that Rinaldi and Turilli expressions were derived using linear theory of magnetoelastic effects described in Section ‘‘Linear magnetoelastic effects’’, so that they do not include higher order corrections like morphic and rotational-magnetostrictive effects mentioned in Appendix ‘‘Higher order magnetoelastic effects and magnetic domain propagation’’ [11]. Moreover, Rinaldi and Turilli did not include the exchange stiffness term 𝐸𝑒𝑥, Eq. (2), in the total magnetic energy density [33]. We perform these calculations using the program VelCrys [34], where this approach has been recently implemented. If one is interested in the fractional change in frequency instead of group velocity, then Eq. (19) can be converted to (𝑓−𝑓0)∕𝑓0 using Eq. (A.1). Results Influence of magnetic field From the analysis of Eqs. (14), it follows that the case of longitudinal wave propagation along 𝐤∥ [100] with polarization 𝐮∥𝐇∥ [100] gives separate elastic wave with longitudinal sound velocity 𝑣𝑒∥=√𝐶11 𝜌= 5556.5 (m∕s),(20) and magnetic wave with frequency 𝜔𝑚=𝛾(2𝐴𝑘2 𝑀𝑠 +2𝐾1 𝑀𝑠 +2𝑏2 1 𝑀𝑠(𝐶11 −𝐶12)+𝐻)= =𝛾(𝐻𝑒𝑥 +𝐻𝐾+𝐻). (21) In other words, a longitudinal acoustic wave is not affected by oscillations of magnetic moments and, accordingly, the Simon effect is absent for such vectors 𝐤, 𝐮 and 𝐇 directions. Here we used the notation 𝐻𝑒𝑥 =2𝐴𝑘2 𝑀𝑠 and 𝐻𝐾=2𝐾1 𝑀𝑠 +2𝑏2 1 𝑀𝑠(𝐶11−𝐶12) to make Eq. (21) more compact and convenient. Further analysis shows that in the case of transverse wave propagation with 𝐤⟂𝐮∥𝜶 or 𝐤∥𝜶⟂𝐮, on the contrary, the magnetoelastic terms connecting the magnetic and elastic degrees of freedom in Eqs. (14) do not vanish and thus coupled magnetoelastic waves arise. This happens, for example, in the case of wave propagation along 𝐤∥ [001] with polarization 𝐮∥𝐇∥ [100] giving the following equations of coupled waves: 𝑖𝜔𝑀𝑠 𝛾𝛼𝑦− (2𝐴𝑘2+ 2𝐾1+2𝑏2 1 𝐶11 −𝐶12 +𝐻𝑀𝑠)𝛼𝑧−𝑏2𝑖𝑘𝑢𝑥= 0, 𝑖𝜔𝑀𝑠 𝛾𝛼𝑧+ (2𝐴𝑘2+ 2𝐾1+2𝑏2 1 𝐶11 −𝐶12 +𝐻𝑀𝑠)𝛼𝑦= 0, (𝜌𝜔2−𝐶44𝑘2)𝑢𝑥+𝑏2𝑖𝑘𝛼𝑧= 0. (22) The dispersion relation corresponding to these equations takes the form (𝜔2−𝜔2 𝑒)(𝜔2−𝜔2 𝑚) − 𝛾𝜔𝑚𝑏2 2 𝑀𝑠𝜌𝑘2= 0,(23) and has two magnetoelastic wave modes 𝜔2 ±=𝜔𝑒2+𝜔𝑚2 2±1 2√(𝜔𝑒2−𝜔𝑚2)2+ 4 𝛾𝜔𝑚𝑏2 2 𝑀𝑠𝜌𝑘2,(24) where 𝜔𝑒 is pure transverse elastic wave frequency 𝜔𝑒=𝑘√𝐶44 𝜌,(25) with sound velocity 𝑣𝑒⟂=√𝐶44 𝜌= 3615.9 (m/s) and 𝜔𝑚 is pure spin wave frequency given by Eq. (21). On Fig. 2 we show phonon–magnon dispersion along the [001] (𝛤-H) line of 𝑘-points in first Brillouin zone for the dimensionless Results in Physics 73 (2025) 108264 4 I. Korniienko et al. Fig. 2. Dispersion relation for magnetoelastic waves in a ferromagnet calculated using Eq. (24) and parameters from Table 1. The blue dotted and dashed lines denote the unperturbed magnonic (𝜔𝑚) and phononic (𝜔𝑒) dispersion relations, respectively, while the solid lines represent the coupled magnetoelastic waves. The gray dotted line corresponds to value 𝑘𝑎0∕2𝜋= 0.008333 used for further investigation of Simon effect in the FEM and spin-lattice simulations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) parameter 𝑘𝑎0∕2𝜋 in ranges from 0 to 0.05, where the value 0 of the parameter corresponds to the 𝛤-point and the value 1 corresponds to the H-point. Now lets mark 4𝛾𝜔𝑚𝑏2 2 𝑀𝑠𝜌𝑘2 in Eq. (24) as 𝛿2 and suggested that 𝛿2< (𝜔𝑒2−𝜔𝑚2)2. Here we note that although such assumption is relevant for the most cases, it fails in close proximity to resonance where the frequencies 𝜔𝑒 and 𝜔𝑚 almost coincide. Then using Maclaurin series and approximation √1 + 𝑥= 1 + 𝑥 2−𝑥2 8+⋯≈ 1 + 𝑥 2 (for |𝑥|<1) one can get 𝜔2 ±=1 2(𝜔𝑒2+𝜔𝑚2± (𝜔𝑒2−𝜔𝑚2)) ± 1 4 𝛿2 (𝜔𝑒2−𝜔𝑚2).(26) The ‘‘−’’ sign gives us frequency 𝜔2 −=𝜔𝑚2−1 4 𝛿2 (𝜔𝑒2−𝜔𝑚2),(27) which after further simplification takes a form 𝜔−≈𝜔𝑚(1 − 1 8 𝛿2 𝜔𝑚2(𝜔𝑒2−𝜔𝑚2)).(28) It is easy to see that it follows to quasi-magnetic wave mode and therefore it is beyond the scope of the interest of this paper. Next, let us consider the ‘‘+’’ sign case: 𝜔2 +=𝜔𝑒2+1 4 𝛿2 (𝜔𝑒2−𝜔𝑚2),(29) giving us the approximate frequency of the quasi-sound wave mode 𝜔+≈𝜔𝑒(1 + 1 8 𝛿2 𝜔𝑒2(𝜔𝑒2−𝜔𝑚2)).(30) Since the first-order approximation linear-magnetoelastic theory does not allow to include contributions from the fractional changes of frequency arising from higher order correction (i.e. morphic, rotational and magnetostrictive effects), Eq. (A.2), the Simon effect can be simply extracted from angular frequencies as 𝑆(𝐻) = 𝜔+−𝜔𝑒 𝜔𝑒 =𝑏2 2𝛾𝜔𝑚𝑘2 2𝑀𝑠𝜌𝜔𝑒2(𝜔𝑒2−𝜔𝑚2).(31) Thus, the Simon effect formula taking into account the dispersion contributions (which were lost in the original Simon effect formula [9,10] Fig. 3. Simon effect 𝑆(𝐻) of a transverse wave propagated along 𝐤∥ [001] with 𝑘𝑎0∕2𝜋= 0.0083, polarization 𝐮∥ [100] and applied magnetic field 𝐇∥ [100]. Comparison between theory (obtained using 𝜔± from Eq. (24)) with simplified formula Eq. (32) (‘‘New formula S(H)’’), Simon effect formula S(H) Eq. (35), and simulation using FEM, atomistic spin-lattice simulation and calculations using effective corrections to elastic constants (Eqs. (18) and (19)). In the atomistic spin-lattice simulations, we originally computed the fractional change in group velocity 𝛥𝑣∕𝑣0 [11], so that to plot here the net contribution of Simon effect 𝑆(𝐻) we subtracted the contributions responsible for other effects (𝑙−𝑙0 𝑙0 +𝐺(𝑚) + 𝑅(𝜆)=8.27 × 10−5) from the obtained in simulation data, according to Eqs. (A.1) and (A.2). Similarly, we subtracted the magnetostrictive contribution 𝑙−𝑙0 𝑙0 = −1.3×10−5 from the data obtained by correcting the elastic constants since we originally computed the fractional change in group velocity Eq. (19), and in this approach higher order corrections are not included [33]. as expected to have no noticeable effect) takes the form 𝑆(𝐻) = 𝑏2 2 2𝑀𝑠𝐶44 𝛾2(𝐻𝐾+𝐻𝑒𝑥 +𝐻) (𝑘2𝐶44 𝜌−𝛾2(𝐻𝐾+𝐻𝑒𝑥 +𝐻)2).(32) The advantage of the refined formula is greater accuracy in the lowfield regime. In particular, the Eq. (32) shows the more correct place of decline of the curve 𝑆(𝐻) when the field decreases thus giving the value of the resonance field as 𝐻𝑟𝑒𝑠 =𝑘 𝛾√𝐶44 𝜌−1 𝑀𝑠(2𝑏2 1 𝐶11 −𝐶12 + 2𝐾1+ 2𝐴𝑘2),(33) revealing the magnitude of the applied magnetic field where an anomaly may be observed in the fractional change in frequency. In case of significantly large applied external magnetic field 𝐻 for which 𝜔𝑚2≫ 𝜔𝑒2, after neglecting of 𝜔𝑒2 in the denominator we can get 𝑆(𝐻) ≈ − 𝛾𝑏2 2 2𝐶44𝑀𝑠𝜔𝑚 = − 𝑏2 2 2𝐶44𝑀𝑠(𝐻𝐾+𝐻𝑒𝑥 +𝐻),(34) which is very similar to the already known Simon effect formula [9,10] for proposed directions of 𝐇,𝐤 and 𝐮 vectors: 𝑆(𝐻) = − 𝑏2 2 2𝐶44𝑀𝑠(2𝐾1 𝑀𝑠 +𝐻) .(35) In this way, the refined formula Eq. (32) turns into Simon’s formula Eq. (35) for fields 𝐻 significantly larger than 𝐻𝑟𝑒𝑠. In Fig. 3 we compare the analytical methods to determine the Simon effect and the three simulation approaches described in Section ‘‘Simulations’’. It should be noted that spin-lattice simulations takes into account terms of higher orders in the energy density expression and thus, in addition to the Simon effect for the considered directions of vectors (𝐤∥ [001], 𝐮∥𝐇∥ [100]), it also takes into account Eq. (A.3) 𝐺(𝑚) = −𝑚2,𝛾 3∕12𝐶44 with 𝑚2,𝛾 3= −14.3 MPa, 𝑅(𝜆) = 3(𝜆001 −𝜆111)∕2, as well as magnetostrictive effect 𝛥𝑙∕𝑙0= −𝜆001∕2 [11]. Thus, in order to get data consistent with data obtained by other methods based on linear Results in Physics 73 (2025) 108264 5 I. Korniienko et al. Fig. 4. The relative difference between S(H) from Eq. (32) (𝑆𝑛𝑒𝑤) and Eq. (35) (𝑆𝑆𝑖𝑚𝑜𝑛) in bcc Fe for a transverse wave propagated along 𝐤∥ [001] with 𝑘𝑎0∕2𝜋= 0.008333, polarization 𝐮∥ [100] and applied magnetic field 𝐇∥ [100]. theory of magnetoelasticity, the difference given by these effects should be subtracted from the data in the case of the atomistic spin-lattice simulation, according to Eqs. (A.1) and (A.2). As can be seen from Fig. 3, in the region of large fields (𝐻 ≫ 𝐻𝑟𝑒𝑠 ≈ 29.4 kOe), all methods, both analytical and based on numerical simulations, give quite similar result. However, with a decrease in the magnitude of the applied magnetic field, the situation begins to change dramatically. Thus, the lowest accuracy is shown by the approximate formula obtained by Simon (purple dashed line) given by Eq. (35), which was unable to correctly represent the region with a sharp change in the speed of the sound wave. A slightly better result is shown by the method based on the magnetic corrections of elastic constants used in VelCrys (orange line with stars). The specified formula from Eq. (32) (blue dotted line), despite all the simplifications, has practically the same accuracy with direct calculation of the effect from the frequencies of coupled magnetoelastic oscillations as (𝜔±−𝜔𝑒)∕𝜔𝑒 with Eqs. (24), (21) and (25) (lilac and green broad solid lines). Thus, the consideration of dispersion effects associated with exchange stiffness significantly improves the accuracy of the 𝑆(𝐻) formula describing the Simon effect in the low field regime. The advantage of the new formula is that it (as well as the direct calculation from the frequency expressions) gives good agreement with numerical simulations not only in the region of large (𝐻 ≫ 𝐻𝑟𝑒𝑠) but also in region of small fields (𝐻 ≪ 𝐻𝑟𝑒𝑠). And we draw attention to the fact that for the calculation of the resonance value of the applied field 𝐻𝑟𝑒𝑠, it is also important to take into account dispersion contributions (see Eq. (33)), the neglect of which can lead to a significant deviation from the actual value, as it happens in the presented example for bcc Fe. For greater clarity of the role of taking into account terms related to dispersion in Fig. 4 we show how the relative difference between S(H) from Eqs. (32) and (35) varies with the magnetic field. In the vicinity of resonance (𝐻→𝐻𝑟𝑒𝑠), all approximate analytical methods give untrustworthy and questionable result with extremely large value of |𝛥𝑓∕𝑓0| (asymptotic behavior). A possible explanation could be that in the region closer to resonance the Simon effect cannot be considered only in the approximation of linearized Eq. (14) and should contain corrections corresponding to higher orders of 𝛼𝑗 and 𝑢𝑗 in them. In contrast to approximate analytical expressions, the simulations show quite plausible results with smaller and finite fractional change in frequency, but results still vary depending on the program in which the simulation was carried out (teal line with circles and violet line with squares). The observed difference between atomistic simulations and FEM might be related to numerical effects or/and physical reasons. For example, FEM approach becomes more accurate as wave vector 𝑘 is decreased where the continuum approximation holds, see Fig. 5. Fig. 5. The role of dispersion in the Simon effect of a transverse wave propagated along 𝐤∥ [001] with polarization 𝐮∥𝐇∥ [100] at constant field 𝐻= 29383 Oe. Comparison of the analytical solution 𝑆(𝑘) given by the Eq. (32), with analytical expressions based on the Eqs. (24), and FEM, atomistic spin-lattice numerical simulation. In the atomistic spin-lattice simulations, to plot the net contribution of Simon effect we subtracted the contributions responsible for other effects (𝑙−𝑙0 𝑙0 +𝐺(𝑚)+𝑅(𝜆) = 8.27×10−5), as in Fig. 3.. In order to further verify the proposed approach and more fully highlight the role of dispersion effect, in Appendix ‘‘Comparison with experimental data for Co-Pt alloy’’ we present a comparison of the proposed theoretical model with the available experimental data for the Co-Pt alloy. Influence of wave vector One more interesting aspect, which naturally follows from the proposed model, is that resonance-related anomalies in the speed of sound waves (analogous to those occurring at 𝐻→𝐻𝑟𝑒𝑠) can also appear for a constant field 𝐻=𝑐𝑜𝑛𝑠𝑡 at certain values of the wave vectors 𝑘→𝑘𝑟𝑒𝑠 (see Figs. 5and 6). This is explained by the fact that the resonance conditions, namely 𝜔𝑒=𝜔𝑚, in Eqs. (31) and (32) can be achieved not only due to the choice of the external magnetic field, but also due to the choice of the wave vector 𝑘. In other words 𝑆(𝐻) in Eq. (32) should be read as 𝑆(𝑘, 𝐻) since the influence of 𝑘 can be significant under certain conditions. In Fig. 5, we show an example of the expected dependence of the value 𝑆(𝑘, 𝐻) on the wave vector 𝑘 from the analytical formulas Eqs. (24) and (32) in the presence of a constant field 𝐻=29383 Oe. Comparison of this dependence with data from FEM and atomistic spin-lattice simulations shows a good agreement in the region of small values of 𝑘 and the correct position of the first resonance at 𝑘𝑎0∕2𝜋≈ 0.0083. As 𝑘 increases, the FEM simulations give null Simon effect (𝑆≈ 0), which means that it cannot account for the magnetoelastic effects on sound waves due to the small size of the required FEM geometries beyond the continuum approximation. In contrast to FEM, the data of atomistic simulations are in better agreement with the analytical results and, in particular, show the presence of a second resonance in an approximately expected location. Small shift of the resonance location in the atomistic simulation can be explained by the fact that its position is mainly determined by the ratio of the saturation magnetization to the exchange constant 𝑀𝑠∕𝐴, which may slightly differ in the simulation and the analytical formula. Discussion These results reveal the existence of anomalies in Simon effect, its wave vector dependency. A more complete picture of possible variations is presented in Fig. 6 which shows the dependence of the Simon effect value 𝑆(𝑘, 𝐻) given by Eq. (32) both on the external magnetic Results in Physics 73 (2025) 108264 6 I. Korniienko et al. Fig. 6. Dispersion dependence of the Simon effect for bcc Fe. The influence on the 𝑆 value for a transverse wave with 𝐤∥ [001] and 𝐮∥𝐇∥ [100] both the magnetic field and the wave vector are shown. The yellow line corresponds to the resonant condition 𝜔𝑒=𝜔𝑚 (written in the form of Eq. (33)). Near the resonance (dark regions) anomalies on the fractional change in sound wave frequency are expected. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) field and on the parameter associated with the wave vector. The dark green and dark violet areas on the graph correspond to the areas near the resonance (the condition for the appearance of the coincidence of frequencies 𝜔𝑒=𝜔𝑚 is shown by the yellow line). As expected, these areas are places where resonance-induced anomalies in the propagation speed of sound waves in ferromagnets can occur. Thus, the largest field at which it is possible to fulfill the condition of resonance is 𝐻𝑚𝑎𝑥 =𝑀𝑠𝐶44 8𝐴𝜌𝛾2−1 𝑀𝑠(2𝑏2 1 𝐶11 −𝐶12 + 2𝐾1),(36) which in the case of bcc Fe has a value of 𝐻𝑚𝑎𝑥 = 44527 Oe. The corresponding to 𝐻𝑚𝑎𝑥 value of the wave vector at the point of resonance is 𝑘𝑚𝑎𝑥 =𝑀𝑠 4𝐴𝛾 √𝐶44 𝜌,(37) and for case of iron 𝑘𝑚𝑎𝑥𝑎0∕2𝜋≈ 0.0198. Certain anomalies in the frequency of the propagating sound wave are also possible in fields slightly larger than 𝐻𝑚𝑎𝑥 due to proximity to resonance (see Fig. 6), but fully realized resonance-induced acoustical anomalies can occur only in fields 𝐻≤𝐻𝑚𝑎𝑥. We note that here we considered an example of coupled magnetoelastic oscillations and analyzed their role in the Simon effect. However, it is known about the strong dependence of the magnitude of the Simon effect on the mutual orientation of the vectors 𝐤, 𝐮 and 𝐌, so the effect we have studied may vary in case of arbitrary orientations of these vectors. Conclusions Although the role of the exchange interaction in the dispersion dependencies of magnetoelastic waves is known and well studied, it is common practice to neglect it when describing the effect of a magnetic field on the fractional change in sound wave velocity. Such approach leads to the disappearance of dispersion terms (i.e. dependence on the wave vector 𝐤) in the resulting analytical expressions and simplifies them. Convenient for ultrasonic measurements in high magnetic fields, mentioned approach can, however, fail in the low-field regime. In our work, we identified the analytical formula to describe the Simon effect in cubic ferromagnets by introducing in it all the necessary dispersion terms. This modification improved the accuracy of Simon’s formula almost to the level of direct calculation of the fractional change in frequency from the analytical expressions of the coupled magnetoelastic oscillations frequencies. As a result, it well describes the change in the speed of the sound wave in the crystal both for field values significantly higher and lower than the resonant value. However, the description of the change in the sound velocity under conditions close to resonance still remains an open question and a goal for further research. The presented analysis provides a plausible explanation of the observed peculiarities in Simon effect at low magnetic fields reported in previous works [10,11]. We hope that our results, in particular, will be useful to scientists dealing with the study of surface acoustic waves in ferromagnetic materials for fast and simple estimations of magnetic field effect. CRediT authorship contribution statement I. Korniienko: Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Data curation. P. Nieves: Writing – review & editing, Supervision, Investigation, Conceptualization. D. Legut: Writing – review & editing, Validation, Supervision, Resources, Funding acquisition, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This work was supported by projects e-INFRA CZ (ID:90254) and QM4ST (CZ.02.01.01/00/22_008/0004572) by The Ministry of Education, Youth and Sports of the Czech Republic and also by Czech Science Foundation by grant No. 22-35410K. P. N. acknowledges support by grant MU-23-BG22/00168 funded by The Ministry of Universities of Spain. Appendix A. Higher order magnetoelastic effects and magnetic domain propagation The magnetoelastic effects on ultrasound wave propagation are typically studied through the fractional change in group velocity (𝑣−𝑣0)∕𝑣0 or wave frequency (𝑓−𝑓0)∕𝑓0, which are related to each other via the following expression [10,11,35] 𝑣−𝑣0 𝑣0 =𝑓−𝑓0 𝑓0 +𝑙−𝑙0 𝑙0 ,(A.1) where (𝑙−𝑙0)∕𝑙0 is the fractional change in length along wave propagation direction. In addition to Simon effect 𝑆(𝐻), which is field dependent and arises from linear magnetoelastic effects (Section ‘‘Linear magnetoelastic effects’’), there are higher order magnetoelastic effects influencing the fractional change in frequency [10,11,35] 𝑓−𝑓0 𝑓0 =𝑆(𝐻) + 𝐺(𝑚) + 𝑅(𝜆) + 𝛥𝑜𝑡ℎ𝑒𝑟𝑠,(A.2) where 𝐺(𝑚) accounts for contribution from the morphic coefficients 𝑚, arising from second order magnetoelastic coupling combined with anharmonic third order elastic constants, and is not explicitly dependent on the magnetic field strength but depends on the magnetization orientation with respect to crystallographic axes. The term 𝑅(𝜆) comes from rotational and magnetostrictive effects that could not be predicted by the linear theory of elasticity but are derived from the rotational invariance and finite strain theory. Lastly, 𝛥𝑜𝑡ℎ𝑒𝑟𝑠 corresponds to other possible effects like 𝛥𝑠-effect due to domain wall propagation which Results in Physics 73 (2025) 108264 7 I. Korniienko et al. is field dependent below technical saturation and difficult to characterize since depends on the zero-field demagnetized state or remanent state [10]. In this work, we focus on transverse sound waves propagating along 𝐤∥ [001] with polarization 𝐮∥ [100] and applied field 𝐇∥ [100], where elastic and magnetic parts are strongly coupled, and anomalies have been observed on the fractional change in the frequency at low applied magnetic field [10,11]. For this case, the high order corrections in Eq. (A.2) for cubic crystals take the form [10,11,35] 𝐺(𝑚)||||| 𝒌∥[001] 𝒖∥𝑴∥[100] = − 𝑚𝛾,2 3 12𝐶44 , 𝑅(𝜆)||||| 𝒌∥[001] 𝒖∥𝑴∥[100] =3 2(𝜆001 −𝜆111),(A.3) where 𝑚𝛾,2 3 is a morphic coefficient as defined in Refs. [10,11,35], while 𝜆001(111) are anisotropic magnetostrictive coefficients. Note that these high order corrections are included in atomistic spin-lattice simulations due to the atomic nature of this approach [11]. Finally, in this case the fractional change in length reads [11] 𝑙−𝑙0 𝑙0||||| 𝒌∥[001] 𝒖∥𝑴∥[100] =𝑙−𝑙0 𝑙0||||| 𝜷∥[001] 𝜶∥[100] = − 𝜆001 2,(A.4) where 𝜷 is the length direction. Appendix B. Magnetoelastic coupling in FEM The FEM simulations in present work are performed using COMSOL Multiphysics software [23] where the magnetoelastic effects can be calculated by combining the Solid Mechanic module with the Micromagnetic module developed by Weichao Yu et al. [25–27]. Thus, the elastic effect is taken into account in the Micromagnetic module by adding an additional term 𝐇𝑚𝑒 eff = − 𝜕𝐸𝑚𝑒 𝜕𝐌,(B.1) in the effective field 𝐇eff of the Landau–Lifshitz–Gilbert equation Eq. (8), while magnetic effect is included in the Solid Mechanic module by an additional term in the applied force to the solid 𝐟𝑚𝑒 given by Eq. (11), as body load. In the Solid Mechanic module, the dynamics of lattice displacement 𝐮 is governed by equation of motion for elastic waves given by Eq. (9). The mechanical force at the boundaries of system is set 𝐟𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦 𝑚𝑒 = −𝐧⋅⎛⎜⎜⎝ 2𝑏1𝛼2 𝑥𝑏2𝛼𝑥𝛼𝑦𝑏2𝛼𝑥𝛼𝑧 𝑏2𝛼𝑥𝛼𝑦2𝑏1𝛼2 𝑦𝑏2𝛼𝑦𝛼𝑧 𝑏2𝛼𝑥𝛼𝑧𝑏2𝛼𝑧𝛼𝑦2𝑏1𝛼2 𝑧⎞⎟⎟⎠ ,(B.2) where 𝐧 is the norm vector on the boundary, which is derived from Lagrangian equations. Appendix C. Atomistic spin-lattice model of bcc Fe For the atomistic spin-lattice simulations we consider the following Hamiltonian 𝑠𝑙(𝐫,𝐩,𝐬) = 𝑚𝑎𝑔 (𝐫,𝐬) + 𝑁 ∑ 𝑖=1 |𝐩𝑖|2 2𝑚𝑖 + 𝑁 ∑ 𝑖,𝑗=1 (𝑟𝑖𝑗 ),(C.1) where 𝐫𝑖, 𝐩𝑖, 𝐬𝑖, and 𝑚𝑖 stand for the position, momentum, normalized magnetic moment and mass for each atom 𝑖 in the system, respectively, (𝑟𝑖𝑗 ) = (|𝐫𝑖−𝐫𝑗|) is the interatomic potential energy and 𝑁 is the total number of atoms in the system with total volume 𝑉. For the classical interatomic potential (𝑟𝑖𝑗 ), we use the spectral neighbor analysis potential (SNAP) [58] for bcc Fe developed by Nikolov et al. [36]. The magnetic energy includes the Zeeman term, exchange interaction and the Néel interaction: 𝑚𝑎𝑔(𝐫,𝐬) = −𝜇0 𝑁 ∑ 𝑖=1 𝜇𝑖𝐇𝐬𝑖−1 2 𝑁 ∑ 𝑖,𝑗=1,𝑖≠𝑗 𝐽(𝑟𝑖𝑗 )𝐬𝑖⋅𝐬𝑗+𝑁 𝑒𝑒𝑙(𝐫,𝐬),(C.2) where 𝜇𝑖 is the atomic magnetic moment, 𝜇0 is the vacuum permeability, 𝐇 is the external magnetic field, 𝐽(𝑟𝑖𝑗 ) is the exchange parameter. The Néel interaction term reads [1,11,37] 𝑁 𝑒𝑒𝑙 = − 1 2 𝑁 ∑ 𝑖,𝑗=1 {𝑔(𝑟𝑖𝑗 ) + 𝑙1(𝑟𝑖𝑗 )[(𝒆𝑖𝑗 ⋅𝒔𝑖)(𝒆𝑖𝑗 ⋅𝒔𝑗) − 𝒔𝑖⋅𝒔𝑗 3] +𝑞1(𝑟𝑖𝑗 )[(𝒆𝑖𝑗 ⋅𝒔𝑖)2−𝒔𝑖⋅𝒔𝑗 3][(𝒆𝑖𝑗 ⋅𝒔𝑗)2−𝒔𝑖⋅𝒔𝑗 3] +𝑞2(𝑟𝑖𝑗 )[(𝒆𝑖𝑗 ⋅𝒔𝑖)(𝒆𝑖𝑗 ⋅𝒔𝑗)3+ (𝒆𝑖𝑗 ⋅𝒔𝑗)(𝒆𝑖𝑗 ⋅𝒔𝑖)3]}, (C.3) where 𝒆𝑖𝑗 =𝒓𝑖𝑗 ∕𝑟𝑖𝑗 , and 𝑙1(𝑟𝑖𝑗 ) = 𝑙(𝑟𝑖𝑗 ) + 12 35 𝑞(𝑟𝑖𝑗 ), 𝑞1(𝑟𝑖𝑗 ) = 9 5𝑞(𝑟𝑖𝑗 ), 𝑞2(𝑟𝑖𝑗 ) = − 2 5𝑞(𝑟𝑖𝑗 ). (C.4) Following the methodology proposed in [37] we set 𝑔(𝑟𝑖𝑗 ) = −𝐽(𝑟𝑖𝑗 ), which leads to an offset of the exchange energy. The spatial dependence of 𝐽(𝑟𝑖𝑗 ), 𝑙(𝑟𝑖𝑗 ) and 𝑞(𝑟𝑖𝑗 ) is described using the Bethe–Slater curve, as implemented in the SPIN package of LAMMPS [31] 𝐽(𝑟𝑖𝑗 ) = 4𝛼𝐽(𝑟𝑖𝑗 𝛿𝐽)2[1 − 𝛾𝐽(𝑟𝑖𝑗 𝛿𝐽)2]𝑒−(𝑟𝑖𝑗 𝛿𝐽)2 𝛩(𝑅𝑐,𝐽 −𝑟𝑖𝑗 ), 𝑙(𝑟𝑖𝑗 ) = 4𝛼𝑙(𝑟𝑖𝑗 𝛿𝑙)2[1 − 𝛾𝑙(𝑟𝑖𝑗 𝛿𝑙)2]𝑒−(𝑟𝑖𝑗 𝛿𝑙)2 𝛩(𝑅𝑐,𝑙 −𝑟𝑖𝑗 ), 𝑞(𝑟𝑖𝑗 ) = 4𝛼𝑞(𝑟𝑖𝑗 𝛿𝑞)2[1 − 𝛾𝑞(𝑟𝑖𝑗 𝛿𝑞)2]𝑒−(𝑟𝑖𝑗 𝛿𝑞)2 𝛩(𝑅𝑐,𝑞 −𝑟𝑖𝑗 ), (C.5) where 𝛩(𝑅𝑐,𝑛−𝑟𝑖𝑗 ) is the Heaviside step function and the 𝑅𝑐,𝑛 (𝑛=𝐽, 𝑙, 𝑞) are the cut-off radii. The parameters 𝛼𝑛, 𝛾𝑛, 𝛿𝑛 must be determined in order to reproduce the Curie temperature, magnetostriction and magnetocrystalline anisotropy. Analytical expressions and the method of calculating these parameters are explained in detail in [37]. Here for bcc Fe we follow the parameters proposed in Ref. [11], thus taking 𝑅𝑐,𝐽 =𝑅𝑐,𝑙 =𝑅𝑐,𝑞 = 2.6Å, 𝛿𝐽=𝛿𝑙=𝛿𝑞= 2.45105 Å, 𝛼𝐽= −14.2048 meV/atom, 𝛼𝑙= 377.32 μeV/atom, 𝛼𝑞= 29.965 μeV/atom, 𝛾𝐽= 2.6125, 𝛾𝑙= 0.78979, 𝛾𝑞= 1.0496. The dynamics of the coupled spins and atoms can be obtained by integrating the following equations of motion [31]: 𝑑𝐫𝑖 𝑑𝑡 =𝐩𝑖 𝑚𝑖 , 𝑑𝐩𝑖 𝑑𝑡 = − 𝜕𝑠𝑙 𝜕𝐫𝑖 , 𝑑𝐬𝑖 𝑑𝑡 = − 1 ℏ 𝜕𝑚𝑎𝑔 𝜕𝐬𝑖 ×𝐬𝑖, (C.6) with ℏ= 6.582 × 10−7 eV/(rad⋅THz) the reduced Planck constant. Appendix D. Comparison with experimental data for Co-Pt alloy In order to further verify the analytical expression proposed in this paper, we will use it to compare with the available experimental data for external field dependence of the ultrasonic echoes for the transverse mode (𝐤∥ [001], 𝐮∥ [100]) in Co-Pt alloy [10]. We note that in the mentioned experimental paper only the direction of the wave vector 𝐤 is indicated and there are no data on its magnitude, therefore we take the value that best fits the data, i.e. 𝑘𝑎0∕2𝜋≈ 0.00131 for both sets of experimental data (with 𝐇∥ [001] and 𝐇∥ [100]), although this parameter could vary and be different in individual experiments. It should also be noted that in order to avoid the cumbersomeness of the formulas, we left the minimum number of interactions in Eq. (1), neglecting the demagnetizing effect and dipolar field associated with the stress, therefore the resulting Eq. (32) does not distinguish between the [100] and [001] external field orientations for mentioned 𝐤 and 𝐮 Results in Physics 73 (2025) 108264 8 I. Korniienko et al. directions. Thus, to properly describe an experiment where such effects have a significant impact, corresponding to them interactions need to be taken into consideration. They can be included in the effective magnetic field easily as an additional constant fields since the system is assumed to be saturated and lead to the appearance of corresponding terms in 𝜔𝑚 Eq. (21) and, as result, in 𝑆(𝐻). The fractional change in frequency (𝑓−𝑓0)∕𝑓0 that observed in experiment contain higher order magnetoelastic effects and other possible effects like 𝛥𝑠-effect as it mentioned in Eq. (A.2) and to obtain the field-dependent contributions responsible for Simon effect, these higher order effects must be subtracted from the experimental data. Thus, for the case of 𝐮∥ [100], 𝐤∥ [001], 𝐇∥ [001] the higher order corrections are 𝐺(𝑚) = −0.000517, 𝑅(𝜆) = 0.0003945 and other contributions including such as 𝛥𝑠-effect, that we estimate from experimental curve, 𝛥𝑜𝑡ℎ𝑒𝑟𝑠 ≈ 0.00437. The analytical expression for describing the phenomenon of interest to us is given by Rouchy et al. [10] as 𝑆(𝐻) = − 𝑏2 2 2𝐶44𝑀𝑠(2𝐾1 𝑀𝑠 +𝐻𝐷+𝐻),(D.1) where 𝐻𝐷= −𝑁𝑀𝑠 is demagnetizing field with demagnetizing factor 𝑁=4𝜋 3 for a sphere of 4.94 mm in diameter as it was used in the experiment. The proposed in this paper refined formula including dispersion effects Eq. (32) for such case will take the form 𝑆(𝐻) = 𝑏2 2 2𝑀𝑠𝐶44 𝛾2(𝐻𝐾+2𝐴𝑘2 𝑀𝑠 +𝐻𝐷+𝐻) (𝑘2𝐶44 𝜌−𝛾2(𝐻𝐾+2𝐴𝑘2 𝑀𝑠 +𝐻𝐷+𝐻)2) ,(D.2) where 𝐻𝐾=2𝐾1 𝑀𝑠 +2𝑏2 1 𝑀𝑠(𝐶11−𝐶12). In Fig. 7(a) we show a comparison of the Simon formula Eq. (D.1), the refined Simon formula Eq. (D.2), and experiment. As can be seen, both formulas give equally good agreement with experiment in the high-field region, but in the low-field region near resonance, both are not sufficiently accurate. More experimental data in this regime will be highly desirable in order to clarify the behavior of Simon effect. The better result we get in case of 𝐮∥ [100], 𝐤∥ [001], 𝐇∥ [100] (see Fig. 7(b)) where practically the entire experimental curve can be explained by a refined analytical expression Eq. (D.4). Since for transverse modes, the field dependence is shifted when the magnetic field is perpendicular to the propagation and this arises because of the dynamical dipolar field associated with the stress [35], the analytical formula for Simon effect given by Rouchy et al. contain additional 4𝜋𝑀𝑠 term [10]: 𝑆(𝐻) = − 𝑏2 2 2𝐶44𝑀𝑠(2𝐾1 𝑀𝑠 + 4𝜋𝑀𝑠+𝐻𝐷+𝐻),(D.3) which can also be included in our refined formula as 𝑆(𝐻) = 𝑏2 2 2𝑀𝑠𝐶44 𝛾2(𝐻𝐾+2𝐴𝑘2 𝑀𝑠 + 4𝜋𝑀𝑠+𝐻𝐷+𝐻) (𝑘2𝐶44 𝜌−𝛾2(𝐻𝐾+2𝐴𝑘2 𝑀𝑠 + 4𝜋𝑀𝑠+𝐻𝐷+𝐻)2) .(D.4) The corrections arising from higher-order magnetoelastic effects and which should be subtracted from the experimental data in this case are 𝐺(𝑚) = −0.000517, 𝑅(𝜆) = −0.0003945 and other contributions 𝛥𝑜𝑡ℎ𝑒𝑟𝑠 ≈ 0.00465. For calculations we used following material parameters of Co-Pt alloy: 𝑎0=3.749 Å [38], 𝐴= 10−6 erg cm−1 [39], 𝐶11 = 2.897 × 1012 erg cm−3, 𝐶12 = 1.785 × 1012 erg cm−3, 𝐶44 = 1.241 × 1012 erg cm−3, 𝑏1= −3.5 × 108 erg cm−3, 𝑏2= 1.19 × 108 erg cm−3, 𝐾1= −4 × 105 erg cm−3, 𝑀𝑠= 749.6 emu cm−3, 𝑚𝛾,2 3= 7.7 × 109 erg cm−3 [10], 𝜆001 = 210 × 10−6, 𝜆111 = −53 × 10−6 [34]. Data availability We will make the data publicly available. Fig. 7. The field-dependent contributions to fractional change in frequency (𝑓−𝑓0)∕𝑓0. (a) Comparison with the experiment for the case 𝐮∥ [100], 𝐤∥ [001], 𝐇∥ [001]. ‘‘Simon formula’’ is given by Eq. (D.1) and ‘‘New formula’’ is Eq. (D.2). The vertical dotted lines show the values of the resonant fields, that follow from the Eqs. (D.1) and (D.2), respectively. (b) Case of 𝐮∥ [100], 𝐤∥ [001], 𝐇∥ [100]. ‘‘Simon formula’’ is given by Eq. (D.3) and ‘‘New formula’’ is given by Eq. (D.4). Experimental curves are taken from Ref. [10]. References [1] Du Trémolet de Lacheisserie E. Magnetostriction theory and applications of magnetoelasticity. CRC-Press; 1993. [2] Reid AH, Shen X, Maldonado P, Chase T, Jal E, Granitzka PW, Carva K, Li RK, Li J, Wu L, Vecchione T, Liu T, Chen Z, Higley DJ, Hartmann N, Coffee R, Wu J, Dakovski GL, Schlotter WF, Ohldag H, Takahashi YK, Mehta V, Hellwig O, Fry A, Zhu Y, Cao J, Fullerton EE, Stöhr J, Oppeneer PM, Wang XJ, Dürr HA. Beyond a phenomenological description of magnetostriction. Nat Commun 2018;9(1):388. http://dx.doi.org/10.1038/s41467-017-02730-7. [3] Danilevich A, L’vov V. Strong influence of ferromagnetic ordering and internal pressure on the elastic modulus of shape memory alloy. J Magn Magn Mater 2013;333:108–13. http://dx.doi.org/10.1016/j.jmmm.2012.12.021, URL https:// www.sciencedirect.com/science/article/pii/S030488531200995X. [4] Li X, Labanowski D, Salahuddin S, Lynch CS. Spin wave generation by surface acoustic waves. J Appl Phys 2017;122(4):043904. http://dx.doi.org/10.1063/1. 4996102. [5] Puebla J, Xu M, Rana B, Yamamoto K, Maekawa S, Otani Y. Acoustic ferromagnetic resonance and spin pumping induced by surface acoustic waves. J Phys D: Appl Phys 2020;53(26):264002. http://dx.doi.org/10.1088/1361-6463/ab7efe. [6] Gomonay HV, Kornienko EG, Loktev VM. Magnetoelastic coupling as a source of shape dependence of afmr spectra. Ukr J Phys 2005;50(8):816–22, URL http://archive.ujp.bitp.kiev.ua/files/journals/50/8/500810p.pdf. [7] du Trémolet de Lacheisserie E. Magnetoelastic properties of amorphous alloys. J Magn Magn Mater 1982;25(3):251–70. http://dx.doi.org/10.1016/ 0304-8853(82)90072-5, URL https://www.sciencedirect.com/science/article/pii/ 0304885382900725. [8] Del Moral A. Handbook of magnetostriction and magnetostrictive materials. Zaragoza: Del Moral Zaragoza; 2008. [9] Simon G. Das verhalten von ebenen ultraschallwellen in homogen magnetisierten einkristallen. Z Für Naturforschung A 1958;13(2):84–9. http://dx.doi.org/10. 1515/zna-1958-0204. [10] Rouchy J, de Lacheisserie ET, Genna J, Waintal A. A 6 T ultrasonic study of the magnetoelastic coupling in ferromagnetic Co—Pt alloy. J Magn Magn Mater 1980;21(1):69–79. http://dx.doi.org/10.1016/0304-8853(80)904898, URL https://www.sciencedirect.com/science/article/pii/0304885380904898. [11] Nieves P, Tranchida J, Nikolov S, Fraile A, Legut D. Atomistic simulations of magnetoelastic effects on sound velocity. Phys Rev B 2022;105:134430. http://dx.doi.org/10.1103/PhysRevB.105.134430, URL https://link.aps.org/doi/ 10.1103/PhysRevB.105.134430. Results in Physics 73 (2025) 108264 9