
\section*{Introduction}

\subsection*{Motivation}
Optical pump-probe experiments provide a powerful avenue to explore collective dynamics in correlated quantum materials. Several experiments have demonstrated dramatic changes in the terahertz and mid-IR reflectivity following photo-excitation. In experiments on  $\rm{YBa_2Cu_3O_{6.5}}$ cuprate superconductors, parametric resonances caused by internal Josephson plasma oscillations can induce extra "edge" features in the optical reflectivity \cite{Kaiser14,vonHoegen19,Marios20}. On the other hand, the phonon-polariton system $\rm SiC$\cite{Limonov17} and the bulk superconductor $\rm K_3 C_{60}$\cite{budden20, dolgirev21} provide examples where oscillations inside these materials result in significant reflectivity enhancement, possibly exceeding unity. A unified interpretation of these seemingly disparate pump induced features can be formulated in terms of Floquet theory under the assumption that collective excitations create a system with time-periodic properties. Previous analysis has demonstrated the existence of four types of drive induced features in reflectivity, depending on the relative strength of parametric driving and dissipation\cite{FloquetFresnel21}. Edge-like features occur from interference between different Floquet components of the transmission channels when dissipation is strong compared to the oscillation amplitude. In the opposite regime, strong amplification of reflectivity occurs when the parametric drive is not compensated by dissipation and the material exhibits a lasing instability. In fact, such reflectivity features can serve as reporters of a lasing instability, indicating that the effective Floquet medium can be used as a gain medium in a laser.

\begin{figure*}[t!]
    \centering
    \includegraphics[width=1.\textwidth]{Figure1.png}
    \caption{Schematic mechanism that leads to amplification of THz optical reflectivity following high frequency pumping. a) An ultrafast laser pulse (0.5eV, 150 fs) photoexcites electrons through direct dipole transitions between the valence and conduction bands. The photo-excited electrons generate pairs of phonons through the electron-phonon interaction squeezing the phonon field (see equation~(\ref{eq:electronphononInt})). b) Reflectivity amplification of pumped $\rm Ta_2NiSe_5$. Coherent phonon field fluctuations oscillate at twice the phonon frequency, $\omega_{ph}$, creating an effective THz Floquet medium. Parametric driving from the phonon oscillations can create pairs of photons at the signal and idler frequency once stimulated by the probe pulse. This enhances the reflectivity and also scatters counter-propagating light oscillating at the idler frequency. c) Relative change in the reflectivity as a function of frequency subsequent to photoexcitation.  The experimental results are shown at two different temperatures together with the theoretical fit which considers parametric amplification by a 9.4 THz oscillating field. }
    \label{fig:process}
\end{figure*}

Here, we develop a theory to explain the experimentally observed pump-induced terahertz amplification of the reflectivity in layered quasi-1D semiconductor $\rm{Ta_2 Ni Se_5}$, a potential excitonic insulator material. An excitonic insulator phase can emerge when bound electron-hole pairs form a condensate in a semi-metallic phase opening up a gap and causing a phase transition to an insulating state\cite{Jerome1967}. In $\rm{Ta_2 Ni Se_5}$, alternating chains of Ta and Ni are aligned along the axis, forming a sheet in the $\textit{ac}$-plane, and multiple layers are stacked along the $\textit{b}$-axis (see Fig. 1A). An exchange of electron between adjacent Ta and Ni chains along the $\textit{a}$-axis creates an exciton across the chain (along the $\textit{c}$-axis) \cite{kaneko2013,Werdehausen2018,Sugimoto2018,Mazza2020,Lu2017,Larkin2017,Larkin2018,Mor2018,Seki2014}. Previous studies have reported a second order monoclinic-orthorhombic phase transition at $\rm{T_c}$ = 326 K. This leads to an insulating gap below the critical temperature which grows to be 160 meV at low temperatures. The pump-probe experiment demonstrates that above-gap frequency optical excitation induces reflectivity amplification which is maximal at $4.7$ THz, corresponding to an IR-active $\rm{B_{3u}}$ phonon mode confirmed by previous infrared measurements \cite{Larkin2018,Larkin2017}. This turns out to be a surprising result given that the high-frequency pump cannot directly excite the collective THz modes. We show phenomenologically and using ab-initio calculations that this dramatic down conversion is caused by a nontrivial interplay of electron-phonon interactions and phonon non-linearities.

The multi-step process described in this article, converting high frequency pumping to THz amplification of reflectivity, is schematically illustrated in figure~\ref{fig:process}(a) and is outlined as follows: a) The pump excites electrons from the valence bands to the conduction bands through direct dipole transitions. b) The photo-excited electrons exhibit an unusual quadratic coupling to the $4.7$ THz phonons (see Results): 
\begin{equation}
    H_{el - ph} = \sum_k g_k n_{el,k} Q^2,
\label{eq:electronphononInt}
\end{equation}
where $Q$ is the IR-phonon displacement, $n_{el,k}$ the photo-excited electron occupation and $g_k$ the effective electron-phonon coupling. 
The pump is non-resonant with the IR active terahertz phonons and thus does not directly initiate the enhanced reflectivity dynamics. Instead, IR-phonon pair generation occurs via a Raman process caused by pump induced changes in the electronic occupation. The result is that the expectation value of the phonon displacement is zero, $\expe{Q} = 0$, but the fluctuations are squeezed and coherently oscillate at twice the phonon frequency, $\expe{Q^2(t)} = \expe{Q^2}_0 + A \cos(2\omega_{\rm ph} t) $. c) The squeezed phonon oscillations and phonon nonlinearities create a Floquet material oscillating at $2 \omega_{\rm ph}  = 9.4 $ THz. d) Parametric resonances due to the oscillating field occur primarily around $\omega_{\rm ph}$, parametrically amplifying the reflectivity as depicted in Figure~\ref{fig:process} (b). 

The theory developed in this article addresses two inter-connected questions: 1) Can a Floquet material oscillating at 9.4 THz provide parametric amplification similar to experimental observations? 2) What is the origin of the 9.4 THz oscillation?.


\subsection*{Summary of results} 
We begin our analysis by showing that the parametrically amplified reflectivity at $4.7$ THz is consistent with the presence of a coherently oscillating Raman mode at twice the frequency, $9.4$ THz. The theoretical fit captures the experimental amplification profile, plotted in figure~\ref{fig:process}(c). A schematic Hamiltonian describing parametric amplification of the optical reflectivity in the presence of a coherently oscillating mode at $\omega_d = 9.4$ THz is given by:
\begin{equation}
    H_{ampl.} = X(t) a_{1,k}^\dag a_{2,-k}^\dag + h.c.,
