# The first superconducting final focus quadrupole prototype of the FCC-ee
study

A. Thabuis  M. Koratzinos  G. Kirby  M. Liebsch  C. Petrone
European Organization for Nuclear Research (CERN)  Geneva  Switzerland

###### Abstract

A single aperture Canted-Cosine-Theta (CCT) quadrupole magnet, made of NbTi
superconductors, has been developed for the final focus region of the FCC-ee
study. The conductor layout is optimised to mitigate edge effects on one of
the two sides of the magnet that typically lead to undesired higher-order
multipoles. Experimental results of a prototype, including paraffin wax
impregnation and cryogenic temperature measurements, are presented. The magnet
exhibits no training behaviour, surpassing the nominal current during the
initial ramp. Field quality is excellent, with higher-order multipoles below
$1\times 10^{-4}$ units, consistent with simulations and room temperature
tests. These findings confirm the potential of superconducting CCT magnets to
offer compact solutions for applications demanding stringent field quality.

##  1 INTRODUCTION

The Future Circular Collider (FCC) study plans on building a Higgs and
electroweak factory through collisions of positrons and electrons (FCC-ee)
[1]. Once these particles have been accelerated to their desired energy,
strong so-called final focus quadrupole magnets, focus the two individual
beams around the so-called interaction region where collisions occur.

A Canted-Cosine-Theta (CCT) magnet topology has been proposed for the
realization of the final-focus quadrupole, following a similar work done at
CERN [2]. The conductor layout has been designed such that the edge effects
are reduced on one side and kept unchanged on the other side for comparison. A
prototype of shorter axial length was designed and built in 2019. Its
specifications are summarized in Table 1. The integrated quadrupole component
$B_{2}$ is $2.2$ times less than the full-size final magnet of length $0.7$ m
(initially $1.2$ m-long but recently adjusted in the project baseline). The
magnet is wound with $8$ individually insulated superconducting wires. The
wires are then connected in series in a joint box positioned at one end of the
magnet.

Superconducting magnets need impregnation for better strand-to-strand thermal
management and better mechanical integrity. This consists of filling the gaps
between windings with an electrically insulating material, improving the
overall thermal conductivity of the pack, and better distributing the stress
across it. Traditional epoxy-based impregnation often exhibits a training
phenomenon, requiring multiple cycles before optimal performance. Recent
studies demonstrate that paraffin wax impregnation eliminates this issue [3,
4, 5], solving a bothersome problem of superconducting magnets. Here, we
describe a dedicated impregnation setup developed for this purpose. Magnetic
measurements of the demonstrator were initially conducted at room temperature
(i.e. warm) [6] and are now presented at cryogenic temperatures (i.e. cold) to
assess superconducting magnet performance accurately. The cold tests were
conducted in Nov. 2023.

Table 1: Specification Final-focus Quadrupole Demonstrator

Aperture $\diameter$ $40$ mm Inductance $6.3$ mH $l_{z}$ $430$ mm Joint R
@$600$ A $130$ n$\Omega$ Conductor NbTi $I_{\text{nom}}$ per wire $764.5$ A
Winding pack $8$ wires $J_{e}$ $728$ A/mm2 Wire type LHC cable polyimid insul.
peak $B$ $2.95$ T Wire $\diameter$ w/ insul. $0.970$ mm peak $\nabla$ $B$
$100$ T/m Wire $\diameter$ w/o insul. $0.825$ mm $\int B_{2}\,dz$ at
r=$r^{\star}$ $0.315$ T.m s/c to Cu ratio $0.526$ $T_{\text{nom}}$ $1.9$ K
$l_{\text{conductor}}$ $276$ m $T_{\text{margin}}$ $2.5$ K

Figure 1: Wax impregnation station inside an oven (top); Cross-section view of
the top part of the magnet mounted within the set-up (bottom).

##  2 WAX IMPREGNATION

The magnet impregnation process involves injecting hot melted wax into the
empty gaps between windings, which then solidifies to ensure mechanical
integrity. Paraffin wax undergoes a liquid-solid phase transition at around 55
°C, contracting by approximately $15\%$ upon solidification. To prevent the
formation of trapped empty spaces within the winding pack, temperature
gradients are applied across the magnet during progressive cooling, following
an approach similar to [3].

Figure 2: Measured temperature at different positions across the magnet during
the wax solidification process.

Figure 1 shows the impregnation station and a cross-section of the top part of
the magnet. Internal wax reservoirs at the top of the magnet provide
additional impregnation volume to compensate for crystallization shrinkage.
The system is placed within an oven to achieve the temperature necessary for
the wax to be melted. The magnet is air-tight and a weak vacuum is created
within the magnet to avoid trapping air bubbles which would degrade the
impregnation quality. Wax flow is facilitated by a peristaltic pump, filling
the magnet from the bottom. A copper shaft is inserted inside the magnet and
connected to a copper plate lying at the bottom of it. This induces two
temperature gradients across the magnet: one from the inner layer to the outer
layer, and another from the top to the bottom part of the magnet. A water
cooling circuit is added below the magnet enhancing temperature control of the
heat exchanger and accelerating wax cooling.

Figure 2 shows the temperature profiles measured at various positions across
the magnet and downsampled by averaging over ten points. An average 3 °C and 5
°C temperature difference was observed between the inner side (in contact with
the heat exchanger shaft) and outer side (in contact with the magnet housing);
and between the top and bottom part of the magnet respectively.

##  3 TESTING CAMPAIGN AT COLD

Once impregnated, the magnet was tested at cryogenic temperatures for
characterization. The prototype was placed inside a cryostat with rotating
coils, as depicted on the right side of Fig. 3. The cryostat is filled with
helium and cooled to $1.9$ K and later to $4.5$ K. The test took place in
November 2023 at the Siegtag cryostat of the SM18 cryogenic facility of CERN.

Figure 3: (Left) Inner former of the CCT magnet. Note the deviation of the
first two turns from a pure sinusoidal shape [7]. (right) Open cryostat set-up
for cold measurements. Figure 4: Quenching campaign of the magnet with various
current ramping rates.

###  3.1 Critical Current Density

Multiple taps within the joint box monitor voltage between wires to detect
quenches and protect the magnet by discharging energy into a dump resistor. A
$50$ $\Omega$ crowbar connected to a $75$ $\Omega$ dump resistor is employed
for energy extraction. A $100$ mV threshold with $10$ ms validation time is
used for quench detection. Inductive voltages are canceled by comparing the
first half of the magnet with the second half, and quarters of the magnets are
also compared for redundancy and in the case of symmetric quenches.

Various ramping tests determine the critical current of the magnet, as
illustrated in Fig. 4. Notably, no training is observed. The short-sample
critical current already measured on previous occasions is reached on the
first attempt. A margin exceeding $20\%$ at $1.9$ K is validated for the
nominal operating current. Additional tests have been performed at $4.5$ K for
demonstration, also reaching the short sample limit associated with this
operating temperature.

Table 2: Comparison of the field quality between the warm and cold ($1.9$ K)
measurements of the CCT quadrupole powered with a current of $5$ A and $765$ A
respectively. The normal $b_{n}$ and skew $a_{n}$ components of different
higher order multipoles $n$ are evaluated at a $R_{\text{ref}}=10$ mm. They
are expressed in $10^{-4}$ units.

$\bm{b_{n}}$ $\bm{a_{n}}$ $\bm{n}$ Case Center Non Corrected Side Corrected
Side Center Non Corrected Side Corrected Side 3 Warm $-0.22$ $-35.98$ $1.83$
$0.32$ $-78.07$ $-2.24$ Cold $0.32$ $-32.04$ $0.71$ $0.11$ $-57.36$ $-0.78$ 4
Warm $0.53$ $11.80$ $3.38$ $0.54$ $-46.56$ $-1.15$ Cold $0.52$ $13.81$ $3.36$
$0.79$ $-37.39$ $-1.44$ 5 Warm $-0.16$ $-8.66$ $-0.30$ $-0.03$ $-5.83$ $0.02$
Cold $-0.18$ $-7.31$ $-0.30$ $-0.06$ $-4.13$ $0.18$ 6 Warm $0.64$ $4.44$
$1.02$ $-0.12$ $1.74$ $-0.11$ Cold $0.79$ $4.11$ $1.22$ $-0.07$ $1.80$ $-0.12$
7 Warm $0.03$ $-1.62$ $0.00$ $-0.00$ $-0.79$ $-0.01$ Cold $0.01$ $-1.47$
$-0.03$ $-0.01$ $-0.91$ $0.01$ 8 Warm $0.00$ $0.65$ $-0.05$ $-0.03$ $0.27$
$-0.02$ Cold $0.00$ $0.65$ $-0.01$ $-0.02$ $0.35$ $-0.04$

###  3.2 Field Quality Measurements

The field quality is determined by the normalized coefficients resulting from
the Fourier analysis of the flux density distribution across the magnet cross-
section [8, Chapter 6]. These coefficients are determined by inductive passive
PCB coils rotating around the axis inside the magnet. Three coils cover the
entire axial length: one at the center for nominal magnetic strength
characterization, and two others to investigate edge effects. A distinctive
conductor layout on one side corrects undesired higher-order multipole
content. This is illustrated in Fig. 3 with the grooves hosting the conductor
diverging from the usual sinusoidal profile that leads to undesired edge
effects. The correction is confirmed by comparing coefficients from the coils
positioned on the corrected and non-corrected sides as shown in TABLE 2. The
field quality analysis results, although worse than the warm tests, are still
below $1$ unit in the central region. We have not repeated the extra step done
at warm which reduced multipole errors and consisted of rotating the magnet by
about $45^{\circ}$ with respect to the measuring system (before repeating the
measurements) [6].

The magnitude of undesired higher-order harmonics is significantly reduced on
the corrected side compared to the non-corrected one, reaching values below
$5$ units (the normalization uses the quadrupole field of that side and not of
the whole magnet). One should note that when evaluating the overall
performance of the magnet, the total length would be considered, resulting in
a larger normalizing quadrupole component which mainly comes from the central
part. This dilutes the contribution of the edges to the field quality of the
prototype, and even more for the full-scale magnet.

This validates the possibility of tuning the conductor layout of CCT magnets
in the edge to remove undesired higher-order multipole content. This approach
can also get rid of crosstalk from the adjacent magnet in a twin aperture
arrangement as envisaged in FCC-ee.

##  4 CONCLUSION

