# Lamellar fluctuations melt ferroelectricity

G. G<EMAIL_ADDRESS>Centro de Investigación en Ciencia
e Ingeniería de Materiales, Universidad de Costa Rica, San José, Costa Rica
11501,  Escuela de Física, Universidad de Costa Rica, San José, Costa Rica
11501,  Department of Materials Science and Metallurgy, University of
Cambridge, Cambridge, UK CB3 0DS.  C. H. Liang  James Franck Institute,
University of Chicago, Chicago, Illinois, USA 60637,  Pritzker School of
Molecular Engineering, University of Chicago, Chicago, Illinois, USA 60637.
P. B<EMAIL_ADDRESS>James Franck Institute, University
of Chicago, Chicago, Illinois, USA 60637.

###### Abstract

We consider a standard Ginzburg-Landau model of a ferroelectric whose
electrical polarization is coupled to gradients of elastic strain [1, 2, 3,
4]. At the harmonic level, such flexoelectric interaction is known to
hybridize acoustic and optic phonon modes and lead to phases with modulated
lattice structures that precede the symmetry broken state for sufficiently
large couplings [5, 6, 7, 8, 9, 10]. Here, we use the we use the self-
consistent phonon approximation to calculate the effects of thermal and
quantum polarization fluctuations on the bare modes to show that such long-
range modulated order is unstable at all temperatures. We discuss the
implications for the nearly ferroelectric SrTiO3 and KTaO3, and propose that
these systems are melted versions of an underlying modulated state which is
dominated by non-zero momentum thermal fluctuations except at the very lowest
temperatures.

A ferroelectric (FE) material is defined by spontaneously broken inversion
symmetry, usually due to a small and coherent displacement of atoms within the
unit cell and hence a spontaneous electrical polarization ${\bm{P}}$ [11].
Because lattice displacements generate the order, the homogeneity of the
polarization depends on coupling between neighboring unit cells, and over
longer distances by a continuous elastic strain $\epsilon$. In a cubic (or
other high symmetry FE) the principal order parameter (OP) ${\bm{P}}$ belongs
to a different irreducible representation from the strain fields $\epsilon$ at
zero wavevector (${\bm{q}}=0$). This means that the leading order
electrostrictive coupling in a uniform Landau theory for the free energy is
bilinear in polarization and linear in strain (i.e. $\mathcal{O}\left(\epsilon
P^{2}\right)$. However, also there will be in general a flexoelectric coupling
[1, 2, 3, 4] (i.e. $\mathcal{O}\left(P\nabla\epsilon\right)$), and it is well
known that such couplings may give rise to modulated incommensurate phases
within a harmonic theory, provided that the coupling is big enough [5, 6, 8,
7, 9, 10]. Note that the flexoelectric coupling - being harmonic - is in
principle more relevant than the nonlinear electrostrictive coupling;
flexoelectricity will be the focus of this paper.

The materials SrTiO3 (STO) and KTaO3 (KTO) have long been cited as examples of
quantum paraelectrics (QPEs). In both cases the dielectric constant follows a
Curie-Weiss (CW) law divergence with a finite-temperature intercept, yet at
low temperatures saturates at extremely large values [12] and there is no FE
phase transition. The idea of quantum tunnelling of domains that in a mean-
field like (and massive) system have already grown to be very large before
quantum effects take over is disconcerting and has, over the years, generated
considerable interest in the general theory of a FE quantum critical point
[13, 14, 15, 16, 17, 18, 19, 20, 21].

Recent experiments, however, have revealed phase diagrams and lattice dynamics
in STO and KTO which cannot be explained with our current understanding of
these materials [22]. Specifically, high precision measurements of the
dielectric response have found several crossovers (specially seen under
pressure) [23, 24, 25, 26], and, very significantly, inelastic neutron
scattering experiments observe softening of a transverse acoustic (TA) mode at
a small but finite wavector in STO ($\simeq 0.025\,$rlu) [27] and at a larger
one in KTO ($\simeq 0.1\,$rlu) [6]. The latter results are consistent with
momentum-resolved electron energy-loss spectroscopy [28] and a previous
neutron scattering study [29] which show an unusual softening of the TA
branch, and with a first-principle simulation 333See Fig. 2 in Supplementary
of Ref. [58] which found a dip in the acoustic dispersion.

Soft acoustic phonons with a minimum at a finite wavector are characteristic
of systems approaching a structural instability corresponding to the onset of
a long-wavelength modulation [31, 32], and can be understood as the
consequence of level-repelling hybrid modes associated with a sufficiently
large coupling between a primary OP (e.g. polarization) and gradients of a
secondary OP (e.g. strain) [33]. However, harmonic and mean-field solutions of
such models [5, 6, 8, 7, 9, 10] predict the condensation of the low energy
branch at a finite wavevector ${\bm{q}}_{mod}$, thus leading to a modulated
phase which is not observed in STO and KTO. Nor can these solutions generate
crossovers.

In this letter, we invoke a well-known model of flexoelectricity [1, 2, 3, 4]
and, crucially, solve it within the self-consistent phonon approximation
(SCPA) in order to account for thermal and quantum fluctuations of
polarization of the disordered phase. We show that in the presence of purely
thermal fluctuations, the phase space for orientational fluctuations is large
enough that the modulated transition is suppressed to zero temperature, and
that by including zero-point quantum fluctuations, the correlation length will
remain finite even at absolute zero.

Figure 1: Phase boundaries in the (a) harmonic approximation and in the SCPA
(b) ignoring and (c) including quantum fluctuations of polarization. (d) SCPA
schematics of the temperature dependence of the zone-centre hybrid TO mode and
the (e) hybrid TA and TO phonon dispersions in the disordered phase. (f)
Schematic correlation function of polarization in the lamellar-like liquid
with correlation length $\xi$ (see Eq. (2)). QPE=quantum paraelectric.

We estimate the parameters in comparison to KTO and STO, and find reasonable
agreement with a model where the flexo-electric coupling is close to the
classical critical value for the onset of a modulated phase. Our analysis
suggests that the putative QPE phase of STO and KTO is largely explained by
classical thermal modulated fluctuations, and that there is a second crossover
to a regime dominated by quantum fluctuations at much lower temperature.
However, the quantum regime in our theory does not require macroscopic quantum
tunnelling, and instead arises more trivially due to zero-point occupation of
orientational zero-modes. In the light of these results, we propose that STO
and KTO are incipient modulated dielectrics with incommensurate order caused
by the coupling of polarization and strain inhomogeneity. In the underlying
ordered phase, the rotational symmetry of the crystal is broken as well as its
translational symmetry, but only in one direction (${\bm{q}}_{mod}$). This
would lead to intertwined stripe order parameters where the polarization and
elastic strain are periodic with a wavevector ${\bm{q}}_{mod}$ [31]. Such
stratified ordering resembles the structural organisation of the “lamellar” or
“smectic” phases in liquid crystals: a nearly periodic array of parallel
layers of polarizable molecules which do not possess long-range positional
order within a layer [34]. In a smectic the director vector ${\bm{n}}$ points
orthogonal to the planar lamellae, as does here the wavevector
${\bm{q}}_{mod}$. The underlying materials science is of course quite
different.

Lastly, we comment on lacunae in the model, the most important of which is the
neglect of crystal anisotropy and compositional disorder. We expect anisotropy
to lead to an ordered (striped) phase at low enough temperatures; disorder to
a glassy frozen version of the same. Both of these effects can suppress the
quantum regime. We also suggest experiments to further test our proposed
picture.

We consider an isotropic polarizable, elastic medium with a polarization order
parameter $P_{\alpha}\left({\bm{x}}\right)$ and a linear strain
$\epsilon_{\alpha\beta}({\bm{x}})=(1/2)\left(\partial u_{\alpha}/\partial
x_{\beta}+\partial u_{\beta}/\partial x_{\alpha}\right)$ as a secondary order
parameter in which $u_{\alpha}({\bm{x}})$ is the displacement field due to
long-wavelength acoustic phonons ($\alpha,\beta=1,2,3$).
$P_{\alpha}\left({\bm{x}}\right)$ is associated with a soft transverse optic
(TO) mode, the condensation of which leads to the FE transition. This leads us
to consider the GL Hamiltonian that is a sum of three terms: for the
polarization $H_{\text{pol}}$; the elastic modes $H_{\text{elastic}}$; and the
static flexoelectric coupling $H_{\text{flexo}}$ [2],

$\displaystyle H=H_{\text{pol}}+H_{\text{elastic}}+H_{\text{flexo}},$ (1)

where,

$\displaystyle H_{\text{pol}}$ $\displaystyle=\frac{1}{2}\int
d^{3}xd^{3}x^{\prime}\sum_{\alpha,\beta}U_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})P_{\alpha}\left({\bm{x}}\right)P_{\beta}\left({\bm{x}}^{\prime}\right)+u\int
d^{3}x\sum_{\alpha,\beta}P_{\alpha}^{2}\left({\bm{x}}\right)P_{\beta}^{2}\left({\bm{x}}\right),$
$\displaystyle H_{\text{elastic}}$ $\displaystyle=\frac{1}{2}\int
d^{3}x\sum_{\alpha,\beta,\gamma,\lambda}C_{\alpha\beta\gamma\lambda}\,\epsilon_{\alpha\beta}({\bm{x}})\epsilon_{\gamma\lambda}({\bm{x}}),$

and,

$\displaystyle H_{\text{flexo}}=\frac{1}{2}\int
d^{3}x\sum_{\alpha\beta\gamma\lambda}f_{\alpha\beta\gamma\lambda}\,\left[\epsilon_{\beta\gamma}({\bm{x}})\frac{\partial
P_{\alpha}({\bm{x}})}{\partial
x_{\lambda}}-P_{\alpha}({\bm{x}})\frac{\partial\epsilon_{\beta\gamma}({\bm{x}})}{\partial
x_{\lambda}}\right].$

Here,
$U_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})=r\delta_{\alpha\beta}\delta({\bm{x}}-{\bm{x}}^{\prime})+F_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})$
is the bare propagator of the optic phonons,
$F_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})$ is the dipole tensor with
Fourier transform
$F_{\alpha\beta}\left({\bm{q}}\right)=cq^{2}\delta_{\alpha\beta}+g_{0}\left(q_{\alpha}q_{\beta}/q^{2}\right)-h_{0}q_{\alpha}q_{\beta}$
[35]. $C_{\alpha\beta\gamma\delta}$ is the elastic constant tensor in units of
energy per unit volume, $f_{\alpha\beta\gamma\lambda}$ is the flexocoupling
tensor, and $r=r_{0}(T-T_{0})$ where $T_{0}$ is the mean field transition
temperature to a homogeneous FE state at $f_{\alpha\beta\gamma\lambda}=0$.

