Electrodynamics of ferroelectric films with negative capacitance

TL;DR

This paper develops a comprehensive electrodynamic theory for ferroelectric ultra-thin films with periodic domain textures, revealing negative capacitance phenomena, resonance modes in the THz range, and implications for THz optics and metamaterials.

Contribution

It introduces a novel theoretical framework linking domain textures to negative capacitance and THz resonances in ferroelectric films, enabling new optical and metamaterial applications.

Findings

Identification of THz resonance modes of domain oscillations

Discovery of singular phase points related to negative capacitance

Potential for THz negative capacitance-based optical devices

Abstract

We construct a comprehensive theory of the electrodynamic response of ferroelectric ultra-thin films containing periodic domain textures (PDT) with 180{\deg} polarization-oriented domains. The focal point of the theory is the negative-capacitance phenomenon which naturally arises from the depolarization field induced by PDT. We derive frequency-dependent dielectric permittivity related to the PDT dynamics across the entire frequency range. We find the resonance mode of domain oscillations in the THz spectral band and the singular points in the phase of the reflected THz beam that are intimately related to the negative capacitance. Our findings provide a platform for the THz negative capacitance-based optics of ferroelectric films and for engineering the epsilon-near-zero plasmonic THz metamaterials.

Figures

Equations

$$\mathbf{D}={\epsilon }_0{\epsilon }_{\Vert }\mathbf{E}+{\mathbf{P}}_{\mathrm{dw}}.$$

$${\epsilon }_{\mathrm{f}}=\frac{D}{{\epsilon }_{0}E},$$

$${\epsilon }_{\mathrm{f}}=\frac{{\epsilon }_{\parallel }D}{D-P_{\mathrm{dw}}}.$$

$${\epsilon }_{\mathrm{f}}={\epsilon }_{\parallel }-\frac{\pi \varsigma }{4\ln 2}{\left(\frac{{\epsilon }_{\perp }}{{\epsilon }_{\parallel }}\right)}^{1/2}\frac{2a_{\mathrm{f}}}{d}{\epsilon }_{\parallel },$$

$$d\simeq \sqrt{3.53{\left({\epsilon }_{\perp }/{\epsilon }_{\Vert }\right)}^{1/2}\varsigma \delta \cdot 2a_{\mathrm{f}}}.$$

$$\mu \stackrel{..}{x}(t)+\eta \stackrel{.}{x}(t)+kx(t)=2P_{\mathrm{s}}{E}_{\mathrm{D}}(t),$$

$$k=\frac{4P_{\mathrm{s}}^2}{{\epsilon }_0{\epsilon }_{\Vert }d}g,g\simeq 1-\frac{4\ln 2}{\pi \varsigma }{\left(\frac{{\epsilon }_{\parallel }}{{\epsilon }_{\perp }}\right)}^{1/2}\frac{d}{2a_{\mathrm{f}}}.$$

$$\gamma (\omega )=\frac{g^{-1}{\omega }_0^2}{{\omega }_0^2-{\omega }^2-i\Gamma \omega },{\omega }_0={\left(\frac{4P_{\mathrm{s}}^{2}g}{\mu d{\epsilon }_0{\epsilon }_{\Vert }}\right)}^{1/2},$$

$${\epsilon }_{\mathrm{f}}(\omega )=\frac{{\epsilon }_{\Vert }}{1-\gamma (\omega )}={\epsilon }_{\Vert }\frac{{\omega }_0^2-{\omega }^2-i\Gamma \omega }{\left(1-g^{-1}\right){\omega }_0^2-{\omega }^2-i\Gamma \omega },$$

$${\epsilon }_{\mathrm{tot}}(\omega )=\frac{{\omega }_{\mathrm{LO}}^2-{\omega }^2-i\Gamma \omega }{{\omega }_{\mathrm{TO}}^2-{\omega }^2-i\Gamma \omega }{\epsilon }_{\mathrm{tot}}(\infty ),$$

$$r_{123}^p=\frac{r_{12}^p+r_{23}^p{e}^{2i\beta }}{1+r_{12}^p{r}_{23}^p{e}^{2i\beta }},$$

$$r_{ij}^p=\frac{{\epsilon }_j{\xi }_i-{\epsilon }_i{\xi }_j}{{\epsilon }_j{\xi }_i+{\epsilon }_i{\xi }_j},i,j=1,2,3$$

$${\xi }_1=\sqrt{{\epsilon }_1}\cos {\theta }_1,{\xi }_2={\left({\epsilon }_{\perp }-\frac{{\epsilon }_{\perp }}{{\epsilon }_2}{\epsilon }_1{\sin }^2{\theta }_1\right)}^{1/2}$$

$${\xi }_3={\left({\epsilon }_3-{\epsilon }_1{\sin }^2{\theta }_1\right)}^{1/2}.$$

$${\epsilon }_{\mathrm{f}}=\frac{{\epsilon }_{\parallel }}{1-\left(\eta (D)/D\right)P_{\mathrm{s}}}.$$

$$({\epsilon }_{\parallel }{\partial }_z^2+{\epsilon }_{\perp }{\partial }_x^2){\phi }^{(f)}=0,$$

$${\epsilon }_{\mathrm{p}}({\partial }_z^2+{\partial }_x^2){\phi }^{(p)}=0.$$

$${\epsilon }_0{\epsilon }_{\parallel }{\partial }_z{\phi }^{(f)}-{\epsilon }_0{\epsilon }_{\mathrm{p}}{\partial }_z{\phi }^{(p)}=\pm {\sigma }_{\pm }(x),$$

$${\phi }^{(f)}={\phi }^{(p)}.$$

$${\sigma }_{\pm }(x)=\pm \vartheta (x)P_{\mathrm{s}}=\pm \left(\eta +\sum _{n=1}^{\infty }{p}_n\cos q_{n}x\right)P_{\mathrm{s}}$$

$$p_n=\frac{4}{\pi n}\sin \frac{1+\eta }{2}\pi n$$

$${\phi }^{(f)}$$