This work presents the first prototype built for the final focus quadrupole of
the FCC-ee study. A CCT superconducting magnet impregnated with paraffin wax
is built and tested. The magnet exhibits no training behavior and reaches the
short-sample critical current ($1000$ A) for various current rates up to $100$
A/s. The field quality is excellent. Additionally, this prototype validates
the possibility of tuning the conductor layout of a CCT magnet to greatly
reduce edge effects for improved field quality. Future work will present a
similar approach for crosstalk compensation of a twin-aperture arrangement.
This will be demonstrated with a CCT sextupole wound with High-Temperature-
Superconducting (HTS) tapes.

##  5 ACKNOWLEDGMENTS

We would like to acknowledge the help of Michael Daly, Colin Müller, and Jaap
Kosse from the Paul Scherrer Institute (PSI), as well as Daniel Barna from the
Wigner Research Centre for Physics for their support with the impregnation
station. Additionally, many thanks to the team at SM18 Arnaud Devred, Gerard
Willering, Jerome Feuvrier, and Franco Julio Mangiarotti for their help with
the cold tests. Finally, we would like to thank Austin Ball, Günther
Dissertori, Markus Klute, and Maf Alidra for their invaluable help during
various stages of the project.

## References

  * [1] A. Abada et al., “FCC-ee: The Lepton Collider”, Eur. Phys. J. Spec. Top., vol. 228, pp. 261–623, 2019.
  * [2] G. A. Kirby et al., “Hi-Lumi LHC Twin Aperture Orbit Correctors 0.5-m Model Magnet Development and Cold Test”, IEEE Trans. Appl. Supercond., vol. 28, pp. 1–5, 2018.
  * [3] D. Barna et al., “The Superconducting Shield (SuShi) Septum Magnet Prototype”, IEEE Trans. Appl. Supercond., vol. 32, pp. 1–5, 2022.
  * [4] M. Daly et al., “Improved training in paraffin-wax impregnated Nb3Sn Rutherford cables demonstrated in BOX samples”, Sup. Sci. and Tech., vol. 35, p. 055014, 2022.
  * [5] D. M. Araujo et al., “Assessment of Training Performance, Degradation and Robustness of ParaffinWax Impregnated Nb 3 Sn Demonstrator under High Magnetic Field”, IEEE Trans. Appl. Supercond., vol. 34, pp. 1–8, 2024.
  * [6] M. Koratzinos et al., “Magnetic Measurements at Warm of the First FCC-ee Final Focus Quadrupole Prototype”, presented at Proc. IPAC’21, (virtual) Campinas, Brazil, May 2021.
  * [7] M. Koratzinos et al., “A Method for Greatly Reduced Edge Effects and Crosstalk in CCT Magnets”, IEEE Trans. Appl. Supercond., vol. 28, pp. 1–5, 2018.
  * [8] S. Russenschuck, Field Computation for Accelerator Magnets Analytical and Numerical Methods for Electromagnetic Design and Optimization, Wiley-VCH Verlag GmbH, 2010.


allow us to control the expected squared KSD between $\hat{\mu}^{(n)}$ and
$\hat{\nu}^{(n)}$.

To this end, let $\mathcal{E}$ denote the event that each random batch
$\mathcal{K}_{t}$ of GB-SVGD is disjoint and contains unique elements for
every $t\in(T)$. Subsequently, let $\mathcal{K}$ denote the set of all indices
that were chosen to be part of some random batch $\mathcal{K}_{t}$. Let
$\Lambda$ be a uniformly random permutation over $\\{0,\dots,n-1\\}$. We note
that, conditioned on $\mathcal{E}$, the distribution of the random set
$\mathcal{K}$ is the same as the distribution of
$\\{\Lambda(0),\dots,\Lambda(KT-1)\\}$. We can couple a uniformly random
permutation $\Lambda$ and $\mathcal{K}_{t}$ for $0\leq t\leq T$ such that
under the event $\mathcal{E}$,
$\mathcal{K}=\\{\Lambda(0),\dots,\Lambda(KT-1)\\}$ and
$\\{\Lambda(tK),\dots,\Lambda((t+1)K-1)\\}$ is the random batch
$\mathcal{K}_{t}$. Thus, under the event $\mathcal{E}$, one can couple a
uniformly random permutation $\Lambda$ and $\mathcal{K}_{t}$ for $t\in(T)$
such that $\mathcal{K}=\\{\Lambda(0),\dots,\Lambda(KT-1)\\}$ and
$\mathcal{K}_{t}=\\{\Lambda(tK),\dots,\Lambda((t+1)K-1)\\}$

With this insight, we couple VP-SVGD and GB-SVGDas follows. We note that, the
random batch $\mathcal{K}_{t}$ in GB-SVGD is sampled independently of the
initial particles. To this end, let
$\bar{\mathbf{x}}_{0}^{(0)},\dots,\bar{\mathbf{x}}_{0}^{(n-1)}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{0}$, and let the random batches
$\mathcal{K}_{t}$ and permutation $\Lambda$ be jointly distributed as
described above, independently of
$\bar{\mathbf{x}}_{0}^{(0)},\dots,\bar{\mathbf{x}}_{0}^{(n-1)}$, i.e.

$\displaystyle\Lambda\sim\mathsf{Uniform}(\mathbb{S}_{(n)}),\quad\mathcal{K}_{t}=\\{\Lambda(tK),\dots,\Lambda((t+1)K-1)\\},\quad
t\in(T)$

We now define $\mathbf{x}_{0}^{(0)},\dots,\mathbf{x}_{0}^{(KT+n-1)}$ as:

$\mathbf{x}_{0}^{(l)}\coloneqq\begin{cases}=\bar{\mathbf{x}}_{0}^{(\Lambda(l))}&\textrm{
for }0\leq l\leq n-1\\\ \sim\mu_{0}\text{ independent of everything
else}&\textrm{ for }n\leq l\leq KT+n-1\end{cases}$ (31)

Let $\bar{\mathbf{x}}_{0}^{(0)},\dots,\bar{\mathbf{x}}_{0}^{(n-1)}$ and
$\mathcal{K}_{t}$ as the initialization and random batches for GB-SVGD, and
let $\mathbf{x}_{0}^{(0)},\dots,\mathbf{x}_{0}^{(KT+n-1)}$ be the
initialization for GB-SVGD. We first show that this construction is indeed a
valid coupling between VP-SVGD and GB-SVGD.

######  Claim 1.

Conditioned on $\mathcal{E}$, the inputs to VP-SVGD and GB-SVGD, as
constructed above is a valid coupling, i.e., the marginal distribution of
$\mathbf{x}_{0}^{(0)},\dots,\mathbf{x}_{0}^{(KT+n-1)}$ is equal to the
distribution of initial particles in VP-SVGD, and the marginal distribution of
$\bar{\mathbf{x}}^{(0)}_{0},\dots,\bar{\mathbf{x}}^{(n-1)}_{0},(\mathcal{K}_{t})_{t\in(T)}$
is the same as the distribution of initial particles and random batches in
$\mathcal{K}_{t}$

###### Proof.

By construction
$\bar{\mathbf{x}}^{(0)}_{0},\dots,\bar{\mathbf{x}}^{(n-1)}_{0}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{0}$. Moreover, conditioned on
$\mathcal{E}$, the distribution of
$\mathcal{K}_{t}=\\{\Lambda(tK),\dots,\Lambda((t+1)K-1)\\}$, has the
distribution of a uniform random batch of size $K$ since
$\Lambda\sim\mathsf{Uniform}(\mathbb{S}_{n})$. Furthermore, since $\Lambda$ is
sampled independently of
$\bar{\mathbf{x}}^{(0)}_{0},\dots,\bar{\mathbf{x}}^{(n-1)}_{0}$,
$\mathcal{K}_{t}$ is independent of
$\bar{\mathbf{x}}^{(0)}_{0},\dots,\bar{\mathbf{x}}^{(n-1)}_{0}$ for any
$t\in(T)$. Thus, the coupling constructed above has the correct marginal with
respect to GB-SVGD.

To establish the same for VP-SVGD, we note that by (31),
$\mathbf{x}^{(n)}_{0},\dots,\mathbf{x}^{(KT+n-1)}_{0}\ \stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\mu_{0}$, _sampled independently of everything else_.
Moreover, since
$\bar{\mathbf{x}}^{(0)}_{0},\dots,\bar{\mathbf{x}}^{(n-1)}_{0}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{0}$, we infer that
$\bar{\mathbf{x}}^{(\Lambda(0))}_{0},\dots,\bar{\mathbf{x}}^{(\Lambda(n-1))}_{0}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{0}$ for any arbitrary
permutation $\Lambda\in\mathbb{S}_{n}$. From this, and (31), we conclude that
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT+n-1)}_{0}\ \stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\ \mu_{0}$. Hence, the coupling constructed above has the
correct marginal with respect to VP-SVGD. ∎

We now show that, under the constructed coupling, the time-evolution of the
particles of VP-SVGD and GB-SVGD satisfy
$\bar{\mathbf{x}}_{t}^{(\Lambda(l))}=\mathbf{x}_{t}^{(l)},\ KT\leq l\leq
n-1,t\in(T+1)$, when conditioned on the event $\mathcal{E}$.

######  Claim 2.

Let the inputs to VP-SVGD and GB-SVGD be coupled as per the construction
above. Then, conditioned on the event $\mathcal{E}$, the particles
$\mathbf{x}^{(s)}_{t}$ and $\bar{\mathbf{x}}^{(s)}_{t}$ of VP-SVGD and GB-SVGD
respectively, satisfy
$\bar{\mathbf{x}}_{t}^{(\Lambda(l))}=\mathbf{x}_{t}^{(l)}$ for every $KT\leq
l\leq n-1$ and $0\leq t\leq T$

###### Proof.

We prove this by an inductive argument. Clearly, the claim holds for $t=0$ by
the construction of our coupling. Assume it holds for some arbitrary
$t\in(T)$. Now, writing the update equation for GB-SVGD for $KT\leq l\leq
n-1$,

$\displaystyle\bar{\mathbf{x}}^{(\Lambda(l))}_{t+1}$
$\displaystyle=\bar{\mathbf{x}}^{(\Lambda(l))}_{t}-\frac{\gamma}{K}\sum_{r\in\mathcal{K}_{t}}h(\bar{\mathbf{x}}^{(\Lambda(l))}_{t},\bar{\mathbf{x}}^{(r)}_{t})$
$\displaystyle=\bar{\mathbf{x}}^{(\Lambda(l))}_{t}-\frac{\gamma}{K}\sum_{l=0}^{K-1}h(\bar{\mathbf{x}}^{(\Lambda(l))}_{t},\bar{\mathbf{x}}^{(\Lambda(Kt+l))}_{t})$
$\displaystyle=\mathbf{x}^{(l)}_{t}-\frac{\gamma}{K}\sum_{l=0}^{K-1}h(\mathbf{x}^{(l)}_{t},\mathbf{x}^{(Kt+l)}_{t})=\mathbf{x}^{(l+1)}_{t}$

