Quantum Electrodynamic Control of Matter: Cavity-Enhanced Ferroelectric Phase Transition

The interaction between light and matter can be utilized to qualitatively alter physical properties of materials. Recent theoretical and experimental studies have explored this possibility of controlling matter by light based on driving many-body systems via strong classical electromagnetic radiation, leading to a time-dependent Hamiltonian for electronic or lattice degrees of freedom. To avoid inevitable heating effects in driven interacting systems, pump-probe setups with ultrashort laser pulses have so far been used to study transient light-induced modifications of material properties. Here, we pursue yet another direction of controlling quantum matter by modifying quantum fluctuations of its electromagnetic environment. In contrast to earlier proposals focusing on light-enhanced electron-electron interactions, we consider a dipolar quantum many-body system embedded in a cavity composed of metal mirrors, and formulate a theoretical framework to manipulate its equilibrium properties on the basis of quantum light-matter interaction. Owing to the strong coupling of polarizable material excitations to the cavity electric field, the cavity confinement leads to the hybridization of different types of the fundamental excitations, including dipolar phonons, cavity photons, and plasmons in metal mirrors. We show that this hybridization qualitatively alters the nature of the collective excitations and can be used to control energy-level structures in a wide range of systems. Most notably, in quantum paraelectric materials, we show that the cavity-induced softening of infrared optical phonons enhances the ferroelectric phase. Our findings suggest an intriguing possibility of inducing a superradiant-type quantum phase transition via the light-matter coupling without external pumping. We also discuss possible applications to collective excitations in molecular materials and excitonic devices.


