Higgs Oscillations in time-resolved Optical Conductivity

TL;DR

This paper calculates the time-resolved optical conductivity of a pumped superconductor using a non-equilibrium approach, revealing Higgs oscillations and order parameter dynamics in the transient state.

Contribution

It introduces a method to compute full vertex corrected optical conductivity in non-equilibrium superconductors, highlighting Higgs mode oscillations.

Findings

Transient suppression of superconductivity observed.

Higgs amplitude mode oscillations detected in the spectrum.

Recovery dynamics of the order parameter characterized.

Abstract

Driving superconductors out of equilibrium is a promising avenue to study their equilibrium properties as well as to control the superconducting state. Non-equilibrium superconductors are often studied using time resolved optical conductivity measurements. Thus, the characterization of a superconducting state in a pump driven non-equilibrium state requires careful attention in the time domain. We calculate time-resolved optical conductivity of a pumped superconducting state using a non-equilibrium Keldysh approach. Through functional derivation, the optical conductivity is obtained with full vertex corrections and used to characterize the transient superconducting state. The transient optical conductivity shows the suppression of the superconducting order parameter in the time domain. The subsequent recovery of the order parameter exhibits oscillatory behavior that corresponds to the Higgs amplitude mode, and may be seen in several parts of the spectrum.

Tables

Table 1. Table 1: Parameters used in the simulation

Phonon frequency ($\Omega$)	$0.2$ $\mathrm{eV}$
Phonon coupling ($g^2$)	$0.12$ $\mathrm{eV}$
Impurity coupling ($\left|V_i\right|$)	0.01 $\mathrm{eV}$
Band parameters	$V_{nn}=0.25\mathrm{eV}$, $\mu =0.0\mathrm{eV}$
Temperature	$\approx 83K$
Pump pulse	${\omega }_p=1.5\mathrm{eV}$ , ${\sigma }_p=8\left[1/\mathrm{eV}\right]$
Probe pulse	$\omega =0.1,\sigma =2\left[1/\mathrm{eV}\right]$

Equations

$$\mathcal{H}$$

$$+\frac{g}{\sqrt{N}}\sum _{\begin{array}{c}\sigma \\ \mathbf{k},\mathbf{q}\end{array}}{c}_{\mathbf{k}+\mathbf{q},\sigma }^{†}{c}_{\mathbf{k},\sigma }^{}\left(b_{\mathbf{q}}^{}+b_{-\mathbf{q}}^{†}\right)$$

$$+\sum _i{V}_i{c}_i^{†}{c}_i$$

$$\mathbf{j}(\mathrm{t})={\mathrm{N}}_{\mathbf{k}}^{-1}\sum _{\mathbf{k}}\nabla \xi (\mathbf{k}-\vec{\mathrm{A}})\mathrm{Im}{\mathrm{G}}_{\mathbf{k}}^{<}(\mathrm{t},\mathrm{t})$$

Full text

I Introduction

Advancements of time-resolved spectroscopic techniques have enabled the measurement Graf et al. (2011) and control Fausti et al. (2011); Hu et al. (2014); Mitrano et al. (2016); Yang et al. (2018a, b, c) of hidden properties of ground states and low energy excitations in correlated materials, which are not easily accessible in equilibrium. In the energy spectrum, relatively low energy scale – terahertz and mid-infrared – frequencies have a special place in such techniques because of their non-invasive nature (to suppress the heat production that may destroy the ordered state) and suitability for observation of low energy excitations of the ground states. This advantage enables the characterization of the ordered ground state and study of the relaxation of the low energy excitations.

Within the context of superconductivity, the excitations of superconducting condensate are of great interest. Experimentally, it is difficult to excite a superconducting order because the order parameter is a scalar field (spinless and chargeless), and does not couple to a vector electromagnetic (EM) field in the linear regime Pekker and Varma (2015). Because a superconducting state exhibits broken $U(1)$ gauge symmetry, when the order is perturbed, $U(1)$ symmetry breaks spontaneously and results in two oscillating bosonic modes: first, the massless phase mode through Goldstone-Nambu mechanism and second, the massive amplitude mode through Higgs-Anderson mechanism. The phase mode gets pushed to the plasma frequency range because of coupling to Coulomb interactions. The other excitation, the Higgs mode resides at a frequency of twice the superconducting gap ($2\Delta$). Recently, it has been shown that the amplitude mode can be excited, and observed, by an EM field using a coexisting order e.g. charge density wave Méasson et al. (2014), by non-linear coupling to EM field (generation of third harmonic) Matsunaga et al. (2014), by time-resolved conductivity Matsunaga et al. (2013) or by using the presence of a supercurrent Nakamura et al. (2018); Moor et al. (2017).

