On the Evaluation of Higher-Harmonic-Current Responses for High-Field Spectroscopies in Disordered Superconductors

Abstract

We discuss a formalism that allows for the calculation of a higher-harmonic-current response to a strong applied electric field for disordered superconducting systems described on the basis of tight-binding models with on- and/or intersite interactions. The theory is based on an expansion of the density matrix in powers of the field amplitudes, where we solve the equation of motion for the individual components. This allows the evaluation of higher-order response functions on significantly larger lattices than one can achieve with a previously used approach, which is based on a direct temporal integration of the equation of motion for the complete density matrix. In the case of small lattices, where both methods can be applied by including also the contribution of collective modes, we demonstrate the agreement of the corresponding results.

1. Introduction

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, the 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 advanced 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 was advanced in [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 SCs, 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] was 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, was accomplished in [15]. In our previous work [19,26], we 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 [19] and [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 considering the influence of collective modes (amplitude, phase, charge) and, in principle, also allows incorporating 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 [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 was demonstrated in [26].

Here, we compared 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 analyse the two different approaches for the computation of the THG in a disordered tight-binding lattice with attractive on-site interaction (“attractive Hubbard model”) 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.

2. Model

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]):

$$H=\sum _{ij\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }-\left|U\right|\sum _i{n}_{i\uparrow }{n}_{i\downarrow }+\sum _{i\sigma }{V}_i{n}_{i\sigma }$$

(1)

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 the mean-field using the Bogoljubov-de-Gennes (BdG) transformation:

$$c_{i\sigma }=\sum _k\left[u_i\left(k\right){\gamma }_{k,\sigma }-\sigma v_i^{*}\left(k\right){\gamma }_{k,-\sigma }^{†}\right]$$

which yields the eigenvalue equations:

$$\begin{array}{cc}\hfill {\omega }_k{u}_n\left(k\right)& ={\sum }_j{t}_{nj}{u}_j\left(k\right)+[V_n-\frac{\left|U\right|}{2}\langle n_n\rangle -\mu ]u_n\left(k\right)+{\Delta }_n{v}_n\left(k\right)\end{array}$$

(2)

$$\begin{array}{cc}\hfill {\omega }_k{v}_n\left(k\right)& =-{\sum }_j{t}_{nj}^{*}{v}_j\left(k\right)-[V_n-\frac{\left|U\right|}{2}\langle n_n\rangle -\mu ]u_n\left(k\right)+{\Delta }_n^{*}{u}_n\left(k\right)\end{array}$$

(3)

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) and (3), one can iteratively determine the ground state density matrix $\mathcal{R}$ with the elements:

$$\begin{array}{ccc}\hfill {\rho }_{ij}& =& \langle c_{i,\uparrow }^{†}{c}_{j,\uparrow }\rangle =\sum _k\left[v_i\left(k\right)v_j^{*}\left(k\right)(1-f\left(E_k\right))+u_i^{*}\left(k\right)u_j\left(k\right)f\left(E_k\right)\right]\hfill \\ \hfill {\overline{\rho }}_{ij}& =& \langle c_{i,\downarrow }{c}_{j,\downarrow }^{†}\rangle =\sum _k\left[u_i\left(k\right)u_j^{*}\left(k\right)(1-f\left(E_k\right))+v_i^{*}\left(k\right)v_j\left(k\right)f\left(E_k\right)\right]\hfill \\ \hfill {\kappa }_{ij}& =& \langle c_{i,\downarrow }{c}_{j,\uparrow }\rangle =\sum _k\left[-u_i\left(k\right)v_j^{*}\left(k\right)(1-f\left(E_k\right))+v_i^{*}\left(k\right)u_j\left(k\right)f\left(E_k\right)\right]\hfill \end{array}$$

which in compact notation can be written as

$$\mathcal{R}=\left(\begin{array}{cc}\rho & {\kappa }^{†}\\ \kappa & \overline{\rho }\end{array}\right).$$

The BdG approximated energy can then be expressed via the density matrix as

$$E^{BdG}=\sum _{ij}{t}_{ij}\left({\rho }_{ij}-{\overline{\rho }}_{ij}\right)+U\sum _i\left({\rho }_{ii}(1-{\overline{\rho }}_{ii})+{\kappa }_{ii}^{*}{\kappa }_{ii}\right)+\sum _i{V}_i\left[{\rho }_{ii}-{\overline{\rho }}_{ii}+1\right],$$

and the BdG-Hamiltonian matrix is defined as

$${\mathcal{H}}_{ij}^{BdG}=\frac{\partial E^{BdG}}{\partial {\mathcal{R}}_{ji}}.$$

(4)

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 constants 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}$:

$$\begin{array}{ccc}\hfill T^{BdG}& =& -t\left\{e^{i{A}_x}{\rho }_{i+x,i}+e^{-i{A}_x}{\rho }_{i-x,i}-e^{-i{A}_x}{\overline{\rho }}_{i+x,i}-e^{i{A}_x}{\overline{\rho }}_{i+x,i}\right\}\hfill \\ & -& t^{\prime }\left\{e^{i{A}_x}{\rho }_{i+x,i+y}+e^{-i{A}_x}{\rho }_{i-x-y,i}--e^{-i{A}_x}{\overline{\rho }}_{i+x+y,i}-e^{i{A}_x}{\overline{\rho }}_{i+x+y,i}\right.\hfill \\ & +& \left.e^{i{A}_x}{\rho }_{i+x,i-y}+e^{-i{A}_x}{\rho }_{i-x+y,i}-e^{-i{A}_x}{\overline{\rho }}_{i+x-y,i}-e^{i{A}_x}{\overline{\rho }}_{i+x-y,i}\right\}\hfill \end{array}$$

