# Nonadiabatic dissociation of molecular Bose-Einstein condensates:
competition between chemical reactions

Rajesh K. Malla  Center for Nonlinear Studies and Theoretical Division, Los
Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

We provide a framework to solve generic models describing the dissociation of
multiple molecular Bose-Einstein condensates in a nonadiabatic regime. The
competition between individual chemical reactions can lead to non-trivial
dependence on critical components such as path interference and symmetries,
thus, affecting the final distribution of atomic population. We find an
analytical solution for an illustrative example model involving four atomic
modes. When the system parameters satisfy $CPT$ symmetry, where $C$ is charge
conjugation, $P$ is parity, and $T$ is time-reversal symmetry, our solution
predicts a population imbalance between atomic modes that is exponentially
sensitive to system parameters. However, a weakly broken symmetry alters the
population in each atomic mode and can reverse the population imbalance. Our
solution also demonstrates a strong quantum correlation between atomic modes
that leads to the spontaneous production of atoms in a multi-mode squeezed
state. Moreover, in our framework, a time-dependent non-Hermitian quantum
mechanics naturally manifests which can alternatively be realized
experimentally in photonic systems.

##  I Introduction

Recent experimental advances in ultracold atomic and molecular gases enable us
to investigate many-body precision physics, with potential application in
quantum control of chemical reactions, precision measurements, quantum
simulation and quantum information processingBohn  _et al._ [2017a]; Carr  _et
al._ [2009]; QuÃ©mÃ©ner and Julienne [2012]. Ultracold platforms allow control
of interaction parameters that has lead to some fascinating experimental
observations including, reaction between atomic and molecular Bose Einstein
condensates (BEC) Zhang  _et al._ [2021], realization of Unruh radiation Hu
_et al._ [2019], and Bose FireworksClark  _et al._ [2017]. Over the past two
decades, there has been a tremendous surge in the theoretical Yurovsky  _et
al._ [2002]; Altland  _et al._ [2009]; Sun and Sinitsyn [2016]; Itin and Törmä
[2009]; Malla  _et al._ [2021] and experimental Herbig  _et al._ [2003]; Xu
_et al._ [2003]; Bohn  _et al._ [2017b]; Matyjaśkiewicz  _et al._ [2008]; Duda
_et al._ [2021] studies on reaction between ultracold molecules and atoms near
Feshbach resonance. There have also been several studies on atom-molecule
conversion by Feshbach resonance due to coupled-channels Chin  _et al._
[2010]. However, the time-dependent sweeping across the Feshbach resonance to
trigger conversion between atomic and molecular BEC has been rarely explored
for systems involving multiple atomic modes.

The goal of this paper is to provide a frame work to solve a dissociation
mechanism in which different types of diatomic molecules undergo dissociation
simultaneously, so that the individual reactions are coupled. The dynamics now
involves multiple parameters and their interplay leads to various non-trivial
effects such as interference and symmetry breaking, which may emerge due to
competition between multiple modes corresponding to different types of atoms.
Multiple modes can also describe higher degrees of freedom, i.e., spin angular
momentum, rotational and vibrational modes.

The general Hamiltonian describing reaction between atomic and molecular
condensates is given by

${\hat{H}}=\sum_{n}\epsilon_{n}^{\Psi}{\hat{\Psi}}_{n}^{\dagger}{\hat{\Psi}}_{n}+\sum_{k,n}\epsilon_{k,n}^{a}{\hat{a}}_{k,n}^{\dagger}{\hat{a}}_{k,n}+\epsilon_{k,n}^{b}{\hat{b}}_{k,n}^{\dagger}{\hat{b}}_{k,n}\\\
+\sum_{k,n,m}g_{k,n,m}{\hat{\Psi}}_{n}^{\dagger}{\hat{a}}_{k,m}{\hat{b}}_{k,n-m}+g_{k,n,m}^{*}{\hat{\Psi}}_{n}{\hat{a}}_{k,m}^{\dagger}{\hat{b}}_{k,n-m}^{\dagger},$
(1)

where $k$ is the reaction channel Malla  _et al._ [2022],
${\hat{\psi}}({\hat{\psi}}^{\dagger})$, ${\hat{a}}({\hat{a}}^{\dagger})$, and
${\hat{b}}({\hat{b}}^{\dagger})$ represent annihilation (creation) of the
molecular field operators and atomic field operators, $n$ corresponding to a
particular molecular mode. Hamiltonian (1), in general, is not analytically
solvable. We propose to solve this model in the mean-field approximation where
the molecular operators can be replaced by their expectation values Yurovsky
_et al._ [2002]. For this approximation to hold, we assume that the number of
atoms produced is small compared to the number of molecules. We propose to
solve these type of models in the Heisenberg picture, where the atomic field
operators satisfy a non-Hermitian evolution equation.

There are two types of dissociation processes, stimulated dissociation and
spontaneous dissociation. In the stimulated process an initially populated
atomic mode (or modes) stimulates the dissociation of molecules. In the
spontaneous process the atomic modes are initially empty and the molecules
spontaneously dissociate into atoms. The spontaneous process could also
describe cosmological production of particles and anti-particles due to vacuum
fluctuations in the early stage of the Universe, which may have been a
consequence of fundamental parameters such as mass or interaction coefficients
becoming time-dependent and the evolution passing through resonance Green and
Porto [2020].

After laying out the framework, we solve an illustrative example containing
diatomic molecules with two atoms and each atom having two degrees of freedom;
thus, the model includes four atomic modes. Let us denote the two degrees of
freedom as spin-up and spin-down, however, note that they are bosons. We
predict several key results for parameters satisfying the $CPT$ symmetry.

First, the analytical result predicts an exponential sensitive imbalance of
atomic population in each mode. For example, an atomic mode initially occupied
with an spin-up atom stimulates an exponentially suppressed atomic production
in spin-down states. This suppression is an artifact of destructive
interference which arises due to the $CPT$ symmetry. However, a weakly broken
symmetry alters the atomic production in each mode and can reverse the spin
imbalance for a finite interference.

Our second key result appears in the spontaneous process, where the chemical
reactions are coupled, and the atoms are produced in a multi-mode squeezed
state. Multi-mode squeezed states promise applications in quantum metrology
Gessner  _et al._ [2018], quantum communication, and quantum imaging Coelho
_et al._ [2009]; Simon  _et al._ [2007]; Duan  _et al._ [2001].

In addition to these two central results, our solution could open pathways for
new investigations in time-dependent non-Hermitian systems.

##  II Theory of single dissociation process in the nonadiabatic regime

Let us start with the dissociation process of a molecular BEC into two-mode
atomic BEC, $AB\rightarrow A+B$. The dissociation process is described by the
Hamiltonian in the curve crossing approximation,Yurovsky  _et al._ [2002];
Kayali and Sinitsyn [2003]

${\hat{H}}_{2}=\mu_{1}(t){\hat{a}}^{\dagger}{\hat{a}}+\mu_{2}(t){\hat{b}}^{\dagger}{\hat{b}}+J{\hat{\psi}}^{\dagger}{\hat{a}}{\hat{b}}+J^{*}{\hat{\psi}}{\hat{a}}^{\dagger}{\hat{b}}^{\dagger}.$
(2)

where the chemical potentials, $\mu_{i}(t)$, of two atomic modes are
considered to time-dependent. In the nonadiabtic limit, only a small fraction
of molecules get converted into atoms Altland  _et al._ [2009]. Assuming the
number of molecules to be much larger than the number of atoms produced, we
can replace the molecular field operator with the expectation value
$\langle{\hat{\psi}}\rangle$Yurovsky  _et al._ [2002]; Kayali and Sinitsyn
[2003]. The system then reduces to interacting two atomic modes and the
corresponding Hamiltonian reads

${\hat{H}}_{eff}(t)=\mu_{1}(t){\hat{a}}^{\dagger}{\hat{a}}+\mu_{2}(t){\hat{b}}^{\dagger}{\hat{b}}+g{\hat{a}}{\hat{b}}+g^{*}{\hat{a}}^{\dagger}{\hat{b}}^{\dagger}.$
(3)

where $g=J\langle\psi^{\dagger}\rangle$.

The model (3) can be solved in the Heisenberg picture where the operators
${\hat{a}}(t)$ and ${\hat{b}}(t)$ satisfy

$\displaystyle i\frac{d}{dt}\begin{pmatrix}{\hat{a}}(t)\\\
{\hat{b}}^{\dagger}(t)\end{pmatrix}=\begin{pmatrix}\mu_{1}(t)&g^{*}\\\
-g&\mu_{2}(t)\end{pmatrix}\begin{pmatrix}{\hat{a}}(t)\\\
{\hat{b}}^{\dagger}(t)\end{pmatrix}.$ (4)

The molecular dissociation occurs near the crossing of two chemical potentials
where the chemical potentials can be linearlized. Also, note that equation (4)
is only applicable for bosonic atomic modes.

The number of atoms in $A$ mode at time $t\rightarrow\infty$ is then given by

$\langle{\hat{a}}^{\dagger}(t){\hat{a}}(t)\rangle_{t\rightarrow\infty}=n_{sp}^{(A)}+n_{st}^{(A)},$
(5)

where $n_{sp}$ corresponds to spontaneous dissociation while $n_{st}$
represents stimulated dissociation. For a linear sweep, $\mu_{1}(t)=\beta
t,\mu_{2}(t)=-\beta t$, equation (4) resembles a non-Hermitian Landau-Zener
(nLZ) model, and can be solved exactly (see appendix A). The general solution
of ${\hat{a}}(t)$ has the form
${\hat{a}}(t)=\phi_{A}(t){\hat{a}}(t_{0})+\phi_{B}(t){\hat{b}}^{\dagger}(t_{0})$,
where $t_{0}$ represents the initial time and is assumed to be far away from
the level crossing. The asymptotic solution of $\phi_{A}(t\rightarrow\infty)$
and $\phi_{B}(t\rightarrow\infty)$, with initial condition
$\phi_{A}(t_{0}\rightarrow-\infty)=1$, $\phi_{B}(t_{0}\rightarrow-\infty)=0$,
leads to the expressions

$\displaystyle|\phi_{A}(\infty)|^{2}=|\phi_{B}(\infty)|^{2}+1=\exp{\left(\frac{\pi|g|^{2}}{\beta}\right)}.$
(6)

In contrast to the Hermitian LZ model where
$|\phi_{A}(\infty)|^{2}+|\phi_{B}(\infty)|^{2}=1$, the non-Hermiticity
conserves the difference $|\phi_{A}(\infty)|^{2}-|\phi_{B}(\infty)|^{2}=1$.

The number of atoms in the spontaneous dissociation is
$n_{sp}^{(A)}=|\phi_{B}(\infty)|^{2}$. In the spontaneous dissociation
process, unpairing is quantum correlated. The atoms are produced in a two-mode
squeezed state, where the maximum uncertainty along one direction is
exponentially suppressed with the exponent proportional to $|g|^{2}/|\beta|$.

The number of stimulated atoms in mode $A$ depends on the initial population
of $A$ and $B$ modes and the phases of the solutions $\phi_{A}(t)$ and
$\phi_{B}(t)$ at large times Kayali and Sinitsyn [2003]. If only the $A$ mode
is initially populated with a single atom, then the number of atoms in the $A$
mode is $|\phi_{A}(\infty)|^{2}$ and in the $B$ mode is
$|\phi_{B}(\infty)|^{2}$. Since the atoms are always produced in pairs the
population difference between $n_{st}^{(A)}$ and $n_{st}^{(B)}$ is conserved.

Figure 1: (a) The schematic diagram of dissociation of diatomic molecules with
four atomic modes. (b) The instantaneous eigenvalues of (9) are shown as a
function of time for parameters
$\beta_{1}=-1,\beta_{2}=0.5,E_{1}=5,E_{2}=1,g=0.5,\gamma=1$. The black curves
are real values and the red curves represent imaginary values. (c) The number
of atoms produced in a particular mode when the system was stimulated by a
single atom in $A_{\uparrow}$ mode. The scattered points are obtained from
numerical evolution of (9) while the solid curves are obtained from the
analytical formulas in (12).

##  III Generalization to dissociation of multiple molecular BECs

