Finite difference magnetoelastic simulator

We describe an extension of the micromagnetic finite difference simulation software MuMax3 to solve elasto-magneto-dynamical problems. The new module allows for numerical simulations of magnetization and displacement dynamics in magnetostrictive materials and structures, including both direct and inverse magnetostriction. The theoretical background is introduced, and the implementation of the extension is discussed. The magnetoelastic extension of MuMax3 is freely available under the GNU General Public License v3.

Magnetoelectrics are a class of spintronic materials that have received particular attention due to their potential to manipulate magnets at very low power and the possibility for new functionalities [20][21][22][23][24] .Besides multiferroic materials [25][26][27][28][29] , which possess simultaneous magnetic and ferroelectric polarization, magnetoelectric composites have been at the center of research on magnetoelectricity for decades 20,23,[29][30][31][32][33][34] .Such composites consist of piezoelectric and magnetostrictive materials and rely on the interaction between the magnetization and mechanical degrees of freedom.Applying an electric field to a piezoelectric material leads to the generation of a mechanical deformation (strain) that can be transferred to an adjacent magnetostrictive material.The strain inside the (ferromagnetic) magnetostrictive material leads then to a magnetoelastic effective magnetic field via the Villari effect (inverse magnetostriction) that can exert a torque on the magnetization.Conversely, the rotation of the magnetization in a magnetostrictive material generates strain that can lead to charge separation and an electric polarization in an adjacent piezoelectric material.This indirect coupling between electric fields and magnetization is schematically represented in Figure 1.This paper (and the software described therein) addresses the magnetoelastic interaction between strain ε and magnetization M. For many applications, in particular in microelectronics, spintronic devices must be miniaturized to the nanoscale to be competitive.The integration of magneto-electric composites into scaled spintronic devices requires the understanding of the piezoelectric as well as the magnetoelastic interactions at the nanoscale.While the fundamental magnetoelectric and magnetoelastic coupling is well understood [35][36][37][38][39][40] , the detailed behavior of nanoscale structures has only very recently been addressed [41][42][43][44][45][46][47][48] .Nanoscale spintronic devices may possess complex structures with different functional areas and materials, leading to spatially varying properties and highly convoluted boundary conditions.It is clear that such devices cannot be represented well by (often approximate) analytical models and that numerical simulations are typically necessary.
While many commercial and free solutions exist to simulate nanoscale piezoelectric devices due to the importance of microelectromechanical systems (MEMS), no comprehensive free software package exists to date to simulate magnetoelastic interactions at the nanoscale.At microwave frequencies in the GHz range, which are often relevant for spintronic applications, the magnetization dynamics become nontrivial and micromagnetic simulators are required to describe and understand the device performance.While it is possible to include magnetoelastic effective fields in micromagnetic simulations using the open-source simulator OOMMF [49][50][51] , MuMax3 52 , or Boris 53 , such approaches cannot describe the backaction of the magnetodynamics on the elastodynamics of the system.Hence, a full description of the magnetoelastic dynamics of such a system requires solve simultaneously both the elastodynamic and magnetodynamic equations of motion.Whereas full magnetoelastic simulations of nanoscale structures have been published using Comsol ® , [54][55][56][57][58][59][60][61][62] the software modules developed for these simulations are not freely available for the scientific community.
In this paper, we present a magnetoelastic software module that allows for the simultaneous numerical solution of magneto-and elastodynamics in complex geometries, including both direct and inverse magnetostrictive interactions.The module is written as an extension of the widely used micromagnetic solver MuMax3 52 .Both the solver and the extension are available under the GNU General Public License v3.This approach allows for the usage of all pre-defined functionalities of MuMax3 in combination with the magnetoelastic module.Furthermore, simulations are accelerated using GPUs with respect to standard simulators that run on CPUs.

