Limitations to Electrical Probing of Spontaneous Polarization in Ferroelectric-Dielectric Heterostructures

An accurate estimate of the ferroelectric polarization in ferroelectric-dielectric stacks is important from a materials science perspective, and it is also crucial for the development of ferroelectric based electron devices. This paper revisits the theory and application of the PUND technique in Metal-Ferroelectric-Dielectric-Metal (MFDM) structures by using analytical derivations and numerical simulations. In an MFDM structure the results of the PUND technique may largely differ from the polarization actually switched in the stack, which in turn is different from the remnant polarization of the underlying ferroelectric. The main hindrances that prevent PUND measurements from providing a good estimate of the polarization switching in MFDM stacks are thus discussed. The inspection of the involved physical quantities, not always accessible in experiments, provides a useful insight about the main sources of the errors in the PUND technique, and clarifies the delicate interplay between the depolarization field and the charge injection and trapping in MFDM stacks with a thin dielectric layer.


I. INTRODUCTION
Thanks to the discovery of a robust ferroelectricity in hafnium oxide thin films, [1] several intriguing applications of ferroelectricity in CMOS electron devices have been proposed and are being presently scrutinized. Device concepts include nanoscale CMOS FETs exploitating the effective negative capacitance to improve the subthresholdswing, [2]- [7] as well as Ferroelectric Tunnelling Junctions (FTJs), [8], [9] and ferroelectric FETs, [10], [11] which may be used as non-volatile memories or as memristors for neuromorphic computing applications. [12] In most of the above material systems and electron devices, one of the materials adjacent to the ferroelectric is a dielectric (see also Figure 1), and the working principle of the devices relies on the influence that the ferroelectric polarization, P, exerts on the band bending inside the device. Quite understandably, a dependable determination of P is of primary importance in ferroelectric materials and ferroelectric based electron devices. To this purpose the Positive-Up-Negative-Down (PUND) measurement technique was originally conceived for Metal-Ferroelectric-Metal (MFM) structures [13], [14] (see also PUND waveforms in Figure 2), and it is still routinely used also in Metal-Ferroelectric-Dielectric-Metal (MFDM) device structures. [15]- [18] The main goal of the PUND technique for an MFM structure is an accurate determination of the difference, 2P r , between the polarization at zero ferroelectric field for the arXiv:2201.12103v2 [cond-mat.mtrl-sci] 12 Apr 2022 positive and negative polarization state. In an MFM stack with ideal metal electrodes the field can be zeroed by using a zero external voltage, and the PUND measurements essentially intend to minimize the contributions to 2P r due to the background ferroelectric polarization and to possible leakage currents (see also the discussion in Section II). In an MFDM structure, however, a zero external voltage cannot force a zero ferroelectric field due to the depolarization field, and the application of the PUND technique to MFDM devices is not straightforward in several respects. In fact the interpretation of PUND results in an MFDM stack can lead to artifacts and to a misleading information about the polarization of the underlying ferroelectric layer. We here revisit the theory and application of PUND measurements in MFDM structures, and to this purpose we use both analytical derivations and a comprehensive modelling framework, that has been previously validated and calibrated against experiments. Our results show that: a) the determination of the spontaneous ferroelectric polarization is challenging in MFDM structures even if the charge injection and trapping in the dielectric stack is negligible; b) in the presence of a non negligible charge injection and trapping, the variations of spontaneous polarization and trapped charge are inextricably entangled, which further complicates the extraction of the switched polarization; c) the t D dependence of the extracted polarization can be an artifact due to the (t D dependent) interplay between the depolarization field and charge trapping. The paper is organized as follows. In Section II we present the theoretical background behind the extraction of the spontaneous polarization from the terminal currents measured by the PUND technique. In Section III we provide an overview of our in house modelling framework and then in Secs. III-B and III-C we report simulation results about the ability of the PUND technique to determine the switched polarization in MFDM structures, and offer an insight about the main sources of errors. In Section IV we propose a few concluding remarks and in the Appendix provide additional details about some aspects of the modelling framework.

