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

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.


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 firstprincipia 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.

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 Ω of R 3 , and its fully embedded layer by an open subset Ω f , Ω f ⊂ Ω, as well as its paraelectric environment by Ω p = Ω \ Ω f . Also, let us denote by S, S = ∂Ω, the boundary of Ω, and by S f , S f = Ω f ∩ Ω p , the interface between Ω f and Ω p . The region Ω is also such that two subparts of its boundary are included in two planes, (z = −γ) and (z = γ), 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, P, on the electric field is nonlinear through only one of its components (here, the third one). More precisely, if we represent by 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: with P a scalar field, P : Ω f −→ R, and ε the electric permittivity.  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, E, is such that: whereF(P, E) is expressed as follows, with F(P, 0) the Ginzburg-Landau potential. The term "− Ω 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 "− 1 8π Ω E · (ε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 F(P, 0) can be written with respect to powers of P ,

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; • κ , the displacive parameter (namely positive); • P 0 , P 0 = 0, the spontaneous polarization at low temperatures; • ξ, the coherence length tensor.
In this expansion, the odd powers of P do not occur because of symmetry reasons. The term "∇P · (ξ∇P )" penalizes spatial inhomogeneities: it becomes important when ferroelectric domains are considered. Furthermore, if Ω is simply connected, the electric field E can be expressed with the help of a scalar electric potential, ϕ, as: E = −∇ϕ.
Finally, the energy of the device can be written from (2), with the help of P and ϕ, also with κ > 0, as follows, It arises from the Euler-Lagrange variation of the energy, δE = 0, that: and as also reported with details in [17]. The relations in (3) and (4) involve the Ginzburg-Landau formalism and Electrostatics equations. In particular, ∇ · D = 0 in Ω f and Ω p , where D is the electric displacement field, in Ω p .
These equations must be completed by conditions at the interface S f = Ω f ∩ Ω p and on the boundary S = ∂Ω. Namely, at the ferroelectric/paraelectric interface, the scalar field ϕ satisfies the continuity condition, and the electric displacement field is such that: where ν = (ν x , ν y , ν z ) T is the outward unit normal to Ω f , defined on S f . It follows that: Also, the scalar field P is subject to the condition: Finally, the electrodes affixed on the device require the potentials ∓φ 0 on the planes (z = ±γ), γ > 0, generating a voltage potential U = 2φ 0 (see Figure 1); that can be expressed through a boundary condition as follows, with ϕ S a given function, taking the values ∓φ 0 on (z = ±γ) in particular. Throughout the paper, the set of equations in (3)- (8), 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 ∈ (−∞, 0]. Also, we will take into account the extended consideration where κ ≥ 0 in the (EGL) system.
Let us mention furthermore that the volume of any considered ferroelectric layer Ω 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 Ω, 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.

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

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 ∈ N, we define the n-level truncation function, T n , as follows: x ∈ R −→ max{min{x, n}, −n} ∈ R, also expressed as The relations below derive then from this definition, Moreover, T n satisfies the property expressed in the following lemma.
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 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, with of course t and P 0 considered as previously.
Proposition 1. Let n ∈ N and (P n , ϕ n ) ∈ W . Furthermore, let κ > 0 and g be defined as aforementioned, as well as w, by assuming also that ξ and ε satisfy (9)- (12). Then, there exists a unique (P n , ϕ n ) ∈ W satisfying, for all (Q, ψ) ∈ W , Proof. If of course there exists (P n , ϕ n ) satisfying (18), then by adding the first relation of (18) to the second one multiplied by κ 16π 2 , we obtain: where, The bilinear map a(·, ·) and the linear map L(·) satisfy the Lax-Milgram Theorem hypotheses. In fact, we have, Furthermore, for (Q, ϕ) ∈ W , we have, thanks to (15), Thus, Finally, we obtain from the Lax-Milgram Theorem that there exists a unique (P n , ϕ n ) ∈ W satisfying (19), for all (Q, ψ) ∈ W . The first equation of (18) can be obtained from (19) by taking ψ = 0, and the second one by taking Q = 0 also from (19).
Proposition 2. Let n ∈ N. Furthermore, let κ > 0 and g be defined as aforementioned, as well as w, by assuming also that ξ and ε satisfy (9)- (12). Then, there exists (P n , ϕ n ) ∈ W such that, for all (Q, ψ) ∈ W , Proof. For each (P n , ϕ n ) ∈ W , we denote by (P n , ϕ n ), where (P n , ϕ n ) ∈ W , the solution of (18), which allows us to introduce f n : (P n , ϕ n ) ∈ W −→ (P n , ϕ n ) ∈ 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 , ϕ n ) ∈ f n (W ). There exists (P n , ϕ n ) ∈ W such that f n (P n , ϕ n ) = (P n , ϕ n ). We have, by choosing (P n , ϕ n ) as a test-function in (18), that: where C L is a continuity constant of the linear map L(·) and α is an ellipticity constant of the bilinear map a(·, ·), resulting from the proof of Proposition 1. Then, f n (W ) is bounded and, if B ⊂ W is the ball of center 0 W and radius C L /α, it follows that f n | B : B −→ B is correctly defined; • let us now note that f n (B) is relatively compact in W . We consider (P We denote by (P n , ϕ n ) its limit. Let (P n , ϕ n ) = f n (P n , ϕ n ). For k ∈ N, let us set We deduce, by using Lemma 1, that: ) k∈N converges to P n in L 2 (Ω f ) and the previous inequality leads to the convergence of (P • it remains to establish the continuity of f n . Let (P (m) n , ϕ (m) n ) m∈N be a sequence which converges to (P n , ϕ n ) in W . By repeating the above calculation and similar notations, we note that: Thus, (f n (P n )) m∈N converges to f n (P n , ϕ 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 , ϕ n ) ∈ W subject to (20).
We have now all ingredients to establish the main result here, namely the following statement.
Proof. The case where κ = 0 is immediate; it is even related to a determination of (P, ϕ) satisfying (13)-(14), where P can be expressed explicitly (as for instance in the frame of Remark 1 below). Let us set κ > 0 in the sequel. Let (P n , ϕ n ) n∈N ∈ 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 , ϕ n k ) k∈N converges weakly in W and almost everywhere in Ω f × Ω to an element (P, ϕ) of W . Next, we establish that ( Finally, we show that (P, ϕ) is a solution of (13)- (14). Before starting, let us note that, for all (Q, ψ) ∈ W , by adding the first equation of (20) to the second one multiplied by κ 16π 2 , it follows that: Then, • since g(x) → ∓∞, when x → ±∞, there exists β > 0 such that g(x) ≤ 0 for all x ≥ β and g(x) ≥ 0 for all x ≤ −β. Let n ∈ N. By choosing (P n , ϕ n ) as a test-function in (21), we have, thanks to the Green formula, that: In (22), the left hand-side terms are all positive. Furthermore, with (16), We deduce: where α is an ellipticity constant of the bilinear map a(·, ·) resulting from the proof of Proposition 1. Thus (P n , ϕ n ) n∈N is bounded in W and we can extract a subsequence (P nm , ϕ nm ) m∈N which converges weakly to (P, ϕ) in W . In particular, (P nm ) m∈N converges weakly to P in H 1 (Ω f ). Furthermore, as the injection of H 1 (Ω f ) in L 2 (Ω f ) is compact, the sequence (P nm ) m∈N converges to P in L 2 (Ω f ). Thus, there exists (cf., e.g., [5]) a subsequence (P nm p ) p∈N which converges almost everywhere to P in Ω f . In the same way, we can show that (ϕ nm p ) p∈N converges to ϕ in L 2 (Ω). It follows that there is a subsequence (ϕ nm pq ) q∈N which converges to ϕ almost everywhere in Ω. To simplify the notation, we will denote by (P n k , ϕ n k ) k∈N the built subsequence; • we show now that the sequence (T n k • g • P n k ) k∈N converges in L 1 (Ω f ) to g •P . Since (P n k ) k∈N converges to P a.e. in Ω f , we deduce, by continuity of g and thanks to (17), that (T n k • g • P n k ) k∈N converges to g(P ) a.e. in Ω f . If we establish that (T n k • g • P n k ) k∈N is uniformly integrable, we can conclude with the help of the Vitali convergence Theorem. Let us start by noting that ((T n k • g • P n k )P n k ) k∈N is bounded in L 1 (Ω f ). By reproducing the previous calculations, we obtain from (22) that, for k ∈ N, Thus, with (23), it derives that: Furthermore, with (16), it follows that: On the one hand, and, by choosing in particular M = 2C η , we obtain: On the other hand, we have, by using (16), that: So, for all E ⊂ Ω f satisfying |E| ≤ η 2 sup |x|≤M |g(x)| , it follows that: Thus, (T n k • g • P n k ) k∈N is uniformly integrable and the Vitali convergence Theorem allows to show that (T n k • g • P n k ) k∈N converges in L 1 (Ω f ) to g • P ; • in order to conclude, let us note that (P, ϕ) satisfies (13)- (14). Let (Q, ψ) ∈ C ∞ (Ω f ) × D(Ω). At first, by using the Green formula, we have that: Furthermore, as S f is the Lipschitz-continuous boundary of Ω f , there exists (cf., e.g., [4]) a constant K 2 , K 2 > 0, such that: L 2 (Ω f ) , and (P n k | S f ) k∈N converges to P | S f in L 2 (S f ). We deduce that: Ω f ∂ a (P n k − P )∂ b Qdx → 0, when k → +∞, and therefore, Similarly, we get that: Furthermore, with the Green formula, we have, With the similar arguments as before, it follows that: Moreover, with the Green formula, we get: and therefore, In the same way, we note that: Finally, with the Hölder inequality, we obtain that: Since it has been established that (T n k • g • P n k ) k∈N converges to g • P in L 1 (Ω f ), it follows that: 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 , ϕ n k ) k∈N , we get, for all (Q, ψ) ∈ C ∞ (Ω f ) × D(Ω), In particular, by choosing ψ = 0 and Q = 0 successively, in the above relation, we deduce that: Finally, the relations of the obtained system are valid for (Q, ψ) ∈ W , due to the density of C ∞ (Ω f ) in H 1 (Ω f ) and to the one of D(Ω) in H 1 0 (Ω). Then, (P, ϕ) is a solution of (13)- (14).
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 = 0, κ = 0, and ϕ S = 0 (as well as w = 0), then there are two distinct pairs (P, ϕ), (−P, −ϕ), with P := |t||P 0 |, satisfying the weak formulation (13)-(14).

