section0em3em subsection2em2em

# Supplemental material: Two-dimensional THz spectroscopy of nonlinear
phononics in the topological insulator MnBi2Te4

T.G.H. Blank  Radboud University, Institute for Molecules and Materials, 6525
AJ Nijmegen, The Netherlands.  K.A. Grishunin  Radboud University, Institute
for Molecules and Materials, 6525 AJ Nijmegen, The Netherlands.  K.A. Zvezdin
Prokhorov General Physics Institute of the Russian Academy of Sciences, 119991
Moscow, Russia.  N.T. Hai  Department of Physics, National Changhua University
of Education, Changhua 500, Taiwan.  J.C. Wu  Department of Physics, National
Changhua University of Education, Changhua 500, Taiwan.  S.-H. Su  Department
of Physics, National Cheng Kung University, Tainan 701, Taiwan.  J.-C. A.
Huang  Department of Physics, National Cheng Kung University, Tainan 701,
Taiwan.  A.K. Zvezdin  Prokhorov General Physics Institute of the Russian
Academy of Sciences, 119991 Moscow, Russia.  A.V. Kimel  Radboud University,
Institute for Molecules and Materials, 6525 AJ Nijmegen, The Netherlands.

###### Contents

  1. S1 THz electric field calibration by electro-optic sampling in GaP
  2. S2 Supplementary results of single THz pump-probe measurements
    1. I Temperature dependence
    2. II THz electric field strength dependence
    3. III Polarization dependencies of the THz pump ($\alpha$) and probe ($\beta$)
  3. S3 Group theory of the point group $D_{3d}$
    1. I Detection mechanism
    2. II Symmetry-allowed coupling terms
  4. S4 Details of 2D THz spectroscopy
  5. S5 Analytical solution to the equations of motion
  6. S6 Numerical simulations of phonon dynamics
    1. I THz-mediated coupling between the modes $\propto Q_{\mathrm{R}}Q_{\mathrm{IR}}E$
    2. II Photonic impact $\propto E^{2}$
    3. III Anharmonic lattice dynamics $\propto Q_{\mathrm{R}}Q_{\mathrm{IR}}^{2}$
    4. IV Distinguishing phonon and complementing peaks

##  S1 THz electric field calibration by electro-optic sampling in GaP

Figure S1: THz calibration using electro-optical sampling in GaP. The inset
shows the associated THz spectrum.

After the generation using tilted pulse front optical rectification in a
LiNbO3 crystal, the THz beam was expanded, collimated, and focused using three
off-axis parabolic mirrors onto a ($110$) GaP crystal of $50$ $\mu$m thickness
that was grown on a ($100$)-cut GaP substrate with $2$ mm thickness. The THz
electric field induces birefringence in the ($110$) GaP [1, 2], which is
mapped as a function of time by measuring the changes of ellipticity of a
linearly polarized probe pulse a central wavelength of $800$ nm in a balanced
detection scheme. From the electro-optic coefficients, we calculated the THz
electric field strength as in Fig. S1 resulting in a peak electric field
strength of $780$ kV/cm.

##  S2 Supplementary results of single THz pump-probe measurements

The THz-induced dynamics shown in main the article have been characterized as
a function of temperature $T$, the strength of the THz electric field, and the
polarization of the pump and probe. Some of these dependencies are already
shown in the main article. In this section, all the supplementary results of
these dependencies are presented.

###  I Temperature dependence

Figure S2 shows the spectra of the waveforms for several temperatures,
revealing a minor red-shift of the phonon frequency.

Figure S2: FFT spectrum of the THz-induced dynamics for several temperatures.
The $E_{g}$ phonon mode frequency slightly decreases for increasing
temperatures, in agreement with previous results [3]. Moreover, the signal at
$T=78$ K is strong, making it suitable for the measurement using 2D THz
spectroscopy.

###  II THz electric field strength dependence

Most of the following results have already been plotted in the main article.
However, we also noted a small contribution in the signal spectrum centered
around $2.07$ THz, which we attribute to the THz-induced birefringence at the
overlap (i.e. Pockels effect). Indeed, this contribution scaled exactly linear
with the applied field (see Figs. S3(a) and S4(a)). Secondly, we got a strong
DC component in the spectrum due to THz-induced offset. This offset scales
with a rather unusual power-law dependence $\propto E_{\mathrm{THz}}^{2.6}$.

Adding to this, besides measuring changes in the rotation of the probe, we
also tracked changes in transmission. In all measurements, a step-like change
in transmission was observed right when the THz pulse arrived. The relative
change in transmission $\Delta T$ scales with exactly the same anomalous
power-law dependence $\propto E_{\mathrm{THz}}^{2.6}$ as the THz-induced DC
offset (see Figs. S3(b) and S4(b)), indicating a relation between the two.