$${\phi }^{(p)}$$

$$\pm P_{\mathrm{s}}\sum _{n=1}^{\infty }{\psi }_n\exp q_n(a_{\mathrm{f}}-\left|z\right|)\cos q_{n}x$$

$${\psi }_n=\frac{p_n}{q_n{\epsilon }_0\sqrt{{\epsilon }_{\perp }{\epsilon }_{\parallel }}\coth \left(\sqrt{{\epsilon }_{\perp }/{\epsilon }_{\parallel }}{q}_n{a}_{\mathrm{f}}\right)+q_n{\epsilon }_0{\epsilon }_{\mathrm{p}}}.$$

$$F_{\mathrm{el}}$$

$$F_{\mathrm{el}}=\frac{2a_{\mathrm{f}}}{{\epsilon }_0{\epsilon }_{\parallel }}\left(\frac{1}{2}{P}_{\mathrm{s}}^2{\eta }^2-D{P}_{\mathrm{s}}\eta \right)+\frac{8}{{\pi }^3{\epsilon }_0\varsigma \sqrt{{\epsilon }_{\perp }{\epsilon }_{\parallel }}}df(\eta )P_{\mathrm{s}}^2,$$

$$f(\eta )=\sum _{n=1}^{\infty }{n}^{-3}{\sin }^2(1+\eta )\frac{\pi n}{2}\stackrel{\eta \to 0}{\simeq }\frac{7}{8}\zeta \left(3\right)-\frac{\ln 2}{4}{\left(\pi \eta \right)}^2.$$

$$\eta =\frac{D}{P_{\mathrm{s}}}{\left(1-\frac{4\ln 2}{\pi \varsigma }{\left(\frac{{\epsilon }_{\parallel }}{{\epsilon }_{\perp }}\right)}^{1/2}\frac{d}{2a_{\mathrm{f}}}\right)}^{-1}.$$

Full text

Introduction

Functional properties of ferroelectric films and superlattices are drastically different from those of the bulk materials. Extensive experimental Streiffer et al. (2002); Fong et al. (2004); Zubko et al. (2010, 2012); Hruszkewycz et al. (2013); Yadav et al. (2016) and theoretical Bratkovsky and Levanyuk (2000, 2001); Kornev et al. (2004); Stephanovich et al. (2005); De Guerville et al. (2005); Aguado-Puente and Junquera (2008); Luk’yanchuk et al. (2009) studies revealed new physics that emerges from the periodic domain textures (PDT) stabilized by the strains and the depolarization fields. One of the most striking features of the PDT is the manifest of the negative capacitance Zubko et al. (2016), the phenomenon that has been a subject of an intense recent attention Salahuddin and Datta (2008); Cano and Jiménez (2010); Zhirnov and Cavin (2008); Appleby et al. (2014); Khan et al. (2015); Catalan et al. (2015). The interest is motivated by both, its fundamental importance, and by its high potential for technological applications, in particular, as a platform for novel low-dissipation field-effect transistors.

Another facet of emergent functionalities of ferroelectric films arises due to the strain-tunability of their multiscale spontaneous polarization dynamics that provides a wide range of operational frequencies from a few kilohertz to tens terahertz ($1\,\mathrm{THz}$ $=10^{12}\,\mathrm{Hz}$). The spectral band $\ 10^{4}\,\mathrm{Hz}-0.3\,\mathrm{THz}$ is covered by the relaxation dynamics of the domain walls (DWs) and polar clusters. The far-infrared spectral region of $3-30\,\mathrm{THz}$ is governed by the soft-mode vibrations of the polar ions. Recent studies suggested that oscillations of domain walls (DW) in PDT occur in the frequency window of $0.3-3\,\mathrm{THz}$ Zhang et al. (2011); Luk’yanchuk et al. (2014); Hlinka et al. (2017), i.e within the least studied frequency range, referred to as a THz gap.

Our work steps into the breach and demonstrates that the resonant behavior of oscillations of PDT is the collective effect similar to the standard plasmonic excitations in metals. We find that this effect is governed by the stiffness of depolarization field and is a consequence of the negative capacitance. We calculate the electrodynamic response of PDT and find the frequency dependence of the effective dielectric permittivity. We investigate the THz optics of the PDT structure and demonstrate the existence of the topologically protected phase-singular points of the absolute darkness in the Fresnel reflection coefficient. These points are a fingerprint of the negative capacitance phenomenon.

Formation of the PDT lowers the energy of the depolarization field induced by surface depolarization charges that appear at polarization termination points. Had the polarization maintained the same direction throughout a slab, see Fig. 1a, the electrostatic energy stored by the depolarization field would be proportional to the volume of the whole system hence huge. Splitting the system into the PDT implies that the surface depolarization charges form stripes of alternating signs. As a result, the depolarization field becomes restricted to the near-surface shell, see Fig. 1b, hence drastically diminishing the related electrostatic energy. Although predicted already in earlier works by Landau Landau and Lifshitz (1935, 2013) and Kittel Kittel (1946) in the context of ferromagnetic systems, the PDT in ferroelectric films has long been viewed as unlikely until the recent direct experimental evidences for equilibrium $180^{\circ}$ stripe domains in strained ferroelectric thin films of PbTiO3 (PTO) deposited on the SrTiO3 (STO) substrate Streiffer et al. (2002); Hruszkewycz et al. (2013) (as sketched in Fig. 1b) and in PTO/STO superlattices Zubko et al. (2010, 2012) (Fig. 1c). The observed PDT behaviors appeared to well follow the theoretical predictions Bratkovsky and Levanyuk (2000, 2001); Kornev et al. (2004); Stephanovich et al. (2005); De Guerville et al. (2005); Aguado-Puente and Junquera (2008); Luk’yanchuk et al. (2009). In particular, it was found that the mono-domain $z$-oriented state in the unscreened films is always unstable with respect the PDT formation.

