Ultrafast terahertz-field-induced dynamics of superconducting bulk and quasi-1D samples

Within a density-matrix formalism based on the Bardeen–Cooper–Schrieffer (BCS) model and the Bogoliubov–de Gennes equations we provide a description of the dynamics of the non-equilibrium superconducting pairing induced by a terahertz (THz) laser pulse in bulk and quasi-one-dimensional (1D) samples of conventional (BCS-type) superconductors. A cross-over from an adiabatic to a non-adiabatic regime takes place for short and intense THz pulses. In the non-adiabatic regime, the order parameter performs a damped oscillation. We discuss how the parameters of the THz pulse influence the amplitude and the mean value of the oscillation in bulk samples. It is demonstrated that for high intensities the non-adiabatic regime can be reached even for pulses longer than the oscillation period. For the 1D samples we find that the oscillation may attenuate with a different power law. This is analysed by comparing the THz-induced dynamics with the dynamics induced by a sudden switching of the pairing strength, which exhibits essentially the same behaviour. The numerical calculations show that the exponent of the power law depends critically on the density of states in the Debye window and therefore changes in an oscillatory way with the confinement strength. Irregularities in the decay of the oscillation are predicted when the 1D quantum wire is cut short to an elongated zero-dimensional quantum dot structure.

The behaviour of superconducting systems far from equilibrium can provide detailed information on carrier relaxation and pairing interaction dynamics and hence has the potential to give new insights into the origin of superconductivity. This has propelled a wealth of experimental and theoretical works in this area starting from the late 1960s [1][2][3][4][5][6][7][8][9]. In order to prepare a non-equilibrium state, many techniques have been used for a long time like, e.g. tunnel injection and optical irradiation. Nowadays, new state-of-the-art experimental probes have been developed, such as scanning tunneling microscopy, ultrafast lasers, spin-polarized tunneling injection and terahertz (THz) spectroscopy.
Of particular interest are time-resolved studies of the dynamics which provide more powerful means to obtain detailed information on the interactions responsible for electron pairing than static spectral characteristics [10][11][12][13][14][15][16]. In a typical experiment, high-energy quasiparticles are quickly excited out of equilibrium by an intense femtosecond optical laser pulse and then relax by carrier-carrier and carrier-phonon scattering, eventually creating an excess population of low-energy quasiparticles and phonons. The resulting changes of the reflectivity and/or the transmittivity are detected by a probe beam. A major driving force for a renewed interest in this subject is the recent progress in developing new powerful sources for intense ultrashort THz pulses [17][18][19][20][21]. These sources give access to the so far unexplored regime of coherent dynamics of superconducting correlations by generating pulses both shorter than the intrinsic time for superconductor dynamics τ ∼h/ | |, where is the pairing gap, and with central frequencies of the order of the superconducting gap. It thus became possible to study the low energy electromagnetic response of conventional superconductors by monitoring the resonant excitation dynamics of quasiparticles in the vicinity of the superconducting gap without heating the phonon bath. Recently an experimental study of the non-equilibrium Bardeen-Cooper-Schrieffer (BCS) state dynamics by intense THz pulses with pulse duration of 90 fs in a superconducting NbN film with 24 nm thickness has been reported [22] which tests the dynamics in the coherent regime.
Understanding the dynamics of superconductors in the coherent regime also provides a challenge for the theory. The most widely used theoretical approaches for the non-equilibrium dynamics of superconducting systems are based either on the time-dependent Ginzburg-Landau theory (TDGL) [23,24], which is a phenomenological description based on a single variable, the order parameter (t), or on the Boltzmann equation [25,26], which describes the dynamics in terms of a kinetic equation for the quasiparticle distribution coupled to a self-consistent equation for (t). The TDGL procedure is applicable when the quasiparticles are able to reach a local equilibrium compatible with the instantaneous value of (t) on a time scale τ (the quasiparticle energy relaxation time at temperature T , τ ∼h E F / max 2 , k 2 B T 2 with E F being the Fermi energy) much shorter than the typical dynamical time scale of the superconducting order parameter τ ∼h/ | |, i.e. τ τ . This requirement restricts the validity of the TDGL to situations where mechanisms destroying Cooper pairs are effective, for example gapless superconductors. In the opposite limit, τ t 0 τ with t 0 being the duration of the external pulse, the Boltzmann kinetic approach goes beyond the TDGL by keeping all quasiparticle distributions as dynamical variables, but coherences are left out. It usually also involves a gradient expansion for the spatial and momentum dependences of the distribution functions and is, thus, applicable only when the latter variations are slow. Theoretical attempts to model experiments have, on the one hand, used quasi-equilibrium models, like the model in [6,7], which uses an equilibrium distribution function at an effective electron temperature relative to the bath temperature, originally proposed by Kaganov and co-workers [27,28]; and the model proposed by Owen and Scalapino [5], where the system is described in terms of a new chemical potential for the excited quasiparticles. So in the latter case the thermal equilibrium is assumed although chemical equilibrium is not reached for the paired and unpaired electrons. On the other hand, rate equation approaches based on the phenomenological Rothwarf-Taylor (RT) model [2] have been used to describe the recovery dynamics of the superconducting state. In this model, a pair of rate equations describes a system of superconducting quasiparticles coupled to phonons. In the RT phenomenology the excited state is characterized by number densities rather than non-equilibrium energy distribution functions, which drastically simplifies the calculations compared with the approach based on the Boltzmann equations. It is assumed that quasiparticles rapidly thermalize to a narrow range of energy just above the gap and the dynamics are governed by the creation of high-energy phonons due to Cooper-pair recombination and subsequent phonon decay.
In the experiments of [22] the limit t 0 τ τ has been reached, where the quasiparticle spectrum loses its physical meaning and the system evolves non-adiabatically. In this case, the evolution of the system is collisionless in the time interval t < τ and the above-mentioned methods are not applicable. Instead, the dynamics have to be calculated using the BCS model directly. For a fermionic condensate, also described by the BCS model, this leads to the prediction of different regimes of oscillatory behaviour of the pairing potential after a sudden perturbation [29][30][31]. Very similar behaviour has been observed in simulations of a bulk superconductor after short THz excitation [32,33,35].
In the present work we investigate the collisionless dynamics of the BCS paired state induced by an ultrashort THz pulse. Our aim is two-fold: firstly, we examine the conditions necessary to reach the oscillatory regime in bulk superconductors; and secondly, we investigate the effects of spatial confinement present in a nanoscale wire. The nanoscale system is modelled by the Bogoliubov-de Gennes (BdG) equation, which provides a direct extension of the BCS model to spatially non-homogeneous systems [36,37]. Recent technological advances resulted in the fabrication of high-quality superconducting nanoscale systems, such as singlecrystalline Pb and In nanofilms with thicknesses down to a few monolayers [38][39][40][41], Sn and Al nanowires (both single-crystalline and made of strongly coupled grains) with diameters down to 8-10 nm [42][43][44]. Structural imperfections were minimized so that such nanowires and nanofilms did not show signatures of suppression of superconductivity due to disorder [39]. The interplay between superconductivity and quantum confinement has been under significant attention for the last decade due to the discovery of several novel phenomena in nanoscale samples [45][46][47][48][49][50][51] like an increase of the critical temperature in nanowires which originates from quantum-size oscillations of superconducting properties with thickness [52]. It can be expected that the confinement influences the coherent dynamics as well.