Figure S3: (a) The FFT peak amplitude of the THz-induced rotation as a
function of THz peak electric field strength (maximally $780$ kV/cm), measured
at a temperature of $T=12$ K. (b) Power-law dependence of the transient change
of transmission, which is shown to have the same unusual power-law dependence
$\gamma=2.6$ as the DC component measured in the rotation data. Figure S4:
Same measurement as described in Fig. S3, but measured $T=78$ K. The power-law
dependencies w.r.t. $T=12$ K are unchanged.

###  III Polarization dependencies of the THz pump ($\alpha$) and probe
($\beta$)

The orientation of the THz pulse was controlled by a pair of wire-grid
polarizers, and the probe polarization was rotated using a half-wave plate.
The THz electric field orientation is characterized by the angle $\alpha$ with
respect to the experimental $y$-axis
$\mathbf{E}_{\mathrm{THz}}=E_{0}(\sin\alpha,\cos\alpha,0)$, while the angle of
the probe pulse electric field $\beta$ is defined with respect to the
experimental $x$-axis
$\mathbf{E}_{\mathrm{probe}}=E_{\mathrm{probe}}(\cos\alpha,\sin\alpha)$. We
have no information about how these experimental axes compare to the crystal
$x$ and $y$-axes. However, we do know that the pulses travel along the
crystallographic $z$-axis, so the experimental $x-y$ plane does coincide with
the crystallographic $x-y$ plane apart from an in-plane rotation. In all
measurements presented in the main article, both $\alpha$ and $\beta$ were put
to $0^{\circ}$, as these parameters gave strong results (see Figs. S5, S6 and
S7) without needing to rotate the THz pulse which is accompanied by a loss in
power as discussed below.

Rotating the THz polarization using two wire-grid polarizers reduces the peak
THz electric field. In particular, the maximum reduction occurs when rotating
to $\alpha=90^{\circ}$, where the field is halved. Of course, during the scans
as a function of THz polarization $\alpha$, we kept the field constant by
picking the right angles of the wire-grid polarizers.

(a) $\beta=0^{\circ}$.

(b) $\beta=45^{\circ}$.

(c) $\beta=90^{\circ}$.

Figure S5: Pump polarization scans for several probe polarizations
$\beta=0,45$ and $90$ degrees, measured at $T=78$ K. The resemblance of (a)
and (c) underlines a $\pi/2$ periodicity of the probe polarization, which will
be confirmed in Fig. S7.

Figure S6: Real FFT amplitudes of the data shown in Fig. S5(a), exhibiting the
$\pi$ periodicity with respect to the pump polarization. It shows that the
signal at $\alpha=0^{\circ}$ is relatively strong, therefore it was used
throughout in the main article to avoid rotating the THz polarization (which
is accompanied by a loss in the field).

Figure S7: (a) Different waveforms for several $\beta$ with $\alpha=0^{\circ}$
fixed and $T=78$ K. (b) Real FFT amplitude at $3.14$ THz as a function of
probe polarization $\beta$.

##  S3 Group theory of the point group $D_{3d}$

Table 1: Character table for point-group $D_{3d}$ (rhombohedral) [4]. | $E$ | $2C_{3}$ | $3C^{\prime}_{2}$ | $i$ | $2iC_{3}$ | $3iC^{\prime}_{2}$ |  |
---|---|---|---|---|---|---|---|---
$A_{1g}$ | 1 | 1 | 1 | 1 | 1 | 1 |  | $x^{2}+y^{2},z^{2}$
$A_{2g}$ | 1 | 1 | -1 | 1 | 1 | -1 | $R_{z}$ |
$E_{g}$ | 2 | -1 | 0 | 2 | -1 | 0 | $(R_{x},R_{y})$ | $(xz,yz),(x^{2}-y^{2},xy)$
$A_{1u}$ | 1 | 1 | 1 | -1 | -1 | -1 |  |
$A_{2u}$ | 1 | 1 | -1 | -1 | -1 | 1 | $z$ |
$E_{u}$ | 2 | -1 | 0 | -2 | 1 | 0 | $(x,y)$ |

From the basis functions in the character table associated with each
irreducible representation, it can be seen that $A_{1g}$ and $E_{g}$ modes are
Raman-active, while $A_{2u}$ and $E_{u}$ are IR-active. Given that MnBi2Te4
has $7$ atoms in the unit cell, group theory predicts the existence of nine
Raman-active $3A_{1g}\oplus 3E_{g}$ and nine IR-active optical phonons
$3A_{2u}\oplus 3E_{u}$ [5, 6].

###  I Detection mechanism