In the absence of flexoelectric coupling and non-linearities $(u=0)$, the
phonon excitations of $H$ are well-known: there are two doubly-fold degenerate
TO and TA modes with (squared) frequencies
$\omega_{\text{TO}}^{2}({\bm{q}})=r+cq^{2}$ and
$\omega_{\text{TA}}^{2}({\bm{q}})=(\hbar^{2}C_{s}/\rho)q^{2}$, respectively,
and two singly-degenerate longitudinal optic and acoustic modes with (squared)
frequencies $\omega_{\text{LO}}^{2}({\bm{q}})=r+g_{0}+(c+h_{0})q^{2}$ and
$\omega_{\text{LA}}^{2}({\bm{q}})=(\hbar^{2}(C_{v}+(4/3)C_{s})/\rho)q^{2}$,
respectively. $C_{v}$ and $C_{s}$ are the bulk and shear elastic moduli,
respectively, and $\rho$ is the density of the solid.

In this paper, we focus on the effects on the hybrid transverse modes only as
the level repulsion is significantly weaker in the longitudinal branches due
to the characteristically large depolarizing fields of FEs which gap the
longitudinal optic excitations. Moreover, it is shown in the Supplementary
Note 1 (SN1) that in isotropic media the flexoelectric interaction does not
mix the transverse and longitudinal excitations.

We now discuss our results. Figure 1 (a) shows schematics of the phase
boundary (PB) computed from the condensation of the bare modes (i.e., $u=0$)
in the parent phase (see SN1). For small shears of the flexocoupling tensor
(i.e., $f_{s}<f_{cr}=\sqrt{c\,C_{s}}$), the soft zone-center TO mode condenses
at a transition temperature which is independent of $f_{s}$ and equal to
$T_{0}$, thus leading to a FE transition. For large flexoelectric shears
(i.e., $f_{s}>f_{cr}$), a soft minimum develops at a non-zero wave-vector
${\bm{q}}_{mod}$ of the TA branch, which upon condensation at a temperature
$T_{mod}>T_{0}$, would lead to the modulated phase with coupled strain and
polarization order parameters that are periodic with a wavevector
${\bm{q}}_{mod}$ [31]. The zone-center TO mode remains gaped with energy
squared $\tau^{0}_{mod}\equiv r_{0}(T_{mod}-T_{0})$. At $f_{s}=f_{cr}$,
${\bm{q}}_{mod}=0$ and both the TO and TA modes condense $T_{mod}=T_{0}$.

We now assess the stability of the bare PB. Near above $T_{0}$ and $T_{mod}$
and in the classical limit, the local correlation functions of polarization
are approximately given as follows (see SN2),

$\displaystyle\left<{\bm{P}}^{2}\right>_{0}\propto\begin{cases}1-b\left(T-T_{0}\right)^{1/2},&f_{s}<f_{cr},\\\
\left(T-T_{0}\right)^{-1/4},&f_{s}=f_{cr},\\\
\left(T-T_{mod}\right)^{-1/2},&f_{s}>f_{cr},\end{cases}$

where $\left<...\right>_{0}$ denotes thermal average at the Gaussian level and
$b$ is a constant. The PB is thus strongly modified by the inclusion of
fluctuations: while the FE transition prevails at $T_{0}$, the long-range
modulated phase change is unstable, except at absolute zero where
$\left<{\bm{P}}^{2}\right>_{0}=0$.

We now consider non-linearities (i.e., $u>0$) in the SCPA (see SN3). We first
discuss the impact of thermal fluctuations alone and then introduce quantum
effects. Fig. 1 (b) shows the purely classical PB. For $f_{s}<f_{cr}$, the
coupling to finite momentum elastic fluctuations suppresses the FE transition
temperature and leads to a thermal regime where the squared of the soft zone-
center TO phonon frequency $\tau$ (or equivalently the inverse dielectric
constant) separates from the CW law behavior, as shown in Fig. 1 (d). For
$f_{s}>f_{cr}$, the bare PB turns into a crossover. This is illustrated best
by the SCPA dispersions shown in Fig. 1 (d). As mentioned above, the bare
coupled acoustic mode softens to zero at a non-zero momentum ${\bm{q}}_{mod}$.
However, in an isotropic theory, such as ours, only the modulus of
${\bm{q}}_{mod}$ is determined but not the direction. The phase space for
fluctuations in the direction of ${\bm{q}}_{mod}$ lies on the surface of a
sphere of one dimension less than the physical dimension. This all but
guarantees that thermal fluctuations at any finite temperature will disorder
the direction, so there will be no long-range stripe order. This is very
analogous to fluctuating laminar or nematic soft matter systems [34].
Consequently, the harmonic transition in Fig. 1 (a) becomes a crossover in
Fig. 1 (b), and the acoustic excitations in the vicinity of ${\bm{q}}_{mod}$
acquire a gap, as seen in Fig. 1 (e). The polarization correlations of the
melted phase are sinusoidal functions attenuated by an exponentially decaying
envelope with a correlation length
$\xi\propto\left(\tau-\tau^{0}_{mod}\right)^{-1/4}$ longer than the period of
the modulation, i.e.,

$\displaystyle\left<{\bm{P}}({\bm{x}})\cdot{\bm{P}}({\bm{0}})\right>\sim
e^{-x/(2\,q_{mod}\,\xi^{2})}\,\sin\left(q_{mod}\,x\right)/\left(q_{mod}\,x\right),~{}~{}~{}~{}~{}~{}(q_{mod}\,\xi\gg
1),$ (2)

where $x=\left|{\bm{x}}\right|$ and $q_{mod}=\left|{{\bm{q}}_{mod}}\right|$
(see SN3). This is schematically shown in Fig. 1 (f). Upon cooling, $\xi$
would grow as $\tau$ approaches $\tau_{mod}^{0}$ (Fig. 1 (d)) and diverge at
absolute zero, thus restoring the long-range stripe order, as shown in Fig. 1
(b).

We now asses the effects of quantum polarization fluctuations. In our model,
the nonlinearity of the theory is controlled by the quartic coefficient $u$,
as shown in Fig. 2. For $f_{s}=0$, the FE transition temperature drops
approximately linearly in $u$, vanishing at a quantum critical point. For
$0<f_{s}<f_{cr}$, the critical nonlinearity is reduced, and there is a much
more substantial parameter regime where quantum fluctuations alone destroy the
transition and lead to the standard QPE phase, as it is shown in Fig. 1 (c).
This is consistent with the much larger phase space for orientational
fluctuations when the long-range modulated order would be at finite $q$.

Figure 2: (a) Dependence of the ferroelectric transition temperature $T_{c}$
on the flexoelectric shear coupling constant $f_{s}$ and the quartic
anharmonic coefficient $u$ of the polarization. (b)-(c) Cross-sections of (a)
in the $T_{c}-u$ and $T_{c}-f_{s}$ plane, respectively. $u_{0}$ is the value
of $u$ at $T_{c}=0\,$K and $f_{s}=0$.

For $f_{s}>f_{cr}$, the zero-point fluctuations keep the energy gap at
${\bm{q}}_{mod}$ open even at absolute zero, thus destabilising the classical
$T=0\,$K transition shown in Fig. 1 (b). In addition, the temperature at which
the quantum fluctuations become relevant marks another crossover in the
spectrum, where the correlation length will saturate. This occurs when
$k_{B}T$ is comparable to the energy of the excitations around
${\bm{q}}_{mod}$, as shown in Fig. 1 (e). So, for $f_{s}>f_{cr}$ we predict
that there will be two crossovers as $T$ is lowered, and no long-range stripe
order, as shown in Fig. 1 (c). We stress that when zero-point quantum
fluctuations are included, the only ordered phase that we find is
ferroelectricity. But outside the envelope shown in Figure 2 (a) and for
$f_{s}>f_{cr}$, we expect the lamellar-like fluctuations to have a very long-
range correlation.

Figure 3: (a),(b) Calculated temperature dependence of the zone center hybrid TO mode compared to neutron [36, 37, 38] and hyper-Raman [39] scattering experiments . (c),(d) Calculated dispersion of the hybrid TA and TO phonons (solid lines) compared to neutron scattering experiments (circles and diamonds) [27, 38]. Dashed line in (d) is the bare linear TA dispersion for $f_{s}=0$ with a sound velocity of 4800 m/s (as shown in Ref. [27] for comparison). Model parameters are given in Table 1. Table 1: Model parameters. |  KTO |  STO
---|---|---
$r_{0}\,$[meV2 K-1] |  $0.38$222From Ref. [40]. |  $0.47$555From Ref. [42].
$T_{0}\,$[K] |  $4.54$222From Ref. [40]. |  $35.5$555From Ref. [42].
$u\,$ [$10^{-4}$ meV3 rlu-3] |  $7.19$333From Ref. [41]. |  $7.61$555From Ref. [42].
$c\,$[$10^{3}$ meV2 rlu-2] |  $4.828$444From Ref. [38]. | $13.5$
$a_{t}$111$a_{t}=\hbar^{2}C_{s}/\rho$ and $v_{t}=\hbar f_{s}/\sqrt{\rho}$. [$10^{3}\,$ meV2 rlu-2] |  $1.553$444From Ref. [38]. |  $2.739$666From Ref. [43].
$v_{t}$111$a_{t}=\hbar^{2}C_{s}/\rho$ and $v_{t}=\hbar f_{s}/\sqrt{\rho}$. [$10^{3}\,$ meV2 rlu-2] | $3.067$ | $6.737$
$h\,$ [$10^{4}$ meV4 rlu-4] | $6.15$ | $406$
$v_{t}/v_{cr}=f_{s}/f_{cr}$ | $1.12$ | $1.17$