Theoretical background
The magnetoelastic software module described below extends the well-known micromagnetic solver MuMax3 and allows for the simulation of the magnetoelastodynamic behavior in magnetostrictive materials.In this section, the fundamental underlying equations are introduced, which are implemented in the solver.
The magnetization dynamics are described by with m = M / M s the normalized magnetization, M s the saturation magnetization, and T the magnetic torque.The dot denotes a time derivative.In the absence of spin currents, the torque is equal to the Landau-Lifshitz torque, given by with γ 0 = µ 0 γ, γ the absolute value of the gyromagnetic ratio, α the phenomenological damping constant, and H eff the effective magnetic field.This effective field contains all effects that influence the magnetization dynamics, including exchange and dipolar interactions.Equation (1) and Equation (2) form the base of the micromagnetic solver already implemented in MuMax3.
Effects of inverse magnetostriction can be inserted into the above framework by adding a magnetoelastic field H mel to the effective field H eff in Equation (2).For materials with cubic (or higher) crystal symmetry, the magnetoelastic effective field is given by 63-65 with B 1 and B 2 the first and second magnetoelastic coupling constants, respectively, and ε ij the components of the mechanical strain tensor � ε .
The elastodynamics in a solid are described by the second order differential equation 66,67 eff , or equivalently by two first order differential equations with u the displacement, v the velocity, a the acceleration, ρ the mass density, η a phenomenological damping parameter, and f eff the effective body force.This body force contains all effects that influence the displacement dynamics, such as elastic or magnetoelastic contributions.These two contributions are the ones considered in this extension.The elastic body force is given by f el = ∇σ , with σ the mechanical stress tensor.For small displacement values, Hooke's law ˆ= σ c ε is valid.Here, ĉ is the fourth order stiffness tensor.Together with the definition of the strain ( ) the elastic body force for materials with cubic crystal structure (or higher symmetry) can be written as @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ 12 44 x y y x z y z @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ Here, c ij are the components of the stiffness tensor in reduced dimensionality, i.e. in Voigt notation.
In magnetostrictive materials, there is also a body force contribution due to the magnetostriction effect.For a magnetostrictive material with cubic crystal structure (or higher symmetry), the magnetoelastic body force is given by x y f @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ (8) The mechanical boundary conditions at the surface are with F s the traction force per unit surface and n the surface normal.

Methods
In this section, we describe the MuMax3 framework as well as the implementation of the magnetoelastic module.The integration of the module into MuMax3 has the benefit that all predefined MuMax3 functions can also be applied to the magnetoelastic extension.The first part of this section explains the numerical implementation of the elastodynamic equations as well as the elastic body force.This is followed by an overview of all novel functions and parameters defined in the magnetoelastic extension together with instructions on how to use them.

MuMax3 framework
MuMax3 is an established finite difference micromagnetic solver.In the following, several important aspects of the MuMax3 framework as well as the procedures for running simulations and extracting data are briefly introduced.This introduction is meant to provide the necessary basis to properly use the magnetoelastic module in a later stage.
In MuMax3, simulations are often defined by a script with the .mx3extension that contains all information necessary to describe the system under study.The script syntax is based on a subset of the Go language and is comprehensively explained in the MuMax3 documentation 52,68 .The most important aspects are, however, briefly explained here.
The first part of the script contains typically information about the dimensions of the simulated system and the mesh size.The user sets the mesh and the dimensions by defining the cell size together with the number of cells in the three orthogonal directions.The cell shape is always a cuboid.Within the simulated system, it is possible to define different regions with different geometries, in which the material parameters can vary, and which can therefore represent different materials.
Next, material parameters need to be assigned to every region.Based on the provided parameters, MuMax3 automatically determines which interactions need to be considered in the simulation.Therefore, it is straightforward to include a wide range of different interactions by just specifying appropriate material parameters.Once the mesh with its regions and all material parameters are provided, the simulation can be run in two ways: the first option consists of a numerical energy minimization to determine the equilibrium ground state of the system.The second option consists of the time-dependent integration of the magnetization dynamics using Equation (1).This option allows for the assessment of the temporal evolution of the dynamic variables, in particular the magnetization dynamics.In this paper, we will focus on the second option.Several numerical integration schemes are readily implemented in MuMax3 and can be selected by the SetSolver(x) command, with x an integer between 1 and 6 68,69 .Once the desired solver is chosen, the simulation can be started by the command run(t_tot), with t_tot the total simulation time set by the user.
MuMax3 can also perform several mathematical operations on physical variables.For example, this can be useful when extracting and postprocessing the output data.Note that all mathematical operations are performed on GPUs, which typically results in short computational times.All variables can be extracted at any time during the simulation and can be converted to different output formats, such as .ovf,.vtk,or .csv.