Now let us focus on a generic scenario where multiple molecular condensates
undergo dissociation process. Multiple molecular condensates can consist of
different types of molecules. These complex dynamics is difficult to handle.
Here we propose a way to solve such a complex dynamics in the curve crossing
approximation. The Hamiltonian describing dissociation of multiple molecular
condensates consisting of diatomic molecules are given by

${\hat{H}}_{n}(t)=\sum\limits_{i=1}^{n}\mu_{1}^{(i)}(t){\hat{a}}_{i}^{\dagger}{\hat{a}}_{i}+\sum\limits_{i=1}^{m}\mu_{2}^{(j)}(t){\hat{b}}_{i}^{\dagger}{\hat{b}}_{i}j\\\
+\sum\limits_{ij}\left[J_{ij}{\hat{\psi}}_{ij}^{\dagger}{\hat{a}}_{i}{\hat{b}}_{j}+h.c.\right],$
(7)

where $\mu_{1}^{(i)}(t)$ and $\mu_{2}^{(i)}(t)$ are the time-dependent
chemical potentials of two-atomic modes for $i$-th type molecule, parameters
$g_{i}$ are the coupling constants of the dissociation process
$AB_{ij}\rightarrow A_{i}+B_{i}$. The simultaneous dissociation of different
types of molecules are correlated.

Hamiltonian (7) can describe dissociation of $n\times m$ types of different
molecules. Employing similar procedures as before, if one replaces the
molecular operators with their complex numbers then the system reduces many
quadratically interacting atomic modes. We can solve model (7) in the
Heisenberg picture in a nonadiabatic regime. The solution requires solving an
equation similar to (4), and the corresponding matrix satisfies non-Hermitian
multi-state Landau-Zener (nMLZ). The Hermitian multi-state Landau-Zener (MLZ)
models have been explored extensively, while the non-Hermitian counterpart has
not been studied at all. In the next section we solve an illustrative example
model of Hamiltonian (7) that includes four atomic modes.

###  III.1 Example: dissociation process including four atomic modes

Now let us model a dissociation process, when different molecular condensates
have same two atoms and each atom can be in different rotational or
vibrational states. For simplicity, we restrict the atomic modes to two
levels, and denote the levels as spin-up ($\uparrow$) and spin-down
($\downarrow$) states, see Fig. 1a. Then the Hamiltonian describing the
molecular dissociation process reads

$\displaystyle{\hat{H}}_{4}=\mu_{1}(t)({\hat{a}}_{\uparrow}^{\dagger}{\hat{a}}_{\uparrow}-{\hat{a}}_{\downarrow}^{\dagger}{\hat{a}}_{\downarrow})+\mu_{2}(t)({\hat{b}}_{\uparrow}^{\dagger}{\hat{b}}_{\uparrow}-{\hat{b}}_{\downarrow}^{\dagger}{\hat{b}}_{\downarrow})$
$\displaystyle+\Big{[}J_{1}{\hat{\psi}}_{\uparrow\uparrow}^{\dagger}{\hat{a}}_{\uparrow}{\hat{b}}_{\uparrow}+J_{1}{\hat{\psi}}_{\downarrow\downarrow}^{\dagger}{\hat{a}}_{\downarrow}{\hat{b}}_{\downarrow}+J_{2}{\hat{\psi}}_{\uparrow\downarrow}^{\dagger}{\hat{a}}_{\uparrow}{\hat{b}}_{\downarrow}-J_{2}{\hat{\psi}}_{\downarrow\uparrow}^{\dagger}{\hat{a}}_{\downarrow}{\hat{b}}_{\uparrow}+h.c.\Big{]},$
(8)

where the coupling parameters describing dissociation of triplet
($\uparrow\uparrow$) or ($\downarrow\downarrow$) molecules are symmetric while
the coupling parameters describing dissociation of singlet
($\uparrow\downarrow$) or ($\downarrow\uparrow$) molecules are anti-symmetric.
This anti-symmetric nature is similar to a spin-orbit coupling Lin  _et al._
[2011].

We solve (8) in the mean-field approximation and replace all the molecular
operators with their expectation values. For a linear sweep, the evolution
equation of the atomic field operators is given by
$i\frac{d}{dt}{\hat{\Phi}}(t)=H_{4}^{(o)}(t){\hat{\Phi}}(t)$, where the non-
Hermitian matrix $H_{4}^{(o)}(t)$, in the basis
${\hat{\Phi}}(t)=\\{{\hat{a}}_{\uparrow}(t),{\hat{a}}_{\downarrow}(t),{\hat{b}}_{\uparrow}^{\dagger}(t),{\hat{b}}_{\downarrow}^{\dagger}(t)\\}$,
has the form

$H_{4}^{(o)}(t)=\begin{pmatrix}b_{1}t+E_{1}&0&g^{*}&-\gamma^{*}\\\
0&-b_{1}t+E_{1}&\gamma^{*}&g^{*}\\\ -g&-\gamma&b_{2}t+E_{2}&0\\\
\gamma&-g&0&-b_{2}t+E_{2}\end{pmatrix},$ (9)

where the coupling parameters are
$g=J_{1}\langle{\hat{\psi}}_{\uparrow\uparrow}^{\dagger}\rangle=J_{1}\langle{\hat{\psi}}_{\downarrow\downarrow}^{\dagger}\rangle$
and
$\gamma=-J_{2}\langle{\hat{\psi}}_{\uparrow\downarrow}^{\dagger}\rangle=-J_{2}\langle{\hat{\psi}}_{\downarrow\uparrow}^{\dagger}\rangle$,
where $J_{1}$ and $J_{2}$ are real valued. The parameters $b_{1}$ and $b_{2}$
are the slopes of the chemical potentials and $E_{i}$ correspond to level
separations. Note that changing the slopes, level spacings, and coupling
parameters one can explore various emergent mechanisms.

Figure 2: The pairwise squeezing of $A$ and $B$ modes is shown from (18) and
(19) and the constant shift is ignored in both directions. The parameters are
$|g|=0.5$ and $|\beta_{1}|=|\beta_{2}|=1$. The expectation values are obtained
from dimensionless operators ${\hat{X}}$.

The non-Hermitian model (9) has a form similar to multi-state LZ (MLZ) model.
The non-hermiticity leads to the emergence of complex eigenvalues near
diabatic (diagonal) level crossings. The anti-Hermitian couplings ensure that
the eigenvalues near level crossing are complex conjugate of each other. This
regime is known as PT-symmetry broken phase in the non-Hermitian literature.
Note that in the two-mode level crossing, there is one PT-broken phase
Longstaff and Graefe [2019]. The model (9), however, has four PT broken
phases. In addition, the dynamics include two paths, purple and green colors
in Fig. 1b, that can amount to a finite path interference in the system. The
interplay between the PT broken phase and the path interference strongly
affects the final occupation of the atomic modes at large times.

The solution of the evolution equation for $H_{4}^{(o)}(t)$ requires finding
the matrix $S$, which satisfies
$|\hat{\Phi}(t\rightarrow-\infty)\rangle=S|\hat{\Phi}(t\rightarrow+\infty)\rangle,$
where $S\equiv U(T,-T)_{T\rightarrow\infty}$ is a non-unitary matrix. In
Hermitian dynamics, $S$ is commonly known as the scattering matrix. The number
of atoms in each mode depends on the matrix elements $S_{ij}$. These elements,
however, cannot be obtained analytically for any arbitrary complex parameters
$g$ and $\gamma$ Malla and Raikh [2017]; Malla  _et al._ [2021]; Yan  _et al._
[2021].

In Ref. Sinitsyn, 2015 it was shown that the $CPT$ symmetry was partially
responsible for the solvability of a Hermitian Hamiltonian similar to (9) with
real coupling parameters. The three simultaneous operations in $CPT$ are
defined as follows:
(i) time-reversal (${T}$): change $t\rightarrow-t$,
${\hat{a}}\rightarrow{\hat{a}}^{\dagger}$, and
${\hat{b}}^{\dagger}\rightarrow{\hat{b}}$;
(ii) complex conjugation ($C$): change the imaginary number $i\rightarrow-i$,
${\hat{a}}\rightarrow{\hat{a}}^{\dagger}$, and
${\hat{b}}^{\dagger}\rightarrow{\hat{b}}$;
(iii) parity (${P}$): rename amplitudes
${\hat{a}}_{\uparrow}\rightarrow-{\hat{a}}_{\downarrow}$,
${\hat{a}}_{\downarrow}\rightarrow{\hat{a}}_{\uparrow}$,
${\hat{b}}_{\uparrow}\rightarrow-{\hat{b}}_{\downarrow}$, and
${\hat{b}}_{\downarrow}\rightarrow{\hat{b}}_{\uparrow}$.
It is straightforward to show that the evolution equation corresponding to the
non-Hermitian model (9) is invariant under $CPT$ transformation if the
coupling parameters are real-valued, $g=g^{*}$ and $\gamma=\gamma^{*}$.

The presence of $CPT$ symmetry imposes the following relation between the
elements of the $S$ matrix:

${\hat{S}}={\hat{S}}^{\prime}\equiv{\hat{C}}{\hat{P}}{\hat{T}}{\hat{S}}=\begin{pmatrix}S_{22}&-S_{12}&S_{42}&-S_{32}\\\
-S_{21}&S_{11}&-S_{41}&S_{31}\\\ S_{24}&-S_{14}&S_{44}&-S_{34}\\\
-S_{23}&S_{13}&-S_{43}&S_{33},\end{pmatrix}$ (10)

where $S_{ij}$ are the matrix elements of ${\hat{S}}$ and the levels {1,2,3,4}
refer to the operators
{${\hat{a}}_{\uparrow},{\hat{a}}_{\downarrow},{\hat{b}}_{\uparrow}^{\dagger},{\hat{b}}_{\downarrow}^{\dagger}$},
respectively. From (10) we obtain the following relations:

$\displaystyle S_{11}=S_{22},S_{33}=S_{44},S_{12}=S_{21}=S_{34}=S_{43}=0.$
(11)

The amplitude of the $S$ matrix elements $S_{ii}$, $S_{13}$, and $S_{24}$ can
be obtained in the independent crossing approximation Demkov and Ostrovsky
[2000]. The matrix elements $S_{14}$ and $S_{23}$ can be obtained under
solvability, when one can substitute the level spacing parameters $E_{i}$ to
be zero, and the model (9) is equivalent to models with four modes crossing at
a single point Li  _et al._ [2017].

Our model has similarity with a particular subclass of Hermitian multi-state
Landau-Zener (MLZ) problems, and the properties of scattering matrices for the
corresponding Hermitian model are listed in Ref. Sinitsyn, 2015, where the
parameters $g$ and $\gamma$ are considered real. In the Hermitian model, the
transition probability is defined as $P_{ij}=|S_{ij}|^{2}$, and the
probabilities satisfy the conservation of probability law, e.g.,
$\sum_{j}P_{ij}=1$. Here, $P_{ij}$ refers to the transition probability to
find an electron in the $i$-th state at large times $t\rightarrow\infty$
starting from $j$-th state at initial time $t\rightarrow-\infty$. The exact
solution refers to the analytical formula equation (15) in Ref. Sinitsyn, 2015
for all the transition matrix elements.

Here, the non-Hermiticity destroys the unitarity and the transition
probabilities do not add up to 1. Therefore, discussion of transition
probability becomes irrelevant. Nevertheless, $|S_{ij}|^{2}$ still describes a
physically relevant quantity and we define $n_{ij}=|S_{ij}|^{2}$, where
$n_{ij}$ refers to the number of atoms produced in the $i$-th mode due to the
dissociation process triggered by a single atom in the $j$-th mode.

To obtain matrix elements $n_{ij}$, we apply similar principles of MLZ theory
from Hermitian systems and treat the dynamics near each PT-broken phase as an
individual nLZ transition. There are four individual PT-broken phases, with
two of them controlled by $g$ while the other two are controlled by $\gamma$.
Now we can replace the survival probability
$p_{g}\equiv\exp(-\pi|g|^{2}/\beta_{2})$ with
${\tilde{n}}_{g}\equiv\exp(\pi|g|^{2}/\beta_{2})$ and the transition
probability $q_{g}=1-p_{g}$ with ${\tilde{n}}_{g}-1$. With this substitution
we find all the matrix elements of ${\hat{n}}$,

