# Supplementary materials for ”Time dynamics of multi-photon scattering in a
two-level mixer”

Suggested keywords

###### Contents

  1. I Sample fabrication
  2. II Measurement setup
  3. III Qubit Spectroscopy
  4. IV Rabi oscillations
  5. V Mixing results for 5 and 6 driving pulses
  6. VI Analytical expressions for mixing maps
  7. VII Numerical simulation of mixing

##  I Sample fabrication

The chip with individual qubits coupled to a coplanar waveguide is fabricated
with typical three-lithography process. The high-resistance Si wafer was
cleaned with the IPA and RCA-1 solution. Surface silicon oxide layer contains
many unwanted defects and thereby it is removed by HF after cleaning. Next,
the aluminum film of 100 nm is evaporated onto clean wafer with e-beam
evaporation machine Plassys MEB-550 in a vacuum better than $10^{-7}$ mBar.
The mask creating the pattern for rough structures was created with two-layer
stack of resists: 800 nm-thick PMGI resist was covered with 200 nm of positive
resist AZ-1505. The stack was baked on hotplate at 160 ∘C for 5 minutes and
then the mask pattern is exposed by laser maskless lithographer Heidelberg
MLA100. The exposed resist is developed with chlorine-based ICP etching, and
thereby waveguide and contact pads are formed. Following, the stack for e-beam
lithography is formed with Copolymer and ARP-6200.04 resists, which is exposed
by electron beam lithographer Crestec-CABL9000 in order to form josephson
junctions. After the mask is developed, we perform shadow-angle two-layer
evaporation of junctions with the thicknesses of 25 and 40 nm respectively,
and then connect junctions to qubit’s capacitor with bandages. This requires
another lithography and evaporation with removing aluminum oxide formed on top
of previous layer by in situ Ar ion milling. The test junctions are fabricated
along with qubits, and their resistance is tested with the probe station.

##  II Measurement setup

Figure S1: The experimental scheme for measurements of fields scattered by a
single qubit in the waveguide.

The scheme of measurements is presented in Fig. S1. The room-temperature part
includes up-conversion and down-conversion schemes built from commercially
available AWGs, mixers, ADCs and rf-generators. For spectroscopic measurements
and calibration we also use VNA and SA. The low-temperature part includes
input cryogenic coaxial line properly thermalized with attenuators and
filters, and output line with amplification chain. Typical solution in our
area is to use pair of cryogenic isolators to prevent thermal noise
propagation from 4K stage down to sample via output cabling. First of all, our
intent is to generate two pulse sequences with close carrier frequencies and
different delays, and to do this we use the following scheme of up-conversion.
We choose intermediate frequency $\omega_{\text{IF}}$ to be equal to 50-100
MHz, and digitally prepare input waveform, which is the sum of two train
sequences of non-overlapping pulses with carriers
$\omega_{\text{IF}}+\delta\omega$ and $\omega_{\text{IF}}-\delta\omega$. This
waveform is characterized by the period of beats due to small detunings
$T_{b}=2\pi/2\delta\omega$. Therefore, the duration of waveforms is
conveniently chosen to be the multiple of $T_{b}$. The individual pulses are
much shorter than $T_{b}$. The time between rising edge of two consecutive
pulses with a certain carrier frequency is the fraction of $T_{b}$, but
nevertheless it is much greater than the relaxation lifetime of a qubit. The
waveforms are loaded into AWG, the output channels of which are connected to I
and Q ports of quadrature mixer. The continuous wave at frequency
$\omega_{\text{LO}}$ is given by the rf-generator, output of which is
connected to L port of mixer. We calibrate up-conversion by adjusting the
phases and offsets of waveforms with carrier frequency of $\omega_{\text{IF}}$
in order to maximize the signal at $\omega_{\text{LO}}+\omega_{\text{IF}}$ and
to minimize all other harmonics.

As for down-conversion part, we exploit the same quadrature mixer. The
scattered field undergoes amplification and then comes into $R$-port of the
output mixer. The LO signal goes from the same generator as used for up-
conversion. Down-converted quadratures from I and Q ports are digitized with
the two synchronized channels of ADC, and the acquisition is triggered by AWG.
After moderate averaging (typically millions of waveform periods), we take
complex Fourier transform of signal and get the results presented in Fig. 2
and 3 of main text.

