A coupled magneto-structural continuum model for multiferroic $\mathrm{BiFeO}_{3}$

We define the free energy density corresponding to the structural distortions of the lattice as $f_{\mathrm{latt}}$,

The energy expansion of $f_{P},f_{A}$ and $f_{AP}$ contains only the terms allowed by symmetry to the fourth orderFedorova et al. (2022),

and

Additionally, the elastic, electrostrictive ($\mathbf{P}$-$\bm{\varepsilon}$), and rotostrictive ($\mathbf{A}$-$\bm{\varepsilon}$) energy is included as

and

respectively. Finally, to evaluate inhomogeneous phases (i.e DWs), we include the lowest-order Lifshitz invariants Cao and Barsch (1990); Li et al. (2001); Hlinka and Marton (2006) for the structural distortions to Eq. (3),

and

for both the $\mathbf{P}$ and $\mathbf{A}$ order parameters respectively. A comma in the subscript denotes a partial derivative with respect to the specified spatial directions. The bulk homogeneous contribution to the energy (i.e. the terms not involving $f_{\nabla P}$ and $f_{\nabla A}$) has been previously parameterized with DFT calculationsFedorova et al. (2022); we refer the reader to this publication for the relevant coefficients.

However, in the case of the gradient energy, the set of coefficients $\{G_{ij}$, $H_{ij}\}$ are difficult to obtain directly from DFT (see for example the approach outlined in Refs. [(48; 49; 50)]) - so we employ a fitting procedure in Sec. II.4 to evaluate them. We should emphasize that if a different bulk homogeneous phenomenological potential is used (i.e. Refs. [(51; 52; 36)]), then the gradient coefficients obtained would be different since they depend strongly on the energetics of the order parameters in the vicinity of the DW.

To find the polar ground states, we evolve the coupled time dependent Landau-Ginzburg (TDLG) equations,

and

along with satisfying the stress-divergence equation for mechanical equilibrium,

where $\sigma_{ij}=\sigma_{ji}=\partial f_{\mathrm{latt}}/\partial\varepsilon_{ij}$ is the elastic stress of the material. We write the components of $\sigma_{ij}$ as

where $\varepsilon_{kl}$ is the elastic strain from Eq. (2) and the eigenstrain is related to the spontaneous strain via,

where $Q_{ijkl}$ and $R_{ijkl}$ are the electrostrictive and rotostrictive coefficients. These are related to our free energy density coefficients $q_{ijkl}$ and $r_{ijkl}$ as (in Voight notation),

and

with similar definitions for the quantities involving $R_{ijkl}$. We also investigate electrostatic phenomena in our model through the Poisson equation,

where $\Phi_{\mathrm{E}}$ is the electrostatic potential which defines the electric field $\mathbf{E}=-\nabla\Phi_{\mathrm{E}}$ in the usual way. The parameter $\epsilon_{b}=30\,\epsilon_{0}$ is the relative background dielectric constant Graf et al. (2015). Eq. (18) is solved at every time step of the evolution of Eq. (10) and (11). In Sec. II.4 we are searching for the local minima due to the relaxation dynamics of Eq. (10) and (11) and as such the time relaxation constants $\Gamma_{P}$ and $\Gamma_{A}$ are set to unity.

To enforce periodicity on the strain tensor components in our representative volume element that includes DWs, we separate the strain fields calculated from Eq. (2) and (12) into homogeneous (global) and inhomogeneous (local) parts. This is done utilizing the method formulated by Biswas and co-workers in Ref. [(54)] which relaxes the stress components along the periodic directions and thus allows corresponding deformation to occur. Here, the homogeneous contribution of the total strain obeys the following integrated quantity at every time step of the relaxation,

where $V$ is the volume of our simulation containing the DW profiles. The total stress tensor, $\sigma_{ij}^{\mathrm{total}}$, is calculated from the sum of homogeneous, inhomogeneous, and eigenstrain components $\varepsilon_{ij}^{\mathrm{total}}=\varepsilon_{ij}^{\mathrm{inhom}}+\varepsilon_{ij}^{\mathrm{hom}}+\varepsilon_{ij}^{\mathrm{eig}}$ for all periodic directions $i$ and corresponding periodic component $j$ at every time step of the simulation.

