# Quenched dynamics and pattern formation in clean and disordered Bogoliubov-
de Gennes superconductors

Bo Fan<EMAIL_ADDRESS>Shanghai Center for Complex Physics, School of
Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China  Antonio M. García-
García<EMAIL_ADDRESS>Shanghai Center for Complex Physics, School of
Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China

(August 28, 2024)

###### Abstract

We study the quench dynamics of a two dimensional superconductor in a lattice
of size up to $200\times 200$ employing the self-consistent time dependent
Bogoliubov-de Gennes (BdG) formalism. In the clean limit, the dynamics of the
order parameter for short times, characterized by a fast exponential growth
and an oscillatory pattern, agrees with the Bardeen-Cooper-Schrieffer (BCS)
prediction. However, unlike BCS, we observe for longer times an universal
exponential decay of these time oscillations that we show explicitly to be
induced by the full emergence of spatial inhomogeneities of the order
parameter, even in the clean limit, characterized by the exponential growth of
its variance. The addition of a weak disorder does not alter these results
qualitatively. In this region, the spatial inhomogeneities rapidly develops
into an intricate spatial structure consisting of ordered fragmented stripes
in perpendicular directions where the order parameter is heavily suppressed
especially in the central region. As the disorder strength increases, the
fragmented stripes gradually turn into a square lattice of approximately
circular spatial regions where the condensate is heavily suppressed. A further
increase of disorder leads to the deformation and ultimate destruction of this
lattice. We explore suitable settings for the experimental confirmation of
these findings.

##  I Introduction

The spontaneous formation of patterns, defects and other spatial structures is
a fascinating phenomenon rather ubiquitous in nature. It can be observed in
various contexts, from the vortex pattern in superconducting thin films
induced by a magnetic field [1] to the formation of spatial structures in a
driven Bose-Hubbard model [2]. A particularly interesting phenomenon of defect
formation, termed Kibble-Zurek mechanism [3, 4], is the spontaneous generation
of vortices as a result of the quench through a second-order phase transition.
It has been observed both experimentally [5, 6, 7] and numerically [8, 9, 10]
in a wide variety of physical systems.

In many of these situations, far from equilibrium dynamics triggers spatial
instabilities that eventually lead to pattern formation. In the context of
superconductivity and superfluidity, although the study of nonequilibrium
dynamics has received a lot of attention, the spontaneous generation of
spatial structures, either due to the Kibble-Zurek mechanism or of different
origin, have been modeled by employing phenomenological approaches such as the
time dependent Ginzburg-Landau [11, 1, 12, 13], the Gross-Pitaevskii equation
[8, 2] or applied holography [14, 9] where the dynamics of the superconductor
is mapped onto that of a gravitational system.

It is expected that these phenomenological approaches will be qualitatively
correct close to a second order phase transition. However, the full non-linear
structure of the time dependent BdG equations [15, 16], due to the self-
consistent condition verified by the order parameter, is fully necessary for
the quantitative description of the out of equilibrium dynamics of a
superconductor. A simpler problem, the quench dynamics of a Bardeen-Cooper-
Schrieffer (BCS) superconductor [17], first investigated in Ref. [18], has
received a lot of recent attention [19, 20, 21, 22, 23, 24, 25]. The
conclusion of these studies is that details of the dynamics of the order
parameter amplitude depend on both the initial state and the quench protocol,
though a generic feature is the existence of oscillations in time. For an
initial state above the critical temperature, it has been argued [19, 26, 21]
that these oscillations do not decay with time, unless collision effects
beyond BCS are taken into account, while quenches in the coupling constant
within the superconducting state [21, 22, 23] lead to oscillations whose
amplitude typically decays either exponentially or as a power-law.

A perturbative analytic treatment of the dynamics using the full time
dependent BdG formalism, that accounts for spatial inhomogeneities, showed
that eventually the order parameter develops a simple oscillating spatial
structure with a typical length directly related to the superconducting
coherence length [27]. This perturbative treatment cannot take into account
the full non-linear nature of the time evolution so the validity of these
results is only assured for relatively short time scales where these non-
linear effects are small. A very similar pattern of spatial oscillations has
been reported in the quench dynamics of the order parameter of a holographic
superconductor [14].

Despite this considerable progress, the dynamics of a superconductor after a
quench, especially the nature of the emergent spatial structures, is still an
open problem. We note that this is mostly due to technical challenges
resulting from the combination of the non-linearity induced by the self-
consistent condition and the requirement of sufficiently large system sizes in
order to account for the emergent spatial structure from the quench dynamics.
Moreover, it is necessary to consider at least a two dimensional
superconductor because fluctuations in one dimension, even at low temperature,
are too large to employ a mean field formalism.

In this paper, we address this problem by studying the quench dynamics of a
two dimensional superconductor by the full self-consistent time-dependent BdG
formalism [28, 29, 30] in a $200\times 200$ square lattice that enables us to
investigate in detail complex spatial patterns. We have found that for no
disorder and short times, in agreement with the BCS results, the order
parameter first grows exponentially and then has an oscillatory behavior.
However, for longer times, time oscillations in the order parameter are
suppressed exponentially independently of the quench protocol. The precursor
of this behavior is the emergence of spatial inhomogeneities, beyond the reach
of the BCS formalism, characterized by an exponential growth of the variance
in space of the order parameter which ultimately results in the appearance of
short stripes in the horizontal and vertical directions where the order
parameter is heavily suppressed. We believe our results are largely universal
as these spatial patterns occur well after the quench ends. The addition of a
weak disordered potential, modeling impurities which are ubiquitous in
experiments, does not change the above results qualitatively. As disorder
increases, the mentioned fragmented stripes gradually morphs into a square
lattice of fake vortices, namely, approximately circular regions where the
amplitude of the order parameter is very small but with a trivial phase.
Finally, as the insulating transition is approached, the lattice symmetry is
eventually lost though the repulsion between fake vortices persists. In the
next section, we introduce the model and the computation scheme.

##  II The model

In order to study the time evolution of a two dimensional superconductor after
a temperature quench, we employ the mean field time dependent BdG equations
[31, 28, 32, 33, 29, 34], which are given by

$\left(\begin{matrix}\hat{K}&\hat{\Delta}({\bf r}_{i},t)\\\
\hat{\Delta}^{*}({\bf r}_{i},t)&-\hat{K}^{*}\\\
\end{matrix}\right)\left(\begin{matrix}u_{n}({\bf r}_{i},t)\\\ v_{n}({\bf
r}_{i},t)\\\ \end{matrix}\right)=i\hbar{\partial\over\partial
t}\left(\begin{matrix}u_{n}({\bf r}_{i},t)\\\ v_{n}({\bf r}_{i},t)\\\
\end{matrix}\right)$ (1)

where

$\hat{K}u_{n}({\bf r}_{i})=-t_{i,i+\delta}\sum_{\delta}u_{n}({\bf
r}_{i+\delta})+(V_{i}-\mu)u_{n}({\bf r}_{i}),$ (2)

$\delta$ stands for the nearest neighbors sites, and $t_{i,i+\delta}$ is the
hopping energy between the nearest neighbors sites and we set
$t_{i,i+\delta}=1$ for simplicity. The onsite random potential $V_{i}$ is
uniformly distributed between $[-V,V]$, where $V$ is the strength of the
disordered potential and $\mu$ is the chemical potential.

In order to simulate the dynamical evolution, we solve the above equations by
using the fourth order Runge Kutta algorithm [35, 29]. The occupation number
is assumed to satisfy the Fermi-Dirac distribution
$f(E_{n})=[1+\exp(E_{n}/k_{B}T)]^{-1}$ at each time-step during the dynamical
process, where $T$ is temperature and $k_{B}=1$ is the Boltzmann constant. The
time dependent order parameter is then defined as,

$\Delta({\bf r}_{i},t)=|U|\sum_{n}u_{n}({\bf r}_{i},t)v_{n}^{*}({\bf
r}_{i},t)[1-2f(E_{n})],$ (3)

where $U$ is the strength of the on-site, phonon induced, attractive electron-
electron interaction, leading to the superconducting state. The time dependent
local density is given by

$n({\bf r}_{i},t)=2\sum_{n}[|u_{n}({\bf r}_{i},t)|^{2}f(E_{n})+|v_{n}({\bf
r}_{i},t)|^{2}(1-f(E_{n}))].$ (4)

For numerical convenience, we start from the equilibrium state at temperature
$T_{i}>T_{c}$, namely, a vanishing order parameter which can be obtained from
the exact solution of the BdG equations [36, 37, 38]. Time evolution is
induced by lowering the temperature $T(t)$ from $T_{i}>T_{c}$ to
$T_{f}=0.1T_{c}$ using the following linear quench protocol [9, 39],