${\hat{n}}=\begin{pmatrix}\tilde{n}_{g}\tilde{n}_{\gamma}&0&\tilde{n}_{\gamma}(\tilde{n}_{g}-1)&\tilde{n}_{\gamma}-1\\\
0&\tilde{n}_{g}\tilde{n}_{\gamma}&\tilde{n}_{\gamma}-1&\tilde{n}_{\gamma}(\tilde{n}_{g}-1)\\\
\tilde{n}_{\gamma}(\tilde{n}_{g}-1)&\tilde{n}_{\gamma}-1&\tilde{n}_{g}\tilde{n}_{\gamma}&0\\\
\tilde{n}_{\gamma}-1&\tilde{n}_{\gamma}(\tilde{n}_{g}-1)&0&\tilde{n}_{g}\tilde{n}_{\gamma}\end{pmatrix},$
(12)

where matrix elements in each column, with index $i$, represent the population
in all the atomic levels when the reaction is triggered by a single atom in
the $i$-th level. The matrix elements in a column satisfy the conservation
law,

$n_{jj}-\sum_{i}n_{ij}=1,$ (13)

where $j$-th mode is occupied with a single atom. This conservation law
essentially implies that the atoms are produced in pairs and the population
difference between $A$ and $B$ is an invariant. If the initial mode is
occupied by $n_{0}$ number of atoms then the number of atoms in each level is
increased by $n_{0}$ times.

The number of stimulated atoms in each mode is shown in Fig. 1c when
$A_{\uparrow}$ is initially populated with a single atom. We find that our
analytical prediction (12) perfectly agrees with the results obtained from
numerical simulation of the evolution equation. The number of atoms in
$A_{\downarrow}$ mode is zero, as predicted by the $CPT$ symmetry. This is an
artifact of destructive interference due to the complete cancellation of two
competing reaction channels, as shown in Fig. 1b. The total number of down
spins is purely produced in the $B_{\downarrow}$ mode. The ratio between the
total number of down spin atoms and the total number of up spin atoms is then
given by the formula

$\frac{n_{\downarrow}}{n_{\uparrow}}=\frac{\tilde{n}_{\gamma}-1}{2\tilde{n}_{g}\tilde{n}_{\gamma}-\tilde{n}_{\gamma}}.$
(14)

For a small value of $|\gamma|^{2}/\beta_{1}$, the expression can be reduced
to
$n_{\downarrow}/n_{\uparrow}=e^{-\frac{\pi|g|^{2}}{2|\beta_{2}|}}(|\gamma|^{2}/2|\beta_{1}|)$,
which is exponentially suppressed with parameter $|g|^{2}/|\beta_{2}|$ and
depends linearly on parameter $|\gamma|^{2}/|\beta_{1}|$. This is one of our
key result. This peculiar dependence is shown in the parametric plot in the
inset of Fig. 1c.

We also observe that the number of atoms in any mode only depends on
$|S_{ij}|^{2}$. Our model (9) remains solvable even for the complex values of
$g$ and $\gamma$ when the phase difference between them is $\theta=n\pi$. In
addition, the number of atoms in each model can be obtained from the same
matrix ${\hat{n}}$ in (12), which was obtained using the $CPT$ symmetry for
real values of $g$ and $\gamma$. This is possible since we can transform the
evolution equation corresponding to the matrix (9) with complex $g$ and
$\gamma$ to a new evolution equation corresponding to a new matrix with real
coupling parameters that satisfy $CPT$ symmetry when the phase difference
between $g$ and $\gamma$ is $n\pi$. To understand this let us start with
parameters of the form $g=|g|e^{i\phi}$ and $\gamma=|\gamma|e^{i\chi}$. Now
the evolution equation reads

$i\frac{d{\hat{\Phi}}}{dt}=\begin{pmatrix}b_{1}t+E_{1}&0&|g|e^{-i\phi}&-|\gamma|e^{-i\chi}\\\
0&-b_{1}t+E_{1}&|\gamma|e^{-i\chi}&|g|e^{-i\phi}\\\
-|g|e^{i\phi}&-|\gamma|e^{i\chi}&b_{2}t+E_{2}&0\\\
|\gamma|e^{i\chi}&-|g|e^{i\phi}&0&-b_{2}t+E_{2}\end{pmatrix}{\hat{\Phi}}.$
(15)

Next let us use the following transformation,
${\hat{\Psi}}=\\{{\hat{a}}_{\uparrow}(t),{\hat{a}}_{\downarrow}(t),{\hat{b}}_{\uparrow}^{\dagger}(t)e^{-i\phi},{\hat{b}}_{\downarrow}^{\dagger}(t)e^{-i\phi}\\}$.
Then ${\hat{\Psi}}$ satisfies the equation of motion

$i\frac{d{\hat{\Psi}}}{dt}=\begin{pmatrix}b_{1}t+E_{1}&0&|g|&-|\gamma|e^{-i(\chi-\phi)}\\\
0&-b_{1}t+E_{1}&|\gamma|e^{-i(\chi-\phi)}&|g|\\\
-|g|&-|\gamma|e^{i(\chi-\phi)}&b_{2}t+E_{2}&0\\\
|\gamma|e^{i(\chi-\phi)}&-|g|&0&-b_{2}t+E_{2}\end{pmatrix}{\hat{\Psi}}.$ (16)

Now if the phase difference $\chi-\phi=n\pi$, then the coupling parameters
become real. This equation (16) is now invariant under $CPT$. Since the number
of atoms corresponds to $\langle{\hat{b}}^{\dagger}{\hat{b}}\rangle$, the
phase $e^{-i\phi}$ does not matter and the final atomic population can be
given by (12).

Figure 3: (a) The instantaneous eigenvalues of (9) are shown as a function of
time for parameters
$\beta_{1}=-1,\beta_{2}=0.5,E_{1}=5,E_{2}=1,g=0.5,\gamma=i$. The black curves
are real values and the red curves represent imaginary values. (b) The
stimulated population of atomic modes from initially prepared single atom in
$A_{\uparrow}$ is shown as a function of $\theta$. The colors blue, green,
red, and black correspond to atomic modes $A_{\uparrow}$, $B_{\uparrow}$,
$A_{\downarrow}$, and $B_{\downarrow}$, respectively. (c) The Schematic setup
to realize the dynamics of our model (9). The red and green circles correspond
to optical waveguides and the nearest coupling between waveguides are
imaginary.

###  III.2 Emergence of multi-mode squeezed states due to spontaneous emission

The spontaneous dissociation process in the four-mode model (9) is quantum
correlated. The number of atoms in each mode is equal due to the $CPT$
symmetry and can be obtained from (12),

$n_{\alpha,\sigma}={\tilde{n}}_{g}{\tilde{n}}_{\gamma}-1,$ (17)

where $\alpha$ refers to $A$ or $B$ and $\sigma$ refers to the spin
orientation. Another key quantum phenomenon emerging in the spontaneous
dissociation process is that atoms are produced in a multi-mode squeezed
state. The solvability of our model allows us to evaluate the amount of
squeezing, and we calculate pair-wise squeezing of atomic modes. First, we
observe that the modes $A_{\uparrow}$ and $A_{\downarrow}$ are not quantum
correlated and are not in a two-mode squeezed state. The maximum uncertainty
along both $X_{+}$ and $X_{-}$ is equal, where $X$ is the position operator of
the two-mode squeezed state. The same applies to modes $B_{\uparrow}$ and
$B_{\downarrow}$.

Therefore, the atoms produced in the spontaneous process are always in a
three-mode squeezed state Rojas-Rojas  _et al._ [2019]. To observe squeezing
between $A_{\uparrow}$ and $B_{\uparrow}$ ($B_{\downarrow}$), we evaluate
$\langle X_{\pm}^{2}\rangle$ (see appendix B) and find the maximum values
along each direction $X_{\pm}$ is given by

$\displaystyle\langle
X_{\pm}^{2}\rangle_{A_{\uparrow}B_{\uparrow}}=c_{\gamma}+\frac{1}{2}|(\tilde{n}_{g}\tilde{n}_{\gamma})^{1/2}\pm((\tilde{n}_{g}-1)\tilde{n}_{\gamma})^{1/2}|^{2},\leavevmode\nobreak\
$ (18) $\displaystyle\langle
X_{\pm}^{2}\rangle_{A_{\uparrow}B_{\downarrow}}=c_{g\gamma}+\frac{1}{2}|(\tilde{n}_{g}\tilde{n}_{\gamma})^{1/2}\pm(\tilde{n}_{\gamma}-1)^{1/2}|^{2},\leavevmode\nobreak\
$ (19)

where $c_{\gamma}=\frac{{\tilde{n}}_{\gamma}-1}{2}$ and
$c_{g\gamma}=\frac{{\tilde{n}}_{\gamma}{\tilde{n}}_{g}-1}{2}$ are constant
shifts. In Fig. 2 we plot $\langle X_{\pm}^{2}\rangle$ for different values of
$|\gamma|$. The atoms produced in $A_{\uparrow}$ and $B_{\uparrow}$ mode are
always squeezed along the $X_{-}$ direction, and this phenomenon has the same
origin as squeezing in a simple two-mode atomic dissociation. The atoms in
$A_{\uparrow}$ and $B_{\downarrow}$ modes are squeezed along $X_{+}$ for small
values of $|\gamma|$. This squeezing is suppressed with an increase in
$|\gamma|$ and eventually changes to $X_{-}$ direction. Note that the
direction of the squeezing is rotated by the phase $\phi$, since $\langle
X^{2}\rangle$ includes the correlator $\langle{\hat{a}}{\hat{b}}\rangle$. From
the symmetry of couplings, we find that the squeezing for up-spin $A$ and $B$
modes and down-spin $A$ and $B$ modes is equal, while the squeezing between
up-spin $A$ and down-spin $B$ modes is the same as down-spin $A$ and up-spin
$B$ modes.

###  III.3 Phase difference between $g$ and $\gamma$ and broken solvability.–

Now we discuss the scenario when there is a finite phase difference between
$g$ and $\gamma$. This finite phase difference, $\theta$, is responsible for
violation of CPT symmetry, and gaps emerge near the exact crossing of levels
with opposite slopes, see Fig. 3a. This violation of exact crossing leads to
asymmetric phase accumulations along the two interfering paths, leading to a
finite interference in the system. Due to the finite interference, initial
atoms in $A_{\uparrow}$ mode can stimulate the production of atoms in
$A_{\downarrow}$ mode. However, an exact analytical solution is impossible,
and we turn to a numerical approach.

In Fig. 3b, we show the $\theta$ dependence of the number of atoms in each
mode stimulated by a single atom in $A_{\uparrow}$ mode. For $\theta=0$ or
$\pi$ we recover the solvability condition and the $CPT$ symmetry of the
evolution equation, which indeed agrees with formula (14), and the down-spin
atoms are exponentially suppressed. The production of down-spin atoms
increases with increasing $\theta$, which corresponds to widening of the gaps
and stronger constructive interference, see Fig. 3b. Note, the imaginary parts
of the eigenvalues, however, remain unaffected. At $\theta=(2m+1)\pi/2$, the
interference becomes constructive and stimulates a maximum down spin.

Concerning spontaneous dissociation, the total number atoms in each mode is no
longer equal, however they still satisfy the condition
$n_{A_{\uparrow}}+n_{A_{\downarrow}}=n_{B_{\uparrow}}+n_{B_{\downarrow}}$. The
condition that atoms are produced in pairs ensures the quantum correlation
between modes, and here the atoms will be produced in four-mode squeezing due
to an effective coupling between $A_{\uparrow}$ and $A_{\downarrow}$.

##  IV Discussion

(i.) We have solved a mechanism of dissociation of molecular condensates when
there are competing chemical reactions due to multiple atomic modes. Our model
can distinguish between the types of atoms produced in comparison to prior
studies. The validity of our results relies on the fact that the molecular
operators can be replaced with their expectation values. This assumption only
describes a strong nonadiabatic conversion mechanism where only a small
fraction of molecules can be converted to atoms. In our previous work Ref.
Malla  _et al._ , 2022 we analytically solved a model describing two-atomic
modes where the exact solution allowed us to investigate both the nonadiabatic
as well as the adiabatic regime. This work was motivated by the experiments in
Ref. Zhang  _et al._ , 2021 where $Cs$ atoms were converted into $Cs_{2}$
molecules.