Equations (10), (11), (12), (18), and (19) are cast into their weak formulation sufficient for the finite element analysis. Our method uses linear Lagrange shape functions for the coupled variable system. The finite element mesh spacing is selected to be $\Delta x\approx 0.1$ nm for all calculations in this work. This small mesh spacing helps resolve the thin DWs in BFO to smoothness which are discussed extensively in Section II.4 and II.7. We implement Newmark-beta time integration Newmark (1959) with convergence between time steps achieved when the nonlinear residuals calculated during the Newton-Raphson iteration (with block Jacobi preconditioning) have been reduced by $10^{-8}$ relative tolerance. If convergence is not obtained, we use adaptive time stepping with reduction factor of $0.5$. The finite element method (FEM) implementation of this work is available within Ferret Mangeri et al. (2017) which is an add-on module for the open source Multiphysics Object Oriented Simulation Environment (MOOSE) framework Permann et al. (2020).

In the absence of order parameter gradients, the homogeneous FE states of $\mathbf{P}$ parallel to $\mathbf{A}$ which we denote as $\mathbf{P}\uparrow\uparrow\mathbf{A}$ can be obtained numerically. To perform this calculation, we evolve Eq. (10) and (11) simultaneously solving Eq. (12) (at every time step) until the relative change in total volume integrated energy density $F$ between adjacent time steps is less than $5\times 10^{-7}$ eV/s. The bulk potential predicts the spontaneous values of the order parameters upon minimization that are $P_{s}=|\mathbf{P}|=0.945\,\,\mathrm{C}/\mathrm{m}^{2}$ and $A_{s}=|\mathbf{A}|=13.398^{\circ}$. The spontaneous normal and shear strains that correspond to these values are $\varepsilon_{n}=\varepsilon_{ii}=1.308\times 10^{-2}$ and $\varepsilon_{s}=\varepsilon_{ij}=2.95\times 10^{-3}$ for $i\neq j$ in agreement with Ref. [(44)]. The free energy density of the ground state given by Eq. (3) is -15.5653 $\mathrm{eV}\cdot\mathrm{nm}^{-3}$. The energy functional used also describes identical energy minima when $\mathbf{P}\uparrow\downarrow\mathbf{A}$ (which is equivalent to a $180^{\circ}$ phase reversal of the tilt field). Since the rotostrictive strains defined in Eq. (14) are invariant upon full reversal of $\mathbf{A}$, then these numbers are left unchanged. In the next Section II.4 we evaluate the inhomogeneous textures of the DWs and parameterize the gradient coefficients $\{G_{ij},H_{ij}\}$ used in our model.

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},...)$ of Eq. (8) and Eq. (9) is needed. To achieve this, we consult DFT calculations reported by Diéguez and co-workers in Ref. [(40)]. It was shown that an assortment of metastable states are allowed in BFO and that this zoology of different DW types forms an energy hierarchy. Due to electrostatic compatibility, this collection of states has specific requirements on the components of the order parameters that modulate across the domain boundary. For example, the lowest energy configurations which we denote (see Table 1) as $2/1(100)$ and $3/0(110)$ are the $109^{\circ}$ and $180^{\circ}$ DWs respectively. In this notation, it is indicated that, for the $2/1$ DW, two components of $\mathbf{P}$ and one component of $\mathbf{A}$ switch sign across the boundary whose plane normal is (100), whereas for the $3/0$ DW, $\mathbf{P}$ undergoes a full reversal where $\mathbf{A}$ is 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 this table. This determines which terms in Eq. (8) and (9) are primary contributions to the DW energy. This is particularly advantageous as it has allowed us to separate the computation of specific DWs in the analysis of fitting the gradient coefficients to the DFT results.

