# Transfer learning for accelerating phase-field modeling of ferroelectric domain formation in large-scale 3D systems

Kévin Alhada-Lahbabi, Damien Deleruyelle, Brice Gautier

► To cite this version:

Kévin Alhada-Lahbabi, Damien Deleruyelle, Brice Gautier. Transfer learning for accelerating phase-field modeling of ferroelectric domain formation in large-scale 3D systems. Computer Methods in Applied Mechanics and Engineering, 2024, 429, pp.117167. 10.1016/j.cma.2024.117167. hal-04814150

HAL Id: hal-04814150

https://hal.science/hal-04814150v1

Submitted on 2 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Transfer learning for accelerating phase-field modeling of ferroelectric domain formation in large-scale 3D systems

Kévin Alhada-Lahbabi *, Damien Deleruyelle, Brice Gautier

INSA Lyon, CNRS, Ecole Centrale de Lyon, Université Claude Bernard Lyon 1, CPE Lyon, INL, UMR5270, Villeurbanne, 69621, France

## ARTICLE INFO

**Keywords:**

- Ferroelectrics
- Phase-field
- Machine learning
- Transfer learning
- 3D large-scale systems
- Surrogate model
- Domain formation

## ABSTRACT

High-throughput phase-field simulations emerge as a compelling technique to predict the evolution of domain structures in ferroelectric materials. Despite their potential, their widespread use is impeded by significant computational costs, particularly for large 3D systems, limiting applicability in real-size scenarios and extensive parameter exploration. Here, we present a machine-learning approach tailored for accelerating 3D ferroelectric phase-field modeling. Our model exhibits remarkable accuracy in forecasting domain formation trajectories, providing stable predictions throughout entire simulation paths. Furthermore, we introduce a transfer learning strategy that facilitates the model extension to larger systems in a computationally tractable way, broadening its applicability to realistic scenarios. With a relative error below 10% when compared to high-fidelity methods, our surrogate serves as an efficient replacement for ferroelectric phase-field. By achieving speed-ups of up to 3300-fold compared to direct numerical solvers, the presented approach contributes to significantly accelerate the discovery and design of ferroelectric materials.

## 1. Introduction

Achieving a precise control of polarization in ferroelectric materials is crucial for their potential applications in modern nanoelectronic devices [1], ranging from nonvolatile memories [2,3] to negative capacitance devices [4] and neuromorphic computing [4–7]. Among ferroelectrics, PbZr<sub>x</sub>Ti<sub>1-x</sub>O<sub>3</sub>(PZT)-based ABX<sub>3</sub> perovskite ceramics have attracted substantial attention due their remarkable properties [8,9]. While substantial research has focused on controlling ferroelectric polarization [10–13], the evolution of domain states remains a phenomenon subject to nuanced interplay among external parameters [2,3]. These intricate microstructural paths underscore the necessity for a comprehensive understanding of ferroelectric domain dynamics.

Phase-field modeling has emerged as a prevalent computational technique for understanding ferroelectric materials [14–17]. Unfortunately, its broader adoption is hindered by the significant computational burden associated with iteratively solving an intricate system of partial differential equations (PDEs). Despite the emergence of robust alternatives to substitute the conventional finite element method (FEM) [18], such as isogeometric analysis (IGA) [19] and spectral methods [14], computational costs continue to escalate notably as the dimensions of the system expand. Particularly, this requirements surge when modeling large 3D systems or simulating long time scales, emphasizing the need for fast substitution models.

Nowadays, machine learning has become a compelling pathway for accelerating ferroelectrics material design and discovery. No-table contributions encompass the discovery of materials with tailored properties [20–23], refinement of microscopy characterization tools [24–26], and advancements in first-principles [27,28] and molecular dynamics computations [28]. Lately, machine-learning

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: kevin.alhada-lahbabi@insa-lyon.fr (K. Alhada-Lahbabi).</footnote>
surrogate models have gained significant attention for expediting phase-field simulation, attributed to their capacity to swiftly infer solutions for complex systems of PDEs [29–43]. Anticipating the microstructural evolution poses a significant challenge due to the co-evolution of material phases in these dynamic systems, characterized by simultaneous slow and fast dynamics, along with high spatiotemporal gradients. A common strategy is to design these models as explicit time-steppers, forecasting the subsequent state by relying on information from the current input state. These frameworks are often trained in a data-driven fashion, utilizing high-fidelity solvers to generate a training dataset consisting of pairs of features and labels.

These methodologies often employ dimensionality reduction techniques, such as Principal Component Analysis (PCA) or autoencoders (AE), to project data into a low-dimensional space. This process significantly enhances the efficiency of forecasting trajectory dynamics [30,35,38,40]. For example, Montes de Oca Zapiain et al. proposed a framework based on PCA and recurrent neural networks, which demonstrates high accuracy and remarkable speedups in learning the dynamics of the two-phase mixture problem in the latent space [30].

Alternatively, several research groups have utilized convolutional neural networks (CNNs) as surrogate models [29,38,40, 41]. Exploiting the image-based structure of phase-field data, they directly use the microstructure as 2D or 3D input arrays. Recently, Peivaste et al. introduced a U-net-based architecture for replacing the grain growth phase-field model. Their model accurately forecasts the order parameter in a two-grain phase-field, highlighting the efficacy of CNNs in emulating general phase-field microstructure evolution [38].

In contrast, a distinct strategy for developing machine-learning emulators with limited or no training data involves the use of Physics-Informed Neural Networks (PINNs) [33,44]. Due to their capacity to incorporate system-specific physical knowledge during training, PINNs have shown remarkable efficacy in addressing PDE-based problems [45,46]. In a recent breakthrough, Lu et al. introduced Deep Operator Networks (DeepONet) learning frameworks [32]. Inspired by the universal operator approximation theorem [47], DeepONet has demonstrated its ability to directly learn the underlying nonlinear operator directly from the data. Building upon this architecture, notable emulators for the phase-field problem have been developed, exploiting the free energy as a physics-informed loss function [31,48].

So far, only few studies have explored the acceleration of ferroelectric phase-field simulations. In a recent work, we introduced a CNN-based framework for inferring domain formation in PZT ferroelectric [29]. The model effectively controlled error accumulation, maintaining a remarkable 5% relative error over entire trajectory pathways unrolled from an initial random polarization state, and achieving speedups exceeding 600. At present, this methodology is confined to 2D scenarios with a relatively small $N_x \times N_y = 32 \times 32$ grid size, thereby limiting its ability to accommodate realistic 3D configurations. The exploration of these larger structures inherently involves time-consuming operations when employing direct numerical solvers. Those stringent requirements underscore the need for 3D surrogates and a methodology to scale up predictions for real-size structures in a computationally feasible manner.

In this work, we introduce a machine-learning approach to accelerate 3D phase-field modeling in PZT ferroelectric. Our model demonstrates remarkable accuracy in inferring entire ferroelectric pattern formation, yielding stable rollout predictions spanning a hundred of frames. In addition, we introduce a transfer learning strategy to effectively scale our approach to larger systems in a computationally tractable manner. To achieve this, the network undergoes initial training on a relatively extensive dataset of small-dimensional domains, followed by fine-tuning on a significantly smaller dataset featuring real-size structures. This approach reduced the necessity for time-intensive data generation of large training data, ensuring both efficiency and high accuracy in inferring full-scale systems.

With a relative error of under 10% while substituting for 90% of the traditional solver in unrolled predictions and achieving an acceleration of up to 3300, our model emerges as a fast real-time surrogate for the 3D ferroelectric phase-field. These outcomes broaden the application scope of the surrogate model to encompass more realistic scenarios and enable exploration of a larger parameter space using phase-field modeling.

## 2. Learning ferroelectric domain formation with machine-learning

In this endeavor, our principal objective is the development of a surrogate model for inferring pattern-formation across 3D ferroelectric domain formation. With classical phase-field methods, the microstructure evolution is given by the time-dependent Ginzburg–Landau (TDGL) equation [15],

$$ \frac{\partial P(r,t)}{\partial t} = -L \frac{\delta \psi}{\delta P(r,t)} \qquad (1) $$

where $P(r,t)$ is the spontaneous polarization, $L$ is a kinetic coefficient and $\psi$ is the total free energy (refer to Appendix A for detailed information regarding the phase-field methodology).

### 2.1. Surrogate model operation

In this study, the surrogate model operates as an explicit time-stepper, forecasting the subsequent state of the microstructure $X^{t_{k+1}}$ at time $t_{k+1}$, based on the current state $X^t$ at time $t_k$ through the operation :

$$ X^{t_{k+1}} = S(X^t), \qquad (2) $$

where $S$ denotes the neural network's forward computation.
We illustrate our approach using a PZT thin film with vertical (001) orientation, aligning with a well-established configuration documented in the literature [49–51]. In this system, we aim to predict the dynamic of vertical c ferroelectric domains. To achieve this, we adopt the polarization variable $P_x^{ik}$ as a succinct representation of the microstructure at any given time $t_k$:

$$X^{ik} = P_x^{ik} \qquad (3)$$