Imagine, however, that surfaces of the slab get short-circuited, for example, by the metallic electrodes. Then charges would acquire the capability of arbitrary sliding along the surfaces and flowing between the electrodes, thus compensating the depolarization field inside the ferroelectric. As a result, the depolarization energy is nullified, the domain walls would creep away, and the mono-domain state stabilizes. This tendency of a ferrolecric with the PDT to self-generating the charge at its surfaces can be expressed via attributing to the ferroelectric slab a negative capacitance. One should bear in mind, however, that in reality, the slab that is not linked to the electrodes, maintains the PDT since charges cannot propagate across the dielectric and the charge redistribution described above would not occur.

It is the overscreeing effect due to the depolarization field that lies in the origin of the negative capacitance. The component of the depolarization field penetrating the bulk of the ferroelectric slab from the surface is opposite to the applied field and exceeds it by the absolute value. As a result, the total field inside the slab is counter-directed to the induced polarization. Hence the dielectric permittivity is negative. Writing down the coupled PDT equations of motion and Maxwell equations, we derive the dynamic frequency-dependent dielectric permittivity of the ferroelectric films and superlattices containing PDT across the entire frequency range and find the corresponding resonance mode. We show that its emergence is intimately related to the low frequency negative dielectric permittivity.

Results

Negative capacitance

We start the description of the domain electrodynamics with the static limit, and show that PDT in response to the applied field develops the fascinating negative dielectric permittivity $\varepsilon_{\mathrm{f}}$ and, accordingly, the negative capacitance, $C=\varepsilon_{0}\varepsilon_{\mathrm{f}}{S}/(2a_{\mathrm{f}})$, of the ferroelectric layer of surface area $S$ and thickness $2a_{\mathrm{f}}$. We further demonstrate that the negative static $\varepsilon_{\mathrm{f}}$ of the PDT is a key element for striking arising of the resonance mode of DWs oscillations at THz frequencies.

The applied field is characterized by the induction ${\mathbf{D}}$ that conserves all across the system, see Fig. 1d, and generates the net polarization of the ferroelectric layer, ${\mathbf{P}}\parallel\mathbf{z}$ (the electrostatic fields in the ferroelectric slab with domains are macroscopically averaged). In our ferroelectric film the induction is related to the electric field through

[TABLE]

The contribution $\mathbf{P_{\mathrm{dw}}}$ stems from the motion of the domain walls. In the absence of the motion, Eq. (1) reduces to standard electrostatic relation between the electric field and induction. The dielectric constant $\varepsilon_{\|}$ is the intrinsic permittivity of the ferroelectric material along the polarization direction. The ${\mathbf{P}}_{\mathrm{dw}}$ term reflects the change in the overall polarization due to the displacement of the DWs altering thus the relative contribution from the “up” and “down” -oriented domains. This implies the imbalance of the related depolarization charges of the opposite sign leading to the nonzero average charge density ${\sigma}=P_{\mathrm{dw}}$. Thus the depolarization field ${\mathbf{E}}_{\mathrm{dep}}=-(\varepsilon_{0}\varepsilon_{\|})^{-1}\mathbf{P}_{\mathrm{dw}}$ penetrates the bulk of the slab and is directed opposite to the applied field $\mathbf{D}$. This field, together with the driving field $\mathbf{E}_{\mathrm{D}}=(\varepsilon_{0}\varepsilon_{\|})^{-1}\mathbf{D}$ which would have settled down if DWs were immobile, contributes to the total field, $\mathbf{E}=\mathbf{E}_{\mathrm{D}}+\mathbf{E}_{\mathrm{dep}}$. Hence the sign of the system’s dielectric permittivity, which is defined as

[TABLE]

depends on the relative magnitude of the oppositely oriented $\mathbf{E}_{\mathrm{D}}$ and $\mathbf{E}_{\mathrm{dep}}$.

Making use of the relation (1), one obtains

[TABLE]

To derive $\varepsilon_{f}$, we have to know how the polarization related to the DW displacement depends on the applied field $D$. The calculations presented in the Appendix result in

[TABLE]

where the first term stems from the positive intrinsic contribution of $\mathbf{E}_{\mathrm{D}}$, whereas the second, the negative one, is generated by the depolarizing field $\mathbf{E}_{\mathrm{dep}}$, reflecting the effect of the moving DW. The depolarizing term outweights the intrinsic term, hence the negative dielectric permittivity settles, provided the film thickness, $2a_{\mathrm{f}}$ exceeds the domain width, $d$. The equilibrium value of $d$ can be either taken from the experiment Streiffer et al. (2002); Zubko et al. (2012) or estimated by the Landau-Kittel square root law Landau and Lifshitz (1935); Kittel (1946); Landau and Lifshitz (2013) specifically adapted for ferroelectrics in Bratkovsky and Levanyuk (2000); De Guerville et al. (2005); Catalan et al. (2012),

[TABLE]

The resulting typical domain configuration for PTO films is shown in Fig. 1b. The DW thickness, $\delta$, is about 1 nm Meyer and Vanderbilt (2002). In the above formulas $\varsigma=1+\varepsilon_{\mathrm{p}}/\left(\varepsilon_{\|}\varepsilon_{\perp}\right)^{1/2}$, and the values of intrinsic permittivities along and across the polarization direction, $\varepsilon_{\|}$, $\varepsilon_{\perp}$, and of paraelectric layer, $\varepsilon_{\mathrm{p}}$, are specified below; in the experimental range $1\leqslant\varsigma\leqslant 4$. If the sandwiching paraelectric layers possess different permittivities, $\varepsilon_{\mathrm{p}}^{+}$ and $\varepsilon_{\mathrm{p}}^{-1}$, the effective parameter $\varsigma_{\mathrm{eff}}$ is defined by the relation $\varsigma_{\mathrm{eff}}^{-1}=\frac{1}{2}\left[\varsigma^{-1}(\varepsilon_{\mathrm{p}}^{+})+\varsigma^{-1}(\varepsilon_{\mathrm{p}}^{-})\right]$. In particular, if the ferroelectric film is deposited on the metallic substrate so that $\varepsilon_{\mathrm{p}}^{-}=1$ and $\varepsilon_{\mathrm{p}}^{+}=\infty$, one has $\varsigma_{\mathrm{eff}}\simeq 2$. The polarization profile in this case can be obtained by the image method.

