Functional renormalization group study of the Anderson–Holstein model

We present a comprehensive study of the spectral and transport properties in the Anderson–Holstein model both in and out of equilibrium using the functional renormalization group (fRG). We show how the previously established machinery of Matsubara and Keldysh fRG can be extended to include the local phonon mode. Based on the analysis of spectral properties in equilibrium we identify different regimes depending on the strength of the electron–phonon interaction and the frequency of the phonon mode. We supplement these considerations with analytical results from the Kondo model. We also calculate the nonlinear differential conductance through the Anderson–Holstein quantum dot and find clear signatures of the presence of the phonon mode.


Introduction
The Anderson-Holstein model is widely used to describe electronic transport through individual molecules contacted between two leads: The lowest unoccupied molecular orbital is described as a single, spin-degenerate (in the absence of a magnetic field) level of a quantum dot with a Coulomb repulsion between electrons of opposite spin.In addition, the molecule can vibrate at some characteristic frequency ω 0 , and the electrons interact with the polarization field generated by this vibration.
Equilibrium properties, such as linear conductance and spectral density, of the Anderson-Holstein model have been studied perturbatively [1,2], and nonperturbatively using numerical renormalization group (NRG) techniques [3,4].In the limit of weak coupling to the leads, a Lang-Firsov transformation together with a generalized Schrieffer-Wolff transformation can be used to obtain an effective Kondo model, which allows for a traditional renormalization group treatment [5].
Using the Meir-Wingreen formula [6] it is possible to obtain a simple expression for the non-linear current-voltage characteristics of the model, under the assumption that the spectral function does not change when a bias voltage is applied.This means essentially that the molecule is kept well-grounded by a good contact with one lead whereas the other lead at a finite voltage is very weakly coupled to the molecule.This can be realized in a scanning tunneling microscope (STM) setup where the differential conductance through the STM tip probes the spectral density of the molecule.Due to the above assumption the many-body calculation can be performed in equilibrium as in [4,5].
Going beyond the STM-like setup the Anderson-Holstein model in a bias voltage driven non-equilibrium steady-state has been studied using real-time diagrammatics [7], rate equations [8,1,9], and slave-boson techniques with the non-crossing approximation [10].In these cases the double occupation of the dot has been forbidden by considering the limit of an infinitely strong Coulomb repulsion, suppressing all charge fluctuations.Some numerical progress in the case of finite Coulomb repulsion has also been made using the imaginary-time non-equilibrium formalism [11,12].
In addition, a related resonant level model with a local phonon mode has been studied using real-time diagrammatic Monte Carlo [13], variational Lang-Firsov transformation [14], numerical path integral approaches [15,16], NRG [17], and scattering states NRG [18].In this model also the time evolution towards the nonequilibrium steady state has been studied [19].The lack of spin excludes all Kondorelated effects, however.This model admits a simple description of the interplay between electron tunneling and phonon dynamics in terms of the ratio ω 0 /Γ, where Γ is the tunneling rate of the electrons, both limits ω 0 /Γ 1 and ω 0 /Γ 1 being theoretically well controlled [17].
An experimental realization of an Anderson-Holstein quantum dot can be for example a suspended carbon nanotube, where the principal vibration mode of the tube plays the role of a local phonon mode [20].The reported strength of the electron-phonon coupling in this case is relatively high -the ratio between the coupling constant and the phonon frequency, λ/ω 0 , reaching values in excess of 5 [20].Moreover, due to the length of the nanotubes the ratio ω 0 /Γ can be of the order of one, a regime which is theoretically challenging.Accessing this calls for a nonperturbative treatment.
In this article we study the spectral and transport properties in the Anderson-Holstein model using the functional renormalization group (FRG) [21].We show that in equilibrium FRG can be successfully used to approach the regime of intermediate electron-phonon coupling, resulting in a Kondo-like behavior.We also use the nonequilibrium Keldysh formulation of FRG [22,21] to study the non-linear electron transport through the Anderson-Holstein quantum dot and compare our results to first order perturbation theory.
This Article is structured in the following way: We introduce the model and express it in terms of the Keldysh Green functions, and briefly describe the FRG in section 2. We present the equilibrium properties, i.e., spectral densities and linear conductance in section 3, and non-equilibrium spectral densities and differential conductance in section 4. We summarize our findings and describe avenues for future research in section 5. Throughout this paper we work with natural units in which = k B = 1.

