**I. Physical modeling of low dimensional ferroelectrics**

Determine your model system. You should determine the dimensionality of the system (i.e., 0-D nanoplatelet 1-D nanowire and 2-D nanofilm), material, and electric and mechanical boundary conditions, which are recommended to be written as optional tags in your program. In this demo, we consider a PbTiO3 ferroelectric nanoplatelet under surface traction and open-circuit condition as shown in **Figure 1**.

Construct the system’s free energy. In general, the free energy of a ferroelectric is a functional of polarization, and should incorporate the effects of polarization inhomogeneity, elastic and electric fields. It can be written as a sum of the bulk free energy (including a Landau-type potential, gradient energy, elastic energy and electric energy) and surface free energy.

**Tip**: The construction of free energy is very important in phase-field model of ferroelectrics. The governing equation, i.e., the Time Dependent Ginzburg–Landau (TDGL) equation, describes that the evolution of domain structure toward its equilibrium is driven by the decrease of the system’s free energy.

**Caution**: The free energy should be generally constructed under the thermodynamic framework.

- Design your experiment. In our case, the experiment is to apply mechanical loads to the nanoplatelet to affect its domain structure.

**II. Programming**

The flow diagram of phase-field simulation at a given condition (e.g., fixed temperature, open-circuited and mechanical load) is as shown in **Figure 2**. Followings are the main functional parts of the corresponding program.

- Initialization. In this part of program, the parameters and variables needed for simulation would be defined and initialized. It is recommended to write this part in a way that it can either use the default values of the parameters and variables or read them from input files.

**Tip**: For numerical accuracy and convenience, it is recommended to make the parameters and variables dimensionless so that their values fall into a suitable range.

- Grid division of the system. This process would generate some tables, from which you can easily find out the information of the nodes and elements, such as its coordinates, adjacent nodes and elements. In our case, we make the following tables,
Tnode ----Each row stores a node’s number, coordinates and adjacent elements,

Tnodes---- Each row stores a node’s number and its adjacent nodes,

Telement---- Each row stores an element’s number, coordinates and its nodes.

**Tip**: These tables can largely simplify the calculation of element stiffness matrices, element node displacement/potential vectors, and the assemblage of global stiffness matrices and displacement/potential vectors.

- Calculation of element stiffness matrices and the assemblage of global stiffness matrices. In most cases, the stiffness matrices (including elastic and electric) can be considered unchanged during the process, therefore the calculation only need to be done by once.

**Tip**: According to the boundary conditions, the computation size and accuracy requirement, the electric and elastic fields can be solved by different methods, such as conjugate gradient method (CGM), finite-element method (FEM), and fast Fourier transformation (FFT). For our case of a ferroelectric nanoplatelet with a moderate size, finite-element method is suitable to solve the electric and elastic fields.

Calculation of the elastic field. Calculate the element node displacement vectors, and assemble them into the global node displacement vector. Solve the equation [Ku]{U}={Fu}. This should be done at each step of polarization evolution.

Calculation of the electric field. Calculate the element node potential vectors, and assemble them into the global node potential vector. Solve the equation [Kφ]{Ф}={Fφ}. This should be done at each step of polarization evolution.

**Tip**: Iterative methods, e.g. the Gauss-Seidel iteration method, would be suitable to equations [Ku]{U}={Fu} and [Kφ]{Ф}={Fφ}.

Calculation of the evolution force and polarization field at next step. Simple explicit difference methods or semi-implicit Fourier-spectral algorithms8 can be used to solve the governing equation, i.e., the TDGL equation.

Error analysis. Calculate the error between the polarization field at this step and next step. If the error is small enough, end program; otherwise, repeat steps 3-6.

Input and output. In our case, the input includes a parameter file storing the value of parameters and a configuration file storing the initial polarization field. The output includes a log file recording the monitoring information during the simulation and files storing the polarization field at selected steps.

**III. Simulation**

We would like to simulate how the mechanical load affects the formation of domain structure in the nanoplatelet.

Prepare a set of parameter files. In our case, the parameter files differ only in the value of surface traction.

For each parameter files, compile and run the program.

**Tip**: You can also make some modifications on the original program, so that the whole simulation can be fulfilled by running the program by once.

**IV. Result analysis & Post processing**

At this stage, we assume that there is an output file storing the polarization field, e.g., p_final.txt. The file stores the node’s number, its coordinates and the polarization components. Now we will show how to use MATLAB (http://www.mathworks.com) to visualize the domain structures.

Start MATLAB.

Import file p_final.txt into the Workspace as shown in Figure 3.

Run the following commands to obtain a vector plot the domain structure as shown in **Figure 4a**,

quiver3(p_final(:,3),p_final(:,2),p_final(:,4),p_final(:,6),p_final(:,5),p_final(:,7));

axis equal;

axis off;

- Run the following commands to add a color plot of the polarization magnitude to
**Figure 4a** as shown in **Figure 4b**,

meshing=[10 3 30];

nx=meshing(1)+1;

ny=meshing(2)+1;

nz=meshing(3)+1;

m=0;

for i=1:nx

```
for j=1:ny
for k=1:nz
m=m+1;
P\(i,j,k)=\(p_final\(m,5)**p_final\(m,5)+p_final\(m,6)**p_final\(m,6)+p_final\(m,7)*p_final\(m,7))^0.5;
end
end
```

end

[x,y,z]=meshgrid([0:ny-1],[0:nx-1],[0:nz-1]);

Xslice=[0,ny-1];Yslice=[0,nx-1];Zslice=[0,nz-1];

hold on;h =slice(x,y,z,P,Xslice,Yslice,Zslice);

set(h,'LineStyle','none');

- Change the colormap to your preference. For example, run the following commands to as shown in
**Figure 4c**,

color_round=[1 0 1];color_floor=[0 1 0];

n_color=64;

for i=1:n_color

```
map\(i,1)=color_round\(1)+\(i-1)*\(color_floor\(1)-color_round\(1))/\(n_color-1);
map\(i,2)=color_round\(2)+\(i-1)*\(color_floor\(2)-color_round\(2))/\(n_color-1);
map\(i,3)=color_round\(3)+\(i-1)*\(color_floor\(3)-color_round\(3))/\(n_color-1);
```

end

colormap(map);