I. INTRODUCTION
One of the central goals in both quantum optics and condensed matter physics is to understand macroscopic phenomena emerging from the fundamental quantum degrees of freedom. Historically, quantum optics strived to study lightmatter interactions mainly in few-body regimes [1]. On another front, condensed matter physics investigated collective phenomena of many-body systems, and largely focused on optimizing structural and electronic phases of solids. In particular, a number of techniques, including applied strain, electronic doping, and isotope substitution, have been used to enhance a macroscopic order such as ferroelectricity or superconductivity. As an interdisciplinary frontier at the intersection of optics and condensed matter physics, there have recently been remarkable developments in using classical electromagnetic radiation to control transient states of matter [2]. The aim of this paper is to explore yet another route towards controlling the phase of matter by quantum light, namely, by modifying quantum fluctuations of the electromagnetic field in equilibrium -in the absence of an external drive.
Studies of light-induced modifications in material properties date back to experiments by Dayem and Wyatt [3,4], who observed that coherent microwave radiation leads to an increase of the critical current in superconductors. This phenomenon has been understood from the Eliashberg electronphoton theory [5][6][7], in which photons with frequencies below the quasiparticle threshold are shown to modify the distribution of quasiparticles in a way that the superconducting order parameter is enhanced. Other seminal works studying control of matter by light include Floquet engineering of electronic band structures [8][9][10][11][12], and photo-induced pumping by ultrashort laser pulses for creating metastable nonequilibrium states [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. Nonlinear optical phenomena can also be understood from the perspective of light-induced changes in optical properties of matter [29][30][31]. For instance, the lightinduced transparency can be considered as strong modification of the refractive index and absorption coefficient by a control optical beam [32], and the optical phase conjugation phenomena arise from the resonant parametric scattering of signal photons and matter excitations [33,34]. The essential element of all these developments is an external pump that provides classical light field acting on quantum matter.
An alternative approach to controlling matter by light is to utilize quantum fluctuations of the vacuum electromagnetic field. A prototypical model to understand such quantum lightmatter interaction is the Dicke model, which describes an ensemble of two-level atoms coupled to a single electromagnetic mode of a cavity [35][36][37][38][39][40]. For a sufficiently strong electric dipole, Dicke has shown that one should find instability into the superradiant phase characterized by spontaneous atomic polarization and macroscopic photon occupation [35]. This transition originates from the competition between the decrease of the total energy due to the light-matter coupling and the energy cost of adding photons to the cavity and admixing the excited state of atoms. The Dicke model is exactly solvable via the Bethe-ansatz equations originally discussed in the central spin model [41][42][43], and the existence of the superradiant transition has been rigorously proven in the thermodynamic limit [37,38].
Subsequently to Dicke's original work, it has been pointed out that the superradiant transition requires such strong interaction that one needs to include theÂ 2 term in the Hamiltonian, whereÂ denotes the vector potential of electromagnetic fields. In turn, the latter term has been shown to prevent the transition [44,45]. Even though ideas have been suggested to circumvent this no-go theorem by using microwave resonators [46,47] or including additional inter-particle interactions [48,49], the problem is still a subject of debate [50]. An alternative proposal using multilevel systems with external pump fields [51,52] has been experimentally realized [53,54]. Analogues of the Dicke transition have also been explored in ultracold atoms, where matter excitations correspond to different momentum states of atoms [55,56] or hyperfine states in spin-orbit-coupled Bose gases [57]. In all these realizations, external driving is essential to obtain the superradiant-type phases.
The present work in this context suggests a promising route towards realizing a genuine equilibrium superradiant transition under no external drive, which has so far remained elusive. Our consideration is motivated by an observation that transition into a superradiant state should be easier to attain in a system naturally close to a phase with spontaneous electric polarization. We investigate cavity-induced changes in a quantum paraelectric material at the verge of the ferroelectric order (see Fig. 1(a-c)); such materials can remain paraelectrics even at zero temperature while the quantum phase transition into the ferroelectric phase can be induced by, e.g., isotope substitution. We thus discuss a regime of a thermodynamically large system, in which multiple continuum modes of both electromagnetic fields and matter excitations must be included. This makes a sharp contrast to the previous studies on the superradiant transitions, which have solely focused on a simplified case of many-body systems coupled to only a few modes of a cavity. We formulate a simple yet general theoretical framework to analyze hybridization of continuum matter excitations and light modes confined in metallic cladding layers acting as cavity mirrors (see Fig. 2). Resulting polariton excitations consist of infrared active phonons in the quantum paraelectric, electromagnetic fields in the cavity, and plasmons in the metallic electrodes. We show that the spectrum of The horizontal axis represents the strength of a tuning parameter such as pressure or substitution, which governs the transition between ferroelectric and paraelectric phases. The vertical axis is the temperature while the red dashed line in the paraelectric phase indicates the fact that the energy gap of an optical phonon decreases and goes down to zero at the paraelectric-to-ferroelectric transition point. The main aim of this paper is to show that this structural quantum phase transition can be induced by softening the phonon mode via the light-matter coupling in the cavity with no external pump. (b,c) Ferroelectric and paraelectric phases in the case of ionic crystals. Green (blue) circles indicate the ions with the positive (negative) charges. The ferroelectric phase in (b) is characterized by a nonzero uniform displacement causing an electric polarization (red arrows) while the mean-phonon displacement is absent in the paraelectric phase in (c). The variables φ and Ω represent the phonon field and the bare phonon frequency, respectively (see Sec. III). such cavity polaritons is qualitatively different from the bulk spectrum, and can exhibit significant mode-softening, indicating enhancement of the ferroelectric instability.
Building upon the analysis of collective hybridized modes, we investigate the quantum phase transition into the ferroelectric phase of the heterostructure cavity configuration (cf. Fig. 2). An important aspect of our analysis is the inclusion of polariton interactions originating from nonlinearities of the dipolar phonons. In the bulk geometry, the importance of such nonlinearities in understanding the quantum phase transition has been discussed in Ref. [58]. In contrast, we focus on the cavity geometry and present a variational analysis to study effects of polariton interactions in such inhomogeneous setups. One of the most notable results is that one can find the cavityinduced transition into a ferroelectric phase even when the corresponding bulk system resides in the paraelectric phase. This finding suggests an intriguing possibility of inducing the superradiant quantum phase transition by utilizing the lightmatter interaction with the vacuum electromagnetic field. (Top) Setup of the heterostructure cavity configuration considered in the present manuscript. A slab of paraelectric material with a finite thickness is sandwiched between two semi-infinite metallic electrodes. Right panel represents spatial profiles of typical energy modes along the spatial direction perpendicular to the metalinsulator interfaces. (Bottom) Schematic figure illustrating three fundamental excitations relevant to the theory. A phonon excitation describes an optical collective mode induced by, e.g., the deformation of ionic crystals around the equilibrium configuration. More generally, it can also represent dipole moments of molecular solids or assembled organic molecules. The interaction between phonons and cavity photons is mediated by the dipole coupling, whose strength is represented by g. A plasmon, a collective excitation of electronic gases in metal mirrors, leads to the electromagnetic coupling characterized by the plasma frequency ωp.
In this context, one of the main novelties in the present work is an accurate analysis taking into account the dispersive nature of polariton excitations in the dielectric of finite thickness as well as the surface plasmonic effects occurring at the metal-insulator interfaces, i.e., a finite plasma frequency of the metallic cavity mirrors. We demonstrate that the proper consideration of these aspects is important to accurately understand several key features in the strong light-matter coupling regime, which are not readily attainable in the conventional bulk settings. We point out that the cavity-induced changes of the polariton dispersions provide a general route towards controlling energy-level structures not only in paraelectric materials, but also in many other systems possessing electric dipole moments, such as molecular solids, organic light-emitting devices, and excitonic systems. Our analysis can further be extended to include fabricated surface nanostructures [120,121]. In view of the versatility of the present formulation, we anticipate that our results will be useful to advance our understanding of a variety of physical systems lying at the intersection of condensed matter physics, quantum optics, and quantum chemistry.
The remainder of the paper is organized as follows. In Sec. II, we summarize our main results at a nontechnical level. In Sec. III, we present a theoretical framework for describing light-matter hybrid collective modes in both bulk and cavity settings. As a first step, we neglect phonon nonlinearities and determine dispersion relations of elementary excitations in the heterostructure cavity configuration. In Sec. IV, we include phonon nonlinearities based on the variational approach and analyze their effect on the paraelectric-to-ferroelectric phase transition in the cavity. We show that the mode confinement in the cavity geometry induces significant softening of polariton modes, indicating the enhanced propensity for ferroelectricity. We present a phase diagram of the paraelectric-to-ferroelectric transition in the cavity geometry. In Sec. V, we present a summary of results and suggest several interesting directions for future investigations.

II. SUMMARY OF MAIN RESULTS
We first summarize the main findings of this work at a nontechnical level before proceeding to the detailed theoretical formulation in the subsequent sections. The main focus of this paper is on many-body systems that are naturally close to the superradiant-type transition, i.e., paraelectric materials in the vicinity of the quantum phase transition into the ferroelectric phase (see Fig. 1(a-c)). For instance, SrTiO 3 and KTaO 3 are paraelectric at ambient conditions, but isotope substitution of 18 O for 16 O or applying pressure can bring the materials into a phase with spontaneous electric polarization [122].
The key property of such materials relevant to our analysis is the existence of a low-energy infrared optical phonon mode. A conventional way to approach the quantum critical point is to vary a tuning parameter (such as pressure or chemical composition) in order to decrease energy of this phonon mode. Transition into the broken symmetry phase occurs when the minimum energy of the phonon dispersion goes down to zero as indicated by the red dashed line in Fig. 1(a).
Our resonator setting facilitates hybridization of the electric dipole moment of the optical phonon with the quantized electromagnetic modes confined in the cavity; the resulting polariton modes further couple to plasmons in metal mirrors (see Fig. 2). A salient feature of our analysis is the emphasis on quantum nature of the light-matter coupling, which arises from strong zero-point fluctuations of the vacuum electromagnetic fields in the cavity, while a similar argument also applies to an equilibrium state at finite temperature. Placing the focus on equilibrium systems makes the present work distinct from the previous studies on materials subject to classical light waves, in which the use of strong external drive is essential. It is such quantum aspect of light-matter interaction that motivates us to ask two questions as follows: (A) How does the light-matter coupling modify polariton modes of paraelectric materials in the cavity geometry? (B) Is it possible to enhance the ferroelectric instability by modifying quantum electromagnetic environment of a many-body system under no external pump?
To address these questions, we formulate a simple yet generic theoretical framework including three essential ingredients: infrared active phonons, cavity photons, and plasmons in metal mirrors. To be concrete, we focus on the prototypical heterostructure configuration that consists of a paraelectric material sandwiched between large metallic electrodes (see the top panel in Fig. 2). One of the notable features of this system is that energies can be confined in volumes much smaller than the free-space wavelengths of photons, leading to strong enhancement of light-matter interaction [123]. Strong hybridization between the continuum of electromagnetic fields and matter excitations, including both phonons and plasmons, necessitates analyzing new collective modes that are absent in the bulk geometry. We assume that the size of the cavity is large enough for collective eigenmodes to be characterized by a conserved in-plane momentum. In this respect, the present consideration is distinct from most of the previous studies on the superradiant phases, where only one or at most a few modes in a cavity have been included.
In our theory, phonon modes possess electric dipole moments and are described by the real-valued vector quantum fieldφ. Physically, this field can correspond to the amplitude of the optical phonon, which describes a distortion of an ionic crystal leading to electric polarization. We note that the same formalism can be applied to systems with other polarizable modes such as collective dipole moments in assembled organic molecules or molecular solids. Interaction between polarizable modes and cavity electric fields is characterized by the light-matter dipole coupling whose strength is represented by g (see Sec. III for more detailed discussions). This coupling results in the formation of the collective excitations known as the phonon polaritons. Another class of elementary excitations included in our analysis is plasmons, i.e., collective excitations of electronic density in metal mirrors. Since plasmons describe fluctuations of electric charges, they are inherently coupled to the electromagnetic oscillations, which in the bulk system result in plasma resonance at frequency ω p . An important feature of plasmons relevant to our analysis is that they can form surface modes localized at the interfaces that are known as the surface plasmon polaritons [91]. This feature originates from the fact that electric fields can be reflected by the metal interfaces in such a way that the fields take maximum values at the boundaries. As schematically illustrated in Fig. 3, the resulting hybridized dispersions exhibit qualitative changes in the nature of cavity polariton excitations compared to the bulk material. Altogether, the richness of the cavity energy modes arises from the intrinsic hybridization of different types of excitations, including infrared phonons in the paraelectric, electromagnetic fields in the cavity, and plasmons in the metallic electrodes.
To answer the first question (A) in a concrete manner, we begin our analysis with a detailed discussion of the quadratic theory, which allows us to identify all the hybridized eigenmodes in the cavity geometry. Interestingly, even when the system is far from the transition point, the cavity confinement can qualitatively alter the nature of collective excitations (see e.g., Fig. 6 in Sec. III). We emphasize that, while our theory is formulated for a heterostructure cavity configuration, the present approach of controlling material excitations via modifying quantum electromagnetic environment can readily be generalized to a wide range of platforms. More specifically, besides quantum paraelectrics, we also discuss several potential applications of such cavity-based control of polariton modes to molecular materials and excitonic devices, such as improving the efficiency of organic light-emitting devices and realizing photon-exciton converters (see Sec. V and App. D). The inclusion of full complexities appropriate for specific experimental systems will require further analyses beyond the present work.
To address the second question (B), we demonstrate that the resonator setup in Fig. 2 provides yet another route towards controlling the phase of matter via its coupling to quantum fluctuations of electromagnetic fields. Building upon the identification of collective modes detailed in Sec. III, we show that the hybridization with the cavity electromagnetic fields can lead to significant softening of the dipolar phonon modes, and ultimately induce the paraelectric-to-ferroelectric phase transition when the softened phonon mode goes down to the zero frequency (see Fig. 3(b)). In particular, we find that this transition to the ferroelectric phase is possible even when the corresponding bulk system is still within the paraelectric phase. In free space, phonon nonlinearities render the formation of a spontaneous dipole moment energetically unfavorable. The role of the cavity here is to alleviate this adverse effect of phonon-phonon interactions, thus allowing one to induce ferroelectricity in quantum paraelectric materials that otherwise exhibit no spontaneous polarizations even at zero temperature. Our results thus indicate an intriguing possibility of inducing the genuine superradiant-type quantum phase transition via the light-matter coupling with no external pump.
We note that this phenomenon is not present at the level of a quadratic model, and can only be found when phonon nonlinearities are included into the analysis. In the case of bulk quantum paraelectrics with no cavity confinement, the crucial role of nonlinear effects at the verge of the quantum phase transition has previously been discussed in Ref. [58]. We take into account the nonlinear effects in the cavity geometry on the basis of the variational approach, leading to renormalization of the hybridized collective excitations. The shift of the paraelectric-to-ferroelectric transition point can directly be related to this renormalization of energy dispersions. We demonstrate that the phenomenon of cavity-enhanced ferroelectricity is more pronounced when confinement effects are enhanced by using a cavity geometry with a thinner thickness or a higher plasma frequency (see Fig. 9(b) in Sec. IV).

III. PHONON-PHOTON-PLASMON HYBRIDIZATION
Our main goal is to understand the qualitative physics of a dipolar insulator coupled to quantized electromagnetic fields in a cavity rather than to make predictions specific to particular experimental setups. In this section, we outline the general theory of the light-matter interaction and formulate a framework to reveal generic features in the hybridization of collective dipolar modes, cavity photons, and plasmons in metal mirrors. While the theory does not feature full complexities of an experimental system, it applies to a broad range of materials supporting a polarizable vector field as matter excitations. In the next section, we will apply the framework to an ionic crystal confined in the cavity as a concrete case and analyze its paraelectric-to-ferroelectric phase transition.

A. Bulk
Before moving on to the case of the cavity setting, we first discuss the three-dimensional bulk setup with no boundary effects and determine the bulk dispersion relations for the phonon-photon hybridized modes, i.e., the phonon polaritons. Later, we generalize the theory to the cavity configuration by taking into account confinement effects and also contributions from plasmons in metal mirrors.
The total system of the bulk setting is governed by the HamiltonianĤ where the first (second) term describes the dynamics of the electromagnetic field (dipolar matter field), and the last term represents the light-matter coupling between these two fields.
We specify each of those terms below.

Electromagnetic field
The free evolution of the electromagnetic field is described by the Hamiltonian where 0 is the vacuum permittivity, c is the vacuum light speed,Â(r) is the vector potential of photons, andΠ(r) is its canonically conjugate variable. The integral is taken over a cubic whose volume is V = L 3 and we choose the Coulomb gauge Assuming the periodic boundary conditions and introducing discrete wavenumbers k i = 2πn i /L with i = x, y, z and n i being integer numbers, we represent the vector fields in terms of their Fourier components aŝ where ω k = c |k| is the vacuum photon dispersion, andâ kλ (â † kλ ) is the annihilation (creation) operator of photons with momentum k and polarization λ. The Coulomb gauge (3) leads to the transversality condition of the orthonormal polarization vectors kλ : The Hamiltonian can be diagonalized in the basis of momentum k and polarization λ aŝ where we omit the vacuum energy.

Matter field
We consider a dipolar insulator whose excitation is represented by a real-valued vector fieldφ(r). For instance, it can describe an infrared optical phonon mode induced by the collective deformation around the equilibrium configuration in ionic crystals such as MgO, SiC and NaCl [124]. In other examples, it can represent collective dipolar modes associated with assembled organic molecules, molecular crystals, and excitonic systems. We focus on a nondispersive matter excitation with frequency Ω as appropriate for, e.g., an optical phonon in dipolar insulators or a collective dipolar mode in molecular materials. The effective Hamiltonian of the matter field is thus given bŷ whereĤ 0 is the quadratic part andĤ int is the most relevant interaction Here,π(r) is a conjugate field of the matter fieldφ(r), v is the unit-cell volume in the insulator, M is the effective mass, which corresponds to the reduced mass of two ions in the case of ionic crystals, and U ≥ 0 characterizes the interaction strength. Physically, the interaction term can originate from, e.g., phonon anharmonicity. When we will discuss in Sec IV the transition between the paraelectric phase with φ = 0 and the ferroelectric phase with φ = 0, the phonon frequency Ω 2 should be considered as the tuning parameter governing the transition in the effective theory (cf. Fig. 1). We note that in the case of ionic crystals the present model neglects the coupling between electrons and soft phonons, which is typically considered as weak [125]. The matter fieldφ(r) and its conjugate variableπ(r) can be quantized through the canonical quantization. In the Coulomb gauge, it is useful to decompose them in terms of the transverse and longitudinal components asφ(r) = φ (r) +φ ⊥ (r) andπ(r) =π (r) +π ⊥ (r), leading to the expressionŝ where N = V /v is the number of modes,φ ⊥ kλ (φ ⊥ † kλ ) is the annihilation (creation) operator of transverse phonons with momentum k and polarization λ, andφ k (φ † k ) is the annihilation (creation) operator of longitudinal phonons. We recall that the transverse polarization vectors kλ satisfy Eqs. (6) and (7). At the quadratic level (U = 0), only the transverse components couple to the dynamical electromagnetic field while the longitudinal components couple to the static longitudinal electric field as discussed below. In the similar manner as in the electromagnetic HamiltonianĤ light , the quadratic partĤ 0 of the matter HamiltonianĤ matter can be diagonalized aŝ

Light-matter coupling
The transverse dipolar matter fieldφ ⊥ can directly couple to the electromagnetic field in a bilinear form. For instance, we recall that the fieldφ ⊥ can represent transverse displacements of ionic charges with opposite signs and thus can couple to the electric field. In terms of the vector potentialÂ, the light-matter interaction is given bŷ where Z * e is the effective charge of a dipolar mode, which corresponds to the charge of ions in the case of ionic crystals, and V C is the Coulomb potential term whose explicit form is specified below. We note that, aside a constant factor, the conjugate fieldπ physically corresponds to the current operator j ∝ ∂ t φ associated with the dipolar fieldφ. We here invoke neither the dipole approximation nor the rotating wave approximation that are often used in literature, but can break down especially in strong-coupling regimes [81,126]. The light-matter interaction (17) including botĥ π ·Â and A 2 terms allows one to explicitly retain the gauge invariance of the theory. It is also noteworthy that the lightmatter coupling (17) is nonvanishing even in the absence of photons. This interaction with the vacuum electromagnetic field originates from the zero-point fluctuations of the electric field and can be viewed as the light-matter interaction mediated by virtual photons. In this paper, we focus on such a low photon, quantum regime with no external pump.
Our choice of the Coulomb gauge leads to the nondynamical contribution V C resulting from the constraint relating the longitudinal component of the electric field to that of the matter field. To see this explicitly, we first decompose the electric field in terms of the transverse and longitudinal components asÊ which satisfy ∇ ·Ê ⊥ = 0 and ∇ ×Ê = 0. On one hand, the transverse part can be directly related to the vector potential On the other hand, the longitudinal part is subject to the constraint resulting from the Coulomb gaugê where we recall thatφ is the longitudinal part of the matter field satisfying ∇ ×φ = 0. This longitudinal component leads to the expression of the Coulomb potential term V C as follows [127]:

Bulk dispersions
One can diagonalize the quadratic part of the total Hamiltonian,Ĥ which allows one to obtain an analytical expression of the bulk dispersions. The presence of light-matter interaction leads to the formation of the hybridized modes consisting of photons and phonons, known as the phonon polaritons. This hybridization strongly modifies the underlying dispersion relations such that the avoided crossing occurs between two polariton dispersions. At the quadratic level, we can decompose the quadratic HamiltonianĤ tot,0 in Eq. (22) into two parts that solely include either transverse or longitudinal components of the field operators as follows: To diagonalize the transverse part of the quadratic Hamilto-nianĤ ⊥ tot,0 , it is useful to introduce the conjugate pairs of variables for the field operatorsπ ⊥ andÂ in the Fourier space: where the superscript P indicates that an operator is a conjugate variable of the corresponding operator, i.e., each pair of conjugate variables satisfies the canonical commutation relations. The transverse part of the quadratic Hamiltonian can then be simplified aŝ We here introduce the coupling strength g as The dispersion relations of the transverse modes can now be obtained by diagonalizing the matrix in the second term in Eq. (28). The resulting Hamiltonian for the transverse components isĤ where s = ± denotes an upper or a lower bulk dispersion with eigenfrequencies ω [Ω/c] k andγ ⊥ ksλ (γ ⊥ † ksλ ) represents an annihilation (creation) operator of a transverse phonon-photon hybridized mode with momentum k, polarization λ, and mode s. The solid curves in Fig. 4 show these two branches of the bulk phonon polaritons. The splitting between the upper and lower dispersive modes is characterized by the coupling strength g. The upper branch exhibits a quadratic dispersion at low k = |k| and is gapped, i.e., it takes a finite value Ω 2 + g 2 at k = 0, while the lower one linearly grows from zero and saturates at Ω even in the limit of k → ∞.
Meanwhile, the longitudinal part of the total quadratic HamiltonianĤ tot,0 only contains the longitudinal matter field φ and can be readily diagonalized aŝ where Ω = Ω 2 + g 2 is the longitudinal phonon frequency [128]. We note that the longitudinal mode is independent of the momentum k and is thus nondispersive as shown in the solid horizontal line in Fig. 4.

B. Cavity configuration
We next consider the cavity setting and take into account confinement effects as well as the hybridization of electromagnetic fields with plasmons in metal mirrors in addition to collective dipolar modes. We focus on the prototypical heterostructure configuration consisting of a dipolar insulator sandwiched between two semi-infinite ideal metals (see Fig. 5). We consider eigenmodes propagating with the twodimensional in-plane momentum q while the formalism developed in this section can be extended to deal with more complex setups such as fabricated metal surface structures.
Thin-film configurations considered here have been previously applied to light-emitting devices [129,130], the electron-tunneling emission [131], and the transmitting waveguide [132]. One of the main novelties introduced by this paper is to reveal the great potential of heterostructure FIG. 5. Schematic figure illustrating the cavity configuration considered in this paper. A dipolar insulator whose center is positioned at z = 0 is sandwiched between two semi-infinite metals with interfaces at z = ±d/2. configurations as a quantum electrodynamic setting towards controlling collective matter excitations and, in particular, to point out a possibility of inducing the superradiant-type quantum phase transition via the vacuum electromagnetic environment. More specifically, analyzing the hybridization of cavity electromagnetic fields and matter excitations, we demonstrate that the present cavity architecture enriches the underlying dispersion relations in a way qualitatively different from the corresponding bulk case. Including phonon nonlinearities, we will further show that the heterostructure configuration enables one to significantly soften the phonon modes, which ultimately induces the structural phase transition.

Hamiltonian
To reveal essential features of the present cavity setting, we consider an ideal metal with the Drude property whose plasma frequency is denoted as ω p = n e e 2 /(m e 0 ), where e is the elementary electric charge, n e is the electron density in the metal, and m e is the electronic mass. We assume that all the materials are nonmagnetic such that the relative magnetic permeability µ is always taken to be µ = 1. As illustrated in Fig. 5, we consider a setup with the center of a dipolar insulator being positioned at z = 0 while the insulator-metal interfaces being positioned at z = ±d/2. We again choose the Coulomb gauge ∇ ·Â = 0 and thus the HamiltonianĤ light of the free electromagnetic field is the same as in the bulk case (cf. Eq. (2)) while the matter Hamiltonian is modified aŝ where we define the integral over the insulator region as The light-matter coupling term is also modified aŝ where the first term describes the plasmon coupling in the metal regions defined as dz, and the second term is the coupling between the insulator and the cavity electromagnetic field. We recall that the Coulomb gauge leads to the additional constraint on the longitudinal components of the electric and matter fields in the insulator region |z| < d/2:Ê As in the bulk case, this constraint results in the Coulomb potential term V C originating from the longitudinal dipolar modes as follows: Summarizing, the present cavity light-matter system consists of two canonically conjugate pairs of variables (Π,Â) and (π,φ); the former describes the dynamical electromagnetic degrees of freedom and only have transverse components while the latter describes matter excitations and have both transverse and longitudinal parts. Their time evolution is governed by the Hamiltonian consisting of Eqs. (2), (33) and (36), and is subject to the constraint (37).
We note that the use of the Coulomb gauge remains to be a convenient choice even in the present case of the inhomogeneous geometry. An alternative naive choice would be the condition ∇ · [ (r, ω)Â(r, ω)] = 0 with being the relative permittivity in the medium, which might be viewed as a quantum counterpart of the gauge constraint in the macroscopic Maxwell equations. However, it has been pointed out that in inhomogeneous media this gauge choice suffers from the difficulty of performing the canonical quantization of the variables due to the frequency-domain formulation of the gauge constraint [127]. In contrast, the Coulomb gauge chosen here can be still imposed on a general inhomogeneous system without such difficulties; we thus employ it throughout this paper.

Elementary excitations
The quadratic part of the total cavity Hamiltonian can still be diagonalized, which allows us to identify the hybridized eigenmodes and the corresponding elementary excitation energies. To this end, we first expand the electromagnetic vector potential in terms of the transverse elementary modes. Owing to the spatial invariance on the xy plane, we obtain where A = L 2 is the area of interest on the xy plane, q is the two-dimensional in-plane wavevector q = (q x , q y ) T , and ρ = (x, y) T is the position vector on the xy plane. The transverse mode function U ⊥ qnλ satisfies ∇ · U ⊥ qnλ = 0. Because the spatial invariance along the z direction is lost in the heterostructure configuration, an elementary excitation is now labeled by a discrete number n ∈ N (instead of q z in the bulk case) as well as the two-dimensional vector q. We also introduce a discrete variable λ = (Λ, P ) that includes the polarization label Λ ∈ {TM, TE} and the parity label P = ±. Here, the polarization label Λ in the discrete subscript λ indicates either the transverse magnetic (TM) or electric (TE) polarization associated with the zero magnetic or electric field in the z direction: Meanwhile, the parity label P indicates the symmetry of the mode function with respect to the spatial inversion along the z direction. From the parity symmetry of the Maxwell equations, it follows that the z component of the mode function must have the opposite parity to that of the x, y components (see Appendix B). We thus assign the parity label P to each mode function according to Finally, we can impose the orthonormal conditions on the mode functions as We can expand the matter degrees of freedom in the similar manner as in the electromagnetic field. Because the matter field resides only in the insulator region, it is natural to expand the transverse components of the dipolar matter field aŝ where N d = A/v = N/d is the number of modes per thickness, and {f ⊥ qnλ } is a complete set of transverse orthonormal functions in the insulator region satisfying ∇ · (f ⊥ qnλ e iq·ρ ) = 0. We note that a position variable r in the matter fieldφ should be interpreted as the parameters restricted to the region |z| < d/2. In the discrete variable λ = (Λ, P ), the polarization label takes either TM or TE modes, Λ ∈ {TM, TE}, and the parity label P = ± is defined in the same manner as in Eqs. (43) and (44). Another possible matter excitations are the longitudinal modes, for which we can expand the matter field aŝ where {f qnλ } is a complete set of longitudinal orthonormal functions in the insulator region, which satisfies ∇ × (f qnλ e iq·ρ ) = 0. We denote this mode by the label λ = (L, P ), where the polarization label Λ is fixed to be the longitudinal one Λ = L while the parity P = ± of the mode function is defined in the same way as in the transverse case. The longitudinal mode associates with an excitation of the longitudinal electric field via Eq. (37). An explicit form of the transverse mode functions {u ⊥ qnλ } and {f ⊥ qnλ }, the longitudinal mode functions {f qnλ }, and the corresponding excitation energies {ω qnλ } can be determined from solving the eigenvalue problem with imposing suitable boundary conditions at the interfaces as detailed in Appendix B. In the transverse modes, we note that both matter and light fields as well as the hybridized eigenmodes share the same spatial profile represented by u ⊥ qnλ owing to the ω  , the TM1,2 modes represent the two highest modes in the TM sector, the SS (=TM0) mode represents the symmetric surface (SS) mode, the AS1,2 modes represent the antisymmetric surface (AS) modes whose dipolar moments point out of the plane at low q. The longitudinal L± modes have the same nondispersive energy independent of the momentum q and the parity P = ±. The assignment of the label TM0 to the SS mode is in accordance with the convention used in studies of plasmonics [133]. In (d) and (e), the RTMn modes with n = 0, 1, 2 . . . represent the resonant TM modes close to the dipolar-phonon frequency Ω. The RTM n≥1 modes extend over the insulator region while the most-softened mode RTM0 is essentially the surface mode localized at the interfaces analogous to the SS mode. In (c), the TE0,1 modes represent the two highest modes in the TE sector. In (f), the RTEn modes with n = 0, 1, 2 . . . represent the resonant TE modes; all of them extend over the insulator region. The parameters are d = c/Ω, g = 3.5 Ω, and ωp = 5 Ω. For the sake of visibility, the plasma frequency is set to be a rather low value while its specific choice will not qualitatively affect the influence of cavity on the ferroelectricity as discussed in Sec. IV.
bilinear form of the light-matter coupling. Meanwhile, the longitudinal modes consist of only the matter field and their spatial profiles are characterized in terms of the mode functions f qnλ . It is also noteworthy that an excitation energy ω qnλ of an elementary eigenmode now explicitly depends on the polarization Λ ∈ {TM, TE, L} in contrast to the bulk case in Eq. (31). In terms of the obtained elementary eigenmodes, the quadratic part of the total light-matter Hamiltonian,Ĥ tot,0 =Ĥ light +Ĥ 0 +Ĥ l−m , can be diagonalized aŝ whereγ qnλ (γ † qnλ ) represents an annihilation (creation) operator of a hybridized elementary excitation with in-plane wavevector q and discrete labels n and λ. The spatial profile of each excitation is characterized by an eigenmode function u ⊥ qnλ or f qnλ depending on the polarization, whose explicit functional forms are given in Appendix B. Figure 6 shows the elementary excitations in the cavity setting for each sector of λ = (Λ, P ), and Table I Fig. 6. The second and third columns represent the parity P = ± and the polarization Λ ∈ {TM, TE, L} of each mode, respectively. The fourth column indicates the direction of the dipole moments in the insulator at low in-plane wavevector q with respect to the plane parallel to the interfaces. The fifth column shows whether or not the divergence of the electric field is vanishing. The sixth and seventh columns represent whether the longitudinal wavenumber κη takes a real value or a purely imaginary value at low and high q, respectively. The subscript η ∈ {q, n, λ} with λ = (Λ, P ) represents a set of all the variables specifying each eigenmode. The final two columns indicate the physical nature of each hybridized band at low and high q, where we abbreviate elementary excitations as plasmons (PL), photons (PT), phonons (PN), polaritons (P), and specify surface modes by S. First of all, a salient feature of the elementary excitations in the cavity setting is the emergence of a large number of soft-phonon modes lying below the bare phonon frequency Ω, which we term as the resonant TM (RTM) modes (see Fig. 6(d,e)). The RTM modes are labeled by an integer number n = 0, 1, 2 . . . as RTM n . The parity P of the corresponding mode function u ⊥ qnλ coincides with the parity of n. The RTM n modes with n ≥ 1 lie between Ω and the lower bulk dispersion (plotted as the lower black dotted curve in Fig. 6), and share the common structures. For instance, all of the RTM n≥1 modes start from frequencies slightly below Ω and saturate to it at high q ≡ |q|. Thus, these modes predominantly consist of the hybridization of dipolar phonons and cavity electromagnetic fields at low q while the phonon contribution eventually becomes dominant at high q. The field amplitudes of these RTM n≥1 modes extend over the insulator region owing to the real-valuedness of the longitudinal wavenumber κ η over the entire plane of the in-plane wavevector q. Here, the subscript η is used to summarize all the indices as η = {q, n, λ}. The origin of the emergence of many RTM modes below the bare phonon frequency can be understood from simple physical arguments based on the frequencydependent dielectric functions (see Appendix C).
Remarkably, the lowest resonant mode RTM 0 , which is the most-softened phonon mode among all the elementary excitations in all the sectors, has a qualitatively different physical nature compared to the RTM n≥1 modes discussed above. Its longitudinal wavenumber κ η changes from a real value to a purely imaginary one as the dispersion crosses the lower bulk dispersion in the bulk (see Fig. 6(d)). Moreover, the saturated frequency value of the RTM 0 mode at q → ∞ remains substantially below Ω with a nonvanishing gap. Accordingly, the RTM 0 mode can be interpreted as the low-energy analogue of the surface polariton mode labeled as the SS mode in Fig 6(a); both modes result from the intrinsic hybridization of dipolar phonons, cavity electromagnetic fields, and plasmons while the major contribution at high q comes from the phonon (plasmon) part in the case of the RTM 0 (SS) mode.
To gain further insights, in Fig. 7 we plot the spatial profiles of the electric field E of the RTM 0 mode. Figures 7(a) and (b) show the electrical flux lines of the RTM 0 mode at low and high in-plane momenta q, respectively, with the color indicating the magnitude of the out-of-plane component E z of the electric field. We note that our convention of denoting the z-component as the out-of-plane one is in accordance with the planar cavity setting whose metal-insulator interfaces are parallel to the xy plane. Without loss of generality, here we focus on the mode propagating in the x direction, i.e. q = q x , leading to the condition E y = 0 (see Appendix B). Figures 7(c) and (d) show the low-q and high-q spatial profiles of the inplane (red solid curve) and the out-of-plane (blue solid curve) components of the electric field at x = 0. At low q before crossing the lower bulk dispersion, the field predominantly points to the in-plane direction and extends over the insulator region |z| ≤ d/2 (Fig. 7(c)). In contrast, after crossing the bulk dispersion, the longitudinal wavenumber κ η turns out to take a purely imaginary value and the field profiles qualitatively change. Specifically, at high q the fields are localized around the interfaces and the z component of the electric field exhibits significant discontinuities due to the surface charges concentrated on the metal interfaces ( Fig. 7(d)); these features are characteristic of surface polariton modes. As shown later, it is the softening of the RTM 0 mode discussed here that will trigger the structural instability of materials, thus ultimately inducing the paraelectric-to-ferroelectric phase transition. When the plasma frequency is close to the phonon frequency and the hybridization with plasmons becomes particularly prominent, the above crossover behavior of the RTM 0 mode also manifests itself as the emergence of roton-type excitations. To see this explicitly, in Fig. 8 we show the energy dispersion of the RTM 0 mode at different light-matter coupling strengths g with a low plasma frequency. It is evident that the energy dispersion takes the minimum value at finite in-plane momentum q * > 0, indicating the formation of roton excitations. Physically, this minimum roughly corresponds to the vanishing point of the longitudinal wavenumber κ η = 0, at which the crossover from low to high in-plane momentum behavior occurs (compare Figs. 7(c) with (d)). While a specific value of the finite in-plane momentum q * depends on system parameters such as the thickness of the insulator, it is a general feature that roton excitations eventually disappear (i.e., q * → 0) as the plasmon component in the eigenmode becomes less significant by, for instance, increasing either the light-matter coupling strength g or the plasma frequency ω p .
We now proceed to discuss physical properties of the other elementary modes. The highest band (denoted as TM 2 in Fig. 6(a)) starting from the plasma frequency ω p at q = 0 is essentially the upper branch of the standard plasmon polaritons. This hybridized mode reduces to the bare plasmon (photon) excitation in the limit of q → 0 (q → ∞). The next highest band (denoted as TM 1 in Fig. 6(b)) initiates from a frequency slightly below the plasma frequency and has a physical origin similar to the TM 2 mode, but there still exists a difference that the TM 1 mode sustains the nonvanishing hybridization with the electromagnetic fields at low q. The mode functions u ⊥ η of the TM 1,2 modes and the corresponding dipole moments predominantly point to the directions parallel to the interfaces at low q. The longitudinal wavenumber κ η remains real over the entire in-plane momentum space, and thus the amplitudes of the TM 1,2 modes extend over the insulator region (see also Fig. 11 in Appendix B).
We next discuss the symmetric and antisymmetric surface modes lying in the frequency regime Ω 2 + g 2 < ω < ω p , which are denoted as SS(=TM 0 ) and AS 2 in Figs. 6(a) and (b), respectively. They essentially arise from the split of the coupled surface phonon polaritons at two metal-insulator interfaces [134]. We remark that the assignment of the label TM 0 to the SS mode is in accordance with the convention of studies in plasmonics [133]. At a high in-plane wavevector q, the SS and AS 2 modes become bound modes that are localized around the interfaces; the field amplitudes take maxima at the interfaces and exponentially decay away from them. In the limit of q → ∞, the energy dispersions for both of the SS and AS 2 modes become flat and saturate to the same, finite frequency below ω p . Meanwhile, at q = 0 the AS 2 mode starts from a frequency lower than that of the SS mode, and thus the former exhibits a larger dispersion. The SS mode can have a real longitudinal wavenumber κ η at low q and thus it can be radiative, i.e., it can be coupled to light modes outside the cavity when metal mirrors have finite thicknesses. As a result, the SS mode can contribute to optical properties of the system. With increasing q, the longitudinal wavenumber κ η changes from a real value to a purely imaginary value and the SS mode changes into the nonradiative mode. In contrast, the AS 2 mode is nonradiative over the entire in-plane momentum and can potentially be applied to nanoscale plasmonic waveguides owing to its localized and dispersive nature. The difference between two surface modes at low q also results in their qualitatively different low-energy polarization behavior; the dipole moment of the SS mode points to the in-plane direction while that of the AS 2 mode is purely out-of-plane. The lowest mode in the TM polarization, which is denoted as AS 1 , lies slightly below the lower bulk dispersion (black dotted curve) and converges to the linear photon dispersion at q → 0 in contrast to the other modes featuring quadratic asymptotes. This mode essentially shares the similar properties with the AS 2 mode; both of them are localized modes over the entire in-plane momentum q and point out of the plane at low q due to the surface charges. The AS 1 mode asymptotically becomes a purely photon mode at q → 0.
We note that these surface modes result from the intrinsic hybridization of the collective dipolar mode, electromagnetic fields in the cavity, and plasmons in metal mirrors, and can emerge only in the case of the TM polarization, in which the continuity of the z component of the electric field can be mitigated. Thus, the surface modes including the SS(=TM 0 ), AS 1,2 , and RTM 0 modes can be excited only by p-polarized light. It is worthwhile to mention that in the context of plasmonics the nanoscale localization of such surface modes has previously found applications to manipulate light waves below the diffraction limit [135,136]. The longitudinal modes L ± in Fig. 6(a) have the parity P = ± and lie at the constant frequency Ω = Ω 2 + g 2 . These nondispersive modes arise from the vanishing insulator permittivity and is characterized by the nonzero divergence of the mode function ∇ · (e iq·ρ f η ) = 0. Both the in-plane and out-of-plane amplitudes of the longitudinal electric field extend over the insulator region while the amplitudes vanish in the metal regions due to the boundary conditions. At low q, the out-of-plane component becomes dominant in the insulator region due to the surface charges.
Finally, Fig. 6(c) plots the TE-polarized bands of both parities P = ±. The corresponding magnified plot around the phonon frequency is also shown in Fig. 6(f). These TE modes have simpler structures than the TM ones, i.e., each of all the modes lies slightly above the higher or lower branches of the bulk dispersions (black dotted curves) due to energy enhancements by the cavity confinement. We label the higher set of the TE modes as TE 0,1 while the lower set of a large number of resonant TE modes as RTE n=0,1,2... . Because the mode functions (and thus the electric fields) for all the modes are continuous even at the interfaces, the origin of the enhanced energies in these TE modes can simply be understood as the acquisition of nonzero longitudinal momentum associated with the cavity confinement. Indeed, the corresponding longitudinal wavenumber κ η remains real for all the TE modes over the entire region of the in-plane momentum q. This real-valuedness of κ η also indicates that the TE modes are essentially radiative photonic modes that can directly couple to light modes outside the system when the thicknesses of the metal mirrors are taken to be finite.
Several remarks are in order. Firstly, we here emphasize the importance of explicitly including plasmon contributions into the analysis. This treatment ensures that the tails of the evanescent fields properly extend into the surrounding metal mirrors. In many of the previous studies on materials or assembled molecules confined in the quantum electromagnetic environment, the mirrors are assumed to be perfectly reflective such that the electric fields must exactly vanish at the interfaces and are absent in the metal regions. In such analyses, the physical nature of a multitude of the localized surface modes discussed here cannot be addressed appropriately. For instance, the most-softened phonon mode RTM 0 is one important example of such localized modes; an analysis assuming the perfect reflectivity cannot capture its unique features, including the qualitative changes associated with increasing inplane momentum q (cf. Fig. 7), the significant nonvanishing softening from the bare phonon frequency Ω even in the limit q → ∞ (cf. Fig. 6(d)), as well as the emergence of roton-type excitations (cf. Fig. 8). Importantly, the correct treatment of this RTM 0 mode turns out to be crucial especially when we discuss the superradiant-type quantum phase transition as detailed in the next section. Meanwhile, we remark that the reflectivity does not qualitatively change physical properties of the TE modes, which are analogous to the modes supported in Fabry-Perot cavities [137].
Secondly, we note that the present results can readily be generalized to include further complexities in practical experimental situations. In particular, effects from losses due to electron damping can be included by introducing phenomenological imaginary parts in the permittivity of the metal mirrors. One can also use a mathematically equivalent, yet more rigorous approach by, e.g., considering coupling to bosonic reservoir modes or relying on the macroscopic Green's function formalism and the microscopic Hopfield approach. Yet, the effects found here can in principle outweigh the losses owing to the strong cavity confinement.
Finally, we remark that the rich structures of the elementary excitations demonstrated above can further be tuned by varying the thickness d of the dipolar insulator. The increase of d generally leads to a generation of more delocalized bulk modes extending over the insulator. Thus, aside from the localized surface modes, the eigenmodes become denser and eventually converge to the corresponding bulk dispersions in the limit of d → ∞. In contrast, the decrease of d leads to the filtering of the delocalized modes; for instance, a degenerate frequency of the TE 0 and SS modes at zero-in-plane momentum shifts up for a smaller d. Thus, with decreasing d, the localized nature of the modes plays more and more crucial roles in determining the physical properties of the present heterostructure configuration. We assess the effects of changing the thickness in more detail in the next section by including the most relevant nonlinear effect in a self-consistent manner. The nonlinearity renormalizes the value of the bare phonon frequency Ω and thus effectively modifies the underlying elementary excitations. We will see that a thinner insulator is more preferable to realize the cavity-enhanced ferroelectricity. The reason is that a thin heterostructure configuration can alleviate the adverse effect of the interaction among the energetically dense phonon excitations and is thus advantageous for inducing the ultimate softening of the RTM 0 mode, which causes the ferroelectric instability or, namely, the superradiant-type phase transition.

IV. CASE STUDY OF CAVITY-ENHANCED FERROELECTRIC PHASE TRANSITION
As a concrete case study of the theoretical formulation in the previous section, we here discuss the influence of the quantum electrodynamic environment on the paraelectric-toferroelectric phase transition, which is the structural quantum phase transition commonly observed in ionic crystals. Most notably, we find that in the cavity setting the transition point can be shifted in favor of the ferroelectricity owing to the softening of the dipolar phonon mode. This softening results from the interplay between the modified dispersions and feedback effects through nonlinear interactions between the energetically dense elementary excitations.

A. Variational principle
We are interested in how the coupling between the dipolar phonon modes and the light modes confined in the cavity affects physical properties of quantum paraelectrics on the border of ferroelectric order. The paraelectric-to-ferroelectric phase transition is essentially a structural transition of the lattice structure in, e.g., ionic crystals. The phonon fieldφ plays the role of the order parameter, which spontaneously breaks the symmetry in the ordered, ferroelectric phase. In ionic crystals, this order is supported by the imbalanced charge configuration with the relative displacement, which accompanies a uniform and nonzero electric polarization (cf. Fig. 1(b)). Thus, the transition is necessarily associated with the softening of an optical phonon mode while the acoustic phonon mode is irrelevant.
We thus consider an effective model described by the total HamiltonianĤ whereĤ matter now includes both the quadratic terms (34) and the most relevant interaction term (35),Ĥ light is the Hamiltonian for the free electromagnetic field (2), andĤ l−m describes the light-matter coupling in the cavity setting as given by Eq. (36). We note that here it is crucial to go beyond the linear regime and to include nonlinear effects because the system is at the verge of the phase transition point and fluctuations can be significant.
To qualitatively understand physical consequences of the cavity confinement, we perform the variational analysis on the present nonlinear model. Specifically, we first replace the bare phonon frequency Ω 2 inĤ 0 (cf. Eq. (34)) by a system tuning parameter r, where r ≥ 0 (r < 0) corresponds to the paraelectric (ferroelectric) phase at the level of the noninteracting linear model. In practice, this tuning parameter can be varied in various manners such as by changing densities, chemical isotopic substitutions, or applying hydrostatic pressure. Accordingly, we denote the quadratic part of the total Hamiltonian asĤ tot,0 (r) and the full interacting Hamiltonian asĤ tot (r) =Ĥ tot,0 (r) +Ĥ int . We then represent the corresponding exact thermal equilibrium state aŝ where Z(r) ≡ Tr[e −βĤtot(r) ]. We then find the best approximation for this interacting equilibrium state by using a tuning variable as a variational parameter r eff : where D KL [·|·] is the Kullback-Leibler divergence [138], and an effective Gaussian state is defined aŝ with Z 0 (r eff ) ≡ Tr[e −βĤtot,0(r eff ) ].
We introduce the variational free energy as where we denote an expectation value with respect to the effective density matrix as · · · eff = Tr[ρ 0 (r eff ) · · · ], and F 0 (r eff ) = − 1 β ln Z 0 (r eff ) is the Gaussian variational free energy. The variational condition (51) can then be simplified as [139] min r eff where F (r) = − 1 β ln Z(r) is the exact free energy. The optimal effective tuning parameter r * eff , which is defined by gives the best approximation of the interacting model in terms of the linear theory. Within the present variational formalism, the phase transition point is determined from the condition r * eff = 0 that can eventually be met as the bare tuning parameter r is decreased from a value corresponding to the paraelectric phase. In other words, the interaction feedback is introduced as the renormalization of the effective phonon frequency r * eff = Ω 2 eff in the variational theoryĤ tot,0 (r * eff ), and the transition occurs when this renormalized frequency of paraelectric materials goes down to zero. We note that in the present analysis the expectation value can be decoupled as in the standard mean-field analysis, and the theoretical procedure can be exact in the spherical model, where the number of components in the order parameter is taken to be infinite.

B. Results
We show the obtained variational free energy F v against the variational parameter r eff at different insulator thicknesses d in Fig. 9(a). From top to bottom panels, we decrease the bare tuning parameter r. The thermal equilibrium value r * eff of the effective phonon frequency is determined from identifying the minimum of each variational free-energy landscape. In all the panels, a thinner insulator thickness d in the cavity setting leads to a smaller equilibrium phonon frequency r * eff , indicating the phonon softening owing to the coupling with light modes strongly confined in the cavity. With further decreasing the thickness d or the tuning parameter r, the equilibrium phonon frequency eventually goes down to zero from above, r * eff → 0 + , at which the paraelectric-to-ferroelectric phase transition occurs (see e.g., the red solid curve in the bottom panel of Fig. 9(a)).
The corresponding phase diagram of the cavity setting with varying thickness d is plotted in Fig. 9(b), where the cavityinduced softening results in favoring of the ferroelectric phase compared to the bulk transition point indicated by the black dashed line. This favoring of the ferroelectric phase in the quantum electromagnetic environment qualitatively persists almost independently of specific choices of the parameters FIG. 9. (a) Variational free energies Fv plotted against the variational parameter r eff with decreasing tuning parameter r from top to bottom panels. The red and blue solid curves correspond to the results for the cavity setting with the insulator thicknesses d = 1 and d = 3, respectively, while the green solid curves correspond to the bulk results. The minima of the variational energies determine the equilibrium values r * eff = Ω 2 eff of the effective phonon frequency including the nonlinear feedback. (b) Phase diagram in the cavity setting. As the tuning parameter r is decreased, the transition from the paraelectric phase (PE) to the ferroelectric phase (FE) occurs when the equilibrium effective phonon frequency r * eff goes down to zero. The transition point in the cavity setting shifts in favor of FE compared to the bulk transition point indicated by the black dashed line. The inset shows the dependence on the plasma frequency ωp at the insulator thickness d = 1. We set c = = β = 1 (corresponding to, e.g., c/(kBT ) 38 µm for T = 60 K), and use ωp = 5, U = 1, and g = 3.5. We set the momentum cutoff to be Λc = 10 2 . and should be a generic feature of various dipolar insulators confined in the cavity. Most notably, this robust effect occurs even when the underlying phase in the bulk setting is paraelectric, thus indicating a possibility of inducing the structural quantum phase transition via the cavity confinement as indicated by the vertical black downarrow in Fig. 9(b). The higher plasma frequency leads to faster spatial decay of low-energy modes inside the metallic regions and thus corresponds to tighter mode confinements. Hence, increasing the plasma frequency has qualitatively similar effect to decreasing the paraelectric slab thickness (see the inset in Fig. 9(b)).
To gain further insights into the phonon softening, in Fig. 10 we plot the most-softened phonon mode (i.e., the RTM 0 mode) with varying the insulator thickness d but at the fixed bare tuning parameter r. It is this softening that destabilizes the paraelectric phase when it goes down to zero energy, triggering the phase transition. Thus, the unstable dipolar phonon mode at the verge of the transition is precisely the transverse one and points to the in-plane direction (cf. Fig. 7(a,c)). We emphasize that the results in Fig. 10 are obtained by using the renormalized effective phonon frequency Ω * eff = r * eff for each thickness d, which is determined from identifying the minimum of the variational free energy including the nonlinear interaction in a self-consistent manner (cf. Fig. 9(a)). We note that all the contributions from the 'spectator' modes, i.e., the modes other than the RTM 0 , are also included in the variational analysis to determine the renormalized phonon frequency. As shown in Fig. 10, to realize the ultimate softening of the RTM 0 mode, a sparser band structure as realized with a thinner layer is preferable since the energy cost due to the nonlinear repulsive interaction between the energetically dense elementary modes can be mitigated.
Finally, towards experimentally testing the present results, we suggest that it is promising to consider paraelectric materials such as SrTiO 3 and KTaO 3 [122,140], as well as SnTe and Pb 1−z Sn z Te [141], where a tuning parameter can be varied by using pressure, chemical composition, and isotope substitution. This flexibility allows one to naturally pre-  10. Softening of the phonon mode relevant to the paraelectricto-ferroelectric phase transition, which is labeled as the RTM0 mode in Fig. 6. The solid colored curves represent the in-plane dispersions of the softened phonon mode with decreasing the insulator thickness d from top to bottom curves. Each dispersion is obtained by setting the phonon-frequency parameter Ω 2 equal to the effective value r * eff determined from the variational analysis for each thickness (cf. Fig. 9(a)). The black dashed line represents the effective phonon frequency in the bulk case. We set c = = β = 1, and the parameters are chosen to be ωp = 5, U = 1, and g = 3.5. We set the momentum cutoff Λc = 10 2 . The tuning parameter is fixed to be r/Λ 2 c = −45. From top to the bottom solid curves, the thickness is varied from d = 4 to d = 1 with step δd = 1.
pare materials at the verge of the ferroelectric phase. Strong nonlinearity is also available in these materials. To be concrete, if we consider SrTiO 3 , one can use the transverse optical phonon mode, for which the phonon nonlinear coupling is U/M 2 2.2 × 10 3 THz 2Å−2 kg −1 [142]. In the unit of c = = β = 1 with, e.g., the temperature T = 60 K, this nonlinear coupling strength can be read as U 1.5, which is close to the value used in the present numerical calculations [143]. We remark that, in this choice, the thickness d = 1 corresponds to an order of several tens of micrometers.

V. SUMMARY AND DISCUSSIONS
Our work demonstrated the great potential of a heterostructure configuration as a platform towards controlling the quantum phase of matter by modifying quantum electrodynamic environments in the absence of external pump. In particular, the present consideration revealed that the strong hybridization between the continuum of matter excitations and cavity electromagnetic fields enriches the spectrum of elementary excitations in a way qualitatively different from the bulk material. More specifically, the main results of this paper can be summarized as follows: (i) Analysis of hybridized collective modes in a cavity.
We analyzed the collective modes in the heterostructure configuration shown in Fig. 2, i.e., the dielectric slab having an infrared active phonon at frequency Ω surrounded by two metallic electrodes. Hybridization of the phonon modes of the dielectric, the cavity electromagnetic modes, and the plasmons of the surrounding metallic mirrors gives rise to a rich structure of the eigenmodes; in particular, we found the emergence of a number of energy modes at frequencies slightly below the bare phonon frequency Ω. The origin of these eigenmodes can be understood from the strong frequency dependence of the dielectric function close to the phonon frequency (see Appendix C).
(ii) Cavity-enhanced ferroelectric phase transition. To go beyond the simple quadratic model of the photonphonon-plasmon hybridization, we included phonon nonlinearities in the dielectric material through the variational method. We found that the collective modes are softened as the thickness of the dielectric slab is decreased. The origin of this softening is the screening of the repulsive component of electric dipoles in the sideby-side configuration by metallic mirrors. The softening leads to the enhancement of the ferroelectric instability, which allows the cavity system to change into ferroelectric even when the bulk material remains paraelectric. We suggest that such cavity-enhanced transition could be experimentally observed in several quantum paraelectric materials at the verge of the ferroelectric phase.
The mechanism proposed in this paper is not restricted to specific setups, but should be applicable to a wide class of sys-tems supporting excitations that possess electric dipole moments. For instance, our theory can be used to analyze the coupling of cavity fields to excitonic molecules [144,145] and also to artificial quantum systems such as ultracold polar molecules [146], trapped ions [147], and Rydberg atoms [148]. The present formulation should also be applicable to analyze polariton modes in the hexagonal boron nitride (hBN) [149][150][151], which is a dielectric insulator material associated with a large optical phonon frequency. We point out again that the predicted changes of the elementary excitations can occur without external driving in contrast to, e.g., the earlier studies of exciton-poalritons in microcavities [152][153][154][155][156].
As by-product of the analyses given in this paper, we envision several potential applications of the cavity-induced changes of the dispersive excitation modes. For instance, one can consider a molecular crystal sandwiched between metallic mirrors, where polaritons originate from the hybridization of singlet excitons and cavity electromagnetic modes. The corresponding triplet excitons have negligible interaction with light and are thus unaffected by the cavity confinement. The heterostructure cavity configuration then allows one to change relative energies of singlet and triplet excitations, which can be applied to improve the efficiency of organic light-emitting devices and realize photon-exciton converters. We expect that the modification of the hybridized dispersion in the proposed cavity setup can also be used to control photochemical reactions in molecular systems. The similar physics may be observed in the presence of a single metal-dielectric interface. Further discussions can be found in Appendix D.
Several open questions remain for future studies. First, we recall that the whole system considered in this paper is assumed to be in thermal equilibrium, i.e., the effective temperature is identical for both of the matter and the cavity. This assumption could be violated if, for example, the metal mirrors and the insulator are further coupled to external bath modes at different temperatures. Moreover, driving the cavity into a high-intensity regime or coherent pumping of atoms should induce the cavity-mediated interactions among the underlying elementary excitations [56,[157][158][159][160]. It remains an intriguing question how our formalism can be generalized to such nonequilibrium setups.
Second, a quantitative amount of the predicted modification of the elementary excitations and also that of the transitionpoint shift can in practice depend on several specifics such as the quality factor of a cavity, the phonon-scattering rate in matter, and the plasma frequency and the conductivity in metals. It merits further study to explore the impact of those experimental complexities on the present results. Meanwhile, on the border of the phase transition, there are also interesting open questions on how quantum critical behavior is modified due to cavity confinement and how it can be affected by dissipative nature of a cavity [122,[161][162][163][164][165][166][167]. It would also be interesting to explore combining cavity-enhanced phase transitions with topological aspects of open systems [168][169][170][171][172][173]. Third, it will be intriguing to consider the situation in which the primary coupling between matter and light is mediated by the magnetic field. A concrete example can be Damon-Eschbach modes in ferromagnetic films [174]. The intrinsi-cally chiral nature of these modes may qualitatively alter the character of collective modes in the cavity configuration. One can also consider a situation with both electric and magnetic couplings being important. This can be relevant to multiferroic systems [175], in which ferroelectric and magnetic orders are intertwined.
Finally, it is also intriguing to consider further features and functionalities that can arise from fabricating additional structures at the metal-insulator interfaces. For instance, periodic texturing or distortion of the surfaces can lead to the formation of gaps in the in-plane dispersion as realized in photonic crystals. This will provide additional flexibility for manipulation of the hybridized light-matter excitations. It also merits further study to consider how changes in the phonon spectrum in the metal-paraelectric heterostructure can affect many-body states in the metallic mirrors, such as, superconducting states. This should be particularly relevant to low-density superconductivity mediated by plasmons as discussed in the case of doped strontium titanate, where the soft transverse optical phonon plays a crucial role [176]. We hope that our work stimulates further studies in these directions. In the linear regime, elementary excitations of the electromagnetic fields in nonabsorbing matter can in general be obtained by solving the macroscopic Maxwell equations with using boundary conditions appropriate for a physical setting of interest. Specifically, the dynamics of the electromagnetic fields is governed by where E is the electric field, B is the magnetic field, µ 0 is the vacuum permeability, and D is the macroscopic electric field in matter. Transforming to the frequency basis, the macroscopic electric field D is related to E via with 0 being the vacuum permittivity and (ω, r) being the frequency-and position-dependent relative permittivity reflecting material properties. We choose the Coulomb gauge In terms of the vector potential, the equation of motion can be simplified as where c = 1/ √ 0 µ 0 is the vacuum light speed. In the bulk case, one can simply solve this eigenvalue problem by employing the spatial invariance and expressing the eigensolutions in terms of the plane waves (cf. Eq. (4)).

Appendix B: Eigenmodes in the cavity configuration
From now on, we focus on the heterostructure configuration as illustrated in Fig. 5, which is symmetric under the spatial parity transformation along the direction perpendicular to the plane, i.e., under the reflection (x, y, z) → (x, y, −z).

Transverse modes
We first consider the transverse modes. In this case, the vector field can be quantized and expanded aŝ where A is the area of the system, η ∈ {q, n, λ} labels each eigenmode in the cavity setting,â η is its annihilation operator, ω η is an eigenfrequency, and U η is the corresponding transverse eigenmode function satisfying ∇ · U ⊥ η = 0. Here, q is the two-dimensional in-plane wavevector, n ∈ N labels a discrete eigenmode in each sector λ = (Λ, P ) with Λ ∈ {TM, TE} specifying the transverse magnetic (TM) or the transverse electric (TE) polarization defined by the conditions, and P = ± specifying the parity symmetry of the eigenmode function under the reflection, Because of the spatial invariance in the x and y directions, we can decompose the eigenfunction as where ρ = (x, y) T . The eigenvalue problem (A7) in the present case can then be given by with the relative permittivity being defined by Here, we use the Heaviside step function and assume the relative permittivity functions in the metals (|z| > d/2) and the insulator (|z| < d/2) as We recall that ω p is the plasma frequency in the metal mirrors, Ω is the frequency of the nondispersive matter excitation such as an optical phonon mode in ionic crystals, and g characterizes the strength of the light-matter coupling between the matter dipolar mode and cavity electromagnetic fields. It is also useful to express the electric and magnetic fields in terms of the eigenmodes as follows:  Fig. 6. The polarization of each mode is specified by Λ ∈ {TM, TE, L}, where TM and TE indicate the transverse magnetic and electric polarizations, respectively, and L denotes the longitudinal polarization. The parity of the mode functions u ⊥ η and f η are specified by P = ± (see also the main text). We also indicate whether or not the divergence of the electric field ∇ · Eη vanishes. The last column shows the explicit functional form of the mode function for each eigenmode, where ei indicates the unit vector in the direction i = x, y, z.

Bands
Λ P ∇ · Eη Mode function TM2, RTMn=0,2,4..., SS (=TM0) The boundary conditions require the continuity of n × E η and n · [ 0 E η ] across the interfaces with n being the unit vector perpendicular to the surfaces. The latter condition is equivalent to the continuity of the magnetic field orthogonal to the surface normal n × B η .

a. TM modes
Without loss of generality, we here consider the TMpolarized eigenmodes satisfying The Maxwell equation ∇ ·B = 0 then leads to q y = 0, i.e., the eigenmodes propagate along the x axis and thus we denote q = q x . The parity conditions (B4) and (B5) in the present case can be read as In terms of the mode function u ⊥ η , the TM modes considered here satisfy We first consider the case of the even parity P = +. Solving the eigenvalue equation (B7) for the mode function u ⊥ η with the transversality condition, i.e., ∇ · (e iq·ρ u ⊥ η ) = 0, we obtain the explicit functional form of u ⊥ η for this class of eigenmodes as shown in the first line of Table II. This set includes the bands labeled as TM 2 , RTM n=0,2,4... , SS (=TM 0 ) that are shown in Fig. 6(a,d). In these modes, the longitudinal wavenumber κ η in the insulator region is defined by We emphasize that κ η takes either a real value or a purely imaginary value. In the latter case, the eigenmodes are localized at the interfaces and represent the surface modes. We define the skin wavenumber ν η , which is real and positive, as follows: The coefficients u M η and u I η in the metals (M) and in the insulator (I), respectively, satisfy the following relations imposed by the boundary conditions: Aside an irrelevant phase factor, these coefficients can be fixed by further imposing the normalization condition dz u ⊥ η (z) 2 = 1.
In the same way as done above, we can also obtain the corresponding eigenmodes in the case of the odd parity P = −. Their functional forms are given in the second line in Table II. The corresponding coefficients satisfy the boundary conditions  Fig. 6 while the in-plane momentum q for low-and high-momentum results is chosen to be q = 0.1 Ω/c and q = 7 Ω/c, respectively. and also the normalization condition dz u ⊥ η (z) 2 = 1.
These modes are labeled as TM 1 , RTM n=1,3,5... , AS 1,2 that are shown in Fig. 6(b,e). We show typical spatial profiles of the electric fields at low and high in-plane momenta q for each of these TM modes in Fig. 11.

b. TE modes
In the case of the TE polarization, without loss of generality we can assume the following conditions: From the transversality condition, and the Maxwell equation ∇ ·D = 0, we obtain the condition q y = 0, i.e., the TE modes considered here also propagate along the x axis. The electromagnetic fields of the TE modes transform under the parity transformation as These conditions can be read in terms of the mode function u ⊥ η as The eigenvalue equation then leads to the functional form in the third (fourth) line in Table II in the case of the even parity P = + (the odd parity P = −). The boundary conditions require that the value of the mode function u ⊥ η and its first spatial derivative must be continuous across the interfaces. More specifically, in the case of the even parity they lead to the conditions on the coefficients u M,I η , while in the odd-parity case we have to impose the conditions u M η /u I η = e νηd/2 sin (κ η d/2) = ν η κ η e νηd/2 cos (κ η d/2) .
The normalization condition dz u ⊥ η (z) 2 = 1 can be also imposed on the coefficients. In Fig. 6(c,f), the TE modes with the even parity are labeled as TE 1 , RTE n=0,2,4... while the ones with the odd parity are labeled as TE 0 , RTE n=1,3,5... . In contrast to some of the TM modes (such as the AS 1,2 and the RTM 0 modes), the spatial profiles of the electric fields are qualitatively insensitive to the in-plane momenta for all of the TE modes; their typical results are shown in Fig. 12.

Longitudinal mode
The other type of solutions is the longitudinal mode that is characterized by the longitudinal electric field Because of our choice of the Coulomb gaugeÂ = 0, the longitudinal field is directly linked with the longitudinal mode of the matter fieldφ as discussed in the main text (cf. Eq. (37)). From the Maxwell equation ∇ ·D = 0, we can see that this mode oscillates in time at a constant frequency Ω corresponding to the vanishing permittivity in the insulator We can expand the longitudinal matter field in terms of the mode function aŝ where the longitudinal mode function satisfies Meanwhile, the Maxwell equation ∇×Ê = − ∂B ∂t ensures the vanishing magnetic fieldB = 0. Thus, it suffices to only consider the longitudinal electric field and, without loss of generality, we here consider the modes satisfying The corresponding in-plane momenta point to the x direction, i.e., q y = 0 and thus q = q x . The boundary conditions at the interfaces lead to the vanishing electromagnetic fields in the metal mirrors. The resulting explicit functional form of the mode function is given in the last line in Table II. Depending on the parity P = ± defined via the relations there are two classes of solutions, which are labeled as the L ± modes in Fig. 6(a). The corresponding longitudinal wavenumber κ η can be simply obtained by the boundary conditions Appendix C: Origin of the resonant soft-phonon modes We here provide simple explanations to elucidate the origin of the resonant softened phonon modes that are denoted by the labels RTM and RTE in Fig. 6. To this end, we consider a frequency regime in which the metallic permittivity is negative p (ω) < 0 while the insulator dielectric function is positive ξ(ω) > 0 (cf. Eqs. (B10) and (B11)). To be specific, we focus on the RTM modes at the zero in-plane momentum q = 0 The mode functions f η of all the TE modes have only the in-plane components, whose spatial profiles are shown by the red solid curve (corresponding to the y component of the electric field).
The parameters are the same as in Fig. 6. We use q = 0.1 Ω/c for the sake of concreteness while the spatial profiles of the TE modes are qualitatively insensitive to a specific choice of the in-plane momentum q.
while similar arguments can be given for the corresponding RTE modes as well.
Inside the insulator region, the RTM modes with q = 0 have the spatial profiles of the magnetic field either sin(κz) or cos(κz), which we denote by the labels P = + or P = −, respectively. We here recall that the parity of the eigenmodes is defined by that of the transverse components of the electric field. In these cases, the longitudinal wavevector κ should satisfy the following equations: ξ(ω) | p (ω)| = cot κd 2 (P = +), − tan κd 2 (P = −). (C1) The Maxwell equations also lead to the additional constraint In the case of the frequency regime slightly below the bare phonon frequency Ω of paraelectric materials, the dielectric function ξ(ω) takes a singularly large value compared with | p (ω)|. Thus, the conditions (C1) and (C2) lead to the approximate eigenmode equation ξ(ω)ω n πc d (n = 1, 2, . . .).
Owing to the singular nature of ξ(ω) in the limit of ω → Ω−0, this condition can in principle be attained for an arbitrary integer n. It is this resonant structure that allows for generating a large number of soft-phonon modes in the heterostructure configuration. In practice, when the damping of the phonon mode is included, there should exist an effective cutoff on possible values of n, but still many eigensolutions should be allowed as long as the damping rate is small enough in comparison with the phonon frequency. Only the singlet state S couples to confined light modes, resulting in the split into the multiple branches of the in-plane dispersions. A lower dispersion (red solid curve denoted as S ) can potentially be applied to enhance the efficiency of organic light-emitting devices. Tuning the in-plane momentum q or the thickness of the layer, a higher dispersion (blue solid curve denoted as S ) can be used to realize the resonance condition ω S = 2ωT, leading to a conversion of a photon into a pair of triplet excitons as denoted by two dashed black arrows.
Appendix D: Applications of cavity-based control of polariton dispersions Substantial modifications of the polariton dispersions in the cavity geometry (see e.g., Fig. 6 in Sec. III) are generic features for collective dipolar modes coupled to confined light modes. In this respect, the cavity-based control of energylevel structures can potentially be applied to a wide range of systems such as ionic crystals, molecular solids, assembled molecules, and excitonic systems. Organic molecules are particularly promising for such applications since their excitations have large dipole moments and crystals can achieve high-molecular densities [177]. Recent experiments have demonstrated the strong coupling between a single molecule and a cavity mode even for the modest cavity quality factors [118,178]. This means that the strength of the coupling between the molecular optical excitation and the cavity mode exceeds the decoherence rate of the excitations [61,82,98,179]. Besides a case study of the cavity-enhanced ferroelectricity in Sec. IV, we discuss below that the cavity-induced changes in the hybridized dispersion can be utilized to improve the efficiency of organic light-emitting devices (LEDs) and to realize exciton-photon converters.
We consider a crystal of polar molecules surrounded by metallic electrodes as shown in Fig. 13. The relevant infrared active mode can be either an electronic excitation, such as the E1 electronic transition, or a vibrational excitation of the molecular crystal. While these two types of excitations have very different frequencies (i.e., optical vs mid-infrared), the general formalism presented in Sec. III should be applicable to both cases. Either one of these modes can be introduced as a nondispersive mode having the frequency Ω in Eq. (16). To be concrete, we focus on the case of electronic excitations, and argue that the proposed cavity configuration can be applied to improve efficiencies in molecular LEDs.
A major challenge in organic LEDs is that the lowest excitonic state is a spin triplet, which is optically dark (i.e., it recombines through nonradiative channels). The reason for this is that Hund's rule for electrons in a molecule favors the spin alignment and thus the spin triplet excitonic states are lower in energy than the singlet ones. Having a practical method of reversing the order of excitonic states by making the singlet state of an exciton lower in energy will have important implications for improving the efficiency of organic LEDs [114]. In the cavity configuration, since triplet excitonic states are dark, they do not couple to light and are not strongly modified. In contrast, the singlet excitonic state hybridizes with confined light modes and forms the softened collective modes such as the RTM 0 mode shown in Fig. 6(a). This should make it possible to have the lowest excitonic state as a spin singlet, thus realizing efficient organic LEDs (see the red solid curve in Fig. 13(b)). While this type of hybridizations have been discussed in the literature before [114], the novelty of the present consideration is the inclusion of the continuum of both collective matter excitations and cavity light modes in a slab of a finite thickness as well as taking into account the dispersive nature of the hybridized excitations there. In general, the dispersive nature could be utilized to further tune the selective control of chemical reactivity landscapes in order to enhance or suppress certain types of photochemical reactions.
We next discuss a complementary idea of using a cavity to facilitate the conversion of photons into pairs of triplet excitons. This process could be useful for generating entangled pairs of excitations. At ambient conditions, such conversion process is challenging to realize because the resonant absorption of light occurs at the energy of singlet excitons and this energy is only slightly above the energy of the triplet excitons; typical energy spacing is characterized by the Hund's splitting and is an order of a few meV.
This limitation can be overcome by utilizing the proposed cavity setup. As discussed before, when a molecular solid is confined in the cavity, the singlet exciton hybridizes with light modes while the triplet one does not. This can lead to the formation of a hybridized singlet mode above the original exciton energy (see e.g., the blue solid curve S' above the black dashed horizontal line S in Fig. 13(b)). One may, for example, use the SS(=TM 0 ) and AS 2 modes for this purpose. Owing to the nontrivial dispersion of such a hybridized mode, one can tune the in-plane momentum q (and the thickness d or the plasma frequency ω p if necessary) such that a hybridized singlet energy ω S becomes equal to twice of the triplet energy 2ω T and thus the conversion process of a single photon into two triplet excitons becomes resonant (cf. the dashed black arrows in Fig. 13(b)). One can create singlet exciton-polaritons at finite in-plane momenta by, e.g., having incident light at nonzero angle or constructing a grating on the surface of the