We only investigate the molecular dissociation process in the mean-field
approximation and assume that molecules and atoms do not interact with each
other. The experimental feasibility of such a system will depend on the
preparation of molecular condensates in four modes. To prepare a simple multi-
component condensate, one approach could be to prepare the molecules in
different quantum states prior to condensationZhang  _et al._ [2017]. This can
be further simplified by assuming diatomic molecules consist of the same type
of atoms, which will require us to prepare only 3-component condensates. When
the molecules dissociate due to the time-varying chemical potentials, they
will form atomic condensates. The number of atoms prepared is small compared
to molecules are small and therefore the interaction between atoms can be
neglected. A complete understanding requires understanding the stability of
such a complex system which is beyond the scope of this article.

(ii) Photonic realization of non-Hermitian multi-state Landau-Zener models\-
The similarity between the paraxial Helmoltz equation and the time-dependent
Schrödinger equation allows to implement various quantum mechanical models in
photonic waveguides. To test various non-Hermitian dynamics, photonic
waveguides provide a desirable platform due to the ability to control and
manipulate the gain and loss of the media. The light propagation in a network
of $N$ arbitrarily coupled waveguides satisfies coupled-mode equations in the
form of the paraxial Helmholtz equation. When the modes are orthogonal to each
other, the paraxial Helmholtz equation resembles the time-dependent
Schrödinger equation $i\partial\psi(z)/\partial z=H(z)\psi(z),$ where the
“Hamiltonian” $H(z)$ is position dependent. To realize nMLZ models the
refractive index in each waveguide can be made to be linearly dependent on the
position ($z$) in order to design a particular setup of diabatic level
crossings, similar to the one shown in Fig. 3c. In order to introduce the nMLZ
dynamics, the coupling between the waveguides must be anti-Hermitian. The
simplest way to satisfy anti-hermiticity is to consider the coupling to be
imaginary. Recently, an imaginary coupling was realized in experiment Ke  _et
al._ [2018], where the authors achieved the imaginary coupling between two
waveguides by introducing an anchila waveguide in between them. This will
allow for any possible arrangement to test a general nMLZ model. In Fig. 3c,
we show a design to test our model (9). This design allows for imaginary
coupling between two nearest neighbors and forbids coupling between diagonal
waveguides. This design should be able to test the model (9) for any set of
parameters, and the interfering paths are shown by red and blue arrows. For
our solvable model, one should expect the following observations. First,
staring light at waveguide 1, one should observe zero intensity in the
waveguide 2 at large distances. Second, in the large coupling limit one should
observe equally strong intense light in 1 and 3. The intensity of light in the
fourth waveguide will not be zero, however in comparison with 1 and 3, the
intensity will be exponentially suppressed.

## Acknowledgements

R.K.M thanks Nikolai Sinitsyn, Avadh Saxena, Julia Cen, Mikhail Raikh, Wilton
Kort-Kamp, and Saikat Banerjee for useful discussions. This work was supported
by the U.S. Department of Energy, Office of Science, Basic Energy Sciences,
Materials Sciences and Engineering Division, Condensed Matter Theory Program,
and the LANL Center for Nonlinear Studies.

##  Appendix A Non-Hermitian Landau-Zener model

The non-Hermitian Landau-Zener (nLZ) model can be described by the matrix

${\hat{H}}_{nLZ}=\begin{pmatrix}\beta t&&g\\\ -g^{*}&&-\beta t\end{pmatrix},$
(20)

where $\pm\beta$ are the two slopes corresponding to two levels and $g$ is the
level coupling, which in general is complex. Unlike the Hermitian model the
level couplings $g$ and $-g^{*}$ are negative complex conjugate of each other.
The eigenvalues of ${\hat{H}}_{nLZ}$ are given by

$E_{nLZ}(t)=\pm\sqrt{\beta^{2}t^{2}-|g|^{2}}.$ (21)

In standard LZ model the eigenvalues are given by
$E_{LZ}=\pm\sqrt{\beta^{2}t^{2}+|g|^{2}}$.

Figure 4: The (a) eigen-energies of a standard LZ model and (b) the
eigenvalues of a matrix describing the nLZ model ($E_{nLZ}(t)$) are
schematically shown as a function of time. The blue corresponds to real part
of eigenvalues while the red corresponds to the imaginary part.

The eigenvalues $E_{nLZ}(t)$ are shown as a function of time in Fig. 4b
together with the eigen-energies for a standard LZ model with parameter
$\beta$ and $g$.

The solution of the evolution equation for matrix (20) has the form of a
$2\times 1$ column vector,

$|\phi(t)\rangle=\begin{pmatrix}a(t)\\\ b(t)\end{pmatrix},$

where $a(t)$ satisfies a second order differential equation

$\ddot{a}(t)+(\beta^{2}t^{2}-|g|^{2}+i\beta)a(t)=0,$ (22)

whose solutions are given by parabolic cylinder functions par .

The solution can now be expressed as follows Malla  _et al._ [2017]

$|\phi(t)\rangle=\phi_{1}\begin{pmatrix}D_{\nu}(z)\\\
-i\sqrt{\nu}D_{\nu-1}(z)\end{pmatrix}+\phi_{2}\begin{pmatrix}D_{\nu}(-z)\\\
-i\sqrt{\nu}D_{\nu-1}(-z)\end{pmatrix},$ (23)

where $D_{\nu}(z)$ is the parabolic cylinder function, with
$\nu=i|g|^{2}/2\beta$ and $z=\sqrt{2\beta}e^{i\pi/4}t$. Since we are only
interested in the dynamics at large times $t\rightarrow\pm\infty$, the
asymptotic behavior of parabolic cylinder functions is good enough for our
analysis. The critical difference between the Hermitian and the non-Hermitian
dynamics comes from the phase of $\nu$, which is $-\pi/2$ for the Hermitian
case and $\pi/2$ for the non-Hermitian case. This crucial difference leads to
the following distinction between LZ and nLZ dynamics. Assuming the system
starts with the upper state, $|a(t\rightarrow-\infty)|^{2}=1$, the asymptotic
solution of $a(t)$ at large positive times is given by
$|a(t\rightarrow\infty)|^{2}=e^{-\pi|g|^{2}/\beta}$ for the Hermitian model
and $|a(t\rightarrow\infty)|^{2}=e^{\pi|g|^{2}/\beta}$ for the non-Hermitian
model. Similarly, the solution $b(t)$ at large positive times is given by
$|b(t\rightarrow\infty)|^{2}=1-e^{-\pi|g|^{2}/\beta}$ for the Hermitian model
and $|b(t\rightarrow\infty)|^{2}=e^{\pi|g|^{2}/\beta}-1$ for the non-Hermitian
model.

This exponential dependence in the non-Hermitian model appears due to the
anti-Hermitian complex couplings. The conservation law in the Hermitian model
describes the common probability law,

$|a(t)|^{2}+|b(t)|^{2}=1.$

The probability of survival is given by $|a(t\rightarrow\infty)|^{2}$ while
the probability of adiabatic transition is given by
$|b(t\rightarrow\infty)|^{2}$. In the non-Hermitian model, however, there are
no such terms, and the exponentially growing nature makes it difficult to call
it probability. Nevertheless, in our model, $|a(t\rightarrow\infty)|^{2}$ and
$|b(t\rightarrow\infty)|^{2}$ correspond to the number of atoms produced in
the respective atomic modes. The conservation law here results from the
asymptotic expansion of parabolic cylinder functions,

$|a(t)|^{2}-|b(t)|^{2}=1,$

and agrees with our dissociation mechanism in which atoms are produced in
pairs.

##  Appendix B Pairwise squeezing of atomic modes

The spontaneous dissociation process is quantum correlated and the atoms are
produced in multi-mode quantum squeezed states. To quantify squeezing, we
define position operators for the two modes corresponding to different types
atoms $A$ and $B$. The position operator reads

${\hat{X}}_{A_{\sigma}B_{\sigma^{\prime}}}^{(\phi)}=\frac{1}{2}\left[e^{-i\phi}({\hat{a}}_{\sigma}+{\hat{b}}_{\sigma^{\prime}})+e^{i\phi}({\hat{a}}_{\sigma}^{\dagger}+{\hat{b}}_{\sigma^{\prime}}^{\dagger})\right],$
(24)

where $\sigma$, $\sigma^{\prime}$ refer to the spins while $\phi$ is the angle
that defines the measured quadrature. It is straightforward to see that the
expectation value of operator $X$ is zero, while the expectation value of
$X^{2}$ can be non-zero. We expand the expression for the $X^{2}$ operator
from (24)

${\hat{X}^{(\theta)^{2}}}_{\pm,A_{\sigma}B_{\sigma^{\prime}}}=\frac{1}{4}\Big{[}e^{-2i\phi}({\hat{a}}_{\sigma}^{2}+{\hat{b}}_{\sigma^{\prime}}^{2}+{\hat{a}}_{\sigma}{\hat{b}}_{\sigma^{\prime}}+{\hat{b}}_{\sigma^{\prime}}{\hat{a}}_{\sigma})+\\\
{\hat{a}^{\dagger}}_{\sigma}{\hat{a}}_{\sigma}+{\hat{a}^{\dagger}}_{\sigma}{\hat{b}}_{\sigma^{\prime}}+{\hat{b}^{\dagger}}_{\sigma^{\prime}}{\hat{a}}_{\sigma}+{\hat{b}^{\dagger}}_{\sigma^{\prime}}{\hat{b}}_{\sigma^{\prime}}\Big{]}+h.c.$
(25)

At large times $t\rightarrow\infty$, the operators satisfy

$\displaystyle{\hat{a}}_{\uparrow}(t)=S_{11}{\hat{a}}_{\uparrow}(t_{0})+S_{12}{\hat{a}}_{\downarrow}(t_{0})+S_{13}{\hat{b}}_{\uparrow}^{\dagger}(t_{0})+S_{14}{\hat{b}}_{\downarrow}^{\dagger}(t_{0}),$
$\displaystyle{\hat{a}}_{\downarrow}(t)=S_{21}{\hat{a}}_{\uparrow}(t_{0})+S_{22}{\hat{a}}_{\downarrow}(t_{0})+S_{23}{\hat{b}}_{\uparrow}^{\dagger}(t_{0})+S_{24}{\hat{b}}_{\downarrow}^{\dagger}(t_{0}),$
$\displaystyle{\hat{b}}_{\uparrow}^{\dagger}(t)=S_{31}{\hat{a}}_{\uparrow}(t_{0})+S_{32}{\hat{a}}_{\downarrow}(t_{0})+S_{33}{\hat{b}}_{\uparrow}^{\dagger}(t_{0})+S_{34}{\hat{b}}_{\downarrow}^{\dagger}(t_{0}),$
$\displaystyle{\hat{b}}_{\downarrow}^{\dagger}(t)=S_{41}{\hat{a}}_{\uparrow}(t_{0})+S_{42}{\hat{a}}_{\downarrow}(t_{0})+S_{43}{\hat{b}}_{\uparrow}^{\dagger}(t_{0})+S_{44}{\hat{b}}_{\downarrow}^{\dagger}(t_{0}).$
(26)

Here, $t_{0}$ corresponds to initial time. The non-zero contribution to
$\langle X^{2}\rangle$ comes from terms proportional to
$\langle{\hat{a}}(t_{0}){\hat{a}}^{\dagger}(t_{0})\rangle$ and
$\langle{\hat{b}}(t_{0}){\hat{b}}^{\dagger}(t_{0})\rangle$. The term
$\langle{\hat{a}}(t_{0}){\hat{a}}^{\dagger}(t_{0})\rangle=1+\langle{\hat{a}}^{\dagger}(t_{0}){\hat{a}}(t_{0})\rangle$
and the second term vanishes since the initial atomic population is zero.

Now we can express the expectation value of $X^{2}$ with only non-zero terms

$\langle{\hat{X}^{(\theta)^{2}}}_{A_{\sigma}B_{\sigma^{\prime}}}\rangle=\frac{1}{4}\left[e^{-2i\phi}({\hat{a}}_{\sigma}{\hat{b}}_{\sigma^{\prime}}+{\hat{b}}_{\sigma^{\prime}}{\hat{a}}_{\sigma})+{\hat{a}^{\dagger}}_{\sigma}{\hat{a}}_{\sigma}+{\hat{b}^{\dagger}}_{\sigma^{\prime}}{\hat{b}}_{\sigma^{\prime}}\right]+h.c.$
(27)