DW oscillations

To quantify the electromagnetic behavior of PDT, we derive the dynamic response function $\gamma(\omega)$, defined by the relation $\mathbf{P}_{\mathrm{dw}\,\omega}=\gamma(\omega)\,\mathbf{D}_{\omega}$, for the periodic structure of the domain stripes aligned along the $\mathbf{y}$-axis which forms in the ferroelectric layer. The polarization axis $\mathbf{z}$ is perpendicular to the film plane, while the in-plane DW motion occurs along the $\mathbf{x}$-axis, the index $\omega$ indicates the Fourrier transforms of the respective quantities. We neglect the longitudinal DW fluctuations which can broad the resonance peak Brierley and Littlewood (2014).

In the harmonic oscillator approximation Sidorkin (2006) the driven dynamics of DWs is given by

[TABLE]

where $x$ is the coordinate of alternating DW displacements (Fig. 1d), the coefficients $k$, $\mu$ and $\eta$ are calculated per unit of DW area, and $2P_{\mathrm{s}}E_{\mathrm{D}}(t)$ is the pressure due to the induction-induced driving field forcing DW to move to flip the surrounding spontaneous polarization from $-P_{\mathrm{s}}$ to $+P_{\mathrm{s}}$.

Equation (6) presumes that the DW is rigid and atomically narrow. However, in the tens-nanometer-thin films the DW broadens significantly due to depolarization effects in the large range of temperatures (yet the effective thickness of the DWs remains smaller then the distance between them as long and the model holds). The broadening is especially pronounced on approach to the film surface Luk’yanchuk et al. (2009) and upon the DW displacement the polarization reverses first at the surface and then propagates into the interior of the film Zhang et al. (2011). As a result of the DW broadening, the underlying pinning potential for the DW due to the periodic atomic structure gets essentially reduced and becomes relevant only at low temperatures.

The stiffness $k$ is of a purely electrostatic origin and, therefore, is not too sensitive to the explicit shape of DW. We calculate it as a coefficient at the restoring force that pushes the domain walls back in order to reduce the depolarization field, $E_{\mathrm{dep}}=-\Delta P_{\mathrm{s}}/\varepsilon_{0}=-2(x/\varepsilon_{0}d)P_{\mathrm{s}}$, arising upon the displacment of DW from the equilibrium. This field is induced by the extra depolarization surface charges, $\sigma=\pm\Delta P_{\mathrm{s}}$, appearing due to the spontaneous polarization excess, $\Delta P_{\mathrm{s}}=2(x/d)P_{\mathrm{s}}$, and is directed antiparallel to $\Delta P_{\mathrm{s}}$. We evaluate the depolarization energy of the system as the electrostatic energy of the dielectric slab with permittivity $\varepsilon_{\|}$ and with surface charges $\sigma=2({x}/{d})P_{\mathrm{s}}$ as $W=\left(2a_{\mathrm{f}}\sigma^{2}/2\varepsilon_{0}\varepsilon_{\|}\right)S$, where $S$ is the surface area. The corresponding energy per unit area of the displaced DW is $w=(d/2Sa_{\mathrm{f}})W$. Relating it with the harmonic oscillator stiffness energy $kx^{2}/2$, we find

[TABLE]

The factor $g\leq 1$, calculated in the Appendix, makes the above qualitative consideration more precise, taking into account the non-uniform part of the depolarization field near the surface caused by the stepwise distribution of the depolarization surface charges. As follows from (7) in very thick samples, these corrections vanish and $g\rightarrow 1$. In realistic cases that we considered here, namely, the PTO films with $2a_{\mathrm{f}}\simeq 10-30$ nm, it varies from 0.4 to 0.9.

The effective DW mass, $\mu$, and the viscosity, $\eta$, in Eq. (6) are related to the motion of the material-constituent polar ions during the displacement of the DWs Kittel (1951) and with the piezoelectric effect of the time-varying depolarization field Sidorkin (2006). The magnitude of $\mu$ was calculated in ab-initio simulations of DW dynamics in the sub-THz range Zhang et al. (2011). We propose an interpolation formula $\mu\,\mathrm{[kg/m}^{2}\mathrm{]}\simeq 1.3\sqrt{2a_{\mathrm{f}}\,\mathrm{[nm]}}\times 10^{-9}$, which fits well the numerical result. Viscosity, $\eta$, is expressed via damping factor of DW motion, $\Gamma=\eta/\mu$. Since the consistent theory for $\Gamma$ is not available, we chose the one that appears in the measurements of the soft-mode relaxation of polar ions.

THZ resonance

In what follows we calculate the complex dynamic permittivity of PDT, $\varepsilon_{\mathrm{f}}(\omega)$. The resulting plots of the real, $\Re\/\varepsilon_{\mathrm{f}}(\omega)$, and imaginary, {\Im}$$\varepsilon_{\mathrm{f}}(\omega), parts of the dispersion $\varepsilon_{\mathrm{f}}(\omega)$ for 30 nm thick film of PTO are shown in Fig. 2a.

The reaction of the ferroelectric layer to the time-dependent applied field $\mathbf{D}(t)=\int d\omega\,\mathbf{D}_{\omega}e^{-i\omega t}$ is characterized by the linear response function, $\gamma(\omega)$. Solving Eq. (6) by Fourier transformation, we find:

[TABLE]

where the characteristic oscillation frequency is the usual $\omega_{0}=\left(k/\mu\right)^{1/2}$ with $k$ and $\mu$ related to the system parameters as discussed above. When deriving Eq. (8), we took into account that $\mathbf{P}_{\mathrm{dw}\,\omega}=2P_{\mathrm{s}}\,\left(x_{\omega}/d\right)\mathbf{z}$ and $\mathbf{D}_{\omega}=\varepsilon_{0}\varepsilon_{\|}\mathbf{E}_{\mathrm{D}\,\omega}$.