(5)

where we included a nearest ($\sim t$) and next-nearest ($\sim t^{\prime }$) neighbour hopping into the Hamiltonian.

3. Computation of the Dynamics

The equation of motion for the density matrix reads

$$i\frac{d}{dt}\mathcal{R}=\left[\mathcal{R},{\mathcal{H}}^{BdG}\right]$$

(6)

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

$$j_x\left(t\right)=-\frac{1}{N}\frac{\partial E^{BdG}}{\partial A_x}=-\frac{1}{N}\frac{\partial T^{BdG}\left(t\right)}{\partial A_x}$$

(7)

with N denoting the number of sites. The task is now to evaluate the current response 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

$$j_x=\left(1-\frac{1}{2}{A}_x^2\right)j_{para}^x+A_x\left(1-\frac{1}{6}{A}_x^2\right)j_{dia}^x$$

(8)

with

$$\begin{array}{ccc}\hfill j_{para}^x& =& it\sum _n\left[{\rho }_{n+x,n}-{\overline{\rho }}_{n-x,n}-{\rho }_{n-x,n}+{\overline{\rho }}_{n+x,n}\right]\hfill \\ & +& i{t}^{\prime }\sum _n\left[{\rho }_{n+x+y,n}-{\overline{\rho }}_{n-x-y,n}-{\rho }_{n-x-y,n}+{\overline{\rho }}_{n+x+y,n}\right]\hfill \\ & +& i{t}^{\prime }\sum _n\left[{\rho }_{n+x-y,n}-{\overline{\rho }}_{n-x+y,n}-{\rho }_{n-x+y,n}+{\overline{\rho }}_{n+x-y,n}\right]\hfill \\ \hfill j_{dia}^x& =& -t\sum _n\left[{\rho }_{n+x,n}-{\overline{\rho }}_{n-x,n}+{\rho }_{n-x,n}-{\overline{\rho }}_{n+x,n}\right]\hfill \\ & -& t^{\prime }\sum _n\left[{\rho }_{n+x+y,n}-{\overline{\rho }}_{n-x-y,n}+{\rho }_{n-x-y,n}-{\overline{\rho }}_{n+x+y,n}\right]\hfill \\ & -& t^{\prime }\sum _n\left[{\rho }_{n+x-y,n}-{\overline{\rho }}_{n-x+y,n}+{\rho }_{n-x+y,n}-{\overline{\rho }}_{n+x-y,n}\right].\hfill \end{array}$$

Here, the subscripts $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$, which leads to the standard diamagnetic contribution to the current in the 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$:

$$\begin{array}{ccc}\hfill j_{para}^x& =& \sum _n{A}_0^n{j}_{para}^{x,\left(n\right)}\hfill \\ \hfill j_{dia}^x& =& \sum _n{A}_0^n{j}_{dia}^{x,\left(n\right)}\hfill \end{array}$$

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 third harmonic contribution to the current density reads

$$j_x^{\left(3\right)}\left(t\right)=j_{para}^{x,\left(3\right)}\left(t\right)-\frac{1}{2}{f}^2\left(t\right)j_{para}^{x,\left(1\right)}\left(t\right)+f\left(t\right)j_{dia}^{x,\left(2\right)}\left(t\right)-\frac{1}{6}{f}^3\left(t\right)j_{dia}^{x,\left(0\right)}$$

(9)

where we find that the dominant paramagnetic 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 [19], we numerically integrated Equation (6) using an Adams–Bashforth algorithm and an initialising fourth-order Runge–Kutta method. The resulting time-dependent currents $j_{para}^x$ and $j_{dia}^x$ then were separated numerically into the individual components $j_{para}^{x,\left(n\right)}$ and $j_{dia}^{x,\left(n\right)}$ from which, after the Fourier transformation, the frequency-dependent third harmonic response Equation (9) was 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 compared this approach with a different strategy, where from the beginning, we expanded the density matrix in powers of the applied vector potential:

$$\mathcal{R}=\sum _{n=0}{A}_0^n{\mathcal{R}}^{\left(n\right)}.$$

(10)

Here, ${\mathcal{R}}^{\left(0\right)}$ is the equilibrium density matrix for which

$$\left[{\mathcal{R}}^{\left(0\right)},{\mathcal{H}}^{BdG}\right]=0,$$

(11)

as we already emphasised 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 the frequency space. The next subsections will address in detail the evaluation of the current responses up to third order, including the contribution from collective mode via the random phase approximation (RPA).

3.1. First Order

The first-order current contribution, relevant for the evaluation of the optical conductivity, is given by

$$j_x^{\left(1\right)}=j_{para}^{x,\left(1\right)}+A_0{j}_{dia}^{x,\left(0\right)}$$

(12)

which requires the evaluation of the density matrix up to order $n=1$.

By selecting all terms $\sim A_0$ in the equation of motion, Equation (6), one obtains

$$i{\underline{\underline{\stackrel{·}{R}}}}^{\left(1\right)}=\left[{\underline{\underline{R}}}^{\left(1\right)},{\underline{\underline{H}}}^{\left(0\right)}\right]-\left|U\right|\left[{\underline{\underline{R}}}^{\left(0\right)},{\underline{\underline{D}}}^{\left(1\right)}\right]+f\left(t\right)\left[{\underline{\underline{R}}}^{\left(0\right)},\underline{\underline{V}}\right]$$