Let us first consider the expectation value of $X^{2}$ for $A_{\uparrow}$ and
$B_{\uparrow}$ modes. Substituting equation (26) in equation (25) we find

$\langle{\hat{X}^{(\theta)^{2}}}_{\pm,A_{\uparrow}B_{\uparrow}}\rangle=\frac{1}{2}\left[Re[S_{11}S_{13}e^{2i\phi}]+|S_{11}|^{2}+|S_{13}|^{2}+|S_{14}|^{2}\right],$
(28)

where

$-2|S_{11}||S_{13}|\leq Re[S_{11}S_{13}e^{2i\phi}]\leq 2|S_{11}||S_{13}|.$

We can obtain the maximum and the minimum value of
$\langle{\hat{X}^{(\theta)^{2}}}_{\pm,A_{\uparrow}B_{\uparrow}}\rangle$, which
are given by

$\langle{\hat{X}^{(\theta)^{2}}}_{\pm,A_{\uparrow}B_{\uparrow}}\rangle=\frac{(\tilde{n}_{\gamma}-1)}{2}+\frac{1}{2}|(\tilde{n}_{g}\tilde{n}_{\gamma})^{1/2}\pm((\tilde{n}_{g}-1)\tilde{n}_{\gamma})^{1/2}|^{2},$
(29)

and we obtain our equation (20) from the main text. For $\gamma=0$ we recover
the two-mode squeezing in the spin-independent atomic dissociation in Ref.
Yurovsky  _et al._ [2002].

Similarly we can obtain equation (21) in our main text by switching $S_{13}$
and $S_{14}$ in equation (28).

## References

  * Bohn  _et al._ [2017a] J. L. Bohn, A. M. Rey, and J. Ye, Science 357, 1002 (2017a), https://www.science.org/doi/pdf/10.1126/science.aam6299 .
  * Carr  _et al._ [2009] L. D. Carr, D. DeMille, R. V. Krems, and J. Ye, New Journal of Physics 11, 055049 (2009).
  * QuÃ©mÃ©ner and Julienne [2012] G. QuÃ©mÃ©ner and P. S. Julienne, Chemical Reviews 112, 4949 (2012), pMID: 22921011, https://doi.org/10.1021/cr300092g .
  * Zhang  _et al._ [2021] Z. Zhang, L. Chen, K.-X. Yao, and C. Chin, Nature 592, 708 (2021).
  * Hu  _et al._ [2019] J. Hu, L. Feng, Z. Zhang, and C. Chin, Nature Physics 15, 785 (2019).
  * Clark  _et al._ [2017] L. W. Clark, A. Gaj, L. Feng, and C. Chin, Nature 551, 356 (2017).
  * Yurovsky  _et al._ [2002] V. A. Yurovsky, A. Ben-Reuven, and P. S. Julienne, Phys. Rev. A 65, 043607 (2002).
  * Altland  _et al._ [2009] A. Altland, V. Gurarie, T. Kriecherbauer, and A. Polkovnikov, Phys. Rev. A 79, 042703 (2009).
  * Sun and Sinitsyn [2016] C. Sun and N. A. Sinitsyn, Phys. Rev. A 94, 033808 (2016).
  * Itin and Törmä [2009] A. P. Itin and P. Törmä, Phys. Rev. A 79, 055602 (2009).
  * Malla  _et al._ [2021] R. K. Malla, V. Y. Chernyak, and N. A. Sinitsyn, Phys. Rev. B 103, 144301 (2021).
  * Herbig  _et al._ [2003] J. Herbig, T. Kraemer, M. Mark, T. Weber, C. Chin, H.-C. Nägerl, and R. Grimm, Science 301, 1510 (2003).
  * Xu  _et al._ [2003] K. Xu, T. Mukaiyama, J. R. Abo-Shaeer, J. K. Chin, D. E. Miller, and W. Ketterle, Phys. Rev. Lett. 91, 210402 (2003).
  * Bohn  _et al._ [2017b] J. L. Bohn, A. M. Rey, and J. Ye, Science 357, 1002 (2017b).
  * Matyjaśkiewicz  _et al._ [2008] S. Matyjaśkiewicz, M. H. Szymańska, and K. Góral, Phys. Rev. Lett. 101, 150410 (2008).
  * Duda  _et al._ [2021] M. Duda, X.-Y. Chen, A. Schindewolf, R. Bause, J. Milczewski, I. Schmidt, R. Bloch, and X.-Y. Luo, Preprit arXiv:2111.04301 (2021).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Reviews of Modern Physics 82, 1225 (2010).
  * Malla  _et al._ [2022] R. K. Malla, V. Y. Chernyak, C. Sun, and N. A. Sinitsyn, Phys. Rev. Lett. 129, 033201 (2022).
  * Green and Porto [2020] D. Green and R. A. Porto, Phys. Rev. Lett. 124, 251302 (2020).
  * Gessner  _et al._ [2018] M. Gessner, L. Pezzè, and A. Smerzi, Phys. Rev. Lett. 121, 130503 (2018).
  * Coelho  _et al._ [2009] A. S. Coelho, F. A. S. Barbosa, K. N. Cassemiro, A. S. Villar, M. Martinelli, and P. Nussenzveig, Science 326, 823 (2009), https://www.science.org/doi/pdf/10.1126/science.1178683 .
  * Simon  _et al._ [2007] C. Simon, H. de Riedmatten, M. Afzelius, N. Sangouard, H. Zbinden, and N. Gisin, Phys. Rev. Lett. 98, 190503 (2007).
  * Duan  _et al._ [2001] L. M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller, Nature 414, 413 (2001).
  * Kayali and Sinitsyn [2003] M. A. Kayali and N. A. Sinitsyn, Phys. Rev. A 67, 045603 (2003).
  * Lin  _et al._ [2011] Y. J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Longstaff and Graefe [2019] B. Longstaff and E.-M. Graefe, Phys. Rev. A 100, 052119 (2019).
  * Malla and Raikh [2017] R. K. Malla and M. E. Raikh, Phys. Rev. B 96, 115437 (2017).
  * Yan  _et al._ [2021] B. Yan, V. Y. Chernyak, W. H. Zurek, and N. A. Sinitsyn, Phys. Rev. Lett. 126, 070602 (2021).
  * Sinitsyn [2015] N. A. Sinitsyn, Phys. Rev. B 92, 205431 (2015).
  * Demkov and Ostrovsky [2000] Y. N. Demkov and V. N. Ostrovsky, Phys. Rev. A 61, 032705 (2000).
  * Li  _et al._ [2017] F. Li, C. Sun, V. Y. Chernyak, and N. A. Sinitsyn, Phys. Rev. A 96, 022107 (2017).
  * Rojas-Rojas  _et al._ [2019] S. Rojas-Rojas, E. Barriga, C. Muñoz, P. Solano, and C. Hermann-Avigliano, Phys. Rev. A 100, 023841 (2019).
  * Zhang  _et al._ [2017] Y.-C. Zhang, S. Ding, and S. Zhang, Phys. Rev. A 95, 041603 (2017).
  * Ke  _et al._ [2018] S. Ke, D. Zhao, Q. Liu, S. Wu, B. Wang, and P. Lu, Journal of Lightwave Technology 36, 2510 (2018).
  * [35] “Higher transcendental functions, edited by a. erdelyi, vol. 2 (mcgraw-hill, new york, 1953).” .
  * Malla  _et al._ [2017] R. K. Malla, E. G. Mishchenko, and M. E. Raikh, Phys. Rev. B 96, 075419 (2017).


# Critical fluorescence of a transmon at the Schmid transition

M. Houzet  Univ. Grenoble Alpes, CEA, IRIG, Pheliqs, F-38000 Grenoble, France
L. I. Glazman  Departments of Physics and Applied Physics, Yale University,
New Haven, CT 06520, USA

###### Abstract

We investigate inelastic microwave photon scattering by a transmon qubit
embedded in a high-impedance circuit. The transmon undergoes a charge-
localization (Schmid) transition upon the impedance reaching the critical
value. Due to the unique transmon level structure, the fluorescence spectrum
carries a signature of the transition point. At higher circuit impedance,
quasielastic photon scattering may account for the main part of the inelastic
scattering cross-section; we find its dependence on the qubit and circuit
parameters.

Introduction.– A quantum-mechanical degree of freedom can be severely affected
by its coupling to a dissipative environment. In a pioneering work Schmid1983
, Schmid predicted that a superconducting circuit as elementary as a Josephson
junction is insulating when it is ohmically shunted by a resistance larger
than the resistance quantum, $R_{Q}=\pi\hbar/2e^{2}$. This result, which was
further studied in Refs. Bulgadaev1984 ; Guinea1985 ; Fisher1985 , reflects a
charge localization transition, which is associated with the breaking of the
ground state degeneracy. Remarkably the prediction holds at any ratio between
the Josephson and charging energies of the junction, $\tilde{E}_{J}$ and
$\tilde{E}_{C}=e^{2}/2\tilde{C}$, respectively, where $\tilde{C}$ is the
junction capacitance Schon1990 . So far evidence for the charge localization
transition by dc Yagi1997 ; Pentilla2001 ; Watanabe2003  or low-frequency
Murani2019  measurements remains elusive.

As a quantum many-body effect, the Schmid transition should not only affect
the ground state, but the excited states as well. Here we find a spectroscopic
signature of the transition in the fluorescence spectrum footnote0  of a
weakly nonlinear Josephson junction with $\tilde{E}_{C}\ll\tilde{E}_{J}$,
a.k.a. a transmon qubit Koch2007 ; Schreier2008 , coupled to a Josephson-
junction chain. The chain realizes a transmission line with an adjustable
impedance Masluk2012 ; Bell2012 , in which plasmon waves or, equivalently,
microwave photons propagate freely. This setup attracted recent experimental
interest Puertas2019 ; Kuzmin2019 ; Leger2019  as a way to emulate quantum
impurity problems with superconducting quantum circuits Garcia-Ripoll2008 ;
LeHur2012 ; Goldstein2013 ; Snyman2015 ; Gheeraert2018 ; Yamamoto2019 . Our
theory predicts a characteristic dependence of the inelastic scattering cross-
sections on the parameters of the setup.

Model.– Our consideration starts with the superconducting circuit of Fig. 1.
The Josephson-junction chain to which the transmon is coupled is characterized
by the Josephson energy $E_{J}$, and charging energies $E_{C}=e^{2}/2C$ and
$E_{g}=e^{2}/2C_{g}$, where $C$ and $C_{g}$ are the chain’s junction and
ground capacitances. Under the conditions

$E_{C}E_{J}\gg{\tilde{E}}_{C}{\tilde{E}}_{J}\gg{\tilde{E}}^{2}_{C},$ (1)

the classical analysis, in which the Josephson junctions are approximated as
linear inductances, yields a narrow transmon resonance that lies deep in the
linear part of the waves’ dispersion. The transmon resonance frequency and
half-width are, respectively,
$\omega_{0}=\sqrt{8{\tilde{E}}_{C}{\tilde{E}}_{J}}/\hbar$ and
$\Gamma=1/2Z\tilde{C}$; the waves’ dispersion is
$\omega_{p}=vp/\sqrt{(vp/\omega_{\mathrm{B}})^{2}+1}$. Here $Z=R_{Q}/2K$ with
$K=\pi\sqrt{E_{J}/8E_{g}}$ is the chain’s impedance at low frequency,
$\omega_{\mathrm{B}}=\sqrt{8{E}_{C}{E}_{J}}/\hbar$ is the photon bandwidth,
$v=a\sqrt{8E_{J}E_{g}}/\hbar$ is the photon velocity and $a$ is the chain’s
unit cell length. Noting that $\Gamma=(4/\pi\hbar)K{\tilde{E}}_{C}$ and
assuming $K\lesssim 1$ in the relevant range of parameters, we find from Eq.
(1) that, indeed, $\Gamma\ll\omega_{0}\ll\omega_{\mathrm{B}}$. Within the
classical analysis, photon scattering is purely elastic.

Figure 1:  We study microwave photon scattering off a transmon coupled to a
Josephson-junction chain. The transmon has Josephson energy $\tilde{E}_{J}$
and capacitance $\tilde{C}$. The junctions in the chain have Josephson energy
$E_{J}$ and capacitance $C$, each superconducting island has ground
capacitance $C_{g}$.

