# Fundamental limits of superconducting quantum computers

Michele Vischi1,2, Luca Ferialdi1,2, Andrea Trombettoni1,2,3, Angelo
Bassi1,2,4 1 Department of Physics, University of Trieste, Strada Costiera 11,
34151 Trieste, Italy  2Istituto Nazionale di Fisica Nucleare, Trieste Section,
Via Valerio 2, 34127 Trieste, Italy  3CNR-IOM DEMOCRITOS Simulation Center,
Via Bonomea 265, I-34136 Trieste, Italy  4National Institute of Optics - CNR -
Research Unit of Trieste, Strada Statale 14 - 34149 Trieste

###### Abstract

The Continuous Spontaneous Localization (CSL) model is an alternative
formulation of quantum mechanics which introduces a noise coupled non linearly
to the wave function to account for its collapse. We consider CSL effects on
quantum computers made of superconducting transmon qubits. As a direct effect
CSL reduces quantum superpositions of the computational basis states of the
qubits: we show the reduction rate to be negligibly small. However, an
indirect effect of CSL, dissipation induced by the noise, also leads transmon
qubits to decohere, by generating additional quasiparticles. Since the
decoherence rate of transmon qubits depends on the quasiparticle density, by
computing their generation rate induced by CSL, we can estimate the
corresponding quasiparticle density and thus the limit set by CSL on the
performances of transmon quantum computers. We show that CSL could spoil the
quantum computation of practical algorithms on large devices. We further
explore the possibility of testing CSL effects on superconducting devices.

††preprint: APS/123-QED

##  I Introduction

Quantum computers bring with themselves the promise to allow for a significant
speedup in the resolution of many relevant complex problems with respect to
current classical computers Shor (1999); Brassard  _et al._ (2002); Montanaro
(2015). When large enough quantum computers will be available, they are
expected to impact many fields such as cybersecurity Shor (1999), drug
synthesis Cao  _et al._ (2018), simulation of quantum systems Childs  _et al._
(2018), to name a few. Such computers will need a large processor comprising
many qubits, the fundamental units of quantum computation (the analog of bits
in classical computers). Different physical realizations of qubits are
currently under development, e.g. superconducting qubits Kjaergaard  _et al._
(2020), trapped ions Bruzewicz  _et al._ (2019), photonic chips Bombin  _et
al._ (2021) and spin qubits Recher and Trauzettel (2010). Superconducting
qubits, which are electrical circuits made of superconducting materials, are
among the most promising for scaling up quantum processors.

Quantum chips made of so called transmon superconducting qubits Koch  _et al._
(2007) are currently fabricated by the major companies investing in quantum
computing, such as Google Arute  _et al._ (2019), IBM Gambetta (2020) and
Rigetti Sete  _et al._ (2016). State of the art transmon devices contain $\sim
10^{2}$ qubits, and more than once were used to reach quantum supremacy Arute
_et al._ (2019); Wu  _et al._ (2021), i.e. the resolution of a problem faster
than what possible with any classical computer. Although this is an
outstanding result, the problems solved so far do not have any practical
application. Many technical hurdles must be overcome for such quantum machines
to be practical. The main challenge to address is the fragility of quantum
states that are stored in quantum computers: even a single qubit suffers
decoherence, i.e. the undesired loss of its quantum properties over time.

Several sources of noise lead to decoherence. We can divide them into two
categories: environmental noise sources, whose effects can be eventually
mitigated with technological development, and possible fundamental noise
sources, that are instead unavoidable. Most of the research focuses on how to
protect qubits from environmental noises Córcoles  _et al._ (2011); Wang  _et
al._ (2015); Klimov  _et al._ (2018); Catelani  _et al._ (2011a); Glazman and
Catelani (2021); Martinis  _et al._ (2009), and less attention has been paid
to identifying fundamental ones Lloyd (2000); Gambini  _et al._ (2005). While
technology develops, we may reach a level of control and accuracy at which
fundamental noises can not be overlooked anymore. Among these, of particular
relevance are those related to the spontaneous collapse of the quantum states,
as predicted by the continuous spontaneous localization (CSL) model Bassi and
Ghirardi (2003).

The CSL model (more generally, collapse models Bassi and Ghirardi (2003)) is
an alternative formulation of quantum mechanics developed in order to solve
the tension between the quantum superposition principle and the wave packet
reduction postulate. At the core of the model there is a classical noise,
suitably coupled to the wave function of particles accounting for its
collapse, which now becomes part of the dynamics, not a separate postulate.
The CSL state vector reduction of quantum states becomes more effective as the
difference in mass density of the states in superposition increases Ferialdi
and Bassi (2020). In this way the model is consistent with quantum mechanics
in the microscopic regime, where the standard theory gives extremely accurate
predictions, at the same time justifying why macroscopic objects are always
localized in space.

Collapse models set an intrinsic limit to the stability of quantum systems
over time. As a consequence, they are expected to set a limit to the
scalability of quantum computers, which is the subject of this work where we
will consider quantum computers based on superconducting technology. As we
will see, the CSL state vector reduction does not limit significantly the
performance of transmon qubits. However an indirect effect of the collapse,
dissipation, is more relevant, given the extreme sensibility of
superconducting devices to perturbations. We will compute this effect and will
show that, when neglecting other sources of noise is (or will be) possible,
dissipation can be detected at current experimental temperatures. Moreover we
will see that CSL could spoil the quantum computation of complex practical
algorithms on large devices.

The paper is organized as follows. In section II we review for convenience the
key elements of the Bardeen-Cooper-Schrieffer (BCS) theory of superconductors
needed for the subsequent analysis, highlighting the role played by the
excited states of a superconductor, the quasiparticles. In section III we
discuss how to treat CSL noise in the BCS formalism. The results provide the
basis for the study of the effect of CSL terms in superconducting devices and
computers. We will then outline the key elements of the CSL model that are
needed to compute CSL effects on transmon qubits. In section IV we focus on
transmon qubits and we compute the reduction rate of CSL, i.e. the rate at
which CSL localizes a superposition of transmon qubits. We estimate that this
effect is at date negligibly small for practical purposes. We then
characterize the effects of the CSL dissipation on transmon devices. In
section V we show how the CSL noise perturbs the superconducting materials in
transmon qubits, by generating additional quasiparticles. These excited states
accumulate over time in the devices’ volume leading to a steady state
quasiparticle density larger than the thermal one. This excess limits the
coherence time of transmon qubits. We perform the computation of the CSL
quasiparticle density in section VI. In section VII, we explore the
possibility of testing collapse models with superconducting devices, showing
the experimental conditions to be met in order to detect the CSL excess of
quasiparticles. In section VIII we estimate the fundamental limit due to CSL
on the performances of transmon quantum computers given the quasiparticle
density computed in section VI. In section IX we conclude with final remarks
and outlook.

##  II BCS theory

Transmon qubits exploit the peculiar properties of superconducting materials.
This section is devoted to introducing the key concepts and quantities of the
microscopic theory of superconductors that will be used in the following.

Conventional superconductors are described by the BCS theory. Below a critical
temperature $T_{c}$, it is energetically convenient for electrons in some
metals to bind in pairs, called Cooper pairs. The attraction is mediated by
phonons, and the two electrons in a Cooper pair have opposite momenta and
spins. The total spin is zero, which allows different pairs to behave
coherently similarly to (but not exactly as) boson condensates.

The BCS Hamiltonian of the system of electrons is given by Schrieffer (2018);
Grosso and Parravicini (2013):

$\hat{H}_{BCS}=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}}\hat{c}_{\mathbf{k}\sigma}^{\dagger}\hat{c}_{\mathbf{k}\sigma}+\sum_{\mathbf{k}\mathbf{k}^{\prime}}U_{\mathbf{k}\mathbf{k}^{\prime}}\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}\hat{c}_{-\mathbf{k}^{\prime}\downarrow}\hat{c}_{\mathbf{k}^{\prime}\uparrow}$
(1)

where $\xi_{\mathbf{k}}=\hbar^{2}k^{2}/2m-\epsilon_{F}$ is the energy measured
with respect to the Fermi energy $\epsilon_{F}$, and
$U_{\mathbf{k}\mathbf{k}^{\prime}}$ are the matrix elements of the interaction
potential. The first term of the Hamiltonian is the kinetic energy while the
second potential term couples pairs of different momenta $k$ and $k^{\prime}$.
The ground state of a superconductor is given by the BCS ground state:

$\ket{\psi_{S}}=\prod_{\mathbf{k}}(u_{\mathbf{k}}+e^{i\phi}v_{\mathbf{k}}\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger})\ket{0}\,,$
(2)

where $\ket{0}$ is the vacuum state of the electrons and the operators
$\hat{c}^{\dagger}_{\mathbf{k}\sigma}$ ($\hat{c}_{\mathbf{k}\sigma}$) create
(destroy) an electron of momentum $k$ and spin $\sigma$. The real coefficients
$u_{\mathbf{k}}$ and $v_{\mathbf{k}}$ satisfy the normalization conditions
$u_{\mathbf{k}}^{2}+v_{\mathbf{k}}^{2}=1$. $v_{\mathbf{k}}^{2}$
($u_{\mathbf{k}}^{2}$) gives the probability that the Cooper pair of momentum
$k$ is occupied (unoccupied).

Since finding the excited states of a superconductor in terms of electron
operators is not easy, it is convenient to perform a Bogoliubov
transformation. Such a transformation diagonalizes the Hamiltonian (1) by
introducing new canonical fermionic operators
$\hat{\gamma}_{\mathbf{k}\sigma}$ (see appendix A):

$\hat{H}_{B}=\sum_{\mathbf{k}\sigma}E_{\mathbf{k}}\hat{\gamma}^{\dagger}_{\mathbf{k}\sigma}\hat{\gamma}_{\mathbf{k}\sigma}\,,$
(3)

where $E_{\mathbf{k}}$ is the energy of a an excited state associated to
momentum $k$:

$E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+\Delta_{\mathbf{k}}^{2}}\,,$ (4)

and $\Delta_{\mathbf{k}}$ are the so called superconducting gap parameters.
The operators $\hat{\gamma}_{\mathbf{k}\sigma}$ are such that the BCS ground
state is their vacuum state, i.e.
$\hat{\gamma}_{\mathbf{k}\sigma}\ket{\Psi_{S}}=0$. Acting with
$\hat{\gamma}_{\mathbf{k}\sigma}^{\dagger}$ on the BCS ground state gives an
excited state, called quasiparticle. The excited states in terms of
$\hat{c}^{\dagger}_{\mathbf{k}\sigma}$ read:

$\displaystyle\ket{\psi_{1}}=\hat{\gamma}_{\mathbf{k}\uparrow}^{\dagger}\ket{\psi_{S}}=\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\prod_{\mathbf{l}\neq\mathbf{k}}(u_{\mathbf{l}}+v_{\mathbf{l}}\hat{c}_{\mathbf{l}\uparrow}^{\dagger}\hat{c}_{-\mathbf{l}\downarrow}^{\dagger})\ket{0}\,$
(5)
$\displaystyle\ket{\psi_{2}}=\hat{\gamma}_{-\mathbf{k}\downarrow}^{\dagger}\ket{\psi_{S}}=\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}\prod_{\mathbf{l}\neq\mathbf{k}}(u_{\mathbf{l}}+v_{\mathbf{l}}c_{\mathbf{l}\uparrow}^{\dagger}c_{-\mathbf{l}\downarrow}^{\dagger})\ket{0}\,,$
(6)

which means that for given momentum $k$ there is an electron with probability
$1$ and the other state of the pair is empty. Quasiparticles can be then
interpreted as fermions created by $\hat{\gamma}^{\dagger}_{\mathbf{k}\sigma}$
which are in one-to-one correspondence with the
$\hat{c}^{\dagger}_{\mathbf{k}\sigma}$.

The coefficients $u_{\mathbf{k}}$ and $v_{\mathbf{k}}$ are equal to:

$\displaystyle
u_{\mathbf{k}}^{2}=\frac{1}{2}\bigg{(}1+\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\bigg{)}$
$\displaystyle
v_{\mathbf{k}}^{2}=\frac{1}{2}\bigg{(}1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\bigg{)}\,.$
(7)

One can see that, as $k$ varies from well below the Fermi surface to well
above it, $v_{\mathbf{k}}^{2}$ goes from $1$ to $0$ (and analogously
$u_{\mathbf{k}}^{2}$ goes from $0$ to $1$), i.e. Cooper pairs of momentum $k$
well below the Fermi surface are occupied in the ground state with probability
$1$. As the momentum of the Cooper pairs increases above the Fermi surface,
the occupation probability decreases to $0$.

In general, the superconducting gap parameters $\Delta_{\mathbf{k}}$ have
different values for different $k$. As well known, BCS simplifies this (see
appendix A) by setting $\Delta_{\mathbf{k}}=\Delta$ for values of $k$ such
that $|\xi_{\mathbf{k}}|<\hbar\omega_{D}$ (with $\omega_{D}$ the Debye
frequency of phonons), and $\Delta_{\mathbf{k}}=0$ otherwise. $\Delta$ is the
superconducting gap, and its meaning becomes clearer by looking at the
quasiparticle energies $E_{\mathbf{k}}$ as function of $\xi_{\mathbf{k}}$ in
Eq.(4). They have a minimum at $\xi_{\mathbf{k}}=0$ and the value of the
excitation energy at this point is $\Delta$. This implies that to have an
excited states at least an amount of energy $\Delta$ is required.

The superconducting gap is in general a function of temperature, but at
sufficiently low temperatures (we will make this assumption hereafter) it is
$\Delta(T)\approx\Delta(0)=1.76k_{B}T_{c}$, where $T_{c}$ is the critical
temperature and $k_{B}$ is the Boltzmann constant.

In the realistic case of a small, but non vanishing $T$, some quasiparticle
will be thermally excited, and they will be distributed according to some
occupation function $f(E,T)$. For a superconductor at thermal equilibrium the
quasiparticles have a Fermi-Dirac occupation function:

$f_{FD}(E,T)=\frac{1}{e^{E/k_{B}T}+1}\,.$ (8)

The normalized quasiparticle density $x_{qp}$, i.e. the ratio between the
number of quasiparticles and the number of Cooper pairs inside a
superconductor, gives an estimate of how many quasiparticles there are in a
superconducting device. Its expression is given by:

$x_{qp}=\int_{\Delta}^{\infty}f(E)\rho(E)dE\,,$ (9)

where $\rho(E)$ is the normalized superconducting density of states:

$\rho(E)=\frac{E}{\sqrt{E^{2}-\Delta^{2}}}\,.$ (10)