II. CHARGE AND CURRENT AT THE ELECTRODES IN THE MFDM STRUCTURE
Let us consider the MFDM structure sketched in Figure 1(a), where r=(x, y) and z are the coordinate in the plane of the ferroelectric-dielectric interface and in the direction normal to the interface, while Q MF , Q MD are the charges per unit area respectively at the MF and MD electrodes, respectively (see Figure 1). Assuming a perfect screening in the metal electrodes, the electrostatic problem is linear and Q MF , Q MD can be written by using appropriate Green's functions for the charges in the structure. More specifically, for Q MF we have where A is the device area, P is the ferroelectric spontaneous polarization, E FT (r,t) denotes the z component of the electric field at the position r of the MF-FE interface (namely at z=−t F ). simulation results reported in Section III (so t U = 2t P ), if not otherwise stated. A preset pulse at V T = −5 V and for 125 µs is used to set an initial negative polarization state; b) Charge waveform sketch during the PUND simulation. Some key points are defined in order to simplify the notation of the paper.
At any time t, the E FT (r,t) is determined by the external bias V T and by the charges in the dielectric stack, that are charges located at coordinate r 0 defined as r at z = 0. In this latter respect, we here define the depolarization field, E DP (r,t), at the MF-FE interface as the field produced by the distribution of the total charge [P(r 0 ,t) + Q S (r 0 ,t)] at the FE-DE interface, where Q S (r 0 ,t) is an interface charge due to fixed Coulomb centers or to interface traps. We also similarly introduce E ρ (r,t) as the field produced by the remaining charge densities ρ(r 0 , z 0 ,t) in the dielectric stack at z 0 = 0. The E FT (r,t) can thus be written as where ε D being the thickness and relative permittivity of the dielectric, and t F , ε F being the thickness and background permittivity of the ferroelectric (see Figure 1(b)). By substituting Equation 2 in Equation 1 we obtain where P AV , E DP,AV and E ρ,AV denote respectively the average polarization and average fields at the MF-FE interface. From A of the Appendix it can be inferred that the average E DP and E ρ can be written as: where G MF (r 0 , z 0 ) is the Green's function defined as with E FT [r 0 , z 0 ] (r) being the field E FT (r) produced by a point charge e located at (r 0 , z 0 ). In B of Appendix we also demonstrate that, under realistic assumptions, the G MF (r 0 , z 0 ) for the MFDM structure in Figure 1(b) is independent of r 0 and it can be evaluated analytically. For a charge located at the FE-DE interface, for example, we have G MF (r 0 , 0) −(C F /C 0 ), which allows to rewrite Equation 4a as At any z 0 = 0 the G MF (r 0 , z 0 ) can be similarly expressed with a z 0 dependent capacitance ratio. We now recall that PUND measurements are based on the integral of the transient current at the electrodes, hence by definition the experiments can probe only the variations of the polarization and charges in the device stack. Here below the discussion is carried out in terms of the current I MF at the MF terminal; in B of Appendix we report the corresponding expression for I MD . In the presence of a trapping distributed throughout the device, it is thus difficult to express the influence of ρ(r 0 , z 0 ,t) on I MF , because Equation 4b shows that the information about the distribution along z 0 is required. Consequently, hereafter we simplify the picture and assume that the time derivative of the charge trapped in the dielectric stack is dominated by the (∂ Q S /∂t) term due to traps at the FE-DE interface, which implies ∂ E ρ,AV /∂t ∂ E DP,AV /∂t. We also assume that Q S (t) can change only through the terminal currents I QS,MF , I QS,MD shown in Figure 3, and we let I lkg denote a possible leakage current, not contributing to trapping. The current at the MF electrode can thus be written as where Q MF has been expressed via Equation 3 assuming ∂ E ρ,AV /∂t ∂ E DP,AV /∂t, and in the last equality we have used Equation 6.
If we now consider the Positive (P) pulse of a PUND experiment starting at t = 0 s with V T (0) = 0 V (see the waveform in Figure 2(a)), we can evaluate the charge Q P (t) (with 0 ≤ t ≤ t P ) by integrating the  Table I). The energy position of acceptor and donor type traps is also depicted. The tunnelling coefficients T MD and T MF depend on the energy, σ E , and geometric, σ T , cross sections of the traps, as well as on the traps energy E T .