The surrogate must then adeptly learn to predict the microstructure one time-step $\Delta t$ ahead:

$$P_x^{ik+1} = S(P_x^{ik}) \qquad (4)$$

By iteratively using its predictions, the network can generate rollout trajectories from the initial state at $t_0$ to the final time $t_N$ across times $t = \{t_0, \dots, t_N\}$, formally expressed as:

$$P_x^{iN} = S^N(...S(P_x^{i0})) \qquad (5)$$

To emulate the incremental update of the polarization field governed by the TDGL equation,

$$P_x^{i(k+1)} = P_x^{i(k)} + \frac{\partial P_x^{i(k)}}{\partial t} \Delta t = P_x^{i(k)} + dP_x^{i(k)}, \qquad (6)$$

we compelled the surrogate model to predict the polarization residual $dP_x^{ik}$ [29]. The model outputs were then determined by residual learning as,

$$P_x^{i(k+1)} = S(P_x^{i(k)}) = P_x^{i(k)} + dP_x^{i(k)} \qquad (7)$$

where $P_x^{ik}$ represents the current microstructural input state, and $dP_x^{ik}$ is the output obtained through the network's forward operation.

## 2.2. Surrogate model architecture

For the surrogate representation, we employed a 3D CNN based on encoder-decoder, similar to a prior work [29]. Specifically, we used a 3D Unet with skip-connections, a well-established architecture in computer vision [52]. The network is designed to receive the current microstructural state as inputs, denoted $X^{ik} = P_x^{ik}$, with a system shape of $[N_x, N_y, N_z]$. Subsequently, the framework predicts the ensuing state $X^{i(k+1)} = P_x^{i(k+1)}$ at the following time step. Further details regarding the network architecture and hyperparameters can be found in Appendix B.

## 2.3. Training loss error

To learn how to replicate the evolution of domain morphology given by 3D ferroelectric phase-field, the surrogate model undergoes supervised training utilizing the $L_2$ error loss function. The training loss is defined as:

$$\mathcal{L} = L_2(P_{xPF}^{i(k+1)}, S(P_x^{ik})) = L_2(P_{xPF}^{i(k+1)}, P_x^{i(k+1)}) \qquad (8)$$

In this context, $P_{xPF}^{i(k+1)}$ signifies the microstructure labels derived from high-fidelity phase-field simulations, while $S(P_x^{ik}) = P_x^{i(k+1)}$ represents the model outputs. Details on the computation of the $L_2$ error loss function can be found in Appendix C.

# 3. Results

In this section, we investigate the prediction of 3D ferroelectric domain formation through a transfer learning approach enabling efficient scalability to larger systems. We achieve this by first training the machine learning framework on relatively small systems ($N_x \times N_y \times N_z = 16 \times 32 \times 32$). Then, we employ transfer learning to extend predictions to significantly larger systems (ranging from $16 \times 64 \times 64$ to $16 \times 256 \times 256$). Results are first presented for the initial training on the $16 \times 32 \times 32$ systems, followed by those obtained through transfer learning on progressively larger systems. Finally, the study elucidates the model's capability to directly extrapolate predictions for even larger systems ($N_x \times N_y \times N_z = 16 \times 512 \times 512$) without any additional fine-tuning, thereby enabling prediction in real-size systems without additional training costs.

## 3.1. Domain formation prediction in $N_x \times N_y \times N_z = 16 \times 32 \times 32$ systems

### 3.1.1. Dataset

The training of data-driven surrogate models typically necessitates the generation of expansive feature-label datasets derived from computationally expensive high-fidelity simulations. To mitigate this challenge, the initial model training is carried out on relatively small $N_x \times N_y \times N_z = 16 \times 32 \times 32$ systems, thereby significantly reducing the time required to generate the initial substantial training dataset.

A total of 1400 phase-field simulations of ferroelectric domain formation were carried out using the methodology detailed in Appendix A. The simulations were initiated from randomly domain states, with the polarization $P_x(r)$ being initialized at each grid point through random sampling from a uniform distribution with values ranging from $-P_{c0}$ to $P_{c0}$. The morphological state evolution
<figure><img src="image_4.png" /><figcaption>Fig. 1. (a) Illustration of the evolution of the polarization field $P_x$ throughout the evolution of domain morphology. (b) Visualization of selected domain trajectory examples in the PCA lower-dimensional space built on the first three principal components. (c) Distribution of ferroelectric polarization in the training dataset. (d) Machine-learning encoder-decoder framework used to learn domain formation forecasting. The surrogate model takes the microstructure at time $t_k$ as input and predicts the subsequent microstructural state at time $t_{k+1}$, computing the polarization through residual learning. The model is trained using a data-driven approach, minimizing the $L_2$ training loss between the model outputs $S(P_x^{t_k})$ and the corresponding high-fidelity phase-field labels $P_{s_{PF}}^{t_{k+1}}$.</figcaption></figure>

was continued until reaching domain equilibrium, which occurred approximately after a duration of 1000Δt. Here, Δt denotes the time-step employed in the phase-field simulations. The microstructure state was then recorded at uniformly spaced time intervals of 20Δt, resulting in trajectories comprising 50 frames per simulation ({$t_0, \dots, t_{50}$}). Each trajectory is represented as a tensor of shape (50, 16, 32, 32). An overview of a typical polarization configuration trajectory drawn from the training dataset is presented in Fig. 1a. The dataset was then partitioned into a training dataset (1000 simulations), a validation dataset (200 simulations), and a test dataset (200 simulations).

Phase-field simulations inherently produce highly nonlinear trajectories, contributing to the complexity of inferring the domain dynamics. To ensure the diversity of the generated data for effective model generalization, Principal Component Analysis (PCA) was employed, following similar methodologies in analogous studies [29,30]. Further details on PCA are outlined in the Appendix D. Fig. 1b showcases the representation of five random selected polarization trajectories from the training data set in the low-dimensional space defined by the first three principal components. The visualization underscores the intricate microstructural pathways from their initial states, emphasizing the diversity encapsulated in the training examples.

Additionally, an analysis of the distribution of ferroelectric polarization is depicted in Fig. 1c. This distribution reveals two primary components within the training dataset, corresponding to the values of spontaneous polarization $P_{c0}$ and $-P_{c0}$. Finally, a schematic representation of the machine learning framework is presented in Fig. 1c. This illustration elucidates the model operational framework as an explicit time-stepper, predicting the next domain state $P_x^{t_{k+1}} = S(P_x^{t_k})$ through residual learning, relying on the current input state $P_x^{t_k}$.

The model was trained on the 1200 structures comprising the training/validation dataset over 100 epochs, ensuring convergence of both training and validation losses. The corresponding training history is illustrated in Fig. 2a. Afterwards, the model's performance was assessed on direct one-step and rollout predictions for the 200 simulations within the test dataset.

### 3.1.2. Evaluation of one-step prediction

For the evaluation of one-step predictions, we present the $L_2$ mean squared error (MSE) and $L_1$ mean absolute error (MAE) metrics in Table 1 (see Appendix C for detailed metric computation procedures). The model exhibits notable accuracy in inferring
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** (a) Progression of the $L_2$ training error in the training and validation datasets during the initial training on the 16 × 32 × 32 systems. (b) Evolution of the $L_2$ training and validation errors for the model fine-tuned on the 16 × 64 × 64, 16 × 128 × 128, and 16 × 256 × 256 systems using transfer learning. The training and validation errors are depicted with solid and dashed lines, respectively.</figcaption></figure>

<table><caption>**Table 1**
Performance evaluation of the model for vertical domain formation with one-step prediction.</caption><thead><tr><th>Field</th><th><i>L</i><sub>2</sub> MSE</th><th><i>L</i><sub>1</sub> MAE (C/m<sup>2</sup>)</th></tr></thead><tbody><tr><td><i>P</i><sub>x</sub></td><td>1.06 × 10<sup>-5</sup></td><td>4.06 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** (a) Evolution of the $L_2$ error and (b) MARE during the unrolled prediction of domain formation. For each starting frame, complete domain formation trajectories are predicted until the final frame at time $t_{50}$. Results are averaged over 200 test dataset simulations, with the solid line indicating the mean and the shaded region representing the interquartile range between the 25th and 75th percentiles.</figcaption></figure>

the subsequent microstructural state, achieving a MSE of $1.06 \times 10^{-5}$ and a quantitative MAE of $4.06 \times 10^{-3}$ C/m$^2$ for the polarization field. However, it is essential to underscore that the one-step prediction error, while indicative of immediate forecasting precision, does not comprehensively depict the model's performance in forecasting over long time scales, where the model must unfold entire simulations. A thorough assessment is therefore crucial to ascertain the model's ability to prevent the gradual accumulation of errors over time and ensure stable predictions.

### 3.1.3. Evaluation of unrolled prediction

For the rollout trajectories evaluation, we initiated simulations from various initial frames to assess the model's robustness in potential error-accumulation scenarios. Our primary goal is to replace the directed numerical solver as early as possible in the process to accelerate simulations. Therefore, we restricted our analysis to initial frames corresponding to time $t_5$ to $t_{25}$ and assigned the model with the prediction of the entire simulation trajectory up to time $t_{50}$. The mean, 25th and 75th quartiles of the MSE and macro average relative error (MARE) computed over the unrolled trajectories are reported for the polarization field in Fig. 3 for each starting frame.

