# Many-body quantum chaos in mixtures of multiple species

Vijay Kumar and Dibyendu Roy  Raman Research Institute, Bangalore 560080,
India

###### Abstract

We study spectral correlations in many-body quantum mixtures of fermions,
bosons, and qubits with periodically kicked spreading and mixing of species.
We take two types of mixing, namely, Jaynes-Cummings and Rabi, respectively,
satisfying and breaking the conservation of a total number of species. We
analytically derive the generating Hamiltonians whose spectral properties
determine the spectral form factor in the leading order. We further analyze
the system-size $(L)$ scaling of Thouless time $t^{*}$, beyond which the
spectral form factor follows the prediction of random matrix theory. The
$L$-dependence of $t^{*}$ crosses over from $\log L$ to $L^{2}$ with an
increasing Jaynes-Cummings mixing between qubits and fermions or bosons in a
finite-sized chain, and it finally settles to $t^{*}\propto\mathcal{O}(L^{2})$
in the thermodynamic limit for any mixing strength. The Rabi mixing between
qubits and fermions leads to $t^{*}\propto\mathcal{O}(\log L)$, previously
predicted for single species of qubits or fermions without total number
conservation.

A series of recent microscopic studies has explored quantum chaos and spectral
correlations in periodically driven (Floquet) many-body systems [1, 2, 3, 4,
5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18] to show the emergence of
universal random matrix theory (RMT) description of the spectral form factor
(SFF) in these models by going beyond the semiclassical periodic-orbit
approaches [19, 20]. These investigations have further strengthened our
understanding of the quantum chaos conjecture [21, 22, 23, 24, 25, 26, 27, 28,
29] for describing the spectral fluctuations of many-body nonintegrable
quantum systems by RMT. Till now, such microscopic derivation of SFF in many-
body quantum models has been restricted to systems with single
components/species, e.g., fermions, bosons, and qubits. Nature, however, is
full of systems consisting of multiple species, such as the crystalline solids
of electrons and phonons and the black-body radiation comprising thermal
electromagnetic radiation within or surrounding a matter in thermodynamic
equilibrium. Inspired by these examples, we derive the leading order
contributions to SFF in various mixed many-body quantum systems with two
different types of species, e.g., qubits and bosons or fermions [30].

We consider many-body quantum mixtures where a base Hamiltonian with the
entries diagonal in the Fock space basis of two different species is kicked
periodically by another Hamiltonian with terms consisting of mixing between
two species and nearest-neighbor hopping of any one species. The diagonal
entries in the base Hamiltonian include random chemical potentials and
transition frequencies along with pairwise long-range interactions of one
species. We consider two forms of the mixing Hamiltonian: (a) Jaynes-Cummings
$(JC)$ [31, 32, 33, 34, 35] and (b) Rabi $(R)$ [36, 37] interaction between
different species. While the $JC$ preserves the total number of excitations of
both species, the $R$ does not. Thus, we have $U(1)$) symmetry in the $JC$
mixing system, which is absent for the $R$ mixing. Our models’ two different
components are either qubit and spinless boson or qubit and spinless fermion.
Since spinless fermions are related to spin-1/2s or qubits, our results here
are valid for many different types of mixture, e.g., the results for a
compound model of qubits and spinless bosons are also helpful for a mix
between spinless fermions and bosons. Similarly, the results for a mixture of
qubits and spinless fermions apply to a mixture of spin-1/2s of different
species, e.g., electrons and atomic nuclei in solids.

First, we rewrite the spectral form factor of the compound systems in terms of
a bi-stochastic many-body process [7, 13] generated by an effective
Hamiltonian. The effective Hamiltonian describes the leading order
contributions of SFF within the random phase approximation (RPA) in the
Trotter regime of small perturbation parameters. We identify symmetries of the
effective Hamiltonian controlling dynamical processes for the emergence of RMT
behavior in these models [13, 38]. These symmetries are important in
determining system-size $(L)$ scaling of the Thouless timescales $t^{*}$
beyond which the SFF has a universal RMT/COE form for our time-reversal
invariant models of a circular orthogonal ensemble (COE). For $JC$ mixing, we
find, $t^{*}\propto L^{2}$ when $L\to\infty$, which is a characteristics of
$U(1)$-symmetric model [39, 5, 7, 13]. However, we show an exciting
competition between the hopping and mixing of the driving Hamiltonian, leading
to a crossover behavior in the $L$-dependence of $t^{*}$ when a finite-size
system is considered. For a finite system, $t^{*}\propto\log L$ when the
mixing strength is smaller than the hopping, and $t^{*}\propto L^{2}$ for a
higher mixing strength compared to hopping. The above crossover in $L$ scaling
of $t^{*}$ is inevitable in many experimental studies with highly controlled
laboratory settings of finite size [35, 40, 41, 42, 43, 44, 45, 18]. For $R$
mixing between fermions and qubits, $t^{*}\propto\log L$ or $L^{0}$ for large
$L$, which is similar to the single species of fermion or spin-1/2 models in
the absence of $U(1)$ symmetry. In contrast to fermions or spin-1/2s, the only
boson model lacking $U(1)$ symmetry shows an algebraic $L$-dependence of
$t^{*}$ [13]. We offer numerical evidence that the $L$-dependence of $t^{*}$
for $R$ mixing between bosons and qubits seems to behave similarly to $R$
mixing between fermions and qubits.

The base (kicked) Hamiltonian $\hat{H}_{0}$ of our systems denotes a one-
dimensional lattice of length $L$ consisting of spinless fermions or bosons
and qubits with no coupling between these two entities/species.

$\displaystyle\hat{H}_{0}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{L}(\omega_{i}\hat{n}_{i}+\Omega_{i}\hat{\sigma}^{\dagger}_{i}\hat{\sigma}_{i})+\sum_{i<j}U_{ij}\hat{n}_{i}\hat{n}_{j},$
(1)

where $\hat{n}_{i}=\hat{a}_{i}^{\dagger}\hat{a}_{i}$ is the number operator
with $\hat{a}_{i}^{\dagger}$ being a fermion or boson creation operator at
site $i$. The raising and lowering operators
$\hat{\sigma}^{\dagger}_{j}\equiv(\hat{\sigma}_{j}^{x}+i\hat{\sigma}^{y}_{j})/2,\hat{\sigma}_{j}\equiv(\hat{\sigma}_{j}^{x}-i\hat{\sigma}^{y}_{j})/2$
are for the qubit at site $j$. Here, $\omega_{i}$ and $\Omega_{i}$ are,
respectively onsite energy/frequency of the fermion/boson and the transition
frequency of the qubit at site $i$. We choose one or both of $\omega_{i}$ and
$\Omega_{i}$ random as Gaussian iid variables of zero mean and finite standard
deviation. We further take long-range interaction between fermions or bosons
at sites $i$ and $j$ given by $U_{ij}=U_{0}/|i-j|^{\alpha}$ with an exponent
in the interval $1<\alpha<2$. The form of $\hat{H}_{0}$ is fixed by minimal
requirements for analytical calculation as well as physical relevance. Our
analytical calculation requires the RPA and integration out of the parameters
of $\hat{H}_{0}$, and both are met by the above choice of $\hat{H}_{0}$. The
model with bosons and qubits and its close variants can physically represent
light-matter interactions in real systems and engineered meta-materials [32,
35, 37] and electron-phonon interactions in crystalline solids.

The driving/kicking Hamiltonian consists of a term denoting the mixing between
fermions/bosons and qubits locally and another term indicating nearest-
neighbor hopping of fermions/bosons. The driving Hamiltonian with $JC$ and $R$
interactions are, respectively,

$\displaystyle\hat{H}_{\text{J}C}=\sum_{i=1}^{L}g(\hat{a}^{\dagger}_{i}\hat{\sigma}_{i}+\hat{\sigma}^{\dagger}_{i}\hat{a}_{i})+\sum_{i=1}^{L}(-J\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+{\rm
h.c.}),$ (2)
$\displaystyle\hat{H}_{\text{R}}=\sum_{i=1}^{L}g(\hat{a}^{\dagger}_{i}+\hat{a}_{i})(\hat{\sigma}_{i}+\hat{\sigma}^{\dagger}_{i})+\sum_{i=1}^{L}(-J\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+{\rm
h.c.}),$ (3)

where $g$ and $J$ are the strength of mixing and hopping. The total excitation
number operator,
$\hat{N}=\sum_{i=1}^{L}(\hat{n}_{i}+\hat{\sigma}^{\dagger}_{i}\hat{\sigma}_{i})$,
commutes with both $\hat{H}_{0}$ and $\hat{H}_{\text{J}C}$, but not with
$\hat{H}_{\text{R}}$. Thus, the time-dependent total Hamiltonian,
$\hat{H}(t)=\hat{H}_{0}+\hat{H}_{\text{J}C/R}\sum_{m\in\mathbb{Z}}\delta(t-m)$,
commutes with $\hat{N}$ for $JC$ interaction but not for $R$ interaction
showing the presence or absence of a $U(1)$ symmetry, which corresponds
respectively to conservation or violation of the total excitation number in
our models. We here use periodic boundary condition (PBC) in real space, i.e.,
$\hat{a}_{i}\equiv\hat{a}_{i+L},\hat{\sigma}_{i}\equiv\hat{\sigma}_{i+L}$.

The SFF, $K(t)$, is defined as a time $(t)$ Fourier transformation of the two-
point correlation of the spectral density of quasienergies, which are
eigenvalues of the unitary one-cycle Floquet propagator $\hat{U}$ of our
periodically driven systems. $K(t)$ can be written as [1, 7]

$\displaystyle
K(t)=\langle(\text{tr}\hat{U}^{t})(\text{tr}\hat{U}^{-t})\rangle-(\mathcal{N}_{\zeta}^{\beta})^{2}\delta_{t,0},$
(4)

where $\mathcal{N}_{\zeta}^{\beta}$ is the dimension of the Hilbert space of
the system with $\zeta=JC,R$ mixing for fermions $(\beta=F)$ and bosons
$(\beta=B)$. Here, $\langle\dots\rangle$ denotes an average over the quench
disorders $\\{\Omega_{i}\\}$ and/or $\\{\omega_{i}\\}$. The one-cycle time-
evolution operator $\hat{U}$ can be expressed as

$\displaystyle\hat{U}=\hat{V}\hat{W},\quad\hat{W}=e^{-i\hat{H}_{0}}\text{ and
}\hat{V}=e^{-i\hat{H}_{\text{J}C/R}}.$ (5)

We consider the basis states
$|\underline{n\sigma}\rangle\equiv|n_{1},\ldots,n_{L}\rangle\otimes|\sigma_{1},\ldots,\sigma_{L}\rangle,$
where the occupation number of spinless fermion/boson and qubit at the lattice
site $j$ are respectively given $n_{j}=0,1~{}(F)~{}{\rm
and}~{}0,1,2,\ldots(B)$, and $\sigma_{j}=0,1$. The total number of excitations
$N\equiv\langle\underline{n\sigma}|\hat{N}|\underline{n\sigma}\rangle=\sum_{j=1}^{L}(n_{j}+\sigma_{j})$
is conserved in the whole system only for $JC$ mixing.

For $JC$ mixing of fermions and qubits, we can distribute total excitations
$N(<2L)$ among $2L$ states consisting of $L$ spatially localized qubit
excitations and another $L$ spatially delocalized fermionic excitations. Thus,
the dimension of the Hilbert space for this system with $N$ excitations
$\mathcal{N}_{\text{J}C}^{\text{F}}=(2L)!/((2L-N)!N!)$. We further have,
$\mathcal{N}_{\text{R}}^{\text{F}}=\sum_{N=0,2..}^{2L}\mathcal{N}_{\rm
JC}^{\rm F}=2^{2L-1}$, which is the dimension of the even sector of Hilbert
space for $R$ mixing of fermions and qubits.

For $JC$ mixing between bosons and qubits, the number of qubit excitations $M$
($\equiv\sum_{j=1}^{L}\sigma_{j}$) can be $0\leq M\leq\min(N,L)$. The total
number of bosons there would be $N-M$. We can find the dimension
$\mathcal{N}_{\rm JC}^{\rm B}$ of the Hilbert space by summing over allowed
$M$. Thus, we get

$\mathcal{N}_{\text{J}C}^{\text{B}}=\sum_{M=0}^{\min(N,L)}\frac{L!}{M!(L-M)!}\frac{(N-M+L-1)!}{(N-M)!(L-1)!}.$
(6)

The Hilbert space dimension $\mathcal{N}_{\text{R}}^{\text{B}}$ becomes
infinite for $R$ mixing of bosons and qubits as $N$ is not conserved and has
no upper bound. However, as discussed later, it is possible to introduce a
truncation for a maximum number of total excitation $N_{max}$ in the lattice
for numerical calculation.

Figure 1: Spectral form factor $K(t)$ using Eq. 7 for different system sizes
$L$ of the kicked chain with $JC$ mixing between fermions and qubits for
$g=0.1,J=0.4$ in (a,b), and $g=0.4,J=0.1$ in (c,d). We take half-filling
$N/L=1/2$. In (b) and (d), we show data collapse in scaled time $t/\log L$ and
$t/L^{1.85}$, respectively.

Both for fermionic and bosonic models, these basis states
$|\underline{n\sigma}\rangle$ are eigenstates of $\hat{H}_{0}$ and $\hat{W}$,
which allows us to integrate out $\hat{H}_{0}$ from $\hat{U}$ and $K(t)$
through the RPA by disorder averaging over different realizations. We further
make an identity permutation approximation to achieve the following simple
form of the SFF [7, 13] by including the leading order contributions at times
$t\ll t_{\rm H}\equiv\mathcal{N}_{\zeta}^{\beta}$:
$K(t)=2t\>\text{tr}\mathcal{M}^{t}$, where $\mathcal{M}$ is a
$\mathcal{N}_{\zeta}^{\beta}\times\mathcal{N}_{\zeta}^{\beta}$ double
stochastic square matrix whose elements are
$\mathcal{M}_{\underline{n\sigma},\underline{n\sigma}^{\prime}}=|\langle\underline{n\sigma}|\hat{V}|\underline{n\sigma}^{\prime}\rangle|^{2}$.
The largest eigenvalue of $\mathcal{M}$ is one due to the unitarity of
$\hat{V}$. Thus, we can write the eigenvalues of $\mathcal{M}$ as
$1,\lambda_{1},\lambda_{2},\lambda_{3},\dots$ with
$1\geq|\lambda_{i}|\geq|\lambda_{i+1}|$. Using these eigenvalues, we express
the SFF as (see Sec. I of SM[46] for a derivation):