expression for I MF in Equation 7 and obtain
where Q QS,MF and Q lkg are the integral of I QS,MF and I lkg , respectively. Moreover it is understood that P AV , E DP,AV and Q QS,MF in Equation 8 denote the variations from the corresponding values at t=0 or, equivalently, that Equation 8 conventionally assumes P AV =E DP,AV =Q S,MF =Q MF =0 at t=0. The charges Q U (t), Q N (t), Q D (t) during respectively the Up (U), Negative (N) and Down (D) pulses of the PUND technique have expressions equivalent to Equation 8. The first and last term at the right hand side of Equation 8 are the contributions due to, respectively, the linear polarization of the dielectrics and the leakage. Even in an MFM structure these contributions complicate the extraction of the remnant polarization 2P r =P AV (t = t P ), and PUND measurements address this issue by subtracting from Q P the charge Q U during the U pulse at the same external bias V T . In our theoretical framework such an approach results in the charge Q PU =(Q P −Q U ), that can be written as where the apices (P), (U) identify the P and U pulse and all charges are evaluated at times corresponding to the same V T value during either a rising or a falling V T ramp. We recall that, as already mentioned about Equation8, the P AV , E DP,AV and Q QS,MF in Equation 9 denote the variations from the corresponding values at the beginning of either the P or the U pulse. The second and third term in the right hand side of Equation 9 are due respectively to the depolarization field and the current at the MF electrode contributing to trapping at the FE-DE interface. In an MFM structure both these terms are negligible and, moreover, it is typically assumed that the polarization can be stabilized after the P pulse, so that P (U) AV is much smaller than P (P) AV . Furthermore, it is also usually assumed that the leakage affects the measurements to a similar extent during the P and U pulse, leading to Q lkg . [19] Under these circumstances Equation 9 shows that the Q PU in an MFM stack can be interpreted as the P (P) AV that we wish to determine. In an MFDM structure, instead, the terms in Equation 9 due to the depolarization field can be comparable to P AV appear much more delicate that in the MFM counterpart. This is systematically investigated in SectionIII by using numerical simulations.

III. NUMERICAL SIMULATIONS AND DISCUSSION
The goal of the PUND measurements is to accurately determine the spontaneous polarization switched by the external bias, which is very challenging in an MFDM structure due to the depolarization field and the possible charge trapping. In this section, we use numerical simulations to investigate the possible errors and artifacts produced by PUND measurements in MFDM structures, and to provide some useful physical insights.

