Topological properties of multi-terminal superconducting nanostructures: effect of a continuous spectrum

TL;DR

This paper explores how the continuous spectrum in multi-terminal superconducting nanostructures influences their topological properties, revealing non-quantized contributions to transconductance and conditions under which these effects vanish.

Contribution

It introduces a redefined Berry curvature accounting for the continuous spectrum and derives explicit expressions for weak energy dependence and Weyl singularities.

Findings

Continuous spectrum adds non-quantized parts to transconductance.

Time-reversible scattering matrices eliminate non-topological contributions.

Derived compact formulas for Berry curvature near Weyl points.

Abstract

Recently, it has been shown that multi-terminal superconducting nanostructures may possess topological properties that involve Berry curvatures in the parametric space of the superconducting phases of the terminals, and associated Chern numbers that are manifested in quantized transconductances of the nanostructure. In this Article, we investigate how the continuous spectrum that is intrinsically present in superconductors, affects these properties. We model the nanostructure within scattering formalism deriving the action and the response function that permits a re-definition of Berry curvature for continuous spectrum. We have found that the re-defined Berry curvature may have a non-topological phase-independent contribution that adds a non-quantized part to the transconductances. This contribution vanishes for a time-reversible scattering matrix. We have found compact expressions for the redefined Berry curvature for the cases of weak energy dependence of the scattering matrix and investigated the vicinity of Weyl singularities in the spectrum.

Figures

Equations

$$\text{det}(e^{2i\chi }-S_{\epsilon }{e}^{i\varphi }{\sigma }_y(S_{-\epsilon }^T{)}^{-1}{\sigma }_y{e}^{-i\varphi })=0,\chi =\arccos (\frac{\epsilon }{\Delta })$$

$$2L=-\mathrm{Tr}\log [{\Pi }_{+}+{\Pi }_{-}{\hat{S}}_{ϵ}],{\Pi }_{\pm }=\frac{1\pm g}{2}$$

$${\hat{S}}_{ϵ}=\left(\begin{array}{cc}{S}_{ϵ}& 0\\ 0& S_{-ϵ}^T\end{array}\right)$$

$$g=U^{†}{\tau }_{z}U,U^{†}=\left(\begin{array}{cc}{e}^{\frac{i\varphi (\tau )}{2}}& 0\\ 0& e^{\frac{-i\varphi (\tau )}{2}}\end{array}\right)\left(\begin{array}{cc}{A}_{-ϵ}& A_{ϵ}\\ A_{ϵ}& A_{-ϵ}\end{array}\right)$$

$$A_{ϵ}=\sqrt{\frac{E+ϵ}{2E}},E=\sqrt{{ϵ}^2+|\Delta {|}^2},$$

$$-2L=2_S\mathrm{Tr}\log [A_{ϵ}{e}^{\frac{-i\varphi (\tau )}{2}}{S}_{ϵ}{e}^{\frac{i\varphi (\tau )}{2}}{A}_{ϵ}+$$

$$+A_{-ϵ}{e}^{\frac{i\varphi (\tau )}{2}}{S}_{-ϵ}^T{e}^{\frac{-i\varphi (\tau )}{2}}{A}_{-ϵ}]$$

$$E_g=-\frac{2_S}{2}\int \frac{dϵ}{2\pi }\mathrm{Tr}\log Q_{ϵ}$$

$$Q_{ϵ}=A_{ϵ}^2{S}_{ϵ}+A_{-ϵ}^2{S}_{-ϵ}^T$$

$$E_g=-\frac{2_S}{2}\int \frac{dϵ}{2\pi }\mathrm{Tr}\log \left(\frac{E+ϵ}{2E}+\frac{E-ϵ}{2E}S{S}^{\ast }\right)+$$

$$+\frac{2_S}{2}\int \frac{dϵ}{2\pi }\log \det (S^T)$$

$$E_g=-\frac{2_S}{2}\sum _{k>0}\int \frac{dϵ}{2\pi }\log [\frac{(E+ϵ{)}^2+(E-ϵ{)}^2+2\cos 2{\chi }_k}{4({ϵ}^2+|\Delta {|}^2)}]$$

$$E_g=-\frac{2_S}{2}\sum _{{ϵ}_k>0}{ϵ}_k$$

$$I_{\alpha }=2e\frac{\partial E_g}{\partial {\varphi }_{\alpha }^{(0)}}.$$

$$\delta L=\sum _{\alpha ,\beta }\int \frac{d\omega }{2\pi }\frac{\delta {\varphi }_{\omega }^{\alpha }\delta {\varphi }_{-\omega }^{\beta }}{2}{R}_{\omega }^{\alpha \beta },$$

$$R_{\omega }^{\alpha \beta }=$$

$$\displaystyle-2_{S}\int\frac{d\epsilon}{2\pi}{\rm Tr}\Big{\{}Q_{\epsilon}^{-1}A_{\epsilon}^{2}[\frac{I_{\alpha}}{2}(S_{\epsilon-\omega}-S_{\epsilon})\frac{I_{\beta}}{2}+$$

$$+\frac{I_{\beta }}{2}(S_{ϵ+\omega }-S_{ϵ})\frac{I_{\alpha }}{2}]+$$

$$+\frac{1}{2}{Q}_{ϵ}^{-1}\frac{{\partial }^2{Q}_{ϵ}}{\partial \alpha \partial \beta }-$$

$$-\frac{1}{2}{Q}_{ϵ+\omega }^{-1}(A_{-(ϵ+\omega )}(\frac{i{I}_{\alpha }}{2}{S}_{-ϵ}^T-S_{-(ϵ+\omega )}^T\frac{i{I}_{\alpha }}{2})A_{-ϵ}-$$

$$-A_{ϵ+\omega }(\frac{i{I}_{\alpha }}{2}{S}_{ϵ}-S_{ϵ+\omega }\frac{i{I}_{\alpha }}{2})A_{\omega })\times$$

$$\times Q_{ϵ}^{-1}(A_{-ϵ}(\frac{i{I}_{\beta }}{2}{S}_{-(ϵ+\omega )}^T-S_{-ϵ}^T\frac{i{I}_{\beta }}{2})A_{-(ϵ+\omega )}-$$

$$\displaystyle-A_{\epsilon}(\frac{iI_{\beta}}{2}S_{\epsilon+\omega}-S_{\epsilon}\frac{iI_{\beta}}{2})A_{\epsilon+\omega})\Big{\}}$$

$$\underset{\omega \to 0}{\lim }{R}_{\omega }^{\alpha \beta }=-\frac{2_S}{2}\frac{{\partial }^2}{\partial \alpha \partial \beta }\int \frac{dϵ}{2\pi }\mathrm{Tr}\log Q_{ϵ}=\frac{{\partial }^2{E}_g}{\partial \alpha \partial \beta }$$

$$R_{\omega }^{\alpha \beta }=\frac{{\partial }^2{E}_g}{\partial \alpha \partial \beta }+\omega B_{\alpha \beta }+O({\omega }^2)$$

$$B_{\alpha \beta }=-\frac{2_S}{2}\int \frac{\mathrm{d}ϵ}{2\pi }(\frac{1}{2}\mathrm{Tr}\left[Q_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial ϵ}{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \alpha }{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \beta }\right]+$$

$$+\frac{\partial }{\partial \beta }\mathrm{Tr}\left[Q_{ϵ}^{-1}{A}^2(ϵ)\{\frac{\partial S_{ϵ}}{\partial ϵ},\frac{i{I}_{\alpha }}{2}\}\right])-(\alpha \leftrightarrow \beta )$$

$$\delta \{\int \frac{dϵ}{2\pi }\mathrm{Tr}\left[Q_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial ϵ}{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \alpha }{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \beta }{e}^{\alpha \beta }\right]\}=$$

$$=\int \frac{dϵ}{2\pi }{\partial }_{ϵ}\mathrm{Tr}\left[Q_{ϵ}^{-1}\delta Q_{ϵ}{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \alpha }{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \beta }\right]e^{\alpha \beta }+$$

$$+\int \frac{dϵ}{2\pi }{\partial }_{\alpha }\mathrm{Tr}[Q_{ϵ}^{-1}\delta Q_{ϵ}{Q}_{ϵ}^{-1}(\frac{\partial Q_{ϵ}}{\partial \beta }{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial ϵ}-$$

$$\frac{\partial Q_{ϵ}}{\partial ϵ}{Q}_{ϵ}^{-1}\frac{\partial Q_{ϵ}}{\partial \beta })]e^{\alpha \beta }$$

$$\frac{1}{2\pi }\delta \{\mathrm{Tr}[S_{-\infty }\frac{I_{\alpha }}{2}{S}_{+\infty }^{†}\frac{I_{\beta }}{2}]\}{e}^{\alpha \beta }$$

$$2\pi q=\mathrm{div}\vec{E},E^{\gamma }\equiv \frac{1}{2}{e}^{\gamma \alpha \beta }{B}_{\alpha \beta }$$

$$2\pi C^{12}={\int }_0^{2\pi }{\int }_0^{2\pi }d{\varphi }_{1}d{\varphi }_2\frac{B_{\alpha \beta }{e}^{\alpha \beta }}{2}=\int (d\vec{s},\vec{E})$$

$${\bar{E}}^{\gamma }=\frac{1}{2\pi }\{\mathrm{Tr}[S_{-\infty }\frac{I_{\alpha }}{2}{S}_{+\infty }^{†}\frac{I_{\beta }}{2}]\}{e}^{\gamma \alpha \beta }$$