We now turn to comparisons to STO and KTO. Table 1 gives the values for the
model parameters, which for the most part were taken from the literature while
others were obtained from experimental data, e.g., $f_{s}$ was fitted to the
reported TA anomalies ($0.025\,$rlu for STO [27], and $0.1\,$rlu for KTO [6]).
The details can be found in SN4. For both materials, we find agreement between
the calculated and measured temperature dependence of the TO frequency, as
shown in Figs. 3 (a) and 3 (b). We note the apparent kink in the calculated
frequencies is a consequence of the small values of $u$ in our
paremetrization. As shown in SN3, $\tau$ is a smooth function of $T$ at all
temperatures. Figures 3 (c) and 3 (d) show the calculated and observed phonon
dispersion curves. While our model and parametrization clearly overestimate
the repulsion between the TA and TO branches, it generates a minimum at a
finite $q$ in the hybrid acoustic dispersion while preventing symmetry
breaking down to the lowest temperatures, as observed in experiments [27, 6].
We stress that no such behavior occurs for $f_{s}\leq f_{cr}$, nor could we
obtained a physically reasonable parametrization in this regime, as explained
in SN4. We attribute the quantitative discrepancies to ignoring crystal
anisotropy (e.g. low temperature STO is tetragonal while KTO is cubic),
broadening, non-zero lifetimes, electrostriction, and coupling to other
lattice modes (e.g., rotations of the TiO6 octahedra). Our fits place both
materials in the liquid phase of Fig. 1 (c) ($f_{s}/f_{cr}=1.17$ for STO and
$f_{s}/f_{cr}=1.12$ for KTO). Of course, these values should be taken as rough
estimates, as our simplified model, approximate solution, and currently
available experimental data prevent us from determining $f_{s}$ precisely.
Nonetheless, they show that both materials are near a modulated instability.

This point of course connects to a long debate about the nature of the FE
phase transition as “displacive” or “order-disorder” [44, 45]. Polar
nanodomains are clearly seen at room temperature in strained (and therefore
low-temperature ferrolectric) films of STO though their presence has been
attributed to long-range Coulomb forces rather than strain [46]. Nuclear
magnetic resonance (NMR) is a technique that is sensitive to local charge
disproportionation and electric field gradients, and was useful in identifying
charge-density waves in transition metal dichalchogenides, even when not fully
ordered [47, 48]. NMR measurements on STO [49] have been interpreted in terms
of dynamic fluctuations of the Ti-ion into off-center sites coupled in a
biased way to elongations of the unit cell; however the length scale of this
correlation could not be determined and the measurements did not go below
$25\,$K. Data at a lower temperature would be helpful. Ab initio calculations
of STO [50] concluded that quantum lattice fluctuations are necessary to
stabilise the paraelectric phase - but spatially modulated phases were not
considered in that work. Modulated phases could also be in principle seen in
ab initio calculations, albeit using very large supercells. In that context we
note recent calculations [51, 52] which predict that STO and cubic BaTiO3
should be polymorphic at $T=0\,$K, with a disordered arrangement of off-site
Ti distortions. A disordered ground state at absolute zero violates the third
law, and perhaps is associated instead with the failure to match the supercell
to the natural incommensurate period.

We now discuss our simplifying approximations. We have ignored crystal
anisotropy, in which case the orientational modes will have a preferred
direction. Once the temperature becomes considerably smaller than the
anisotropy in the spectrum, the phase space volume of critical fluctuations
will be reduced, which may lead to the stabilization of an ordered modulated
phase which could precede the establishment of an ordered FE phase (this is of
course the standard situation seen, for example, in quartz [53]). Similarly,
disorder will break the $\mathcal{O}(3)$ symmetry and should be expected to
lead to a glassy phase where orientational fluctuations are frozen. This is a
familiar situation in, e.g. charge-density-wave systems, where impurities and
disorder pin its order parameter and generate a strongly nonlinear response to
applied electric fields [54, 55]. Our conclusion that STO should be near a
modulated instability is in stark contrast with ab initio calculations [56],
which suggest that $f_{s}/f_{cr}\simeq 0.7$ at most. We note, however, that we
are using a phenomenological continuum theory where the parameters are
determined from experiment (but nonetheless can make explicit predictions
because the number of parameters is small). We cannot assert that the
parameters can be directly calculated from harmonic microscopic theory because
we have already coarse-grained to a scale of many unit cells, so that, for
example, twin boundaries and defect structures have been integrated over [57].
Therefore, the difference between parameters obtained from macroscopic fits
and those determined by microscopic calculations hints at some mesoscale
physics to be unearthed.

Our self-consistent harmonic solution of the model implicitly assumes that in
the melted phases the amplitudes of fluctuations are small. Beyond this
regime, electrostrictive couplings (effectively quartic in the primary order
parameter) will be important. These of course enhance the anisotropy (at least
at very long length scales) but also effectively soften the sound velocity and
hence enhance the flexoelectric effect. In a strong coupling theory, the model
would become one of quasi-periodic FE domains separated by fluctuating, sharp
domain walls. In a better theory, the modes will acquire a linewidth, as seen
in experiment. This is an interesting matter for future exploration, but we
expect that the qualitative features of our theory should survive.

We conclude that a strong enough flexoelectric coupling and purely thermal
fluctuations of polarization can mimic quantum paraelectricity by saturating
the low-temperature dielectric constant. We explain that in an isotropic
system, the divergent orientational fluctuations will suppress ordering even
at low temperatures so that one should expect a series of crossovers upon
cooling: from a high-temperature CW law; to a classical fluctuating lamellar
phase; to a quantum lamellar phase; and finally to a striped phase when
crystal anisotropy becomes important. Such a sequence of crossovers is loosely
consistent with experiments in STO [23, 24] and the soft finite $q$ modes have
also been observed [27, 28]. Coak et al.[23] ventured that a loss peak seen in
their data might be ascribed to quantum fluctuations of domain walls. We also
stress that the quantum regime of our theory does not require tunneling of
heavy atoms through potential barriers, but is simply zero-point physics in a
quasi-harmonic well. With an albeit simplified theory, our parameters are
physically reasonable and we believe our model provides a fair description of
the disordered state in STO and KTO.

In addition to the NMR experiments suggested above to further test our
proposed scenario, we point out that measurements of the dynamical structure
factor $S({\bm{q}},\omega)$ would detect fluctuations associated with the
melting of the modulated phase. $S({\bm{q}},\omega)$ has of course been
measured by neutrons and X-rays in STO [37, 29, 58] and KTO [6, 38], but more
data is needed in order to clearly resolve any dips along their acoustic
branches. Finally, we point out that our model could be adapted to investigate
flexoelectricity in other systems such as recently synthesized complex oxide
membranes [59], has a relation to structural Leggett modes in pyrochlore
relaxors [60], and could be relevant to metallicity and superconductivity in
doped STO [61].

###### Acknowledgements.

We acknowledge helpful feedback from Benoît Fauqué, Massimiliano Stengel,
Peter Abbamonte, Caitlin Kengle, Michael Norman, Xing He, and Siddharth
Saxena. C.H.L. acknowledges support from NSF GRFP DGE-1746045, NSF INTERN, and
DOE SCGSR. G.G.G.-V. acknowledges support from the Vice-rectory for Research
at the University of Costa Rica under project No. 816-C1-601.