Regarding the detection via the depolarized scattering mechanism, a Raman-
active phonon may scatter light with a rotated polarization, if it transforms
as any of the off-diagonal components of the polarizability tensor [7]. Given
that the polarizability tensor is a symmetric tensor, these components must
transform as $xy,xz$, and $yz$. Therefore, from the character table, it can be
seen that only the $E_{g}$ mode contributes to the depolarized scattering
mechanism. The scattering efficiency $S_{ij}$, and thereby the degree of
ellipticity induced on the light, of a Raman-active mode with symmetry
$\Delta$ can be expressed in terms of its Raman tensor [8]:

$S_{ij}\propto(\mathbf{e}_{i}\cdot\mathcal{R}(\Delta)\cdot\mathbf{e}_{j})^{2}$
(S1)

where $\mathcal{R}(\Delta)$ the Raman tensor of the $\Delta$ phonon and
$\mathbf{e}_{i}$ and $\mathbf{e}_{j}$ are the unit polarization vectors of the
incident and scattered light, respectively. For the two-fold degenerate
$E_{g}$ mode there are two Raman tensors, and the total intensity can be found
by adding the scattering efficiencies of them together [8]:

$\mathcal{R}(E_{g}^{(I)})=\begin{pmatrix}c&0&0\\\ 0&-c&d\\\
0&d&0\end{pmatrix}\ \ \ \ \ \ \ \ \ \
\mathcal{R}(E_{g}^{(II)})=\begin{pmatrix}0&-c&-d\\\ -c&0&0\\\
-d&0&0\end{pmatrix}$ (S2)

where the $c,d$ are constants. On the contrary, the Raman tensor of $A_{1g}$
contains only diagonal components [8].

However, in the experiment, we measured rotation and not ellipticity. We
believe that the ellipticity induced by the depolarized scattering mechanisms
is converted to a rotation due to static birefringence in the sample and
substrate.

###  II Symmetry-allowed coupling terms

Table 2: Direct product table of the irreducible representations [4]. $\bigotimes$ | $A_{1g}$ | $A_{1u}$ | $A_{2g}$ | $A_{2u}$ | $E_{u}$ | $E_{g}$
---|---|---|---|---|---|---
$A_{1g}$ | $A_{1g}$ | $A_{1u}$ | $A_{2g}$ | $A_{2u}$ | $E_{u}$ | $E_{g}$
$A_{1u}$ | - | $A_{1g}$ | $A_{2u}$ | $A_{2g}$ | $E_{g}$ | $E_{u}$
$A_{2g}$ | - | - | $A_{1g}$ | $A_{1u}$ | $E_{u}$ | $E_{g}$
$A_{2u}$ | - | - | - | $A_{1g}$ | $E_{g}$ | $E_{u}$
$E_{u}$ | - | - | - | - | $A_{1g}\oplus A_{2g}\oplus E_{g}$ | $A_{1u}\oplus A_{2u}\oplus E_{u}$
$E_{g}$ | - | - | - | - | - | $A_{1g}\oplus A_{2g}\oplus E_{g}$

When considering which of the coupling terms in the interaction potential were
allowed, we used the direct product table of the point group $D_{3d}$ (Table
2). As the potential energy $U$ is scalar and invariant under point-group
operations, it transforms as the totally symmetric representation $A_{1g}$.
Therefore, the formal selection rule for coupling terms to exist, is that the
product of their irreducible representations contains the $A_{1g}$
representation [7, 9]. Here, it is important to note that the THz electric
field lies in the $xy$ plane and therefore it transforms as $E_{u}$. For
example, the nonlinear coupling term described by the interaction potential
$\propto Q_{\mathrm{R}}Q_{\mathrm{IR}}^{2}$ involving the Raman-active $E_{g}$
phonon and an IR-active $A_{2u}$ is not allowed, since $E_{g}\otimes
A_{2u}\otimes A_{2u}=E_{g}$ does not contain the totally symmetric
representation. Moreover, we note that the linear excitation (electric dipole
interaction) of the IR-active $A_{2u}$ described by the potential $\propto
Q_{\mathrm{IR}}E$ is not allowed for fields in the $x-y$ plane (in-plane
fields transform as $E_{u}$), and therefore we do not consider phonons of this
species in the main article. It is only allowed when the field is along the
$z$-axis (in which case the field transforms as $A_{2u}$). This fact is also
intuitive, as the IR-active $A_{2u}$ mode corresponds to out-of-plane
vibrations [10, 5] and can therefore only be excited by an out-of-plane field.

Regarding the terms involving the Raman-active $E_{g}$ phonon and the IR-
active $E_{u}$ phonon, all terms discussed in the main article are allowed by
symmetry.

##  S4 Details of 2D THz spectroscopy

Figure S8: (a) Calibrated waveforms of the two THz pump pulses used for 2D
spectroscopy, measured by electro-optical sampling in GaP. (b) Associated
normalized Fourier spectra of the two THz pulses. Although the two THz
waveforms in the time domain are a bit different, their absolute Fourier
spectra are nearly identical.