$\displaystyle
K(t)=2t\big{(}1+\sum_{i=1}^{\mathcal{N}_{\zeta}^{\beta}-1}\lambda_{i}^{t}\big{)},$
(7)

where $K(t)\simeq 2t$ is a leading order in $t/t_{\rm H}$ result of RMT/COE.
The RMT/COE form of $K(t)\simeq 2t$ in a leading order appears beyond the
Thouless timescales $t^{*}(L)$ when the contribution from the second term in
Eq. 7 becomes negligible. The contribution from the second term depends on the
properties of $\lambda_{i}$ for $i\geq 1$. We next try to understand the
features of $\mathcal{M}$ and its eigenvalues. We can find Hermitian quantum
Hamiltonians generating $\mathcal{M}$ in the Trotter regime of small $g,J$ for
fermionic and bosonic models with $JC$ and $R$ mixing. The Hamiltonians are
derived by writing $\mathcal{M}$ using an element-wise commutative product
(also known as the Hadamard product) of $\hat{V}$ with $\hat{V}^{*}$ in the
basis $|\underline{n\sigma}\rangle$, and then expanding $\hat{V}$ in the
Trotter regime of small parameters of $\hat{H}_{\rm JC}~{}(\hat{H}_{\rm R})$
up to second order in $\hat{H}_{\rm JC}~{}(\hat{H}_{\rm R})$. The emergent
symmetries of these generating Hamiltonians control the dynamical properties
of these models, such as $t^{*}(L)$, and they can be significantly different
from the symmetries of $\hat{H}(t)$.

We first analyze $\mathcal{M}$ for the fermionic model with $JC$ mixing. The
generating Hamiltonian for PBC is (see Sec. II of SM[46] for a derivation)

$\displaystyle\mathcal{M}_{\text{J}C}^{\text{F}}$ $\displaystyle=$
$\displaystyle\big{(}1-\frac{(g^{2}+J^{2})L}{2}\big{)}\mathds{1}_{\mathcal{N}_{\text{J}C}^{\text{F}}}+\sum_{i=1}^{L}\sum_{\nu}\Big{(}\frac{J^{2}}{2}\hat{\tau}_{i}^{\nu}\hat{\tau}_{i+1}^{\nu}$
(8) $\displaystyle+$
$\displaystyle\frac{g^{2}}{2}\hat{\tau}_{i}^{\nu}\hat{\sigma}_{i}^{\nu}\Big{)}+\mathcal{O}(J^{4},g^{4}),$

where $\hat{\tau}^{\nu}_{i}$ and $\hat{\sigma}^{\nu}_{i}$ are the $\nu$th
component of Pauli matrix at site $i$ and $\nu\in\\{x,y,z\\}$. Here,
$\hat{\tau}^{\nu}_{i}$ and $\hat{\sigma}^{\nu}_{i}$ represent, respectively,
the spinless fermions and qubits. The largest eigenvalue one of
$\mathcal{M}_{\text{J}C}^{\text{F}}$ corresponds to a state in which all
$\tau$ and $\sigma$ spins are polarized in one particular direction, say along
$z$ axis. $\mathcal{M}_{\text{J}C}^{\text{F}}$ commutes with the operators,
$\sum_{j=1}^{L}(\hat{\tau}^{\nu}_{j}+\hat{\sigma}^{\nu}_{j})/2$ for
$\nu\in\\{x,y,z\\}$, which satisfy $SU(2)$ algebra. Thus,
$\mathcal{M}_{\text{J}C}^{\text{F}}$ has $SU(2)$ symmetry, which implies that
there would be degenerate symmetry multiplets of the subleading eigenvalues of
$\mathcal{M}_{\text{J}C}^{\text{F}}$ for different $N~{}(=1,2,3,\dots 2L-1)$.
Nevertheless, other energy eigenvalues can also appear between different
descendent states for higher $N$. Since we are interested in $L$-dependence of
$t^{*}$ at finite filling fractions $(N/L)$, the ordering of descendant states
in the full spectrum of $\mathcal{M}_{\text{J}C}^{\text{F}}$ for $N>1$ is
important. It can be shown that the value of $\lambda_{1}$ is the same for all
$N$, including $N=1$ at any value of $g,J$. The eigenvalues of
$\mathcal{M}_{\text{J}C}^{\text{F}}$ excluding the largest eigenvalue one for
$N=1$ are

$\displaystyle\lambda_{i}=1-g^{2}-J^{2}\big{(}1-\cos\frac{2i\pi}{L}\big{)}+\sqrt{J^{4}\big{(}1-\cos\frac{2i\pi}{L}\big{)}^{2}+g^{4}},$
(9)

for $i=1,2,3\dots,L-1$. In the thermodynamic limit of $L\to\infty$, we find
from Eq. 9, $\lambda_{1}\approx 1-(2\pi^{2}J^{2})/L^{2}$ for any value of
$g,J$. However, such approximation of Eq. 9 is also applicable for large
$L~{}(>{\it l}_{c}\equiv\pi/\sin^{-1}(g/\sqrt{2}J)$, a critical length-scale
depending on $g,J$) when $(1-\cos\frac{2\pi}{L}\big{)}\ll(g/J)^{2}$. We then
further approximate $K(t)$ at long time $t$, $1\ll
t\ll\mathcal{N}_{\text{J}C}^{\text{F}}$, by keeping up to the second largest
eigenvalue $\lambda_{1}$ of $\mathcal{M}_{\text{J}C}^{\text{F}}$. Thus, we get
for SFF

$\displaystyle K(t)\simeq 2t(1+\lambda_{1}^{t})\simeq 2t(1+e^{-t/t^{*}(L)}),$
(10)

where we take the scaling of $\lambda_{1}$ with system size $L$ as $1-1/t^{*}$
and $t^{*}=L^{2}/(2\pi^{2}J^{2})$ [7, 13]. The above $L$-dependence of $t^{*}$
is similar to our earlier result in Ref. [7] for a $U(1)$ symmetric fermionic
model without the qubits.

However, there is another interesting parameter regime when
$(1-\cos\frac{2i\pi}{L}\big{)}\gg(g/J)^{2}$ for a finite $L$, which is the
case in many proposed highly controlled laboratory settings to test
predictions for SFF [47, 43, 44, 45, 18]. For
$(1-\cos\frac{2i\pi}{L}\big{)}\gg(g/J)^{2}$ at finite $L~{}(<{\it l}_{c})$, we
have $\lambda_{i}\approx 1-g^{2}+(g^{2}/2J)^{2}\csc^{2}(i\pi/L)$ for
$i=1,2,3\dots,L-1$. Therefore, the second largest eigenvalues for small $g/J$
are approximately $L-1$ fold degenerate. These features of $\lambda_{i}$ give
a different form of SFF, and the system-size scaling of $t^{*}$: $K(t)\simeq
2t(1+\sum_{i=1}^{L-1}\lambda_{i}^{t})$, which leads to
$t^{*}\approx\mathcal{O}(\log L)$ (check Sec. II of SM[46]). Such logarithmic
system size dependence of $t^{*}$ has been previously reported for $U(1)$
symmetry-broken models in the absence of total particle number conservation
[1]. We can also increase the value of $g/J$ at a fixed finite $L$ to access
the other condition, $(1-\cos\frac{2i\pi}{L}\big{)}\ll(g/J)^{2}$, to get the
SFF in Eq. 10, and $t^{*}\propto L^{2}$. Thus, we find a crossover in system
size scaling of $t^{*}$ with a varying scaled mixing strength $g/J$ at finite
lengths in our model with $JC$ mixing between fermions and qubits.

To demonstrate two different $L$ scaling of $t^{*}$, we plot $K(t)$ with $t$
using Eq. 7, which is obtained by applying the RPA and identity permutation
for leading order contributions. In Figs. 1(a,c), we show $K(t)$ with $t$ for
$g=0.1,J=0.4~{}({\it l}_{c}\sim 17)$ and $g=0.4,J=0.1~{}({\it l}_{c}\sim 0)$,
respectively. We take the half-filled case with $N/L=1/2$. We can understand
the $L$ dependence of $t^{*}$ for these two parameter sets by scaling $t$ and
$K(t)$ by predicted $L$ dependence. For this, we plot $K(t)/\log L$ against
$t/\log L$ in Fig. 1(b) and $K(t)/L^{1.85}$ against $t/L^{1.85}$ in Fig. 1(d).
Figs. 1(b,d) display a nice data collapse for different $L$ at a time above
$t^{*}$ for the universal RMT behavior of the SFF. Such data collapse confirms
our above-predicted crossover of the $L$ dependence of $t^{*}$ with an
increasing $g/J$. We could not get $t^{*}$ growing exactly as $L^{2}$ for a
large $g/J$ in our numerics with limited $L$. Still, our obtained exponent
$(\sim 1.85)$ in this region is close to the predicted value of $2$.

The generating Hamiltonian for $JC$ mixing between bosons and qubits in the
Trotter regime reads as (see Sec. III of SM[46] for a derivation)

$\displaystyle\mathcal{M}_{\text{J}C}^{\text{B}}$ $\displaystyle=$
$\displaystyle\big{(}1+\frac{(g^{2}+J^{2})L}{2}\big{)}\mathds{1}_{\mathcal{N}_{\text{J}C}^{\text{B}}}+\sum_{i=1}^{L}\Big{(}2J^{2}\big{(}\hat{K}^{1}_{i}\hat{K}^{1}_{i+1}$
(11) $\displaystyle+$
$\displaystyle\hat{K}^{2}_{i}\hat{K}^{2}_{i+1}-\hat{K}^{0}_{i}\hat{K}^{0}_{i+1}\big{)}+g^{2}(\hat{K}^{1}_{i}\hat{\sigma}^{x}_{i}-\hat{K}^{2}_{i}\hat{\sigma}^{y}_{i}$
$\displaystyle+$
$\displaystyle\hat{K}^{0}_{i}-\hat{\sigma}^{\dagger}_{i}\hat{\sigma}_{i})\Big{)}+\mathcal{O}(J^{4},g^{4}),$

where
$\hat{K}^{1}_{j}=(\hat{K}^{+}_{j}+\hat{K}^{-}_{j})/2,\hat{K}^{2}_{j}=(\hat{K}^{+}_{j}-\hat{K}^{-}_{j})/2i$.
We define a set of local operators
$\hat{K}^{0}_{j}=-(\hat{n}_{j}+1/2),~{}\hat{K}^{+}_{j}=\hat{a}_{j}\sqrt{\hat{n}_{j}},~{}\hat{K}^{-}_{j}=\sqrt{\hat{n}_{j}}\hat{a}^{\dagger}_{j}$,
which satisfy the commutation relations of $SU(1,1)$ algebra at the same site,
and commute otherwise:
$[\hat{K}^{+}_{i},\hat{K}^{-}_{j}]=-2\hat{K}^{0}_{i}\delta_{ij},~{}[\hat{K}^{0}_{i},\hat{K}^{\pm}_{j}]=\pm\hat{K}^{\pm}_{i}\delta_{ij}$.
However, $\mathcal{M}_{\text{J}C}^{\text{B}}$ in Eq. 11 does not commute with
$\hat{K}^{\alpha}=\sum_{i=1}^{L}\hat{K}^{\alpha}_{i}$, $\alpha\in\\{+,-,0\\}$
for a non-zero $g$. Thus, $\mathcal{M}_{\text{J}C}^{\text{B}}$ does not
possess $SU(1,1)$ symmetry unlike the only boson model investigated in Roy
_et al._ [13]. Nevertheless, we find
$[\mathcal{M}_{\text{J}C}^{\text{B}},\sum_{i}(\hat{K}^{0}_{i}+\hat{\sigma}^{\dagger}_{i}\hat{\sigma}_{i})]=0$,
which indicates a $U(1)$ symmetry of $\mathcal{M}_{\text{J}C}^{\text{B}}$. As
shown in SM[46], the $L$-dependence of $t^{*}$ for this model is similar to
that of fermions and qubits. For a finite $L$, there is a crossover in the
$L$-dependence of $t^{*}$ from $\log L$ to $L^{2}$ with an increasing $g/J$
for $JC$ mixing between bosons and qubits. The eigenvalues of
$\mathcal{M}_{\text{J}C}^{\text{B}}$ are identical to those of
$\mathcal{M}_{\text{J}C}^{\text{F}}$ for $N=1$. The largest eigenvalues of
$\mathcal{M}_{\text{J}C}^{\text{B}}$ for any finite $N$ become degenerate with
those for $N=1$ with an increasing $L$ due to an emergent approximate symmetry
of $\mathcal{M}_{\text{J}C}^{\text{B}}$. The above features lead to the
similarity between the fermionic and bosonic models with $JC$ mixing.

Next, we consider $R$ mixing between fermions or bosons and qubits. We start
with the fermionic case having a finite-dimensional Hilbert space. The
generating Hamiltonian in this case is (see Sec. IV of SM[46] for a
derivation)

$\displaystyle\mathcal{M}_{\text{R}}^{\text{F}}$ $\displaystyle=$
$\displaystyle\big{(}1-\frac{(2g^{2}+J^{2})L}{2}\big{)}\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{F}}}+\sum_{i=1}^{L}\big{(}\sum_{\nu}\frac{J^{2}}{2}\hat{\tau}_{i}^{\nu}\hat{\tau}_{i+1}^{\nu}$
(12) $\displaystyle+$ $\displaystyle
g^{2}\hat{\tau}_{i}^{z}\hat{\sigma}_{i}^{z}\big{)}+\mathcal{O}(J^{4},g^{4}),$