Finite Element approximations and discrete nonlinear system
We discretize the formulation (13)- (14) by making use of a mesh of Ω 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 Ω f and Ω that are here polyhedral. In a correlative way, the discrete region associated with the boundary of Ω 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 Ω f and to another one of the discrete region associated with Ω p .
Let us introduce the discrete space W h , where, 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 , ϕ h ) ∈ W h such that, for all We are thus in particular concerned here with conforming discretizations of the regions Ω f , Ω, as well as with a conforming discretization of the space W ; namely, W h ⊂ 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 , ϕ h ) of this formulation leads to the determination of a discrete electric polarization associated with P, by using (1), wherein P is replaced by P h and E by −∇(ϕ h + w). Let us denote by {x j , 1 ≤ j ≤ N } the set of vertices of T h , numbered in such a way that x 1 , ..., x N f ∈ Ω f and x 1 , ..., x N 0 ∈ Ω, where N f , N 0 and N are respectively the number of vertices in Ω f , Ω and Ω; namely, N f < N 0 < N . Let v j be the shape function associated with x j , for j ∈ 1, N . We have thus, for (P h , ϕ h ) ∈ W h , where P j ∈ R for all j ∈ 1, N f and ϕ j ∈ R for all j ∈ 1, N 0 . The discrete formulation (24) can then be written as a nonlinear matrix system having two vector unknowns. It consists of finding P = (P 1 , ..., P N f ) T ∈ R N f and Φ = (ϕ 1 , ..., ϕ N 0 ) T ∈ R N 0 satisfying : where, 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 ε and ξ 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 ϕ 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.