$$C=n+2\pi (\overrightarrow{\stackrel{ˉ}{E}},\vec{n})$$

Full text

We have found that the re-defined Berry curvature may have a non-topological phase-independent contribution that adds a non-quantized part to the transconductances. This contribution vanishes for a time-reversible scattering matrix. We have found compact expressions for the redefined Berry curvature for the cases of weak energy dependence of the scattering matrix and investigated the vicinity of Weyl singularities in the spectrum.

I Introduction

The study of topological materials has been on the front edge of the modern research in condensed matter physics for the past decade Xu and Balents (2018); Yao and Wang (2018); Pacholski et al. (2018); Hossain et al. (2018); Tan et al. (2018). These materials are appealing from fundamental point of view and for possible applications (TI-based PhotodetectorBrems et al. (2018); Tang et al. (2017), spintronicsGötte et al. (2014), field-effect transistorMaciejko et al. (2010), catalystChen et al. (2011) and quantum computingNayak et al. (2008); Aasen et al. (2016)). The basis for applications is the topological protection of quantum states, which makes the states robust against small perturbations and leads to many unusual phenomena, e.g. topologically protected edge statesKane and Mele (2005); Wu et al. (2006); Fu et al. (2007). The topological superconductorsQi et al. (2009, 2010); Das et al. (2012); Fu and Berg (2010) and Chern insulatorsHaldane (1988); Regnault and Bernevig (2011); Zhang and Qi (2014); Thonhauser and Vanderbilt (2006) are the classes of topological materials that are relevant for the present paper. In the case of the Chern insulator the topological characteristic is an integer Chern numberQi and Zhang (2011); Moore and Balents (2007) computed with the Green’s function of electrons occupying the bands in a Brillouin zone of a material - WZW formWitten (1983); Wang et al. (2010a, b); Essin and Gurarie (2011). The first Chern number reduces to the sum of first Chern numbers of the filled bands. For each band, the first Chern number is defined as an integral of the Berry curvature over the Brillouin zoneNiu et al. (1985); Thouless et al. (1982). The Berry curvature is commonly definedM. V. Berry (1984) as $B_{\alpha\beta}=2{\rm Im}\langle\partial_{\alpha}k|\partial_{\beta}k\rangle$ with $|k\rangle$ being the wavefunction in this band and $\alpha,\beta$ being the parameters: in this case two components of a wavevector. If the Chern number of a crystal is not zero, the edge states necessarily appear at the interface between the crystal and the vacuum (since the Chern number of the vacuum is zero). The dimensionality of topological materials in real space is restricted by three from above, which significantly limits possible topological phases.

However, there is a way to circumvent this fundamental limitation. Recently, the multi-terminal superconducting nanostructures with conventional superconductors were proposed to realize the topological solids in higher dimensionsRiwar et al. (2016). Such nanostructures host discrete spectrum of so called Andreev bound statesAndreev (1964); de Gennes and Saint-James (1963); Beenakker and van Houten (1991). The energies and wavefunctions of these states depend periodically on the phases of superconducting terminals. This sets an analogy with a bandstructure that depends periodically on the wavevectors. The dimensionality of this bandstructure is the number of terminals minus one. Also, as it was notedRiwar et al. (2016), the multi-terminal superconducting nanostructures cannot be classified as the high-dimensional topological superconductors from the standard periodic table of topological phasesKitaev (2009). The authors of Riwar et al. (2016) have considered in detail 4-terminal superconducting nanostructures and proved the existence of Weyl singularitiesLu et al. (2015); Soluyanov et al. (2015) in the spectrum. The Weyl singularity is manifested as level crossing of Andreev bound states at a certain point in 3-dimensional phase space. Each Weyl singularity can be regarded as a point-like source of Berry curvature. Owing to this, a nonzero two-dimensional Chern number can be realized and is manifested as a quantized transconductance of the nanostructure. This transconductance is the response of the current in one of the terminals on the voltage applied to the other terminal in the limit of small voltage, this signifies an adiabatic regime.

The peculiarity of the system under consideration is the presence of a continuous spectrum next to the discrete one. These states are the extended states in the terminals with energies above the superconducting gap. Were a spectrum discrete, the adiabaticity condition would imply the level spacing being much larger than the driving frequency. The level-spacing is zero for a continuous spectrum, so this complicates the adiabaticity conditions. This has been pointed out already in Ref.Riwar et al. (2016) but was not investigated in detail. We note the generality of the situation: a generic gapped system might have a continuous spectrum above the certain threshold, and the adiabaticity condition required for the manifestations of topology needs to be revisited in this situation.

The aim of the present article is to investigate this question in detail for a generic model of a superconducting nanostructure. We have studied the linear response of currents on the changes of superconducting phases in the terminals. We model a multi-terminal superconducting nanostructure within the scattering approachNazarov and Blanter (2009). In this approach the terminals of the nanostructure are described with semiclassical Green’s functions and the scatterer coupled to the terminals is described by a unitary (in real time) S-matrix. Although it is not crucial, we made use of Matsubara formalism which conveniently allows us to concentrate on the ground state of the system and the limit of zero temperature is formally achieved by considering continuous Matsubara frequencies. So we do the calculations in imaginary time formalismSchön and Zaikin (1990). At the first step, we obtain the general effective action describing the nanostructure in terms of the S-matrix and time-dependent semiclassical Green’s functions of the terminals. At the second step, we expand the action to the second power in time-dependent phases of the terminals. At the third step, we concentrate on the limit of small voltage and driving frequency, to obtain the response function relevant for topological properties.

We can use the properly anti-symmetrized response function as a generalized definition of the Berry curvature that is suitable for the systems with and without a continuous spectrum. The main result of the present article is that so-defined Berry curvature is contributed to by a continuous spectrum as well as discrete one even in the case of energy-independent S-matrix. We derive an explicit formula for it. This solves the paradox mentioned in Riwar et al. (2016): the Berry curvature associated with discrete Andreev bands is discontinuous when the highest Andreev bound state merges with the continuum, which indicates that the integral of the Berry curvature defined only for discrete spectrum will not reduce to an integer. The redefined Berry curvature that we find is continuous. It gives rise to integer Chern numbers if the S-matrix is time-reversible. If it does not we reveal a specific additional non-topological contribution that does not depend on the superconducting phases. We note the the importance of the energy scales much larger than superconducting gap $|\Delta|$ in this context. This is why we also discuss in detail the case of an energy-dependent S-matrix the energy scale of variation of which may be in any relation with superconducting gap. We find that the non-topological contribution depends on the regularization of the S-matirx at large energies. In particular, it vanishes if the S-matrix is regularized as $S_{\pm\infty}=1$, this corresponds to no conduction between the terminals.

The paper is organized as follows. In Sec. II we introduce the details of a model of a multi-terminal superconducting nanostructure and review the main aspects of a scattering matrix approach formalism in this case. The derivation and discussion of the response function are given in Sec. IV. In Sec. VI we discuss the specific behaviour near the Weyl singularities, in the absence and presence of a weak spin-orbit coupling. In Sec. V we apply the general formulae to the case of a scattering matrix that varies only slightly on the scale of the superconducting gap $|\Delta|$. In Sec. VII we address the energy-dependent S-matrices at arbitrary energy scale for a specific model of an energy dependence. We conclude the paper with the discussion of our results (Sec. VIII). The technical details of the derivations are presented in Appendices.

II Multi-terminal superconducting nanostructure

Generally a multi-terminal superconducting nanostructure (Fig. 1) is a small conducting structure that connects $n$ superconducting leads. The leads are macroscopic and are characterized by the phases of the superconducting order parameter. Each lead labeled by $\alpha\in\{0,1,\cdots,n-1\}$ has its own superconducting phase $\phi_{\alpha}$ and one of the leads’ phase can be set to zero value $\phi_{0}=0$, according to the overall gauge invariance. The nanostructure design and these phases determine the superconducting currents $I_{\alpha}$ in each lead, that are the most relevant quantities to observe experimentally.

We aim to describe a general situation without specifying the nanostructure design. To this end, we opt to describe the system within the scattering approach pioneered by Beenakker Beenakker (1997). The superconducting leads are treated as terminals: they are regarded as reservoirs which contain macroscopic amount of electrons and are in thermal equilibrium. A common assumption that we also make in this article is that all terminals are made from the same material and thus have the same modulus of the superconducting order parameter $|\Delta|$. At sufficiently low temperatures and applied voltages one can disregard possible inelastic processes in the nanostructure and concentrate on elastic scattering only. Following the basics of the scattering approachNazarov and Blanter (2009), we assume $N_{\alpha}$ spin-degenerate transport channels in terminal $\alpha$. The conducting structure connecting the terminals is a scattering region and is completely characterized by a scattering matrix $S$ which generally depends on energy $\varepsilon$ and is a unitary matrix at any $\varepsilon$. In Matsubara formalism we use imaginary energy $\epsilon$ and the matrix $S$ satisfies the condition $S_{\epsilon}S_{-\epsilon}^{\dagger}=1$. All the details of the nanostructure design are incorporated into the scattering matrix.

The electrons and holes in the superconducting transport channels involved in the scattering process may be described as plane waves that scatter in the region of the nanostructure and then return to the corresponding terminals. Amplitudes of incoming and outgoing waves are linearly related by the $S$-matrix. The numbers of transport channels in the terminal $\alpha$ denoted as $N_{\alpha}$ determines the dimension of the scattering matrix: $\textrm{dim}S=M\times M$, where $M=2_{S}\sum_{\alpha}N_{\alpha}$ and $2_{s}$ counts for the spin.