Theory
The theory is presented for the general case of spatially inhomogeneous systems. More details about the special case of bulk systems can be found in [32].

BCS model and the Bogoliubov-de Gennes equation
We consider a system described by the usual BCS-Hamiltonian where ψ † α (r) and ψ α (r) are the field operators for an electron with spin α. H e is the singleparticle electron Hamiltonian, which also includes the confinement potential. A point-like form of the electron-electron interaction characterized by the coefficient g is assumed, and an appropriate momentum cut-off confining the interaction to a narrow layer near the Fermi surface (Debye window) is implied.
In the mean field description one obtains the Hartree-Fock-Bogoliubov (HFB) Hamiltonian with the order parameter (r, t) = g ψ ↑ (r) ψ ↓ (r) . This Hamiltonian can be diagonalized with the help of the canonical Bogoliubov transformation, which expresses the electron field operators in terms of new Fermi operators The amplitudes u p (r), v p (r) satisfy at the initial (equilibrium) time the BdG equations As long as the current value of (r) is equal to its initial value (0) (r) (i.e. strictly speaking only in the ground state), the mean-field Hamiltonian is reduced to the diagonal form where the constant E g is the ground state energy of the superconductor and E p can be understood as the excitation energy of one bogolon.
In this work we make use of Anderson's approximate solution to the BdG equation, i.e. we assume that u p (r) = u p ϕ p (r), v p (r) = v p ϕ p (r), where ϕ p (r) are the single-electron wave functions with energy ξ p . In this approximation the equilibrium state of the superconductor is characterized by the relations u 2 . A detailed discussion of the range of validity of the Anderson approximation in equilibrium can be found in [53], and we will here only argue that it can be extended to the nonequilibrium calculations. The approximation implies that the relation p,q (t) p (t) δ p,q needs to be fulfilled approximately, which becomes exact when the order parameter is spatially homogeneous. Therefore, if it is valid in equilibrium, it should also hold in the driven system provided that the external field does not introduce inhomogeneities that significantly modify those already present due to the confinement. We shall assume in the following that the spatial variation of the external field is sufficiently weak in order to justify the approximation p,q (t) = p (t)δ p,q for all times.
Laser driving influences the order parameter and therefore the actual Hamiltonian will in general not be diagonal in the initial bogolon basis. If the order parameter (r, t) is not equal to the initial value (0) (r), the Hamiltonian is given by where we have used the Anderson approximation and omitted the constant term. An index with a bar indicates a time-reversed state, i.e. r| p = p|r . R p and C p depend on the current value of the order parameter and are defined as The order parameter itself is given in the bogolon basis by with the interaction matrix element