Given that the optical conductivity is the primary probe of the Higgs mode, a proper theoretical description of the conductivity is a necessity. Numerous work has been done to calculate the conductivity (and other response functions) of correlated electrons in an equilibrium Mattis and Bardeen (1958); Zimmermann et al. (1991); Chen (1993); Chubukov and Maslov (2017) and non-equilibrium stateEckstein and Kollar (2008); Eckstein et al. (2010); Orenstein and Dodge (2015); Krull et al. (2014). In a non-equilibrium state, such as induced by a pump-probe setup, a calculation of the conductivity through the Bethe-Salpeter (BS) equation, which is necessary to capture effects beyond the bare-bubble susceptibility, becomes computationally prohibitive because the Hamiltonian loses time-translational symmetry, resulting in 4 separate time variables in the BS equation. Previous solutions to this problem includes exact diagonalization of the Hamiltonian (which is limited by system size) and mean-field analyses of the BCS Hamiltonian (which is a priori not gauge invariant, and neglects inelastic collisions and dynamics of the interactions) Nambu (1960). A recent proposal – dotted tdDMFTTsuji et al. (2016) – has been used to calculate the pair field susceptibility of a superconductor. However, the tdDMFT method is only fully justified for an infinite dimensional case.

In this work, we go beyond these limitations and calculate the time-dependent optical conductivity using a functional derivation approach based on non-equilibrium Green’s functions. We solve the Nambu-Gor’kov equations for electron-phonon mediated superconductivity self-consistently on the Keldysh contour and evaluate the non-equilibrium interacting Green’s functions in the time domain. We consider the Holstein model with impurity scattering as a particular instance to study the transient optical conductivity of a superconductor. The optical conductivity is calculated by a functional derivative of the current with respect to the applied field. One of the advantages of this particular method is that it naturally includes vertex corrections Säkkinen et al. (2012a, b); Stefanucci (2013), but bypasses the calculation of the BS equation in the time domain. In equilibrium, our results reproduce several features of the known conductivity of dirty superconductorsMattis and Bardeen (1958); Zimmermann et al. (1991) such as an upturn towards low frequencies inside the gap. In a pump driven, non-equilibrium case, the conductivity reflects the temporal dynamics of superconducting order including suppression, recovery, and the Higgs oscillations. These are clearly present in the features of the conductivity that are commonly associated with the superconducting order in equilibrium, i.e. the energy location of the gap, the coherence peak, and the phonon features in the real part of the conductivity, as well as in the inductive $1/\omega$ low-frequency response in the imaginary part. We quantify and characterize the transient superconducting state using these features as well a purely time domain feature (the probe current). All the quantities show excellent correlation to the gap dynamics which are known from the underlying simulationsKemper et al. (2015).

The structure of the paper is as follows: first, we describe the Hamiltonian and the method to calculate optical conductivity in equilibrium and non-equilibrium. Then, we discuss the results for the equilibrium state and the dynamics of the Cooper pairs in a pump-driven non-equilibrium state.

II Methods

We use a minimal Hamiltonian to simulate a phonon-mediated, $s$-wave superconductor. We consider the Holstein Hamiltonian on 2D square lattice

[TABLE]

Here, $\xi(\mathbf{k})$ ($=-2\mathrm{V_{nn}}\left[\cos(k_{x})+\cos(k_{y})\right]-\mu$) is the nearest neighbor tight-binding energy dispersion measured relative to the chemical potential, $c^{\dagger}_{\mathbf{k}},c^{\phantom{\dagger}}_{\mathbf{k}}$ ($b^{\dagger}_{\mathbf{q}},b_{\mathbf{q}})$ are the standard creation and annihilation operators for an electron (phonon), $g$ is the momentum-independent $e$-ph coupling constant, and $\Omega$ is the frequency for the Einstein phonon. $V_{i}$ is the coupling between electrons and impurities which are distributed randomly on lattice sites.

