{'fraction_non_alphanumeric': 0.044141770641745785, 'fraction_numerical': 0.017746184818809962, 'mean_word_length': 4.3383309008889475, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Full-counting statistics for molecular junctions: Fluctuation theorem and singularities

We study the full-counting statistics of charges transmitted through a single-level quantum dot weakly coupled to a local Einstein phonon which causes fluctuations in the dot energy. An analytic expression for the cumulant generating function, accurate up to second order in the electron-phonon coupling and valid for finite voltages and temperatures, is obtained in the extended wide-band limit. The result accounts for nonequilibrium phonon distributions induced by the source-drain bias voltage, and concomitantly satisfies the fluctuation theorem. Extending the counting field to the complex plane, we investigate the locations of possible singularities of the cumulant generating function, and exploit them to identify regimes in which the electron transfer is affected differently by the coupling to the phonons. Within a large-deviation analysis, we find a kink in the probability distribution, analogous to a first-order phase transition in thermodynamics, which would be a unique hallmark of the electron-phonon correlations. This kink reflects the fact that although inelastic scattering by the phonons once the voltage exceeds their frequency can scatter electrons opposite to the bias, this will never generate current flowing against the bias at zero temperature, in accordance with the fluctuation theorem.


I. INTRODUCTION
It has been revealed in electric transport measurements that the electron-phonon interaction induces unique features in the nonequilibrium current through molecular junctions 1,2 and atomic wires. [3][4][5] In particular, when the phonon energy is small compared to the resonance width on the junction, the inelastic phonon scattering increases/decreases the current for small/large transmission probabilities as the source-drain bias voltage V exceeds the local phonon frequency ω 0 . 1,3-5 Quite a number of theoretical microscopic models, e.g. Ref. 6,7, including density-functional theories, 8 have been devoted to such junctions (see Ref. 9 and references therein). They revealed that a simplified model, 10 of a single-level quantum dot coupled to a local Einstein phonon mode causing fluctuations of the dot energy level, 11 seems to suffice to capture this ubiquitous feature. Recently, the current noise of an atomic wire has been measured, 5 and it was observed that the electron-phonon interaction can enhance or reduce the noise, depending on the value of the transmission probability. Based on the theory of Avriller and Levy Yeyati, 12 the negative correction is understood as resulting from the anti-bunching of two electrons: An electron cannot be inelastically scattered by a phonon when the final state is already occupied by another electron.
Avriller and Levy Yeyati considerations follow from the theory of full-counting statistics (FCS), 13,14 which is most convenient for analyzing nonequilibrium electric transport. Indeed, considerable effort has been invested in recent years in exploiting FCS to study various aspects of nonequilibrium quantum transport [e.g. Refs. 12-28 and references therein]. FCS refers to the probability distribution P τ (q) of the charge q to be transmitted through a quantum conductor during a certain measurement time τ at out-of-equilibrium conditions (we set e = = 1). The effect of coupling to a vibrational mode on the electric transport has also been investigated in this context, both for a weak 12,21-25 and a strong 26-28 electron-phonon coupling.
In general, it is a rather formidable task to calculate the FCS of interacting electrons. For this reason, most of the investigators have taken advantage of the Keldysh field-theory technique. There, the characteristic function (CF), or the scaled cumulant generating-function (CGF) 29 pertaining to the steady state, can be formally written as a 'partition function' or a 'free energy', respectively, defined along the Keldysh timecontour. The λ appearing in these formulae is termed the counting field or the counting parameter. However, even if one calculates the CGF successfully, one still needs to find a way to characterize the electronic correlations in it.
One promising approach would be to utilize the location distribution of the zeros of the CF, or equivalently that of the singularities of the CGF, by allowing the counting field λ in Eq. (1) to attain complex values, [30][31][32][33] similarly to the Yang-Lee theory of phase transitions. 34 This idea is based on the recent observation that, upon transforming λ into u, the singularities of the CGF of noninteracting electrons transported between two terminals are all on the negative real axis of the u−plane. 30,35,36 It suggests that singularities off the negative real axis would characterize electronic correlations. From this aspect, molecular junctions are rather advantageous since second-order perturbation theory in the electron-phonon coupling would capture most relevant features of the electron-phonon correlations in them, allowing for obtaining the location distribution of the singularities. Another recent ingredient is the fluctuation theorem (FT). [37][38][39][40][41][42][43][44][45][46][47][48][49] The FT is a consequence of micro-reversibility and can be understood as a microscopic extension of the second law of thermodynamics. Despite its simple appearance, a detailed-balance like relation [see Eq. (6) below], the FT reproduces the linear-response results, i.e. it ensures the fluctuation-dissipation theorem and Onsager's reciprocal relations close to equilibrium, [38][39][40][41][42][43] while conveying invaluable information at nonequilibrium conditions. For molecular junctions, the FT has been addressed using the master-equation approach of FCS for incoherent electron transport. 26,27 The FT is considered to be a basic symmetry, such as gauge invariance, which the CF should fulfill.
In the present paper we investigate the FCS of electrons coupled to phonons under out-of-equilibrium conditions. The quantity to be calculated and analyzed is the cumulant generating-function. 12,[21][22][23][24][25] Employing the Luttinger-Ward functional, 50-52 , we obtain its full analytic expression, or equivalently all cumulants, in the wide-band limit, treating systematically the nonequilibrium phonon distribution via a self-consistent condition accurate up to second order in the electron-phonon coupling. There are other attempts in the literature 23,25 to account for the nonequilibrium phonon effect avoiding the Luttinger-Ward functional; these produce contradicting results; our first two cumulants agree with those of Ref. 23, which has not gone beyond the second cumulant. (Our results disagree with those of Ref. 25.) The structure of the paper is as follows. We begin in Sec. II with brief general explanations of the FT, the large-deviation analysis and singularities of the CGF. We then introduce in Sec. III the model Hamiltonian and present analytical results for the CGF and detailed explanations of the calculations and the approximations involved. A self-consistent calculation based on the nonequilibrium Luttinger-Ward functional, 50-52 is relegated to Appendices A and B. In Sec. IV we analyze the singularities of the CGF and demonstrate the prob-ability distribution within the large-deviation analysis. Section V summarizes our results. Technical details of the calculations are given in Appendices C and D.

II. THE FLUCTUATION THEOREM AND SINGULARITIES
The definition of the probability distribution P τ (q) in quantum systems requires special care. 13,14 Full-counting statistics theory begins with by introducing the CF [see Eq. (1) and Eq. (14) below for the definition employed in the Keldysh technique] and then defining the quasiprobability distribution by the inverse Fourier transform of the CF, More details are given in Sec. III. The n-th cumulant, in steady state, is given by the n-th derivative of the CGF, Eq. (2), For example, the average current is the first cumulant, and the current noise is the second one. The FT in the context of quantum electric transport relates the probability of the current to flow oppositely to the bias voltage because of thermal agitations, P τ (−q), [38][39][40][41][42][43][44][45][46][47][48] to the distribution P τ (q), where β is the inverse temperature. The FT can be equivalently written in terms of the CGF F , Eq. (2), as The relation (7) restricts the possible locations of the singularities of the CGF in the λ−plane. As an example, we depict in Fig. 1 (a) the branch cuts corresponding to the continuous singularities of the CGF pertaining to two-terminal transport of noninteracting electrons, Eq. (25) below. This CGF is 2π−periodic along the real axis of λ, which guarantees integer values of charge. 13 The branch cuts, depicted by thick lines, are at Re λ = (2n − 1)π, where n is an integer. The FT ensures that the branch cuts are symmetrically distributed around λ = iβV /2. (In Fig. 1 (a), the upper left thick line is identical to the lower right thick line, etc.) The 2π−periodicity is removed by the conformal transformation Eq. (3). Then the branch cuts are on the negative real axis of the u-plane [ Fig. 1 (b)]. 30,31,35,36 The steady-state probability distribution, beyond the central-limit theorem, is derived within the theory of large deviations. 29 At steady state, realized in the τ → ∞ limit, we scale q = Iτ and Z τ ≈ e τ F . Then the integral of the inverse Fourier transform Eq. (4) can be estimated by the saddle-point approximation 53 and the result is written with the rate function 29 I as P τ ≈ e −τ I . Since P τ is real and positive, the saddle point is expected to reside on the imaginary axis of complex λ−plane. Then the rate function is written in the form of a Legendre-Fenchel transform, Here λ is a purely imaginary number. In most cases, the CGF is real, i.e. the imaginary part of the exponent of the Fourier integral (4), ln Z τ (λ) − iλq, is zero on the imaginary axis of λ−plane. Then the imaginary axis is expected to be the steepest contour of the integral (4). 53 There are few exceptions where singularities of the CGF are on the imaginary axis 17,18 as we will also find below. The relation between the CGF and the rate function is analogous to that of thermodynamic potentials. It suggests that singularities on the imaginary axis of the λ−plane would also cause characteristic features in the rate function.

III. MODEL HAMILTONIAN AND CUMULANT-GENERATING FUNCTION
Our explicit calculations are carried out for a simple model [ Fig. 2], a single-level quantum dot coupled to a local Einstein phonon, which induces fluctuations in that level energy. [9][10][11]54 The model Hamiltonian is where the ± subscript is due to the presence of counting fields (see below). The molecular bridge or the atomic wire is described by the Hamiltonian H mol , in which c 0 (c † 0 ) destroys (creates) an electron on the localized level representing the molecule, of energy ǫ 0 , b and b † are the destruction and creation operators of the vibrations to which the electron is coupled while residing on the dot, and γ is the strength of the electron-phonon coupling. The vibrational modes obey the harmonic Hamiltonian where ω 0 is the frequency of the Einstein phonon. The leads are represented by free electron gases, of creation and destruction operators c † rk and c rk , and eigen energies ǫ rk , Finally there is the tunneling Hamiltonian, coupling the leads to the molecule. This part of the Hamiltonian is augmented by the counting fields. Those appear as phase factors on the operators c rk and c † rk , where J L and J R are the tunneling amplitudes between the left and the right lead, and the dot. In terms of the tunneling Hamiltonian Eq. (13), the characteristic function is 14 where K denotes the Keldysh contour, which runs from t = −τ /2 to t = τ /2 on the upper branch and returns to t = −τ /2 on the lower one (see Fig. 3) and T K is the timeordering operator along that contour. The subscript I indicates time dependence in the interaction picture. The ± notation indicates the branch of the Keldysh contour on which the tunneling Hamiltonian is effective, + for the upper branch in Fig. 3 and − for the lower one. In the long-time limit, the scaled CGF depends only on the the difference of the two counting fields As λ L +λ R counts the number of electrons flowing into the dot, the fact that the CGF depends solely on λ implies current conservation.

FIG. 3: The Keldysh contour.
Since in our case the electron-phonon coupling γ is weak, we carry out the calculation up to second-order accuracy, O(γ 2 ). There is a subtle point in this expansion. A naive second-order perturbation theory is not capable of producing the correct nonequilibrium phonon distribution. 6,8,[21][22][23]55,56 We therefore have to perform a re-summation of infinite diagrams by adopting the linked cluster expansion, see e.g. Ref. 19, or a more advanced method, the nonequilibrium Luttinger-Ward functional, Φ. [50][51][52]57 The first approach had been adopted in Refs. 24 and 25, while the second had been employed implicitly by Gogolin and Komnik,16 and hence by later studies 12,21-23 based on Ref. 16. Although the second scheme is physically transparent since it relies on the self energy of the electron Green function, the price to pay is that a self-consistent calculation is required in order to ensure various conservation laws. Here we follow the second approach based explicitly on the Luttinger-Ward functional, which enables us to perform the self-consistent calculations of the CGF in a transparent manner. Details, are given in Appendices A and B.
In terms of the Luttinger-Ward potential, the cumulant generating-function is where F 0 is the CGF pertaining to noninteracting electrons, and given in terms of the Keldysh Green function g λ [see Eq. (19)] of the electronic part of the Hamiltonian Eq. (9), with the trace being performed over the 2 × 2 Keldysh space. The CGF satisfies the normalization condition F 0 (0) = 0. The effect of the electron-phonon interaction is included in the scaled Luttinger-Ward potentialΦ (2) , which consists of linked diagrams up to O(γ 2 ) accuracy, A. The electronic part The inverse of the Keldysh Green function g λ reads where Γ L and Γ R are the partial widths of the localized level induced by the coupling with the leads, Here ν L and ν R are the densities of states in the left and right leads. Each of the leads is specified by its chemical potential µ r , such that µ L − µ R = V , and their electron/hole distribution is accordingly given by Our calculation is confined to the extended wideband limit, [21][22][23] in which either the level broadening Γ = r=L,R Γ r is larger than the all other energy scales except for the dot level |ǫ 0 | or the dot level |ǫ 0 | is larger than all other energy scales, i. e. or In both cases, we can neglect the energy dependence of the dot density of states. The normalized dot density of states is then replaced by its value at the Fermi energy [see Eq. (C3)] The calculation detailed in Appendix C yields 13 where Here T is the transmission of the localized level, written in terms of the normalized density of states Eq. (24) and the transmission probability at resonance The counting field dependence is only through X λ . This result obeys the FT since when λ → −λ + iβV , X λ is unchanged. This CGF possesses continuous lines of singularities for X λ ∈ [1, ∞). As discussed in Sec. II, in the u−plane these become branch cuts (−∞, u − ] and [u + , 0), where the branch points are at u ± = e −βV /2 (x ± √ The zeroth-order CGF (25) takes a particular simple form for a symmetric bridge Γ L = Γ R at resonance ǫ 0 = 0, for which the transmission is perfect T = 1 [see Eqs. (27), (28), and (24)] which describes Gaussian thermal fluctuations. On the other hand, when the bridge is extremely askew, Γ R ≪ Γ L or |ǫ 0 | ≫ Γ and thus T ≪ 1 one may expand the CGF to obtain which is the bi-directional Poisson form.

B. The phonon-induced part
The nonequilibrium Luttinger-Ward functional Φ in Eq. (18) results from the coupling of the charge carriers to the vibrational modes. We expand it diagrammatically in powers of the small parameter (see Appendix B for details) The relevant diagrams are shown in Fig This approximation accounts for the relaxation of the phonon mode by the particle-hole excitations in the electrodes. The RPA is expected to provide results accurate up to second order in γ, The dressed phonon Green function D λ in Eq. (32) is given by [see Appendix D] In Eq. (34) there appears the Keldysh particle-hole propagator Π, whose lesser/greater components are given by It is convenient to express these propagators in the forms whereΠ ±∓ rr ′ λ describe the particle-hole excitations created in the r and r ′ leads. Explicitly, where s r = ±1 for r = L/R and α rr ′ = 4Γ r Γ r ′ /Γ 2 . Here The lesser components are obtained from the symmetry relations For λ = 0, Eq. (37) reduces to the well-known form Having determined the dressed phonon Green function Eq. (34), we now use it to obtainΦ RPA , Eq. (32). By using Eqs. (D10) and (D16), the determinant of the dressed phonon Green function is found to be where the retarded component Π R is expressed by the lesser and greater components upon exploiting the Kramers-Kronig relation, Eq. (D17). The counting-field dependent phonon life time broadening, √ A λ is a crucial ingredient in our analysis. It is given by for the condition (22) and zero for the condition (23).
Collecting the results, we find that the derivative of Eq (32) is The ω−integrations is carried out to leading order in g. The integrand has four poles, which within that accuracy are given by where we have used A λ (ω) = A λ (−ω). The real part of the retarded component in Eq. (45) also shifts the argument of A, which can be neglected since it is already proportional to the small parameter g, Eq. (31). The locations of poles of the integrand in Eq. (44) are complex functions of λ. Consider them first in the limit λ → 0. Since A 0 (ω 0 ) = (2gρ 2 0 ω 0 ) 2 is a positive real number, then at small enough λ, the real part of the squire root satisfies Re A λ (ω 0 ) > 0 and thus the integration yields up to O(g). It follows that which satisfies the normalization conditionΦ (2) (0) = 0. Away from the origin λ = 0, Re A λ (ω 0 ) may be negative. When this is the case we obtain Hence a branch cut appears at λ satisfying Re A λ (ω 0 ) = 0.
In order to carry out a large-deviation analysis it suffices to consider the imaginary axis of the λ−plane, since there the rate function (8) is real. On the imaginary axis, A λ (ω 0 ) is a real function. When A λ (ω 0 ) > 0 the imaginary part ofΦ (2) is zero and thus the imaginary axis can serve as the steepest contour (the steepest ascent path) 53 of the integral Eq. (4). When A λ (ω 0 ) ≤ 0, there appears a branch cut on the imaginary axis. One might have thought that the branch cut would be detrimental to the saddle-point approximation of the integral in Eq. (4). Here we point out that a complex integration along a path encircling the branch cut oscillates rapidly in the τ → ∞ limit and thus would be averaged out. Therefore we will neglect the contribution from the branch cut.
Equations (42), (47), and (48), supplemented by Eqs. (37) and (43), are the main results of this paper. Closed expressions, but confined to the first two cumulants, have been derived in Ref. 23. In contrast, we obtain analytic expressions for the entire CGF, which enable us to examine its singularities and to fully analyze the rate function itself, as will be detailed in the Sec. IV.
To conclude this section we verify that our results obey the FT. Exploiting the extended detailed-balance relation, and the relations we find Therefore Eqs. (47) and (48) satisfy the FT,

IV. RESULTS AND DISCUSSION
In the following we confine ourselves to a symmetric junction, α = 1 (Γ L = Γ R ), and thus the normalized density of states on the dot at the Fermi level dominates the transmission probability, T = ρ 0 [see Eq. (27)]. The numerical results are all obtained for the electron-phonon coupling constant, Eq. (31), g = 0.1, unless otherwise specified.
A. Average current and noise Figure 5 (a) shows the source-drain bias voltage dependence of the current I at a finite temperature βω 0 = 10 (solid lines) and at zero temperature (dashed lines) for a perfect, T = 1, and a relatively weak, T = 0.5, transmission probabilities. At perfect transmission the current is suppressed above the threshold |V | > ω 0 , because electrons are inelastically backscattered by phonons. When the transmission is weak, T = 0.5, the current is slightly enhanced above the threshold. These results are consistent with previous ones. 21,22 A finite temperature tends to smear the kink structure of the T = 1 curve; it affects far less the average current at weak transmission, T = 0.5, where the solid and dashed lines almost overlap. Figure 5 (b) depicts the current noise I 2 . At perfect transmission the noise is absent below the threshold |V | < ω 0 at zero temperature. Thermal fluctuations which arise at finite temperatures induce additional noise below the threshold. Although the current is suppressed above the threshold |V | = ω 0 by the inelastic phonon scattering, the noise is significantly enhanced. This indicates that inelastic phonon scattering broadens the probability distribution of the current. In the case of a weak transmission, T = 0.5, the temperature effect is less dramatic-the noise is simply enhanced.

B. Singularities and the rate function
At zero temperature, it is possible to obtain a simpler form for the scaled CGF, which is useful for finding its singularities. For positive voltages V ≥ 0, the explicit form of the electronic part of the CGF, Eq. (25), reads where u is defined in Eq. (3). The function A λ which determines the phonon part is given by for 0 < V < ω 0 and for V ≥ ω 0 . Below, we investigate the analytic properties of the phonon-induced part of the CGF.

Elastic phonon scattering
As at zero temperature phonons cannot be excited when 0 < V < ω 0 , and electron transport at such voltages is hence affected only by elastic phonon scattering. Figure 6 (a) shows schematically the square-root branch cut ofΦ (2) , Eqs. (47) and (48). We find that the u−plane is separated into two by the the brach cut. This branch cut intersects the real axis at u 0 = 1 − 1/T [u 0 is indicated by empty dots in Fig. 6 (a)] and u 1 = (T − 1)(V + ω 0 )/[V (T −1)+T ω 0 ]. Around these points, A λ can be expanded as and Upon sweeping the transmission T from 0 to 1, u 0 increases from −∞ to 0, while u 1 increases but at T C = V /(V + ω 0 ) jumps from +∞ to −∞. Accordingly, we may define two regimes, or phases, 32 I (T < T C ) and II (T > T C ) as indicated in Fig. 6 (a). This classification roughly captures the behavior of the average current and the current noise. Figures 6 (b) and (c) depict the corrections induced by the electron-phonon interaction in the average current and in the current noise, respectively, We find that the electron-phonon interaction always increases the average current under the conditions adopted here [ Fig. 6 (b)]. On the other hand, the noise can be either enhanced or suppressed, depending on which regime the transmission is in [Fig. 6 (c)].
In regime I, one of the intersection points, u 1 , is on the positive real axis, outside the unit circle |u| = 1 [the left panel of Fig. 6 (a)]. Therefore, in the λ−plane, there is a nonanalytic point on the positive iλ−axis, which induces a weak non-convexity of the CGF as shown in Fig. 7 (a).
[The non-convex region is indicated by an arrow there. For comparison, we also plot the g = 0 case (the dotted line) for which the CGF is convex.] Figure 7  Then, similar to the way a firstorder phase transition manifests itself in thermodynamics, a kink appears in the rate function. We note that the location of the kink does not coincide with that of the peak; the peak of the rate function is at I = I where λ = 0 and I = 0. The physical consequence is that the elastic phonon scattering broadens the distribution by enhancing the probability of currents larger than the average value. It is important to note that the kink is a consequence of the non-convexity of the CGF on the positive real axis of the u−plane, a feature which is absent in the noninteracting-electrons case.
In regime II the two intersection points are on the negative real axis of the u−plane [the right panel of Fig. 6 (a)]. The CGF and the rate function pertaining to this case are plotted in Figs. 7 (c) and (d). The CGF is convex and the corresponding rate function is concave. In this regime the peak position is shifted from that for g = 0 [dashed line]. The elastic scattering by the phonons can either broaden or shrink the width of the rate function depending on the transmission probability as we deduce from Fig. 6 (c).

Inelastic phonon scattering
Above threshold V ≥ ω 0 , inelastic phonon scattering becomes possible. The analytic properties of the scaled CGF in this regime depend on the bias voltage. The branch cuts in the u−plane are schematically shown in Fig. 8 (a). Three branch points, u 1 , u 2 , and u 3 can be obtained by searching for the roots of A λ (u) = 0 [the filled dots in Fig. 8 (a)]. In addition there is another point, u 0 = 1−1/T [the empty dots in Fig. 8 (a)], around which A λ can be expanded The branch point u 1 is always on the negative real axis such that u 1 ≤ u 0 . From the positions of the other two branch points, u 2 and u 3 , we identify three regimes, see Fig. 8 (a). In regime I the two branch points are on the positive real axis outside the unit circle. As T increases u 2 and u 3 approach one another and meet at T = T − . Then in regime II, the two branch points are located symmetrically off the real axis. Upon further increasing T , u 2 and u 3 move in the complex plane and at T = T + they meet on the real axis again. In regime III, the two branch points are on the positive real axis inside the unit circle. In Figs. 9 (a) and (b) we plot the CGF and the rate function pertaining to regime I. The overall tendencies are similar to those found in regime I below threshold, 0 < V < ω 0 . In the shaded area of Fig. 9 (a) the CGF is nonanalytic and non-convex. This gives rise to a stronger kink structure in the rate function [ Fig. 9 (b)]. Figures  9 (c) and (d) present the CGF and the rate function in regime II. The overall tendency is again similar to those of regime II below threshold. Because of the nonanalyticities off the real axis, the statistics would not be reduced to that of noninteracting electrons.
We note that our identification of the three regimes roughly captures the behavior of the corrections to the current and the noise induced by the electron-phonon interaction [Figs. 8 (b) and (c)], which oscillate as a function of the transmission probability.  The bias voltage is V /ω0 = 1.5. The boundaries between I and II and between II and III are at T− ≈ 0.551 and T+ ≈ 0.908, respectively. Regime II roughly corresponds to the negativenoise correction region.

C. Fluctuation theorem
The analysis of the CGF and the rate function in regime III [see Fig. 8] is rather subtle. At zero temperature, there develops a non-convex region in the CGF, and the origin λ = 0 might enter it. When this happens, second-order perturbation fails since the rate function does not satisfy the relation I(I = I ) = 0 [this relation is related to the normalization condition F (λ = 0) = 0]. We therefore study regime III at finite temperatures, taking as an example βω 0 = 10, and show that the FT is crucial for obtaining a physically-reasonable result. Note that when the symmetry (7) holds, the FT (6) is also preserved within the large-deviation analysis, 58 (60) Figure 10 exhibits the CGF and the rate function at perfect transmission. For comparison, we plot the corresponding curves for noninteracting electrons, Eq. (29). As we have already noted when discussing the current noise, Fig. 5 (b), the width of the rate function is enhanced by inelastic phonon scattering [ Fig. 10 (b)]. The CGF obeys the FT, Eq. (7), and the curves are symmetric around the dot-dashed vertical line at iλ = −βV /2 [ Fig. 10 (a)]. The peak of the probability distribution is shifted in the negative direction and the probability to find large current fluctuations is suppressed as compared with the noninteracting case [ Fig. 10 (b)]. In the shaded area of Fig. 10 (a), the CGF is non-analytic and non-convex. Correspondingly, the rate function has a non-differentiable point at I = 0, see Fig. 10 (b). As a result, although the probability to observe currents smaller than the average value is enhanced by the inelastic phonon scattering, the probability to find negative currents I < 0 is strongly suppressed. This is consistent with the FT (6), which states that although thermal agitations generate current flowing in the opposite direction to the source-drain bias, that probability is exponentially suppressed at low temperatures. Note that previous studies report on a finite current flowing oppositely to the bias at zero temperature [see e.g. Eq. (13) in Ref. 12] in disagreement with the FT, although it may be quantitatively small. This can be easily seen by calculating the probability distribution of the transmitted charge q using the CGF as given by Eq. (13) of Ref. 12 and the inverse Fourier transform Eq. (4), where the parametersq 0 andq 1 are defined in Ref. 12. In the limit of zero temperature β → ∞ this probability distribution remains finite at q <q 0 including negative the q regime, which violates the FT (6) .

D. Discussion
Recently, Kumar et al 5 have explored experimentally the possibility to identify different regimes, as implied by the sign of the correction to the noise induced by the coupling with the phonons, 12,24 as depicted in Fig. 8  (c). Above the threshold, our classification predicts three regimes, similarly to Ref. 5 It is worthwhile to expound upon this point. Previous studies 5,12,23,24 have considered the changes of slope of the nth cumulant at threshold, as guidelines for the classifications of the various regimes. If we adopt this scheme, we find [see Eqs. (47), (54) and (55)] which upon expanding becomes We see that the slope difference of the n = 1 cumulant changes its sign at T = 0.635, that of n = 2 at T = 0.434, 0.816, etc. In the n = 2 case, the result is compatible with that of Ref. 23. For n = 3, we obtain three zeros (at T = 0.449, 0.576, 0.885) similarly to Refs. 12,24, though the positions are different. In general, the coefficient of λ n in the expansion is an n + 2thdegree polynomial function of T , yielding oscillations as a function of T , which is one example of the universal oscillations. 59,62 From the higher cumulants, we obtain more zeros, but it is unclear what useful information can be extracted of them. In contrast, the topology of the singularity distribution is distinct and, in our opinion, provides a better way to classify regimes where electron transport is affected differently by the coupling with the phonons. However, as a direct connection between the singularity distribution and the phonon-induced noise seems to be absent at the moment, it is not surprising that we obtain just poor quantitative agreements between our critical points and the positions of the zeros for the phonon-induced noise in Figs. 6 (c) and 8 (c). Another comment related to previous studies has to do with the high-bias limit. Urban et. al. 25 have suggested that in that limit, i.e., for ω 0 ≪ V ≪ Γ, the correction induced by the electron-phonon coupling to the cumulants scales as I n+1 ph / I n ph ∼ V /ω 0 . We can use our result Eq. (55) to obtain in this regime the expansion which implies a stronger scaling, Our result thus extends that of Ref. 23, where I 2 ph ∝ V 4 was reported. We note that our phonon-induced part of the CGF, Eq. (32), is already different form that of Urban et. al.,Eq. (11) in Ref. 25.
A comment on the validity of second-order perturbation theory is called for. Flindt et. al. 61 have analyzed the FCS of sequential transport through a quantum dot containing two levels. 63 They have assumed identical incoming rates but different outgoing ones for the two levels and analyzed the approximate CGF, derived by expanding in the ratio of the two outgoing rates, taken to be a small parameter. They have noticed a peculiar behavior: The agreement between the cumulants obtained by differentiating the approximate CGF (with respect to the counting field) and those derived by differentiating the CGF before expanding it was good for the low-order cumulants, but failed completely for the higher-order ones. In our case, the CGF is obtained as an expansion in the electron-phonon coupling and thus it is plausible that the unphysical results which we encountered around regime III for V > ω 0 at zero temperature, may have a similar origin to the apparent discrepancy reported in Ref. 61. This fault may be resolved by accounting for all orders in the electron-phonon coupling. However, getting analytical results seems to be technically complicated and almost inevitably requires numerical methods. 57 Another example of the FCS for interacting electrons is found in transport through a quantum dot in the Kondo regime. Recently, Sakano et. al. have calculated the CGF for the SU (N )−impurity Anderson model, 20 using the renormalized perturbation theory. Their result is exact up to cubic order in the source-drain bias voltage O(V 3 ) for a particle-hole symmetric case, and is quadratic in u. This is interpreted as the sum of the CGF for single-particle transfer and that for two-particle transfer. In our case, the phonon part of the CGF at zero temperature, to leading order in V , can be obtained from Eqs. (47) and (54) in the form where u 0 = 1 − 1/T . One may interpret the first righthand side of Eq. (67) as a sum of independent n-electron transfers. However, the second line indicates that the CGF is non-analytic on the negative real axis of the u−plane at u 0 ≤ 0. Therefore, within the first order expansion in V , the electron transfer statistics may be reduced to that of noninteracting electrons. 30 This example suggests that there exist certain subtleties in interpreting the CGF for interacting electrons.

V. SUMMARY
We have investigated the full-counting statistics of currents mediated by elastic and inelastic electron-phonon scattering. In the extended wide-band limit, we obtained analytic expressions for the cumulant generatingfunction, accurate up to second order in the electronphonon coupling. Our results are applicable for finite temperatures and bias voltages and satisfy the fluctuation theorem. Using those we analyzed the locations of singularities of the CGF. The singularities are symmetrically distributed in the λ−plane as to obey the fluctuation theorem.
The singularities in u− plane, which appear because of the electron-phonon interaction, classify specific regimes in which the dependence of the electron transfer statistics on the bare transmission is distinct. For small transmission probabilities we find singularities of the CGF on the positive real axis satisfying u > 1. Around the singularities, the CGF is non-convex, which results in a kink of the rate function. Such a kink, derived within the largedeviation analysis resembles a first-order phase transition in thermodynamics. It signifies the tendency of the phonon scattering to enhance the probability to find currents larger than the average value.
When the bias voltage is larger than the phonon frequency, V > ω 0 , we find singularities in 0 < u < 1 around perfect transmission. This results in a kink at I = 0 and a strong reduction for I < 0 in the rate function. This behavior can be understood in the following way: In this regime, phonons scatter electrons inelastically opposite to the bias. This broadens the rate function but the probability for current flowing opposite to the bias voltage is suppressed exponentially at low temperatures in accordance with the fluctuation theorem.
The kink structures in the rate function characterizes the electron-phonon interactions. Although measurements of the rate function of molecular junctions would be technically demanding, the FCS can be in principle monitored experimentally, 64 as is proven by existing measurements of higher cumulants 59,65 and of the FCS itself, 46,66 for metallic and semiconducting nanostructures. This gives hope that our predictions could be put to test. tion (R2214) from the JSPS, and by the MEXT kakenhi "Quantum Cybernetics". This paper was written while the authors were members of a research group on molecular electronics at the Institute for Advanced Studies, Jerusalem.
Appendix A: The Luttinger-Ward potential The way to construct a partition function based on the self energy is to exploit the Luttinger-Ward functional approach, 51 or the self-consistent Φ-derivable approximation. 50 This method can be straightforwardly extended to a nonequilibrium situation. 52 The saddlepoint approximation for the CGF can also be constructed by this approach. 48 The underlying idea is the observation that by introducing the Luttinger-Ward functional Φ, which includes all skeleton diagrams, the total generating functional can be formally written as where the trace and the product mean integrations over time along the Keldysh contour. The Green function G, also defined on the Keldysh contour, is where the arguments 1, 2, · · · stand for t 1 , t 2 , · · · . Here g λ is the Green function of the noninteracting electrons .
Both functions G and Σ depend implicitly on the counting field only through g λ . By differentiating Eq. (A1) with respect to the counting field λ and using Eq. (A5) one obtains where we have used the relation dΦ dλ = Tr δΦ δG ∂G ∂λ .
It should be emphasized that in this formulation the self energy has to be determined self-consistently in order to satisfy conservation laws. 52 To second order in the electron-phonon coupling γ, Eq. (A6) becomes where we have replaced Σ(G) by Σ(g λ ), since it is already O(γ 2 ). Upon using the identity and consequently Comparing this expression with the original one, Eq. (A1), one notes that the term depending explicitly on the self energy has disappeared. The lowest-order scaled Luttinger-ward potentialΦ (2) (λ) is obtained by expand-ingΦ(g λ ) (which may depend on g λ and not on G) up to O(γ 2 ) [see Eq. (18)]. For brevity, the calculation above is presented in the time domain. One may also switch to the frequency representation in which Eq. (A6) reads The Pauli matrix τ 3 , (A13) appears once we project the time from the Keldysh contour on the real time axis, where t ± ∈ K ± (see Fig. 3). This leads to Eqs. (16) and (17) in the main text. Note that Eq. (A6) or (A12) corresponds to the "generalized current expression" given by Eq. The Fock term depends on the equilibrium phonon distribution n + , which suggests that a re-summation of an infinite series of diagrams is needed in order to account for the nonequilibrium phonon distribution. We carry out this summation within the random-phase approximation, by summing over all ring diagrams [see diagrams (c) and (d) in Fig. 4]. The RPA is also known to be relevant for the AC conductance. 69 In this way, we obtain the functional where we have also included the zeroth-order term, Eq. (B1) and the full phonon propagator, D λ , Equation (B10) yields Eq. (32) in Fourier space. This RPA can be also formulated using the Keldysh path-integral approach within the saddle-point approximation 48 and further accounting for the Gaussianfluctuation correction around it. A comment on the FT and current conservation in the diagrammatic expansion is called for. The FT, including current conservation as represented by Eq. (15), has been proved using perturbation expansion in the interaction. 40 Although the proof has been constructed for the Coulomb interaction, it may be extended to the electron-phonon interaction case as well.