Making the assumptions that $f(E)$ is a Fermi-Dirac distribution and that
$x_{qp}$ is small, so that most of the states are close to the gap, one can
find the following expression for $x_{qp}$ Catelani  _et al._ (2011b):

$x_{qp}=\sqrt{2\pi k_{B}T/\Delta}e^{-\Delta/k_{B}T}\,.$ (11)

Inserting in this equation temperatures close to the experimental regime of
$20$ mK, and the parameters of aluminum (a typical superconductor used in
transmon qubits), one can see that the quasiparticle density should be
exponentially suppressed: $x_{qp}\sim 10^{-52}$ de Visser  _et al._ (2014).
For the purposes of this paper, it is very important to stress that however
experiments on superconducting qubits and superconducting resonators show
higher density values, $x^{\text{\tiny}{exp}}_{qp}\sim 10^{-9}-10^{-6}$
Serniak  _et al._ (2019). Thermal equilibrium seems not able to explain such
an excess of quasiparticles and for this reason the latter are usually called
non-equilibrium quasiparticles. Appendix A contains further technical details
about the BCS theory.

##  III CSL model in the BCS framework

The CSL model is an alternative formulation of quantum mechanics devised to
solve the problem of the quantum-to-classical transition in quantum theory.
CSL unifies the Schrödinger evolution, which is linear and deterministic, with
the non-linear and stochastic dynamics giving wave packet reduction. This is
done by modifying the Schrödinger equation, adding stochastic and non-linear
terms that implement the collapse of the wave function. This dynamical
modification is consistent with quantum mechanics in the microscopic regime,
where the standard theory gives extremely accurate predictions, at the same
time justifying why macroscopic objects are always localized in space. The
strength and spatial extension of the collapse (which is white in time) are
dictated respectively by two parameters of the theory, $\lambda$ and $r_{c}$.
Theoretical arguments suggest that $r_{c}\sim 10^{-7}$ m Bassi and Ghirardi
(2003); Carlesso and Donadi (2019) and $\lambda$ in the range $\sim
10^{-8}-10^{-16}$s-1Bassi and Ghirardi (2003); Adler (2007). From now on, to
fix a possible value we set $\lambda=10^{-10}$s-1.

The full non-linear and stochastic CSL dynamics is not easy to work with.
Since we are interested in prediction of observable effects, one can
equivalently use the following simplified linear dynamics (see appendix B for
further details):

$i\hbar\frac{d\ket{\psi(t)}}{dt}=(\hat{H}+\hat{H}_{\text{\tiny
CSL}})\ket{\psi(t)}\,,$ (12)

where $\hat{H}$ is the Hamiltonian of the system ($\hat{H}_{BCS}$ in the
present case) and $\hat{H}_{\text{\tiny CSL}}$ is the CSL contribution. By
performing a Fourier transform to momentum space in a normalization box of
volume $V$ to avoid potential divergences (see appendix B for the details),
$\hat{H}_{\text{\tiny CSL}}$ takes the following expression:

$\begin{split}&\hat{H}_{\text{\tiny
CSL}}=-\frac{m\hbar\sqrt{\lambda}}{m_{0}V}\sum_{\mathbf{k_{1}}\mathbf{k_{2}},s}\widetilde{W}_{\mathbf{k_{1}}-\mathbf{k_{2}}}(t)\,\widetilde{G}_{\mathbf{k_{1}}-\mathbf{k_{2}}}\,\hat{c}^{\dagger}_{\mathbf{k}_{1}s}\hat{c}_{\mathbf{k}_{2}s}\,,\end{split}$
(13)

where $m_{0}$ is the nucleon mass and $m$ is the mass of the particle
considered (in our case, the electrons). The stochastic processes
$\widetilde{W}_{\mathbf{k}}(t)$ have expectation value and two-point
correlator given by:

$\displaystyle\mathbb{E}[\widetilde{W}_{\mathbf{k_{1}}}(t)]=0$ (14)
$\displaystyle\mathbb{E}[\widetilde{W}_{\mathbf{k_{1}}}(t)\widetilde{W}_{\mathbf{k_{2}}}(s)]=V\delta_{(\mathbf{k_{1}}+\mathbf{k_{2}})}\delta(t-s)\,,$
(15)

and they are weighted by the Gaussian function:

$\widetilde{G}_{\mathbf{k}}=(4\pi
r_{c}^{2})^{3/4}e^{-\frac{r_{c}^{2}k^{2}}{2}}\,.$ (16)

Equation (13) shows that the CSL noise scatters electrons, effectively acting
as a kick which adds energy to the system. The goal of the present analysis is
to quantify the main effect of those kicks on the BCS ground state of a
superconductor, which, as we will see, amounts to the breaking of Cooper pairs
into quasiparticles. More precisely, we will compute the transition
probability from a superconductor initially in its ground state to states
containing quasiparticles. We will use perturbation theory, thus considering
the CSL Hamiltonian (13) as an interacting term with the diagonal Bogoliubov
Hamiltonian $\hat{H}_{B}$ as the free term.

In this setup, the transition probability from the initial state $\ket{i}$ to
a final state $\ket{f}$ is given by:

$P_{fi}=\mathbb{E}[|T_{fi}|^{2}]=\mathbb{E}[|\bra{f}\hat{U}_{I}(t,t_{i})\ket{i}|^{2}]\,,$
(17)

where $T_{fi}$ is the transition amplitude from $\ket{i}$ to a state $\ket{f}$
and $\hat{U}_{I}(t,t_{i})$ is the time evolution operator in the interaction
picture:

$\hat{U}_{I}(t,t_{i})=e^{\frac{i}{\hbar}\hat{H}_{B}t}\hat{U}(t,t_{i})e^{-\frac{i}{\hbar}\hat{H}_{B}t}\,.$
(18)

The operator $\hat{U}_{I}(t,t_{i})$ can be expanded in the Dyson series, which
at first order reads:

$\hat{U}(t,t_{i})=1+\int_{t_{i}}^{t}ds\,\hat{H}^{(I)}_{\text{\tiny
CSL}}(s)\,.$ (19)

where $\hat{H}^{(I)}_{\text{\tiny
CSL}}(t)=e^{\frac{i}{\hbar}\hat{H}_{B}t}\hat{H}_{\text{\tiny
CSL}}e^{-\frac{i}{\hbar}\hat{H}_{B}t}$ is the Hamiltonian (13) in the
interaction picture.

The diagonal form of the Bogoliubov Hamiltonian in terms of
$\hat{\gamma}_{\mathbf{k}\sigma}$ and the properties of the BCS ground state
for these operators make it easier to work with creation and annihilation
operators for quasiparticles $\hat{\gamma}^{\dagger}_{\mathbf{k}\sigma}$ and
$\hat{\gamma}_{\mathbf{k}\sigma}$. Using the inverse Bogoliubov transformation
[see Eq.(47) in appendix A], we can transform the CSL Hamiltonian in terms of
the operators $\hat{\gamma}_{\mathbf{k}s}$ and then find its expression in the
interaction picture, reading:

$\begin{split}&\hat{H}^{(I)}_{\text{\tiny
CSL}}=-\frac{\hbar\sqrt{\lambda}m}{m_{0}V}\sum_{\mathbf{k_{1}}\mathbf{k_{2}}}\widetilde{G}_{\mathbf{k_{1}}-\mathbf{k_{2}}}\times\\\
&\bigg{[}\widetilde{W}_{\mathbf{k_{1}}-\mathbf{k_{2}}}(t)L(k_{1},k_{2})e^{\frac{i}{\hbar}(E_{\mathbf{k}_{1}}-E_{\mathbf{k}_{2}})t}\hat{\gamma}^{\dagger}_{\mathbf{k}_{1}\uparrow}\hat{\gamma}_{\mathbf{k}_{2}\uparrow}\\\
&+\widetilde{W}_{\mathbf{k_{2}}-\mathbf{k_{1}}}(t)L(k_{1},k_{2})e^{\frac{i}{\hbar}(E_{\mathbf{k}_{2}}-E_{\mathbf{k}_{1}})t}\hat{\gamma}^{\dagger}_{-\mathbf{k}_{2}\downarrow}\hat{\gamma}_{-\mathbf{k}_{1}\downarrow}\\\
&+\widetilde{W}_{\mathbf{k_{2}}-\mathbf{k_{1}}}(t)M(k_{1},k_{2})e^{-i\phi}e^{-\frac{i}{\hbar}(E_{\mathbf{k}_{1}}+E_{\mathbf{k}_{2}})t}\hat{\gamma}_{\mathbf{k}_{2}\uparrow}\hat{\gamma}_{-\mathbf{k}_{1}\downarrow}\\\
&+\widetilde{W}_{\mathbf{k_{1}}-\mathbf{k_{2}}}(t)M(k_{1},k_{2})e^{i\phi}e^{\frac{i}{\hbar}(E_{\mathbf{k}_{1}}+E_{\mathbf{k}_{2}})t}\hat{\gamma}^{\dagger}_{\mathbf{k}_{1}\uparrow}\hat{\gamma}^{\dagger}_{-\mathbf{k}_{2}\downarrow}\bigg{]}\,.\end{split}$
(20)

The first two terms in the square brackets are associated to quasiparticle
scattering: a quasiparticle of momentum $k_{1}$ is annihilated and another one
of momentum $k_{2}$ is created. The third term is associated to quasiparticle
recombination: two quasiparticles of different momenta are annihilated. The
fourth term is the inverse process, called quasiparticle generation: two
quasiparticle of different momenta are created. The functions $M$ and $L$,
usually called coherence factors, are given by:

$\displaystyle
L(k_{1},k_{2})=(u_{\mathbf{k}_{1}}u_{\mathbf{k}_{2}}-v_{\mathbf{k}_{1}}v_{\mathbf{k}_{2}})$
(21) $\displaystyle
M(k_{1},k_{2})=(u_{\mathbf{k}_{1}}v_{\mathbf{k}_{2}}+v_{\mathbf{k}_{1}}u_{\mathbf{k}_{2}})\,.$
(22)

Substituting Eqs. (7) we have that:

$\displaystyle
L^{2}(E_{1},E_{2})=\frac{1}{2}\bigg{(}1-\frac{\Delta^{2}-\xi_{1}\xi_{2}}{E_{1}E_{2}}\bigg{)}$
(23) $\displaystyle
M^{2}(E_{1},E_{2})=\frac{1}{2}\bigg{(}1+\frac{\Delta^{2}-\xi_{1}\xi_{2}}{E_{1}E_{2}}\bigg{)}\,,$
(24)

where according to Eq.(4),
$\xi_{\mathbf{k}}=\sqrt{E_{\mathbf{k}}^{2}-\Delta_{\mathbf{k}}^{2}}$.

##  IV CSL reduction on transmon qubits

As anticipated, the direct effect of CSL on superconducting quantum computers
is to destroy superposition states. We now give an estimate of the rate at
which CSL reduces superpositions of the computational basis states (labeled by
$\ket{0}$ and $\ket{1}$) in a transmon qubit. To estimate this rate, we first
need a description of the physical form of the computational basis states of
transmon qubits, keeping in mind that the CSL collapse mechanism is sensitive
only to space superpositions of different masses.

A transmon qubit can be described as an electrical circuit in which a
superconducting island (usually made of aluminum) of volume $V\sim
10^{2}\mu$m3 is linked to a superconducting reservoir (of approximately the
same size) through an insulating barrier of width $d\sim 1-10^{2}$ nm Wang
_et al._ (2014), forming a Josephson junction. By applying a gate voltage,
Cooper pairs can tunnel through the Josephson junction from the reservoir into
the island. The dimensions and the components of the transmon circuit are such
that the number of excess Cooper pairs in the island becomes a quantum number,
and the tunneling of a single Cooper pair can be controlled by manipulating
the gate voltage.

The computational basis states of the transmon qubit are characterized by the
number of excess Cooper pairs that have tunneled from the reservoir into the
island. The transmon qubit works in a regime for which $\ket{0}$ and $\ket{1}$
are not exactly eigenstates of the number of extra Cooper pairs operator, but
rather both $\ket{0}$ and $\ket{1}$ are a superposition of states with a
different number of extra Cooper pairs on the island. Nevertheless, for our
estimate we can assume that, for a typical device, the difference in the
number of pairs of the two computational basis states is of the order of $4$
Koch  _et al._ (2007). For our purposes, we can then effectively identify
$\ket{0}$ as the state with $4$ Cooper pairs on one side of the Josephson
junction, the reservoir, and $\ket{1}$ as the state with $4$ Cooper pairs on
the other side of the junction, the island. We can think of these four Cooper
pairs as two groups of four electrons. The two groups are separated by the BCS
coherence length $\xi_{c}$ (a measure of the average distance between the two
electrons in a Cooper pair) which for aluminum, the superconducting material
used inside transmon circuits, is $\sim 10^{-6}$ m Grosso and Parravicini
(2013). Given these assumptions, we can compute the reduction rate with the
formula Adler (2007):

$\Gamma_{R}=\lambda n^{2}N\bigg{(}\frac{m_{e}}{m_{0}}\bigg{)}^{2}$ (25)

for $n$ particles within a radius smaller than the correlation length $r_{c}$,
$N$ groups of particles separated by more than the correlation length $r_{c}$,
and with $m_{e}$ the mass of electrons. In our effective model, if we
substitute $n=4$ and $N=2$ into Eq.(25), considering an optimal configuration
with four electrons packed within a distance smaller than $r_{c}$ in each of
the $N=2$ groups of electrons separated by a distance $\xi_{c}>r_{c}$, we find
the reduction rate $\Gamma=32\lambda(m_{e}/m_{0})^{2}\approx 10^{-16}$s-1. We
chose the value for $n$ that gives the highest value for $\Gamma_{R}$, to
estimate the strongest theoretical CSL effect. This shows that CSL reduction
is negligibly weak, as the lifetime of a single qubit would be of the orders
of billions of years. Even a large quantum computer composed by millions or
even billions of such qubits, would be safe against the localization of
superpositions dictated by collapse models.

However, besides collapsing superposition states, the CSL noise induces also
dissipation as any other environmental noise: it perturbs the transmon qubits
materials, also leading to decoherence effects. In fact, dissipation generates
excited quasiparticle states that accumulate over time and their presence
destroys quantum superpositions of transmon qubits.

##  V Transition probability due to CSL dissipation

CSL is ineffective in directly suppressing superposition states of
superconducting qubits, mainly because too few electrons are involved in the
superposition, which moreover have a very small mass. However, CSL impacts
superpositions of superconducting qubits also indirectly. As discussed in
section III, the CSL noise couples to the Cooper pairs inside a
superconductor, generating its excited states, the quasiparticles. Since the
rate of decoherence of transmon qubits is proportional to the density of
quasiparticles Catelani  _et al._ (2011b), these limit their performance. We
now quantify the quasiparticle density due to CSL, to infer a limit on the
coherence time of transmon qubits. To do so we compute the CSL transition
probability to quasiparticle states, that will be used in the next section to
obtain the CSL generation rate of quasiparticles and the evolution of the
occupation function of quasiparticles $f(E)$.