## References

  * Zubko  _et al._ [2013] P. Zubko, G. Catalan, and A. K. Tagantsev, Flexoelectric effect in solids, Ann. Rev. Mater. Res. 43, 387 (2013).
  * Yudin and Tagantsev [2013] P. V. Yudin and A. K. Tagantsev, Fundamentals of flexoelectricity in solids, Nanotechnology 24, 432001 (2013).
  * Tagantsev and Yudin [2016] A. K. Tagantsev and P. V. Yudin, eds.,  _Flexoelectricity in Solids: From Theory to Applications_ (World Scientific, 2016).
  * Wang  _et al._ [2019] B. Wang, Y. Gu, S. Zhang, and L.-Q. Chen, Flexoelectricity in solids: Progress, challenges, and perspectives, Prog. Mater. Sci. 106, 100570 (2019).
  * Vaks [1968] V. G. Vaks, Phase transitions of the displacement type in ferroelectrics, J. Exp. Th. Phys. 27, 486 (1968).
  * Axe  _et al._ [1970] J. D. Axe, J. Harada, and G. Shirane, Anomalous acoustic dispersion in centrosymmetric crystals with soft optic phonons, Phys. Rev. B 1, 1227 (1970).
  * Aslanyan and Levanyuk [1980] T. A. Aslanyan and A. P. Levanyuk, Possibility of formation if an incommensurate structure due to fluctuation effects, Sov. Phys. Solid State 22, 1262 (1980).
  * Heine and McConnell [1981] V. Heine and J. D. C. McConnell, Origin of modulated incommensurate phases in insulators, Phys. Rev. Lett. 46, 1092 (1981).
  * Yudin  _et al._ [2014] P. V. Yudin, R. Ahluwalia, and A. K. Tagantsev, Upper bounds for flexoelectric coefficients in ferroelectrics, Appl. Phys. Lett. 104, 082913 (2014).
  * Morozovska  _et al._ [2016] A. N. Morozovska, E. A. Eliseev, C. M. Scherbakov, and Y. M. Vysochanskii, Influence of elastic strain gradient on the upper limit of flexocoupling strength, spatially modulated phases, and soft phonon dispersion in ferroics, Phys. Rev. B 94, 174112 (2016).
  * Strukov and Levanyuk [1998] B. A. Strukov and A. P. Levanyuk,  _Ferroelectric Phenomena in Crystals: Physical Foundations_ (Springer, 1998).
  * Müller and Burkard [1979] K. A. Müller and H. Burkard, SrTiO3: An intrinsic quantum paraelectric below 4 K, Phys. Rev. B 19, 3593 (1979).
  * Müller  _et al._ [1991] K. Müller, W. Berlinger, and E. Tosatti, Indication for a novel phase in the quantum paraelectric regime of SrTiO3, Z. Physik B - Condensed Matter 84, 277 (1991).
  * Martoňák [1991] R. Martoňák,  _New Phenomena in the Quantum Paraelectric State of SrTiO 3_, Master’s thesis, International School for Advanced Studies (1991).
  * Tosatti and Martoňák [1994] E. Tosatti and R. Martoňák, Rotational melting in displacive quantum paraelectrics, Solid State Communications 92, 167 (1994).
  * Martoňák and Tosatti [1994] R. Martoňák and E. Tosatti, Path-integral monte carlo study of a model two-dimensional quantum paraelectric, Phys. Rev. B 49, 12596 (1994).
  * Zhong and Vanderbilt [1996] W. Zhong and D. Vanderbilt, Effect of quantum fluctuations on structural phase transitions in SrTiO3 and BaTiO3, Phys. Rev. B 53, 5047 (1996).
  * Roussev and Millis [2003] R. Roussev and A. J. Millis, Theory of the quantum paraelectric-ferroelectric transition, Phys. Rev. B 67, 014105 (2003).
  * Bussmann-Holder  _et al._ [2007] A. Bussmann-Holder, H. Büttner, and A. R. Bishop, Polar-soft-mode-driven structural phase transition in SrTiO3, Phys. Rev. Lett. 99, 167603 (2007).
  * Pálová  _et al._ [2009] L. Pálová, P. Chandra, and P. Coleman, Quantum critical paraelectrics and the casimir effect in time, Phys. Rev. B 79, 075101 (2009).
  * Conduit and Simons [2010] G. J. Conduit and B. D. Simons, Theory of quantum paraelectrics and the metaelectric transition, Phys. Rev. B 81, 024102 (2010).
  * Chandra  _et al._ [2017] P. Chandra, G. G. Lonzarich, S. E. Rowley, and J. F. Scott, Prospects and applications near ferroelectric quantum phase transitions: a key issues review, Rep. Prog. Phys. 80, 112502 (2017).
  * Coak  _et al._ [2018] M. J. Coak, C. R. S. Haines, C. Liu, D. M. Jarvis, P. B. Littlewood, and S. S. Saxena, Dielectric response of quantum critical ferroelectric as a function of pressure, Sci. Rep. 8, 14936 (2018).
  * Coak  _et al._ [2020] M. J. Coak, C. R. S. Haines, C. Liu, S. E. Rowley, G. G. Lonzarich, and S. S. Saxena, Quantum critical phenomena in a compressible displacive ferroelectric, Proc. Natl. Acad. Sci. 117, 12707 (2020).
  * Coak  _et al._ [2019] M. J. Coak, C. R. S. Haines, C. Liu, G. G. Guzmán-Verri, and S. S. Saxena, Pressure dependence of ferroelectric quantum critical fluctuations, Phys. Rev. B 100, 214111 (2019).
  * Rowley  _et al._ [2014] S. E. Rowley, L. J. Spalek, R. P. Smith, M. P. M. Dean, M. Itoh, J. F. Scott, G. G. Lonzarich, and S. S. Saxena, Ferroelectric quantum criticality, Nature Phys. 10, 367 (2014).
  * Fauqué  _et al._ [2022] B. Fauqué, P. Bourges, A. Subedi, K. Behnia, B. Baptiste, B. Roessli, T. Fennell, S. Raymond, and P. Steffens, Mesoscopic fluctuating domains in strontium titanate, Phys. Rev. B 106, L140301 (2022).
  * Kengle  _et al._ [2022] C. S. Kengle, S. I. Rubeck, M. Rak, J. Chen, F. Hoveyda, S. Bettler, A. Husain, M. Mitrano, A. Edelman, P. Littlewood, T.-C. Chiang, F. Mahmood, and P. Abbamonte, Anharmonic multiphonon origin of the valence plasmon in SrTi1-xNbxO3 (2022), arXiv:2210.15026 .
  * Vacher  _et al._ [1992] R. Vacher, J. Pelous, B. Hennion, G. Coddens, E. Courtens, and K. A. Müller, Anomalies in the TA-Phonon Branches of SrTiO3 in the Quantum Paraelectric Regime, Europhys. Lett. 17, 45 (1992).
  * Note [1] See Fig. 2 in Supplementary of Ref. [58].
  * Blinc and Levanyuk [1986] R. Blinc and A. P. Levanyuk, eds.,  _Incommesurate phases in dielectrics_, Vol. 1 & 2 (Elsevier, 1986).
  * Scott [1983] J. F. Scott, Acoustic-phonon dispersion at incommensurate phase transitions, Ferroelectrics 47, 33 (1983).
  * Heine and McConnell [1984] V. Heine and J. D. C. McConnell, The origin of incommensurate structures in insulators, J. Phys. C: Solid State Physics 17, 1199 (1984).
  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Aharony and Fisher [1973] A. Aharony and M. E. Fisher, Critical behavior of magnets with dipolar interactions. i. renormalization group near four dimensions, Phys. Rev. B 8, 3323 (1973).
  * Shirane  _et al._ [1967] G. Shirane, R. Nathans, and V. J. Minkiewicz, Temperature dependence of the soft ferroelectric mode in KTaO3, Phys. Rev. 157, 396 (1967).
  * Yamada and Shirane [1969] Y. Yamada and G. Shirane, Neutron scattering and nature of the soft optical phonon in SrTiO3, J. Phys. Soc. Jpn. 26, 396 (1969).
  * Farhi  _et al._ [2000] E. Farhi, A. K. Tagantsev, R. Currat, B. Hehlen, E. Courtens, and L. A. Boatner, Low energy phonon spectrum and its parameterization in pure KTaO3 below 80 K, Eur. Phys. J. B. 15, 615 (2000).
  * Vogt [1995] H. Vogt, Refined treatment of the model of linearly coupled anharmonic oscillators and its application to the temperature dependence of the zone-center soft-mode frequencies of ${\mathrm{KTaO}}_{3}$ and ${\mathrm{SrTiO}}_{3}$, Phys. Rev. B 51, 8046 (1995).
  * Fujishita  _et al._ [2016] H. Fujishita, S. Kitazawa, M. Saito, R. Ishisaka, H. Okamoto, and T. Yamaguchi, Quantum paraelectric states in SrTiO3 and KTaO3: Barrett model, Vendik model, and quantum criticality, J. Phys. Soc. Jpn. 85, 074703 (2016).
  * Skoromets  _et al._ [2016] V. Skoromets, C. Kadlec, H. Němec, D. Fattakhova-Rohlfing, and P. Kužel, Tunable dielectric properties of KTaO3 single crystals in the terahertz range, J. Phys. D: Appl. Phys. 49, 065306 (2016).
  * Rabe  _et al._ [2007] K. M. Rabe, C. H. Ahn, and J.-M. Triscone, eds.,  _Physics of Ferroelectrics A Modern Perspective_ (Springer-Verlag, 2007).
  * Carpenter [2007] M. A. Carpenter, Elastic anomalies accompanying phase transitions in (Ca,Sr)TiO3 perovskites: Part I. Landau theory and a calibration for SrTiO3, Am. Mineral 92, 309 (2007).
  * Jona and Shirane [1993] F. Jona and G. Shirane,  _Ferroelectric Crystals_ (Dover, 1993).
  * Comès  _et al._ [1970] R. Comès, M. Lambert, and A. Guinier, Désordre linéaire dans les cristaux (cas du silicium, du quartz, et des pérovskites ferroélectriques), Acta Crys. A26, 244 (1970).
  * Salmani-Rezaie  _et al._ [2020] S. Salmani-Rezaie, K. Ahadi, W. M. Strickland, and S. Stemmer, Order-disorder ferroelectric transition of strained SrTiO3, Phys. Rev. Lett. 125, 087601 (2020).
  * Stiles and Williams [1976] J. A. R. Stiles and D. L. Williams, NMR measurements of phenomena related to charge density waves in NbSe2, J. Phys. C: Solid State Phys. 9, 3941 (1976).
  * Pfeiffer  _et al._ [1982] L. Pfeiffer, R. E. Walstedt, R. F. Bell, and T. Kovacs, Temperature dependence of the orthorhombic charge-density-wave parameters in $2H$-TaSe2 by 77Se nuclear magnetic resonance, Phys. Rev. Lett. 49, 1162 (1982).
  * Zalar  _et al._ [2005] B. Zalar, A. Lebar, J. Seliger, R. Blinc, V. V. Laguta, and M. Itoh, NMR study of disorder in BaTiO3 and SrTiO3, Phys. Rev. B 71, 064107 (2005).
  * Shin  _et al._ [2021] D. Shin, S. Latini, C. Schäfer, S. A. Sato, U. De Giovannini, H. Hübener, and A. Rubio, Quantum paraelectric phase of SrTiO3 from first principles, Phys. Rev. B 104, L060103 (2021).
  * Zhao  _et al._ [2021] X.-G. Zhao, Z. Wang, O. I. Malyi, and A. Zunger, Effect of static local distortions vs. dynamic motions on the stability and band gaps of cubic oxide and halide perovskites, Materials Today 49, 107 (2021).
  * Zhao  _et al._ [2022] X.-G. Zhao, O. I. Malyi, S. J. L. Billinge, and A. Zunger, Intrinsic local symmetry breaking in nominally cubic paraelectric BaTiO3, Phys. Rev. B 105, 224108 (2022).
  * Aslanyan and Levanyuk [1979] T. A. Aslanyan and A. P. Levanyuk, On the possibility of the incommensurate phase near $\alpha\leftrightarrow\beta$ transition point in quartz, Solid State Commun. 31, 547 (1979).
  * Fukuyama and Lee [1978] H. Fukuyama and P. A. Lee, Dynamics of the charge-density wave. I. Impurity pinning in a single chain, Phys. Rev. B 17, 535 (1978).
  * Lee and Rice [1979] P. A. Lee and T. M. Rice, Electric field depinning of charge density waves, Phys. Rev. B 19, 3970 (1979).
  * Stengel [2016] M. Stengel, Unified ab initio formulation of flexoelectricity and strain-gradient elasticity, Phys. Rev. B 93, 245107 (2016).
  * Pesquera  _et al._ [2018] D. Pesquera, M. A. Carpenter, and E. K. H. Salje, Glasslike dynamics of polar domain walls in cryogenic SrTiO3, Phys. Rev. Lett. 121, 235701 (2018).
  * He  _et al._ [2020] X. He, D. Bansal, B. Winn, S. Chi, L. Boatner, and O. Delaire, Anharmonic eigenvectors and acoustic phonon disappearance in quantum paraelectric SrTiO3, Phys. Rev. Lett. 124, 145901 (2020).
  * Harbola  _et al._ [2021] V. Harbola, S. Crossley, S. S. Hong, D. Lu, Y. A. Birkhölzer, Y. Hikita, and H. Y. Hwang, Strain gradient elasticity in SrTiO3 membranes: bending versus stretching, Nano Lett. 21, 2470 (2021).
  * Meier  _et al._ [2022] Q. N. Meier, D. Hickox-Young, G. Laurita, N. A. Spaldin, J. M. Rondinelli, and M. R. Norman, Leggett modes accompanying crystallographic phase transitions, Phys. Rev. X 12, 011024 (2022).
  * Rischau  _et al._ [2017] C. W. Rischau, X. Lin, C. P. Grams, F. Finck, S. Harms, J. Engelmayer, T. Lorenz, T. Gallais, B. Fauqué, J. Hemberger, and K. Behnia, A ferroelectric quantum phase transition inside the superconducting dome of Sr1-xCaxTiO3-δ, Nat. Phys. 13, 643 (2017).

