Bi-directional coupling in strain-mediated multiferroic heterostructures with magnetic domains and domain wall motion

Strain-coupled multiferroic heterostructures provide a path to energy-efficient, voltage-controlled magnetic nanoscale devices, a region where current-based methods of magnetic control suffer from Ohmic dissipation. Growing interest in highly magnetoelastic materials, such as Terfenol-D, prompts a more accurate understanding of their magnetization behavior. To address this need, we simulate the strain-induced magnetization change with two modeling methods: the commonly used unidirectional model and the recently developed bidirectional model. Unidirectional models account for magnetoelastic effects only, while bidirectional models account for both magnetoelastic and magnetostrictive effects. We found unidirectional models are on par with bidirectional models when describing the magnetic behavior in weakly magnetoelastic materials (e.g., Nickel), but the two models deviate when highly magnetoelastic materials (e.g., Terfenol-D) are introduced. These results suggest that magnetostrictive feedback is critical for modeling highly magnetoelastic materials, as opposed to weaker magnetoelastic materials, where we observe only minor differences between the two methods’ outputs. To our best knowledge, this work represents the first comparison of unidirectional and bidirectional modeling in composite multiferroic systems, demonstrating that back-coupling of magnetization to strain can inhibit formation and rotation of magnetic states, highlighting the need to revisit the assumption that unidirectional modeling always captures the necessary physics in strain-mediated multiferroics.


Note S1 Initial magnetization classifications
For nanoscale magnetic structures, it is mainly the competition between magnetic anisotropies and the geometry of the structures that determines the initialization state at the equilibrium, which relaxes from the saturated initial state, where all individual magnetic moments throughout the ring point towards +x-direction ( Figure 1b). 1 During relaxation from the saturated initial state, the competition between exchange energy and magnetostatic energy (i.e., demagnetization energy) plays a pivotal role for total energy minimization in the magnetic structure. The micromagnetic simulation model can be simplified for characterizing the dependence of the initial mapping of domain states on the ring dimensions. For obtaining initial mapping of domain states, we only consider the dominant energy terms, i.e., exchange energy and magnetostatic energy, from the Landau-Lifshitz-Gilbert (LLG) equation. 2 In this way, such initial remanent states for various ring geometries can be approximated using micromagnetics alone.
To find the magnetic ground state in the magnetoelastic rings, we use OOMMF micromagnetic simulator to identify domain states attainable in rings with outer diameter (OD) on micron scale, width (w) and thickness (t) on submicron scale. According to previous investigation on Ni rings, at relaxation states after initialization, these submicron size scales rings are capable of achieving an "onion state" with either transverse DWs or vortex DWs on opposite sides of the ring, and a vortex state with flux closure domains around the ring. 3 DW states at equilibrium in Terfenol-D rings with different dimensions have been modeled using the OOMMF micromagnetic simulation with an initial magnetization all aligned in +x direction. 4 In this work, we simulated nanoscale rings with a variety of dimensions, including rings with OD of 1 µ m, w of 50 nm, 150 nm, 300 nm, and 400 nm, and t of 15 nm, 30 nm and 45 nm, along with rings with OD of 2 µ m, w of 150 nm, 300 nm, 400 nm, 600 nm, and t of 15 nm, 30 nm and 45 nm. Correlation between energy density and magnetic configurations in those ring structures is further examined to produce phase diagrams with design specifications for rings with geometries that are energetically favorable for the desired onion states. 2,5 The relevant length scales for the ring structure are the OD, w, t and the magnetostatic exchange length ! "# expressed as ! "# = %& ' ( ) * + , where A is the exchange stiffness, and ! " is the saturation magnetization. 5 For Terfenol-D, the calculation is performed based on the following parameters, A = 1.0 × 10 -11 J/m and ! " = 8.0 × 10 5 A/m, giving ! "# = 5.0 nm. 6 For the micromagnetic model, the maximum dm/dt is set as 0.02 for a more precise simulation of domain states in Terfenol-D ring structures with OD of 1 µ m and 2 µ m. Due to the competition between demagnetization and exchange energies, the magnetic ring states vary from onion states to vortex state with a flux closure configuration when the ring's width is increased.
Onion states with transverse DWs are of interest due to their large energy density and stray field.
The energy flux out of the ring can have practical applications in trapping and/or interacting with nano-and micro-scale particles in the surroundings of the ring via localized magnetostatic interaction. 7 Depending on the ring dimension and w/OD ratio, the magnetic energy density of rings with the same OD versus thickness of the ring are plotted in Figure S1a