Shown in Fig. 2b are the real and imaginary parts of $\gamma(\omega)$ for our sample exhibiting a characteristic Lorentzian-type form. The latter is a legacy of a singularity at $\omega_{0}$, which would have been exhibited by a system without the dissipation, illustrated by the dashed lines. The emergence of the resonance response to the applied field is clearly demonstrated by the frequency dependencies of the absolute value of the response function, $|\gamma(\omega)|$, for different film thicknesses, see Fig. 3a, and reaching its maximum at the resonance frequency $\omega_{\mathrm{r}}^{2}=\omega_{0}^{2}-\Gamma^{2}/2$. The dependence of $\omega_{\mathrm{r}}$ upon the film thickness, $2a_{\mathrm{f}}$, is shown in Fig. 3b. Shown by the blue line in the same panel is the Landau-Kittel dependence of the domain width, $d$, on the film thickness given by Eq. (5) .

Equation (8) enables optimization of the materials parameters and the film thickness to make sure that $\omega_{\mathrm{r}}$ falls within the desired THz frequency range. In particular, for the strained films of PTO with the high spontaneous polarization, $P_{s}\simeq 0.65\,\mathrm{C\,m}^{-2}$, see Pertsev et al. (2003), relatively low permittivities $\varepsilon_{\|}\simeq 100$, $\varepsilon_{\perp}\simeq 30$, soft mode damping factor $\Gamma\simeq 20\,\mathrm{cm}^{-1}$ ($0.6\,\mathrm{THz}$) Hlinka et al. (2011) and $\mu$ defined as above, the resonance frequency $\omega_{r}$ decreases and spans the range from $1.5$ to $0.75\,\mathrm{THz}$ when $2a_{\mathrm{f}}$ increases from $10$ to $40\,\mathrm{nm}$. We find the damping frequency, $\omega_{\mathrm{d}}^{2}=\omega_{0}^{2}-\Gamma^{2}/4$, of the attenuated oscillations of domains in PTO, $x(t)=x_{0}e^{-(\Gamma/2)\,t}\sin\omega_{\mathrm{d}}t$, which is slightly larger then $\omega_{\mathrm{r}}$, see Fig. 3b. Remarkably, our formulas perfectly describe the results of ab-initio simulations of DWs oscillations in PZT free-standing ultrathin films Zhang et al. (2011). The calculated damping frequency, $\omega_{\mathrm{d}}$, is shown by the dashed line in the Fig. 3b, the symbols display the results of simulations. Here we used the following parameters for PZT: $P_{\mathrm{s}}\simeq 0.40\,\mathrm{C\,m}^{-2}$ Pertsev et al. (2003), $\varepsilon_{\|},\varepsilon_{\perp}\simeq 350$, $\Gamma\simeq 27\mathrm{cm}^{-1}$ ($0.8\,\mathrm{THz}$) Buixaderas et al. (2010) and the same value for DW mass, $\mu$. At the same time, the resonance frequency of PZT, $\omega_{\mathrm{r}}$, drops rapidly with the thickness and vanishes above 4 nm. The data in this paragraph are given for ceramic samples. The values for strained films can be slightly different.

Making use of the relation $\gamma(\omega)=1-\varepsilon_{\|}/\varepsilon_{\mathrm{f}}(\omega)$ that follows directly from Eq. (1), we obtain the dynamic permittivity of a ferroelectric layer

[TABLE]

with $1-g^{-1}<0$. At small $\omega$ this yields the negative static permittivity described by Eq. (4). We plot the real, $\Re\/\varepsilon_{\mathrm{f}}(\omega)$, and imaginary, {\Im}$$\varepsilon_{\mathrm{f}}(\omega), parts of complex dispersion, $\varepsilon_{\mathrm{f}}(\omega)$, calculated for 30 nm thick film of PTO in paraelectric environement, see Fig. 2a. At low frequencies $\Re\varepsilon_{\mathrm{f}}(\omega)$ remains negative. At high frequencies the oscilations of DWs freeze out, and therefore $\varepsilon_{\mathrm{f}}(\omega)$ should become positive and equal to intrinsic permettivity, $\Re\/\varepsilon_{\mathrm{f}}(\infty)=\varepsilon_{\|}>0$. At frequency $\omega=\omega_{0}$, the real part of permittivity turns zero, $\Re\/\varepsilon_{\mathrm{f}}(\omega_{0})=0$. This behavior of PDT permittivity resembles that of a metal, having the negative real part of permittivity below the plasma frequency at which it becomes equal to zero. Accordingly, this will lead to the plasma resonance at $\omega=\omega_{0}$ in a PDT as a response to the driving field $\mathbf{D}$.

The discovered resonance mode gives rise to the Drude-Lorentz frequency response of a biased capacitor consisting of the ferroelectric layer (of thickness $2a_{\mathrm{f}}$) sandwiched between paraelectric buffer layers (of thickness $a_{\mathrm{p}}$ and permittivity $\varepsilon_{\mathrm{p}}$ each). The same response appears in $2a_{\mathrm{f}}/2a_{\mathrm{p}}$ ferroelectric-paraelectric superlattice. The effective permittivity, calculated for the system of in-series capacitors, $\varepsilon_{\mathrm{tot}}^{-1}(\omega)=\alpha_{\mathrm{p}}\varepsilon_{\mathrm{p}}^{-1}+\alpha_{\mathrm{f}}\varepsilon_{\mathrm{f}}^{-1}(\omega)$, results from the interplay of the positive, paralectric, and negative, ferroelectric, contributions and remains positive at small $\omega$. Here $\alpha_{\mathrm{p}}=a_{\mathrm{p}}/(a_{\mathrm{p}}+a_{\mathrm{f}})$ and $\alpha_{\mathrm{f}}=a_{\mathrm{f}}/(a_{\mathrm{p}}+a_{\mathrm{f}})$ are the relative weights of the respective layers. That it is the paraelectric contribution that dominates the behavior of the effective permittivity at $\omega=0$, is inherently related to the very emergence of the domains. Namely, it is related to the fact that the thickness of the paraelectric buffer, $a_{\mathrm{p}}$, is larger then the domain width, $d$, and thus domain depolarization stray fields do not interact with electrodes Mokrỳ et al. (2004). The frequency dependence of the $\varepsilon_{\mathrm{tot}}(\omega)$ for PTO/STO $30/50$ $\mathrm{nm}$ system is shown in Fig. 2c. One sees (the detail of calculations are given further) that $\Re\,\varepsilon_{\mathrm{tot}}(\omega)$ maintains a positive value for low frequencies and reflects the presence of the oscillatory mode at $\omega$ slightly less then $\omega_{0}$. As before, the dashed line depicts the behavior of $\Re\,\varepsilon_{\mathrm{tot}}(\omega)$ in the absence of dissipation. We thus introduce a new metamaterial operating in the THz frequency-range $\varepsilon_{\mathrm{tot}}(\omega)$ tunable by the variation of ratio $2a_{\mathrm{f}}/2a_{\mathrm{p}}$.

