doc changes #252

doc/content/MSCA_EU_Horizon2020_Results/horizon2020_ex1.md

@@ -119,9 +119,9 @@ and `AuxVariables` blocks,
          link=False
          language=python
 
-where we use (i.e. for vector component $P_z$) `initial_from_file_var = polar_z` and `initial_from_file_timestep = 'LATEST'` as code flags to load in the values of `polar_z` and others at the latest timestep from the Exodus output `BFO_P0A0.e`. In this sense, we have changed $\{\mathbf{P},\mathbf{A}\}$ to be an AuxVariable class instead of a Variable since they will be fixed for the duration of the evolution of the spins. Our variables to solve for are the set of sublattice magnetizations $\{\mathbf{m}_1, \mathbf{m}_2\}$ with postprocessed outputs $\mathbf{L} = \mathbf{m}_1 - \mathbf{m}_2$ and $\mathbf{m} = \mathbf{m}_1 + \mathbf{m}_2$. Note that the conventional factor of 2 is missing from the output of these calculations.
+where we use (i.e. for vector component $P_z$) `initial_from_file_var = polar_z` and `initial_from_file_timestep = 'LATEST'` as code flags to load in the values of `polar_z` and others at the latest timestep from the Exodus output `BFO_P0A0.e`. In this sense, we have changed $\{\mathbf{P},\mathbf{A}\}$ to be an AuxVariable class instead of a Variable since they will be fixed for the duration of the evolution of the spins. Our variables to solve for are the set of sublattice magnetizations $\{\mathbf{m}_1, \mathbf{m}_2\}$ with postprocessed outputs of the Néel vector $\mathbf{L} = \mathbf{m}_1 - \mathbf{m}_2$ and net (weak) magnetization $\mathbf{m} = \mathbf{m}_1 + \mathbf{m}_2$. Note that the conventional factor of 2 is missing from the output of these calculations.
 
-This problem considers the evolution of the Landau-Lifshitz-Gilbert (Bloch) equation for the two sublattices,
+This problem considers the evolution of the normalized Landau-Lifshitz-Gilbert (Bloch) equation for the two sublattices,
 
 \begin{equation}\label{eqn:LLG_LLB}
   \begin{aligned}

doc/content/MSCA_EU_Horizon2020_Results/horizon2020_ex4.md

@@ -4,6 +4,52 @@
 
 This page details how to obtain ME switching trajectories in FERRET. We use an electric field to switch $\mathbf{P}$ explicitly in a fully-time dependent simulation. The magnetization follows $\mathbf{P}$ since any change in $\mathbf{P}$ (and $\mathbf{A}$) corresponds to a drastic change in the magnetic energy, $f_\mathrm{MP}$, that couples to the structural distortions.
 
