We introduce a new theoretical approach for analyzing pump and probe experiments in non-linear systems of optical phonons. In our approach, the effect of coherently pumped polaritons is modeled as providing time-periodic modulation of the system parameters. Within this framework, propagation of the probe pulse is described by the Floquet version of Maxwell’s equations and leads to phenomena such as frequency mixing and resonant parametric production of polariton pairs. We analyze light reflection from a slab of insulating material with a strongly excited phonon-polariton mode and obtain analytic expressions for the frequency-dependent reflection coefficient for the probe pulse. Our results are in agreement with recent experiments by Cartella *et al.* [Proc. Natl. Acad. Sci. U. S. A. **115**, 12148 (2018)], which demonstrated light amplification in a resonantly excited SiC insulator. We show that, beyond a critical pumping strength, such systems should exhibit Floquet parametric instability, which corresponds to resonant scattering of pump polaritons into pairs of finite momentum polaritons. We find that the parametric instability should be achievable in SiC using current experimental techniques and discuss its signatures, including the non-analytic frequency dependence of the reflection coefficient and the probe pulse afterglow. We discuss possible applications of the parametric instability phenomenon and suggest that similar types of instabilities can be present in other photoexcited non-linear systems.

## I. INTRODUCTION

### A. Motivation and overview

Control of light by collective excitations in solids, such as phonons, opens exciting possibilities for realizing materials and devices with new functionalities.^{2,3} Recent experiments by Cartella *et al.*^{1} demonstrated terahertz optical amplification in a SiC insulator following a strong mid-IR pump that resonantly excited the SiC stretching mode. Motivated by these experiments, we study the problem of light reflection from a slab of insulating material with a strongly excited polariton mode. We analyze this system using the Floquet formalism in which polariton oscillations provide periodic time modulation of the material properties. We present an analytic solution of the Fresnel light reflection problem in which we include mixing of two frequencies analogous to mixing between signal and idler modes in parametric amplifiers. We demonstrate discontinuous dependence of the reflection coefficient on the incoming photon frequency, which we attribute to the existence of unstable polariton modes. We interpret these instabilities as resonantly enhanced polariton parametric wave mixing. More broadly, the interaction between light and active media has recently been of theoretical interest for its unusual properties. Parametric amplification on time-periodic media is studied intensively.^{4–7} We find that such parametric amplification is possible in pumped SiC media and is mediated by phonon-polaritons.

A simple physical picture of light amplification by a coherently oscillating polariton emerging from coupling of light to phonons is shown in Fig. 1(a). Non-linear coupling between an IR-active phonon *Q* and an electric field *E* leads to the interaction term of the form *Q*^{2} *E*^{2}. When the polariton mode is strongly excited by the pump pulse, we find coherent oscillations of the phonon field ⟨*Q*^{2}(*t*)⟩ = *A* cos^{2}*ω*_{p}*t*, which can be interpreted as Floquet driving of the material. Such driving can result in the production of pairs of photons at frequencies that add up to twice the drive frequency *ω*_{1} + *ω*_{2} = 2*ω*_{p}. This is analogous to four-wave mixing in non-linear optics, with two of the waves corresponding to the collective mode of the material. When there is a probe pulse at frequency *ω*_{s}, we find not only an increase in the number of photons at the same frequency upon reflection but also generation of the complementary idler photons with frequency *ω*_{i} = 2*ω*_{p} − *ω*_{s}. If we neglect feedback of the signal pulse on the coherently oscillating polariton mode, we can analyze the light amplification effect by solving Maxwell’s equations using the two-mode mixing Floquet formalism. In this case, reflection and transmission coefficients are described by the Floquet eigenmodes in which the components of the signal and idler frequencies are mixed. We note that nonlinearity arising from hybridization of a pair of photons with a pair of polaritons is expected to be considerably stronger than the usual optical four-wave mixing since it comes from a resonant process.

In principle, in Floquet systems, one has mixing between more than two frequencies. We discuss generalization of our analysis beyond the two-mode approximation and show that it does not change results appreciably, except at discrete frequencies where higher-order resonances occur for strong pumping. Moreover, the Floquet formalism also allows us to study unstable polariton modes induced by coherent oscillations of the phonon field. We consider the light reflection coefficient as a function of complex frequency and look for its poles in the complex plane. Our analysis treats the weak probe pulse only to linear order; hence, we can apply the standard linear response argument that poles of the response functions correspond to collective modes in the system. Collective modes with a positive imaginary part imply instabilities, corresponding to exponential growth in time of small initial fluctuations, which will ultimately be limited by non-linear terms. In the case of static gain media, the analysis of unstable polariton modes through poles of the reflection coefficients has been discussed in Refs. 8 and 9.

A simple physical picture of the Floquet parametric instability comes from the perspective of resonant polariton wave mixing. Linear coupling between light and phonon *Q* leads to the formation of upper and lower polariton modes, as shown in Fig. 1(b). An initial pump pulse creates a condensate of upper polaritons at zero momentum. Non-linear interactions between phonons and the electric field, such as *Q*^{4}, *Q*^{3}*E*, and *Q*^{2} *E*^{2}, result in scattering of polaritons. In particular, there is a process in which two polariton modes scatter in a way that conserves energy and momentum [see Fig. 1(b)]. As we discuss below, this argument correctly predicts frequencies of the unstable modes, which are naturally connected by the complementarity condition *ω*_{1} + *ω*_{2} = 2*ω*_{p}.^{10} As is clear from this physical picture, large amplitude of phonon-polaritons is crucial for realizing our Floquet parametric instability. Without the pump field, no amplification is observed.

Our results provide a simple physical interpretation of light amplification observed in recent experiments by Cartella *et al.*^{1} Moreover, we argue that Floquet parametric instability could already be present in these experiments since the critical pumping field strength is smaller than the one used by the Hamburg group. Our approach utilizing the Floquet formalism is applicable to a broad class of systems,^{11} and thus, it provides a general framework for investigating instabilities in pumped systems. Finally, our analysis reveals a general mechanism that could be useful for realization of practical sources of THz radiation, relevant for many potential applications.

### B. Relation to earlier work on optical phase conjugation

Before proceeding, it is useful to put our results in the context of earlier works on optical phase conjugation phenomena.^{12,13} While phase conjugation phenomena can arise from a variety of optical nonlinearities, including four-wave mixing, photon echoes, and stimulated Raman scattering,^{14–17} we will focus on a comparison to the Brillouin scattering mechanism.^{13} In this case, one considers a system with coupling of optical fields to excitations of the matter. Examples include phonon excitations and molecular, vibrational, and rotational excitations. The corresponding coupling can be represented as a scattering process,

where *E*_{L} is the pump laser field, *E*_{S} is the stimulated scattered light, and *Q*_{B} describes a matter excitation amplitude, such as a phonon or a molecular vibration. Stimulated scattering leading to optical wave front reversal can be understood as follows:^{13} after the system is pumped by laser light at *ω*_{L}, it can generate a smaller frequency photon *ω*_{S} and excitation Ω_{B}, provided that the condition *ω*_{L} = *ω*_{S} + *ω*_{B} is satisfied. Another way of describing this phenomenon is to observe that after the pump, the system becomes an amplifying medium for *E*_{S} photons. In the context of optical systems, front reversal phenomena correspond to optical frequencies *ω*_{L} and *ω*_{S} being much larger than the matter excitation frequency *ω*_{B}.

The parametric instability of phonon-polaritons that we discuss in our paper does not distinguish between optical fields and matter excitations. Due to strong linear coupling, all relevant fields are collective excitations of the medium: upper and lower polaritons. Parametric instability of a pumped SiC insulator that we describe below corresponds to polaritons pumped to the bottom of the upper polariton branch undergoing stimulated scattering into a pair of polaritons at different energies such that one polariton goes to the lower branch and the other becomes a higher energy excitation in the upper branch (see Fig. 1). We do not have separation of energy scales since all modes are in the mid-IR range of tens of THz. Note that due to strong non-linearities, the relative composition of polaritons as photon/phonon hybrid excitations depends on the amplitude of the driving (see, e.g., Ref. 18). Hence, in order to account for the changing character of the polaritons, we keep track of both the electric field and the phonon amplitudes in our treatment.