where the second equality uses the fact that
$\mathcal{K}_{t}=\\{\Lambda(tK),\dots,\Lambda((t+1)K-1)\\}$ when conditioned
on $\mathcal{E}$ and the third equality uses the induction hypothesis
$\bar{\mathbf{x}}^{(\Lambda(l))}_{t}=\mathbf{x}^{(l)}_{t}$ for $KT\leq l\leq
n-1$. Hence, the claim is proven true by induction. ∎

Equipped with the above coupling between the inputs of VP-SVGD and GB-SVGD,
one can now couple their outputs by sampling an $S\sim\mathsf{Uniform}((n))$
and using this sampled $S$ as the random timestep chosen by both VP-SVGD (Step
6 in Algorithm 1) and GB-SVGD (Step 7 in Algorithm 2) that are run with the
coupled input constructed above. It is easy to see that this results in a
coupling of the outputs $\mathbf{Y}$ and $\bar{\mathbf{Y}}$ of VP-SVGD and GB-
SVGD respectively. Furthermore, by Claim 2, we note that, conditioned on the
event $\mathcal{E}$, $\mathbf{y}^{(l-TK)}=\bar{\mathbf{y}}^{(\Lambda(l))}$ for
every $KT\leq l\leq n-1$. We now define the permutation
$\tau\in\mathbb{S}_{(n)}$ as follows:

$\tau(\Lambda(l))=\begin{cases}l+n-KT&\text{ for }0\leq l\leq KT-1\\\
l-KT&\text{ for }KT\leq l\leq n-1\end{cases}$ (32)

It follows that $\bar{\mathbf{y}}^{\tau(l)}=\mathbf{y}^{(l)}$ for $KT\leq
l\leq n-1$. Thus, by definition of Kernel Stein Discrepancy (Definition 1), we
can infer that the following holds when conditioned on the event $\mathcal{E}$

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})|\mathcal{E}]$
$\displaystyle=\mathbb{E}\left[\|h_{\hat{\nu}^{(n)}}-h_{\hat{\mu}^{(n)}}\|^{2}_{\mathcal{H}}\
|\ \mathcal{E}\right]$
$\displaystyle=\mathbb{E}\left[\|\frac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\bar{\mathbf{y}}^{(l)})-\frac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{y}^{(l)})\|^{2}_{\mathcal{H}}\
|\ \mathcal{E}\right]$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}[\|\sum_{l=0}^{n-1}h(\cdot,\bar{\mathbf{y}}^{(\tau(l))})-h(\cdot,\mathbf{y}^{(l)})\|^{2}_{\mathcal{H}}\
|\ \mathcal{E}]$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}[\|\sum_{l=0}^{KT-1}h(\cdot,\bar{\mathbf{y}}^{(\tau(l))})-h(\cdot,\mathbf{y}^{(l)})\|^{2}_{\mathcal{H}}\
|\ \mathcal{E}]$
$\displaystyle\leq\frac{KT}{n^{2}}\sum_{l=0}^{KT-1}\mathbb{E}[\|h(\cdot,\bar{\mathbf{y}}^{(\tau(l))})-h(\cdot,\mathbf{y}^{(l)})\|_{\mathcal{H}}^{2}\
|\ \mathcal{E}]$ $\displaystyle\leq\frac{2K^{2}T^{2}\xi^{2}}{n^{2}}$ (33)

where the second step uses the permutation invariance of summation, the third
step uses the fact that $\bar{\mathbf{y}}^{\tau(l)}=\bar{\mathbf{y}}^{(l)}$
for $KT\leq l\leq n-1$, the fourth step uses the convexity of
$\|\cdot\|^{2}_{\mathcal{H}}$ and the last step uses the almost-sure iterate
bounds of Lemma 9 and 10

Under the event $\mathcal{E}^{c}$, we directly apply the almost-sure iterate
bounds of Lemma 9 and 10 to obtain the following:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})|\mathcal{E}^{c}]$
$\displaystyle=\mathbb{E}\left[\|h_{\hat{\mu}^{(n)}}-h_{\hat{\nu}^{(n)}}\|^{2}_{\mathcal{H}}\
|\ \mathcal{E}^{c}\right]$
$\displaystyle=\frac{1}{n^{2}}\mathbb{E}[\|\sum_{l=0}^{n-1}h(\cdot,\bar{\mathbf{y}}^{(l)})-h(\cdot,\mathbf{y}^{(l)})\|^{2}_{\mathcal{H}}\
|\ \mathcal{E}^{c}]$ $\displaystyle\leq 2\xi^{2}$ (34)

From Equations (D.1) and (D.1), it follows that:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})]$
$\displaystyle=\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})|\mathcal{E}]\mathbb{P}(\mathcal{E})+\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})|\mathcal{E}^{c}]\mathbb{P}(\mathcal{E}^{c})$
$\displaystyle\leq\frac{2K^{2}T^{2}\xi^{2}}{n^{2}}\mathbb{P}(\mathcal{E})+2\xi^{2}\mathbb{P}(\mathcal{E}^{c})$

Recall that $P(\mathcal{E})=1$ under sampling without replacement and
$P(\mathcal{E})=1-\tfrac{K^{2}T^{2}}{n}$ under sampling with replacement.
Thus, we conclude that the following holds under the constructed coupling of
$\mathbf{Y}$ and $\bar{\mathbf{Y}}$

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})]\leq\begin{cases}\frac{2K^{2}T^{2}\xi^{2}}{n^{2}}\quad&\text{(without
replacement sampling)}\\\
\frac{2K^{2}T^{2}\xi^{2}}{n^{2}}\left(1-\frac{K^{2}T^{2}}{n}\right)+\frac{2K^{2}T^{2}\xi^{2}}{n}\quad&\text{(with
replacement sampling)}\end{cases}$

∎

##  Appendix E Finite-Particle Convergence Guarantees for VP-SVGD and GB-SVGD

In this section, we proove the finite-particle convergence rates of VP-SVGD
and GB-SVGD. To this end, present the proof of Corollary 1 in Appendix E.1 and
the proof of Corollary 2 in Appendix E.2. Finally, we compare the oracle
complexity (i.e., the number of evaluations of $\nabla F$) of VP-SVGD and GB-
SVGD to that of SVGD in Appendix E.3

###  E.1 VP-SVGD : Proof of Corollary 1

###### Proof.

Recall from Algorithm 1 that the outputs of VP-SVGD are
$\mathbf{x}^{(KT)}_{S},\dots,\mathbf{x}^{(KT+n-1)}_{S}$ where
$S\sim\mathsf{Uniform}(\\{0,\dots,T-1\\})$. Hence, their empirical measure
$\hat{\mu}^{(n)}$ is given by
$\hat{\mu}^{(n)}=\tfrac{1}{n}\sum_{l=0}^{n-1}\delta_{\mathbf{x}^{(KT+l)}_{S}}$.
From the definition of the Kernel Stein Discrepancy (Definition 1), it follows
that,

$\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})=\|h_{\hat{\mu}^{(n)}}\|^{2}_{\mathcal{H}}=\|\tfrac{1}{n}\sum_{l=1}^{N}h(\cdot,\mathbf{x}^{(KT+l)}_{S})\|^{2}_{\mathcal{H}}$
(35)

For the sake of clarity, only in this proof, we use $\mathcal{C}$ to denote
the conditioning on the virtual particles
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(KT-1)}_{0}$. Now, consider any
arbitrary $t\in\\{0,\dots,T-1\\}$. Taking conditional expectations on both
sides of Equation (35) by conditioning on $\mathcal{C}$ and the event
$\\{S=t\\}$, we obtain the following:

$\displaystyle\mathbb{E}\left[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})\
|\ \mathcal{C},S=t\right]$
$\displaystyle=\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{S})\|^{2}_{\mathcal{H}}\
|\ \mathcal{C},S=t\right]$
$\displaystyle=\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\ (\mathbf{x}^{(s)}_{0})_{0\leq s\leq KT-1}\right]$ (36)

Recall from Equation (14) in Appendix C.1 that for any $l\in\\{0,\dots,n-1\\}$
$\mathbf{x}^{(KT+l)}_{t}$ depends only on
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}$ and
$\mathbf{x}^{(KT+l)}_{0}$. Furthermore, from Appendix C.1, we recall that the
filtration $\mathcal{F}_{t}$ is defined as
$\mathcal{F}_{t}=\sigma(\\{\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}\\})$.
It follows that,

$\displaystyle\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\ (\mathbf{x}^{(s)}_{0})_{0\leq s\leq KT-1}\right]$
$\displaystyle=\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=1}^{N}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\ (\mathbf{x}^{(s)}_{0})_{0\leq s\leq Kt-1}\right]$
$\displaystyle=\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\ \mathcal{F}_{t}\right]$ (37)

To control
$\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\ \mathcal{F}_{t}\right]$, we apply the arguments used in the proof of Lemma
3. To this end, note that when conditioned on the virtual particles
$\mathbf{x}^{(0)}_{0},\dots,\mathbf{x}^{(Kt-1)}_{0}$, the particles
$\mathbf{x}^{(KT)}_{t},\dots,\mathbf{x}^{(KT+n-1)}_{t}\
\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\ \mu_{t}|\mathcal{F}_{t}$.
Furthermore, since $\gamma\leq\nicefrac{{1}}{{2A_{1}L}}$ (as per the parameter
settings of Theorem 1),
$\|h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|_{\mathcal{H}}\leq\xi\ \forall\ l\in(n)$
by Lemma 6. Finally,
$\mathbb{E}[h(\mathbf{x},\mathbf{x}^{(KT+l)}_{t})|\mathcal{F}_{t}]=h_{\mu_{t}|\mathcal{F}_{t}}(\mathbf{x})\
\forall\ l\in(n),\ \mathbf{x}\in\mathbb{R}^{d}$. Hence, from Lemma 8, we
conclude that $h_{\mu_{t}|\mathcal{F}_{t}}$ is the Gelfand-Pettis integral of
the map $\mathbf{x}\to h(\mathbf{x},\mathbf{x}^{(KT+l)}_{t})$ with respect to
the measure $\mu_{t}|\mathcal{F}_{t}$, i.e.,

