# Numerical modeling and investigations of 3D devices with ferroelectric layer fully embedded in a paraelectric environment

Pierre-William Martelli, Séraphin Mefire

► To cite this version:

Pierre-William Martelli, Séraphin Mefire. Numerical modeling and investigations of 3D devices with ferroelectric layer fully embedded in a paraelectric environment. Computers & Mathematics with Applications, 2018, 76 (3), pp.451-470. 10.1016/j.camwa.2018.04.027. hal-02937221

HAL Id: hal-02937221

https://hal.science/hal-02937221

Submitted on 12 Sep 2020

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
Numerical modeling and investigations of 3D devices
with ferroelectric layer fully embedded in a
paraelectric environment

P.-W. Martelli* and S.M. Mefire

Université de Lorraine, CNRS, IECL, F-54000 Nancy, France.

*Corresponding author. Email address: pierre-william.martelli@univ-lorraine.fr

**Abstract.** We investigate three-dimensional devices made up of a ferroelectric layer that is fully embedded in a paraelectric environment, by modeling based on the Ginzburg-Landau formalism as well as on the Electrostatics equations, and boundary conditions that are suitable for applications. From finite element approximations and inexact Newton techniques for solving numerically the resulting nonlinear system, we develop two numerical protocols. The first protocol concerns a determination of states related to the system and incorporates a process of heating as well as of cooling of devices, in terms of the temperature, whereas the second one is devoted to the existence study of hysteresis loops. The efficiency of these protocols is particularly emphasized from numerical simulations.

**Keywords.** Electroactive media, Ferroelectricity, Electrostatics, Ginzburg-Landau systems, Finite Elements, Inexact Newton methods, Numerical Simulations.

**MS Codes.** 82D45, 82D80, 78A25, 35Q56, 47H10, 65N30, 49M15

# 1 Introduction

In the frame of the development of new memory concepts, the main desired qualities are the non-volatility of the stored information, the low power consumption and the low writing/reading times. Among the emerging concepts, the ferroelectric memory (FeRAM) has a significant lead over its competitors (see, e.g., [22]). The FeRAM is based on ferroelectric devices which have a hysteretic behavior and thus allow to store bits. The miniaturization required for associated applications led to the discovery of new nanoscopic properties of devices, due to long-range Coulomb interactions which are no longer negligible (see, e.g., [14]). A notable property is that the ferroelectric body can be organized into a finite number of distinct regions, the *domains*, separated by *domain walls*, in which the electric polarization field is in opposite directions. A pertinent numerical investigation of such devices requires thus a description of the polarization texture as accurate as possible, in terms of the physical and geometrical parameters.

Usually, two numerical approaches are considered for studying such devices: the approach based on first-principia simulations (see, e.g., [13, 24]), related to a good accuracy but having the particularity to be however subject to expensive CPU costs, and the approach that consists of solving Partial Differential Equations, namely very suited to a description of the parameters dependence as regards the results.

Two phases are in association with a ferroelectric body. At high-temperature, we have the paraelectric phase, in which the polarization field depends linearly on the electric field. Below a certain temperature, called the Curie temperature (see, e.g., [14]), this dependence becomes nonlinear: it is the ferroelectric phase. In the 1940s, V.L. Ginzburg applied the Landau phase transition theory to ferroelectricity (see, e.g., [11]). E.V. Chenskii and V.V. Tarasenko [6] had the idea, in the 1980s, to complete the original model in order to take into account the long-range Coulomb interactions, in particular by making also use of the Electrostatics equations.

Recently, I.A. Luk'yanchuck et al. [15] have established the universal properties of ferroelectric domains from this model together with periodic conditions, for the context of 2D devices made up of a ferroelectric film partly embedded in a paraelectric environment. In this work, we are concerned with 3D devices made
up of a ferroelectric layer that is fully embedded in a paraelectric environment, and where the geometrical configurations of nanofilms, nanorods or nanodots can be considered. We are namely interested in an extended frame of investigations in the sense that the considerations of the physical and geometrical parameters can be also the ones that are not systematically in correspondence with particular concrete situations, in contrast with [16] where the considerations of the parameters are essentially related to applications by involving furthermore the geometrical context of nanodots.

This paper is subdivided in six sections. We start by describing a model based on the Ginzburg-Landau formalism and Electrostatics equations, together with boundary conditions that are suitable for applications, for studying such 3D devices. In Section 3, we associate a weak formulation with the nonlinear system that represents this model, and establish the existence of weak solutions. We are concerned, in Section 4, with a numerical approach using Finite Elements for the discretization of the formulation and inexact Newton techniques for solving the associated discrete nonlinear system. In Section 5, we develop two numerical protocols and present the obtained results. The first protocol is devoted to the determination of states, and incorporates also a process of heating as well as of cooling of devices, in terms of the temperature. The second one aims at investigating the existence of hysteresis loops. The performed numerical simulations deal with various considerations of geometrical and physical parameters. We report in Section 6 the concluding remarks.

## 2 The ferroelectric model

We are here interested in a model for the study of 3D devices made up of a ferroelectric layer that is fully embedded in a paraelectric environment. Let us represent geometrically such a device by an open bounded subset $\Omega$ of $\mathbf{R}^3$, and its fully embedded layer by an open subset $\Omega_f$, $\overline{\Omega_f} \subset \Omega$, as well as its paraelectric environment by $\Omega_p = \Omega \setminus \overline{\Omega_f}$. Also, let us denote by $S$, $S = \partial\Omega$, the boundary of $\Omega$, and by $S_f$, $S_f = \overline{\Omega_f} \cap \overline{\Omega_p}$, the interface between $\Omega_f$ and $\Omega_p$. The region $\Omega$ is also such that two subparts of its boundary are included in two planes, ($z = -\gamma$) and ($z = \gamma$), so it is possible to affix there electrodes, in order to prescribe a voltage potential (see, e.g., Figure 1). It is assumed that the ferroelectric body is uniaxial, namely that the dependence of the polarization field, $\mathbf{P}$, on the electric field is nonlinear through only one of its components (here, the third one). More precisely, if we represent by $\mathbf{E} = (E_x, E_y, E_z)^T$ the electric field, where $E_x, E_y, E_z$ are of course scalar fields, and the superscript "T" denotes here as well as in the next sections the transpose, we then have (see, e.g., [14]) that:

$$ \mathbf{P} = \frac{1}{4\pi}(\varepsilon\mathbf{E} - \mathbf{E}) + \begin{pmatrix} 0 \\ 0 \\ P \end{pmatrix} \quad \text{in } \Omega_f, \quad (1) $$

with $P$ a scalar field, $P: \Omega_f \to \mathbf{R}$, and $\varepsilon$ the electric permittivity.

<figure><img src="image_2.png" /><figcaption>Figure 1: An example of geometrical configuration of the device.</figcaption></figure>

By making use of the general theory presented in [14] or of the same arguments as those described in [15],
we have here that the generating functional energy, $\mathcal{E}$, is such that:

$$ \mathcal{E}(\mathbf{P}, \mathbf{E}) = \int_{\Omega} \tilde{\mathcal{F}}(\mathbf{P}, \mathbf{E})dx, \qquad (2) $$

where $\tilde{\mathcal{F}}(\mathbf{P}, \mathbf{E})$ is expressed as follows,