Hamiltonian and effective action
The system is described by the Anderson-Holstein Hamiltonian where The fermionic operators d σ and c krσ refer to electrons on the dot and in reservoir r = R, L, respectively, while the bosonic operator b refers to phonons with frequency ω 0 and coupling constant λ.Furthermore, eV G is the gate voltage and U the strength of the repulsive interaction between electrons on the dot.With our definition of , eV G = 0 corresponds to half filling of the dot.
Since we consider only equilibrium, as well as steady state non-equilibrium, the system becomes translationally invariant in time and it is sufficient to work in frequency space.We treat the reservoirs in the wide band limit with a constant density of states D r and assume that the coupling between the reservoir and the dot states does not depend on momentum and spin, After integrating out the reservoirs, held locally in equilibrium with chemical potential µ r and temperature T r , the noninteracting but reservoir dressed electron propagator on the dot reads and n r (ω) = 1 e (ω−µr)/Tr + 1 .
We present the derivation in terms of the non-equilibrium Keldysh Green functions.Analogous but somewhat simpler expressions hold for the equilibrium Matsubara Green functions [23].Finally, we integrate out the phonons to obtain an effective action for the interacting electrons on the dot: where the multi-indices 1 = (l 1 , ω 1 , σ 1 ) contain the Keldysh index §, frequency, and spin, and repeated indices are summed or integrated over: l σ dω 2π .The direct ‡ We note that sometimes Γ is defined in the literature with an additional prefactor of 2. § We use the convention where the 2 × 2 matrix Green function is given by antisymmetrized electron-electron interaction vertex is given by νee and the retarded phonon-mediated electron-electron interaction vertex (electronphonon-electron vertex) by with the phonon propagators Here, T p is the phonon temperature and η a positive infinitesimal, which also defines the width of the delta peaks in (20).In the numerical calculations we assign it a finite value which is small enough not to affect the results.We have also defined the bosonic frequencies where the second equalities hold due to frequency conservation.

Functional renormalization group
After integrating out exactly both the leads and the phonons we tackle the remaining problem of the interacting electrons on the dot with the functional renormalization group [21].For the non-equilibrium calculations we use the Keldysh formulation developed for the single-impurity Anderson model (SIAM) in [22].We use the hybridization to a structureless auxiliary lead, Λ, as the flow parameter and truncate the hierarchy of FRG equations at the second order.This includes all interactions at least to fourth order in λ as well as contributions beyond plain perturbation theory by virtue of the RG flow.We also use the approximate mixing of the channels and the static self-energy feedback described in sections VI B and VI D of [22], respectively.Since the Keldysh structure of the bare electron-phonon-electron vertex is compatible with the structure used in the flow equations of [22], only the initial conditions of the flow need to be adjusted to account for the phonon-mediated interaction.Compared to the flow without phonons, we need to add λ 2 2 D(∆) to the initial value of (a d Λ ) σσ (∆) and of (a d Λ ) σσ (∆), and to add λ 2 2 K(∆) to the initial value of (b d Λ ) σσ (∆) and of (b d Λ ) σσ (∆); for the definitions of the functions (a d Λ ) σσ (∆) and (b d Λ ) σσ (∆), see Appendix A of [22].The initial values of U d Λ and W d Λσ change according to (87b) and (87c) of [22].For the self-energy we obtain the initial conditions where the difference to (39c) in [22] results from the Hartree-Fock diagram with the interaction vertex νepe .The resulting system of differential equations is integrated from Λ = ∞ to Λ = 0 numerically.For the calculations of equilibrium quantities, namely effective mass (section 3.2) and linear conductance (section 3.3), we also use the Matsubara formulation of FRG.To include the local phonon mode we make analogous changes as described above to the initial conditions presented in [23].In addition, with Matsubara FRG we use frequency cut-off as the flow parameter to simplify the calculations.
Once both the self-energy and the vertex function are determined at the end of the FRG flow, we obtain the full interacting spectral function via and the current by [6] The (differential) conductance is then obtained by numerically differentiating the current with respect to the bias voltage.