To obtain the (100)- or (110)-oriented DWs within our phase field scheme, we choose an initial condition for the components of the order parameters to be a sin(x) or sin(x+y) profile respectively. We then relax Eq. (10), and (11) until convergence along with satisfying the conditions of mechanical equilibrium of Eq. (12) at every time step. The periodic boundary conditions on the components of $\mathbf{P},\mathbf{A},$ and $\mathbf{u}$ for (100)- or (110)-oriented domain walls are enforced along the [100] and [110] directions respectively. We compute the DW energy with

where $F_{0}$ is the corresponding monodomain energy from Eq. (3) integrated over the computational volume $V$. The energy $F$ is computed from the solution that contains the DW profile. The number of DWs in the simulation box is $N$ and $S$ the surface area of the DW plane. We find convergence on the computed energies within 1 $\mathrm{mJ}/\mathrm{m}^{2}$ provided that the DW-DW distances are greater than 30 nm due to long-range strain interactions. For fourth-order thermodynamic potentials, a fit function of the form $W_{k}\tanh{\left[\left(r-r_{0}\right)/t_{k}\right]}$ is sufficient to fit the evolution of order parameters that switch across the DWMarton et al. (2010) where $W_{k}$ is the value of the switched spontaneous order parameters far from a DW plane localized at $r_{0}$ and $t_{k}$ corresponds to the thickness of the polar or octahedral tilt parameters for $k=P,A$ respectively.

As a first example, consider the lowest energy DW predicted by DFT, the so-called $109^{\circ}$ 2/1 (100) DW which is indeed frequently observed in thin film samples of BFO Zhang et al. (2018); Parsonet et al. (2022). The primary gradient coefficients governing the energy of the wall are the $H_{11}$ and $G_{44}$ coefficients owing to the fact that $A_{x,x}$, $P_{y,x}$, and $P_{z,x}$ are nonzero (see Table 1). The resulting DW profile for the 2/1 (100) wall is presented in Fig. 1(a). The profile is a smooth rotation of both $A_{x}$ and $P_{y}=P_{z}$ across the wall region. The inset on the left reveals that the non-switching component $P_{x}$ experiences a slight decrease ($\approx-3\%$) at the wall. The quantitative value of the modulation of the non-switched component is consistent with DFT results of the same DW typeKörbel and Sanvito (2020).

The small change of $P_{x}$ corresponds with a built-in $\Phi_{\mathrm{E}}$ shown in the right inset panel which is of comparable order ($\approx 10$ mV) to those estimated from DFTKörbel and Sanvito (2020). Fitting the $\mathbf{P}$-$\mathbf{A}$ profile shows that the DW is quite thin (thickness $2t_{P}\approx 0.08$ nm). Hence, we obtain DWs with marked Ising character. We provide an energy profile scan across the primary coefficients $H_{11}$ and $G_{44}$ in (b). The dashed white line outlines the predictions from DFT results in Ref. [(40)]. We also should mention that the dependence on other coefficients is quite weak due to the relatively small gradients in non-switching components. These calculations (and others not shown here) reveal that the choice of $\{G_{ij},H_{ij}\}$ is not unique, i.e. one can find the same DW energies (with very similar profiles) for different combinations of the primary coefficients. Therefore, it is necessary to visit other DW configurations to constrain the values of the entire set.