##  III Qubit Spectroscopy

The single-tone spectroscopy of the qubit is shown in Fig. S2. We measure the
transmission of the signal through the waveguide with the use of VNA. The
frequency of stimulus is swept in a wide range, and when it is resonant with
the qubit, the dip in transmission is observed, see Fig. S2(a). We normalize
the transmission data and fit the qubit line with a circle-fit method. The
optimal fit parameters include internal and external (coupling) quality
factors $Q_{i}$ and $Q_{c}$, respectively, and the frequency of the qubit
$f_{r}$. Sweeping current values and repeating transmission measurement gives
us the spectroscopy of the qubit, Fig. S2(b).

Figure S2: Single-tone transmission spectroscopy of the qubit in the
waveguide. (a) Resonant line of the qubit tuned close to upper sweetspot.
Points are data measured with VNA, and line is the best fit to the data. The
important fitting parameters are written in the figure title. (b) The current-
dependent single tone spectroscopy of the qubit.

##  IV Rabi oscillations

With the presented experimental setup we are able to measure the dynamics of
the qubit under strong resonant microwave pulses. For this, we direct into a
waveguide a sequence of short pulses divided by long intervals. The sequence
is formed by our up-conversion scheme. Being resonant with the qubit, the
pulse causes the rotational dynamics of the latter, also known as Rabi
oscillations. Let us consider the rotation around $x$-axis on the Bloch
sphere, then the qubit state evolves in time as:

$\displaystyle\braket{\sigma_{x}}$ $\displaystyle=0,$ (S1)
$\displaystyle\braket{\sigma_{y}}$
$\displaystyle=\frac{2\Omega}{2\Omega^{2}+\Gamma_{1}^{2}}\left[\Gamma_{1}-e^{-3\Gamma_{1}t/4}\left(\Gamma_{1}\cos\Omega^{\prime}t+\frac{(4\Omega^{2}-\Gamma_{1}^{2})}{4\Omega^{\prime}}\sin\Omega^{\prime}t\right)\right],$
(S2) $\displaystyle\braket{\sigma_{z}}$
$\displaystyle=-\frac{\Gamma_{1}}{2\Omega^{2}+\Gamma_{1}^{2}}\left[\Gamma_{1}+e^{-{3}\Gamma_{1}t/4}\left(\frac{3\Omega^{2}}{2\Omega^{\prime}}\sin\Omega^{\prime}t+\frac{2\Omega^{2}}{\Gamma_{1}}\cos\Omega^{\prime}t\right)\right].$
(S3)

Here we suppose that the pure dephasing and non-radiative relaxation are
negligible, and that the drive is exactly in resonance with the transition.
Also we introduced the generalized Rabi frequency
$\Omega^{\prime}=\sqrt{\Omega^{2}-\Gamma_{1}^{2}/16}$. For measuring the
amplitude of the emitted field (which is proportional to
$\braket{\sigma_{y}}$), we down-convert the field emitted right after the
driving pulse is gone. After that, we digitize both quadratures the signal at
$\omega_{\text{IF}}$ and average over many single traces. Then we apply fast
Fourier transform (FFT) in order to get the complex amplitude of field. The
single quadrature of field measured for different lengths of driving pulse is
presented at Fig. S3(a).

Figure S3: Rabi oscillations. (a) A single trace representing qubit dynamics
under resonance microwave pulses. Blue crosses show raw data, and the line is
optimal fit according to simplified version of Eq. (S2). (b) The Rabi
frequencies versus voltage amplitude of driving pulses with a linear fit.

The fitting of experimental points is done with the help of Eq. (S2) and
allows to extract the Rabi frequency and the decay time. We repeat this for
various driving amplitudes and finally obtain the dependence of $\Omega$ on
the voltage amplitude, see Fig. S3(b).

##  V Mixing results for 5 and 6 driving pulses

Here we present the extension of results presented in Fig. 3 of main text. The
angle dependencies of field amplitudes of side peaks measured for case of
five- and six-pulse configuration are presented in Fig. S4. Comparing to the
driving sequences with less applied pulses, we can note two features of maps.
First, the number of elastic side peaks is increasing: each additional pulse
causes the emergence of two extra components.