First order perturbation theory
For U = 0 and small U epe /Γ we compare FRG results with a perturbation theory to first order in U epe .The two addends (Hartree and Fock) of νepe in (18) yield two contributions to the retarded self-energy in first order perturbation theory.The first addend is frequency independent, where nσ denotes the average occupation of the spin state σ in the noninteracting dot.We consider it reasonable to evaluate this addend self-consistently in order to maintain particle-hole symmetry.Then it satisfies where we have suppressed the spin index and used In the special case Consider this equation for Then it has a unique solution for U epe /Γ ≤ π/2.

In particular at eV
For U epe /Γ > π/2 and V G sufficiently close to µ the equation has three solutions and leads to a hysteresis of Σ R H (V G ). Therefore we use the self-consistent approach only for U epe /Γ ≤ π/2.The second addend in (18) yields in first order perturbation theory the contribution Here we choose the free propagation to include the self-consistent Hartree contribution determined above.For T L = T R = 0 the integral can be evaluated analytically, yielding where we have again suppressed the spin index and set The total retarded self-energy in first order perturbation theory is We note that Σ R F (ω) diverges logarithmically at ω = ±ω 0 + µ r = ±ω 0 ± eV /2.Consequently, the spectral density in first order perturbation theory vanishes at these frequencies.An exception is the particle-hole symmetric situation where the eV = 2ω 0 dip at zero energy is missing (cf.figure 4 and figure 6).This is due to a cancellation of these logarithmic divergencies.

Choice of the parameters
The Anderson-Holstein model has a large number of physical parameters which can be varied, thus making a thorough study of the model complicated.To keep things manageable, we restrict our study in this Article to the case of a vanishing magnetic field, and zero temperature, T L = T R = T p = 0.Moreover, to enhance the visibility of phonon-related effects, we will focus on the case where the electron-phonon mediated interaction dominates over the bare Coulomb interaction and therefore mainly set U = 0. We will also restrict our study to a left-right symmetric structure, Γ L = Γ R .Using the tunneling rate to the leads, Γ, as a characteristic energy scale we are left with the freedom to vary the strength of the electron-phonon coupling, λ/Γ, and the frequency of the phonon mode, ω 0 /Γ.Furthermore, it turns out (see section 3.1) that the physically most relevant combinations are λ/ω 0 and 2λ 2 /(ω 0 Γ) ≡ U epe /Γ.We will also study how the results change when the average occupation of the dot is changed by tuning the gate voltage eV G .