Coupling to the external THz laser field
The Hamiltonian H HFB describes a superconductor without external perturbations. In order to model the driving by a THz laser we have to add the coupling to an external electromagnetic field. The latter shall be represented by a vector potential A (r, t), where we adopt the Coulomb gauge ∇A (r, t) = 0. We consider a single Gaussian pulse with full-width at half-maximum (FWHM) of 2 √ ln 2τ such that the corresponding vector potential reads where ω is the central frequency of the pulse and Q = ω/c the wave vector. Here, we have neglected the finite spread of wave vectors that is necessarily present in a pulse of finite width.
In the case of the quantum wire we assume for simplicity that the laser field propagates along the wire, i.e. in z-direction. The amplitude A 0 is perpendicular to the propagation direction. The Hamiltonian for the interaction with the electromagnetic field can be written as Substituting the Bogoliubov transformation equation (3) we can write Altogether, the full Hamiltonian of our model is then given by

Equations of motion
The goal of the present paper is the analysis of the pairing dynamics of a superconductor as manifested in the time evolution of the order parameter (r, t), or, equivalently, p (t).
According to equation (9), its dynamics are determined by the four expectations values . The equations of motion for these dynamical variables are set up via the Heisenberg equation. For the special case of a quantum wire, the quantum number p = ( j, m, k) consists of the radial quantum number j, orbital angular momentum number m and the continuous wave vector k for the free motion along the wire. The numerically necessitated cut-off is performed similarly to the bulk case [32]. Dynamical variables whose pair of indices ( j, m, k) and ( j , m , k ) are too far apart (| j − j | > n j , etc) are neglected. We chose n j = 0, n m = 2 and n k = 4Q, where Q is the photon wave vector. Because only states with j = 0 lie within the Debye window of the thin wire under consideration, setting n j to zero only excludes correlations between states of which at least one lies outside of the Debye window.
We assume that before the arrival of the laser pulse the system is in its ground state, which in our case is the quasiparticle vacuum. This means that all four correlators are zero and the dynamics starts with p = (0) p . With these initial values the subsequent dynamics is then determined by numerically solving our equations of motion.
Most of the results presented in the following sections have been obtained by using material parameters for Pb. Calculations using material parameters for Al [34] as well as for Sn (see section 4) show that the general dynamical behaviour does not qualitatively depend on the choice of material parameters. Pb parameters lead to a considerably smaller number of states in the Debye window than, e.g. Sn parameters, which considerably shortens the very time-consuming numerical calculations in particular in the case of laser excitation. Therefore, despite the known deficiencies of the mean-field theory in quantitatively reproducing the static properties of Pb, we have chosen Pb parameters as a representative example for the laser driven dynamics predicted on the mean-field level. Comparisons in section 4 with a few simulations that we have carried out with Sn parameters confirm that within the mean-field theory the qualitative aspects of the dynamics are insensitive to the choice of parameters.