which commutes with $\sum_{i=1}^{L}\hat{\tau}_{i}^{z}$ and
$\hat{\sigma}_{j}^{z}$ for $j\in\\{1,2,3,\dots,L\\}$, indicating a global
$U(1)$ symmetry for fermions and local $U(1)$ symmetry for each qubit.
Interestingly, $\hat{H}(t)$ does not have a global $U(1)$ symmetry for $R$
mixing. The generating Hamiltonian for $R$ mixing does not have $SU(2)$
symmetry due to magnetic anisotropy created by coupling to the qubits in
contrast to that in Eq. 8 for $JC$ mixing between fermions and qubits. The
eigenvalues $\lambda_{i}$ of $\mathcal{M}_{\text{R}}^{\text{F}}$ can be
determined by fixing $\sum_{i=1}^{L}\hat{\tau}_{i}^{z}$ and
$\hat{\sigma}_{j}^{z}$ for $j\in\\{1,2,3,\dots,L\\}$ as these are good quantum
numbers. The eigenvalues are doubly degenerate since
$\mathcal{M}_{\text{R}}^{\text{F}}$ is invariant under
$\prod_{i=1}^{L}\hat{\tau}_{i}^{x}\hat{\sigma}_{i}^{x}$, which implies a state
obtained by flipping all the $\tau$ and $\sigma$ spins of an eigenstate of
$\mathcal{M}_{\text{R}}^{\text{F}}$ is also an eigenstate with the same
eigenvalue. The largest eigenvalue one of $\mathcal{M}_{\text{R}}^{\text{F}}$
is a state $|\lambda_{0}\rangle$ in which all $\tau$ and $\sigma$ spins are
polarized in $z$ direction. The second largest eigenvalues of
$\mathcal{M}_{\text{R}}^{\text{F}}$ are $L+1$ and $L$ fold degenerate,
respectively, for $(g/J)^{2}<2/3$ and $(g/J)^{2}>2/3$ (see SM[46] for
details). For $(g/J)^{2}<2/3$, the second largest eigenvalues are $1-2g^{2}$,
which consist of $L$ eigenstates with anyone $\sigma$ spin being flipped in
$|\lambda_{0}\rangle$ and another superposition state with a single $\tau$
spin flipping in $|\lambda_{0}\rangle$. For $(g/J)^{2}>2/3$, the second
largest eigenvalues are $1-4g^{2}-2J^{2}(1-\sqrt{1+4(g/J)^{4}})$, which are
$L$ eigenstates with one $\tau$ spin flipping and one $\sigma$ spin being
flipped in $|\lambda_{0}\rangle$. Thus, the second largest eigenvalues for any
$g/J$ are $L$ independent. So we get $t^{*}\propto\log L$ or $\log(L+1)$ for
$R$ mixing between fermions and qubits. Such $L$-dependence of $t^{*}$ is
similar to that in a periodically kicked transverse-field Ising model in Kos
_et al._ [1] with local kicking terms. Interestingly, a similar $L$-scaling of
$t^{*}$ can also be obtained for the $U(1)$-symmetry broken model explored in
Ref. [7] when the pairing $\Delta$ and tunneling $J$ strengths are the same.
We can also get $t^{*}\propto L^{0}$ of Roy and Prosen [7] for arbitrary
$\Delta$ and $J$ when $g$ is different/random for different qubits to lift the
degeneracy in the second largest eigenvalues.

Finally, we consider the mixture of bosons and qubits with $R$ mixing between
them. The generating Hamiltonian for this case in the Trotter regime reads as
(see Sec. V of SM[46] for a derivation)

$\displaystyle\mathcal{M}_{\text{R}}^{\text{B}}$ $\displaystyle=$
$\displaystyle\big{(}1+\frac{J^{2}L}{2}\big{)}\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{B}}}+\sum_{i=1}^{L}\Big{(}2J^{2}\big{(}\hat{K}^{1}_{i}\hat{K}^{1}_{i+1}+\hat{K}^{2}_{i}\hat{K}^{2}_{i+1}$
(13) $\displaystyle-$
$\displaystyle\hat{K}^{0}_{i}\hat{K}^{0}_{i+1}\big{)}+2g^{2}(\hat{K}^{1}_{i}\hat{\sigma}^{x}_{i}+\hat{K}^{0}_{i})\Big{)}+\mathcal{O}(J^{4},g^{4}),$

which commutes with $\hat{\sigma}_{j}^{x}$ for $j\in\\{1,2,3,\dots,L\\}$. We
could not calculate the spectrum of $\mathcal{M}_{\text{R}}^{\text{B}}$
analytically. Instead, we determine it numerically by varying $N_{max}$ for a
fixed $L$ to get an estimate of $\lambda_{1}$ in the large $N_{max}$ limit. We
use linear extrapolations in $1/N_{max}$ towards $1/N_{max}=0$ to find
evidence for a gap between the largest and second largest eigenvalues, as
shown in Sec. V of SM[46]. The second largest eigenvalues are also $L$-fold
degenerate, suggesting a scaling of $t^{*}\propto\log L$. We remind here that
a periodically kicked boson model without particle number conservation shows
$t^{*}\propto\mathcal{O}(L^{\gamma}),\gamma=0.7\pm 0.1$ [13], which is sharply
different from the present case of bosons and qubits without total number
conservation.

We have analytically calculated the SFF in many-body quantum mixtures of
fermions, bosons and qubits with periodically kicked spreading and mixing of
species. Different types of mixing between species can drastically alter the
timescale for the emergence of RMT behavior of $K(t)$ in quantum mixtures. We
show how competition between mixing and hopping/spreading of species in
$U(1)$-symmetric finite-size systems can lead to a logarithmic $L$ scaling of
$t^{*}$, which has been predicted before only for $U(1)$-symmetry broken
single-species models [1, 4]. This finding is practical and vital as quantum
mixtures of different species are abundant in nature as well as controlled
experimental set-ups of cold atoms and photonic systems, and many of these
systems are finite-sized. We further show the $t^{*}$ scaling for $R$ mixing
of fermions and qubits is similar to those obtained for a single species of
spin-1/2s or fermions. Finally, our results indicate that the $R$ mixing of
species with different statistics (e.g., bosons and qubits) can lead to
completely new features for the main species (e.g., bosons) with individual
hopping.

Acknowledgements: We thank Prof. Toma$\check{\rm z}$ Prosen for many useful
discussions.

## References

  * Kos  _et al._ [2018] P. Kos, M. Ljubotina, and T. Prosen, Many-body quantum chaos: Analytic connection to random matrix theory, Phys. Rev. X 8, 021062 (2018).
  * Bertini  _et al._ [2018] B. Bertini, P. Kos, and T. Prosen, Exact spectral form factor in a minimal model of many-body quantum chaos, Phys. Rev. Lett. 121, 264101 (2018).
  * Chan  _et al._ [2018a] A. Chan, A. De Luca, and J. T. Chalker, Solution of a minimal model for many-body quantum chaos, Phys. Rev. X 8, 041019 (2018a).
  * Chan  _et al._ [2018b] A. Chan, A. De Luca, and J. T. Chalker, Spectral statistics in spatially extended chaotic quantum many-body systems, Phys. Rev. Lett. 121, 060601 (2018b).
  * Friedman  _et al._ [2019] A. J. Friedman, A. Chan, A. De Luca, and J. T. Chalker, Spectral statistics and many-body quantum chaos with conserved charge, Phys. Rev. Lett. 123, 210603 (2019).
  * Bertini  _et al._ [2019] B. Bertini, P. Kos, and T. Prosen, Exact correlation functions for dual-unitary lattice models in $1+1$ dimensions, Phys. Rev. Lett. 123, 210601 (2019).
  * Roy and Prosen [2020] D. Roy and T. Prosen, Random matrix spectral form factor in kicked interacting fermionic chains, Phys. Rev. E 102, 060202(R) (2020).
  * Bertini  _et al._ [2021] B. Bertini, P. Kos, and T. Prosen, Random matrix spectral form factor of dual-unitary quantum circuits, Commun. Math. Phys. 387, 597 (2021).
  * Garratt and Chalker [2021] S. J. Garratt and J. T. Chalker, Local pairing of feynman histories in many-body floquet models, Phys. Rev. X 11, 021051 (2021).
  * Li  _et al._ [2021] J. Li, T. Prosen, and A. Chan, Spectral statistics of non-Hermitian matrices and dissipative quantum chaos, Phys. Rev. Lett. 127, 170602 (2021).
  * Kos  _et al._ [2021] P. Kos, B. Bertini, and T. Prosen, Chaos and ergodicity in extended quantum systems with noisy driving, Phys. Rev. Lett. 126, 190601 (2021).
  * Moudgalya  _et al._ [2021] S. Moudgalya, A. Prem, D. A. Huse, and A. Chan, Spectral statistics in constrained many-body quantum chaotic systems, Phys. Rev. Research 3, 023176 (2021).
  * Roy  _et al._ [2022] D. Roy, D. Mishra, and T. Prosen, Spectral form factor in a minimal bosonic model of many-body quantum chaos, Phys. Rev. E 106, 024208 (2022).
  * Winer and Swingle [2022a] M. Winer and B. Swingle, The Loschmidt spectral form factor, J. High Energ. Phys. 2022, 137 (2022).
  * Winer and Swingle [2022b] M. Winer and B. Swingle, Hydrodynamic theory of the connected spectral form factor, Phys. Rev. X 12, 021009 (2022b).
  * Liao and Galitski [2022a] Y. Liao and V. Galitski, Emergence of many-body quantum chaos via spontaneous breaking of unitarity, Phys. Rev. B 105, L140202 (2022a).
  * Liao and Galitski [2022b] Y. Liao and V. Galitski, Emergence of many-body quantum chaos via spontaneous breaking of unitarity, Phys. Rev. B 105, L140202 (2022b).
  * Dag  _et al._ [2023] C. Dag, S. Mistakidis, A. Chan, and H. R. Sadeghpour, Many-body quantum chaos in stroboscopically-driven cold atoms, Commun. Phys. 6, 136 (2023).
  * Haake [2001] F. Haake,  _Quantum Signatures of Chaos, 2nd ed._ (Springer, New York, 2001).
  * Berry [1985] M. V. Berry, Semiclassical theory of spectral rigidity, Proc. R. Soc. Lond. A 400, 229 (1985).
  * Bohigas  _et al._ [1984] O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. Lett. 52, 1 (1984).
  * McDonald and Kaufman [1979] S. W. McDonald and A. N. Kaufman, Spectrum and eigenfunctions for a Hamiltonian with stochastic trajectories, Phys. Rev. Lett. 42, 1189 (1979).
  * Casati  _et al._ [1980] G. Casati, F. Valz-Gris, and I. Guarnieri, On the connection between quantization of nonintegrable systems and statistical theory of spectra, Lett. Nuovo Cimento 28, 279 (1980).
  * Berry and Tabor [1977] M. V. Berry and M. Tabor, Level clustering in the regular spectrum, Proc. R. Soc. Lond. A 356, 375 (1977).
  * Berry [1981] M. V. Berry, Quantizing a classically ergodic system: Sinai’s billiard and the KKR method, Ann. Phys. 131, 163 (1981).
  * Sieber and Richter [2001] M. Sieber and K. Richter, Correlations between periodic orbits and their role in spectral statistics, Phys. Scr. T90, 128 (2001).
  * Sieber [2002] M. Sieber, Leading off-diagonal approximation for the spectral form factor for uniformly hyperbolic systems, J. Phys. A 35, L613 (2002).
  * Müller  _et al._ [2004] S. Müller, S. Heusler, P. Braun, F. Haake, and A. Altland, Semiclassical foundation of universality in quantum chaos, Phys. Rev. Lett. 93, 014103 (2004).
  * Müller  _et al._ [2005] S. Müller, S. Heusler, P. Braun, F. Haake, and A. Altland, Periodic-orbit theory of universality in quantum chaos, Phys. Rev. E 72, 046207 (2005).
  * Chávez-Carlos  _et al._ [2019] J. Chávez-Carlos, B. López-del Carpio, M. A. Bastarrachea-Magnani, P. Stránský, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Quantum and classical Lyapunov exponents in atom-field interaction systems, Phys. Rev. Lett. 122, 024101 (2019).
  * Jaynes and Cummings [1963] E. Jaynes and F. Cummings, Chaos and ergodicity in extended quantum systems with noisy driving, Proc. IEEE. 51, 89 (1963).
  * Scully and Zubairy [1997] M. Scully and M. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Hartmann  _et al._ [2006] M. J. Hartmann, F. G. S. L. Brandão, and M. B. Plenio, Strongly interacting polaritons in coupled arrays of cavities, Nature Physics 2, 849 (2006).
  * Angelakis  _et al._ [2007] D. G. Angelakis, M. F. Santos, and S. Bose, Photon-blockade-induced Mott transitions and $XY$ spin models in coupled cavity arrays, Phys. Rev. A 76, 031805 (2007).
  * Roy  _et al._ [2017] D. Roy, C. M. Wilson, and O. Firstenberg, Colloquium: Strongly interacting photons in one-dimensional continuum, Rev. Mod. Phys. 89, 021001 (2017).
  * Rabi [1937] I. I. Rabi, Space quantization in a gyrating magnetic field, Phys. Rev. 51, 652 (1937).
  * Xie  _et al._ [2017] Q. Xie, H. Zhong, M. T. Batchelor, and C. Lee, The quantum Rabi model: solution and dynamics, J. Phys. A 50, 113001 (2017).
  * Agarwal  _et al._ [2023] L. Agarwal, S. Sahu, and S. Xu, Charge transport, information scrambling and quantum operator-coherence in a many-body system with U(1) symmetry, J. High Energ. Phys. 2023, 37 (2023).
  * Gharibyan  _et al._ [2018] H. Gharibyan, M. Hanada, S. H. Shenker, and M. Tezuka, Onset of random matrix behavior in scrambling systems, J. High Energ. Phys. 2018, 124 (2018).
  * Islam  _et al._ [2015] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Cronenberger  _et al._ [2019] S. Cronenberger, C. Abbas, D. Scalbert, and H. Boukari, Spatiotemporal spin noise spectroscopy, Phys. Rev. Lett. 123, 017401 (2019).
  * Swar  _et al._ [2021] M. Swar, D. Roy, S. Bhar, S. Roy, and S. Chaudhuri, Detection of spin coherence in cold atoms via Faraday rotation fluctuations, Phys. Rev. Res. 3, 043171 (2021).
  * Joshi  _et al._ [2022] L. K. Joshi, A. Elben, A. Vikram, B. Vermersch, V. Galitski, and P. Zoller, Probing many-body quantum chaos with quantum simulators, Phys. Rev. X 12, 011018 (2022).
  * [46] See Supplemental Materials (SM) for the derivation of spectral form factor and generating Hamiltonians for different mixing, and the details of spectral properties of these Hamiltonians in the Trotter regime.
  * Roy  _et al._ [2015] D. Roy, R. Singh, and R. Moessner, Probing many-body localization by spin noise spectroscopy, Phys. Rev. B 92, 180205(R) (2015).
  * Economou [2006] E. N. Economou,  _Green’s functions in quantum physics_ , Vol. 7 (Springer Science & Business Media, 2006).