Identification of the parameter regimes
The spinless resonant level model with a local phonon mode is commonly characterized in terms of the ratio ω 0 /Γ [17]: In the adiabatic regime, ω 0 /Γ 1, the local phonon is too slow to effectively respond to the rapidly tunneling electrons and its presence has little effect on their motion.In the anti-adiabatic regime, ω 0 /Γ 1, an electron on the dot can lower its energy by forming a polaron with the phonon mode, the energy shift being E p = −2λ 2 /ω 0 = −U epe .In the limit λ/ω 0 1 the tunneling rate becomes strongly renormalized to Γ eff = Γe −(λ/ω 0 ) 2 Γ.For this reason it is expected that for strong electron-phonon coupling the crossover from anti-adiabatic to the adiabatic regime happens at ω 0 /Γ eff ≈ 1 instead of ω 0 /Γ ≈ 1.For the resonant level model this has been confirmed using NRG in [17].
In the Anderson-Holstein model the presence of the spin leads to richer physics than what can be simply described by an effective tunneling rate.However, the width of the central Kondo resonance in the spectral function is similarly an important quantity.The renormalization of the width can be described by Γ ren = Γ/m * with the effective mass, m * = 1 − ∂ ω Σ R σ ω=0 .To leading order in Γ, analytical estimates can be derived from an effective Kondo model [4].
In the limit U epe U , the Kondo coupling constant is given by [4] and the effective mass by i.e., the electron-phonon coupling reduces the sharpening of the Kondo resonance.This is especially intuitive in the limit ω 0 U , where the effective electron-electron interaction is given simply by U eff = U − U epe .
We are more interested in the limit U epe U .In this case the effective model is the anisotropic Kondo model with the coupling constants where For J ⊥ J , and for When λ ω 0 , U we find and For λ ω 0 we have J /J ⊥ ∝ exp[2(λ/ω 0 ) 2 ], and On the other hand, in the limit ω 0 λ U , J ⊥ D = J D = 8Γ/(πU epe ), and which is equal to the effective mass in the regular Anderson model with U replaced by U epe .
Let us now consider the case where U = 0.The effective anisotropic Kondo model is still valid when U epe /Γ 1.For the opposite limit U epe /Γ 1 and weak electronphonon coupling, λ/ω 0 1, the problem is perturbative and no Kondo-like exponential behavior can be expected.For λ/ω 0 1 we find that the effective mass behaves exponentially, but it is given by (46) instead of (45) as for λ/ω 0 1, U epe /Γ 1 (see section 3.2).
For λ/ω 0 1, ω 0 ≈ Γ, it follows that U epe /Γ 1.This lies in the regime of strong renormalization, regime (II) in figure 1, and is therefore similar to the extended antiadiabatic regime of the resonant level model [17].Decreasing ω 0 /Γ traces a vertical line towards the adiabatic regime (IV), with the crossover happening roughly at U epe /Γ ≈ 1 where ω 0 Γ.Regime (I) we identify with the negative-U , i.e., attractive SIAM with U = −U epe .Finally, in regime (III) the electron-phonon interaction is weak, and the model can be studied perturbatively.

Effective mass from FRG
The dependence of the effective mass on λ/ω 0 with ω 0 /Γ fixed, obtained numerically from the equilibrium Matsubara FRG, is shown in figure 2. A few data points from Keldysh FRG are included for comparison.Varying λ/ω 0 in this way traces a line across the four different regimes described above (cf.figure 1).
In the regime (I), the effective mass should be described by ( 46), indicated as a dashed line in figure 2(a).The numerically calculated effective mass for (λ/ω 0 ) 2 ∈ [0.02, 0.05] grows with a slightly smaller prefactor in the exponential function than the expected πω 0 /(4Γ), however.A similar discrepancy is seen in the FRG studies of the SIAM [23,22].For very large U epe /Γ, (λ/ω 0 ) 2 0.05, the effective mass bends away from the exponentially increasing trend.This effect is also known from the SIAM, signaling the limit of the validity of our FRG method.
In the regime (IV), the effective mass is well described by (46) as is evident from figure 2(d) for (λ/ω 0 ) 2 ∈ [4,50].As U epe /Γ approaches unity at (λ/ω 0 ) 2 = 125 the effective mass starts to increase rapidly.We are unable to reach far into regime (II), with m * described by (45), since the renormalized interaction vertex grows too large , again signaling the limit of the validity of our FRG method.This is even more clear from figure 2(b) in which ln(m * ) is plotted against (λ/ω 0 ) 4 instead of (λ/ω 0 ) 2 as in (a), (c), and (d).The behavior described by (45), indicated as a dashed line, is not observed even for the largest values of λ/ω 0 .Here, U epe /Γ reaches unity at (λ/ω 0 ) 4 = 25.In addition, as we approach regime (II) so that U epe /Γ 1, the effective mass from We use a criterion where the FRG result is deemed unreliable if the renormalized value of the vertex function grows to more than ten times its original value at any energy during the flow.
Keldysh FRG starts to deviate from that of Matsubara FRG.This can be clearly seen from the rightmost Keldysh data point in figure 2(c).
In the regime (III), c. f. figure 2(c), the effective mass from FRG and first order perturbation theory agree up to λ/ω 0 ≈ 0.3, or, equivalently, U epe /Γ ≈ 0.45.The linear conductance as a function of the gate voltage is shown in figure 3.In general the linear conductance exhibits a peak at a vanishing gate voltage, i.e., at the particle-hole symmetric point.In the case considered here, U epe U , the low-energy excitations are described by charge fluctuations instead of the spin fluctuations of the ordinary SIAM.As a result the roles of the charge and spin are essentially switched.Application of a gate voltage then has a similar effect as the application of a magnetic field in the ordinary SIAM, leading to a collapse of the conductance plateau.We observe numerically that the width of the resulting conductance peak decreases linearly with increasing U epe /Γ except when approaching the regime (II) of strong electronphonon coupling, where the width of the peak decreases significantly faster.In all cases G = 2e 2 /h for eV G = 0.