\end{equation}
where $X(t)$ is a Raman mode and $ \{ a^\dag_{1,k} , a^\dag_{2,-k} \} $ are the creation operators of photons with opposite momentum which may be associated with different phonon-polariton bands $\{1, 2\}$ at frequencies parametrically resonant with the drive, $\omega_{1} (k) + \omega_2 ( - k)  = \omega_d$. As shown in figure~\ref{fig:process}(b), Floquet drive can generate pairs of photons at the signal frequency of the incoming probe, $\omega_s$, and a photon oscillating at the idler frequency, $\omega_{\rm id} = \omega_d - \omega_s$. Reflectivity amplification is then caused by stimulated emission of photon pairs by the driven material. To compute the reflection coefficient, we use degenerate Floquet perturbation theory to construct the eigenmodes and solve the Fresnel equations at the boundary.

We then proceed to investigate the origin of the $9.4$ THz oscillation. We show that electron-phonon coupling between electronic bands involved in the photo-excitation and the 4.7 THz IR-phonon naturally leads to phonon squeezing. Phonon squeezing oscillations have twice the frequency of the phonon and act as a Raman mode. We find that this mechanism is generic and should apply to all IR-phonons present in the system. However, for phonon squeezing oscillations to create an effective Floquet material that can significantly enhance the reflectivity, the phonon needs to be both strongly coupled to the photo-excited bands and also have a strong coupling to the electric field. 

We hence perform ab-intio calculations to determine both the magnitude of the electron phonon coupling and the IR activity of the different phonon modes. We find that, indeed, the $4.7$ THz IR-phonon is special, exhibiting strong electron-phonon coupling and appreciable IR activity as compared to other IR-phonons in the same frequency range. Moreover, we confirm that phonon fluctuations are responsible for the THz amplification. The theory shows that the predicted reflectivity amplification is a reporter of a lasing instability, implying that pumped $\rm Ta_2 Ni Se_5$ could, in principle, serve as a Floquet gain medium for THz lasing.

As a final intriguing aspect, we find that the $4.7$ THz mode is strongly coupled to the elusive and hotly debated order parameter of $\rm Ta_2 Ni Se_{5}$\cite{Baldini20,Sugimoto2018,Mazza2020,Golez2020,Sun2020,Subedi2020,kaneko2013}, which is thought to have some excitonic insulating character. A number of state-of-the-art experiments have been performed to shed light on the origin of the phase transition at $326$ K as well as what role the excitonic and lattice instabilities play in this phase transition process. However, the results show contradictory results \cite{Werdehausen2018,Werdehausenagain2021,Lu2017,Larkin2017,Larkin2018,Mor2017,Mor2018,Seki2014,Kim2020,Kim2021,Bretscher2021,Bretscheragain2021,Baldini20,Volkov2021,Pal2020,Tang2020,Chen2020,Seo2018,Wakisaka2009,Nakano2018,Saha2021,Lee2019,Suzuki2021} and the question whether $\rm Ta_2 Ni Se_{5}$ hosts an excitonic insulator phase remains an open question. Using DFT frozen phonon calculations, we find that the electron-phonon coupling for the $4.7$ THz phonon effectively vanishes in the high temperature orthorhombic phase. This indicates that the THz parametric amplification of reflectivity is mediated by the low temperature phase and is expected to be sensitive to the phase transition. This interplay between the order parameter and phase transition has been observed in pump-probe \cite{Werdehausen2018,Werdehausenagain2021}, near infrared \cite{Bretscher2021,Bretscheragain2021}, and time-resolved ARPES experiments \cite{Baldini20,Tang2020,Suzuki2021,Denis2021} , and has unraveled key ingredients in the physics of $\rm{Ta_2NiSe_5}$ in addition to providing a clearer understanding of the nature of the order parameter.  The link between the 4.7 THz phonon and the order parameter opens new avenues for using parametric amplification as a novel method to track order parameter dynamics.







\section*{Results}

\subsection*{Parametric amplification of reflectivity in pumped $\rm{Ta_2 Ni Se_5}$}
\label{sec:paraAmpl}

\subsubsection{Equilibrium reflectivity} 
The reflectivity of a material is captured by the frequency dependent refractive index appearing in the Maxwell equations:
\begin{equation}
    \left( \frac{n(\omega)^2\omega^2}{c^2} - k^2 \right) E = 0
\end{equation}
where $E$ is the electric field and all information about phonons and other IR active modes is encoded in $n(\omega)$. We assume that the probe corresponds to an electromagnetic wave reflected from the sample at normal incidence. The propagation direction is along the $b$-axis of the crystal, which we refer to as the $y$-direction, whereas the electric field points in the $a$-direction which we refer to as the $z$-axis. The refractive index can be directly extracted experimentally from the complex reflection coefficient at normal incidence in an equilibrium system\cite{WOOTEN72}:
\begin{equation}
    n(\omega) = \frac{1 - r (\omega)}{1 + r (\omega)}.
\label{eq:indexRefl}
\end{equation}

\subsubsection{Eigenstates in Floquet state}
Once we have obtained the refractive index from the equilibrium reflectivity, we model the Floquet material as experiencing a parametric drive oscillating at frequency, $\omega_d$, which mixes signal and idler frequencies:
\begin{equation}
    E(t) =  e^{ i k y} \left( E_s e^{ - i \omega_s t} + E_{id}e^{i \omega_{id} t} \right).
\label{eq:ansatzTNS}
\end{equation}
where $\omega_s$ is the frequency of the incoming probe and $\omega_{id} = \omega_d - \omega_s$. Such mixing corresponds to degenerate perturbation theory in Floquet systems and the idler component is the nearest Floquet band contribution to $\omega_s$ which is responsible for parametric instabilities \cite{Marios20, Sho19, Eckardt15}. The oscillating mode is included phenomenologically through a time-periodic contribution to the electric permittivity:
\begin{equation}
    \delta \epsilon(t) = 2 A_{\rm drive} \cos(\omega_d t).
\end{equation}
Using the ansatz in equation~(\ref{eq:ansatzTNS}), the equations of motion in the Floquet state for the different oscillating components of the electric field become:
\begin{subequations}
\begin{align}
\left( \frac{n(\omega_s)^2\omega_s^2}{c^2} - k^2 \right) E_s + A_{dr} E_{id} =& 0 \\ 
\left( \frac{n(\omega_{id})^2\omega_{id}^2}{c^2} - k^2 \right) E_{id} + A_{dr} E_s  =& 0.
\end{align}
\label{eq:FloquetEigen}
\end{subequations}
To compute the reflectivity at normal incidence, we first find the allowed $k$ values for a given $\omega_s$. Due to the coupling of signal and idler components, two such $k$ values exist, associated with two transmission channels, both of which oscillate at signal and idler frequencies. The transmission channels correspond to eigenvectors of the Floquet equations of motion inside the material,
\begin{equation}
    E_{i} = t_i E_0 e^{i k_i y} \left(e^{  - i \omega_s t}  + \alpha_i e^{ i \omega_{id} t} \right),
