# Solving 2D Landau-Khalatnikov equation & Poisson equation using finite element approach

I would like to solve Landau-Khalatnikov (LK) equation and Poisson equation using finite element approach to simulate ferroelectric hysteresis loop, electric potential distribution and electric polarization distribution. Here, I will provide detail as much as possible. LK equation, which is also called time-dependent Ginzburg-Landau (TDGL) equation, can be simply expressed as follow:

#1: $$-\frac{1}{\Gamma}\frac{\partial P}{\partial t}= \alpha P + \beta P^{3} + \gamma P^{5}-E$$

Where Γ is kinetic coefficient, α, β, γ are Landau coefficients, P is polarization and E is electric field (of course, t is time). The picture is much like a ferroelectric material sandwiched by electrode, one end has fixed potential and other end is grounded to generate electric field. First, to simply test feasibility of MMA to solve such an equation, I wrote very simple code to solve LK equation in 1D:

Clear[\[Alpha], \[Beta], \[Gamma], Efieldamp, \[Omega]];
\[Alpha] = -1.35*10^9;
\[Beta] = -2.64*10^10;
\[Gamma] = 2.5*10^11;
Efieldamp = 2*10^9;
\[Omega] = 10^6*Pi;
\[CapitalGamma] = 1;

Efield[t_?NumericQ] := Efieldamp*Sin[\[Omega]*t];

SimpleLKSolver =
First@NDSolve[{-1/\[CapitalGamma]*D[P[t], t] == \[Alpha]*(P[t]*0.01) + \[Beta]*(P[t]*0.01)^3 +
\[Gamma]*(P[t]*0.01)^5 - Efield[t], P == 0}, P, {t, 0, 2.5*1/\[Omega]*Pi}, Method -> Automatic, MaxSteps -> Infinity]



As shown in the code, I put E as Efieldamp * Sin[ω * t] where Efieldamp and ω are defined as above. The multiplication of 0.01 to P in rhs is for unit conversion (from C/m2 to uC/cm2). The code works nicely, and plotting P vs. V (V was calculated assuming that the distance between voltage source and ground is 10 nm (10-8 m)) gives nice hysteresis curve:

LK =
ParametricPlot[{Efield[t]*10^-8, P[t] /. SimpleLKSolver}, {t, 0, 2.5*1/\[Omega]*Pi}, PlotRange -> All, AspectRatio -> 2/3,
PlotStyle -> {Bold, Black}, Frame -> True,
FrameStyle -> {{Black, Bold}, {Black, Bold}}] Now, I would like to solve the equation in 2D using finite element method. Actually, I would also like to simultaneously solve Poisson’s equation to see electric potential. In 2D, the equations are:

#2: $$-\frac{1}{\Gamma}\frac{\partial P(x,z)}{\partial t}= \alpha P(x,z) + \beta P(x,z)^{3} + \gamma P(x,z)^{5}-g_{11}\frac{\partial^2 P(x,z)}{\partial z^2}-g_{44}\frac{\partial^2 P(x,z)}{\partial x^2}+\frac{\mathrm{d} \phi(x,z)}{\mathrm{d} x}$$

#3: $$-\varepsilon _{0}(\varepsilon _{x}\frac{\partial^2 \phi }{\partial x^2}+\varepsilon _{z}\frac{\partial^2 \phi }{\partial z^2})=-\frac{\mathrm{d} P}{\mathrm{d} z}$$

In 2D space, I defined x as horizontal direction and z as vertical direction. The spatial second derivative terms are now added to LK (1st) equation in order to calculate spatial gradient of polarization; g11 (g44) are gradient coefficient, ϕ is electric potential, εx and εz are anisotropic dielectric permittivity in x and z direction. dϕ(x,z)/dz is electric field term in 2D, considering electric field only in Z-direction. In Poisson equation (2nd equation), -dP(x,z)/dz gives volume charge density considering polarization in ferroelectric materials. So, the workflow is basically this: first solve LK equation to obtain polarization distribution in 2D and obtain PV hysteresis curve, then obtain electric potential distribution by solving Poisson equation with obtained P value. My try is the following:

Needs["NDSolveFEM"]

region = Rectangle[{0, 0}, {100*10^-9, 10*10^-9}];
mesh = ToElementMesh[region, MaxCellMeasure -> 2*10^-18]