To analyze inelastic scattering, we first notice that $K=1/2$ is the critical
value for the Schmid transition Schmid1983 . Below that value, the environment
induces charge localization. The elementary processes responsible for it are
phase slips Hekking1997 . In an isolated transmon, the phase slip amplitude
for the first excited band is Koch2007

$\lambda_{1}=\frac{64}{\sqrt{\pi}}\frac{{\tilde{E}}_{C}}{\hbar}\left(\frac{2\tilde{E}_{J}}{\tilde{E}_{C}}\right)^{5/4}e^{-\sqrt{8{\tilde{E}_{J}}/{\tilde{E}_{C}}}}$
(2)

(parametrically larger than the phase slip amplitude for the ground state,
$\lambda_{1}/\lambda_{0}=-8\sqrt{2{\tilde{E}}_{J}/{\tilde{E}}_{C}}$).
Therefore, $\Gamma$ may become parametrically smaller than $\lambda_{1}$ only
deep in the localized regime, $K\ll 1$. Barring that, we treat phase slips
perturbatively. In most of the discussion below we assume
$\lambda_{1}\ll\Gamma$. At the same time, the proliferation of phase slips in
the chain’s junctions, with an amplitude $\lambda_{\rm chain}\propto
e^{-\sqrt{8E_{J}/E_{C}}}$, is known to drive a superfluid-to-insulating
transition below the critical value $K=2$ Bradley1984 ; Korshunov1989 ;
Basko2020 . That physics can be disregarded either if the bandwidth $\Gamma$
that limits frequency exchange in a quasielastic process (see below) is larger
than the insulating gap, $\Gamma\gg\lambda_{\mathrm{chain}}$, or if the chain
is short enough to evade the thermodynamic limit, $L<v/\lambda_{\rm chain}$
(but still long enough to ignore effects related to a finite level spacing,
$\Delta=\pi v/L\to 0$).

Ignoring the chain’s nonlinearity we use Hamiltonian

$H=4\tilde{E}_{C}(\hat{N}-\hat{n})^{2}+\tilde{E}_{J}[1-\cos{\hat{\varphi}}]+\sum_{p}\hbar\omega_{p}a^{\dagger}_{p}a_{p}$
(3)

to describe the setup. Here $\hat{\varphi}$ and $\hat{N}$ are the transmon’s
conjugate phase and charge operators, $a_{p}$ and $a^{\dagger}_{p}$ are the
annihilation and creation operators of a linearly-dispersing photon of wave
vector $p$, and

$\hat{n}=\frac{1}{\pi}\sum_{p>0}f_{p}(a_{p}+a_{p}^{\dagger})\quad\mathrm{with}\quad
f_{p}=\sqrt{{K\Delta}/{\omega_{p}}}$ (4)

is the charge displacement operator at the qubit. The sum in Eq. (4) includes
only the dynamical variables of the chain. The static mode ($p=0$) compensates
for an eventual charge offset, which is effectively attenuated by the total
capacitance of the array, $\sim C_{g}\\!\cdot\\!(L/a)$ vlad .

The low-energy properties of Eq. (3) are described by the boundary sine-Gordon
model, $H_{\rm
sG}=\sum_{p}\hbar\omega_{p}a^{\dagger}_{p}a_{p}-\lambda_{0}\cos(2\pi\hat{n})$
Gogolin-book . Keeping the ratio $\tilde{E}_{C}/\tilde{E}_{J}$ small but
finite generalizes it to a new quantum impurity problem, in which the transmon
resonance brings in a nontrivial structure of the high-frequency spectrum, as
we discuss now.

Phase-slip vs quartic anharmonicities.– We first show, by perturbative
analysis in coupling $f_{p}$, that phase slips play the dominant role in
quasielastic (or soft) photon scattering, i.e., inelastic scattering with a
small energy transfer between the incoming and one of the outgoing photons.
For this we write the Hamiltonian (3) in the transmon basis of discrete
eigenstates (we set $\hbar=1$ hereinafter),

$\displaystyle H$ $\displaystyle=$ $\displaystyle H_{0}+V,\quad
H_{0}=\sum_{s}\varepsilon_{s}|s\rangle\langle
s|+\sum_{p}\omega_{p}a^{\dagger}_{p}a_{p}+4E_{C}{\hat{n}}^{2},$ $\displaystyle
V$ $\displaystyle=$
$\displaystyle-\hat{n}\sum_{ss^{\prime}}W_{ss^{\prime}}|s\rangle\langle
s^{\prime}|\quad{\rm with}\quad W_{ss^{\prime}}=8E_{C}\langle
s|\hat{N}|s^{\prime}\rangle.$ (5)

Here we ordered the transmon eigenenergies,
$\varepsilon_{s}>\varepsilon_{s^{\prime}}$ if $s>s^{\prime}\geq 0$, and set
$\varepsilon_{0}=0$. The partial inelastic cross-section for a photon with
frequency $\omega$ to be converted into three photons with frequencies
$\omega_{1},\omega_{2},\omega_{3}$, such that
$\omega=\omega_{1}+\omega_{2}+\omega_{3}$, is obtained with Fermi’s Golden
Rule,

$\gamma(\omega_{1},\omega_{2},\omega_{3}|\omega)=\frac{2\pi^{2}}{3!}\frac{\left|{\cal
A}_{p;p_{1}p_{2}p_{3}}\right|^{2}}{\Delta^{4}}$ (6)

($3!$ accounts for permutations of momenta that describe the same final state)
with $\omega=\omega_{p}$, $\omega_{i}=\omega_{p_{i}}$, and a matrix element
obtained perturbatively in ${f_{p}}$,

$\displaystyle{\cal A}_{p;p_{1}p_{2}p_{3}}$ $\displaystyle=$
$\displaystyle\langle
0|a_{p_{1}}a_{p_{2}}a_{p_{3}}V\left(\frac{1}{\omega_{{p}}-H_{0}}V\right)^{3}a^{\dagger}_{p}|0\rangle$
(7) $\displaystyle=$
$\displaystyle-\frac{K^{2}\Delta^{2}}{\pi^{4}\sqrt{\omega\omega_{1}\omega_{2}\omega_{3}}}\sum_{srt}W_{0s}W_{sr}W_{rt}W_{t0}$
$\displaystyle\times\left\\{\frac{1}{(\varepsilon_{s}+\omega_{1}+\omega_{2}-\omega)(\varepsilon_{r}+\omega_{1}-\omega)(\varepsilon_{t}-\omega)}\right.$
$\displaystyle+\left.\frac{1}{(\varepsilon_{s}+\omega_{1}+\omega_{2}-\omega)(\varepsilon_{r}+\omega_{1}-\omega)(\varepsilon_{t}+\omega_{1})}\right.$
$\displaystyle+\left.\frac{1}{(\varepsilon_{s}+\omega_{1}+\omega_{2}-\omega)(\varepsilon_{r}+\omega_{1}+\omega_{2})(\varepsilon_{t}+\omega_{1})}\right.$
$\displaystyle+\left.\frac{1}{(\varepsilon_{s}+\omega_{1}+\omega_{2}+\omega_{3})(\varepsilon_{r}+\omega_{1}+\omega_{2})(\varepsilon_{t}+\omega_{1})}\right.$
$\displaystyle+\left.\text{permutations
of}\,\,\omega_{1},\omega_{2},\omega_{3}\,\right\\}.$

Here we ignored the last term in $H_{0}$ (it vanishes in the thermodynamic
limit), and the summation is over the transmon levels. The most divergent
terms correspond to the following sequences of transmon virtual states,
$0\shortrightarrow 1\shortrightarrow 0\shortrightarrow t\shortrightarrow 0$,
$0\shortrightarrow 1\shortrightarrow r\shortrightarrow 1\shortrightarrow 0$,
or $0\shortrightarrow s\shortrightarrow 0\shortrightarrow 1\shortrightarrow
0$, and to frequencies such that

$|\omega-\varepsilon_{1}|,|\omega_{1}-\varepsilon_{1}|,\omega_{2},\omega_{3}\ll\varepsilon_{1}$
(8)

(as well as two other inequalities after permutation of
$\omega_{1},\omega_{2},\omega_{3}$). The contribution to ${\cal
A}_{p;p_{1}p_{2}p_{3}}$ that corresponds to the on-shell condition,
$\omega=\omega_{1}+\omega_{2}+\omega_{3}$, and satisfies the inequalities (8),
is

$-\frac{2K^{2}\Delta^{2}}{\pi^{4}\varepsilon_{1}\sqrt{\omega_{2}\omega_{3}}}\frac{|W_{01}|^{2}}{(\omega-\varepsilon_{1})(\omega_{1}-\varepsilon_{1})}\sum_{s}\left(\frac{|W_{1s}|^{2}}{\varepsilon_{s}-\varepsilon_{1}}-\frac{|W_{0s}|^{2}}{\varepsilon_{s}}\right).$
(9)

By substituting the operator $\hat{n}$ with a gate charge $\cal N$ in Eq. (5),
one can calculate the gate sensitivity of the energy levels
$\varepsilon_{1,0}({\cal N})$ perturbatively in $\cal N$, and identify
footnote7

$\sum_{s}\left(\frac{|W_{1s}|^{2}}{\varepsilon_{s}-\varepsilon_{1}}-\frac{|W_{0s}|^{2}}{\varepsilon_{s}}\right)=\frac{1}{2}\left.\frac{\partial^{2}[\varepsilon_{1}({\cal
N})-\varepsilon_{0}({\cal N})]}{\partial{\cal N}^{2}}\right|_{{\cal N}=0}.$
(10)

In the transmon limit [cf. Eqs. (1) and (2)], $\varepsilon_{s}({\cal
N})=\varepsilon_{s}(1/4)+\lambda_{s}\cos(2\pi{\cal N})$, allowing us to
replace the right-hand side of the exact Eq. (10) with
$2\pi^{2}(\lambda_{1}-\lambda_{0})\approx 2\pi^{2}\lambda_{1}$. Evaluating all
remaining factors in Eq. (9) within the transmon’s harmonic approximation, in
which $\varepsilon^{(0)}_{s}\approx s\omega_{0}$ and
$W_{sr}^{(0)}=8\tilde{E}_{C}(\tilde{E}_{J}/32\tilde{E}_{C})^{1/4}[\sqrt{s}\delta_{s,r+1}+\sqrt{s+1}\delta_{s,r-1}]$,
we find the leading contribution to Eq. (6) for frequencies satisfying the
conditions (8),

$\gamma_{\rm
si}(\omega_{1},\omega_{2},\omega_{3}|\omega)=\frac{16}{3}\frac{K^{2}\Gamma^{2}\lambda_{1}^{2}}{\omega_{2}\omega_{3}(\omega_{1}-\omega_{0})^{2}(\omega-\omega_{0})^{2}}.$
(11)

The partial cross-section of quasielastic scattering, which characterizes the
transmon fluorescence,

$\gamma_{\rm
si}(\omega^{\prime}|\omega)=\int_{0}^{\omega-\omega^{\prime}}d\omega_{2}\gamma_{\rm
si}(\omega^{\prime},\omega_{2},\omega-\omega^{\prime}-\omega_{2}|\omega)\,,$
(12)

is obtained at $|\omega-\omega_{0}|,|\omega^{\prime}-\omega_{0}|\ll\omega_{0}$
from Eq. (11) and two other equations with permutations of
$\omega_{1},\omega_{2},\omega_{3}$, as

$\gamma_{\rm
si}(\omega^{\prime}|\omega)=\frac{32K^{2}\Gamma^{2}\lambda_{1}^{2}}{(\omega-\omega_{0})^{2}(\omega^{\prime}-\omega_{0})^{2}(\omega-\omega^{\prime})}\ln\left(\frac{\omega-\omega^{\prime}}{\Delta}\right);$
(13)

here we used $\Delta$ as a low-frequency cut-off.

We contrast Eq. (13) with the one obtained for a weak anharmonic oscillator,
after approximating $1-\cos\varphi\approx\varphi^{2}/2-\varphi^{4}/24$ in the
Josephson term of Eq. (3) and neglecting the phase slips. In that case, the
gate sensitivity is absent, the right-hand side of the identity (10) is zero,
rendering $\gamma_{\rm si}(\omega^{\prime}|\omega)=0$. In harmonic
approximation, Eq. (7) vanishes identically. Treating the anharmonic
corrections to $\varepsilon_{s}$ and $W_{sr}$ appearing in Eq. (7)
perturbatively, and assuming the incident photon to be close to resonance,
$|\omega-\omega_{0}|\ll\omega_{0}$, yields