A. Modelling framework and validation
Our in house developed simulation framework comprises models for the ferroelectric dynamics, a dynamic equation for the traps at the FE-DE interface and a description of a tunnelling injection from the MF and MD electrodes to the traps. The dynamics of the ferroelectric domains is described by a formulation of the multi-domain Landau-Ginzburg-Devonshire (LGD) model more thoroughly discussed in [20] where α, β , γ are the anisotropy constants, ρ denotes the resistivity that sets a time scale t ρ =ρ/(2|α|) for the ferroelectric switching. Moreover, the parameters 1/C i, j describe the depolarization energy and depolarization field in the MFDM structure, while k and w are the coupling constant and the domain wall width involved in the formulation of the domain wall energy. The mean values for α, β , γ used in simulations are reported in Table I and are referred to Hafnium-Zirconium-Oxide (HZO) ferroelectric which is the current state-of-the art ferrolectric used in FTJs. In all simulations the domain wall coupling k was set to zero, by following recent first principles calculations for HfO 2 . [21] Moreover, we assumed a resistivity ρ=115 Ωm, which is consistent with recently reported values for Hafnium-Zirconium (HZO) based capacitors, [22], [23] and results in a time scale for the ferroelectric dynamics t ρ ≈ 119.8 ns. All simulations include a n D = 1024 and a domain size d=5nm; we verified that results are insensitive to a further increase of n D . The modelling for the ferroelectric dynamics has been extensively compared to transient negative capacitance measurements, demonstrating a good agreement with experiments in asymmetric MFDM dielectric stacks, [20], [24] and also in symmetric MFDFM structures. [25] The charge trapping model follows a first-order dynamic equation for the occupation f T of either acceptor or donor corresponding to a ratio σ EC = 10% between the standard deviation and the mean value of the coercive field E C . m D and m F are the effective tunnelling masses for respectively Al 2 O 3 and HZO, while σ E , σ T denote respectively the energy and geometric cross section of the traps, that are used in the tunnelling model. The maximum energy values for the traps are 0.6 and 1.3 eV below the conduction band minimum at the FE-DE interface for acceptor and donor traps, respectively. Both traps type extend in energy for 2 eV below their maximum. The electron affinity was set to χ D = 1.4 eV for Al 2 O 3 [27] and to χ F = 2.4 eV for HZO, [28] while the workfunction for both TiN metal electrodes was taken as Φ M = 4.5 eV. [29] α 8 1.46 · 10 9 3.14 · 10 10 34 10 0.4 0.18 1 · 10 −19 7 · 10 −3 type traps at the FE-DE interface. By denoting with c MD0 , c MF0 the capture rate from the metal MD and MF electrodes the equation governing f T can be written as In the derivation of Equation 11 we used a detailed balance condition, ensuring that the steady state f T value at the equilibrium (i.e. for V T =0 V) is given by the Fermi function. In this work, the capture rates were attributed to tunnelling from and to the electrodes, and the tunnelling transmission described according to a WKB approximation that involves the tunnelling effective mass in the two dielectrics m D , m F , and the area and energy cross sections σ T [m 2 ], σ E [eV] (see also Figure  3). More details about the trapping and tunnelling models may be found in [9], [26] reporting also a good agreement with the polarization versus voltage curves in FTJs with an MFDM structure and a thin dielectric layers. The values of m D , m F , σ T , σ E used in this work are reported in Table I.
From the occupations f T one can readily calculate the charges Q acc and Q don respectively in acceptor and donor type traps, and finally the overall interface trapped charge Q S =(Q acc +Q don ). The knowledge of the time dependent polarization and trapped charge, in turn, allows one to numerically calculate all the quantities discussed in Section II, such as I MF , P AV , Q S,AV , and also I MF,QS , Q MF,QS .