(13)

with

$$\underline{\underline{V}}=\left(\begin{array}{cc}\underline{\underline{v}}& \underline{\underline{0}}\\ \underline{\underline{0}}& \underline{\underline{v}}\end{array}\right)$$

(14)

and

$$v_{nm}=-it\left[{\delta }_{m,n+x}-{\delta }_{m,n-x}\right]-i{t}^{\prime }\left[{\delta }_{m,n+x+y}-{\delta }_{m,n-x-y}\right]-i{t}^{\prime }\left[{\delta }_{m,n+x-y}-{\delta }_{m,n-x+y}\right].$$

The matrix ${\underline{\underline{D}}}^{\left(1\right)}$ is defined as

$${\underline{\underline{D}}}^{\left(1\right)}=\left(\begin{array}{cc}-{\overline{\rho }}_D^{\left(1\right)}& {\kappa }_D^{†,\left(1\right)}\\ {\kappa }_D^{\left(1\right)}& -{\rho }_D^{\left(1\right)}\end{array}\right)$$

(15)

and the subscript D indicates that it only contains the diagonal elements of the respective matrices, e.g., ${\left[{\rho }_D^{\left(1\right)}\right]}_{nm}\equiv {\left[{\rho }^{\left(1\right)}\right]}_{nm}{\delta }_{nm}$, which are part of ${\underline{\underline{R}}}^{\left(1\right)}$.

The non-perturbed Hamiltonian ${\underline{\underline{H}}}^{\left(0\right)}$ (i.e., for $A_0=0$) can be diagonalised:

$$\begin{array}{ccc}\hfill {\underline{\underline{\tilde{H}}}}^{\left(0\right)}& =& {\underline{\underline{T}}}^{-1}{\underline{\underline{H}}}^{\left(0\right)}\underline{\underline{T}}=\left(\begin{array}{cccccc}-E_N& \cdots & 0& 0& \cdots & 0\\ ⋮& \ddots & 0& ⋮& \cdots & ⋮\\ 0& \cdots & -E_1& 0& \cdots & 0\\ 0& \cdots & 0& E_1& \cdots & 0\\ ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& 0& \cdots & E_N\end{array}\right)\hfill \end{array}$$

and the same transformation also diagonalises the non-perturbed density matrix:

$$\begin{array}{ccc}\hfill {\underline{\underline{\tilde{R}}}}^{\left(0\right)}& =& {\underline{\underline{T}}}^{-1}{\underline{\underline{R}}}^{\left(0\right)}\underline{\underline{T}}=\left(\begin{array}{cccccc}1& \cdots & 0& 0& \cdots & 0\\ ⋮& \ddots & 0& ⋮& \cdots & ⋮\\ 0& \cdots & 1& 0& \cdots & 0\\ 0& \cdots & 0& 0& \cdots & 0\\ ⋮& & ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& 0& \cdots & 0\end{array}\right).\hfill \end{array}$$

With this transformation, Equation (13) can be written as

$$\begin{array}{ccc}\hfill i{\stackrel{·}{\tilde{R}}}_{nm}^{\left(1\right)}& =& (E_{mm}-E_{nn}){\tilde{R}}_{nm}^1-\left|U\right|({\tilde{R}}_{nn}^0-{\tilde{R}}_{mm}^0){\tilde{D}}_{nm}^{\left(1\right)}+({\tilde{R}}_{nn}^0-{\tilde{R}}_{mm}^0){\tilde{V}}_{nm}f\left(t\right)\hfill \end{array}$$

(16)

where $\underline{\underline{\tilde{V}}}$ and ${\underline{\underline{\tilde{D}}}}^{\left(1\right)}$ denote the transformed matrices, Equations (14) and (15).

We now perform a Fourier transformation:

$$\begin{array}{ccc}\hfill {\tilde{R}}_{nm}^1\left(\omega \right)& =& \int dt{e}^{i\omega t}{\tilde{R}}_{nm}^1\left(t\right)\hfill \\ \hfill f\left(\omega \right)& =& \int dt{e}^{i\omega t}f\left(t\right)\hfill \end{array}$$

so that Equation (16) reads

$$\begin{array}{ccc}\hfill {\tilde{R}}_{nm}^1\left(\omega \right)& =& \frac{{\tilde{R}}_{nn}^0-{\tilde{R}}_{mm}^0}{\omega -E_{mm}+E_{nn}}{\tilde{V}}_{nm}f\left(\omega \right)-\left|U\right|\frac{{\tilde{R}}_{nn}^0-{\tilde{R}}_{mm}^0}{\omega -E_{mm}+E_{nn}}{\tilde{D}}_{nm}^{\left(1\right)}\hfill \\ & \equiv & {\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{V}}_{nm}f\left(\omega \right)-\left|U\right|{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{D}}_{nm}^{\left(1\right)}.\hfill \end{array}$$

(17)

On the BCS level ($U=0$), the density matrix is now obtained by transforming back to $R_{ij}^{\left(1\right)}$ in the original site representation. For the practical computation, ${\tilde{\chi }}_{nm}(\omega \to \omega -iϵ)$ should be shifted into the complex plane in order to avoid singularities.

Including fluctuations means including the corrections due to the matrix ${\underline{\underline{D}}}^{\left(1\right)}$. In the original site representation and in the case of local interactions (as in the present case of the attractive Hubbard model), ${\underline{\underline{D}}}^{\left(1\right)}$ has only diagonal elements in $\rho$ and $\kappa$, which in the following, we denote by Greek letters, i.e., $D_{\alpha \beta }^{\left(1\right)}$ refers to a non-zero element of the matrix ${\underline{\underline{D}}}^{\left(1\right)}$. The case of intersite interactions, as, e.g., relevant for the description of d-wave superconductivity, requires a corresponding modification of the following discussion.