Adiabatic and non-adiabatic regime
The energy gap of superconducting lead at zero temperature is 2 0 = 2.7 meV. Starting from this equilibrium value, figure 1 shows how laser pulses with two different pulse widths influence the order parameter (t). Both pulses have the same time-integrated intensity and a photon energy ofhω = 2.96 meV slightly above the gap. The pulses can therefore create quasiparticles, which leads to a reduction of the modulus of the order parameter or, equivalently, of the energy gap. While | | is constant after the longer pulse with FWHM 10 ps, a fast oscillation occurs after excitation by the shorter pulse with FWHM 1 ps. The oscillation decays proportional to 1/ √ t as is demonstrated by the dotted lines. This result is general within the BCS formalism and does not depend on the material [34]. In the quasiparticle picture the behaviour of the order parameter can be understood as follows: both pulses excite quasiparticle occupations γ † γ and coherences γ γ , which lowers the value of the order parameter. In the case of the longer pulse, the generation of the quasiparticles occurs adiabatically; when the pulse has subsided, the coherences vanish in the basis that belongs to the instantaneous value of (t), and a stationary state is reached. The shorter pulse, on the other hand, affects the system so fast that it can no longer adiabatically follow and coherences remain after the pulse. The result is a non-stationary state far from equilibrium which exhibits oscillations of the order parameter.
The time dependence of the oscillation closely follows the equation which has been derived in [30] for a BCS system without external driving for a certain set of non-equilibrium initial states. ∞ is the long-time value of | | to which the oscillation eventually decays and a defines the amplitude of the oscillation. The parameters φ and t 0 allow for a shift of the phase of the oscillation and of the origin of the time scale. Obviously the oscillation frequency is connected to the mean value. This leads to the interesting consequence that the oscillation cannot be observed in probe spectra due to the trade-off between temporal resolution and energy resolution. In order to energetically resolve the energy gap, a time interval larger thanh/ ∞ would be required, but during this interval the oscillation averages out [32]. Further simulations suggest that the oscillation should become observable in a twopulse experiment [33].