B. Simulations results
We used numerical simulations to emulate PUND measurements with a 1 kHz waveform and a 5 V peak voltage in an MFDM structure with different dielectric thicknesses t D and trap densities N acc and N don . The ferroelectric HZO layer is 10 nm thick in all simulations.
We do not account for a possible leakage current flowing directly from MD to MF. Here we did not attempt to model the leakage current because we think that leakage is strongly technology dependent, frequently governed by Poole-Frenkel and hopping mechanisms in the the HZO and, as such, difficult to describe in simulations. [30] Moreover, while it is understood that a failure of the condition Q can induce artifacts in PUND experiments even in an MFM stack, this issue goes beyond the scope of the present paper, that is focused on the influence that the depolarization field and charge trapping have on the results of PUND measurements in MFDM structures.
In Figure 4 we report simulation results for t D = 1.5 nm and for different trap densities. The Q PUND is here defined as either Q PU =(Q P − Q U ) or Q ND =(Q N − Q D ), respectively for the positive and negative V T values. The sign of the interface charge Q S,AV is typically opposite to the sign of the polarization, and Figure 4 reports −Q S,AV , which together with P AV determines E DP,AV according to Equation 6. All charges in Figure 4 are referred to the corresponding value at the beginning of the P pulse, namely at t=0 s and V T =0 V in Figure 2 (see also the discussion about Equations 8 and 9). The Q PUND in Figure 4: Simulated charges corresponding to a 1 kHz PUND waveform applied to an MFDM structure. The Q PUND is either  Figure 4(a) shows a hysteresis loop that is much more tilted and stretched than in the corresponding MFM curves (filled triangles). The features for an MFDM are similar to those experimentally observed in the P-V curves for an HZO capacitor serially connected to a discrete ceramic capacitor ensuring a negligible charge injection, [31] or to measurements in MFDM structures with thicker Al 2 O 3 layers. [25] In fact, the relatively low density of traps in Figure 4(a) results in an interface charge Q S (green diamonds) that is practically negligible compared to the ferroelectric polarization (red squares). The lack of any compensation of the polarization results in a large depolarization field E DP,AV (see Equation 6), which in turn leads to a vast discrepancy between Q PUND and P AV . In fact, Figure 4(a) shows that for |V T | above about 4 V a complete polarization switching occurs. Nevertheless the corresponding Q PUND is much smaller than P AV , mainly because the E DP,AV term in Equation 9 subtracts from the P AV term due to the opposite sign. The results for Q PUND are quite different in Figure 4(b), because the Q S,AV can now compensate P AV to a large extent, thus drastically reducing the depolarization field. The hysteresis loop of the Q PUND curve in Figure 4(b) is qualitatively similar to the experimental behaviour observed in FTJ structures with a thin tunnel oxide, [9], [32] and the discrepancy between Q PUND and P AV is much smaller than in Figure 4(a). Figure 5 illustrates the same analysis as in Figure 4, but for a larger dielectric thickness t D = 2.5 nm. In this case the results of both small and large interface traps densities have a qualitative behaviour similar to Figure 4(a), namely the Q S is very small (green diamonds) and the compensation of the ferroelectric polarization is minimal. For both traps densities, the depolarization field results in a Q PUND much smaller than P AV . The different behaviour in Figure 5(b) compared to Figure 4(b) is due to the fact that, according to the tunneling effective masses and traps cross-sections reported in Table I, the trapping and de-trapping dynamics cannot follow the 1 kHz V T waveform for t D =2.5 nm or larger. In fact, because the HZO layer is 10 nm thick, in the simulations of this work the trapping dynamics is essentially set by the tunnelling through the much thinner dielectric layer. The lack of Q S modulation in Figure 5(b) is thus a dynamic effect. This emphasizes that the trapping induced compensation of the ferroelectric polarization requires both a large enough trap density at the FE-DE interface, and a trapping dynamics fast enough to respond to the V T waveform. This latter observation has been crucial in transient negative capacitance experiments, where thick dielectrics and fast bias waveforms were used to avoid the undesired compensation of the ferroelectric polarization and to achieve a hysteresis free behaviour. [20], [24], [33] While the t D values at which traps can no longer respond to a given V T waveform depend on the tunnelling model and the corresponding parameters in Table I, the qualitative trend is expected to be independent of the modelling details.