One notable difference between our analysis and standard discussion of the optical phase conjugation phenomena is that we focus on non-linear processes taking place near the material boundary rather than in the bulk. Hence, in our case, momentum perpendicular to the interface is not conserved and the system is less constrained by the phase matching condition. This argument also implies that we need to carefully account for evanescent solutions present near the boundary. This requires us to go beyond the standard Slowly Varying Envelope Approximation (SVEA)^{12,19} since at the interface between the medium and air, the amplitude changes on the length scale smaller than the wavelength. This aspect of the problem has a direct analog in the scattering theory of electrons in condensed matter physics. When describing electron states in metals, one is usually interested in single electron states with energies close to the Fermi energy. Hence, one can often represent electronic wave functions as

where *A*(*x*) is a slowly varying amplitude (for simplicity, we discuss the 1D case). This approximation works well when external parameters vary smoothly on the scale of the Fermi momentum *λ*_{F}. Approximation (2) reduces the Schrödinger equation for the electron to a first-order differential equation in *x*, and as the boundary conditions, one requires the continuity of *ψ*(*x*) but imposes no conditions on $d\psi dx$ since $d\psi dx\u2248kF(x)\psi (x)$. Form (2) is insufficient when one considers a problem of electron transport at the interface of two materials that have different values of *k*_{F}. The mismatch of *k*_{F} gives rise to scattering of electrons,^{20} which one misses when imposing a single boundary condition of the continuity of *ψ*(*x*).

In our discussion of the pump and probe experiments in SiC, we compute the reflection coefficient of the probe beam using the full solution of Maxwell’s equations and not resorting to the SVEA. In the latter approximation when looking at light propagation in a non-uniform system, one only imposes the continuity of one of the fields, for example, *E*, and assumes that continuity of *B* is satisfied automatically since it can be related to the spatial derivative of *E* using the dielectric function. This is accurate only when parameters vary smoothly on the scale of a wavelength. At sharp interfaces, material properties change abruptly on the scale shorter than the wavelength; hence, the full solution of Maxwell’s equations is necessary. One of the most surprising consequences of going beyond the SVEA is the discontinuities in the frequency dependence of the reflection coefficient (see Fig. 5).

To capture these instabilities, we treat the material as a Floquet medium whose optical and elastic properties are modulated periodically in time by the strongly excited polariton mode. We truncate the resulting Floquet equations, taking into account the hybridization between resonant processes non-perturbatively while neglecting non-resonant processes. The difference of this approach from conventional non-linear optics approaches^{21,22} is analogous to doing degenerate perturbation theory in quantum mechanics: mixing between degenerate states should be treated exactly, while non-degenerate processes can be neglected. Hybridization, then, leads to a reorganization of normal modes, which reflects the mutual feedback of signal and idler modes that are brought into resonance by the periodic modulation of the medium. Then, when the frequencies of the signal and idler simultaneously match those of natural excitations of the system [e.g., the finite momentum upper and lower polariton pair shown in Fig. 1(b)], this leads to a parametric instability in which these modes get spontaneously amplified. Hence, a key element of our approach is the exact treatment of resonant processes.

### C. Organization of this paper

This paper is organized as follows. In Sec. II, we introduce both linear and non-linear models of light–matter coupling in SiC. We first solve the linear model and discuss the upper and lower polariton excitations. We then consider the non-linear model and discuss propagation of the probe pulse following strong pumping resonant with the bottom of the upper polariton branch. In Sec. III, we solve the Fresnel reflection problem of the medium with the large amplitude polariton mode. We obtain an analytic solution and find that the reflectivity has a discontinuous dependence on the probe light frequency. Furthermore, we numerically show that the parametric instability coincides with these discontinuities. Up to this point, we perform full analysis of our problem with respect to each frequency mode of the probe light. In Sec. IV, we apply our result to analyzing temporally dependent probe light. We study the optical properties arising from the parametric instability. By analyzing the temporal problem of the propagation of the probe light, we show that the parametric instability induces the probe pulse afterglow. We also show that for the pump field strength and the angle of incidence in the current experiments, the parametric instability should already be present. In Sec. V, we discuss possible applications of the Floquet parametric instability phenomenon and future directions.

## II. EFFECTIVE MODEL AND POLARITONS

### A. Linear theory and polaritons

Polaritons in solids can be understood as collective excitations arising from hybridization of IR-active phonons with light. While this phenomenon is a standard textbook material (see, e.g., Ref. 23), we present a brief derivation here to establish notations for subsequent discussions. We limit our discussion to the transverse optical (TO) polaritons in which electric field and polarization are perpendicular to the wave vector.

In the linear response regime, the interaction between an IR-active phonon mode and the electric field is described by the effective potential energy density

where **Q** is the phonon amplitude, **E** is the electric field, *Q* ≡ |**Q**|, *E* ≡ |**E**|, Ω_{TO} is the phonon frequency, *ϵ*_{0} is the permittivity of vacuum, and *ϵ*_{∞} is the dielectric response at high frequency. The parameter *η* characterizes the electric dipole moment of the phonon mode. Equation (3) can be used to obtain the medium polarization density,

and the equation of motion for the phonon amplitude,

where *γ* is the damping coefficient of the phonon. Collective modes can be found by combining the last two equations with the Maxwell equations,

It is convenient to combine the two Maxwell equations by taking a curl of (7), using (8), and substituting **P** from (5),

Equations (6) and (9) provide two coupled equations describing propagation of small amplitude electromagnetic waves inside the insulator. In a bulk material, we take solutions to be plane waves,

where *ω*_{k} is the collective mode frequency and **k** is its wave vector. We assume that **Q**_{p} is parallel to **E**_{p}, and both are perpendicular to **k**. We find a linear set of equations for the Fourier amplitudes,

where *k* ≡ |**k**| and $g(\omega )\u2261\u2212\omega 2\u2212i\omega \gamma +\Omega TO2$. To obtain eigenmodes, we require the determinant to be zero and obtain

Solutions of the last equation describe lower and upper polaritons shown in Fig. 1. At *k* = 0 and assuming *η* ≠ 0, the frequency of the lower polariton is at *ω* = 0; that of the upper polariton is given by

For the upper polariton at *k* = 0, we find a simple relation between the field and phonon amplitudes,

### B. Nonlinear interactions

When the amplitude of phonon oscillations and electromagnetic fields becomes large, we need to include non-linear terms in the phonon–light interaction. We include terms up to fourth order in potential (3), consistent with the symmetry of the system,

We point out that Eq. (16) is valid when an incident electric field is polarized in the ab plane of the 4H-SiC, as was the case in experiments by Cartella *et al.*^{1} The space group of 4H-SiC is P63mc (aka “Wurtzite crystal structure”). The corresponding point group is 6 mm, which has six mirror planes,^{24} all of them parallel to the *c* axis (i.e., they are all “vertical” mirror planes). Analysis presented in Ref. 1 shows that the *k* = 0 phonon participating in the formation of the polariton relevant to their experiments should have polarization in the ab plane and transform according to the two-dimensional irreducible representation of 6 mm called E2.^{25} The square of E2 contains the trivial representation; hence, the term *Q*^{2} can be part of the Hamiltonian. On the other hand, the cube of E2 does not contain the trivial representation (to see this, one can raise the characters of the E2 row in the table to the third power and observe that this gives a set of characters that has zero overlap with the trivial representation). In fact, one can easily check that only even powers of E2 contain the trivial representation. We note that our analysis will also apply to the phonon that obeys the two-dimensional irreducible representation E1 that transforms as the vectors x and y. In writing (16), we have omitted terms of the form *E*^{3}*Q* and *Q*^{4}, as their coefficients are known to be small in SiC.^{1} Note also the absence of *E*^{4}, which allows us to preserve the standard form of the Maxwell equations. Nonlinearities are absorbed into the electric field dependence of the polarization. Note that the phonon–light interaction presented here is specifically tailored using the symmetry considerations to describe SiC. However, our Floquet approach is not restricted to this type of symmetry. Different symmetry allows different terms, e.g., cubic terms, which leads to instabilities at different frequencies. We show the analysis with more terms in Sec. IV B.