\end{equation}
where $\alpha_i$ is the relative amplitude of the signal and idler component in the eigenvector of wave-vector $k_i$, $t_i$ is the transmission coefficient of the $i$-th channel and $E_0$ the amplitude of the incoming field. 

\subsubsection{Floquet-Fresnel equations}
In a reflection problem, the eigenvalue equation~(\ref{eq:FloquetEigen}) enables computation of the transmitted wavevectors as a function of a fixed frequency set by the incoming light. However, the answer is given in terms of $k_i^2$ rather than $k_i$ and the correct root is chosen such that the field vanishes at infinity, $\mbox{Im} \{k_i \} > 0  $. The reflectivity is computed by solving the Floquet-Fresnel equations, matching electric and magnetic fields parallel to the surface both at frequency $\omega_s$ and at frequency $\omega_{id}$. Inside the material we have the electric field:
\begin{equation}
    E_{mat} = E_0 \sum_i t_i e^{i k_i y} \left( e^{- i \omega_s t} + \alpha_i e^{ i \omega_{id} t} \right) ,
\end{equation}
and for vacuum we have:
\begin{equation}
\begin{split}
    E_{vac} = &E_0 \left(e^{ i \omega_{s}/c y - i \omega_s t} + r_s e^{- i \omega_s/c y - i \omega_s t} \right) +\\ 
    &E_0 r_{id} e^{  i \omega_{id}t. + i \omega_{id} y }.
\end{split}
\end{equation}
Using the homogeneous Maxwell equations, $ \nabla \times E = - \partial_t B$ to compute the magnetic field and matching boundary conditions at $y = 0$, we obtain the Fresnel equations for the driven system:
\begin{subequations}
\begin{align}
1 + r_s =& t_1 + t_2, \\
1 - r_s =& \frac{k_1}{\omega_s} t_1 + \frac{k_2}{\omega_s} t_2, \\
r_{id} =& t_1 \alpha_1 + t_2 \alpha_2, \\
r_{id} =& \frac{ k_1}{\omega_{id}} t_1 \alpha_1 + \frac{k_2}{ \omega_{id} } t_2 \alpha_2.
\end{align}
\label{eq:FloquetFresnel}
\end{subequations}
To fit the data, we choose a drive at $9.4$ THz and allow for small changes in the static properties of the system such as a photo-induced conductivity stemming from photo-excited charge carriers. The fitted parameters are given in the Materials and Methods section.

\subsection*{ Phonon squeezing initiated by photoexcited electrons} 
\label{sec:eltoPhon}

\subsubsection{Electron phonon coupling}
In this section, we develop a microscopic theory of the coupling of electronic bands to IR active phonons. In the dipolar gauge\cite{Jiajun2020}, the effective dipole of the phonon is linearly coupled to the electronic dipole. For two electronic bands with an allowed dipole transition, this coupling takes the form:
\begin{equation}
    H_{\rm el- ph} = Q \sum_k \lambda_k \left( \hat{c}^\dag_{1,k} \hat{c}_{2,k} + \hat{c}^\dag_{2,k} \hat{c}_{1,k} \right),
\end{equation}
where $Q$ is the phonon coordinate, $\hat{c}_1$ and $\hat{c}_2$ the annihilation operators of the two electronic bands and $\lambda_k$ the electron-phonon coupling matrix element as a function of momentum $k$. Using a Schrieffer-Wolff transformation, we "integrate out" the linear electron-phonon coupling, which is non-resonant due to the different energy scales between the IR-phonon and the electronic transition, to reveal the resonant non-linear coupling between the electron occupation number and the phonon squeezing operator, $Q^2$. In the materials and methods material, we show that this procedure leads to the effective coupling:
\begin{equation}
    H_{\rm el-ph, eff.} = \sum_k \frac{\lambda_k^2}{ \Delta_k}  \left( n_{1,k} - n_{2,k} \right) Q^2,
\label{eq:Squeezingcoupling}
\end{equation}
where $\Delta_k = E_{1,k} - E_{2,k}$ is the energy difference between the electronic states and $\hat{n}_{i,k} = \hat{c}^\dag_{i,k} \hat{c}_{i,k } $ the number operator.  The above equation applies to phonons that are coupled to independent pairs of electronic bands. If three or more bands are simultaneously coupled to a specific phonon the effective electron-phonon Hamiltonian is more complicated but with similar qualitative features, such as the  coupling of electron density to the square of the phonon coordinate.

\subsubsection{Phonon squeezing}

Once electrons have been photo-excited, the finite number of optically excited electrons quenches the frequency of the phonon:
\begin{equation}
    H_{\rm el-ph, eff.} = \frac{M f(t) Q^2}{2},
\label{eq:ParametricPhonon}
\end{equation}
where $M$ is the mass of the phonon. In equation~(\ref{eq:ParametricPhonon}), the parametric driving $f(t)$ comes from the photo-excited distribution of electrons that is strongly coupled to the phonon:
\begin{equation}
    f(t) = \sum_k \frac{2 \lambda_k^2}{M \Delta_k} \left( \expe{n_{1,k}}(t) - \expe{n_{2,k}(t) }\right),
\end{equation} 
The photo-excited electron dynamics are fast and can be approximated as a delta function
in time. More generally we can approximate the photo-excited distribution as having a characteristic life time, $f(t) = f_0 \theta(t)e^{-t/t_{decay}} $, but this does not change our conclusions. 

The above electron induced drive describes a Raman process that does not excite the phonon directly (i.e. $\expe{Q} = 0$). However, the squeezing operator starts oscillating at frequency equal to twice the phonon frequency, as shown in the materials and methods:
\begin{equation}
   \left( \partial_t^2 - \left( 2 \omega_{\rm ph } \right)^2  \right) \expe{Q^2} = - 4 f(t) \expe{Q^2}_0, 
\end{equation}
where $\expe{Q^2}_0$ is the equilibrium fluctuations and $\omega_{\rm ph}$ is phonon frequency at zero momentum which is renormalized by the coulomb force, $\omega_{\rm ph}^2 = \omega_{\rm ph,0}^2 + \frac{Z^2}{\epsilon \epsilon_0 M }$. 

To show that this phenomenon is related to squeezing, we expand $Q^2$ in terms of creation and annihilation operators, using $Q = \frac{a + a^\dag}{\sqrt{2 M \omega_{ph}}}$:
\begin{equation}
    Q^2(t) = \frac{1}{2 M \omega_{ph}} \left( a^\dag(t) a^\dag(t) + a(t) a(t) + a^\dag(t) a(t) + a(t) a^\dag(t) \right).