Initially, we can observe that the prediction error rises when simulations commence from earlier states. This phenomenon can be attributed to two primary factors. Firstly, the less defined initial polarization configuration at the outset of domain formation complicates the anticipation of the domain trajectory, given the higher spatiotemporal gradients characterizing the domain dynamics at this stage. Secondly, the gradual error accumulation is a consequence of the increased number of inference steps required, a phenomenon commonly encountered in surrogate networks, where incremental prediction errors may accumulate over consecutive steps [37]. Nonetheless, the model demonstrates noteworthy robustness and stability, effectively preventing the accumulation of errors, even in the early stages of simulation initiation. In particular, the MARE falls below 10% when the surrogate starts from the frame at time $t_6$, and dips below 5% from an initial time $t_8$. These error levels fall below an acceptable threshold, especially for the challenging task of replacing about 90% of the high-fidelity solver, showcasing the model's effectiveness in balancing accuracy and computational efficiency.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison between unrolled predictions and corresponding high-fidelity ground truth, initiated from the frame at time $t_5$, for a test trajectory. Microstructural polarization states ($P_x$) are visualized at (a) the initial ($t_5$) and (b-c) final ($t_{50}$) times. (d) PCA representation of the prediction and high-fidelity trajectories, using the first three principal components.</figcaption></figure>

### 3.1.4. Illustration of forecasted trajectories

Fig. 4 depicts a forecasting example starting at time $t_5$, comparing model predictions with the ground truth of a test trajectory. Despite the intricate initial configuration, the model displays accurate roll-outs without significant error growth. While minor errors persist throughout the inferred trajectory, their impact on the final state remains negligible. Furthermore, to enhance the visualization of the model's accuracy along the prediction trajectory, both the ground truth and model predictions at each time step ($t_5, ..., t_{50}$) are depicted in Fig. 4d in the PCA low-dimensional space based on the first three principal components. This PCA representation offers valuable insights into the morphology evolution pathways taken by the surrogate. Specifically, the framework accurately replicates the long-term evolution trajectories of the high-fidelity solver, demonstrating similar inference dynamics. Interestingly, the main deviation in trajectories occurs in the first predicted frames, primarily due to the intricacies of the domain structure at this early stage. Despite these initial deviations, the model consistently preserves the overall trend inherent to domain formation, highlighting its intrinsic understanding of the underlying phase-field evolution equations.

It is noteworthy that trajectory accuracy could be enhanced by initiating predictions from a later frame. However, this improvement would come at a higher computational, as a larger portion of the simulations would involve the time-consuming drawbacks of high-fidelity solvers. For instance, starting from frames beyond $t_{10}$ yields a MARE below 2%, indicating very high precision in forecasting more advanced trajectories. Overall, these findings underscore the machine learning emulator's capability to substitute the direct numerical solver early in the simulation, providing a flexible trade-off between computation speed and accuracy.

## 3.2. Transfer learning prediction scaling to larger systems

### 3.2.1. Dataset

To faithfully model real-life setups, the phase-field method often necessitates operation in large systems. However, a notable challenge in developing surrogates for such extensive systems arises from the prohibitive computational burden associated with generating a substantial training dataset using classical methods. Here, we introduce a transfer learning-based approach to bridge this computational gap, aiming to efficiently scale up predictions within a reduced and reasonable timeframe.

Transfer learning is a widely employed technique for efficiently addressing material discovery and design tasks [53–55]. In the context of domain formation forecasting, the core idea involves utilizing a pre-trained surrogate model initially trained on small 16 × 32 × 32 structures and subsequently fine-tuning its parameters with a reduced quantity of larger systems. To facilitate this, datasets are generated for domains of sizes 16 × 64 × 64, 16 × 128 × 128, and 16 × 256 × 256, each comprising 300 phase-field simulations. The datasets are then partitioned into a small training set (50 simulations), a validation set (50 simulations), and a test set (200 simulations). Subsequently, for each system size, the pre-trained model undergoes a fine-tuning process on the training dataset, which includes a relatively low number of large structures (refer to Appendix E for training details). The respective training histories for various system sizes are presented in Fig. 2b. In each scenario, both the training and validation losses exhibit a swift convergence, underscoring the surrogate's proficiency in scaling predictions to the dimensions of the targeted computational regions. Despite variations in training errors depending on the target system size, notably with the lowest error for the 16 × 64 × 64 domains, it is essential to note that the training error may not fully reflect unrolled simulation accuracy. Such variations may stem from statistical fluctuations within the dataset due to the relatively low number of training examples or potentially from a more effective scaling of the pre-trained model to domains with similar dimensions.

### 3.2.2. Evaluation of unrolled prediction

We present the results from the roll-out evaluation after training on the 200 test simulations, highlighting the dependence on the initial starting frame for each system size in Fig. 5. Overall, transfer learning demonstrates a comparable level of accuracy to the initial training on the 16 × 32 × 32 domains, even as the system dimensions increase. Specifically, the MARE remains consistently low, even when initiated from early frames. Across all scenarios, it is contained below or very close to 10% when starting from the frame at time $t_5$, and rapidly decreases below 5% when commencing from time $t_8$.
<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) Illustration of the $L_2$ error evolution and (b) MARE in the unrolled prediction for models trained using transfer learning on larger system sizes, specifically 16 × 64 × 64, 16 × 128 × 128, and 16 × 256 × 256. Results are averaged over 200 examples from the respective test datasets, with the model tasked to complete the simulation until $t_{50}$.</figcaption></figure>

### 3.2.3. Illustration of forecasted trajectories

An illustration of the predicted and high-fidelity polarization trajectories, obtained after fine-tuning the models, is provided in Fig. 6 for a 16 × 64 × 64 test system. The model was provided with the initial microstructural state at time $t_5$ and tasked with unrolling the entire simulation up to $t_{50}$. The results demonstrate high precision throughout the domain pathways, indicating the model's effective learning to scale predictions to larger inputs and capturing the dependency of long and short-range interactions with system sizes. Nevertheless, the surrogate may produce errors concerning domain wall dynamics, as illustrated in the MAE prediction error representation (Fig. 6c). Here, it can be noticed that the model occasionally fails to accurately forecast the precise location of domain walls, leading to error accumulation in these regions. However, the surrogate adeptly predicts the sharp gradients in the domain wall regions, as evidenced in Fig. 6d, where the predicted and reference polarization evolutions are depicted in a one-dimensional fashion across a domain wall region during domain formation.

We further elucidated the evolution of prediction errors at the domain wall depending on the initial state for domain prediction. To illustrate this, the final domain state at the wall region shown in Fig. 6d is depicted in Fig. 7 at time $t_{50}$ for unrolled predictions initiated from $t_0, t_5, t_{10}$, and $t_{15}$.

As the model starts from a more advanced initial microstructure, we notice a progressively closer match in domain wall location and profile to the traditional phase-field, especially noticeable between $t_0$ and $t_{10}$. This convergence further refines until achieving precise agreement at $t_{15}$, offering a comprehensive perspective on the trade-off between accuracy and computational speed-up provided by such surrogate models.

Illustrations of the microstructural evolution on 16 × 128 × 128 and 16 × 256 × 256 test systems are provided in Fig. 8a and b, respectively, following fine-tuning of the model on these respective domain sizes. While the surrogate delivers faithful and accurate predictions of the overall domain structure, it is crucial to recognize the existence of minor errors in the domain shape, particularly noticeable in the 16 × 256 × 256 system. This underscores the inherent challenge associated with predicting the dynamical response of more complex systems. Nevertheless, even when initiating predictions from such an early state, the surrogate model consistently yields final states closely resembling the ground truth microstructure.

While leveraging the high-fidelity solver for additional initial steps could mitigate model deviation, the accompanying computational cost, particularly in those larger systems, would be significant. Overall, these findings emphasize the feasibility of efficiently scaling up predictions using a limited amount of extensive training data while maintaining a satisfactory level of accuracy, enabling the early substitution of high-fidelity solvers.

### 3.3. Extrapolation to larger system sizes

While the previously discussed fine-tuning based transfer approach could be extended to even larger systems, the associated costs of data generation rapidly become prohibitive, even with a reduced training dataset. Therefore, this section explores the direct application of pre-trained models on 16 × 256 × 256 systems to extrapolate domain formation in larger 16 × 512 × 512 domains without additional fine-tuning.

To achieve this, a test dataset comprising 100 phase-field simulations of domain formation was generated on a grid with dimensions 16 × 512 × 512. Subsequently, the previously fine-tuned model on the 16 × 256 × 256 systems was directly assessed in a roll-out context on the 16 × 512 × 512 domains. The corresponding $L_2$ error and MARE are reported in Fig. 9 when the model carries out the domain formation trajectory until $t_{50}$, with starting frames varying from $t_5$ to $t_{20}$.