Optics

The effective permittivity of a ferroelectric/paraelectric layered system is calculated as $\varepsilon_{\mathrm{tot}}^{-1}(\omega)=\alpha_{\mathrm{p}}\varepsilon_{\mathrm{p}}^{-1}+\alpha_{\mathrm{f}}\varepsilon_{\mathrm{f}}^{-1}(\omega)$ and can be conveniently written in a standard form:

[TABLE]

where the characteristic frequencies $\omega_{\mathrm{\scriptscriptstyle LO}}$, $\omega_{\mathrm{\scriptscriptstyle TO}}$ are the analogues of the longitudinal and traversal oscillation frequencies in the polar mode spectroscopy. The former one, is the DW oscillation frequency, $\omega_{\mathrm{\scriptscriptstyle LO}}=\omega_{0}$. The latter frequency can be obtained from the relation $\omega_{\mathrm{\mathrm{\scriptscriptstyle TO}}}^{2}/\omega_{\mathrm{\mathrm{\scriptscriptstyle LO}}}^{2}=\varepsilon_{\mathrm{tot}}(\infty)/\varepsilon_{\mathrm{tot}}(0)$ that is the analog of the Lyddane-Sachs-Teller relation. Here $\varepsilon_{\mathrm{tot}}^{-1}(0)=\alpha_{\mathrm{p}}\varepsilon_{\mathrm{p}}^{-1}+\alpha_{\mathrm{f}}\varepsilon_{\mathrm{f}}^{-1}$ and $\varepsilon_{\mathrm{tot}}^{-1}(\infty)=\alpha_{\mathrm{p}}\varepsilon_{\mathrm{p}}^{-1}+\alpha_{\mathrm{f}}\varepsilon_{\|}^{-1}$. Equation (10) describes the dispersion curve in Fig. 2c. It clearly shows that the obtained PDT resonance mode freezes out when the frequency exceeds $\omega_{0}$ and what is left are the molecular vibrations including the polar soft mode oscillations.

The derived dispersion relation for the dielectric permittivity provides a foundation for the description of the THz optical properties of a ferroelectric film with PDT, in particular, for finding its complex refraction index, $n(\omega)=\sqrt{\varepsilon(\omega)}$. This enables to devise an approach for inferring the resonance properties and negative permittivity effects from the optical measurements in the sub(THz) frequency band. Importantly, in the systems that we address, the typical thicknesses of the ferroelectric layer, $2a_{\mathrm{f}}\simeq 10\div 40\,\mathrm{nm}$ and the PDT period, $2d$, are both less then the wavelength of THz radiation, $\lambda\simeq 1\div 0.15\,\mathrm{mm}$ for $0.3\div 2\,\mathrm{THz}$. Hence the THZ wave will interact with the system hosting the PDT as with an ultra-thin dielectric film endowed with the anisotropic dielectric permittivity, the latter having the in- and out of plane components $\varepsilon_{\perp}$ and $\varepsilon_{\mathrm{f}}(\omega)$ respectively.

Reflection

We consider the situation where the incident THz beam gets reflected from the dielectric substrate, STO, with dielectric permittivity $\varepsilon_{\mathrm{p}}=300+10i$ in the sub-THz frequency range, passing twice through the PTO ferroelectric film with PDT. To ensure the interaction of the beam with the PDT polarization, the electric field of the light should be confined to the plane of incidence, hence the p-wave polarization geometry should be used. The experimental setup is shown in the inset in Fig. 4a. The calculated reflectance, $R$, of such a system is shown in the main panel of Fig. 4a as function of the angle of incidence, $\theta$, by the solid blue line. The plot illustrates a typical p-wave reflectance with the amplitude that nearly vanishes in the vicinity of the Brewster angle of the substrate, $\theta_{\mathrm{\scriptscriptstyle{B}}}\approx\arctan\left(\Re\mathrm{\,}\varepsilon_{\mathrm{p}}\right)^{1/2}\simeq 86.7^{\circ}$. For comparison, shown in Fig. 4a by the blue dashed line is also the net contribution of the substrate. Note that because of the very small thickness of the film, it is practically indistinguishable from the reflectance of the whole system. The difference between them, $\Delta R=R_{0}-R\simeq 10^{-3}$, is seen only under a magnification. Thus, to ensure the efficient filtering of the resonance PDT signal from the background substrate contribution, one has to focus instead on the relative variance of the signals difference, $\Delta R/R_{0}$, usually referred to as to reflectively and used to detect the tiny molecular mono-layers adsorbed on the substrate Tolstoy et al. (2003), see Fig. 4.