Linear conductance
Looking at the inset of figure 3, results from first order perturbation theory agree perfectly with the results from FRG for U epe /Γ = 0.2, whereas for U epe /Γ = 1.0 deviations start to appear.
It is important to note that the Keldysh formulation of FRG produces markedly different results than the equilibrium FRG when approaching the regime (II).This is another signature of reaching the validity limit of our truncated FRG equations.

Spectral density
To understand the nonlinear current-voltage characteristics of the model, let us first look at how the spectral density changes upon applying a bias voltage.
In general, features in the spectral density appear at various energies.For eV G = 0, the main resonance lies at ω = 0 with phonon steps at ω = ±ω 0 .The phonon steps are accompanied by shoulders at higher energies with a width that depends on U epe /Γ.At finite bias voltages these steps split to ω = ±ω 0 ± eV /2, each with an accompanying shoulder.As shown in section 2.3, the spectral density in first order perturbation theory goes to zero at these energies.The most pronounced effect of the higher order terms included in FRG is the smoothening of these features.
For U epe /Γ 1 (see figure 4) the height of the central resonance decreases as the bias voltage is increased, whereas the height of the shoulders increase.When eV = 2ω 0 , the steps at ω = ±(ω 0 − eV /2) merge as can be seen from the double peak structure at zero energy.For eV > 2ω 0 these features vanish completely, leaving only the kinks at ω = ±(ω 0 + eV /2) (barely visible in the main panel of figure 4).
For U epe /Γ 1 (see figure 5) the spectral weight is significantly redistributed: The height of the central resonance decreases quickly as the bias voltage is increased.The maxima of the shoulders approach ω = ±(ω 0 + eV /2) and they increase in height.The behavior of the steps at ω = ±ω 0 ± eV /2 is similar as in the case U epe /Γ 1.
We note that the validity of FRG when approaching regime (II) may be somewhat questionable in view of the discussion in section 3.2 and section 3.3.However, since there are no comparable results from any other method available in the literature, we have decided to present non-equilibrium results at the border towards this regime as well.We hope that this motivates further study of the topic by experts using different, suitable strong coupling methods.Moreover, for the ordinary SIAM the differential conductance obtained by FRG agrees with the differential conductance obtained by other numerical methods [24] despite deviations in the spectral function.This fact lends credence to the U epe /Γ = 5 results presented in this section.
Finally, for a finite gate voltage (see figure 6) the central resonance and the phonon shoulders shift approximately by eV G , while the phonon steps stay at the same energies.