Our findings demonstrate that the surrogate model effectively extrapolates the microstructure of previously unobserved larger systems while maintaining high accuracy, achieving an MARE below 5% even when initiating from the early stages of domain formation. An example of the forecasted microstructure, through unrolled prediction from an initial state at $t_5$, is presented in Fig. 10 along with its corresponding ground truth. The extrapolated prediction exhibits a close alignment with high-fidelity phase-field computations, faithfully reproducing the complex domain and domain wall patterns across the ferroelectric thin films. Thus, it becomes feasible to make predictions for real-world system sizes without incurring additional computational costs associated with generating training examples of extensive dimensions.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Illustration of domain prediction scaled to a larger system of 16 × 64 × 64 using transfer learning. Depiction of unrolled predictions initiated from the frame at time $t_5$ (a), corresponding high-fidelity ground truth (GT) (b), MAE prediction error (c), and (d) one-dimensional cross-section across a domain wall region over the simulation time.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Evolution of the predicted and reference polarization, along with the MAE, across a one-dimensional cross-section of a domain wall region as a function of the starting frame. For each starting frame, the simulation is unrolled until $t_{50}$, where the final domain wall configuration is depicted.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Illustration of domain prediction scaled to larger systems of (a) 16 × 128 × 128 and (b) 16 × 256 × 256 using transfer learning. In each case, the comparison between the unrolled prediction and the corresponding high-fidelity ground truth, initiated from the frame at time $t_5$ and run until the final timestep $t_{50}$, is shown for a trajectory selected from the test dataset.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** (a) Illustration of the $L_2$ error evolution and (b) MARE in the unrolled prediction for models trained using transfer learning on a larger system size of 16 × 512 × 512, extrapolated without fine-tuning from a pre-trained model on 16 × 256 × 256 systems. Results are averaged over 100 test examples, where the model is tasked with completing the simulation up to $t_{50}$.</figcaption></figure>

### 3.4. Computation efficiency

Accurately quantifying the computational speed-ups facilitated by a surrogate model poses a non-trivial challenge, considering the inherent computational demands involved in the surrogate creation process, including dataset generation and model training. In this context, we initiate our analysis by presenting the inference times of both the high-fidelity phase field and the surrogate for
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Illustration of domain prediction extrapolated to a larger 16 × 512 × 512 system using a pre-trained model on structures with dimensions 16 × 256 × 256. The comparison between the unrolled prediction and the high-fidelity phase-field is shown at the final timestep $t_{50}$, with the simulation initiated from the frame at time $t_5$.</figcaption></figure>

<caption>**Table 2**</caption>

Computational efficiency of the machine learning framework across system sizes. Reported times include Phase-Field (PF) and Machine Learning (ML) inference time per frame, acceleration (CPU-GPU), dataset generation, and model training.

<table><thead><tr><th rowspan="2">System dimensions<br>N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub></th><th rowspan="2">PF Inference (s/frame)<br>Phase-Field</th><th colspan="2">ML Inference (s/frame)</th><th colspan="2">Acceleration</th><th rowspan="2">Dataset (days)</th><th rowspan="2">Training (h)</th></tr><tr><th>CPU</th><th>GPU</th><th>CPU</th><th>GPU</th></tr></thead><tbody><tr><td>16 × 32 × 32</td><td>11.4</td><td>9.7 × 10<sup>-3</sup></td><td>3.6 × 10<sup>-3</sup></td><td>1170</td><td>3170</td><td>6.6</td><td>16.5</td></tr><tr><td>16 × 64 × 64</td><td>44.3</td><td>2.6 × 10<sup>-2</sup></td><td>1.3 × 10<sup>-2</sup></td><td>1710</td><td>3290</td><td>1.3</td><td>0.9</td></tr><tr><td>16 × 128 × 128</td><td>170.4</td><td>1.7 × 10<sup>-1</sup></td><td>5.1 × 10<sup>-2</sup></td><td>1010</td><td>3340</td><td>4.9</td><td>1.7</td></tr><tr><td>16 × 256 × 256</td><td>618</td><td>8.4 × 10<sup>-1</sup></td><td>2.0 × 10<sup>-1</sup></td><td>732</td><td>3030</td><td>17.8</td><td>5.3</td></tr></tbody></table>

each system size, as depicted in Table 2. To ensure representative outcomes, the complete test dataset, comprising 200 simulations in each scenario, was regenerated using the surrogate model.

Assuming the surrogate is utilized to replace the phase-field from time $t_5$, we further elucidate the acceleration factor for each system size in the case of an unrolled simulation. The evaluation encompassed both CPU and GPU assessments, ensuring a fair comparison with simulations conducted using traditional solvers. The results, reported in Table 2, underscore the substantial computational speed-ups provided by the surrogate, ranging from approximately 1000 to 3300, contingent on the material employed.

Although our approach demonstrates remarkably fast inference times, it is imperative to acknowledge that the creation of such a surrogate entails an initial computational investment. To provide a comprehensive overview of the costs associated with model creation, we present the dataset generation times and model training times in Table 2. Leveraging the transfer learning approach, we strive to minimize the time required for generating datasets, particularly for large system sizes.

The outcomes presented in this section are intricately linked to the specifics of our computational material and the phase-field methodology. In this study, training and GPU assessments employed an NVIDIA GeForce RTX 3080 with 10 GB of RAM, while dataset generation and CPU assessments were conducted using an Intel i9 CPU clocked at 5.1 GHz. Different computational configurations, resources, or numerical methodologies, such as choosing finite-element methods over spectral methods to solve the PDEs set, may result in divergent execution times.

## 4. Discussion

This study highlights the surrogate model's potential as a highly efficient and accurate alternative to conventional phase-field modeling, demonstrating remarkable precision in inferring entire unrolled simulations of ferroelectric 3D domain formation. The framework accurately captures intricate domain dynamics, forecasting morphological state evolution and polarization across diverse initial states and pathways.

With the ability to mitigate trajectory deviation, the framework consistently keeps the relative error below 10% when initiating the simulation from time $t_5$, replacing 90% of the direct numerical solver. Additionally, prediction relative errors quickly drop below 5% when initiating the simulation from time $t_8$ and beyond, providing even more accurate results for 85% of the trajectory. This comprehensive overview guides the trade-off between accuracy and computational speedup, based on the desired application.

In terms of computational efficiency, the surrogate consistently achieves substantial speed-ups, ranging from approximately 1000 to 3300. The framework effectively alleviates the computational burden associated with 3D phase-field modeling, representing a significant advancement for handling realistic and resource-intensive scenarios. However, it is crucial to note that developing such surrogates requires a substantial time investment, particularly during the dataset generation phase using the high-fidelity solver. To address this challenge, strategies involve incorporating data augmentation techniques, including physically plausible transformations, to expedite the generation of necessary training simulations. Despite the significant initial time commitment for
dataset generation, it is crucial to note that this phase remains a one-time expenditure. Subsequent applications of the surrogate to efficiently address direct and inverse problems incur minimal time costs, highlighting the long-term efficiency of our approach.

As an illustrative application example of our framework in accelerating direct problem-solving, one could consider the analysis of domain stability following tip-induced mechanical switching. In these scenarios, high-throughput phase-field is commonly leveraged to elucidate polarization state stability after a mechanical tip loading [12,56]. Depending on energetic considerations, the switched domain can either stabilize or undergo a reversal back to its initial state. The presented framework could significantly expedites this relaxation phase, enabling a broader exploration of the influence of initial morphological configurations.

Beyond its application in direct simulations, the surrogate model can be effectively employed to address inverse problems, such as determining optimal defect arrangements to achieve desired material properties. In this context, the framework can be extended to encompass 3D defects or charge arrangements by incorporating them as additional input tensors. Leveraging the differentiability and speed enhancements of the surrogate model, it becomes feasible to expeditiously sample diverse defect configurations across numerous simulations, aligning with target properties through back-propagation within a condensed timeframe.

Further research efforts could potentially encompass the development of surrogates for electrical and mechanical switching phenomena in ferroelectric materials. Such models must accurately incorporate electrical and elastic boundary conditions at each time step to ensure faithful reproduction of domain switching trajectories. Given the prevalence of phase-field methods in analyzing experimental setups involving electrical or mechanical stimuli, this advancement would significantly broaden the applicability of the surrogate. Moreover, an extension to materials exhibiting a more intricate polarization arrangement is conceivable. For example, scenarios involving polarization evolving along the three directions ($P_x, P_y, P_z$) and resulting in $a_1/a_2/c$ domains could be addressed by augmenting the number of input and output channels of the model, as demonstrated in a 2D prior study [29].

## 5. Conclusion

In conclusion, we have introduced a machine learning framework that accurately forecasts domain formation in 3D ferroelectric phase-field. Achieving less than 10% relative error when replacing 90% of the high-fidelity solver, the model proves itself as a reliable surrogate for ferroelectric modeling. The use of transfer learning further enables scaling the model to real-size systems in an efficient computational manner. By achieving substantial computational acceleration, the presented approach holds promise in accelerating the design and discovery of ferroelectric materials.