\end{equation}
Expectation values of $a(t) a(t)^\dag$ do not oscillate rapidly while the anomalous pairs $a^\dag(t) a^\dag(t)$ and $a(t) a(t)$ oscillate at twice the phonon frequency. As a result a state with phonon fluctuations, $Q^2$, that oscillate at twice the phonon frequency implies the existence of a condensate of phonon pairs $\expe{a^\dag(t) a^\dag(t)} \neq 0$. 






\subsubsection{Floquet matter}

Coherent oscillations of long lived modes, such as the phonon squeezing oscillations, turn matter into Floquet matter with time-periodic properties through non-linear interactions with the rest of the material. Following the same steps as in references \cite{Sho19,Michael21}, we show in the materials and methods how a lattice anharmonicity  proportional to $Q^4$ and squeezing dynamics lead to the signal and idler coupling in equation~(\ref{eq:FloquetEigen}). For this model, we find that the signal idler coupling is proportional to the IR activity of the mode, $Z$, the phonon anharmonicity, $u$, and the amplitude of the oscillations, $B$ : $A_{\rm drive} \propto Z^2 u B$.



 \subsection*{Ab-initio calculations : the 4.7 THz IR-Phonon}
 \label{sec:DFT}
 \begin{figure*}
    \centering
    \includegraphics[trim = 0pt 0pt 0pt 0pt,clip, scale = 0.4]{finalPlotv3.pdf}
    \caption{Summary of ab-initio calculations. a) Optical matrix elements along the \textit{a}-axis calculated for the frequency region between 0.4 to 0.9 eV. The size of the circles is proportional to the magnitude of the matrix elements and indicates which electronic valence and conduction bands are involved in photo-absorption. b) IR activity of the most IR active phonons; the phonon at 4.7 THz which dominates electron-phonon interactions in that region is highlighted with the mode character shown in the inset. A complete list of IR active phonons in this region for TNS is given in materials and methods. c) and d) show the effect of electron-phonon coupling on the band structure in the monoclinic and orthorhombic phase respectively. The electron-phonon coupling is captured by calculating the electronic band-structures using DFT in the presence of either negative or positive phonon displacement (frozen phonon approximation). The phonon 4.7 THz in the monoclinic phase is strongly coupled to the electronic bands 2 and 3 which are highlighted in the figure. In the high temperature orthorhombic phase the effective electron-phonon coupling is reduced by an order of magnitude. This indicates that the electron-phonon interaction is very sensitive to the order parameter.}
    \label{fig:DFTsummary}
\end{figure*}
The above discussion is generic and in principle applies to every IR-phonon inside a given material. In this section we use DFT calculations to determine which phonons make the dominant contribution to the parametric amplification of reflectivity in $\rm Ta_2 Ni Se_5$. 
We start our discussion by identifying which valence and conduction bands are involved in photo-absorption by evaluating the optical dipole transitions matrix elements. In Fig.~\ref{fig:DFTsummary}(a), we plot the momentum resolved optical contribution, for field polarization along the \textit{a}-axis, of each band to dipole allowed transitions within the experimentally relevant energy window between 0.33 and 0.9 eV set by the pump parameters. For a given valence(conduction) band and k-point, the optical contribution is defined by summing the square of the dipole transition matrix elements associated with the transition. We find that it is the first three conduction bands that are predominantly excited by the pump. 

Turning our attention to phonons, we use ab-initio calculations to identify all IR active phonons. In particular, we find a number of phonons in the low temperature monoclinic phase between 4 and 5 THz shown in the supplementary material. In Fig.~\ref{fig:DFTsummary} (b) we plot the IR activity of phonons as a function of frequency and identify the 4.7 THz phonon which, as shown below, turns out to be the most strongly electron coupled phonon. 

We compute the electron-phonon interaction between IR active phonons and electrons. This is to identify the phonons that dominate the interaction with the photo-excited electronic bands. To accomplish this we use the method of frozen phonons. In this approach, the electronic bands are recalculated with the lattice shifted along a phonon eigendisplacement. To quantify the coupling strength of a specific phonon to a particular electronic band we integrate the energy changes of the band in the presence of the phonon over the Brillouin zone as outlined in the supplementary material. We find that in the vicinity of $4.5$ THz, relevant to the experimental observables, the $4.7$ THz mode has roughly an order of magnitude larger electron-phonon coupling strength compared to the nearby phonons. In particular, as shown in Fig.~\ref{fig:DFTsummary}(c), the 2nd and 3rd conduction bands (in this paper we number the conduction bands from lowest to highest in energy) are significantly renormalised in the presence of the 4.7 THz phonon. Since the phonon is mostly coupled to two electronic bands, the frozen phonon calculations are consistent with equation~(\ref{eq:Squeezingcoupling}) leading to the two bands shifting in energy by an equal and opposite amount given by:
\begin{equation}
\begin{split}
   \left.E_{k,3}\right|_{\expe{Q} \neq 0} -  \left.E_{k,3}\right|_{\expe{Q} = 0} =&  - \left(\left.E_{k,2}\right|_{\expe{Q} \neq 0} -  \left.E_{k,2}\right|_{\expe{Q} = 0}\right)\\
   =& \expe{Q}^2 \frac{\lambda_k^2}{\Delta_k \left( 1- \frac{\omega_{ph,0}^2}{\Delta^2_k} \right)}
\end{split}
\end{equation}
where $E_{k,2}$ and $E_{k,3}$ is energy of conduction band two and three at momentum $k$, and $\Delta_k$ is the energy difference of the two bands in equilibrium. Therefore, in this case, calculating the energy shifts as a function of momentum in the frozen phonon approximation allows for direct computation of the electron-phonon interaction.

To summarize, we use ab-initio calculations to identify the electron bands excited by 
pumping and to establish the IR-active phonons in $\rm Ta_2 Ni Se_5$. Subsequently, frozen phonon calculations allowed for the identification of the phonon with the dominant electron-phonon coupling to the excited electronic bands. This leads to the identification of the $4.7$ THz mode as responsible for the parametric amplification observed in experiments through the mechanism outlined in previous subsections. We note that we do not exclude the possibility of subdominant amplification in reflectivity spectra arising from other phonons at different frequencies. In the supplementary material, we provide more details of our DFT analysis of IR-phonons in $\rm Ta_2 Ni Se_5$, the frozen phonon calculations and the calculation of optical matrix elements.

\subsubsection{Connection to the order parameter}


Finally, we discuss the effects of the phase transition on the phonon squeezing process. In figure~\ref{fig:DFTsummary}, we compute the electron-phonon interaction of the $4.7$ THz phonon in the low temperature monoclinic phase while for the high temperature orthorhombic phase we compute the electron-phonon interaction of the phonon adiabatically connected to the $4.7$ THz phonon eigenstate (for details on the identification of the adiabatically connected phonon see supplementary material). The electron-phonon interaction effectively disappears in the high temperature phase providing evidence that electron-phonon coupling is mediated by the order parameter. In particular, we suggest that close to the phase transition the electron-phonon coupling is proportional to the order parameter, $ \lambda_k = \Phi B_k$, where $\Phi$ is the order parameter and $B_k$ a constant. This, in turn, suggests that parametric amplification can be used as a nontrivial probe to investigate order parameter dynamics. The microscopic reason for the strong dependence of the electron-phonon coupling to the order parameter is a very interesting question that could reveal new insights about the nature of the order parameter in $\rm Ta_2 Ni Se_5$ and will be addressed in a subsequent publication. 





