Theory of time-resolved optical conductivity of superconductors:comparing two methods for its evaluation

Time-resolved optical conductivity is an oft-used tool to interrogate quantum materials driven out of equilibrium. Theoretically calculating this observable is a complex topic with several approaches discussed in the literature. Using a non-equilibrium Keldysh formalism and a functional derivative approach to the conductivity, we present a comparison of two particular approaches to the calculation of the optical conductivity, and their distinguishing features, as applied to a pumped superconductor. The two methods are distinguished by the relative motion of the probe and gate times; either the probe or gate time is kept fixed while the other is swept. We find that both the methods result in same qualitative features of the time-resolved conductivity after pump is over. However, calculating the conductivity by keeping the gate fixed removes artifacts inherent to the other method. We provide software that, based on data for the first method, is able to construct the second approach.


Introduction
Time-resolved optical conductivity is one of the workhorse experiments for studying quantum materials driven out of equilibrium. Recent advances in THz technology have enabled the time-resolved measurement of the conductivity at low frequencies, and this approach has been applied to a variety of systems, including superconductors driven out of equilibrium, where several novel features have been observed. These include a low-frequency upturn in the inductive response [1], which indicates a potential enhancement of superconductivity and oscillations at a frequency of twice the superconducting gap (2∆) that has been attributed to the Higgs amplitude mode of the superconductor [2-10], although the contribution from light-induced excitation of the Cooper pairs is also shown to be important [11][12][13].
From the theoretical side, the calculation of time-resolved optical conductivity has been limited to few cases, or evaluated [6,11,[14][15][16] using simple models for the electronic states and the time evolution. Notable exceptions are Eckstein et al. [17], Tsuji et al. [18] and Kumar et al. [19] who used a non-equilibrium Green's function approach for the driven electronic states and in one case a numerical functional derivative approach to calculate the optical conductivity. The important advance of the latter is that it includes the vertex corrections due to the included interactions automatically. Kumar et al. [19] studied the time-resolved optical conductivity of a driven superconductor, and observed signatures of the Higgs oscillations across the spectrum.
Fundamentally, the conductivity σ is the linear proportionality between the applied electric field E(t) and the resulting current J(t). In the time domain, this is expressed as arXiv:1909.04669v1 [cond-mat.supr-con] 10 Sep 2019 where we have suppressed vector indices for clarity. That is, we apply a field at some time, and observe the resulting current at some time later. Given this relation, we may take a functional derivative to obtain the conductivity or equivalently, a ratio in the frequency domain Out of equilibrium, the situation becomes more complex. There are now three separate time points: the pump time, the probe time, and the time at which the current is measured (the gate time). These three are illustrated in Fig. 1. The presence of the pump pulse which induces system dynamics breaks the time-translation invariance, which implies the equation for the conductivity now reads and an ambiguity arises for the evaluation in the frequency domain. That is, given that there are now three fixed points in time rather than two, which temporal axis should be Fourier transformed? And, which ones (or which differences) correspond to t and t ? The three time points, t pump , t probe and t gate are shown, as well as the variable time to be Fourier transformed over: τ. The green curve shows a sample probe current (J(t)) obtained with the non-equilibrium Keldysh method. The pump time defines t = 0. (a) Method I sweeps t gate to later times and the delay-time t delay is set by the pump-probe spacing. (b) Method II sweeps the probe pulse att probe to earlier times and t delay is set by the pump-gate spacing. (c) Schematic of the observed currents in the (t gate , t probe ) plane, with the relative time points from (a) and (b) indicated. Methods I and II correspond to horizontal/vertical cuts in this plane, respectively. (d) Probe currents obtained from the non-equilibrium Keldysh formalism after interpolation. Inset: horizontal and vertical cuts of the current at t gate = t pp = 100. Beneath is a plot of the difference (Method II -Method I) between the two curves, scaled by a factor of 10 to increase visibility.
The measurement or calculation of the optical conductivity is typically performed with a pump and a probe at times t pump and t probe , respectively. The simplest approach is to measure the emitted field after the probe as a function of sampling time (by a gate) and take a Fourier transform along this axis, using the pump-probe separation t probe − t pump as the time delay axis in σ(ω, t delay ). This is schematically shown as "Method I" in Fig. 1(a). However, while the signal is collected, the dynamics induced by the pump are still occurring in the system, which are in a sense averaged over the time during which the emitted field is measured. To remedy this, another approach is used where the pump and gate time are kept at a fixed separation, and the probe pulse is swept backwards-this is shown as "Method II" in the figure. This second method has the advantage that the system is always in the same state after the pump whenever the measurement occurs (the probe is assumed to be small and to not affect the dynamics). As was pointed out here [20][21][22], this remedies issues such as the appearance of dynamics before the pump occurs (termed "perturbed free induction"). In this work, we will apply both Methods to obtain the conductivity of a pumped superconductor, and contrast the approaches.

Methods
The conductivities, regardless of which Method, are determined by a functional derivative of the current J(τ). The current is obtained by a non-equilibrium Keldysh Green's function formalism: the self-consistent solution of the Dyson equation on the Keldysh contour [23]. The equations of motion are solved in the superconducting state using a the Nambu Green's functions [24], with strong-electron phonon interactions mediating the pairing interactions. In addition to electron-phonon interactions (which also scatter in addition to providing the pairing glue), we include impurity scattering to consider the dirty limit of BCS and its resulting signal below the energy of the pairing boson [25,26]. The parameters for the calculation are listed in Tab. 1. These parameters were chosen for simplicity of calculation and do not represent any specific material; they may be adjusted to simulate real materials. Samples of the resulting currents are shown in Fig. 1.
We calculated the conductivity using two Methods, as illustrated in Fig. 1, and as explained here. In both cases, a current is measured as a function of a time τ, and a functional derivative is performed numerically by Fourier transforming J(τ) and the electric field, and taking the ratio. The difference arises in which time is kept fixed, and which is swept to evaluate the current J(τ). There are three time points. t pump is the arrival time of the pump pulse, which is used as time zero. t probe is when the probe pulse hits the sample. Finally, t gate is the time when the generated current is measured. The relative time between the pump and the probe is t pp ≡ t pump − t probe .
Method I: we compute for fixed values of t pp the ratio σ(ω, t delay ) = J(ω,t pp ) E(ω,t pp ) , taking the Fourier transform along the t gate axis (in the horizontal direction in Fig. 1(d)). This Method is from a computational perspective straightforward since it simply involves the application of two pulses, and calculating the resulting current.
Method II: t gate is kept at a fixed distance from the pump, and the probe is swept backwards to generate the current J(τ). Then, for fixed values of t gate we compute the ratio σ(ω, t delay ) = J(ω,t gate ) E(ω,t gate ) , taking the Fourier transform along the τ axis [in the vertical direction in Fig. 1(d)]. This Method is computationally more complex since a large number of pump-probe delay sets need to be generated. Here, we have taken data generated from method I and performed Akima spline interpolation [27] to be able to take vertical cuts. The interpolation results are shown in Fig. 1(d).
The advantage of Method II is that the system dynamics, which are driven by the pump, are always in the same state when the current is measured (at t gate ). Since the goal of the measurement is to determine the pump-induced dynamics, this Method may be able to more selectively observe these and provide better time resolution. In Method I, the effective time resolution for the dynamics is set by the decay time of the current (signal-length), which may be long and is generally not known in advance. This long decay time (and thus long effective resolution) produces an averaging over the system dynamics, which may obscure them or in extreme cases hide them if the pump-induced dynamics is comparable to the signal-length.
From an experimental viewpoint, there appears to be a preference towards Method I, potentially due to its simplicity and the long relaxation of the excitations compared to the probe width (see e.g. [1, 3,28]. A notable exception is Ref. [22], which discusses the perturbed free induction decay in some detail. As we will demonstrate below, Method II has potential upsides which may prove useful within the experimental context.

Results
To demonstrate the difference between the two approaches, we consider a driven system that has interesting dynamics after the pump-a pumped superconductor. This system shows non-trivial changes in the conductivity, most notably oscillations of the superconducting order. Oscillations are complex when it comes to evaluating them in the optical conductivity since this requires a Fourier transform; in principle this could simply average over the oscillations and result in a peak in the conductivity rather than any time-dependent behavior, and thus we expect the two methods to show marked differences here. To increase visibility of each curve there is a fixed offset between each conductivity. A video is included as a supplement. Fig. 2 displays the evolution of the real (σ 1 ) and imaginary (σ 2 ) parts of the conductivity computed with both Methods as a function of their respective t delay values. The bottom curves corresponding to t delay = −50.0 indicate the conductivity components in the equilibrium state, before the pump is applied. The conductivity shows the expected features for a strong-coupling BCS superconductor in the presence of impurity scattering: in σ 1 an upturn at low frequencies in the real part, a step near 2∆, and a minimum at Ω + 2∆, and in σ 2 a divergence towards low frequencies. As the pump is applied, the superconductor is partially melted and the features who positions involve ∆ red-shift.
In addition to a reduction in the order parameter, the system exhibits Higgs (or Anderson-Higgs) oscillations, which are an oscillatory decay in the relaxation of the excited population of the Cooper pairs in superconducting condensates subject to perturbation by ultrafast pump fields-these were previously discussed based on similar calculations using the non-equilibrium Keldysh formalism [5,19,24,29]. Higgs oscillations arise here due to a time-dependence of the superconducting order parameter and we observe them in the time-dependent conductivity σ(ω, t delay ) as time-dependent oscillations of the gap edge and minimum around the phonon energy. It is important to note that a critical aspect of the method for observing the Higgs oscillations with Method I is that the probe current decays. If this were not the case, Method I would effectively have no time resolution, and only oscillations in the peak height would be visible [14].
Following the analysis of Kumar et al. [19], we further investigate the time-evolution of the superconducting order parameter in the aftermath of the pump by considering the dynamics of four time-resolved quantities as functions of t delay for the two Methods: 1. The probe current minimum, which was demonstrated to be a measure of the order parameter[2,19]. 2. The location of the gap edge in σ 1 (ω, t delay ), which we define to the the point ω edge on the frequency axis where the mean of σ sc 1 /σ ns 1 max and σ sc 1 /σ ns 1 min is located within the range from ω = 0 to ω ≈ 2∆ equilibrium . 3. The location of the σ 1 (ω, t delay ) minimum about the phonon frequency Ω (measured with respect to Ω). 4. The conductivity at a fixed frequency: σ 1 (ω = 0.083, t delay ).
These four quantities, obtained using both Methods, are shown in Fig. 3, panels (b)-(e), in the order in which they were discussed. The time axis t delay is provided by t pp and t gate for Methods I and II, respectively. For reference, we also show the anomalous "density": F < (t, t). This quantity is an instantaneous measure of the superconductivity in the system; in equilibrium it is equivalent to the right hand side of the gap equation, summed over momenta: where E k = ξ 2 k + ∆ 2 eq , and T is the temperature. Although this is an imperfect measure of the amount of superconductivity in the system because it only captures the amplitude of the order and not the phase [30], it has shown to be correlated with the spectral gap and its dynamics [24]. σ 1 (ω 0 , t delay ) at a fixed frequency around the coherence peak ω 0 ≈ 2∆ eq . There are several differences between the two Methods, which we will now discuss. In all cases, the curves show a suppression of the gap, followed by the characteristic Higgs oscillations, which were discussed in detail in Ref. [19].
However, the most striking difference is an apparent horizontal shift between the curves. For the Fourier transformed quantities obtained from the conductivity (panels (c)-(e)), the dynamics using Method I occur earlier by approximately 50eV −1 . These shifts are due to a mechanism termed "perturbed free induction decay," where for negative times the pump arrives while the probe-induced current is still decaying, causing an earlier than expected observation of changes in the conductivity due to the pump [see e.g. the third current in Fig. 1(c)]. This is in particular relevant for low-frequency features in the conductivity, which by nature require long time signals. This effect produces both the earlier appearance of the maximal change and following shift of the Higgs oscillations, but also the appearance of changes before the pump is active. In contrast, Method II does not have this artifact -the maximal change occurs more closely to where F < (t, t) reaches its maximal change.
In the case of the minimum current, the reverse occurs: the current from Method II reaches its largest change earlier than Method I. Here, this is due to the delay in when the minimum current occurs: it appears some fixed time after the pump, which in Method I shifts the time to higher values. In comparing the currents, Method II more correctly identifies when the change in the gap occurs-it stops decreasing when the pump is off.

Summary
We have evaluated the conductivity of a pumped superconductor with two experimentally accessible methods. The methods involve two arrangements of the three times involved in an optical pump-probe experiment: the pump, probe, and gate times. The conductivities obtained from the two Methods are qualitatively similar, but differ in some key details. The pumped superconductor typically exhibits a suppression, followed by recovery with Higgs oscillations. Both methods discussed here exhibit these features, but with notable differences. First, the (simpler) Method I observes "perturbed free induction" changes before the pump arrives. More generally, due to an effective averaging over the current decay time, Method I's features are somewhat smeared. In contrast, Method II has sharp changes, and exhibits no changes before the pump arrives. Our results suggest that Method II offers improved time resolution over Method I, although Method I does reproduce some of the observed effects.

Materials and Methods
The simulations were performed with the self-consistent non-equilibrium Keldysh formalism described in Ref. [23]. The Holstein Hamiltonian is used to simulate a phonon-mediated, s-wave superconductor on 2D square lattice Here, ξ(k) (= −2V nn cos(k x ) + cos(k y ) − µ) is the nearest neighbor tight-binding energy dispersion with hopping parameter V nn measured relative to the chemical potential µ, c † k , c k (b † q , b q ) are the standard creation and annihilation operators for an electron (phonon), g is the momentum-independent e-ph coupling constant, and Ω is the frequency for the Einstein phonon. V i is the coupling between electrons and impurities which are distributed randomly on lattice sites.
We used the parameters listed in Tab. 1. An oscillating Gaussian pump pulse with a width σ p and a central frequency ω p was applied, followed by a probe pulse of similar shape but with σ and ω as width and central frequency, respectively. As illustrated in Fig. 1, the pump-probe delay time was varied. The generated data was interpolated with Akima splines [27] before taking Fourier transforms in the two directions indicated in the figure.
The software and data used in this manuscript is publicly available at [31].