An advantage of using the reflectivity is that the value of the substrate reflectance, $R_{0}$, almost vanishes at $\theta_{\mathrm{\scriptscriptstyle{B}}}$ and the use of the ratio $\Delta R/R_{0}$ unravels the net contribution of the ultra-thin film per se. Shown by the solid red line, the reflectivity of our system does experience the drastic variation of about unity and clearly displays the Lorentz-like spike in vicinity of $\theta_{\mathrm{\scriptscriptstyle B}}$. This spike is a fingerprint of the PDT resonance and can be used for its experimental identification. The complete qualitative behavior of the resonance feature is illuminated by the three-dimensional color plot of $\Delta R/R_{0}$, as function of the frequency and the angle of incidence (Fig. 4b). The Lorentz-spike sharp behaviour gets smoothed upon the deviation from the resonance frequency.

The reflectancies $R$ and $R_{0}$ are calculated as a square of the magnitude of the corresponding Fresnel reflection coefficients, $R=|r^{p}_{123}|^{2}$ and $R_{0}=|r^{p}_{13}|^{2}$. Here index $p$ refers to the beam polarization and indices 1, 2, and 3 correspond to the ordering of the media: the air, the film, and the substrate, respectively.

The Fresnel reflection coefficient of the substrate-deposited film is given by the expression Tolstoy et al. (2003)

[TABLE]

where

[TABLE]

are the generalized Fresnel reflection coefficients between two adjacent media and $\beta=2\pi\left(2a_{f}/\lambda\right)\xi_{2}$ is the phase shift of the electromagnetic wave after a single pass through the film. The generalized indices of refraction, $\xi_{i}$, are defined via the angle of incidence, $\theta_{1}$, and dielectric permittivities of the media, $\varepsilon_{i}$, as

[TABLE]

While the refraction indices $\xi_{1}$ and $\xi_{3}$ represent the customary properties of the air, $\varepsilon_{1}=1$, and substrate, $\varepsilon_{3}=\varepsilon_{\mathrm{p}}$ endowed with the isotropic permittivities, the targeted electrodynamic response of PDT with negative capacitance is encoded in the anisotropic refraction index $\xi_{2}$ that is deduced from Schopper (1952). The latter grasps both the transversal permittivity of the film, $\varepsilon_{\perp}$, and the effective longitudinal permittivity of PDS, $\varepsilon_{2}=\varepsilon_{\mathrm{f}}(\omega)$.

The described approach identifies the PDT resonance behavior via detecting the characteristic Lorentz-like spike features in the angular and frequency dependencies of the reflectance. However, the developed theory of the electromagnetic response equips us with the extraordinary more advanced technique. Namely, the resonance point $\omega_{0}$ where $\varepsilon_{\mathrm{f}}(\omega)$ changes the sign, is found, with the great precision from the phase map of the Fresnel coefficient $r_{123}^{p}$. The latter (unlike the real-value reflectance $R=|r_{123}^{p}|^{2}$, where the phase information is lost) is a complex quantity which is measured, for example, by the phase-resolved ellipsometry technique Marsik et al. (2016).

To develop a feasible protocol for observation of the frequency domain endowed with negative $\varepsilon_{\mathrm{f}}$ let us describe the dependence $r_{123}^{p}(\theta)$ near the Brewster angle in the system with the ideally transparent substrate with $\Im\varepsilon_{\mathrm{p}}=0$. If the film were absent, the Fresnel coefficient would have been real and changed the sign upon passing zero at $\theta=\theta_{\mathrm{\scriptscriptstyle B}}$. Placing the film on the substrate, gives the contribution, albeit very small, that pulls $r_{123}^{p}$ into a complex plane, where the trajectory $r_{123}^{p}(\theta)$ circumflexes the zero from the positive- to negative semi-axis acquiring the phase $\pm\pi$. Equations (11)-(14) offering a full description of the electromagnetic response show that the direction (clock- vs. anticlock-wise) of the phase rotation is determined by the sign of $\varepsilon_{\mathrm{f}}$. The frequency $\omega_{0}$ at which $\varepsilon_{\mathrm{f}}$ changes the sign corresponds, therefore, to a particular point in the $(\theta,\omega)$-plane where phase acquires a factor of $2\pi$ when making a close loop around it, and the Fresnel coefficient, $r_{123}^{p}$, vanishes. These zeros of the complex function, $r_{123}^{p}(\theta,\omega)$, where the intensity of the reflected light is zero are often called ‘points of absolute darkness’ Kravets et al. (2013). They are stable with respect to small variations of parameters hence topologically protected, and we refer to them as to topological darkness points.

Figure 5a displays the 3D color plot of the phase and amplitude of the Fresnel reflection coefficient, $r_{123}^{p}$, of the PDT in the 30nm PTO film as a function of an incident angle and frequency for the case of transparent substrate with $\varepsilon_{\mathrm{p}}=300$, derived from Eqs. (11)-(14). Two topological dark points resulting from the negative permittivity are clearly visible at the manifold depicting $r_{123}^{p}(\theta,\omega)$. The higher-frequency point A corresponds to the discussed above sign change of $\varepsilon_{\mathrm{f}}$. The second, lower-frequency one, appears when $\varepsilon_{\mathrm{f}}\simeq-\sqrt{\varepsilon_{\mathrm{p}}\varepsilon_{\perp}}$. This manifestly illustrates the topological protection effect: the topological darkness points are not smoothed but merely shifted by dissipation. Figure 5b shows the the same phase plot for the PTO on a real substrate STO with $\varepsilon_{\mathrm{p}}=300+10i$. One sees that if the dissipation is sufficiently large, the darkness points merge and cease to exist, i.e. the sharp Brewster-angle singularities in the phase plot get smoothened.

The obtained results and recent measurements of the negative capacitance in PDT Zubko et al. (2016) open a new area of research in the field of the THz optics in ferroelectric materials. Exciting opportunities open in the area of the plasmonic epsilon-near-zero (ENZ) metamaterials having the unique property that the electromagnetic wave propagates with almost no phase advance. Although such materials have been made artificially in the microwave, visible and far-infrared spectral ranges Maas et al. (2013), the engineered ENZ structures in the THz spectral range have not been explored so far.

Acknowledgments

The work was supported by the U.S. Department of Energy, Office of Science, Materials Sciences and Engineering Division (VV) and European FP7-MC-ITN-NOTEDEV and H2020-MC-RISE-ENGIMA actions (IL).

