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

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.

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

where

$\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,

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

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],

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.

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.

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.

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.

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.

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,

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.

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.

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

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.

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.

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.

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

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.

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.

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.

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.