$$ \left\{ \begin{array}{ll} \tilde{\mathcal{F}}(\mathbf{P}, \mathbf{E}) = \mathcal{F}(P, 0) - E_z P - \frac{\mathbf{E} \cdot (\varepsilon \mathbf{E})}{8\pi} & \text{in } \Omega_f, \\ \tilde{\mathcal{F}}(\mathbf{P}, \mathbf{E}) = -\frac{\mathbf{E} \cdot (\varepsilon \mathbf{E})}{8\pi} & \text{in } \Omega_p, \end{array} \right. $$

with $\mathcal{F}(P, 0)$ the Ginzburg-Landau potential. The term $-\int_{\Omega_f} E_z P dx$ penalizes the energy if $E_z$ and $P$ are of opposed signs, and thus traduces the effect of depolarizing field. The energy $-\frac{1}{8\pi} \int_{\Omega} \mathbf{E} \cdot (\varepsilon \mathbf{E}) dx$ represents the effect of the domain walls. Thanks to the Ginzburg-Landau theory of phase transitions (cf., e.g., [14, 15]), the term $\mathcal{F}(P, 0)$ can be written with respect to powers of $P$,

$$ \mathcal{F}(P, 0) = \frac{4\pi}{\varkappa_{\parallel}} \left( \frac{t}{2} P^2 + \frac{1}{4} \frac{1}{P_0^2} P^4 + \frac{1}{2} \nabla P \cdot (\xi \nabla P) \right), $$

with

* $t$, the reduced temperature, which can be expressed via the Curie temperature $T_0$ and the temperature $T$ of the device as: $t = T/T_0 - 1$;

* $\varkappa_{\parallel}$, the displacive parameter (namely positive);

* $P_0, P_0 \neq 0$, the spontaneous polarization at low temperatures;

* $\xi$, the coherence length tensor.

In this expansion, the odd powers of $P$ do not occur because of symmetry reasons. The term $\nabla P \cdot (\xi \nabla P)$ penalizes spatial inhomogeneities: it becomes important when ferroelectric domains are considered.

Furthermore, if $\Omega$ is simply connected, the electric field $\mathbf{E}$ can be expressed with the help of a scalar electric potential, $\varphi$, as: $\mathbf{E} = -\nabla \varphi$.

Finally, the energy of the device can be written from (2), with the help of $P$ and $\varphi$, also with $\varkappa_{\parallel} > 0$, as follows,

$$ \mathcal{E}(\mathbf{P}, \mathbf{E}) = \int_{\Omega_f} \frac{4\pi}{\varkappa_{\parallel}} \left( \frac{t}{2} P(x)^2 + \frac{1}{4} \frac{1}{P_0^2} P(x)^4 + \frac{1}{2} \nabla P(x) \cdot (\xi(x) \nabla P(x)) \right) dx + \int_{\Omega_f} \partial_z \varphi(x) P(x) dx - \int_{\Omega} \frac{1}{8\pi} \nabla \varphi(x) \cdot (\varepsilon(x) \nabla \varphi(x)) dx. $$

It arises from the Euler-Lagrange variation of the energy, $\delta\mathcal{E} = 0$, that:

$$ tP + \left(\frac{P}{P_0}\right)^2 P - \nabla \cdot (\xi \nabla P) + \frac{\varkappa_{\parallel}}{4\pi} \partial_z \varphi = 0 \quad \text{in } \Omega_f, \qquad (3) $$

and

$$ \left\{ \begin{array}{ll} -\nabla \cdot (\varepsilon \nabla \varphi) + 4\pi \partial_z P = 0 & \text{in } \Omega_f, \\ -\nabla \cdot (\xi \nabla \varphi) = 0 & \text{in } \Omega_p, \end{array} \right. \qquad (4) $$

as also reported with details in [17]. The relations in (3) and (4) involve the Ginzburg-Landau formalism and Electrostatics equations. In particular, $\nabla \cdot \mathbf{D} = 0$ in $\Omega_f$ and $\Omega_p$, where $\mathbf{D}$ is the electric displacement field,

$$ \mathbf{D} = \begin{cases} \varepsilon \mathbf{E} + \begin{pmatrix} 0 \\ 0 \\ 4\pi P \end{pmatrix} & \text{in } \Omega_f, \\ \varepsilon \mathbf{E} & \text{in } \Omega_p. \end{cases} $$

These equations must be completed by conditions at the interface $S_f = \overline{\Omega_f} \cap \overline{\Omega_p}$ and on the boundary $S = \partial \Omega$. Namely, at the ferroelectric/paraelectric interface, the scalar field $\varphi$ satisfies the continuity condition,

$$ \varphi|_{\Omega_f} = \varphi|_{\Omega_p} \quad \text{on } S_f, \qquad (5) $$
and the electric displacement field is such that:

$$ \nu \cdot (\mathbf{D}|_{\Omega_f} - \mathbf{D}|_{\Omega_p}) = 0 \quad \text{on } S_f, $$

where $\nu = (\nu_x, \nu_y, \nu_z)^T$ is the outward unit normal to $\Omega_f$, defined on $S_f$. It follows that:

$$ \nu \cdot ((\varepsilon \nabla \varphi)|_{\Omega_f} - (\varepsilon \nabla \varphi)|_{\Omega_p}) - 4\pi\nu_z P = 0 \quad \text{on } S_f. \qquad (6) $$

Also, the scalar field $P$ is subject to the condition:

$$ \nu \cdot (\xi \nabla P) = 0 \quad \text{on } S_f. \qquad (7) $$

Finally, the electrodes affixed on the device require the potentials $\mp\phi_0$ on the planes ($z = \pm\gamma$), $\gamma > 0$, generating a voltage potential $U = 2\phi_0$ (see Figure 1); that can be expressed through a boundary condition as follows,

$$ \varphi = \varphi_S \quad \text{on } S, \qquad (8) $$

with $\varphi_S$ a given function, taking the values $\mp\phi_0$ on ($z = \pm\gamma$) in particular.

Throughout the paper, the set of equations in (3)-(8),

$$ \text{(EGL)} \quad \left\{ \begin{array}{ll} tP + \left(\frac{P}{P_0}\right)^2 P - \nabla \cdot (\xi \nabla P) + \frac{\kappa_{\parallel}}{4\pi} \partial_z \varphi = 0 & \text{in } \Omega_f, \\ -\nabla \cdot (\varepsilon \nabla \varphi) + 4\pi \partial_z P = 0 & \text{in } \Omega_f, \\ -\nabla \cdot (\varepsilon \nabla \varphi) = 0 & \text{in } \Omega_p, \\ \varphi|_{\Omega_f} = \varphi|_{\Omega_p} & \text{on } S_f, \\ \nu \cdot ((\varepsilon \nabla \varphi)|_{\Omega_f} - (\varepsilon \nabla \varphi)|_{\Omega_p}) - 4\pi\nu_z P = 0 & \text{on } S_f, \\ \nu \cdot (\xi \nabla P) = 0 & \text{on } S_f, \\ \varphi = \varphi_S & \text{on } S, \end{array} \right. $$

is called the *Electrostatics Ginzburg-Landau (EGL)* system.

Although, in the applications, the temperature $T$ must be considered in such a way that the reduced temperature $t$ takes values in $[-1, 0)$, we will be interested in a more large context of values for this parameter, namely by investigating with $t \in (-\infty, 0]$. Also, we will take into account the extended consideration where $\kappa_{\parallel} \ge 0$ in the (EGL) system.

Let us mention furthermore that the volume of any considered ferroelectric layer $\Omega_f$, in order to avoid to disturb the emergent fringing field (see, e.g., [16]), will be of course sufficiently more small than the one of the medium $\Omega$, but not too small. In particular also, the related considerations will not then be subject to numerical instabilities in experiments, although the asymptotic framework developed by H. Ammari et al. in [1] may remedy this issue.

Of course, in the frame of the study of ferroelectric materials or devices, modeling approaches deal also with the contexts of ferroelectric layers that are not embedded or that are only partly embedded, and involve the time dependence (see, e.g., [2, 17, 22, 23] and the references therein). In particular in [2], H. Ammari and K. Hamdache develop a theoretical analysis for studying a model where the electric polarization is time dependent and a ferroelectric material represents alone the 3D physical region of interest.

# 3 Weak formulation and existence of solutions

We associate a weak formulation with the (EGL) system, and study the existence of solutions of this formulation.

## 3.1 Weak formulation
Let $t \le 0$, $\varkappa_\parallel \ge 0$ and $P_0 \ne 0$ the parameters considered as previously. We assume that each of the bounded open subsets of $\mathbf{R}^3$ that are $\Omega_f$ and $\Omega$ has a Lipschitz-continuous boundary. Also, we consider the tensors $\xi$ and $\varepsilon$, where $\xi = (\xi_{i,j})_{1\le i,j\le 3}$, $\varepsilon = (\varepsilon_{i,j})_{1\le i,j\le 3}$, in such a way that:

$$ \xi_{i,j} \in L^{\infty}(\Omega_f), \forall 1 \le i, j \le 3, \qquad (9) $$

$$ \exists \xi_0 > 0; \forall d = (d_1, d_2, d_3)^T \in \mathbf{R}^3, \sum_{i,j=1}^{3} \xi_{i,j} d_j d_i \ge \xi_0 |d|^2 \text{ a.e. in } \Omega_f, \qquad (10) $$

$$ \varepsilon_{i,j} \in L^{\infty}(\Omega), \forall 1 \le i, j \le 3, \qquad (11) $$

$$ \exists \varepsilon_0 > 0; \forall d = (d_1, d_2, d_3)^T \in \mathbf{R}^3, \sum_{i,j=1}^{3} \varepsilon_{i,j} d_j d_i \ge \varepsilon_0 |d|^2 \text{ a.e. in } \Omega, \qquad (12) $$

with $|\cdot|$ denoting here the euclidian norm on $\mathbf{R}^3$.

We set:

$$ W = H^1(\Omega_f) \times H_0^1(\Omega), $$

and represent by $(\cdot, \cdot)_W$ the scalar product defined on $W \times W$, where for $(P, \varphi), (Q, \psi) \in W$,

$$ ((P, \varphi), (Q, \psi))_W = \int_{\Omega_f} PQdx + \int_{\Omega_f} \nabla P \cdot \nabla Qdx + \int_{\Omega} \nabla \varphi \cdot \nabla \psi dx. $$

The associated norm is denoted by $\|\cdot\|_W$, where for $(P, \varphi) \in W$,

$$ \|(P, \varphi)\|_W = \left( \|P\|_{H^1(\Omega_f)}^2 + |\varphi|_{H^1(\Omega)}^2 \right)^{1/2}. $$

Assume that we have $(P, \tilde{\varphi}) \in H^2(\Omega_f) \times H^2(\Omega)$ satisfying (EGL) with $\xi$, $\varepsilon$ and $\varphi_S$ sufficiently regular, by defining $\varphi = \tilde{\varphi} - w$, where $w$ is a scalar field sufficiently regular and such that $w|_S = \varphi_S$; namely, in particular, $\varphi|_S = 0$. Let $(Q, \psi) \in W$. From (3) and with (7), it results thus that:

$$ \int_{\Omega_f} tPQdx + \int_{\Omega_f} \frac{P^3}{P_0^2}Qdx + \int_{\Omega_f} (\xi\nabla P) \cdot \nabla Qdx + \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi}(\partial_z\varphi)Qdx = - \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi}(\partial_z w)Qdx. $$

By only considering $P$ such that $P \in H^1(\Omega_f)$, it derives in particular (see, e.g., [12]) that $P \in L^6(\Omega_f)$, and therefore the term $\int_{\Omega_f} \frac{1}{P_0^2}P^3Qdx$ is also well defined.

With the same assumptions as above, it follows from (4), and by making use of the transmission condition (6), that:

$$ \int_{\Omega} (\varepsilon \nabla \varphi) \cdot \nabla \psi dx - \int_{S_f} 4\pi \nu_z P \psi d\sigma + \int_{\Omega_f} 4\pi (\partial_z P) \psi dx = - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla \psi dx. $$

Finally, the weak formulation associated with (EGL) consists of finding $(P, \varphi) \in W$ such that, for $(Q, \psi) \in W$,

$$ \left\{ \int_{\Omega_f} tPQdx + \int_{\Omega_f} \frac{P^3}{P_0^2}Qdx + \int_{\Omega_f} (\xi\nabla P) \cdot \nabla Qdx + \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi}(\partial_z\varphi)Qdx = - \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi}(\partial_z w)Qdx, \right. \qquad (13) $$

$$ \left. \int_{\Omega} (\varepsilon \nabla \varphi) \cdot \nabla \psi dx - \int_{S_f} 4\pi \nu_z P \psi d\sigma + \int_{\Omega_f} 4\pi (\partial_z P) \psi dx = - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla \psi dx, \right. \qquad (14) $$

with $\xi$ and $\varepsilon$ satisfying (9)-(12), $w$ considered in such a way that $w \in H^1(\Omega)$ and $w|_S = \varphi_S$, where $\varphi_S \in H^{1/2}(S)$, is given.
## 3.2 Existence of weak solutions

In order to prove the existence of solutions of the weak formulation, we proceed in two stages, by making use also of an idea of [8]. We first consider a sequence of weak problems, also related to (13)-(14), but in which a truncation of certain terms at "the n-level" (made explicit below) is done. We show, with the help of the Schauder Theorem (see, e.g., [10]), that each of these problems has a solution (see Propositions 1 and 2). Next, we extract from this sequence of solutions a subsequence (see Theorem 1) for which we prove, by using the Vitali Theorem (see, e.g., [9]), the convergence to a solution of the formulation (13)-(14).

For all $n, n \in \mathbf{N}$, we define the n-level truncation function, $T_n$, as follows: $x \in \mathbf{R} \mapsto \max\{\min\{x, n\}, -n\} \in \mathbf{R}$, also expressed as

$$T_n(x) = \begin{cases} -n & \text{if } x \leq -n, \\ x & \text{if } -n \leq x \leq n, \\ n & \text{if } x \geq n. \end{cases}$$

The relations below derive then from this definition,

$$\left\{ \begin{array}{l} |T_n(x)| \leq n, \quad \forall (n, x) \in \mathbf{N} \times \mathbf{R}, \\ |T_n(x)| \leq |x|, \quad \forall (n, x) \in \mathbf{N} \times \mathbf{R}, \\ \forall x \in \mathbf{R}, \exists N_x \in \mathbf{N}; \forall n \geq N_x, T_n(x) = x. \end{array} \right. \tag{15}$$

Moreover, $T_n$ satisfies the property expressed in the following lemma.

**Lemma 1.** Let $f: \mathbf{R} \to \mathbf{R}$ a continuously differentiable function, satisfying $\lim_{x \to -\infty} f(x) = +\infty$ and $\lim_{x \to +\infty} f(x) = -\infty$. Then, for $n \in \mathbf{N}$, $T_n \circ f$ is a continuous Lipschitzian function.

This other property of $T_n$, considering such a function $f$, can be established as done in [17] from some explicit calculations or by noting directly that the derivative of $T_n \circ f$ in $\mathbf{R} \setminus \{-n, n\}$ is bounded.

Let us introduce the two preliminary statements (see Propositions 1 and 2) of the present study of the existence of weak solutions. Also, in order to simplify the notation hereafter, let us denote by $g$ the function defined as follows,

$$g: x \in \mathbf{R} \mapsto (1-t)x - \frac{x^3}{P_0^2} \in \mathbf{R},$$

with of course $t$ and $P_0$ considered as previously.

**Proposition 1.** Let $n \in \mathbf{N}$ and $(\bar{P}_n, \bar{\varphi}_n) \in W$. Furthermore, let $\varkappa_\parallel > 0$ and $g$ be defined as aforementioned, as well as $w$, by assuming also that $\xi$ and $\varepsilon$ satisfy (9)-(12). Then, there exists a unique $(P_n, \varphi_n) \in W$ satisfying, for all $(Q, \psi) \in W$,

$$\left\{ \begin{aligned} & \int_{\Omega_f} P_n Q dx + \int_{\Omega_f} (\xi \nabla P_n) \cdot \nabla Q dx + \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi} (\partial_z \varphi_n) Q dx \\ &= \int_{\Omega_f} (T_n \circ g \circ \bar{P}_n) Q dx - \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi} (\partial_z w) Q dx, \\ & \int_{\Omega} (\varepsilon \nabla \varphi_n) \cdot \nabla \psi dx - \int_{S_f} 4\pi \nu_z P_n \psi d\sigma + \int_{\Omega_f} 4\pi (\partial_z P_n) \psi dx \\ &= - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla \psi dx. \end{aligned} \right. \tag{18}$$