\section*{Discussion}

We have investigated the microscopic mechanism of amplification of THz optical reflectivity in $\rm{Ta_2NiSe_5}$ arising from high frequency optical pumping. We showed that strong electron-phonon coupling opens new pathways towards realizing THz parametric amplification through high frequency pumping. Ab-initio calculations highlight the importance of the $4.7$ THz IR-active phonon which is strongly coupled to electrons allowing for the amplification to manifest in the THz reflection spectrum.

Our theory indicates that choosing which electronic band to photoexcite selects the IR-phonon that is most strongly coupled to that electronic band. As a result, we can use different pumping frequencies in the same material to tune the frequency of the THz parametric amplification through mode-selective phonon squeezing.

Finally, we showed that the electron-phonon coupling is strongly dependent on the order parameter and becomes suppressed in the high temperature orthorhombic phase. This suggests that THz amplification of reflectivity can be used as a new probe to order parameter dynamics.

\section*{Materials and Methods}

\subsection*{Theoretical fit of the parametric reflectivity amplification in pumped $\rm Ta_2 Ni Se_5$} 
\label{app:TNSfit}

As mentioned in the main text in equation~(\ref{eq:indexRefl}), in principle the complex refractive index can be directly computed by the complex reflectivity amplitude, $r(\omega)$. However,  small phase errors upon experimental extraction of $r(\omega)$ could lead to unphysical behaviour of the reflectivity. To overcome this complication, we fit the data by assuming that for an insulator like $\rm Ta_2NiSe_5$ the refractive index is real. As a result, we can instead use the absolute value of $|r| = \sqrt{R}$ and express the refractive index as:
\begin{equation}
    n(\omega) = \frac{1 + \sqrt{R}}{1- \sqrt{R}}.
\end{equation}
Upon parametric resonance and solving equation~(\ref{eq:FloquetEigen}) together with the Floquet-Fresnel boundary conditions in equation~(\ref{eq:FloquetFresnel}), assuming a purely real spectrum with no dissipation leads to highly divergent behavior on parametric resonance. As a result, we put back dissipation by including an imaginary component to the refractive index. In the driven case we argue that this is physical, since $n_{\rm driven} (\omega)^2 = n^2(\omega) + i \sigma_{\rm driven}/\omega$ and dissipation can arise through through a transient contribution to the conductivity by electrons excited across the gap of the insulator. To fit the data we choose a parametric drive at $\omega_d = 9.4$ THz, drive amplitude $A_{\rm drive} = 7.5 \frac{THz^2}{c^2}$, overall constant renormalization of the refractive, $n^2_{\rm drive} = n^2(\omega) + \delta n^2 $, with $\delta n^2 = 0.1 + 0.01 i $ and an overall Gaussian broadening function with standard deviation of $0.1$ THz.

\subsection*{Electron - phonon interaction }
\label{app:SchriefferWolff}

The Hamiltonian of two electronic bands coupled by an allowed direct dipole transition that is, in turn, coupled to a phonon is given by:
\begin{equation}
    H_{\rm el}  = \frac{\Delta}{2} \left(\hat{c}_1^\dag \hat{c}_1 - \hat{c}_2^\dag \hat{c}_2 \right) + \lambda \left( \hat{c}_1^\dag \hat{c}_2+ \hat{c}_2^\dag \hat{c}_1 \right) Q,
\label{eq:elecHamil}
\end{equation}
where $\Delta = E_1 - E_2$ is the difference in energy between the two bands, $\lambda$ is the coupling constant coming from the dipole-dipole interaction in the dipole gauge, $Q$ is the phonon coordinate and $\{\hat{c}_1, \hat{c}_2\} $ are the annihilation operators of the two electron bands. The IR-phonon quadratic hamiltonian is given by:
\begin{equation*}
    H_{\rm ph} = Z E Q + M \omega_{\rm ph,0}^2 \frac{Q^2}{2} + \frac{\Pi^2}{2 M },
\label{eq:phononHamil}
\end{equation*}
where $\Pi$ is the conjugate momentum of the phonon coordinate $Q$, $Z$ the effective coupling to the electromagnetic field, E, and $\omega_{\rm ph,0} $ the phonon frequency. Since the transition itself is not resonant with the phonon mode, we decouple the linear electron phonon interaction perturbatively using a Schrieffer-Wolff transformation. This generates an interaction between the electron-hole pair density and the fluctuations of the phonon field which can be a resonant process. To perform this transformation, it is convenient to note that bilinear combinations of $\{ \hat{c}_1^\dag, \hat{c}_1, \hat{c}_2^\dag, \hat{c}_2 \}$ appearing in the Hamiltonian, obey SU(2) commutation relations by making the following identification:
\begin{subequations}
\begin{align}
    S^x = \frac{\hat{c}_1^\dag \hat{c}_2 + \hat{c}_2^\dag \hat{c}_1}{2}, \\
    S^y = -i \frac{\hat{c}_1^\dag \hat{c}_2 - \hat{c}_2^\dag \hat{c}_1}{2}, \\
    S^z = \frac{\hat{c}_1^\dag \hat{c}_1 - \hat{c}_2^\dag \hat{c}_2}{2},
\end{align}
\label{eq:spin}
\end{subequations}
where their commutators are $\left[S^x, S^y\right]  = i S^z$ and its cyclic permutations. In terms of the spin operators, we separate the Hamiltonian into a non-interacting and an  interacting part:
\begin{align}
    H_0 =& \Delta S^z + M \omega_{ph,0}^2 \frac{ Q^2}{2} + \frac{\Pi^2}{2 M }, \\
    V =& 2 \lambda S^x Q
\end{align}
To remove the interacting part $V$ to linear order in $\lambda$, we consider a unitary transformation of the type:
\begin{align}
    U =& {\rm Exp} \{ i A\}, \\
    A =& \alpha Q S^y + \beta \Pi S^x 
\end{align}
The Schrieffer-Wolff expansion is given by:
\begin{equation}
    U H_{total}U^\dag = H_0 + V - i \left[H_0, A \right] - i \left[V , A \right] - \frac{1}{2} \left[\left[H_0, A\right],A \right].
\end{equation}
The parameters $\alpha$ and $\beta$ are found such that:
\begin{equation}
    V = i \left[ H_0 , A\right]