Implementation of the magnetoelastic extension
To simulate the magnetoelastic dynamics, both the magnetic and elastic equations of motion, i.e.Equation (1) and Equation (4), respectively, must be solved simultaneously.For the implementation, the two first order differential equations in Equation ( 5) were chosen instead of the the single second order differential equation in Equation ( 4) for stability reasons.Hence, the system of equations that must be solved becomes which can be written in a compact manner as Note that the magnetic torque T and the effective body force f eff include the magnetoelastic field and the magnetoelastic body force, respectively, to account for the mutual interactions between the elastic and magnetic domains.Hence, instead of solving a single differential equation, i.e.Equation ( 1), the magnetoelastic module must solve a set of differential equations, including elastic and magnetoelastic interaction terms, as a function of time.
Time integration.The time integration of Equation ( 13) has been numerically implemented via the fourth order Runge-Kutta (RK4) method given by ), 6 where the subscript i corresponds to the number of the time step with length ∆t.The k-parameters are given by To use this new solver instead of the regular MuMax3 solvers that integrate only Equation (1), the user must specify SetSolver(9) in the .mx3script.The regular MuMax3 solvers run by default with adaptive time stepping.By contrast, the novel solver only runs with a predefined fixed time step ∆t, which also needs to be defined by the user in the .mx3script.Typical time steps that give accurate results and lead to reasonable simulation times are of the order of ∆t ≈ 0.1 ps.The user can also always verify whether the chosen time step is accurate enough by running an additional simulation with a much smaller time step.If both simulations lead to identical results, the chosen time steps are sufficiently small.If the two results differ, the time step should be shortened further.

Calculation of R.
For every time step i, the right hand side of Equation ( 13), R(t, w), needs to be calculated several times.The vector R comprises the magnetic torque T, the velocity v , and the acceleration a.The numerical calculation of T is already implemented and optimized in the regular MuMax3 framework.Hence, the calculation of T is kept untouched in the magnetoelastic extension.The velocity vector v is obtained by the Runge-Kutta integration of Equation (12).Finally, the acceleration a is determined by calculating a = (f eff − ηv) /ρ using the previously found velocity v in combination with Equation (7) and Equation ( 8), whose sum is the effective body force.
In the following, the numerical implementation of the elastic body force will be explained in more detail since it is nontrivial in presence of free elastic boundary conditions.The elastic body force depends on second order spatial derivatives of the displacement, as indicated by Equation ( 6).These spatial derivatives are approximated using the finite difference method.In the bulk of the system, the derivatives are calculated based on the sum of a second order forward and second order backward finite difference using averaged stiffness constants.This is done to properly account for parameters that may vary between different positions (regions).
To illustrate this approach, the implemented equations to calculate the second order derivative of u (i, j) at position (i, j) with respect to x and the mixed derivative with respect to x and y are ( 1, 1) ( 1, 1) ( 1, 1) ( 1, 1) .4 The second derivative in the y-direction as well as the mixed derivative with respect to y and x directions are defined in an straightforward analogous way.For uniform stiffness constants, i.e. c(i, j) ≡ c, these equations reduce to the regular central finite difference scheme, i.e.