*Proof.* If of course there exists $(P_n, \varphi_n)$ satisfying (18), then by adding the first relation of (18) to the second one multiplied by $\frac{\varkappa_\parallel}{16\pi^2}$, we obtain:

$$a((P_n, \varphi_n), (Q, \psi)) = L(Q, \psi) \quad \forall (Q, \psi) \in W, \tag{19}$$

where,

$$a((P_n, \varphi_n), (Q, \psi)) = \int_{\Omega_f} P_n Q dx + \int_{\Omega_f} (\xi \nabla P_n) \cdot \nabla Q dx + \int_{\Omega} \frac{\varkappa_\parallel}{16\pi^2} (\varepsilon \nabla \varphi_n) \cdot \nabla \psi dx \\ + \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi} (\partial_z \varphi_n) Q dx - \int_{S_f} \frac{\varkappa_\parallel}{4\pi} \nu_z P_n \psi d\sigma + \int_{\Omega_f} \frac{\varkappa_\parallel}{4\pi} (\partial_z P_n) \psi dx,$$
and

$$L(Q, \psi) = \int_{\Omega_f} (T_n \circ g \circ \bar{P}_n) Q dx - \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z w) Q dx - \int_{\Omega} \frac{\kappa_{\parallel}}{16\pi^2} (\varepsilon \nabla w) \cdot \nabla \psi dx.$$

The bilinear map $a(\cdot, \cdot)$ and the linear map $L(\cdot)$ satisfy the Lax-Milgram Theorem hypotheses. In fact, we have, for $(P, \varphi), (Q, \psi) \in W$,

$$\begin{aligned} |a((P, \varphi), (Q, \psi))| & \le \left|\int_{\Omega_f} P Q dx\right| + \left|\int_{\Omega_f} (\xi \nabla P) \cdot \nabla Q dx\right| + \left|\int_{\Omega} \frac{\kappa_{\parallel}}{16\pi^2} (\varepsilon \nabla \varphi) \cdot \nabla \psi dx\right| \\ & \qquad + \left|\int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z \varphi) Q dx\right| + \left|\int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} P (\partial_z \psi) dx\right| \\ & \le C \|(P, \varphi)\|_W \|(Q, \psi)\|_W, \end{aligned}$$

where $C = 1 + 9 \max_{1 \le i,j \le 3} \|\xi_{i,j}\|_{L^\infty(\Omega_f)} + \frac{9\kappa_{\parallel}}{16\pi^2} \max_{1 \le i,j \le 3} \|\varepsilon_{i,j}\|_{L^\infty(\Omega)} + \frac{\kappa_{\parallel}}{2\pi}$. Also, we note with $(P, \varphi) \in W$ that:

$$\begin{aligned} a((P, \varphi), (P, \varphi)) & \ge \|P\|_{L^2(\Omega_f)}^2 + \xi_0 |P|_{H^1(\Omega_f)}^2 + \frac{\kappa_{\parallel}}{16\pi^2} \varepsilon_0 |\varphi|_{H^1(\Omega)}^2 \\ & \ge \min \left\{1, \xi_0, \frac{\kappa_{\parallel}}{16\pi^2} \varepsilon_0\right\} \|(P, \varphi)\|_W^2. \end{aligned}$$

Furthermore, for $(Q, \varphi) \in W$, we have, thanks to (15),

$$\left|\int_{\Omega_f} (T_n \circ g \circ \bar{P}_n) Q dx\right| \le n \int_{\Omega_f} |Q| dx \le n \sqrt{|\Omega_f|} \|Q\|_{L^2(\Omega_f)} \le n \sqrt{|\Omega_f|} \|(Q, \psi)\|_W.$$

Thus,

$$|L(Q, \psi)| \le \left( n \sqrt{|\Omega_f|} + \frac{\kappa_{\parallel}}{4\pi} \|\partial_z w\|_{L^2(\Omega_f)} + \frac{9\kappa_{\parallel}}{16\pi^2} \max_{1 \le i,j \le 3} \|\varepsilon_{i,j}\|_{L^\infty(\Omega)} |w|_{H^1(\Omega)} \right) \|(Q, \psi)\|_W.$$

Finally, we obtain from the Lax-Milgram Theorem that there exists a unique $(P_n, \varphi_n) \in W$ satisfying (19), for all $(Q, \psi) \in W$. The first equation of (18) can be obtained from (19) by taking $\psi = 0$, and the second one by taking $Q = 0$ also from (19). □

**Proposition 2.** Let $n \in \mathbb{N}$. Furthermore, let $\kappa_{\parallel} > 0$ and $g$ be defined as aforementioned, as well as $w$, by assuming also that $\xi$ and $\varepsilon$ satisfy (9)-(12). Then, there exists $(P_n, \varphi_n) \in W$ such that, for all $(Q, \psi) \in W$,

$$\left\{ \begin{aligned} & \int_{\Omega_f} P_n Q dx - \int_{\Omega_f} (T_n \circ g \circ P_n) Q dx + \int_{\Omega_f} (\xi \nabla P_n) \cdot \nabla Q dx \\ & \qquad + \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z \varphi_n) Q dx = - \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z w) Q dx, \\ & \int_{\Omega} (\varepsilon \nabla \varphi_n) \cdot \nabla \psi dx - \int_{S_f} 4\pi \nu_z P_n \psi d\sigma + \int_{\Omega_f} 4\pi (\partial_z P_n) \psi dx \\ & \qquad = - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla \psi dx. \end{aligned} \right. \quad (20)$$

*Proof.* For each $(\bar{P}_n, \bar{\varphi}_n) \in W$, we denote by $(P_n, \varphi_n)$, where $(P_n, \varphi_n) \in W$, the solution of (18), which allows us to introduce $f_n : (\bar{P}_n, \bar{\varphi}_n) \in W \mapsto (P_n, \varphi_n) \in W$. The idea of the proof is to check that $f_n$ satisfies the hypotheses of the Schauder fixed point Theorem. Each fixed point of $f_n$ will be a solution of (20). We organize the proof into three parts. Firstly, we establish that $f_n$ is bounded and so maps a closed ball $B$ into itself. Next, we prove that $f_n(B)$ is relatively compact in $W$, and finally that $f_n$ is continuous. Thus,

* let us start by showing that the image of $f_n$ is bounded. Let $(P_n, \varphi_n) \in f_n(W)$. There exists $(\bar{P}_n, \bar{\varphi}_n) \in W$ such that $f_n(\bar{P}_n, \bar{\varphi}_n) = (P_n, \varphi_n)$. We have, by choosing $(P_n, \varphi_n)$ as a test-function in (18), that: $\|(P_n, \varphi_n)\|_W \le \frac{C_L}{\alpha}$, where $C_L$ is a continuity constant of the linear map $L(\cdot)$ and $\alpha$ is an ellipticity constant of the bilinear map $a(\cdot, \cdot)$, resulting from the proof of Proposition 1. Then, $f_n(W)$ is bounded and, if $B \subset W$ is the ball of center $0_W$ and radius $C_L/\alpha$, it follows that $f_n|_B : B \to B$ is correctly defined;
• let us now note that $f_n(B)$ is relatively compact in $W$. We consider $(P_n^{(m)}, \varphi_n^{(m)})_{m \in \mathbb{N}}$ a sequence of elements of $f_n(B)$. For all $m \in \mathbb{N}$, there exists $(\bar{P}_n^{(m)}, \bar{\varphi}_n^{(m)}) \in B$ satisfying $f_n(\bar{P}_n^{(m)}, \bar{\varphi}_n^{(m)}) = (P_n^{(m)}, \varphi_n^{(m)})$. Since $(\bar{P}_n^{(m)}, \bar{\varphi}_n^{(m)})_{m \in \mathbb{N}}$ is bounded, we can extract a weakly convergent subsequence in $W$, represented by $(\bar{P}_n^{(m_k)}, \bar{\varphi}_n^{(m_k)})_{k \in \mathbb{N}}$. We denote by $(\bar{P}_n, \bar{\varphi}_n)$ its limit. Let $(P_n, \varphi_n) = f_n(\bar{P}_n, \bar{\varphi}_n)$. For $k \in \mathbb{N}$, let us set $\alpha_{m_k} = P_n^{(m_k)} - P_n$ and $\beta_{m_k} = \varphi_n^{(m_k)} - \varphi_n$. By choosing $(\alpha_{m_k}, \beta_{m_k})$ as a test-function in (18), we get:

$$ 
\begin{aligned}
& \int_{\Omega_f} \alpha_{m_k}^2 dx + \int_{\Omega_f} (\xi \nabla \alpha_{m_k}) \cdot \nabla \alpha_{m_k} dx + \int_{\Omega} \frac{\varkappa_{\parallel}}{16\pi^2} (\varepsilon \nabla \beta_{m_k}) \cdot \nabla \beta_{m_k} dx \\
& + \int_{\Omega_f} \frac{\varkappa_{\parallel}}{4\pi} (\partial_z \beta_{m_k}) \alpha_{m_k} dx + \int_{\Omega_f} \frac{\varkappa_{\parallel}}{4\pi} (\partial_z \alpha_{m_k}) \beta_{m_k} dx - \int_{S_f} \frac{\varkappa_{\parallel}}{4\pi} \nu_z \alpha_{m_k} \beta_{m_k} d\sigma \\
& = \int_{\Omega_f} (T_n \circ g \circ \bar{P}_n^{(m_k)} - T_n \circ g \circ \bar{P}_n) \alpha_{m_k} dx.
\end{aligned}
$$

We deduce, by using Lemma 1, that: $\|(\alpha_{m_k}, \beta_{m_k})\|_W^2 \le \frac{K_n}{\alpha} \|\bar{P}_n^{(m_k)} - \bar{P}_n\|_{L^2(\Omega_f)} \|(\alpha_{m_k}, \beta_{m_k})\|_W$, where $K_n$ is a Lipschitz constant of $T_n \circ g$. Thus, $\|(\alpha_{m_k}, \beta_{m_k})\|_W \le \frac{K_n}{\alpha} \|\bar{P}_n^{(m_k)} - \bar{P}_n\|_{L^2(\Omega_f)}$. As the injection of $H^1(\Omega_f)$ in $L^2(\Omega_f)$ is compact and since $(\bar{P}_n^{(m_k)})_{k \in \mathbb{N}}$ converges weakly to $\bar{P}_n$ in $H^1(\Omega_f)$, we have that $(\bar{P}_n^{(m_k)}, \varphi_n^{(m_k)})_{k \in \mathbb{N}}$ converges to $\bar{P}_n$ in $L^2(\Omega_f)$ and the previous inequality leads to the convergence of $(P_n^{(m_k)}, \varphi_n^{(m_k)})_{k \in \mathbb{N}}$ to $(P_n, \varphi_n)$ in $W$. So $f_n(B)$ is relatively compact in $W$;