After including nonlinear terms, we find modified expressions for polarization and the equation of motion of the phonon,

Nonlinear polaritons can now be obtained by solving the phonon equation [Eq. (19)] together with the Maxwell equation,

One important consequence of including nonlinear terms is the amplitude dependence of the polariton frequency. Another consequence is frequency mixing that we discuss in Subsection II C.

As the first step in analyzing the consequences of non-linearities, we neglect higher harmonics but consider renormalization of the polariton frequency by the finite excitation amplitude. From now on, when we refer to the upper polariton frequency at *k* = 0 as *ω*_{p}, we imply that its renormalization by the finite drive amplitude has been included. The ratio between **E**_{p} and **Q**_{p} is also modified by the non-linear terms, and we denote it as *ν*,

In the context of pump and probe experiments, we are interested in the situation where an upper polariton is excited by a strong and short pump pulse, which is then followed by a much weaker probe pulse. For simplicity, we assume that the pump excites the upper polariton at the frequency *ω*_{p} only. Realistically, a finite range of frequencies surrounding *ω*_{p} is excited. In analyzing reflection of the probe pulse, we will neglect its feedback on the pumped polariton. Furthermore, we will assume that on timescales of duration of the probe pulse, we can neglect the decay of the pumped polariton. Hence, everywhere in the analysis below, we assume that the pumped polariton is characterized by time-independent parameters **Q**_{p}, *ω*_{p}, and *ν*. Note that the dependence of *ω*_{p} and *ν* on the polariton amplitude **Q**_{p} is always included in our numerical calculations.

### C. Frequency mixing in the bulk

The general setting of pump and probe experiments has a strong pump pulse followed by a weaker probe. In this spirit, we consider a problem of a small amplitude electro-magnetic pulse propagating in the background of a large amplitude polariton oscillation. Our philosophy is to treat polariton oscillations at *ω*_{p} as providing periodic time modulation of the medium properties, which can be interpreted as Floquet driving. The Floquet modulation frequency is 2*ω*_{p} since nonlinear corrections to the linear equations of motion are always of the form *Q*^{2}**E** or *E*^{2}**Q**. The Floquet driving provides mixing between different frequencies, and in general, one needs to include an infinite number of harmonics. We limit our discussion to the lowest order mixing and consider small amplitude oscillations at two frequencies only, which satisfy the relation *ω*_{s} + *ω*_{i} = 2*ω*_{p}. In the spirit of the general formalism of parametric resonance, we will refer to the components of the frequencies *ω*_{s} and *ω*_{i}(≡2*ω*_{p} − *ω*_{s}) as the signal and the idler, respectively. The approximation of only including mixing between two modes is justified when nonlinearities are not too strong. As a consistency check, we generalized the analysis presented in this paper to include more modes and verified that including them does not change any results qualitatively. We discuss the case of mixing between three modes in Sec. IV B.

We introduce three frequency components for the electric field and the phonon amplitude as

where, for concreteness, we assume a semi-infinite medium in the *z* < 0 half-plane. Here, **E**_{p} and **Q**_{p} are large fields, created by the strong pump pulse that drives the polariton, whereas **E**_{κ} and **Q**_{κ} are small fields, created by the weak probe pulse. The latter fields are of the form

so that **E**_{1} and **E**_{2} should be interpreted as amplitudes of the signal and idler components, respectively. Our goal is to solve Maxwell’s equations linearized with respect to the small *ω*_{s} and *ω*_{i} components. This allows us to neglect their feedback on the strong polariton component at *ω*_{p}. In writing Eqs. (22) and (23), we changed the notation for the spatial dependence of the fields to the evanescent wave form *e*^{κz}, where *κ* ≡ *ik*. This is done for convenience of subsequent discussions. Note that the frequency components *ω*_{s} and *ω*_{i} have the same spatial structure because they correspond to different Fourier components of the same Floquet eigenmode and the pumped polariton is at *k* = 0. Having the same spatial structure for different frequency components should be contrasted to solutions of the Maxwell equations in stationary media, where each frequency *ω* has its own wave vector *κ*(*ω*). In the Floquet formalism, a single eigenmode is constructed as a superposition of several frequency components, all sharing the same wave vector.

After straightforward calculations, which are presented in Appendix A, we obtain linear equations for the imaginary wavevectors *κ* of Floquet eigenmodes,

where

Equations (27)–(32) are arranged to provide a simpler physical interpretation. The function $\u03f5\u0303(\omega )$ has been defined as an effective dielectric constant when solving the linearized phonon-Maxwell equations in the presence of strong oscillations at *ω*_{p} but without including the frequency mixing between the signal and idler Fourier components. This dielectric constant function is then used to define the “imaginary wave vectors” $\kappa \u0303s$ and $\kappa \u0303i$ corresponding to *ω*_{s} and *ω*_{i}. Figure 2(a) shows the dependence of $Im[\kappa \u0303s]$ and $Im[\kappa \u0303i]$ on *ω*_{s} (in the latter case, *ω*_{i} = 2*ω*_{p} − *ω*_{s}). Material parameters used in our analysis have been taken from Ref. 1 and are summarized in Table I. The function $Im[\kappa \u0303s(\omega s)]$ is qualitatively similar to the dispersion for linear polaritons. We can identify the lower and upper polariton branches below and above the bare phonon frequency Ω_{TO}.

Ω_{TO} | 2π · 24.0 (THz) |

γ | 0.2 (THz) |

ϵ_{∞} | 5.91 |

E_{p} | 8.00 × 10^{2} (MV/m) |

η | 8.39 × 10^{8} ($As/kg12m32$) |

α | 7.11 × 10^{25} ($As\u2009/kg32m12$) |

β | 1.69 × 10^{7} (A^{2}s^{4}/kg^{2}m^{2}) |

ν | 9.82 × 10^{11} (MV/$kg12m12$) |

Ω_{TO} | 2π · 24.0 (THz) |

γ | 0.2 (THz) |

ϵ_{∞} | 5.91 |

E_{p} | 8.00 × 10^{2} (MV/m) |

η | 8.39 × 10^{8} ($As/kg12m32$) |

α | 7.11 × 10^{25} ($As\u2009/kg32m12$) |

β | 1.69 × 10^{7} (A^{2}s^{4}/kg^{2}m^{2}) |

ν | 9.82 × 10^{11} (MV/$kg12m12$) |

From the solution of (26), we find the Floquet eigenmodes

The eigensolutions of $(E1,E2*)$ satisfy

where $(E1+,E2+*)$ and $(E1\u2212,E2\u2212*)$ are the eigensolutions of *κ*_{+} and *κ*_{−}, respectively, and

To summarize the discussion in this section, the Floquet eigensolutions are expressed as superpositions of two waves of frequencies *ω*_{s} and *ω*_{i},

The results of numerical calculations of *κ*_{±} are shown in Figs. 2(c) and 2(d). Note that to compute the renormalization of *ω*_{p} due to the finite amplitude of the pumped polariton, we need to solve $Re[\u03f5\u0303(\omega p)]=0$. For example, we find *ω*_{p} = 29.4 THz when |**E**_{p}| = 8 MV/cm. This should be compared to the bare polariton frequency of 30.5 THz.