We start by computing the transition probability in Eq.(17) with
$\ket{i}=\ket{\psi_{S}}$, using the time evolution operator expanded at first
order as in Eq.(19), with the CSL Hamiltonian in Eq.(20), expressed in terms
of quasiparticle operators. By recalling that $\ket{\psi_{S}}$ is the vacuum
state for the operators $\hat{\gamma}_{\mathbf{k}s}$, only the fourth term of
the Hamiltonian in Eq.(20) gives a non vanishing contribution when acting on
$\ket{\psi_{S}}$. This produces a transition to a final state different from
$\ket{\psi_{S}}$, which contains quasiparticles. From these considerations,
the main effect of the CSL noise on superconductors at first order in
perturbation theory is the generation of quasiparticles. In particular, the
only final states which give non zero contribution to the expectation value in
(17) are those of the form:

$\ket{f}=\hat{\gamma}^{\dagger}_{\mathbf{q}\uparrow}\hat{\gamma}^{\dagger}_{-\mathbf{p}\downarrow}\ket{\psi_{S}}\,,$
(26)

where $q$ and $p$ are fixed. Choosing $\ket{f}$ as in Eq.(26), the zero order
term of the Dyson series gives a zero contribution when inserted into Eq.(17).
The first order contribution of the Dyson series is:

$T^{(1)}_{qp}=-\frac{i}{\hbar}\int_{t_{i}}^{t}dt_{1}\bra{\psi_{S}}\hat{\gamma}_{-\mathbf{p}\downarrow}\hat{\gamma}_{\mathbf{q}\uparrow}\hat{H}^{(I)}_{\text{\tiny
CSL}}(t_{1})\ket{\psi_{S}}\,,$ (27)

and by substituting Eq.(20) for the interacting Hamiltonian we find:

$\begin{split}T^{(1)}_{qp}=&\frac{i\sqrt{\lambda}m}{m_{0}V}\int_{t_{i}}^{t}dt_{1}\,\sum_{\mathbf{k_{1}}\mathbf{k_{2}}}\widetilde{G}_{\mathbf{k_{1}}-\mathbf{k_{2}}}\widetilde{W}_{\mathbf{k_{1}}-\mathbf{k_{2}}}(t_{1})M(k_{1},k_{2})\times\\\
&e^{i\phi}e^{\frac{i}{\hbar}(E_{\mathbf{k}_{1}}+E_{\mathbf{k}_{2}})t_{1}}\bra{\psi_{S}}\hat{\gamma}_{-p\downarrow}\hat{\gamma}_{q\uparrow}\hat{\gamma}^{\dagger}_{\mathbf{k}_{1}\uparrow}\hat{\gamma}^{\dagger}_{-\mathbf{k}_{2}\downarrow}\ket{\psi_{S}}\,.\end{split}$
(28)

The expectation value in the last line is simplified using the anticommutation
rules:

$\bra{\psi_{S}}\hat{\gamma}_{-p\downarrow}\hat{\gamma}_{q\uparrow}\hat{\gamma}^{\dagger}_{\mathbf{k}_{1}\uparrow}\hat{\gamma}^{\dagger}_{-\mathbf{k}_{2}\downarrow}\ket{\psi_{S}}=\delta_{\mathbf{k_{1}},\mathbf{q}}\delta_{\mathbf{k_{2}},\mathbf{p}}\,,$
(29)

thus leading to:

$T^{(1)}_{qp}=\frac{i\sqrt{\gamma}m}{m_{0}V}\int_{t_{i}}^{t}dt_{1}\,\widetilde{G}_{\mathbf{q}-\mathbf{p}}\widetilde{W}_{\mathbf{q}-\mathbf{p}}(t_{1})M(q,p)e^{i\phi}e^{\frac{i}{\hbar}(E_{q}+E_{p})t_{1}}\,.$
(30)

This is the transition amplitude to a specific final state with fixed momenta
$p$ and $q$. The transition probability from the BCS ground state to a state
as in Eq.(26) is computed by taking the expectation value of the square
modulus of the transition amplitude, according to Eq.(17). In doing so, one
has to compute the two-point correlator between the noise and its complex
conjugate. Using the fact that
$\widetilde{W}^{*}_{\mathbf{k}}(t)=\widetilde{W}_{-\mathbf{k}}(t)$ and
Eq.(15), one has that:

$P_{qp}=\frac{\lambda
m^{2}}{m_{0}^{2}V}\widetilde{G}^{2}_{\mathbf{q}-\mathbf{p}}\,M^{2}(q,p)\,t\,,$
(31)

where we assumed that $t_{0}=0$ and we performed integration over $t_{1}$.
This probability grows linearly with time. By dividing Eq.(31) by $t$, we
obtain the transition rate, i.e. the rate at with which two quasiparticles
with given momentum $p$ and $q$ are generated by the CSL noise. This
transition rate will be the starting point to compute the generation rate in
the next section, through which we can obtain the evolution of the
quasiparticle occupation function $f(E)$, and thus the quasiparticle density.
We present further calculations of the total generation rate of quasiparticles
per unit volume in appendix C.

##  VI Quasiparticle density due the CSL noise

In the following, we consider a superconductor at thermal equilibrium and we
neglect all sources of environmental noise. To account for thermal effects, we
should consider the interaction of electrons with other electrons and with
phonons, but since typically the electron-phonon interaction is dominant, we
will consider only this one. The electron-phonon interaction comprises three
main physical processes: quasiparticle scattering (both by emission and
absorption of a phonon), quasiparticle recombination (by the emission of a
phonon) and quasiparticle generation (by the absorption of a phonon).

The kinetic equation for the quasiparticle occupation function describes how
$f(E)$ redistributes over time because of the above processes. It contains the
rates of the electron-phonon processes, and a generation rate given by an
external source Martinis  _et al._ (2009); Goldie and Withington (2012):

$\begin{split}&\frac{df(E)}{dt}=\gamma_{g}^{\text{\tiny
ext}}(E)+\frac{\gamma_{0}}{\Delta^{3}}\int_{E}^{\infty}dE^{\prime}S(E,E^{\prime})\times\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,[(\bar{f}(E))f(E^{\prime})(N(E^{\prime}-E)+1)\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(E)\bar{f}(E^{\prime})N(E^{\prime}-E)]\\\
&+\frac{\gamma_{0}}{\Delta^{3}}\int_{\Delta}^{E}dE^{\prime}S(E,E^{\prime})[(\bar{f}(E))f(E^{\prime})N(E-E^{\prime})\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(E)\bar{f}(E^{\prime})(N(E-E^{\prime})+1)]\\\
&+\frac{\gamma_{0}}{\Delta^{3}}\int_{\Delta}^{\infty}dE^{\prime}G(E,E^{\prime})[(\bar{f}(E))\bar{f}(E^{\prime})N(E+E^{\prime})\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-f(E)f(E^{\prime})(N(E+E^{\prime})+1)]\,,\end{split}$
(32)

where $S(E,E^{\prime})=(E-E^{\prime})^{2}\rho(E^{\prime})L^{2}(E,E^{\prime})$,
$G(E,E^{\prime})=(E+E^{\prime})^{2}\rho(E^{\prime})M^{2}(E,E^{\prime})$ and
the rate $\gamma_{0}$ is a characteristic electron-phonon rate and it is a
constant for a given material (for aluminum $1/\gamma_{0}=\tau_{0}=438$ns
Martinis  _et al._ (2009)). The factor $N(\Omega)$ is the occupation function
of phonons, which is taken to follow a Bose-Einstein distribution:

$N(\Omega)=\frac{1}{e^{\Omega/k_{B}T_{ph}}-1}\,,$ (33)

for a bath of phonons at temperature $T_{ph}$. Phonons are supposed to be in
equilibrium at the refrigerator temperature and thus $N(\Omega)$ does not
change in time. The rate $\gamma_{g}^{\text{\tiny ext}}(E)$ is a generation
rate per unit time due to external sources. The second and third terms in
Eq.(32) are associated to thermal quasiparticle scattering, and the fourth
term to thermal quasiparticle recombination and generation.

One can check that when quasiparticles and phonons are in thermal equilibrium,
i.e. when $f(E)$ is a Fermi-Dirac distribution and $N(\Omega)$ is a Bose-
Einstein distribution with $T_{ph}=T$, the terms in every square brackets of
Eq.(32) cancel. This is not surprising since, if $\gamma_{g}^{\text{\tiny
ext}}=0$, i.e. there is no source of quasiparticles, the superconductor
remains at equilibrium and the occupation function $f(E)$ does not change in
time.

Similarly to the phonon case, the CSL Hamiltonian (20) contains terms
associated to scattering, recombination and generation of quasiparticles
because of the interaction with the CSL noise. For temperatures below $100$mK,
we checked that the generation rate is the only significant process. Let us
then compute this rate when two quasiparticles, one at a fixed energy $E$ and
the other at any energy $E^{\prime}$, are created by the interaction with the
CSL noise.

To do so, we start from the transition probability in Eq.(31), integrating
over one momentum:

$\gamma_{g}^{\text{\tiny CSL}}(E_{\mathbf{q}})=\frac{\lambda(4\pi
r_{c})^{3/2}m^{2}}{m_{0}^{2}(2\pi)^{3}}\int
d^{3}pe^{-r_{c}^{2}(p-q)^{2}}M^{2}(p,q)\bar{f}(E_{\mathbf{p}})\,.$ (34)

Defining $x=E/\Delta$ and following the calculations in appendix D we find:

$\begin{split}&\gamma_{g}^{\text{\tiny CSL}}(x)=\frac{m^{2}\lambda
r_{c}}{2\sqrt{\pi}m_{0}^{2}}\frac{\sqrt{2m\Delta}}{\hbar}\frac{1}{\sqrt{s(x)+\beta}}e^{-\frac{\Delta}{k_{B}T_{\text{\tiny
CSL}}}s(x)}\\\ &e^{-\frac{2T_{F}}{T_{\text{\tiny
CSL}}}}\int_{1}^{\infty}dy\,e^{-\frac{\Delta}{k_{B}T_{\text{\tiny
CSL}}}s(y)}\,e^{\frac{2\Delta}{k_{B}T_{\text{\tiny
CSL}}}(\sqrt{(s(x)+\beta)(s(y)+\beta)}}\\\
&\times\rho(y)\frac{1}{2}\bigg{(}1-\frac{\sqrt{x^{2}-1}\sqrt{y^{2}-1}}{xy}+\frac{1}{xy}\bigg{)}\bar{f}(\Delta
y)\,,\end{split}$ (35)

where $s(x)=\sqrt{x^{2}-1}$, $\beta=\epsilon_{F}/\Delta$ and
$k_{B}T_{\text{\tiny CSL}}=\hbar^{2}/(2mr_{c}^{2})$.

We solved numerically Eq.(32) with $\gamma_{g}^{\text{\tiny
ext}}=\gamma_{g}^{\text{\tiny CSL}}$ and $\lambda=10^{-10}$s-1. Figure 1 shows
the evolved occupation function obtained numerically (blue solid line) for a
starting equilibrium temperature of $20$ mK, which is a typical operational
temperature of transmon qubits. The steady state deviates strongly from the
starting Fermi Dirac distribution, given by Eq.(8) with $T=20$ mK. The evolved
occupation function can not be approximated by a thermal distribution at some
effective temperature. More on the evolution of the quasiparticle occupation
function can be found in appendix E.

The quasiparticle density $x^{\text{\tiny}{CSL}}_{qp}$ generated by CSL is
obtained easily by using Eq.(9) with $f(E)$ given by the evolved quasiparticle
occupation function. $x^{\text{\tiny}{CSL}}_{qp}$ turns out to be $\sim
10^{-18}$. This quasiparticle density is $\sim 9$ orders of magnitude lower
than the lowest reported experimental value of $x_{qp}=10^{-9}$ Serniak  _et
al._ (2019). This means that other environmental noise sources are yet
dominant with the current techological implmentation of the transmon qubits.
Nevertheless we can use the CSL quasiparticle density to get a limit on the
coherence time of transmon qubits.

Figure 1: Log plot of the quasiparticle occupation function obtained by
solving Eq.(32) numerically for a starting equilibrium temperature of $20$ mK.
The evolved occupation function deviates significantly from the initial
thermal distribution $f_{FD}$, here not shown because too small.

##  VII Possibility of testing CSL with superconducting devices

Our analysis not only leads to possible limits on the performances of
superconducting devices imposed by CSL, but also suggests that CSL dissipation
could be tested with such devices.

We compare the electron-phonon generation and recombination rates with the CSL
generation rate computed in section VI [see Eq.(35)], to find the temperatures
for which the CSL generation rate dominates over the electron-phonon
generation and/or recombination rates.

The recombination and generation rates enter the last term of Eq.(32). The
electron-phonon recombination rate $\gamma_{r}^{e-ph}(E)$ is the inverse
lifetime for a quasiparticle at some given energy $E$ to recombine with
another quasiparticle of any energy $E^{\prime}$ by emitting a phonon of
energy $E+E^{\prime}$. Its expression is given byChang and Scalapino (1977);
Martinis  _et al._ (2009):

$\begin{split}\gamma_{r}^{e-ph}(E)=\frac{\gamma_{0}}{\Delta^{3}}\int_{\Delta}^{\infty}dE^{\prime}&(E+E^{\prime})^{2}\rho(E^{\prime})\bigg{(}1+\frac{\Delta^{2}}{EE^{\prime}}\bigg{)}\\\
&(N(E+E^{\prime})+1)f(E^{\prime})\,.\end{split}$ (36)

We now briefly outline the different contributions appearing in the above
integral. The term $(E+E^{\prime})^{2}$ is the square energy of the emitted
phonon during the recombination process. The normalized superconducting
density of states $\rho(E^{\prime})$ (see Eq.(10)) appears since we are
integrating over $E^{\prime}$. The term $(1+\Delta^{2}/(EE^{\prime}))$ is the
coherence factor squared $M^{2}(E,E^{\prime})$ in Eq.(21). The factor
$N(\Omega)$ is the occupation function of phonons, a Bose-Einstein
distribution Eq.(33) at fixed $T=T_{ph}$. The recombination rate is
proportional to the quasiparticle density because of the factor
$f(E^{\prime})$ in the integrand Martinis  _et al._ (2009). This implies that
recombination mediated by phonons is slow when the density of quasiparticles
is small.

The electron-phonon generation rate $\gamma_{g}^{e-ph}(E)$ is the inverse
lifetime for a quasiparticle at some given energy $E$ to be generated with
another quasiparticle of any energy $E^{\prime}$ by absorbing a phonon of
energy $E+E^{\prime}$ Chang and Scalapino (1977); Martinis  _et al._ (2009).
Its expression is given by:

$\begin{split}\gamma_{g}^{e-ph}(E)=\frac{\gamma_{0}}{\Delta^{3}}\int_{\Delta}^{\infty}dE^{\prime}&(E+E^{\prime})^{2}\rho(E^{\prime})\bigg{(}1+\frac{\Delta^{2}}{EE^{\prime}}\bigg{)}\\\
&\bar{f}(E^{\prime})N(E+E^{\prime})\,.\end{split}$ (37)

where $\bar{f}(E)=1-f(E)$.

In figure 2 we fix the energy at $\Delta$ and plot the two differences
$D_{1}=\gamma_{g}^{\text{\tiny CSL}}(\Delta,T)-\gamma_{g}^{e-ph}(\Delta,T)$
and $D_{2}=\gamma_{g}^{\text{\tiny
CSL}}(\Delta,T)-\gamma_{r}^{e-ph}(\Delta,T)$ as $T$ varies from $20$ mK to
$80$mK. We assume that quasiparticles and phonons are in equilibrium so that
$T_{ph}=T$ for every temperature $T$.

We see that the CSL generation rate is dominant over the electron-phonon
generation rate ($D_{1}>0$) for temperatures lower than $\sim 70$mK, and
similarly $D_{2}>0$ for temperatures lower than $\sim 35$mK. In particular, at
temperatures lower than $35$mK the CSL generation rate is many orders of
magnitude larger than the electron-phonon generation rate. Since a
superconducting sample can be easily cooled down to $\sim 20$mK, this regime
is already accessible by experiments. Of course, in order to perform a
meaningful test of CSL, one should isolate the system from any other external
source of quasiparticles. Note however that the CSL generation rate is itself
very small $\sim 10^{-18}$ s-1 (computed through Eq.(35) at $E=\Delta$), thus
making this kind of experiment challenging.

A straightforward way to perform such an experiment is by measuring the
quasiparticle subgap current through a Josephson junction, which as we saw is
a key component in superconducting qubits. When two identical superconductors
at absolute zero are linked together with an insulator, at voltages smaller
than $V<2\Delta/e$ no quasiparticle current flows, since no quasiparticle
state is excited. At finite temperature and $V<2\Delta/e$, there will be a
finite current of quasiparticles, called quasiparticle subgap current Milliken
_et al._ (2004). The formula for the subgap current at the voltage difference
$V$ is given by:

$\begin{split}I_{qp}&=\frac{1}{eR_{N}}\int
dE\,\rho(E)\rho(E+eV)(f(E)-f(E+eV))\\\ &=\frac{I_{c}}{\Delta}\int
dE\,\rho(E)\rho(E+eV)(f(E)-f(E+eV))\\\ &=I_{c}\int
dx\,\rho(x)\rho(x+eV/\Delta)(f(x)-f(x+eV/\Delta))\end{split}$ (38)

where $R_{N}$ and $I_{c}$ are, respectively, the normal state resistance and
the critical current of the junction. Note that we used the Ambegaokar-
Baratoff relation $R_{N}e=\pi\Delta/2I_{c}\approx\Delta/I_{c}$ Barone and
Paterno (1982), with $e$ the electron charge, in the second line, and in the
third line we performed the substitution $x=E/\Delta$.

Figure 2: Difference between the CSL generation rate and the electron-phonon
generation rate $D_{1}$ (blue solid line), and difference between the CSL
generation rate and the electron-phonon recombination rate $D_{2}$ (red dashed
line), as functions of temperature. $D_{1}$ is positive for temperatures lower
than $\sim 70$mK and $D_{2}$ is positive for temperatures lower than $\sim
35$mK. This implies that for temperatures lower than $\sim 35$mK the CSL
generation rate is dominant over both the electron-phonon generation and
recombination rates. This regime is accessible by experiments, since current
refrigerators reach temperatures down to $\sim 20$ mK. Isolating a
superconductor at these temperatures could lead to the detection of
quasiparticles generated by the CSL noise.

Plugging our computed occupation function in this equation, for a junction
with critical current $I_{c}\sim 10^{-4}A$, we find a quasiparticle current of
the order of $I^{\text{\tiny}{CSL}}_{qp}\sim 10^{-24}A$. This value has to be
compared to the experimental values which are of the order of
$I^{\text{\tiny}{exp}}_{qp}\sim 10^{-12}A$ Milliken  _et al._ (2004). A CSL
quasiparticle current $I^{\text{\tiny}{CSL}}_{qp}$ of such intensity is
extremely difficult to detect, as experiments that measure the quasiparticle
current have a sensitivity of the order of pA.

The value of the CSL quasiparticle subgap current depends on the dimension of
the Josephson junction and on the volumes of the two linked superconductors.
Superconducting qubits have small superconducting volumes and thin Josephson
junctions to work in the quantum regime and achieve a coherent control of the
tunneling of single Cooper pairs, but larger volumes imply a larger number of
quasiparticles.

Also, one can consider superconducting materials different from aluminum. The
critical temperature $T_{c}$ of the superconductor might play a role for two
reasons. First, a superconductor with lower $T_{c}$ may have a lower
temperature at which the subgap current saturates because of external sources
de Visser  _et al._ (2011); Milliken  _et al._ (2004). Second, superconductors
with higher $T_{c}$ increase the CSL quasiparticle generation as Eq.(69)
depends on $\Delta=1.76k_{B}T_{c}$.

It is worthwhile noticing that the perturbative approach developed in this
work can be easily extended to higher orders. Second order effects could be
relevant for CSL induced tunneling events of quasiparticles through Josephson
junctions.

##  VIII Limits set by CSL in superconducting quantum computers

We use the quasiparticle density obtained in section VI, to estimate the
ultimate limits set by CSL dissipation on the performances of quantum
computers based on transmon qubits Koch  _et al._ (2007). First, we estimate
the coherence time allowed by CSL for a single qubit. Next, we estimate how
this would limit the performance of a quantum computer made of $N$ of these
qubits.

State of the art transmon qubits have a relaxation time $T_{1}$ (the inverse
of the relaxation rate $\Gamma_{1}$) of the order of $10-100$ $\mu$s
Kjaergaard  _et al._ (2020). A future goal is to achieve greater $T_{1}$ in
order to have more reliable qubits, but this is experimentally challenging.
Relaxation can be driven by many loss channels, such as radiation losses
Córcoles  _et al._ (2011), dielectric losses Wang  _et al._ (2015), two-level
fluctuators in the junction materials Klimov  _et al._ (2018) and by the
excess of quasiparticles in the superconducting materials Catelani  _et al._
(2011a); Glazman and Catelani (2021); Martinis  _et al._ (2009). The
relaxation rate depends on every of these loss channel affecting the qubit
Vepsäläinen  _et al._ (2020).

For a transmon qubit, the contribution to the relaxation rate due to
quasiparticles depends linearly on the normalized quasiparticle density
$x_{qp}$ Vepsäläinen  _et al._ (2020):

$\Gamma_{1}=\sqrt{\frac{2\omega_{q}\Delta}{\pi^{2}\hbar}}x_{qp}$ (39)

where $\omega_{q}$ is the frequency of the given qubit (in Vepsäläinen  _et
al._ (2020) $\omega_{q}=2\pi\times 3.48$ GHz).

Given $x^{\text{\tiny}{CSL}}_{qp}\sim 10^{-18}$ that we found in the previous
section, we have $\Gamma^{\text{\tiny}{CSL}}_{1}\approx 10^{-6}$s-1. This is
$\sim 10$ orders of magnitude larger than the CSL reduction rate estimated in
section IV, showing that CSL dissipation is way more effective than the direct
collapse process in corrupting superpositions of transmon qubits basis states.
Note that $T^{\text{\tiny}{CSL}}_{1}$ is anyhow $10$ orders of magnitude
larger than the $T_{1}=100$ $\mu$s of current transmon qubits, implying that
CSL dissipation would not influence effectively the performance of a single
qubit.

By knowing the limit on the coherence time of a single transmon qubit, we can
estimate the limit on the performances of a quantum computer with $N$ of these
transmon qubits. Decoherence of each qubit accumulates during the operational
time of the quantum computer, eventually spoiling the quantum computation.

Before proceeding, we mention that schemes to recover from errors during a
quantum computation, so called quantum error correction schemes Gottesman
(2010), are planned to be implemented in future devices. However scaling a
quantum computer with an implemented error correction scheme is not an easy
task.

In the near term, so called Noisy Intermediate Scale Quantum (NISQ) Preskill
(2018) devices could become a viable tool, whose qubits are subject to noise,
without any implemented error correction scheme. We then consider the effect
of CSL on the performance of a NISQ transmon quantum computer composed of $N$
qubits, each with a decay rate $\Gamma^{\text{\tiny}{CSL}}_{1}$. Quantum
algorithms on such a device might require the storage of a maximally entangled
state of all these $N$ qubits Mooney  _et al._ (2019), whose total decoherence
rate $\Gamma_{tot}$ scales with $N$ and depends on the individual decay rate
$\Gamma_{1}$ of the qubits Palma  _et al._ (1996); Ischi  _et al._ (2005).

Since the volume of a typical transmon qubit is hundreds of microns cube
Martinis  _et al._ (2009); Wang  _et al._ (2014), we roughly and
conservatively estimate that the distance among such qubits is greater than
$r_{c}$. Recalling Eq.(25), one thus has that the CSL rate scales linearly
with $N$:
$\Gamma^{\text{\tiny}{CSL}}_{tot}=N\times\Gamma^{\text{\tiny}{CSL}}_{1}$.
Thus, the CSL limit on the decoherence time of a single qubit,
$T^{\text{\tiny}{CSL}}_{1}\approx 10^{6}$ s, gives a limit on
$T^{\text{\tiny}{CSL}}_{tot}=1/\Gamma^{\text{\tiny}{CSL}}_{tot}$, the
decoherence time of a quantum computer of $N$ qubits. In order to gain a good
fidelity of the output of the quantum computation, i.e. an output as close as
possible to the desired result, the operational time of the quantum algorithm
should be substantially smaller than the total decoherence time. The
operational time of a quantum algorithm can be naively defined as
$T_{op}=n_{g}\times t_{g}$, where $n_{g}$ is the number of quantum gates and
$t_{g}$ is the time to implement one of these gate operations. State of the
art devices have a $t_{g}\approx 10-100$ ns. Then, what one requires is that:

$T_{op}\ll T^{\text{\tiny}{CSL}}_{tot}\Rightarrow
n_{g}\ll\frac{T^{\text{\tiny}{CSL}}_{1}}{t_{g}\,N}$ (40)

Note the resemblance of this condition to the so called ”rule of thumb” metric
for the performance of a quantum computer found in Salm  _et al._ (2020), if
one defines $\epsilon=T^{\text{\tiny}{CSL}}_{1}/t_{g}$.

As outlined in Salm  _et al._ (2020), it is not easy to estimate how smaller
$T_{op}$ should be with respect to $T^{\text{\tiny}{CSL}}_{tot}$. Having a
$T_{op}\approx T^{\text{\tiny}{CSL}}_{tot}$ implies that the algorithm is
highly probable to fail Moll  _et al._ (2018). Note moreover that the naive
estimate we gave for $T_{op}$ is optimistic, as many factors can contribute to
increase it. The most important examples are the restricted connectivity
between qubits in a quantum chip and the possibility of implementing only the
native gates allowed by the hardware Salm  _et al._ (2020). The first factor
introduces additional swap gates so that two-qubit gates are performed only
between physically connected qubits, the second factor requires the
decomposition of the gates of the algorithm into those belonging to the native
gate set, adding extra gate operations. Given these considerations, we will
conservatively require that $T_{op}\approx 10^{-3}\times
T^{\text{\tiny}{CSL}}_{tot}$.

State-of-the-art transmon quantum computers with $N\approx 10^{2}$ would have
a total decoherence time induced by CSL of $T^{\text{\tiny}{CSL}}_{tot}\sim
10^{4}$ s. Thus the total operational time to reach a sufficiently accurate
result is $T_{op}\approx 10$s, and by having $t_{g}=10^{-7}$ s, this allows a
maximum of $n_{g}\approx 10^{8}$ operations. A $N=10^{3}$ transmon quantum
computer, targeted by IBM in 2023, allows for $n_{g}=10^{7}$. Many important
quantum algorithms such as molecular simulation, Shor algorithm for prime
number factorization Reiher  _et al._ (2017); Kutin (2006), which require
respectively at least $N\sim 10^{2}$, $n_{g}\sim 10^{14}$, $N\sim 10^{3}$,
$n_{g}\sim 10^{9}$ Childs  _et al._ (2018) will be corrupted by the CSL noise.
NISQ transmon quantum computers will scale up, eventually reaching the
milestone of $N=10^{6}$, which is thought to be the number of qubits necessary
to apply error correction schemes. Note that in this case the total time
allowed by CSL without error correction will decrease to
$T^{\text{\tiny}{CSL}}_{tot}\sim 10^{-3}$s.

Figure 3: Log-log plot of the number $n_{g}$ of allowed gates as a function of
the number $N$ of qubits (blue dotted line) as dictated by Eq.(40) with
$T_{1}^{\text{\tiny}{CSL}}=10^{6}$s and fixed $t_{g}=10^{-7}$s. The colored
area under this curve corresponds to algorithms that use resources (i.e.
number of qubits and number of gates) that allow to reach a good fidelity of
the output. The three red dot correspond to the resources necessary to
complete algorithms such as Shor factorization, molecular simulation and
derivative pricing.

This would allow for a maximum number of $n_{g}\approx 10^{4}$ gate
operations, which is very far for what is needed to complete quantum
algorithms for real life application Shor (1999); Childs  _et al._ (2018).

Figure 3 summarizes these results: we plot the number of allowed gates as a
function of the number of qubits (blue dotted line) as dictated by Eq.(40).
Quantum algorithms that exploit a number of qubits and a number of gates that
stay under this curve (colored blue area) can reach a good fidelity of the
output under the influence of CSL. The three red dots show the resources,
estimated in Childs  _et al._ (2018) and Chakrabarti  _et al._ (2021), in
order to complete important quantum algorithms such as prime number
factorization, molecular simulation and derivative pricing Chakrabarti  _et
al._ (2021). The points corresponding to Shor factorization and molecular
simulation lie outside of the colored area: this means that CSL may spoil
quantum computation with transmon NISQ quantum computers and also stresses the
need of scaling up the devices with the possibility of performing quantum
error correction or other error mitigation techniques. In fact even if the CSL
limit on the coherence time of a single qubit is of the order of $10^{6}$s,
which seems an extraordinary long time, the performances of a quantum computer
could be spoiled without any scheme to recover from errors.

We conclude by pointing out that the values we obtained should be taken as
rough estimates of the fundamental limitation imposed by collapse models on
the performances of transmon quantum computers. To get a more accurate result,
one needs to focus on a specific algorithm and calculate in detail the
resources, in terms of number of qubits and of quantum gates, needed to
complete it with a sufficiently good fidelity of the output. Indeed the number
of quantum gates may vary a lot depending on different aspects Salm  _et al._
(2020), such as the specific physical hardware used and the algorithm to be
solved.

##  IX Conclusions and Outlook

We showed how the CSL model affects superconducting quantum computers. The
intrinsic localization of superpositions dictated by collapse models leaves
the superposition of basis states of transmon qubits intact for very long
times. However CSL contributes to decoherence also in an indirect way:
dissipation induced by the CSL noise perturbs the superconducting material and
leads to the generation of quasiparticles. These accumulate over time inside
the volume of the device leading to relaxation at a rate proportional to their
density. We estimated the quasiparticle density due to CSL by adding the CSL
generation rate of quasiparticles to the kinetic equation for the
quasiparticle occupation function. We solved the kinetic equation numerically
to find its steady state solution. With this calculation we obtained a lower
quasiparticle density than the experimental one, so we conclude that other
environmental noise sources are currently giving the dominant contribution to
the experimental excess of quasiparticles.

However, assuming one can eliminate environmental noises, the CSL excess
quasiparticle density still limits the coherence time of a transmon quantum
computer. CSL dissipation does not influence significantly a single qubit, as
the coherence time allowed by CSL dissipation is of the order of $10^{6}$s,
but it is relevant for a NISQ quantum computer composed of many qubits in
which a complex quantum algorithm is run. Indeed the total decoherence time
$T^{\text{\tiny}{CSL}}_{tot}$ of a quantum computer is inversely proportional
to the number of qubits $N$ stacked together in its processor. This implies
that, as the technology scales up to a larger $N$, less gate operations can be
applied before the state of the quantum computer is corrupted by noise, as
shown in figure 3. We showed that important algorithms such as prime number
factorization and molecular simulation could be spoiled by CSL. Our analysis
is performed by assuming that there is no quantum error correction scheme
implemented in the devices. More accurate results could be obtained by
focusing on a specific algorithm to find the resources needed to solve it,
possibly including quantum error correction.

We further explored the possibility of testing CSL models with superconducting
devices. The fact that the experimental values for the density of
quasiparticle and for the subgap quasiparticle current are bigger than the CSL
ones, implies that the detection of CSL effects is currently beyond the
experimental sensitivity of superconducting devices, for which other
environmental sources are dominant. We do not exclude that testing CSL models
may be possible in the future as the technology develops Cardani  _et al._
(2021), given the importance that superconducting devices have for quantum
computing. Our result show that when a superconducting sample is sufficiently
shielded against environmental noises, CSL quasiparticles could be detected at
the current refrigerators temperature.

###### Acknowledgements.

Discussions with M. Paternostro and A. Varlamov are gratefully ackowledged. AB
and LF acknowledge financial support from the H2020 FET Project TEQ (Grant No.
766900). AT and AB acknowledge support form the CNR/RS (London) project
”Testing fundamental theories with ultracold atoms”. AB acknowledges the
Foundational Questions Institute and Fetzer Franklin Fund, a donor advised
fund of Silicon Valley Community Foundation (Grant No. FQXi-RFP-CPW- 2002),
INFN and the University of Trieste.

##  Appendix A BCS theory

The BCS Hamiltonian of the system of electrons is given by:

$\hat{H}_{BCS}=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}}\hat{c}_{\mathbf{k}\sigma}^{\dagger}\hat{c}_{\mathbf{k}\sigma}+\sum_{\mathbf{k}k^{\prime}}U_{\mathbf{k}k^{\prime}}\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}\hat{c}_{-\mathbf{k}^{\prime}\downarrow}\hat{c}_{\mathbf{k}^{\prime}\uparrow}$
(41)