Supplementary Material for “Many-Body Quantum Chaos in Mixtures of Multiple
Species”

Vijay Kumar and Dibyendu Roy

Raman Research Institute, Bangalore 560080, India

##  I Spectral Form Factor for Periodically Kicked Systems

The spectral form factor (SFF) is defined as the Fourier transform of the two-
point correlation function of spectral density. For a periodically kicked
system, the spectral density is,
$\rho(\varphi)=\left(2\pi/\mathcal{N}\right)\sum_{n=1}^{\mathcal{N}}\delta(\varphi-\varphi_{n})$,
where $\varphi_{n}$’s are the eigenphases of the Floquet propagator,
$\hat{U}$, and $\mathcal{N}$ is the dimension of the system’s Hilbert space.
The prefactor, $2\pi/\mathcal{N}$, is chosen such that the averaged spectral
density is normalized to unity, i.e.,

$\displaystyle\langle\rho(\varphi)\rangle_{\varphi}=\int_{0}^{2\pi}d\varphi\>\rho(\varphi)=1.$
(S1)

We define a two-point correlation function of $\rho(\varphi)$ as

$\displaystyle
R(\vartheta)=\langle\rho(\varphi+\vartheta/2)\rho(\varphi-\vartheta/2)\rangle_{\varphi}-\langle\rho(\varphi)\rangle_{\varphi}^{2}.$
(S2)

We can then write the SFF as

$\displaystyle K(t)=\frac{\mathcal{N}^{2}}{2\pi}\int_{0}^{2\pi}d\vartheta
R(\vartheta)e^{-i\vartheta
t}=(\text{tr}\hat{U}^{t})(\text{tr}\hat{U}^{-t})-\mathcal{N}^{2}\delta_{t,0}.$
(S3)

Since $K(t)$ in Eq. S3 is not self-averaging over disorder in onsite energy
and transition frequency, we average it over different disorder realizations
of them.

$\displaystyle
K(t)=\langle(\text{tr}\hat{U}^{t})(\text{tr}\hat{U}^{-t})\rangle-\mathcal{N}^{2}\delta_{t,0},$
(S4)

where $\langle...\rangle$ represent averaging over disorder realizations. In
the main paper, we consider periodically kicked quantum mixtures of multiple
species whose Hamiltonian reads

$\displaystyle\hat{H}(t)$
$\displaystyle=\hat{H}_{0}+\hat{H}_{\text{JC/R}}\sum_{m\in\mathcal{Z}}\delta(t-m),$
(S5)

where $\hat{H}_{0}$ is the base Hamiltonian and $\hat{H}_{\text{JC/R}}$ is the
driving Hamiltonian. These Hamiltonians are given in the main paper. The
Floquet propagator over the unit cycle of kicking is

$\displaystyle\hat{U}$
$\displaystyle=\mathcal{T}e^{-i\int_{0}^{1}dt\hat{H}(t)}=\lim_{\epsilon\rightarrow
0}e^{-i\int_{\epsilon}^{1+\epsilon}dt\hat{H}(t)}=\hat{V}\hat{W},$ (S6)

where $\hat{W}=e^{-i\hat{H}_{0}},\hat{V}=e^{-i\hat{H}_{\text{JC/R}}}$. To
proceed with the calculation of the SFF in Eq. S4, we choose the occupation
number basis,
$|\underline{n\sigma}\rangle\equiv|n_{1},...,n_{L}\rangle\otimes|\sigma_{1},...,\sigma_{L}\rangle$,
which are the eigenbasis of $H_{0}$ [7, 13]. Here, $n_{j},\sigma_{j}$ are
respectively the numbers of fermions/bosons and qubit excitation at site $j$.
Thus, we have

$\displaystyle\hat{W}|\underline{n\sigma}\rangle$
$\displaystyle=e^{-i\theta_{\underline{n\sigma}}}|\underline{n\sigma}\rangle,$
(S7) $\displaystyle\theta_{\underline{n\sigma}}$
$\displaystyle=\sum_{i=1}^{L}(\omega_{i}n_{i}+\Omega_{i}\sigma_{i})+\sum_{i<j}U_{ij}n_{i}n_{j}.$
(S8)

We derive $\text{tr}\hat{U}^{t}$ by inserting the identity operator
$\sum_{\underline{n\sigma}_{\mu}}|\underline{n\sigma}_{\mu}\rangle\langle\underline{n\sigma}_{\mu}|=\mathds{1}_{\mathcal{N}}$,
at different time steps $\mu=1,2,...,t$:

$\displaystyle\text{tr}\hat{U}^{t}$
$\displaystyle=\sum_{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}}\langle\underline{n\sigma}_{1}|\hat{V}\hat{W}|\underline{n\sigma}_{2}\rangle\langle\underline{n\sigma}_{2}|\hat{V}\hat{W}...|\underline{n\sigma}_{t}\rangle$
$\displaystyle\qquad\qquad\qquad\qquad\times\langle\underline{n\sigma}_{t}|\hat{V}\hat{W}|\underline{n\sigma}_{1}\rangle,$
(S9)

where trace requires periodic boundary condition (PBC) in time,
$\underline{n\sigma}_{t+1}=\underline{n\sigma}_{1}$. Using Eq. S7 in Eq. I, we
get for the SFF:

$\displaystyle K(t)$
$\displaystyle=\sum_{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}}\sum_{\underline{n\sigma}_{1}^{\prime},...,\underline{n\sigma}_{t}^{\prime}}\langle
e^{-i\sum_{\mu=1}^{t}\left(\theta_{\underline{n\sigma}_{\mu}}-\theta_{\underline{n\sigma}_{\mu}^{\prime}}\right)}\rangle$
$\displaystyle\qquad\times\prod_{\mu=1}^{t}V_{\underline{n\sigma}_{\mu},\underline{n\sigma}_{\mu+1}}V^{*}_{\underline{n\sigma}_{\mu}^{\prime},\underline{n\sigma}_{\mu+1}^{\prime}}.$
(S10)

For the randomness in the parameters of $\hat{H}_{0}$ and the long range
interaction, we take $\theta_{\underline{n\sigma}_{\mu}}$(mod $2\pi$) as
uniform iid’s over $[0,2\pi)$. The last consideration is a random phase
approximation (RPA). Within the RPA, we approximate $\langle
e^{-i\sum_{\mu=1}^{t}\left(\theta_{\underline{n\sigma}_{\mu}}-\theta_{\underline{n\sigma}_{\mu}^{\prime}}\right)}\rangle\simeq\delta_{\\{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}\\},\\{\underline{n\sigma}_{1}^{\prime},...,\underline{n\sigma}_{t}^{\prime}\\}}$,
where

$\displaystyle\delta_{\\{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}\\},\\{\underline{n\sigma}_{1}^{\prime},...,\underline{n\sigma}_{t}^{\prime}\\}}$
$\displaystyle=\begin{cases}1,\text{ if
$\\{\underline{n\sigma}_{1}^{\prime},...,\underline{n\sigma}_{t}^{\prime}\\}$
is a}\\\ \;\;\;\;\;\text{permutation of}\\\
\;\;\;\;\;\text{$\\{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}\\}$}\\\
0,\text{ otherwise}\end{cases}.$ (S11)

For $t\ll t_{\text{H}}$ (where $t_{\text{H}}$ is the Heisenberg time) when a
fraction ($\approx
t^{2}/t_{\text{H}}+\mathcal{O}\left(1/t_{\text{H}}^{2}\right)$) of
configurations in $\\{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}\\}$
with repeated basis vectors is negligible, Kos  _et al._ [1] have also shown
that such configurations with repeated basis do not contribute to leading
order term in the SFF. We also perform exact numerical computations of the SFF
using Eq. S4 to compare it to that obtained using RPA. We show these
comparisons in Figs. S2,S5 to validate the RPA analysis for certain parameter
regimes of $\hat{H}_{0}$ in finite systems. Therefore, we can safely assume
that the relevant permutations are, $\pi\in S_{t}$, where, $S_{t}$, is the
permutation group of $t$ distinct objects. Thus, the Eq. I reduces to

$\displaystyle K(t)$
$\displaystyle=\sum_{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}}\sum_{\pi\in
S_{t}}\prod_{\mu=1}^{t}V_{\underline{n\sigma}_{\mu},\underline{n\sigma}_{\mu+1}}V^{*}_{\underline{n\sigma}_{\pi(\mu)},\underline{n\sigma}_{\pi(\mu+1)}}.$
(S12)

Considering $t$ cyclic and $t$ anti-cyclic permutations, the SFF can be
written as

$\displaystyle
K(t)=2t\sum_{\underline{n\sigma}_{1},...,\underline{n\sigma}_{t}}\prod_{\mu=1}^{t}|V_{\underline{n\sigma}_{\mu},\underline{n\sigma}_{\mu+1}}|^{2}=2t\text{
tr}\mathcal{M}^{t},$ (S13)

where, the elements of $\mathcal{M}$ are
$\mathcal{M}_{\underline{n\sigma},\underline{n\sigma}^{\prime}}=|\langle\underline{n\sigma}|\hat{V}|\underline{n\sigma}^{\prime}\rangle|^{2}$.
Here, $\mathcal{M}$ is a doubly stochastic matrix due to unitarity of
$\hat{V}$, therefore, the largest eigenvalue of $\mathcal{M}$ is 1 and other
eigenvalues ($\lambda_{1},...,\lambda_{\mathcal{N}-1}$) satisfy,
$1\geq|\lambda_{1}|\geq|\lambda_{2}|\geq...\geq|\lambda_{\mathcal{N}-1}|$.
Now, the SFF in Eq. S13 in terms of eigenvalues of $\mathcal{M}$ reads

$\displaystyle
K(t)=2t\left(1+\sum_{i=1}^{\mathcal{N}-1}\lambda_{i}^{t}\right),$ (S14)

which is given in the main paper. For circular orthogonal ensemble (COE), the
leading order behavior of universal RMT form $2t$, which appears at long time
beyond the Thouless timescale $t^{*}(L)$ when the second term,
$\sum_{i=1}^{\mathcal{N}-1}\lambda_{i}^{t}$, diminishes as $|\lambda_{i}|<1$.
Since, $\lambda_{i}^{t}$ vanishes faster for smaller value of $|\lambda_{i}|$,
the nature of Thouless time is mainly determined by the largest eigenvalues.
To find Thouless time, $t^{*}(L)$, we study the eigenspectrum of $\mathcal{M}$
in the Trotter regime of small parameters of the driving Hamiltonian. We write
$\mathcal{M}$ in terms of the Hadamard product of $\hat{H}_{\text{JC/R}}$ by
Taylor expanding $\hat{V}$:

$\displaystyle\mathcal{M}=\sum_{k=0}^{\infty}\sum_{l=-k}^{k}(-1)^{l}\frac{\left(\hat{H}_{\text{JC/R}}\right)^{k+l}\bullet\left(\hat{H}_{\text{JC/R}}\right)^{k-l}}{(k+l)!(k-l)!}.$
(S15)

We approximate $\mathcal{M}$ by keeping the terms in Eq. S15 upto two leading
orders ($k=0,1$) in the Trotter regime:

$\displaystyle\mathcal{M}=\mathds{1}+\hat{H}_{\text{JC/R}}\bullet\hat{H}_{\text{JC/R}}-(\hat{H}_{\text{JC/R}})^{2}\bullet\mathds{1}+\mathcal{O}(\hat{H}_{\text{JC/R}}^{4}).$
(S16)

We explicitly calculate $\mathcal{M}$ in the Trotter regime for different
types of mixing between fermions or bosons and qubits investigated in the main
paper.

##  II Jaynes-Cummings mixing between fermions and qubits

The driving Hamiltonian $\hat{H}_{\text{JC}}$ for Jaynes-Cummings $(JC)$
mixing between fermions and qubits is given in Eq. (2) of the main paper.

$\displaystyle\hat{H}_{\text{JC}}=\sum_{i=1}^{L}\left(g\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}-J\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+{\rm
H.c.}\right).$ (S17)

We rewrite the fermion creation and annihilation operators in terms of
spin-1/2 operators using the Jordan-Wigner transformation.

$\displaystyle\hat{a}_{i}=\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right)\hat{\tau}_{i}\>,~{}~{}\hat{a}_{i}^{\dagger}=\hat{\tau}_{i}^{\dagger}\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right),$
(S18)