Note S2 Equivalent coupled model setup
To reduce simulation time while maintaining accurate results, the multiphysics finite element model can also be set up with a Terfenol-D or Ni ring on top of a SiO 2 substrate, as shown in Figure S2. This separates the computation required for piezoelectric strain to a separate modeling step. Tensile strain is then induced along 45 o in the counter clockwise direction away from the +x axis, and compressive strain along -45 o to the +x axis. Therefore, for Terfenol-D with positive magnetostriction, magnetization and DWs tend to orient along 45 o away from the +x axis (i.e., along tensile strain direction), same DW behavior as on PMN-PT substrate.
The purpose of choosing a thickness of 500 nm for the substrate rather than a few hundred of microns to match PMN-PT substrate thickness is to make the computations feasible while not affecting the results.

Figure S2
Geometry of the setup in COMSOL Multiphysics, where strain is applied to SiO 2 substrate to achieve similar effect as applying voltage to the piezoelectric substrate.

Note S3 Influence of ramping speed on the domain wall dynamics
The speed of ramping electric field/ strain can affect the domain wall rotation dynamics. Strain of different amplitudes are applied to the substrate with the following ramping slopes k = 10 10 s -1 , 10 9 s -1 , and 10 8 s -1 , and the corresponding time for ramping t is 0.1 ns, 1 ns and 10 ns. Here, we show the average magnetization angle rotation dynamics predicted by bidirectional model for Terfenol-D ring with an applied strain of 500 ppm (as shown in Figure S3). By comparison, we conclude that the ramping slope of the electric field/strain affects the dynamics of the system, but can lead to similar magnetization states in the ring at equilibrium following the application of strain for this case.

Figure S3
Domain wall angle as a function of strain application time in Terfenol-D rings calculated by bidirectional models, with a strain of 500 ppm generated in the substrate.

Note S4 Influence of mesh size and time step on convergence
To compare the effect of mesh size and time step on convergence, we used the BD model to predict magnetization variation in the Terfenol-D ring when a 1000 ppm strain is applied. The model setup that produces the result in the main text adopts a mesh size of 10 nm for the ring, and 40 nm for the substrate. A model with much finer mesh element setup has a mesh size of 5 nm (close to the exchange length in Terfenol-D) for the ring, and 10 nm for the substrate. The time step is taken as 5 ps for the first model, and 1ps for the second model. As shown in Figure   S4, magnetization rotation predicted by the BD model with larger mesh size (S4, left) predicts the same tendency as modeled by the BD model with smaller mesh size (S4, right). With a dramatically increased number of mesh elements, the latter model obtains better convergence during simulation (convergence error below 10 -$ ), while still following the same domain rotation trend as predicted by the one adopted in the letter. In addition, Figure S5 shows the magnetization configurations in Terfenol-D ring predicted by both models at equilibrium. We thus conclude that the finest mesh size used here is not necessary as it is much more timeconsuming and differs little from the ones using larger mesh size. The model used for the letter is sufficient to describe the overall magnetization dynamics in the systems.

List of Supporting Videos
Supporting Video S1 Domain wall rotation in Terfenol-D ring predicted by a unidirectional model, with an applied strain of 500 ppm (avi format, 1400 frames, 1280*1280px) Supporting Video S2 Domain wall rotation in Terfenol-D ring predicted by a bidirectional model, with an applied strain of 500 ppm (avi format, 1400 frames, 1280*1280px)