Both, linear and non-linear response spectroscopies provide valuable and complementary information on the excitations of high-temperature superconductors. Since the discovery of these materials optical conductivity measurements have been central in advancing our understanding of the unusual electronic properties, including e.g. the superconducting gap, pseudogap or the transition from a Mott-insulating state to a (non-) Fermi liquid (for a review see e.g. [1]).

On the other hand, the development of non-equilibrium spectroscopies has significantly put forward our understanding of complex materials, due to the possibility of disentangling different dynamical processes via their different relaxation times. [2] With regard to superconducting materials, measurements of the non-linear current response have been recently been applied in order to elucidate the order parameter dynamics which, as a scalar quantity, is not visible in the equilibrium response. Corresponding experiments have been conducted in NbN [3,4,5], MgB₂ [6,7], pnictides [8] and cuprate superconductors [9,10,11,12] whereas the theoretical understanding of these studies has been advanced in Refs. [13,14,15,16,17,18,19,20,21,22,23,24]. Basically, the current density in response to an applied vector potential $\mathbf{A}\left(t\right)$ can be expanded up to third order as $j_{\alpha }={\chi }_{\alpha \beta }^{\left(1\right)}{A}_{\beta }+{\chi }_{\alpha \beta \gamma \delta }^{\left(3\right)}{A}_{\beta }{A}_{\gamma }{A}_{\delta }$ where ${\chi }^{\left(1\right)}$ is the linear response which is related to the optical conductivity. On the other hand, ${\chi }^{\left(3\right)}$ incorporates processes where the (scalar) order parameter fluctuations $\delta \Delta$ are driven by terms quadratic in $\mathbf{A}\left(t\right)$ so that the third harmonic generation (THG) is expected to be enhanced at twice the frequency of the incoming field $2\omega$ corresponding to the spectral gap $2\Delta$ of the superconductor (SC). However, it has been shown [13] that for clean single-band s-wave SC’s these amplitude (’Higgs’) excitations yield only a minor contribution to the THG which instead is dominated by the BCS quasiparticle excitations across the SC spectral gap. For a square lattice, the amplitude excitations only become visible when the THG is measured at an angle of $\pi /4$ with respect to the bond direction which suppresses the QP contribution. For a clean system the response is only due to the diamagnetic current while disorder induces also a paramagnetic contribution [4,15,16,19]. It has been shown [19] that at moderate disorder the response is still dominated by the BCS part while collective modes may yield comparable contributions only in the limit of strong disorder.

In this context various approximations for the theoretical description of disorder within the BCS theory have been considered. In the weakly disordered limit $k_{F}l\gg 1$ previous work [16,21] is either based on the Mattis-Bardeen model [25] or on the self-consistent Born approximation [4]. The summation of diagrams with impurity ladders, equivalent to the solution of quasiclassical Eilenberger equations and formally valid for arbitrary scattering rate, has been accomplished in Ref. [15]. In our previous work [19,26] we have treated the influence of disorder on the THG exactly by solving the time-dependent Bogoljubov-de Gennes equations on finite lattices with local Anderson-type disorder. Formally this has been achieved by adding a time-dependent vector potential to the hamiltonian and by computing the resulting dynamics from the equation of motion for the time-dependent density matrix. This formalism can be accomplished in two different ways which have been used in Ref. [19] and Ref. [26], respectively. First, the dynamics of the full density matrix can be computed from the equation of motion and at the end the various harmonic contributions, proportional to the corresponding power in the amplitude of the applied vector potential $\sim A_0^n$, have to be extracted numerically, see [19]. This is a rather flexible approach, which allows to consider the influence of collective modes (amplitude, phase, charge) and in principle also allows to incorporate different pump-probe protocols. However, for a lattice with N sites, the density matrix for a superconductor has dimensions $2N\times 2N$ so that the formalism is restricted to small lattices only. Second, as outlined in Ref. [26], the density matrix can be expanded from the beginning in powers of the applied vector potential. The equations of motion can be immediately formulated in frequency space for the individual components and allow for the computation of the current response at the respective order. This approach is less flexible in the time domain but can be applied to much larger lattices as relevant for the investigation of d-wave superconductors, at least on the BCS level, as has been demonstrated in Ref. [26].

Here we compare in detail both formalisms and also discuss how collective modes can be incorporated into the second approach. Section 2 introduces the model and we will analyze the two different approaches for the computation of the THG in disordered tight-binding lattices in Section 3. In the same section we also compare the outcome of both procedures for the case of a disordered s-wave system. We conclude our discussion in Section 4.