The electrons and holes experience Andreev reflection in the superconducting terminals: the electrons are converted into holes and turn back, the same happens to holes. The Andreev reflection is complete at the energies smaller than the superconducting gap $\Delta$. Therefore, electron-hole waves may be confined in the nanostructure giving rise to discrete energy levels called Andreev bound states (ABS). The amplitudes and phases of these confined states are determined by the scattering matrix and Andreev reflection phases that involve the superconducting phases of the corresponding terminals. One can find the energies of the ABS $\varepsilon$ through Beenakker’s determinant equationBeenakker and van Houten (1991):

[TABLE]

where $S_{\varepsilon}$ is the S-matrix at the real energy $\varepsilon$, $\sigma_{y}=\begin{pmatrix}0&-i\\ i&0\end{pmatrix}$ is a Pauli matrix acting in the spin space and $e^{i\phi}$ is the diagonal matrix in channel space ascribing the stationary superconducting phases of the terminals to the corresponding channels, $e^{i\phi}\to\delta_{ab}e^{i\phi_{\alpha}}$ where $a,b$ label the channels and $\alpha$ is the terminal corresponding to the channel $a$. The ABS energies and the corresponding eigenvectors in the space of the channels depend parametrically on $n-1$ independent phases $\phi_{\alpha}\in[0,2\pi]$ and thus can be viewed as a bandstructure defined in a ”Brilluoin zone” of phases. It was notedRiwar et al. (2016) that (without spin-orbit interaction) three independent parameters are needed to tune the $n-1$ dimensional band structure of energy levels of ABS to reach the Weyl singularity at zero energy. It was also notedRiwar et al. (2016) that only one parameter is required to satisfy the condition for the highest ABS to touch the continuum above the gap ($\varepsilon=|\Delta|$). The ABS merges the continuum in this case and this implies that one cannot change this level adiabatically even for arbitrarily slow change of the parameters. When the incommensurate small voltages are applied to two terminals to sweep the phasesRiwar et al. (2016), the system passes the points where the highest level merges with the continuum. This makes it questionable to apply the adiabaticity reasoning in this case. This makes it necessary to consider the contribution of the continuous spectrum to the response function of the currents in the limit of slow change of the parameters.

III Action

The most general way to describe the nanostructure under consideration is to use an action method. This method has been pioneered in the context of a simple Josephson junction in Schön and Zaikin (1990). In this method one deals with an action of the nanostructure that depends on the time-dependent superconducting phases $\phi_{\alpha}(\tau)$. The transport properties of the nanostructure as well as quantum fluctuations of the phases in case the nanostructure is embedded in the external circuit Schön and Zaikin (1990), can be derived from this action.

One of the advances of this Article is the derivation of such action for multi-terminal nanostructure and arbitrary S-matrix in Matsubara formalism. The details of the derivation are given in IX. Here we give the answer:

[TABLE]

here $\Pi_{\pm}$ and $\hat{S}_{\epsilon}$ are matrices in a space that is a direct product of the space of channels, the imaginary-time space, spin and Nambu space. The matrix $\hat{S}_{\epsilon}$ is diagonal in the corresponding energy representation, therefore it depends on the difference of the imaginary time indices only. Its Nambu structure is given by

[TABLE]

where $S_{\epsilon}$ is the electron energy-dependent S-matrix (see App. IX). The matrix $g$ is composed of the matrices diagonal in energy and diagonal in time in the following way:

[TABLE]

where

[TABLE]

where $\tau_{z}$ is the 3rd Pauli matrix acting in Nambu space and the Nambu structure has been made explicit in $U^{\dagger}$. This form assumes that $|\Delta|$ is the same in all the terminals. If it is not so, the matrix $A_{\epsilon}$ also acquires the dependence on the channel index. It is worth noting that $g^{2}=1$ so that $\Pi_{\pm}$ are projectors. The matrix $g$ can be associated with the semiclassical Green’s function in a terminalEilenberger (1968); Nazarov and Blanter (2009): $e^{i\phi(\tau)}$ is the diagonal matrix in channel space ascribing the time-dependent superconducting phases of the terminals to the corresponding channels, $e^{i\phi(\tau)}\to\delta_{ab}e^{i\phi_{\alpha}(\tau)}$ where $a,b$ label the channels and $\alpha$ is the terminal corresponding to the channel $a$. We note the gauge invariance of the action: due to the invariance of the trace under unitary transformations, the superconducting phases can be ascribed to the terminal Green’s functions $g$ as well as to the scattering matrix. Let us assume that the matrix $S_{\epsilon}$ does not depend on spin. Then the trace over spin is trivial. It is convenient to apply the unitary transformation $U^{\dagger}$ as in $\eqref{g}$ to all the matrices in $\eqref{Bigaction}$. This transforms the matrix $g$ to $\tau_{z}$. Then the projectors take a simple form $\Pi_{\pm}\to\frac{1\pm\sigma_{z}}{2}$ and the matrix in $\eqref{Bigaction}$ reduces to the lower block-triangular form in Nambu space. The determinant is then equal to the determinant of the lower right block of the transformed matrix $\bar{S}_{\epsilon}$. Then the action takes the form

[TABLE]

[TABLE]

the S-matrix in Matsubara formalism is subject to the unitarity constraint,

[TABLE]

In what follows we concentrate on the zero-temperature limit $k_{B}T\ll|\Delta|$, so the summations over discrete frequencies are replaced with integrations $\int\frac{d\epsilon}{2\pi}$.

III.1 Stationary phases

In the stationary case $\phi(\tau)=\phi+\delta\phi(\tau)$ with constant $\phi$ and $\delta\phi(\tau)\equiv 0$ the value of the action gives the stationary phase-dependent ground state energy of the nanostructure $E_{g}=\lim_{k_{B}T\to 0}TL_{0}$.

[TABLE]

[TABLE]

where Trace is now over the channel space and the Trace over spin space is taken explicitly as a factor of $2_{S}$ unless specifically addressed. The operator $Q_{\epsilon}$ introduced here has the properties of the inverse of the Green’s function although it is not related to an operator average: its determinant as function of complex $\epsilon$ vanishes, ${\rm det}Q_{\epsilon}=0$, at imaginary values $\epsilon=\pm i\varepsilon_{k}$ corresponding to the ABS energies (compare with $\eqref{beenakker}$). In addition to these singularities the operator $Q_{\epsilon}$ has two cuts in the plane of complex $\epsilon$ corresponding to the presence of a continuous spectrum in the terminals above the gap $|\Delta|$. We choose the cuts as shown in Fig. 2. The expression $\eqref{GSE}$ can be simplified in the case when the S-matrix does not depend on energy

[TABLE]

the second (divergent) contribution here does not depend on the superconducting phases so we omit it. To compute the integral it is convenient to choose the basis in which the unitary matrix $\Lambda=SS^{*}$ is diagonal. This is a unitary matrix, so the eigenvalues are unimodular complex numbers. The phases of the eigenvalues are related to the energies of ABS: $\Lambda_{k}=e^{2i\chi_{k}},\quad\chi_{k}=\arccos[\epsilon_{k}/|\Delta|],\quad\chi\in[-\pi/2;\pi/2]$. The eigenvalue $\Lambda_{k}=1$ is doubly degenerate and corresponds to the values $\epsilon_{k}=\pm|\Delta|$. The eigenvalues come in complex conjugated pairs $\Lambda_{k}^{*}=\Lambda_{-k}$, where $(-k)$ corresponds to the Nambu-counterpart of the $k-$th eigenvector. So only the eigenvalues ${\rm Im}\Lambda_{k}>0$ correspond to the quasiparticle states with positive energies. We will label them with positive indices $k$. In what follows we define a ”bar” operation that links these pairs $|\bar{k}\rangle=S|k^{\star}\rangle=|-k\rangle$ where $|k\rangle$ is some eigenvector of $\Lambda$. We note, however, that this operation is not a convolution, since $|\bar{\bar{k}}\rangle=\Lambda_{k}|k\rangle$.

In this basis we can rewrite the integral as

[TABLE]

Evaluation of the integral brings to the known result

[TABLE]

where $\epsilon_{k}$ are the stationary phase-dependent ABS energies, as discussed above. The derivative of the ground state energy with respect to a stationary phase in terminal $\alpha$ gives the stationary current in the corresponding terminal,

[TABLE]

We expect this relation to hold in the adiabatic limit. In the following Section, we will access the time-dependent currents concentrating on the next order correction in the limit of small frequencies.

IV Response function of the currents

To compute the response function of the currents we assume small nonstationary phase addition to the stationary phases $\phi$, $\phi(\tau)=\phi+\delta\phi(\tau)$, $\delta\phi(\tau)\ll 2\pi$ and expand the action to the second order in $\delta\phi(\tau)$ (first order vanishes automatically since $\delta\phi(\tau)$ is nonstationary $\int_{0}^{\beta}d\tau\delta\phi(\tau)=0$). We give the details in Append. X. The total contribution to the action reads

[TABLE]

$\delta\phi_{\omega}$ being the Fourier transform of $\delta\phi(\tau)$. The frequency-dependent response function of the current $R_{\omega}^{\alpha\beta}$ is given by

[TABLE]