APPENDIX

Calculations of the electrostatic properties of PDT in the applied field $\mathbf{D}$ generalize the zero-field Landau-Kittel calculations of PDT given in the textbook Landau and Lifshitz (2013) in application to magnetic domains. The geometry of the system is shown in Fig. 6. The $x$-axis is directed along the PDT texture, the $y$-axis is parallel to the domains and the $z$-axis is directed across the film. Since the polarization distribution does not change along the $y$-axis, the problem is reduced to the 2D one in the $xz$-plane. In the external field the up- and down-oriented polarization domains have different widths, $d_{+}$ and $d_{-}$, the field-induced polarization due to DW motion being expressed via the field dependence of the asymmetry factor $\eta=\left(d_{+}-d_{-}\right)/2d$ as $P_{\mathrm{dw}}=\eta(D)P_{\mathrm{s}}$, with $2d=d_{+}+d_{-}$ being the period of the structure. Following the relation (3) we express the dielectric permittivity of the PDT as:

[TABLE]

The linear dependence $\eta$ on $D$ is found by minimization of the electrostatic energy of the system, considering $\eta$ as variation parameter.

We start with the derivation of the spacial distribution of the electric fields inside the ferroelectric ($f$) and paraelectric ($p$) layers, $\mathbf{E}^{(p,f)}=(E_{x}^{(p,f)},E_{z}^{(p,f)})$, induced by the depolarization surface charges arising at the interfaces between $f$- and $p$- layers at $z=\pm a_{\mathrm{f}}$. These charges, appearing due to termination of the alternating polarization in PDT, are conventionally presented as $\sigma_{\pm}(x)=\pm\vartheta(x)P_{\mathrm{s}}$ where the sign-alternating function, $\vartheta$, is defined as $\vartheta(x)=\pm 1$ if $x\in d_{\pm}$. The corresponding electrostatic potentials in the $f$- and $p$-layers $\nabla\varphi^{(p,f)}=-\mathbf{E}^{(p,f)}$ satisfy the Laplace equations:

[TABLE]

The boundary conditions for the potential at the interfaces are fixed by the charges $\sigma_{\pm}(x)$:

[TABLE]

System (16) with boundary conditions (17)** ** is solved by the Fourier-series expansion method, taking into account that

[TABLE]

with $q_{n}=\pi n/d$ and

[TABLE]

After straightforward calculations we obtain:

[TABLE]

where $\varphi_{0}=\mp\left(\varepsilon_{0}\varepsilon_{\parallel}\right)^{-1}\left(D-\eta P_{\mathrm{s}}\right)a_{\mathrm{f}}$ is the average value of the potential at the upper/lower surface of the ferroelectric slab, and

[TABLE]

The electrostatic energy of the system can be calculated as the integral $\sum_{\pm}\int dx\int\left[\varphi(x)\right]_{z=\pm a_{\mathrm{f}}}d\sigma_{\pm}(x)\,$ over the surfaces of the ferroelectric layer. The integration over $\sigma_{\pm}$ reflects the self-consistent adiabatic charging of the surfaces, driven by the self-polarization of the system from the paraelectric state, until the state with the equilibrium spontaneous polarization, $\pm P_{\mathrm{s}}$, is achieved inside domains.

Taking into account the dependicies given by (20) of $\varphi(x)$ and $\sigma_{\pm}(x)$ on $P_{\mathrm{s}}$ we calculate the energy per unit length of PDT ( here the factor 2 is because of the two 2 sides of the slab):

[TABLE]

In the experimentally relevant limit $\sqrt{\varepsilon_{\perp}/\varepsilon_{\parallel}}\,a_{\mathrm{f}}/d\gg 1$, Eq. (APPENDIX) is simplified as:

[TABLE]

where

[TABLE]

Minimization of (23) with respect $\eta$ gives:

[TABLE]

Finally, making use (15) we obtain the dielectric permittivity of the PDT (4).

The DW displacements in PDT shown in Fig. 1d are given by $\pm x=\pm\eta d/2$. Equating the stiffness energy of the DW displacement, $2a_{\mathrm{f}}\cdot kx^{2}/2$, to (23) (taken in harmonic approximation at $D=0$ and multiplying by $d$, i.e. using the energy per a single DW), we obtain the expression for the stiffness coefficient, $k$, (7).

Note that the expression (4) can also be derived from the formulas given in Bratkovsky and Levanyuk (2001) by expressing the ratio $P_{a}/U$ in Eq. (31) of Bratkovsky and Levanyuk (2001) via the ratio $R/S$ defined in their Eq. (25). The Landau-Kittel formula (5) for the domain width, $d$, is obtained by minimization of $F_{\mathrm{el}}+F_{\mathrm{dw}}$ energy with respect $d$ at $\eta=0$. Here $F_{\mathrm{dw}}=(a_{f}/d)w_{\mathrm{dw}}$ is the  domain wall energy, calculated per unit length of PDT and $w_{\mathrm{dw}}$ is the surface energy of DW. The spatial distribution of the electric and polarization fields in the $p$- and $f$- slabs are found from (20) as: $\mathbf{E}^{(p,f)}=(E_{x}^{(f)},E_{z}^{(f)})=-\nabla\varphi^{(p,f)}$, $\mathbf{P}^{(p)}=\varepsilon_{0}\left(\varepsilon_{\mathrm{p}}-1\right)\mathbf{E}^{(p)}$, $\mathbf{P}^{(f)}=\left(P_{x}^{(f)},P_{z}^{(f)}\right)$ where $P_{x}^{(f)}=\varepsilon_{0}\left(\varepsilon_{\perp}-1\right)E_{x}^{(f)},P_{z}^{(f)}=\varepsilon_{0}\left(\varepsilon_{\parallel}-1\right)E_{z}^{(f)}+\vartheta(x)P_{\mathrm{s}}$. For the very thin films the field lines of $\mathbf{P}((r)$ form the closed loops, resembling the alternating vortex-antivortex patterns Yadav et al. (2016)