Next, we present $\mathbf{P}$-$\mathbf{A}$ profiles of three higher energy (100)-oriented domain walls (1/1, 0/3, and 2/2) in Fig 2 (a), (b) and (c) respectively. These three calculations correspond to those using our best estimates of the gradient coefficients $\{G_{ij},H_{ij}\}$ in Table 2. In all three cases, we find the presence of a small changes in the non-switching components of the order parameters shown in circles for $A_{k}$ and diamonds for $P_{k}$. For example, in the $71^{\circ}$ 1/1 shown in DW Fig 2 (a), $P_{y}$ (in red) which does not change sign, grows at the DW by about 15$\%$. This is in contrast to the $P_{x}$ component (in blue) which only grows by 2.5$\%$ demonstrating the influence of the weak built-in field which reduces the magnitude of this component to keep this wall neutral. Similar changes on the order of about 10$\%$ are also seen in $A_{x}=A_{y}$ components shown in blue. This DW-induced change in P seems to be the largest in the 0/3-type DW shown in (b). Due to the influence of built in electric fields from the solution of the Poisson equation (and our best estimates of the anisotropic gradient coefficients), the value of $P_{x}$ component grows by about 5$\%$ whereas the $P_{y}=P_{z}$ components diminish by almost -35$\%$ (shown in black). Again, we also find changes in the non-switching components in the $109^{\circ}$ 2/2 wall, with $P_{x}$ (blue diamonds) growing by about 2$\%$; by contrast and $A_{x}$ decreases by $-6.4^{\circ}$ (blue circles).

In panels (d), (e), and (f) of Fig. 2, we depict the corresponding spontaneous strain profiles corresponding to the cases in panels (a), (b), and (c) respectively. Importantly, far from the DW plane, the spontaneous values of the normal (triangles) and shear (squares) components of the strain converge to their respective values of the single domain state. However, the strained state of the DW causes various components of $\varepsilon_{ij}$ to grow or depress by large percentages to accommodate the electro- and rotostrictive coupling intrinsic in this structure. In the case of the 1/1 DW in (d), the value of the $\varepsilon_{zz}$ (in black) shrinks until eventually changing sign (smoothly) at the domain boundary. For the 2/2 DW, there is a large tensile strain in $\varepsilon_{xx}$ (in blue) growing by about a factor of three across the wall.

Also presented in Table 1 are the DW thicknesses associated to the corresponding order parameters, which differ between $\mathbf{P}$ and $\mathbf{A}$. We should note that the thicknesses of the DW corresponding to $\mathbf{P}$ and $\mathbf{A}$ differ. This arises because our resulting fit parameters are anisotropic (i.e. $H_{11}\ll H_{44}$) and also the presence of growth/decrease in non-switching components of $\mathbf{P}$ and $\mathbf{A}$ due to the roto- and electrostrictive coupling. Nevertheless, as seen in the table, the domain walls are quite thin ($2t_{k}\approx 0.05-0.5$ nm) which agrees quite-well with the available literature on BFO suggesting atomistically thin DWs Diéguez et al. (2013); Hlinka et al. (2017); Körbel and Sanvito (2020); Borisevich et al. (2010). The smaller value of the DW thickness in the $\mathbf{A}$ (as compared to $\mathbf{P}$) also shows good qualitative agreement with measurements from experiments using Z-contrast scanning transmission electron microscopyBorisevich et al. (2010).

We extend this type of analysis iteratively for the possible DWs listed in Table 1 so that we can converge our set of coefficients yielding reasonable $F_{\mathrm{DW}}$ values comparable to DFT; importantly, capturing the energy hierarchyDiéguez et al. (2013); Körbel and Sanvito (2020); Xue et al. (2014) predicted for the collection of walls. Our best estimates of the gradient coefficients found through our fitting procedure are presented in Table 2. We find that $H_{11}\ll-H_{12}<H_{44}$ in agreement with similar studies on BFOXue et al. (2014, 2021). This is an important relationship that results from harmonic models of antiferrodistortive cubic perovskite materials which has been connected to an asymmetry in the phonon bands at the R pointGesi et al. (1972); Stirling (1972); Cao and Barsch (1990). Another result from our fits is that the energy hierarchy yields $F_{\mathrm{DW}}(109^{\circ})<F_{\mathrm{DW}}(180^{\circ})<F_{\mathrm{DW}}(71^{\circ})$ for the lowest energy wallsLubk et al. (2009); Diéguez et al. (2013); Xue et al. (2014, 2021).