However, in the present case, the elements $D_{\alpha \beta }^{\left(1\right)}$ are related to the diagonal elements of the density matrix, which we obtain by back-transforming Equation (17):

$$\begin{array}{ccc}\hfill R_{\alpha \beta }^{\left(1\right)}& =& T_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{V}}_{nm}{T}_{m\beta }^{-1}f\left(\omega \right)-\left|U\right|T_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{D}}_{nm}^{\left(1\right)}{T}_{m\beta }^{-1}\hfill \\ & \equiv & -{\tau }_{\alpha \nu }^y{D}_{\nu \mu }^{\left(1\right),†}{\tau }_{\mu \beta }^y=-{\tau }_{\beta \nu }^y{D}_{\nu \mu }^{\left(1\right)}{\tau }_{\mu \alpha }^y\hfill \end{array}$$

(18)

where we used the following identity for the diagonal elements of the density matrix:

$${\underline{\underline{R}}}_D^{\left(1\right)}=-{\underline{\underline{\tau }}}_y{\underline{\underline{D}}}^{\left(1\right),†}{\underline{\underline{\tau }}}_y$$

(19)

with

$${\underline{\underline{\tau }}}_y=\left(\begin{array}{cc}\underline{\underline{0}}& -i\underline{\underline{1}}\\ i\underline{\underline{1}}& \underline{\underline{0}}.\end{array}\right)$$

(20)

Equation (18) can be solved for $D_{\nu \mu }^{\left(1\right)}$ as

$$\begin{array}{ccc}\hfill D_{\nu \mu }^{\left(1\right)}& =& -{\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{V}}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^{y}f\left(\omega \right)+\left|U\right|{\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{D}}_{nm}^{\left(1\right)}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\hfill \\ & =& -{\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{V}}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^{y}f\left(\omega \right)+\left|U\right|{\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right)T_{n\rho }^{-1}{D}_{\rho \sigma }^{\left(1\right)}{T}_{\sigma m}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\hfill \end{array}$$

where, in the last step, we transformed ${\tilde{D}}^{\left(1\right)}$ back into the original site representation.

We now define

$$\begin{array}{ccc}\hfill K_{\nu \mu }& =& -{\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{V}}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill W_{\nu \mu ,\rho \sigma }& =& {\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right)T_{n\rho }^{-1}{T}_{\sigma m}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\hfill \end{array}$$

(22)

so that the equation for ${\underline{\underline{D}}}^{\left(1\right)}$ is given by

$$D_{\nu \mu }^{\left(1\right)}=K_{\nu \mu }\left(\omega \right)f\left(\omega \right)+\left|U\right|W_{\nu \mu ,\rho \sigma }{D}_{\rho \sigma }^{\left(1\right)}$$

(23)

or

$$\underset{M_{\nu \mu ,\rho \sigma }}{\underbrace{\left[{\delta }_{\nu \mu ,\rho \sigma }-\left|U\right|W_{\nu \mu ,\rho \sigma }\right]}}{D}_{\rho \sigma }^{\left(1\right)}=K_{\nu \mu }\left(\omega \right)f\left(\omega \right)$$

(24)

and therefore,

$${\underline{\underline{D}}}^{\left(1\right)}\left(\omega \right)={\underline{\underline{M}}}^{-1}\left(\omega \right)\underline{\underline{K}}\left(\omega \right)f\left(\omega \right).$$

(25)

Inserting the transformed solution of Equation (25) into Equation (17) yields the transformed solution for the density matrix.

Figure 1 shows the magnitude of the first harmonic response for a particular disorder configuration ($V/t=1$) on a $8\times 8$ square lattice. We compared the current, obtained from the direct time integration of Equation (6), with the result from Equation (17). For the BCS result, we neglected the time evolution of local charge densities and anomalous correlations in the BdG Hamiltonian Equation (4). This amounts to neglecting the contribution of $\tilde{D}$ in Equation (17), which instead is relevant for the inclusion of collective modes within the RPA. In particular, the phase modes are responsible for the excitations (peaky structures in Figure 1b,d) below the optical gap $2\Delta$; cf. [19]. Note that Figure 1 reports the magnitude of the first-order current response so that the finite BCS response below $2\Delta$ is due to the real part of the current–current correlations. Obviously, the energy resolution in the direct time integration (blue dotted lines) depends on the time interval over which the time integration is performed. In the expansion approach, Equation (17), this resolution can be mimicked by using different values for the parameter $ϵ$, which shifts the energy into the complex plane. However, a finite $ϵ$ describes an exponential damping of the time-dependent density matrix over a time scale $\sim 1/ϵ$. On the other hand, there is no damping in the time integration method, but the integration is simply performed over a fixed time interval. In Figure 1, we use 10 (Panels a, b) and 50 (Panels c, d) periods of the applied vector potential as the time interval for the integration. Note that, for each frequency point, a separate time integration is required.

3.2. Second Order

We proceed by evaluating the diamagnetic contribution to the third harmonic current $A_0{j}_{dia}^{x,\left(2\right)}$; cf. Equation (9). Collecting all terms $\sim A_0^2$, we find for the correction to the density matrix in second order:

$$\begin{array}{ccc}\hfill i{\underline{\underline{\stackrel{·}{R}}}}^{\left(2\right)}\left(t\right)& =& \left[{\underline{\underline{R}}}^{\left(2\right)}\left(t\right),{\underline{\underline{H}}}^{\left(0\right)}\right]+\left[{\underline{\underline{R}}}^{\left(1\right)}\left(t\right),\underline{\underline{V}}\right]f\left(t\right)+\frac{1}{2}\left[{\underline{\underline{R}}}^{\left(0\right)},\underline{\underline{C}}\right]f^2\left(t\right)\hfill \\ & -& \left|U\right|\left[{\underline{\underline{R}}}^{\left(1\right)},{\underline{\underline{D}}}^{\left(1\right)}\right]-\left|U\right|\left[{\underline{\underline{R}}}^{\left(0\right)},{\underline{\underline{D}}}^{\left(2\right)}\right]\hfill \end{array}$$

(26)

where we defined the matrix:

$$\underline{\underline{C}}=\left(\begin{array}{cc}\underline{\underline{c}}& \underline{\underline{0}}\\ \underline{\underline{0}}& -\underline{\underline{c}}\end{array}\right)$$

(27)

and

$$\begin{array}{ccc}\hfill c_{nm}& =& t\left[{\delta }_{m,n+x}+{\delta }_{m,n-x}\right]+t^{\prime }\left[{\delta }_{m,n+x+y}+{\delta }_{m,n-x-y}\right]+t^{\prime }\left[{\delta }_{m,n+x-y}+{\delta }_{m,n-x+y}\right].\hfill \end{array}$$

The Fourier transformation yields

$$\begin{array}{ccc}\hfill \omega {\underline{\underline{R}}}^{\left(2\right)}\left(\omega \right)& =& \left[{\underline{\underline{R}}}^{\left(2\right)}\left(\omega \right),{\underline{\underline{H}}}^{\left(0\right)}\right]+\int d\nu \left[{\underline{\underline{R}}}^{\left(1\right)}\left(\nu \right),\underline{\underline{V}}\right]f(\omega -\nu )\hfill \\ & +& \frac{1}{2}\left[{\underline{\underline{R}}}^{\left(0\right)},\underline{\underline{C}}\right]\int d\nu f\left(\nu \right)f(\omega -\nu )-\left|U\right|\int d\nu \left[{\underline{\underline{R}}}^{\left(1\right)}\left(\nu \right),{\underline{\underline{D}}}^{\left(1\right)}(\omega -\nu )\right]\hfill \\ & -& \left|U\right|\left[{\underline{\underline{R}}}^{\left(0\right)},{\underline{\underline{D}}}^{\left(2\right)}\left(\omega \right)\right]\hfill \end{array}$$

(28)

which, upon applying the transformation to diagonal states, can be written as

$$\begin{array}{ccc}\hfill {\tilde{R}}_{nm}^{\left(2\right)}\left(\omega \right)& =& \frac{1}{2}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{C}}_{nm}\int d\nu f\left(\nu \right)f(\omega -\nu )-\left|U\right|{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{D}}_{nm}^{\left(2\right)}\left(\omega \right)\hfill \\ & +& \frac{1}{\omega +E_{nn}-E_{mm}}\int d\nu {\left[{\tilde{r}}^{\left(1\right)}\left(\nu \right),\tilde{\vartheta }(\omega -\nu )\right]}_{nm}f\left(\nu \right)f(\omega -\nu ).\hfill \end{array}$$

(29)

Here, we defined

$$\begin{array}{ccc}\hfill {\tilde{R}}_{nm}^{\left(1\right)}\left(\omega \right)& \equiv & {\tilde{r}}_{nm}^{\left(1\right)}\left(\omega \right)f\left(\omega \right)={\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\vartheta }_{nm}\left(\omega \right)f\left(\omega \right)\hfill \\ \hfill {\vartheta }_{nm}\left(\omega \right)& \equiv & {\tilde{V}}_{nm}-\left|U\right|{\left[{\tilde{d}}^{\left(1\right)}\left(\omega \right)\right]}_{nm}\hfill \\ \hfill {\tilde{D}}_{nm}^{\left(1\right)}\left(\omega \right)& \equiv & {\tilde{d}}_{nm}^{\left(1\right)}\left(\omega \right)f\left(\omega \right).\hfill \end{array}$$

We can now follow the same procedure as in the case of the first-order RPA calculation. By transforming to the real space representation, where $D_{nm}^{\left(2\right)}$ is again diagonal (similar to Equation (15)), one obtains

$$\begin{array}{ccc}& & D_{\rho \sigma }^{\left(2\right)}\left(\omega \right)=M_{\rho \sigma ,\nu \mu }^{-1}\left(\omega \right)G_{\nu \mu }\left(\omega \right)\int d\nu f\left(\nu \right)f(\omega -\nu )\hfill \\ & -& M_{\rho \sigma ,\nu \mu }^{-1}\left(\omega \right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{\omega +E_{nn}-E_{mm}}\int d\nu {\left[{\tilde{r}}^{\left(1\right)}\left(\nu \right),\tilde{\vartheta }(\omega -\nu )\right]}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^{y}f\left(\nu \right)f(\omega -\nu )\hfill \end{array}$$

(30)

where the matrix $\underline{\underline{M}}$ is the same as in Equation (24), and we defined

$$G_{\nu \mu }=-\frac{1}{2}{\tau }_{\mu \alpha }^y{T}_{\alpha n}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{C}}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y.$$

(31)

Then, by solving Equation (30) and plugging the transformed result into Equation (29), one obtains the second-order frequency-dependent contribution to the density matrix in response to an external field $f\left(\omega \right)$.