where $\xi_{\mathbf{k}}=\hbar k^{2}/2m-\epsilon_{F}$ is the energy measured
with respect to the Fermi energy $\epsilon_{F}$ and $U_{\mathbf{k}k^{\prime}}$
is the interaction potential. The first term of the Hamiltonian is the kinetic
energy while the second potential term couples pairs of different momenta $k$
and $k^{\prime}$. The ground state of the Hamiltonian is the BCS ground state
(2) in the main text. The form of the BCS Hamiltonian is involved since the
potential term contains four fermionic operators. A simplified form is found
through the mean field procedure. One defines

$a_{\mathbf{k}}=\langle\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}\rangle$
(42)

and assumes that the fluctuations
$(\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}-a_{\mathbf{k}})$
are negligible. Then the following substitution:

$\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}=a_{\mathbf{k}}+\big{(}\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}-a_{\mathbf{k}}\big{)}$
(43)

(and its conjugate) is performed in the BCS Hamiltonian. By keeping terms up
to first order in the fluctuations, one obtains the Bogoliubov Hamiltonian:

$\begin{split}\hat{H}_{B}&=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}}\hat{c}_{\mathbf{k}\sigma}^{\dagger}\hat{c}_{\mathbf{k}\sigma}\\\
&+\sum_{\mathbf{k}k^{\prime}}U_{\mathbf{k}k^{\prime}}[a_{\mathbf{k}^{\prime}}\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}+a_{\mathbf{k}}\hat{c}_{-\mathbf{k}^{\prime}\downarrow}\hat{c}_{\mathbf{k}^{\prime}\uparrow}-a_{\mathbf{k}}a_{\mathbf{k}^{\prime}}]\,,\end{split}$
(44)

and defining
$\Delta_{\mathbf{k}}=\sum_{\mathbf{k}^{\prime}}U_{\mathbf{k}k^{\prime}}a_{\mathbf{k}^{\prime}}$
the Hamiltionian becomes:

$\hat{H}_{B}=\sum_{\mathbf{k}\sigma}\xi_{\mathbf{k}}\hat{c}_{\mathbf{k}\sigma}^{\dagger}\hat{c}_{\mathbf{k}\sigma}-\sum_{\mathbf{k}}\Delta_{\mathbf{k}}[\hat{c}_{\mathbf{k}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}\downarrow}^{\dagger}+\hat{c}_{-\mathbf{k}\downarrow}\hat{c}_{\mathbf{k}\uparrow}+a_{\mathbf{k}}]\,.$
(45)

At this point one performs the Bogoliubov transformation:

$\begin{split}&\hat{\gamma}_{\mathbf{k}\uparrow}=u_{\mathbf{k}}\hat{c}_{\mathbf{k}\uparrow}-v_{\mathbf{k}}e^{i\phi}\hat{c}^{\dagger}_{-\mathbf{k}\downarrow}\,\,\,\,\,\,\hat{\gamma}_{-\mathbf{k}\downarrow}=v_{\mathbf{k}}e^{i\phi}\hat{c}^{\dagger}_{\mathbf{k}\uparrow}+u_{\mathbf{k}}\hat{c}_{-\mathbf{k}\downarrow}\\\
&\hat{\gamma}^{\dagger}_{\mathbf{k}\uparrow}=u_{\mathbf{k}}\hat{c}^{\dagger}_{\mathbf{k}\uparrow}-v_{\mathbf{k}}e^{-i\phi}\hat{c}_{-\mathbf{k}\downarrow}\,\,\,\,\,\,\hat{\gamma}^{\dagger}_{-\mathbf{k}\downarrow}=v_{\mathbf{k}}e^{-i\phi}\hat{c}_{\mathbf{k}\uparrow}+u_{\mathbf{k}}\hat{c}^{\dagger}_{-\mathbf{k}\downarrow}\,,\end{split}$
(46)

and the inverse Bogoliubov transformation is given by:

$\begin{split}&\hat{c}_{\mathbf{k}\uparrow}=u_{\mathbf{k}}\hat{\gamma}_{\mathbf{k}\uparrow}+v_{\mathbf{k}}e^{i\phi}\hat{\gamma}^{\dagger}_{-\mathbf{k}\downarrow}\,\,\,\,\,\,\hat{c}_{-\mathbf{k}\downarrow}=u_{\mathbf{k}}\hat{\gamma}_{-\mathbf{k}\downarrow}-v_{\mathbf{k}}e^{i\phi}\hat{\gamma}^{\dagger}_{\mathbf{k}\uparrow}\\\
&\hat{c}^{\dagger}_{\mathbf{k}\uparrow}=u_{\mathbf{k}}\hat{\gamma}^{\dagger}_{\mathbf{k}\uparrow}+v_{\mathbf{k}}e^{-i\phi}\hat{\gamma}_{-\mathbf{k}\downarrow}\,\,\,\,\,\,\hat{c}^{\dagger}_{-\mathbf{k}\downarrow}=u_{\mathbf{k}}\hat{\gamma}^{\dagger}_{-\mathbf{k}\downarrow}-v_{\mathbf{k}}e^{-i\phi}\hat{\gamma}_{\mathbf{k}\uparrow}\,.\end{split}$
(47)

Substituting Eqs.(47) in the Bogoliubov Hamiltonian one finds:

$\begin{split}\hat{H}_{B}&=\sum_{\mathbf{k}}\big{[}\xi_{\mathbf{k}}(u_{\mathbf{k}}^{2}-v_{\mathbf{k}}^{2})+2\Delta_{\mathbf{k}}u_{\mathbf{k}}v_{\mathbf{k}}\big{]}\big{[}\hat{\gamma}^{\dagger}_{\mathbf{k}\uparrow}\hat{\gamma}_{\mathbf{k}\uparrow}+\hat{\gamma}^{\dagger}_{-\mathbf{k}\downarrow}\hat{\gamma}_{-\mathbf{k}\downarrow}\big{]}\\\
&+\sum_{\mathbf{k}}\big{[}2\xi_{\mathbf{k}}u_{\mathbf{k}}v_{\mathbf{k}}-\Delta_{\mathbf{k}}(u_{\mathbf{k}}^{2}-v_{\mathbf{k}}^{2})\big{]}\big{[}\hat{\gamma}^{\dagger}_{\mathbf{k}\uparrow}\hat{\gamma}^{\dagger}_{-\mathbf{k}\downarrow}+\hat{\gamma}_{-\mathbf{k}\downarrow}\hat{\gamma}_{\mathbf{k}\uparrow}\big{]}\\\
&+\sum_{\mathbf{k}}\big{[}2\xi_{\mathbf{k}}v_{\mathbf{k}}^{2}-2\Delta_{\mathbf{k}}u_{\mathbf{k}}v_{\mathbf{k}}+\Delta_{\mathbf{k}}a_{\mathbf{k}}\big{]}\,.\end{split}$
(48)

This expression contains undesired terms of the type
$\hat{\gamma}\hat{\gamma}$ and $\hat{\gamma}^{\dagger}\hat{\gamma}^{\dagger}$
so the coefficients of these terms are set to zero Grosso and Parravicini
(2013); Tinkham (2004):

$2\xi_{\mathbf{k}}u_{\mathbf{k}}v_{\mathbf{k}}-\Delta_{\mathbf{k}}(u_{\mathbf{k}}^{2}-v_{\mathbf{k}}^{2})=0\,.$
(49)

This condition together with the normalization condition
$u_{\mathbf{k}}^{2}+v_{\mathbf{k}}^{2}=1$ gives:

$\displaystyle
u_{\mathbf{k}}^{2}=\frac{1}{2}\bigg{(}1+\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\bigg{)}$
$\displaystyle
v_{\mathbf{k}}^{2}=\frac{1}{2}\bigg{(}1-\frac{\xi_{\mathbf{k}}}{E_{\mathbf{k}}}\bigg{)}$
(50)

The quasiparticle energies are given by:

$\xi_{\mathbf{k}}(u_{\mathbf{k}}^{2}-v_{\mathbf{k}}^{2})+2\Delta_{\mathbf{k}}u_{\mathbf{k}}v_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+\Delta_{\mathbf{k}}^{2}}\,,$
(51)

and the Bogoliubov Hamiltionian can be expressed in the form:

$\hat{H}_{B}=\sum_{\mathbf{k}}E_{\mathbf{k}}\big{[}\hat{\gamma}^{\dagger}_{\mathbf{k}\uparrow}\hat{\gamma}_{\mathbf{k}\uparrow}+\hat{\gamma}^{\dagger}_{-\mathbf{k}\downarrow}\hat{\gamma}_{-\mathbf{k}\downarrow}\big{]}+W_{S}\,,$
(52)

with $E_{\mathbf{k}}=\sqrt{\xi_{\mathbf{k}}^{2}+\Delta_{\mathbf{k}}^{2}}$, and
$W_{S}=\sum_{\mathbf{k}}\big{[}2\xi_{\mathbf{k}}v_{\mathbf{k}}^{2}-2\Delta_{\mathbf{k}}u_{\mathbf{k}}v_{\mathbf{k}}+\Delta_{\mathbf{k}}a_{\mathbf{k}}\big{]}$.
The superconducting gap parameters are given by:

$\begin{split}\Delta_{\mathbf{k}}=&\sum_{\mathbf{k}^{\prime}}U_{\mathbf{k}k^{\prime}}a_{\mathbf{k}^{\prime}}=\sum_{\mathbf{k}^{\prime}}U_{\mathbf{k}k^{\prime}}\langle\hat{c}_{\mathbf{k}^{\prime}\uparrow}^{\dagger}\hat{c}_{-\mathbf{k}^{\prime}\downarrow}^{\dagger}\rangle\\\
&=\sum_{\mathbf{k}^{\prime}}U_{\mathbf{k}k^{\prime}}u_{\mathbf{k}^{\prime}}v_{\mathbf{k}^{\prime}}\langle
1-\hat{\gamma}^{\dagger}_{\mathbf{k}^{\prime}\uparrow}\hat{\gamma}_{\mathbf{k}^{\prime}\uparrow}-\hat{\gamma}^{\dagger}_{-\mathbf{k}^{\prime}\downarrow}\hat{\gamma}_{-\mathbf{k}^{\prime}\downarrow}\rangle\\\
&=\sum_{\mathbf{k}^{\prime}}U_{\mathbf{k}k^{\prime}}u_{\mathbf{k}^{\prime}}v_{\mathbf{k}^{\prime}}(1-f(E_{\mathbf{k}^{\prime}}))\end{split}$
(53)