\end{equation}
Matching linear terms in $\Pi$ and $Q$, leads to the parameters:
\begin{subequations}
\begin{align}
    &\beta =  \frac{\alpha}{M \Delta }, \\
    -&\Delta \alpha + 2 \lambda + M \omega^2_{\rm ph,0} \beta = 0, \\
    \Rightarrow &\alpha = \frac{2 \lambda}{\Delta \left( 1 - \frac{\omega_{\rm ph,0}^2}{\Delta^2} \right)},
\end{align}
\end{subequations}
which confirms that this perturbation theory can be carried out as long as $\omega_{\rm ph, 0}$ is off-resonant with the transition energy $\Delta$.  The effective electron-phonon interaction after the Schrieffer-Wolff transformation is given by:
\begin{equation}
\begin{split}
    H_{\rm eff} =& - i\left[V, A \right] - \frac{1}{2} \left[\left[H_0, A\right] ,A\right] = - \frac{i}{2} \left[V, A\right],\\
    =& \frac{2 \lambda^2}{\Delta \left(1 - \frac{\omega^2_{\rm ph,0}}{\Delta^2} \right)}Q^2 S^z + \frac{2 \lambda^2}{\Delta \left( 1 - \frac{\omega^2_{\rm ph,0} }{\Delta^2} \right)} \left(S^x\right)^2,
\end{split}
\end{equation}
where the second term does not depend on the phonons. The residual electron phonon interaction is thus give by:
\begin{equation}
    H_{\rm el-ph, eff} = \frac{\lambda^2}{\Delta \left( 1 - \frac{\omega_{\rm ph, 0}^2}{\Delta^2} \right)} Q^2 \left( n_1 - n_2 \right)
\end{equation}
where $n_1 = \hat{c}_1^\dag \hat{c}_1$ is the occupation number of mode 1. This implies that if either mode 1 or mode 2 are photo-excited the phonon can be squeezed. 

The linear term in the electric field, $Z E Q$, appearing in equation~\ref{eq:phononHamil}, is also transformed by the Schrieffer-Wolff transformation. It gives rise to a term $\beta Z E S^x$ which provides a small renormalization of the dipole transition amplitude between the bands and does not affect our discussion. 


The above result can be generalized to an arbitrary number of independent pairs of electronic states. For example, for two electronic states the above result in the limit of $\frac{\omega_{ph,0}^2}{\Delta_k^2} \rightarrow 0 $ is generalized to:
\begin{equation}
     H_{\rm el - ph, eff} = \sum_k \frac{\lambda_k^2}{\Delta_k } Q^2 \left( n_{1,k} - n_{2,k}\right),
\label{eq:elphonEff}
\end{equation}
This expression is the one used in the main text. 

\subsection*{Phonon squeezing}
\label{app:phononSqueeze}
Using a Hartree-Fock type approximation on the effective electron-phonon Hamiltonian in equation~(\ref{eq:elphonEff}), we derive an effective Hamiltonian for the phonon system only:
\begin{equation}
    H_{\rm ph} = Z E Q + M \left( \omega_{\rm ph, 0 }^2 + f(t) \right) \frac{Q^2}{2} + \frac{\Pi^2}{2 M } 
\label{eq:phononEff}
\end{equation}
where $Z$ is the IR activity of the phonon mode, $\omega_{\rm ph,0}$ the bare phonon frequency and the effective parametric drive, $f(t)$, is given by photo-excited electron density coupled to the phonon mode:
\begin{equation}
    f(t) = \sum_k \frac{2 \lambda_k^2}{M \Delta_k } \left( \expe{n_{1,k}} - \expe{n_{2,k} }\right),
\end{equation}
Due to the fast dynamics of electrons, $f(t)$ acts as an impulsive delta function like parametric drive. Such a drive is not periodic but it can linearly excite phonon fluctuations $\expe{Q^2}$ which will oscillate in time. To show this, we compute the equations of motion for fluctuations of the phonon field:
\begin{subequations}
\begin{align}
\partial_t \expe{Q^2} =& \frac{\expe{ \Pi Q+ Q \Pi}}{M} , \\
\begin{split}
    \partial_t \expe{ \Pi Q + Q \Pi}=& -2 M \left( \omega^2_{\rm ph,0} + f(t) \right) \expe{ Q^2} + \\& 2\frac{\expe{ \Pi^2}}{M} - 2 Z \expe{E Q} ,
\end{split}\\
\begin{split}
\partial_t \expe{\Pi ^2} =& - M \left( \omega^2_{\rm ph, 0} + f(t) \right)\expe{ \left( \Pi Q+ Q \Pi \right) } - \\&Z \expe{E \Pi + \Pi E} . 
\end{split}
\end{align}
\end{subequations}
At $k = 0$, we can use Maxwell's equations to remove the electric field dependence, $E = \frac{Z}{\epsilon \epsilon_0 } Q $. Performing this substitution simplifies the equations of motion,
\begin{subequations}
\begin{align}
\partial_t \expe{Q^2} =& \frac{\expe{ \Pi Q + Q \Pi}}{M} , \\
    \partial_t \expe{ \Pi Q + Q \Pi}=& -2 M \left( \omega^2_{\rm ph} + f(t) \right) \expe{ Q^2} +  2\frac{\expe{ \Pi^2}}{M} ,\\
\partial_t \expe{\Pi ^2} =& - M \left( \omega^2_{\rm ph} + f(t) \right) \expe{  \left( \Pi Q+ Q \Pi \right) }. 
\end{align}
\label{eq:fluceom}
\end{subequations}
where $\omega_{\rm ph}^2 = \omega^2_{\rm ph,0} + \omega_{\rm pl, phonon}^2$ is the frequency of the phonon at $k = 0$ which differs from the bare frequency by the phonon plasma frequency given by $\omega^2_{\rm pl, phonon} = \frac{Z^2}{\epsilon \epsilon_0 M }$. Being perturbative in the drive $f(t)$, we expand the phonon fluctuations,
\begin{equation}
    \expe{Q^2} = \expe{Q^2}_0 + \expe{Q^2}_1,
\end{equation}
where $\expe{Q^2}_0$ is the thermal expectation value and $\expe{Q^2}_1 \propto f(t)$. To linear order in $f(t)$, the equations of motion imply, $\expe{\Pi^2}_1/ M = M \omega^2_{\rm ph} \expe{Q^2}_1 + M f(t) \expe{Q^2}_0 + \mathcal{O}(f^2) $. Finally, combining equations in (\ref{eq:fluceom}) we find that:
\begin{equation}
   \left( \partial_t^2 + 4 \omega_{\rm ph}^2 \right) \expe{Q^2}_1 = - 4 f(t) \expe{Q^2}_0.
\end{equation}
This result shows that phonon fluctuations are linearly driven by photo-excitation and oscillate at twice the phonon frequency, $2 \omega_{\rm ph}$. These coherent oscillations of phonon fluctuations behave as a Raman mode.