## CRediT authorship contribution statement

**Kévin Alhada-Lahbabi:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Damien Deleruyelle:** Writing – review & editing, Validation, Supervision, Funding acquisition. **Brice Gautier:** Writing – review & editing, Validation, Supervision, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

This work has received the financial support of IPCEI France 2030 Programs POI and eFerroNVM together with ANR-23-CE24-0015-01 ECHOES project.

## Appendix A. Phase-field modeling

In traditional ferroelectric phase-field modeling, the evolution of ferroelectric polarization microstructure during 3D domain formations is described by the TDGL equation [14–16]:

$$ \frac{\partial P_i(r, t)}{\partial t} = -L \frac{\delta \psi}{\delta P_i(r, t)}, \quad (i = 1, 2, 3) \qquad (A.1) $$

where $P_i(r, t)$ is the spontaneous polarization, $L$ is a kinetic coefficient, $\psi$ is the total free energy and $r = (x, y, z)$ denotes the spatial vector in 3D. The total free energy includes the bulk, gradient, electric, and elastic free energy density:

$$ \psi = \int_V (\psi_{bulk} + \psi_{grad} + \psi_{elec} + \psi_{elas}) \delta V \qquad (A.2) $$

The polarization was updated at each time step following the explicit scheme:

$$ P_i(r, t + \Delta t) = P_i(r, t) + \frac{\partial P_i(r, t)}{\partial t} \Delta t \qquad (A.3) $$
where $\Delta t$ is the time step for integration.

The bulk energy is described by:

$$ \psi(P_i) = \alpha_i P_i^2 + \alpha_{ij} P_i^2 P_j^2 + \alpha_{ijk} P_i^2 P_j^2 P_k^2, \qquad (A.4) $$

where $\alpha_i$, $\alpha_{ij}$ and $\alpha_{ijk}$ are the second-, fourth- and sixth-order PZT Landau coefficients, which are taken from the literature [57], and $\alpha_1 = \frac{T-10}{2\epsilon C}$ refers to the dielectric permittivity $\epsilon$, the Curie temperature $T_0$ and the Curie constant $C$.

The energy caused by the domain walls is described by the gradient energy, which in a cubic system is calculated by:

$$ \psi_{grad}(\nabla P_i) = \frac{1}{2} G_{ijkl} \left( \frac{\partial P_i}{\partial x_j} \frac{\partial P_k}{\partial x_l} \right), \qquad (A.5) $$

where $G$ represents the gradient energy coefficient tensor.

The electric energy is given by:

$$ \psi_{elec}(P_i, E_i) = -P_i E_i - \frac{1}{2} \epsilon_0 \epsilon_r E_i E_j, \qquad (A.6) $$

where $E_i = -\nabla_i V$ is the electric field inside the cell obtained by solving the electrostatic equilibrium for the electrostatic potential $V(r)$:

$$ \Delta V(r) = - \frac{\rho(r) - \nabla \cdot \mathbf{P}(r)}{\epsilon_0 \epsilon_r} \qquad (A.7) $$

where $\rho$ is the electric charge and in which a depolarization field is induced by the polarization charges term $-\nabla \cdot \mathbf{P}(r)$. The electrostatic equilibrium is determined using the fast Fourier transform spectral method, as detailed in prior studies [29,58]. Periodic boundary conditions are applied along the in-plane axes (y and z). Along the vertical x axis, a short-circuit electrostatic condition is imposed, employing the discrete sine transform (DST), assuming complete charge screening of both the bottom electrode and the top surface:

$$ V(0, y, z) = V(x_{top}, y, z) = 0 \quad V \qquad (A.8) $$

The elastic energy density is described by:

$$ \psi_{elas}(P_i, \epsilon_{ij}) = \frac{1}{2} C_{ijkl} (\epsilon_{ij}(\mathbf{r}) - \epsilon_{ij}^0(\mathbf{r})) (\epsilon_{kl}(\mathbf{r}) - \epsilon_{kl}^0(\mathbf{r})) \qquad (A.9) $$

where $C$ is the elastic stiffness tensor, $\epsilon$ is the total strain and $\epsilon^0$ is the electrostrictive strain caused by the polarization as

$$ \epsilon_{ij}^{0}(\mathbf{r}) = Q_{ijkl} P_{k}(\mathbf{r}) P_{l}(\mathbf{r}) \qquad (A.10) $$

where $Q$ is the electrostrictive tensor. The total strain contains the homogeneous and heterogeneous strains:

$$ \epsilon_{ij}(\mathbf{r}) = \epsilon_{ij}^{s} + \delta\epsilon_{ij}(\mathbf{r}) \qquad (A.11) $$

which is linked to the mechanical displacement $u_i$ by:

$$ \delta\epsilon_{ij}(\mathbf{r}) = \frac{1}{2} \left( \frac{\partial u_j(\mathbf{r})}{\partial x_j} + \frac{\partial u_i(\mathbf{r})}{\partial x_i} \right) \qquad (A.12) $$

The mechanical equilibrium equation $\sigma_{ij,j} = 0$ is given by (using Einstein notation):

$$ C_{ijkl} \frac{\partial^2 u_k(\mathbf{r})}{\partial x_j \partial x_l} = C_{ijkl} \frac{\partial \epsilon_{ij}^0(\mathbf{r})}{\partial x_j} \qquad (A.13) $$

which is solved in Fourier space, looking for a periodic solution, as detailed in previous study [29].

As an exemplary model, and with the objective of maintaining a general setting, we chose to illustrate our approach using PZT thin film vertically oriented, a configuration commonly encountered in reference studies [49–51]. The PZT parameters were employed as established in prior literature [49,59,60]. Details regarding these parameters and their normalization [60] can be found in Tables A.1 and A.2.

Simulations were performed on 3D grids with varying dimensions ($N_x \times N_y \times N_z = 16 \times 32 \times 32 \sim 16 \times 256 \times 256$), as outlined through the manuscript. All grids used uniform spacing $\Delta x/l_0 = \Delta y/l_0 = \Delta z/l_0 = 1$, where $l_0 \approx 1$ nm (based on $G_{110}$ and $\alpha_0$ at 25 °C, $\alpha_0$ from Ref. [59]). Gradient energy coefficients were employed based on Ref. [59]: $G_{11}/G_{110} = 0.6$, $G_{12}/G_{110} = 0$, $G_{44}/G_{110} = 0.3$. A typical –1% in-plane misfit strain was applied to the films, as commonly encountered in vertically oriented PZT [49]. The time step was $\Delta t = 0.02t_0$, where $t_0 = 1/(\alpha_0 L_0)$.

## Appendix B. Model architecture and hyperparameters details

The machine learning model employed in this study is a CNN-based encoder–decoder, specifically following the U-Net architecture.

The input data has a shape of $[N_x, N_y, N_z, 1]$, corresponding to the microstructural state input. The initial phase of the network involves a CNN encoder comprising five downsampling blocks. Each of the first four downsampling blocks (Down 1–4) consists of a
<table><caption>**Table A.1**
Energy coefficients and input parameters in the phase-field simulation.</caption><thead><tr><th>Landau coefficients</th><th></th><th>Elastic/Electrostrictive</th><th></th><th>Simulation parameters</th><th></th></tr></thead><tbody><tr><td>α<sub>1</sub> (C<sup>-2</sup> m<sup>2</sup> N)</td><td>-1.7 × 10<sup>8</sup></td><td>C<sub>11</sub> (GPa)</td><td>174.6</td><td>P<sub>c0</sub> (Cm<sup>-2</sup>)</td><td>0.757</td></tr><tr><td>α<sub>11</sub> (C<sup>-4</sup> m<sup>6</sup> N)</td><td>-7.3 × 10<sup>7</sup></td><td>C<sub>12</sub> (GPa)</td><td>79.4</td><td>E<sub>0</sub> (V m<sup>-1</sup>)</td><td>1.3 × 10<sup>8</sup></td></tr><tr><td>α<sub>12</sub> (C<sup>-4</sup> m<sup>6</sup> N)</td><td>+7.5 × 10<sup>8</sup></td><td>C<sub>44</sub> (GPa)</td><td>1111.1</td><td>ε<sub>r</sub> (a.u.)</td><td>50</td></tr><tr><td>α<sub>111</sub> (C<sup>-6</sup> m<sup>10</sup> N)</td><td>+2.6 × 10<sup>8</sup></td><td>Q<sub>11</sub> (C<sup>-2</sup> m<sup>4</sup>)</td><td>0.089</td><td>T (°C)</td><td>25</td></tr><tr><td>α<sub>112</sub> (C<sup>-6</sup> m<sup>10</sup> N)</td><td>+6.1 × 10<sup>8</sup></td><td>Q<sub>12</sub> (C<sup>-2</sup> m<sup>4</sup>)</td><td>-0.026</td><td>N<sub>x</sub> (a.u.)</td><td>16</td></tr><tr><td>α<sub>123</sub> (C<sup>-6</sup> m<sup>10</sup> N)</td><td>-3.7 × 10<sup>9</sup></td><td>Q<sub>44</sub> (C<sup>-2</sup> m<sup>4</sup>)</td><td>0.034</td><td>N<sub>y</sub>, N<sub>z</sub> (a.u.)</td><td>32–256</td></tr></tbody></table>