We exemplify the result for a harmonic external field with $f\left(\omega \right)=\delta (\omega -\mathsf{\Omega })+\delta (\omega +\mathsf{\Omega })$. Then, from Equation (9) it turns out that the diamagnetic contribution to $j_x^{\left(3\right)}\left(t\right)$ is given by $f\left(t\right)j_{dia}^{x,\left(2\right)}\left(t\right)$, which, upon Fourier transformation, implies that $j_x^{\left(3\right)}\left(\omega \right)$ is given by $j_{dia}^{x,\left(2\right)}(\omega -\mathsf{\Omega })$. Thus, the diamagnetic response at $\omega =3\mathsf{\Omega }$ is determined by the density matrix ${\tilde{R}}_{nm}^{\left(2\right)}(\omega -\mathsf{\Omega })$. From Equations (29) and (30), one finds

$$\begin{array}{ccc}\hfill {\tilde{R}}_{nm}^{\left(2\right)}(\omega -\mathsf{\Omega })& =& \frac{1}{2}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(2\mathsf{\Omega }\right){\tilde{C}}_{nm}\delta (\omega -3\mathsf{\Omega })-\left|U\right|{\tilde{\chi }}_{nm}^{\left(0\right)}\left(2\mathsf{\Omega }\right){\tilde{D}}_{nm}^{\left(2\right)}(\omega -\mathsf{\Omega })\hfill \\ & +& \frac{1}{2\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(1\right)}\left(\mathsf{\Omega }\right),\tilde{\vartheta }\left(\mathsf{\Omega }\right)\right]}_{nm}\delta (\omega -3\mathsf{\Omega }),\hfill \end{array}$$

(32)

with

$$\begin{array}{ccc}\hfill D_{\rho \sigma }^{\left(2\right)}(\omega -\mathsf{\Omega })& =& M_{\rho \sigma ,\nu \mu }^{-1}\left(2\mathsf{\Omega }\right)G_{\nu \mu }\left(2\mathsf{\Omega }\right)\delta (\omega -3\mathsf{\Omega })\hfill \\ & -& M_{\rho \sigma ,\nu \mu }^{-1}\left(2\mathsf{\Omega }\right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{2\mathsf{\Omega }+E_{nn}-E_{mm}}\hfill \\ & \times & {\left[{\tilde{r}}^{\left(1\right)}\left(\mathsf{\Omega }\right),\tilde{\vartheta }\left(\mathsf{\Omega }\right)\right]}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\delta (\omega -3\mathsf{\Omega }).\hfill \end{array}$$

(33)

Figure 2 compares the second harmonic response from the direct time integration of Equation (6) with the expansion from Equations (32) and (33) for a particular disorder realisation. As discussed in [19], disorder washes out the resonance at $\omega =\Delta$, and collective modes only slightly increase the intensity of the diamagnetic response. One can also observe that a single parameter $ϵ$ allows adjusting the response, evaluated from the expansion (red line) to the time-integrated result (blue dotted line) at low energy; however, the agreement in intensity is lost at larger values of $\omega$. For larger time integration intervals (cf., Panels c, d), the agreement is pushed to higher energies.

3.3. Third Order

Finally, we evaluated the paramagnetic contribution to the third harmonic current $j_{para}^{x,\left(3\right)}$. Collecting all terms $\sim A_0^3$ in the equation of motion, Equation (6), results in the following equation for the third-order correction to the density matrix

$$\begin{array}{ccc}\hfill i{\underline{\underline{\stackrel{·}{R}}}}^{\left(3\right)}\left(t\right)& =& \left[{\underline{\underline{R}}}^{\left(3\right)}\left(t\right),{\underline{\underline{H}}}^{\left(0\right)}\right]+\left[{\underline{\underline{R}}}^{\left(2\right)}\left(t\right),\underline{\underline{V}}\right]f\left(t\right)+\frac{1}{2}\left[{\underline{\underline{R}}}^{\left(1\right)},\underline{\underline{C}}\right]f^2\left(t\right)\hfill \\ & -& \frac{1}{6}\left[{\underline{\underline{R}}}^{\left(0\right)}\left(t\right),\underline{\underline{V}}\right]f^3\left(t\right)-\left|U\right|\left[{\underline{\underline{R}}}^{\left(0\right)},{\underline{\underline{D}}}^{\left(3\right)}\right]-\left|U\right|\left[{\underline{\underline{R}}}^{\left(1\right)},{\underline{\underline{D}}}^{\left(2\right)}\right]\hfill \\ & -& \left|U\right|\left[{\underline{\underline{R}}}^{\left(2\right)},{\underline{\underline{D}}}^{\left(1\right)}\right].\hfill \end{array}$$

(34)

The Fourier transformation yields

$$\begin{array}{ccc}\hfill \omega \underline{\underline{R^3}}\left(\omega \right)& =& \left[{\underline{\underline{R}}}^{\left(3\right)}\left(\omega \right),{\underline{\underline{H}}}^{\left(0\right)}\right]+\int d{\omega }_1\left[{\underline{\underline{R}}}^{\left(2\right)}\left({\omega }_1\right),\underline{\underline{V}}\right]f(\omega -{\omega }_1)\hfill \\ & +& \frac{1}{2}\int d{\omega }_1\int d{\omega }_2\left[{\underline{\underline{R}}}^{\left(1\right)}\left({\omega }_1\right),\underline{\underline{C}}\right]f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)\hfill \\ & -& \frac{1}{6}\left[{\underline{\underline{R}}}^{\left(0\right)},\underline{\underline{V}}\right]\int d{\omega }_1\int d{\omega }_{2}f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)\hfill \\ & -& \left|U\right|\left[{\underline{\underline{R}}}^{\left(0\right)},{\underline{\underline{D}}}^{\left(3\right)}\left(\omega \right)\right]-\left|U\right|\int d{\omega }_1\left[{\underline{\underline{R}}}^{\left(1\right)}\left({\omega }_1\right),{\underline{\underline{D}}}^{\left(2\right)}(\omega -{\omega }_1)\right]\hfill \\ & -& \left|U\right|\int d{\omega }_1\left[{\underline{\underline{R}}}^{\left(2\right)}\left({\omega }_1\right),{\underline{\underline{D}}}^{\left(1\right)}(\omega -{\omega }_1)\right].\hfill \end{array}$$