Differential conductance
According to (28), the electric current at T = 0 and µ L,R = ±eV /2 is proportional to the integrated spectral weight in the transport window from ω = −eV /2 to eV /2.For the differential conductance follows The two addends in square brackets describe two mechanisms how the spectral weight in the transport window changes when V is increased.The first is due to the fact that the transport window gets larger when voltage is increased; it is always positive.The second stems from the fact that the spectral function changes its form when voltage is increased.
The differential conductance as a function of bias voltage is shown in figure 7 for G = 0. a vanishing bias voltage only the first addend contributes, and ρ σ (0) = 1/(πΓ) yields G(V = 0) = 2e 2 /h (with h = 2π in our units) as already noted for the linear conductance in section 3.3.When voltage is increased, both contributions lead to a decrease in differential conductance: the first addend decreases, since the value of the spectral function decreases with increasing |ω| from its maximal value at ω = 0; the second addend becomes negative, since spectral weight is moved to higher energies outside the transport window with increasing voltage.For small U epe this second effect is not very pronounced (compare figure 4).The width of the peak in the differential conductance at V = 0 is then essentially determined by the first addend and hence given by twice the width of the central peak in the spectral function.For larger U epe , the second effect is of considerable importance (compare figure 5).The width of the peak in the differential conductance at V = 0 is then significantly smaller than twice the width of the central peak in the spectral function.For large U epe and large bias voltage, the negative second contribution may even dominate over the first one.A negative differential conductance results, as seen in figure 7 for U epe /Γ = 5, ω 0 /Γ = 2.5 and eV /Γ around 5. At eV = ω 0 , the phonon steps of the spectral function situated at ω = ±(ω 0 − eV /2) enter the transport window.This leads to the kinks of the differential conductance at eV = ω 0 visible in figure 7.
Let us for completeness look at the effect of the gate voltage on the differential conductance.This is shown in figure 8.The interplay of the gate voltage and bias voltage induced collapse of the main resonance is complicated, but the low-bias behavior is reminiscent of the differential conductance in the SIAM at finite magnetic field [25,22]: The peak in the differential conductance splits as the gate voltage is increased.The kink at eV = ω 0 appears also at finite gate voltages, but there is an additional step close to eV = 2ω 0 , which can be related to the merge and subsequent vanishing of the features in the spectral density at ω = ±(ω 0 − eV /2).
In first order perturbation theory the kinks at eV = ω 0 are far more pronounced due to the fact that the first term in (47) vanishes, but the steps at eV = 2ω 0 are absent.The latter is then clearly a higher order effect.Overall, FRG and perturbation theory agree over the whole range of bias voltages as long as U epe /Γ 1.

Conclusions
We have studied the spectral and transport properties of the Anderson-Holstein model both in and out of equilibrium using the functional renormalization group.We have identified the parameter regimes in which our approach is reliable by studying the renormalization of the quantum dot level width and comparing to analytical estimates.We have found clear signatures of the phonon mode in the transport properties: The width of the linear conductance peak as a function of gate voltage is decreased by increasing the coupling strength to the local phonon mode.The differential conductance at finite bias voltages exhibits a kink when the bias voltage equals the phonon frequency, and, away from the particle-hole symmetric point, also when the bias voltage equals twice the phonon frequency.The thermoelectric properties of the SIAM with an attractive interaction between the dot electrons have been predicted in an NRG study to differ significantly from those of the model with a repulsive interaction [26].Since a sufficiently strong interaction with a local phonon mode gives rise to an overall effective attractive interaction between the electrons, it is a natural extension to our work to next study the thermoelectric properties of the Anderson-Holstein model.

Figure 1 .
Figure 1.Illustration of the division into four different parameter regimes: attractive single-impurity Anderson model (I), extended anti-adiabatic (II), perturbative (III), and adiabatic (IV).The colored lines and data points correspond to those of figure 2.

Figure 3 .
Figure 3. Linear conductance through the Anderson-Holstein quantum dot as a function of the gate voltage for different values of U epe /Γ and λ/ω 0 .The lines fall roughly into three groups based on the value of U epe /Γ.Solid lines are obtained with the equiliubrium Matsubara formalism and dashed lines with the non-equilibrium Keldysh formalism.Inset shows a comparison to first order perturbation theory (dashdotted lines).

Figure 7 .
Figure 7. Differential conductance through the Anderson-Holstein quantum dot as a function of bias voltage for V G = 0, and different values of U epe /Γ and ω 0 /Γ, calculated with the Keldysh FRG.Dash-dotted lines are obtained with first order perturbation theory.Inset shows a detailed view of the differential conductance near eV = ω 0 for ω 0 /Γ = 0.1.

Figure 8 .
Figure 8. Differential conductance through the Anderson-Holstein quantum dot as a function of bias voltage for different values of gate voltage and U epe /Γ = 1, ω 0 /Γ = 2.5, calculated with the Keldysh FRG.Inset shows a comparison to first order perturbation theory (dash-dotted lines).