In the double-pump experiment, we generated two THz pulses with peak electric
fields of approximately $350$ and $330$ kV/cm (see Fig. S8) in a single LiNbO3
crystal from two separate optical pulses with a central wavelength of $800$
nm, whose mutual delay (the excitation time $\tau$) can be adjusted. We
modulated the THz pump $1$ and pump $2$ at unequal repetition rates of $500$
and $250$ Hz, respectively, by mechanical choppers. Thereby, we could isolate
the signal from pump $1$ ($\theta_{1}(t,\tau)$) and pump $2$ ($\theta_{2}(t))$
separately with a data-acquisition card. However, due to the non-linear nature
of the generation process by optical rectification, the simultaneous presence
of the two pulses inside the LiNbO3 crystal ($\tau\approx 0$) results in a
parasitic mutual interaction between the pulses, which affects the generation
process. This causes anomalous results around the pump-pump overlap $\tau=0$,
as superposition is broken by the generation process itself. Therefore, during
the analysis of the results (when taking the 2D FFT of the data) presented in
the main article, we excluded a small region around the overlap of the two THz
pulses $|\tau|<1$ ps, in order to exclude any artificial signals.

##  S5 Analytical solution to the equations of motion

The THz electric field can be modeled accurately by the Gaussian derivative
function:

$E(t)=t_{0}\frac{d}{dt}\mathcal{G}(t),\qquad\mathcal{G}(t)=E_{0}\sin\omega_{0}t\exp\left(-\frac{t^{2}}{t_{0}^{2}}\right)$
(S3)

where $\omega_{0}t_{0}\approx 1$ and $E_{0}$ is the peak electric field
strength. This electric field is the direct driving force of the IR-active
phonon modes described by the normal coordinate $Q_{\eta}$, whose equations of
motion in absence of dissipation are of the form:

$\ddot{Q}_{\eta}+\omega_{\eta}^{2}Q_{\eta}=\frac{Z_{\eta}^{*}E(t)}{M_{\eta}}$
(S4)

In order to solve this equation, we first substitute
$\xi_{\eta}=\dot{Q}_{\eta}+i\omega_{\eta}Q_{\eta}$ and solve the resulting
first order differential equation for $\xi_{\eta}$. Then we can obtain our
solution $Q_{\eta}=\frac{1}{\omega_{\eta}}\Im{\xi_{\eta}}$, which gives:

$\displaystyle Q_{\eta}(t)\rvert_{t\gg t_{0}}$ $\displaystyle=$
$\displaystyle\frac{Z_{\eta}^{*}}{\omega_{\eta}M_{\eta}}\Im{\exp(i\omega_{\eta}t)\int_{-\infty}^{\infty}E(t)\exp(-i\omega_{\eta}t)\,\mathrm{d}t}$
(S5)

after integration by parts:

$\displaystyle Q_{\eta}(t)\rvert_{t\gg t_{0}}$ $\displaystyle=$
$\displaystyle\frac{Z_{\eta}^{*}E_{0}t_{0}}{M_{\eta}}\Im{i\exp(i\omega_{\eta}t)\int_{-\infty}^{\infty}\sin\omega_{0}t\exp(-t^{2}/t^{2}_{0})\exp(-i\omega_{\eta}t)\,\mathrm{d}t}$
$\displaystyle=$
$\displaystyle\frac{Z_{\eta}^{*}E_{0}t_{0}}{2M_{\eta}}\Im{\exp(i\omega_{\eta}t)\int_{-\infty}^{\infty}\left[\cos(\omega_{0}-\omega_{\eta})t-\cos(\omega_{0}+\omega_{\eta})t\right]\exp(-t^{2}/t^{2}_{0})\,\mathrm{d}t}$

where we used that thy symmetric integration over an odd function around zero
is zero. After the substitution $t^{\prime}\mapsto\sqrt{2}t/t_{0}$, we can
apply the identity
$\int_{-\infty}^{\infty}\cos(ax)\exp(-\tfrac{1}{2}x^{2})\,\mathrm{d}x=\sqrt{2\pi}\exp(-\tfrac{1}{2}a^{2})$
and obtain after rearranging the terms:

$\displaystyle Q_{\eta}(t)\rvert_{t\gg t_{0}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2\pi}Z_{\eta}^{*}E_{0}t^{2}_{0}}{2M_{\eta}}\left[\exp(-\tfrac{1}{4}(\omega_{0}-\omega_{\eta})^{2}t_{0}^{2})-\exp(-\tfrac{1}{4}(\omega_{0}+\omega_{\eta})^{2}t_{0}^{2})\right]\sin\omega_{\eta}t$
(S7) $\displaystyle=$
$\displaystyle\frac{Z_{\eta}^{*}E_{0}t^{2}_{0}}{M_{\eta}}\Xi(\omega_{\eta})\sin\omega_{\eta}t$