$T(t)=\left\\{\begin{array}[]{lr}T_{i}-\tau_{Q}t,&t_{i}\leq t\leq t_{f}\\\
T_{f},&t>t_{f}\end{array}\right.$ (5)

where $\tau_{Q}$ is the slope of the quench, which characterizes the quench
duration. $t_{i}=0$ is the starting time of the quench, and
$t_{f}=(T_{i}-T_{f})/\tau_{Q}$ is the quench ending time corresponding to the
final temperature $T_{f}$. In our study, we let $T_{i}=1.2T_{c}$. Our main
focuses are fast quenches leading to a non-adiabatic time evolution, so we set
$\tau_{Q}=50$. We did not find qualitative differences using smaller values of
$\tau_{Q}$. We are mostly interested in the generation of stable spatial
patterns by the dynamics which occurs for relatively long time scales so we
expect our results to be largely independent on the quench protocol provided
that it is fast enough. Indeed, for zero disorder, we have checked, see
Appendix B, that a quench in the coupling constant leads to qualitatively
similar results for sufficiently long times. Moreover, since we aim to compare
with the BCS dynamics for short times, our quench results in a superconducting
state which for no disorder is still spatially homogeneous in the $T\lesssim
T_{c}$ region.

##  III Quench dynamics of the order parameter: initial exponential growth,
time oscillations and its eventual suppression

We now proceed to study the dynamics of the condensate amplitude triggered by
lowering the temperature of the system from $T_{i}>T_{c}$ to $T_{f}<T_{c}$
using the quench protocol Eq. (5). Since one of our main goals is the modeling
of stable spatial patterns of the order parameter in the long time limit [27],
we use the self-consistent time dependent BdG equations introduced earlier
which results in an initial spatially homogeneous evolution of the order
parameter. We shall see that a weak disordered potential does not change this
picture substantially. This is indeed a welcome feature as another aim of the
paper is to compare for short times our results with previous theoretical
predictions using the simpler BCS approach that cannot account for spatial
inhomogeneities.

Figure 1: The dynamics of the spatially averaged order parameter
$\langle\Delta(\bm{r})\rangle$ (black dot) together with an exponential
fitting $f(t)=a_{f}\exp(b_{f}t)$ (red line) at short times corresponding to
temperatures slightly below the critical one that marks the transition into
the superconducting state.

We first study the initial growth of the condensate as the system enters the
superconducting phase. We observe that the amplitude of the order parameter
increases rapidly from zero as the system enters the superconducting phase by
lowering the temperature $T<T_{c}$. A careful fitting of the numerical
results, see Appendix A, indicates that the growth of the order parameter is
exponential in this region. This confirms that the quench is fast enough to
induce a highly non-adiabatic time evolution. An exponential growth is also
observed, see Figure 1, in a relatively broad range of disorder strengths with
a growth rate, $0.86\sim 0.89$, that is not very sensitive to disorder and it
is also similar ($0.86$) to that found in the clean limit.

Figure 2: The real time evolution, after spatial average, of the order
parameter $\langle\Delta(\bm{r})\rangle$ in the presence of a random potential
for different values of the disorder strength $V$. The system size is
$N=200\times 200$, the coupling constant $U=-3$ and the chemical potential
$\mu=-0.34$ corresponding to a mean charge density $\langle n\rangle\simeq
0.875$. $\Delta_{0}$ is the value of the order parameter in the clean limit at
zero temperature.

The initial exponential growth of the spatial averaged order parameter
amplitude, obtained from the solution of the BdG equations Eq. (1), is
followed, see Figure 2, by relatively simple oscillations in time. At this
early stage, we did not yet observe spatial inhomogeneity in the clean limit,
or even in the presence of a weak disordered potential $V\leq 0.1$. The time
evolution of the order parameter within the simple BCS formalism [21, 19, 26,
22, 23], that by design neglects spatial inhomogeneities, also shows
oscillations in time whose details depend on both the initial state and the
quench protocol. More specifically, for an initial state of the order
parameter characterized by uncorrelated phases in momentum space [26] or with
a very small initial value [21], which simulates the system above the critical
temperature $T_{c}$, the lowering of the temperature below $T_{c}$ induces
undamped periodic oscillations in the amplitude of the order parameter. By
contrast, for a quench in the coupling constant at zero temperature, and
therefore inside the superconducting phase [22, 21, 23], the amplitude of the
order parameter oscillations can decay either as a power-law or exponential
way, or not decay at all, depending on the values of the initial and final
coupling constants. As expected, we have found excellent agreement with the
BCS prediction for the protocols that we have tested explicitly, see Appendix
B. This is not surprising as BCS theory and BdG theory should agree in the
limit in which the order parameter is spatially homogeneous.

The dynamics becomes more interesting for longer times. Results depicted in
Figure 2 indicate that this simple dynamical regime ends rather abruptly due
to the sharp exponential suppression of the amplitude of these oscillations,
which does not occur in the BCS dynamics [21, 19, 26, 22, 23], see also
Appendix B. The time scale of this suppression is sensitive to the addition of
a weak random potential. For a stronger disorder still deep in the metallic
phase, the oscillations are almost fully suppressed after a few periods. We
show next, by employing the mentioned time dependent BdG formalism, that the
origin of that exponential suppression lies in the development, even for no
disorder, of spatial inhomogeneities in the order parameter.

In order to carry out a more quantitative analysis of the decay of the
amplitude of the order parameter, we define
$\delta\Delta=(\langle\Delta(\bm{r})\rangle_{\rm
peak}-\langle\Delta(\bm{r})\rangle_{\rm valley})/\Delta_{0}$, see Figure 3,
where peak and valley refer to consecutive local maxima and minima of the
oscillations in time of the order parameter. For sufficiently short times,
where the order parameter is spatially homogeneous, the reduction of the
amplitude $\delta\Delta$ is consistent with a power-law decay. A fit with a
power-law decaying function $\sim t^{-\gamma_{p}}$ yields $\gamma_{p}\sim 1.4$
for weak or no disorder which illustrates that this slow decay is not related
to the presence of a random potential.

The power-law decay is followed by a much faster exponential decay
$\sim\exp(-\gamma_{e}t)$ even in the clean limit. As is expected, $\gamma_{e}$
becomes larger as disorder increases. Even for a relatively weak disorder,
$V=0.1$, it is already about two times larger that in the clean limit.
Qualitatively, the dependence of the crossover time between power-law and
exponential decay on the disorder strength $V\ll 1$ seems to be rather weak
which reinforces the idea that disorder does not play a leading role in this
phenomenon. We investigate in more detail this region of exponential
suppression in the next section.

Figure 3: Left y axis: the amplitude of the time oscillation $\delta\Delta$
(black dots) of the order parameter in the presence of different disorder
strengths $V$ and the corresponding power-law $y\propto t^{-\gamma_{p}}$ (red
line) and exponential fittings $y\propto\exp(-\gamma_{e}t)$ (blue line). Right
y axis: the variance of the order parameter Var$[\Delta(\bm{r})]$. We combine
both results using a double y axis plot in order to show explicitly that the
exponential suppression of time oscillations is induced by the exponential
growth of spatial inhomogeneities. Only when the spatial inhomogeneities
become sufficiently large, because of the exponential growth of the variance,
the exponential suppression of the time oscillations occurs.

##  IV Exponential growth of emergent spatial inhomogeneities and exponential
suppression of time oscillations

We now carry out a comprehensive study of the full form of the exponential
decay of the order parameter oscillations in time aiming to relate this
exponential suppression with the emergence of spatial inhomogeneities in the
order parameter even in the absence of disorder.

Figure 4:  The time dependence of the spatially averaged order parameter
$\langle\Delta(\bm{r})\rangle$ (blue line) in the clean $V=0$ and weak
disorder regions $V=0.1$. The red dash-dot line is the best fitting, by Eq.
(6), in the region of exponential decay of the amplitude of the time
oscillations of the order parameter that follows the region of slower power-
law decay. As is observed in the inset plots of the spatial distribution of
the order parameter at times $t_{1},t_{2},t_{3}$, this time scale is related
to the formation of spatial inhomogeneities. As was expected, the fitting
parameter $\gamma$ which characterizes the exponential decay increases with
disorder. The system size is $N=200\times 200$, the coupling constant $U=-3$,
and the chemical potential $\mu=-0.34$.

As a first step, we employ a simple oscillatory function with an exponential
decay of the amplitude [27, 19] to fit the time dependent spatial averaged
order parameter obtained from the BdG formalism,

$\Delta(t)=\Delta_{f}-A(\Delta_{f}+C+\cos(\omega(t-t_{0})\Delta_{0}))/\exp(\gamma(t-t_{0})\Delta_{0}),$
(6)

where $\Delta_{0}\sim 0.83$ is the value of the order parameter in the clean
limit at zero temperature, and $\Delta_{f}$ is the order parameter in the
final equilibrium state after the quench. The four fitting parameters are
$\gamma$, the decay ratio of the amplitude, $A$, $C$ and $\omega$. As is shown
in Figure 4, for times right after the crossover to an exponential decay in
$\delta\Delta$ (see Figure 3) , we find a very good agreement with the BdG
results in the clean and weak disorder limit. We note that in the BCS
approach, the value of $\omega$ is sensitive to the initial state and the
quench protocol [21, 19, 40]. In our case, the fitting yields $\omega\sim
0.3\Delta_{0}$ which is in the same ballpark as the BCS prediction [40] for a
initial state characterized by a very small order parameter with respect to
$\Delta_{0}$. In any case, we do not expect quantitative agreement because
this frequency may also be affected by the emergence of spatial
inhomogeneities.

For earlier times, as expected, the decay of time oscillations is much slower
so the fitting is much worse. Moreover, we find that this exponential decay,
even for no disorder, seems to be closely related with the emergence of
spatial inhomogeneities in the order parameter. This can be seen from the
similarity between the time scale $t_{2}$ in which the exponential suppression
of oscillations occurs and the time scale $t_{3}$ in which spatial
inhomogeneities, already existent for $t_{2}$, become substantial, see the
insets of Figure 4. In other words, the emergence of spatial inhomogeneities
eventually, namely, when they are large enough, triggers the exponential decay
of oscillations in time that terminates approximately when the spatial
inhomogeneities are fully formed.

Figure 5: Variance of the spatial inhomogeneities of the amplitude of the
order parameter $\Delta(\bm{r},t)$ for different disorder strengths $V$. The
starting of the exponential growth observed for intermediate times shows much
earlier than the beginning of the exponential suppression of the oscillations
of the averaged order parameter.

In order to establish a more quantitative relation between spatial
inhomogeneities and the exponential decay of time oscillations, we compute the
variance of the order parameter
Var$|\Delta(\bm{r})|=\langle\Delta^{2}(\bm{r})\rangle-\langle\Delta(\bm{r})\rangle^{2}$
as a function of time in the clean limit and in the absence of weak disorder
$V\leq 0.1$. Results depicted in Figure 5 show a region of intermediate times
where the variance grows exponentially. We define $t_{e}$ as time in which
this exponential growth terminates because the spatial patterns are completely
developed, see Figure 5.

Another interesting feature of the time dependence of the spatial variance is
the observation of a period of no growth, only for no disorder $V=0$, right
after the initial growth of the condensate. This feature strongly suggests
that the later exponential growth is independent of the quench protocol. For
$V\neq 0$, the situation is different. The early exponential growth in time of
the amplitude of the order parameter is followed by the exponential growth of
the variance which indicates that the seeds of spatial inhomogeneities due to
disorder are amplified exponentially by the dynamics of the BdG
superconductor. The flat behavior of the variance for later times, after the
exponential growth, confirms the earlier claim that the end of the time
oscillations is related to the approach to a quasi equilibrium state where the
change with time of the order parameter is heavily suppressed.

We have now all ingredients to compare explicitly to what extent the
exponential growth of the spatial inhomogeneities and the crossover between
power-law and exponential decay of oscillations in time are closely related.
For that purpose, we depict in Figure 3, back to back, the plots of
$\delta\Delta(t)$, that characterizes the amplitude of time oscillations, and
Var$[\Delta(r)]$, representing the exponential growth of the variance of the
spatial inhomogeneities. For both, the clean and the weak disordered case,
$V\leq 0.1$, the crossover to an exponential decay of the oscillations
precisely occurs when the exponential growth of the variance of spatial
fluctuations is close to its termination, namely, when it has reached a value
sufficiently large to have an impact on the quench dynamics. Therefore,
$t_{d}$ and $t_{e}$ has a very similar dependence on disorder and $t_{d}$ is a
bit larger than $t_{e}$ in all cases as the exponential suppression does not
require a full development of spatial inhomogeneities described by $t_{e}$.
These results fully confirm that the oscillations in time of the order
parameter are eventually suppressed exponentially by the emergence of spatial
inhomogeneities characterized by a variance that grows exponentially. This is
therefore a quite robust feature of the non-adiabatic dynamics of BdG
superconductors.

We turn to a quantitative investigation of the dependence on disorder of these
results. For that purpose, we show in Figure 6, the dependence on the strength
of disorder $V$ of $t_{d}$, the typical crossover time between power-law and
exponential decay of $\delta\Delta(t)$, and $t_{e}$ the time at which the
exponential growth of the variance terminates. We observe that in the weak
disorder limit $V\leq 0.1$ that is of interest, both typical times have not
only similar values but also a simple logarithmic dependence with the disorder
strength, with a finite value for $V=0$, confirming, at least for weak
disorder, that the exponential decay of the order parameter time oscillations
has its origin in spatial inhomogeneities induced by the quench dynamics with
the disorder potential playing the secondary role, at least in this region, of
shifting the development of these spatial pattern to earlier times. We note
that the slightly larger value of $t_{e}$ is expected because the effect of
the spatial inhomogeneities must be felt earlier, but not much earlier because
the growth is exponential, than the time at which the spatial inhomogeneities
are fully formed.

Figure 6: Typical crossover time $t_{d}$ (blue rhombs), see Figure 3, between
the power-law and exponential decay of the amplitude of the time oscillations
of the order parameter as a function of the disorder strength $V$. Typical
time $t_{e}$ (red rhombs) signaling the end of the exponential increase of the
variance of the order parameter spatial inhomogeneities depicted in Figure 5.
The solid lines are fits using the function ${\rm time}=a_{d}\log(b_{d}\times
V+c_{d})$. We not only observe a similar logarithmic dependence in both cases,
with a finite value at $V=0$, but $t_{d}$ and $t_{e}$ have similar values.
Therefore, the exponential growth of the spatial inhomogeneities induces the
late exponential decay of the amplitude of the time oscillations.

We now move to the quantitative description of the spatial patterns for
sufficiently long times when the system reaches the final quasi equilibrium
state.

##  V Long times: spontaneous formation of spatial patterns resulting from the
quench dynamics

The conclusion of the previous section is that, within the time dependent BdG
formalism, time oscillations of the order parameter are exponentially
suppressed for sufficiently long times. This is markedly different from the
BCS result in which this specific exponential suppression does not occur,
because the order parameter is spatially homogeneous.

Figure 7: Left column: the time dependence of the spatial averaged order
parameter $\langle\Delta(\bm{r})\rangle$ for different disorder strengths $V$.
Central column: spatial distribution of the order parameter at time $t_{4}$,
see left column, corresponding with a quasi-equilibrium state with intricate
spatial patterns and a very weak time dependence. Right column: structure
factor Eq. (7) that reveals patterns in spatial distribution of the order
parameter at time $t_{4}$. The parameters are: system size $N=200\times 200$,
coupling constant $U=-3$ and the chemical potential $\mu=-0.34$.

In this section, we present a detailed description of those stable,
spontaneously formed, spatial patterns, see Figure 7, together with its
formation process, see Figure 8, as a function of the disorder strength.
Videos of the full time evolution are available here [41]. We recall that for
no disorder, spatial inhomogeneities start to appear when the time
oscillations are significantly suppressed. Initially, see left column of
Figure 8, they resemble soft periodic or quasi-periodic domains where the
order parameter amplitude is substantially smaller than in the surroundings.
For longer times, these domains become more pronounced and adopt the shape of
relatively thin finite-size stripes, most times organized in a perpendicular
fashion. These broken stripes have well defined centers where the suppression
of the order parameter is even stronger. This seems to be an equilibrium or
quasi-equilibrium state because we do not appreciate further changes even for
much larger time scales. A small but finite disorder does not change much the
emergence of these spatial patterns.

Figure 8: The spatial distribution of the order parameter $\Delta(\bm{r})$ at
times $t_{1},t_{2}$ and $t_{3}$, defined in Figure 7, corresponding to
different stages of the time evolution before the full quasi-equilibrium state
($t_{4}$), and for different values of the disorder strength $V$.

In order to have a more quantitative description of this spatial structure, we
have computed the averaged gap correlation function
$\langle\Delta(\bm{r})\Delta(0)\rangle$. Results depicted in Figure 9 show
oscillations with a typical length $\ell_{p}\sim 12.5$ in the clean limit,
which is much larger than the superconducting coherence length at thermal
equilibrium. This is a strong indication that the pattern of spatial
inhomogeneities, especially its periodicity, is due to the quench dynamics.

Figure 9:  The static order parameter correlation function
$\langle\Delta(\bm{r})\Delta(0)\rangle$ at time $t_{4}$, defined in Figure 7,
normalized by $\langle\Delta(\bm{r}=0)\Delta(0)\rangle$, the spatial average
of the square of the order parameter.

As a further probe of the periodic nature of the spatial structure, we compute
the structure factor,

$S(\bm{q})=\sum_{ij}\Delta(\bm{r}_{i})\Delta(\bm{r}_{j})\exp(i\bm{q}(\bm{r}_{i}-\bm{r}_{j}))/\sum_{i}\Delta^{2}(\bm{r}_{i}),$
(7)

of the distribution of the order parameter for long times where the system
seems to have reached a quasi equilibrium state. We termed quasi-equilibrium
state for two reasons, there is still some residual time dependence and also
the resulting state is very different from that corresponding to the solution
of the static BdG equations. As is shown in Figure 7, the spatial patterns
resulting from the out of equilibrium dynamics have a square crystal-like
structure. Moreover, according to the Bragg pattern depicted in the right
column of Figure 7, the distance from the center to the peaks is around $0.25$
in momentum space which in real space corresponds to a typical length of the
strips, or the fake vortex lattice, of about $\ell_{p}\sim\pi/0.25=12.56$.
This is consistent with the previous findings for the order parameter spatial
correlation function in Figure 9. The addition of a weak disorder potential
does not modify substantially these results.

Figure 10: The variation of the charge density $\delta
n(\bm{r})=n(\bm{r},t)-n(\bm{r},t=0)$ at time $t_{4}$, defined in Figure 7. The
spatial pattern of suppressed charge density matches that of the
superconducting order parameter, especially when disorder is not very strong.

Results depicted in Figure 10 for the charge density variation $\delta
n(\bm{r})$ further confirm the existence of spatial patterns consistent with
the one found for the spatial structure of the order parameter. Therefore, the
observed stripe-like domains where the order parameter is suppressed cannot be
attributed to quantum coherence effects [42, 37, 43], but are rather related
to modulations of the charge density caused by the strong suppression of time
oscillations.

###  V.1 From fragmented stripes to fake vortices

As disorder strength increases, but still deep in the metallic region, we
observe a similar periodic pattern that does not yet seem much affected by
disorder. Although the periodicity is robust, see Figure 7, the shape of the
domains where the order parameter is suppressed undergoes a gradual change.
The stripes becomes more fragmented until they become completely disconnected.
However, a rather strict periodicity of these patterns is still observed even
for substantially larger disorder $V\leq 0.5$. The regions where the order
parameter is suppressed are now circular and with a typical length that is
much larger than the superconducting coherence length. As we said earlier, the
phase of the order parameter is fixed in this region, so the periodic
suppression of the order parameter is not related to non-trivial topological
properties. For that reason, we have termed this spatial pattern fake vortex
lattice. Both, the correlation function
$\langle\Delta(\bm{r})\Delta(0)\rangle$, see Figure 9, and the structure
factor of the spatial distribution of the superconducting order parameter, see
Figure 7, confirm this lattice structure of the spatial patterns.

For $0.5<V<2$, disorder effects become gradually more pronounced. The periodic
pattern of the fake vortex is gradually deformed though clear vortex repulsion
is still observed even when $V=1.0$. The latter is confirmed by an explicit
calculation of the structure form factor, see Figure 7, that still shows a
clear circular pattern. This is also confirmed in the order parameter
correlation function $\langle\Delta(\bm{r})\Delta(0)\rangle$, see Figure 9,
where the vestiges of a lattice structure resulting in vortex repulsion show
up as a valley at $|\bm{r}|=8$ and a peak at $|\bm{r}|=16$ for a relatively
large strength $V=1.0$ of the random potential.

In the strong disorder limit, close or around the transition $V\sim 2$, the
spatial distribution of the order parameter is mostly controlled by disorder.
It becomes highly inhomogeneous, has no specific periodic pattern but it seems
that some vortex-like structures remain. For instance, the gap correlation
function $\langle\Delta(\bm{r})\Delta(0)\rangle$, see Figure 9, only has short
correlations in space and the Fourier transform has no Bragg peaks signaling
no periodic pattern in space.

##  VI Discussion

A natural question to ask is the origin of the crossover from stripes to fake
vortices that we observe as the strength of disorder is increased. Since
disorder suppresses finite size effects, we cannot rule out that the so called
fake vortices, that arise for stronger disorder, should still be observed for
weak or no disorder, instead of the fragmented stripes, if it were possible to
reach much larger lattice sizes. Unfortunately, substantially larger sizes are
currently beyond the reach of our computing capabilities.

Another important issue is the dependence of the results to the strength of
the electron-phonon coupling. We have set the coupling constant in the strong
coupling region $U=-3$, and set the Debye energy to the full band of the
energy spectrum, in order to be able to explore quantitatively the rich
spatial structure of the order parameter resulting from the out of equilibrium
dynamics whose typical length is around $\ell_{p}\sim 12.56$, which is much
larger than the superconducting coherence length. As is expected, this is not
fully consistent with the weak-coupling analytic prediction [27] that this
typical length should be the superconducting coherence length. Results for an
even stronger coupling constant $U=-5$, presented in the Appendix D, show
similar spatial patterns but, as expected, the typical length of the patterns
is shorter. A more quantitative understanding on the precise nature of this
length scale would require a more systematic, and therefore numerically
costly, analysis of its dependence on the quenched dynamics, for instance by
varying the coupling constant, which is beyond the scope of the paper.

More specifically, it may be interesting to explore the weak coupling region
$|U|\leq 1$ where quantum coherence effects induce multifractal-like features
for intermediate to strong disorder. However, we anticipate that typically
strong quantum coherence effects such as multifractality will suppress the
spatial patterns induced by the quench dynamics even for a relatively weak
disorder strength $V\leq 0.5$. This bring us to the issue of the experimental
confirmation of these results. The strong coupling region we have explored in
this paper is more amenable for experiments with Bose-Einstein condensate at
very low temperature where the dynamics is triggered by a change in the
coupling constant which is feasible to carry out in this setting. Disorder-
like effects in this problem can be modeled by quasi-periodic optical lattice
configurations.

We note that in superconducting materials, a quench protocol based on the
change of coupling constant is not in principle possible but even a controlled
quench in temperature is challenging. This is why some experiments opted to
induce out of equilibrium dynamics by bombarding the sample with photons
leading to heating and subsequent cooling [44]. However, the theoretical
modeling of such systems is well beyond the mean-field approach, and
relatively simple quench protocols, that we are considering here.

Finally, we address the dependence of the results on the quench speed. Results
depicted in Appendix C for a much slower quench, points to a more nuanced
picture. For sufficiently long times, spatial inhomogeneities are clearly
observed but they do not have a periodic pattern. Large domains with similar
values of the order parameter are separated by filamentary domain walls where
the order parameter is highly suppressed. The domains walls becomes thinner
for longer times but they persist in the range of times we can explore
numerically. We therefore expect that a sufficiently slow quench will in
principle lead to an essentially adiabatic dynamics where early time
oscillations are suppressed and spatial patterns may not develop at all. We
also note that sufficiently fast temperature quenches inducing far out of
equilibrium effects are beyond the finite temperature mean-field formalism
that we employ in the paper.

##  VII Conclusions

We have investigated the quenched dynamics of a BdG superconductor. In order
to compare with previous results, which use the simpler, spatially
homogeneous, BCS formalism, we have employed two quench protocols: an abrupt
change in the coupling constant, see Appendix B, and a smooth linear drop in
temperature starting above the critical temperature. For zero disorder and
sufficiently long times, where the BCS approach ceases to be applicable, we
have obtained similar results so the study of the role of disorder was carried
out only for the second quench protocol.

For short times, we observe similar results as in the simpler BCS approach,
the order parameter remains homogeneous, it first grows exponentially and then
oscillates in time with a pattern that depends on the quench protocol and the
initial state. However, in contrast with previous BCS findings, the amplitude
of these time oscillations eventually decreases exponentially in time because
of the emergence of spatial inhomogeneities of the order parameter. We have
characterized the emergence of these spatial inhomogeneities by the
exponential growth of the variance in space of the order parameter and showed
that this exponential growth in space causes the exponential suppression of
the oscillations in time. This feature cannot be accounted for in the BCS
formalism and it is observed even in the limit of no disorder. A weak
disordered potential does not change qualitatively these results.

For longer times, these spatial instabilities turn into rich spatial patterns.
In the clean limit, or for sufficiently weak disorder, we observe ordered
filamentary structures resembling finite size stripes in closed perpendicular
directions, where the order parameter is heavily suppressed. The suppression
of the order parameter in the central region of these broken stripes is much
more pronounced. For stronger disorder, but still deep in the metallic region,
the stripe-like structures morph into a square lattice of fake vortices,
namely, the amplitude of the order parameter is heavily suppressed in circular
regions whose typical length is much larger than the superconducting coherence
length, but the phase has no topological properties so no real vortex is
formed. Larger sizes, beyond our current numerical capabilities, are needed to
clarify whether the observed broken stripes for weak or no disorder become
also a vortex lattice in the thermodynamic limit. A further increase of
disorder leads to a gradual deformation of the lattice, though isolated fake
vortices that repel each other are still clearly observed. Finally, close to
the insulating transition, spatial inhomogeneities induced by disorder become
also relevant leading to a quite intricate spatial structure that is
characterized by a lack of vortex lattice symmetry but persistence of vortex
repulsion. Although the time scale we can simulate numerically is limited, our
results suggest that, neglecting the effect of collisions beyond the time
dependent BdG formalism, these spatial structures correspond with a quasi-
equilibrium state of the superconductor which is still very different from
that corresponding with the solution of the static BdG equations. The observed
spatial patterns are confirmed by a careful analysis in Fourier space based on
the calculation of the structure factor.

It would be interesting to gain a more comprehensive understanding of the
precise nature of the spatial patterns, especially in the so called stripe
region, either by numerical or analytic techniques, in order to determine
whether, for sufficiently large sizes, the spatial structures that emerge for
intermediate times lead to a full checkerboard-like shape instead of the
observed broken stripes. It would also be worthwhile to extend these results
to p-wave and d-wave superconductors in order to explore its potential
relevance in topological superconductivity and the physics of cuprates
superconductors.

Acknowledgments

We acknowledge financial support from a Shanghai talent program and from the
National Key R&D Program of China (Project ID: 2019YFA0308603), China
Postdoctoral Science Foundation (Grant numbers 2023M732256, 2023T160409,
GZB20230420). We acknowledge Prof. Zi Cai for helpful and illuminating
discussions.

##  Appendix A The initial growth of the condensate

In this appendix, we study the initial growth of the superconducting order
parameter when the temperature is just lowered below the critical temperature.
For that purpose, we fit numerical results from the solution of the BdG
equation with Gaussian, exponential and power-law test functions. The outcome
of the fitting, see Figure 11, is that the order parameter increases
exponentially during the early stages of the time evolution in the
superconducting phase. This is consistent with previous results in
phenomenological models for sufficiently fast quenches [9].

Figure 11: The dynamics of the spatially averaged order parameter
$\langle\Delta(\bm{r})\rangle$ (black dot) for times right after the system
has entered in the superconducting phase. The results of the different
fittings using Gaussian, power-law and exponential functions clearly indicate
that the latter is the one closer to the numerical results.

##  Appendix B  Comparison between the dynamic BCS and BdG results in the
sudden coupling quench case

In this appendix, we compare the quench dynamics using the BCS and BdG
formalism. We cannot use the quench in temperature of the main manuscript
because it is difficult to model it in the BCS formalism which cannot, at
least in a fully self-consistent way, account for spatial inhomogeneities.
Therefore, we compare the quench dynamics in BCS and BdG by using a quench
protocol with an abrupt change in the coupling constant at zero temperature
where the initial state is spatially homogeneous and therefore it is possible
to model it with both approaches. Results depicted in Figure 12 show that in
the beginning, when the spatial dependence is not yet observable, the BCS and
BdG quench dynamics is quantitatively very similar. However, at time
$t\Delta_{f}$ around 400, we start observing the spatial structures in the BdG
results which triggers a sharp departure from the BCS prediction. Those
results, together with those presented in the main text corresponding the
dynamics after lowering the temperature into the superconducting phase, show
that the emergent spatial structure as a consequence of the quench dynamics is
rather universal and independent on the quench protocol.

Figure 12: The time evolution of the order parameter $\Delta(t)/\Delta_{0}$ in
the clean limit obtained with the BCS (blue) and BdG (red) formalisms. The
insets show (left) the time range when the BdG result start deviating from the
BCS one and (right) the spatial distribution of the order parameter at the
time marked by the black arrow. The system is prepared in a initial state with
coupling constant $|U|=1$. The quench protocol consists in an abrupt change of
the coupling constant from $|U|=1$ to $|U|=3$. The other parameters are system
size $N=40\times 40$ and chemical potential $\mu=0$.

##  Appendix C  Dynamics resulting from slow quenches

The main text provides a detailed analysis of the time evolution of the system
under a fast quench protocol. It is expected that the results for a
sufficiently slow quench would be somehow less interesting as the dynamics
would be essentially adiabatic corresponding to a slow, and largely
homogeneous, at least initially, growth of the order parameter as temperature
is gradually lowered below the critical temperature. For the sake of
completeness, in this appendix, we present results for the dynamics, and
spatial pattern formation, for a slow quench characterized by $\tau_{Q}=500$.
We find striking similarities with the fast quench case but also important
differences. Results depicted in Figure 13 in the very weak disorder limit
$V=0.001$ show that the observation of time oscillations is greatly delayed
with respect to the fast quench limit but, once they occur, they are
qualitatively similar with a fast suppression once the filamentary spatial
patterns are fully developed. Interestingly, the crossover from the broken
stripe to the fake vortex phase is not clearly observed though it may be due
to the limited time scale we have explored numerically or simply to the weak
strength of the random potential. In the presence of a stronger disorder
$V=0.5$, see Figure 14, there is no visible time oscillations. The shape of
spatial domains become irregular, and domain walls become thinner over time,
but persist within the range of times that we can explore numerically. This
rather different, with respect to the fast quench case, spatial structure
deserves further exploration.

(a) $t_{1}$

(b) $t_{2}$

(c) $t_{3}$

(d) $t_{4}$

Figure 13: Top: The quenched dynamics of the spatially averaged order
parameter $\langle\Delta(\bm{r})\rangle$ in the presence of a random potential
with disorder strength $V=0.001$. Bottom: the spatial distribution of the
order parameter at times $t_{1},t_{2},t_{3}$ and $t_{4}$ defined in the top
plot.

(a) $t_{1}$

(b) $t_{2}$

(c) $t_{3}$

(d) $t_{4}$

Figure 14: Top: The slow quench dynamics of the spatially averaged order
parameter $\langle\Delta(\bm{r})\rangle$ in the presence of a weak disorder of
strength $V=0.5$. Bottom: the spatial distribution of the order parameter at
times $t_{1},t_{2},t_{3}$ and $t_{4}$ defined in the top plot.

##  Appendix D  Dynamics in the fast quench, strong coupling limit $U=-5$

In this appendix, we present results for the quench dynamics in the clean
limit ($V=0$) of the order parameter, using the protocol of the main text
(fast quench), in the region of stronger coupling constant $U=-5$. We observe,
see Figure 15, qualitatively similar features as for the $U=-3$ case studied
in the main text. The spatial averaged order parameter increases exponentially
when temperature is below $T_{c}$, and then exhibits damped oscillations in
time. The development of spatial patterns in the order parameter also followed
a similar path: the sharp growth of the spatial inhomogeneities occurs around
the time in which time oscillations are fully suppressed and a broken stripe
phase is followed by the fake vortex phase though the size of these structures
is substantially smaller than for $U=-3$. For a more quantitative
understanding of this size difference, we study the order parameter
correlation $\langle\Delta(\bm{r})\Delta(0)\rangle$ function and the structure
factor at time $t_{4}$ in Figure 15(d), when the order parameter is almost at
equilibrium, namely, it only experiences very small, non harmonic,
oscillations due to the residual collective behavior of Cooper pairs. The
results, depicted in Figure 16, show that
$\langle\Delta(\bm{r})\Delta(0)\rangle$ undergoes oscillations in space.
Likewise, the structure factor reveals the existence of a square lattice of
fake vortices with typical length $\ell_{p}\sim 5.6$ in real space. This
typical length is much shorter than that for $U=-3$, but still much larger
than the superconducting coherence length which, for this value of the
coupling, is of the order of the lattice spacing.

(a) $t_{1}$

(b) $t_{2}$

(c) $t_{3}$

(d) $t_{4}$

Figure 15: Top: The time evolution of the spatial averaged order parameter
$\langle\Delta(\bm{r})\rangle$ in the clean limit, $V=0$. Bottom; The spatial
distribution of the order parameter at the corresponding times
$t_{1},t_{2},t_{3}$ and $t_{4}$ shown in the top figure. The system size is
$N=200\times 200$, the coupling constant $U=-5$ and chemical potential
$\mu=-0.45$.

Figure 16:  Left: The order parameter correlation function
$\langle\Delta(\bm{r})\Delta(0)\rangle$ at the equilibrium time, $t_{4}$ in
Figure 15(d) which is normalized by $\langle\Delta(\bm{r}=0)\Delta(0)\rangle$.
Right: The structure factor at quasi equilibrium time $t_{4}$, in Figure
15(d). The Bragg’s pattern shows a distance to the peaks of around $\pi/5.6$
in momentum space, which corresponds to a typical length of the lattice of
about $\ell_{p}\sim 5.6$ while the superconducting coherence length is much
smaller, about the lattice spacing.

## References

  * [1] WY Córdoba-Camacho, RM Da Silva, AA Shanenko, Alexei Vagov, AS Vasenko, BG Lvov, and J Albino Aguiar.  Spontaneous pattern formation in superconducting films.  Journal of Physics: Condensed Matter, 32(7):075403, 2019.
  * [2] Zijian Wang, Carlos Navarrete-Benlloch, and Zi Cai.  Pattern formation and exotic order in driven-dissipative bose-hubbard systems.  Physical Review Letters, 125(11):115301, 2020.
  * [3] Tom WB Kibble.  Some implications of a cosmological phase transition.  Physics Reports, 67(1):183–199, 1980.
  * [4] Wojciech H Zurek.  Cosmological experiments in superfluid helium?  Nature, 317(6037):505–508, 1985.
  * [5] Giacomo Lamporesi, Simone Donadello, Simone Serafini, Franco Dalfovo, and Gabriele Ferrari.  Spontaneous creation of kibble–zurek solitons in a bose–einstein condensate.  Nature Physics, 9(10):656–660, 2013.
  * [6] S Ulm, J Roßnagel, G Jacob, C Degünther, ST Dawkins, UG Poschinger, R Nigmatullin, A Retzker, MB Plenio, F Schmidt-Kaler, et al.  Observation of the kibble–zurek scaling law for defect formation in ion crystals.  Nature communications, 4(1):2290, 2013.
  * [7] Guillaume Labeyrie and Robin Kaiser.  Kibble-zurek mechanism in the self-organization of a cold atomic cloud.  Physical Review Letters, 117(27):275701, 2016.
  * [8] Timothy Chi Hin Liew, OA Egorov, M Matuszewski, Oleksandr Kyriienko, X Ma, and EA Ostrovskaya.  Instability-induced formation and nonequilibrium dynamics of phase defects in polariton condensates.  Physical Review B, 91(8):085413, 2015.
  * [9] Paul M Chesler, Antonio M García-García, and Hong Liu.  Defect formation beyond kibble-zurek mechanism and holography.  Physical Review X, 5(2):021015, 2015.
  * [10] Han-Chuan Kou and Peng Li.  Interferometry based on the quantum kibble-zurek mechanism.  Physical Review B, 106(18):184301, 2022.
  * [11] Ze Jing, Huadong Yong, and Youhe Zhou.  Thermal coupling effect on the vortex dynamics of superconducting thin films: time-dependent ginzburg–landau simulations.  Superconductor Science and Technology, 31(5):055007, 2018.
  * [12] Antonio Lara, César González-Ruano, and Farkhad G Aliev.  Time-dependent ginzburg-landau simulations of superconducting vortices in three dimensions.  Low Temperature Physics, 46(4):316–324, 2020.
  * [13] Bakhrom Oripov and Steven M Anlage.  Time-dependent ginzburg-landau treatment of rf magnetic vortices in superconductors: Vortex semiloops in a spatially nonuniform magnetic field.  Physical Review E, 101(3):033306, 2020.
  * [14] Antonio M. García-García, Hua Bi Zeng, and Hai-Qing Zhang.  A thermal quench induces spatial inhomogeneities in a holographic superconductor.  Journal of High Energy Physics, 2014(7), jul 2014.
  * [15] Rupert L Frank, Christian Hainzl, Benjamin Schlein, and Robert Seiringer.  Incompatibility of time-dependent bogoliubov–de-gennes and ginzburg–landau equations.  Letters in Mathematical Physics, 106:913–923, 2016.
  * [16] Christian Hainzl and Jonathan Seyrich.  Comparing the full time-dependent bogoliubov-de-gennes equations to their linear approximation: a numerical investigation.  The European Physical Journal B, 89:1–10, 2016.
  * [17] John Bardeen, Leon N Cooper, and John Robert Schrieffer.  Theory of superconductivity.  Physical review, 108(5):1175, 1957.
  * [18] AF Volkov and Sh M Kogan.  Collisionless relaxation of the energy gap in superconductors.  Soviet Journal of Experimental and Theoretical Physics, 38:1018, 1974.
  * [19] RA Barankov, LS Levitov, and BZ Spivak.  Collective Rabi oscillations and solitons in a time-dependent BCS pairing problem.  Physical Review Letters, 93(16):160401, 2004.
  * [20] ShengQuan Zhou.  BCS Pairing Dynamics, 2006.  Physics 569 Emergent States of Matter, Fall 2006, at the University of Illinois at Urbana-Champaign.
  * [21] Emil A Yuzbashyan, Oleksandr Tsyplyatyev, and Bori s L Altshuler.  Relaxation and persistent oscillations of the order parameter in ferm ionic condensates.  Physical review letters, 96(9):097005, 2006.
  * [22] Emil A. Yuzbashyan and Maxim Dzero.  Dynamical vanishing of the order parameter in a fermionic condensate.  Phys. Rev. Lett., 96:230404, Jun 2006.
  * [23] Roman A Barankov and Leonid S Levitov.  Synchronization in the BCS pairing dynamics as a critical phenomenon.  Physical Review Letters, 96(23):230403, 2006.
  * [24] HP Ojeda Collado, Gonzalo Usaj, CA Balseiro, Damián H Zanette, and José Lorenzana.  Emergent parametric resonances and time-crystal phases in driven bardeen-cooper-schrieffer systems.  Physical Review Research, 3(4):L042023, 2021.
  * [25] HP Collado, Gonzalo Usaj, CA Balseiro, Damián H Zanette, and José Lorenzana.  Dynamical phases transitions in periodically driven bardeen-cooper-schrieffer systems.  arXiv preprint arXiv:2210.15693, 2022.
  * [26] R. A. Barankov and L. S. Levitov.  Dynamical selection in developing fermionic pairing.  Phys. Rev. A, 73:033614, Mar 2006.
  * [27] M Dzero, EA Yuzbashyan, and BL Altshuler.  Cooper pair turbulence in atomic fermi gases.  Europhysics Letters, 85(2):20004, 2009.
  * [28] KJ Challis, RJ Ballagh, and CW Gardiner.  Bragg scattering of cooper pairs in an ultracold fermi gas.  Physical review letters, 98(9):093002, 2007.
  * [29] Peng Zou and Franco Dalfovo.  Josephson oscillations and self-trapping of superfluid fermions in a double-well potential.  Journal of Low Temperature Physics, 177(5):240–256, 2014.
  * [30] Gia-Wei Chern and Kipton Barros.  Nonequilibrium dynamics of superconductivity in the attractive hubbard model.  Physical Review B, 99(3):035162, 2019.
  * [31] PGf de Gennes.  Boundary effects in superconductors.  Reviews of Modern Physics, 36(1):225, 1964.
  * [32] Renyuan Liao and Joachim Brand.  Traveling dark solitons in superfluid fermi gases.  Physical Review A, 83(4):041604, 2011.
  * [33] RG Scott, F Dalfovo, LP Pitaevskii, and S Stringari.  Rapid ramps across the bec-bcs crossover: A route to measuring t he superfluid gap.  Physical Review A, 86(5):053604, 2012.
  * [34] Jun Tokimoto, Shunji Tsuchiya, and Tetsuro Nikuni.  Josephson oscillation and self-trapping in a fermi superfluid gas a cross the bcs-bec crossover.  Journal of Low Temperature Physics, 208(5):372–378, 2022.
  * [35] H Pu, CK Law, S Raghavan, JH Eberly, and NP Bigelow.  Spin-mixing dynamics of a spinor bose-einstein condensate.  Physical Review A, 60(2):1463, 1999.
  * [36] Amit Ghosal, Mohit Randeria, and Nandini Trivedi.  Inhomogeneous pairing in highly disordered s-wave superconductors.  Physical Review B, 65(1):014501, 2001.
  * [37] Bo Fan and Antonio M. García-García.  Enhanced phase-coherent multifractal two-dimensional superconductivity.  Phys. Rev. B, 101:104509, Mar 2020.
  * [38] Bo Fan, Abhisek Samanta, and Antonio M García-García.  Tuning superinductors by quantum coherence effects for enhancing quantum computing.  Physical Review Letters, 130(4):047001, 2023.
  * [39] Julian Sonner, Adolfo Del Campo, and Wojciech H Zurek.  Universal far-from-equilibrium dynamics of a holographic superconductor.  Nature communications, 6(1):7406, 2015.
  * [40] G. L. Warner and A. J. Leggett.  Quench dynamics of a superfluid fermi gas.  Phys. Rev. B, 71:134514, Apr 2005.
  * [41] Time evolution of the order parameter after temperature quench for different disorder strengths v:.
$V=0:$
https://www.physics.sjtu.edu.cn/amgg/fakevortex_video/L200_quench_Ti1p2_Tf0p1_dT50_U3_V0p0_dt0p1_subplot.avi
$V=0.001:$
https://www.physics.sjtu.edu.cn/amgg/fakevortex_video/L200_quench_Ti1p2_Tf0p1_dT50_U3_V0p001_dt0p1_subplot.avi
$V=0.1:$
https://www.physics.sjtu.edu.cn/amgg/fakevortex_video/L200_quench_Ti1p2_Tf0p1_dT50_U3_V0p1_dt0p1_subplot.avi
$V=0.5:$
https://www.physics.sjtu.edu.cn/amgg/fakevortex_video/L200_quench_Ti1p2_Tf0p1_dT50_U3_V0p5_dt0p1_subplot.avi
$V=1:$
https://www.physics.sjtu.edu.cn/amgg/fakevortex_video/L200_quench_Ti1p2_Tf0p1_dT50_U3_V1p0_dt0p1_subplot.avi
$V=2:$
https://www.physics.sjtu.edu.cn/amgg/fakevortex_video/L200_quench_Ti1p2_Tf0p1_dT50_U3_V2p0_dt0p1_subplot.avi.

  * [42] Carmen Rubio-Verdu, Antonio M Garcia-Garcia, Hyejin Ryu, Deung-Jang Choi, Javier Zaldivar, Shujie Tang, Bo Fan, Zhi-Xun Shen, Sung-Kwan Mo, Jose Ignacio Pascual, et al.  Visualization of multifractal superconductivity in a two-dimensional transition metal dichalcogenide in the weak-disorder regime.  Nano letters, 20(7):5111–5118, 2020.
  * [43] Kun Zhao, Haicheng Lin, Xiao Xiao, Wantong Huang, Wei Yao, Mingzhe Yan, Ying Xing, Qinghua Zhang, Zi-Xiang Li, Shintaro Hoshino, et al.  Disorder-induced multifractal superconductivity in monolayer niobium dichalcogenides.  Nature Physics, 15(9):904–910, 2019.
  * [44] V. M. H. Ruutu, Vladimir Eltsov, A. J. Gill, A. J. Gill, Tom W. B. Kibble, Matti Krusius, Yuriy Makhlin, Yuriy Makhlin, Bernard Plaçais, Grigory E. Volovik, Grigory E. Volovik, and Wen-Bin Xu.  Vortex formation in neutron-irradiated superfluid 3he as an analogue of cosmological defect formation.  Nature, 382:334–336, 1996.


$\displaystyle y_{\text{pred}}=\text{sign}\;f(\bm{\theta},\mathbf{x}).$ (29)

The following is the quantum circuit used to evaluate overlaps in the
QuantumMetricLearner for two inputs
$\mathbf{x},\mathbf{x^{\prime}}\in\mathbb{R}^{4}$:

The number of circuits run for each training step and prediction scales
linearly with the size of the batch of example pairs used – in this paper we
fixed $32$ for both. Larger batch sizes allow a more reliable estimate of the
cost and predicted label.

hyperparameter | values
---|---
_learning_rate_ | [0.001, 0.01, 0.1]
_n_layers_ (embedding layers) | [1, 3, 4]

###  C.8 TreeTensorClassifier [19]

This model was designed to avoid the phenomenon of barren plateaus. We
implement the ‘tree tensor’ structure shown in Figure 1 of [19]. The
variational circuit in this model has a tree-like structure and therefore
requires a number of qubits that is a power of 2. In [19] one first optimizes
a variational circuit that finds a state that approximates an amplitude-
encoded data state. The reason for this is to make the algorithm more
efficient; here, to avoid an additional variational optimisation, we assume
direct access to the exact amplitude encoded state $V(\bm{x})\ket{0}$ (padding
with constant values $1/2^{n}$ with $n$ the number of qubits when necessary).

This state is then fed into the variational circuit $U(\bm{\theta})$
consisting of trainable single-qubit $Y$ rotations and CNOTs. The tree
structure means that there are few parameters; for a circuit with $n$ qubits
one has only $2n-1$ parameters.

Prediction is given by measurement of $Z$ on the first qubit,

$\displaystyle y_{\text{pred}}=\text{sign}f(\bm{\theta},\bm{x}),$ (30)

where

$\displaystyle
f(\bm{\theta},\bm{x})=\bra{0}V^{\dagger}(\bm{x})U^{\dagger}(\bm{\theta})Z_{1}U(\bm{\theta})V(\bm{x})\ket{0},$
(31)

and training is via the square loss:

$\displaystyle\ell(\bm{\theta},\bm{x},y)=(f(\bm{\theta},\bm{x})-y)^{2}.$ (32)

The following is an example of the quantum circuit used in the
TreeTensorClassifier for a $4$-dimensional input encoded into state
$\ket{\psi_{\mathbf{x}}}$:

This model has few suggested hyperparameters, so we vary only the learning
rate.

hyperparameter | values
---|---
_learning_rate_ | [0.001, 0.01, 0.1]

###  C.9 IQPKernelClassifier [52]

This model is a kernel equivalent of the IQP variational model. The embedding
$V(\bm{x})\ket{0}$ is the same IQP-inspired embedding given by PennyLane’s
IQPEmbedding template. This defines a kernel

$\displaystyle
k(\bm{x},\bm{x}^{\prime})=\text{tr}[\rho(\bm{x})\rho(\bm{x}^{\prime})]=|\bra{0}V^{\dagger}(\bm{x})V(\bm{x}^{\prime})\ket{0}|^{2},$
(33)

which we evaluate by applying the unitary
$V^{\dagger}(\bm{x})V(\bm{x}^{\prime})$ to an input state $\ket{0}$ and
calculating the probability to measure $\ket{0}$.

The following is an example of the quantum circuit used in the
IQPKernelClassifier for two inputs $\mathbf{x},\mathbf{x}\in\mathbb{R}^{4}$:

The kernel matrix $K$ is fed to scikit-learn’s SVC class, which trains a
support vector machine classifier. Prediction is given by

$\displaystyle
y_{\text{pred}}=\text{sign}\;\sum_{i}\alpha_{i}k(\bm{x}_{i},\bm{x}),$ (34)

with $\alpha_{i}$ the weights of the support vector machine.

The hyperparameter search values are as follows:

hyperparameter | values
---|---
_repeats_ (embedding layers) | [1, 5, 10]
_C_ (SVC regularization) | [0.1,1,10,100]

###  C.10 ProjectedQuantumKernel [54]

This model is a kernel method that uses a Hamiltonian-inspired data embedding
that resembles a Trotter evolution of a 1D-Heisenberg model with random
couplings. This consists of applying random single-qubit rotations to an input
state $\ket{0}$ of $n=d+1$ qubits, followed by $L$ layers of two-qubit
rotations with generators $XX$, $YY$, $ZZ$ to pairs of adjacent qubits, with
angles given by the elements of $\bm{x}$:

$\prod_{j=1}^{n}\exp\left(-\mathrm{i}\frac{t}{L}x_{j}\left(X_{j}X_{j+1}+Y_{j}Y_{j+1}+Z_{j}Z_{j+1}\right)\right),$
(35)

where $t$ is a hyperparameter of the model. Writing the embedded states as
$\rho(\bm{x}_{i})$, the kernel function $k(\bm{x}_{i},\bm{x}_{j})$ is

$\displaystyle\exp\left(-\gamma\sum_{k=1}^{n}\sum_{P\in\\{X,Y,Z\\}}\left(\text{tr}[P\rho_{k}(\bm{x}_{i})]-\text{tr}[P\rho_{k}(\bm{x}_{j})]\right)^{2}\right)$
(36)

where $\rho_{k}$ is the reduced state of $k^{th}$ qubit of $\rho$. This is
simply the RBF kernel with bandwidth $\gamma$ applied to feature mapped
vectors $\bm{\phi}(\bm{x})$ with elements $(\text{tr}[P\rho_{k}(\bm{x})])$.

The following is an example of the quantum circuit used to compute feature
vectors from measurements in ProjectedQuantumKernel for the input
$\mathbf{x}=(0.1,0.2,0.3,0.4)^{T}$:

According to [54], the default $\gamma$ value is set as

$\displaystyle\gamma_{0}=\frac{1}{\text{Var}(\bm{\phi})d}$ (37)

where $\text{Var}(\bm{\phi})$ is the variance of the elements of all the
vectors $\bm{\phi}(\bm{x}_{i})$. We include another hyperparameter
‘gamma_factor’ that scales this default value.

A support vector machine is then trained using scikit-learn’s SVC class, and
prediction is given by

$\displaystyle
y_{\text{pred}}=\text{sign}\;\sum_{i}\alpha_{i}k(\bm{x}_{i},\bm{x}).$ (38)

hyperparameter | values
---|---
_trotter_steps_(embedding layers) $L$ | [1, 3, 5]
_C_ (SVC regularization) | [0.1, 1, 10, 100]
t (time) | [0.01, 0.1, 1.0]
_gamma_factor_ | [0.1,1,10]

###  C.11 QuantumKitchenSinks [55]

This model uses a quantum circuit to define a feature map given by the
concatenation of its output bit-strings. These features are then used to train
a linear classifier. The feature map procedure works as follows:

  * •

Linearly transform the input feature vector $\bm{x}$ as
$\bm{x}^{\prime}_{k}=W_{k}\bm{x}+\bm{b}_{k}$ using randomly sampled $W_{k}$,
$\bm{b}_{k}$ for $k=1,\cdots,k_{\text{max}}$. Here $W_{k}$, $\bm{b}_{k}$ are
such that $\bm{x}^{\prime}$ has dimension $n$ which may be different from $d$.

  * •

Feed each $\bm{x}^{\prime}_{k}$ into a circuit (described below) that returns
a single measurement sample $\bm{z}_{k}\in\\{0,1\\}^{n}$. The concatenated
vector $\bm{z}_{1}\oplus\cdots\oplus\bm{z}_{k_{\text{max}}}$ is the feature
mapped vector of size $n\cdot k_{\text{max}}$.

The circuit used in the second step above consists of angle encoding the
feature mapped vectors with $X$ rotations on individual qubits, and applying
two layers of CNOT gates: the first between adjacent qubits $(j,j+1)$; the
second between qubits $(j,j+2)$ a distance 2 apart. This choice is a natural
generalisation of the example present in [55], which does not present a
specific circuit structure for circuits beyond 4 qubits.

The following is an example of the quantum circuit used to compute feature
vectors QuantumKitchenSinks for the input $\mathbf{x}=(0.1,0.2,0.3,0.4)^{T}$
(here without first applying the classical linear transformation):

The feature mapped vectors are then passed to a linear classifier, which we
implement using scikit-learn’s LogisticRegression class (note logistic
regression can be used for linear classification via the cross entropy loss).

For prediction, a new feature-mapped vector $\bm{z}$ is created via the same
process, and

$\displaystyle y_{\text{pred}}=\text{sign}\;\bm{w}\cdot\bm{z}$ (39)

where $\bm{w}$ is the linear classifier found during training.

The hyperparameter search values are as follows:

hyperparameter | values
---|---
_n_qfeatures_ (circuit size) | [$d$, $\lfloor d/2\rfloor$]
_n_episodes_ (number of circuits $k_{\text{max}}$) | [10, 100, 500, 2000]

###  C.12 QuanvolutionalNeuralNetwork [15]

This model consists of a fixed quantum feature map followed by a trainable
convolutional neural network. The data is first scaled to lie in the range
$[-1,1]$ and a binary threshold function (with threshold zero) is applied and
the data scaled by $\pi$ so that it takes values in $\\{0,\pi\\}$. One then
applies a convolutional layer to the data, where the convolutional filters are
given by a $n_{q}^{2}$ qubit quantum circuits that take as input $n_{q}\times
n_{q}$ sized windows of the input data, in an analogous manner to a classical
convolutional filter. The quantum circuits consists of random gates that we
implement via PennyLane’s RandomLayers class, and the output of the circuits
is the number of ones in the bitstring that has the highest probability to be
sampled from the circuit. As with convolutional neural networks, we allow for
more than one channel in this layer, controlled by the hyperparameter
_n_qchannels_.

The following is an example of the quantum circuit used in the
QuanvolutionalNeuralNetwork for a $2\times 2$-dimensional input window, where
the thresholded and rescaled pixel values form a pre-processed input vector
$(0,0,\pi,\pi)^{T}$:

The resulting data is then fed into a convolutional neural network with the
same specifications as for the ConvolutionalNeuralNetwork.

hyperparameter | values
---|---
_learning_rate_ | [0.0001, 0.001, 0.01]
_n_qchannels_ | [1, 5, 10 ]
_qkernel_shape_ ($n_{q}$) | [2, 3]
_kernel_shape_ (CNN filter size) | [2, 3, 5]

###  C.13 WeiNet [20]

This model (that we call WeiNet following the first author of the paper)
implements a convolutional layer as a unitary operation that acts on input
data that is amplitude encoded into a quantum state. The model has two
registers: the ancilliary register and the work register. The ancilliary
register is used to parameterise a 4 qubit state which in turn controls a
number of unitaries that act on the work register, where the data is encoded
via amplitude encoding. Note that in figure 2 of [20], the Hadamard gates on
the ancilla register have no effect since we trace this register out. The
effect of this register is then to simply perform a classical mixture of the
unitaries $Q_{i}$ defined therein on the work register. For simplicity (and to
save qubit numbers), we parameterise this distribution via $16$ real trainable
parameters.

Two of the qubits are then traced out, which is equivalent to a type of
pooling. All single and double correlators $\langle Z\rangle$ and $\langle
ZZ\rangle$ are measured, and a linear model on these values is used for
classification.

The following is an example of one of the quantum circuits used in the WeiNet
model for a $4\times 4$-dimensional input window encoded into state
$\ket{\psi_{\mathbf{x}}}$ acted on by “filter unitary” $U$:

The paper does not specify the training loss; like the other convolutional
models we use the binary cross entropy with a sigmoid (logistic) activation.
The corresponding probabilities from the sigmoid activation are used for
prediction.

hyperparameter | values
---|---
_learning_rate_ | [0.0001, 0.001, 0.01]
_filter_type_ | [edge_detect, smooth, sharpen ]

###  C.14 MLPClassifier

This is a multiplayer perception model implemented via scikit-learn’s
MLPClassifier class. An input feature vector $\bm{x}$ is transformed using a
sequence of linear transformations $W_{l}$ and element wise-activation
functions $a$ as

$\displaystyle f(\bm{x})=a(W_{L}\cdots a(W_{2}(a(W_{1}\bm{x})))\cdots).$ (40)

‘We use the rectified linear unit activation $a(x)=\max(0,x)$. The trainable
parameters of the model are the weights of the matrices $W_{k}$. We vary the
number of layers $L$ and the shapes of the weight matrices $W_{k}$ via the
model parameter _hidden_layer_size_ , which we set to be one of
$[(100,),(10,10,10,10),(50,10,5)]$. Here, each element of the tuple
corresponds to a different layer and the values give the output dimensions of
the corresponding matrices. The last layer is not included here since it is
always of dimension 1. Training is done via gradient descent with the binary
cross entropy loss using the adam update and the default class fit method. We
vary the initial learning rate and regularisation strength _alpha_. All other
parameters are set to the class defaults, except the maximum number of
iterations, _max_iter_ , that we set to 3000.

hyperparameter | values
---|---
_learning_rate_ | [0.001, 0.01, 0.1]
_hidden_layer_sizes_ | [(100,), (10,10,10,10), (50,10,5) ]
_alpha_ (regularisation) | [0.01, 0.001, 0.0001]

###  C.15 SVC

This model is a support vector machine classifier implemented via scikit-
learn’s SVC class. We use the radial basis function kernel:

$\displaystyle
k(\bm{x},\bm{x}^{\prime})=\exp(-\emph{gamma}||x-x^{\prime}||^{2}).$ (41)

During hyperparameter search, we vary the bandwidth parameter gamma and the
regularisation strength C.

hyperparameter | values
---|---
_C_ (SVC regularization) | [ 0.1, 1, 10, 100]
_gamma_ | [0.001, 0.01, 0.1, 1]

###  C.16 ConvolutionalNeuralNetwork

This model is a vanilla implementation of a convolutional neural network
(CNN), written in flax. The structure of the network is as follows

  * •

a 2D convolutional layer with 32 output channels

  * •

a max pool layer

  * •

a 2D convolutional layer with 64 output channels

  * •

a max pool layer

  * •

a two layer fully connected feedforward neural network with 128 hidden neurons
and one output neuron

The probability of class 1 is given by

$\displaystyle P(+1|\bm{w},\bm{x})=\sigma(f(\bm{w}),\bm{x})$ (42)

where $\bm{w}$ are the weights of the model, $f(\bm{w})$ is the value of the
final neuron, and $\sigma$ is the logistic function. These probabilities are
fed to binary cross entropy loss for training.

hyperparameter | values
---|---
_learning_rate_ | [0.0001, 0.001, 0.01 ]
_kernel_shape_ | [2, 3, 5]

###  C.17 SeparableVariationalClassifier

This is a simple quantum neural network model that does not use entanglement.
The data encoding $V($ consists of $L$ layers, where in each layer arbitrary
trainable single-qubit rotations are performed followed by a product angle
embedding of the data via Pauli Y rotation gates. The encoding is proceeded by
another layer of trainable single-qubit rotations, and prediction is given by
measurement of $O=\frac{1}{n}(Z_{1}+\cdots+Z_{n})$, i.e.

$\displaystyle y_{\text{pred}}=\text{sign}f(\bm{\theta},\bm{x})$ (43)

where

$\displaystyle
f(\bm{\theta},\bm{x})=\langle\frac{1}{n}(Z_{1}+\cdots+Z_{n})\rangle,$ (44)

where $\bm{\theta}$ represents all trainable parameters in the $L$ layers.

The following is an example of the single-qubit quantum circuit used in the
SeparableVariationalClassifier model to process the first value of the input
vector $\mathbf{x}=(0.1,0.2,0.3,0.4)^{T}$ (whereas the remaining features are
processed by similar circuits):

To train the model, we define class probabilities via the logistic function

$\displaystyle P(+1|\bm{\theta},\bm{x})=\sigma(6\langle O\rangle)$ (45)

where we multiply the observable value by $6$ since the sigmoid function
varies significantly over the range $[-6,6]$. These probabilities are then
used in a binary cross entropy loss function.

hyperparameter | values
---|---
_learning_rate_ | [0.001, 0.01, 0.1]
_encoding_layers_ ($L$) | [1,3,5,10]

###  C.18 SeparableKernelClassifier

This model is the kernel equivalent of the above. The data encoding consists
$L$ layers, where in each layer an $X$ rotation with angle $\pi/4$ is applied
to each qubit followed by a product of $Y$ rotations that encode each element
of $\bm{x}$ into individual qubits (so $n=d)$.

The kernel is given by (6):

$\displaystyle
k(\bm{x},\bm{x}^{\prime})=\text{tr}[\rho(\bm{x})\rho(\bm{x}^{\prime})],$ (46)

and the model is trained using scikit-learn’s SVC class.

The following is an example of the single-qubit quantum circuit used in the
SeparableKernelClassifier model to process the first value of the inputs
$\mathbf{x},\mathbf{x^{\prime}}=(0.1,0.2,0.3,0.4)^{T}$, using 2 layers:

hyperparameter | values
---|---
_encoding layers_ ($L$) | [1,3,5,10]
SVM regularization | [0.1, 1, 10, 100]

##  Appendix D Convergence criteria of variational models

Here we describe the criterion used to decide convergence during training of
variational models. This was used for all quantum neural network and quantum
convolutional neural network models, as well as the vanilla convolutional
neural network model. Convergence is decided as follows:

  * •

During training, the previous 400 values of the loss are stored, and at each
step after step 400, the sample mean $\mu_{1}$ of the loss values for the
first 200 and the sample mean $\mu_{2}$ of the loss values of last 200 steps
is calculated, as well as the standard deviation $\sigma_{2}$ of the last 200
loss values.

  * •

If the model has converged, the statistics of the loss values in the two
intervals should be approximately the same, so that the statistics of the
sample means $\mu_{1}$ and $\mu_{2}$ are the same. Under an assumption of
normality, the standard deviation of $\mu_{2}$ is $\sigma_{2}/\sqrt{200}$, and
so for a converged model one expects

$\displaystyle|\mu_{1}-\mu_{2}|<\frac{\sigma_{2}}{\sqrt{200}}$ (47)

with reasonable probability, whereas a model that is still training and whose
loss decreases significantly in 200 steps is unlikely to satisfy the
criterion.

  * •

We adopted a slightly smaller interval than this to decrease the likelyhood of
accidental early stopping. Namely, we used the criterion

$\displaystyle|\mu_{1}-\mu_{2}|<\frac{\sigma_{2}}{2\sqrt{200}}$ (48)

to decide convergence. From visual inspection of loss plots, we found that
this always gave a reasonable stopping criterion that corresponded to an
approximately flat loss curve over the last 200 steps.

##  Appendix E Details on the datasets used in this study

We used 6 types of data generation procedures, five of which are purely
synthetic, while one (MNIST) uses post-processing to reduce the dimension of
the feature vectors. The data generation procedures provide different
variables that can be tuned to create individual datasets. We refer to
datasets for which one such variable is tuned as “benchmarks” and denote them
by capital spelling.

This section describes the procedures in more detail, and gives the settings
for generating the benchmarks. The code can be found in the github repository.

Note that unless otherwise stated, data was split into training and test sets
using a ratio of test/train=$0.2$.

###  E.1 Linearly separable

The linearly separable data generation procedure creates data by a perceptron
with fixed weights that labels input vectors sampled from a hypercube. The
data is hence guaranteed to be linearly separable, and the dataset can be
considered as the easiest or “fruit-fly” task in classification.

Procedure:

  1. 1.

Sample input vectors $\mathbf{x}\in\mathbb{R}^{d}$ uniformly from a
$d$-dimensional hypercube spanning the interval $[-1,1]$ in each dimension.

  2. 2.

Pick a weight vector for some $\mathbf{w}\in\mathbb{R}^{d}$ to define the
linear decision boundary at $\mathbf{w}\mathbf{x}=0$. Only retain the first
$N$ input vectors that do not lie within a margin of size $0.02d$ around the
decision boundary, or $\ \lVert\mathbf{w}\mathbf{x}\rVert>0.02d$.

  3. 3.

Generate continuous-valued labels $y^{\prime}=\mathbf{wx}$.

  4. 4.

Binarise the labels via

$y=\begin{cases}\phantom{-}1\text{ if }y^{\prime}-y_{\text{med}}>0\\\ -1\text{
else},\end{cases}$

where $y_{\text{med}}$ is the median of all continuous labels. This
standarisation procedure ensures that the classes are balanced.

Settings for the linearly separable benchmark:

  * •

Number of features $d\in\\{2,...,20\\}$

  * •

Number of samples $N=300$

  * •

Perceptron weights $\mathbf{w}=(1,...,1)^{T}$

###  E.2 Bars and stripes

The bars and stripes generation procedure is intended to be a simple task for
the three convolutional models. It creates gray-scale images of either
vertical bars or horizontal stripe on a 2D pixel grid.

_Procedure_ :

  1. 1.

Sample $N$ labels $y_{i}=-1,1$ uniformly at random.

  2. 2.

For each $y_{i}$ create a pixel grid of shape $d\times d$ that will store the
data $\bm{x}_{i}$. If $y_{i}=-1$, for each column, sample a random variable
taking values $\pm 1$ with equal probability, and fill the column of
$\bm{x}_{i}$ with this value. If $y_{i}=1$, do the same for the rows.

  3. 3.

For each image $\bm{x}_{i}$ add independent Gaussian noise with standard
deviation $\sigma$ and mean $0$.

Settings for the bars & stripes benchmark:

  * •

Image width $d=4,8,16,32.$

  * •

noise standard deviation $\sigma=0.5$

  * •

Number of data samples $N=1000$

###  E.3 Downscaled MNIST

The MNIST datasets are based on the famous handwritten digits data [43] using
digits $3$ and $5$ which are amongst the hardest to distinguish. The ratio
between test and training set for this procedure is test/train=$0.17$. The
original data was processed by different methods for dimensionality reduction.

The mnist pca, mnist pca- benchmarks use principal component analysis (PCA) to
reduce dimensions.

Procedure:

  1. 1.

Flatten and standarise the inputs images, which is important for PCA to work
well. The standarisation parameters are derived from the training set, and
then used to standarise the test set.222222This best practice takes into
account that in applications one does not necessarily have access to the test
set at training time.

  2. 2.

Compute the $d$ largest principal components of the pre-processed training set
inputs via Principal Component Analysis.

  3. 3.

Project the training and test set inputs onto those components to gain new
input vectors of dimension $d$. This is the mnist pca dataset.

  4. 4.

For the mnist pca- dataset, we sampled a subset of $250$ data points for each
of the training and test set from mnist pca.

Settings for the mnist pca benchmark:

  * •

Number of features $d\in\\{2,...,20\\}$

Settings for the mnist pca- benchmark:

  * •

Number of features $d\in\\{2,...,20\\}$

  * •

Number of samples $N=250$ for training and test set each.

The mnist cg benchmark coarse-grains the pixels of the original images to try
and preserve the correlation structure used by convolutional neural networks.
As before, we use the digits 3 and 5 only.

Procedure:

  1. 1.

Resize the original 28x28 pixel data to a pixel grid of size $H\times H$,
using bilinear interpolation.

  2. 2.

Flatten and standardize the images.

Settings for the mnist cg benchmark:

  * •

Pixel grid height/width $H~{}\in~{}\\{4,8,16,32\\}$.

###  E.4 Hidden manifold model

This data generation procedure is based on Goldt  _et al._ [79], who classify
data sampled from a $m$-dimensional manifold by a neural network, and then
embed the data into a $d$-dimensional space. The properties of this data
generation process allow the authors to compute analytical generalisation
error dynamics using tools from statistical physics. The structure intends to
mimic datasets used in image recognition (such as MNIST), which have been
shown to effectively use low-dimensional manifolds.

Procedure:

  1. 1.

Randomly sample $N$ feature vectors $\mathbf{c}^{m}\in\mathbb{R}^{m}$ with
entries from a standard normal distribution. These vectors lie on the “hidden
manifold”.

  2. 2.

Create an embedding matrix $F\in\mathbf{R}^{d\times m}$.

  3. 3.

Embed the feature vectors via

$\mathbf{x}=\mathbf{\phi}(\mathbf{Fc}/\sqrt{m})$

where $\mathbf{\phi}_{i}(\mathbf{x})=\tanh(x_{i}-b_{i})$

  4. 4.

Generate continuous-valued labels using a neural network applied to the
vectors on the manifold

$y^{\prime}=\mathbf{v}^{T}\mathbf{\varphi}(\mathbf{Wc}/\sqrt{m})$

using component-wise $tanh$ functions as the activation $\varphi$, and the
entries in $W\in\mathbb{R}^{m,m},\mathbf{v}\in\mathbb{R}^{m}$ sampled from a
standard distribution.

  5. 5.

In order to get balanced classes, we rescale the data by subtracting the
median $y_{\text{med}}$ of all labels and then apply a thresholding function

$y=\begin{cases}\phantom{-}1\text{ if }y^{\prime}-y_{\text{med}}>0\\\ -1\text{
else}\end{cases}$

Settings for the hidden manifold benchmark:

  * •

Number of features $d\in\\{2,...,20\\}$

  * •

Number of samples $N=300$

  * •

Manifold dimension $m=6$

  * •

Entries of feature matrix $F$ sampled from a standard distribution

Settings for the hidden manifold diff benchmark:

  * •

Number of features $d=10$

  * •

Number of samples $N=300$

  * •

Manifold dimension $m\in\\{2,...,20\\}$

  * •

Entries of feature matrix $F$ sampled from a standard distribution

###  E.5 Two curves

This data generation procedure is inspired by Buchanan  _et al._ [80], who
consider data sampled from two curves – one for each class – embedded into a
$d$-dimensional space to prove that the maximum curvature and minimum distance
of these curves determine the resources required by a neural network to
generalise well. We can hence understand curvature and distance as two
variables that influence the difficulty of the data.

To control the curvature we use a one-dimensional Fourier series of a maximum
degree $D$ in each dimension. To control the average distance of the curves we
use the same embedding, but shift one curve by some constant.

Procedure:

  1. 1.

Sample $N$ values $t\in\mathbb{R}$ uniformly at random from the interval
$[0,1]$. (This value defines the position of a data point on the curve we
embed.)

  2. 2.

To create the inputs for class $1$, embed half of the $t$-values into a
$d$-dimensional space via a Fourier series defined in every dimension,

$x_{i}=\sum_{n=0}^{D}\alpha^{i}_{n}\cos(nt)+\beta^{i}_{n}\sin(nt)+\epsilon,$

where $D$ is the maximum degree of the Fourier series and
$\\{\alpha^{i}_{n}\\},\\{\beta^{i}_{n}\\}$ are real-valued Fourier
coefficients that we sample uniformly from the interval $[0,1]$. The noise
factor $\epsilon$ determines the variance of a random “spread” added to the
curves around their trajectory in the high-dimensional space.

  3. 3.

To create the inputs for class $-1$, embed the other half of the $t$-values
using the same procedure and Fourier coefficients, but adding an offset of
$\Delta$ to each dimension.

Settings for the two curves benchmark:

  * •

Number of features $d\in\\{2,...,20\\}$

  * •

Number of samples $N=300$

  * •

Noise factor $\epsilon=0.01$

  * •

Maximum degree $D=5$

  * •

Curve offset $\Delta=0.1$

Settings for the two curves diff benchmark:

  * •

Number of features $d=10$

  * •

Number of samples $N=300$

  * •

Noise factor $\epsilon=0.01$

  * •

Maximum degree $D\in\\{2,...,20\\}$

  * •

Curve offset $\Delta=\frac{1}{2D}$

###  E.6 Hyperplanes and parity

We created an artificial dataset that classifies low-dimensional feature
vectors by whether they lie on the “positive” side of an even or odd number of
a set of $k$ hyperplanes. The feature vectors are then embedded into a higher-
dimensional space via a linear transform. The result is a division of the
space into regions of different classes that are delineated by hyperplane
intersections. The parity operation makes sure that a model implicitly has to
learn all hyperplane positions to guess the right label. The difficulty of the
classification problem is expected to increase with the number of hyperplanes.

Procedure:

  1. 1.

Sample $N$ feature vectors $\mathbf{c}\in\mathbb{R}^{m}$ from a standard
normal distribution.

  2. 2.

Embed each feature vector into $\mathbb{R}^{d}$ by multiplying the feature
vectors with a matrix $\mathbf{M}\in\mathbb{R}^{d\times m}$,

$\mathbf{x}=\mathbf{M}\mathbf{c}.$

  3. 3.

Compute $k$ predictions for the $m$-dimensional feature vectors via

$p^{(j)}=\begin{cases}1\text{ if }\mathbf{w}^{(j)}\mathbf{c}+b^{(j)}>0\\\
-1\text{ else}\end{cases}j=1,...,k$

using uniformly sampled weight vectors
$\\{\mathbf{w}^{(j)}\in\mathbb{R}^{m}\\}$ and biases $\\{b^{(j)}\\}$.

  4. 4.

The final label is defined as the parity of these predictions, or whether the
number of 1-predictions is even:

$y=\begin{cases}\phantom{-}1\text{ if
}\sum\limits_{j=1}^{k}\frac{p^{(j)}+1}{2}\text{ even}\\\ -1\text{
else}\end{cases}$

  5. 5.

To ensure balanced classes we initially sample a larger number of datapoints
from which we subsample the desired number for each class.

  6. 6.

Standarise the inputs.

Settings for the hyperplanes diff benchmark:

  * •

Dimension $d=10$

  * •

Number of hyperplanes $k\in\\{2,...,20\\}$

  * •

Number of data samples $N=1000$

  * •

Dimension of hyperplane and initial feature vectors $m=3$

  * •

Entries of $\mathbf{M}$ are uniformly sampled from $[0,1]$

##  Appendix F Collection of detailed results

We add the ranking and accuracy plots of all benchmarks here for readers who
are interested in the detailed results.

Figure 20: Ranking plots like shown in Figure 10 for selected benchmarks.

Figure 21: Ranking plots like shown in Figure 10 for selected benchmarks
(continued).

Figure 22: Detailed training and test accuracies for the benchmarks not shown
in Figure 11.

Figure 23: Results for separable models for the benchmarks not shown in Figure
12. Figure 24: Total ranking results as reported in Figure 10, but with the
two separable models included.