The phonon subsystem is treated as a heat reservoir whose properties do not change in time as we drive the electronic subsystem, which is valid for the small pump fluences considered here Abdurazakov et al. (2018); Murakami et al. (2015). The interactions are treated within a self-consistent perturbative framework, which sums the diagrams to all orders. The self energy is calculated using from a Luttinger-Ward functional to ensure conservation-laws Baym and Kadanoff (1961). We use the first order diagram for the electron-phonon interaction, and the self-consistent Born approximation for the impurity scattering. The superconducting state is treated within a self-consistent Migdal-Eliashberg formalism, and the time evolution is done by solving the Gor’kov equations self-consistently on the Keldysh contour Kemper et al. (2015).

Table 1 lists the parameters used in the simulation; these parameters result in an equilibrium superconducting gap $\Delta\approx 46$ meV. The choice of parameters does not represent a specific material. Rather, the parameters were chosen for numerical tractability. Real systems can be simulated by appropriately rescaled parameters.

The pump field, which is applied in the $(11)$ direction, is included to all orders via Peierls’ substitution. In addition to the pump pulse, a secondary (probe) pulse is included in the same way as the pump. However, the probe amplitude and frequency are optimized to ensure that the probe is in the linear-response regime and is able to probe the conductivity within the $2\Delta$ range. The pump and probe pulse envelopes are taken to be Gaussian curves $\mathrm{|\vec{A}(t)|=A_{max}cos(\omega t)\exp(\frac{-t^{2}}{2\sigma})}$.

To calculate transient conductivity we have used the algorithm proposed in Refs. Stefanucci (2013); Shao et al. (2016). The central idea is that first we calculate non-equilibrium current $\mathrm{\vec{J}_{pump}}$ for the pumped state without a probe as follows:

[TABLE]

where the derivative is taken along the field $(11)$ direction. Then, for each pump-probe delay time ($\mathrm{t}$) we calculate change in the current as function of time ($\delta\vec{J}=\mathrm{\vec{J}_{pump+probe}-\vec{J}_{pump}}$). The current and the probe time profiles are then used to calculate time-dependent conductivity as $\mathrm{\sigma(t,\omega)=\frac{\delta\vec{J}(t,\omega)}{\vec{E}_{probe}(\omega)}}$. Here, we have taken the Fourier transform along average-time axis ($t=t^{\prime}$). However, depending on the experimental settings other time axes can also be used to take the Fourier transform, as described in Refs. Hamm (1995); Němec et al. (2002); Averitt and Taylor (2002); Nienhuys and Sundström (2005). For this particular choice, length of the time signal average outs the amplitude of the Higgs mode. We limits the time signal length to $0.01$ [1/eV] to capture the Higgs oscillations within the numerical accuracy.

III Results

III.1 Equilibrium

First, we calculate the conductivity of the system in equilibrium state, i.e. without a pump field. The results are shown in Fig. 1 as a function of temperature $T$. For reference, we have also labeled the curves by their equilibrium superconducting gap $2\Delta$ as determined from the static component of the anomalous retarded self-energy $\Sigma^{F}_{R}(\omega=0)$. In the normal state ($T>T_{c}$) we observe the Drude features in the conductivity near zero frequency, and the effect of the Einstein phonon at the phonon frequency $\Omega$. The presence of a phonon lowers the optical spectral weight in the vicinity of the phonon frequency ($\Omega$); this may be observed as a flattening of the spectral weight in $\sigma_{1}$ at $\Omega$. It is important to note that the minimum of the real part of the conductivity lie at the phonon frequency $\Omega$ in the normal state and shifts by $2\Delta$ in the superconducting state. This particular feature will be used to study the dynamics of the superconducting edge (the gap) in the superconducting state.