<table><caption>**Table A.2**
Normalization of the parameters in the phase-field simulation.</caption><thead><tr><th>Landau coefficients</th><th></th><th>Elastic/Electrostrictive</th><th></th><th>Simulation parameters</th><th></th></tr></thead><tbody><tr><td>α<sup>*</sup><sub>1</sub> = α<sub>1</sub>/|α<sub>1</sub>|</td><td>-1</td><td>C<sup>*</sup><sub>11</sub> = C<sub>11</sub>/(|α<sub>1</sub>|P<sup>2</sup><sub>c0</sub>)</td><td>1766</td><td>P<sup>*</sup></td><td>P<sub>c0</sub></td></tr><tr><td>α<sup>*</sup><sub>11</sub> = α<sub>11</sub>P<sup>2</sup><sub>c0</sub>/|α<sub>1</sub>|</td><td>-0.24</td><td>C<sup>*</sup><sub>12</sub> = C<sub>12</sub>/(|α<sub>1</sub>|P<sup>2</sup><sub>c0</sub>)</td><td>802</td><td>G<sup>*</sup><sub>11</sub> = G<sub>11</sub>/G<sub>110</sub></td><td>0.6</td></tr><tr><td>α<sup>*</sup><sub>12</sub> = α<sub>12</sub>P<sup>2</sup><sub>c0</sub>/|α<sub>1</sub>|</td><td>2.49</td><td>C<sup>*</sup><sub>44</sub> = C<sub>44</sub>/(|α<sub>1</sub>|P<sup>2</sup><sub>c0</sub>)</td><td>1123</td><td>G<sup>*</sup><sub>44</sub> = G<sub>44</sub>/G<sub>110</sub></td><td>0.3</td></tr><tr><td>α<sup>*</sup><sub>111</sub> = α<sub>111</sub>P<sup>4</sup><sub>c0</sub>/|α<sub>1</sub>|</td><td>0.49</td><td>Q<sup>*</sup><sub>11</sub> = Q<sub>11</sub>P<sup>2</sup><sub>c0</sub></td><td>0.05</td><td>r<sup>*</sup></td><td>r/l<sub>0</sub></td></tr><tr><td>α<sup>*</sup><sub>112</sub> = α<sub>112</sub>P<sup>4</sup><sub>c0</sub>/|α<sub>1</sub>|</td><td>1.16</td><td>Q<sup>*</sup><sub>12</sub> = Q<sub>12</sub>P<sup>2</sup><sub>c0</sub></td><td>-0.014</td><td>t<sup>*</sup></td><td>t/l<sub>0</sub></td></tr><tr><td>α<sup>*</sup><sub>123</sub> = α<sub>123</sub>P<sup>4</sup><sub>c0</sub>/|α<sub>1</sub>|</td><td>-7.04</td><td>Q<sup>*</sup><sub>44</sub> = Q<sub>44</sub>P<sup>2</sup><sub>c0</sub></td><td>0.019</td><td>E<sup>*</sup></td><td>E/E<sub>0</sub></td></tr></tbody></table>

<table><caption>**Table B.1**
Details of the surrogate model architecture for a $N_x \times N_y \times N_z$ system size.</caption><thead><tr><th colspan="3">Encoder</th><th colspan="3">Decoder</th></tr><tr><th>Layer</th><th>Filters</th><th>Shape</th><th>Layer</th><th>Filters</th><th>Shape</th></tr></thead><tbody><tr><td>Input</td><td>None</td><td>[N<sub>x</sub>, N<sub>y</sub>, N<sub>z</sub>, 1]</td><td></td><td></td><td></td></tr><tr><td>Down 1<sup>a</sup></td><td>32</td><td>[N<sub>x</sub>/2, N<sub>y</sub>/2, N<sub>z</sub>/2, 32]</td><td>Up 1<sup>c</sup></td><td>256</td><td>[N<sub>x</sub>/16, N<sub>y</sub>/16, N<sub>z</sub>/16, 256 + 256<sup>e</sup>]</td></tr><tr><td>Down 2<sup>a</sup></td><td>64</td><td>[N<sub>x</sub>/4, N<sub>y</sub>/4, N<sub>z</sub>/4, 64]</td><td>Up 2<sup>d</sup></td><td>128</td><td>[N<sub>x</sub>/8, N<sub>y</sub>/8, N<sub>z</sub>/8, 128 + 128<sup>e</sup>]</td></tr><tr><td>Down 3<sup>a</sup></td><td>128</td><td>[N<sub>x</sub>/8, N<sub>y</sub>/8, N<sub>z</sub>/8, 128]</td><td>Up 3<sup>d</sup></td><td>64</td><td>[N<sub>x</sub>/4, N<sub>y</sub>/4, N<sub>z</sub>/4, 64 + 64<sup>e</sup>]</td></tr><tr><td>Down 4<sup>a</sup></td><td>256</td><td>[N<sub>x</sub>/16, N<sub>y</sub>/16, N<sub>z</sub>/16, 256]</td><td>Up 4<sup>d</sup></td><td>32</td><td>[N<sub>x</sub>/2, N<sub>y</sub>/2, N<sub>z</sub>/2, 32 + 32<sup>e</sup>]</td></tr><tr><td>Down 5<sup>b</sup></td><td>256</td><td>[N<sub>x</sub>/16, N<sub>y</sub>/32, N<sub>z</sub>/32, 256]</td><td>Output</td><td>N<sub>fields</sub></td><td>[N<sub>x</sub>, N<sub>y</sub>, N<sub>z</sub>, 1]</td></tr></tbody></table>

3D convolutional layer (Conv3D) with a 4 × 4 kernel, padding equal to same, and a stride of 2. These blocks successively reduce the image size by a factor of 2 along the (x, y, z) axes and double the number of features. Batch Normalization (BN) and Leaky Rectified Linear Units (LeakyReLU) with a slope of $\alpha = 0.2$ follow each convolutional layer in every downsampling block, except for the first block. Specifics for each layer, including the number of filters, are provided in Table B.1.

The last downsampling block (Down 5) involves a Conv3D layer with a 4 × 4 kernel, padding equal to same, and a stride of 1. Subsequently, a BN layer and LeakyReLU activation are applied, followed by a MaxPooling3D layer with same padding and strides equal to 2. This 3D pooling layer employs a pool size of (1, 2, 2), halving the image size along the (y, z) axes.

The decoder operates in reverse, employing four upsampling blocks and a final output layer. The first upsampling block (Up 1) features a 3D convolutional transpose layer (Conv3DTranspose) with a 4 × 4 kernel, padding equal to same, and a stride of 1. This is succeeded by a BN layer and a Dropout layer with a rate of 0.5 to mitigate overfitting, followed by a Rectified Linear Units (ReLU) activation layer. Subsequently, an UpSampling3D layer with a pool size of (1, 2, 2) is applied, doubling the image size along the (y, z) axes. The subsequent upsampling blocks (Up 2–4) each employ a Conv3DTranspose layer with a 4 × 4 kernel, same padding, and a stride of 2, followed by BN layers and ReLU activation. Each block upsamples the image size by a factor of 2 along the (x, y, z) axes.

A final Conv3DTranspose layer with a 4 × 4 kernel, same padding, and a stride of 2 generates the full-size $N_x \times N_y \times N_z$ network predictions with hyperbolic tangent activation. Skip connections, link encoder and decoder layers to improve the coherence of decoder reconstructions with the input image [52].

The choice of this architecture is driven by the dimensions of the selected microstructure input. Specifically, for a size of $N_x \times N_y \times N_z = 16 \times 32 \times 32$, the final encoding compresses the vector along all (x, y, z) axes to a fully compressed shape of [1, 1, 1, 256] in the latent space. Furthermore, this architectural framework enables its straightforward application to the larger systems considered in this study.

<footnote>a Downsampling block.</footnote>

<footnote>b Downsampling Pool-Based block.</footnote>

<footnote>c Upsampling Pool-Based block.</footnote>

<footnote>d Upsampling block.</footnote>

<footnote>e Skip-connection features.</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. D.1. Representation of the cumulative explained variance evolution obtained through PCA analysis on the training dataset. The solid line depicts the cumulative explained variance for the $P_x$ field, while the dotted line represents the threshold corresponding to 95% of the cumulative explained variance.</figcaption></figure>

## Appendix C. Error metrics

The mean squared error (MSE) loss function $\mathcal{L}_2$ can be computed over the $\{(Y_i, X_i)\}_{i=1}^N$ training samples as :

$$ \mathcal{L}_2(Y^{t_{k+1}}, S(X^{t_k})) = \frac{1}{N} \sum_{i=1}^{N} (Y_i^{t_{k+1}} - S(X_i^{t_k}))^2 \quad (C.1) $$

where $Y^{t_{k+1}}$ are the microstructure labels obtained by real phase-field simulations and $S(X^t)$ are the model outputs.