(35)

Defining

$$\begin{array}{ccc}\hfill {\underline{\underline{D}}}^{\left(2\right)}\left(\omega \right)& \equiv & \int d\nu {\underline{\underline{d}}}^{\left(2\right)}(\omega ,\nu )f\left(\nu \right)f(\omega -\nu )\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill {\underline{\underline{R}}}^{\left(2\right)}\left(\omega \right)& \equiv & \int d\nu {\underline{\underline{r}}}^{\left(2\right)}(\omega ,\nu )f\left(\nu \right)f(\omega -\nu )\hfill \end{array}$$

(37)

and transforming to diagonal states, Equation (35) becomes

$$\begin{array}{ccc}& & {\tilde{R}}_{nm}^{\left(3\right)}\left(\omega \right)=-\frac{1}{6}{\tilde{\chi }}^{\left(0\right)}\left(\omega \right){\tilde{V}}_{nm}\int d{\omega }_1\int d{\omega }_{2}f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)\hfill \\ & +& \frac{1}{\omega +E_{nn}-E_{mm}}\int d{\omega }_1\int d{\omega }_2{\left[{\tilde{r}}^{\left(2\right)}({\omega }_1+{\omega }_2,{\omega }_2),\underset{=\tilde{\vartheta }(\omega -{\omega }_1-{\omega }_2)}{\underbrace{\tilde{V}-\left|U\right|{\tilde{d}}^{\left(1\right)}(\omega -{\omega }_1-{\omega }_2)}}\right]}_{nm}\hfill \\ & \times & f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)\hfill \\ & +& \frac{1}{\omega +E_{nn}-E_{mm}}\int d{\omega }_1\int d{\omega }_2{\left[{\tilde{r}}^{\left(1\right)}\left({\omega }_1\right),\frac{1}{2}\tilde{C}-\left|U\right|{\tilde{d}}^{\left(2\right)}(\omega -{\omega }_1,{\omega }_2)\right]}_{nm}\hfill \\ & \times & f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)-\left|U\right|{\tilde{\chi }}_{nm}^{\left(0\right)}\left(\omega \right){\tilde{D}}_{nm}^{\left(3\right)}\left(\omega \right).\hfill \end{array}$$

(38)

Now, we follow the usual procedure and write Equation (38) in terms of the diagonal elements, i.e., $D_{\nu \rho }^{\left(3\right)}$, which yields

$$\begin{array}{ccc}& & D_{\rho \sigma }^{\left(3\right)}\left(\omega \right)=-\frac{1}{6}{M}_{\rho \sigma ,\nu \rho }^{-1}{K}_{\nu \rho }\left(\omega \right)\int d{\omega }_1\int d{\omega }_{2}f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)\hfill \\ & -& M_{\rho \sigma ,\nu \rho }^{-1}\left(\omega \right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{\omega +E_{nn}-E_{mm}}\int d{\omega }_1\int d{\omega }_2{\left[{\tilde{r}}^{\left(2\right)}({\omega }_1+{\omega }_2,{\omega }_2),\tilde{\vartheta }(\omega -{\omega }_1-{\omega }_2)\right]}_{nm}\hfill \\ & \times & T_{m\beta }^{-1}{\tau }_{\beta \nu }^{y}f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2)\hfill \\ & -& M_{\rho \sigma ,\nu \rho }^{-1}\left(\omega \right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{\omega +E_{nn}-E_{mm}}\int d{\omega }_1\int d{\omega }_2{\left[{\tilde{r}}^{\left(1\right)}\left({\omega }_1\right),\frac{1}{2}\tilde{C}-\left|U\right|{\tilde{d}}^{\left(2\right)}(\omega -{\omega }_1,{\omega }_2)\right]}_{nm}\hfill \\ & \times & T_{m\beta }^{-1}{\tau }_{\beta \nu }^{y}f\left({\omega }_1\right)f\left({\omega }_2\right)f(\omega -{\omega }_1-{\omega }_2).\hfill \end{array}$$

(39)

which, upon inserting into Equation (38), yields the third-order correction to the density matrix.

We considered again a harmonic external field with $f\left(\omega \right)=\delta (\omega -\mathsf{\Omega })+\delta (\omega +\mathsf{\Omega })$. The contribution of ${\tilde{R}}_{nm}^{\left(3\right)}\left(\omega \right)\sim \delta (\omega -3\mathsf{\Omega })$ is then given by

$$\begin{array}{ccc}\hfill {\tilde{R}}_{nm}^{\left(3\right)}\left(\omega \right)& =& -\frac{1}{6}{\tilde{\chi }}^{\left(0\right)}\left(3\mathsf{\Omega }\right){\tilde{V}}_{nm}\delta (\omega -3\mathsf{\Omega })+\frac{1}{3\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega }),\vartheta \left(\mathsf{\Omega }\right)\right]}_{nm}\delta (\omega -3\mathsf{\Omega })\hfill \\ & +& \frac{1}{3\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(1\right)}\left({\mathsf{\Omega }}_1\right),\frac{1}{2}\tilde{C}-\left|U\right|{\tilde{d}}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega })\right]}_{nm}\delta (\omega -3\mathsf{\Omega })\hfill \\ & -& \left|U\right|{\tilde{\chi }}_{nm}^{\left(0\right)}\left(3\mathsf{\Omega }\right){\tilde{d}}_{nm}^{\left(3\right)}\left(3\mathsf{\Omega }\right)\delta (\omega -3\mathsf{\Omega }).\hfill \end{array}$$