Supplementary Information

G. G. Guzmán-Verri1,2,3, C. H. Liang4,5, P. B. Littlewood4

1Centro de Investigación en Ciencia e Ingeniería de Materiales (CICIMA),
Universidad de Costa Rica, San José, Costa Rica 11501,
2Escuela de Física, Universidad de Costa Rica, San José, Costa Rica 11501,

3Department of Materials Science and Metallurgy, University of Cambridge,
Cambridge, UK CB3 0DS,

4James Franck Institute, University of Chicago, 929 E 57 St, Chicago,
Illinois, USA 60637,
5Pritzker School of Molecular Engineering, University of Chicago, 5640 S Ellis
Ave, Chicago, Illinois, USA 60637.

Contents

  * Note 1. Bare hybrid modes.

  * Note 2. Phase stability.

  * Note 3. Self-consistent phonon approximation.

  * Note 4. Fits to experiments in STO and KTO.

  * References.

##  Supplementary Note 1. Bare hybrid modes

In the absence of non-linearities ($u=0$), we rewrite the GL Hamiltonian of
Eq. (1) in momentum space and in terms of the displacement field,

$H=\fint\frac{d^{3}q}{(2\pi)^{3}}\,\sum_{\alpha\beta}U_{\alpha\beta}({\bm{q}})P_{\alpha}({\bm{q}})P^{*}_{\beta}({\bm{q}})+\sum_{\alpha\beta}M_{\alpha\beta}({\bm{q}})u_{\alpha}({\bm{q}})u^{*}_{\beta}({\bm{q}})\\\
+\sum_{\alpha\beta}V_{\alpha\beta}({\bm{q}})\left[u_{\alpha}({\bm{q}})P_{\beta}^{*}({\bm{q}})+u_{\alpha}^{*}({\bm{q}})P_{\beta}({\bm{q}})\right],$
(S1)

where
$M_{\alpha\beta}({\bm{q}})=\sum_{\gamma\lambda}C_{\gamma\alpha\lambda\beta}q_{\gamma}q_{\lambda},$
is the bare dynamical matrix of the acoustic phonons, and
$V_{\alpha\beta}({\bm{q}})=\sum_{\gamma\lambda}f_{\alpha\beta\gamma\lambda}q_{\gamma}q_{\lambda}$
is the flexoelectric interaction. $\fint d^{3}q$ indicates to integrate over
half reciprocal space only.

For isotropic media [1, 2],

$\displaystyle
U_{\alpha\beta}\left({\bm{q}}\right)=\left[r+cq^{2}\right]\delta_{\alpha\beta}+g_{0}\left(q_{\alpha}q_{\beta}/q^{2}\right)-h_{0}q_{\alpha}q_{\beta},$
$\displaystyle
C_{\alpha\beta\gamma\lambda}=\left(C_{v}-2C_{s}/3\right)\delta_{\alpha\beta}\delta_{\gamma\lambda}+C_{s}\left(\delta_{\alpha\gamma}\delta_{\beta\lambda}+\delta_{\alpha\lambda}\delta_{\beta\gamma}\right),$

and [3],

$\displaystyle
f_{\alpha\beta\gamma\lambda}=\left(f_{v}-2f_{s}/3\right)\delta_{\alpha\lambda}\delta_{\beta\gamma}+f_{s}\left(\delta_{\alpha\gamma}\delta_{\beta\lambda}+\delta_{\alpha\beta}\delta_{\gamma\lambda}\right),$

where $C_{v}=(C_{11}+2C_{12})/3$ and $C_{s}=(C_{11}-C_{12})/2$ are the bulk
and shear moduli and $f_{v}=(f_{11}+2f_{12})/3$ and $f_{s}=(f_{11}-f_{12})/2$
are bulk and shear components of the flexocoupling tensor.

Thus,

$\displaystyle
M_{\alpha\beta}({\bm{q}})=\left(C_{v}+C_{s}/3\right)q_{\alpha}q_{\beta}+C_{s}q^{2}\delta_{\alpha\beta},$

and,

$\displaystyle
V_{\alpha\beta}({\bm{q}})=\left(f_{v}+f_{s}/3\right)q_{\alpha}q_{\beta}+f_{s}q^{2}\delta_{\alpha\beta},$

In the normal coordinates of the bare modes without flexoelectricity, the
inverse propagator ${\bm{G}}_{0}^{-1}({\bm{q}})$ of Eq. (S1) is block-
diagonal,

$\displaystyle{\bm{G}}_{0}^{-1}({\bm{q}})=\begin{pmatrix}\omega_{\text{TO}}^{2}({\bm{q}})&V_{t}({\bm{q}})&0&0&0&0\\\
V_{t}({\bm{q}})&\omega_{\text{TA}}^{2}({\bm{q}})&0&0&0&0\\\
0&0&\omega_{\text{TO}}^{2}({\bm{q}})&V_{t}({\bm{q}})&0&0\\\
0&0&V_{t}({\bm{q}})&\omega_{\text{TA}}^{2}({\bm{q}})&0&0\\\
0&0&0&0&\omega_{\text{LO}}^{2}({\bm{q}})&V_{l}({\bm{q}})\\\
0&0&0&0&V_{l}({\bm{q}})&\omega_{\text{LA}}^{2}({\bm{q}})\end{pmatrix},$ (S2)

where,

$\displaystyle\omega_{\text{TO}}^{2}({\bm{q}})$
$\displaystyle=r+cq^{2}+hq^{4},$
$\displaystyle\omega_{\text{LO}}^{2}({\bm{q}})$
$\displaystyle=r+g_{0}+(c+h_{0})q^{2},$
$\displaystyle\omega_{\text{TA}}^{2}({\bm{q}})$ $\displaystyle=a_{t}q^{2},$
$\displaystyle\omega_{\text{LA}}^{2}({\bm{q}})$ $\displaystyle=a_{l}q^{2},$
$\displaystyle V_{t}({\bm{q}})$ $\displaystyle=v_{t}q^{2},$ $\displaystyle
V_{l}({\bm{q}})$ $\displaystyle=v_{l}q^{2},$

with $a_{t}\equiv\hbar^{2}C_{s}/\rho$,
$a_{l}\equiv\hbar^{2}(C_{v}+(4/3)C_{s})/\rho$, $v_{t}\equiv\hbar
f_{s}/\sqrt{\rho}$, and $v_{l}\equiv\hbar(f_{v}+(4/3)f_{s})/\sqrt{\rho}$,
where $\rho$ is the density in units of mass per unit volume. In writing Eq.
(S2), we have chosen the following ordering for the entries of
${\bm{G}}_{0}^{-1}({\bm{q}})$:
$\left\\{P_{\text{T1}}({\bm{q}}),u_{\text{T1}}({\bm{q}}),P_{\text{T2}}({\bm{q}}),u_{\text{T2}}({\bm{q}}),P_{\text{L}}({\bm{q}}),u_{\text{L}}({\bm{q}})\right\\}$,
where T and L refer to transverse and longitudinal components relative to
${\bm{q}}$, respectively. Note that $V_{\alpha\beta}({\bm{q}})$ does not mix
the transverse and longitudinal modes and that the transverse blocks are
identical. These are consequences of assuming an isotropic medium and do not
generally occur in crystals. Note we have added a $\mathcal{O}(q^{4})$ term in
$\omega_{\text{TO}}^{2}({\bm{q}})$ by hand. This is done to produce a well-
conditioned hybrid TA branch for $f_{s}>f_{cr}$ that would allow us to capture
the modulated wavevector in the hybrid modes (otherwise it would be negatively
unbounded, which is unphysical). We note such terms can result from higher-
order gradients in the effective Hamiltonian (e.g.,
$\left|{\nabla{\bm{P}}({\bm{x}})}\right|^{4}$). We do not explicitly write
them down in Eq. (1) for simplicity, as it is sometimes done in the continuum
theory of incommensurate transitions [4]. We also note that a similar term
could be added to $\omega_{\text{TA}}^{2}({\bm{q}})$, however, we do not do so
as it that would only add an extra fitting parameter with no qualitative
difference in the resulting hybrid modes.

Since $G_{0}^{-1}({\bm{q}})$ is a dynamical matrix, the bare frequencies are
given by its eigenvalues,

$\displaystyle\Omega_{\text{TO}}^{0}({\bm{q}})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\omega_{\text{TA}}^{2}({\bm{q}})+\omega_{\text{TO}}^{2}({\bm{q}})+\sqrt{4V_{t}^{2}({\bm{q}})+\left[\omega_{\text{TO}}^{2}({\bm{q}})-\omega_{\text{TA}}^{2}({\bm{q}})\right]^{2}}\right)^{1/2},$
(S3a) $\displaystyle\Omega_{\text{TA}}^{0}({\bm{q}})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\omega_{\text{TA}}^{2}({\bm{q}})+\omega_{\text{TO}}^{2}({\bm{q}})-\sqrt{4V_{t}^{2}({\bm{q}})+\left[\omega_{\text{TO}}^{2}({\bm{q}})-\omega_{\text{TA}}^{2}({\bm{q}})\right]^{2}}\right)^{1/2},$
(S3b) $\displaystyle\Omega_{\text{LO}}^{0}({\bm{q}})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\omega_{\text{LA}}^{2}({\bm{q}})+\omega_{\text{LO}}^{2}({\bm{q}})+\sqrt{4V_{l}^{2}({\bm{q}})+\left[\omega_{\text{LO}}^{2}({\bm{q}})-\omega_{\text{LA}}^{2}({\bm{q}})\right]^{2}}\right)^{1/2},$
(S3c) $\displaystyle\Omega_{\text{LA}}^{0}({\bm{q}})$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\omega_{\text{LA}}^{2}({\bm{q}})+\omega_{\text{LO}}^{2}({\bm{q}})-\sqrt{4V_{l}^{2}({\bm{q}})+\left[\omega_{\text{LO}}^{2}({\bm{q}})-\omega_{\text{LA}}^{2}({\bm{q}})\right]^{2}}\right)^{1/2},$
(S3d)