To assess rollout simulations, the score is averaged across each trajectory, such as:

$$ \mathcal{L}_2 = \frac{1}{MN} \sum_{m=1}^{M} \sum_{n=1}^{N} (Y_m^{t_{n+1}} - S_m(X^{t_n}))^2 \quad (C.2) $$

where $M$ denotes the number of trajectories in the validation dataset and $N$ is the number of frames per simulation.

The Mean Absolute Error (MAE) can be determined by:

$$ \text{MAE} = \frac{1}{N} \sum_{n=1}^{N} |Y_m^{t_{k+1}} - S_m(X^{t_k})| \quad (C.3) $$

Finally, the Macro Average Relative Error (MARE) can be computed in the context of rollout evaluations by:

$$ \text{MARE} = \frac{1}{M} \sum_{m=1}^{M} \frac{\sum_{n=1}^{N} |Y_m^{t_{n+1}} - S_m(X^{t_n})|}{\sum_{n=1}^{N} |Y_m^{t_{n+1}}|} \quad (C.4) $$

## Appendix D. Principal component analysis

Principal Component Analysis (PCA), following methodologies established in previous Refs. [30,35], was executed to represent our dataset in a condensed, low-dimensional space. PCA, as a dimensionality reduction technique, transforms the dataset into a reduced basis that optimally captures the majority of data variance. The orthonormal vectors of this new basis, denoted as $\{\phi_i\}_{i=1}^N$, are referred to as principal components (PC), designed to maximize variance while minimizing the number of dimensions.

In the resultant low-dimensional space, the diverse trajectories originating from initial states underscore the variety present in the simulated datasets. The nonlinear behavior inherent in TDLG equations contributes to distinct domain formation trajectories based on varying initial conditions. Fig. D.1 elucidates the evolution of cumulative variance as a function of the number of principal components dedicated to the $P_x$ field. In particular, the analysis emphasizes the requirement for a significant number of principal components to capture 95% of the cumulative explained variance. This underscores the intricate and diverse nature of simulation trajectories encapsulated in the dataset.

## Appendix E. Training details

To improve stability in rollout predictions, we introduced Gaussian noise into the input training features. Employing a progressive random walk scheme, the noise was incrementally augmented along simulation trajectories to emulate typical error accumulation data [29,61]. Importantly, no noise was introduced into the label targets obtained from high-fidelity phase-field simulations. The magnitudes of noise for the polarization field was set at $\sigma_p = 10^{-3}$, aligning with the methodology outlined in Ref. [29].
The model parameters were optimized during training using the Adam optimizer [62] with a batch size of 32. For the training on the $N_x \times N_y \times N_z = 16 \times 32 \times 32$ structures, the initial learning rate was set at $10^{-3}$ and gradually reduced to $10^{-6}$, following an exponential decay over during 100 epochs. For the fine-tuning of the network parameters to larger structures of $N_x \times N_y \times N_z = 16 \times 64 \times 64$ to $N_x \times N_y \times N_z = 16 \times 256 \times 256$ grid dimensions, the model was trained during 20 epochs with an initial learning rate of $10^{-5}$ and gradually reduced to $10^{-6}$, following an exponential decay.

## References

[1] J.F. Scott, C.A.P. de Araujo, Ferroelectric memories, Science 246 (4936) (1989) 1400-1405.

[2] A. Crassous, S. Sluka, A.K. Tagantsev, N. Setter, Polarization charge as a reconfigurable quasi-dopant in ferroelectric thin films, Nature Nanotechnol. 10 (7) (2015) 614-618.

[3] P. Sharma, D. Sando, Q. Zhang, X. Cheng, S. Prosandeev, R. Bulanadi, S. Prokhorenko, L. Bellaiche, L.-Q. Chen, V. Nagarajan, J. Seidel, Conformational domain wall switch, Adv. Funct. Mater. 29 (18) (2019) 1807523.

[4] A.I. Khan, K. Chatterjee, B. Wang, S. Drapcho, L. You, C. Serrao, S.R. Bakaul, R. Ramesh, S. Salahuddin, Negative capacitance in a ferroelectric capacitor, Nat. Mater. 14 (2) (2015) 182-186.

[5] T. Mikolajick, S. Slesazeck, H. Mulaosmanovic, M.H. Park, S. Fichtner, P.D. Lomenzo, M. Hoffmann, U. Schroeder, Next generation ferroelectric materials for semiconductor process integration and their applications, J. Appl. Phys. 129 (10) (2021) 100901.

[6] M.H. Park, Y.H. Lee, H.J. Kim, Y.J. Kim, T. Moon, K.D. Kim, J. Müller, A. Kersch, U. Schroeder, T. Mikolajick, C.S. Hwang, Ferroelectricity and antiferroelectricity of doped thin $HfO_2$-based films, Adv. Mater. 27 (11) (2015) 1811-1831.

[7] R. Khosla, S.K. Sharma, Integration of ferroelectric materials: An ultimate solution for next-generation computing and storage devices, ACS Appl. Electron. Mater. 3 (7) (2021) 2862-2897.

[8] N. Setter, D. Damjanovic, L. Eng, G. Fox, S. Gevorgian, S. Hong, A. Kingon, H. Kohlstedt, N.Y. Park, G.B. Stephenson, I. Stolitchnov, A.K. Taganstev, D.V. Taylor, T. Yamada, S. Streiffer, Ferroelectric thin films: Review of materials, properties, and applications, J. Appl. Phys. 100 (5) (2006) 051606.

[9] N. Izyumskaya, Y.-I. Alivov, S.-J. Cho, H. Morkoç, H. Lee, Y.-S. Kang, Processing, structure, properties, and applications of PZT thin films, Crit. Rev. Solid State Mater. Sci. 32 (3-4) (2007) 111-202.

[10] L. McGilly, P. Yudin, L. Feigl, A. Tagantsev, N. Setter, Controlling domain wall motion in ferroelectric thin films, Nature Nanotechnol. 10 (2) (2015) 145-150.

[11] N. Setter, D. Damjanovic, L. Eng, G. Fox, S. Gevorgian, S. Hong, A. Kingon, H. Kohlstedt, N.Y. Park, G.B. Stephenson, I. Stolitchnov, A.K. Taganstev, D.V. Taylor, T. Yamada, S. Streiffer, Ferroelectric thin films: Review of materials, properties, and applications, J. Appl. Phys. 100 (5) (2006) 051606.

[12] S. Gonzalez Casal, X. Bai, K. Alhada-Lahbabi, B. Canut, B. Vilquin, P. Rojo Romeo, S. Brottet, D. Albertini, D. Deleruyelle, M. Bugnet, I. Canero Infante, B. Gautier, Mechanical switching of ferroelectric domains in 33-200 nm-thick sol-gel-grown $PbZr_{0.2}Ti_{0.8}O_3$ films assisted by nanocavities, Adv. Electron. Mater. 8 (9) (2022) 1-9.

[13] E.J. Guo, R. Roth, S. Das, K. Dörr, Strain induced low mechanical switching force in ultrathin $PbZr_{0.2}Ti_{0.8}O_3$ films, Appl. Phys. Lett. 105 (1) (2014) 012903.

[14] Y.L. Li, S.Y. Hu, Z.K. Liu, L.Q. Chen, Effect of substrate constraint on the stability and evolution of ferroelectric domain structures in thin films, Acta Mater. 50 (2) (2002) 395-411.

[15] L.-Q. Chen, Phase-field models for microstructure evolution, Annu. Rev. Mater. Res. 32 (1) (2002) 113-140.

[16] J.J. Wang, B. Wang, L.Q. Chen, Understanding, predicting, and designing ferroelectric domain structures and switching guided by the phase-field method, Annu. Rev. Mater. Res. 49 (2019) 127-152.

[17] Y. Zhao, Understanding and design of metallic alloys guided by phase-field simulations, njp Comput. Mater. 9 (2023) 94.

[18] O.C. Zienkiewicz, R.L. Taylor, R.L. Taylor, The finite element method: solid mechanics, vol. 2, Butterworth-heinemann, 2000.

[19] J.A. Cottrell, T.J. Hughes, Y. Bazilev, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009.

[20] P.V. Balachandran, B. Kowalski, A. Sehirlioglu, T. Lookman, Experimental search for high-temperature ferroelectric perovskites guided by two-step machine learning, Nature Commun. 9 (1) (2018) 1668.

[21] P. Xu, D. Chang, T. Lu, L. Li, M. Li, W. Lu, Search for ABO3 type ferroelectric perovskites with targeted multi-properties by machine learning strategies, J. Chem. Inf. Model. 62 (21) (2022) 5038-5049.

[22] X.-Y. Ma, H.-Y. Lyu, K.-R. Hao, Y.-M. Zhao, X. Qian, Q.-B. Yan, G. Su, Large family of two-dimensional ferroelectric metals discovered via machine learning, Sci. Bull. 66 (3) (2021) 233-242.

[23] Y. Liu, K.P. Kelley, R.K. Vasudevan, H. Funakubo, M.A. Ziatdinov, S.V. Kalinin, Experimental discovery of structure-property relationships in ferroelectric materials via active learning, Nat. Mach. Intell. 4 (4) (2022) 341-350.