+Our full approach involves solving the couple dynamic equation system,
+
+\begin{equation}\label{eqn:TDLG}
+  \begin{aligned}
+    \frac{\partial \mathbf{P} }{\partial t} &= -\Gamma_P \frac{\delta f}{\delta \mathbf{P}}, \\
+    \frac{\partial \mathbf{A} }{\partial t} &= -\Gamma_A \frac{\delta f}{\delta \mathbf{A}},
+  \end{aligned}
+\end{equation}
+
+for the structural order ($\{\mathbf{P},\mathbf{A}\}$) along with,
+
+\begin{equation}
+  \begin{aligned}
+    \frac{\partial \sigma_{ij}}{\partial x_j} = \frac{\partial}{\partial x_j}\left[C_{ijkl} \left(\varepsilon_{kl} + Q_{klmn} P_m P_n + R_{klmn} A_m A_n \right)\right] = 0,
+  \end{aligned}
+\end{equation}
+
+at every time step, where $Q_{klmn}$ and $R_{klmn}$ are the electro- and rotostrictive coefficients that couple $\mathbf{P}$ and $\mathbf{A}$ to the elastic strain respectively. For the spins, we solve,
+
+\begin{equation}\label{eqn:LLG}
+  \begin{aligned}
+    \frac{\partial \mathbf{M}_\eta}{\partial t} = -\left(\frac{\gamma}{1+\alpha^2}\right)\mathbf{M}_\eta\times \mathbf{H}_\eta - \frac{1}{M_s} \left(\frac{\gamma \alpha}{1+\alpha^2}\right) \mathbf{M}_\eta \times \left(\mathbf{M}_\eta \times \mathbf{H}_\eta\right).
+  \end{aligned}
+\end{equation}
+
+Here, $\Gamma_P, \Gamma_A$ are a relaxation coefficients related to the time scales involved in the structural phase transition. The parameter $\gamma$ is the electron gyromagnetic ratio and $\mathbf{H}_\eta$ is the effective field acting on sublattice $\eta$. The coeffiicent $\alpha$ is a phenomenological damping constant which if made nonzero (and positive) drives the magnetic system to the ground state. We select a time-dependent electric field $\mathbf{E}$ of the form,
+
+\begin{equation}\label{eqn:E}
+  \begin{aligned}
+    \mathbf{E}(\omega) = \langle 0,0, E_0 \rangle \sin{\left(\omega t\right)}.
+  \end{aligned}
+\end{equation}
+
+which facilitates a reversal of the $P_z$ component. This means we expect the transition of $\mathbf{P}$ to go from $[111]$ to $[11\bar{1}]$ orientation. For this problem, we use the input file `BFO_P111_TO_P111b_switch_m1_a1.i` located in the tutorials subdirectory. The `Exodus` input that we use corresponds to that of the second simulation in Example 1 (a fully relaxed polar-magnetic solution). We load this via the `Mesh` block,
+
+INSERT
+
+The `Kernels` are both the ones due to structural evolution (TDLGD) along with the ones from micromagnetic evolution (LLG-LLB). They are listed in the following lengthy block,
+
+INSERT
+
+A possible visualization of the output using ParaView is provided below,
+
+!media media/tut_ME_sw.png style=display:block;margin:auto;width:67%; caption=Homogeneous switching at $\alpha = 0.003$ with Top: Neel vector switching and Bottom: net magnetization switching. These dynamics are acquired using a time-dependent electric field at frequency $\omega = 600$ MHz.  id=fig_ME_sw
+
+
 -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 
 This project [SCALES - 897614](https://cordis.europa.eu/project/id/897614) was funded for 2021-2023 at the [Luxembourg Institute of Science and Technology](https://www.list.lu/) under principle investigator [Jorge Íñiguez](https://sites.google.com/site/jorgeiniguezresearch/). The research was carried out within the framework of the [Marie Skłodowska-Curie Action (H2020-MSCA-IF-2019)](https://ec.europa.eu/info/funding-tenders/opportunities/portal/screen/opportunities/topic-details/msca-if-2020) fellowship.

doc/content/MSCA_EU_Horizon2020_Results/horizon2020_results1.md

@@ -48,7 +48,7 @@ coupled to the direction of the antiphase tilts $\mathbf{A}$ which acts to cant
 
 A schematic of the spin order is shown in Fig. 2 (left) easy-plane anisotropy of $\mathbf{m}_\eta$ relative to $\mathbf{P}$ direction and (right) DMI-induced canting in the easy-plane (shown as cyan circles). The quantities $\mathbf{L}$ and $\mathbf{m}$ are defined such that $|\mathbf{L}| + |\mathbf{m}| = 1$ with, in general, $|\mathbf{L}| \gg |\mathbf{m}|$ reflecting the presence of a strong AFM coupling between the sublattices but with a weak noncollinearity in $\mathbf{m}_1$ and $\mathbf{m}_2$. The total weak magnetization is $\mathbf{M} = M_s \mathbf{m}$ where $M_s$ is the saturation magnetization density of the Fe sublattice ($\approx 4.0 \mu$ B/Fe) [!cite](Dixit2015).
 
-We solve the couple dynamic equation system,
+Our full approach involves solving the couple dynamic equation system,
 
 \begin{equation}\label{eqn:TDLG}
   \begin{aligned}
@@ -57,17 +57,25 @@ We solve the couple dynamic equation system,
   \end{aligned}
 \end{equation}
 
-for the structural order ($\{\mathbf{P},\mathbf{A}\}$) and
+for the structural order ($\{\mathbf{P},\mathbf{A}\}$) along with,
+
+\begin{equation}
+  \begin{aligned}
+    \frac{\partial \sigma_{ij}}{\partial x_j} = \frac{\partial}{\partial x_j}\left[C_{ijkl} \left(\varepsilon_{kl} + Q_{klmn} P_m P_n + R_{klmn} A_m A_n \right)\right] = 0,
+  \end{aligned}
+\end{equation}
+
+at every time step, where $Q_{klmn}$ and $R_{klmn}$ are the electro- and rotostrictive coefficients that couple $\mathbf{P}$ and $\mathbf{A}$ to the elastic strain respectively. For the spins, we solve,
 
 \begin{equation}\label{eqn:LLG}
   \begin{aligned}
-    \frac{\partial \mathbf{M}_\eta}{\partial t} = -\gamma \mathbf{M}_\eta\times \mathbf{H} - \frac{1}{M_s} \frac{\gamma \alpha}{1+\alpha^2} \mathbf{M}_\eta \times \mathbf{M}_\eta \times \mathbf{H},
+    \frac{\partial \mathbf{M}_\eta}{\partial t} = -\left(\frac{\gamma}{1+\alpha^2}\right)\mathbf{M}_\eta\times \mathbf{H}_\eta - \frac{1}{M_s} \left(\frac{\gamma \alpha}{1+\alpha^2}\right) \mathbf{M}_\eta \times \left(\mathbf{M}_\eta \times \mathbf{H}_\eta\right).
   \end{aligned}
 \end{equation}
 
-for the spin system. Here, $\Gamma_P, \Gamma_A$ are a relaxation coefficients related to the time scales involved in the structural phase transition. The coeffiicent $\alpha$ is a phenomenological damping constant which if made nonzero (and positive) drives the magnetic system to the ground state.
+Here, $\Gamma_P, \Gamma_A$ are a relaxation coefficients related to the time scales involved in the structural phase transition. The parameter $\gamma$ is the electron gyromagnetic ratio and $\mathbf{H}_\eta$ is the effective field acting on sublattice $\eta$. The coeffiicent $\alpha$ is a phenomenological damping constant which if made nonzero (and positive) drives the magnetic system to the ground state.
 
-A detailed description of our model is shared in the preprint on arXiv at [!cite](Mangeri2023).
+A detailed description of our model is shared in the preprint on arXiv at [!cite](Mangeri2023) and in the following sections we step through our key results and provide representative example files and documentation to reproduce the calculations.
 
 -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 

doc/content/MSCA_EU_Horizon2020_Results/horizon2020_results2.md

@@ -44,27 +44,27 @@ where $Q_{ijkl}$ and $R_{ijkl}$ are the electrostrictive and rotostrictive coeff
 
 Note that $(\mathbf{P} \uparrow \downarrow \mathbf{A})$ is also a possible minimized energy solution of the thermodynamic potential $f$. Due to the symmetry of the electrostrictive and rotostrictive coupling terms, the table is left invariant under full reversal of $\mathbf{A}$. The free energy density of the eight-fold domain possibilities is -15.5653 $\mathrm{eV}\cdot\mathrm{nm}^{-3}$.
 
-To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md) for more details on the input file.
+To reproduce these results, we refer the reader [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md) for more details on the input file and the documentation to reproduce these results.
 
 # Antiferromagnetic ringdown (homogeneous spin state)
 