here the stationary phases are ascribed to the S-matrix. We use a shorthand notation $\partial/\partial_{\alpha}=\partial/\partial_{\phi_{\alpha}}$ and define a set of matrices that project channel space onto the space of the channels in the terminal $\alpha$, $(I^{\alpha})^{ab}=\delta^{ab}$ if $a$ is a channel in terminal $\alpha$ and $(I^{\alpha})^{ab}=0$ otherwise. The first term in $\eqref{general}$ vanishes at zero frequency and in the case of the energy-independent S-matrix. The second term does not depend on frequency $\omega$. In the limit of zero frequency the second and the third terms reproduce the stationary response function of the currents

[TABLE]

Let us consider the limit of small $\omega\ll|\Delta|$ and concentrate on the first order correction to the adiabatic limit

[TABLE]

We note that the response function is analytic in the vicinity of $\omega=0$. This is guaranteed by the gap in the density of states, which is given by the energy of the lowest ABS. Away from the zero-energy Weyl singularity it can be estimated as $|\Delta|/N$ with $N$ being the total number of ABS in the nanostructure. The vicinity of a Weyl singularity has to be treated more carefully as we discuss in Sec. VI. Let us note that for any system with a discrete spectrum the quantity $B_{\alpha\beta}$ can be related to the Berry curvatureM. V. Berry (1984); Niu et al. (1985); Thouless et al. (1982). For any state in the discrete spectrum the Berry curvature corresponding to this state is given by $B_{\alpha\beta}^{(i)}=2{\rm Im}\langle\partial_{\alpha}i|\partial_{\beta}i\rangle$ with $i$ labeling discrete states and $|i\rangle$ being the wavefunction of the corresponding state. In our case we are interested in the total Berry curvature of the superconducting ground state defined as $B_{\alpha\beta}=-\frac{1}{2}\sum_{i}B_{\alpha\beta}^{(i)}$ where $i$ labels the (spin-degenerate) wavefunctions of the BdG equation with positive eigenvaluesRiwar et al. (2016). However, the adiabaticity condition which justifies the expansion in $\eqref{linear}$ for the case of discrete spectrum requires the frequency to be much smaller than the smallest energy spacing between the levels.

In our system, the continuous spectrum above the superconducting gap is present. In principle, any continuous spectrum can be approximated with a discrete spectrum with a vanishing level spacing $\overline{\delta}\to 0$. By doing this we can utilize the previous expression for the response function $B_{\alpha\beta}$ since it is valid for the discrete spectrum. However, the adiabaticity condition which is necessary for this expression to be valid would reduce to $\omega\ll\overline{\delta}\to 0$. This condition contains an artificially introduced $\overline{\delta}$ and is by construction very restrictive in $\omega$. On the other hand, the expansion in Eq. $\eqref{linear}$ is valid under a physically meaningful and less restrictive condition $\omega\ll|\Delta|/N$. Taken all that into account, we conclude that the response function $B_{\alpha\beta}$ defined in Eq $\eqref{linear}$ does not have to reduce to the expression for a total Berry curvature of a superconducting ground state of a system discussed above. The topological properties of this quantity also have to be investigated separately.

One may conjecture that the resulting response function in Eq. $\eqref{linear}$ reduces to the sum of the Berry curvatures of the discrete ABS spectrum, so that it is not contributed to by the continuous spectrum. This conjecture relies on the analogy between the expressions for the total Berry curvature and the superconducting ground state energy. In the case when the S-matrix is energy-independent, only the discrete states contribute to the ground state energy. Thus motivated, in the following we investigate the response function $B_{\alpha\beta}$ defined by means of Eq $\eqref{linear}$ in detail. We find that there is a contribution from the continuous spectrum to this quantity as well as from the discrete one. We also find that in general the integral of $B_{\alpha\beta}$ over the phases $\phi_{\alpha},\phi_{\beta}$ that would normally define an integer Chern number, is not integer. Therefore, $B_{\alpha\beta}$ contains a non-topological contribution. This non-topological part is contributed by the continuous as well as the discrete part of the spectrum.

The tensor $B_{\alpha\beta}$ defined in Eq. $\eqref{linear}$ is antisymmetric (since $R^{\alpha\beta}_{\omega}=R^{\beta\alpha}_{-\omega}$). The concrete expression for $B_{\alpha\beta}$ reads:

[TABLE]

The first term here resembles the usual WZW formEssin and Gurarie (2011) for a Chern number. Usually, the form contains the matrix Green’s functionsEssin and Gurarie (2011), in our case the form utilizes the matrix $Q_{\epsilon}$ defined by Eq. $\eqref{qmatrix}$. We note however that in distinction from common applications of WZW forms here one cannot regard $Q_{\epsilon}$ as a smooth function of parameters $\phi_{\alpha},\phi_{\beta},\epsilon$ defined on a compact manifold without a boundary. This is because in general this matrix has different limits at positive and negative infinite energies $S_{-\infty}$ for $\epsilon\to-\infty$ and $S_{-\infty}^{T}$ for $\epsilon\to+\infty$ that also depend on the phases. Due to this reason the integral of the first term over a compact surface without a boundary in a space of phases does not have to reduce to an integer $\cdot(2\pi)^{-1}$. The second term in Eq. $\eqref{result1}$ has a form of a total derivative with respect to a phase of a periodic and smooth function, so the integral of this one over a compact surface will give zero.

In order to obtain the value of this integral let us consider first the variation of this value upon the small smooth variation of the matrix $Q_{\epsilon}\to Q_{\epsilon}+\delta Q_{\epsilon}$ that comes from the small variation of the S-matrix $\delta S_{\epsilon}$, so $\delta Q_{\epsilon}=A_{\epsilon}^{2}\delta S_{\epsilon}+A_{-\epsilon}^{2}\delta S_{-\epsilon}^{T}$. The value of the integral of the second contribution in Eq. $\eqref{result1}$ does not contribute to the integral over a compact submanifold in phase space, so we needn’t consider its variation. It is known Qi et al. (2008) that the variation of the first contribution to $B_{\alpha\beta}$ reduces to the total derivatives

[TABLE]

The value of the integral of second term in $\eqref{varder}$ over a compact submanifold in phase space vanishes if the submanifold does not pass Weyl singularities corresponding to ${\rm det}Q_{\epsilon}^{-1}\to\infty$, because it has a form of a total derivative of a smooth function. Evaluation of the integral in $\eqref{varcond}$ yields the following contribution to the variation of $B_{\alpha\beta}$

[TABLE]

We note that this contribution is generally nonzero and does not depend on phases.

Let us turn to the evaluation of the topological charge that is proven to be very useful in the field Wang et al. (2010b). The value of the topological charge is defined in a usual way with the divergence of the topological field $\vec{E}$

[TABLE]

To compute the topological charge we need to consider a special variation of the S-matrix that just corresponds to the stationary phase derivative $\delta S_{\epsilon}=[S_{\epsilon},\frac{iI_{\gamma}}{2}]\delta\phi_{\gamma}$. Since the expression under the trace in $\eqref{landvar1}$ does not depend on phases, the topological charge vanishes at any point where the field $\vec{E}$ is well-defined, or alternatively ${\rm det}Q_{\epsilon}^{-1}$ is finite. The Weyl singularities give rise to the point-like integer charges being the sources of the field $\vec{E}$. We consider this in detail in Sec. VI. This situation is in complete analogy with that of the standard Berry curvature of a discrete spectrum where Weyl singularities correspond to band crossings. However, we have computed the topological charge for the particular phase-dependence of the S-matrix on phases ($e^{-\frac{i\phi}{2}}Se^{-\frac{i\phi}{2}}$). We have not considered the topological charge in the space of 2 phases $\phi_{\alpha},\phi_{\beta}$ and some other parameter characterizing the scattering matrix, this charge could be nonzero and have a continuous distribution. The investigation of the general parametric dependence of the S-matrix is beyond the scope of the present article.

We separate the field $\vec{E}$ into three parts: a part produced by the point-like charges, divergenceless field that is zero in average, and a constant part $\vec{\bar{E}}$. The value of the integral

[TABLE]

is given by the flux of the topological field through the corresponding surface. This flux reduces to the integer for the first contribution to $\vec{E}$, vanishes for the second divergenceless contribution and may result in some value for the constant part of the field. We stress that the last contribution being present is the main distinction from the common case. The value of this constant field is then given by the integration of the variation $\eqref{landvar1}$:

[TABLE]

This constant field can contribute to the flux through any plane in the phase space.

[TABLE]

where $\vec{n}$ is the normal vector to this plane. As it has been shown in Ref. Riwar et al. (2016) the value of $C^{12}$ is directly related to the observable transconductance between the leads $\alpha$ and $\beta$. Therefore, in contrast to the conclusions of Ref. Riwar et al. (2016) the value of transconductance does not always quantize although the change of transconductance with a phase can be quantized.

So, in principle a nonzero non-topological contribution to $\eqref{flux}$ can be present. This contribution is nonzero if the S-matrix is not regularized at infinite energy such that $[S_{-\infty},I_{\alpha}]=0$. If the S-matrix is regularized in this way, then the $Q_{\epsilon}$ matrix is defined on a compact space of parameters $(\epsilon,\alpha,\beta)$, so the first contribution to Eq. $\eqref{result1}$ would reduce to an integer $n$ (with proper normalization). If it is not regularized this way, then this boundary term leads to the presence of a non-topological contribution to the response function, that comes due to the presence of a continuous spectrum and, formally, from the fact that the matrix $Q_{\epsilon}$ is not defined on a compact space, as discussed above. In the limit of energy-independent S-matrix, this contribution reduces to the antisymmetric part of the Landauer conductanceNazarov and Blanter (2009); Landauer (1957). In this case, if the bare S-matrix (without the stationary phases of terminals ascribed) is non-symmetric (which means the breaking the time-reversibility condition) we obtain a nonzero value of $\eqref{landauer}$. If the S-matrix is time-reversible, the non-topological contribution is zero and the integer quantization of transconductance is restored.