Influence of pulse intensity and pulse width
We have seen that a sufficiently short THz pulse can put a superconductor into the non-adiabatic regime where the order parameter oscillates. We will now analyse how width and intensity of the pulse influence the amplitude of this oscillation. Figure 2(a) shows the influence of the time-integrated intensity on the oscillation for three different pulse widths. I 0 denotes the same intensity as has been applied in figure 1. The amplitude a (shown in the main part) and the long-time value 2 ∞ (shown in the inset) have been obtained by fitting equation (16) to the numerically calculated data | (t)| for times after the pump pulse. The amplitude of the oscillation rises monotonically with increasing intensity. The increase is almost linear for the shortest pulse, whereas for the longer pulses it starts with a slope close to zero which increases at higher intensities. This is due to the fact that the shortest pulse is shorter than the oscillation period of | (t)|, which is about 1.5 ps. Therefore this pulse always initiates non-adiabatic dynamics. The other pulses are longer than the oscillation period and do not create oscillations in the linear regime. Interestingly, we find, however, that by increasing the intensity we can again excite an oscillation, even though the regime of ultrafast excitation with pulses shorter than the typical time scale of the coherent evolution of the superconductor has not been reached. At these higher intensities we have indeed entered the regime of nonlinear driving, as can also be seen in the intensity dependence of the long-time value shown in the inset. For the longer pulses, the curves start linearly and then flatten at higher intensities. This is due to Pauli blocking and consequently is more pronounced for longer pulses. The flattening occurs in the same intensity range where the oscillations start to appear. For the shortest pulse, the intensity dependence of the long-time value has a quadratic component visible at small intensities. This is because for short pulses the peak intensity is higher and the lowering of the gap does no longer depend linearly on the intensity but gains a quadratic component. The curve for the amplitude in the case of the shortest pulse shows a kink at high intensities. At this point, as can be seen in the inset, | | is substantially decreased to a very low value and equation (16) no longer fits the numerical results very well.
The influence of the pulse width is detailed in figure 2(b) for two values of the timeintegrated intensity I . The main part again shows the amplitude a and the inset the long-time value of the gap, 2 ∞ . For both intensities we see that the amplitude of the oscillation tends to zero for long pulse widths and rises steeply with decreasing widths below a value of about 1.5 ps, which roughly coincides with the oscillation period of the order parameter. This is in agreement with our previous findings that pulses shorter than the oscillation period always lead to non-adiabatic dynamics. While the curve for the lower intensity exhibits a monotonic decrease with increasing pulse width, a flattening in the curve for I = 2I 0 can be seen around 2 ps and the amplitude vanishes only for considerably longer pulses. This means that pulses with a sufficiently high intensity may initiate a non-adiabatic behaviour even if their width is longer than the oscillation period. The mean value plotted in the inset also shows a change of the behaviour around 2 ps. The reason is that for pulses shorter than this value an increasing part of the spectrum falls into the energy gap leading to a reduced generation of quasiparticles and thus to a reduced shrinkage of the gap. For very short pulse widths the mean value ∞ falls again. As already discussed, this is because of the very high peak intensity of the short pulses and the resulting importance of quadratic contributions.

Results: nanowire
As was shown in [54], the superconducting gap in ultrathin slabs and wires exhibits thicknessdependent oscillations accompanied by pronounced resonant enhancements (see section 4.2). The physics behind such oscillations can be outlined as follows. Due to the transverse quantization of the electron motion, the conduction band in nanowires and nanofilms splits up into a series of subbands, and superconductivity is supported by a set of quantum channels. Such single-electron subbands move in energy while changing the nanowire/nanofilm thickness. When the bottom of a subband passes through the Fermi surface, we have an abrupt increase in the density of single-electron states at the Fermi level and, in turn, an enhancement of the superconducting correlations, resulting in an enlargement of the order parameter and other basic superconducting characteristics. This results in quantum-size oscillations of superconducting properties with thickness. It is of high interest to investigate the influence of such resonance effects on the dynamics of the pairing correlations.