where the transverse and longitudinal branches are two- and one-fold
degenerate, respectively. The TA, TO, LA, and LO labels refer to optic- and
acoustic-like mode behavior in the long-wavelength limit.

We now discuss the behavior of the hybrid modes. For
$0<v_{t}<v_{cr}\equiv\sqrt{c\,a_{t}}$, the TA an TO phonons repel at finite
wavevectors and soften as $T\to T_{0}$. At $T_{0}$,
$\Omega_{\text{TO}}^{0}({\bm{q}}=0)$ condenses, which leads to a FE
transition. For $v_{t}\geq v_{cr}$, the level repulsion generates a minimum at
a finite wavevector $q_{mod}$ in the acoustic branch as temperature decreases.
At a temperature
$T_{mod}=T_{0}+(v_{t}^{2}-v_{cr}^{2})^{2}/(4hr_{0}a_{t}^{2})$,
$q_{mod}=\sqrt{(v_{t}^{2}-v_{cr}^{2})/(2ha_{t})}$ and
$\Omega_{\text{TA}}^{0}({\bm{q}}_{mod})$ condenses, while the zone-center TO
phonon remains gaped with energy
$\tau^{0}_{mod}\equiv\Omega_{\text{TO}}^{0}({\bm{q}}=0)=\sqrt{r_{0}(T_{mod}-T_{0})}=\left(v_{cr}^{2}-v_{t}^{2}\right)/(2a_{t}\sqrt{h})$.
We note that $v_{t}=v_{cr}$ is a triple point between the paraelectric, FE and
modulated phases. N.B.: since $f_{s}/f_{cr}=v_{t}/v_{cr}$, where
$f_{cr}=\sqrt{c\,C_{s}}$, the above results can be rewritten as follows:
$T_{mod}=T_{0}+(f_{s}^{2}-f_{cr}^{2})^{2}/(4hr_{0}C_{s}^{2})$,
$q_{mod}=\sqrt{(f_{s}^{2}-f_{cr}^{2})/(2hC_{s})}$,
$\Omega_{\text{TO}}^{0}({\bm{q}}=0)=\left(f_{cr}^{2}-f_{s}^{2}\right)/(2C_{s}\sqrt{h})$.
These resulting phase diagram is shown in Fig. 1 (a).

For the longitudinal modes, the relevant parameter is $v_{l}$: for $0\leq
v_{l}<v_{l,cr}\equiv\sqrt{\left(c+h_{0}\right)\,a_{l}}$, the zone center LO
phonon does not condense due to the depolarizing field $g_{0}$. For $v_{l}\geq
v_{l,cr}$, $\Omega_{\text{LA}}^{0}({\bm{q}}_{l,mod})$ condenses at a
temperature
$T_{l,mod}=T_{0}-(g_{0}/r_{0})+\left(v_{l}^{2}-{v_{l,cr}}^{2}\right)^{2}/(4ha_{l}^{2}r_{0})$,
where $q_{l,mod}=\sqrt{\left(v_{l}^{2}-v_{l,cr}^{2}\right)/(2a_{l}h)}$. The
zone-center LO phonon remains gaped with (squared) energy
$\Omega_{\text{LO}}^{0}({\bm{q}}=0)=\left(v_{l,cr}^{2}-{v_{l}}^{2}\right)/(2a_{l}\sqrt{h})$.
We note that $T_{mod}$ and $T_{l,mod}$ depend for the most part on independent
model parameters, therefore one cannot conclude in general which modulated
transition precedes the other. Nonetheless, $g_{0}$ is typically very large in
pervoskite FEs, thus we expect $T_{mod}>T_{l,mod}$, which will be assumed
hereafter.

##  Supplementary Note 2. Phase stability

We now discuss the stability of the PB. We will show that in the presence of
purely thermal fluctuations the FE PB is stable, while the modulated solid and
triple point in the phase diagram are unstable.

At the Gaussian level and in the classical limit, the local correlation
functions $\left<\left|{{\bm{P}}^{2}}\right|\right>$ of the disordered phase
of our isotropic model are given in terms of its transverse
($\left<P_{\perp}^{2}({\bm{q}})\right>_{0}\equiv\left<P_{\text{
T1}}^{2}({\bm{q}})\right>_{0}=\left<P_{\text{T2}}^{2}({\bm{q}})\right>_{0}$)
and longitudinal fluctuations
($\left<P_{\parallel}^{2}({\bm{q}})\right>_{0}\equiv\left<P_{\text{L}}^{2}({\bm{q}})\right>_{0}$)
as follows,

$\displaystyle\left<\left|{{\bm{P}}^{2}}\right|\right>_{0}=\int\frac{d^{3}q}{(2\pi)^{3}}\left[2\left<P_{\perp}^{2}({\bm{q}})\right>_{0}+\left<P_{\parallel}^{2}({\bm{q}})\right>_{0}\right],$
(S4)

where,

$\displaystyle\left<P_{\perp}^{2}({\bm{q}})\right>_{0}=\frac{k_{B}T}{r_{0}(T-T_{0})+(1-v_{t}^{2}/v_{cr}^{2})cq^{2}+hq^{4}},$

and,

$\displaystyle\left<P_{\parallel}^{2}({\bm{q}})\right>_{0}=\frac{k_{B}T}{r_{0}(T-T_{0})+g_{0}+(1-v_{l}^{2}/v_{l,cr}^{2})(c+h_{0})q^{2}+hq^{4}}.$

At the onset of the FE transition ($T\to T_{0}$ for $0\leq v_{t}<v_{cr}$),
$\left<P_{\perp}^{2}({\bm{q}})\right>_{0}\to\infty$ as $q\to 0$, while
$\left<P_{\parallel}^{2}({\bm{q}})\right>_{0}$ remains finite due to the
depolarizing field term $g_{0}$ and therefore we will ignored it. By
integrating $\left<P_{\perp}^{2}({\bm{q}})\right>_{0}$ over all of
${\bm{q}}$-space, we obtain the following result for $T\to T_{0}$,

$\displaystyle(k_{B}T)^{-1}\int
d^{3}q\left<P_{\perp}^{2}({\bm{q}})\right>_{0}$
$\displaystyle\propto\frac{\pi}{2\sqrt{\left(1-v_{t}^{2}/v_{cr}^{2}\right)hc}}\left[1-\sqrt{\frac{h}{c}\frac{\left(T-T_{0}\right)}{\left(1-v_{t}^{2}/v_{cr}^{2}\right)}}\right].$

The thermal fluctuations are finite at $T_{0}$, thus the long-range FE order
is stable in the purely classical limit.

We now look at the fluctuations near above the modulated transition. For $T\to
T_{mod}$ and $v_{t}\geq v_{cr}$,
$\left<P_{\perp}^{2}({\bm{q}})\right>_{0}\to\infty$ as $q\to q_{mod}$, whereas
$\left<P_{\parallel}^{2}({\bm{q}})\right>_{0}$ remains finite since
$T_{mod}>T_{l,mod}$ by assumption and thus we ignore it. By integrating the
transverse fluctuations over momentum, we obtain the following result for
$T\to T_{mod}$,

$\displaystyle(k_{B}T)^{-1}\int
d^{3}q\left<P_{\perp}^{2}({\bm{q}})\right>_{0}\propto\left(\frac{T_{mod}-T_{0}}{T-T_{mod}}\right)^{1/2}.$
(S5)

Equations (S5) diverges as $T\to T_{mod}$, thus the modulated solid is
unstable.

At the triple point ($v_{t}=v_{cr}$), $T_{mod}=T_{0}$ and $q_{mod}=0$, thus,

$\displaystyle(k_{B}T)^{-1}\int
d^{3}q\left<P_{\perp}^{2}({\bm{q}})\right>_{0}\propto\left(\frac{1}{T-T_{0}}\right)^{1/4}.$
(S6)

Equation (S6) diverges as $T\to T_{0}$, therefore the coexistence is unstable.

##  Supplementary Note 3. Self-consistent phonon approximation

We follow a standard procedure [5] to solve the statistical mechanical problem
posed by our model Hamiltonian in disordered phase. At the Gaussian level, the
quartic anharmonicities in Eq. (1) are approximated by,

$\displaystyle\int
d^{3}x\sum_{\alpha,\beta}uP^{2}_{\alpha}({\bm{x}})P^{2}_{\beta}({\bm{x}})\simeq\int
d^{3}x\sum_{\alpha,\beta}\left[2u\left<\left|{{\bm{P}}({\bm{x}})}\right|^{2}\right>\delta_{\alpha\beta}+4u\left<P_{\alpha}({\bm{x}})P_{\beta}({\bm{x}})\right>\right]P_{\alpha}({\bm{x}})P_{\beta}({\bm{x}}),$

where $\left<P_{\alpha}({\bm{x}})P_{\beta}({\bm{x}})\right>$ are correlation
functions of polarization and
$\left<\left|{{\bm{P}}({\bm{x}})}\right|^{2}\right>=\sum_{\gamma=1}^{3}\left<P^{2}_{\gamma}({\bm{x}})\right>$.
We thus have an effective Hamiltonian identical to Eq. (1), but with
$U_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})$ now given as follows,

$\displaystyle
U_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})=\left[\left(r+4u\left<\left|{{\bm{P}}({\bm{x}})}\right|^{2}\right>\right)\delta_{\alpha\beta}+8u\left<P_{\alpha}({\bm{x}})P_{\beta}({\bm{x}})\right>\right]\delta({\bm{x}}-{\bm{x}}^{\prime})+F_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime}).$

By symmetry,
$\left<P_{1}({\bm{x}})P_{2}({\bm{x}})\right>=\left<P_{1}({\bm{x}})P_{3}({\bm{x}})\right>=\left<P_{2}({\bm{x}})P_{3}({\bm{x}})\right>=0$
and
$\left<P_{1}^{2}({\bm{x}})\right>=\left<P_{2}^{2}({\bm{x}})\right>=\left<P_{3}^{2}({\bm{x}})\right>=(1/3)\left<\left|{\bm{P}}\right|^{2}\right>$,
thus,

$\displaystyle
U_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})=\left[r+(20u/3)\left<\left|{\bm{P}}\right|^{2}\right>\right]\delta_{\alpha\beta}\delta({\bm{x}}-{\bm{x}}^{\prime})+F_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime}).$
(S7)