V Weak energy dependence of the $S-$matrix

In the description of the realistic nanostructure a reasonable approximation is to consider the S-matrix to be constant on the scale of $|\Delta|$. It corresponds to the case of a short nanostructure (smaller than the superconducting coherence length). So a logical approximation would be to describe the nanostructure with a constant S-matrix at all energies. The response function $B_{\alpha\beta}$ is given by an integral over energy in Eq. $\eqref{result1}$. Would this integral accumulate in the region $\epsilon\sim|\Delta|$, then the approximation of a constant S-matrix at all energies would be accurate. However, there can be a significant contribution from the energy scales $\epsilon\gg|\Delta|$ to the integral yielding $B_{\alpha\beta}$. In this case the energy dependence of the S-matrix at the large energies becomes important. To investigate this we consider the contributions from the small scales $\epsilon\gtrsim|\Delta|$ and from the large scales $\epsilon\gg|\Delta|$ in the Subsections V.1 and V.2 respectively.

V.1 Energy-independent $S-$matrix:

In this Subsection we analyze the small-scale ($\epsilon\sim|\Delta|$) contribution to $\eqref{result1}$. For this we approximate the S-matrix to be constant at all energies and extend the integration limits to infinity. The second term in $\eqref{result1}$ vanishes since $\frac{\partial S_{\epsilon}}{\partial\epsilon}=0$. The integral in the first term in $\eqref{result1}$ converges on the scale $\epsilon\gtrsim|\Delta|$. This statement only necessarily holds if the S-matrix is energy-independent. Otherwise, the contribution from the larger scales can be present and we investigate it in V.2. Similarly to $\eqref{gseint}$, the result of integration under consideration can be expressed in terms of the eigenvalues and eigenvectors of the unitary matrix $\Lambda=SS^{*}$. We use the same notations $\left|k\right\rangle$ and $\left|\bar{k}\right\rangle$ for the eigenvectors related to the complex conjugated eigenvalues pair $\Lambda_{k}$ and $\Lambda_{k}^{*}$ correspondingly as described after Eq.$\eqref{holycow}$. We remind that the phase of the eigenvalue $\Lambda_{k}=e^{2i\chi_{k}}$ with $k>0$ is related to the energy of ABS as $\chi_{k}=\arccos[\epsilon_{k}/|\Delta|]$. We also remind that $\Lambda_{k}=1$ is degenerate and corresponds to the energy of one of the ABS $\epsilon_{k}=|\Delta|$. Upon crossing this point in phase space, this ABS state exchanges the wave function with its Nambu counterpart with the eigenvalue $\epsilon_{k^{\prime}}=-|\Delta|$. Due to this we call such points gap touching singularities.

Evaluating the integral yields

[TABLE]

where $k,j$ label the eigenvalues of $\Lambda$, and the summation goes over indices with both signs. If the number of channels is odd, there is an eigenvector of $\Lambda$ corresponding precisely to the eigenvalue $\Lambda_{k}=1$. Then the index $k=0$ corresponds to this state. If the number of channels is even, the indices in Eq.$\eqref{result2}$ do not take the zero value. In the following we consider the number of channels to be even. The logarithm here has a branch cut along the real axis as $[0,+\infty]$ (see Fig. 3) to avoid the gap touching singularity ambiguity $\Lambda_{k}=1$. Let us consider the behaviour of $B_{\alpha\beta}$ in the vicinity of the gap touching singularity. Since the wave function corresponding to $\Lambda_{k}\to 1+i0$ is discontinuous upon crossing this singularity, it is not obvious that $B_{\alpha\beta}$ is continuous. However, one can observe that the first term is a sum of Berry curvatures of individual levels multiplied by the eigenvalue-dependent prefactors $\log\Lambda_{k}$. This prefactors vanish for the discontinuous wavefunctions at the gap touching degeneracy and guarantee the continuity of the first term. Also, one can show that the second term in Eq.$\eqref{result2}$ is continuous. Consequently, $B_{\alpha\beta}$ is continuous at this point (see Fig. 4). The only possibility for $B_{\alpha\beta}$ to be ill-defined at some points in phase space is the zero-energy Weyl singularity where $\det Q_{\epsilon}^{-1}$ diverges (see Sec.VI).

The response function $B_{\alpha\beta}$ is expressed in terms of eigenvalues and eigenvectors of the matrix $\Lambda$. So is the ABS contribution to the ground state Berry curvature, which was conjectured as a result for $B_{\alpha\beta}$ (see Sec. IV). It was shownRiwar et al. (2016) that this ABS contribution is given by $B^{\rm ABS}_{\alpha\beta}=-\frac{2_{S}}{2}\sum_{k>0}B^{(k)}_{\alpha\beta}$, $B^{(k)}_{\alpha\beta}=2{\rm Im}\langle\partial_{\alpha}k|\partial_{\beta}k\rangle$. Since one of the wavefunctions contributing to this sum is discontinuous at the gap touching singularity, we conclude that $B^{\rm ABS}_{\alpha\beta}$ is discontinuous contrary to $B_{\alpha\beta}$. One can understand the difference between $B_{\alpha\beta}$ and $B^{\rm ABS}_{\alpha\beta}$ by considering the computation of the integral in the first term in Eq. $\eqref{result1}$ by means of complex analysis (in the plane of complex $\epsilon$). By shifting the integration contour to the upper half-plane, one can see that the integral is contributed to by the poles, corresponding to ABS and the cut above the gap (see Fig. 2). The contribution from the poles results in $B_{\rm ABS}$, but the contribution from the cut, $B_{\alpha\beta}^{\rm cut}=B_{\alpha\beta}-B_{\alpha\beta}^{\rm ABS}\neq 0$, is equally important (see Fig. 4).

For the integrated $B_{\alpha\beta}$ we obtain in accordance with Eq. $\eqref{flux}$

[TABLE]

so the value of transconductance is not necessarily quantized in the approximation of the energy-independent S-matrix.

V.2 Contribution from the large scales

In the previous Section we have shown that the non-topological contribution to the transconductance comes from the boundary terms at $\epsilon=\pm\infty$ (see Eq.$\eqref{varcond}$). This means that, contrary to intuition, there is an essential contribution to $B_{\alpha\beta}$ coming from the energy scales much larger than the energy gap. In order to investigate the large energy contribution we assume the regularization of the S-matrix at large energies. So, in this Subsection we consider $B_{\alpha\beta}$ for a particular energy-dependence of the S-matrix. It is chosen such that the S-matrix is regularized at infinity such that it varies slowly on the scale of a superconducting gap $|\Delta|$ and $S_{\pm\infty}=1$. This S-matrix corresponds to a complete isolation of the terminals at the largest energies. With this regularization, the matrix $Q_{\epsilon}$ is defined on a compact parameter space $(\alpha,\beta,\epsilon)$ and the first contribution in $\eqref{result1}$ must reduce to an integer. Due to the scale separation, there are two contributions to $B_{\alpha\beta}$. One comes from the scales $\epsilon\sim|\Delta|$ and is given by the same result $\eqref{result2}$. Another one comes from the scales $\epsilon\gg|\Delta|$.

For negative energies, the large scale contribution with asymptotic accuracy equals

[TABLE]

with the notation $S=S_{\epsilon=0}$.

For positive ones:

[TABLE]

So, the both contributions give the following addition to the response function

[TABLE]

Both terms here do not depend on phases. The first one is exactly equal to the constant part of the topological field defined previously with an opposite sign (computed for an energy-independent S-matrix case). So after integration over two phases, it cancels the non-topological contribution from small scales in $\eqref{Smallint1}$. Since we assume a regularization $S_{\pm\infty}=1$, the second term is zero (${\rm Tr}[S_{+\infty}^{\dagger}\frac{I_{\alpha}}{2}S_{-\infty}\frac{I_{\beta}}{2}]=0$), so the total mean value of the transconductance is quantized in correspondence with the theory of characteristic classes.

The second contribution to $B_{\alpha\beta}$ in Eq. $\eqref{result1}$ contains the energy-derivative of the S-matrix under the integral. Due to this the energy scale of its dependence drops out from the integral. So, one may expect that it contributes to the large scale contribution to $B_{\alpha\beta}$. However, with asymptotic accuracy it vanishes in the limit when the S-matrix varies slowly on the scale $|\Delta|$. Indeed, in the limit $|\epsilon|\gg|\Delta|$

[TABLE]

In this limit for $\epsilon<0$, the integrand equals

[TABLE]

with asymptotic accuracy, since the expression under the trace does not depend on phases. For $\epsilon>0$ the integrand vanishes since $A_{\epsilon}^{2}\to 0$ for $\epsilon\gg|\Delta|$.

VI The vicinity of a Weyl point