-To find the possible magnetic states in this material, we evolve the Landau-Lifshitz-Bloch (LLB) equation,
+To find the possible magnetic states in this material, we evolve the normalized two-sublattice Landau-Lifshitz-Bloch (LLB) equation,
 
 \begin{equation}\label{eqn:LLG_LLB}
   \begin{aligned}
     \frac{d \mathbf{m}_\eta}{dt} = -\frac{\gamma}{1+\alpha^2} \left(\mathbf{m}_\eta \times \mathbf{H}_\eta\right) - \frac{\gamma\alpha}{1+\alpha^2}\mathbf{m}_\eta\times\left(\mathbf{m}_\eta\times\mathbf{H}_\eta\right) + \frac{\gamma\tilde\alpha_\parallel}{(1+\alpha^2)} m_{\eta}^2 \left[ m_{\eta}^2 - 1 \right]\mathbf{m}_\eta,
   \end{aligned}
 \end{equation}
 
-where $\mathbf{P}$ and $\mathbf{A}$ are held fixed corresponding to one of the possible eight-fold degenerate domain directions. This approximation for evaluation of ground states is suitable since the magnetic system does not appreciably influence the relaxation dynamics of the polarization since the work required to move ions (i.e. those relative displacements corresponding to $\mathbf{P}$ and $\mathbf{A}$) is much larger than those to reorient the spins on the Fe sites. The constant $\gamma$ corresponds to the electron gyromagnetic factor equal to $2.2101\times10^5$ rad. m A${}^{-1}$ s ${}^{-1}$, $\alpha$ the Gilbert damping, and $\tilde\alpha_\parallel$ the longitudinal damping constant from the LLB approximation. The effective fields are defined as $\mathbf{H}_\eta = - \mu_0^{-1} M_s^{-1} \delta f / \delta \mathbf{m}_\eta$ with $\mu_0$ the permeability of vacuum. The saturation magnetization density of the BFO sublattices is $M_s = 4.0$ $\mu$B/Fe [!cite](Dixit2015). We set $\alpha = 0.01$ and $\tilde\alpha_\parallel = 10^2$ and ringdown the magnetic system to an energy minimum.
+where $\mathbf{P}$ and $\mathbf{A}$ are held fixed corresponding to one of the possible eight-fold degenerate domain directions. This approximation for evaluation of ground states is suitable since the magnetic system does not appreciably influence the relaxation dynamics of the polarization since the work required to move ions (i.e. those relative displacements corresponding to $\mathbf{P}$ and $\mathbf{A}$) is much larger than those to reorient the spins on the Fe sites. The constant $\gamma$ corresponds to the electron gyromagnetic factor equal to $2.2101\times10^5$ rad. m A${}^{-1}$ s ${}^{-1}$, $\alpha$ the Gilbert damping, and $\tilde\alpha_\parallel$ the longitudinal damping constant from the LLB approximation. The effective fields are defined as $\mathbf{H}_\eta = - \mu_0^{-1} M_s^{-1} \delta f / \delta \mathbf{m}_\eta$ with $\mu_0$ the permeability of vacuum. The saturation magnetization density of the BFO sublattices ($\eta = 1,2$) is $M_s = 4.0$ $\mu$B/Fe [!cite](Dixit2015). We set $\alpha = 0.01$ and $\tilde\alpha_\parallel = 10^2$ and ringdown the magnetic system to an energy minimum.
 