The dressed propagator ${\bm{G}}^{-1}({\bm{q}})$ and frequencies
$\Omega_{\lambda}({\bm{q}}),~{}(\lambda=\text{TO, TA, LO, LA})$ are therefore
identical to their corresponding bare forms given in Eqs. (S2) and (S3) with
$r$ replaced by,

$\displaystyle\tau\equiv\Omega_{\text{TO}}^{2}({\bm{q}}=0)=r+(20u/3)\left<\left|{\bm{P}}\right|^{2}\right>,$
(S8)

where $\Omega_{\text{TO}}(0)$ is the zone-center TO mode frequency. The
correlation functions are calculated from the fluctuation-dissipation theorem
[6],

$\displaystyle\left<P_{\alpha}({\bm{x}})P_{\beta}({\bm{x}}^{\prime})\right>=\int\frac{d^{3}q}{(2\pi)^{3}}e^{i{\bm{q}}\cdot\left({\bm{x}}-{\bm{x}}^{\prime}\right)}\sum_{\lambda\lambda^{\prime}}b_{\alpha\lambda}({\bm{q}})\psi_{\lambda\lambda^{\prime}}({\bm{q}})b_{\lambda^{\prime}\beta}^{\dagger}({\bm{q}}),$
(S9)

where,

$\displaystyle\psi_{\lambda\lambda^{\prime}}({\bm{q}})=\frac{1}{2\Omega_{\lambda}({\bm{q}})}\coth\left(\frac{\Omega_{\lambda}({\bm{q}})}{2k_{B}T}\right)\,\delta_{\lambda\lambda^{\prime}},$

and $b_{\alpha\lambda}({\bm{q}})$ is a unitary transformation that
diagonalizes ${\bm{G}}^{-1}({\bm{q}})$. Eq. (S9) provides the self-consistency
equation for $\Omega_{\lambda}({\bm{q}})$.

###  .1 Dependence of $\tau$ on $u$ for $f_{s}>f_{cr}$

Figs. S1 (a) and (b) above show the temperature dependence of
$\tau\equiv\Omega_{\text{TO}}^{2}(0)$ and its first derivative $d\tau/dT$ for
several values of $u$. We first point out the dramatic change in the behavior
of $\tau$ between $u=0$ and $u>0$ is expected as infinitesimally small
fluctuations prevent long-range modulated order. We then note that for finite
$u$, $d\tau/dT$ seems to become multi-valued at $T_{mod}$ as $u\to 0$ (but
remains finite). However, a closer inspection reveals that it is well-defined.

Figure S1: (a) Temperature dependence of the zone-center TO mode $\tau$ and
(b) its first derivative for several values of the dimensionaless quartic
coefficient $\tilde{u}\equiv(20/3\pi)(uq_{mod}^{3}/\sqrt{\tau_{mod}})$. All
other parameters are taken from Table I of the main text for STO.

For large flexocouplings ($f_{s}>f_{cr}$), we rewrite Eq. (S8) as follows,

$\displaystyle\tau-\tau_{mod}^{0}=r_{0}\left(T-T_{mod}\right)+\left(20u/3\right)\left<\left|{{\bm{P}}(0)}\right|^{2}\right>,$
(S10)

where $\left<\left|{{\bm{P}}(0)}\right|^{2}\right>$ are the local correlation
functions of polarization, and $\tau_{mod}^{0}\equiv
r_{0}\left(T_{mod}-T_{0}\right)$ is the harmonic ($u=0$) zone-center TO
frequency at $T_{mod}$.

In the classical limit and $\tau-\tau_{mod}^{0}\to 0$,

$\displaystyle\left<\left|{{\bm{P}}(0)}\right|^{2}\right>=\int\frac{d^{3}q}{(2\pi)^{3}}\frac{k_{B}T}{\left(\tau-\tau_{mod}^{0}\right)+h\left(q^{2}-q_{mod}^{2}\right)^{2}}\propto\frac{k_{B}T}{\left(\tau-\tau_{mod}^{0}\right)^{1/2}}.$

We now substitute this into Eq. (S10),

$\displaystyle\tau-\tau_{mod}^{0}\simeq
r_{0}\left(T-T_{mod}\right)+\frac{u^{\prime}T}{\left(\tau-\tau_{mod}^{0}\right)^{1/2}},$

with $u^{\prime}\propto u$. By taking the derivative with respect to $T$ on
both sides of this equation and solve for $d\tau/dT$, we obtain the following
result,

$\displaystyle\frac{d\tau}{dT}=\frac{r_{0}+u^{\prime}/\left(\tau-\tau_{mod}^{0}\right)^{1/2}}{1+u^{\prime}T/2\left(\tau-\tau_{mod}^{0}\right)^{3/2}}.$

For the harmonic theory (i.e. $u^{\prime}=u=0$), we get $d\tau/dT=r_{0}$ at
all temperatures, as expected. For $u^{\prime}>0$ and $T=T_{mod}$,
$\tau-\tau_{mod}^{0}\simeq\left(u^{\prime}T_{mod}\right)^{2/3}$ and,

$\displaystyle\left.\frac{d\tau}{dT}\right|_{T=T_{mod}}=\frac{2}{3}\left(r_{0}+\frac{{u^{\prime}}^{2/3}}{T_{mod}^{1/3}}\right).$

We thus find $d\tau/dT$ is single-valued at $T_{mod}$ for $u^{\prime}>0$ and
thus there is no kink.

###  .2 Correlation function of polarization in the lamellar-like liquid phase

The purpose of this subsection is the derive the correlation function of
polarization given in Eq. (2) of the main text. We consider the correlation
function of polarization given as follows,

$\displaystyle\left<{\bm{P}}({\bm{x}})\cdot{\bm{P}}({\bm{0}})\right>$
$\displaystyle=\sum_{\alpha}\left<P_{\alpha}({\bm{x}})P_{\alpha}({\bm{0}})\right>=\int\frac{d^{3}q}{(2\pi)^{3}}e^{i{\bm{q}}\cdot{\bm{x}}}\left[2\left<P_{\perp}^{2}({\bm{q}})\right>+\left<P_{\parallel}^{2}({\bm{q}})\right>\right],$
(S11)

where $\left<P_{\perp}^{2}({\bm{q}})\right>$ and
$\left<P_{\parallel}^{2}({\bm{q}})\right>$ are transverse and longitudinal
polarization correlations. In the classical limit,

$\displaystyle\left<P_{\perp}^{2}({\bm{q}})\right>$
$\displaystyle=\frac{k_{B}T}{\tau+(1-f_{s}^{2}/f_{cr}^{2})cq^{2}+hq^{4}},$
(S12)

and,

$\displaystyle\left<P_{\parallel}^{2}({\bm{q}})\right>=\frac{k_{B}T}{\tau+g_{0}+(1-f_{l}^{2}/f_{l,cr}^{2})(c+h_{0})q^{2}+hq^{4}}.$

where $\tau$ is given by Eq. (S8) and $q_{mod}$ and $\tau_{mod}^{0}$ are given
in SN1. Since the longitudinal fluctuations are much weaker than the
transverse ones due to the large depolarizing term $g_{0}$, we will ignore
them.

For $f_{s}>f_{cr}$, we rewrite Eq. (S12) as follows,

$\displaystyle\left<P_{\perp}^{2}({\bm{q}})\right>=\frac{q_{mod}^{4}}{\tau_{mod}^{0}}\frac{k_{B}T}{\xi^{-4}+(q^{2}-q_{mod}^{2})^{2}},$

where $\xi$ is the correlation length,

$\displaystyle
q_{mod}\,\xi\equiv\left(\frac{\tau^{0}_{mod}}{\tau-\tau^{0}_{mod}}\right)^{1/4}.$

Since the fluctuations are isotropic, the angular integrals in Eq. (S11) can
be readily calculated,

$\displaystyle\left<{\bm{P}}({\bm{x}})\cdot{\bm{P}}({\bm{0}})\right>=\frac{k_{B}T}{\pi^{2}}\frac{q_{mod}^{4}}{\tau_{mod}^{0}}\int_{0}^{\infty}dq\,\frac{q^{2}}{\xi^{-4}+(q^{2}-q_{mod}^{2})^{2}}\frac{\sin
qx}{qx},$

where $x=\left|{\bm{x}}\right|$. After evaluating this integral by countour
integration and the residue theorem, we obtain the following result,