### D. Discussion of Floquet eigenmodes

#### 1. Discontinuities and phase matching condition

The most surprising feature of the computed eigenmodes, shown in Fig. 2(c), is the discontinuous dependence of Im[*κ*_{α}] at the frequencies 20, 30, and 39 THz (here, *α* is either + or −). From a mathematical point of view, these discontinuities arise because Eq. (33) is an equation for $\kappa \alpha 2$, and *κ*_{α} itself has an overall sign ambiguity. For a semi-infinite sample, we fix the sign using the condition Re[*κ*_{α}] ≥ 0 so that in the light reflection problem, we can construct solutions decaying from the interface into the bulk of the material. In Fig. 3, we show that when $\kappa \alpha 2$ crosses the negative real axis, Im[*κ*_{α}] jumps discontinuously, whereas Re[*κ*_{α}] has a kink as it touches zero while always remaining non-negative. Figure 2(d) shows $1Re[\kappa \u0303\alpha ]$, which can be understood as the localization length of the Floquet eigenmodes. In agreement with our discussion, the points where $Im[\kappa \u0303\alpha ]$ is discontinuous coincide with divergences in $1/Re[\kappa \u0303\alpha ]$, implying that the Floquet solutions at these frequencies are perfectly traveling waves with no decay. As we will show in the following, for a semi-infinite medium, the discontinuities in $Im[\kappa \u0303\alpha ]$ result in discontinuous jumps in the reflection coefficient *R*_{s}(*ω*).

For a finite slab geometry, one may object that fixing the sign of *κ*_{α} is not physical since solutions with both signs (corresponding to both exponential decay and growth as a function of *z*) must be kept in the reflection problem. However, as the slab becomes thicker, the coefficient of the growing solution becomes exponentially small, and for a thick enough slab, it plays a negligible role in the reflection problem. In other words, we expect that for a thick sample, the discontinuities in *R*_{s}(*ω*) are replaced by rapid crossovers near the frequencies of the discontinuities in Im[*κ*_{α}]. The crossover regions are given the frequency windows over which the localization lengths in Fig. 2(d) exceed the sample thickness. The widths of these windows decay as the inverse slab thickness. We verified that this is, indeed, the case through explicit computation of reflection from finite slab geometries.

Let us now discuss the physical conditions required for the discontinuities in Im[*κ*_{α}] to occur. A comparison of Figs. 2(a) and 2(c) shows that these occur near frequencies where

Comparing with Eq. (33), we observe that when condition (37) is satisfied, Floquet frequency mixing set by the coefficients Γ_{1} and Γ_{2} becomes strong and significantly shifts the mode vectors *κ*_{±}, which become very different from *κ*_{{s,i}}. Condition (37) can be understood as the phase matching condition: the two frequency components have the same spatial structure and can be efficiently mixed by the spatially uniform but time varying modulation of system parameters at frequency *ω*_{p}. We also observe that the discontinuities occur in the frequency regime where $|Im[\kappa \u0303{s,i}]|\u226b|Re[\kappa \u0303{s,i}]|$ so that the signal and idler modes are primarily traveling waves with relatively weak damping.

In addition to the discontinuities discussed above, note that Im[*κ*_{α}] in Fig. 2(c) also displays sharp crossovers accompanied by peaks at the frequencies 24 and 35 THz. The crossover at 24 THz occurs at the original phonon mode frequency Ω_{TO}, whereas that at 35 THz is an “idler satellite” peak at frequency 2*ω*_{p} − Ω_{TO}. These crossovers occur because $1g\u0303(\omega )$ has a sharp maximum at Ω_{TO}. The phonon dissipation *γ* controls the sharpness of these crossovers and their width.

#### 2. Eigenmodes away from frequency matching condition

Far from the special frequencies satisfying the phase matching condition, we have

and hybridization between the signal, *ω*_{s}, and idler, *ω*_{i}, modes is weak. We find that for $|\kappa \u0303s|<|\kappa \u0303i|$,

whereas for $|\kappa \u0303s|>|\kappa \u0303i|$, the roles of + and − are reversed.

## III. FRESNEL–FLOQUET REFLECTION PROBLEM

### A. Solution of Fresnel–Floquet problem

We now consider the problem of light reflection from an insulator with an excited polariton mode. Here, we present the analysis for a semi-infinite slab of material and normal angle of incidence. Generalizations to cases of a slab of finite thickness and arbitrary angle of incidence are discussed in Appendix C. The main difference between our system and the canonical Fresnel reflection problem is frequency mixing. An incident beam at frequency *ω*_{s} excites an electromagnetic response in the insulator at two frequencies: *ω*_{s} and *ω*_{i} = 2*ω*_{p} − *ω*_{s}. Therefore, for both the transmitted and reflected beams, we need to include both frequency components (see Fig. 4). When describing light penetrating into the nonequilibrium insulator, we need to use eigenmodes of the system. Hence, we will decompose the transmitted wave into Floquet eigenmodes *κ*_{+} and *κ*_{−}. We then have the following decomposition of the fields:

In the air, *z* > 0, for linearly polarized light,

where

Here,

Note that we take the opposite signs for *κ*_{s} and *κ*_{i} in Eqs. (45) and (46) so that both describe waves traveling upward, i.e., reflected from the surface of the insulator. The complex coefficient *r*_{s} appearing in Eq. (44) is the reflection coefficient for the signal, whereas *r*_{i} describes the amplitude of the reflected idler, both of which are determined below.

In the insulator, for *z* < 0,

Expressions for *E*_{±}(*t*) should be taken from (36). By substituting *E*_{±}(*t*) into (7), we find that the magnetic field in (47) is given by

The coefficients *r*_{s/i}, *t*_{±} can be determined using the boundary conditions for the electric and magnetic fields. We require continuity of the electric and magnetic fields at the boundary (note that, for normal incidence, the fields are parallel to the interface), and we find

### B. Numerical results for the reflectivity *r*_{s}

The results of numerical calculations of the reflectivity

for three values of the pump field amplitude, |**E**_{p}| = 0, 5, and 8 MV/cm, are represented in Fig. 5. Due to the non-linearities, the frequency of the polariton mode differs somewhat in the three cases, taking values *ω*_{p} = 30.5, 29.9, and 29.4 THz, respectively. System parameters are taken from Table I. The green line corresponds to the situation without the pump field, i.e., |**E**_{p}| = 0. Note that in this case, the incident light is reflected almost perfectly in the frequency gap between the two branches of polaritons. This spectral region Ω_{TO} < *ω*_{s} < *ω*_{p} is called the reststrahlen band. Orange and blue lines correspond to *R*_{s} for |**E**_{p}| = 5 and 8 MV/cm, respectively. For finite |*E*_{p}|, the reststrahlen band is still visible, although its width is reduced due to the renormalization of *ω*_{p} by non-linear effects. At finite pump fields, *R*_{s} starts to show non-analytic dependence on frequency. At frequencies *ω*_{s} = 20, 29, 39 THz, we observe abrupt jumps in the reflection coefficient. At *ω*_{s} = 35 THz, we observe a peak in the reflection coefficient but no discontinuity [see Figs. 7 and 8 for a more detailed information about the frequency dependence of *R*_{s}(*ω*) around these points]. The jumps and the peak of *R*_{s} at these frequencies are obviously related to the discontinuities and the crossover in *κ*_{±}, which, in turn, originates from the fact that phase matching conditions are satisfied at these frequencies. We point out that for different pump intensities, singularities at 35 and 39 THz occur at slightly different frequencies because *ω*_{i} depends on the value of *ω*_{p}.

Far from the special frequencies that satisfy the phase matching condition (37), *r*_{s} is essentially the same as in the linear model regardless of the value of |**E**_{p}|. In this case, we can use (39) and approximate (50) as