In the superconducting state ($T<T_{c}$) we observe the opening of the gap in the conductivity i.e. the lowering of the optical spectral weight inside the $2\Delta$ window (marked by the first shaded region in the figure near zero frequency). In addition, we note the shift of the minimum around the phonon frequency from $\Omega$ to $\Omega+2\Delta$. The normalized conductivity ($\sigma_{1}(T<T_{c})/\sigma_{1}(T=25)$) is plotted in inset (c), which clearly shows the opening of the gap and the development of the coherence peak. The imaginary part of the conductivity (Fig. 1(b)) can also be analyzed in a similar way to characterize the superconducting state; it shows $\frac{1}{\omega}$ behavior inside the $2\Delta$ window. For comparison we include the conductivity calculated from the bare-bubble susceptibility $\chi(\mathbf{q}=0,\omega)=\int d(t^{\prime}-t^{\prime\prime})\sum_{\mathbf{k}}|\vec{v}_{\mathbf{k}}|^{2}[G^{>}_{\mathbf{k}}(t^{\prime},t^{\prime\prime})G^{<}_{\mathbf{k}}(t^{\prime\prime},t^{\prime})-G^{<}_{\mathbf{k}}(t^{\prime},t^{\prime\prime})G^{>}_{\mathbf{k}}(t^{\prime\prime},t^{\prime})]e^{-i\omega(t^{\prime}-t^{\prime\prime})}$ in the figure (dashed line). We observe noticeable qualitative differences because of the vertex corrections in the conductivity calculated from the functional derivation of probe current, mainly near the gap edge and for low energies.

As has been shown in THz pump-probe experiments Matsunaga et al. (2013), the probe current maximum/minimum can also be used as an indicator to study changes in the physical state of a system, e.g. phase transitions or pump induced changes in the system response. We analyze our data in a similar way and plot the equilibrium probe current at different temperatures ranging from normal state to the superconducting state in Fig. 2. As the temperature is reduced, rapid oscillations appear in the probe current, and the minimum is reduced (see Fig. 2(b)). The minimum of the probe current directly correlates with the superconducting order parameter, which is shown in panel (c) where we plot the first minimum of the probe current and the superconducting order parameter. This particular feature can also be used to characterize the transient conductivity in a pump driven superconductor which is shown in the following sections.

III.2 Non-equilibrium

In this section, we discuss the dynamics of electrons in a non-equilibrium state. We calculate the conductivity of a pump-driven non-equilibrium state. Figure 3 shows the real part ($\sigma_{1}$) of the transient conductivity as a function of pump-probe delay time for two pump fluences $\mathrm{A_{max}=0.2,0.4\;[1/a_{0}]}$. We observe noticeable changes in the conductivity from the equilibrium state (at $\mathrm{t=-40\;[1/eV]}$). The suppression of superconducting order can be observed as the edge of the gap (indicated by red markers in region I) moving towards zero during early delay times. As noted above, the location of the minimum in the real part of the conductivity near the phonon frequency $\Omega$ can also be used as an indicator of superconductivity. We observe that the minimum location also shifts on the frequency axis. These minima are located within the shaded region II in the figure. Such suppression of superconductivity is expected in the transient state of the system when the pump drives the system because, intuitively, the pump injects energy in the system and creates excitations. These excitations raise the effective temperature of the system and result in the observed superconducting order suppression. It is important to notice that the effective-temperature picture does not imply a local equilibrium in the transient state as shown previously Kemper et al. (2015); Nosarzewski et al. (2017); Kumar et al. (2019). Rather, the system is in a dynamic non-thermal state where oscillation of the superfluid condensate is observed (this will be discussed in the following sections). For the larger fluence $\mathrm{A_{max}=0.4\;[1/a_{0}]}$, the melting of the superconducting order is stronger and the gap is closed further. Furthermore, it is important to notice that for larger fluences the conductivity ($\sigma_{1}$) does not show a ”clear” spectral gap in the spectrum as shown in Figure 3 panel (b) at delay time $\mathrm{t=24[1/eV]}$. However, the superconductivity remains in the system (as evidenced by a finite off-diagonal order) as shown here Kemper et al. (2015). Out of equilibrium, the direct connection between a gap in the optical spectra and finite off-diagonal order is broken.

III.3 Gap suppression and the Higgs oscillation

It has been proposed and observed Matsunaga et al. (2013, 2014) that when a superconducting condensate is perturbed by an ultrafast pump field, the relaxation of the excited population of the Cooper pairs follows an oscillatory decay mirroring the dynamics predicted for the order parameterVolkov and Kogan (1974); Littlewood and Varma (1982); Krull et al. (2014); Cea et al. (2015); Kemper et al. (2015); Kumar et al. (2019); Murakami et al. (2016). Such amplitude oscillations are called Higgs (or Anderson-Higgs) oscillations. Here they arise from a time dependence of the underlying order parameter, which is reflected in the time-dependent conductivity.