$\displaystyle\left<{\bm{P}}({\bm{x}})\cdot{\bm{P}}({\bm{0}})\right>$
$\displaystyle=\frac{k_{B}T}{2\pi\xi^{-2}}\frac{q_{mod}^{5}}{\tau_{mod}^{0}}e^{-q_{mod}\,r\sqrt{\frac{\sqrt{1+(q_{mod}\xi)^{-4}}-1}{2}}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\sin\left(q_{mod}\,r\sqrt{\frac{\sqrt{1+(q_{mod}\xi)^{-4}}+1}{2}}\right)/(q_{mod}x).$

For large correlation lengths ($q_{mod}\,\xi\gg 1$), we obtain,

$\displaystyle\left<{\bm{P}}({\bm{x}})\cdot{\bm{P}}({\bm{0}})\right>\simeq\frac{k_{B}T}{2\pi\xi^{-2}}\frac{q_{mod}^{5}}{\tau_{mod}^{0}}\,e^{-x/(2\,q_{mod}\,\xi^{2})}\,\sin\left(q_{mod}\,x\right)/\left(q_{mod}\,x\right),$

which is the desired result. We note the calculation of the correlation
function in the Gaussian approximation is identical to the one presented here
except that one must replace $\tau\to r_{0}(T-T_{0})$ in Eq. (S12) and define
the correlation length as
$q_{mod}\,\xi=\left(T_{mod}/(T-T_{mod})\right)^{1/4}$.

##  I  Supplementary Note 4. Fits to experiments in STO and KTO

In this note, we describe our fits to the phonon spectra of STO and KTO. Our
model has seven parameters $r_{0},T_{0},u,c,a_{t},v_{t}$ and $h$, which are
given in Table I of the main text. For both materials, we have attempted fits
for $f_{s}\geq f_{cr}$ and $0\leq f_{s}<f_{cr}$.

###  I.1 Large couplings: $f_{s}\geq f_{cr}$ (i.e. $v_{t}\geq v_{cr}$)

STO: We take $r_{0},T_{0}$ and $u$ from Landau theory [7]. Values for
$c,a_{t}$, and $v_{t}$ are available in Ref. [8], however, they do not provide
a good fit to the recent neutron scattering experiments of Ref. [9]. We thus
provide a new parametrization.

From the Supplementary Note 2, we have $a_{t}=\hbar^{2}C_{s}/\rho$. We take
$C_{s}=1.24\times 10^{11}$ N m-2 from the observed low-temperature shear
modulus and $\rho=5.104\times 10^{3}\,$kg m-3 [10]. To obtain $c$, we expand
Eq. (S3a) to leading order in $q$, i.e.,
$\Omega_{\text{TO}}^{2}({\bm{q}})\simeq\tau+cq^{2}$ and fit it to the observed
TO mode dispersion at 10 K [9]. $v_{t}$ and $h$ were obtained by fitting
$\tau_{mod}^{1/2}$ and $q_{mod}$ to the zone center TO energy at low
temperatures ($\simeq 1.70\,\text{meV}$) [11, 12] and to the wavevector
minimum in the TA branch ($\simeq 0.024\,$rlu) [9], respectively.

KTO: To determine $r_{0}$, $T_{0}$, and $u$, we first take the Landau limit of
our model with a homogeneous polarization
(${\bm{P}}({\bm{x}})=\left<{\bm{P}}\right>_{\text{MFT}}$) and
$f_{\alpha\beta\gamma\lambda}=F_{\alpha\beta}({\bm{x}},{\bm{x}}^{\prime})=0$.
We then fit $r_{0}$, $T_{0}$, and $u$ to the observed temperature and applied
electric field dependence of the dielectric constant [13, 14]. The values for
$c$ and $a_{t}$ are known from neutron scattering experiments Ref. [15].
$v_{t}$ and $h$ are obtained by fitting $\tau_{mod}^{1/2}$ and $q_{mod}$ to
the observed zone center TO energy at low temperatures ($\simeq
2.5\,\text{meV}$) [16, 12, 15] and the wavevector minimum in the TA branch
$(\simeq 0.1\,$rlu) [17], respectively.

###  I.2 Small couplings: $0\leq f_{s}<f_{cr}$ (i.e. $0\leq v_{t}<v_{cr}$)

STO and KTO: We assume there is no long-range polar order in STO and KTO down
to absolute zero, i.e., the energy gap of the zone-center TO phonon is finite
at $0\,$K (see Eq. (S8)),

$\displaystyle\tau=\Omega_{\text{TO}}^{2}({\bm{q}}=0)=-r_{0}T_{0}+\left(20u/3\right)\left<\left|{\bm{P}}\right|^{2}\right>>0,$
(S13)

where $\left<\left|{\bm{P}}\right|^{2}\right>$ are the fluctuations of
polarization calculated within the SCPA (see Supplementary Note 3). By taking
$r_{0},T_{0},u,c,a_{t}$ and $h$ from the previous section and setting
$\tau=1.7\,\text{meV}\,(2.5\,\text{meV})$ for STO (KTO), we find that the
inequality (S13) cannot be satisfied for $0\leq v_{t}<v_{cr}$ in either
material. We find the same result if we attempt to fit $r_{0},T_{0}$, and $u$
while keeping the polarization physical, which we did by constraining
$r_{0},T_{0}$, and $u$ to produce the Landau value at $0\,$K, i.e.,
$\left<{\bm{P}}\right>_{\text{MFT}}\,=\left(\sqrt{2\epsilon_{0}}A\times
10^{2}\right)\sqrt{r_{0}T_{0}/(2u)}=8.6\,$$\mu$C /cm${}^{2}\,(2.3\,$$\mu$C
/cm${}^{2})$ for STO (KTO), where $\epsilon_{0}$ is the vacuum permittivity in
units of C${}^{2}\,$N${}^{-1}\,$m-2, and $A=194.4\,$meV ($158.7\,$meV) is a
proportionality factor between the experimental zone-center TO frequency
($\Omega_{\text{TO}}^{\text{exp}}$) and the dielectric constant
$(\epsilon^{\text{exp}})$:
$\Omega_{\text{TO}}^{\text{exp}}=A/\sqrt{\epsilon^{\text{exp}}}$ [11].
$\sqrt{2\epsilon_{0}}A\times 10^{2}$ is a conversion factor from our model
units to $\mu$C /cm2.

By relaxing the polarization constraint and fitting $r_{0},T_{0}$ and $u$ to
the observed zone-center TO mode, we find a fair agreement with experiments,
as shown in Fig. S2. However, the resulting parametrization yields an
unphysical $\left<{\bm{P}}\right>_{\text{MFT}}\,$ (see Table S1).

Figure S2: (a),(b) Calculated temperature dependence of the zone center hybrid TO mode compared to neutron [16, 11, 15] and hyper-Raman [12] scattering experiments . (c),(d) Calculated dispersion of the hybrid TA and TO phonons (solid lines) compared to neutron scattering experiments (circles and diamonds) [9, 15]. Dashed line in (d) is the bare linear TA dispersion for $f_{s}=0$ with a sound velocity of 4800 m/s (as shown in Ref. [9] for comparison). Model parameters are given in Table S1. Table S1: Model parameters for small flexocoplings and no polarization constrain. For STO, we took $c$ from Table I, $v_{t}$ from Ref. [8], and obtained $r_{0},T_{0}$ and $u$ by fitting Eq. (S8) to $\tau^{1/2}=1.7,6.7$ and $11.1\,$meV at $T=0,127$, and $292\,$K, respectively, with $h=0$. For KTO, we took $c,a_{t}$ and $v_{t}$ from Ref. [15], and obtained $r_{0},T_{0}$ and $u$ by fitting Eq. (S8) to $\tau^{1/2}=2.5,7.1$ and $10.2\,$meV at $T=0,131$, and $297\,$K, respectively, with $h=0$. We could not fit the data for either material for $h>0$. |  KTO |  STO
---|---|---
$r_{0}\,$[meV2 K-1] | $0.16$ | $0.20$
$T_{0}\,$[K] | $100$ | $473$
$u\,$ [$10^{3}$ meV3 rlu-3] | $0.12$ | $1.29$
$c\,$[$10^{3}$ meV2 rlu-2] |  $4.828$111From Ref. [15]. | $13.5$
$a_{t}$ [$10^{3}\,$ meV2 rlu-2] |  $1.553$111From Ref. [15]. |  $2.739$222From Ref. [10].
$v_{t}$ [$10^{3}\,$ meV2 rlu-2] |  $2.451$111From Ref. [15]. |  $2.802$333From Ref. [8].
$h\,$ [meV4 rlu-4] | $0.0$ | $0.0$
$v_{t}/v_{cr}=f_{s}/f_{cr}$ | $0.90$ | $0.50$
$\left<{\bm{P}}\right>_{\text{MFT}}\,$ [$\mu$C/cm2] | $0.02$ | $0.02$

## References

  * Aharony and Fisher [1973] A. Aharony and M. E. Fisher, Critical behavior of magnets with dipolar interactions. i. renormalization group near four dimensions, Phys. Rev. B 8, 3323 (1973).
  * Oliver and Agelet de Saracibar [2017] X. Oliver and C. Agelet de Saracibar,  _Continuum Mechanics for Engineers. Theory and Problems_ (2nd edition, 2017).
  * Le Quang and He [2011] H. Le Quang and Q.-C. He, The number and types of all possible rotational symmetries for flexoelectric tensors, Proc. Royal Soc. A 467, 2369 (2011).
  * Strukov and Levanyuk [1998] B. A. Strukov and A. P. Levanyuk,  _Ferroelectric Phenomena in Crystals: Physical Foundations_ (Springer, 1998).
  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Pytte [1972] E. Pytte, Theory of perovskite ferroelectrics, Phys. Rev. B 5, 3758 (1972).
  * Rabe  _et al._ [2007] K. M. Rabe, C. H. Ahn, and J.-M. Triscone, eds.,  _Physics of Ferroelectrics A Modern Perspective_ (Springer-Verlag, 2007).
  * Vaks [2007] V. G. Vaks, ed.,  _Introduction to the Microscopic Theory of Ferroelectrics_ (Nakuna, 2007).
  * Fauqué  _et al._ [2022] B. Fauqué, P. Bourges, A. Subedi, K. Behnia, B. Baptiste, B. Roessli, T. Fennell, S. Raymond, and P. Steffens, Mesoscopic fluctuating domains in strontium titanate, Phys. Rev. B 106, L140301 (2022).
  * Carpenter [2007] M. A. Carpenter, Elastic anomalies accompanying phase transitions in (Ca,Sr)TiO3 perovskites: Part I. Landau theory and a calibration for SrTiO3, Am. Mineral 92, 309 (2007).
  * Yamada and Shirane [1969] Y. Yamada and G. Shirane, Neutron Scattering and Nature of the Soft Optical Phonon in SrTiO3, J. Phys. Soc. Jpn. 26, 396 (1969).
  * Vogt [1995] H. Vogt, Refined treatment of the model of linearly coupled anharmonic oscillators and its application to the temperature dependence of the zone-center soft-mode frequencies of ${\mathrm{KTaO}}_{3}$ and ${\mathrm{SrTiO}}_{3}$, Phys. Rev. B 51, 8046 (1995).
  * Fujishita  _et al._ [2016] H. Fujishita, S. Kitazawa, M. Saito, R. Ishisaka, H. Okamoto, and T. Yamaguchi, Quantum paraelectric states in SrTiO3 and KTaO3: Barrett model, Vendik model, and quantum criticality, J. Phys. Soc. Jpn. 85, 074703 (2016).
  * Skoromets  _et al._ [2016] V. Skoromets, C. Kadlec, H. Němec, D. Fattakhova-Rohlfing, and P. Kužel, Tunable dielectric properties of KTaO3 single crystals in the terahertz range, J. Phys. D: Appl. Phys. 49, 065306 (2016).
  * Farhi  _et al._ [2000] E. Farhi, A. K. Tagantsev, R. Currat, B. Hehlen, E. Courtens, and L. A. Boatner, Low energy phonon spectrum and its parameterization in pure KTaO3 below 80 K, Eur. Phys. J. B. 15, 615 (2000).
  * Shirane  _et al._ [1967] G. Shirane, R. Nathans, and V. J. Minkiewicz, Temperature dependence of the soft ferroelectric mode in KTaO3, Phys. Rev. 157, 396 (1967).
  * Axe  _et al._ [1970] J. D. Axe, J. Harada, and G. Shirane, Anomalous acoustic dispersion in centrosymmetric crystals with soft optic phonons, Phys. Rev. B 1, 1227 (1970).