This is the familiar Fresnel’s formula. Since Im[*κ*_{s}] and $Im[\kappa \u0303s]$ have the same sign and Re[*κ*_{s}] = 0, Eq. (55) shows that in this case, *R*_{s} is continuous and smaller than one.

We note that the value of the pump electric field used in experiments by Cartella *et al.* was as high as 8.68 MV/cm. Hence, we expect the Floquet parametric instability to be easily achievable in experiments with a SiC insulator. Since our analysis is based on a simplifying assumption that the pump pulse is monochromatic at the frequency of the *k* = 0 upper polariton, a detailed comparison between Fig. 5 and the experimental result in Ref. 1 is not possible. Rather, we aim to propose that the novel instability is present at the frequencies which Cartella *et al.* did not focus on.

For a direct quantitative comparison with the experimental results, we need to take into account the following effects. In reality, the finite duration of the pump pulse broadens its spectrum such that a finite window of polariton energies is initially excited. This is expected to lead to broadening of the non-analyticities of the reflection coefficient that we computed. In addition, the pulse used in Ref. 1 for probing was sufficiently broadband so that it is not possible to separate signal and idler components by frequency. The directions of outgoing signal and idler components are different; thus, to understand what fraction of the idler is being captured by the detectors, one needs to know the range of angles of the incident probe pulse and the range of angles captured by the detector. Furthermore, under realistic experimental conditions, two effects are expected to limit the non-analytic frequency dependence of the reflection coefficient and turn discontinuities into sharp crossovers. First, the pump pulse has finite penetration depth into the material (see, e.g., Ref. 26 for a discussion of the effect of Floquet driving with a finite penetration depth in superconductors, although we emphasize that in insulators, such as SiC, the penetration depth is expected to be much larger). Hence, the reflection of the probe pulse will have the character of reflection from a finite slab discussed in Appendix B. Second, at frequencies close to the parametric instability, the reflection coefficient becomes large and non-linear terms for the probe pulse should be taken into account. These are expected to suppress the reflection coefficient, as was recently discussed for Josephson microwave parametric amplifiers by Sivak *et al.*^{27} This latter effect can be neglected for weak enough probe pulses.

### C. Understanding *r*_{s}: Non-analytic behavior and lasing modes

The most surprising features of the results presented in Fig. 5 are the discontinuities in the frequency dependence of the reflectivity *R*_{s}. There are two close pairs of discontinuities, with the first pair at 19/20.5 THz and the second pair at 38/39.5 THz. There is also an isolated spike at frequency 29 THz. These discontinuities are a direct consequence of the discontinuous dependence of *κ*_{±} on *ω* discussed in Sec. II D. Previously, we related such discontinuities to the existence of special frequencies that satisfy phase matching conditions for the signal and idler modes, at which Floquet modulation gives rise to strong mixing between the modes. We will now discuss that behind these discontinuities lies an even more dramatic physical effect: parametric instabilities of the Floquet medium. In the case of a close pair of discontinuities, e.g., the 38/39.5 THz pair, unstable modes appear in the whole range of frequencies between the two. In the case of a single spike, parametric instability occurs at the same frequency. As we discuss below, these modes are spatially extended inside the slab.

We now review the arguments that allow us to identify parametric instability as the origin of singular behavior of *R*_{s}. Let us consider reflectivity from the perspective of a response function of a physical system. It is useful to extend the definition of reflection coefficient from real to complex frequencies. We can look for poles of *R*_{s/i}(*ω*) and identify them with collective modes of the system. In a static dissipative medium, such poles should have a negative imaginary part of the frequency, indicating a finite decay rate. For gain media, one can find poles with a positive imaginary part, which correspond to parametric instability instabilities.^{9} The relation of these poles and the lasing phenomena are presented in Sec. IV.

In Fig. 6, we present results of a more detailed analysis of the reflectivity in a slab of finite width of pumped SiC at complex frequencies close to the frequency *ω*_{s} = 39 THz. This corresponds to one of the frequencies where we observed non-analytic behavior of *R*_{s} in the case of a semi-infinite slab. We observe that discrete poles are distributed on the circumference of a circle, which extends into both upper and lower half-planes. The circumference passes near the curves where Re[*κ*] = 0 for either *κ*_{+} or *κ*_{−}, and in the limit of a semi-infinite sample, the circumference coincides exactly with these curves. The poles on the upper-half of the circumference have positive imaginary frequency, i.e., they correspond to the lasing modes. Discrete character of eigenmodes comes from using a slab of finite thickness, which gives an analog of discrete modes in a cavity. As we take the slab thickness to infinity, the poles become denser and their residues are reduced. In the limit of a semi-infinite medium, discrete poles merge into a continuous line of poles with finite residue density per length. Then, the reflectivity *R*_{s} develops a step-like discontinuity across the line of poles.

Plots of *R*_{s} for an infinitely thick slab of SiC as a function of complex frequency *ω*_{s} are represented in Fig. 7. We also see that some poles are located near the real axis. However, they are not important in our analysis because closer examination shows that the line that satisfies Re[*κ*] = 0 near the real axis is in the lower half-plane, indicating that the poles around this line have a finite lifetime.

We observe a similar pattern of discontinuities in the imaginary part of *R*_{s} around frequencies *ω*_{s} = 20 THz. This reflects the symmetry between *ω*_{s} and 2*ω*_{p} − *ω*_{s} in *κ*_{±}. On the other hand, at *ω*_{s} = *ω*_{p}, we find poles not on a circle but on a straight line [see Fig. 7(c)]. This can be understood by observing that symmetry *κ*_{+} = *κ*_{−} holds at this frequency.

### D. Enhanced reflection around 24 and 35 THz

Examination of the reflection coefficient in the case without pumping (green line in Fig. 5) shows that it approaches one for frequencies between 24 and 30 THz. This frequency range corresponds to the reststrahlen band of SiC for which we do not find bulk polariton modes. In systems with pumping, we observe that in the same frequency range, Re[*κ*_{+}] is large [see Fig. 2(d)]. This indicates that light incident at this frequency cannot penetrate deeply into the material. This is the remnant of the reststrahlen band in the static case. We also observe that Re[*κ*_{+}] goes through the maximum at 24 THz. This is the frequency of the original phonon mode and corresponds to the strongest losses inside the medium. Due to the Floquet nature of the system, Re[*κ*_{+}] is also large for the “Floquet mirror frequencies” 24 THz $<\omega s<$ 35 THz. The largest value of Re[*κ*_{+}] occurs close to 35 THz, which is a Floquet partner of the phonon frequency of 24 THz. Large values of Re[*κ*_{+}] suggest that light cannot penetrate the Floquet medium, and we should find an increase in the reflection coefficient. This is, indeed, the case, as can be seen from Fig. 5.

It is interesting to point out that extending the reflection coefficient to the complex *ω*_{s}, we find isolated poles in the lower half-plane with Re[*ω*_{s}] close to 24 and 35 THz and an imaginary part slightly below −*γ*/2, where *γ* is the phonon dumping coefficient, as shown in Fig. 8. However, we do not expect these poles to have any consequences for the reflection coefficient at real frequencies. These poles are separated from the real axis by the branch cut in the reflection coefficient at Im[*ω*_{s}] = −*γ*/2. Even in the static case, one finds an isolated pole in the reflection coefficient in the lower half-plane with the real part of frequency matching the original phonon frequency and imaginary part slightly below −*γ*/2. In the static case, one also has a branch cut in the reflection coefficient at Im[*ω*_{s}] = −*γ*/2 that separates these poles from the physical real axis. In Table II, we summarize the collective modes discussed in Sec. II C and D.

20 THz | Lasing mode |

24 THz | Re[κ_{±}] is large, strong reflection |

29 THz | Lasing mode |

35 THz | Re[κ_{±}] is large, strong reflection |

39 THz | Lasing mode |