• it remains to establish the continuity of $f_n$. Let $(\bar{P}_n^{(m)}, \bar{\varphi}_n^{(m)})_{m \in \mathbb{N}}$ be a sequence which converges to $(\bar{P}_n, \bar{\varphi}_n)$ in $W$. By repeating the above calculation and similar notations, we note that:

$$ \|(\alpha_m, \beta_m)\|_W \le \frac{K_n}{\alpha} \|\bar{P}_n^{(m)} - \bar{P}_n\|_{L^2(\Omega_f)} \le \frac{K_n}{\alpha} \|(\bar{P}_n^{(m)} - \bar{P}_n, \bar{\varphi}_n^{(m)} - \bar{\varphi}_n)\|_W. $$

Thus, $(f_n(\bar{P}_n^{(m)}, \bar{\varphi}_n^{(m)}))_{m \in \mathbb{N}}$ converges to $f_n(\bar{P}_n, \bar{\varphi}_n)$, and the continuity of $f_n$ is established.

Finally, we can apply the Schauder Theorem and deduce that $f_n$ has a fixed point; there exists $(P_n, \varphi_n) \in W$ subject to (20). □

We have now all ingredients to establish the main result here, namely the following statement.

**Theorem 1.** Let us consider $t \le 0$, $P_0 \ne 0$, $\varkappa_{\parallel} \ge 0$, by assuming also that $\xi$ and $\varepsilon$ satisfy (9)-(12). Furthermore, let $\varphi_S \in H^{1/2}(S)$ and $w \in H^1(\Omega)$ such that $w|_S = \varphi_S$. Then, there exists a solution of the formulation (13)-(14).

*Proof.* The case where $\varkappa_{\parallel} = 0$ is immediate; it is even related to a determination of $(P, \varphi)$ satisfying (13)-(14), where $P$ can be expressed explicitly (as for instance in the frame of Remark 1 below). Let us set $\varkappa_{\parallel} > 0$ in the sequel. Let $(P_n, \varphi_n)_{n \in \mathbb{N}} \in W^N$ a sequence of solutions of (20), built from Proposition 2. We organize this proof in three parts. Firstly, we show that a subsequence $(P_{n_k}, \varphi_{n_k})_{k \in \mathbb{N}}$ converges weakly in $W$ and almost everywhere in $\Omega_f \times \Omega$ to an element $(P, \varphi)$ of $W$. Next, we establish that $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ converges in $L^1(\Omega_f)$ to $g \circ P$. Finally, we show that $(P, \varphi)$ is a solution of (13)-(14).

Before starting, let us note that, for all $(Q, \psi) \in W$, by adding the first equation of (20) to the second one multiplied by $\frac{\varkappa_{\parallel}}{16\pi^2}$, it follows that:

$$ 
\begin{aligned}
& \int_{\Omega_f} P_n Q dx - \int_{\Omega_f} (T_n \circ g \circ P_n) Q dx + \int_{\Omega_f} (\xi \nabla P_n) \cdot \nabla Q dx \\
& + \int_{\Omega} \frac{\varkappa_{\parallel}}{16\pi^2} (\varepsilon \nabla \varphi_n) \cdot \nabla \psi dx + \int_{\Omega_f} \frac{\varkappa_{\parallel}}{4\pi} (\partial_z \varphi_n) Q dx + \int_{\Omega_f} \frac{\varkappa_{\parallel}}{4\pi} (\partial_z P_n) \psi dx \\
& - \int_{S_f} \frac{\varkappa_{\parallel}}{4\pi} \nu_z P_n \psi d\sigma = - \int_{\Omega_f} \frac{\varkappa_{\parallel}}{4\pi} (\partial_z w) Q dx - \int_{\Omega} \frac{\varkappa_{\parallel}}{16\pi^2} (\varepsilon \nabla w) \cdot \nabla \psi dx.
\end{aligned}
\quad (21) $$

Then,

• since $g(x) \to \mp\infty$, when $x \to \pm\infty$, there exists $\beta > 0$ such that $g(x) \le 0$ for all $x \ge \beta$ and $g(x) \ge 0$ for all $x \le -\beta$. Let $n \in \mathbb{N}$. By choosing $(P_n, \varphi_n)$ as a test-function in (21), we have, thanks to the Green
formula, that:

$$
\begin{align}
& \int_{\Omega_f} P_n^2 \mathrm{d}x - \int_{\Omega_f} (T_n \circ g \circ P_n) P_n \chi_{\{|P_n| > \beta\}} \mathrm{d}x + \int_{\Omega_f} (\xi \nabla P_n) \cdot \nabla P_n \mathrm{d}x \nonumber \\
& + \int_{\Omega} \frac{\varkappa_{\|}}{16\pi^2} (\varepsilon \nabla \varphi_n) \cdot \nabla \varphi_n \mathrm{d}x = \int_{\Omega_f} (T_n \circ g \circ P_n) P_n \chi_{\{|P_n| \le \beta\}} \mathrm{d}x \nonumber \\
& - \int_{\Omega_f} \frac{\varkappa_{\|}}{4\pi} (\partial_z w) P_n \mathrm{d}x - \int_{\Omega} \frac{\varkappa_{\|}}{16\pi^2} (\varepsilon \nabla w) \cdot \nabla \varphi_n \mathrm{d}x. \tag{22}
\end{align}
$$

In (22), the left hand-side terms are all positive. Furthermore, with (16),

$$
\left| \int_{\Omega_f} (T_n \circ g \circ P_n) P_n \chi_{\{|P_n| \le \beta\}} \mathrm{d}x \right| \le \left( \left| \Omega_f \right| \sup_{x \in [-\beta, \beta]} (g(x))^2 \right)^{1/2} \|P_n\|_{L^2(\Omega_f)}.
$$

We deduce:

$$
\|(P_n, \varphi_n)\|_W^2 \le \frac{C}{\alpha} \|(P_n, \varphi_n)\|_W, \quad (23)
$$

with $C = (|\Omega_f| \sup_{x \in [-\beta, \beta]} (g(x)^2))^{1/2} + \frac{\varkappa_{\|}}{4\pi} \|\partial_z w\|_{L^2(\Omega_f)} + \frac{9\varkappa_{\|}}{16\pi^2} \max_{1 \le i,j \le 3} \|\varepsilon_{i,j}\|_{L^\infty(\Omega)} |w|_{H^1(\Omega)}$, where
$\alpha$ is an ellipticity constant of the bilinear map $a(\cdot, \cdot)$ resulting from the proof of Proposition 1. Thus
$(P_n, \varphi_n)_{n \in \mathbb{N}}$ is bounded in $W$ and we can extract a subsequence $(P_{n_m}, \varphi_{n_m})_{m \in \mathbb{N}}$ which converges weakly to $(P, \varphi)$ in $W$. In particular, $(P_{n_m})_{m \in \mathbb{N}}$ converges weakly to $P$ in $H^1(\Omega_f)$. Furthermore, as the injection of $H^1(\Omega_f)$ in $L^2(\Omega_f)$ is compact, the sequence $(P_{n_m})_{m \in \mathbb{N}}$ converges to $P$ in $L^2(\Omega_f)$. Thus, there exists (cf., e.g., [5]) a subsequence $(P_{n_{m_p}})_{p \in \mathbb{N}}$ which converges almost everywhere to $P$ in $\Omega_f$. In the same way, we can show that $(\varphi_{n_{m_p}})_{p \in \mathbb{N}}$ converges to $\varphi$ in $L^2(\Omega)$. It follows that there is a subsequence $(\varphi_{n_{m_{p_q}}})_{q \in \mathbb{N}}$ which converges to $\varphi$ almost everywhere in $\Omega$. To simplify the notation, we will denote by $(P_{n_k}, \varphi_{n_k})_{k \in \mathbb{N}}$ the built subsequence;

* we show now that the sequence $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ converges in $L^1(\Omega_f)$ to $g \circ P$. Since $(P_{n_k})_{k \in \mathbb{N}}$ converges to $P$ a.e. in $\Omega_f$, we deduce, by continuity of $g$ and thanks to (17), that $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ converges to $g(P)$ a.e. in $\Omega_f$. If we establish that $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ is uniformly integrable, we can conclude with the help of the Vitali convergence Theorem. Let us start by noting that $((T_{n_k} \circ g \circ P_{n_k}) P_{n_k})_{k \in \mathbb{N}}$ is bounded in $L^1(\Omega_f)$. By reproducing the previous calculations, we obtain from (22) that, for $k \in \mathbb{N}$,

$$
- \int_{\Omega_f} (T_{n_k} \circ g \circ P_{n_k}) P_{n_k} \chi_{\{|P_{n_k}| > \beta\}} \mathrm{d}x \le C \|(P_{n_k}, \varphi_{n_k})\|_W.
$$

Thus, with (23), it derives that:

$$
- \int_{\Omega_f} (T_{n_k} \circ g \circ P_{n_k}) P_{n_k} \chi_{\{|P_{n_k}| > \beta\}} \mathrm{d}x \le \frac{C^2}{\alpha}.
$$

Furthermore, with (16), it follows that:

$$
\int_{\Omega_f} |(T_{n_k} \circ g \circ P_{n_k}) P_{n_k}| \chi_{\{|P_{n_k}| \le \beta\}} \mathrm{d}x \le |\Omega_f| \sup_{x \in [-\beta, \beta]} |g(x)x|.
$$

Therefore, the sequence $((T_{n_k} \circ g \circ P_{n_k}) P_{n_k})_{k \in \mathbb{N}}$ is bounded in $L^1(\Omega_f)$, by the constant $C'$, $C' = \frac{c^2}{\alpha} + |\Omega_f| \sup_{x \in [-\beta, \beta]} |g(x)x|$. Let $\eta > 0$, $E \subset \Omega_f$ and $M > 0$. We have:

$$
\int_E |T_{n_k} \circ g \circ P_{n_k}| \mathrm{d}x = \int_{E \cap \{|P_{n_k}| > M\}} |T_{n_k} \circ g \circ P_{n_k}| \mathrm{d}x + \int_{E \cap \{|P_{n_k}| \le M\}} |T_{n_k} \circ g \circ P_{n_k}| \mathrm{d}x.
$$

On the one hand,