In this Section, we investigate the Berry curvature in the vicinity of a Weyl singularity, that occurs at some point $\vec{\phi}_{0}$ in the 3-dimensional phase space. Such Weyl points have been analyzed in Riwar et al. (2016) assuming spin symmetry, in Yokoyama and Nazarov (2015) the analysis has been extended to cover weak spin-orbit interaction. Without spin-orbit coupling, the Weyl points are situated at zero energy and $\det Q^{-1}_{\epsilon=0}$ diverges near the point. A conical spectrum of ABS is found in the vicinity of the point Riwar et al. (2016). A weak spin-orbit coupling splits the energy cones in spin and shifts the Weyl point to a finite energy Yokoyama and Nazarov (2015). Further, we discuss separately the cases of vanishing and weak spin-orbit coupling.

VI.1 Vanishing spin-orbit coupling

When the spin-orbit (SO) coupling is absent, the Weyl singularities are located at some points in the phase space ${\vec{\phi}_{0}}$ and occur at zero energy $\epsilon_{\pm}=0$. To consider the vicinity of the singularity, we assume a small phase deviation $\delta\hat{\phi}=\hat{\phi}-\hat{\phi}_{0}\ll 1$ from the singularity point and assign it to each channel via the diagonal matrix $e^{\delta\hat{\phi}}$. In the vicinity, $B_{\alpha\beta}$ defined by Eq. $\eqref{result1}$ only has non-zero contributions from the first term of quasi-WZW term. The second term vanishes asymptotically when the energy approaches zero, as shown in Eq. 34. Conform to these approximations, we extend the domain of the integration over the phases to infinity since $B^{\alpha\beta}$ is concentrated near the singularity point.

To compute $B^{\alpha\beta}$, we approximate the $Q$ matrix near the Weyl point with the expression that keeps the first orders in $\epsilon$ and of the variation: $Q=(\epsilon+\frac{1}{2}\delta\Lambda)S^{T}=MS^{T}$, $S$ being the scattering matrix in the singularity point at $\epsilon=0$. Conveniently, we can replace $Q$ with $M$ in Eq.(21). We find the variation $\delta\Lambda$ by expanding the S-matrix in $\delta\vec{\phi}$:

[TABLE]

We can contract the dimension of $M$ projecting it to two eigenvectors of $\Lambda$ that achieve singular values at the Weyl point. Following Riwar et al. (2016), we separate the singular part of $M$ and write in the basis of ABS eigenvectors $\left|+\right\rangle$ and $\left|-\right\rangle$ satisfying $S\left|\pm\right\rangle=\pm\left|\mp\right\rangle^{*}$, $\Lambda|\pm\rangle=-|\pm\rangle$:

[TABLE]

where $\vec{\tau}$ are the Pauli matrices in the space of these two eigenvectors, and the components of $\vec{h}$ are proportional to the components of $\vec{\phi}$: $h_{x}+ih_{y}=2\left\langle-\right|\delta\hat{\phi}\left|+\right\rangle$, $h_{z}=\left\langle+\right|\delta\hat{\phi}\left|+\right\rangle-\left\langle-\right|\delta\hat{\phi}\left|-\right\rangle$.

The form of $M$ is similar to the generic form of Green’s function of a two-level system. We expect that the two poles of $M^{-1}$ should be positioned symmetrically on the imaginary axis $\epsilon$ due to BdG particle-hole symmetry. Indeed, we find these poles at $\epsilon_{\pm}=\pm i\frac{|\vec{h}|}{2}$. Using the trace relations of Pauli matrices, we reduce in the leading order $B_{\alpha\beta}$ to the Berry curvature of the corresponding levels :

[TABLE]

We note that in this section all the matrices have the spin index. For an $N$ dimensional space of superconducting phases, the singularities are concentrated in the $N-3$ dimensions and the relevant space is reduced to a $3$-dimensional subspace $\{\delta\phi_{1},\delta\phi_{2},\delta\phi_{3}\}$. For certainty, we set the indices $\alpha,\beta=1,2$, and consider the curvature defined in the $\phi_{1}-\phi_{2}$ plane at a fixed phase $\phi_{3}$.

The $\phi_{3}$ dependence of the integral of the curvature with respect to superconducting phases $\phi_{1}$, $\phi_{2}$ witnesses the change of first Chern number $C^{12}$ when the integration plane passes the singularity point. Since we only concentrate on the vicinity of the Weyl singularity, the integral under the approximations made can only indicate the change of the Chern number, rather than its total value that can be determined by integration over the regions far from the singularity point. To compute the integrated $B_{\alpha\beta}$, we notice from Eq.(39) that the energy spectrum is linear in $\delta\phi$, and introduce a linear relation $h_{i}=\sum_{\alpha}\delta\phi_{\alpha}T_{\alpha i}$ with $T_{\alpha i}=\partial_{\alpha}h_{i}$ being a real invertible matrix. The integrated $B_{12}$ is then obtained as:

[TABLE]

$\textrm{sgn}(\delta\phi_{3})$ determining the orientation of the $\delta\phi_{3}$ deviation.

This implies that whenever the integration plane passes the Weyl point, the first Chern number is changed by $\Delta C^{12}=\frac{1}{2}\textrm{sgn}(\delta\phi_{3}\det T)-\frac{1}{2}\textrm{sgn}(-\delta\phi_{3}\det T)=\pm 1$. This manifest the the integer values of the topological charge. The integrated $B_{\alpha\beta}$ in Eq.(41) specifies the flux of the Berry field penetrating the plane which is either above or below the singularity point. This flux, owing to symmetry, is a half of the total flux, this explains the half-integer values. Therefore, the main contribution to Eq.$\eqref{result2}$ in the vicinity the Weyl point is given by the Berry curvatures of the two levels that are close to zero energy, and can be presented as

[TABLE]

VI.2 Weak Spin-Orbit Coupling

Let us turn on a weak spin-orbit interaction and take it into account perturbatively giving a small spin-dependent change to the scattering matrix that preserves its unitarity, as is done in Yokoyama and Nazarov (2015). The first order variation thus reads

[TABLE]

where the last equality sign implies the commutation relation $\sigma_{y}\sigma_{i}^{*}\sigma_{y}=-\sigma_{i}$. Here, $\vec{\sigma}$ are the Pauli matrices in spin space and $\vec{K}$ being the corresponding Hermitian matrix in the channel space characterizing the spin-orbit effects. Owing to the time reversibility, $\vec{K}(\vec{\phi})=-\vec{K}(-\vec{\phi})$, yet in the vicinity of the singularity we may disregard its dependence on superconducting phases.

As in the previous Subsection, we project the matrix $Q$ onto singular subspace that has now dimension of $4$ to account for spin, and replace it with the matrix $M$. Writing the latter in the basis of eigenvectors $\left|\pm\right\rangle\left|\uparrow(\downarrow)\right\rangle$:

[TABLE]

$\vec{K}^{\prime}=\left\langle+\right|\vec{K}^{*}\left|+\right\rangle+\left\langle-\right|\vec{K}^{*}\left|-\right\rangle$. We can conveniently choose the spin quantization axis in the direction of $\vec{K}^{\prime}$ replacing the operator $\vec{\sigma}\cdot\vec{K^{\prime}}$ with its eigenvalues $\pm|K_{0}|=\pm\sqrt{|\vec{\sigma}\cdot\vec{K}^{\prime}|}$ for spin up and down, respectively.

The spin-orbit coupling lifts the spin degeneracy of the ABS in the vicinity of a Weyl point. The poles at imaginary energies become $\epsilon_{\uparrow}=i(\pm\frac{|\vec{h}|}{2}+\frac{|K_{0}|}{2})$ for spin up and $\epsilon_{\downarrow}=i(\pm\frac{|\vec{h}|}{2}-\frac{|K_{0}|}{2})$ for spin down. Contrary to the spin-degenerate case, the singularities at $|\vec{h}|=0$ are no longer at zero energy. Instead, they are shifted to $\pm i|K_{0}|$, see Fig. 7. The conical singularity of the spectrum remains and the topology is still protected, as we will explain below in detail.

The ABS energies cross zero energy when

[TABLE]