C. Discrepancies between Q PUND and P AV
The discrepancies between Q PUND and P AV shown in Figures 4 and 5 correspond to errors in the outcome of the PUND method. Hence in this section we evaluate the relative error E PU =|Q PU − P where from hereon all the quantities are evaluated at the end of the P or the U pulse, namely when V T is zero. Similar definitions apply to the N and D pulses, and the simulation results are also completely similar (not shown). Figure 6 shows the Q PU of the MFDM structures. It can be seen that a combination of a large t D and low concentrations of traps lead to low simulated Q PU values, because the corresponding ε 0 ε F E (P) DP,AV term is comparable to P (P) AV (as later shown in Figure 8). Figure 7 reports the evaluation of the error E PU for different dielectric thicknesses and trap densities. As it can be seen, the error tends to decrease for increasing trap densities, due to the corresponding reduction of the depolarization field E DP . For the same reason the error increases for thicker dielectrics. This latter behavior results in a t D dependence of the Figure 6: Q PU extracted from PUND simulations in MFDM structures for V T = 0 V, for different dielectric thickness t D and different traps density N acc =N don (in units of cm −2 eV −1 ). The remnant polarization 2P r extracted for an MFM structure is also reported for comparison. P AV estimated by the PUND method, which is an artifact of the method when it is applied to an MFDM structure. For t D = 2.5 nm the error is fairly insensitive to the trap density, because the Q S in the traps cannot respond to the V T waveform according to our tunnelling model. To gain an insight about the main causes of the errors shown in Figure 7, we first rewrite Equation 9 as where we have used Equation 6 to express the depolarization field E (P) DP,AV ; here we have omitted the leakage part because the leakage current is not included in our simulations, and all the quantities in Equation 12 are evaluated at the end of the P or the U pulse. Then we report in Figure 8(a) the quantities in the right hand side of Equations 6, 9 and 12, for a dielectric thickness t D = 1.5 nm and evaluated in the same condition used to evaluate the PUND error in Figure 7 (i.e. V T = 0 V). QS,MF (diamonds) related to the trapping and de-trapping current at the MF electrode are very small even for large trap densities, hence they do not appreciably influence Q PU in Equations 9 and 12. This is not surprising because, in the MFDM structures at study, traps exchange electrons primarily with the MD electrode as the dielectric is much thinner than the ferroelectric layer. Moreover, at large trap densities the Q (P) S,AV (filled circles) in the P pulse is comparable to P AV of PUND measurements in an MFDM structure for different thicknesses t D and different trap densities N acc = N don . The error is calculated for the P and U pulses. a) Error evaluated at the end of the P pulse (see inset); b) Error evaluated at the peak of the P pulse pulse.
always negligible compared to P (P) AV . This is because, for the case at study in Figure 8, the band bending in the dielectric at the end of the P pulse is such that the energy levels of both acceptor and donor traps fall below the Fermi level of the MD contact (see Figure 3(b)). Hence, essentially all traps have been filled at the end of the P pulse, and their occupation is not appreciably changed during the following U pulse. Figure 8(a) shows that also P (U) AV in the U pulse is much smaller than P (P) AV . This is because the P AV in the P and U pulse is a measure of the non reversible switching, whereas most of the switching in the U pulse is reversible in nature because it is the switching of those domains that have back switched after the P pulse. As mentioned above, Figure 8(a) shows that at low trap densities we have |Q AV , as it is confirmed by Figure 8(b). These are the conditions that in Figure 7 correspond to the maximum discrepancy between Q PU and P AV | term, which can be observed in Figure 8(b). The error in Figure 7 is correspondingly reduced at large trap densities; in fact Equation 9 suggests that Q PU tends to P

IV. CONCLUSIONS
We have revisited the theory and application of the PUND technique in MFDM structures by using analytical derivations and numerical simulations. The interplay between the depolarization field and charge trapping in an MFDM stack makes it difficult to obtain from the terminal currents alone an accurate estimate of the spontaneous polarization switched in the P or in the N pulse.
The discrepancies between Q PU and P AV that can be identified as an error of the PUND technique, it should be understood that neither the Q PU nor the P (P) AV of an MFDM structure are a good estimate of the 2P r of the underlying ferroelectric. This is because the depolarization field can be large at zero external bias, so that the MFDM structure at V T =0 is not at all representative of the ferroelectric material at zero ferroelectric field. More precisely the Q PU of the PUND technique tends to underestimate the nonreversible switched polarization P (P) AV , which in turn is an underestimate of the 2P r of the ferroelectric. The differences between these quantities depend on t D and on the density of traps, which may lead to artifacts in the characterization of a possible t D dependence of the properties of the underlying ferroelectric layer. Of course, we acknowledge that it would be very useful to suggest corrections to the PUND technique or to propose a novel technique for MFDM structures in order to duly account for the depolarizing field, and maybe even separate the switched polarization from the trapped charge. At the time of writing, however, we are not able to suggest a clear way of achieving such targets in an MFDM structure, and by relying exclusively on quantities accessible in experiments. In this respect, we cannot but conclude that more work is needed to improve the electrical probing of spontaneous polarization in ferroelectric-dielectric heterostructures.