$\mathbb{E}[\left\langle
h(\cdot,\mathbf{x}^{(KT+l)}_{t}),f\right\rangle_{\mathcal{H}}|\mathcal{F}_{t}]=\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},f\right\rangle\ \forall\ f\in\mathcal{H}$ (38)

To control
$\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\ \mathcal{F}_{t}\right]$, we proceed as follows:

$\displaystyle\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}$
$\displaystyle=\frac{1}{n^{2}}\sum_{l_{1},l_{2}=0}^{n-1}\left\langle
h(\cdot,\mathbf{x}^{(Kt+l_{1})}_{t}),h(\cdot,\mathbf{x}^{(Kt+l_{2})}_{t})\right\rangle_{\mathcal{H}}$
$\displaystyle=\frac{1}{n^{2}}\sum_{l=0}^{m-1}\|h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\|^{2}_{\mathcal{H}}+\frac{1}{n^{2}}\sum_{0\leq
l_{1}\neq l_{2}\leq n-1}\left\langle
h(\cdot,\mathbf{x}^{(Kt+l_{1})}_{t}),h(\cdot,\mathbf{x}^{(Kt+l_{2})}_{t})\right\rangle_{\mathcal{H}}$
$\displaystyle\leq\frac{\xi^{2}}{n}+\frac{1}{n^{2}}\sum_{0\leq l_{1}\neq
l_{2}\leq n-1}\left\langle
h(\cdot,\mathbf{x}^{(Kt+l_{1})}_{t}),h(\cdot,\mathbf{x}^{(Kt+l_{2})}_{t})\right\rangle_{\mathcal{H}}$

where the last inequality uses the fact that
$\|h(\cdot,\mathbf{x}^{(Kt+l)}_{t})\|_{\mathcal{H}}\leq\xi$ almost surely as
per Lemma 9.

To control the conditional expectation of the off-diagonal terms, let
$i=Kt+l_{1}$ and $j=Kt+l_{2}$ for any arbitrary $l_{1},l_{2}$ with $0\leq
l_{1}\neq l_{2}\leq n-1$. Conditioned on $\mathcal{F}_{t}$,
$\mathbf{x}^{(i)}_{t}$ and $\mathbf{x}^{(j)}_{t}$ are i.i.d samples from
$\mu_{t}|\mathcal{F}_{t}$. Thus, by Equation (38) and Fubini’s Theorem,

$\displaystyle\mathbb{E}_{\mathbf{x}^{(i)}_{t},\mathbf{x}^{(j)}_{t}}\left[\left\langle
h(\cdot,\mathbf{x}^{(i)}_{t}),h(\cdot,\mathbf{x}^{(j)}_{t})\right\rangle_{\mathcal{H}}\bigr{|}\mathcal{F}_{t}\right]$
$\displaystyle=\mathbb{E}_{\mathbf{x}^{(i)}_{t}}\left[\mathbb{E}_{\mathbf{x}^{(j)}_{t}}\left[\left\langle
h(\cdot,\mathbf{x}^{(i)}_{t}),h(\cdot,\mathbf{x}^{(j)}_{t})\right\rangle_{\mathcal{H}}\bigr{|}\mathcal{F}_{t}\right]\right]$
$\displaystyle=\mathbb{E}_{\mathbf{x}^{(i)}_{t}}\left[\left\langle
h_{\mu_{t}|\mathcal{F}_{t}},h(\cdot,\mathbf{x}^{(i)}_{t})\right\rangle_{\mathcal{H}}\bigr{|}\mathcal{F}_{t}\right]$
$\displaystyle=\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}$

It follows that,

$\displaystyle\mathbb{E}\left[\|\tfrac{1}{n}\sum_{l=0}^{n-1}h(\cdot,\mathbf{x}^{(KT+l)}_{t})\|^{2}_{\mathcal{H}}\
|\
\mathcal{F}_{t}\right]\leq\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}+\frac{\xi^{2}}{n}$

Substituting the above into equation 36 and equation 37, we obtain the
following:

$\displaystyle\mathbb{E}\left[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})\
|\ \mathcal{C},S=t\right]$
$\displaystyle\leq\frac{\xi^{2}}{n}+\|h_{\mu_{t}|\mathcal{F}_{t}}\|^{2}_{\mathcal{H}}=\frac{\xi^{2}}{n}+\mathsf{KSD}^{2}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{*})$

where the second step applies Definition 1. Finally, taking expectations with
respect to $\mathcal{C}$ and $S\sim\mathsf{Uniform}(\\{0,\dots,T-1\\})$ on
both sides of the above inequality, we get:

$\displaystyle\mathbb{E}\left[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})\right]\leq\frac{\xi^{2}}{n}+\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\mu_{t}|\mathcal{F}_{t}||\pi^{\star})]$

Substituting the bound from Theorem 1 into the above inequality, we conclude
that:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq\frac{\xi^{2}}{n}+\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{K}$

We note that for $\gamma=O(\tfrac{(Kd)^{\eta}}{T^{1-\eta}})$ and
$R=\sqrt{\nicefrac{{d}}{{L}}}$,
$\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)=O(d)$
by Lemma 4 and

$\displaystyle\xi^{2}\leq 4\zeta_{0}+4\zeta_{1}(\gamma
T)^{\nicefrac{{2}}{{\alpha}}}+4\zeta_{2}(\gamma^{2}T)^{\nicefrac{{2}}{{\alpha}}}+4\zeta_{3}R^{\nicefrac{{4}}{{\alpha}}}\leq
O\left(\left(KdT\right)^{\tfrac{1}{1+\alpha}}+d^{\nicefrac{{2}}{{\alpha}}}\right)$

Furthermore,

$\displaystyle\frac{2\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{K}$ $\displaystyle\leq O\left(\frac{d}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{2K}\right)\leq
O\left(\frac{d^{1-\eta}}{(KT)^{\eta}}\right)$ $\displaystyle\leq
O\left(\frac{d^{\tfrac{2+\alpha}{2(1+\alpha)}}}{(KT)^{\tfrac{\alpha}{2(1+\alpha)}}}\right)$

It follows that,

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq
O\left(\frac{d^{\nicefrac{{2}}{{\alpha}}}}{n}+\frac{(KTd)^{\tfrac{1}{1+\alpha}}}{n}+\frac{d^{\tfrac{2+\alpha}{2(1+\alpha)}}}{(KT)^{\tfrac{\alpha}{2(1+\alpha)}}}\right)$

$KT=d^{\tfrac{\alpha}{2+\alpha}}n^{\tfrac{2(1+\alpha)}{2+\alpha}}$, we
conclude:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq
O\left(\frac{d^{\tfrac{2}{2+\alpha}}}{n^{\tfrac{\alpha}{2+\alpha}}}+\frac{d^{\nicefrac{{2}}{{\alpha}}}}{n}\right)$

∎

###  E.2 GB-SVGD : Proof of Corollary 2

###### Proof.

Let $\bar{\mathbf{Y}}=(\bar{\mathbf{y}}^{(0)},\dots,\bar{\mathbf{y}}^{(n-1)})$
denote the $n$ particles output by GB-SVGD and let
$\hat{\nu}^{(n)}=\tfrac{1}{n}\sum_{l=0}^{n-1}\delta_{\bar{\mathbf{y}}^{(l)}}$
denote their empirical measure. Let $\mathcal{E}$ denote the event that each
random batch $\mathcal{K}_{t}$ of GB-SVGD is disjoint and contains unique
elements for every $t\in(T)$. Moreover, let
$\mathbf{Y}=(\mathbf{y}^{(0)},\dots,\mathbf{y}^{(n-1)})$ denote the $n$
particles output by VP-SVGD, run with the parameter settings stated above, and
coupled with $\bar{\mathbf{Y}}$ as per the coupling constructed in the proof
of Theorem 3 in Appendix D.1. Let
$\hat{\mu}^{(n)}=\tfrac{1}{n}\sum_{l=0}^{n-1}\delta_{\mathbf{y}^{(l)}}$ denote
their empirical measure. By definition of Kernel Stein Discrepancy (Definition
1) and the convexity $\|\cdot\|^{2}_{\mathcal{H}}$, it follows that:

$\displaystyle\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]$
$\displaystyle=\mathbb{E}[\|h_{\hat{\nu}^{(n)}}\|^{2}_{\mathcal{H}}]$
$\displaystyle=\mathbb{E}[\|h_{\hat{\nu}^{(n)}}-h_{\hat{\mu}^{(n)}}+h_{\hat{\mu}^{(n)}}\|^{2}_{\mathcal{H}}]$
$\displaystyle\leq
2\mathbb{E}[\|h_{\hat{\nu}^{(n)}}-h_{\hat{\mu}^{(n)}}\|^{2}_{\mathcal{H}}]+2\mathbb{E}[\|h_{\hat{\mu}^{(n)}}\|^{2}_{\mathcal{H}}]$
$\displaystyle=2\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\hat{\mu}^{(n)})]+2\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]$