\subsection*{Floquet matter from squeezing dynamics of phonons}
\label{app:IRphonon}

Material properties such as the electric permittivity, become time-periodic in the presence of oscillating fields through interactions. Here, we demonstrate how lattice potential anharmonicities lead to a time-periodic index of refraction:

We consider a phonon system with a $Q^4$ anharmonicity for the IR-phonon with a Hamiltonian:
\begin{equation}
     H_{\rm ph} = Z E Q + M \left( \omega_{\rm ph, 0 }^2 + f(t) \right) \frac{Q^2}{2} + \frac{\Pi^2}{2 M } + u Q^4
\label{eq:NonlinHamil}
\end{equation}
The equations of motion for the phonon given by the Hamiltonian in equation~(\ref{eq:NonlinHamil}) is
\begin{equation}
    \left( \partial_t^2 + \gamma \partial_t + \omega_{\rm ph,0}^2 + 4 uQ^2 \right) Q = \frac{Z}{M} E.
\end{equation}
Using a gaussian ansatz for the phonons, we can linearize the above equation as:
\begin{equation}
     \left( \partial_t^2 + \gamma \partial_t + \omega_{\rm ph,0}^2 + 12 u \expe{Q^2}(t) \right) Q = \frac{Z}{M} E,
\label{eq:phononeom}
\end{equation}
where the fluctuations $\expe{Q^2} = \expe{Q^2}_0 + A \cos( 2 \omega_{ph} t )$.
The phonon mode appears in the Maxwells equations as:
\begin{equation}
    \left( \frac{1}{c^2}\partial_t^2 - k^2 \right) E = - Z \partial^2_t Q.
\end{equation}
To find the effective signal idler mixing presented in equation~(\ref{eq:FloquetEigen}), we expand the equations of motion in signal and idler contributions:
\begin{equation}
    Q = Q_{s} e^{-i \omega_s t} + Q_{\rm id} e^{i \omega_{\rm id} t} 
\end{equation}
Equation~(\ref{eq:phononeom}) becomes:
\begin{equation}
\begin{split}
&\begin{pmatrix} Q_s \\ Q_{id} \end{pmatrix} = \begin{pmatrix} \frac{Z}{ \omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 } && 0 \\ 0 &&\frac{Z}{ \omega_{\rm id}^2 +i \gamma \omega_{\rm id} - \omega_{\rm ph}^2}  \end{pmatrix} \cdot \begin{pmatrix} E_{s} \\ E_{\rm id}  \end{pmatrix} + \\
&\frac{Z u A}{\left(\omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 \right)\left(\omega_{\rm id}^2 +i \gamma \omega_{\rm id} - \omega_{\rm ph}^2 \right) } \begin{pmatrix} E_{\rm id} \\ E_s \end{pmatrix} .
\end{split}
\label{eq:Floqphonon}
\end{equation}
Substituting equation~\ref{eq:Floqphonon} in Maxwells equation we find the equations of motion for the signal and idler component of the electric field to be:
\newline
\begin{subequations}
\begin{align}
\left( \frac{n^2_{eq.} (\omega_s) }{c^2} \omega_s^2 - k^2 \right) E_s+  A_{\rm drive, s} (\omega_s, \omega_{\rm id} ) E_{\rm id}=&0 ,\\
\left( \frac{n^2_{eq.} (\omega_{\rm id}) }{c^2} \omega_s^2 - k^2 \right) E_s+  A_{\rm drive, id} (\omega_s, \omega_{\rm id} ) E_{\rm id}=&0
\end{align}
\end{subequations}
where the signal and idler driving amplitudes, $A_{\rm drive, s}$ and $A_{\rm drive, id}$ are given by:
\begin{subequations}
\begin{align}
\begin{split}   
&A_{\rm drive,s} =\\
&\frac{Z^2 u  A \omega_s^2}{\left(\omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 \right)\left(\omega_{\rm id}^2 +i \gamma \omega_{\rm id} - \omega_{\rm ph}^2 \right) } ,
\end{split} \\
\begin{split} 
     &A_{\rm drive,id} = \\
     &\frac{Z^2 u  A \omega_{\rm id}^2}{\left(\omega_s^2 + i \gamma \omega_s - \omega_{\rm ph}^2 \right)\left(\omega_{\rm id}^2 +i \gamma \omega_{\rm id}- \omega_{\rm ph}^2 \right) } .
\end{split}
\end{align}
\end{subequations}


\subsection*{DFT calculations of IR-phonons in $\rm Ta_2 Ni Se_5$}
The phonon spectrum at $\Gamma$ has been computed using the Density Functional Perturbation Theory routines of the VASP code \cite{Vasp1,Vasp2,Vasp3,Vasp4} using the vdW-opt88 functional \cite{Klime2011,Klimes2010} on a 48x4x6 mesh. The Born effective charges have been computed on a similar k-mesh, the IR-activity has been computed using phonopy \cite{phonopy} and the phonon spectroscopy code by Skelton et. al \cite{IR_activity}. The electronic bandstructure plots have been obtained using a 24x4x6 k-mesh and a 320 eV cutoff using the standard PBE functional. 

\subsubsection{Phonons}
For the monoclinic geometry there are 21 symmetry allowed IR-active modes, 11 B$_\text{u}$ and 10 A$_\text{u}$ modes.Of these modes, only four are in a frequency range near 4.5THz, listed in table~\ref{tab: IR_phonons}.

\begin{table}[]
    \centering
    \begin{tabular}{c|c|c}
    Phonon & Symmetry & Frequency (THz) \\ \hline
    21  &  B$_\text{u}$  &  4.238  \\
    22  &  A$_\text{u}$  &  4.314  \\
    25  &  A$_\text{u}$  &  4.634  \\
    26  &  B$_\text{u}$  &  4.699  \\
    \end{tabular}
    \caption{IR-active phonon modes with a frequency close to 4.5 THz. }
    \label{tab: IR_phonons}
\end{table}

In the main text we identify phonon 26 which has a frequency of 4.699 THz as the relevant for phonon for parametric amplification. Therefore, we have investigated how this phonon evolves across the structural phase transition. To do this we have linearly interpolated between the monoclinic and high temperature orthorhombic atomic and lattice geometry using an interpolation parameter d. 

\begin{equation}
\begin{split}
    \vec{v}_\text{i}(d) &= \vec{v}_\text{o} + d \cdot \vec{v}_\text{t} \\
    \vec{v}_\text{lat,i}(d) &= \vec{v}_\text{lat,o} + d \cdot \vec{v}_\text{lat,t} \quad , \\