where
$\Xi(\omega)\equiv\sqrt{\pi}\exp(-\tfrac{1}{4}(1+t_{0}^{2}\omega^{2}))\sinh(\tfrac{1}{2}t_{0}\omega)$.
For later purpose we define
$K_{\eta}=\frac{Z_{\eta}^{*}}{M_{\eta}}\Xi(\omega_{\eta})$. Now, we are
interested in the nonlinear excitation of the $E_{g}$ Raman-active phonon
through these IR-active modes. When an IR-active phonon with normal coordinate
$Q_{\eta}$ is excited by the first THz pulse (electric field $E(t)$), the
second THz pulse (electric field $E^{\prime}(t,\tau)$) delayed at a time
$\tau$ may excite the Raman-active phonon through nonlinear coupling term
$\propto Q_{\mathrm{R}}Q_{\mathrm{IR}}E$. The equation of motion for the
Raman-active phonon (with normal coordinate $Q_{R}$), taking into account this
nonlinear coupling term, is given by:

$\ddot{Q}_{\mathrm{R}}+\omega_{\mathrm{R}}^{2}Q_{\mathrm{R}}=\frac{1}{M_{\mathrm{R}}}\sum_{\eta=1}^{4}b_{\eta}Q_{\mathrm{IR}\eta}(t)E^{\prime}(t,\tau)$
(S8)

where the $b_{\eta}$ coupling constants. The electric field of the second THz
pulse with a peak amplitude of $E_{0}^{\prime}$ is similar to the first THz
pulse, but delayed at a time $\tau>0$:

$\frac{E^{\prime}(t,\tau)}{E^{\prime}_{0}}=\frac{E(t-\tau)}{E_{0}}$ (S9)

This second pulse will efficiently excite the Raman-active mode via the
already excited IR-active mode of the first pulse. Note that for simplicity,
we ignore the occurrence of this process within a single THz pulse. The
solution for $t\gg\tau$ of the Raman-active mode normal coordinate in Eq. S8
is given by:

$\displaystyle Q_{\mathrm{R}}(t)\rvert_{t\gg\tau}$
$\displaystyle=\frac{1}{\omega_{\mathrm{R}}M_{\mathrm{R}}}\sum_{\eta=1}^{4}b_{\eta}\Im{e^{i\omega_{\mathrm{R}}t}\int_{-\infty}^{\infty}Q_{\mathrm{IR}\eta}(t)E^{\prime}(t)e^{-i\omega_{\mathrm{R}}t}\,\mathrm{d}t}$
(S10)
$\displaystyle=\frac{E_{0}E^{\prime}_{0}t_{0}^{3}}{\omega_{\mathrm{R}}M_{\mathrm{R}}}\sum_{\eta=1}^{4}b_{\eta}K_{\eta}\Im{e^{i\omega_{\mathrm{R}}t}\int_{-\infty}^{\infty}\sin\omega_{\mu}t\frac{\mathrm{d}}{\mathrm{d}t}\left[\sin\omega_{0}(t-\tau)\exp(-\frac{(t-\tau)^{2}}{t_{0}^{2}})\right]e^{-i\omega_{\mathrm{R}}t}\,\mathrm{d}t}$
$\displaystyle=\frac{E_{0}E^{\prime}_{0}t_{0}^{3}}{2\omega_{\mathrm{R}}M_{\mathrm{R}}}\sum_{\eta=1}^{4}b_{\eta}K_{\eta}\Bigg{[}\Im{-(\omega_{\mathrm{R}}+\omega_{\eta})e^{i\omega_{\mathrm{R}}t}e^{-i(\omega_{\mathrm{R}}+\omega_{\eta})\tau}\int_{-\infty}^{\infty}\sin\omega_{0}te^{-\frac{t^{2}}{t_{0}^{2}}}e^{-i(\omega_{\mathrm{R}}+\omega_{\eta})t}\,\mathrm{d}t}$
$\displaystyle\qquad\qquad\qquad\qquad+\Im{(\omega_{\mathrm{R}}-\omega_{\eta})e^{i\omega_{\mathrm{R}}t}e^{-i(\omega_{\mathrm{R}}-\omega_{\eta})\tau}\int_{-\infty}^{\infty}\sin\omega_{0}te^{-\frac{t^{2}}{t_{0}^{2}}}e^{-i(\omega_{\mathrm{R}}-\omega_{\eta})t}\,\mathrm{d}t}\Bigg{]}$
$\displaystyle=E_{0}E^{\prime}_{0}t_{0}^{4}\sum_{\eta=1}^{4}b_{\eta}K_{\eta}\Bigg{(}L^{+}_{R\eta}\cos(\omega_{\mathrm{R}}t-(\omega_{\mathrm{R}}+\omega_{\eta})\tau)-L^{-}_{R\eta}\cos(\omega_{\mathrm{R}}t-(\omega_{\mathrm{R}}-\omega_{\eta})\tau)\Bigg{)}$
(S11)