and by substituting the expression of $u_{\mathbf{k}}$ and $v_{\mathbf{k}}$
one finds the usual self consistent equations:

$\Delta_{\mathbf{k}}=-\frac{1}{2}\sum_{\mathbf{k}^{\prime}}U_{\mathbf{k}k^{\prime}}\frac{\Delta_{\mathbf{k}^{\prime}}}{E_{\mathbf{k}^{\prime}}}(1-f(E_{\mathbf{k}^{\prime}}))$
(54)

The BCS assumption is that $U_{\mathbf{k}k^{\prime}}=-V$, a negative constant,
for $k$ such that $|\xi_{\mathbf{k}}|<\hbar\omega_{D}$, and
$U_{\mathbf{k}k^{\prime}}=0$ otherwise. In this way one has:

$\Delta_{\mathbf{k}}=\begin{cases}\Delta\,\,for\,\,|\xi_{\mathbf{k}}|<\hbar\omega_{D}\\\
0\,\,for\,\,|\xi_{\mathbf{k}}|>\hbar\omega_{D}\end{cases}$ (55)

The equation for the gap becomes then:

$\Delta=\frac{V}{2}\sum_{\mathbf{k}^{\prime}}\frac{\Delta}{E_{\mathbf{k}^{\prime}}}(1-2f(E_{\mathbf{k}^{\prime}}))$
(56)

Simplifying the common $\Delta$ factor we are left with:

$\begin{split}&1=\frac{V}{2}\sum_{\mathbf{k}^{\prime}}\frac{1}{E_{\mathbf{k}^{\prime}}}(1-2f(E_{\mathbf{k}^{\prime}}))\\\
&1=g(\epsilon_{F})V\bigg{(}\int_{\Delta}^{\hbar\omega_{D}}\frac{dE}{\sqrt{E^{2}-\Delta^{2}}}-\int_{\Delta}^{\hbar\omega_{D}}dE\rho(E)\frac{1}{E}2f(E)\bigg{)}\\\
&1\approx
g(\epsilon_{F})V\bigg{(}\ln{\frac{2\hbar\omega_{D}}{\Delta}}-x_{qp}\bigg{)}\end{split}$
(57)

where in the third line we switched to an integration over the quasiparticle
energies $E$. The latter equation can be solved for the gap:

$\Delta=2\hbar\omega_{D}e^{-1/g(\epsilon_{F})V-x_{qp}}=\Delta(0)e^{-x_{qp}}\approx\Delta(0)(1-x_{qp})$
(58)

where $\Delta_{0}=2\hbar\omega_{D}e^{-1/g(\epsilon_{F})}=1.76k_{B}T_{c}$ is
the superconducting gap with no quasiparticles, so at absolute zero. Notice
that the gap depends on the normalized quasiparticle density and therefore on
the occupation function of quasiparticles. However for small occupation
function, and for small enough temperatures, it is reasonable to approximate
$\Delta=\Delta(0)$.

##  Appendix B CSL Model

The CSL model is usually formulated in position space. The collapse of the
wave function is described by a non-linear and stochastic interaction with a
classical noise through the Itô equation:

$\begin{split}d\ket{\psi}=\bigg{[}&-\frac{i}{\hbar}\hat{H}dt+\frac{\sqrt{\lambda}}{m_{0}}\int
d^{3}x\big{(}\hat{M}(\mathbf{x})-\langle\hat{M}(\mathbf{x})\rangle\big{)}dW_{t}(\mathbf{x})\\\
&-\frac{\lambda}{2m_{0}^{2}}\int
d^{3}x\,d^{3}y\big{(}\hat{M}(\mathbf{x})-\langle\hat{M}(\mathbf{x})\rangle\big{)}G(\mathbf{x}-\mathbf{y})\times\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\big{(}\hat{M}(\mathbf{y})-\langle\hat{M}(\mathbf{y})\rangle\big{)}dt\bigg{]}\ket{\psi}\end{split}$
(59)

where $m_{0}$ is the nucleon mass, $\langle\,\cdot\,\rangle$ denotes the
expectation value on the state $\ket{\psi}$ and $\hat{M}(\mathbf{x})$ is the
mass density operator defined by:

$\hat{M}(\mathbf{x})=\sum_{j}m_{j}\hat{a}_{j}^{\dagger}(\mathbf{x})\hat{a}_{j}(\mathbf{x})\,.$
(60)

The operators $\hat{a}_{j}^{\dagger}(\mathbf{x})$ and
$\hat{a}_{j}(\mathbf{x})$ are the creation and annihilation operators at
position $\mathbf{x}$ of a particle of type $j$ with mass $m_{j}$ (in our case
we will have a single $m$ given by the mass of electrons). The
$G(\mathbf{x}-\mathbf{y})$ in (59) are Gaussian functions of the form:

$G(\mathbf{x}-\mathbf{y})=\frac{1}{(4\pi
r_{c}^{2})^{3/2}}e^{-\frac{1}{4r_{c}^{2}}(x-y)^{2}}\,,$ (61)

that characterize the statistical properties of the noise $W_{t}(\mathbf{x})$.
Indeed, by calling $\xi_{t}(\mathbf{x})=dW_{t}(\mathbf{x})/dt$, one has that
$\mathbb{E}[\xi_{t}(\mathbf{x})]=0$, and the two point correlator
$\mathbb{E}[\xi_{t}(\mathbf{x}),\xi_{s}(\mathbf{y})]=G(\mathbf{x}-\mathbf{y})\delta(t-s)$
where $\mathbb{E}[\cdot]$ denotes the stochastic average. These properties
together with the mass proportionality of $\hat{M}$ guarantee respectively
localization in space and the amplification mechanism: the collapse rate of a
body of $N$ constituents gets amplified linearly in $N$. It is generally
difficult to work directly with Eq. (59), mainly because of its non-linearity.
Since we are interested in expectation values we can use the simplified
linear, but still stochastic, dynamic given by Eq.(12) in the main text. The
CSL term in position space and in the Stratonovich form is given by:

$\hat{H}_{\text{\tiny CSL}}=-\frac{\hbar\sqrt{\lambda}}{m_{0}}\int
d^{3}x\,\xi_{t}(\mathbf{x})\hat{M}(\mathbf{x})\,.$ (62)

This term is related to the second term in equation (59), but now it is linear
because it does not contain $\langle\hat{M}(\mathbf{x})\rangle$ anymore.
Moreover the third term of Eq.(59) is not present. These simplifications are
possible because of the equivalence of Eqs.(59) and (62) at the statistical
level: non-linearity effects are washed away when expectation values are
computed. The term (62) is Fourier transformed to obtain Eq.(13) of the main
text: we work in the normalization volume $V$ to avoid any divergences and the
position representation of the field operators is related to the momentum
operators via:

$\displaystyle\hat{a}(\mathbf{x},s)$
$\displaystyle=\frac{1}{\sqrt{V}}\sum_{\mathbf{k}}e^{i\mathbf{k}\mathbf{x}}\hat{c}_{\mathbf{k}s}$
(63) $\displaystyle\hat{a}^{\dagger}(\mathbf{x},s)$
$\displaystyle=\frac{1}{\sqrt{V}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\mathbf{x}}\hat{c}^{\dagger}_{\mathbf{k}s}\,.$
(64)

##  Appendix C Calculation of the total generation rate of quasiparticles per
unit time and unit volume

We can have a first estimate of the CSL effects by computing the total rate
$\Gamma$ of generation of quasiparticles per unit time and unit volume. This
corresponds to performing a sum over momenta $p$ and $q$ on the transition
probability (31):

$\begin{split}&\Gamma=\frac{\lambda
m^{2}}{m_{0}^{2}V}\sum_{q,p}\widetilde{G}^{2}_{\mathbf{q}-\mathbf{p}}\,M^{2}(q,p)\\\
&=\frac{\lambda m^{2}Vr_{c}}{4m_{0}^{2}\pi^{5/2}}\int dqq\int
dpp\,\big{[}e^{-r_{c}^{2}(p-q)^{2}}-e^{-r_{c}^{2}(p+q)^{2}}\big{]}M^{2}(q,p)\,,\end{split}$
(65)

where in the second line we expressed the summations as integrals in spherical
coordinates.

To make the integral adimensional, we perform the following substitutions:
$\hbar p/\sqrt{2m\Delta}=x$ and $\hbar q/\sqrt{2m\Delta}=y$. We recall that
$\Delta_{\mathbf{k}}=\Delta$ for $|\xi_{\mathbf{k}}|<\hbar\omega_{D}$ and zero
otherwise, which gives the following constraints for the modulus of $x$:

$A_{-}<x<A_{+}$ (66)

where $A_{\pm}=\sqrt{\frac{\epsilon_{F}\pm\hbar\omega_{D}}{\Delta}}$. The same
applies to the modulus of $y$. Performing these substitutions we find:

$\begin{split}\Gamma=&\frac{\lambda
m^{2}Vr_{c}}{4m_{0}^{2}\pi^{5/2}}\bigg{(}\frac{2m\Delta}{\hbar^{2}}\bigg{)}^{2}\int_{A_{-}}^{A_{+}}dy\int_{A_{-}}^{A_{+}}dx\,xy\\\
&\big{(}e^{-\frac{2m\Delta r_{c}^{2}}{\hbar^{2}}(x-y)^{2}}-e^{-\frac{2m\Delta
r_{c}^{2}}{\hbar^{2}}(x+y)^{2}}\big{)}\times\\\
&\bigg{(}1-\frac{(x^{2}-\beta)(y^{2}-\beta)}{[((x^{2}-\beta)^{2}+1)((y^{2}-\beta)^{2}+1)]^{\frac{1}{2}}}\\\
&+\frac{1}{[((x^{2}-\beta)^{2}+1)((y^{2}-\beta)^{2}+1)]^{\frac{1}{2}}}\bigg{)}\\\
\end{split}$ (67)

where $\epsilon_{F}/\Delta=\beta$. We can compare the different values of the
adimensional parameters appearing in the above integral. We have that
$\beta\sim 10^{4}$, $\frac{2m\Delta r_{c}^{2}}{\hbar^{2}}\sim 10^{2}$,
$A_{-}^{2}=\beta-\frac{\hbar\omega_{D}}{\Delta}\sim 10^{4}-10$ and
$A_{+}^{2}=\beta+\frac{\hbar\omega_{D}}{\Delta}\sim 10^{4}+10$. Then the
Gaussian functions can be considered as Dirac deltas, and since we are
integrating in two intervals for $x$ and $y$ where $x$ has the same sign of
$y$, the second Gaussian gives no contribution. By exploiting the Dirac delta
representation:

$\delta(t)=\lim_{\epsilon\to
0}\frac{1}{\epsilon\sqrt{2\pi}}e^{-\frac{1}{2\epsilon^{2}}t^{2}}\,,$ (68)

where in our case $\frac{1}{2\epsilon^{2}}=\frac{2m\Delta
r_{c}^{2}}{\hbar^{2}}$, we are then left with:

$\begin{split}&\Gamma=\frac{\lambda
m^{2}V}{4\pi^{2}m_{0}^{2}}\bigg{(}\frac{\sqrt{2m}}{\hbar}\bigg{)}^{3}\Delta^{2}\bigg{(}\sqrt{\epsilon_{F}}\int_{0}^{+\hbar\omega_{D}}d\xi\,\frac{1}{(\xi^{2}+\Delta^{2})}\bigg{)}\\\
&=\frac{\lambda
m^{2}V}{8m_{0}^{2}\pi}\bigg{(}\frac{\sqrt{2m}}{\hbar}\bigg{)}^{3}\sqrt{\epsilon_{F}}\Delta\,,\end{split}$
(69)

where in the first line we already performed the change of variables from $x$
to $\xi$, and in the second line we approximated the integral to $\pi/2$ since
it is equal to $\tan^{-1}(\hbar\omega_{D}/\Delta)$ with
$\hbar\omega_{D}/\Delta>>1$. Assuming that each quasiparticle is generated
with energy $\Delta$, the total power density supplied by the CSL to a
superconducting sample of volume $V$ is Eq.(69) multiplied by $\Delta$.
Plugging in these equations the parameters for aluminum ($\epsilon_{F}=11.6$eV
and $\Delta=3.4\times 10^{-4}$eV), and a value $\lambda=10^{-10}$s-1, the
total CSL generation rate is $\Gamma\approx 3\times 10^{-11}$s${}^{-1}\mu$m-3
and the power density $P_{tot}=1\times 10^{-33}$W$\mu$m-3. Other works
Martinis  _et al._ (2009) estimated the total generation rate per unit volume
and the power density that would account for the experimental quasiparticle
density measured in transmon qubits. The values that we obtained are $\sim
14-16$ orders of magnitude smaller than the values estimated in Martinis  _et
al._ (2009), $\Gamma_{g}^{ext}=2.4\times 10^{3}$s-1 $\mu$m-3 and
$P_{tot}=6\times 10^{-14}$ W$\mu$m-3, thus showing that we can not attribute
the current experimental excess of quasiparticles to the CSL noise, which is
due to other sources. However we can neglect them, and compute the steady
state quasiparticle density due to CSL, as we do in the main text. In order to
do so, the information that the total generation rate gives is incomplete.
Thermal processes and the CSL generation of quasiparticles contribute to the
evolution of the occupation function $f(E)$ of quasiparticles. In this regime
$f(E)$ redistributes over time to a steady state different from a thermal
state, that enters Eq.(9) to give the quasiparticle density due to CSL.

##  Appendix D Calculation of generation rate of quasiparticles per unit time

We start from Eq.(34) of the main text:

$\begin{split}\gamma_{g}^{\text{\tiny
CSL}}&(E_{\mathbf{q}})=\frac{\lambda(4\pi
r_{c})^{3/2}m^{2}}{m_{0}^{2}(2\pi)^{3}}\int
d^{3}pe^{-r_{c}^{2}(p-q)^{2}}M^{2}(p,q)\bar{f}(E_{\mathbf{p}})\\\
&=\frac{\lambda m^{2}(4\pi)^{3/2}r_{c}}{m_{0}^{2}2(2\pi)^{2}}\int
dp\frac{p}{q}\bigg{(}e^{-r_{c}^{2}(p-q)^{2}}-e^{-r_{c}^{2}(p+q)^{2}}\bigg{)}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,M^{2}(p,q)\bar{f}(E_{\mathbf{p}})\,,\end{split}$
(70)

where in the second line we expressed the integral in polar coordinates. We
make the following substitutions:
$q=\sqrt{2m}/\hbar(\sqrt{E^{2}-\Delta^{2}}+\epsilon_{F})^{1/2}$ and
$p=\sqrt{2m}/\hbar(\sqrt{E^{\prime 2}-\Delta^{2}}+\epsilon_{F})^{1/2}$ to
obtain:

$\begin{split}&\gamma_{g}^{\text{\tiny CSL}}(E)=\frac{m^{2}\lambda
r_{c}}{2\sqrt{\pi}m_{0}^{2}}\frac{\sqrt{2m}}{\hbar}\frac{1}{\sqrt{\sqrt{E-\Delta^{2}}+\epsilon_{F}}}\\\
&\int_{\Delta}^{\infty}dE^{\prime}\bigg{(}e^{-\frac{2mr_{c}^{2}}{\hbar^{2}}(\sqrt{\sqrt{E-\Delta^{2}}+\epsilon_{F})}-\sqrt{\sqrt{E^{\prime}-\Delta^{2}}+\epsilon_{F})})^{2}}\\\
&-e^{-\frac{2mr_{c}^{2}}{\hbar^{2}}(\sqrt{\sqrt{E-\Delta^{2}}+\epsilon_{F})}+\sqrt{\sqrt{E^{\prime}-\Delta^{2}}+\epsilon_{F})})^{2}}\bigg{)}\\\
&\times\rho(E^{\prime})M^{2}(E,E^{\prime})\bar{f}(E^{\prime})\,.\end{split}$
(71)

The above expression is simplified by expanding the squares in the exponential
and collecting the common factors:

$\begin{split}\gamma_{g}^{\text{\tiny CSL}}(E)=&\frac{m^{2}\lambda
r_{c}}{2\sqrt{\pi}m_{0}^{2}}\frac{\sqrt{2m}}{\hbar}\frac{1}{\sqrt{\sqrt{E^{2}-\Delta^{2}}+\epsilon_{F}}}e^{-\frac{2mr_{c}^{2}}{\hbar^{2}}\sqrt{E^{2}-\Delta^{2}}}\\\
&e^{-\frac{4mr_{c}^{2}\epsilon_{F}}{\hbar^{2}}}\int_{\Delta}^{\infty}dE^{\prime}e^{-\frac{2mr_{c}^{2}}{\hbar^{2}}\sqrt{E^{2}_{1}-\Delta^{2}}}\\\
&\bigg{(}e^{\frac{2mr_{c}^{2}}{\hbar^{2}}2(\sqrt{(\sqrt{E^{2}-\Delta^{2}}+\epsilon_{F})(\sqrt{E^{\prime
2}-\Delta^{2}}+\epsilon_{F})}}\\\
&-e^{-\frac{2mr_{c}^{2}}{\hbar^{2}}2(\sqrt{(\sqrt{E^{2}-\Delta^{2}}+\epsilon_{F})(\sqrt{E^{\prime
2}-\Delta^{2}}+\epsilon_{F})}}\bigg{)}\\\
&\times\rho(E^{\prime})M^{2}(E,E^{\prime})\bar{f}(E^{\prime})\,.\end{split}$
(72)

By defining $T_{\text{\tiny CSL}}$ through $k_{B}T_{\text{\tiny
CSL}}=\hbar^{2}/(2mr_{c}^{2})$, and making the substitution $x=E/\Delta$ and
$y=E^{\prime}/\Delta$, we get:

$\begin{split}&\gamma_{g}^{\text{\tiny CSL}}(x)=\frac{m^{2}\lambda
r_{c}}{2\sqrt{\pi}m_{0}^{2}}\frac{\sqrt{2m\Delta}}{\hbar}\frac{1}{\sqrt{s(x)+\beta}}e^{-\frac{\Delta}{k_{B}T_{\text{\tiny
CSL}}}s(x)}\\\ &e^{-\frac{2T_{F}}{T_{\text{\tiny
CSL}}}}\int_{1}^{\infty}dy\,e^{-\frac{\Delta}{k_{B}T_{\text{\tiny
CSL}}}s(y)}\rho(y)M^{2}(x,y)\bar{f}(\Delta y))\\\
&\bigg{(}e^{\frac{2\Delta}{k_{B}T_{\text{\tiny
CSL}}}(\sqrt{(s(x)+\beta)(s(y)+\beta)}}-e^{-\frac{2\Delta}{k_{B}T_{\text{\tiny
CSL}}}(\sqrt{(s(x)+\beta)(s(y)+\beta})}\bigg{)}\,,\end{split}$ (73)

where we called $s(x)=\sqrt{x^{2}-1}$. We can further simplify this expression
by neglecting the negative term in the last line because it is exponentially
suppressed. We can finally write:

$\begin{split}&\gamma_{g}^{\text{\tiny CSL}}(x)=\frac{m^{2}\lambda
r_{c}}{2\sqrt{\pi}m_{0}^{2}}\frac{\sqrt{2m\Delta}}{\hbar}\frac{1}{\sqrt{s(x)+\beta}}e^{-\frac{\Delta}{k_{B}T_{\text{\tiny
CSL}}}s(x)}\\\ &e^{-\frac{2T_{F}}{T_{\text{\tiny
CSL}}}}\int_{1}^{\infty}dy\,e^{-\frac{\Delta}{k_{B}T_{\text{\tiny
CSL}}}s(y)}\,e^{\frac{2\Delta}{k_{B}T_{\text{\tiny
CSL}}}(\sqrt{(s(x)+\beta)(s(y)+\beta)}}\\\
&\times\rho(y)\frac{1}{2}\bigg{(}1-\frac{\sqrt{x^{2}-1}\sqrt{y^{2}-1}}{xy}+\frac{1}{xy}\bigg{)}\bar{f}(\Delta
y)\,.\end{split}$ (74)

This is Eq.(35) of the main text.

##  Appendix E Steady state solution

In the main text, we have seen that the steady state solution to Eq.(32) when
$\gamma_{g}^{\text{\tiny ext}}=0$, is the Fermi Dirac distribution function.
This observation, together with the fact that the CSL injection rate
$\gamma_{g}^{\text{\tiny CSL}}(E)$ is small, suggests the following procedure
to approximate the steady state solution when $\gamma_{g}^{\text{\tiny
ext}}=\gamma_{g}^{\text{\tiny CSL}}$: since the steady state solution without
injection rate is a Fermi-Dirac distribution $f_{FD}(E,T)$, at some
temperature $T$, when we add the CSL injection term, the modified steady state
solution can be written as $f_{SS}(E)=f_{FD}(E,T)+\delta f(E)$ with $\delta
f(E)$ a small perturbation.

We start from Eq.(32) and we neglect the rate terms associated to
recombination and generation, since they are exponentially small. We plug
$f_{SS}=f_{FD}+\delta f$ into Eq.(32), and we perform the substitutions
$x=E/\Delta$, $y=E^{\prime}/\Delta$ to get:

$\begin{split}&\gamma_{g}^{\text{\tiny
CSL}}(x)+\gamma_{0}\int_{x}^{4}dyS(x,y)\big{[}(1-f_{FD}(x)-\delta f(x))\\\
&(f_{FD}(y)+\delta f(y))(N(y-x)+1)-(f_{FD}(x)+\delta f(x))\\\
&(1-f_{FD}(y))-(f_{FD}(x)+\delta f(x))(-\delta f(y))(N(y-x))\big{]}\\\
&+\gamma_{0}\int_{1}^{x}dyS(x,y)\big{[}(1-f_{FD}(x)-\delta
f(x))(f_{FD}(y)+\delta f(y))\\\ &(N(x-y))-(f_{FD}(x)+\delta
f(x))(1-f_{FD}(y)\\\ &-\delta f(y))(N(x-y)+1)\big{]}=0\,.\end{split}$ (75)

The integration limit is $x=4$ since we are interested in the low energy
behaviour of the occupation function.

Developing the product of all these factors, the terms coming from the Fermi-
Dirac distribution cancel out. Then, keeping terms up to first order in
$\delta f$ one has:

$\begin{split}&\gamma_{g}^{\text{\tiny
CSL}}(x)+\gamma_{0}\int_{x}^{4}dyS(x,y)\big{[}\bar{f}_{FD}(x)\delta
f(y)(N(y-x)+1)\\\ &-\delta f(x)f_{FD}(N(y-x)+1)-\delta
f(x)\bar{f}_{FD}(y)N(y-x)\\\ &+\delta
f(x)f_{FD}(y)N(y-x)\big{]}+\gamma_{0}\int_{1}^{x}dyS(x,y)\big{[}\bar{f}_{FD}(x)\\\
&\delta f(y)N(x-y)-f_{FD}(y)\delta f(x)(N(x-y))-\delta f(x)\bar{f}_{FD}(y)\\\
&(N(x-y)+1)+f_{FD}(x)\delta f(y)(N(x-y)+1)\big{]}=0\,.\end{split}$ (76)

where $\bar{f}_{FD}(x)=1-f_{FD}(x)$. This is further simplified by neglecting
terms proportional to $N(\Omega)$ (see Eq.(33) of the main text), since they
are exponentially small for the temperatures that we consider. In this way one
ends up with:

$\begin{split}&\gamma_{g}^{\text{\tiny
CSL}}(x)+\gamma_{0}\int_{x}^{4}dyS(x,y)[(\bar{f}_{FD}(x))\delta
f(y)(N(y-x)+1)\\\ &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-\delta
f(x)f_{FD}(y)(N(y-x)+1)]\\\ &+\gamma_{0}\int_{1}^{x}dyS(x,y)[(f_{FD}(x))\delta
f(y)(N(x-y)+1)\\\ &\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,-\delta
f(x)\bar{f}_{FD}(y)(N(x-y)+1)]=0\,\end{split}$ (77)

where $\bar{f}_{FD}=(1-f_{FD})$.

From now on, we will approximate $\bar{f}_{FD}(x)=1-f_{FD}(x)\approx 1$ and
$N(\Omega)+1\approx 1$, because at milliKelvin temperatures and in the
considered energy interval $[\Delta,4\Delta]$, $f_{FD}(E)\sim
e^{-E/k_{b}T}\approx 0$ and $N(\Omega)\sim e^{-\Omega/k_{b}T}\approx 0$. The
above equation thus simplifies to:

$\begin{split}&\gamma_{g}^{\text{\tiny
CSL}}(x)+\gamma_{0}\int_{x}^{4}dyS(x,y)[\delta f(y)-\delta f(x)f_{FD}(y)]\\\
&+\gamma_{0}\int_{1}^{x}dyS(x,y)[(f_{FD}(x))\delta f(y)-\delta
f(x)]=0\,.\end{split}$ (78)

To further simplify this expression, we note that, for milliKelvin
temperatures, we can neglect the addends of the form $\delta f(x)f_{FD}(y)$
and $f_{FD}(x)\delta f(y)$, since we computed the corresponding integrals and
they turned out to be small quantities with respect to the other terms in
Eq.(78). Finally one has that:

$\begin{split}&\gamma_{g}^{\text{\tiny
CSL}}(x)+\gamma_{0}\int_{x}^{4}dyS(x,y)\delta f(y)\\\
&-\gamma_{0}\int_{1}^{x}dyS(x,y)\delta f(x)=0\,,\end{split}$ (79)

which can be inverted to find the following equation for $\delta f(x)$:

$\begin{split}\delta f(x)&=\frac{\gamma_{g}^{\text{\tiny
CSL}}(x)}{\gamma_{0}\int_{1}^{x}dyS(x,y)}+\frac{\int_{x}^{4}dyS(x,y)\delta
f(y)}{\gamma_{0}\int_{1}^{x}dyS(x,y)}\\\ &\approx\frac{\gamma_{g}^{\text{\tiny
CSL}}(x)}{\gamma_{0}\int_{1}^{x}dyS(x,y)}\,.\end{split}$ (80)

This shows that, by knowing the CSL generation rate, we can compute the
correction to the Fermi-Dirac distribution, and thus the steady state solution
of Eq.(32) without the need to solve it. The approximation in the second line
of Eq.(80) is justified because if we substitute $\delta
f(y)=\frac{\gamma_{g}^{\text{\tiny CSL}}(y)}{\gamma_{0}\int_{1}^{y}dzS(y,z)}$,
in the second term of Eq.(80) we obtain a small correction to the first term.
In the next section we check the validity of this approximation scheme.

##  Appendix F Validity of the approximation scheme

To test the validity of this approximation scheme, we have chosen three
representative starting equilibrium temperatures ($65$ mK, $45$ mK and $25$
mK) to check two aspects. First, we compared $f_{SS}=f_{FD}+\delta f$ with
$f_{FD}$, to check if $\delta f$ is indeed a small perturbation to $f_{FD}$.
Second, we compared $f_{SS}$ with the numerical solution of (32) (with
$\gamma_{g}^{\text{\tiny ext}}(E)=\gamma_{g}^{\text{\tiny CSL}}(E)$), to check
if $f_{SS}$ is a good approximation to the steady state solution of Eq.(32).
We summarize the results in figure 4, where we plot the occupation function
obtained with the numerical simulation (blue solid line), the analytical
expression $f_{SS}=f_{FD}+\delta f$ (red dashed line) and the initial Fermi
Dirac distribution $f_{FD}$ (black dotted line), for a starting equilibrium
temperature of $65$ mK (a), $45$ mK (b) and $25$ mK (c).

In figure 4a the perturbation $\delta f$ is small with respect to $f_{FD}$
only for energies lower than $\sim 1.5\Delta$, as the red dashed line
representing $f_{SS}=f_{FD}+\delta f$ is close to the black dotted line
representing $f_{FD}$ only in this energy interval. Nevertheless the
analytical expression $f_{SS}$ (red dashed line) is close to the numerical
result represented by the blue solid line in the whole energy interval. This
means that, despite the fact that we are not allowed to treat $\delta f$ as a
small perturbation to $f_{FD}$, $f_{SS}$ is a good approximation to the steady
state solution of Eq.(32).

In figure 4b the analytical solution (red dashed line) is many orders of
magnitude larger than the initial Fermi-Dirac occupation function (black
dotted line) in the whole energy interval. Thus, in this case $\delta f$ is
not a small perturbation to $f_{FD}$. The same applies to figure 4c. Note
however that, analogously to figure 4a, in both figures 4b and 4c the
numerical result (blue solid line) is well approximated by $f_{SS}$ (red
dashed line).

In conclusion, for the temperatures considered, even when $\delta f$ in
Eq.(80) is not a small perturbation to $f_{FD}$, we can use
$f_{SS}=f_{FD}+\delta f$ as the steady state solution of Eq.(32).

Note moreover that, the fact that for $T=45$ mK and $T=25$ mK, the correction
$\delta f\gg f_{FD}$ implies that $f_{SS}=f_{FD}+\delta f\approx\delta f$.
That is to say, $\delta f$ given by Eq.(80), is itself the good approximation
to the steady state solution of Eq.(32). As $\delta f$ turns out to be almost
independent of temperature for milliKelvin temperatures, this suggests that
the steady state solution for any $T<45$ mK, can always be well approximated
by $\delta f$. This fact was confirmed by additional simulations (here not
shown) with starting equilibrium temperatures at $T=10$ mK and $T=5$ mK, for
which the numerical result was still well approximated by $f_{SS}\approx\delta
f$.

Figure 4: Log plot of the quasiparticle occupation function obtained solving
Eq.(32) numerically (blue solid line), of the the analytical expression
$f_{SS}=f_{FD}+\delta f$ with $\delta f$ given by Eq.(80) (red dashed line)
and of the starting Fermi-Dirac distribution $f_{FD}$ (black dotted line) for
a starting equilibrium temperature of $65$ mK (a), $45$ mK (b) and $25$ mK
(c). In a) $\delta f$ is a small perturbation to $f_{FD}$ only for energies
lower than $\sim 1.5\Delta$, where the red dashed line representing $f_{SS}$
is close to the black dotted line representing $f_{FD}$. Despite this, the red
dashed line is very close to the blue solid line, meaning that $f_{SS}$ is a
good approximation to the occupation function obtained by numerical
simulation. In b) the correction $\delta f\gg f_{FD}$, since the red dashed
line representing $f_{SS}$ is many orders of magnitude larger than the black
dotted line representing $f_{FD}$. The same applies for $T=25$mK in c).
Analogously to a), for both b) and c) the red dashed line representing the
analytical solution $f_{SS}$ is a good approximation to the blue solid line
representing the occupation function obtained by numerical simulation. Note
that in c) we multiplied by $10^{16}$ the initial $f_{FD}$ to have the same
interval in the vertical axis of b).