$\gamma_{\rm{q}}(\omega_{1},\omega_{2},\omega_{3}|\omega)=\frac{({256}/{3\pi^{2}})\Gamma^{4}\tilde{E}_{C}^{2}\omega^{3}_{0}\omega_{1}\omega_{2}\omega_{3}}{[(\omega_{0}-\omega)(\omega_{0}^{2}-\omega_{1}^{2})(\omega_{0}^{2}-\omega_{2}^{2})(\omega_{0}^{2}-\omega_{3}^{2})]^{2}}.$
(14)

Assuming that the outgoing photon is also close to resonance,
$|\omega^{\prime}-\omega_{0}|\ll\omega_{0}$, we get

$\gamma_{\rm{q}}(\omega^{\prime}|\omega)=\frac{32}{9\pi^{2}}\frac{\Gamma^{4}\tilde{E}_{C}^{2}(\omega-\omega^{\prime})^{3}}{\omega_{0}^{6}(\omega-\omega_{0})^{2}(\omega^{\prime}-\omega_{0})^{2}}.$
(15)

The comparison of Eqs. (13) and (15) shows that phase slips are much more
effective in coupling the resonant modes to the low-frequency ones than the
anharmonic corrections to the qubit levels. The low-frequency modes, being far
away from the resonance, do not hybridize well with the qubit. The phase slips
are free from that drawback (at the expense of a potentially small value of
$\lambda_{1}$). Thus phase slips dominate in the inelastic processes at
$\omega^{\prime}\to\omega$.

Differential cross-section.– We proceed further by accounting for higher-order
quasielastic processes at finite $K$. The dichotomy between the high-frequency
photon modes that are in resonance with the transmon and the low-frequency
modes motivates a two-band approximation:

$\displaystyle H_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\omega_{0}|1\rangle\langle
1|+\sum_{p>p_{c}}\omega_{p}|p\rangle\langle
p|+\sum_{p>p_{c}}\left[t|p\rangle\langle 1|+{\rm H.c.}\right]$ (16)
$\displaystyle+\sum_{0<p<p_{c}}\omega_{p}a^{\dagger}_{p}a_{p}+\lambda_{1}|1\rangle\langle
1|\cos(2\pi\tilde{n}).$

Here the first line represents the hybridization of the transmon with high-
frequency photons; $|p\rangle=a^{\dagger}_{p}|0\rangle$, where $|0\rangle$ is
the ground state, $|1\rangle$ is the state in which (only) the transmon is
excited, and $t=\sqrt{\Gamma\Delta/\pi}$ is the hybridization matrix element.
The second line in Eq. (16) accounts for the coupling of low-frequency photons
and the transmon excited state through phase slips footnote-a ; the local
charge operator $\tilde{n}$ differs from Eq. (4) by the restriction of the sum
to low-frequency modes, $0<p<p_{c}$. The separation between low- and high-
frequency photons is set by frequency $\omega_{c}=vp_{c}$ such that
$\Gamma\ll\omega_{0}-\omega_{c}\ll\omega_{0}$.

The first line in the Hamiltonian (16) is equivalent to the Fano-Anderson
model. Its eigenstate $|k\rangle$ with energy $\omega_{k}$ such that
$\omega_{k}-\omega_{0}=-\Gamma\tan({\omega_{k}d}/{v})$ has an overlap with the
transmon state

$\beta^{2}_{k}\equiv|\langle
1|k\rangle|^{2}=\frac{{{\Gamma\Delta/\pi}}}{(\omega_{k}-\omega_{0})^{2}+\Gamma^{2}}\,,$
(17)

assuming $\omega_{k}$ is close to the resonance,
$|\omega_{k}-\omega_{0}|\ll\omega_{0}$ footnote1 . In new variables,
Hamiltonian (16) then reads

$\displaystyle H_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\sum_{k>p_{c}}\omega_{k}|k\rangle\langle
k|+\sum_{0<p<p_{c}}{\omega_{p}}a^{\dagger}_{p}a_{p}+H_{1},$ $\displaystyle
H_{1}$ $\displaystyle=$
$\displaystyle\lambda_{1}\sum_{k,k^{\prime}>p_{c}}\beta_{k}\beta_{k^{\prime}}|k\rangle\langle
k^{\prime}|\cos(2\pi\tilde{n}).$ (18)

It is the “backaction” of the qubit on the dynamic charge $\tilde{n}$ that
leads to the emission of “soft” photon modes by the resonant ones footnote8 .
Using Eq. (Critical fluorescence of a transmon at the Schmid transition), we
apply Fermi’s Golden Rule to calculate the (quasielastic) fluorescence
spectrum perturbatively in $\lambda_{1}$,

$\gamma_{\rm
si}(\omega_{{k^{\prime}}}|\omega_{{k}})=\frac{2\pi^{2}}{\Delta^{2}}\sum_{f}\left|\langle
k^{\prime},f|H_{1}|k,0\rangle\right|^{2}\delta(\omega_{{k}}-\omega_{{k^{\prime}}}-E_{f}).$
(19)

Here $|k,f\rangle=|k\rangle\otimes|f\rangle$ and $|f\rangle$ is a multiphoton
state with energy $E_{f}$ formed out of low-frequency photon modes. By
standard manipulations, we express Eq. (19) in terms of a photon correlation
function for an array disconnected from a transmon,

$\gamma_{\rm
si}(\omega^{\prime}|\omega)=\frac{\lambda_{1}^{2}\Gamma^{2}/\pi}{[(\omega-\omega_{0})^{2}+\Gamma^{2}][(\omega^{\prime}-\omega_{0})^{2}+\Gamma^{2}]}{\cal
C}(\omega-\omega^{\prime})$ (20)

with

$\\!{\cal C}(\Omega)=2{\rm
Re}\int_{0}^{\infty}\\!\\!\\!dte^{i(\Omega+i0^{+})t}\langle\cos
2\pi\tilde{n}(t)\cos 2\pi\tilde{n}(0)\rangle.$ (21)

Here
$\tilde{n}(t)=(1/\pi)\sum_{0<p<p_{c}}f_{p}(a_{p}e^{-i\omega_{p}t}+a_{p}^{\dagger}e^{i\omega_{p}t})$.
We use extensively Baker-Hausdorff formula Baker  to find

${\cal C}(\Omega)=2e^{-4\sum_{p}f_{p}^{2}}{\rm
Re}\int_{0}^{\infty}\\!\\!\\!dte^{i(\Omega+i0^{+})t}\cosh\left(4\sum_{p}f_{p}^{2}e^{-i\omega_{p}t}\right)$
(22)

at zero temperature. The Taylor expansion of the $\cosh$-factor in Eq. (22)
allows interpreting Eq. (20) as a partial cross-section of a high-frequency
photon scattering into another high-frequency photon, while an even number of
low-frequency photons is produced. At $K\to 0$, Eq. (20) reproduces Eq. (13)
upon the renormalization footnote4  of the phase slip amplitude,
$\lambda_{1}\to\lambda_{1}e^{-2\sum_{p}f_{p}^{2}}\approx\lambda_{1}(\Delta/\omega_{0})^{2K}$,
and not too close to the resonance,
$\Gamma\ll|\omega-\omega_{0}|,|\omega^{\prime}-\omega_{0}|\ll\omega_{0}$.

Being proportional to $\lambda_{1}^{2}(\Delta/\omega_{0})^{4K}$, the three-
photon amplitude vanishes in the thermodynamic limit at finite $K$. Thus
higher-order processes should be included. Instead of evaluating and summing
them we note that at $\Delta\to 0$ Eq. (21) can be simplified,

${\cal C}(\Omega)\approx{\rm
Re}\int_{0}^{\infty}\\!\\!\\!dte^{i(\Omega+i0^{+})t}\langle
e^{{i}2\pi\tilde{n}(t)}e^{-{i}2\pi\tilde{n}(0)}\,\rangle\,.$ (23)

This correlator has been much studied Weiss1985 ; Schon1990 ,

${\cal
C}(\Omega)=\frac{\pi}{\Gamma(4K)}\frac{1}{\Omega}\left(\frac{\Omega}{\omega_{0}}\right)^{4K}e^{-\Omega/\omega_{0}}\,,\quad{\Omega>0}\,.$
(24)

Here $\Gamma(4K)$ is the Gamma function. Equation (20) with ${\cal C}(\Omega)$
of Eq. (24) is our main result footnote-b . It relates the fluorescence
spectrum with the dynamical phase-boost susceptibility $\propto{\cal
C}(\omega)$ at $\omega\ll\omega_{0}$ Guinea1985  in the same range of validity
defined by Eq. (1).

Inserting Eq. (24) in Eq. (20) we find that at resonant excitation,
$|\omega-\omega_{0}|\lesssim\Gamma$, and in the frequency range
$\omega_{0}\gg\omega_{0}-\omega^{\prime}\gg\Gamma$ of the emitted photons, the
fluorescence intensity is a power-law footnote5  of
$\omega_{0}-\omega^{\prime}$,

${\gamma_{\rm
si}(\omega^{\prime}|\omega)=\frac{1}{\Gamma(4K)}\frac{\lambda_{1}^{2}}{\omega_{0}^{3}}\left(\frac{\omega_{0}-\omega^{\prime}}{\omega_{0}}\right)^{4K-3}}\,.$
(25)

The perturbative-in-$\lambda_{1}$ result for the differential cross-section of
quasielastic scattering works at any $K$, except its smallest values allowing
for ${\lambda_{1}\gtrsim\Gamma}$. The behavior of Eq. (25) parallels the one
of the dynamical susceptibility $A(\omega)$ of Eq. (10) in Ref. Guinea1985b .
These two quantities are not normalizable at $K<1/2$, which is the signature
of the charge localized phase. At the critical point ($K=1/2$), we find
$\gamma(\omega^{\prime}|\omega_{0})\propto 1/(\omega_{0}-\omega^{\prime})$ as
the dynamical critical signature of the Schmid transition.

Total inelastic cross-section.– Finally, we show that at $K<1/2$ the
quasielastic transitions may yield the main contribution to the total
inelastic cross-section footnote6 ,

$\gamma(\omega)=\int_{0}^{\omega}d\omega^{\prime}\gamma(\omega^{\prime}|\omega)\,.$
(26)

Indeed, at $K<1/2$ the dominant contribution of the partial cross-section (20)
to the integral comes from a vicinity of order $\Gamma$ near its upper bound.
We may thus extend the lower bound in Eq. (26) to $-\infty$ and evaluate the
quasielastic component of the total cross-section as

$\gamma_{\rm si}(\omega_{0})=\frac{\pi}{2\sin(2\pi
K)\Gamma(4K)}\frac{\lambda_{1}^{2}}{\Gamma^{2}}\left(\frac{\Gamma}{\omega_{0}}\right)^{4K}$
(27)

for the incoming photons lying within the width of the resonance. Furthermore,
the inelastic spectral linewidth is asymmetric, see Fig. 2, with asymptotes

$\gamma_{\rm
si}(\omega)=\frac{\pi\lambda_{1}^{2}\Gamma/\omega_{0}^{3}}{\Gamma(4K)}\left(\frac{\omega_{0}}{\omega-\omega_{0}}\right)^{3-4K}$
(28)

at $\omega-\omega_{0}\gg\Gamma$ and

$\gamma_{\rm
si}(\omega)=\frac{\pi(1-4K)\lambda_{1}^{2}\Gamma^{2}/\omega_{0}^{4}}{\sin(4\pi
K)\Gamma(4K)}\left(\frac{\omega_{0}}{\omega_{0}-\omega}\right)^{4-4K}$ (29)

at $\omega_{0}-\omega\gg\Gamma$. Importantly, the found $\gamma_{\rm
si}(\omega)$ is independent of the artificially-introduced partition frequency
$\omega_{c}$ thus justifying the use of Eq. (16).

To assess the contribution to Eq. (26) of deeply-inelastic processes, we use
the differential cross-section (14) stemming from the transmon anharmonicity
and favoring large energy transfer between the incoming and any of the
outgoing photons. The corresponding contribution to the total cross-section is