where
$L^{\pm}_{R\eta}=\frac{\omega_{\mathrm{R}}\pm\omega_{\eta}}{\omega_{\mathrm{R}}M_{\mathrm{R}}}\Xi(\omega_{R}\pm\omega_{\eta})$.
Therefore, we will observe peaks in the 2D FFT spectrum at the complementing
excitation frequencies
$(\omega_{\mathrm{det}},\omega_{\mathrm{ex}})=(\omega_{\mathrm{R}},\omega_{\mathrm{R}}\pm\omega_{\eta})$.

To explain the peak seen at the excitation frequency of the IR-active phonon
frequency, i.e. at
$(\omega_{\mathrm{det}},\omega_{\mathrm{ex}})=(\omega_{\mathrm{R}},\omega_{\eta})$,
we have to look at the region where moving THz pulse $E^{\prime}(t)$ arrives
before the ”first” THz pulse $E(t)$ (when $\tau<0$). Therefore, we first need
to solve the excitation of the IR-active mode with pulse $E^{\prime}(t-\tau)$
with $\tau<0$ as in Eq. (S9), which gives:

$Q_{\eta}(t)\rvert_{t>\tau,t<0}=\frac{Z_{\eta}^{*}E^{\prime}_{0}t^{2}_{0}}{M_{\eta}}\Xi(\omega_{\eta})\sin\omega_{\eta}(t-\tau).$
(S12)

Afterward, we again solve Eq. (S8) with $E(t)$ as in Eq. (S3), in the presence
of an IR-active phonon described by Eq. S12. The solution is given by
substituting $t\mapsto t+\tau$ in the solution of Eq. (S11):

$Q_{\mathrm{R}}(t)\rvert_{t>\tau,t>0}=E_{0}E^{\prime}_{0}t_{0}^{4}\sum_{\eta=1}^{4}b_{\eta}K_{\eta}\Bigg{(}L^{+}_{R\eta}\cos(\omega_{\mathrm{R}}t-\omega_{\eta}\tau)-L^{-}_{R\eta}\cos(\omega_{\mathrm{R}}t+\omega_{\eta}\tau)\Bigg{)}.$
(S13)

It means that depending on whether the moving pulse comes first or second, one
can distinguish between the real and complementing phonon peak. If one takes
the 2D FFT where the moving pulse comes first $\tau<0$, a real phonon peak
will appear in the spectrum. And if the Fourier transform is taken over the
region where the moving pulse comes after THz pulse at $t=0$ (i.e. $\tau>0$),
the complementary phonon peak $\omega_{\mathrm{R}}-\omega_{\eta}$ will be
present in the spectrum. However, due to the internal reflection in our
experimental data, a moving pulse always comes before and after and we cannot
test this hypothesis experimentally. But we can check our hypothesis with the
numerical simulations. This is done in Section IV, and it confirms our
expectation.

##  S6 Numerical simulations of phonon dynamics

In the simulations, we model the two THz electric field pulse shapes by the
Gaussian derivative function (S3), using the parameters $t_{0}=0.36$ ps,
$\omega_{0}=1/t_{0}$ and the electric field $E_{0}$ set to unity. The
simulated waveforms and their associated spectra are shown in Fig. S9.

Figure S9: (a) Simulated THz pulses by the Gaussian derivative function,
including the reflection to resemble the experimental situation. (b) The
simulated THz spectrum, which is highly similar to the actual spectrum shown
in Fig. S8(b).

The two coupled equations of motion that were simulated are:

$\displaystyle\ddot{Q}_{\mathrm{IR}}+2\zeta_{\mathrm{IR}}\omega_{\mathrm{IR}}\dot{Q}_{\mathrm{IR}}+\omega_{\mathrm{IR}}^{2}Q_{\mathrm{IR}}$
$\displaystyle=$
$\displaystyle\frac{Z^{*}}{M_{\mathrm{IR}}}E+\frac{1}{M_{\mathrm{IR}}}bQ_{\mathrm{R}}E+\frac{2}{M_{\mathrm{IR}}}cQ_{\mathrm{R}}Q_{\mathrm{IR}}$
(S14)
$\displaystyle\ddot{Q}_{\mathrm{R}}+2\zeta_{\mathrm{R}}\omega_{\mathrm{R}}\dot{Q}_{\mathrm{R}}+\omega_{\mathrm{R}}^{2}Q_{\mathrm{R}}$
$\displaystyle=$
$\displaystyle\frac{\delta_{\mathrm{R}}}{M_{\mathrm{R}}}E^{2}+\frac{1}{M_{\mathrm{R}}}bQ_{\mathrm{IR}}E+\frac{1}{M_{\mathrm{R}}}Q_{\mathrm{IR}}^{2}$
(S15)