20 THz | Lasing mode |

24 THz | Re[κ_{±}] is large, strong reflection |

29 THz | Lasing mode |

35 THz | Re[κ_{±}] is large, strong reflection |

39 THz | Lasing mode |

## IV. DISCUSSION OF THE PARAMETRIC INSTABILITY

### A. Lasing, causality, and afterglow

The existence of lasing modes modifies how one needs to analyze reflection of the probe pulse in the time domain. The standard approach is to decompose the incoming pulse in the Fourier series as

Then, one uses scattering amplitude for individual Fourier components of *r*_{s}(*ω*) and *r*_{i}(*ω*) to construct the reflected wave as

As long as the poles of the reflection coefficient $rs,\u2009i(\omega )$ are in the lower half-plane, one can show that this procedure satisfies causality: reflected light does not appear before the incident pulse approaches the surface of the material. One can verify this by taking a pulse that only hits the surface at *t* = 0, i.e., for *t* < 0, we have *E*_{0}(*z* = 0, *t*) = 0 and the function *e*_{0}(*z* = 0, *ω*) is analytic in the upper half-plane of complex *ω*. When computing the integral in (57) for *t* < 0, one can close the contour of *ω* integration in the upper half-plane and prove that *E*_{r}(*z* = 0, *t* < 0) = 0. When $rs,\u2009i(\omega )$ has poles in the upper half-plane, this argument does not work and Eq. (57) formally implies nonzero amplitude of the reflected light before the incident beam interacts with the material. To avoid this contradiction, one needs to define the range of frequency integration to be a line in the complex *ω* plane parallel to the real axis but going above any poles of *r*_{{s,i}}(*ω*) (see Fig. 9). By deforming this contour as shown in Fig. 9, we see that this is equivalent to taking a Fourier integral on the real axis and adding contributions from the poles of *r*_{{s,i}}(*ω*). Note that poles of *r*_{s}(*ω*) and *r*_{i}(*ω*) always coincide because the denominators of *r*_{s}(*ω*) and *r*_{i}(*ω*) are the same in (50) and (51). Then, the temporal solution of the electric field is given by

Taking $\omega j=\omega j\u2032+i\omega j\u2033$, we find that for $\omega j\u2033>0$, the contribution coming from pole *j* grows exponentially in time as $e\omega j\u2033t$. At long times, exponential growth of the reflected light should be limited by processes that we did not include in our analysis, such as non-linear interaction of the signal and idler modes and their feedback depleting the pump field. However, we expect that our approach correctly predicts the characteristic time of the “afterglow” that develops following the probe pulse after time

We note that even in the absence of the probe pulse, quantum and thermal noise may be sufficient to seed the parametric instability. The “afterglow” should have Fourier components that match the frequency of unstable modes [see the $ei\omega j\u2032t$ part of the second line of Eq. (58)]. Hence, the “afterglow” signal should carry direct signatures of unstable modes. A similar residual response after a pulse is known as ringing.^{28} Ringing is the time response followed by an overshoot, resulting in decaying oscillations of increasing frequency.^{29} When we assume the phonon decays rapidly, the “afterglow” may exhibit a similar behavior to ringing. However, we also expect that the exponential growth of the reflected light will be observed in the presence of the strong pumped phonon-polariton.

Before concluding this section, we would like to address the question of what makes frequencies of 20, 30, and 39 THz special so that parametric instabilities take place at these specific frequencies. Physical mechanism responsible for parametric instabilities is shown in Fig. 1. Pumping creates a condensate of polaritons at the bottom of the upper branch. Non-linearities present in the system allow a pair of polaritons from the condensate into a pair of excitations with opposite momenta so that the total energy and momentum of the pair are conserved. Energies 20 and 39 THz correspond to exactly such scattering where one of the excitations stays in the upper branch and the other goes to the lower branch. In addition, the 30 THz instability corresponds to a scattering of a pair of polaritons with a very small change of momenta and energy of individual excitations. A similar mechanism of coherent scattering of bosonic excitations from a photoexcited coherent state has been observed earlier in experiments with exciton-polaritons in microcavities.^{30–32}

### B. Lasing threshold and nature of unstable mode

One interesting question related to lasing modes is whether they are present for infinitesimal amplitude of the pumped polariton amplitude or appear beyond a certain threshold. In Fig. 10, we show the maximum value Im[*ω*_{s}] of the unstable modes for different pump field strength *E*_{p}. In equilibrium, i.e., when *E*_{p} = 0, all collective modes are in the lower half-plane. The first unstable mode appears for *E*_{p} = 1 MV/cm. Since solutions for *κ*_{±} are symmetric with respect to *ω*_{s} and *ω*_{i}, unstable modes around frequencies *ω*_{s} = 20 and 39 THz appear simultaneously. Finite threshold for lasing in our system can be understood as a result of competition between parametric driving and losses in the medium. Increasing *E*_{p} beyond the threshold value, we find that the maximum value of Im[*ω*_{s}] increases linearly with amplitude of the pumped polariton. Another family of unstable modes that we discussed before is at frequency *ω*_{p}. We find that the growth rate of this instability, i.e., the corresponding maximum value of Im[*ω*_{s}], is always smaller than the growth rate of the modes at *ω*_{s} = 20 and 39 THz. Note that the pump field strength used in Ref. 1 is far beyond this threshold.

Another important test for the lasing mode is the robustness of discontinuities in the reflection coefficient with respect to changes in the direction of the incident probe pulse. Analysis of the angular dependence of discontinuities in the reflection coefficient also tells us about Floquet parametric instabilities for modes with finite momentum along the surface of the material. In Fig. 11, we compare light reflection coefficients *R*_{s} for different angles of incidence. At frequencies 20 and 39 THz, discontinuities remain albeit with a reduced magnitude. However, the discontinuity at 30 THz disappears and is replaced by a sharp peak. To understand the persistence of the discontinuities at 20 and 39 THz, we observe that Floquet eigenmode (33) does not depend on the direction of light inside the medium. Therefore, for any incidence angle, the phase matching condition in (37) is satisfied at these frequencies and the eigenmode is infinitely extended inside the slab. From the perspective of two polariton scattering shown in Fig. 1, it is possible to find a resonant scattering process that involves two polariton bands when the scattered photons have a fixed momentum along one of the axes. On the other hand, the instability arising from the small-momentum resonant polariton scattering within the upper polariton band does not exist when the scattered phonons have a finite component along one of the axes. Hence, incident light cannot penetrate the surface of the slab and the discontinuity at 30 THz for normal incidence becomes a peak for finite incidence angle. In real experiments, incident beams strike the surface with a finite range of incident angles. In light of the earlier discussion, we expect that discontinuities in the reflection coefficients *r*_{s} and *r*_{i} at frequencies 20 and 39 THz should be more readily observed in experiments. Our result in Fig. 11 is for S-polarized wave incidence. Regarding the P-polarized case, we note that the effective potential (16) has cubic terms, which change the result quantitatively. However, the scattering process at 20 and 39 THz should be robust as far as fourth order terms are present in the effective potential.

Before concluding this section, we will discuss effects of including more modes, arising from other harmonics in the Floquet analysis. In particular, we observe that to linear order in $Ep2$, we have mixing between *ω*_{s} and frequencies *ω*_{i} = 2*ω*_{p} − *ω*_{s} and *ω*_{sum} ≡ 2*ω*_{p} + *ω*_{s}. So far, our discussion focused on the role of *ω*_{i}. We can extend our analysis to include *ω*_{sum}. The results of the numerical solution of the model with three coupled frequencies are shown in Fig. 12. We find few differences with the solution of the two frequencies case. Generally, the wave vector corresponding to *ω*_{sum} does not match the wave vector corresponding to *ω*_{s}. Thus, hybridization between *ω*_{s} and *ω*_{sum} is strongly suppressed. One exception is the points shown in Fig. 13. Closer examination of Fig. 12 shows that, indeed, for the frequency *ω*_{s} = 23.8 THz, analysis that includes three mode mixing exhibits a small discontinuity of *R*_{s}. The origin of this discontinuity is a resonant excitation from LP branch point S to UP point S′ (see Fig. 13) by the Floquet modulation of the medium (coherent oscillations excited by the pump pulse). Generalization of this argument suggests that we can also find higher-order resonances when states in the upper and lower bands at the same momentum differ in energy by 2*nω*_{p}. All of them are expected to exhibit lasing, although they involve higher non-linearities and, therefore, require a larger threshold pump intensity for lasing to occur.