( , )
( 1, ) 2 ( , ) ( 1, ) and Boundary conditions.At the edges of a structure, elastic boundary conditions need to be considered, which alter the differential equations and thus their numerical implementation.Three different elastic boundary conditions are implemented in the magnetoelastic module: 1. Periodic boundary conditions, which connect the left edge to the right edge and the bottom to the top of the structure.These boundary conditions do not represent physical boundaries, and therefore Equation (19) and Equation (20) can be used at the edges as well.
2. Fixed boundary conditions, which also do not require modifications of the bulk finite difference scheme.They can be implemented by defining a region at the boundary with a fixed displacement value (e.g. of 0).
3. Free boundary conditions, which correspond to zero traction force at the surfaces.
The third type of free boundary conditions requires the reformulation of the finite difference scheme at the boundaries to satisfy Equation (9).For simplicity, the simulated structure is assumed to be uniform in the z-direction with boundaries only in the x-and y-directions.Moreover, we consider zero traction force.Then, the Neumann boundary conditions at the free surface are 0 0 0.
These expressions must be inserted into the elastic body force at the edges.The implementation of the equations at the mesh edge perpendicular to the x-direction is described below as an example.The same procedure can be used to derive expressions for the corresponding equations at the other interfaces.
The mesh near a boundary in the x-direction is visualized in Figure 2. The interface itself is located at positions (−1/2, j), which are the exact locations where the boundary conditions in Equation ( 23) are valid.Note that this boundary also describes the limits of the magnetization, i.e. the magnetic material ends at (−1/2, j).By contrast, the coordinates of the points next to the boundary edge can be written as (0, j).
For free boundary conditions, the differential equations need to be modified.In the following, the three components of the elastic body force are rewritten in their finite difference approximations taking into account The x-component of the elastic body force is , (0, ) (0, ) .= + xy xx @ j @ j @x @y σ σ f el x (24)   At location i = 0, the first term of Equation ( 24) is given by ( 1, 1) 2 ( , ) 2 ( 1, ) , 8 whereas the second term becomes Note that the boundary condition formally states that σ xy (−1/2, j) = 0, whereas this is approximated by σ xy (0, j) = 0 in the numerical implementation.A possible solution to circumvent this approximation is to use a staggered grid where one of the grids contains the displacement values and the other one the stress values.Another solution could be to implement an additional raster of mesh points on the outside of the "original" mesh that can be used to set the proper boundary conditions.In this extension, we have chosen to use the approximation σ xy (0, j) ≈ σ xy (−1/2, j) = 0 because the mesh is already defined and fixed for the magnetization.However, using another or an additional mesh would strongly complicate the problem since the magnetoelastic extension is built based on the MuMax3 framework, which is designed to work with one mesh only.
The y-component of the elastic body force is given by Considering again that the surface of the magnetic material is located at i = −1/2 and taking into account the boundary conditions given in Equation ( 23), the first term can be rewritten as x y @ @ @ @ @ @ ( ) and the second term as ( , 1) 2 ( , ) ( , 1) ( ( , 1) ( , 1) 2 ( , )) 4 x y i j i j i j u i j u i j u i j u u i j i j x x i j i j i j u i j u yy y x y y y x x y y y y y @ @ @ @ @ @ @ @ @ @ @ ( ) ( ) Finally, the z-component of the elastic body force is given by el, (0, ) (0, ) .
In this case, the first term can be rewritten as and the second term as .
Within MuMax3, all these mathematical operations are parallelized and performed on a GPU to reduce the computational time.

Operation of the magnetoelastic extension
The instructions to install the software are provided on the MuMax3 github page https://github.com/mumax/The software requires the installation of Go and the CUDA Toolkit as well as the availability of an NVIDIA GPU.
In MuMax3, several different data types exist, such as parameters, fields, excitations, etc., which can be defined, modified, or extracted via various methods 52,68 .This section presents an overview of additional material parameters, vector fields, excitations, energies, boundary conditions, and initial states that are defined in the magnetoelastic extension.Their usage is analogous to the usage of the equivalent elements defined in MuMax3 52,68 and will therefore be only briefly discussed.Note that all regular MuMax3 methods can also be applied to the new elements defined in the MuMax3 extension.