Substituting the bounds of Theorem 3 and Corollary 1 into the above
inequality, we conclude the following:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\frac{4K^{2}T^{2}\xi^{2}}{n^{2}}\mathbb{P}(\mathcal{E})+4\xi^{2}\mathbb{P}(\mathcal{E}^{c})+\frac{2\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$

We recall that, $\mathbb{P}(\mathcal{E})=1$ under without-replacement sampling
of the random batches $\mathcal{K}_{t}$ and
$\mathbb{P}(\mathcal{E})=1-\nicefrac{{K^{2}T^{2}}}{{n}}$ under with-
replacement sampling. Thus, under without-replacement sampling, the following
holds:

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\frac{4K^{2}T^{2}\xi^{2}}{n^{2}}+\frac{2\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$

Moreover, the following holds under with-replacement sampling

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\frac{4K^{2}T^{2}\xi^{2}}{n^{2}}\left(1-\frac{K^{2}T^{2}}{n}\right)+\frac{4K^{2}T^{2}\xi^{2}}{n}+\frac{2\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$

Now, let us consider GB-SVGD without replacement with
$R=\sqrt{\nicefrac{{d}}{{L}}}$, $\gamma=O(\tfrac{(Kd)^{\eta}}{T^{1-\eta}})$
and $KT=n^{\nicefrac{{1}}{{2}}}$ It follows that
$\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)=O(d)$
by Lemma 4 and

$\displaystyle\xi^{2}$ $\displaystyle\leq 4\zeta_{0}+4\zeta_{1}(\gamma
T)^{\nicefrac{{2}}{{\alpha}}}+4\zeta_{2}(\gamma^{2}T)^{\nicefrac{{2}}{{\alpha}}}+4\zeta_{3}R^{\nicefrac{{4}}{{\alpha}}}$
$\displaystyle\leq
O\left(\left(KdT\right)^{\tfrac{1}{1+\alpha}}+d^{\nicefrac{{2}}{{\alpha}}}\right)$
$\displaystyle\leq
O\left(d^{\nicefrac{{2}}{{\alpha}}}+d^{\tfrac{1}{1+\alpha}}n^{\tfrac{1}{2(1+\alpha)}}\right)$

Furthermore,

$\displaystyle\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$ $\displaystyle\leq O\left(\frac{d}{\gamma
T}+\frac{\gamma B(4+L)\xi^{2}}{2K}\right)\leq
O\left(\frac{d^{1-\eta}}{(KT)^{\eta}}\right)$ $\displaystyle\leq
O\left(\frac{d^{\tfrac{2+\alpha}{2(1+\alpha)}}}{n^{\tfrac{\alpha}{4(1+\alpha)}}}\right)$

Hence, we conclude that,

$\displaystyle\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]$
$\displaystyle\leq\frac{4K^{2}T^{2}\xi^{2}}{n^{2}}+\frac{2\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$
$\displaystyle\leq\frac{6\xi^{2}}{n}+\frac{4\mathsf{KL}\left(\mu_{0}|\mathcal{F}_{0}\boldsymbol{||}\pi^{\star}\right)}{\gamma
T}+\frac{2\gamma B(4+L)\xi^{2}}{K}$ $\displaystyle\leq
O\left(\frac{d^{\nicefrac{{2}}{{\alpha}}}}{n}+\frac{d^{\tfrac{1}{1+\alpha}}}{n^{\tfrac{1+2\alpha}{2(1+\alpha)}}}+\frac{d^{\tfrac{2+\alpha}{2(1+\alpha)}}}{n^{\tfrac{\alpha}{4(1+\alpha)}}}\right)$

∎

###  E.3 Oracle Complexity of SVGD, VP-SVGD and GB-SVGD

We now compare the gradient oracle complexity, (i.e., the number of
evaluations of $\nabla F$) of VP-SVGD (as implied by Corollary 1) and GB-SVGD
(as implied by Corollary 2) with that of SVGD as implied by the state-of-the-
art finite particle guarantee of Shi and Mackey [42].

####  E.3.1 SVGD

From Equation (1), We note that $T$ steps of SVGD run with $n$ particles
requires $n^{2}T$ evaluations of $\nabla F$.

##### Subgaussian $\pi^{\star}$

For subgaussian $\pi^{\star}$, the finite-particle convergence rate obtained
by Shi and Mackey [42] is
$\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}_{\mathsf{SVGD}}||\pi^{\star})=\tilde{O}(\tfrac{\mathsf{poly}(d)}{\sqrt{\log\log
n^{\Theta(\nicefrac{{1}}{{d}})}}})$, where $\hat{\mu}^{(n)}_{\mathsf{SVGD}}$
denotes the empirical measure of the $n$ particles output by SVGD. By
carefully following the analysis of Shi and Mackey [42], we infer that, to
achieve
$\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}_{\mathsf{SVGD}}||\pi^{\star})\leq\epsilon$,
SVGD requires $T=\tilde{O}(\tfrac{\mathsf{poly}(d)}{\epsilon^{2}})$ and
$n=\tilde{O}(\exp(\Theta(de^{\tfrac{\mathsf{poly}(d)}{\epsilon^{2}}})))$. Thus
the oracle complexity of SVGD (as implied by Shi and Mackey [42]) for
achieving
$\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}_{\mathsf{SVGD}}||\pi^{\star})$ is
$\tilde{O}(\tfrac{\mathsf{poly}(d)}{\epsilon^{2}}\cdot\exp(\Theta(de^{\tfrac{\mathsf{poly}(d)}{\epsilon^{2}}})))$

####  E.3.2 VP-SVGD

From Algorithm 1, we note that $T$ steps of VP-SVGD run with $n$ particles and
a batch-size of $K$ requires $K^{2}T^{2}+KTn$ evaluations of $\nabla F$.

##### Subgaussian $\pi^{\star}$

For subgaussian $\pi^{\star}$, Corollary 1 implies a finite-particle
convergence rate of
$\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]=O(\tfrac{d^{\nicefrac{{1}}{{2}}}}{n^{\nicefrac{{1}}{{2}}}}+\tfrac{d}{n})$
(where $\hat{\mu}^{(n)}$ denotes the empirical measure of the $n$ particles
output by VP-SVGD) assuming
$KT=d^{\nicefrac{{1}}{{2}}}n^{\nicefrac{{3}}{{2}}}$. Hence, to achieve
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq\epsilon$,
VP-SVGD requires $n=O(\tfrac{d}{\epsilon^{4}})$ and
$KT=d^{\nicefrac{{1}}{{2}}}n^{\nicefrac{{3}}{{2}}}=\tfrac{d^{2}}{\epsilon^{6}}$.
The resulting oracle complexity for achieving
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq\epsilon$
is $O(\tfrac{d^{4}}{\epsilon^{12}})$. Compared to the oracle complexity of
SVGD obtained above, this is a _double exponential improvement in both $d$ and
$\nicefrac{{1}}{{\epsilon}}$_. Notably, the obtained oracle complexity
guarantee _completely eliminates the curse of dimensionality_.

##### Subexponential $\pi^{\star}$

For subexponential $\pi^{\star}$, Corollary 1 implies a finite-particle
convergence rate of
$\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]=O(\tfrac{d^{\nicefrac{{2}}{{3}}}}{n^{\nicefrac{{1}}{{3}}}}+\tfrac{d^{2}}{n})$
(where $\hat{\mu}^{(n)}$ denotes the empirical measure of the $n$ particles
output by VP-SVGD) assuming
$KT=d^{\nicefrac{{1}}{{3}}}n^{\nicefrac{{4}}{{3}}}$. Hence, to achieve
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq\epsilon$,
VP-SVGD requires $n=O(\tfrac{d^{2}}{\epsilon^{6}})$ and
$KT=d^{\nicefrac{{1}}{{3}}}n^{\nicefrac{{4}}{{3}}}=\tfrac{d^{3}}{\epsilon^{8}}$.
The resulting oracle complexity for achieving
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\mu}^{(n)}||\pi^{\star})]\leq\epsilon$
is $O(\tfrac{d^{6}}{\epsilon^{16}})$.

####  E.3.3 GB-SVGD

From Algorithm 2, we note that $T$ steps of GB-SVGD run with $n$ particles and
a batch-size of $K$ requires $KTn$ evaluations of $\nabla F$.

##### Subgaussian $\pi^{\star}$

For subgaussian $\pi^{\star}$, Corollary 2 implies a finite-particle
convergence rate of
$\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]=O(\tfrac{d^{\nicefrac{{2}}{{3}}}}{n^{\nicefrac{{1}}{{6}}}}+\tfrac{d}{n})$
(where $\hat{\nu}^{(n)}$ denotes the empirical measure of the $n$ particles
output by GB-SVGD) assuming $KT=n^{\nicefrac{{1}}{{2}}}$. Hence, to achieve
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\epsilon$,
GB-SVGD requires $n=\tfrac{d^{4}}{\epsilon^{12}}$ and
$KT=\sqrt{n}=\tfrac{d^{2}}{\epsilon^{6}}$. Under this setting, the oracle
complexity of GB-SVGD as implied by Corollary 2 is
$O(\tfrac{d^{6}}{\epsilon^{18}})$. Compared to the oracle complexity of SVGD
obtained above, this is _a double exponential improvement in both $d$ and
$\nicefrac{{1}}{{\epsilon}}$_. Notably, the obtained oracle complexity
guarantee _completely eliminates the curse of dimensionality_

##### Subexponential $\pi^{\star}$

For subexponential $\pi^{\star}$, Corollary 2 implies a finite-particle
convergence rate of
$\mathbb{E}[\mathsf{KSD}^{2}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]=O(\tfrac{d^{\nicefrac{{3}}{{4}}}}{n^{\nicefrac{{1}}{{8}}}}+\tfrac{d^{2}}{n})$
(where $\hat{\nu}^{(n)}$ denotes the empirical measure of the $n$ particles
output by GB-SVGD) assuming $KT=n^{\nicefrac{{1}}{{2}}}$. Hence, to achieve
$\mathbb{E}[\mathsf{KSD}_{\pi^{\star}}(\hat{\nu}^{(n)}||\pi^{\star})]\leq\epsilon$,
GB-SVGD requires $n=\tfrac{d^{6}}{\epsilon^{16}}$ and
$KT=\sqrt{n}=\tfrac{d^{3}}{\epsilon^{8}}$. Under this setting, the oracle
complexity of GB-SVGD as implied by Corollary 2 is
$O(\tfrac{d^{9}}{\epsilon^{24}})$.

##  Appendix F Literature Review

Initial works on the analysis of SVGD such as Liu [29], Lu et al. [34], Duncan
et al. [13], Chewi et al. [6], Nüsken and Renger [37] consider the continuous-
time population limit, i.e., the limit of infinite particles and vanishing
step-sizes. In this regime, Liu [29], Lu et al. [34], Nüsken and Renger [37]
show that the behavior of SVGD is characterized by a Partial Differential
Equation (PDE), and established asymptotic convergence of this PDE to the
target distribution. The work of Duncan et al. [13] proposes the Stein
Logarithmic Sobolev Inequality which ensures exponential convergence of this
PDE to the target distribution. However, characterizing the conditions under
which this inequality holds is an open problem. The work of Chewi et al. [6]
show that the PDE governing SVGD in the continuous-time population limit can
be interpreted as an approximate Wasserstein gradient flow of the Chi-squared
divergence. To this end, Chewi et al. [6] shows that the (exact) Wasserstein
gradient flow of the Chi-squared divergence exhibits exponential convergence
to the target distribution when $\pi^{\star}$ satisfies a Poincare Inequality.
To the best of our knowledge, the first discrete-time non-asymptotic
convergence result for population-limit SVGD was established in Korba et al.
[25], where the authors interpreted population-limit SVGD as projected
Wasserstein gradient descent. Their result relied on the assumption that the
Kernel Stein Discrepancy to the target is uniformly bounded along the
trajectory of SVGD, a condition which is hard to verify apriori. This result
was significantly improved in Salim et al. [41], which established convergence
of population-limit SVGD assuming the potential $F$ is smooth the target
$\pi^{\star}\propto e^{-F}$ satisfies Talagrand’s inequality $\mathsf{T}_{1}$,
an assumption which is equivalent to subgaussianity of $\pi^{\star}$. This
result was extended in Sun et al. [44] to accommodate for potentials $F$ that
satisfy a more general smoothness condition.