$$
\int_{E \cap \{|P_{n_k}| > M\}} |T_{n_k} \circ g \circ P_{n_k}| \mathrm{d}x \le \frac{1}{M} \int_{E \cap \{|P_{n_k}| > M\}} |(T_{n_k} \circ g \circ P_{n_k}) P_{n_k}| \mathrm{d}x \le \frac{C'}{M},
$$
and, by choosing in particular $M = \frac{2C'}{\eta}$, we obtain:

$$ \int_{E \cap \{|P_{n_k}| > M\}} |T_{n_k} \circ g \circ P_{n_k}| \, dx < \frac{\eta}{2}. $$

On the other hand, we have, by using (16), that:

$$ \int_{E \cap \{|P_{n_k}| \le M\}} |T_{n_k} \circ g \circ P_{n_k}| \, dx \le \int_{E \cap \{|P_{n_k}| \le M\}} |g \circ P_{n_k}| \, dx \le |E| \sup_{|x| \le M} |g(x)|. $$

So, for all $E \subset \Omega_f$ satisfying $|E| \le \frac{\eta}{2\sup_{|x| \le M} |g(x)|}$, it follows that:

$$ \int_E |T_{n_k} \circ g \circ P_{n_k}| \, dx < \eta. $$

Thus, $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ is uniformly integrable and the Vitali convergence Theorem allows to show that $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ converges in $L^1(\Omega_f)$ to $g \circ P$;

• in order to conclude, let us note that $(P, \varphi)$ satisfies (13)-(14). Let $(Q, \psi) \in C^\infty(\bar{\Omega}_f) \times \mathcal{D}(\Omega)$. At first, by using the Green formula, we have that:

$$ \left| \int_{\Omega_f} \partial_a (P_{n_k} - P) \partial_b Q dx \right| \le \left| \int_{\Omega_f} (P_{n_k} - P) \partial_a \partial_b Q dx \right| + \left| \int_{S_f} \nu_a (P_{n_k} - P) \partial_b Q d\sigma \right| \\ \le \|P_{n_k} - P\|_{L^2(\Omega_f)} \|\partial_a \partial_b Q\|_{L^2(\Omega_f)} + \|P_{n_k} - P\|_{L^2(S_f)} \|\nu_a \partial_b Q\|_{L^2(S_f)}, $$

where $a, b \in [1, 3]$, $(\partial_1, \partial_2, \partial_3) := (\partial_x, \partial_y, \partial_z)$, $(\nu_1, \nu_2, \nu_3)^T := (\nu_x, \nu_y, \nu_z)^T$. Since $(P_{n_k}, \varphi_{n_k})_{k \in \mathbb{N}}$ converges weakly to $(P, \varphi)$ in $W$, we obtain that $(P_{n_k})_{k \in \mathbb{N}}$ is bounded in $H^1(\Omega_f)$ by a constant $K_1$, $K_1 > 0$. Furthermore, as $S_f$ is the Lipschitz-continuous boundary of $\Omega_f$, there exists (cf., e.g., [4]) a constant $K_2$, $K_2 > 0$, such that:

$$ \|P_{n_k} - P\|_{L^2(S_f)} \le K_2 \|P_{n_k} - P\|_{H^1(\Omega_f)}^{1/2} \|P_{n_k} - P\|_{L^2(\Omega_f)}^{1/2}. $$

Thus, $\|P_{n_k} - P\|_{L^2(S_f)} \le K_2 \sqrt{K_1 + \|P\|_{H^1(\Omega_f)} \|P_{n_k} - P\|_{L^2(\Omega_f)}^{1/2}}$, and $(P_{n_k}|_{S_f})_{k \in \mathbb{N}}$ converges to $P|_{S_f}$ in $L^2(S_f)$. We deduce that: $\left| \int_{\Omega_f} \partial_a (P_{n_k} - P) \partial_b Q dx \right| \to 0$, when $k \to +\infty$, and therefore,

$$ \left| \int_{\Omega_f} (\xi \nabla (P_{n_k} - P)) \cdot \nabla Q dx \right| \to 0, \quad \text{when } k \to +\infty. $$

Similarly, we get that:

$$ \left| \int_{\Omega} (\varepsilon \nabla (\varphi_{n_k} - \varphi)) \cdot \nabla \psi dx \right| \to 0, \quad \text{when } k \to +\infty. $$

Furthermore, with the Green formula, we have,

$$ \left| \int_{\Omega_f} (\partial_z(\varphi_{n_k} - \varphi)) Q dx \right| \le \|\varphi_{n_k} - \varphi\|_{L^2(\Omega_f)} \|\partial_z Q\|_{L^2(\Omega_f)} + \|\varphi_{n_k} - \varphi\|_{L^2(S_f)} \|\nu_z Q\|_{L^2(S_f)}. $$

With the similar arguments as before, it follows that:

$$ \left| \int_{\Omega_f} (\partial_z(\varphi_{n_k} - \varphi)) Q dx \right| \to 0, \quad \text{when } k \to +\infty. $$

Moreover, with the Green formula, we get:

$$ \left| \int_{\Omega_f} (\partial_z(P_{n_k} - P)) \psi dx - \int_{S_f} \nu_z(P_{n_k} - P) \psi d\sigma \right| \le \|P_{n_k} - P\|_{L^2(\Omega_f)} \|\partial_z \psi\|_{L^2(\Omega_f)}, $$
and therefore,

$$ \left| \int_{\Omega_f} (\partial_z(P_{n_k} - P)) \psi dx - \int_{S_f} \nu_z(P_{n_k} - P) \psi d\sigma \right| \to 0, \quad \text{when } k \to +\infty. $$

In the same way, we note that:

$$ \left| \int_{\Omega_f} (P_{n_k} - P) Q dx \right| \to 0, \quad \text{when } k \to +\infty. $$

Finally, with the Hölder inequality, we obtain that:

$$ \left| \int_{\Omega_f} (T_{n_k} \circ g \circ P_{n_k} - g \circ P) Q dx \right| \le \|T_{n_k} \circ g \circ P_{n_k} - g \circ P\|_{L^1(\Omega_f)} \|Q\|_{L^\infty(\Omega_f)}. $$

Since it has been established that $(T_{n_k} \circ g \circ P_{n_k})_{k \in \mathbb{N}}$ converges to $g \circ P$ in $L^1(\Omega_f)$, it follows that:

$$ \left| \int_{\Omega_f} (T_{n_k} \circ g \circ P_{n_k} - g \circ P) Q dx \right| \to 0, \quad \text{when } k \to +\infty. $$

By taking the limit in each term of the weak formulation (21), which is namely satisfied by each element of the sequence $(P_{n_k}, \varphi_{n_k})_{k \in \mathbb{N}}$, we get, for all $(Q, \psi) \in C^\infty(\bar{\Omega}_f) \times \mathcal{D}(\Omega)$,

$$ \int_{\Omega_f} P Q dx - \int_{\Omega_f} (g \circ P) Q dx + \int_{\Omega_f} (\xi \nabla P) \cdot \nabla Q dx + \int_{\Omega} \frac{\kappa_{\|}}{16\pi^2} (\varepsilon \nabla \varphi) \cdot \nabla \psi dx + \int_{\Omega_f} \frac{\kappa_{\|}}{4\pi} (\partial_z \varphi) Q dx \\ + \int_{\Omega_f} \frac{\kappa_{\|}}{4\pi} (\partial_z P) \psi dx - \int_{S_f} \frac{\kappa_{\|}}{4\pi} \nu_z P \psi d\sigma = - \int_{\Omega_f} \frac{\kappa_{\|}}{4\pi} (\partial_z w) Q dx - \int_{\Omega} \frac{\kappa_{\|}}{16\pi^2} (\varepsilon \nabla w) \cdot \nabla \psi dx. $$

In particular, by choosing $\psi = 0$ and $Q = 0$ successively, in the above relation, we deduce that:

$$ \left\{ \begin{aligned} & \int_{\Omega_f} t P Q dx + \int_{\Omega_f} \frac{P^3}{P_0^2} Q dx + \int_{\Omega_f} (\xi \nabla P) \cdot \nabla Q dx + \int_{\Omega_f} \frac{\kappa_{\|}}{4\pi} (\partial_z \varphi) Q dx = - \int_{\Omega_f} \frac{\kappa_{\|}}{4\pi} (\partial_z w) Q dx, \\ & \qquad \int_{\Omega} (\varepsilon \nabla \varphi) \cdot \nabla \psi dx - \int_{S_f} 4\pi \nu_z P \psi d\sigma + \int_{\Omega_f} 4\pi (\partial_z P) \psi dx = - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla \psi dx. \end{aligned} \right. $$

Finally, the relations of the obtained system are valid for $(Q, \psi) \in W$, due to the density of $C^\infty(\bar{\Omega}_f)$ in $H^1(\Omega_f)$ and to the one of $\mathcal{D}(\Omega)$ in $H_0^1(\Omega)$. Then, $(P, \varphi)$ is a solution of (13)-(14).

□

Let us note that if $(P, \varphi)$ verifies the weak formulation (13)-(14), then it follows that $(P, \varphi + w)$ is a weak solution of (EGL), and the (weak) electric polarization $P$ is determined from (1), wherein $E = -\nabla(\varphi + w)$.

As indicates the following remark, the uniqueness of solution of the weak formulation is not ensured for any choice of the parameters.

**Remark 1.** Let us mention that there is not a unique solution for any choice of parameters of the system (EGL). For instance, if $t < 0$, $P_0 \neq 0$, $\kappa_{\|} = 0$, and $\varphi_S = 0$ (as well as $w = 0$), then there are two distinct pairs $(P, \varphi)$, $(-P, -\varphi)$, with $P := \sqrt{|t|} |P_0|$, satisfying the weak formulation (13)-(14).

# 4 Finite Element approximations and discrete nonlinear system

We discretize the formulation (13)-(14) by making use of a mesh of $\bar{\Omega}$ and of the Lagrange finite elements. Such a mesh consists of a collection of tetrahedra, $T_h$, obtained from a usual process of triangulation, and is fine enough near the boundaries of $\Omega_f$ and $\Omega$ that are here polyhedral. In a correlative way, the discrete region associated with the boundary of $\Omega_f$ is entirely made up of faces of tetrahedra. Each of these faces is
common both to a tetrahedron of the discrete region associated with $\overline{\Omega_f}$ and to another one of the discrete region associated with $\overline{\Omega_p}$.

Let us introduce the discrete space $W_h$,

$$W_h = X_h \times Y_h,$$

where,

$X_h = \{P_h \in C^0(\overline{\Omega_f}) ; P_h|_T \in P_{1,T} \forall T \in \mathcal{T}_h, T \subset \overline{\Omega_f}\}, \quad Y_h = \{\varphi_h \in C^0(\overline{\Omega}) ; \varphi_h|_T \in P_{1,T} \forall T \in \mathcal{T}_h, \varphi_h|_S = 0\}$, and $P_{1,T}$ is the space of polynomials of degree less than or equal to 1, defined on $T$.

The discrete formulation associated with (13)-(14) consists of finding $(P_h, \varphi_h) \in W_h$ such that, for all $(Q_h, \psi_h) \in W_h$,

$$ 
\left\{ 
\begin{aligned}
& \int_{\Omega_f} t P_h Q_h dx + \int_{\Omega_f} \frac{P_h^3}{P_0^2} Q_h dx + \int_{\Omega_f} (\xi \nabla P_h) \cdot \nabla Q_h dx + \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z \varphi_h) Q_h dx \\
& = - \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z w) Q_h dx, \\
& \int_{\Omega} (\varepsilon \nabla \varphi_h) \cdot \nabla \psi_h dx - \int_{S_f} 4\pi \nu_z P_h \psi_h d\sigma + \int_{\Omega_f} 4\pi (\partial_z P_h) \psi_h dx \\
& = - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla \psi_h dx.
\end{aligned}
\right. 
\qquad (24) $$

We are thus in particular concerned here with conforming discretizations of the regions $\overline{\Omega_f}$, $\overline{\Omega}$, as well as with a conforming discretization of the space $W$; namely, $W_h \subset W$. By considering the data as in the assumptions of Theorem 1, it derives the existence of solutions of the discrete formulation. We already mention that each solution $(P_h, \varphi_h)$ of this formulation leads to the determination of a discrete electric polarization associated with $\mathbf{P}$, by using (1), wherein $P$ is replaced by $P_h$ and $\mathbf{E}$ by $-\nabla(\varphi_h + w)$.

Let us denote by $\{x_j, 1 \le j \le N\}$ the set of vertices of $\mathcal{T}_h$, numbered in such a way that $x_1, \dots, x_{N_f} \in \overline{\Omega_f}$ and $x_1, \dots, x_{N^0} \in \Omega$, where $N_f, N^0$ and $N$ are respectively the number of vertices in $\overline{\Omega_f}$, $\Omega$ and $\overline{\Omega}$; namely, $N_f < N^0 < N$. Let $v_j$ be the shape function associated with $x_j$, for $j \in [1, N]$. We have thus, for $(P_h, \varphi_h) \in W_h$,

$$P_h = \sum_{j=1}^{N_f} P_j v_j, \quad \varphi_h = \sum_{j=1}^{N^0} \varphi_j v_j,$$

where $P_j \in \mathbf{R}$ for all $j \in [1, N_f]$ and $\varphi_j \in \mathbf{R}$ for all $j \in [1, N^0]$. The discrete formulation (24) can then be written as a nonlinear matrix system having two vector unknowns. It consists of finding $\mathcal{P} = (P_1, \dots, P_{N_f})^T \in \mathbf{R}^{N_f}$ and $\Phi = (\varphi_1, \dots, \varphi_{N^0})^T \in \mathbf{R}^{N^0}$ satisfying:

$$ 
\left\{ 
\begin{array}{l}
\mathcal{A}\mathcal{P} + \mathcal{N}(\mathcal{P}) + \mathcal{B}\Phi = \mathcal{F}_1, \\
\mathcal{C}\mathcal{P} + \mathcal{D}\Phi = \mathcal{F}_2,
\end{array}
\right. 
\qquad (25) $$

where,

$$A = \left( \int_{\Omega_f} (t v_j v_i + (\xi \nabla v_j) \cdot \nabla v_i) dx \right)_{(i,j) \in [1, N_f]^2},$$

$$B = \left( \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z v_j) v_i dx \right)_{(i,j) \in [1, N_f] \times [1, N^0]},$$

$$C = \left( - \int_{S_f} 4\pi \nu_z v_j v_i d\sigma + \int_{\Omega_f} 4\pi (\partial_z v_j) v_i dx \right)_{(i,j) \in [1, N^0] \times [1, N_f]},$$

$$D = \left( \int_{\Omega} (\varepsilon \nabla v_j) \cdot \nabla v_i dx \right)_{(i,j) \in [1, N^0]^2},$$

$$N(\mathcal{P}) = \left( \int_{\Omega_f} \frac{1}{P_0^2} \left( \sum_{j=1}^{N_f} P_j v_j \right)^3 v_i dx \right)_{i \in [1, N_f]},$$
$$F_1 = \left( - \int_{\Omega_f} \frac{\kappa_{\parallel}}{4\pi} (\partial_z w) v_i dx \right)_{i \in [1, N_f]} , \quad F_2 = \left( - \int_{\Omega} (\varepsilon \nabla w) \cdot \nabla v_i dx \right)_{i \in [1, N^0]} .$$

The terms of these blocks are computed with a 5-order quadrature formula (that can be found for instance in [26]). Note that, when the tensors $\varepsilon$ and $\xi$ are constant in the geometrical elements and if $w$ is in particular a polynomial of degree 1 (it will always be the case from the consideration of $\varphi_S$ explicitly done hereafter), this formula is exact regarding the computation of these terms.

In order to solve the nonlinear matrix system, we deal with an iterative algorithm (see [17] for more details) from the combination of two inexact Newton methods: the first one is globalized with a linesearch technique (see, e.g., [17, 19]), and converges slowly, for a large choice of initializations, whereas the second one is classical (see, e.g., [7]), converges faster but for a restricted choice of initializations. The GMRES algorithm (see [21]), using a preconditioner based on an incomplete LU factorization, was incorporated into each iteration.

# 5 Numerical simulations

## 5.1 Computational configurations

We consider three computational configurations. The meshes of $\bar{\Omega}$ associated with these configurations are achieved in accordance with the previous section, and represented here by $\mathcal{T}_h^{(i)}$, $1 \le i \le 3$. On the one hand, geometrically, $\Omega_f$ and $\Omega$ are polyhedral and shaped like cylinders; they share the same center, namely the origin, and the same axis, represented by $(Oz)$. We denote by $R_f$ the radius and $H_f$ the height of the “cylinder” represented by $\Omega_f$. In a similar way, the radius and the height of the “cylinder” represented by $\Omega$ are denoted by $R$ and $H$. On the other hand, the context where $\Omega_f$ and $\Omega$ are parallelepipedal is considered, in particular for the last configuration. Namely,

* for $\mathcal{T}_h^{(1)}$, the corresponding values are fixed such that $(R_f, H_f) = (0.4, 0.1)$ and $(R, H) = (1.2, 1.7)$;

* and for $\mathcal{T}_h^{(2)}$, the associated values are such that $(R_f, H_f) = (4, 1)$ and $(R, H) = (12, 17)$;

* whereas for $\mathcal{T}_h^{(3)}$, we set $\bar{\Omega}_f = [-4, 4] \times [-4, 4] \times [-0.5, 0.5]$ and $\bar{\Omega} = [-12, 12] \times [-12, 12] \times [-8.5, 8.5]$.

Let us already mention that, for each computational configuration, several values have been associated with the mesh size $h$ in the numerical simulations. Most of the results that will be presented in the next subsections derive from simulations based on $\mathcal{T}_h^{(1)}$ with $h = 0.0479$, as well on $\mathcal{T}_h^{(2)}$ with $h = 0.37$, as on $\mathcal{T}_h^{(3)}$ with $h = 0.36$, involving 1499982 degrees of freedom about $\mathcal{T}_h^{(1)}$, 1865401 in the frame of $\mathcal{T}_h^{(2)}$ and 2541899 in the context of $\mathcal{T}_h^{(3)}$; the associated mesh size, as regards $\bar{\Omega}_f$, being equal to 0.0228 about $\mathcal{T}_h^{(1)}$, 0.31 in the frame of $\mathcal{T}_h^{(2)}$ and 0.34 in the context of $\mathcal{T}_h^{(3)}$. More particularly, when we will be concerned with a description of results related to other choices for $h$, the values used there will be specified.

Regarding the physical parameters, different values are in particular associated with the coefficients of $\xi$ and $\varepsilon$ defined respectively in $\Omega_f$ and $\Omega$. In these considerations, the tensor $I_3$ denotes the 3 × 3 identity and $\varepsilon|_{\Omega_f}$ is represented as below,

$$\varepsilon|_{\Omega_f} = \begin{pmatrix} \varepsilon_{\perp x} & 0 & 0 \\ 0 & \varepsilon_{\perp y} & 0 \\ 0 & 0 & \varepsilon_{\parallel z} \end{pmatrix},$$

where $\varepsilon_{\perp x} > 0$, $\varepsilon_{\perp y} > 0$, $\varepsilon_{\parallel z} > 0$. Namely,

* for $\mathcal{T}_h^{(1)}$, the corresponding values are fixed such that $\varepsilon_{\perp x} = \varepsilon_{\perp y} = 10$, $\varepsilon_{\parallel z} = 0.1$, $\xi = 0.1I_3$;

* and for $\mathcal{T}_h^{(2)}$, the associated values are such that $\varepsilon_{\perp x} = 10$, $\varepsilon_{\perp y} = 5$, $\varepsilon_{\parallel z} = 1$, $\xi = I_3$;

* whereas for $\mathcal{T}_h^{(3)}$, the same values as in the context of $\mathcal{T}_h^{(2)}$ are used for $\varepsilon_{\perp x}$, $\varepsilon_{\perp y}$, $\varepsilon_{\parallel z}$ and $\xi$.

Furthermore, we set: $P_0 = 1$, $\varkappa_{\parallel} = 11.5$, $\varepsilon|_{\Omega_p} = I_3$, for each computational configuration. We define $\varphi_s$ as below,

$$(x, y, z) \in S \mapsto \varphi_S(x, y, z) = \frac{\gamma - z}{\gamma} \phi_0 - \phi_0 \in \mathbf{R},$$
in such a way that $\phi_S(x, y, -\gamma) = \phi_0$ and $\phi_S(x, y, \gamma) = -\phi_0$, by recalling that the lower and upper faces of $\bar{\Omega}$ are subparts of the planes $(z = -\gamma)$ and $(z = \gamma)$ respectively. Namely, $\gamma = 0.85$ in the frame of $\mathcal{T}_h^{(1)}$, and 8.5 in the context of $\mathcal{T}_h^{(2)}$ as well as of $\mathcal{T}_h^{(3)}$. Let us denote by $U$ the voltage potential applied to the device,

$$U = 2\phi_0,$$

in association with the use in (13)-(14) and (24) of the function $w$ such that $w(x, y, z) = \frac{\gamma - z}{\gamma}\phi_0 - \phi_0$, for $(x, y, z)^T \in \Omega$. Namely, $w|_S = \phi_S$.

In the next subsections, we will describe the numerical results, in each configuration, in particular through the behavior of the discrete "average polarization", which is represented by $\bar{P}_h$,

$$\bar{P}_h = \frac{1}{|\Omega_f|} \int_{\Omega_f} P_h dx,$$

in correspondence with $(P_h, \phi_h)$; any state $(P, \phi)$, subject to (EGL), having a finite element approximation denoted here by $(P_h, \phi_h)$ in accordance with Section 4.

As it was indicated in Section 2, although, in the applications, the reduced temperature $t$ must take values in $[-1, 0]$, we will be interested in a more large context of values for this parameter, namely by investigating with $t \in (-\infty, 0]$.

Let us already mention that other values for $P_0$ and $\kappa_\parallel$ were considered and the numerical results, obtained with these, were qualitatively similar to the ones that will be described hereafter. We also indicate that comparisons with analytic solutions including the case of data mentioned in Remark 1 were done from simulations, by retrieving successfully also numerical convergences to these solutions.

## 5.2 Numerical investigations of states

We deal here with a protocol devoted to the determination of non obvious states, and which incorporates also a process of heating as well as of cooling of the device. Several challenges face us. At first, we need to get different non obvious ferroelectric states, and this leads us to study how to initialize the iterative algorithm. Although having a globalization step, this algorithm however remains sensitive to the consideration of initial data. The second difficulty is that we must find a way to "heat up" as to "cool down" numerically these different states. We initialize our iterative algorithm with random vectors of size $N_f + N^0$, distributed according to the uniform law of $[-1, 1]^{N_f+N^0}$, as well as with deterministic vectors.

We obtain several profiles; each one allowing us to identify the (discrete) states $(P_h, \phi_h)$ subject to (EGL) for which the polarization structures are similar. Figure 2 presents some of these profiles, resulting from $\mathcal{T}_h^{(2)}$ and $\mathcal{T}_h^{(3)}$, with $t = -10$, $U = 0$, and distinguished by the shape and the number of polarization domains:

* the Monodomain Profile (see, e.g., (MP) in Figure 2(a)-(b)), case of a unique domain of polarization;

* the Concentric-cylindrical Profile (see, e.g., (CP) in Figures 2(c)-(d)), context of a central domain shaped like a cylinder or distorted cylinder, surrounded by another domain;

* the Band Profiles (see, e.g., (4BP) and (3BP) in Figures 2(e) and (f)), case of several domains, organized into bands.

These profiles arise also from simulations based on $\mathcal{T}_h^{(1)}$, namely for instance with $t = -20$. In addition to these, the Obvious Profile, denoted by (OP) and corresponding to the case where $P_h \equiv 0$, occurs from simulations performed here. Some of the obtained profiles have already been observed in other works (see, e.g., [16, 25]); for instance, a profile analogous to (3BP) as well as to (4BP), namely the 2-bands profile, has been found experimentally in the case of a nanofilm of Lead Titanate (see [25]).

Now, we specify how to achieve numerically the heating or cooling of the device. Let $(P^{(0)}, \phi^{(0)})$ be a non obvious state at the reduced temperature $t_0$, associated with the discrete representation $(\mathcal{P}^{(0)}, \Phi^{(0)})^T \in \mathbb{R}^{N_f+N^0}$ resulting from the nonlinear system (25). At the reduced temperature $t = t_0 + \delta t$ (where $\delta t > 0$ for heating and $\delta t < 0$ for cooling), we initialize with this vector our algorithm for solving (25), and obtain a new solution, $(\mathcal{P}^{(1)}, \Phi^{(1)})^T$, of this system. By proceeding in this way, step by step, we get a sequence $((P^{(k)}, \phi^{(k)})^)_k$ of states, namely a sequence $((\mathcal{P}^{(k)}, \Phi^{(k)})^)_k$ of discrete representations, where the solution $(\mathcal{P}^{(k)}, \Phi^{(k)})^T$ is provided at the reduced temperature $t = t_0 + k\delta t$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Profiles obtained from $\mathcal{T}_h^{(2)}$ (see (a),(c),(e)) and $\mathcal{T}_h^{(3)}$ (see (b),(d),(f)), with $t = -10$. At left, a cross-section of $P_h$ following ($z=0$) is represented, whereas at middle and at right, the cross-sections of the discrete electric polarization, associated with $\mathbf{P}$, following ($y=0$) and ($z=0$) are respectively represented.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

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

<figure><img src="image_6.png" /><figcaption>Figure 3: Variations of $\bar{P}_h$ with respect to the reduced temperature, from $\mathcal{T}_h^{(1)}$ (see (a)), as well as from $\mathcal{T}_h^{(2)}$ (see (b)) and $\mathcal{T}_h^{(3)}$ (see (c)).</figcaption></figure>
Figure 3(a) presents from $\mathcal{T}_h^{(1)}$ the variations of the discrete “average polarization” marked by solid lines, when we “heat up” (with $\delta t = 0.1$) or “cool down” (with $\delta t = -0.1$) several states, for $U = 0$; namely for instance (3BP). While $\bar{P}_h$ depends continuously on $t$, the distributions of the polarization stay similar. A discontinuity corresponds to a change of profile. For instance, (3BP) exists for a range of small values taken by the reduced temperature and, through heating, it jumps to the obvious profile at $t = -16.9$.

This behavior is characteristic of the ferroelectric devices (see, e.g., [18]), and appears also in the case of simulations from $\mathcal{T}_h^{(2)}$, as indicates Figure 3(b). For instance, the jumps schematized as below occur from $\mathcal{T}_h^{(2)}$:

$$ (\text{MP}) \underset{t=-13.5}{\longrightarrow} (\text{CP}) \underset{t=-6.4}{\longrightarrow} (\text{OP}), \quad (\text{4BP}) \underset{t=-4.8}{\longrightarrow} (\text{OP}). $$

As indicates Figure 3(c), the behavior of $\bar{P}_h$ with respect to the reduced temperature, in the context of $\mathcal{T}_h^{(3)}$, is similar to the one observed from the cylindrical configurations. In particular, similar jumps, schematized as below,

$$ (\text{MP}) \underset{t=-13.7}{\longrightarrow} (\text{CP}) \underset{t=-7.6}{\longrightarrow} (\text{4BP}) \underset{t=-5.1}{\longrightarrow} (\text{OP}), $$

occur from $\mathcal{T}_h^{(3)}$.

The most enduring profile with nonzero polarization is for $\mathcal{T}_h^{(1)}$ distinguished by (3BP), whereas for $\mathcal{T}_h^{(2)}$, as well as for $\mathcal{T}_h^{(3)}$, it is the one corresponding to (4BP). Each of the states associated with (4BP) has a constant “average polarization”; namely, $\bar{P}_h = 0$.

At each temperature, we identified the most stable state by evaluating the energy $\mathcal{E}$ for each of the obtained states. For instance in the context of $\mathcal{T}_h^{(3)}$, the states related to (4BP) remain the most stable on cooling, namely till $t = -7.7$ and then, the avail in energy passes to the ones involving the concentric-cylindrical profile (CP). This temperature value is slightly lower than the one of superheating presented above and in correlation with (CP), namely $t = -7.6$. Finally, the states related to (CP) remain the most stable for the considerations of $t$ such that $-15 \le t \le -7.8$.

Let us specify, as mentioned otherwise in Section 2, that the ferroelectric layer must be, in real applications, sufficiently small when compared to the medium $\Omega$, in order to avoid to disturb the emergent fringing field (see, e.g., [16] and the references therein). It is then particularly suited to deal with devices that take into account this aspect; moreover, numerical instabilities will not be involved as long as the volume of the ferroelectric layer does not appear as too small (although the asymptotic framework developed in [1] may remedy the contrary situation). We already mention that we deal with such a device with a ferroelectric layer for which the volume does not appear as too small, in each of the three computational configurations considered here. For instance in the context of $\mathcal{T}_h^{(3)}$, as for an inspection, let us embed successively the associated ferroelectric layer in two larger paraelectric environments, by defining thus: $\Omega^{(1)} = \bar{\Omega} = [-12, 12] \times [-12, 12] \times [-8.5, 8.5]$, $\Omega^{(2)} = [-15, 15] \times [-15, 15] \times [-11, 11]$, $\Omega^{(3)} = [-18, 18] \times [-18, 18] \times [-13.5, 13.5]$. Let us consider, in accordance with the previous section also, three meshes corresponding to triangulations of $\Omega^{(1)}$, $\Omega^{(2)}$ and $\Omega^{(3)}$ respectively, and having similar sizes approximately equal to 0.69, with correlative, as regards the ferroelectric layer, similar mesh sizes approximately equal to 0.60. A unique ferroelectric/paraelectric interface is here involved of course. Let us represent by $\{x_k^{(i)}\}_{1 \le k \le N^{(i)}}$ the set of internal vertices of the mesh of $\Omega^{(i)}$, where $i$ is fixed, $1 \le i \le 3$. For $u: \bar{\Omega} \to \mathbb{R}$ and $1 \le i \le 3$, let $\|u\|_2^{(i)} = (\frac{1}{N^{(i)}} \sum_{k=1}^{N^{(i)}} |u(x_k^{(i)})|^2)^{1/2}$. For a state $(P^{(i)}, \varphi^{(i)})$ subject to (EGL) and obtained by involving $\Omega^{(i)}$, with $\varphi_S = 0$, we consider the following error, $err_2^{(i)} = \frac{\|\varphi^{(3)} - \varphi^{(i)}\|_2^{(i)}}{\|\varphi^{(3)}\|_2^{(i)}}$, evaluated relatively with respect to the state $(P^{(3)}, \varphi^{(3)})$ which is a priori the more accurate one. By inspecting the behavior of the profile (MP), with respect to the reduced temperature, namely from $t = -15$ and with $\delta t = 0.1$, it follows that the jumps to (CP), (4BP) and (OP) occur again at $t = -13.7, -7.6$ and $-5.1$ respectively, for each of these meshes. The errors $err_2^{(i)}$ ($i = 1, 2$) associated with the profile related to (MP), or to (CP), or again to (4BP), as well to (OP), are, in percent, less than 6%. The similar evaluations of errors about $P^{(i)}$ provide values which, in percent, remain also smaller than 6%. So the paraelectric environment related to $\bar{\Omega}$ appears numerically as sufficiently large when compared to the ferroelectric layer; the influence of its extension being not particularly significant. Let us mention that the CPU time corresponding to one simulation (computation of $(P^{(i)}, \varphi^{(i)})$ obtained from $\Omega^{(i)}$, $1 \le i \le 3$) is approximately equal to 15mn, 1h30 or 3h on an “Intel Xeon Processor E5520 with a frequency of 2.27 GHz”, in the frame of $\Omega^{(1)}, \Omega^{(2)}$ or $\Omega^{(3)}$ respectively.
<figure><img src="image_7.png" /><figcaption>Figure 4: Variations of $\beta(P_h)$ with respect to $U$.</figcaption></figure>

## 5.3 Numerical investigations of hysteresis loops

The protocol we adopt for the numerical study of the existence of hysteresis loops is similar to the one that was defined earlier when we investigated the behavior of the polarization with respect to the reduced temperature $t$. Suppose that $(\mathcal{P}^{(0)}, \Phi^{(0)})^T$ is the discrete representation of the state $(P^{(0)}, \varphi^{(0)})$ subject to (EGL) and associated with the voltage potential $U = U^{(0)}$. For each fixed value $k$ ($k \in \mathbf{N}$), we denote by $(\mathcal{P}^{(k+1)}, \Phi^{(k+1)})^T$ the discrete representation that derives from (25) and corresponds to the voltage potential $U = U^{(0)} + (k+1)\delta U$, $\delta U \in \mathbf{R}^*$, after initializing our iterative algorithm with $(\mathcal{P}^{(k)}, \Phi^{(k)})^T$, the discrete representation of the state $(P^{(k)}, \varphi^{(k)})$ subject to (EGL). In this way, we deal with a sequence $((\mathcal{P}^{(k)}, \Phi^{(k)})^T)_k$ of discrete representations of states $(P^{(k)}, \varphi^{(k)})_k$ subject to (EGL). Namely, when $U^{(0)} = 0$, each of the discrete representations of the states discovered in the previous subsection can be used to initialize our iterative algorithm. Another alternative for the consideration of $U^{(0)}$ consists of choosing $U^{(0)}$ sufficiently large in order to deal with only one state; namely, based on physical arguments (see, e.g., [20]), this state must be quasi-uniform. Hence, we introduce a new parameter $\beta$,

$$ \beta(P_h) = \frac{\frac{1}{|\Omega_f|} \int_{\Omega_f} P_h^4 \mathrm{d}x}{\left(\frac{1}{|\Omega_f|} \int_{\Omega_f} P_h^2 \mathrm{d}x\right)^2}, $$

which reports in fact that $\beta(P_h) \approx 1$ when the state is quasi-uniform. For instance, as indicated in Figure 4, from simulations based on $\mathcal{T}_h^{(1)}$, with $t = -20$, we obtain $|\beta(P_h) - 1| < 10^{-1}$ for $|U| \ge 150$ and so we can consider $U^{(0)} = 150$ as the starting voltage potential for an investigation of hysteresis loops. Typically, the variation of states is inspected by increasing $(\delta U > 0)$ the voltage potential and by decreasing it $(\delta U < 0)$.

Figure 5 presents the numerical results of such an investigation based on each of the three computational configurations and involving different considerations of the reduced temperature. These results correspond to hysteresis loops. In the context of $\mathcal{T}_h^{(1)}$, it has been considered $t = -20$ and $\delta U = \pm 0.1$. For $U = -20$, there exists a unique profile, namely (CP). By increasing the voltage potential, (CP) jumps to (3BP) abruptly when $U = -1.3$. Then, (3BP) jumps for $U = 14.9$ to another profile, (CP), where the polarization domains are inverted. These two transitions allow to describe a branch of the hysteresis loop resulting here. The other branch of this loop occurs in a symmetric way, with the jumps schematized as below,

$$ (\mathrm{CP}) \xrightarrow[U=1.3]{} (\mathrm{3BP}) \xrightarrow[U=-14.9]{} (\mathrm{CP}). $$
<figure><img src="image_8.png" /></figure>

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

<figure><img src="image_10.png" /><figcaption>Figure 5: Hysteresis loops from $\mathcal{T}_h^{(1)}$ with $t = -20$ (see (a)), from $\mathcal{T}_h^{(2)}$ with $t = -15$ (—) and $t = -10$ (---) (see (b)), as well as from $\mathcal{T}_h^{(3)}$ with $t = -10$ (see (c)).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 6: Semi-log representation of the relative errors $err_2^{(i)}$ with respect to $U$, from $\mathcal{T}_{h_i}^{(1)}$ ($i=1,2$).</figcaption></figure>

These abrupt switchings are irreversible: for instance, in the previous case, the profile (3BP) does not jump to (CP) by a decrease of the voltage potential from $U = -1.3$. This leads, in particular, to the existence of several distinct states when $U = 0$; a phenomenon suggesting the applications for which the multibit storage would represent the main tool. Also, let us mention that each jump, here and in the following results, corresponds to a decrease of the energy (2), which is expected from Physics (see, e.g., [17]).

Let us indicate that the numerical results obtained here do not deteriorate from the considerations of some larger mesh sizes. As for an inspection, let us deal for instance with the context of the first computational configuration for which the corresponding mesh, $\mathcal{T}_h^{(1)}$, is also denoted here by $\mathcal{T}_{h_3}^{(1)}$, with $h_3 \equiv h = 0.0479$, recalling that the mesh size of the embedded layer was, correlatively, equal to 0.0228. By making use of two other meshes, that are coarser, we thus distinguish here three meshes of the same polyhedral region $\bar{\Omega}$, containing the same embedded polyhedral layer, also shaped like a cylinder. Namely, in addition to $\mathcal{T}_{h_3}^{(1)}$, we have a coarser discretization of $\bar{\Omega}$, denoted by $\mathcal{T}_{h_2}^{(1)}$, where the mesh size is $h_2 = 0.0872$ while for the embedded layer, the mesh size is, correlatively, equal to 0.0444. We also have a more coarser discretization of $\bar{\Omega}$, denoted by $\mathcal{T}_{h_1}^{(1)}$, where the mesh size is $h_1 = 0.166$, while for the embedded layer, the mesh size is, correlatively, equal to 0.0731. A unique ferroelectric/paraelectric interface is here involved of course. For $1 \le i \le 3$ and $u : \bar{\Omega} \to \mathbf{R}$, let us set: $\|u\|_2^{(i)} = (\frac{1}{N^{0(i)}} \sum_{k=1}^{N^{0(i)}} |u(x_k^{(i)})|^2)^{1/2}$, where $N^{0(i)}$ is the number of internal vertices of $\mathcal{T}_{h_i}^{(1)}$, which are represented here by $x_k^{(i)}$, $1 \le k \le N^{0(i)}$. For a state $(P^{(i)}, \varphi^{(i)})$, subject to (EGL) and resulting from $\mathcal{T}_{h_i}^{(1)}$, $1 \le i \le 3$, let us consider the following error, $err_2^{(i)} = \frac{\|\varphi^{(3)} - \varphi^{(i)}\|_2^{(i)}}{\|\varphi^{(3)}\|_2^{(i)}}$, for $i = 1,2$, evaluated relatively with respect to the state $(P^{(3)}, \varphi^{(3)})$. Figure 6 shows the behavior of $err_2^{(i)}$, with respect to the voltage potential from $U = -20$, with $\delta U = 1$ and for $t = -15$, $i = 1,2$. These relative errors are, in percent, less than 2%. The similar evaluations of relative errors about $P^{(i)}$ provide values which, in percent, remain also smaller than 2%. The switchings of states that occur here arise from the same values of $U$. Similar observations arise from the same inspection based on the two other computational configurations $\mathcal{T}_h^{(2)}$ and $\mathcal{T}_h^{(3)}$.

As indicates Figure 5(b), in the context of $\mathcal{T}_h^{(2)}$, by considering $t = -15$, with $\delta U = 1$, the jumps, schematized as below, occur for instance,

$$ (MP) \xrightarrow{U=26} (CP) \xrightarrow{U=49} (MP), $$

and for $t = -10$, with $\delta U = 1$, we obtain abrupt transitions again, represented as below,

$$ (MP) \xrightarrow{U=-39} (CP) \xrightarrow{U=27} (3BP) \xrightarrow{U=32} (MP). $$
As it can be observed from Figure 5(c), in the context of $T_h^{(3)}$, by considering $t = -10$, with $\delta U = 1$, we obtain a very similar behavior of the states, with particular jumps schematized below,

$$ (MP) \xrightarrow[\ U=-30 \ ]{ } (CP) \xrightarrow[\ U=16 \ ]{ } (3BP) \xrightarrow[\ U=35 \ ]{ } (CP) \xrightarrow[\ U=61 \ ]{ } (MP). $$

# 6 Concluding remarks

We have introduced a model based on the Ginzburg-Landau formalism and Electrostatics equations for the study of 3D devices made up of a ferroelectric layer fully embedded in a paraelectric environment, and involving a prescription of boundary voltage potentials. The existence of weak solutions of the model has been established from a weak formulation for which the unknowns are related to the polarization field and to the electric field.

In view of numerical simulations, we have discretized the formulation with the help of finite elements and dealt with inexact Newton techniques for solving the associated discrete nonlinear system. More specifically, we have developed two numerical protocols for our investigations. The first one has been devoted to the determination of non obvious states related to the system, by incorporating also a heating/cooling process of devices, in terms of the temperature, while the second one has been devoted to the study of the existence of hysteresis loops.

We have been interested in three computational configurations, with various considerations of geometrical and physical parameters. In particular, the contexts of devices having a cylindrical or parallelepipedic geometry, as the fully embedded ferroelectric layer, have been considered. In each computational configuration, similar polarization profiles have been obtained; namely, a monodomain profile, a concentric-cylindrical profile and several bands profiles. The same behavior of the states, with respect to the temperature, was founded; each state jumps into another at a certain temperature. Moreover, in each configuration, hysteresis loops, with abrupt transitions of states, occur. Remarkably, the same types of transitions have been obtained, as for instance the one reported by the transition of the concentric-cylindrical profile to the 3-bands profile.

As an extension of this work, we will be concerned with investigations in the context of 3D devices subject again to a prescription of boundary voltage potentials, but made up of a ferroelectric layer which is only partly embedded in a paraelectric environment. This context, also related to applications, requires of course a different model but which can be introduced as the present one from the Ginzburg-Landau formalism and Electrostatics equations. Beyond the investigations of non obvious states to perform in this context, one of the aims shall concern the numerical study of the so-called *effective permittivity* (cf, e.g., [3]), by associating also concrete values with parameters in the simulations.

# References

[1] H. Ammari, L. Halpern and K. Hamdache, Asymptotic behaviour of thin ferromagnetic films, *Asymptot. Anal.*, **24** (2000), 277-294.

[2] H. Ammari and K. Hamdache, Global existence and regularity of solutions to a system of nonlinear Maxwell equations, *J. Math. Anal. Appl.*, **286** (2003), 51-63.

[3] A.M. Bratkovsky and A.P. Levanyuk, Depolarizing field and “real” hysteresis loops in nanometer-scale ferroelectric films, *App. Phys. Lett.*, **89** (2006), 253108.

[4] S.C. Brenner and L.R. Scott, The mathematical theory of finite element methods, Springer-Verlag, New York, 2nd edition, 2002.

[5] H. Brézis, Analyse fonctionnelle: théorie et applications, Dunod, Paris, 1999.

[6] E.V. Chenskii and V.V. Tarasenko, Theory of phase transitions into inhomogeneous state in organic ferro-electrics in an external electric field, *Zh. Eksp. Teor. Fiz.*, **83** (1982), 1089. [Sov. Phys. JETP, **56** (1982), 618.]

[7] R.S. Dembo, S.C. Eisenstat and T. Steihaug, Inexact Newton methods, *SIAM Journal on Numerical Analysis*, **19** (1982), 400-408.

[8] J. Droniou, Degrés topologiques et applications, *Notes de Cours*, Université de Montpellier, http://www-gm3.univ-mrs.fr/polys/gm3-06/gm3-06.pdf, 2006.
[9] G.B. Folland, Real analysis, John Wiley & Sons, Inc., New York, 2nd edition, 1999.

[10] D. Gilbarg and N.S. Trudinger, Elliptic partial differential equations of second order, Grundlehren der mathematischen Wissenschaften, Springer, 2nd edition, 1998.

[11] V.L. Ginzburg, Polarization and piezoelectric effect in BaTiO3 near the ferroelectric transition point, Zh. Eksp. Teor. Fiz., **19** (1949), 36.

[12] P. Grisvard, Elliptic problems in nonsmooth domains, Pitman Advanced Pub. Program Boston, 1985.

[13] B.K. Lai, I. Ponomareva, I.I. Naumov et al., Electric-field-induced domain evolution in ferroelectric ultra-thin films, Phys. Rev. Lett., **96** (2006), 137602.

[14] L.D. Landau and E.M. Lifshitz, Course of theoretical physics - Electrodynamics of continuous media, volume 8, Elsevier, New-York, 1985.

[15] I.A. Luk'yanchuk, L. Lahoche and A. Sené, Universal properties of ferroelectric domains, Phys. Rev. Lett. **102** (2009).

[16] P.-W. Martelli, S.M. Mefire and I.A. Luk'yanchuk, Multidomain switching in the ferroelectric nanodots, Europhysics Lett. (EPL), **111** (2015), 50001.

[17] P.-W. Martelli, Modélisation et simulations numériques de la formation de domaines ferroélectriques dans des nanostructures 3D, Ph.D. Thesis, Université de Lorraine, 2016.

[18] V.N. Nechaev and A.V. Viskovatykh, On the change in the mechanism of the phase transition in a ferroelectric inclusion of a ferroelectric-dielectric nanocomposite as a function of its size, Phys. Solid State, **57** (2015), 722.

[19] J. Nocedal and S.J. Wright, Numerical optimization, Springer, New York, 2nd edition, 2006.

[20] K. Rabe, C.H. Ahn and J.-M. Triscone, Physics of ferroelectrics, a modern perspective, Springer-Verlag, Berlin, 2007.

[21] Y. Saad and M.H. Schultz, GMRES: A Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput., **7** (1986), 856-869.

[22] J.F. Scott, Ferroelectric memories, Advanced microelectronics. Springer, 2000.

[23] A. Sené, Théorie des domaines et des textures non uniformes dans les ferroélectriques, Ph.D. Thesis, Université de Picardie, 2010.

[24] M. Sepliarsky and R.E. Cohen, First-principles based atomistic modeling of phase stability in PMN-xPT, Journal of Physics: Condensed Matter, **23** (2011).

[25] S.K. Streiffer, J.A. Eastman, D.D. Fong, C. Thompson, A. Munkholm, M.V. Ramana Murty, O. Auciello, G.R. Bai and G.B. Stephenson, Observation of nanoscale 180 stripe domains in ferroelectric PbTiO3 thin films, Phys. Rev. Lett., **89** (2002).

[26] A.H. Stroud, Fifth degree integration formula for the n-simplex, SIAM J. Numer Anal., **6** (1969), 90-98.