where $\hat{\tau}_{j}=\left(\hat{\tau}_{j}^{x}-i\hat{\tau}_{j}^{y}\right)/2$
and
$\hat{\tau}_{j}^{\dagger}=\left(\hat{\tau}_{j}^{x}+i\hat{\tau}_{j}^{y}\right)/2$
are the spin-1/2 lowering and raising operators, respectively. Here,
$\hat{\tau}_{j}^{x,y,z}$ are the Pauli matrices at $j$th site. Substituting
Eq. S18 in Eq. S17, we get with PBC:

$\displaystyle\hat{H}_{\text{JC}}=g\sum_{i=1}^{L}\left(\hat{\tau}_{i}^{\dagger}\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right)\hat{\sigma}_{i}+\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right)\hat{\tau}_{i}\hat{\sigma}_{i}^{\dagger}\right)$
$\displaystyle-J\sum_{i=1}^{L-1}\left(\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i+1}+\hat{\tau}_{i+1}^{\dagger}\hat{\tau}_{i}\right)+J(-1)^{\hat{N}_{f}}\left(\hat{\tau}_{L}^{\dagger}\hat{\tau}_{1}+\hat{\tau}_{1}^{\dagger}\hat{\tau}_{L}\right),$
(S19)

where,
$\hat{N}_{\text{f}}=\sum_{i=1}^{L}\hat{a}_{i}^{\dagger}\hat{a}_{i}=\sum_{i=1}^{L}\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i}$,
is the total number of fermions or corresponding qubit excitations. We notice
that the chosen occupation number basis are equivalent to the eigenbasis of
$\prod_{i=1}^{L}\hat{\tau}_{i}^{z}\hat{\sigma}_{i}^{z}$. To calculate
$\hat{H}_{\text{JC}}\bullet\hat{H}_{\text{JC}}$, we need to find the matrix
elements of $\hat{H}_{\text{JC}}$. In the chosen basis,
$\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right)|\underline{n\sigma}\rangle=\pm|\underline{n\sigma}\rangle$,
and the operators
$\hat{\tau}_{i},\hat{\tau}_{i}^{\dagger},\hat{\sigma}_{i},\hat{\sigma}_{i}^{\dagger}$
flip the related spin states. Thus, the nonzero matrix elements of
$\hat{\tau}_{i}^{\dagger}\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right)\hat{\sigma}_{i}$
are $\pm 1$. Similarly, the nonzero matrix elements of
$\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i+1}$ are 1. Therefore,

$\displaystyle(\hat{\tau}_{i}^{\dagger}\prod_{j=1}^{i-1}\left(-\hat{\tau}_{j}^{z}\right)\hat{\sigma}_{i})\bullet(\hat{\tau}_{i^{\prime}}^{\dagger}\prod_{j=1}^{i^{\prime}-1}\left(-\hat{\tau}_{j}^{z}\right)\hat{\sigma}_{i^{\prime}})=\hat{\tau}_{i}^{\dagger}\hat{\sigma}_{i}\delta_{ii^{\prime}},$
(S20)
$\displaystyle(\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i+1})\bullet(\hat{\tau}_{i^{\prime}}^{\dagger}\hat{\tau}_{i^{\prime}+1})=\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i+1}\delta_{ii^{\prime}}.$
(S21)

The Eqs. S20,S21 are true in the chosen basis but not in an arbitrary basis.
This implies that

$\displaystyle\hat{H}_{\text{JC}}\bullet\hat{H}_{\text{JC}}=\sum_{i=1}^{L}(g^{2}(\hat{\tau}_{i}^{\dagger}\hat{\sigma}_{i}+\hat{\tau}_{i}\hat{\sigma}_{i}^{\dagger})+J^{2}(\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i+1}+\hat{\tau}_{i}\hat{\tau}_{i+1}^{\dagger}))$
$\displaystyle=\sum_{i=1}^{L}\left(\frac{g^{2}}{2}(\hat{\tau}_{i}^{x}\hat{\sigma}_{i}^{x}+\hat{\tau}_{i}^{y}\hat{\sigma}_{i}^{y})+\frac{J^{2}}{2}(\hat{\tau}_{i}^{x}\hat{\tau}_{i+1}^{x}+\hat{\tau}_{i}^{y}\hat{\tau}_{i+1}^{y})\right).$
(S22)

(a)  (b)

Figure S1: The largest eigenvalues $\lambda_{i}$ of
$\mathcal{M}_{\text{JC}}^{\text{F}}$ with $i\in\\{0,1,2,...,49\\}$ for
$g=0.1,J=0.4$ in $(a)$, and $i\in\\{0,1,2,...,20\\}$ for $g=0.4,J=0.1$ in
$(b)$. The plots $(a)$ and $(b)$ show respectively $L$ and $3$ largest
$\lambda_{i}$ being the same for different $N$.

Since $(\hat{H}_{\text{JC}})^{2}\bullet\mathds{1}$ is a matrix of diagonal
elements of $(\hat{H}_{\text{JC}})^{2}$, we find

$\displaystyle(\hat{H}_{\text{JC}})^{2}\bullet\mathds{1}=\sum_{i=1}^{L}\left(g^{2}\frac{1-\hat{\tau}_{i}^{z}\hat{\sigma}_{i}^{z}}{2}+J^{2}\frac{1-\hat{\tau}_{i}^{z}\hat{\tau}_{i+1}^{z}}{2}\right).$
(S23)

Plugging Eqs. II,S23 in Eq. S16, we get

$\displaystyle\mathcal{M}_{\text{JC}}^{\text{F}}=$
$\displaystyle\left(1-\frac{(g^{2}+J^{2})L}{2}\right)\mathds{1}_{\mathcal{N}_{\text{JC}}^{\text{F}}}$
$\displaystyle+\sum_{i=1}^{L}\sum_{\nu}\left(\frac{J^{2}}{2}\hat{\tau}_{i}^{\nu}\hat{\tau}_{i+1}^{\nu}+\frac{g^{2}}{2}\hat{\tau}_{i}^{\nu}\hat{\sigma}_{i}^{\nu}\right)+\mathcal{O}\left(J^{4},g^{4}\right),$
(S24)

where $\nu=x,y,z$ and
$\hat{\sigma}_{j}^{x}=\hat{\sigma}_{j}^{\dagger}+\hat{\sigma}_{j},\hat{\sigma}_{j}^{y}=\left(\hat{\sigma}_{j}^{\dagger}-\hat{\sigma}_{j}\right)/i$.

The generating Hamiltonian $\mathcal{M}_{\text{JC}}^{\text{F}}$ commutes with
the operators
$\sum_{i}\left(\hat{\tau}_{i}^{\nu}+\hat{\sigma}_{i}^{\nu}\right)/2$ for
$\nu=x,y,z$. These operators satisfy $SU(2)$ algebra, which suggests a $SU(2)$
symmetry of $\mathcal{M}_{\text{JC}}^{\text{F}}$. The operator for total
number of excitations in the present model is
$\hat{N}=\sum_{i=1}^{L}\left(\hat{a}_{i}^{\dagger}\hat{a}_{i}+\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i}\right)=\sum_{i=1}^{L}\left(\hat{\tau}_{i}^{\dagger}\hat{\tau}_{i}+\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i}\right)$.
The total number of excitations in a state can be changed by the action of
operators,
$\hat{S}^{+}=\sum_{j}\left(\hat{\tau}_{j}^{\dagger}+\hat{\sigma}_{j}^{\dagger}\right),\hat{S}^{-}=\sum_{j}\left(\hat{\tau}_{j}+\hat{\sigma}_{j}\right)$.
Since, $[\hat{N},\hat{S}^{\pm}]=\pm\hat{S}^{\pm}$, the application of
$\hat{S}^{\pm}$ once on a state leads to a change of the total number of
excitations by $\pm 1$, respectively. We notice $\hat{S}^{\pm}$ commutes with
$\mathcal{M}_{\text{JC}}^{\text{F}}$. Thus, if $|\psi\rangle$ is an eigenstate
of $\mathcal{M}_{\text{JC}}^{\text{F}}$ with an eigenvalue $\lambda$ and total
excitation $N=1$, $\hat{S}^{+}|\psi\rangle$ is another eigenstate with the
same eigenvalue but $N=2$. We can then construct eigenstates of
$\mathcal{M}_{\text{JC}}^{\text{F}}$ with a higher number of total excitations
but with the same eigenvalue $\lambda$ by repeating the applications of
$\hat{S}^{+}$.

We find from the numerics in Fig. S1 that the second largest eigenvalue
$\lambda_{1}$ of $\mathcal{M}_{\text{JC}}^{\text{F}}$ is the same for
$N=1,...,2L-1$. We further notice that $L$ largest eigenvalues of
$\mathcal{M}_{\text{JC}}^{\text{F}}$ are the same for any $N$ in finite-length
chains when $g=0.1,J=0.4$. The $L$ largest eigenvalues can be computed
analytically for $N=1$ sector for any value of $g,J$. They are

$\displaystyle\lambda_{i}$ $\displaystyle=1-g^{2}-J^{2}\left(1-\cos\frac{2\pi
i}{L}\right)$ $\displaystyle\quad+\sqrt{J^{4}\left(1-\cos\frac{2\pi
i}{L}\right)^{2}+g^{4}},$ (S25)

where $i=0,1,...,L-1$. We have $\lambda_{0}=1$ for $i=0$ as expected for a
double stochastic square matrix. The eigenvalues
$\lambda_{1},\lambda_{2},\dots\lambda_{L-1}$ in Eq. II are approximately
degenerate when $g/J\ll 1$ as shown in Fig. S1(a). We try to find how the
Thouless time scales with system size $L$ using these eigenvalues. We
approximate the SFF in Eq. S14 as

$\displaystyle K(t)=2t\left(1+\sum_{i=1}^{L-1}\lambda_{i}^{t}\right).$ (S26)

We can approximate $\lambda_{i}$ in Eq. II for $g/J\ll 1$ and finite $L$ as

$\displaystyle\lambda_{i}\approx
1-g^{2}+\frac{g^{4}}{2J^{2}\left(1-\cos\frac{2\pi i}{L}\right)}.$ (S27)

Plugging $\lambda_{i}$ from Eq. S27 in Eq. S26, the SFF reads

$\displaystyle K(t)$
$\displaystyle=2t\left(1+\sum_{i=1}^{L-1}\left(1-g^{2}+\frac{g^{4}}{2J^{2}\left(1-\cos\frac{2\pi
i}{L}\right)}\right)^{t}\right).$ (S28)

We can further simplify the second part of SFF in Eq. S28 since
$1-g^{2}\gg\frac{g^{4}}{2J^{2}\left(1-\cos\frac{2\pi i}{L}\right)}$ for
$g/J\ll 1$ and finite $L$. We thus write

$\displaystyle\sum_{i=1}^{L-1}\left(1-g^{2}+\frac{g^{4}}{2J^{2}\left(1-\cos\frac{2\pi
i}{L}\right)}\right)^{t}$
$\displaystyle=(1-g^{2})^{t}(L-1)+\frac{tg^{4}(1-g^{2})^{t-1}}{4J^{2}}\sum_{i=1}^{L-1}\frac{1}{\sin^{2}\left(\frac{\pi
i}{L}\right)}$
$\displaystyle=(1-g^{2})^{t}(L-1)+\frac{tg^{4}(1-g^{2})^{t-1}}{4J^{2}}\frac{L^{2}-1}{3},$
(S29)

where we apply the identity $\sum_{i=1}^{L-1}\csc^{2}\left(\frac{\pi
i}{L}\right)=\frac{L^{2}-1}{3}$ in the last line. The Thouless time $t^{*}$ is
defined as the time $t$ when the contribution in Eq. II becomes order of one.
Thus, we have

$\displaystyle(1-g^{2})^{t^{*}-1}\left((L-1)(1-g^{2})+\frac{t^{*}g^{4}}{4J^{2}}\frac{L^{2}-1}{3}\right)\approx
1.$ (S30)

From the first part of left hand side of Eq. S30, we find
$t^{*}(L)\approx-\log(L-1)/\log(1-g^{2})$. To get a better estimate of
$L$-dependence of $t^{*}(L)$ including the second part in Eq. S30, we
numerically solve the above equation to find $t^{*}$ for different
$L\in[6,12]$ when $g=0.1,J=0.4$. Then we fit the data to find $L$-dependence
of $t^{*}(L)$ to get $t^{*}(L)=129.24116\log(L+0.10467)-67.74830$, which
implies that $t^{*}(L)\sim\log(L+0.10467)-0.52420$ explaining the $\log(L)$
scaling of $t^{*}(L)$ observed in our numerical study of the SFF using the
eigenvalues of $\mathcal{M}_{\text{JC}}^{\text{F}}$ in the main paper.

Even for $g/J\ll 1$, the above analysis only works for finite system sizes
because of the requirement,
$g^{4}\ll\text{min}\left(J^{4}\left(1-\cos\frac{2\pi i}{L}\right)\right)$,
which breaks down in the thermodynamic limit of $L\rightarrow\infty$ for any
finite $i$. For $L\rightarrow\infty$, the second-largest eigenvalue,
$\lambda_{1}\approx 1-\frac{2\pi^{2}J^{2}}{L^{2}}$, determines the Thouless-
time scaling for any value of $g,J$ as the largest eigenvalues $\lambda_{i}$
are no longer nearly degenerate. By relating $\lambda_{1}^{t^{*}}\approx 1/e$,
we find

$\displaystyle
t^{*}\approx-1/\log\lambda_{1}\approx\frac{L^{2}}{2\pi^{2}J^{2}}\propto\mathcal{O}(L^{2}).$
(S31)

For finite system sizes, the largest eigenvalues $\lambda_{i}$ are no longer
nearly degenerate when $g/J\gg 1$ as shown in Fig. S1(b). Thus, the scaling of
$t^{*}(L)$ with $L$ is again determined by $\lambda_{1}\approx
1-\frac{2\pi^{2}J^{2}}{L^{2}}$, which again leads to $t^{*}(L)\propto L^{2}$.
Thus, we find a change in the $L$-dependence of $t^{*}(L)$ in finite-sized
systems when the ratio of mixing and hopping parameters are increased from low
to high value. We discuss such a crossover in Fig. 1 of the main paper.