In comparison to prior works on population-limit SVGD, the literature on
finite-particle SVGD is relatively sparse. The works of Liu [29] and Gorham et
al. [17] establish that the dynamics of finite-particle SVGD asymptotically
converge to that of population-limit SVGD in bounded Lipschitz distance and
Wasserstein-1 distance respectively, as the number of particles approaches
infinity. Under the stringent condition of bounded $F$ (which is violated in
various scenarios, e.g. log-strongly concave $\pi^{\star}$), Korba et al. [25]
derived a non-asymptotic bound between the expected Wasserstein-2 distance
between finite-particle SVGD and population-limit SVGD. The work of Liu et al.
[32] considers the special case of Gaussian initialization and bilinear
kernels, and provides a finite-particle convergence guarantee for the first
and second moments whenever the target density is also a Gaussian. To the best
of our knowledge, Shi and Mackey [42] is the only prior work that explicitly
establishes a complete non-asymptotic convergence guarantee of finite-particle
SVGD for subgaussian targets, and shows that the empirical measure of SVGD run
with $n$ particles converges to the target density in KSD at a rate of
$O\left(\sqrt{\tfrac{\mathsf{poly}(d)}{\log\log
n^{\Theta(\nicefrac{{1}}{{d}})}}}\right)$.

##  Appendix G Additional Experimental Details

As discussed in Section 8, we benchmark SVGD and GB-SVGD on the task of
Bayesian Logistic Regression using the Covertype dataset, obtained from the
UCI Machine Learning Repository, whch contains around $580,000$ datapoints
with dimension $d=54$. We use $n=100$ particles for both SVGD and GB-SVGD, and
use batch-size $K=40$ for GB-SVGD. For both SVGD and GB-SVGD, we use AdaGrad
with momentum to set step-sizes as per the procedure used by Liu and Wang
[30]. Our experiments were performed using Python 3 on a 2.20 GHz Intel Xeon
CPU with 13 GB of memory.


[figure]subcapbesideposition=top

# Combined effect of mutually frequency-detuned strong and weak drives on a
two-level system: Envelope of the Rabi oscillations

M. E. Raikh  Department of Physics and Astronomy, University of Utah, Salt
Lake City, UT 84112

###### Abstract

Near-resonant ac-drive acting on a two-level system induces the Rabi
oscillations of the level occupations. It is shown that additional weak drive
properly frequency-detuned from the primary drive causes a resonant response.
This response manifests itself in the emergence of the envelope of the
oscillations. At resonance, the inverse period of the envelope is proportional
to the amplitude of the weak drive. The resonant condition reads: difference
of frequencies between the two drives is equal to the ac-splitting of
quasilevels in the field of the strong drive. Technically, the resonance can
be inferred from the analogy between the equations for the time-evolution of
the spin amplitude and the Mathieu equation, which describes e.g. the
parametric resonance.

##  I Introduction

Dynamics of a spin placed in a magnetic field, ${\bf B}={\bf z_{0}}B$, and
subjected to a sinusoidal ac drive, ${\bf b}=2b{\bf x}_{0}\cos\omega t$, is
known in tiniest details since the celebrated papers [Rabi, ; Bloch, ; Autler-
Townes, ; Shirley, ; Mollow, ]. Most prominent is the situation when $\omega$
is close to the Zeeman splitting, $\Delta_{Z}=\gamma B$, of the spin levels.
Here $\gamma$ is the gyromagnetic ratio. Then the probabilities of $\uparrow$
and $\downarrow$ spin projections oscillate with the generalized Rabi
frequency

$\Theta_{R}=\Bigl{[}\left(\Delta_{Z}-\omega\right)^{2}+\Omega_{R}^{2}\Bigr{]}^{1/2},$
(1)

where $\Omega_{R}=\gamma b$. Concerning the $x$ and $y$ spin projections,
their dynamics comprises three frequencies:Mollow  $\omega$, and
$\omega\pm\Theta_{R}$ (Mollow triplet).

Four underlying frequencies of the spin dynamics suggest the way in which the
driven spin can be manipulated by the secondary ac field. Namely, the
frequency of the secondary field, $\omega_{m}$, can either be low, of the
order of $\Omega_{R}\ll\omega$ (see e.g. Refs. [Fedaruk1, ; Fedaruk2, ;
Fedaruk3, ; Fedaruk4, ; Glenn, ]), or it can be highTHEORY ; WeakAndStrong ;
Experiment1 ; Experiment2 ; Experiment3 ; EXPERIMENTcold ; TwoDrives ; numeric
; numeric1  and detuned from $\omega$ approximately by $\Omega_{R}$. Within
the first technique, the modification of the primary Rabi oscillations by the
secondary drive depends strongly on the relation between the magnitude and
frequency of the secondary drive. In particular, strong secondary drive slows
down the Rabi oscillationsGlenn . When the secondary drive is weak, it can
affect the Rabi oscillations only under the resonance condition
$\omega_{m}\approx\Omega_{R}$. Under this condition, primary oscillations
develop an envelope.Fedaruk3 ; Glenn  In the language of quasienergies this
effect can be interpreted as follows. While the primary drive causes the ac-
splitting of the bare levels into the quasilevels, under the resonant
secondary drive each of these quasilevels gets split additionally.

Within the second technique, the effect of two near-resonant drives on the
two-level system is much more dramatic. Theoretical studiesTHEORY ;
WeakAndStrong ; TwoDrives ; numeric ; numeric1  indicate that adding a weak
second (probe) drive can covert the regime of the Rabi oscillations induced by
the primary drive (pump) into a chaotic behavior of the observables. The
underlying reason for chaos is that, in the presence of two drives, the
Floquet theorem is lifted (if the pump and probe frequencies are
incommensurate). Meanwhile, the behavior of fluorescence observed in
experiments on bichromatically driven two-level systemsExperiment1 ;
Experiment2 ; Experiment3 ; EXPERIMENTcold  indicate that the effect of two
drives is much less dramatic than in theory. In the time domain, the probe
drive leads to the additional modulation of the Rabi oscillations. In the
frequency domain, the pump leads to the splitting of each individual peak of
the Mollow triplet into the daughter triplets.

In the theoretical papers on bichromatic drive the results are presented in
the form of numerical curves calculated for certain sets of parameters. This
is certainly justified for revealing of chaos. But in the regimes when the
effect of secondary drive is modest, there are no analytical formulas
predicting the modified behavior of the Rabi oscillations for given amplitudes
and frequencies of two drives. Such formulas are derived in the present paper.
The reason why these formulas can be derived is that, within the rotating wave
approximation (RWA), the Floquet theorem applies not for the primary drive
frequency but rather for the difference between the frequencies of two drives.
Our analytical treatment yields the resonance condition for which the effect
of a weak secondary drive on the Rabi oscillations is most pronounced.

Figure 1: (Color online) Schematic illustration of a two-level system in a
magnetic field driven by the ac fields with frequencies $\omega_{1}$ and
$\omega_{2}$ close to the level separation, $\Delta_{Z}$. Strong drive with
magnitude, $\lambda_{1}$, and ac-splitting of the Zeeman levels, which
manifests itself in the Rabi oscillations of the level
$\lambda_{2}\ll\lambda_{1}$ properly detuned from the primary drive, see Eq.
14, leads to an additional splitting, $s_{1}$, of the ac-split levels, which,
in turn, results in the envelope of the Rabi oscillations with a period
$\frac{2\pi}{s_{1}}$.

##  II Evolution of the amplitudes of the spin components

Denote with $a_{1}$ and $a_{2}$ the amplitudes of $\uparrow$ and $\downarrow$
spin components. In the presence of two drives these amplitudes evolve with
time as

$\displaystyle
i\frac{da_{1}}{dt}=\frac{\Delta_{Z}}{2}a_{1}+\left[\lambda_{1}e^{-i\omega_{1}t}+\lambda_{2}e^{-i\omega_{2}t}\right]a_{2},$
(2) $\displaystyle
i\frac{da_{2}}{dt}=-\frac{\Delta_{Z}}{2}a_{2}+\left[\lambda_{1}e^{i\omega_{1}t}+\lambda_{2}e^{i\omega_{2}t}\right]a_{1},$
(3)

where $\omega_{1}$ and $\omega_{2}$ are the frequencies of the primary and
secondary drives, respectively, while $\lambda_{1}$ and $\lambda_{2}$ are the
amplitudes of these drives in the frequency units. For simplicity we neglect
the phase shift between the two drives, and, thus, choose $\lambda_{1}$,
$\lambda_{2}$ to be real. We also assume that both drives are weak, namely,
$\lambda_{1},\lambda_{2}\ll\leavevmode\nobreak\ \Delta_{Z}$. This justifies
the replacement of $\cos\omega_{1}t$ and $\cos\omega_{2}t$ in Eqs. 2 and 3 by
the exponents, which is the essence of RWA. Incorporating of counter-rotating
exponents amounts to renormalization of $\Delta_{Z}$ by the Bloch-Siegert
shiftBloch  of the order of
$\frac{\lambda_{1}^{2}}{\Delta_{Z}},\frac{\lambda_{2}^{2}}{\Delta_{Z}}\ll\Delta_{Z}$.
Following the standard procedure for a single-frequency drive, we introduce
new variables

$A_{1}=a_{1}\exp\left(\frac{i\Delta_{Z}}{2}t\right),\leavevmode\nobreak\
\leavevmode\nobreak\ A_{2}=a_{2}\exp\left(-\frac{i\Delta_{Z}}{2}t\right),$ (4)

which allows to exclude high frequency, $\Delta_{Z}$, from the system Eqs. 2
and 3, which takes the form

$\displaystyle
i\frac{dA_{1}}{dt}=\left[\lambda_{1}e^{-i\delta_{1}t}+\lambda_{2}e^{-i\delta_{2}t}\right]A_{2},$
(5) $\displaystyle
i\frac{dA_{2}}{dt}=\left[\lambda_{1}e^{i\delta_{1}t}+\lambda_{2}e^{i\delta_{2}t}\right]A_{1},$
(6)

where the detunings $\delta_{1}$ and $\delta_{2}$ are defined as

$\delta_{1}=\Delta_{Z}-\omega_{1},\leavevmode\nobreak\ \leavevmode\nobreak\
\delta_{2}=\Delta_{Z}-\omega_{2}.$ (7)