[24] M. Ziatdinov, C. Nelson, R.K. Vasudevan, D.Y. Chen, S.V. Kalinin, Building ferroelectric from the bottom up: The machine learning analysis of the atomic-scale ferroelectric distortions, Appl. Phys. Lett. 115 (5) (2019) 052902.

[25] S.M. Neumayer, S. Jesse, G. Velarde, A.L. Kholkin, I. Kravchenko, L.W. Martin, N. Balke, P. Maksymovych, To switch or not to switch – a machine learning approach for ferroelectricity, Nanoscale Adv. 2 (2020) 2063-2072.

[26] Y. Liu, B. Yu, Z. Liu, D. Beck, K. Zeng, High-speed piezoresponse force microscopy and machine learning approaches for dynamic domain growth in ferroelectric materials, ACS Appl. Mater. Interfaces 12 (8) (2020) 9944-9952.

[27] M.M. Rahman, S. Janwari, M. Choi, U.V. Waghmare, J. Lee, Accelerating search for the polar phase stability of ferroelectric oxide by machine learning, Mater. Des. 236 (2021) 112518.

[28] T.M. Linker, K.-i. Nomura, S. Fukushima, R.K. Kalia, A. Krishnamoorthy, A. Nakano, K. Shimamura, F. Shimojo, P. Vashishta, Induction and ferroelectric switching of flux closure domains in strained $PbTiO_3$ with neural network quantum molecular dynamics, Nano Lett. 23 (16) (2023) 7456-7462.

[29] K. Alhada-Lahbabi, D. Deleruyelle, B. Gautier, Machine learning surrogate model for acceleration of ferroelectric phase-field modeling, ACS Applied Electronic Materials 5 (7) (2023) 3894-3907.

[30] D. Montes de Oca Zapiain, J.A. Stewart, R. Dingreville, Accelerating phase-field-based microstructure evolution predictions via surrogate models trained by machine learning methods, njp Comput. Mater. 7 (1) (2021) 1-11.

[31] W. Li, M.Z. Bazant, J. Zhu, Phase-field DeepONet: Physics-informed deep operator neural network for fast simulations of pattern formation governed by gradient flows of free-energy functionals, Comput. Methods Appl. Mech. Engrg. 416 (2023) 116299.

[32] L. Lu, P. Jin, G. Pang, Z. Zhang, G.E. Karniadakis, Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators, Nat. Mach. Intell. 3 (3) (2021) 218-229.

[33] M. Raissi, P. Perdikaris, G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686-707.

[34] A. Hemmasian, F. Ogoke, P. Akbari, J. Malen, J. Beuth, A. Barati Farimani, Surrogate modeling of melt pool temperature field using deep learning, Addit. Manuf. Lett. 5 (2023) 100123.

[35] S. Fetni, T.Q.D. Pham, T.V. Hoang, H.S. Tran, L. Duchêne, X.-V. Tran, A.M. Habraken, Capabilities of auto-encoders and principal component analysis of the reduction of microstructural images; application on the acceleration of phase-field simulations, Comput. Mater. Sci. 216 (2023) 111820.
[36] J.Y. Choi, T. Xue, S. Liao, J. Cao, Accelerating phase-field simulation of three-dimensional microstructure evolution in laser powder bed fusion with composable machine learning predictions, Addit. Manuf. 79 (2024) 103938.

[37] V. Oommen, K. Shukla, S. Desai, R. Dingreville, G.E. Karniadakis, Rethinking materials simulations: Blending direct numerical simulations with neural operators, 2023.

[38] I. Peivaste, N.H. Siboni, G. Alahyarizadeh, R. Ghaderi, B. Svendsen, D. Raabe, J.R. Mianroodi, Machine-learning-based surrogate modeling of microstructure evolution using phase-field, Comput. Mater. Sci. 214 (2022) 111750.

[39] T. Xue, Z. Gan, S. Liao, J. Cao, Physics-embedded graph network for accelerating phase-field simulation of microstructure evolution in additive manufacturing, npj Comput. Mater. 8 (1) (2022) 201.

[40] V. Oommen, K. Shukla, S. Goswami, R. Dingreville, G.E. Karniadakis, Learning two-phase microstructure evolution using neural operators and autoencoder architectures, npj Comput. Mater. 8 (1) (2022) 190.

[41] K. Yang, Y. Cao, Self-supervised learning and prediction of microstructure evolution with convolutional recurrent neural networks, Patterns 2 (5) (2021) 100243.

[42] P. Wu, A.S. Iquebal, K. Ankit, Emulating microstructural evolution during spinodal decomposition using a tensor decomposed convolutional and recurrent neural network, Comput. Mater. Sci. 224 (2023) 112187.

[43] F.P. Kemeth, S. Alonso, B. Echebarria, T. Moldenhawer, C. Beta, I.G. Kevrekidis, Black and gray box learning of amplitude equations: Application to phase field systems, Phys. Rev. E 107 (2023) 025305.

[44] G.E. Karniadakis, I.G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nat. Rev. Phys. 3 (6) (2021) 422–440.

[45] E. Samaniego, C. Anitescu, S. Goswami, V.M. Nguyen-Thanh, H. Guo, K. Hamdia, X. Zhuang, T. Rabczuk, An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications, Comput. Methods Appl. Mech. Engrg. 362 (2020) 112790.

[46] G.H. Teichert, A.R. Natarajan, A. Van der Ven, K. Garikipati, Machine learning materials physics: Integrable deep neural networks enable scale bridging by learning free energy functions, Comput. Methods Appl. Mech. Engrg. 353 (2019) 201–216.

[47] T. Chen, H. Chen, Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems, IEEE Trans. Neural Netw. 6 (4) (1995) 911–917.

[48] S. Goswami, M. Yin, Y. Yu, G.E. Karniadakis, A physics-informed variational DeepONet for predicting crack path in quasi-brittle materials, Comput. Methods Appl. Mech. Engrg. 391 (2022) 114587.

[49] Y. Cao, A. Morozovska, S.V. Kalinin, Pressure-induced switching in ferroelectrics: Phase-field modeling, electrochemistry, flexoelectric effect, and bulk vacancy dynamics, Phys. Rev. B 96 (18) (2017) 184109.

[50] Y.-J. Wang, J. Li, Y.-L. Zhu, X.-L. Ma, Phase-field modeling and electronic structural analysis of flexoelectric effect at 180% domain walls in ferroelectric PbTiO3, J. Appl. Phys. 122 (22) (2017) 224101.

[51] Y. Cao, S.V. Kalinin, Phase-field modeling of chemical control of polarization stability and switching dynamics in ferroelectric thin films, Phys. Rev. B 94 (2016) 235444.

[52] O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation, in: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18, Springer, 2015, pp. 234–241.

[53] L. Jiang, Z. Zhang, H. Hu, X. He, H. Fu, J. Xie, A rapid and effective method for alloy materials design via sample data transfer machine learning, npj Comput. Mater. 9 (1) (2023) 26.

[54] Z. Wang, H. Zhang, J. Ren, X. Lin, T. Han, J. Liu, J. Li, Predicting adsorption ability of adsorbents at arbitrary sites for pollutants using deep transfer learning, npj Comput. Mater. 7 (1) (2021) 19.

[55] Y. Kim, Y. Kim, C. Yang, K. Park, G.X. Gu, S. Ryu, Deep learning framework for material design space exploration using active transfer learning and data augmentation, npj Comput. Mater. 7 (2021) 140.

[56] K. Alhada-Lahbabi, D. Deleruyelle, B. Gautier, Phase-field study of nanocavity-assisted mechanical switching in pbti03 thin films, advanced electronic materials 2300744.

[57] H.-L. Hu, L.-Q. Chen, Three-dimensional computer simulation of ferroelectric domain formation, J. Am. Ceram. Soc. 81 (3) (1998) 492–500.

[58] J.J. Wang, X.Q. Ma, Q. Li, J. Britson, L.Q. Chen, Phase transitions and domain structures of ferroelectric nanoparticles: Phase field model incorporating strong elastic and dielectric inhomogeneity, Acta Mater. 61 (20) (2013) 7591–7603.

[59] Y.L. Li, S.Y. Hu, Z.K. Liu, L.Q. Chen, Effect of substrate constraint on the stability and evolution of ferroelectric domain structures in thin films, Acta Mater. 50 (2) (2002) 395–411.

[60] J. Wang, S.Q. Shi, L.Q. Chen, Y. Li, T.Y. Zhang, Phase-field simulations of ferroelectric/ferroelastic polarization switching, Acta Mater. 52 (3) (2004) 749–764.

[61] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, P. Battaglia, Learning to simulate complex physics with graph networks, in: International Conference on Machine Learning, PMLR, 2020, pp. 8459–8468.

[62] D.P. Kingma, J.L. Ba, Adam: A method for stochastic optimization, in: 3rd Int. Conf. Learn. Represent. ICLR 2015 - Conf. Track Proc., 2015, pp. 1–15.