is satisfied. (Here, we introduce a positively defined matrix $X_{\alpha\beta}=\sum_{i}T_{\alpha i}T_{i\beta}$. Eq.(46) defines an ellipsoidal surface in the 3D superconducting phase space that encloses the singularity at $\hat{\phi}_{0}$ where $|\vec{h}|=0$. Outside the ellipsoid, two positive imaginary poles at $\epsilon_{+\uparrow(\downarrow)}=\frac{i}{2}(|\vec{h}|\pm|K_{0}|)$ hold a half of the residue of the spin degenerate pole $\epsilon_{+}$ each. Two negative imaginary poles $\epsilon_{-\uparrow(\downarrow)}$ at $\epsilon_{-\uparrow(\downarrow)}=\frac{i}{2}(-|\vec{h}|\pm|K_{0}|)$ have the opposite residues. Inside the ellipsoid, poles of $\epsilon_{+\uparrow}$ and $\epsilon_{-\downarrow}$ exchange their values as well as wave functions, thus canceling the contributions from the other two poles. Thus, $B_{\alpha\beta}$ is zero inside the ellipsoid and is the same as in the spin-degenerate case outside the ellipsoid,

[TABLE]

The result of integration of $B^{12}$ over two superconducting phases $\phi_{1}$,$\phi_{2}$ at a fixed $\delta\phi_{3}$ thus reads

[TABLE]

One can understand this result geometrically by presenting Eq. (49) as an integral over the corresponding plane in $\vec{h}$ space,

[TABLE]

where $\hat{n}_{12}^{h}$ is the vector normal of the corresponding plane and $\Omega_{12}$ is eventually the solid angle at which a part of the $\phi_{1}-\phi_{2}$ plane outside the ellipsoid is seen from the Weyl singularity (see Fig.7). Generally, this angle is expressed through elliptic integrals.

The integral can be simplified if we choose the coordinate system in 3D space of the phases in such a way that $T_{13}=T_{31}=T_{23}=T_{32}=0$. With this, the integral can be evaluated as

[TABLE]

We see that in the vicinity of a Weyl point the $C^{12}$ is not a topologically protected quantity confined to the integer values: rather, it changes linearly in an interval of $\delta\phi_{3}$ defined by the strength of the spin-orbit coupling (Fig. 7 )

To explain this, and eventually restore the topological protection of $C_{12}$, let us consider many-body states in the vicinity of the Weyl point. Their energies are given by the eigenvalues of the many-body Hamiltonian $H_{\rm MB}$

[TABLE]

where $E_{\uparrow(\downarrow)}=|\vec{h}|\pm|K_{0}|$ are the energies of quasiparticle excitations with spin up(down), $\hat{n}_{\uparrow(\downarrow)}$ are the number operators of the quasiparticles with the corresponding spin. The energy spectrum $E_{\rm MB}$ for each of the four possible states is given in Fig. 8. As we see from the Figure, the ground state of the superconducting nanostructure corresponds to ${n}_{\uparrow}={n}_{\downarrow}=0$ at $|\vec{h}|>|K_{0}|$ and to ${n}_{\downarrow}=1,{n}_{\uparrow}=0$ within the ellipsoid $|\vec{h}|<|K_{0}|$. These states differ in fermion parity, that is the conserving quantity for the superconducting Hamiltonian. This is why the parity transition that takes place at $|\vec{h}|=|K_{0}|$ is accompanied by the discontinuity of the wave functions, which violates the topological quantization of $C^{12}$. It is evident from Fig. 8 that the states of the odd fermion parity do not depend on phases in the vicinity of the Weyl point therefore corresponding to zero $B^{12}$.

The topological protection is restored if one considers the ground state at fixed parity. Then for the even ground state $C^{12}$ is the same as for the spin-degenerate case and experiences an integer jump when the integration plane passes the singularity point. No change of topological charge occurs for the odd ground state and it remains topologically trivial.

VII Energy-dependent $S-$matrix

In this Section we consider the effect of the energy dependence of the S-matrix on $B_{12}$ given by $\eqref{result1}$ for arbitrary relation between the energy scales of the scattering matrix and the gap $|\Delta|$.

We make use of the following model scattering matrix:

[TABLE]

where $\hat{\Gamma},\hat{H}$ are Hermitian dimensionless matrices with eigenvalues of the order of one. This expression can be regarded as a rather general polar decomposition of an energy-dependend scattering matrix. Since the matrices $\hat{\Gamma},\hat{H}$ can be diagonalized simultaneously, the expession has poles at the complex energies $E=\mu+\mathscr{E}({H}_{n}-i{\Gamma}_{n}/2)$ defined by the corresponding eigenvalues. The poles can be seen as the scattering resonances. The eigenvalues $H_{n}$ set the energies of those resonances and the corresponding eigenvalues $\Gamma_{n}$ give the inverse lifetimes of these resonances, $\Gamma_{n}$ must be positive to assure the correct causal properties of the scatterimg. Real energy scale $\mathscr{E}$ then sets the typical spread of the poles in energy around their average position $\mu$. We note that $S_{\epsilon}\to 1$ as $\epsilon\to\infty$, so the conditions of regularization described in a previous Section are fulfilled and the integral of $B_{12}$ over a compact subspace in phase space that does not cross the Weyl singularities, reduces to an integer. We remind that the limit $S_{\epsilon}\to 1$ corresponds to isolated terminals.

In distinction from the weak energy dependence case, the ABS energies defined by Eq. $\eqref{beenakker}$ can not be readily obtained and the resulting spectrum may be complicated with more ABS per transport channel. It is no more plausible to separate the contributions to $B_{\alpha\beta}$ coming from discrete and continuous spectrum. This, however, does not change the qualitative features of these contributions discussed above.

Let us consider and illustrate the dependence of $B_{12}$ on these two energy scales. We choose random matrices $\hat{H}$, $\hat{\Gamma}$ that satisfy the conditions stated, and compute $B_{12}$ from Eq. 21 at rather arbitrary settings of 3 phases. The integration over the imaginary energy in Eq. 21 permits the evaluation with no regard for the details of a complicated ABS spectrum. We plot the result versus the energy scale $\mathscr{Epsilon}$ at several settings of $\mu$. (Fig. 9)

Let us consider $\mu\neq 0$ first. In this case, at $\mathscr{E}\to 0$ the transmission between the terminals is limited to a small circle of the radius $\simeq\mathscr{E}$ near $\mu$. This suppresses the Andreev scattering that requires good transmission at opposite energies, and all quantities that depend on the phase differences including $B_{\alpha\beta}$. In Fig. 9, this is manifested as almost zero $B_{12}$ at $\mathscr{E}<\mu$. The further increase of $\mathscr{E}$ restores the Andreev scattering bringing $B_{12}$ to its typical values of $\sim(2\pi)^{-2}$. We note a non-monotonous dependence on $\mathscr{E}$ and explain it by the fact that different poles of the scattering matrix contribute to $B_{12}$ with typically different signs, and the magnitude of the contribution depends on the position of the pole with respect to the energy scale $\simeq\Delta$. At $\mathscr{E}\gg\Delta$ the energy dependence of the scattering matrix is weak at $\epsilon\simeq\Delta$ and $B_{12}$ saturates at a value that does not depend on $\mu$ and is given by Eqs. $\eqref{result2}$ and $\eqref{Large}$ (dashed line in the Figure 9).

The case of $\mu=0$ is special at small $\mathscr{E}$ since the concentration of transmission in a small circle of energies does not suppress the Andreev scattering. The ABS in this case are concentrated in this small energy circle (see Golubov et al. (2004)) and depend on all phases. This is why $B_{12}$ does not drop to 0 but rather approaches a finite limit at $\mathscr{E}\to 0$. At $\mathscr{E}\gg\Delta$ $B_{12}$ still saturates at the value corresponding to the weak energy dependence case.

VIII Summary and Conclusions

In this Article, we address the topological properties of multi-terminal superconducting nanostructures. This involves Berry curvatures in the parametric space of the superconducting phases of the terminals and associated Chern numbers that manifest themselves in quantized transconductances Riwar et al. (2016).

The specifics of the superconducting nanostuctures is the presence of continuous spectrum along with the discrete one. The Berry curvature is readily defined for a discrete spectrum. Its generalization for a (partly) continuous spectrum is not straightforward, and is a problem of general interest. It has not been solved in Ref. Riwar et al. (2016).

We perform the calculation in imaginary time, and model the nanostructure with an energy-dependend scattering matrix. We have derived a general action of superconducting nanostructure with time-dependent phases, this is a separate advance. We expand the action near a point in the space of phases to compute the response function at finite frequency. We define the tensor quantity $B_{\alpha\beta}$ (Eq. 21) as a first term in the expansion of the response function at small frequency. This quantity would have been Berry curvature if the spectrum were entirely discrete.

We analyze the topological properties of the computed quantity. Like for Berry curvature, the topological charge associated with divergence of $B_{\alpha\beta}$ is concentrated in the singular points of 3d phase space where ABS cross zero energy — Weyl points. Unlike Berry curvature, the quantity $B_{\alpha\beta}$ has a non-topological contribution that is constant over the space of phases (Eq. 27). This in general adds a non-quantized part to ”Chern” numbers defined as integrals of $B_{\alpha\beta}$ over two superconducting phases, and to the corresponding transconductances. This contribution is determined by the scattering matrix at $\epsilon\to\infty$. It vanishes if the scattering matrix without superconducting phases is time-reversible and if the scattering matrix approaches isolation limit $S_{\epsilon}=1$ at large energies. For an energy-independent scattering matrix, the non-topological term is associated with the anti-symmetrized part of the conductance matrix of the structure in the normal state.

We consider in detail the case of weak energy dependence of the scattering matrix. We separate the contributions of the discrete and continuous spectrum, find them equally important and derive a compact relation for $B_{\alpha\beta}$ (Eq. 29).

We analyze in detail the Berry curvature in the vicinity of Weyl points. We have found a violation of topological protection of ”Chern” number in case of weak spin-orbit coupling. This, however, is rather trivially related to the transition between the ground states of different parity near the Weyl point and associated discontinuity of the wave functions. The topological protection is restored if one considers a ground state of a fixed parity.

We also investigate the properties of $B_{\alpha\beta}$ for the scattering matrices that essentially depend on energy at the energy scale $\simeq\Delta$.

IX Appendix A: derivation of the action

In this Appendix, we derive the effective action for a multi-terminal superconducting junction within the scattering approach. We follow the lines of Ref.Nazarov (2015).

In contrast to Ref.Nazarov (2015) we proceed in Matsubara formalism. Let us start with the formulation of a concrete microscopic model. Since the scattering formalism is universal, there is a great degree of arbitrariness in the choice of the model: all models that are characterized by the same scattering matrix will result in the same action. Properties of the scatterer are to be completely described by an S-matrix, the details of the model that describes the system are not important. So we choose the model in a way we find it convenient (see Fig. 10). We consider a system of independent 1-dimensional channels with pairwise opposite velocities and a linear spectrum. They are defined in the interval $-\infty<x<0$. The total number of channels is $2N$, number $N$ includes the spin doubling. Two channels in a pair with opposite velocities are coupled to the same superconducting reservoir: this is required to assure the time-reversibility of the model at this level. The coupling is a tunnel one, and the coupling strength is characterized by the dwell time scale $\tau$: at this time scale, an electron in a channel would tunnel to a reservoir. The tunneling results in an addition of self-energy to Green’s functions in the channels, which is proportional to the tunneling rate $1/\tau$ and to a matrix Green’s function $g$ characterizing a reservoir (see its concrete definition below). The channels defined in such a way model the electron states coming from and going to the reservoirs that are scattered at the nanoscructure. In the scattering region with a coordinate $y\in[0,l]$, there are $N$ spin-degenerate channels of the same velocity direction. At the boundary $y=0$ the electron amplitudes in the channels match those in the channels of positive velocity at $x=0$ (incoming states), while at $y=l$ the amplitudes match those in the channels with the negative velocity(outgoing states). As we will show, the S-matrix relates the amplitudes at $y=l$ and $y=0$.

To find the action for the nanostructure, we will compute its variation with respect to the variation of $g$. To this end, we require the values of the Green’s functions in the channels $x,x^{\prime}<0$ in close points $x\approx x^{\prime}$. We find the variation in three steps. At the first step, we express the Green’s functions at any $x$ in terms of the Green’s functions at $x\approx 0$. At the second step, we consider the scattering region that provides a boundary condition. With this, we relate these Green’s functions, and solve for them. This permits to find the variation and the action at the third step.

In the channels, we choose the basis in the following form

[TABLE]

where $u_{\pm},v_{\pm}$ are $N$ vectors in the space of the channels associated with the electron and hole amplitudes of the Bogolyubov wave function, and $\pm$ refers to the sign of the velocity in corresponding channels. In this basis, the equation for Green’s function reads

[TABLE]

where $v$ is the velocity that we can set the same for all the channels, $\epsilon$ is the Matsubara frequency, $\tau_{i}$ are Pauli matrices in Nambu space, and $\eta_{3}=\pm$ distinguishes channels with positive and negative velocities. The matrix $g$ is block-diagonal in the channel space. For a given reservoir, it is given by

[TABLE]

$g^{2}=1$, $\Delta$ being the superconducting order parameter in the corresponding reservoir.

We define a block structure

[TABLE]

We are only interested in the diagonal blocks $G_{1;2}$ since the off-diagonal blocks will not contribute to the variation of the action. We integrate the equation assuming $\epsilon\tau\ll 1$ for $G(x,x^{\prime})$ at $x<x^{\prime}$ we obtain

[TABLE]

[TABLE]

where we use special notations for the Green’s functions in the close points

[TABLE]

Since the solution for the Green’s function should not grow $x\to-\infty$, these Green’s functions should satisfy the following conditions

[TABLE]

[TABLE]

These matrices $G_{1;2}^{-}$ can be fixed if we consider the boundary conditions, that can be obtained by solving the equations for the Green’s functions in the the scattering region $y\in[0;l]$. To derive these condition, let us introduce the amplitude vectors $\Psi(y)=G(y,x)$, $X(y)=G(x,y)$ that have Nambu structure $\begin{pmatrix}u(y)\\ v(y)\end{pmatrix}$ and satisfy the equations

[TABLE]

[TABLE]

where $U(y,\epsilon)$ is the $N\times N$ matrix potential acting on electrons inside the scattering region and mixing different channels. The solution of the Eq.$\eqref{ampleq1}$ gives a linear relation on the amplitudes

[TABLE]

where we define the S-matrix for electrons and holes arranged in Nabmu structure

[TABLE]

[TABLE]

The electron scattering matrix is given by

[TABLE]

where $T_{y}$ implies the ordering of the $U(y)$ operators in the exponent according to the values of $y$ in the increasing order. We do not need to specify the energy dependence fo the S-matrix except for the general condition $S_{\epsilon}S_{-\epsilon}^{\dagger}=1$.

The relation on the amplitude $\eqref{bncxi}$ gives the relation between the diagonal and off-diagonal blocks of the Green’s function $\eqref{blocks}$ outside the scattering region but close to it $|x\epsilon/v|\ll 1,|x^{\prime}\epsilon/v|\ll 1$

[TABLE]

The solution of Eq. $\eqref{ampleq}$

[TABLE]

yields another relation between the blocks

[TABLE]

Combining Eq. $\eqref{32}$ and $\eqref{13}$ we obtain the required boundary condition that relates the diagonal sub-blocks

[TABLE]

Combining the equations $\eqref{boundc}$, $\eqref{eq1}$ and $\eqref{eq2}$, and the condition

[TABLE]

that follows directly from $\eqref{Gequation}$ we solve the complete linear system of the equations to obtain the follwing for the diagonal blocks of the general Green’s function $\eqref{blocks}$

[TABLE]

[TABLE]

Next, we employ the formula that expresses the action variation in terms of Green’s functions. We vary the reservoir Green’s function $g$ keeping normalization $g^{2}=1$, so that $\{g,\delta g\}=0$, then the variation of the action $L$ is

[TABLE]

where $\delta\Sigma=\frac{-i}{2\tau}\delta g$ is the variation of self-energy of electrons in channels and $G_{\rm Ch}(x,x)$ is their Green’s function at coinciding points. We note here that indeed only the diagonal blocks $G_{1;2}$ in Eq.$\eqref{blocks}$ contribute since $\Sigma$ is diagonal in this basis. The contribution from the channels corresponding to $G_{1}$ gives

[TABLE]

[TABLE]

The futher calculations is convenient to do in the basis that diagonalizes $g$. In this basis,

[TABLE]

[TABLE]

we find

[TABLE]

[TABLE]

where all the realtions are valid up to the first order in variations. The contribution from the outgoing channels reads

[TABLE]

[TABLE]

Summing both contributions, we obtain

[TABLE]

Hence

[TABLE]

This so-called block-determinant result for the action is similar to the one obtained previously Nazarov (2015) within the Keldysh formalism.

X Appendix B: derivation of the response function

In this Appendix, we present the details of the derivation of the Eq.$\eqref{general}$ and Eq.$\eqref{result1}$. We start with the action as given by Eq. $\eqref{lagrange}$. In order to derive the response function, we assume that the time-dependent deviation ($\delta\phi(\tau)$) from the stationary phase denoted as $\phi$ is small ($\delta\phi(\tau)\ll 2\pi$) so we can expand the action in Taylor series in $\delta\phi(\tau)$. We also note that in time representation the total phase operator is diagonal ($\phi_{\tau\tau^{\prime}}=\delta_{\tau\tau^{\prime}}\phi(\tau)$), which implies that the energy representation of $\phi$ reads

[TABLE]

We consider here the general case of the energy-dependent scattering matrix. The action from Eq.$\eqref{lagrange}$ reads

[TABLE]

$T$ implies the complete operator transposition that includes the reversing of the sign of energy. We remind the definition

[TABLE]

We ascribe the stationary part of the phases to an S-matrix $S_{\epsilon}\to S_{\epsilon}(\phi)$ and expand in small nonstationary deviation $\delta\phi(\tau)$.

[TABLE]

We introduce

[TABLE]

With this,

[TABLE]

We remind the definition of the matrix, that projects on the channels connected to a given terminal $\alpha$:

[TABLE]

where $a,b$ indices are in channels. With the help of this matrix the phase variation can be conveniently expressed as

[TABLE]

For simplicity of the notations, we denote the stationary phase derivatives $\partial_{\phi_{\alpha}}=\partial_{\alpha}$. With all this we consider the expansion of the S-matrix

[TABLE]

Let us we also note the identities for the derivatives with respect to the stationary phases:

[TABLE]

[TABLE]

the first term in the expansion $\eqref{expansion}$ vanishes since $\delta\phi_{\omega=0}=0$. The second term is

[TABLE]

The first term here does not depend on frequency and does not vanish in the limit $\omega\to 0$. The second term up to linear order in $\omega$ can be rewritten as

[TABLE]

The second term in the expansion $\eqref{expansion}$ reads

[TABLE]

where subscripts mean taking the function at the frequency $\epsilon_{1,2}:\quad\epsilon_{1}=\epsilon_{2}+\omega$ and we denoted $\epsilon_{2}=\epsilon$. Summing it with $\eqref{part1R}$ we get the general response function as in Eq. $\eqref{general}$.

To perform the adiabatic expansion in the small parameter $\omega/|\Delta|$ here we keep $\omega$ as an independent parameter. We will use the identities

[TABLE]

where we introduced ”classical” and ”quantum” S-matrices as

[TABLE]

With this, we rewrite the term

[TABLE]

Next, we expand the terms that are taken at $\epsilon_{1}=\epsilon_{2}+\omega$. Thery come from three factors here. The expansion of the first factor $Q^{-1}_{1}\simeq Q_{2}^{-1}+\omega\frac{\partial Q_{\epsilon}^{-1}}{\partial\epsilon}$ gives rise to

[TABLE]

The expansion of the product of the classical parts is symmetric with respect to $\alpha,\beta$, so it vanishes. The product of quantum parts vanishes in linear order in $\omega$. So we only need to consider quantum times classical and expand the quantum one

[TABLE]

it yields

[TABLE]

where the first doubling is due to the same contribution with $\alpha\leftrightarrow\beta$. Summing it with $\eqref{first}$ we obtain the total response function as given by $\eqref{result1}$

[TABLE]

Acknowledgements.

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement №694272).