Next we reduce the system Eqs. 5, 6 to a single second-order differential
equation. Expressing $A_{2}$ from Eq. 5 and substituting into Eq. 6, we get

$\displaystyle-\frac{d^{2}A_{1}}{dt^{2}}-if(t)\frac{dA_{1}}{dt}$
$\displaystyle=\left[\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{1}\lambda_{2}e^{i(\delta_{1}-\delta_{2})t}+\lambda_{2}\lambda_{1}e^{-i(\delta_{1}-\delta_{2})t}\right]A_{1}.$
(8)

where the function $f(t)$ is given by

$\displaystyle
f(t)=\frac{\lambda_{1}\delta_{1}e^{i\delta_{1}t}+\lambda_{2}\delta_{2}e^{i\delta_{2}t}}{\lambda_{1}e^{i\delta_{1}t}+\lambda_{2}e^{i\delta_{2}t}}$
$\displaystyle=\frac{\delta_{1}+\delta_{2}}{2}+\frac{\delta_{1}-\delta_{2}}{2}\Biggl{(}\frac{1-\frac{\lambda_{2}}{\lambda_{1}}e^{i(\delta_{2}-\delta_{1})t}}{1+\frac{\lambda_{2}}{\lambda_{1}}e^{i(\delta_{2}-\delta_{1})t}}\Biggr{)}.$
(9)

It is seen from the expression Eq. II that the function $f(t)$ is periodic
with a period $2\pi/(\delta_{2}-\delta_{1})$. Equally, the right-hand side of
Eq. IV is periodic with the same period. We thus conclude that the Floquet
theorem applies to Eq. (IV). To boost the similarity to the Mathieu equation,
we exclude the first derivative from Eq. (IV) by introducing the new variable

$\displaystyle{\tilde{A}}_{1}(t)=A_{1}(t)\exp\Biggl{[}\frac{i}{2}\int_{0}^{t}dt^{\prime}f(t^{\prime})\Biggr{]}$
$\displaystyle=A_{1}(t)\Biggl{[}1+\frac{\lambda_{2}}{\lambda_{1}}e^{i(\delta_{2}-\delta_{1})t}\Biggr{]}^{1/2}\exp\Biggl{(}{\frac{i\delta_{1}t}{2}}\Biggr{)}.$
(10)

Substituting Eq. II into Eq. IV we arrive to the following equation for
${\tilde{A}}_{1}$

$\displaystyle-\frac{d^{2}{{\tilde{A}}_{1}}}{dt^{2}}-\frac{{\tilde{A}}_{1}}{4}\Biggl{\\{}\delta_{1}^{2}+\delta_{2}^{2}-\Biggl{(}\frac{\lambda_{2}\delta_{1}e^{i\delta_{2}t}+\lambda_{1}\delta_{2}e^{i\delta_{1}t}}{\lambda_{2}e^{i\delta_{2}t}+\lambda_{1}e^{i\delta_{1}t}}\Biggr{)}^{2}\Biggr{\\}}$
$\displaystyle=\left[\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{1}\lambda_{2}e^{i(\delta_{1}-\delta_{2})t}+\lambda_{2}\lambda_{1}e^{-i(\delta_{1}-\delta_{2})t}\right]{\tilde{A}}_{1}.$
(11)

Note that in the limit when the secondary drive is absent, $\lambda_{2}=0$,
Eq. II reproduces the Rabi oscillations in the field of the primary drive,
namely, the solutions of Eq. II have the form
${\tilde{A}}_{1}\propto\exp\left[\pm
i\left(\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}\right)^{1/2}t\right].$

##  III Weak secondary drive

Within the adopted RWA approximation Eq. II applies for arbitrary relation
between two drives $\lambda_{1}$ and $\lambda_{2}$ and arbitrary relation
between the detunings $\delta_{1}$ and $\delta_{2}$ as long as they are much
smaller than $\omega$. Next we take the limit of a weak secondary drive
$\lambda_{2}\ll\lambda_{1}$. In the lowest order in the secondary drive we can
neglect the term $\lambda_{2}^{2}$ in the right-hand side. Concerning the
left-hand side, we expand the ratio in the brackets as follows

$\displaystyle\frac{\lambda_{1}\delta_{2}e^{i\delta_{1}t}+\lambda_{2}\delta_{1}e^{i\delta_{2}t}}{\lambda_{1}e^{i\delta_{1}t}+\lambda_{2}e^{i\delta_{2}t}}$
$\displaystyle\approx\delta_{2}\Biggl{[}1+\frac{\lambda_{2}}{\lambda_{1}}\Bigl{(}\frac{\delta_{1}}{\delta_{2}}-1\Bigr{)}e^{i(\delta_{2}-\delta_{1})t}\Biggr{]}.$
(12)

With the above simplifications Eq. II takes the form

$\displaystyle-\frac{d^{2}{{\tilde{A}}_{1}}}{dt^{2}}-\Biggl{[}\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}+\frac{\lambda_{2}}{2\lambda_{1}}\delta_{2}(\delta_{1}-\delta_{2})e^{-i(\delta_{1}-\delta_{2})t}\Biggr{]}{\tilde{A}}_{1}$
$\displaystyle=\left[\lambda_{1}\lambda_{2}e^{i(\delta_{1}-\delta_{2})t}+\lambda_{2}\lambda_{1}e^{-i(\delta_{1}-\delta_{2})t}\right]{\tilde{A}}_{1}.$
(13)

There are two terms proportional to $\lambda_{2}$ in Eq. III. Assuming that
$\delta_{1}$ and $\delta_{2}$ are of the same order, the term in the left-hand
side can be estimated as $\lambda_{2}\delta_{1}^{2}/\lambda_{1}$, while the
term in the right-hand side is of the order of $\lambda_{2}\lambda_{1}$. If
$\lambda_{1}$ is much greater than $\delta_{1}$, which corresponds to the
limit of developed Rabi oscillations, the term in the right-hand side
dominates. Note now that, if one neglects the term in the left-hand side, Eq.
III reduces to the classical Mathieu equation, which describes e.g. the
electron motion in a weak one-dimensional potential. The time in Eq. III plays
the role of coordinate. The sum $\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}$
plays the role of energy. Despite the similarity, there is a dramatic
difference between Eq. III and the Mathieu equation, since the Mathieu
equation describes the ”forbidden gaps” or the domains of instabilities, while
Eq. III does not. Still, as we will see below, the position of resonance
captured by Eq. III can be found from the same ”Bragg condition” as the
position of the center of the forbidden gap, namely

$2\Biggl{(}\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}\Biggr{)}^{1/2}=|\delta_{1}-\delta_{2}|.$
(14)

To trace how the resonance Eq. 14 emerges, we search for the solution of Eq.
III in the Floquet form

${\tilde{A}}_{1}(t)=f\exp{\\{ist\\}}+g\exp{\Big{\\{}i(s+\delta_{1}-\delta_{2})t\Big{\\}}}.$
(15)

Substituting this form into Eq. III and equating the coefficients in front of
the two exponents, we arrive to the system

$\displaystyle\left(\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}-s^{2}\right)f=-\lambda_{1}\lambda_{2}\Bigg{[}\frac{\delta_{2}(\delta_{1}-\delta_{2})}{2\lambda_{1}^{2}}+1\Bigg{]}g,$
$\displaystyle\left(\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}-(s+\delta_{1}-\delta_{2})^{2}\right)g=-\lambda_{1}\lambda_{2}f.$
(16)

The brackets in the left-hand sides of the two equations coincide under the
condition: $s+\delta_{1}-\delta_{2}=-s$. Thus we set

$s=\frac{\delta_{2}-\delta_{1}}{2}+s_{1},$ (17)

and assume that $s_{1}\ll|\delta_{2}-\delta_{1}|$. Upon neglecting the
$s_{1}^{2}$-term, the system Eq. III takes the form

$\displaystyle\Bigg{[}\lambda_{1}^{2}-\frac{\delta_{2}^{2}-2\delta_{2}\delta_{1}}{4}-s_{1}(\delta_{2}-\delta_{1})\Bigg{]}f=$
$\displaystyle-\lambda_{1}\lambda_{2}\Bigg{[}\frac{\delta_{2}(\delta_{1}-\delta_{2})}{2\lambda_{1}^{2}}+1\Bigg{]}g,$
$\displaystyle\Bigg{[}\lambda_{1}^{2}-\frac{\delta_{2}^{2}-2\delta_{2}\delta_{1}}{4}+s_{1}(\delta_{2}-\delta_{1})\Bigg{]}g=$
$\displaystyle-\lambda_{1}\lambda_{2}f.$ (18)

Finally, upon multiplying the two last equations, we obtain the expression for
the Floquet exponent

$\displaystyle s_{1}^{2}\left(\delta_{2}-\delta_{1}\right)^{2}=$
$\displaystyle\left[\lambda_{1}^{2}-\frac{\delta_{2}^{2}}{4}+\frac{\delta_{1}\delta_{2}}{2}\right]^{2}-\lambda_{2}^{2}\left[\frac{\delta_{1}\delta_{2}-\delta_{2}^{2}}{2}+\lambda_{1}^{2}\right].$
(19)

Note now that the first bracket on the right-hand side of Eq. III turns to
zero when Eq. 14 is satisfied, i.e. at resonance. It is seen from Eq. III that
the resonance manifests itself via the anomalous sensitivity of the Floquet
exponent to the magnitude, $\lambda_{2}$, of the secondary drive. It is also
instructive to rewrite Eq. III in a different form by using the resonant
condition Eq. 14. Substituting
$\lambda_{1}^{2}=\frac{1}{4}\left(\delta_{2}^{2}-2\delta_{1}\delta_{2}\right)$
into Eq. III we get

$s_{1}=\pm\frac{\lambda_{2}\delta_{2}}{2(\delta_{2}-\delta_{1})}=\pm\frac{\lambda_{2}\delta_{2}}{4\left(\lambda_{1}^{2}+\frac{\delta_{2}^{2}}{4}\right)^{1/2}}.$
(20)

In the second identity we have again used the resonant condition Eq. 14.

Yet another way to cast Eq. III is to express $\delta_{2}$ from the resonant
condition. This yields

$s_{1}=\pm\frac{\lambda_{2}}{4\Theta_{R}}\left(\Delta_{Z}-\omega_{1}\pm
2\Theta_{R}\right).$ (21)