We illustrate our formalism within the attractive Hubbard model on a square lattice plus local on-site disorder (cf. e.g., [19,27,28,29,30,31]) where the local potential $V_i$ is taken from a flat distribution $-V_0\le V_i\le +V_0$.

To describe the SC state Equation (1) is solved in mean-field using the BdG transformation which yields the eigenvalue equations and the SC order parameter is defined as ${\Delta }_n=-\left|U\right|\langle c_{n,\downarrow }{c}_{n,\uparrow }\rangle$.

From the eigenvalue problem Equations (2,3) one can iteratively determine the ground state density matrix $\mathcal{R}$ with the elements and which in compact notation can be written as

The BdG approximated energy can then be expressed via the density matrix as and the BdG-hamiltonian matrix is defined as

In the absence of an external field the density matrix $\mathcal{R}$ and the Hamiltonian ${\mathcal{H}}^{BdG}$ commute, so that the density matrix has no time evolution. The dynamics of $\mathcal{R}\left(t\right)$ is induced via the coupling to the electromagnetic field $\overrightarrow{E}\left(t\right)=-\partial \overrightarrow{A}\left(t\right)/\partial t$. Let us consider e.g the case of a (spatially constant) field along the x direction. $A_x\left(t\right)$ is coupled to the system via the Peierls substitution $c_{i+x,\sigma }^{†}{c}_{i,\sigma }\to e^{i{A}_x\left(t\right)}{c}_{i+x,\sigma }^{†}{c}_{i,\sigma }$, where for simplicity we will drop from the equations all the constant by putting the lattice spacing, the electronic charge e, the light velocity c and the Planck constant ℏ equal to one. The Peierls substitution, modifies the kinetic-energy part, leading to the following contribution to $E^{BdG}$: where we have included a nearest ($\sim t$) and next-nearest ($\sim t^{\prime }$) neighbor hopping into the hamiltonian.

The equation of motion for the density matrix reads with the BdG hamiltonian matrix given by Equation (4).

Solving Equation (6) yields a time-dependent BdG energy $E^{BdG}\left(t\right)$ and thus a time-dependent current density which is obtained from with N denoting the number of sites. The task is now to evaluate the current reponse for a given order in the amplitude of the applied vector potential. As a first step we expand Equation (5) up to third order in $A_x$ which yields with Here the subscript $para$ and $dia$ refer to the usual identification of the leading terms coupling the gauge field to the fermionic operators in the Hamiltonian, i.e. the linear coupling between the paramagnetic term and $A_x$, and a quadratic coupling between the electronic density and $A_x^2$, that leads to the standard diamagnetic contribution to the current in linear response. However, both $j_{para}^x$ and $j_{dia}^x$ still contain the vector potential to all orders in $A_x$.

Writing $A_x\left(t\right)=A_{0}f\left(t\right)$ we can expand the currents in a power series in $A_0$ which upon inserting into Equation (8) allows us to extract the various current contributions to order n, $j_x^{\left(n\right)}$. In particular, the 3rd harmonic contribution to the current density reads where we find that the dominant para- and diamagnetic contributions are given by $j_{para}^{x,\left(3\right)}$ and $A_0{j}_{dia}^{x,\left(2\right)}$. On the other hand, $j_{para}^{x,\left(1\right)}$ and $j_{dia}^{x,\left(0\right)}$ also enter the calculation of the optical conductivity in first order.

In Ref. [19] we have numerically integrated Equation (6) using an Adams-Bashforth algorithm and an initializing 4th order Runge-Kutta method. The resulting time-dependent currents $j_{para}^x$ and $j_{dia}^x$ then have been separated numerically into the individual components $j_{para}^{x,\left(n\right)}$ and $j_{dia}^{x,\left(n\right)}$ from which, after Fourier transformation, the frequency dependent third harmonic response Equation (9) has been evaluated. In particular at low energies, this procedure is rather time consuming since the integration has to be performed over several periods of the incoming field.

Here we compare this approach with a different strategy where from the beginning we expand the density matrix in powers of the applied vector potential Here ${\mathcal{R}}^{\left(0\right)}$ is the equilibrium density matrix for which as we already emphasized above.

According to Equation (9) higher order contributions to the density matrix ${\mathcal{R}}^{\left(n\right)}$ allow for the computation of the non-harmonic current responses $j_{para}^{x,\left(n\right)}$ and $j_{dia}^{x,\left(n\right)}$ which, as we will show in the following, can be directly obtained in frequency space. The next subsections will address in detail the evaluation of current responses up to 3rd order including the contribution from collective mode via the random phase approximation (RPA).