In the main paper, we show $K(t)$ vs. $t$ calculated using $\lambda_{i}$ of
$\mathcal{M}_{\text{JC}}^{\text{F}}$ obtained applying the RPA. We numerically
compute the SFF for $JC$ mixing between fermions and qubits using Eq. S4 to
compare it to that obtained within the RPA using the identity permutation
yielding the first-order term in time and also the second-order term of RMT
form [1]. We show a good match between two different computations in Fig. S2
to validate the RPA analysis in finite-length chains for high values of long-
range interaction and random energies to ensure the RPA. We clarify here that
we use the full $\mathcal{M}_{\text{JC}}^{\text{F}}$ instead of its Trotter-
regime form in Eq. II for our numerics in Fig. S2.

Figure S2: Comparison between the exact numerically computed SFF, $K(t)$ vs.
$t$, with that obtained using the RPA for Jaynes-Cummings mixing between
fermions and qubits. The red curve is exact SFF computed numerically using Eq.
S4 and the blue (black dashed) curve is that calculated using the first-order
and second-order term (only the first-order term) in time within the RPA. All,
$\mathcal{N}_{\text{JC}}^{\text{F}}=15504$, eigenvalues of
$\mathcal{M}_{\text{JC}}^{\text{F}}$ are used for the RPA result. For exact
numerical computation, we fix $U_{0}=10,\alpha=1.4$, and
$\omega_{i},\Omega_{i}$ are chosen as Gaussian random variables with a mean
$\langle\omega_{i}\rangle=\langle\Omega_{i}\rangle=1$ and a standard deviation
$\sigma_{\omega_{i}}=\sigma_{\Omega_{i}}=0.3$. Averaging over $520$
realizations of disorder is performed for the direct SFF computation.

##  III Jaynes-Cummings mixing between bosons and qubits

The driving Hamiltonian for $JC$ mixing between bosons and qubits is the same
as $\hat{H}_{\text{JC}}$ in Eq. S17, where $\hat{a}_{i}^{\dagger},\hat{a}_{i}$
are now bosonic creation and annihilation operators. Similar to the fermionic
case, the hopping (and/or mixing) terms in $\hat{H}_{\text{JC}}$ at any two
different bonds (and/or sites) have no simultaneously non-zero matrix elements
in the occupation number basis. Therefore, we get

$\displaystyle\hat{H}_{\text{JC}}\bullet\hat{H}_{\text{JC}}=g^{2}\sum_{i=1}^{L}\left(\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}\bullet\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}+\hat{a}_{i}\hat{\sigma}_{i}^{\dagger}\bullet\hat{a}_{i}\hat{\sigma}_{i}^{\dagger}\right)$
$\displaystyle+J^{2}\sum_{i=1}^{L}\left(\hat{a}_{i}^{\dagger}\hat{a}_{i+1}\bullet\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+\hat{a}_{i}\hat{a}_{i+1}^{\dagger}\bullet\hat{a}_{i}\hat{a}_{i+1}^{\dagger}\right).$
(S32)

Since
$\langle\underline{n\sigma}|\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}\bullet\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}|\underline{n\sigma}^{\prime}\rangle=\langle\underline{n\sigma}|\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}|\underline{n\sigma}^{\prime}\rangle^{2}=\langle\underline{n\sigma}|\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}|\underline{n\sigma}^{\prime}\rangle$,
and
$\langle\underline{n\sigma}|\hat{a}_{i}^{\dagger}\hat{a}_{i+1}\bullet\hat{a}_{i}^{\dagger}\hat{a}_{i+1}|\underline{n\sigma}^{\prime}\rangle=\langle\underline{n\sigma}|\hat{a}_{i}^{\dagger}\hat{a}_{i+1}|\underline{n\sigma}^{\prime}\rangle^{2}=\langle\underline{n\sigma}|\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger}\hat{a}_{i+1}\sqrt{\hat{n}_{i+1}}|\underline{n\sigma}^{\prime}\rangle$,
we can simplify as

$\displaystyle\hat{H}_{\text{JC}}\bullet\hat{H}_{\text{JC}}=g^{2}\sum_{i=1}^{L}(\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}+\hat{a}_{i}\sqrt{\hat{n}_{i}}\hat{\sigma}_{i}^{\dagger})$
$\displaystyle+J^{2}\sum_{i=1}^{L}(\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger}\hat{a}_{i+1}\sqrt{\hat{n}_{i+1}}+\sqrt{\hat{n}_{i+1}}\hat{a}_{i+1}^{\dagger}\hat{a}_{i}\sqrt{\hat{n}_{i}}).$
(S33)

Since $(\hat{H}_{\text{JC}})^{2}\bullet\mathds{1}$ is a matrix of diagonal
elements of $(\hat{H}_{\text{JC}})^{2}$, we get

$\displaystyle(\hat{H}_{\text{JC}})^{2}\bullet\mathds{1}$
$\displaystyle=g^{2}\sum_{i=1}^{L}(\hat{n}_{i}+\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i})$
$\displaystyle+J^{2}\sum_{i=1}^{L}\left(\hat{n}_{i}+\hat{n}_{i+1}+2\hat{n}_{i}\hat{n}_{i+1}\right).$
(S34)

Substituting Eqs. III,III in Eq. S16, we find

$\displaystyle\mathcal{M}_{\text{JC}}^{\text{B}}$
$\displaystyle=\mathds{1}_{\mathcal{N}_{\text{JC}}^{\text{B}}}+g^{2}\sum_{i=1}^{L}\left(\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger}\hat{\sigma}_{i}+\hat{a}_{i}\sqrt{\hat{n}_{i}}\hat{\sigma}_{i}^{\dagger}-\hat{n}_{i}-\hat{\sigma}_{i}\hat{\sigma}_{i}^{\dagger}\right)$
$\displaystyle+J^{2}\sum_{i=1}^{L}\Big{(}\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger}\hat{a}_{i+1}\sqrt{\hat{n}_{i+1}}+\sqrt{\hat{n}_{i+1}}\hat{a}_{i+1}^{\dagger}\hat{a}_{i}\sqrt{\hat{n}_{i}}$
$\displaystyle-\hat{n}_{i}-\hat{n}_{i+1}-2\hat{n}_{i}\hat{n}_{i+1}\Big{)}+\mathcal{O}(J^{4},g^{4}).$
(S35)

The above expression can be written in terms of operators,
$\hat{K}_{i}^{-}=\sqrt{\hat{n}_{i}}\hat{a}_{i}^{\dagger},\hat{K}_{i}^{+}=\hat{a}_{i}\sqrt{\hat{n}_{i}}$,
and $\hat{K}_{i}^{0}=-\left(\hat{n}_{i}+1/2\right)$, satisfying $SU(1,1)$
algebra,
$[\hat{K}_{i}^{+},\hat{K}_{j}^{-}]=-\hat{K}_{i}^{0}\delta_{ij},[\hat{K}_{i}^{0},\hat{K}_{i}^{\pm}]=\pm\hat{K}_{i}^{\pm}\delta_{ij}$,
as

$\displaystyle\mathcal{M}_{\text{JC}}^{\text{B}}$
$\displaystyle=(1+\frac{(g^{2}+J^{2})L}{2})\mathds{1}_{\mathcal{N}_{\text{JC}}^{\text{B}}}+g^{2}\sum_{i=1}^{L}\big{(}\hat{K}_{i}^{-}\hat{\sigma}_{i}+\hat{K}_{i}^{+}\hat{\sigma}_{i}^{\dagger}$
$\displaystyle+\hat{K}_{i}^{0}-\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i}\big{)}+J^{2}\sum_{i=1}^{L}\big{(}\hat{K}_{i}^{+}\hat{K}_{i+1}^{-}+\hat{K}_{i}^{-}\hat{K}_{i+1}^{+}$
$\displaystyle-2\hat{K}_{i}^{0}\hat{K}_{i+1}^{0}\big{)}+\mathcal{O}(J^{4},g^{4}).$
(S36)

We can further define these operators,
$\hat{K}_{j}^{1}=(\hat{K}_{j}^{+}+\hat{K}_{j}^{-})/2,\hat{K}_{j}^{2}=(\hat{K}_{j}^{+}-\hat{K}_{j}^{-})/2i$,
and rewrite Eq. III as

$\displaystyle\mathcal{M}_{\text{JC}}^{\text{B}}$
$\displaystyle=(1+\frac{(g^{2}+J^{2})L}{2})\mathds{1}_{\mathcal{N}_{\text{JC}}^{\text{B}}}+g^{2}\sum_{i=1}^{L}\big{(}\hat{K}_{i}^{1}\hat{\sigma}_{i}^{x}-\hat{K}_{i}^{2}\hat{\sigma}_{i}^{y}$
$\displaystyle+\hat{K}_{i}^{0}-\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i}\big{)}+2J^{2}\sum_{i=1}^{L}\big{(}\hat{K}_{i}^{1}\hat{K}_{i+1}^{1}+\hat{K}_{i}^{2}\hat{K}_{i+1}^{2}$
$\displaystyle-\hat{K}_{i}^{0}\hat{K}_{i+1}^{0}\big{)}+\mathcal{O}(J^{4},g^{4}),$
(S37)

where
$\hat{\sigma}_{i}^{x}=\hat{\sigma}_{i}^{\dagger}+\hat{\sigma}_{i},\hat{\sigma}_{i}^{y}=\left(\hat{\sigma}_{i}^{\dagger}-\hat{\sigma}_{i}\right)/i$.
The generating Hamiltonian $\mathcal{M}_{\text{JC}}^{\text{B}}$ commutes with
$\sum_{i=1}^{L}\hat{K}_{i}^{0}+\hat{\sigma}_{i}^{\dagger}\hat{\sigma}_{i}$,
which suggests a $U(1)$ symmetry of $\mathcal{M}_{\text{JC}}^{\text{B}}$. But
unlike the fermionic case, the lowering and raising operators like
$\sum_{i=1}^{L}(\hat{K}_{i}^{+}+\hat{\sigma}_{i}),\sum_{i=1}^{L}(\hat{K}_{i}^{-}+\hat{\sigma}_{i}^{\dagger})$
do not commute with $\mathcal{M}_{\text{JC}}^{\text{B}}$, which indicates an
absence of $SU(2)$ or $SU(1,1)$ symmetry of
$\mathcal{M}_{\text{JC}}^{\text{B}}$. Thus, the second-largest eigenvalue
$\lambda_{1}$ of $\mathcal{M}_{\text{JC}}^{\text{B}}$ is not required to be
the same for different total number of excitations $N$. Numerical study in
Fig. S3 confirms that $\lambda_{1}$ is different for different $N$ at a fixed
$L$ when $g=0.4,J=0.1$. However, Fig. S3 also shows that $\lambda_{1}$ for
different $N$ seems to converge with an increasing $L$. For system sizes
$L\in[10,26]$, we numerically find $\lambda_{1}(N)\approx 1-c_{N}/L^{2}$ for
different $N$’s, and the difference $(\lambda_{1}(N)-\lambda_{1}(N+1))$ falls
as $c^{\prime}_{N}/L^{3}$, where $c_{N}$ and $c^{\prime}_{N}$ are
$N$-dependent constants. These scaling suggests that $\lambda_{1}$ becomes
independent of $N$ at large $L$ due to an emergent approximate symmetry of
$\mathcal{M}_{\text{JC}}^{\text{B}}$. A similar $N$-independence of $L$
largest eigenvalues of $\mathcal{M}_{\text{JC}}^{\text{B}}$ is also observed
with increasing $L$ when $g/J\ll 1$.

(a)  Figure S3: The second-largest eigenvalue $\lambda_{1}$ of
$\mathcal{M}_{\text{JC}}^{\text{B}}$ vs. $L$ for different total excitations
$N$ when $g=0.4,J=0.1$. $\lambda_{1}$ for different $N$’s approach the same
value with an increasing $L$ suggesting emergence of an approximate symmetry
of $\mathcal{M}_{\text{JC}}^{\text{B}}$.

Similar to the fermionic case in the earlier section, $\lambda_{i}$ of
$\mathcal{M}_{\text{JC}}^{\text{B}}$ can be calculated analytically for $N=1$,
and $\lambda_{i}$ are then identical to those in Eq. II. Incorporating the
above observations for an emergent approximate symmetry of
$\mathcal{M}_{\text{JC}}^{\text{B}}$ with $\lambda_{i}$ for $N=1$, we then
expect a change in the $L$-dependence of $t^{*}(L)$ from $\log L$ to $L^{2}$
with an increasing $g/J$ from much lower than one to much larger than one at
any filling fraction $N/L$ in finite-size systems of bosons and qubits. We
demonstrate such a crossover in the $L$-scaling of $t^{*}(L)$ by studying
$K(t)$ vs. $t$ within the RPA in Fig. S4 at half filling $N/L=1/2$ for finite
system sizes.

Figure S4: Spectral form factor $K(t)$ using Eq. S14 for different system
sizes $L$ of the kicked chain with $JC$ mixing between bosons and qubits for
$g=0.1,J=0.4$ in $(a,b)$, and $g=0.4,J=0.1$ in $(c,d)$. We take half-filling
$N/L=1/2$. In $(b)$ and $(d)$, we show data collapse in scaled time $t/\log L$
and $t/L^{1.86}$, respectively.

##  IV Rabi mixing between fermions and qubits

The driving Hamiltonian $\hat{H}_{\text{R}}$ for Rabi $(R)$ mixing between
fermions and qubits is given in Eq. (3) of the main paper.

$\displaystyle\hat{H}_{\text{R}}$
$\displaystyle=g\sum_{i=1}^{L}(\hat{a}_{i}^{\dagger}+\hat{a}_{i})(\hat{\sigma}_{i}^{\dagger}+\hat{\sigma}_{i})-J\sum_{i=1}^{L}\left(\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+{\text{H.c.}}\right).$
(S38)