Figure S4: The dependence of complex amplitudes of elastically scattered
fields on rotation angle of pulses for each side spectral component.

Second, for components with frequencies close to initial drives maps are
becoming more and more complicated. The reason of this is increasing number of
contributions of different orders, whereas for most-detuned components there
is only one contribution into furthest components appeared only after the very
last pulse.

##  VI Analytical expressions for mixing maps

In this section, we describe in detail how analytical expressions for the
mixing components are calculated. To derive the dynamics of $p$-th side
spectral components when $N$ pulses are applied, we calculate the matrix
element $M^{(N)}_{p}$, defined in Eq. (2) of the main text. It is the
expectation value of operator $\sigma^{-}_{p}$ on the quantum state of the
kind:

$\ket{\phi(t)}=\Lambda_{N}\ket{\phi(0)},$ (S4)

where $\Lambda_{N}$ is total evolution operator of the system defined in Eq.
(3) of the main text. Operators $\Lambda_{N}$ can be directly written by
multiplying all the evolution operators defined in Eq. (4), but the resulting
expressions are very cumbersome. However, the calculations could be simplified
significantly, because we are not interested in the full expression for
$\Lambda_{N}$, but instead looking for $M^{(N)}_{p}$. In the main text, we
introduce operators $a$ and $b$ such that $a^{\dagger}\propto\sigma^{+}_{-}$,
$b^{\dagger}\propto\sigma^{+}_{+}$. Thus, typically for this kind of
calculations,
$a\ket{\phi(0)}=b\ket{\phi(0)}=\bra{\phi(0)}a^{\dagger}=\bra{\phi(0)}b^{\dagger}=0$.
As well, any term which contains any of the following multipliers:
$(\cdot)(\cdot),(\cdot)^{\dagger}(\cdot)^{\dagger}$, where $(\cdot)=\\{a,b\\}$
— will reduce to zero and give no contribution to $M^{(N)}_{p}$. This allows
to introduce the operator $\Lambda^{0}_{N}$ which contains less terms than
$\Lambda_{N}$. Then the simplified expression will be:

$M^{(N)}_{p}=\bra{\phi(0)}(\Lambda^{0}_{N})^{\dagger}\sigma^{-}_{p}\Lambda^{0}_{N}\ket{\phi(0)},$
(S5)

Taking into account Eqs. (2-5) of main text, we do algebraic computations and
express $\Lambda^{0}_{N}$ as:

$\displaystyle\begin{split}\Lambda^{0}_{2}&=1-a^{\dagger}-ba^{\dagger}-b^{\dagger},\quad\Lambda_{3}=1-aa^{\dagger}-ab^{\dagger}-2a^{\dagger}+a^{\dagger}ba^{\dagger}-ba^{\dagger}-b^{\dagger},\end{split}$
$\displaystyle\begin{split}\Lambda^{0}_{4}&=1-aa^{\dagger}-ab^{\dagger}-2a^{\dagger}+a^{\dagger}ba^{\dagger}-3ba^{\dagger}+(ba^{\dagger})^{2}-bb^{\dagger}-2b^{\dagger}+b^{\dagger}aa^{\dagger}+b^{\dagger}ab^{\dagger}+b^{\dagger}ba^{\dagger},\end{split}$
$\displaystyle\begin{split}\Lambda^{0}_{5}&=1-3aa^{\dagger}+aa^{\dagger}ba^{\dagger}-3ab^{\dagger}+ab^{\dagger}aa^{\dagger}+(ab^{\dagger})^{2}+ab^{\dagger}ba^{\dagger}-3a^{\dagger}+a^{\dagger}aa^{\dagger}+a^{\dagger}ab^{\dagger}+4a^{\dagger}ba^{\dagger}-\\\
&-(a^{\dagger}b)^{2}a^{\dagger}+a^{\dagger}bb^{\dagger}-3ba^{\dagger}+(ba^{\dagger})^{2}-bb^{\dagger}-2b^{\dagger}+b^{\dagger}aa^{\dagger}+b^{\dagger}ab^{\dagger}+b^{\dagger}ba^{\dagger},\end{split}$
(S6)
$\displaystyle\begin{split}\Lambda^{0}_{6}&=1-3aa^{\dagger}+aa^{\dagger}ba^{\dagger}-3ab^{\dagger}+ab^{\dagger}aa^{\dagger}+(ab^{\dagger})^{2}+ab^{\dagger}ba^{\dagger}-3a^{\dagger}+a^{\dagger}aa^{\dagger}+a^{\dagger}ab^{\dagger}+4a^{\dagger}ba^{\dagger}-(a^{\dagger}b)^{2}a^{\dagger}+\\\
&+a^{\dagger}bb^{\dagger}-6ba^{\dagger}+ba^{\dagger}aa^{\dagger}+ba^{\dagger}ab^{\dagger}+5(ba^{\dagger})^{2}-(ba^{\dagger})^{3}+ba^{\dagger}bb^{\dagger}-3bb^{\dagger}+bb^{\dagger}aa^{\dagger}+bb^{\dagger}ab^{\dagger}+bb^{\dagger}ba^{\dagger}-\\\
&-3b^{\dagger}+4b^{\dagger}aa^{\dagger}-b^{\dagger}aa^{\dagger}ba^{\dagger}+4b^{\dagger}ab^{\dagger}-(b^{\dagger}a)^{2}a^{\dagger}-b^{\dagger}(ab^{\dagger})^{2}-b^{\dagger}ab^{\dagger}ba^{\dagger}+4b^{\dagger}ba^{\dagger}-b^{\dagger}(ba^{\dagger})^{2}+b^{\dagger}bb^{\dagger}.\end{split}$