Laser induced perturbations
We will first consider a very thin lead cylindrical nanowire with D = 1.18 nm. For the calculations presented here, we have used the following set of parameters which is typical for Pb: g N (0) = 0.39 (N (0) is the bulk energy density of states for one spin projection at the Fermi surface), and Debye temperature T D =hω D /k B = 96 K. Here only subbands with j = 0, m = 0, ±1 have states within the Debye window. The gap equilibrium values are j=0,m=0 = 1.521 meV and 0,±1 = 1.464 meV. The laser energy is set tohω = 3.0 meV, which lies between the subband gaps 2 01 and 2 00 . Because the k-discretization is chosen to match the photon wave vector, the periodic boundary conditions correspond to a wire with length L = 412 µm. Figure 3 shows the dynamics of | 00 (t)| and | 01 (t)| after excitation with a very short THz pulse with FWHM 0.4 ps. The pulse lowers the subband gaps and excites a decaying oscillation very similar to the one we have seen in the bulk case. This is the non-adiabatic regime in a confined superconductor. Just like in the bulk case, long pulses also lower the gap, but after that p is constant (not shown).
For comparison, figure 3 also shows the dynamics induced by an instantaneous change of the coupling constant g from its initial values g i to a value g f . This kind of excitation is computationally much easier to simulate, because the field-dependent terms are responsible for the bulk of the computing time. We obtain a behaviour that is close to the case of field-driven excitation. This will allow us to explore the non-adiabatic regime in confined superconductors in more detail.
For fermionic alkali gases, which are also described by the BCS model, this method of excitation is experimentally accessible [55][56][57]. Of course, an instantaneous change of g is not feasible in an experiment with a conventional superconductor. Our simulations presented in figure 3 nevertheless reveal that it can serve as a reliable approximation to the non-equilibrium state that a THz pulse in the non-adiabatic regime induces in the system. This can be rationalized by recalling that a major effect of the excitation is the generation of quasiparticles which in turn leads to an effective reduction of the coupling g. When the system state right after the pulse is described to a good approximation by calculating the effect of a sudden change of g then also the subsequent dynamics should be the same because after the pulse has vanished the

Quench dynamics
The change in g means that what has formerly been the ground state of the system now becomes an excited state, and because this change happens instantaneously, the non-adiabatic regime with its oscillating order parameter is reached. The influence of the perturbation strength on the dynamics of the spatially averaged order parameter (t) is shown in figure 4. Here we have introduced the spatially averaged order parameter as (t) behaves very similarly to its subband-dependent parts.
One can see that a stronger perturbation, characterized in this particular case by the ratio between the final and the initial value of the coupling strength, g f /g i , results (i) in a smaller stationary value of the order parameter for long times, and (ii) in a reduction of the frequency of the order parameter oscillations. Moreover, after a certain threshold value of the perturbation strength the superconducting correlations are completely destroyed and the dynamics becomes overdamped. This is in accordance with the results of [31]. A more detailed analysis of the influence of the quantum confinement on the threshold value of the pairing strength change will be given elsewhere.
As mentioned before, the exact thickness of the quantum wire has a profound impact on the equilibrium properties; as an example, figure 5 shows the spatial average of the order parameter as a function of the diameter for Pb (left panel) and Sn (right panel) quantum wires. The confinement is also reflected in the dynamical behaviour, as can be seen in figure 6, which   (cf equation (16)) [30]. However, in the resonant cases of D = 6.556 nm (Pb) and D = 5.184 nm (Sn) the oscillation decays much faster, with approximately t −3/4 . The decay exponent therefore strongly depends on the resonance conditions. For wires with a resonant thickness, when the resonant subband gives the largest (main) contribution to the superconducting properties, the decay exponent is found to be close to the value α = 3/4; in the opposite case α = 1/2 is obtained, which coincides with the value found analytically for bulk samples. While making a wire thicker the contribution of the resonant subband becomes less and less important, and the evolution of the pairing condensate exhibits damped oscillations with a power exponent approaching from above its bulk value. In the opposite case of extremely narrow nanowires, when there are only few, but more than one, subbands in the Debye window, and the contribution of the resonant subbands to the superconducting characteristics is dominant, the damping is not a simple power law ∼t −α anymore. Instead the damping is irregular and oscillations can exhibit beating pattern, for example. Nevertheless, even in such cases we can say that the damping of the oscillation occurs faster for resonant wire thicknesses than that for off-resonant thicknesses.

