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

1. Evaluation of the order parameter dynamics in a disordered
superconducting system.

2. Detailed investigation of the interplay between order parameter relaxation and emergence of inhomogeneities.

3. Thorough characterization of the emerging inhomogeneities

1. Discussion of the clean case can be improved, see requested changes
2. Details of the computation should be specified.
3. Some of the figures are too small, see requested changes

In this manuscript authors study the quench dynamics of a
superconducting in real space by solving the time-dependent
Bogoljubov-de Gennes equations. In contrast to the homogeneous case
they find the emergence of inhomogeneities on a longer time scale which
induces an exponential decay of the (averaged) order parameter oscillations.
This is an interesting paper which presents a thorough analysis of the
dynamics on short and longer time scales and provides a detailed discussion
of the emerging inhomogeneities, also as a function of disorder.
I believe that this paper is of significant relevance for the
community and will stimulate further work in this direction.
I therefore recommend publication of the manuscript in SciPost
after the authors have considered the requested changes.

- What induces the emergence of the inhomogeneities in the clean
system and what determined the corresponding time scale?
There should be some "seed" which leads to the development
of the inhomogeneous structures which can be either due to
uncertainties in the preparation of the initial state and (or)
due to the finite time step in the Runge-Kutta integration.
I recommend that authors elaborate a bit more on this point.
In this context authors should also specifiy the accuracy of the
initial iteration and the magnitude of the time step.

- Please enlarge Figs. 1,3,5

- I would guess that most readers (including myself) associate
with "quench" a reduction of the interaction parameter.
While authors specify on page 4 that they consider a temperature
quench (and they also state that quenches in the interaction lead
to quantitatively similar results) it would nevertheless be clearer
if from the beginning (including the abstract) the term "temperature
quench" would be used.

1. An interesting qualitative prediction related with the formation of position dependent fluctuations of the order parameter after a quench.
2. Careful comparison of the initial transient behaviour to other approaches.
3. Careful studying of the effects of disorder

1. The chosen model for the adiabatic time evolution of the electron occupation number with the fast-changing temperature seems unphysical.
2. Many of the quantities in the numerics have been left undefined, making it difficult or impossible for others to repeat the findings (see report).
3. Use of dimensionless units is at some points confusing.
4. Manuscript does not provide information on the possible checks of how discretisation errors or finite-size effects affect the predictions.
5. Manuscript does not explicitly tell how the inhomogeneities are allowed to take place in the numerics in the clean case. In other words, how does the translation symmetry get spontaneously broken?

This manuscript presents an interesting study of what happens in an ordinary two-dimensional superconductor when the temperature is suddenly brought from above the critical temperature to much below the critical temperature. Using a numerical solution of the Bogoliubov-de Gennes equation in a finite lattice, they show how the position averaged amplitude of the order parameter oscillates first in time, but after some initial transient the oscillations decay first in a power law and later exponentially. Then they show how the exponential decay of the oscillations is accompanied by formation of position dependent fluctuations in the order parameter amplitude. They also argue how the oscillations and the power law decay are consistent with earlier treatments of the problems within the position independent BCS approach, but the exponential decay and the formation of position dependent irregularities are naturally beyond them. All this is studied both in the clean and in the disordered limit, showing how the predicted behavior for the clean and weakly disordered cases are quite similar, and in particular the formation of position dependent irregularities does not require disorder.

The results of this study are interesting and I can probably recommend publishing them in SciPost physics after the authors have considered by comments. Namely, the authors do not completely specify what is done in their model (see below). Moreover, the results should be checked against possible issues in discretization. I also suggest them to discuss the results from two points of views, namely fluctuations in two-dimensional systems, and from the models used for the amplitude (Higgs) mode fluctuations in superconductors.