However, we note that $\Lambda^{0}_{N}$ in this form still has terms which do
not contribute into $M^{(N)}_{p}$. These are the terms where leftmost operator
is $a$ or $b$, and they disappear because products $\sigma^{-}_{p}a$ and
$\sigma^{-}_{p}b$ result in zero during the calculation of (S5). However, for
$(\Lambda^{0}_{N})^{\dagger}$ the different terms could be removed — those
where the rightmost operator is $a$ or $b$, similarly, because
$a\sigma^{-}_{p}$ and $b\sigma^{-}_{p}$ result in zero. Therefore, (S5) can be
now rewritten as:

$M^{(N)}_{p}=\bra{\phi(0)}\Lambda^{\leftarrow}_{N}\sigma^{-}_{p}\Lambda^{\rightarrow}_{N}\ket{\phi(0)}.$
(S7)

Here, operators $\Lambda^{\rightarrow}_{N}$ and $\Lambda^{\leftarrow}_{N}$ are
obtained from $\Lambda^{0}_{N}$ after discarding the proper terms. We come to
the following expressions:

$\displaystyle\begin{split}\Lambda^{\rightarrow}_{2}&={}-a^{\dagger}-b^{\dagger},\quad\Lambda^{\rightarrow}_{3}=-2a^{\dagger}+a^{\dagger}ba^{\dagger}-b^{\dagger},\end{split}$
$\displaystyle\begin{split}\Lambda^{\rightarrow}_{4}&={}-2a^{\dagger}+a^{\dagger}ba^{\dagger}-2b^{\dagger}+b^{\dagger}aa^{\dagger}+b^{\dagger}ab^{\dagger}+b^{\dagger}ba^{\dagger},\end{split}$
$\displaystyle\begin{split}\Lambda^{\rightarrow}_{5}&={}-3a^{\dagger}+a^{\dagger}aa^{\dagger}+a^{\dagger}ab^{\dagger}+4a^{\dagger}ba^{\dagger}-(a^{\dagger}b)^{2}a^{\dagger}+a^{\dagger}bb^{\dagger}-2b^{\dagger}+b^{\dagger}aa^{\dagger}+b^{\dagger}ab^{\dagger}+b^{\dagger}ba^{\dagger},\end{split}$
(S8)
$\displaystyle\begin{split}\Lambda^{\rightarrow}_{6}&={}-3a^{\dagger}+a^{\dagger}aa^{\dagger}+a^{\dagger}ab^{\dagger}+4a^{\dagger}ba^{\dagger}-(a^{\dagger}b)^{2}a^{\dagger}+a^{\dagger}bb^{\dagger}-3b^{\dagger}+4b^{\dagger}aa^{\dagger}-b^{\dagger}aa^{\dagger}ba^{\dagger}+4b^{\dagger}ab^{\dagger}-\\\
&-{}(b^{\dagger}a)^{2}a^{\dagger}-b^{\dagger}(ab^{\dagger})^{2}-b^{\dagger}ab^{\dagger}ba^{\dagger}+4b^{\dagger}ba^{\dagger}-b^{\dagger}(ba^{\dagger})^{2}+b^{\dagger}bb^{\dagger}.\end{split}$