New material parameters.
The usage and the assignment of the newly defined material parameters are the same as for traditional MuMax3 material parameters.Details can be found in the MuMax3 documentation 68 .
Parameters can depend on time and-through the definition of regions-on position.The following list provides an overview of novel material parameters and their units, which are defined in the magnetoelastic extension of MuMax3: • • eta: Phenomenological elastic damping constant with unit [kgs −1 m −3 ]; • rho: Mass density with unit [kgm −3 ]; • B1: first magnetoelastic coupling constant with unit [Jm −3 ]; • B2: second magnetoelastic coupling constant with unit [Jm −3 ].
New vector fields.There are two vector field types in MuMax3.The first one is called a variable and can be treated in the same way as the magnetization.The second type is called a quantity and represents a field that can be derived from the fundamental variables.The usage of the new vector fields is the same as for the traditional MuMax3 68 .The novel variables and their units defined within the magnetoelastic extension of MuMax3 are: • u: elastic displacement vector with unit [m]; • du: velocity vector with unit [ms −1 ]; The novel quantities and their units defined within the magnetoelastic extension of MuMax3 are: • normStrain: vector that contains the normal strain components [ε xx , ε yy , ε zz ], calculated according to ( ) . Note that the strain corresponds to the real strain and not the engineering strain; • shearStrain: vector that contains the shear strain components [ε xy , ε yz , ε xz ], calculated according to ( ) . Again, note that the strain corresponds to the real strain and not the engineering strain; • normStress: vector that contains the normal stress components [σ xx , σ yy , σ zz ], calculated according to Hooke's law, ˆ= σ c ε, with unit [Nm −2 ]; • shearStress: vector that contains the shear stress components [σ xy , σ yz , σ xz ], calculated according to Hooke's law, with unit [Nm −2 ]; • poynting: elastic Poynting vector, calculated via • B_mel: magnetoelastic field µ 0 H mel with unit [T] implemented according to Equation (3); • F_mel: magnetoelastic body force with unit [Nm −3 ] implemented according to Equation (3).
The magnetoelastic MuMax3 extension can simulate magnetoelastodynamics for materials with cubic or higher symmetries.Hence, the module only contains the c 11 , c 12 , and c 44 stiffness constants.For the strain, the magnetoelastic field, and the magnetoelastic body force, standard central finite difference schemes are used to calculate the derivatives in the respective equations.
Excitations.Excitations, such as applied fields and current densities, can depend on time and space following the form f(t) × g(x, y, z).Here, g can be a continuously varying spatial profile.This is different from the material parameters, which need to be uniform within each region.The following gives an overview of novel excitations and their units within the magnetoelastic MuMax3 extension.Their usage is analogous to the regular MuMax3 usage for excitations 68 : • force_density: external body force f ext that is added to the effective body force • exx: ε xx strain component; • eyy: ε yy strain component; • ezz: ε zz strain component; • exz: ε xz strain component; • exy: ε xy strain component; • eyz: ε yz strain component.
Energies.Energies can be extracted at different moments in time to study the system.Every energy contribution can be extracted as either a field or a scalar value.The field corresponds to the energy density and the scalar object corresponds to the total energy of that contribution inside the system.The following list gives an overview of novel energies and their units defined in the magnetoelastic MuMax3 extension: • Edens_el: elastic energy density (unit [Jm −3 ]) in the linear regime where Hooke's law is valid, • E_el: the total elastic energy, i.e. el dV ∫ E , with unit [J]; • Edens_mel: magnetoelastic energy density (unit [Jm −3 ]) for a material with cubic (or higher) crystal symmetry, implemented via • E_mel: total magnetoelastic energy, i.e. • Edens_kin: kinetic energy density (unit [Jm −3 ]), implemented via • E_kin: total kinetic energy, i.e. kin dV Initial states.Three new initial states have been implemented for both magnetization and displacement.Their usage is equivalent to the standard MuMax3 usage for the definition of the initial state 68 .
• A Gaussian pulse with spherical distribution of the out-of-plane vector component: GaussianSpheri-cal_outplane(A, pos_x, pos_y, sig_x, sig_y float64), implemented as • with x 0 = pos_x, y 0 = pos_y, σ x = sig_x and σ y = sig_y; • A Gaussian pulse with spherical (symmetric) distributions of the in-plane vector components: GaussianSpherical(A, pos_x, pos_y, sig_x, sig_y, angle float64), implemented as • Uniformly oriented in-plane vector components along a specific direction with a Gaussian distribution along the transverse direction GaussianUniform(A, pos, sig, angle1, angle2 float64), implemented as with θ = angle2 in degrees, x 0 = cos(ϕ)x + sin(ϕ)y, ϕ = angle1, x 0 = pos and σ = sig.Boundary conditions.It is possible to apply several types of elastic boundary conditions in the new magnetoelastic MuMax3 module.All elastic boundary conditions correspond to interfaces in the x-and y-directions of the film since uniform displacement is assumed along the z-direction.
• Free boundary conditions, corresponding to zero forces at the edges of the film.These are the default boundaries in the x-and y-directions of the film.
• Periodic boundary conditions, which can be enabled in the same way as in the traditional MuMax3 framework via the SetPBC() command 68 .
• Absorbing boundaries, which can be implemented by defining regions with gradually increasing elastic damping η, similarly to the gradually increasing damping used to absorb spin waves in a purely magnetic system 70 .This is, for example, useful to prevent wave reflection at surfaces or interfaces.
Fixed boundary conditions for the displacement are also possible.Such boundary conditions can be obtained by defining regions at the edges of the mesh with fixed displacement values.The parameter frozenDispLoc defines a region with fixed displacement, and the parameter FrozenDispVal defines the fixed displacement value.This works in the same way as the regular frozenspins function in MuMax3, which fixes the spins in a specific region.However, spin (magnetic moment) has always the same norm, whereas the norm can vary for the displacement.Therefore, the additional function FrozenDispVal has been defined.Note that this value can be time dependent and thus can also act as an excitation source.