a) Model problems. There are some issues both in the assumptions and in the description of the employed model:
(i) Above Eq. (3), the authors state that the occupation number follows always adiabatically the time-dependent temperature within the fast quench. This seems unphysical to me, because in practice such adiabaticity requires very strong inelastic scattering (scattering rate exceeding the quench time scale, which on the other hand is comparable with the lattice energy scales). It seems to me that within this model the adiabatic quench is possible to reach only with the slow quench discussed in Appendix C.
(ii) The authors set the hopping constant t to unity, implying that all energies should be compared to t and time scales to hbar/t. Then, can the authors clarify if the tau_Q in Eq. (5) has a unit t^2/hbar, or is it perhaps rather related with T_c^2 as implied from the surrounding text.
(iii) Fig. 1 contains the symbol Delta_0 (vertical and horisontal scales), which later in the text is defined as the zero-temperature Delta. It would be better to define it on the first use and also provide its value in terms of t for the chosen parameters (is it 0.83 t as specified below Eq. (6)). Then, can the authors explain why the long-time limit of Delta is so much below Delta_0 (according to Fig. 4 this Delta \approx 0.83 Delta_0, is this a coincidence?). According to the BCS mean field model, Delta(0.1 T_C) is almost the same as Delta_0, the deviation is perhaps one permille.
(iv) The values of U and mu or at least Delta_0 (in units of t) should be provided in Fig. 1 as it depends on them.
(v) I did not find a statement about the boundary conditions employed in the numerics. They should be specified.
(vi) In the clean case, nothing depends on position initially. Therefore, a position independent Delta must be a solution of the dynamics. However, this solution may be unstable to position dependent fluctuations which seems to be the case based on the results. In other words, in order to get a position dependence in the numerics, it has to be put in there either on purpose or accidentally. Can the authors please explain in the text how they did it.
(vii) In Sec. V, the authors claim that the length scale of irregularities is longer than the coherence length, but they do not tell how they define the coherence length, nor its precise value in their model.
(viii) The authors seem to claim that in all of their model, the phase of the order parameter remains essentially fixed. I find it surprising because a 2D system at a finite temperature should exhibit quite large phase fluctuations.

b) Possible issues in discretization
(i) The authors tell that 200x200 lattice is the largest that can be efficiently simulated. However, it would be important to understand the effect of discretization (or finite size effects) on the results. In particular, can the irregularities come from discretization? To check this, perhaps they could show what happens in a slightly smaller lattice (say, 150x150 sites): does the size scale of the irregularities change? If yes, can one argue what might happen in the thermodynamic limit?
(ii) The stripes in Figs. 7 and 8 are either horizontal or vertical. There are two possible reasons for this: they either orient along the lattice directions or along the edge directions, because the two are the same in the simulation. Which one is it? This should be possible to test by changing the relative directions of the edges with respect to the lattice directions (resulting into an irregular edge, but perhaps it does not matter).

c) Physical arguments: could the authors please discuss these
(i) We know from Mermin-Wagner theorem and Berezinskii-Kosterlitz-Thouless physics that in 2D there are no long-range correlations for the phase of the order parameter in superconductors. It is then unclear if there should be a related effect on the amplitude. However, it seems likely that the dimensionality matters for the fluctuations quite much. Can the authors at least emphasize that their results hold for two-dimensional systems, and perhaps speculate what might happen in three dimensions.

(ii) Time-dependent fluctuations of the amplitude of the superconducting order parameter is closely linked with the dynamics of the amplitude or the Higgs mode. It behaves somewhat similarly to what has been discussed in here as it exhibits oscillations and a power-law decay of the order parameter (see, e.g., Eq. (1) in Moor, Volkov, Efetov, PRL 118, 047001 (2017)). Such models typically mostly consider the zero momentum Higgs mode. To my understanding, the results in the present manuscript imply that also the finite-momentum Higgs mode gets excited (as it describes the inhomogeneous state) and even stabilized.

d) Small issues in the paper:
(i) Fig. 4 has the text "increases with disorder". Replace with "increases with increasing disorder"
(ii) On page 9, the authors state "... as a function of time in the clean limit and in the absence of weak disorder". They probably mean "... in the presence of weak disorder".
(iii) Videos on Ref. 41 would require axis labels and some explanation of what the plotted quantities are. Perhaps you could include them in a single site along with a README text. Editors of SciPost might comment on the best way of storing them for long-term use.

1. Define all employed quantities on their first use, and explain how the numerics is done (in particular, how it allows for spontaneous translation symmetry breaking)
2. Check results against discretisation and lattice orientation
3. Discuss results related to the BKT model and Higgs modes