Similarly, for left-acting operators:

$\displaystyle\begin{split}\Lambda^{\leftarrow}_{2}&=1-ab^{\dagger},\quad\Lambda^{\leftarrow}_{3}=1-aa^{\dagger}-ab^{\dagger}-ba^{\dagger},\end{split}$
$\displaystyle\begin{split}\Lambda^{\leftarrow}_{4}&=1-aa^{\dagger}-3ab^{\dagger}+(ab^{\dagger})^{2}-ba^{\dagger}-bb^{\dagger},\end{split}$
$\displaystyle\begin{split}\Lambda^{\leftarrow}_{5}&=1-3aa^{\dagger}+aa^{\dagger}ba^{\dagger}-3ab^{\dagger}+ab^{\dagger}aa^{\dagger}+(ab^{\dagger})^{2}+ab^{\dagger}ba^{\dagger}-3ba^{\dagger}+a^{\dagger}-bb^{\dagger},\end{split}$
(S9)
$\displaystyle\begin{split}\Lambda^{\leftarrow}_{6}&=1-3aa^{\dagger}+aa^{\dagger}ab^{\dagger}+aa^{\dagger}ba^{\dagger}+aa^{\dagger}bb^{\dagger}-6ab^{\dagger}+ab^{\dagger}aa^{\dagger}+5(ab^{\dagger})^{2}-(ab^{\dagger})^{3}+\\\
&+ab^{\dagger}ba^{\dagger}+ab^{\dagger}bb^{\dagger}-3ba^{\dagger}+ba^{\dagger}ab^{\dagger}+(ba^{\dagger})^{2}+ba^{\dagger}bb^{\dagger}-3bb^{\dagger}+bb^{\dagger}ab^{\dagger}.\end{split}$

With the use of (S8) and (S9), after trivial algebra we come to the expression
for $M^{(2)}_{p}$, which is presented as Eq. (6) of the main text. As well, we
get expressions for all studied cases, in particular, for 3 and 4 applied
pulses:

$M^{(3)}_{p}=\cos^{-4}\frac{\theta_{-}}{2}\cos^{-2}\frac{\theta_{+}}{2}\bra{\phi(0)}-\underbrace{ba^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}}_{\exp{(-5i\delta\omega
t)}}+\underbrace{\sigma^{-}_{p}a^{\dagger}ba^{\dagger}-aa^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}+2ba^{\dagger}\sigma^{-}_{p}a^{\dagger}}_{\exp({-3i\delta\omega
t})}\\\
-\underbrace{2\sigma^{-}_{p}a^{\dagger}+2aa^{\dagger}\sigma^{-}_{p}a^{\dagger}-ab^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}+ba^{\dagger}\sigma^{-}_{p}b^{\dagger}}_{\exp{(-i\delta\omega
t)}}\\\
-\underbrace{\sigma^{-}_{p}b^{\dagger}+aa^{\dagger}\sigma^{-}_{p}b^{\dagger}+2ab^{\dagger}\sigma^{-}_{p}a^{\dagger}}_{\exp{(i\delta\omega
t)}}+\underbrace{ab^{\dagger}\sigma^{-}_{p}b^{\dagger}}_{\exp{(3i\delta\omega
t)}}\ket{\phi(0)},$ (S10)