where we took a Raman-active phonon with a frequency of
$\tfrac{1}{2\pi}\omega_{\mathrm{R}}=3.14$ THz and an IR-active phonon with a
frequency of $\tfrac{1}{2\pi}\omega_{\mathrm{IR}}=1.47$ THz. The damping
ratios were set at $\zeta_{\mathrm{R}}=0.005$ and $\zeta_{\mathrm{IR}}=0.015$
and the effective masses were put to unity $M_{\mathrm{IR}}=M_{\mathrm{R}}=1$
as well as the effective charge $Z^{*}=1$. In the following sections, we
display the numerical solutions for each of the three nonlinear coupling terms
in Eq. (S15), separately. In all cases, the relevant parameter
($b,c,\delta_{\mathrm{R}})$ was set to $1$.

The equations were solved numerically using the ordinary differential equation
integrator scipy.integrate.odeint in Python. Just as in the experiment, we can
take the signal with only one of the pulses present ($Q_{\mathrm{R,1}}$ and
$Q_{\mathrm{R,2}}$) or both pulses present ($Q_{\mathrm{R},12}$) and calculate
the by definition nonlinear part:

$Q_{\mathrm{R,NL}}(t,\tau)=Q_{\mathrm{R},12}(t,\tau)-Q_{\mathrm{R},1}(t,\tau)-Q_{\mathrm{R},2}(t,\tau),$
(S16)

###  I THz-mediated coupling between the modes $\propto
Q_{\mathrm{R}}Q_{\mathrm{IR}}E$

Figure S10: Simulated results of $Q_{\mathrm{R}}$ for the driving force
$\propto Q_{\mathrm{IR}}E$ $(b=1)$ in time domain, of both the total signal
$Q_{\mathrm{R},12}$ and the non-linear signal $Q_{\mathrm{R,NL}}$.  Figure
S11: Normalized 2D FFT associated to Fig. S10.

###  II Photonic impact $\propto E^{2}$

The result for the photonic impact term is shown in Fig. S12. Regarding the
nonlinear part, there is only a signal around the pulse overlap $\tau=0$. This
is in stark contrast with the experimental result, as the experiment showed
that the long-lived IR-active phonons can mediate an interaction between two
subsequent THz pulses over pump-pump separation times $\tau$ much longer than
the pulse duration itself. Additionally, when taking the 2D FFT over the
entire time domain, we obtain the result presented in Fig. S13. It shows a
broad line in the excitation spectrum, which clearly does not produce the
experimental result. Moreover, in the analysis of the experimental data, we
exclude the region around the pulse overlap (see Section S4), meaning that we
removed most of the photonic impact contributions to the signal (if present).
Therefore, we exclude photonic impact as the dominant mechanism describing our
data.

Figure S12: Simulated results when taking only photonic impact $\propto E^{2}$
into account $(\delta_{\mathrm{R}}=1)$. Both the total signal
$Q_{\mathrm{R},12}$ and the non-linear signal $Q_{\mathrm{R,NL}}$ are
depicted. Figure S13: Normalized 2D FFT associated to Fig. S12.

###  III Anharmonic lattice dynamics $\propto
Q_{\mathrm{R}}Q_{\mathrm{IR}}^{2}$