-The below figure details the ringdown behavior of the numerical sample. Here, we also look for homogeneous spins (no domain walls or texture in the $\mathbf{m}_\eta$ field).
+The below figure details the numerical ringdown behavior of the sublattices $\mathbf{m}_\eta$. Here, we also look for homogeneous spins (no domain walls or texture in the $\mathbf{m}_\eta$ field).
 
 !media media/ringdown_abc.png style=display:block;margin:auto;width:100%; caption=Ringdown time dependence of Left: components of selected $\mathbf{m}_\eta$. Center: easy-plane angle $\theta_\eta = \cos^{-1}{(\mathbf{m}_\eta\cdot\hat{\mathbf{P}})}$ and Right: canted moment angle $\phi^\mathrm{WFM} = \cos^{-1}{(\mathbf{m}_1\cdot \mathbf{m}_2)}$.  id=fig_ringdown_abc
 
 Importantly, this analysis allowed us to parameterize the DMI coupling at our level of theory, which leads to the non-vanishing moment of $\mathbf{m} = (\mathbf{m}_1 + \mathbf{m}_2)/2$. By tracking the angle of the canted moment $\phi^\mathrm{WFM} = \cos^{-1}{\left(\mathbf{m}_1 \cdot \mathbf{m}_2 \right)}$ as a function of the strength of the DMI coupling, $D_0$, we can identify good agreement ($\phi^\mathrm{WFM} \approx 1.22^\circ$) with the available literature corresponding to $\mathbf{M} = M_s\mathbf{m} = 0.03 \mu$B/f.u..
 