Here we study the dynamics of the superconducting order parameter in the pump-induced non-equilibrium state from the perspective of the time-resolved optical conductivity. For reference, we will use the superconducting gap $\Delta=\Sigma^{F}_{R}(\omega=0)$ as a function of average time. We may estimate the order parameter using the gap edge in $\sigma_{1}(\omega)$, which in equilibrium occurs at $\omega=2\Delta$. We define the edge of the gap as the point $\omega_{edge}$ on the frequency axis where the mean of $(\frac{\sigma_{1}^{sc}}{\sigma_{1}^{ns}})_{max}$ and $(\frac{\sigma_{1}^{sc}}{\sigma_{1}^{ns}})_{min}$ is located within the shaded region I in the figure 3. Similarly, we use the $\sigma_{1}(\omega)$ minimum around the phonon frequency $\Omega$ as a reference. These markers are shown in Fig. 4 panels (a-c), respectively. We observe that the order parameter determined from $\Sigma^{R}_{F}$ is suppressed when the pump is active at early times, and it recovers back to the equilibrium value for later times, exhibiting Higgs oscillations as it recovers. The gap edge and the minimum location exhibit similar behavior, although the relative change is larger at the minimum. All three quantities show Higgs oscillations at approximately the same frequency—in principle the frequency is time dependent as it scales with the local (in time) gapKemper et al. (2015), however in this time range is it approximately constant. Note that there is a small discrepancy in the order parameter value calculated using various pieces of the conductivity data which may arise due to the particular choice of delay axis to Fourier transform, and due to the frequency resolution of the probe signal.

Besides these markers, the Higgs oscillations occur across the response. For example, the probe-current minimum as a function of delay time show the Higgs oscillations as well (c.f. Fig. 4d). Finally, the oscillatory behavior can also be observed when we considering $\sigma_{1}$ at some fixed frequency $\omega_{0}$ as shown in Fig. 4e. Here we have chosen $\omega_{0}$ within $2\Delta$ window, but the oscillation of the conductivity may be seen at all frequencies as a function of time delay (panel (c) Figure 4). A movie of the time-resolved conductivity is available as a supplement.

Panel (f) presents the Fourier transform of the quantities shown in panels (a-e). Although the limited data length leads to wide peaks, the various measurements all oscillate at the same frequency. This is further confirmed by a curve fit to a decaying oscillation (shown on the individual panels), which yields the same frequency for all the measures.

IV Summary

We have presented the time-resolved optical conductivity for a pumped superconductor based on gauge invariant, fully vertex corrected method. The results show that the entire spectrum undergoes changes that reflect the underlying changes in the gap. There are shifts (in energy) of features in the conductivity due to the reduction from $\Delta_{\mathrm{equilibrium}}$ to some reduced $\Delta(t)$, which itself oscillates in time. These “Higgs” oscillations are thus visible in essentially the entire spectrum. We quantify several features that are known to correspond to the gap in equilibrium, e.g. the gap edge and the phonon minimum, and connect them to the underlying gap dynamics which are known from the calculations.

In this work, our analysis is mainly based on the real part of the transient conductivity. For the Higgs mode matter, the imaginary part also exhibits similar oscillations. However, the imaginary part plays an important role in experimental data analysis. Recently Mitrano et al. (2016), the equilibrium-like properties of the imaginary part were utilized to support the claim of the enhanced superconductivity observed in organic superconductor $K_{3}C_{60}$. It is important to understand whether transient optical conductivity exhibits such equilibrium-like properties. This should be analyzed in the future work.

Finally, we stress on the suitability of the method used in this work to calculate transient optical conductivity. The method enables calculation of the response functions beyond the bare-bubble susceptibility. The effect of vertex corrections varies depending on the particulars of the system. For example, they are expected to be minor for an electron-phonon system in the Migdal limit, but not negligible when it comes to impurity scattering in certain regimes (this effect is observed in Fig. 1 for low energies where the impurity scattering is significant). The functional derivative method captures these faithfully and may have broader applicability in the evaluation of equilibrium and non-equilibrium two-particle quantities.

V ACKNOWLEDGMENTS

We acknowledge fruitful discussion with Avinash Rustagi. This work was partly supported by NSF DMR-1752713. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.