Now we turn to the AFM order present in BFO. To encapsulate the magnetic behavior of single crystalline BFO, we propose a continuum-approximation to the magnetic free energy density. We consider the total free energy density of the magnetic subsystem ($f_{\mathrm{mag}}$) to be a sum of the terms responsible for the nominally collinear AFM sublattices ($f_{\mathrm{sp}}$) and those producing the noncollinearity (canted magnetism) by coupling to the structural order ($f_{\mathrm{MP}}$). We first consider the magnetic energy due to the spin subsystem that is not coupled to the structural order,

where $D_{e}<0$ controls the strength of the short-range superexchange energy which favors the spins to have collinear AFM ordering Rezende et al. (2019). At our coarse-grained level of theory, we only consider the first nearest-neighbor exchange coupling which has been calculated from first-principles methodsXu et al. (2019) to be approximately 6 meV/f.u. corresponding to $D_{e}=-23.4\,\mathrm{meV}/\mathrm{nm}^{3}$ in our simulations. The second term describes the AFM non-local exchange stiffness proposed in Ref. [(15)] with $A_{e}=18.7$ meV/nm (or $3\times 10^{-7}$ ergs/cm). The third term corresponds to a weak single-ion anisotropy Weingart et al. (2012) with $K_{1}^{c}=2.2\times 10^{-3}\,\mu$eV/$\mathrm{nm}^{-3}$; this term reflects the cubic symmetry of the lattice and breaks the continuous degeneracy of the magnetic easy-plane into a six-fold symmetry.

The remaining terms are due to the magnetostructual coupling,

where

is due to the antisymmetric DMI which acts to break the collinearity by competing energetically with the first term of Eq. (21). It should be emphasized here that the local oxygen octahedral environments of adjacent Fe atoms underpins the DMI vector Ederer and Spaldin (2005); Fennie (2008); de Sousa and Moore (2009); Rahmedov et al. (2012); Meyer et al. (2022). Therefore, the $\mathbf{A}$ order parameter enables the DMI coupling. Ref. [(13)] provides an estimate of the DMI energy corresponding to $304$ $\mu$eV/f.u. It should be mentioned that the weak canting between $\mathbf{m}_{1}$ and $\mathbf{m}_{2}$ arises from a competition between $D_{e}$ and $D_{0}$ and that different estimates of their values can provide the same degree of canting of the sublattices provided they have the same ratio $D_{e}/D_{0}$. We come back to this in the next section.

BFO is an easy-plane antiferromagnetDixit et al. (2015), in which the magnetic sublattices lie in a plane defined by the direction of $\mathbf{P}$. We include the magnetocrystalline anisotropy termRezende et al. (2019) requisite for easy-plane AFMs as,

with the usual definition of $K_{1}>0$ enforcing the easy-plane condition for $\mathbf{m}_{\eta}$ with $\eta=1,2$. Using DFT methods, Dixit and co-workers Dixit et al. (2015), determined that the relative energy difference between aligning the magnetic sublattices along $\mathbf{P}$ or in the plane normal to $\mathbf{P}$ is $-2.0$ meV/f.u. Therefore, we choose $K_{1}=31.25\mathrm{meV}/\mathrm{nm}^{3}$ for our simulations.

We further couple the magnetic energy surface to the structural order by allowing the weak single-ion anisotropy to also depend on the antiphase tilts $\mathbf{A}$ Weingart et al. (2012) through,

which is in addition to the term in Eq. (21). The choice of $K_{1}^{c}>0$ and $0<a|\textbf{A}|^{2}<K_{1}^{c}$ corresponds to a small energy barrier between the 6-fold possible orientations of the weak magnetization $\mathbf{m}$ thus breaking the continuous degeneracy in the easy-plane. These coefficients can be obtained from DFT calculations as shown in Refs [(13)] and [(41)]. Therefore, we choose our coefficients (see Table 3) such that the relative energy density barrier for the six-fold symmetry is 0.01 meV/$\mathrm{nm}^{3}$ which is a reasonable approximation based on the aforementioned works. We find no influence of this choice of coupling constant on the results presented in this manuscript. The coefficients for $f_{\mathrm{sp}}$ and $f_{\mathrm{MP}}$ are listed in Table 3.