$M^{(4)}_{p}=\cos^{-4}\frac{\theta_{-}}{2}\cos^{-4}\frac{\theta_{+}}{2}\bra{\phi(0)}\underbrace{(ab^{\dagger})^{2}\sigma^{-}_{p}b^{\dagger}ab^{\dagger}}_{\exp{(7\delta\omega
t)}}-\underbrace{3ab^{\dagger}\sigma^{-}_{p}b^{\dagger}ab^{\dagger}-2(ab^{\dagger})^{2}\sigma^{-}_{p}b^{\dagger}+(ab^{\dagger})^{2}\sigma^{-}_{p}b^{\dagger}aa^{\dagger}}_{\exp{(5i\delta\omega
t)}}\\\
+\underbrace{\sigma^{-}_{p}b^{\dagger}ab^{\dagger}-aa^{\dagger}\sigma^{-}_{p}b^{\dagger}ab^{\dagger}+6ab^{\dagger}\sigma^{-}_{p}b^{\dagger}-3ab^{\dagger}\sigma^{-}_{p}b^{\dagger}aa^{\dagger}-2(ab^{\dagger})^{2}\sigma^{-}_{p}a^{\dagger}+(ab^{\dagger})^{2}\sigma^{-}_{p}b^{\dagger}ba^{\dagger}-bb^{\dagger}\sigma^{-}_{p}b^{\dagger}ab^{\dagger}}_{\exp{(3i\delta\omega
t)}}\\\
-\underbrace{2\sigma^{-}_{p}b^{\dagger}\\!+\\!\sigma^{-}_{p}b^{\dagger}aa^{\dagger}\\!+\\!2aa^{\dagger}\sigma^{-}_{p}b^{\dagger}\\!-\\!aa^{\dagger}\sigma^{-}_{p}b^{\dagger}aa^{\dagger}\\!+\\!6ab^{\dagger}\sigma^{-}_{p}a^{\dagger}\\!-\\!3ab^{\dagger}\sigma^{-}_{p}b^{\dagger}ba^{\dagger}\\!+\\!(ab^{\dagger})^{2}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}\\!-\\!ba^{\dagger}\sigma^{-}_{p}b^{\dagger}ab^{\dagger}\\!+\\!2bb^{\dagger}\sigma^{-}_{p}b^{\dagger}\\!-\\!bb^{\dagger}\sigma^{-}_{p}b^{\dagger}aa^{\dagger}}_{\exp{(i\delta\omega
t)}}\\\
-\underbrace{2\sigma^{-}_{p}a^{\dagger}+\sigma^{-}_{p}b^{\dagger}ba^{\dagger}+2aa^{\dagger}\sigma^{-}_{p}a^{\dagger}-aa^{\dagger}\sigma^{-}_{p}b^{\dagger}ba^{\dagger}-3ab^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}+2ba^{\dagger}\sigma^{-}_{p}b^{\dagger}-ba^{\dagger}\sigma^{-}_{p}b^{\dagger}aa^{\dagger}+2bb^{\dagger}\sigma^{-}_{p}a^{\dagger}-bb^{\dagger}\sigma^{-}_{p}b^{\dagger}ba^{\dagger}}_{\exp{(-i\delta\omega
t)}}\\\
+\underbrace{\sigma^{-}_{p}a^{\dagger}ba^{\dagger}-aa^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}+2ba^{\dagger}\sigma^{-}_{p}a^{\dagger}-ba^{\dagger}\sigma^{-}_{p}b^{\dagger}ba^{\dagger}-bb^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}}_{\exp{(-3i\delta\omega
t)}}-\underbrace{ba^{\dagger}\sigma^{-}_{p}a^{\dagger}ba^{\dagger}}_{\exp{(-5i\delta\omega
t)}}\ket{\phi(0)}.$ (S11)

Expressions for $N=5$ and $N=6$ are omitted here, but could be written by
similar procedure. Finally, substituting operators from Eqs. (5) into Eqs.
(S10), (S11), and then using the result in Eq. (7) of the main text, we come
to the final expressions for measured quadrature components of side
frequencies. For $N=3$ we have:

$\displaystyle V^{(3)}_{-5}$
$\displaystyle=-\sin^{3}{\frac{\theta_{-}}{2}}\sin^{2}{\frac{\theta_{+}}{2}}\cos{\frac{\theta_{-}}{2}},\qquad
V^{(3)}_{-3}=\left(2\cos{\theta_{-}}+1\right)\sin^{2}{\frac{\theta_{-}}{2}}\sin{\frac{\theta_{+}}{2}}\cos{\frac{\theta_{+}}{2}},$
$\displaystyle V^{(3)}_{-1}$
$\displaystyle=-\frac{1}{8}\left(3\cos{\theta_{+}}+1\right)\sin{2\theta_{-}},\qquad
V^{(3)}_{1}=\left(1-2\cos{\theta_{-}}\right)\sin{\frac{\theta_{+}}{2}}\cos^{2}{\frac{\theta_{-}}{2}}\cos{\frac{\theta_{+}}{2}},$
(S12) $\displaystyle V^{(3)}_{3}$
$\displaystyle=\sin{\frac{\theta_{-}}{2}}\sin^{2}{\frac{\theta_{+}}{2}}\cos^{3}{\frac{\theta_{-}}{2}}.\qquad\quad$

