{'fraction_non_alphanumeric': 0.0371825715049452, 'fraction_numerical': 0.011921946003742315, 'mean_word_length': 4.526813414093662, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of Multi-Domains in Ferroelectric Tunnel Junction

The Discovery of giant tunnel electroresistance (TER) in Ferroelectric Tunnel Junction (FTJ) paves a futuristic possibility of utilizing the FTJ as a bistable resistive device with an enormously high ON/OFF ratio. In the last 20 years, numerous studies have reported that the formation of multidomain in ferroelectric material is an inevitable process to minimize the total system energy. Recent studies based on phase-field simulations have demonstrated that domain nucleation/motion substantially alters the electrostatics of a ferroelectric material. However, the impact of domain dynamics on quantum transport in FTJ remains elusive. This paper presents a comprehensive study of multidomain dynamics in a ferroelectric tunnel junction. Analysis of this article is twofold; firstly, we study the impact of domain dynamics on electrostatics in an FTJ. Subsequently, the obtained electrostatics is used to study the variations in tunneling current, and TER originated from multidomain dynamics. We show that ON/OFF current density and TER vary locally in the ferroelectric region. Furthermore, the device's electrostatics and quantum transport exhibit an oscillatory nature due to periodic domain texture. ON/OFF current density shows a sine/cosine distribution in ferroelectric, and approximately one-decade local variation in current density is observed. These local fluctuations in current density cause oscillations in the device's ON/OFF ratio. Optimization techniques to achieve a uniform and maximum TER are also discussed. A 2D analytical and explicit model is derived by solving coupled 2D Poisson's equation and Landau-Ginzburg equation. The model incorporates the switching and nucleation of domains by minimizing net ferroelectric energy (depolarization+free+gradient energy density). Furthermore, the impact of the bottom insulator layer on ferroelectric's gradient energy is also studied.


INTRODUCTION
A ferroelectric tunnel junction (FTJ) utilizes the mechanism of quantum tunneling to switch between bistable states 1−34 . To explore the microscopic physics of quantum tunneling in FTJ, first principle-based atomistic/numerical simulations were performed 1,6,9,16,20 . The ON/OFF characteristic of FTJ is defined by modulation in electrical conductance due to polarization switching. Subsequently, by utilizing the ON/OFF ratio of current density, tunnel electroresistance (TER) is calculated 1,2,4,9,13,22 . Moreover, many types of research have demonstrated that the ability of room temperature restive switching in FTJ exhibits a profound candidacy of the FTJ to utilize as a non-volatile memory component 4,5,7,8 . Furthermore, the researchers have studied epitaxial fabricated FTJs to achieve a higher value of the TER 11,12,14,15 . On the other hand, various researches demonstrated that the TER of FTJ can be tuned by altering the potential barrier or by engineering the electrical junctions 17,18,19,21,22 . Recently, the gigantic value of TER ≈ 10 7 is experimentally achieved by the potential barrier modulation due to Van der Waals heterojunctions in FTJ 23 . On top of that, the discovery of ferroelectricity in thin films of Hafnium Oxide (HfO 2 ) and zirconium doped Hafnium Oxide (H 0.5 Zr 0.5 O 2 /HZO) has triggered a vast amount of researches on CMOS compatible FTJ 24−27 . In the last 5-6 years, enormous studies on HZO based FTJ have demonstrated its promising future for non-volatile memories, incorporating excellent compatibility with the CMOS process 28−46 .
Even though plenty of studies is being reported on FTJs, the domain dynamics of FTJs are still not precise yet. Therefore, a comprehensive literature review is divided into four parts to understand the domain dynamics in ferroelectric materiel.

Analytical Models of domain dynamics in FE
The phenomenal paper by Bratkovsky et al. was the fundamental work that studied the impact of the insulator layer on the domain formation in a ferroelectric dielectric stack 47 . By the minimization of total ferroelectric energy, it is shown that enhancement in the insulator layer thickness leads to decrement in domain period (i.e., number of domain increases). A recent paper by Luk'yanchuk et al. has studied the motion of domains with applied voltages 48 . However, these developed mod-els assumed that thickness of the domain wall is zero. Hence, an abrupt transition from the upward domain to the downward domain is assumed . A finite domain wall width is considered by Morozovska et al., and a variational approach is used to derive the expression of polarization profile 49 . Later, Eliseev et al. proposed an analytical model which considers the finite domain wall width, and the impact of surface effect on domain formation is studied 50 . However, only vertical direction gradient is considered in the polarization.
Additionally, numerous researches also studied the domain dynamics in ferroelectric material by analyzing the domain wall velocity. The motion of domain wall by an applied electric field is classically defined by Merz law 51 . An exponential relationship between the wall velocity and the applied electric field is established. Later, a theory of domain wall motion was developed by Miller et al. 52 However, a discrepancy is observed between the developed approach and experiments 53,54 . Therefore, to study the domain wall velocity at the molecular level, density functional theory (DFT) is used 55−58 . The impact of the dead layer on the wall velocity by the polarization screening charges is reported by numerical and analytical methods 59 . Recently, the effect of reversed domain size over the domain wall velocity has been reported 60 . Furthermore, a charge and heat transport-based model of domain wall motion is proposed by Yudin et al. 61

Numerical/Phase-Field Models of domain dynamics in ferroelectric
Over the many years, numerical simulations or Phasefield models have played a vital role in analyzing the physics and motion of domains in the ferroelectric material. Domain texture for a wide range of temperature is captured by numerical solutions of coupled Poisson's equation and Landau-Ginzburg equation 62 . The study of super-lattice (ferroelectric/paraelectric) is reported by Stephanovich et al. 63 It is demonstrated that the polarizing domains in different layers interact with each other, subsequently, calculation of critical temperature was performed. The phase-field model for micro-structural domain growth and its macroscopic effect was proposed by Zhang et al., the nucleation and motion of domain walls are studied under the impact of both electric and mechanical boundary conditions 64 . Size effect on domain pattern growth and domain switching is studied by Wang et al. 65 Besides, the transition from multi-domain state to single domain state is also shown by phase-field simulations. Later, the 2-D polarization profile of multi-domains was proposed by Luk'yanchuk et al. 66 The textures of soft and hard domains w.r.t. ferroelectric thickness were also studied by the authors. The domain texture of a thin film ferroelectric is numerically evaluated by Kontsos et al. 67 Impact of strain and applied electric field on the domain evolution is shown in the paper.
In recent years phase field models have been used for the ferroelectric-dielectric layer to study the nucleation and domain wall motion in ferroelectric with applied voltage. The phase-field model was formulated by solving the time-dependent Ginzburg-Landau (TDGL) equation with the combination of Chensky-Tarasenko (C-T) formalism 68 . Strip domain pattern is observed when the ferroelectric thickness is scaled down to a critical value. Besides, to minimize net system energy, the formation of multi-domains is observed. The Phase-field modeling approach was also applied to study nucleation and domain dynamics in a ferroelectric FET 70 . The 2-D Poisson's equation coupled with Ginzburg-Landau equation is numerically solved to study the electrostatics of Metal-ferroelectric-Insulator-Metal (MFIM) and Metalferroelectric-Insulator-Semiconductor (MFIS) 70 . Additionally, the dynamics of soft and hard domain walls in ferroelectric heterostructures have been reported recently 71 . Note that the basis of these phase field models is the minimization of total ferroelectric energy by nucleation or motion of domain walls. Competition between domain wall energy, free energy, and depolarizing energy leads to two possible scenarios: a moment of the domain wall and the nucleation of a new domain. Whether a new domain will form or the domain wall will move is determined by the interplay between various energy components, which minimizes net ferroelectric energy.
Another phase field modeling class focuses on studying polycrystalline grains/domains in the ferroelectric material. Effect of grain boundary and its orientation on domain switching is studied by Choudhury et al. 72 A review article by Chen et al. shows that the formation of domains in a polycrystalline ferroelectric takes place to reduce the total free energy and domain wall energy of the ferroelectric material 73 . Liu, N et al. studies the impact of grain boundaries on the ferroelectric hysteresis loop. 74 Recently, a time-dependent domain switching model of a polycrystalline ferroelectric is presented by Dasgupta et al. 75 Model was based on Poisson Voronoi Tessellation (PVT) algorithm, which is used to nucleate the grains in the ferroelectric material. Later, the PVT algorithm is also used to study the variability in the P-E loop of a 3-D multi-granular ferroelectric capacitor 76 .

Ferroelectric Tunnel Junction Models
An analytical model based on Tsu-Esaki's current equation reveals that the TER of an FTJ can achieve an enormous value near zero bias by the polarization switching and direct tunneling through the FE barrier 77 . Theoretical calculations of polarization reversal and restive switching in an FTJ were performed by Kohlstedt et al. 78 Subsequently, numerical technique is used to evaluate the Tsu-Esaki current transport with depolarizing fields. A simple 1-D model of TER is presented by Zhuravlev et al. 79 screening effect of the metal electrode is included in the model by a simple expression of 1-D potential profile. It is shown that an increase in FE thickness leads to an enhancement in the TER ratio. An analytical model of tunneling transport through a FE is reported by Pantel  After a detailed literature review of domain physics in the ferroelectric, we notice that the dynamics of domains in an FTJ are still not coherent. Therefore, in the Problem Statement section, the significance/necessity of this work is described. And the methodology to capture the necessary physics is also explained. The Results and Discussion section is focused on the impact of multi-domains on the device's electrostatics and quantum transport. The device optimization section describes the guidelines to extract the highest performance from an FTJ. Furthermore, in the supplementary material, detailed mathematics of developed model and algorithm is given, which is used in this paper.   nucleation/motion of the domain wall due to the applied voltage is given. Therefore, the proposed study lacks to explain the oscillatory nature of electrostatics and TER/current in an FTJ.

PROBLEM STATEMENT & METHODOLOGY
Hence, up to the author's knowledge, no study has analyzed the impact of multi-domain on electrostatics and transport characteristics of an FTJ. Fig. 1 shows a schematic of two-terminal FTJ with a top insulator (dead) layer. The fundamental problem is to capture domain dynamics in the FE region. Firstly, we solve 2-D Poisson's equation with the gradient in polarization. x The 2-D solutions of the above equation are given in supplementary material. Green's function approach is generally used to obtain the electrostatic of a 2-D region bounded by Dirichlet and Neumann boundary conditions 113−120 . The next step is to incorporate the 2-D polarization wave with multi-domain. In section 1.2 of supplementary material, details about polarization modeling are presented. Firstly, we consider the equilibrium condition (zero bias condition), and the equilibrium domain period is calculated. Subsequently, the domain dynamics for non-zero applied bias is evaluated, a comprehensive discussion on methodology is given in section 1 of the supplementary material. Note that the devolved model is explicit and analytical in nature. Therefore, we also provide the first-ever compact model for FTJ with multi-domains. An increase in t f e allows the domain to diffuse further along the y− direction. Therefore, a domain's vertical directional gradient energy decreases ((∂P (x, y)/∂y) ↓) with an enhancement in t f e . Reduction in vertical directional gradient energy corresponds to a larger domain period because a lesser number of domains will be needed to minimize the net energy of the ferroelectric material. In other words, decre- ∂y reduces the gradient energy, which compensates the depolarizing energy, therefore, allows for a reduction in total system energy with a lesser number of domains. The potential profile of t f e = 5 nm and t f e = 10 nm are shown in Fig. 2 (c) and (d) respectively, reduction in the number of domains with increased t f e corresponds to lesser oscillations in potential distribution (frequency to oscillate from +V to -V). Fig. 3 shows the domain dynamics with the alterations in dielectric layer thickness (t ox1 ). In Fig. 3 (a)& (b) domain pattern becomes denser as t ox1 increases which is an opposite trend compared to t f e scaling. Reduction in domain period with increment in the t ox1 can be explained by observing the voltage drop across the FE layer. As t ox1 increases, voltage drop across the DE (V tox1 ) layer rises, leading to an enhancement in potential drop across the FE layer (V tox1 + V t f e = 0, K.V.L.). Depolarizing fields of the FE layer increases as the voltage across FE rises to compensate for increased depolarizing energy, more number of domains nucleates which minimizes the total system energy. Fig. 3 (c) and (d) show the potential profiles for t ox1 = 2 nm and t ox1 = 6 nm, respectively. An enhancement in t ox1 raises the oscillations in the potential profile due to the increased number of domains. Fig. 4 shows the polarization and potential distribution in FE and DE regions with applied bias (V 2 ). Alteration in the magnitude of applied bias triggers the moment of the domain wall. Polarization texture with positive applied bias is shown in Fig. 4 (a)-(c). A positive enhancement in V 2 leads to an expansion in the upward domain period (or reduction in the downward domain period). At V 2 = 0.5 V, the growth in the upward domain period is substantial, and a reduction in the downward domain period is observed simultaneously. Therefore, potential distributions' oscillations (frequency to oscillate from + to -value) decline as V 2 increases (see Fig. 4

Non-equilibrium Electrostatics of domains
On the other hand, when applied bias is shifted along the negative direction, an expansion in the downward domain period is observed (see Fig. 4 (a), (g) & (h)). However, oscillations in the potential profile again decrease as the V 2 swapped along the negative path. Therefore, all the domains will switch upward/downward at a relatively high bias, and the FE layer will attain a homogeneous potential (no oscillations) distribution.

TER and Quantum Transports with Multi-Domain
From Fig. 1 to Fig. 4, it is crystalline that the formation and moment of domains in the FE region remarkably  affect the device's electrostatics. Since quantum transport of electronic charge carries is originated from electrostatics. Therefore, the domain dynamics will significantly alter the device's carrier transport characteristics. This section presents a thorough analysis of quantum transport in FTJ with multi-domain. The T(E) varies enormously with different horizontal positions. Such large T(E) variations can be explained by observing Fig. 2(c). At the location, x a + T D /2 potential is positively maximum. On the other hand, at the x a = 3 + T D /4 position, the potential is minimum (negative), besides the difference between +V to −V is in the order of a few hundred mV. These huge oscillations in potential profiles lead to a spatial dependent T(E) that is extremely sensitive to the x-positions. Due to the periodicity of the domain wave, the T(E) at x a + T D is approximately the same T(E) at x a (see Fig. 2(c), potential along the x-direction is also periodic).
Due to multi-domains, the spatial variations in T(E) will introduce local fluctuations in current density. The absolute magnitude of ON and OFF current densities are shown in Fig. 6 (a)-(d). Fig. 6 (a) shows the 3-D plot of ON current density w.r.t applied bias (V 2 ) and lateral direction (x), respectively. It is fascinating to observe that ON current density oscillates along the lateral direction. The periodicity is originated due to the oscillations in polarization charges of domains (see Fig. 2 Fig. 2(c)).
The domain wall starts to move with the enhancement in applied bias at that moment the switching mechanism starts to kick in (see Fig. 4(a)-(f)). Switching of many domains decreases the oscillations in both polarization charges as well as in potential profile (see Fig. 4(c) & (f)). Fig. 6(b) and Fig. 6(d) show the surface plot of ON current density, as voltage rises, the domain period of the upward domain increases, which corresponds to the same amount of reduction in the downward domain period. Hence, the oscillations in current density decline by switching of domains (moment of domain wall). Depending on the device's physical parameters and value of spontaneous polarization, the ratio of J max /J 0 can approach up to 10, which shows significant variability in the current density via dynamics of multi-domains. The OFF current density (J OF F ) is shown in Fig. 6(c)-(d), a π/2 phase difference is observed between ON and OFF current densities. Therefore, if J ON follows the sine profile, then J OF F will follow cosine distribution and vice-versa.
The next step is to analyze the spatial variations in the ON/OFF ratio of the device. TER is defined as a modulation in conductivity due to the polarization reversal. Here, TER is defined as the ON/OFF current density ratio. Fig. 6 (e) shows 3-D variations in the TER for various combinations of V 2 and x-positions. As anticipated, TER also exhibits periodicity along the spatial x-direction. However, TER sharply shoots up to a peak value at some specific points. These peak points can be evaluated physically by observing minima positions of J OF F at which J ON /J OF F attains its maximum value. Hence, it can be stated that if J ON and J OF F follow sine and cosine profiles, subsequently, the spatial variations in TER can be approximate as a tangent (tan(x)) function. A 2-D surface plot of TER is shown in Fig. 6 (f). The TER is maximum only for a narrow range of x and V 2 combinations. But, oscillations in TER can be diminished by the interplay in domain period with FE and DE layer thickness, and a broader range of peak TER can be achieved. Hence, in the next section, we present device optimization techniques which are useful for achieving a larger TER with a broader peak.
Device optimization for maximum TER Spatial variations in TER are originated due to the periodicity in upward and downward domains (see Fig.  6(e) & (f)). And peak value of TER is obtained only for a narrow range of x-locations which is unwanted for an adequate device. A broader range of peak TER can be achieved by engineering the domain texture of FE. Fig. 7 shows the optimization of TER by tuning the domain density in the FE region. Enhancement in FE thickness reduces the gradient energy (along y-direction), which relaxes the necessity of minimization of depolarization energy. Therefore, system energy can now be minimized by fewer domains. Hence, domain density decreases as FE thickness increases. At more significant t f e , the rate of oscillations in electrostatic device decreases, therefore, the spatial variations in TER start to attenuate, and a wider range of peak TER is obtained. At t f e = 5 nm, the maximum to minimum TER ratio is ≈ 15. On the other hand, t f e = 25 nm, and this ratio approaches unity due to the significantly reduced oscil-lations. The value of TER rises with an enhancement in the t f e . This is due to the increased potential drop across the FE layer at a larger value of t f e . A larger potential drop enhances the asymmetry between tunneling barriers of up and down polarization states, leading to a higher TER. Additionally, domain texture starts to move from multi-domain state to mono-domain state with increasing FE layer thickness. Fig. 7(d)-(f) show the TER plots for various insulator thicknesses. The asymmetry in tunneling barriers for ON and OFF state polarization rises as t ox1 increases, leading to an enhancement in the TER value. Nevertheless, the depolarization energy also rises with increasing t ox1 . Hence, a denser domain pattern nucleates to minimize the depolarizing energy. Therefore, spatial variations in TER begin to shoot up as insulator thickness rises and the range of the x-positions at which TER is maximum begins to shrink. Therefore, a trade-off exists between the broader peak of TER and   Fig. 8 shows the polarization profile and corresponding potential surface plot of an FTJ with a bottom insulator layer (t ox2 ). Polarization is now maximum at the middle of the FE region (y = t ox1 + t f e /2) and minimum at FE/DE interfaces. Hence, FTJ with front and bottom insulator layers will have a larger polarization gradient (along the vertical direction) than FTJ with a single insulator layer. The comparison between gradient energy density and polarization profile of with and without bottom insulator layer is shown in Fig. 9. Note that the gradient energy density of ferroelectric along y-direction (E Dy ) is higher when the bottom dead layer is incorporated (see Fig. 9(b) and (d)). Enhancement in the E Dy can be attributed to the fact that now polarization profile is minimum-maximum-minimum at y = t ox1 , y = t ox1 + t f e /2, and y = t ox1 + t f e + t ox2 , respectively, which raises the gradient in the y-direction. Hence, in Fig. 9(d) E Dy is higher than Fig. 9(b).
More domains are nucleated to compensate for increased energy density due to the higher E Dy decreasing the depolarization energy to obtain the minima of total system energy. But, increment in domain density decreases the domain wall width, which leads to enhancement in a gradient along the x-direction. Therefore, the enhancement in x-directional gradient energy will dominate at a certain level, and further nucleation of domains will not help minimize the total system energy. Thus, the interplay between gradient (both x and y directional), depolarization, and free energy density is crucial in determining the optimum texture of domains in the FE region. Incorporating the bottom dead layer increases the oscillations in the polarization profile (number of domain rises). Therefore, the frequency of oscillations in potential will also go up. Hence, the variations in TER and current density will follow a similar trend as observed in electrostatics.

CONCLUSION
A comprehensive study of domain dynamics in a ferroelectric tunnel junction is presented. Spatial periodicity in polarization charges induces the oscillations in device electrostatics. Since transport characteristics of devices such as tunnel electroresistance (TER), transmission probability (T(E)), and tunneling current originate from electrostatic potential, therefore, electronic transport also exhibits spatial variations. The ON and OFF current density varies ≈ 12 times in the transient region from the upward domain to downward domain. Due to these significant variations in current density along x-direction (horizontal-direction), the device's ON/OFF ratio (TER) also exhibits an oscillatory nature along the x-direction. The TER of the device attains its peak value only for a narrow range of x-positions and is periodically distributed throughout the ferroelectric region. To make TER maximum and uniformly along the x-direction, an optimization in domain texture is required, which can be achieved by altering the DE/FE thicknesses. Furthermore, it is demonstrated that incorporating the bottom DE layer enhances the gradient energy (y-directional). Hence, domain texture becomes denser to compensate for increased energy density. However, the denser domain pattern will enhance the oscillatory nature of the device's electrostatics and its transport characteristics.

Solution of 2-D Poisson's equation
Eq (1)-(3) are a system 2-D partial differential equations (PDE) that possess Dirichlet and Neumann boundary conditions at the interfaces. Here, we use Green's function approach to obtain accurate and analytical solutions of potential functions for ferroelectric (FE) and dielectric (DE) regions. The first step is to derive Green's functions for various FE/DE regions. Subsequently, these functions are plugged in Green's theorem to obtain the expression of potential distribution 1  Our recent article derived an analytical model for metal-ferroelectric-insulator-semiconductor (MFIS) negative capacitance FET. The model provides coupled solutions of 2-D Poisson's equation with the Landau-Ginzburg equation 9 . However, a mono-domain state in the ferroelectric region is assumed. Hence, the model is not capable of capturing the effect of multi-domains. Therefore, in this work, the ∆.P term in (2) is considered to incorporate the multi-domain dynamics in the ferroelectric (FE) region. Note that both horizontal and vertical directional gradient in the polarization is incorporated. The model is derived for a 1/14 multi-domain ferroelectric tunnel junction with top and bottom insulator (dead) layers. The 2-D potential distribution functions for the top insulator (t ox 1 ) , ferroelectric (t f e ), and bottom insulator (t ox 2 ) are as follows.

Top Insulator layer
Ferroelectric region Bottom Insulator layer All the necessary series constants and coefficients are given in the Appendix.
The potential equation in the FE region (5) comprises six different Fourier series terms. The methodology to derive each term and significance of the same is as follows.
Physics and dynamics of multi-domain are primarily incorporated by the φ MD (x, y) and φ n=0 (x) terms. The gradient in the polarization charges is responsible for the φ MD (x, y) expression, which is evaluated by plugging Green's function of the ferroelectric region into the following equation.
Above identity leads to the following potential expression.
Eq. 8 contains two indices: n and j. The polarization distribution of multi-domains is expressed in the Fourier series (see section 1.2), which is responsible for summation over j in (8). Another series over n is arising due to Green's function of the FE region.
At the limit of n → 0, sinh(kx) can be approximate as kx. Hence, Green's function for the n = 0 term is evaluated (given in Appendix) and the corresponding potential function is derived.
The term φ L,R (x, y) is originated due to the left (x = 0) and right (x = L) boundaries and calculated as.
cos k II n t ox 1 + t f e − y B n 1 sinh k II n (L − x) + B n 2 sinh k II n x sinh(k II n L) The interfaces along top-insulator/FE (x, t ox 1 ) and FE/Bottom-insulator (x, t ox 1 + t f e ) contribute to the φ T,B (x, y) term in Eq. 5.
D 1 and D 2 are the Fourier series coefficients and are calculated by the potential continuity conditions.
The expression of D 1 and D 2 is obtained as follows.
Interface Fourier series coefficients (D 1 and D 2 ) are defined in a fashion such that electric displacement field continuity at FE/DE interface is satisfied. Thus, the developed model satisfies both 2-D potential continuity and 2-D electric displacement field continuity at the FE/DE interfaces.
The remaining two potential terms in (5): φ m=0 (x) and φ G x (x, y) are evaluated from m = 0 and G x (x, y; x , y ), respectively.
All the necessary Fourier series constants are given in Appendix. Figure 2. Schematic of polarization wave with different x 1 and x 2 upward and downward domain periods, respectively.

Polarization wave modeling
Polarization in the FE region is assumed to follow a sinusoidal wave with a period of T D and expressed in the form of the Fourier series.
where, a 0 is an average of Fourier wave and A j , B j are Fourier series coefficients of polarization wave. The horizontal direction polarization follows a sinusoidal nature with the upward domain period x 1 and downward domain period x 2 . Fig. 2 shows the schematic of polarization wave plotted at the mid of ferroelectric region (y = t ox 1 + t f e /2). The vertical directional (y) variations in the polarization are incorporated by the function f (y), which is defined by a parabolic function.
f (y) = ay 2 + by + c Parabolic constants a, b and c are given in Appendix. Fig. 3 describes the nature of the function f (y) . The black curve shows an FTJ with front and back dead layers (t ox 2 = 0). The function takes its maximum value in the middle of the ferroelectric region (y = t ox 1 + t f e /2), and it is 0 at the insulator/ferroelectric interfaces. On the other hand, the red curve plotted for t ox 2 → 0, the f (y) is maximum (1) at the ferroelectric/metal interface, and it is minimum (0) at the insulator/ferroelectric interface. The green curve shows a constant f (y) with no gradient along the y-direction.
Impact of f (y) over the multi-domain polarization profile is shown in Fig. 4. In Fig. 4(a), polarization is maximum at the mid of ferroelectric region ( f (y) = 1) and it is zero on the ferroelectric/insulator interfaces ( f (y) = 0). Fig. 4(b) shows that polarization attains its maximum value at the metal/ferroelectric interface, which is also observed in the numerical phase field simulations 10,11 . In the absence of both insulator layers, it is assumed that there is no gradient in the polarization along the vertical direction (see Fig. 4(c)).

Ferroelectric energy dynamics
The following equation defines the net energy density of a ferroelectric material 10−12 .
f net = f f ree + f dep + f grad (20) where, f f ree , f dep , and f grad represent the free energy density, depolarization energy density, and gradient energy density, respectively.

4/14
The free energy density of a ferroelectric material can be expressed in the power series of polarization with associated Taylor series coefficients 13,14 .
f f ree = f 0 + αP 2 (x, y) + β P 4 (x, y) + γP 6 (x, y) (21) where, α, β , and γ are the material-dependent parameters evaluated or taken from the experimental data used to validate/calibrate the developed model 15,16 . The strain is not exclusive as an order parameter in (21). Rather, it is included in the Landau coefficients as reported by Saha et al. 10 Depolarization energy of a ferroelectric material plays a crucial role in the formation of multi-domains. It is a well -established fact that to minimize depolarization energy domain formation occurs 10−14 . The depolarization energy density ( f dep ) is obtained from the 2-D local distributions of electric fields and polarization profile in the FE region 10,11 .
The first term in (22) shows the energy density due to spatial stray fields, and the second term signifies the contribution due to depolarizing fields in the FE region. The 2-D local electric fields are evaluated by differentiating potential function.
Domain polarization exhibits gradients in both x and y directions 10−12 . These gradients in polarization texture cause one more energy density term defined as gradient energy density ( f grad ). where, g 11 and g 22 are the constants and tuned to calibrate the model with experimental results 15,16 , note that both x and y directional gradients have been considered.

Calculation of equilibrium domain period
Both upward and downward domains exhibit an equal period at the zero applied voltage (V 1 = V 2 = 0), 11,17,18 . Hence, net domain period is calculated as. . 5 shows the polarization distribution in an FTJ for the equilibrium condition (V 1 = V 2 = 0). Fig. 5(a) shows the 3-D polarization profile in the ferroelectric region, and Fig. 5(b) shows a surface plot for the same. Fig. 5(c) shows a schematic of polarization wave with equal upward and downward domain periods (x 0 ). The equilibrium domain period (x 0 ) is to be determined by minimizing net ferroelectric energy, which corresponds to a thermodynamically stable state 13,17,18 . where, F net is the net ferroelectric energy and f net is calculated from (20). For the given, physical and material parameters such as t ox 1 , t ox 2 , t f e , L, ε x , ε y , ε ox 1 , ε ox 2 α, β , and γ, the only remaining unknown in (27) is x 0 , calculated by the following equation.  Evaluated x 0 from (28) is a function of physical parameters such as insulator and ferroelectric thicknesses. Discussion related to the impact of various device parameters on the equilibrium domain period is given in the manuscript's main text (see Fig. 2 and Fig. 3 of the main manuscript). Fig. 6 shows the polarization profile for a non-zero applied bias, (V 2 = 0). Due to the applied bias the equilibrium condition is no longer valid. Hence, domain periods of upward and downward domains are not identical (x 1 = x 2 ). For a non-zero value of V 2 domain wall starts to move, and depending on the polarity of voltage, an expansion or reduction in a domain period is observed. Therefore, a ∆x change in the x 1 leads to an equal amount of variation in the x 2 . It is also crystalline from Fig. 6 that expansion in a domain period reduces the adjacent domain period. Hence, for V 2 = 0, and assuming a stationary condition (no transient) the domain period is given as.

Algorithm to calculate the domain period for non zero applied bias
Note that, total domain period evaluated by (29) is the same as calculated by (28). However, due to a non-zero applied bias, an asymmetry in the upward and downward domains is observed. (x 1 = x 2 ). The minimization of net ferroelectric energy calculates the shift in domain wall (∆x).
The above equation is solved for each applied bias, and calculated ∆x is plugged in (5) to derive the 2-D ferroelectric potential distribution incorporating multi-domains effect. The flow chart of the developed algorithm is shown in Fig. 7.

Transition of Domain wall from Soft to Hard
The phase-field simulations demonstrate that as ferroelectric thickness increases, the domain wall gradually shifts from soft wall to hard wall. These secondary effects can be incorporated into the model by molding the polarization wave. A generic expression of polarization profile is derived as follows.
P(x, y) = f (y) a 0 2 + 2P s d w T D ∑ j A j cos (k j x) + B j sin (k j x) k 2 j (31) Fig. 8(a) shows the polarization wave with a hard domain wall. A sudden change in polarization occurs, leading to an extremely high gradient along the lateral direction. Fig. 8(b) shows the polarization wave texture for various ∆ values. Domain wall gradually transits from soft domain wall to hard domain wall as the width of domain wall (∆) increases. Therefore, ∆ can be used as a parameter to incorporate the domain wall transition. However, corresponding Fourier coefficients of polarization wave should be used in the model to capture the domain dynamics.

Quantum Transport Model with Multi-Domains
The transmission probability of the charge carriers is evaluated by the Wentzel-Kramers-Brillouin (WKB) approximation. T (E) = exp(−2γ) (14) where, m * is the electron mass in DE/FE layer 19 , and φ (x s , y) is the DE+FE layer potential profile at a particular x s position. The 3-dimensional current density is calculated by a generic expression of quantum transport which is valid for all ranges of temperature 20 .
Eq. (33) comprises both tunneling and thermionic transport of the charge carrier. The J T →B denotes current from top to bottom (t ox 1 → t f e ), U T is the reference energy point = 0, and µ T is the electron Fermi level (µ T ≈ U T + 0.5q) 20 . The difference between metal work function and electron affinity of the DE layer is defined by ∆φ ox 1 B . Note that the first part of (33) is due to tunneling transport, and the second part is due to thermionic emission. Therefore, all possible current transports have been considered in the model. In similar fashion J B→T can be evaluated by replacing U T , µ T , and ∆φ ox 1 B with U B , µ B , and ∆φ ox 2 B /∆φ FE respectively, and net current density is calculated as.

Metal Screening Effect
The metal electrode does not entirely screen polarization charges. Hence, some residual part of these polarization charges remains unscreened in the ferroelectric layer. To incorporate such effect, the Thomas-Fermi screening approach is used 21,22 .
where, l i is the screening length, Q s is screening charge density, and ε M is the metal permittivity 22 . Note that, screening charge density just changes the effective applied voltage by the ∆φ i amount, hence, effective applied voltage is given as.
where, i = 1, 2 for top and bottom metal electrodes.