$\gamma_{\rm
q}(\omega_{0})=\alpha\frac{\Gamma^{2}\tilde{E}_{C}^{2}}{\omega_{0}^{4}}\quad\mathrm{with}\quad\alpha\approx
0.45.$ (30)

The comparison with Eq. (27) shows that quasielastic processes dominate the
total inelastic cross-section $\gamma(\omega)$ if
$\lambda_{1}/\tilde{E}_{C}\gg(\Gamma/\omega_{0})^{2(1-K)}$. Under the
assumption $\lambda_{1}\ll\Gamma$, this condition is possible to satisfy only
at $K<1/2$: in terms of the Schmid transition, the phase-slip mechanism may
dominate the total inelastic cross-section only in the charge-localized phase.

While Eqs. (20) and (24) remain valid at $K>1/2$, their use in evaluation of
$\gamma(\omega)$ is not justified: the dominant contribution to the integral
in Eq. (26) at $K>1/2$ comes from $\omega-\omega^{\prime}\gtrsim\omega_{c}$
and depends on $\omega_{c}$, rendering the model (16) inapplicable. This is
consistent with our perturbative analysis of Eq. (5): there is no parameter
allowing to single out the phase-slip-induced transitions from other processes
at energy losses comparable to $\omega_{0}$.

Figure 2:  Frequency dependence of the total cross-section at zero temperature
for $K=0.01$ (blue), 0.1 (orange), 0.2 (green), 0.25 (red), 0.3 (purple), and
0.4 (brown).

It is straightforward to generalize Eqs. (20) and (27)-(29) to finite
temperatures $T$, by using the corresponding finite-$T$ generalization
Weiss1985 ; Schon1990  of ${\cal C}(\Omega)$. A low temperature,
$T\ll\omega_{0}$, leaves the power-law spectrum (25) intact at
$\omega_{0}-\omega^{\prime}\gg T$. The total inelastic cross-section in the
scaling region [$\Gamma,|\omega-\omega_{0}|,T\ll\omega_{0}$] is found as

$\frac{\gamma(\omega,T)}{\gamma(\omega_{0},0)}=\frac{\sin(2\pi
K)}{\pi^{2}}\frac{\tau^{4K-1}}{1+\nu^{2}}\int_{0}^{\infty}dx\frac{e^{\pi
x/\tau}|\Gamma(2K+ix/\tau)|^{2}}{(x-\nu)^{2}+1}$ (31)

with $\tau=2\pi T/\Gamma$ and $\nu=(\omega-\omega_{0})/\Gamma$. Its
temperature dependence at resonance, $\omega=\omega_{0}$, is shown in Fig. 3.
The total cross-section increases (decreases) with the temperature at $K>1/4$
($K<1/4$). At low temperature $T\ll\Gamma$,

$\frac{\gamma(\omega_{0},T)}{\gamma(\omega_{0},0)}\approx
1-\frac{1}{2}\left(\frac{2\pi T}{\Gamma}\right)^{4K},$ (32)

where the second term is a significant correction in a wide temperature range
at $K\ll 1$.

Figure 3:  Temperature dependence of the total cross-section at resonance
(same color code as in Fig. 2).

Conclusion.– We believe the spectrum of the fluorescence that we predict, see
Eqs. (20)-(25), charts an interesting direction for future experiments, while
the found total inelastic cross-section, see Eqs. (26)-(32), is directly
related to the ongoing experiments Kuzmin2020  in the spirit of Refs.
Puertas2019 ; Kuzmin2019 ; Leger2019 . Indeed, the internal quality factor
$Q(\omega)$ of a discrete mode in a finite-length array, which is routinely
measured in such experiments, can be expressed in terms of an inelastic decay
rate, $Q(\omega)=\omega/\Gamma_{\rm in}(\omega)$. The latter is related to the
total inelastic cross-section through $\Gamma_{\rm
in}(\omega)=\gamma(\omega)\Delta/\pi$. Lastly, in the weakly nonlinear
transmon regime, which we focussed upon, photon scattering remains mostly
elastic. As the nonlinearity increases, we may anticipate large inelastic
cross-sections, which would manifest a different kind of quantum impurity
problem than the Kondo regime studied in Refs. LeHur2012 ; Goldstein2013 .

###### Acknowledgements.

We acknowledge stimulating discussions with M. Goldstein and with V.
Manucharyan who attracted our attention to the role of phase slips, as well as
B. Beri and N. Roch for useful comments on the manuscript. This work is
supported by the DOE contract DE-FG02-08ER46482 (LG), and by ANR through Grant
No. ANR-16-CE30-0019 and ARO grant W911NF-18-1-0212 (MH).

## References

  * (1) A. Schmid, Phys. Rev. Lett. 51, 1506 (1983).
  * (2) S. A. Bulgadaev, Pisma Zh. Eksp. Teor. Fiz. 39, 264 (1984) [JETP Lett. 39, 315 (1984)].
  * (3) F. Guinea, V. Hakim and A. Muramatsu, Phys. Rev. Lett. 54, 263 (1985).
  * (4) M. P. A. Fisher and W. Zwerger, Phys. Rev. B 32, 6190 (1985).
  * (5) G. Schön and A. D. Zaikin, Phys. Rep. 198, 237 (1990).
  * (6) R. Yagi, S. I. Kobayashi, and Y. Ootuka, J. Phys. Soc. Japan 66, 3722 (1997).
  * (7) J. S. Penttilä, P. J. Hakonen, E. B. Sonin, and M. A. Paalanen, J. Low Temp. Physics, 125, 89 (2001).
  * (8) M. Watanabe and D. B. Haviland, Phys. Rev. B 67, 094505 (2003).
  * (9) A. Murani, N. Bourlet, H. le Sueur, F. Portier, C. Altimiras, D. Esteve, H. Grabert, J. Stockburger, J. Ankerhold, and P. Joyez, Phys. Rev. X 10, 021003 (2020).
  * (10) Let us clarify that the fluorescence considered here can be studied at weak ac drive; it differs from the resonance fluorescence of a strongly driven qubit Astafiev2010  resulting in the appearance of Mollow triplet in the spectrum.
  * (11) J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 76, 042319 (2007).
  * (12) J. A. Schreier, A. A. Houck, J. Koch, D. I. Schuster, B. R. Johnson, J. M. Chow, J. M. Gambetta, J. Majer, L. Frunzio, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. B 77, 180502(R) (2008).
  * (13) N. A. Masluk, I. M. Pop, A. Kamal, Z. K. Minev, and M. H. Devoret, Phys. Rev. Lett. 109, 137002 (2012).
  * (14) M. T. Bell, I. A. Sadovskyy, L. B. Ioffe, A. Yu. Kitaev, and M. E. Gershenson, Phys. Rev. Lett. 109, 137003 (2012).
  * (15) J. Puertas Martínez, S. Léger, N. Gheeraert, R. Dassonneville, L. Planat, F. Foroughi, Y. Krupko, O. Buisson, C. Naud, W. Hasch-Guichard, S. Florens, I. Snyman, and N. Roch, npj Quantum Information 5, 19 (2019).
  * (16) R. Kuzmin, N. Mehta, N. Grabon, R. Mencia, and V. E. Manucharyan, npj Quantum Information 5, 20 (2019).
  * (17) S. Léger, J. Puertas-Martínez, K. Bharadwaj, R. Dassonneville, J. Delaforce, F. Foroughi, V. Milchakov, L. Planat, O. Buisson, C. Naud, W. Hasch-Guichard, S. Florens, I. Snyman, and N. Roch, Nature Commun. 10, 5259 (2019).
  * (18) J. J. García-Ripoll, E. Solano, and M. A. Martin-Delgado, Phys. Rev. B 77, 024522 (2008).
  * (19) K. Le Hur, Phys. Rev. B 85, 140506(R) (2012).
  * (20) M. Goldstein, M. H. Devoret, M. Houzet, and L. I. Glazman, Phys. Rev. Lett. 110, 017002 (2013).
  * (21) I. Snyman and S. Florens, Phys. Rev. B 92, 085131 (2015).
  * (22) N. Gheeraert, X. H. H. Zhang, T. Sépulcre, S. Bera, N. Roch, H. U. Baranger, and S. Florens, Phys. Rev. A 98, 043816 (2018).
  * (23) T. Yamamoto and T. Kato, J. Phys. Soc. Jpn. 88, 094601 (2019).
  * (24) F. W. J. Hekking and L. I. Glazman, Phys. Rev. B 55, 6551 (1997).
  * (25) R. M. Bradley and S. Doniach, Phys. Rev. B 30, 1138 (1984).
  * (26) S. E. Korshunov, Zh. Eksp. Teor. Fiz 95, 1058 (1989) [Sov. Phys. JETP 68, 609 (1989)].
  * (27) D. M. Basko, F. Pfeiffer, P. Adamus, M. Holzmann, and F. W. J. Hekking, Phys. Rev. B 101, 024518 (2020).
  * (28) V. Manucharyan, private communication, 2019.
  * (29) A. O. Gogolin, A. A. Nersesyan, and A. M. Tsvelik, Bosonization and Strongly Correlated Systems (Cambridge University Press, Cambridge, England, 1998).
  * (30) Equation (10) is reminiscent of a famous expression for the effective mass tensor in Bloch band theory, see Eq. (59.5) in [E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics, Part 2: Theory of the Condensed State (Pergamon Press, Oxford, 1980)].
  * (31) In Eq. (16) we dispense with the term $-\lambda_{0}|0\rangle\langle 0|\cos(2\pi\tilde{n})$, which is responsible for charge localization in the ground state ($\lambda_{0}\ll\lambda_{1}$).
  * (32) The density of states (DOS), ${(1/\pi)}dk/d\omega_{k}$, contains a Lorentzian peak corresponding to the resonant level. In a finite chain, it translates into a non-equidistant spectrum $\omega_{k}$, making the DOS width $2\Gamma$ an experimentally accessible parameter Kuzmin2019 .
  * (33) The resonance renders the effect of phase slips on inelastic scattering markedly different from the one considered in homogeneous Josephson-junction chains Bard2018 ; Houzet2019 .
  * (34) $e^{A}e^{B}=e^{A+B}e^{\frac{1}{2}[A,B]}$ when $[A,B]$ is a $c$-number.
  * (35) In this limit, the intermediate scale $\omega_{c}$ can be replaced with $\omega_{0}$, and thus drops from the result.
  * (36) U. Weiss and H. Grabert, Phys. Lett. A 108, 63 (1985).
  * (37) Equation (20) also holds at finite $\lambda_{0}$ provided that the quantum average in ${\cal C}(\Omega)$ is performed with $H_{\rm sG}$. At $K<1/2$, standard methods yield ${\cal C}(\Omega\ll\omega_{\star})\propto\Omega/\omega_{\star}^{2}$ with the soft gap $\omega_{\star}\sim\omega_{0}(\lambda_{0}/\omega_{0})^{1/(1-2K)}$ Gogolin-book . Thus charge localization in the ground state is unveiled by the suppression of $\gamma_{\rm si}$ at small frequency exchange. We safely ignore that structure in our results, as $\omega_{\star}\ll\Gamma$.
  * (38) Again, the chosen frequency range allowed us to replace $\omega_{c}$ that was initially appearing in Eq. (24) by $\omega_{0}$.
  * (39) F. Guinea, Phys. Rev. B 32, 4486 (1985).
  * (40) In Ref. Goldstein2013  we were invoking energy conservation instead,  $\omega\gamma(\omega)=\int_{0}^{\omega}d\omega^{\prime}\omega^{\prime}\gamma(\omega^{\prime}|\omega).$ (33) The two formulas agree with each other when quasielastic transitions dominate.
  * (41) R. Kuzmin, N. Grabon, N. Mehta, A. Burshtein, M. Goldstein, M. Houzet, L. I. Glazman, and V. E. Manucharyan, arXiv:2010.02099.
  * (42) O. Astafiev, A. M. Zagoskin, A. A. Abdumalikov, Y. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai, Science 327, 840 (2010).
  * (43) M. Bard, I. V. Protopopov, and A. D. Mirlin, Phys. Rev. B 98, 224513 (2018).
  * (44) M. Houzet and L. I. Glazman, Phys. Rev. Lett. 122, 237701 (2019).