For $N=4$, the expressions are:

$\displaystyle V^{(4)}_{7}$
$\displaystyle=\frac{\sin^{2}{\frac{\theta_{-}}{2}}\sin^{4}{\frac{\theta_{+}}{2}}\sin{{\theta_{-}}}}{\cos{\theta_{+}}+1},\qquad
V^{(4)}_{5}=-\frac{\left(3\cos{\theta_{-}}+2\right)\sin^{2}{\frac{\theta_{-}}{2}}\sin^{3}{\frac{\theta_{+}}{2}}}{\cos{\frac{\theta_{+}}{2}}},$
$\displaystyle V^{(4)}_{3}$
$\displaystyle=\frac{\left(15\cos{\theta_{-}}\cos{\theta_{+}}+11\cos{\theta_{-}}+\cos{\theta_{+}}+1\right)\sin{\theta_{-}}\sin^{2}{\frac{\theta_{+}}{2}}}{4\left(\cos{\theta_{+}}+1\right)}$
$\displaystyle V^{(4)}_{1}$
$\displaystyle=\left[\sin^{2}{\theta_{-}}\left(5\cos{\theta_{+}}+1\right)-\cos{\theta_{+}}\left(\cos{\theta_{-}}+3\right)\right]\tan{\frac{\theta_{+}}{2}}$
(S13) $\displaystyle V^{(4)}_{-1}$
$\displaystyle=\frac{\sin{\theta_{-}}}{8\left(\cos{\theta_{+}}+1\right)}\left[\left(1-\cos{\theta_{+}}\right)^{2}\left(-15\cos{\theta_{-}}-1\right)+\cos\theta_{-}\left(20-36\cos\theta_{+}\right)\right]$
$\displaystyle V^{(4)}_{-3}$
$\displaystyle=\left(\frac{\cos{\theta_{-}}}{2}\left(3\cos{\theta_{+}}+1\right)+\cos{\theta_{+}}\right)\sin^{2}{\frac{\theta_{-}}{2}}\tan{\frac{\theta_{+}}{2}}$
$\displaystyle V^{(4)}_{-5}$
$\displaystyle=-\sin^{3}{\frac{\theta_{-}}{2}}\sin^{2}{\frac{\theta_{+}}{2}}\cos{\frac{\theta_{-}}{2}}.$

The dependencies (S12) and (S13) are plotted in Fig. 3 of main text, and the
good agreement with experimentally measured angle dependencies is observed.

##  VII Numerical simulation of mixing

In this section we disclose the technical details of numerical calculations
exploited to obtain side spectral components, which emerge as a result of wave
mixing between ultrashort bichromatic pulses.

Figure S5: The numerically generated sequence of driving pulses. (a) The one
of pulse sequences used to generate lines at Fig. 2 of the main text. The
qubit is driven by ultra-short pulses with the length of 8-12 ns, and the
period of pulses is 40 us. The slowly varying envelope is from exponential
multipliers $e^{\pm\delta\omega t}$, with $\delta\omega/2\pi=50$ kHz. (b) the
magnified view of three pulses at $\omega_{-},\omega_{+},\omega_{-}$ with the
durations of 8,12 and 8 ns. The padding is 8 ns. The sequence is updated at
each step, and the relative delay of the middle pulse is swept from -50 to 50
ns with respect to the position plotted in (b).