Another important case is when the pump pulse is short enough as to drive polaritons in a wide window of frequencies. Our analysis, based on Floquet theory, is not easily adapted to treat such a situation. However, this may be amenable to numerical simulations, which we leave for future research.

## V. CONCLUSIONS

In this paper, we developed an analytic theory for analyzing pump and probe experiments in a polariton system with resonant pumping of the upper polariton mode. We computed the response of the system to the probe pulse by solving the coupled set of phonon and Maxwell equations linearized around the coherently oscillating upper polariton mode. Our method utilizes the Floquet approach in which pump-induced coherent polariton oscillations provide time-periodic modulation of the medium properties. The key feature of such Floquet–Fresnel light reflection analysis is frequency mixing analogous to that of parametric amplifiers. The probe light incident at frequency *ω*_{s} results in the appearance of the reflected and transmitted waves at frequency *ω*_{i} = 2*ω*_{p} − *ω*_{s}, where *ω*_{p} is the frequency of the polariton mode excited by the pump. In the context of non-linear optics, such complementary frequency modes are usually referred to as idler modes. Earlier experiments in SiC demonstrated light amplification of the probe beam upon reflection, which in our analysis can be understood as originating from photon pairs produced by the oscillating polariton mode. We find that when the pump-induced polariton amplitude is strong enough, the system should not only exhibit amplification but also undergo a true parametric instability. Such instability can be understood as arising from the two polariton scattering process that satisfies energy and momentum conservation (see Fig. 1). This instability mechanism is analogous to stimulated resonant scattering of exciton-polaritons in the case of pumping at the “magic angle” demonstrated by Baumberg *et al.*^{32} While earlier pump and probe experiments with polaritons in SiC were performed using broadband probe pulses, we showed that using frequency resolved probes should provide new insights into properties of the photoexcited system. In particular, reflection coefficients for both the signal and idler components should exhibit discontinuities in the frequency dependence arising from the parametric instability. We note, however, that true discontinuities are present only in the idealized situation of a uniformly excited infinitely large crystal and infinitesimal intensity of the probe beam. Under realistic experimental conditions of finite penetration depth of the pump pulse and small but finite intensity of the probe pulse, discontinuities in the reflection coefficient should become sharp crossovers. We discuss the possibility of the “afterglow” that can appear after the probe pulse is gone, which can be understood as the parametric instability triggered by the probe pulse. Another interesting feature in the frequency dependence of the reflection coefficients that we discuss is sharp peaks arising from the surface modes that appear near the edge of the reststrahlen band for both the signal and idler frequency components.

In addition, another scattering process may be realized when we pump polaritons at finite momentum in the lower polariton branch. In Fig. 13, we show the situation when the pumped polaritons have a macroscopic occupation at U_{p} and scatter to U′ and U″. This is a direct analog of the exciton resonant scattering observed by Baumberg *et al.*^{32} In this process, the scattering occurs in the lower branch, in contrast to the situation shown in Fig. 1, where the polaritons are scattered into both the upper and lower branches. Our Floquet formalism can be extended to study the optical properties of the system in the presence of polaritons pumped at finite momentum.

We envision several potential technological applications of the system discussed in this paper. Lasing instability allows us to obtain coherent sources at 20 and 40 THz starting from radiation at 30 THz. Finding strong sources at lower frequencies is particularly valuable in the mid-IR domain. Moreover, we expect that the parametric frequency can be tunable. In this paper, we analyzed the simplest situation of pumping at zero wave vector. One can extend Floquet–Fresnel analysis to the situation when the pump excites upper polaritons at finite momentum. This will modify the resonant scattering process shown in Fig. 1 and should result in the modification of lasing frequencies. We also note that resonant polariton scattering discussed in our paper results in the production of entangled THz photons, which may find applications in quantum information. Another promising direction is to consider pumping regime when the system is about to develop parametric instability, so the poles move from the lower to the upper half-plane. At the transition, we should find polariton modes at real frequencies, which indicates that they can propagate in the material for a long time without getting damped and becoming unstable. This is the regime of exact balance between dissipation and Floquet pumping. Such modes can be used for efficient frequency filtering. By changing the momentum of the pump photons, we expect to be able to modify the frequency that can be filtered.

Discussion presented in our paper considered a specific type of optical non-linearity. We expect that mechanism of Floquet parametric instability that we discussed here may be present in a broad class of systems. When a phonon (or two phonons) can scatter into a pair of collective modes such that energy and momentum are satisfied, we expect to find Floquet parametric instabilities.

## ACKNOWLEDGMENTS

The authors thank A. Cartella, A. Cavalleri, J. Faist, B. Halperin, A. Imamoglu, and A. Rosch for useful discussions and valuable comments. The authors thank A. Leitenstorfer for a careful reading of the manuscript and insightful comments. The authors acknowledge support from the Harvard-MIT CUA, Harvard-MPQ Center, AFOSR-MURI: Photonic Quantum Matter (Award No. FA95501610323), the DARPA DRINQS program (Award No. D18AC00014), the Israel Science Foundation (Grant No. 1803/18), and the National Science Foundation through a grant to ITAMP at the Harvard–Smithsonian Center for Astrophysics. S.S. acknowledges support from JSPS Overseas Research Fellowships. D.P. acknowledges the Aspen Center for Physics, where part of this work was done, with support of the NSF under Grant No. PHY-1066293.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: NONLINEAR POLARITON ANALYSIS

The nonlinear dynamics of the coupled phonon amplitude and electric field in SiC is described in Eqs. (19) and (20). We begin by linearizing Eq. (19) with respect to the signal and idler frequency components,

where we used relation (21) to express *E*_{p} = *νQ*_{p}. Separating *Q*_{s} and *E*_{s} terms of the last equation, we obtain

where

We write analogous equations for the idler frequency component and solve for the two phonon amplitudes, *Q*_{s} and $Qi*$ in terms of the electric fields *E*_{s} and $Ei*$. Assuming that the polariton amplitude *Q*_{p} is not too large, we keep terms up to quadratic order in *Q*_{p}. We find

and

We now consider the electric field equation [Eq. (20)] and separate the *ω*_{s} frequency component of this equation,

where

Repeating the same calculation for the *ω*_{i} frequency component, we obtain

### APPENDIX B: ANALYSIS OF LIGHT SCATTERING FROM A SLAB OF FINITE THICKNESS

In this section, we present analysis of light reflection from a slab of finite thickness in the case of normal angle of incidence. We set the boundaries of the material to be at *z* = 0 and *z* = −*L*. In the air, for *z* > 0, we use (43) and (44). In the insulator, for −*L* < *z* < 0,

In the air, for *z* < −*L*,

In Sec. III, we discussed the case of a semi-infinite slab and included only evanescent waves $e\kappa \xb1z$ that vanish as *z* → −∞. In the case of a slab with finite thickness, we need to include solutions corresponding to evanescent waves propagating from the lower interface, which correspond to the $e\u2212\kappa \xb1z$ terms in equation (B2). The reflection and transmission coefficients can now be obtained by analyzing boundary conditions for the Maxwell equations at *z* = 0 and *z* = −*L*. We have for *z* = 0,

and for *z* = −*L*,

where