Conclusion
This paper describes an extension of the established finite difference micromagnetic solver MuMax3, which adds capabilities to calculate elastodynamics including the magnetoelastic coupling between mechanical and magnetic degrees of freedom.It therefore allows for the finite difference simulation of magnetoelastodynamics in 1D, 2D or quasi-3D systems.The implementation in MuMax3 means that all standard MuMax3 functionalities can also be used in the magnetoelastic extension, and that the mathematical operations can be performed on GPUs, resulting in low computational times.Multiple different boundary conditions have been implemented, specifically periodic, fixed, and free bounadry conditions, as described in Sec. .The software is freely available under the GNU General Public License v3.
The fundamental equations in which the code is based, its implementation in MuMax3 and considerations on the treatment of elastic boundary conditions are discussed with clarity and simplicity, making the paper understandable to non-specialists in simulations and mechanics of solids, among which I include myself.It is particularly welcome the section entitled Operation of the magnetoelastic extension, giving useful practical information and details for making it easier for new users to getting started in the new software, which is not always the case in the literature describing new simulation codes.In a similar vein, the authors provide the commented input script of the magnetoelastic module which was used to simulate a nanoscale magnetostrictive CoFeB waveguide for one of their recent publications.
Overall, this new module is poised to become a valuable support for research on a range of systems and applications in which magnetostriction can be relevant and thus the manuscript should be indexed on Open Research Europe.
However, I think that a short discussion on the limitations of the module would improve the manuscript.More specifically, i) The code has been implemented for cubic crystal symmetries.Since many relevant magnetic materials have lower symmetries than the cubic one it would be interesting a comment on the possible extensions to lower symmetries and in particular to systems with uniaxial anisotropy.ii) Since the full 3D elastodynamics description assumes of uniform displacement along the z direction, it would be helpful some discussion on the limitations of this assumption and on how to assess, using the module, how realistic this is for a particular model, if this is possible.
Other additional minor remarks: Page 4, end of line 16: "solve" should be replaced by "solving".At the end of page 7, mention is made to eq. 6.Since it is indeed eq. 7, a simplification of eq.6 considering uniform displacement along z, which is used in the code, to avoid confusion in the discussion of boundary conditions that follows it would be better to refer to eq. 7.
At the end of page 14 the text description accompanying the example of defining a time dependent displacement boundaries refers to "region 3" but it seems that the code is rather for a region 2… The paper repeatedly refers to "cubic (or higher) crystal symmetry".To my understanding this is a bit confusing because the cubic one is the highest of all the crystal symmetries.A material which has a higher symmetry than that of a cubic crystal must be polycrystalline and therefore it is not strictly a crystal but a polycrystalline material.
Is the rationale for developing the new software tool clearly explained?Yes