Following the steps presented for deriving the generating Hamiltonian in the
Trotter regime for $JC$ mixing between fermions and qubits in Sec. II, we find

$\displaystyle\mathcal{M}_{\text{R}}^{\text{F}}$
$\displaystyle=\left(1-\frac{(2g^{2}+J^{2})L}{2}\right)\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{F}}}+\sum_{i=1}^{L}\Big{(}\frac{J^{2}}{2}\sum_{\nu}\hat{\tau}_{i}^{\nu}\hat{\tau}_{i+1}^{\nu}$
$\displaystyle\qquad+g^{2}\hat{\tau}_{i}^{x}\hat{\sigma}_{i}^{x}\Big{)}+\mathcal{O}(J^{4},g^{4}),$
(S39)

where $\nu=x,y,z$. By performing a rotation
$\big{(}\mathcal{R}^{y}_{\theta}\equiv\exp(i\theta\sum_{j=1}^{L}(\hat{\tau}_{j}^{y}+\hat{\sigma}_{j}^{y})/2)\big{)}$
around $y$-axis by $\theta=\pi/2$, we transform the operators as
$\mathcal{R}^{y}_{\theta}\hat{\tau}_{i}^{x}\mathcal{R}_{\theta}^{y\dagger}=\hat{\tilde{\tau}}_{i}^{z},\mathcal{R}^{y}_{\theta}\hat{\tau}_{i}^{y}\mathcal{R}_{\theta}^{y\dagger}=\hat{\tilde{\tau}}_{i}^{y},\mathcal{R}^{y}_{\theta}\hat{\tau}_{i}^{z}\mathcal{R}_{\theta}^{y\dagger}=-\hat{\tilde{\tau}}_{i}^{x},\mathcal{R}^{y}_{\theta}\hat{\sigma}_{i}^{x}\mathcal{R}_{\theta}^{y\dagger}=\hat{\tilde{\sigma}}_{i}^{z},\mathcal{R}^{y}_{\theta}\hat{\sigma}_{i}^{y}\mathcal{R}_{\theta}^{y\dagger}=\hat{\tilde{\sigma}}_{i}^{y},\mathcal{R}^{y}_{\theta}\hat{\sigma}_{i}^{z}\mathcal{R}_{\theta}^{y\dagger}=-\hat{\tilde{\sigma}}_{i}^{x}$,
which lead to
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}=\mathcal{R}^{y}_{\theta}\mathcal{M}_{\text{R}}^{\text{F}}\mathcal{R}_{\theta}^{y\dagger}$:

$\displaystyle\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$
$\displaystyle=\left(1-\frac{(2g^{2}+J^{2})L}{2}\right)\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{F}}}+\sum_{i=1}^{L}\Big{(}\frac{J^{2}}{2}\sum_{\nu}\hat{\tilde{\tau}}_{i}^{\nu}\hat{\tilde{\tau}}_{i+1}^{\nu}$
$\displaystyle\qquad+g^{2}\hat{\tilde{\tau}}_{i}^{z}\hat{\tilde{\sigma}}_{i}^{z}\Big{)}+\mathcal{O}(J^{4},g^{4}).$
(S40)

The generating Hamiltonian $\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$
commutes with $\hat{\tilde{\sigma}}_{i}^{z}$ for each $i\in[1,L]$. We call
this symmetry collectively as $u^{\otimes L}(1)$.
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ also commutes with
$\sum_{i=1}^{L}\hat{\tilde{\tau}}_{i}^{z}$, which is a global $U(1)$ symmetry.
So, $\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ has $U(1)\otimes u^{\otimes
L}(1)$ symmetry. Thus, the eigenvalues and eigenstates of
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ can be labelled by excitation of
individual $\tilde{\sigma}$ qubit,
$\tilde{\sigma}_{i}\left(\equiv\hat{\tilde{\sigma}}_{i}^{\dagger}\hat{\tilde{\sigma}}_{i}\right)$,
and total excitations of $\tilde{\tau}$ qubits,
$\tilde{N}_{\text{f}}\left(\equiv\sum_{i=1}^{L}\hat{\tilde{\tau}}_{i}^{\dagger}\hat{\tilde{\tau}}_{i}\right)$.
The states with the largest eigenvalue $\lambda_{0}=1$ have
$\tilde{N}_{\text{f}}=0,\tilde{\sigma}_{i}=0$ or
$\tilde{N}_{\text{f}}=L,\tilde{\sigma}_{i}=1,$ for all $i\in[1,L]$. The
degeneracy of $2$ for the largest eigenvalue is due to the invariance of
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ under the transformation,
$\prod_{i=1}^{L}\hat{\tilde{\tau}}_{i}^{x}\hat{\tilde{\sigma}}_{i}^{x}$. The
last operator flips all the spins to transform
$\tilde{N}_{\text{f}}=0\rightarrow L,\tilde{\sigma}_{i}=0\rightarrow 1$ for
all $i\in[1,L]$ leading to the degeneracy of $2$. This is equivalent to
$\mathds{Z}_{2}$ symmetry also present in $\hat{H}(t)$ (the $U(1)$ breaking
$\hat{H}_{\text{R}}$ does not mix between the even and odd number of total
excitation sectors). To study chaos, we choose either the even or odd sector
of the Hilbert space of $\hat{H}(t)$ only. Similarly, we choose only one
largest eigenvalue of $\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$, and study
excitations (second-largest eigenvalues) around it. We here choose the largest
eigenvalue state with $\tilde{N}_{\text{f}}=0,\tilde{\sigma}_{i}=0,$ for all
$i\in[1,L]$. The numerical study for different values of $g,J$ shows that
there are two type of $\tilde{\tau}$ and $\tilde{\sigma}$ spin configurations
giving the second-largest eigenvalue of
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ depending on the ratio $g/J$ in
the Trotter regime. We describe them below.

Case 1: For $g/J\ll 1$, the configurations leading to the second-largest
eigenvalues appear by one change in the excitation from the largest eigenvalue
state, e.g.,
$\tilde{N}_{\text{f}}=0,\tilde{\sigma}_{j}=1,\tilde{\sigma}_{i\neq
j}=0,i=1,...,j-1,j+1,...,L$, and $\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{i}=0$
for all $i\in[1,L]$. We find

$\displaystyle\sum_{i=1}^{L}\sum_{\nu}\hat{\tilde{\tau}}_{i}^{\nu}\hat{\tilde{\tau}}_{i+1}^{\nu}|\tilde{N}_{\text{f}}=0,\\{\tilde{\sigma}\\}=\\{0,...,0,1,0,...,0\\}\rangle$
$\displaystyle\qquad\quad=L|\tilde{N}_{\text{f}}=0,\\{\tilde{\sigma}\\}=\\{0,...,0,1,0,...,0\\}\rangle,$
(S41)
$\displaystyle\sum_{i=1}^{L}\hat{\tilde{\tau}}_{i}^{z}\hat{\tilde{\sigma}}_{i}^{z}|\tilde{N}_{\text{f}}=0,\\{\tilde{\sigma}\\}=\\{0,...,0,1,0,...,0\\}\rangle$
$\displaystyle\qquad=(L-2)|\tilde{N}_{\text{f}}=0,\\{\tilde{\sigma}\\}=\\{0,...,0,1,0,...,0\\}\rangle.$
(S42)

Therefore,

$\displaystyle\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}|\tilde{N}_{\text{f}}=0,\\{\tilde{\sigma}\\}=\\{0,...,0,1,0,...,0\\}\rangle$
$\displaystyle\qquad=(1-2g^{2})|\tilde{N}_{\text{f}}=0,\\{\tilde{\sigma}\\}=\\{0,...,0,1,0,...,0\\}\rangle,$
(S43)

and we get $\lambda_{1}=1-2g^{2}$ with $L$-fold degeneracy, which is due to
appearance of $\tilde{\sigma}_{j}=1$ at any $j=1,2,..,L$. For
$\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{i}=0$,
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ becomes the $XXX$ Heisenberg
spin-1/2 chain for $\tilde{\tau}$ spin in a uniform external magnetic field
($2g^{2}$).

$\displaystyle\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}|_{\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{i}=0}=\left(1-\frac{(2g^{2}+J^{2})L}{2}\right)\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{F}}}$
$\displaystyle+\sum_{i=1}^{L}\left(\frac{J^{2}}{2}\sum_{\nu}\hat{\tilde{\tau}}_{i}^{\nu}\hat{\tilde{\tau}}_{i+1}^{\nu}-g^{2}\hat{\tilde{\tau}}_{i}^{z}\right)+\mathcal{O}(J^{4},g^{4}).$
(S44)

For total number of $\tilde{\tau}$ excitations, $\tilde{N}_{\text{f}}=1$, the
eigenvalues of $\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ are

$\displaystyle\lambda_{k_{m}}=1-2g^{2}-2J^{2}(1-\cos k_{m}),$ (S45)

where $k_{m}=2\pi m/L,m=0,1,...,L-1$. For $k_{m}=0$, we get the second-largest
eigenvalue $\lambda_{1}=1-2g^{2}$. Therefore, we get total $(L+1)$-fold
degeneracy in this case. The Thouless-time scaling in the thermodynamic limit
of $L\to\infty$ can be obtained by setting $(L+1)\lambda_{1}^{t^{*}}\approx 1$
as

$\displaystyle
t^{*}(L)\approx-\frac{\log(L+1)}{\log(1-2g^{2})}\approx\mathcal{O}(\log L).$
(S46)

Case 2: For $g/J\gg 1$, the configurations leading to the second-largest
eigenvalues appear by one change in the excitation from the largest eigenvalue
state, e.g.,
$\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{j}=1,\tilde{\sigma}_{i\neq
j}=0,i=1,...,j-1,j+1,...,L$. We can cast
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ in this case as a tight-binding
chain with an impurity in the onsite energy:

$\displaystyle\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}|_{\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{j}=1}=1-4g^{2}-2J^{2}+4g^{2}\hat{\tilde{\tau}}_{j}^{\dagger}\hat{\tilde{\tau}}_{j}$
$\displaystyle+J^{2}\sum_{i=1}^{L}(\hat{\tilde{\tau}}_{i}^{\dagger}\hat{\tilde{\tau}}_{i+1}+\hat{\tilde{\tau}}_{i}\hat{\tilde{\tau}}_{i+1}^{\dagger}).$
(S47)

The spectrum of the Hamiltonian in Eq. IV can be derived using the Dyson
equation to treat the impurity as a perturbation following Economou [48]. We
separate
$\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}|_{\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{j}=1}$
in two parts as

$\displaystyle\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}|_{\tilde{N}_{\text{f}}=1,\tilde{\sigma}_{j}=1}=\mathcal{M}_{0}+\mathcal{M}_{1},$
(S48)

where
$\mathcal{M}_{0}=1-4g^{2}-2J^{2}+J^{2}\sum_{i=1}^{L}(\hat{\tilde{\tau}}_{i}^{\dagger}\hat{\tilde{\tau}}_{i+1}+\hat{\tilde{\tau}}_{i}\hat{\tilde{\tau}}_{i+1}^{\dagger}),\mathcal{M}_{1}=4g^{2}\hat{\tilde{\tau}}_{j}^{\dagger}\hat{\tilde{\tau}}_{j}=4g^{2}|j\rangle\langle
j|$. Here, $|j\rangle=\hat{\tilde{\tau}}_{j}^{\dagger}|\varphi\rangle$, where
$|\varphi\rangle$ is a vacuum state. We consider $G_{0}(z)$ as the Green’s
function for $\mathcal{M}_{0}$, and $G(z)$ as that for
$\mathcal{M}_{0}+\mathcal{M}_{1}$. The Dyson equation [48] gives,

$\displaystyle G$ $\displaystyle=G_{0}$
$\displaystyle~{}+G_{0}\left(\mathcal{M}_{1}+\mathcal{M}_{1}G_{0}\mathcal{M}_{1}+\mathcal{M}_{1}G_{0}\mathcal{M}_{1}G_{0}\mathcal{M}_{1}+...\right)G_{0}$
$\displaystyle=G_{0}+G_{0}\Big{(}(2g)^{2}+(2g)^{4}G_{0}(j,j;z)$
$\displaystyle\quad+(2g)^{6}\left(G_{0}(j,j;z)\right)^{2}+...\Big{)}|j\rangle\langle
j|G_{0}$
$\displaystyle=G_{0}+G_{0}\frac{(2g)^{2}}{1-(2g)^{2}G_{0}(j,j;z)}|j\rangle\langle
j|G_{0},$ (S49)

where $G_{0}(j,j;z)=\langle j|G_{0}(z)|j\rangle$. We find from Eq. IV that $G$
has an isolated pole at

$\displaystyle G_{0}(j,j;E_{b})=\frac{1}{4g^{2}},$ (S50)

which gives the eigenvalue of the corresponding bound state. The free Green’s
function can be determined using the eigenvalues and eigenstates of
$\mathcal{M}_{0}$ as $G_{0}(z)=\sum_{m}|k_{m}\rangle\langle
k_{m}|/(z-E(k_{m}))$, where $|k_{m}\rangle$ are eigenstates of
$\mathcal{M}_{0}$ with eigenvalues $E(k_{m})=1-4g^{2}-2J^{2}+2J^{2}\cos
k_{m},\;k_{m}=2\pi m/L$ for $m=0,1,2,...,L-1$. Therefore, we get

$\displaystyle G_{0}(j_{1},j_{2};z)=\sum_{m}\frac{\langle
j_{1}|k_{m}\rangle\langle
k_{m}|j_{2}\rangle}{z-E(k_{m})}=\frac{1}{L}\sum_{m}\frac{e^{ik_{m}(j_{1}-j_{2})}}{z-E(k_{m})}.$
(S51)

It is easy to evaluate the above sum in the thermodynamic limit of
$L\to\infty$ as