The numerical simulation is done for different configurations of bichromatic
pulsed drive. In previous section we considered the process of analytical
derivation of side spectral components especially for non-overlapping pulses.
In that relatively simple case, the qubit interacts only with one mode at a
time, and never with two modes simultaneously. This allows us to calculate
total evolution by multiplying time-dependent evolution operators of single
mode Jaynes-Cummings model, which come one by one from distinct pulses. The
expression of this evolution operator is commonly known and may be found, for
example in [1], Ch. 15. However, for the case of overlapping red- and blue-
detuned pulses, the similar calculation becomes less straightforward [2],
because there is no analytical solution for the atom driven by two detuned
continuous waves. Particularly, there is a Bessel function dependencies of
side components given in [2], but they are only correct if evolution starts
from ground state, which is not the case for partially overlapped pulses. For
example, if we have two overlapped pulses, then the atom is, first of all, is
rotating by a first single mode driving pulse, and already afterwards it
evolves under two drives, from the start moment of the second pulse.
Therefore, numerical calculations are convenient to demonstrate mixing
spectrum for more complicated sequence of driving pulses.

We numerically simulate the system evolution under time-dependent Hamiltonian
under rotating wave approximation:

$H_{\text{RWA}}(t)=\left[\begin{matrix}0&\Omega_{+}(t)e^{-i\delta\omega
t}+\Omega_{-}(t)e^{i\delta\omega t}\\\ \Omega_{+}(t)e^{i\delta\omega
t}+\Omega_{-}(t)e^{-i\delta\omega t}&0\end{matrix}\right],$ (S15)

where we kept only slowly varying terms from both drives. We employ open-
source `QuTiP` package to construct the short pulses with envelopes
$\Omega_{+}(t)$ and $\Omega_{-}(t)$ which are plotted at Fig. S5 for three
driving pulses. With this, we construct Lindblad master equation with
radiational relaxation and then use `mesolve` function to calculate
$\braket{\sigma_{-}(t)}$.

Figure S6: Fourier transform of $\braket{\sigma_{-}(t)}$

The Fourier transform of expectation value is presented in Fig. S6. We can
distinctly observe three sidebands, two of which are red-shifted and one is
blue- shifted. This corresponds to the qualitative picture of multi-photon
interaction which we discuss in the main text. However, as we described
earlier, the spectrum is more complicated for overlapping pulses. We
illustrate the variety of spectral properties by performing above-described
simulation for various values of the middle pulse delay and in wide range of
Rabi drives $\Omega_{+}=\Omega_{-}=\Omega$. The results are presented at Fig.
S8. The all-dark regions for components $\omega_{-7},\omega_{7}$ and
$\omega_{5}$ are observed for all range of $\Omega$ when the middle pulse (at
$\omega_{+}$) is either very first, very last, or in between other two pulses
(at $\omega_{-}$). This regions correspond to $t_{2}^{\text{center}}<-28$ ns,
$-8<t_{2}^{\text{center}}<8$ ns and $t_{2}^{\text{center}}>28$ ns.

Figure S7: Scattering of three-pulse sequence shown at Fig. S5. Intensities of
elastic components are plotted as a function of effective Rabi amplitude
$\Omega$ and the middle point $t_{2}^{\text{center}}$ of the second pulse at
$\omega_{+}$ shifted with respect to another pair of pulses at $\omega_{-}$.

Figure S8: Scattering of two-pulse sequence. Intensities of elastic components
plotted as a function of effective Rabi amplitude of pulses $\Omega$ and the
middle point $t_{2}^{\text{center}}$ of the second pulse at $\omega_{+}$
shifted with respect to another single pulse at $\omega_{-}$.

Similarly, we make the simulations for a sequence of two pulses, the result of
which also exhibits dark regions where the corresponding mixing processes are
prohibited, see Fig. S8. As expected, the maps are symmetric with respect to
delay and frequency inversion:
$\sigma^{-}_{p}(\Omega,t_{2}^{\text{center}})=\sigma^{-}_{-p}(\Omega,-t_{2}^{\text{center}})$.
These results confirm that the missing peaks at Fig. 1 and 2 of main text are
not only disappear at specific Rabi angles of driving pulses, but are absent
for any amplitudes.

## References

  * Schleich [2011] W. P. Schleich,  _Quantum optics in phase space_ (John Wiley & Sons, 2011).
  * Dmitriev  _et al._ [2017] A. Y. Dmitriev, R. Shaikhaidarov, V. N. Antonov, T. Hönigl-Decrinis, and O. V. Astafiev, Quantum wave mixing and visualisation of coherent and superposed photonic states in a waveguide, Nature Communications 8, 1 (2017).