-By selecting initial conditions corresponding to the six-fold possible orientations of the spin sublattice system, the below table can be obtained yielding 48 different orientations of the the Neel order $\mathbf{L} = \left(\mathbf{m}_1 - \mathbf{m}_2\right) / 2$ and total magnetization $\mathbf{m}$ dependent on the 8 possible orientations of $(\mathbf{P}\uparrow\uparrow\mathbf{A})$.
+By selecting initial conditions corresponding to the six-fold possible orientations of the spin sublattice system, the below table can be obtained yielding 48 different orientations of the the Néel order $\mathbf{L} = \left(\mathbf{m}_1 - \mathbf{m}_2\right) / 2$ and total magnetization $\mathbf{m}$ dependent on the 8 possible orientations of $(\mathbf{P}\uparrow\uparrow\mathbf{A})$.
 
 !equation
 \begin{array}{c|cccccccc}
@@ -90,7 +90,7 @@ By selecting initial conditions corresponding to the six-fold possible orientati
  & [\bar{2}11] & [211] & [\bar{2}\bar{1}\bar{1}] & [2\bar{1}\bar{1}] & [121] & [\bar{1}21] & [1\bar{2}\bar{1}] & [\bar{1}\bar{2}\bar{1}] \\
 \end{array}
 
-To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md) for details on the input file.
+To reproduce these results, we refer the reader to our examples documentation [here](MSCA_EU_Horizon2020_Results/horizon2020_ex1.md) for details on the simulation (input/output).
 
 # Structural domain walls
 
@@ -110,11 +110,11 @@ and
   \end{aligned}
 \end{equation}
 
-A comma in the subscript denotes a partial derivative with respect to the specified spatial directions. In order to study the domain wall topology involving spatial variations of $\mathbf{P}$, $\mathbf{A}$, and strain, a good parameter set estimate of the gradient coefficients $(G_{11}, H_{11}, ...)$ is needed. To achieve this, we consult DFT calculations reported in [!cite](Dieguez2013). It was shown that an assortment of different metastable (along with low energy) domain walls form an energetic hierarchy. This was advantageous for our work, as we were allowed to seperate the computation (fitting) of various $(G_{11}, H_{11}, ...)$ due to only certain terms in the gradient energy density expressions being identified as primary contributions to the DW energy.
+A comma in the subscript denotes a partial derivative with respect to a specified spatial direction. In order to study the domain wall topology involving spatial variations of $\mathbf{P}$, $\mathbf{A}$, and strain, a good parameter set estimate of the gradient coefficients $(G_{11}, H_{11}, ...)$ is needed. To achieve this, we consult DFT calculations reported in [!cite](Dieguez2013). It was shown that an assortment of different metastable (along with low energy) domain walls form an energetic hierarchy. This was advantageous for our work, as we were allowed to seperate the computation (fitting) of various $(G_{11}, H_{11}, ...)$ due to only certain terms in the gradient energy density expressions being identified as primary contributions to the DW energy.
 
 For example, consider the 2/1 (100) DW which is a commonly observed domain boundary observed in experiments. In this notation, it is indicated that, for the $2/1$ DW, two components of $\mathbf{P}$ and one component of $\mathbf{A}$ vary across the boundary whose plane normal is (100) whereis for the $3/0$ DW, $\mathbf{P}$ undergoes a full reversal where $\mathbf{A}$ is approximately unchanged across the (110)-oriented boundary plane. We label the pairs of the domains characterizing the DW as $\mathbf{P}^\mathrm{I}/\mathbf{A}^\mathrm{I}$ and $\mathbf{P}^\mathrm{II}/\mathbf{A}^\mathrm{II}$ in the Table below.
 
-After seeding a sin(x) profile for the initial state of $\mathbf{P}$ and $\mathbf{A}$ corresponding to the expected DW configuration, a typical profile of the order parameters and strain can be obtained as shown below for the $2/1 [100]$ DW when relaxing the time-dependent Landau-Ginzburg equations - Eqs. (\ref{eqn:TDLGD}).
+After seeding a $sin(x)$ profile for the initial state of $\mathbf{P}$ and $\mathbf{A}$ corresponding to the expected DW configuration, a typical profile of the order parameters and strain can be obtained as shown below for the $2/1 [100]$ DW when relaxing the time-dependent Landau-Ginzburg equations - Eqs. (\ref{eqn:TDLGD}).
 
 !media media/tut_21_DW.png style=display:block;margin:auto;width:47%; caption=Different components of $\mathbf{P}$ and $\mathbf{A}$ across the $2/1 [100]$ domain wall.  id=fig_dw_21_prof