We should note that a long-period ($\lambda\approx 64$ nm) cycloidal rotation of the weak magnetization is often observed in BFO samplesBurns et al. (2020); Agbelele et al. (2017); Xu et al. (2021). It is possible to eliminate the cycloidal order by doping Sosnowska et al. (2002), epitaxial strain Bai et al. (2005); Sando et al. (2013); Burns et al. (2020), applied electric fields de Sousa et al. (2013), or by some processing techniques (i.e. via a critical film thickness) Burns et al. (2020) during synthesis. The spin-cycloid could be incorporated into our model by including coupling terms associated with a proposed spin-current mechanism Agbelele et al. (2017); Popkov et al. (2016); Xu et al. (2022). However, in order to provide the simplest model of the ME multiferroic effects, we have neglected them in this work.

In order to find the spin ground states in the presence of an arbitrary structural fields, we consider the Landau-Lifshitz-Bloch (LLB) equation Garanin and Chubykalo-Fesenko (2004) that governs the sublattices $\mathbf{m}_{\eta}$,

where $\alpha$ is the phenomenological Gilbert damping parameter and $\gamma$ is the electronic gyromagnetic coefficient equal to $2.2101\times 10^{5}$ rad. m A^-1 s^-1. 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/FeWeingart et al. (2012); Xu et al. (2019, 2022). The third term arises from the LLB approximation in the zero temperature limit where $\tilde{\alpha}_{\parallel}$ is a damping along the longitudinal direction of $\mathbf{m}_{\eta}$. We implement the LLB equation as a numerical resource in order to provide a restoring force and bind the quantities $\mathbf{m}_{\eta}$ to the unit sphere $(|\mathbf{m}_{\eta}|=1)$. In this context, we consider our spin subsystem to be at $T=0$ K in results presented throughout in this paper. In the Appendix, we provide a short derivation of the LLB torque in the zero temperature limit. We set $\tilde{\alpha}_{\parallel}=10^{3}$ in all results in this work to satisfy the constraint on $\mathbf{m}_{\eta}$.

To look for homogeneous spin ground states, we consider $\alpha=0.05$ and evolve Eq. 26 (utilizing the numerical approach described in Sec. II.3) until the relative change in the total energy computed from the summation of Eq. (21) and Eq. (22) between adjacent time steps is $\Delta F<10^{-8}$ eV/$\mu$s. Also, we stress that the influence of $\tilde{\alpha}_{\parallel}$ is negligible in all results presented in this work provided that its unitless value is around $10^{3}$ or above. To verify that our ground states predict the magnetic ordering consistent with the literature of BFO, we define two angular variables $\phi^{\mathrm{WFM}}=\cos^{-1}{\left(\mathbf{m}_{1}\cdot\mathbf{m}_{2}\right)}$ and $\theta_{\eta}=\cos^{-1}{\left(\mathbf{m}_{\eta}\cdot\hat{\mathbf{P}}\right)}$. The former tracks the degree of canting between the sublattices and the latter tracks the orientation of the magnetization with respect to $\hat{\mathbf{P}}=\mathbf{P}/P_{s}$, the magnetic easy-plane normal. As an example, we first set $\mathbf{P}\uparrow\uparrow\mathbf{A}$ along the [111] direction to be static. The time evolution (ringdown) of Eq. (26) is highlighted in Fig. 3(a) for $\theta_{\eta}$ showing that the sublattices have relaxed into the easy plane defined by $\hat{\mathbf{P}}$ with $\theta_{1}=\theta_{2}=90.0^{\circ}$ In (b) the time dependence of the canting angle $\phi^{\mathrm{WFM}}$ during the relaxation is shown. At the conclusion of the ringdown, $\phi^{\mathrm{WFM}}$ reaches a value of $\approx 1.22^{\circ})$. This demonstrates that the angular quantities $\{\theta_{\eta},\phi^{\mathrm{WFM}}\}$ detail an orthogonal system of the $\{\textbf{P},\textbf{m},\textbf{L}\}$ vectors as often discussed in the literature Heron et al. (2014).