Is the description of the software tool technically sound? Yes
Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?Yes

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes
Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Synthesis of oxides, magnetic and structural characterization of oxides.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
The authors developed an open-source extension to the finite difference simulation software MuMax3 that allows simulating dynamic magneto-elastic interactions.MuMax3 is widely used for simulation of magnetization dynamics due to its open-source availability, ease of use and fast computation times.Magneto-elastic interactions are generally relevant for magnetization dynamics as most magnetic materials are magneto-elastic.Magneto-elasticity can enable magnon-phonon hybridization.This software extension should be a very powerful tool to simulate magneto-elastic interactions and hybrid magnon-phonon dynamics beyond the limits of simple analytical models.Previous open-source software implementations have only been able to simulate the effect of magneto-elastics on magnetization dynamics, but not the back-action of magnetization dynamics onto the elastics.Considering this back-action is however required for accurate description of magneto-elastic phenomena.The contribution is thus of high relevance and novelty for the field.
The description of the underlying physics and parameters used in the code is clear and the authors introduce the relevant background to motivate their work in a suitable manner.I found this to be a very high quality and accessible work that will be important in particular for scientists intending to employ and tailor the extension for their own problems.
As a minor suggestion for clarification, the authors should consider providing a short description of the limitations of the extension.The authors demonstrated in Ref. 71 that the extension can be used to simulate confined magneto-elastic waves in thin waveguides.It would be nice to briefly mention here if there are other geometries where their extension can (or cannot) be used.As an example, many researchers are interested in the interaction of surface acoustic waves (SAWs) and magnetization dynamics in magnetic thin films.The typical system thereby is a substrate / thin magnetic film bilayer.To me it is not clear to what capacity the extension can be used as is or modified to address magneto-elastic phenomena in this configuration.For SAW-spin wave interactions in laterally extended systems, is it possible to work around the restriction of uniform displacement in the z direction by taking a coordinate system with SAW propagation along x and interface normal along y (plane surface waves)?A major challenge might be accurate simulation of SAWs that dominantly propagate in the non-magnetic substrate but still interact with spin waves in the magnetic film.Can the authors comment on this challenge?
As a second point for potential future modifications, it would be interesting to discuss further mechanisms that can lead to a coupling of magnetization dynamics and strain even in the absence of magneto-elasticity.In this regard, Xu and coworkers 1 recently discussed magneto-rotation coupling which would add an additional term to Eq. ( 3) stemming from the rotation of the magnetic anisotropy axis due to shear strain.This coupling can be comparable or even larger than magneto-elastic coupling

Greg Carman
Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, Los Angeles, CA, USA I believe that the authors have done a wonderful job of expanding the availability of software to address the important applications of magnetoelastic problems.To my knowledge the only other approaches are housed in Finite Element based programs such as Comsol which requires the partial differential equations to be input AND the structure is fairly slow and limited in size of problems that can be solved.Below, I provide some comments the authors may wish to consider adding or commenting on in the manuscript.Can authors comment on how this code which I believe is FDTD based compares with FEM based approaches.I am hopeful that the authors know the answer and can provide comparison metrics so that others working in this area can make informed decisions on which numerical platform for this problem is superior. 1.
Authors reduce the full 3D problem to a subset with the assumption of zero strain in the z direction, eq 7. I am not sure why they do this or is it because the full 3D problem is to cumbersome to solve.The previous FEM solutions in the literature do not make this 2.
assumption.Adding a sentence or two describing why this assumption is made will help the community.

Is the rationale for developing the new software tool clearly explained? Yes
Is the description of the software tool technically sound?Yes Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others?Yes Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool?Yes Competing Interests: No competing interests were disclosed.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

Figure 1 .
Figure 1.Schematic of the magnetoelectric coupling in magnetoelectric composite materials, in which the coupling between dielectric polarization P and magnetization M occurs indirectly via strain ε.The software described here allows for the numerical simulation of the magnetoelastic coupling in complex structures.

Figure 2 .
Figure 2. Definition of the grid at a boundary.

©
2021 Carman G.This is an open access peer review report distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Yes Is the description of the software tool technically sound? Yes Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others? Yes Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes Competing Interests:
1,2.Is it in principle possible to extend the code to consider magnetorotation coupling?No competing interests were disclosed.Experimental solid state physics, magnetism, magnetization dynamics.

confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
Reviewer Report 05 May 2021 https://doi.org/10.21956/openreseurope.14372.r26749