Computational configurations
We consider three computational configurations. The meshes of Ω associated with these configurations are achieved in accordance with the previous section, and represented here by T On the one hand, geometrically, Ω f and Ω 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 Ω f . In a similar way, the radius and the height of the "cylinder" represented by Ω are denoted by R and H. On the other hand, the context where Ω f and Ω are parallelepipedal is considered, in particular for the last configuration. Namely, • for T h . 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.
Furthermore, we set: P 0 = 1, κ = 11.5, ε| Ωp = I 3 , for each computational configuration. We define ϕ s as below, in such a way that ϕ S (x, y, −γ) = φ 0 and ϕ S (x, y, γ) = −φ 0 , by recalling that the lower and upper faces of Ω are subparts of the planes (z = −γ) and (z = γ) respectively. Namely, γ = 0.85 in the frame of T (1) h , and 8.5 in the context of T (2) h as well as of T (3) h . Let us denote by U the voltage potential applied to the device, in association with the use in (13)- (14) and (24) of the function w such that w(x, y, z) = γ − z γ φ 0 − φ 0 , for (x, y, z) T ∈ Ω. Namely, w| S = ϕ 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 P h , in correspondence with (P h , ϕ h ); any state (P, ϕ), subject to (EGL), having a finite element approximation denoted here by (P h , ϕ 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 ∈ (−∞, 0].
Let us already mention that other values for P 0 and κ 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.

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 , ϕ h ) subject to (EGL) for which the polarization structures are similar. Figure 2 presents some of these profiles, resulting from T These profiles arise also from simulations based on T h , namely for instance with t = −20. In addition to these, the Obvious Profile, denoted by (OP) and corresponding to the case where P h ≡ 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) , ϕ (0) ) be a non obvious state at the reduced temperature t 0 , associated with the discrete representation (P (0) , Φ (0) ) T ∈ R N f +N 0 resulting from the nonlinear system (25). At the reduced temperature t = t 0 + δt (where δt > 0 for heating and δt < 0 for cooling), we initialize with this vector our algorithm for solving (25), and obtain a new solution, (P (1) , Φ (1) ) T , of this system. By proceeding in this way, step by step, we get a sequence ((P (k) , ϕ (k) )) k of states, namely a sequence ((P (k) , Φ (k) ) T ) k of discrete representations, where the solution (P (k) , Φ (k) ) T is provided at the reduced temperature t = t 0 + kδt.         (1) h the variations of the discrete "average polarization" marked by solid lines, when we "heat up" (with δt = 0.1) or "cool down" (with δt = −0.1) several states, for U = 0; namely for instance (3BP). While 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 T (2) h , as indicates Figure 3(b). For instance, the jumps schematized as below occur from T h . The most enduring profile with nonzero polarization is for T (1) h distinguished by (3BP), whereas for T (2) h , as well as for T (3) h , it is the one corresponding to (4BP). Each of the states associated with (4BP) has a constant "average polarization"; namely, P h ≡ 0.
At each temperature, we identified the most stable state by evaluating the energy E for each of the obtained states. For instance in the context of T (3) h , 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 ≤ t ≤ −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 Ω, 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 T  13.5, 13.5]. Let us consider, in accordance with the previous section also, three meshes corresponding to triangulations of Ω (1) , Ω (2) and Ω (3) respectively, and having similar sizes approximately equal to 0.69, with correlatively, 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 . For a state (P (i) , ϕ (i) ) subject to (EGL) and obtained by involving Ω (i) , with ϕ S = 0, we consider the following error, err , evaluated relatively with respect to the state (P (3) , ϕ (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 δ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

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 (P (0) , Φ (0) ) T is the discrete representation of the state (P (0) , ϕ (0) ) subject to (EGL) and associated with the voltage potential U = U (0) . For each fixed value k (k ∈ N), we denote by (P (k+1) , Φ (k+1) ) T the discrete representation that derives from (25) and corresponds to the voltage potential U = U (0) +(k +1)δU , δU ∈ R , after initializing our iterative algorithm with (P (k) , Φ (k) ) T , the discrete representation of the state (P (k) , ϕ (k) ) subject to (EGL). In this way, we deal with a sequence ((P (k) , Φ (k) ) T ) k of discrete representations of states (P (k) , ϕ (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 β, which reports in fact that β(P h ) ≈ 1 when the state is quasi-uniform. For instance, as indicates Figure 4, from simulations based on T h , with t = −20, we obtain |β(P h ) − 1| < 10 −1 for |U | ≥ 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 (δU > 0) the voltage potential and by decreasing it (δ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 T       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, T (1) h , is also denoted here by T (1) h3 , with h 3 ≡ 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 Ω, containing the same embedded polyhedral layer, also shaped like a cylinder. Namely, in addition to T (1) h3 , we have a coarser discretization of Ω, denoted by T (1) h2 , 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 Ω, denoted by T (1) h1 , 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 ≤ i ≤ 3 and u : Ω −→ R, let us set: u  , for i = 1, 2, evaluated relatively with respect to the state (P (3) , ϕ (3) ). Figure 6 shows the behavior of err (i) 2 , with respect to the voltage potential from U = −20, with δ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 T

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.