(40)

with

$$\begin{array}{ccc}\hfill {\tilde{r}}_{nm}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega })& =& \frac{1}{2}{\tilde{\chi }}_{nm}^{\left(0\right)}\left(2\mathsf{\Omega }\right){\tilde{C}}_{nm}+\frac{1}{2\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(1\right)}\left(\mathsf{\Omega }\right),\tilde{\vartheta }\left(\mathsf{\Omega }\right)\right]}_{nm}\hfill \\ & -& \left|U\right|{\tilde{\chi }}_{nm}^{\left(0\right)}\left(2\mathsf{\Omega }\right){\tilde{d}}_{nm}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega }).\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill {\tilde{d}}_{nm}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega })& =& M_{\rho \sigma ,\nu \mu }^{-1}\left(2\mathsf{\Omega }\right)G_{\nu \mu }\left(2\mathsf{\Omega }\right)\hfill \\ & -& M_{\rho \sigma ,\nu \mu }^{-1}\left(2\mathsf{\Omega }\right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{2\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(1\right)}\left(\mathsf{\Omega }\right),\tilde{\vartheta }\left(\mathsf{\Omega }\right)\right]}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill {\tilde{d}}_{nm}^{\left(3\right)}\left(3\mathsf{\Omega }\right)& =& -\frac{1}{6}{M}_{\rho \sigma ,\nu \rho }^{-1}\left(3\mathsf{\Omega }\right)K_{\nu \rho }\left(3\mathsf{\Omega }\right)\hfill \\ & -& M_{\rho \sigma ,\nu \rho }^{-1}\left(3\mathsf{\Omega }\right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{3\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega }),\tilde{\vartheta }\left(\mathsf{\Omega }\right)\right]}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y\hfill \\ & -& M_{\rho \sigma ,\nu \rho }^{-1}\left(3\mathsf{\Omega }\right){\tau }_{\mu \alpha }^y{T}_{\alpha ,n}\frac{1}{3\mathsf{\Omega }+E_{nn}-E_{mm}}{\left[{\tilde{r}}^{\left(1\right)}\left(\mathsf{\Omega }\right),\frac{1}{2}\tilde{C}-\left|U\right|{\tilde{d}}^{\left(2\right)}(2\mathsf{\Omega },\mathsf{\Omega })\right]}_{nm}{T}_{m\beta }^{-1}{\tau }_{\beta \nu }^y.\hfill \end{array}$$

(43)

For the same disorder configuration as was used for Figure 1 and Figure 2, we show in Figure 3 the third harmonic response from Equations (17) and (25) as compared to the direct time integration of Equation (6). Consistent with our previous results [19], the strongly disordered ordered sample displays a low paramagnetic energy response at $\omega =\Delta$, both in the BCS and RPA, where the latter is enhanced by the contribution from the collective modes. As in case of the diamagnetic contribution (cf. Figure 2), the “expansion result” (red) for a fixed $ϵ$ parameter can be adjusted to the time-integrated spectrum (blue dotted line) at low energies, but with a decreasing number of periods in the time integration, the agreement in intensity is lost at higher energies. This is particularly visible in Figure 3d, where the contribution from band excitations leads to significantly larger intensities for the small $ϵ=0.004$ as compared to the time integration over 50 periods of the applied vector potential.

4. Conclusions

We presented a detailed comparison of two approaches to evaluate the higher-harmonic-current response to an applied electromagnetic field for disordered and superconducting systems on a lattice. The first method is based on the direct time integration of the equation of motion, Equation (6), as was used in [19] for the investigation of the influence of collective modes in disordered s-wave superconductors. Since, in this case, the higher harmonic contribution has to be extracted numerically from the total response, the calculation has to be performed for at least three different magnitudes of the vector potential for each frequency. Together with the fact that, in order to obtain a reasonable frequency resolution, the integration has to be performed over a significant number of periods of the applied vector potential, this method is limited to a small number of lattice sites. On the other hand, it is rather flexible with regard to the simulation of different pump–probe protocols, which can be easily implemented in the formalism.

Alternatively, one can compute the THG from an expansion of the density matrix in powers of the applied vector potential. As we demonstrated in the present paper, the equations of motion for the individual components can be directly solved in the frequency space from which the currents in the various orders are obtained. In [26], this approach was applied to the evaluation of the third harmonic response in d-wave superconductors, where, at least in the BCS limit, one could treat much larger systems than via the direct time integration of the density matrix. In this paper, we showed how RPA corrections can be included in the formalism. An open issue is the problem of how these RPA corrections can be separated into contributions from the amplitude, phase, and charge modes, which, on the other hand, can be easily accomplished within the time integration method.