ACKNOWLEDGEMENTS
This work was supported by European Union through the BeFerroSynaptic project (GA:871737).

APPENDIX A GREEN'S FUNCTION OF A POINT CHARGE IN THE MFDM STACK
In this section we discuss the analytical expression for the Green's function of the point charge defined in Equation 5. The potential ψ(r, z) produced by a point charge located in (r 0 , z 0 ) in a dielectric material having a relative dielectric constant ε r can be obtained by solving the Poisson equation [34] where e is the elementary charge. We now introduce the 2D Fourier transform of ψ(r, z) with respect to the coordinates r=(x,y), and define the Fourier pair Equation 14 assumes that the device area A is large enough that the integral over A is a good approximation of the indefinite integral over the entire (x,y) plane. 1 1 The formalism may be rephrased in terms of a Fourier series by assuming periodic boundary conditions for ψ(r, z) at the edges of the area A, that would however lead to identical results. [34]  If we now recall the identity and substitute Equations 14, 15 into Equation 13, we can readily infer that the unknown potential Ψ(q, z) takes the form [34] ψ(q, where φ (q, z) must satisfy the differential equation Let us now assume that the point charge is located at z 0 =0, namely at the FE-DE interface (see Figure 9). In this case the potential φ F (q, z) in the ferroelectric and φ D (q, z) in the dielectric region can be written as where the four q dependent constants C 1 , C 2 , C 3 and C 4 can be determined by using appropriate boundary conditions. At the interface with metal electrodes we used φ F (−t F )=0, φ D (t D )=0 , whereas at the FE-DE interface we employed the conditions φ F (0)=φ D (0) and [ε 0 ε F (∂ φ F (0)/∂ z) − ε 0 ε D (∂ φ D (0)/∂ z)]=e. By doing so we obtain where we have the notation more compact by introducing e F =exp (−2 qt F ) and e D =exp (2 qt D ). The Green's function G MF (r 0 , z 0 ) that we wish to determine is defined as where E FT (r) denotes the z component of the electric field at the MF-FE interface (i.e. at z=−t F ) produced by a point charge e located at (r 0 , z 0 ). By recalling the definition of the Fourier transform pairs in Equation 14 and then using Equation 16, we have Equations21, 22 finally provide For z 0 =−t F the limit in Equation23 can be readily calculated by using Equations 19a and 20a, so as to obtain G MF (r 0 , 0) = − C F C F +C D The corresponding Green's function G MD at the MD electrode, defined in Equation 32, can be derived with an entirely similar procedure. The result is so that (G MF (r 0 , 0) + G MD (r 0 , 0))=−1.
Similar derivations apply to the case of a point charge located in the ferroelectric (i.e. for −t F < z 0 < 0) or in the dielectric (i.e. for 0 < z 0 < t d ). In the former case we obtain G MD (r 0 , z 0 ) = − whereas in the latter case we have As it can be seen, even for an z 0 =0 we have [G MF (r 0 , z 0 ) + G MD (r 0 , z 0 )] = −1.

APPENDIX B CHARGE AND CURRENT AT THE MD ELECTRODE
The analysis of PUND measurements presented in the main paper is based on the current I MF at the MF electrode. According to the I MF , I MD definitions sketched in Figure 1, we see that I MD must be equal to I MF . For the completeness of definitions and derivations, we here report a concise analysis about the charge Q MD and current I MD at the MD electrode. We start with Q MD written as (see Figure 9) where E DB (r,t) is the z component of the field at the DE-MD interface at z=t D . The term ε 0 ε D E DB (r) can be expressed as ε 0 ε D E DB (r) = C S V T + ε 0 ε D E DI (r) where E DI (r) is the contribution to the field due to the total charge [P(r 0 ) + Q S (r 0 )] at the FE-DE interface. As already discussed in Section II of the main paper, we are here assuming that trapping is dominated by interface traps at the FE-DE interface.
We can now define the Green's function G MD (r 0 , z 0 ) at the MD electrode that allows us to write the average E DI,AV as