\end{split}
\end{equation}
with the transition vector for the atomic configuration and the lattice defined as
\begin{equation}
\begin{split}
\vec{v}_\text{t} &= \vec{v}_{\text{m}} - \vec{v}_{\text{o}} \\
\vec{v}_\text{lat,t} &= \vec{v}_\text{lat,m} - \vec{v}_\text{lat,o} \quad .\\
\end{split}
\end{equation}
The subscripts m denotes the monoclinic geometry and o the orthorhombic geometry.
Because the phonon eigenmodes form an orthonormal basis for the movement of all atoms we can trace how the monoclinic 4.7 THz phonon hybridizes along the phase transition into the orthorhombic phonons, by calculating the overlap of the monoclinic phonon eigenvector with the phonon eigenvectors along the transition. The result is shown in figure \ref{fig: IR_phonon_overlap}. It shows that the 4.7 THz phonon predominantly disperses into two orthorhombic phonons at 4.38 THz and 4.56 THz. 
\begin{figure}
    \centering
    \includegraphics[scale=0.42]{fatbands_IRmode26.pdf}
    \caption{Overlap of the 4.7 THz monoclinic phonon with the phonons calculated along the orthorhombic to monoclinic phase transition. The dot size is proportional to the overlap. The figure shows that the 4.7 THz phonon predominantly hybridizes into two orthorhombic phonons at frequencies of 4.38 THz and 4.56 THz. Note that the orthorhombic structure exhibits two phonon instabilities with complex frequency which are displayed as negative energies here.}
    \label{fig: IR_phonon_overlap}
\end{figure}

\subsubsection{Frozen phonon calculations}
To estimate the electron phonon coupling of the four IR-active phonon modes displayed in table \ref{tab: IR_phonons}, we show frozen phonon band structure calculations for these modes. As discussed in the main text, only the strongly amplified IR mode should display a considerable coupling to the electronic bands.
Upon displacing the atomic positions along the eigendisplacement of the four phonon modes and recalculating the electronic bandstructure, only mode 26 shows a strong modulation of the conduction bands, indicating a strong coupling of this mode to the electronic degrees of freedom upon excitation (see Fig. \ref{fig: IR_phonons}). 


\begin{figure}
    \centering
    \includegraphics[scale=0.42]{PBE_bands.pdf}
    \caption{Frozen phonon bandstructure of all possible IR active active phonons in the monoclinic phase between 4.2 THz and 4.7 THz using the PBE functional. Phonon 26 shows she strongest coupling to the bandstructure upon displacement along its eigenmode. Most importantly, the second and third conduction band are modulated by this mode and are displaced in opposite directions.}
    \label{fig: IR_phonons}
\end{figure}

To quantify the strength of the electron-phonon coupling $\lambda_i$,  we define 

\begin{equation}\label{EQ: e-ph coupling}
    \lambda_i = \sum_k  w_k | \,\, (E_{i,k})\bigr\rvert_{\langle Q \rangle = 0} - (E_{i,k})\bigr\rvert_{\langle Q \rangle \neq 0} \,\, | ,
\end{equation}
as the difference of the Kohn-Sham energies at the atomic equilibrium position $\langle Q \rangle = 0$ and the frozen phonon displaced position $\langle Q \rangle \neq 0$. k runs over all k-Points along the chosen Z-G path near the bandedge and $w_k$ is the symmetry weight of the k-Point. Computing this coupling for the band structures shown in figure \ref{fig: IR_phonons} one obtains for the first four conduction bands the values presented in table \ref{tab: IR_couplings}. It shows that the strongest coupling of the conduction bands is to phonon mode 26. For this mode the second and third conduction band are shifted into opposite directions along the Z-$\Gamma$ near the band edge which is consistent with the above proposed two band model. Therefore, we identify this mode as the relevant mode for the parametric amplification process. 


\begin{table*}[]
    \centering
    \begin{tabular}{c|c|c|c|c}
         & Phonon21 (meV) & Phonon22 (meV) & Phonon25 (meV) & Phonon26 (meV) \\ \hline
           &         &         &        &         \\
cond. Band 1 &  1.197 &       2.762 &       3.321 &       8.171 \\
cond. Band 2 &  1.568 &       2.378 &       1.324 &      24.953 \\
cond. Band 3 &  5.363 &       3.391 &       6.473 &      34.436 \\
cond. Band 4 &  1.329 &       5.657 &      17.629 &      11.814 \\
    \end{tabular}
    \caption{Coupling as defined by equation \eqref{EQ: e-ph coupling} for the bandstructure presented in figure \ref{fig: IR_phonons} . Only states near the bandedge along Z-$\Gamma$ have been considered for the summation in equation \ref{EQ: e-ph coupling}. One can see that conduction bands 2 and 3 shift in opposite directions upon displacing along phonon mode 26 (see table \ref{tab: IR_phonons}). This is consistent with the parametric driving mechanism presented in the main text.}
    \label{tab: IR_couplings}
\end{table*}

\subsubsection{Optical transitions}
\begin{figure}
    \centering
    \includegraphics[scale=0.35]{PBE_fatbands_equ_05ev_allDir.pdf}
    \caption{Calculation of the dipole matrix elements for all three cartesian directions. We filtered only transitions in the experimentally relevant energy range between 0.33 eV and 0.9 eV. The size of the dots shows the size of the dipole matrix element at a given k-Point and band.}
    \label{fig: opticalTransitions}
\end{figure}
In the experiment, parametric amplification has been pumped with a 0.5 eV laser. In the above section we have identified a phonon mode which couples strongly to the second and third conduction band and might be responsible for the parametric amplification process. To fully justify this claim we need to show that using the 0.5 eV pump electrons are indeed excited into these strongly coupled conduction band states. Here we have computed the optical transition matrix elements and plotted them as dots on the bandstructure, where the dot size for the valence and conduction bands corresponds to 
\begin{equation}
\begin{split}
    M_{\lambda, v,k} &= \sum_c M_{\lambda,v,c,k} \\
    M_{\lambda, c,k} &= \sum_v M_{\lambda,v,c,k} \\
\end{split}
\end{equation}
with $M_{\lambda,v,c,k} = |\langle \phi_{v,k} | (-i) \frac{d}{d k_\lambda} | \phi_{c,k} \rangle|^2 $.
$k$ labels the k-point, $v$ the valence, $c$ the conduction band involved in the optical transition and $\lambda$ the Cartesian directions. \newline
The result for all three Cartesian directions is shown in figure \ref{fig: opticalTransitions}. We have have filtered only the experimentally relevant transitions between 0.33 eV and 0.9 eV transition energy. Such a large energy window should compensate for any underestimation of the electronic gap given by the PBE functional. One sees that optical transitions in x-directions are strongest for the first 3 conduction band states. Thus, we can conclude that upon pumping with 0.5 eV there are many electrons excited into the electronic states that couple strongly to phonon 26 and can lead to the parametric amplification process.



\label{app:DFT}