## References

  * Shor (1999) Peter W Shor, “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,” SIAM review 41, 303–332 (1999).
  * Brassard  _et al._ (2002) Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp, “Quantum amplitude amplification and estimation,” Contemporary Mathematics 305, 53–74 (2002).
  * Montanaro (2015) Ashley Montanaro, “Quantum speedup of monte carlo methods,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471, 20150301 (2015).
  * Cao  _et al._ (2018) Yudong Cao, Jhonathan Romero, and Alán Aspuru-Guzik, “Potential of quantum computing for drug discovery,” IBM Journal of Research and Development 62, 6–1 (2018).
  * Childs  _et al._ (2018) Andrew M Childs, Dmitri Maslov, Yunseong Nam, Neil J Ross, and Yuan Su, “Toward the first quantum simulation with quantum speedup,” Proceedings of the National Academy of Sciences 115, 9456–9461 (2018).
  * Kjaergaard  _et al._ (2020) Morten Kjaergaard, Mollie E Schwartz, Jochen Braumüller, Philip Krantz, Joel I-J Wang, Simon Gustavsson, and William D Oliver, “Superconducting qubits: Current state of play,” Annual Review of Condensed Matter Physics 11, 369–395 (2020).
  * Bruzewicz  _et al._ (2019) Colin D Bruzewicz, John Chiaverini, Robert McConnell, and Jeremy M Sage, “Trapped-ion quantum computing: Progress and challenges,” Applied Physics Reviews 6, 021314 (2019).
  * Bombin  _et al._ (2021) Hector Bombin, Isaac H Kim, Daniel Litinski, Naomi Nickerson, Mihir Pant, Fernando Pastawski, Sam Roberts, and Terry Rudolph, “Interleaving: Modular architectures for fault-tolerant photonic quantum computing,” arXiv preprint arXiv:2103.08612 (2021).
  * Recher and Trauzettel (2010) Patrik Recher and Björn Trauzettel, “Quantum dots and spin qubits in graphene,” Nanotechnology 21, 302001 (2010).
  * Koch  _et al._ (2007) Jens Koch, T.M. Yu, Jay Gambetta, A.A. Houck, D.I. Schuster, Johannes Majer, Alexandre Blais, M.H. Devoret, S.M. Girvin, and R.J. Schoelkopf, “Charge-insensitive qubit design derived from the cooper pair box,” Physical Review A 76, 042319 (2007).
  * Arute  _et al._ (2019) Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Rupak Biswas, Sergio Boixo, Fernando GSL Brandao, David A Buell, _et al._ , “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * Gambetta (2020) Jay Gambetta, “Ibm’s roadmap for scaling quantum technology,” IBM Research Blog https://www. ibm. com/blogs/research/2020/09/ibm-quantum-roadmap (2020).
  * Sete  _et al._ (2016) Eyob A Sete, William J Zeng, and Chad T Rigetti, “A functional architecture for scalable quantum computing,” in  _2016 IEEE International Conference on Rebooting Computing (ICRC)_ (IEEE, 2016) pp. 1–6.
  * Wu  _et al._ (2021) Yulin Wu, Wan-Su Bao, Sirui Cao, Fusheng Chen, Ming-Cheng Chen, Xiawei Chen, Tung-Hsun Chung, Hui Deng, Yajie Du, Daojin Fan, _et al._ , “Strong quantum computational advantage using a superconducting quantum processor,” arXiv preprint arXiv:2106.14734 (2021).
  * Córcoles  _et al._ (2011) Antonio D Córcoles, Jerry M Chow, Jay M Gambetta, Chad Rigetti, James R Rozen, George A Keefe, Mary Beth Rothwell, Mark B Ketchen, and Matthias Steffen, “Protecting superconducting qubits from radiation,” Applied Physics Letters 99, 181906 (2011).
  * Wang  _et al._ (2015) Chen Wang, Christopher Axline, Yvonne Y Gao, Teresa Brecht, Yiwen Chu, Luigi Frunzio, MH Devoret, and Robert J Schoelkopf, “Surface participation and dielectric loss in superconducting qubits,” Applied Physics Letters 107, 162601 (2015).
  * Klimov  _et al._ (2018) P.V. Klimov, Julian Kelly, Z Chen, Matthew Neeley, Anthony Megrant, Brian Burkett, Rami Barends, Kunal Arya, Ben Chiaro, Yu Chen, A Dunsworth, A Fowler, B Foxen, C Gidney, M Giustina, R Graff, T Huang, E Jeffrey, E Lucero, J.Y. Mutus, O Naaman, C Neill, C Quintana, P Roushan, D Sank, A Vainsencher, J Wenner, T.C. White, S Boixo, R Babbush, V.N. Smelyanskiy, H Neven, and J.M. Martinis, “Fluctuations of energy-relaxation times in superconducting qubits,” Physical review letters 121, 090502 (2018).
  * Catelani  _et al._ (2011a) Gianluigi Catelani, Robert J Schoelkopf, Michel H Devoret, and Leonid I Glazman, “Relaxation and frequency shifts induced by quasiparticles in superconducting qubits,” Physical Review B 84, 064517 (2011a).
  * Glazman and Catelani (2021) Leonid Glazman and Gianluigi Catelani, “Bogoliubov quasiparticles in superconducting qubits,” SciPost Physics Lecture Notes , 031 (2021).
  * Martinis  _et al._ (2009) John M Martinis, M Ansmann, and J Aumentado, “Energy decay in superconducting josephson-junction qubits from nonequilibrium quasiparticle excitations,” Physical review letters 103, 097002 (2009).
  * Lloyd (2000) Seth Lloyd, “Ultimate physical limits to computation,” Nature 406, 1047–1054 (2000).
  * Gambini  _et al._ (2005) Rodolfo Gambini, Rafael A Porto, and Jorge Pullin, “Fundamental gravitational limitations to quantum computing,” arXiv preprint quant-ph/0507262 (2005).
  * Bassi and Ghirardi (2003) Angelo Bassi and GianCarlo Ghirardi, “Dynamical reduction models,” Physics Reports 379, 257–426 (2003).
  * Ferialdi and Bassi (2020) Luca Ferialdi and Angelo Bassi, “Continuous spontaneous localization reduction rate for rigid bodies,” Physical Review A 102, 042213 (2020).
  * Schrieffer (2018) John Robert Schrieffer,  _Theory of superconductivity_ (CRC press, 2018).
  * Grosso and Parravicini (2013) Giuseppe Grosso and Giuseppe Pastori Parravicini,  _Solid state physics_ (Academic press, 2013).
  * Catelani  _et al._ (2011b) Gianluigi Catelani, Jens Koch, Luigi Frunzio, RJ Schoelkopf, Michel H Devoret, and LI Glazman, “Quasiparticle relaxation of superconducting qubits in the presence of flux,” Physical review letters 106, 077002 (2011b).
  * de Visser  _et al._ (2014) P.J. de Visser, D.J. Goldie, P Diener, S Withington, J.J.A. Baselmans, and T.M. Klapwijk, “Evidence of a nonequilibrium distribution of quasiparticles in the microwave response of a superconducting aluminum resonator,” Physical review letters 112, 047004 (2014).
  * Serniak  _et al._ (2019) K Serniak, S Diamond, M Hays, V Fatemi, S Shankar, L Frunzio, R.J. Schoelkopf, and M.H. Devoret, “Direct dispersive monitoring of charge parity in offset-charge-sensitive transmons,” Physical Review Applied 12, 014052 (2019).
  * Carlesso and Donadi (2019) Matteo Carlesso and Sandro Donadi, “Collapse models: main properties and the state of art of the experimental tests,” in  _Advances in Open Systems and Fundamental Tests of Quantum Mechanics_ (Springer, 2019) pp. 1–13.
  * Adler (2007) Stephen L Adler, “Lower and upper bounds on csl parameters from latent image formation and igm heating,” Journal of Physics A: Mathematical and Theoretical 40, 2935 (2007).
  * Wang  _et al._ (2014) Chen Wang, Yvonne Y Gao, Ioan M Pop, Uri Vool, Chris Axline, Teresa Brecht, Reinier W Heeres, Luigi Frunzio, Michel H Devoret, Gianluigi Catelani, _et al._ , “Measurement and control of quasiparticle dynamics in a superconducting qubit,” Nature communications 5, 1–7 (2014).
  * Goldie and Withington (2012) DJ Goldie and S Withington, “Non-equilibrium superconductivity in quantum-sensing superconducting resonators,” Superconductor Science and Technology 26, 015004 (2012).
  * Chang and Scalapino (1977) Jhy-Jiun Chang and DJ Scalapino, “Kinetic-equation approach to nonequilibrium superconductivity,” Physical Review B 15, 2651 (1977).
  * Milliken  _et al._ (2004) FP Milliken, RH Koch, JR Kirtley, and JR Rozen, “The subgap current in nb/ alox/ nb tunnel junctions,” Applied physics letters 85, 5941–5943 (2004).
  * Barone and Paterno (1982) Antonio Barone and Gianfranco Paterno,  _Physics and applications of the Josephson effect_ , Vol. 1 (Wiley Online Library, 1982).
  * de Visser  _et al._ (2011) P.J. de Visser, J.J.A. Baselmans, P. Diener, S.J.C. Yates, A Endo, and T.M. Klapwijk, “Number fluctuations of sparse quasiparticles in a superconductor,” Physical review letters 106, 167004 (2011).
  * Vepsäläinen  _et al._ (2020) Antti P Vepsäläinen, Amir H Karamlou, John L Orrell, Akshunna S Dogra, Ben Loer, Francisca Vasconcelos, David K Kim, Alexander J Melville, Bethany M Niedzielski, Jonilyn L Yoder, _et al._ , “Impact of ionizing radiation on superconducting qubit coherence,” Nature 584, 551–556 (2020).
  * Gottesman (2010) Daniel Gottesman, “An introduction to quantum error correction and fault-tolerant quantum computation,” in  _Quantum information science and its contributions to mathematics, Proceedings of Symposia in Applied Mathematics_ , Vol. 68 (2010) pp. 13–58.
  * Preskill (2018) John Preskill, “Quantum computing in the nisq era and beyond,” Quantum 2, 79 (2018).
  * Mooney  _et al._ (2019) Gary J Mooney, Charles D Hill, and Lloyd CL Hollenberg, “Entanglement in a 20-qubit superconducting quantum computer,” Scientific reports 9, 1–8 (2019).
  * Palma  _et al._ (1996) G Massimo Palma, Kalle-Antti Suominen, and Artur Ekert, “Quantum computers and dissipation,” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 452, 567–584 (1996).
  * Ischi  _et al._ (2005) Boris Ischi, Michael Hilke, and Martin Dubé, “Decoherence in a n-qubit solid-state quantum register,” Physical Review B 71, 195325 (2005).
  * Salm  _et al._ (2020) Marie Salm, Johanna Barzen, Frank Leymann, and Benjamin Weder, “About a criterion of successfully executing a circuit in the nisq era: what $wd\ll 1/\epsilon_{eff}$ really means,” in  _Proceedings of the 1st ACM SIGSOFT International Workshop on Architectures and Paradigms for Engineering Quantum Software_ (2020) pp. 10–13.
  * Moll  _et al._ (2018) Nikolaj Moll, Panagiotis Barkoutsos, Lev S Bishop, Jerry M Chow, Andrew Cross, Daniel J Egger, Stefan Filipp, Andreas Fuhrer, Jay M Gambetta, Marc Ganzhorn, Abhinav Kandala, Antonio Mezzacapo, Peter Müller, Walter Riess, Gian Salis, John Smolin, Ivano Tavernelli, and Kristan Temme, “Quantum optimization using variational algorithms on near-term quantum devices,” Quantum Science and Technology 3, 030503 (2018).
  * Reiher  _et al._ (2017) Markus Reiher, Nathan Wiebe, Krysta M Svore, Dave Wecker, and Matthias Troyer, “Elucidating reaction mechanisms on quantum computers,” Proceedings of the National Academy of Sciences 114, 7555–7560 (2017).
  * Kutin (2006) Samuel A Kutin, “Shor’s algorithm on a nearest-neighbor machine,” arXiv preprint quant-ph/0609001 (2006).
  * Chakrabarti  _et al._ (2021) Shouvanik Chakrabarti, Rajiv Krishnakumar, Guglielmo Mazzola, Nikitas Stamatopoulos, Stefan Woerner, and William J Zeng, “A threshold for quantum advantage in derivative pricing,” Quantum 5, 463 (2021).
  * Cardani  _et al._ (2021) Laura Cardani, Francesco Valenti, Nicola Casali, Gianluigi Catelani, Thibault Charpentier, Massimiliano Clemenza, Ivan Colantoni, Angelo Cruciani, G D’Imperio, L Gironi, _et al._ , “Reducing the impact of radioactivity on quantum circuits in a deep-underground facility,” Nature communications 12, 1–6 (2021).
  * Tinkham (2004) Michael Tinkham,  _Introduction to superconductivity_ (Courier Corporation, 2004).