Quantization effects in a finite-length quantum wire
So far we have considered wires with large lengths. We will now show that the energy quantization introduced by a finite length of the quantum wire also changes the coherent nonadiabatic dynamics. Figure 7 illustrates the dynamics of the spatially averaged order parameter after a sudden change of the coupling constant with g f /g i = 0.978 for two different lengths of the quantum wire (in this case, parameters for tin have been used). When comparing these plots with those for long wire lengths one can see several differences and similarities. First the time evolution for short wires can be separated into two regions. The evolution during earlier times is characterized by only one main frequency of oscillation like in the case of long wires. The evolution for large times reveals significant differences to the earlier evolution. After a period of decreasing oscillations the evolution starts to exhibit an irregular behaviour of undamped oscillations with many frequencies. Figure 7 shows that the transition time t c increases with increasing L. Up to approximately 18 ps (the t c of the shorter wire) both curves are almost identical, after that the oscillation in the shorter wire becomes irregular. At about 50 ps similar irregularities occur for the longer wire. Since the length is related to the kdiscretization via k = 2π L , one can expect that the finite discretization is responsible for the effect. In the bulk case the oscillatory behaviour of the pairing correlations is due to the superposition of many continuous band frequencies [30], which cannot come into phase again in finite time. It is understandable that even with such large L the discreteness of the k-values becomes important: the Debye window is very narrow so that even with relatively closely spaced energy states corresponding to rather long wires, there are not many relevant states within the Debye window.
The transition time t c between these two regimes of evolution can be estimated as follows. The evolution obviously consists of a sum of oscillating terms, each term oscillating at a particular frequency 2R p . If two neighbouring terms in the sum are oscillating π out of phase with each other we can expect at least an approximate cancellation of these terms. The transition time should be determined by the energy difference between two neighbouring R p . If two oscillating terms with frequency difference E/h are in phase at the beginning, they are out of phase after the time t c = πh E . The energy difference E = R jmk − R j,m,k− k of two neighbouring terms depends on the subband indices j, m as well as the k-value. E is determined by the wire length L, which governs two neighbouring k-values. Since t c is the earliest time, at which one can observe deviations from the single frequency at the beginning, one must find the indices ( j, m, k), for which E is maximal. The subband with the largest energy difference of terms ξ jmk − ξ j,m,k− k is in the lowest one with j = m = 0. Within this subband, the energies for k and k − k discretized mesh points have the largest difference at the upper border of the Debye window. If the subband gaps jm for all subbands are approximately the same, the same conclusion is true for the R jmk , i.e. E is maximal at j = m = 0 and the largest k for which ξ 0,0,k is still within the Debye window, denoted as k max . In general, the E jmk and time averaged R jmk are different. For the off-resonant thicknesses (if jm (t) and (0) jm are not too different from each other) we can write E jmk − E j,m,k− k ∼ R jmk − R j,m,k− k . Thus This estimate gives t c ≈ 11 ps for the shorter wire and t c ≈ 43 ps for the longer wire, which approximately matches the values observed in the numerical calculations.

Conclusions
In this work we have discussed the non-adiabatic dynamics induced by ultrashort THz pulses in bulk and in nanowire superconductors, calculated within in the mean-field BCS (Bogoliubov-Hartree-Fock) approach. The bulk simulations show that the non-adiabatic regime is generally reached when the THz pulse is shorter than the oscillation period h/ (2 ) corresponding to the order parameter ; but nonlinear processes induced by a higher-intensity pulse can allow the pulse length to be increased to more than twice that value. A very similar non-adiabatic behaviour is found in quantum wires excited by a short THz pulse or by a sudden change of the coupling constant. The main difference to the bulk case lies in the exponent of the power law that governs the decay of the non-adiabatic oscillation: the oscillation decays faster when the size of the quantum wire is at a resonant point, i.e. when a bottom of a subband falls into the Debye window, in the vicinity of the Fermi surface, and thereby supplies a much larger number of electronic states. The power law changes from t −1/2 in the bulk and non-resonant case to t −3/4 in the resonant case. For a finite-length quantum wire, the decay of the oscillation is perturbed after a time determined by the length of the wire. Due to the smallness of the Debye window, this effect occurs for relatively large lengths of several tens of micrometres.