The simulation result for the driving term $\propto Q_{\mathrm{IR}}^{2}$ for
the Raman-active phonon is depicted in Fig. S14. In the corresponding 2D FFT
shown in Fig. S15, one can also see the two peaks, similar to the experimental
result. As the driving term $\propto Q_{E_{u}}^{2}$ is not resonant with the
$E_{g}$ mode frequency, i.e. $2f_{E_{u}}\neq f_{E_{g}}$, we can argue it
effectively reduces to the THz-mediated phonon-phonon coupling mechanism as
discussed in Section I. Namely, during the pulse overlap, the $Q_{E_{u}}$
normal coordinate undergoes a forced oscillation in which essentially
$Q_{E_{u}}(t)\propto E_{\mathrm{THz}}(t)$. Therefore, while the THz pulse is
present, one can effectively replace one $Q_{E_{u}}$ by $E$ in the anharmonic
lattice term $Q_{E_{u}}^{2}$ to obtain the driving term $Q_{E_{u}}E$, which
does obey the resonant condition due to the broad spectrum of the THz pulse.
Hence, the same analytical analysis from Sec. S5 holds (where we only
considered the interaction potential $Q_{\mathrm{R}}Q_{\mathrm{IR}}E$.
However, in addition to the two peaks, we see that the anharmonic lattice term
in Fig. S15 also produces a large spectral component away from the detection
frequency of the $E_{g}$ phonon, something which is not seen in the
experiment. This is because the presence of an IR-phonon constantly forces the
oscillation of the Raman-active mode, and exists also when the THz electric
field is not present. Therefore, although this driving term also reproduces
the two peaks, we exclude it as the dominant mechanism.

Figure S14: Simulation of $Q_{\mathrm{R,12}}(t,\tau)$ and
$Q_{\mathrm{R,NL}}(t,\tau)$, taking only the nonlinear coupling term $\propto
Q_{\mathrm{R}}Q_{\mathrm{IR}}^{2}$ into account ($c=1$). Figure S15:
Normalized 2D FFT associated to Fig. S14.

###  IV Distinguishing phonon and complementing peaks

By deriving the analytical solutions of the equations of motion (Section S5),
we hypothesized the possibility to distinguish the actual phonon peak from the
complementary ones when considering the THz-mediated mode coupling term
($\propto Q_{\mathrm{R}}Q_{\mathrm{IR}}E$), by taking the FFT over the
negative and positive $\tau$, respectively. However, this effect is hampered
experimentally by internal reflections. Therefore, we repeated the numerical
simulation as in Sec. I. but without the simulated internal reflections. The
results by taking the FFT over negative and positive $\tau$ confirm our
hypothesis, as can be seen in Fig. S16.

Figure S16: Simulated results as in Fig. S10 (interaction potential $\propto
Q_{\mathrm{R}}Q_{\mathrm{IR}}E$), but without the internal THz reflection.
Taking the 2D FFT in the region where the moving pulse arrives later $\tau>1$
shows the complementary phonon peak in the frequency domain. However, when
making the 2D FFT transformation in the region where the moving pulse arrives
first ($\tau<-1$), we obtain the actual phonon frequency in the spectrum.

## References

  * Leitenstorfer  _et al._ [1999] A. Leitenstorfer, S. Hunsche, J. Shah, M. C. Nuss, and W. H. Knox, Detectors and sources for ultrabroadband electro-optic sampling: Experiment and theory, Applied Physics Letters 74, 1516 (1999).
  * Casalbuoni  _et al._ [2005] S. Casalbuoni, H. Schlarb, B. Schmidt, B. Steffen, P. Schmuser, and A. Winter, Numerical studies on the electro-optic sampling of relativistic electron bunches, in  _Proceedings of the 2005 Particle Accelerator Conference_ (2005) pp. 3070–3072.
  * Mal  _et al._ [2019] P. Mal, G. Bera, G. R. Turpu, S. K. Srivastava, A. Gangan, B. Chakraborty, B. Das, and P. Das, Vibrational spectra of $\mathrm{Pb}_{2}\mathrm{Bi}_{2}\mathrm{Te}_{3}$, $\mathrm{PbBi}_{2}\mathrm{Te}_{4}$, and $\mathrm{PbBi}_{4}\mathrm{Te}_{7}$ topological insulators: temperature-dependent Raman and theoretical insights from DFT simulations, Phys. Chem. Chem. Phys. 21, 15030 (2019).
  * Altmann and Herzig [1994] S. L. Altmann and P. Herzig,  _Point-group theory tables_ (Oxford University Press, 1994).
  * Rodriguez-Vega  _et al._ [2020] M. Rodriguez-Vega, A. Leonardo, and G. A. Fiete, Group theory study of the vibrational modes and magnetic order in the topological antiferromagnet $\mathrm{Mn}\mathrm{Bi}_{2}\mathrm{Te}_{4}$, Phys. Rev. B 102, 104102 (2020).
  * Kobiałka  _et al._ [2022] A. Kobiałka, M. Sternik, and A. Ptok, Dynamical properties of the magnetic topological insulator $T{\mathrm{Bi}}_{2}{\mathrm{Te}}_{4}(T=\mathrm{Mn},\mathrm{Fe})$: Phonons dispersion, Raman active modes, and chiral phonons study, Phys. Rev. B 105, 214304 (2022).
  * Johnson  _et al._ [2019] C. L. Johnson, B. E. Knighton, and J. A. Johnson, Distinguishing nonlinear terahertz excitation pathways with two-dimensional spectroscopy, Phys. Rev. Lett. 122, 073901 (2019).
  * Loudon [1964] R. Loudon, The Raman effect in crystals, Advances in Physics 13, 423 (1964).
  * Juraschek  _et al._ [2017] D. M. Juraschek, M. Fechner, and N. A. Spaldin, Ultrafast structure switching through nonlinear phononics, Phys. Rev. Lett. 118, 054101 (2017).
  * dos Santos Silva  _et al._ [2022] A. dos Santos Silva, É. V. Guimarães, N. F. Cano, and R. S. da Silva, Study of vibrational properties of $\mathrm{Bi}_{2-x}\mathrm{Mn}_{x}\mathrm{Te}_{3}$ nanocrystals in host glass: Effect of $x\mathrm{Mn}$, Journal of Raman Spectroscopy 53, 95 (2022).