## 2 The Phase-Field Model  
The present phase-field model describes ferroelectric domain configuration and evolution by polarization distributions \(\mathbf{P}=(P_{1},P_{2},P_{3})\). According to the Landau-Ginzburg-Devonshire theory, total free energy density for a ferroelectric single crystal is given by [14; 32; 33]  
\[F=F_{\text{Landau}}+F_{\text{grad}}+F_{\text{elas}}+F_{\text{elec}} \tag{1}\]  
in which \(F_{\text{Landau}},\ F_{\text{grad}},\ F_{\text{elas}},\ F_{\text{elec}}\) represent Landau energy, gradient energy, elastic energy, and electrostatic energy, respectively. Landau energy is expanded as a 4 th order polynomial of the polarization components \(P_{i}(i=1,\ 2,\ 3)\):  
\[F_{\text{Landau}}=\alpha_{1}\Big{(}P_{1}^{2}+P_{2}^{2}+P_{3}^{2}\Big{)}+ \alpha_{11}\Big{(}P_{1}^{4}+P_{2}^{4}+P_{3}^{4}\Big{)}+\alpha_{12}\Big{(}P_{ 1}^{2}P_{2}^{2}+P_{1}^{2}P_{3}^{2}+P_{2}^{2}P_{3}^{2}\Big{)} \tag{2}\]  
in which \(\alpha_{1},\alpha_{11},\alpha_{12}\) are Landau coefficients, and only \(\alpha_{1}\) is assumed to be temperature-dependent. Gradient energy, which represents the energy arising from the polarization gradient across domain walls, here is written as  
\[F_{\text{grad}}=\frac{1}{2}G_{11}\Bigg{(}\frac{\partial P_{i}}{\partial x_{j }}\Bigg{)}^{2},\ (i,j=1,2,3) \tag{3}\]  
in which \(G_{11}\) is gradient energy coefficient, and \(x_{i}\) is the \(i\)th component of the Cartesian coordinate system. Elastic energy, which is related to elastic strains originating from phase transitions or lattice accommodations, is written as  
\[F_{\text{elas}}=\frac{1}{2}c_{ijkl}e_{ij}e_{kl} \tag{4}\]  
in which \(c_{ijkl}\) is the elastic stiffness tensor. \(e_{ij}=e_{ij}-e_{ij}^{0}\) is the elastic strain, where \(e_{ij}\) is total strain, and \(e_{ij}^{0}\) is stress-free strain, induced by spontaneous polarization via the electrostrictive effect:  
\[e_{ij}^{0}=Q_{ijkl}P_{k}P_{l},\ (i,j,k,l=1,2,3) \tag{5}\]where \(Q_{ij}\) is the electrostrictive coefficient. Electrostatic energy is obtained from  
\[F_{\text{elec}}=-(E_{i}P_{i}+\frac{1}{2}\varepsilon_{b}E_{i}^{2}) \tag{6}\]  
in which \(E_{i}\) is electric field, and \(\varepsilon_{b}\) is the background dielectric constant. Temporal evolution of the polarization distribution is governed by the time-dependent Ginzburg-Landau equation:  
\[\frac{\partial P_{i}(x,t)}{\partial t}=-L\frac{\partial F}{\partial P_{i}(x,t)} \tag{7}\]  
in which \(L\) is the kinetic coefficient related to the domain-wall mobility. In the model setting, we consider an epitaxial grown BFO film on a rigid substrate in the \(\left(001\right)_{\text{c}}\) direction, as shown in Figure 1. Dimensions of the film are set as \(128\Delta~{}\times~{}128\Delta~{}\times~{}h_{\text{f}}\Delta\), in which \(h_{\text{f}}\) is the grid number for the film thickness, and \(\Delta\) is the grid size in the simulation. In the present work, \(\Delta\) is set as \(1\) nm. A periodical boundary condition is applied along the \(x_{1}\) and \(x_{2}\) axes. A thin film mechanical boundary condition is adopted in this simulation [33]. As for the electric boundary condition of the film, a short circuit condition is applied in this model [34]. The top and bottom electrodes used in this model are shown in Figure 1. It is believed that electrode materials also have effects on the properties of films, because of their differences in chemical and thermal stabilities, especially the Schottky barrier height, which has an influence on leakage current [35; 36; 37]. However, in this work, with the purpose of generally studying the thickness effect in electric-field-induced multi-domain switching, we did not take the effect of electrodes as the primary concern of this study.