$\displaystyle
G_{0}(j_{1},j_{2};z)=\frac{1}{2\pi}\int_{0}^{2\pi}dk\frac{e^{ik(j_{1}-j_{2})}}{z-E(k)}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}dk\frac{e^{ik(j_{1}-j_{2})}}{z-1+4g^{2}+2J^{2}-2J^{2}\cos
k}.$ (S52)

We observe $G_{0}(j_{1},j_{2};z)=G_{0}(j_{2},j_{1};z)$ as $k\rightarrow
2\pi-k$, which implies that $G_{0}(j_{1},j_{2};z)$ is a function of
$|j_{1}-j_{2}|$. We define $w=e^{ik},\epsilon=(z-1+4g^{2}+2J^{2})/(2J^{2})$ to
rewrite as

$\displaystyle G_{0}(j_{1},j_{2};z)=-\frac{1}{2\pi iJ^{2}}\oint
dw\frac{w^{|j_{1}-j_{2}|}}{w^{2}-2\epsilon w+1},$ (S53)

where the contour is a unit circle on complex $w$ plane. The integrand has
poles at $w=w_{\pm}=\epsilon\pm\sqrt{\epsilon^{2}-1}$. We further notice
$w_{+}w_{-}=1$, which implies that $|w_{+}||w_{-}|=1$. Thus, we get

$\displaystyle|w_{+}|=\begin{cases}<1,\quad\frac{\pi}{2}<\arg(\epsilon)<\frac{3\pi}{2}\\\
>1,\quad\text{otherwise}\end{cases},$ (S54)
$\displaystyle|w_{-}|=\begin{cases}>1,\quad\frac{\pi}{2}<\arg(\epsilon)<\frac{3\pi}{2}\\\
<1,\quad\text{otherwise}\end{cases}.$ (S55)

For $\epsilon$ real, and $-1<\epsilon<1$, $|w_{\pm}|=1$, and the poles exist
on the integration contour. This happens for $1-4g^{2}-4J^{2}<z<1-4g^{2}$,
which is a range of continuous spectrum of $\mathcal{M}_{0}$. We have $w_{+}$
or $w_{-}$ inside the integration contour depending on the argument of
$\epsilon$. We finally get

$\displaystyle
G_{0}(j_{1},j_{2};z)=-\frac{1}{J^{2}}\frac{w_{\pm}^{|j_{1}-j_{2}|}}{w_{\pm}-w_{\mp}}.$
(S56)

Equating Eqs. S56 and S50 at $z=E_{b}$, we find

$\displaystyle\epsilon$ $\displaystyle=\pm\sqrt{1+\frac{4g^{4}}{J^{4}}},$
(S57) $\displaystyle\implies E_{b\pm}$ $\displaystyle=1-4g^{2}-2J^{2}\pm
2J^{2}\sqrt{1+\frac{4g^{4}}{J^{4}}}.$ (S58)

We have two bound states in Eq. S58, and the one corresponding to the second-
largest eigenvalue of $\tilde{\mathcal{M}}_{\text{R}}^{\text{F}}$ in this case
is

$\displaystyle\lambda_{1}\equiv
E_{b+}=1-4g^{2}-2J^{2}+2J^{2}\sqrt{1+\frac{4g^{4}}{J^{4}}}.$ (S59)

Here, $\tilde{\sigma}_{j}=1$ can be possible at any one site of $L$ sites,
which leads to an $L$-fold degeneracy. Thus, we again get,
$t^{*}(L)\propto\mathcal{O}(\log L)$ for this case also.

There is a transition in the gap between $\lambda_{0}$ and $\lambda_{1}$ from
case 1 to case 2 as $g/J$ changes. We can find the precise transition point by
equating the eigenvalue obtained for the two cases. The transition point turns
out to be

$\displaystyle\left(\frac{g}{J}\right)_{\text{c}}=\sqrt{\frac{2}{3}}.$ (S60)

We notice that the $L$-fold degeneracy is due to the spatially uniform
coupling $(g)$ between fermions and qubits. Instead, we can consider the
driven Hamiltonian with site-dependent coupling $g_{i}$ at site $i$ for $R$
mixing as

$\displaystyle\hat{H}_{\text{R}}$
$\displaystyle=\sum_{i=1}^{L}g_{i}(\hat{a}_{i}^{\dagger}+\hat{a}_{i})(\hat{\sigma}_{i}^{\dagger}+\hat{\sigma}_{i})-J\sum_{i=1}^{L}\left(\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+{\text{H.c}}.\right),$
(S61)

which would transform the generating Hamiltonian in Eq. IV as

$\displaystyle\mathcal{M}_{\text{R}}^{\text{F}}$
$\displaystyle=\big{(}1-\frac{J^{2}L}{2}-\sum_{i=1}^{L}g_{i}^{2}\big{)}\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{F}}}+\sum_{i=1}^{L}\Big{(}\frac{J^{2}}{2}\sum_{\nu}\hat{\tau}_{i}^{\nu}\hat{\tau}_{i+1}^{\nu}$
$\displaystyle\qquad+g_{i}^{2}\hat{\tau}_{i}^{x}\hat{\sigma}_{i}^{x}\Big{)}+\mathcal{O}(J^{4},g_{i}^{4}).$
(S62)

The second-largest eigenvalues of $\mathcal{M}_{\text{R}}^{\text{F}}$ are
$1-2g_{i}^{2}$. Since $g_{i}$ is different for different $i$, the $L$-fold
degeneracy of the second-largest eigenvalue is lifted for site-dependent $R$
mixing, and we get $t^{*}(L)\sim L^{0}$.

We next numerically compute the SFF for $R$ mixing between fermions and qubits
using Eq. S4 to compare it to that obtained within the RPA using the full
$\mathcal{M}_{\text{R}}^{\text{F}}$ instead of its Trotter-regime form in Eq.
IV. For $R$ mixing violating the total excitation number conservation, we show
a better match between two different computations in Fig. S5 in comparison to
that in Fig. S2 for the $JC$ mixing satisfying the total excitation number
conservation.

Figure S5: Comparison between the exact numerically computed SFF, $K(t)$ vs.
$t$, with that obtained using the RPA for Rabi mixing between fermions and
qubits. The red curve is exact SFF computed numerically using Eq. S4 and the
blue (black dashed) curve is that calculated using the first-order and the
second-order term (only the first-order term) in time within the RPA. All,
$\mathcal{N}_{\text{R}}^{\text{F}}=8192$, eigenvalues of
$\mathcal{M}_{\text{R}}^{\text{F}}$ are used for the RPA result. For exact
numerical computation, we fix $U_{0}=10,\alpha=1.4$, and
$\omega_{i},\Omega_{i}$ are chosen as Gaussian random variables with a mean
$\langle\omega_{i}\rangle=\langle\Omega_{i}\rangle=1$ and a standard deviation
$\sigma_{\omega_{i}}=\sigma_{\Omega_{i}}=0.3$. Averaging over $500$
realizations of disorder is performed for the direct SFF computation.

##  V Rabi mixing between bosons and qubits

The driving Hamiltonian for $R$ mixing between bosons and qubits is the same
as $\hat{H}_{\text{R}}$ in Eq. S38, where $\hat{a}_{i}^{\dagger},\hat{a}_{i}$
are now bosonic creation and annihilation operators. Following the steps
presented for deriving the generating Hamiltonian in the Trotter regime for
$JC$ mixing between bosons and qubits in Sec. III, we get

$\displaystyle\mathcal{M}_{\text{R}}^{\text{B}}$
$\displaystyle=\left(1+\frac{J^{2}L}{2}\right)\mathds{1}_{\mathcal{N}_{\text{R}}^{\text{B}}}+2g^{2}\sum_{i=1}^{L}\left(\hat{K}_{i}^{1}\hat{\sigma}_{i}^{x}+\hat{K}_{i}^{0}\right)$
$\displaystyle\qquad+2J^{2}\sum_{i=1}^{L}\left(\hat{K}_{i}^{1}\hat{K}_{i+1}^{1}+\hat{K}_{i}^{2}\hat{K}_{i+1}^{2}-\hat{K}_{i}^{0}\hat{K}_{i+1}^{0}\right)$
$\displaystyle\qquad+\mathcal{O}(J^{4},g^{4}).$ (S63)

The dimension of the Hilbert space $\mathcal{N}_{\text{R}}^{\text{B}}$ is
infinite for $R$ mixing between bosons and qubits due to the non-conservation
of total excitation $N$ and no bound to number of bosons at any site. We could
not find the spectrum of $\mathcal{N}_{\text{R}}^{\text{B}}$ analytically.
Rather, we introduce a truncation to the maximum number of total excitations
$N_{max}$ for our numerical study of $\lambda_{i}$ of
$\mathcal{M}_{\text{R}}^{\text{B}}$ as in Roy  _et al._ [13]. For a fixed $L$,
we vary $N_{max}$ to find the asymptotic behavior $(N_{max}\to\infty)$ of the
second-largest eigenvalue of $\mathcal{M}_{\text{R}}^{\text{B}}$. Following
$R$ mixing between fermions and qubits, we explore two different parameter
regimes, $g/J\ll 1$ and $g/J\gg 1$ in our numerics as given in Figs. S6,S7.

We further notice that $\mathcal{M}_{\text{R}}^{\text{B}}$ commutes with
$\hat{\sigma}_{i}^{x}$ for all $i\in[1,L]$. Thus, we can label the eigenvalues
and eigenvectors of $\mathcal{M}_{\text{R}}^{\text{B}}$ by
$m_{i}=(\sigma_{i}^{x}+1)/2$, where $\sigma_{i}^{x}$ is an eigenvalue of
$\hat{\sigma}_{i}^{x}$. Nevertheless, we again need to introduce a truncation
to the maximum number of bosons $N^{b}_{max}$ for our numerical study of
$\lambda_{i}$ of $\mathcal{M}_{\text{R}}^{\text{B}}$ by fixing the qubit
excitations. We find from our numerics for $g=0.1,J=0.4$ that the eigenstates
with largest eigenvalue of $\mathcal{M}_{\text{R}}^{\text{B}}$ have $m_{i}=0$
or $m_{i}=1$ for all $i\in[1,L]$. The degeneracy of $2$ in the largest
eigenvalue is due to a $\mathds{Z}_{2}$ symmetry of
$\mathcal{M}_{\text{R}}^{\text{B}}$ as like of
$\mathcal{M}_{\text{R}}^{\text{F}}$. Let us choose the first configuration
($m_{i}=0$) with the largest eigenvalue. The state corresponding to the
second-largest eigenvalue $\lambda_{1}$ of $\mathcal{M}_{\text{R}}^{\text{B}}$
also has the same qubit configuration. The states corresponding to the third-
largest eigenvalue are $L$-fold degenerate and nearly degenerate with the
second-largest eigenvalue state. The qubit configuration for these states with
the third-largest eigenvalue is $m_{j}=1,m_{i\neq j}=0$ for
$i=1,...,j-1,j+1...,L$. This is similar to the results obtained for $R$ mixing
between fermions and qubits.

Interestingly, we find from our numerics that the trend of $\lambda_{1}$ with
increasing $N^{b}_{max}$ matches nicely to that with increasing $N_{max}$ at
large values of $N^{b}_{max}$ and $N_{max}$ when $g/J\ll 1$ as shown in Fig.
S6 for $L=4$. Since larger $N^{b}_{max}$ is accessible with a fixed qubit
configuration in numerics with a fixed $L$, we employ the numerics by changing
$N^{b}_{max}$ to find $\lambda_{1}$ for different $L$’s as
$N^{b}_{max}\to\infty$. We show them in Fig. S6 for $L=4,5,6,7$, which depict
$\lambda_{1}$ at asymptotic $N^{b}_{max}$ slowly drifting towards smaller
values with increasing $L$. The decrease in $\lambda_{1}$ with an increasing
$L$ at $N^{b}_{max}\to\infty$ is probably appearing from the linear
extrapolation of the last few large $N^{b}_{max}$ points, which are not so
large for longer $L$. Nevertheless, since the $L$-dependence of $\lambda_{1}$
is small and there are nearly $(L+1)$-fold degeneracy around the second-
largest eigenvalue, we predict $t^{*}(L)\approx\log L$ when $g/J\ll 1$.

The similarity between the asymptotic trend of $\lambda_{1}$ with increasing
$N_{max}$ and $N_{max}^{b}$ is lost when $g/J\gg 1$. Therefore, we can only
rely on numerics of $\lambda_{1}$ with an increasing $N_{max}$ to understand
the Thouless-time scaling. We show the behavior $\lambda_{1}$ with an
increasing $N_{max}$ for three different lengths in Fig. S7, which seems to
suggest a $L$-independent gap between the first- and second-largest
eigenvalues of $\mathcal{M}_{\text{R}}^{\text{B}}$ for $g/J\gg 1$. We also
observe an $L$-fold degeneracy for the second-largest eigenvalues in our
numerics. Therefore, we expect $t^{*}(L)\approx\log L$ in this regime of
parameters too. Nevertheless, our predictions here for the Thouless-time
scaling for $R$ mixing between bosons and qubits are based on finite-size
numerics with limited data, and there is good scope to improve the current
study with more data for more extended system sizes.

Figure S6: Second-largest eigenvalue $\lambda_{1}$ of
$\mathcal{M}_{\text{R}}^{\text{B}}$ with inverse maximum number of bosons
$(1/N_{max}^{b})$ (also inverse maximum number of total excitations
$1/N_{max}$ for $L=4$) for four different lengths $(L)$. The dashed curves
indicate a linear extrapolation of the last few large $N_{max}^{b}$ (also
$N_{max}$ for $L=4$) points. The parameters are $g=0.1,J=0.4$. The inset shows
a finite gap for $\lambda_{1}$ as $1/N_{max}^{b}\to 0$. Figure S7: Second-
largest eigenvalue $\lambda_{1}$ of $\mathcal{M}_{\text{R}}^{\text{B}}$ with
inverse maximum number of total excitations $(1/N_{max})$ for three different
lengths $(L)$. The dashed curves indicate a linear extrapolation of the last
few large $N_{max}$ points. The parameters are $g=0.4,J=0.1$.