We assume Re[*κ*_{+}] < Re[*κ*_{−}] without loss of generality. Then, the exponents in (B11)–(B14) satisfy

Here, the equation in the middle holds when Re[*κ*_{+}] = 0. From these relations, one finds that $t3\u2032$ and *t*_{4} are exponentially small, which allows us to simplify the analysis by taking $t3\u2032=0$ and *t*_{4} = 0. We can now write equations (B3)–(B10) as

at *z* = 0 and

at *z* = −*L*. Solving these equations, we obtain *r*_{s}. We do not present explicit expression for *r*_{s} because it is somewhat cumbersome and is not needed for our discussion. Instead, we discuss the condition for *r*_{s} to have a pole, which corresponds to a polariton eigenmode of the system. The condition for the pole of *r*_{s} to be zero is given by the following equation:

Here, A is the coefficient determined by the boundary conditions that can be written as

Note that *A* is a function of *κ* but is independent of *L*. In all cases relevant to our discussion, *A* is finite, so we can take the logarithm of equation (B25) and find the eigenmode conditions

where *n* is an arbitrary integer.

Physically, eigenmodes correspond to waves that undergo an infinite number of reflections between the upper and lower boundaries of the slab. Equation (B26) describes the balance between loss/gain during propagation in the slab coming from the dissipation/Floquet pumping (Re[*κ*_{+}]*L*) and loss during the reflection at the interface with air (*A*). Equation (B27) is the usual condition of phase matching after one full cycle. In the limit of *L* → ∞, these equations reduce to

where *a* is an arbitrary real number because $nL\pi $ becomes a dense set. In this case, Eq. (B29) becomes immaterial as there is always some *a* that satisfies (B29).

We also verified the appearance of poles in *r*_{s} by solving the reflection problem (B1) and (B2) numerically. In Fig. 6, we show numerical results for *r*_{s} for slabs with thicknesses *L* = 3 × 10^{−3} and 9 × 10^{−4} m, respectively. Poles correspond to the small white dots in Fig. 6.

A useful way of understanding the solutions in the limit of large *L* is to consider solutions of the equation Re[*κ*] = 0 as defining lines in the complex *ω* plane. By imposing additional condition (B27), we find discrete points corresponding to polariton eigenmodes. The distance between solutions obtained using this procedure scales as *O*(1/*L*) in agreement with our intuition for the mode quantization in a cavity of length *L*. In our numerical calculation, for the solution of (B27), we have an eigenmode close to the real *ω* axis, which correspond to the usual standing wave solutions. We also find solutions that extend into the complex *ω* plane, which are special to the Floquet system. These solutions extend between −1.7 THz < Im[*ω*_{s}] < 1.7*x* THz. As we increase the thickness of the slab, the poles become denser, and in the limit of infinite slab width, we find lines of discontinuities of the reflection coefficient.

### APPENDIX C: ANGULAR DEPENDENCE OF REFLECTIVITY

In this section, we consider the problem of light reflection from the Floquet insulator described in Eq. (26) for a finite angle of incidence. The goal of this analysis is twofold. First, understanding of the angular dependence of the reflection coefficient is important for analyzing results of light reflection the probe part of the experiments. Second, when discussing lasing instabilities in Sec. III C, we only considered modes that have zero momentum in the direction parallel to the surface of the material. We expect that unstable modes with a finite in-plane component are also present in the system. We remind the readers that our method of analyzing lasing instabilities is to study poles of the reflection coefficient. When the incident light has a finite momentum along the surface of the sample, the component of the momentum parallel to the interface is conserved upon reflection. Hence, by analyzing poles of the reflection coefficient at finite angle of incidence, we obtain information about lasing modes with a finite parallel component of momentum. Experimentally, they correspond to modes radiating at finite angle to the normal of the sample.

Following the usual convention for analyzing reflection problems of electro-magnetic waves, we separately discuss the cases of s- and p-polarization of light. In the case of s-polarization, the electric field is parallel to the surface of the material for the incident wave and for the reflected and transmitted waves. In the case of p-polarization, the magnetic field is parallel to the interface for the incident, scattered, and reflected waves. In this section, we present the results for the reflection coefficients and focus on the discussion of their physical implications. Details of the derivation can be found in Appendix D.

To analyze reflection of the S-polarized wave, it is convenient to formulate Maxwell equations using the electric field. We separate the signal and the idler components of the electric field and write

Here, *p*, *q*_{s}, *q*_{i}, *q*_{±} are imaginary wavevectors that satisfy

and *r*_{s}, *r*_{i} and *t*_{s+}, *t*_{s−}, *t*_{i+}, *t*_{i−} are the reflection and transmission coefficients, respectively. These notations are summarized in Fig. 14. The boundary conditions for the Maxwell equations require that the in-plane component of the wave vector to be the same in the air and inside the medium. Note that in (C1) and (C2), *p* is purely imaginary, corresponding to the real in-plane momentum. We find

For the P-polarized wave, we start with the magnetic field

Analysis presented in Appendix D shows that in this case, the reflectivity is given by

To understand the effect of finite *p* in (C6) and (C9), we recall that discontinuities of the reflection coefficient and the corresponding lasing modes appear when Im[*κ*^{2}] crosses zero and Re[*κ*^{2}] < 0. The latter condition is important since for Re[*κ*^{2}] ≥ 0, there is no discontinuity of Im*κ*. If we examine (C5) and recall that *p*^{2} is negative, we realize that Re[*κ*^{2}] − *p*^{2} can become positive when |*p*| is sufficiently large. We also recall that we defined the unstable mode to be *κ*_{+} and that at the lasing transition, Im[*κ*^{2}] becomes equal to zero. Hence, the lasing instability becomes suppressed when

where *ω** is the corresponding instability frequency in the bulk (we remind the readers that instabilities at 20 and 39 THz correspond to the same unstable Floquet eigenmode, and for both of them, we should take *ω** = 20 THz). When condition (C10) is not satisfied, we find that instabilities persist although the jump in *q*_{±} becomes suppressed. This explains why discontinuities in *r*_{s} at 20 and 39 THz become suppressed with the increasing angle of incidence (see Fig. 11). This suggests that lasing instability exists for all modes that have |*p*| smaller than the critical value of $\omega *c\u03f50$, beyond which we can no longer find transverse momentum to match the resonant scattering process shown in Fig. 1. For the unstable mode at *ω** = 29 THz, we observe that Re[*κ*^{2}] is small; hence, this instability gets suppressed at small values of |*p*|.

### APPENDIX D: REFLECTION PROBLEM FOR FINITE ANGLE OF INCIDENCE

We present the details of Appendix C. When the incident wave is not perpendicular to the surface, the reflection and transmission coefficients depend on the polarization of the incident light. A common approach is to decompose incoming waves into S- and P-polarization components. For the S-polarized light, the electric field is given by (C1) and (C2). We take the direction of the electric field parallel to the *y* axis. We can find expressions for the magnetic field using Eq. (7) and relating it to **E**_{s}(**r**) and $Ei*(r)$,

The tangential component of the electric field and the magnetic field should be continuous at *z* = 0. Therefore, from (C1), (C2), (D1), and (D2), we find

Electromagnetic eigenmodes inside the slab are solutions of the Floquet problem discussed in Sec. II C; hence, *t*_{s±} and *t*_{i±} satisfy relation (34). Substituting (C1) and (C2) into (34), we obtain

The reflection problem is solved for the P-polarized incident light in the same way. We can use Eq. (7) to find the electric fields from the expressions for **B**_{s}(**r**) and $Bi*(r)$ in (C7) and (C8),

Boundary conditions for the tangential component of the electric field and the magnetic field at *z* = 0 give

We note that we ignore cubic terms in the effective potential (16) here. In the P-polarized case, the effective potential of 4H-SiC has cubic terms, in general. We should take these terms into account for more accurate analysis.