In the last identity we have returned to the original notations, generalized
Rabi frequency and the detuning of the primary drive (see Eq. 1). Since we
assumed that $s_{1}$ is much smaller than the frequency spacing between the
two drives, the criterion of applicability of Eq. 21 is
$\lambda_{2}\ll|\delta_{2}-\delta_{1}|$.

Eq. 21 is our main result. It suggests that, upon tuning the frequency of the
secondary drive to the resonance Eq. 14, the Floquet exponent, describing the
drive-induced modification of the Rabi oscillations, increases linearly with
the amplitude of the secondary drive.

##  IV Discussion

(i). Upon increasing of $\lambda_{2}$ the spectrum of the secondary Rabi
oscillations becomes progressively richer due to emergence of the higher-order
Bragg resonances for which the detunings satisfy the condition

$2\Biggl{(}\lambda_{1}^{2}+\frac{\delta_{1}^{2}}{4}\Biggr{)}^{1/2}=N|\delta_{1}-\delta_{2}|$
(22)

with higher $N$. To trace the emergence of these resonances, we expand the
ratio Eq. III to the second order in the ratio
$\frac{\lambda_{2}}{\lambda_{1}}$

$\displaystyle\Biggl{(}\frac{\lambda_{2}\delta_{1}e^{i\delta_{2}t}+\lambda_{1}\delta_{2}e^{i\delta_{1}t}}{\lambda_{2}e^{i\delta_{2}t}+\lambda_{1}e^{i\delta_{1}t}}\Biggr{)}^{2}\approx\delta_{2}^{2}\Biggl{[}1+2\frac{\lambda_{2}}{\lambda_{1}}\Biggl{(}\frac{\delta_{1}}{\delta_{2}}-1\Biggr{)}e^{i(\delta_{2}-\delta_{1})t}$
$\displaystyle+\Biggl{(}\frac{\lambda_{2}}{\lambda_{1}}\Biggr{)}^{2}\Biggl{(}\frac{\delta_{1}}{\delta_{2}}-1\Biggr{)}\Biggl{(}\frac{\delta_{1}}{\delta_{2}}-3\Biggr{)}e^{2i(\delta_{2}-\delta_{1})t}\Biggr{]}.$
(23)

We see that the new term in the expansion gives rise to the second harmonics
in the ”effective potential” which causes the resonance with $N=2$.

(ii). Floquet exponent translates into an envelope. The situation addressed in
the present paper is standard and was previously addressed in the literature,
see e.g. Refs. THEORY, ; WeakAndStrong, ; Experiment1, ; Experiment2, ;
Experiment3, ; EXPERIMENTcold, ; TwoDrives, ; numeric, ; numeric1, . We
considered a two-level system driven by two ac fields with strongly different
amplitudes. A new finding reported in the present paper is the ”Bragg
resonance” Eq. 14 which takes place when the difference of frequencies of the
drives is equal to the ac-splitting of quasienergies in the field of the
strong drive. Interplay of the two drives can be interpreted as a modulation
of the amplitude of a primary drive with a frequency
$\omega_{1}-\omega_{2}=\delta_{2}-\delta_{1}$, which couples the ac-split
levels. The result of this coupling is the secondary splitting of
quasienergies. Most importantly, at resonance, the magnitude of this secondary
splitting is linear in the amplitude of the weak drive. Emergence of the Bragg
resonance is illustrated schematically in the figure. Physically, the Bragg
resonance manifests itself in the modulation of the primary Rabi oscillations
with a frequency, $s_{1}$, proportional to the amplitude, $\lambda_{2}$, of
the weak drive, see Eq. 21. It should be emphasized that the depth of
modulation is independent of $\lambda_{2}$. Indeed, this depth is determined
by the interference of the amplitudes $f$ and $g$ in Eq. 15. At resonance, the
ratio of these amplitudes can be expressed from the system Eq. III as follows

$\frac{f}{g}=-\frac{s_{1}(\delta_{2}-\delta_{1})}{\lambda_{1}\lambda_{2}}.$
(24)

Substituting the Floquet exponent from Eq. 20 we find

$\frac{f}{g}=\pm\frac{\delta_{2}}{2\lambda_{1}},$ (25)

i.e. the amplitude $\lambda_{2}$ drops out. In other words, no matter how weak
is the secondary drive, at resonance, it lead s to the modulation of the
primary Rabi oscillations with a depth $\sim 1$. If the detuning,
$\delta_{2}$, of the secondary drive does not satisfy the Bragg condition, Eq.
14 by a small quantity, $\Delta\delta_{2}$, the modulation of the Rabi
oscillations gradually vanishes. As can be seen from the system Eq. III, the
Floquet exponent changes with $\Delta\delta_{2}$ as

$s_{1}(\Delta\delta_{2})=\Bigl{[}s_{1}(0)^{2}+\frac{(\Delta\delta_{2})^{2}}{4}\Bigr{]}^{1/2},$
(26)

i.e. similarly to the primary Rabi oscillations. This suggests that the
dependence of the Floquet exponent on the magnitude to the secondary drive
gradually transforms from linear to quadratic.

Note that the resonance Eq. 14 was not uncovered in the previous theoretical
studies. The reason is that these studies attempted to incorporate finite
lifetimes of the levels from the very beginning. Then, the $2\times 2$ system
Eqs. 2, 3 transforms into the $4\times 4$ system for the elements of the
density matrix, which necessarily complicates the analytical treatment.

(iii). Clearly, the solution of Eq. II exhibits anomalous behavior in a
particular case when the amplitudes of two drives are equal:
$\lambda_{1}=\lambda_{2}=\lambda$, while their frequencies are detuned
symmetrically with respect to the inter-level distance, $\omega$, i.e.
$\delta_{2}=-\delta_{1}=\delta$. In this particular case Eq. II takes the form

$-\frac{d^{2}{{\tilde{A}}_{1}}}{dt^{2}}-\frac{{\tilde{A}}_{1}}{4}\Biggl{(}\delta^{2}+16\lambda^{2}\cos^{2}\delta
t+\frac{\delta^{2}}{\cos^{2}\delta t}\Biggr{)}=0.$ (27)

From the form of Eq. 27 it can be concluded that the solution is singular in
the vicinity of the time moments, $t_{n}=\frac{(2n+1)\pi}{2\delta}$, where it
can be viewed as a Schrödinger equation in the attractive potential
$-\frac{1}{4(t-t_{n})^{2}}$. The form of the singular solution is
${\tilde{A}}_{1}\propto(t-t_{n})^{1/2}$. Note that, within the RWA, Eq. 27 is
exact. The only small-time scale, which can ”cure” the singularity is
$\sim\Delta_{Z}^{-1}\ll\lambda^{-1},\delta^{-1}$ is related to the violation
of the RWA. Certainly, any asymmetry, $(\Delta\lambda)$, in the amplitudes of
the two drives would cut off the singularity at
$(t-t_{n})\lesssim\frac{(\Delta\lambda)}{\lambda}$.

## References

  * (1) I. Rabi, Phys. Rev. 51, 652 (1937).
  * (2) F. Bloch and A. Siegert, ”Magnetic Resonance for Nonrotating Fields,” Phys. Rev. 57, (1940).
  * (3) S. H. Autler and C. H. Townes, ”Stark Effect in Rapidly Varying Fields,” Phys. Rev. 100, 703 (1955).
  * (4) J. H. Shirley, “Solution of the Schrödinger equation with a hamiltonian periodic in time,” Phys. Rev. 138, B979 (1965).
  * (5) B. R. Mollow, ”Power Spectrum of Light Scattered by Two-Level Systems,” Phys. Rev. 188, 1969 (1969).
  * (6) A. P. Saiko, R. Fedaruk, and S. A. Markevich, “Relaxation, decoherence and steady-state population inversion in qubits doubly dressed by microwave and radiofrequency fields,” J. Phys. B: At. Mol. Opt. Phys. 47, 155502 (2014).
  * (7) A. P. Saiko, R. Fedaruk, and A. Kolasa, “Higher-order geometric phase for qubits in a bichromatic field,” J. Phys. B: At. Mol. Opt. Phys. 45, 235501 (2012).
  * (8) A. P. Saiko, R. Fedaruk, “Multiplication of Qubits in a Doubly Resonant Bichromatic Field,” Journal of Experimental and Theoretical Physics Letters 91, 681 (2010).
  * (9) A. P. Saiko, R. Fedaruk, A. Kolasa, and S. A. Markevich, “Dissipative dynamics of qubits driven by a bichromatic field in the dispersive regime,” Phys. Scr. 85, 045301 (2012).
  * (10) R. Glenn, M. E. Limes, B. Pankovich, B. Saam, and M. E. Raikh, “Magnetic resonance in slowly modulated longitudinal field: Modified shape of the Rabi oscillations,” Phys. Rev. B 87, 155128 (2013).
  * (11) M. Wilkens and K. Rzazewski, ”Resonance fluorescence of an arbitrarily driven two-level atom,” Phys. Rev. A 40, 3164 (1989).
  * (12) Z. Ficek and H. S. Freedhoff, ”Fluorescence and absorption by a two-level atom in a bichromatic field with one strong and one weak component,” Phys. Rev. A 53, 4275 (1996).
  * (13) Q. Wu, D. J. Gauthier, and T. W. Mossberg, ”Phase-sensitive dynamics of bichromatically driven two-level atoms,” Phys. Rev. A 49, R1519(R) (1994).
  * (14) Q. Wu, D. J. Gauthier, and T. W. Mossberg, ”Optical resonance and coherent transients in dressed atomic systems,” Phys. Rev. A 50, 1474 (1994).
  * (15) C. C. Yu, J. R. Bochinski, T. M. V. Kordich, T. W. Mossberg, and Z. Ficek, ”Driving the driven atom: Spectral signatures,” Phys. Rev. A 56, R4381(R) (1997).
  * (16) M. Yan, E. G. Rickey, and Y. Zhu, ”Observation of doubly dressed states in cold atoms,” Phys. Rev. A 64, 013412 (2001).
  * (17) H. S. Freedhoff and Z. Ficek, ”Resonance fluorescence and Autler-Townes spectra of a two-level atom driven by two fields of equal frequencies,” Phys. Rev. A 55, 1234 (1997).
  * (18) J. Eidson and F. Fox, ”Quantum chaos in a two-1evel system in a semiclassical radiation field,” Phys. Rev. A 34, 3288 (1986).
  * (19) Y. Pomeau, B. Dorizzi, and B. Grammaticos, ”Chaotic Rabi Oscillations under Quasiperiodic Perturbation,” Phys. Rev. Lett. 56, 681 (1986).