Clear[\[Alpha], \[Beta], \[Gamma], g44, g11, \[Epsilon]x, \[Epsilon]z, \[Epsilon]0, \[Omega], \[CapitalGamma]];
\[Alpha] = -1.35*10^9;
\[Beta] = -2.64*10^10;
\[Gamma] = 2.5*10^11;
g44 = 2*10^-10;
g11 = 0.5*10^-10;
\[Epsilon]x = 25;
\[Epsilon]z = 18;
\[Epsilon]0 = 8.854*10^-12;
\[Omega] = 10^6*Pi;
\[CapitalGamma] = 1;
Vapp[t_] := 20*Sin[\[Omega]*t]

PDESolution =
First@NDSolve[{Inactive[Div][{{\[Epsilon]0*\[Epsilon]x, 0}, {0, \[Epsilon]0*\[Epsilon]z}}.Inactive[Grad][\[Phi][t, x, z], {x, z}], {x, z}] ==
D[P[t, x, z], z], -1/\[CapitalGamma]*
D[P[t, x, z], t] == \[Alpha]*P[t, x, z] + \[Beta]*
P[t, x, z]^3 + \[Gamma]*P[t, x, z]^5 -
g11*D[P[t, x, z], {z, 2}] - g44*D[P[t, x, z], {x, 2}] +
D[\[Phi][t, x, z], z],
DirichletCondition[\[Phi][t, x, z] == 0,
z == 0 && 0 < x < 10*10^-8],
P[0, x, z] == RandomVariate[NormalDistribution[]], \[Phi][0, x, z] == 0,
DirichletCondition[\[Phi][t, x, z] == Vapp[t],
z == 10*10^-9 && 0 < x < 10*10^-8],
PeriodicBoundaryCondition[P[t, x, z], x == 0,
TranslationTransform[{10*10^-8, 0}]],
PeriodicBoundaryCondition[\[Phi][t, x, z], x == 0,
TranslationTransform[{10*10^-8, 0}]]}, {\[Phi],
P}, {x, z} \[Element] mesh, {t, 0, 2.5*1/\[Omega]*Pi},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement"}}}]



For boundary condition, I put ϕ = Vapp (defined above) for top surface (z = 10-8) and ϕ = 0 for bottom surface (z = 0). I also applied periodic boundary conditions (in x direction) for both ϕ and P. For initial condition, I put ϕ(t =0) = 0, whereas I put P= RamdomVariate[NormalDistribution[]] to assign random initial condition spatially (but sum is zero!) This code gives me the following error:

NDSolve::ivcon: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values.

LinearSolve::sing: Matrix <<1>> is singular.

LinearSolve::sing: Matrix SparseArray[<<1>>] is singular.

General::stop: Further output of LinearSolve::sing will be suppressed during this calculation.



With the given error message, I know that something is wrong with my initial condition of P, which is P(0,x,z)= RandomVariate[NormalDistribution[]]; therefore, I put 0 in rhs instead. Then, I obtained repeated error messages of the following lines:

LinearSolve::sing: Matrix SparseArray[Automatic,<<1>>,<<3>>,{1,{{<<1>>},{<<1>>}},{4.*10^9,-4.44089*10^-7,-4.*10^9,1.77636*10^-6,4.*10^9,-4.*10^9,<<51862>>,2.72987*10^-10,-5.8139*10^-10,-7.82473*10^-10,-1.44678*10^-9,-7.82473*10^-10,4.*10^9}}] is singular.



So, I realized that some singularity issue exists in the solution. To remove complication, I also tried removing gradient terms (second-order spatial derivatives), but similar error message came out (singularity error). Having confirmed that 1D solution can easily be obtained using NDSolve command, I assume that 2D solution also exists and that something is wrong with my code for solving 2D problem.

## More explanation on my question

Basically, I'm trying to reproduce the results from the following paper published in arXiv: https://arxiv.org/ftp/arxiv/papers/2105/2105.04647.pdf

The major difference is that the authors here utilized "finite difference method", whereas I'm trying to utilize finite element method. Actually, the equations should be solved in "self-consistent" manner, meaning that, for each time step, polarization calculated from equation #2 (Now I labelled the equations!) are put into equation #3 to calculate electric potential, which is again put into equation #2; this iteration stops until the solutions converge. Though equation #3 is not time-dependent, it is solved statically for each P calculated from each time step. I believe how it is solved is analogous to the case of coupled Poisson-schrodinger equation, where charge distribution is calculated from the wave functions calculated using Schrodinger's equation. and this calculated charge distribution is put into Poisson equation to calculate the electric potential, and so on.

To build the code, I have read some tutorial articles in mathematica website (such as components and data structures) but I have no idea. Any hints or helps will be deeply appreciated!

• I think is going to be tricky unless you can add a time dependent derivative for phi. Apr 13, 2022 at 8:14
• I have deleted my comment/answer such that this question shows as unanswered. Apr 13, 2022 at 15:14
• Could you give a link with your model description? Do you have some data to test code? Apr 14, 2022 at 7:29
• @YSP In your model electric potential equation (Poisson's equation) depends on polarization as D[P[t, x, z], z] , while in the original model there is external electric field only. This term makes model unstable so that electric potential increases up to $Vapp/(\epsilon_0 \epsilon_z)$. Apr 15, 2022 at 2:21
• @YSP In the paper they simulate dielectric layer compose be Al203. I add this layer in my answer. Apr 19, 2022 at 5:15

We can compute this problem by combination static FEM and dynamic method of lines as follows (we use all parameters from the paper "Multi-domain Polarization Switching in Hf0.5Zr0.5O2-Dielectric Stack: The Role of Dielectric Thickness" for dielectric layer of 1 nm)

Needs["NDSolveFEM"];
L = 3 10^-8;(*Al2O3 - 1nm, HZO - 10nm*)dL = 1 10^-9; Lz = 10^-8;
region = Rectangle[{0, 0}, {1, (Lz + dL)/L}];
mesh = ToElementMesh[region, MaxCellMeasure -> .002]
region1 = Rectangle[{0, 0}, {1, Lz/L}];
mesh1 = ToElementMesh[region1, MaxCellMeasure -> .002]

Clear[\[Alpha], \[Beta], \[Gamma], g44,
g11, \[Epsilon]x, \[Epsilon]z, \[Epsilon]0, \[CapitalGamma]];
g11d = 0.1*10^-10; g44d = 1.5*10^-10; lambda = 10^-9; omega =
10^6; t0 = 2 Pi/omega;
\[CapitalGamma] = 1/10; \[Alpha] = -1.625*10^9 t0;
\[Beta] = -6.357*10^10 t0;
\[Gamma] = 6.994*10^11 t0;
g44 = t0/L^2 g44d;
g11 = t0/L^2 g11d;
appro = With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
\[Epsilon]x1 = 22; \[Epsilon]z1 = 18; \[Epsilon]x0 = 10; \[Epsilon]z0 \
= 10;
\[Epsilon]x =
SimplifyPWToUnitStep@
PiecewiseExpand@
If[0 <= x <= 1 && 0 <= z <= Lz/L, \[Epsilon]x1, \[Epsilon]x0] /.
UnitStep -> appro; \[Epsilon]z =
SimplifyPWToUnitStep@
PiecewiseExpand@
If[0 <= x <= 1 && 0 <= z <= Lz/L, \[Epsilon]z1, \[Epsilon]z0] /.
UnitStep -> appro;
\[Epsilon]0 = 8.854*10^-12; dt = 1/100;
conu0 = -0.0001; noise = .00001;
Vapp[t_] := 9 Sin[2 Pi t]; n = Length[mesh["Coordinates"]];
eq1 = Inactivate[
Div[{{\[Epsilon]0*\[Epsilon]x/L,
0}, {0, \[Epsilon]0*\[Epsilon]z/L}} .
Grad[\[Phi][x, z], {x, z}], {x, z}] - pz[i - 1][x, z],
D | Div | Grad]; eq2 =
1/\[CapitalGamma]*(D[P[x, z, t], t])/dt + \[Alpha]*
P[x, z, t] + \[Beta]*P[x, z, t]^3 + \[Gamma]*P[x, z, t]^5 -
g11*D[P[x, z, t], {z, 2}] - g44*D[P[x, z, t], {x, 2}] +
t0/L D[phi[i][x, z], z];

T = Table[dt i, {i, 0, 1000}];
p[x_, z_, t_] := 0; pz[x_, z_] := 0;
Do[phi[i] =
NDSolveValue[{Activate[eq1] == 0,
DirichletCondition[{\[Phi][x, z] == 0}, z == 0 && 0 < x < 1],
DirichletCondition[{\[Phi][x, z] == Vapp[dt i]},
z == (Lz + dL)/L && 0 < x < 1]}, \[Phi], {x, z} \[Element]
mesh] // Quiet;
p0 = p[i - 1][x, z, T[[i + 1]]]; {p[i], pz0[i]} =
NDSolveValue[{eq2 == 0, P[x, z, T[[i]]] == p0,
lambda/L Derivative[0, 1, 0][P][x, Lz/L, t] == P[x, Lz/L, t],
lambda/L Derivative[0, 1, 0][P][x, 0, t] == 0,
Derivative[1, 0, 0][P][0, z, t] == 0,
Derivative[1, 0, 0][P][1, z, t] == 0}, {P,
Derivative[0, 1, 0][P][x, z, T[[i + 1]]]}, {x, 0, 1}, {z, 0,
Lz/L}, {t, T[[i]], T[[i + 1]]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 41, "MaxPoints" -> 81,
"DifferenceOrder" -> 4}}] // Quiet;
pz[i] = Interpolation[
Flatten[Table[{{x, z}, If[z < Lz/L, pz0[i], 0]}, {x, 0,
1, .02}, {z, 0, (Lz + dL)/L, .005 (Lz + dL)/L}], 1]];, {i, 1,
151}]


It takes about 765 s on my notebook, finally we can visualize electric potential and polarization on the last step

{Plot3D[phi[x, z], {x, z} \[Element] mesh,
ColorFunction -> "Rainbow", PlotTheme -> "Marketing",
MeshStyle -> White, AxesLabel -> Automatic],
Plot3D[p[x, z, T[]], {x, z} \[Element] mesh1,
ColorFunction -> "Rainbow", PlotTheme -> "Marketing",
MeshStyle -> White, AxesLabel -> Automatic, PlotRange -> All,
PlotPoints -> 100]} Then we need to compute average potential and polarization in the form

VFE = Table[
NIntegrate[(-Derivative[0, 1][phi[i]][x, z]), {x, 0, 1}, {z, 0,
Lz/L}, PrecisionGoal -> 2, AccuracyGoal -> 2], {i, 1, 151}] //
Quiet;
pav = L/Lz Table[
NIntegrate[
p[i][x, z,
T[[i + 1]]] + \[Epsilon]0*\[Epsilon]z/
L (-Derivative[0, 1][phi[i]][x, z]), {x, 0, 1}, {z, 0, Lz/L},
PrecisionGoal -> 2, AccuracyGoal -> 2], {i, 1, 151}] // Quiet;


To compare with static case from the paper we need to scale polarization and potential with using $$1/\sqrt{2}$$, therefore we have unscaled picture (left) and scaled (right)

{ListPlot[Table[{VFE[[i]], 100 pav[[i]]}, {i, Length[VFE]}],
PlotRange -> All],
ListPlot[Drop[
Table[{VFE[[i]], 100 pav[[i]]}/Sqrt, {i, Length[VFE]}], 7],
Frame -> True, AspectRatio -> 1/1,
FrameLabel -> {"\!$$\*SubscriptBox[\(VFE$$, $$avg$$]\), V",
"\!$$\*SubscriptBox[\(P$$, $$avg$$]\), \
\[Mu]C/\!$$\*SuperscriptBox[\(cm$$, $$2$$]\)"}, PlotStyle -> Blue,
PlotRange -> {{-6.5, 6.5}, {-40, 40}},
PlotMarkers -> "OpenMarkers"]} Note, that scale factor 100 at pav used for units $$\mu C/cm^2$$. For dielectric layer with thickness of 3 nm we need to increase voltage ap to Vapp[t_] := 11 Sin[2 Pi t] then result (red points) together with the previous case 1 nm (blue points) shown below In a case of dielectric layer with thickness 5 nm we increase voltage to Vapp[t_] := 20 Sin[2 Pi t], and have result (green points) together with the case 1 nm (blue points) and 3 nm (red points) • Fabulous!! Thank you very much for your help :) Can I ask you some questions later regarding the code? Currently, I'm trying to understand how it works!
– YSP
Apr 19, 2022 at 7:18
• @YSP Please, see update for the 3 nm dielectric layer case. Apr 19, 2022 at 9:25
• @YSP See also result for 5 nm layer with Vapp[t_] := 20*Sin[\[Omega]*t] Apr 20, 2022 at 11:36
• So appreciative...I'm actually learning a lot from your code. I see you scaled the dimension and parameters to length L. Is that to reduce the computation time?
– YSP
Apr 20, 2022 at 22:50
• @YSP Actually we used scaling in space as x/L and in time as t/t0 to make this code more stable. Apr 21, 2022 at 2:03