Strong coupling in thermoelectric nanojunctions: a reaction coordinate framework

We study a model of a thermoelectric nanojunction driven by vibrationally-assisted tunneling. We apply the reaction coordinate formalism to derive a master equation governing its thermoelectric performance beyond the weak electron-vibrational coupling limit. Employing full counting statistics we calculate the current flow, thermopower, associated noise, and efficiency without resorting to the weak vibrational coupling approximation. We demonstrate intricacies of the power-efficiency-precision trade-off at strong coupling, showing that the three cannot be maximised simultaneously in our model. Finally, we emphasise the importance of capturing non-additivity when considering strong coupling and multiple environments, demonstrating that an additive treatment of the environments can violate the upper bound on thermoelectric efficiency imposed by Carnot.


I. INTRODUCTION
The understanding of quantum transport and quantum thermodynamics is at the forefront of designing practical quantum devices. The study of molecular and solid-state nanojunctions combines these two areas, with the aim of realising functional quantum circuitry [1][2][3][4][5][6][7][8][9][10][11][12][13]. Even minimal nanojunction models often require the inclusion of several different environments, which also presents an interesting theoretical challenge. Typically, these environments represent two key features of a quantum circuit: the electronic leads, which are modelled as fermionic reservoirs; and the vibrational effects of the molecule or solid-state lattice, modelled by a reservoir of bosonic modes.
For both molecular and solid-state nanojunctions understanding the influence of electron-vibrational coupling on transport properties is of crucial importance, and has been studied widely (see for example ). In particular, strong vibrational coupling can result in striking qualitatively changes to transport behaviour, such as the suppression of low-bias current known as Franck-Condon blockade [15,18].
In this work, rather than focussing on electron transport driven by a voltage bias, we study the influence of strong vibrational coupling on the thermoelectric performance of a model nanojunction, where current is driven against a bias by a temperature gradient. The system thus generates a finite power output with an associated efficiency, determined also by the heat flux. We are particularly interested in understanding the interplay between the power, efficiency, and noise, which is of key importance in assessing the performance of thermodynamic devices [39][40][41][42][43][44][45][46].
We consider a two-site nanojunction model where inter-site transport is mediated via phonon-assisted tunneling (PAT). This is a natural setting to study the impact of vibrational coupling in a thermoelectric context as phonons are then directly responsible for promoting electron transport. To model the system at strong electron-phonon coupling we employ the reaction coordinate master equation (RCME) approach [32,37,[47][48][49][50]. This involves applying a unitary Hamiltonian-level mapping to incorporate a single collective mode from the phonon environment into our electronic system, generating what we shall term an augmented system. The remaining phonon modes are then weakly coupled to the augmented system and traced out within the standard Born-Markov (BM) approximations, though retaining the dominant strong-coupling and non-Markovian effects on the original electronic system. Several other theoretical techniques have been applied to related models of vibrationally coupled charge transport, including investigations of thermoelectric performance. Examples include Green's function approaches [22,23,28,31], linear response methods [25], and full counting statistics with the hierarchical equations of motion [38] and quantum master equations [26,31], where strong lead couplings have also been considered. The impact of strong light-matter (rather than vibrational) coupling in thermoelectric devices has also been explored [51].
One significant advantage of the RCME approach is that it allows for a straightforward incorporation of nonadditive effects between the strongly-coupled phonons and the fermionic leads, by deriving the lead influences within the full eigenbasis of the augmented system. This is crucial for studying thermoelectric systems since they operate at finite bias, and an additive treatment in which the leads are insensitive to the strong system-phonon coupling can only be justified in the infinite bias limit [37]. In fact, our non-additive and additive comparisons show directly that the additive approximation can lead to unphysical artefacts in thermoelectric performance no matter how sophisticated the treatment of phonons, such as thermoelectric efficiencies beyond the Carnot bound.
We shall also compare the RCME treatment to a weak coupling master equation (WCME) that treats all environments, including the phonons, within the standard BM approximations. This confirms consistency between β L ,μ L βR,μR λ ε L εR β ph Γ L ΓR FIG. 1. Schematic diagram of the two-site thermoelectric model. Each site is connected to an individual lead and both sites are coupled through a common phonon bath at inverse temperature β ph , with reorganisation energy λ. The site energies are given by i, couplings to the leads are denoted Γi, and each lead is held at a chemical potential µi and inverse temperature βi, where i = L, R.
all methods at weak enough phonon coupling strength, and highlights the limitations of the WCME as the phonon coupling is increased. For all methods we employ counting statistics techniques [15,29,30,[52][53][54] to calculate the transport cumulants necessary to evaluate the thermoelectric power output, noise, and efficiency.
The paper is organised as follows. In Sec. II we outline the theoretical model and techniques used, including the WCME, the reaction-coordinate (RC) mapping, and the implementation of counting statistics. In Sec. III we consider a thermoelectric regime in which a temperature gradient is applied across the leads, with the left lead acting as a resource such that the phonons not only mediate transport but also act as a heat loss mechanism. We contrast this with the case of a phonon resource in Sec. IV. Finally, we summarise our results and discuss future directions in Sec. V.

II. THEORETICAL MODEL
Our two-site model is shown schematically in Fig. 1. Each site is coupled individually to a fermionic lead and both are coupled to a common phonon bath that mediates electron transport. This could describe, for example, a molecule or double quantum dot with two regions of high electron density (the dots), coupled via molecular vibrations or the solid-state lattice. Direct inter-site tunneling can also be included straightforwardly, but to isolate the impact of vibrations we shall not do so here.
We consider the left lead to be held at a lower chemical potential (µ L ) than the right (µ R ). In order for electrons to overcome the chemical potential gradient, leading to a useful current flow, we require a temperature gradient. In the first instance (Sec. III) we set the left lead temperature to be greater than that of the right lead or phonons, implying β L < β R = β ph for the inverse temperatures. Subsequently, in Sec. IV, we shall set the phonon temperature to be greater than that of both leads, which we fix equal (β L = β R = β el ), so that β ph < β el and the phonons become the resource. The right dot is also held at a higher energy than the right lead to promote electron flow from left to right and suppress the back-flow of electrons. As a consequence some of the work done to push the electrons onto the right dot is lost as heat when they are transferred to the right lead. However, as µ R > µ L we still have a positive work output overall.
The Hamiltonian for our set-up is with system Hamiltonian where d L and d R are fermionic annihilation operators for the left and right dot, respectively, L ( R ) is the energy of the left (right) dot, and U is the Coulomb energy. We assume that on-site Coulomb repulsion is large such that each dot is never occupied by more than a single excess electron. The terms corresponding to the fermionic leads are and their interactions with the system are given by Here ki is the energy of the k th fermionic mode in lead i, which is coupled to the system via t ki , and c ki are fermionic annihilation operators for electrons in the leads.
As the system and lead fermionic operators anticommute, the interaction Hamiltonians H IL and H IR do not have the tensor product structure required for our subsequent master equation derivations. However, we can use a Jordan-Wigner transformation to rectify this [30], after which the state space for the central system can be written {|G , |L , |R , |D }. Here |G is the ground state, with neither dot occupied by an excess electron, |L (|R ) denotes an electron present the on the left (right) dot, and |D is the state in which both dots are occupied. Then, the system and interaction Hamiltonians can be written as and respectively, where the negative terms in the left lead couplings preserve the system fermionic anti-commutation relations after the Jordan-Wigner transformation.
We consider a phonon environment coupled to the system in Brownian oscillator form [32,55,56] with system coupling operator Here, x q and p q are phonon position and momentum operators satisfying [x q , p q ] = iδ qq ( is set to 1 throughout), ω q is the frequency of the q th phonon mode, which couples to the system with strength c q . Expressing the phonon position and momentum operators in terms of standard bosonic creation and annihilation operators a † q and a q , the phonon environment Hamiltonian becomes (ignoring a constant term) Likewise, after applying the Jordan-Wigner transformation the interaction Hamiltonian is given by with h q = −c q / 2ω q . The final term in Eq. (12) compensates for energetic renormalisation of the system caused by the phonon coupling [56] and ensures that the sites remain within the lead bias window even for strong phonon interactions.

A. Weak coupling master equation
For sufficiently weak coupling between the system and environments we expect to be able to apply the standard BM approximations to arrive at a Redfield master equation [57]. This treats the environments additively, meaning that each environment influences the system in the same manner as it would if treated in isolation.
We write the lead interaction Hamiltonians as and with Considering H IL first, and following the textbook BM procedure, we arrive at the following Liouvillion for the left lead contribution to the dynamics of the reduced system density operator ρ S (t) [57], Here, the time-dependent system operators and bath correlation functions are given by where is the reduced density operator of the left lead, which we have taken to be in equilibrium. We also define N Ei = k c † ki c ki and the Fermi-Dirac distribution f i (ω) = (e βi(ω−µi) + 1) −1 , where β i is the inverse temperature of fermionic environment i = L, R and µ i its chemical potential.
The right lead contribution follows in complete analogy, where For the phonon contribution we have , and the final term in Eq. (12) can be neglected in the weak-coupling limit as it cancels with an energetic shift obtained within the BM approximations [58]. The resulting phonon Liouvillian reads Here is written in terms of the Bose-Einstein distribution n(ω) = (e β ph ω − 1) −1 for a thermal equilibrium phonon bath ρ ph = e −β ph H E ph /tr(e −β ph H E ph ).
Combining Eqs. (15), (20) and (25) we arrive at the final WCME for our system evolution, written in the Schrödinger picture aṡ Throughout the following we consider a phonon spectral density of Drude-Lorentz form with peak around ω 0 , broadening parameter γ, and reorganisation energy λ. For both leads we take the wideband limit, B. RC mapping for phonon-assisted tunnelling

Non-additive RCME
To move beyond the limitations of the WCME treatment of the phonons we now apply the unitary reaction coordinate mapping to include a collective mode of the phonon environment within an enlarged augmented system [32,37,[47][48][49][50] (see Fig. 2). Following the derivations Schematic of the reaction coordinate (RC) mapping applied to the two-site thermoelectric model. The interactions are treated non-additively, meaning that the leads couple to the augmented system S which is obtained after the mapping. This is comprised of the original system S and the RC of frequency Ω, interacting with strength κ. E represents the residual bosonic environment and other symbols are as defined in Fig. 1.
given in [32,47,50] and using Eqs. (7) and (12) the augmented system Hamiltonian becomes Here, the RC mode is defined such that with system-RC coupling κ 2 = (1/Ω) q ω q h 2 q and mode frequency Ω = ω 0 . Note that we must now explicitly include the final term in Eq. (12) within our mapped Hamiltonian as the phonons are no longer treated perturbatively, leading in Eq. (28) to the shift λ = ∞ 0 dωJ(ω)/ω as the continuum limit of q h 2 q /ω q . As a result of the mapping the remaining bosonic modes are transformed such that they are now decoupled from the original system, but couple directly to the RC via where b k (b † k ) are bosonic annihilation (creation) operators. To compute the augmented system dynamics, it is not necessary to derive precise forms for the interaction strengths f k or frequencies ν k , but rather it suffices to obtain a mapped spectral density J RC (ν) = k f 2 k δ(ν −ν k ). This can be achieved through matching Heisenberg equations of motion for the dynamics before and after the mapping, leading to J RC (ν) = γν/2πω 0 [32,47,50]. As the mapping has no effect on the fermionic leads, the total Hamiltonian now becomes Following the procedure described in [37] we trace out the leads and the residual bosonic bath within the BM approximations to derive a master equation for the augmented system density operator ρ S (t), which describes the original system plus RC, and is given bẏ Here A m for m = 1, 2, 3, 4 are the same system operators as defined in the previous section, A ph = (a † + a) and where i = L, R as before. Defining the system-RC eigenbasis through Eqs. (33) and (34) describe the action of the residual bosonic bath on the mapped system through its coupling to the RC. Eqs. (35) and (36) encode the influence of the leads, and have been derived within the full system-RC eigenbasis so that the leads are "aware" of the strong system-phonon coupling. Thus the treatment is inherently non-additive with respect to the lead and phonon couplings, in contrast to the WCME. Note that, despite the BM approximation applied to the residual bosonic bath, tracing over the RC in ρ S (t) results in non-Markovian dynamics for the reduced system state ρ S (t) that is also valid at large phonon reorganisation energies λ, again in contrast to the WCME [32,47,48].

Additive RCME
As we have just seen, when phonon coupling is strong the RC treatment lends itself naturally to a non-additive description of the fermionic leads in our thermoelectric model. This may not, however, be the case for other techniques that can treat phonons beyond the weak-coupling limit, and in such situations it might only be feasible to include further environments phenomenologically. It is thus important to assess the impact of such a simplifying assumption on thermoelectric performance, which can be done through an additive reaction-coordinate master equation (aRCME).
The aRCME is obtained straightforwardly from the full RCME of Eq. (32) by replacing the left and right lead terms with their respective weak-coupling counterparts given in Eqs. (15) and (20). Specifically, the aRCME readsρ where I RC denotes the identity operator applied within the RC Hilbert space, such that the leads couple only to the original (unmapped) two-site system and are unaffected by the strong phonon coupling.

C. Counting Statistics
We characterise the thermoelectric performance of our setup by tracking the counting statistics of electrons that enter or leave the right lead under steady-state conditions. Following Refs. [53,54] we define the matter current flow (first cumulant or mean I ) and the associated noise (second cumulant or variance I 2 ) as where I = I + − I − , J = I + + I − , |0 is the right eigenvector corresponding to the zero eigenvalue of the Liouvillian under consideration (WCME, RCME, or aR-CME), and 0 | is the identity operator. The current operator I + (I − ) represents the terms in the relevant Liouvillian that add (remove) electrons to (from) the right lead. Defining projection operators P = |0 0 | and Q = I − P we have a pseudo-inverse R = QL −1 Q, which is well defined as the inversion is only performed on the subspace spanned by Q.
Obtaining the thermoelectric power output P is straightforward once we have the current, The thermoelectric efficiency is defined as where Q in is the heat flow entering the system. In Sec. III we have T L > T R = T ph such that the left lead acts as the resource and Q in = Q L . Since the total number of electrons is conserved, the steady-state current from the left lead is given by I . In the following examples we always take U to be by far the largest energy scale, such that we work in the Coulomb blockade regime and the state |D is never occupied. In the WCME and the aRCME this means that each electron has an associated energy L , and the left lead energy current becomes I E L = L I . In the full RCME however, the left lead now couples to a manifold of states within the augmented system and we cannot simply consider a single energy scale L . Instead, we use the fact that in the steady-state tr(H S ρ S ) = 0 to obtain the left lead energy current from Eq. (32) as Here ρ S (∞) is the RCME steady-state obtained from Eq. (32) by settingρ S = 0. Note that we could also have calculated I E L for the WCME and the aRCME in an analogous way, though the result would be the same as setting I E L = L I for these cases. The heat flow can be calculated in all cases from Q L = I E L − µ L I using the appropriate energy current.
In Sec. IV, it is the phonons that act as a resource rather than the left lead, such that Q in = Q ph . As no work is exchanged between the phonon environment and the system, I E ph = Q ph and the phonon heat flow is equal to the phonon energy current. We may then use conservation of energy in the steady-state to calculate the heat flow from I E L + I E R + I E ph = 0, where I E R here is positive for energy entering the system from the right lead.

III. THERMOELECTRIC REGIME I: LEFT LEAD RESOURCE
In our first example we consider the impact of strong electron-phonon coupling when a temperature gradient is introduced between the left and right leads. The phonon bath is held at the same temperature as the right lead (the colder of the two) so that it is only able to act as a transport mediator without being a thermodynamic re-

source.
We begin in Fig. 3(a) by studying the thermoelectric current flow as a function of phonon coupling (reorganisation energy, λ) for a fixed bias voltage. We see that at small reorganisation energies the WCME and RCME agree, as expected, and both predict an increase in current with λ due to stronger phonon-mediated tunneling between the left and right dot. As the phonon reorganisation energy becomes larger, and we move out of the validity range of the BM approximations, the WCME incorrectly predicts a plateau in the current due to saturation of the available transport channels. Specifically, in the WCME transport proceeds via sequential incoherent transitions with electrons tunneling in turn between the left lead and the left dot, then via phonons between the left and right dots, and finally tunneling onto the right lead. As the reorganisation energy is increased, phonon mediated tunneling between the left and right dots becomes increasingly rapid (the WCME single-phonon rate is proportional to λ), until the tunneling rate between the dots and the leads becomes the limiting factor (sat-uration) and the current plateau.
For small reorganisation energies, transport proceeds in the same manner within the RCME and single-phonon processes dominate. In contrast, an enhancement of current is observed in the RCME predictions for λβ R ∼ 1 − 100. In this regime multi-phonon processes that are ignored in the WCME, but are captured through the coupling of the leads to the numerous eigenstates of the augmented system within the RCME, become increasingly important, opening up further transport channels that allow for greater current flow. Furthermore, the energy levels of the dominant transport channels are renormalised by the strong system-phonon coupling; they are increased in comparison to µ L but differences between them are reduced in comparison to the WCME, allowing tunnelling to occur more readily. Nevertheless, for very large reorganisation energies the RCME predicts a complete suppression of current flow due to electron blockade, which again is not captured within the WCME. Similar to Franck-Condon blockade [15,18], in the RCME displacement of the RC becomes large at strong vibrational coupling. Hence overlaps between the augmented system states responsible for transport and the empty state, for which the RC is not displaced, are substantially suppressed, leading in turn to the reduction of current flow.
We now consider variations in the left lead temperature, shown in Fig. 3(b) for an intermediate reorganisation energy of λβ R = 3, where the RCME predicts a moderate increase in current compared to the WCME. Both the WCME and RCME predict current flow decreases as β L increases due to the reduction in temperature gradient between the leads, which suppresses the thermoelectric activity. Both also predict a crossover to negative current flow once the relevant stopping voltage is surpassed and the current then follows the chemical potential (rather than temperature) gradient. In fact, the RCME predicts a slightly more resilient thermoelectric, with the chosen value of V acting as the stopping voltage for a smaller temperature gradient than for the WCME.
This can also be seen in Fig. 3(c), where the RCME predicts a larger current for all bias voltages and an extended stopping voltage compared to the WCME prediction of V S = ( L − µ L )(β R − β L )/β R . Again, this increase in stopping voltage is due to the non-additive interplay between the phonon bath and the fermionic leads, encapsulated through the RC mapping and subsequent master equation derivation, which effectively increases the relevant system energies in comparison to µ L . We see a similar behaviour in the thermoelectric power in Fig. 3(d). For a constant bias, the power as a function of reorganisation energy is of the same form as the current in Fig.  3(a), simply scaled by the value of V .
Having considered the current flow and power output we now turn our attention to the associated thermoelectric efficiency, shown as a function of bias voltage in Fig. 4(a). In the WCME case the efficiency increases linearly with V up to a maximum of the Carnot efficiency η C at the stopping voltage. The RCME prediction is qualitatively different, however, and we see that the increase in power due to strong phonon coupling comes at the expense of a reduced efficiency for heat to work conversion, and the inability of the system to reach the Carnot efficiency even for vanishing power output. This drop in efficiency is caused by the power decreasing towards zero more quickly than the heat flow, with the left lead providing heat energy that is preferentially absorbed by the phonon bath rather than promoting electrons to tunnel against the increasing bias. This can be thought of as a cost imposed by the strong electron-phonon coupling, and the resultant power-efficiency trade off can be seen by comparing the RCME curves in Fig. 4(b), which shows efficiency as a function of phonon coupling strength, and in Fig. 3(a). In contrast, the WCME efficiency is unaffected by increasing the phonon reorganisation energy, another failing of the BM approximations. Instead, we have which is independent of λ, and for V → V S becomes η WCME = (β R − β L )/β R = η C . Likewise, the WCME fails to predict any variation of the stopping voltage with phonon coupling, as shown in Fig. 4(c). On the other hand, the RCME predicts an increase in stopping voltage with λ, consistent again with a larger parameter regime supporting thermoelectric behaviour at strong phonon coupling. These features of strong phonon coupling are also evident in the parametric plots of power and efficiency shown in Fig. 4(d), which are calculated from the RCME by fixing the reorganisation energy and then sweeping through the entire bias window for which the system operates as a thermoelectric. As we increase the bias voltage the power and efficiency both increase until a maximum power is reached, at which point the power begins to decrease back to zero, while the efficiency increases to a maximum and then also returns to zero. Increasing the reorganisation energy results in an increase in the maximum power attainable, however the maximum efficiency is reduced as a result. This trade-off implies that a system with a given level of phonon coupling is suited to different roles (e.g. maximising efficiency or power) depending on how weak or strong the reorganisation energy is. We also note that for increasing λ the parametric plot maps out a larger area due to the increased bias range over which the system acts as a thermoelectric.
We now extend our analysis to consider higher cumulants, specifically the zero frequency noise I 2 , and potential trade-offs between fluctuations and the thermoelectric power and efficiency. We are motivated in part by recent advances in the study of thermodynamic uncertainty relations (TURs) [39][40][41][42][43][44][45][46]. Steady-state TURs provide cost-precision trade-off relations whereby the relative uncertainty in the current is constrained by the entropy production, implying that for a more precise output with less noise, there will be a greater cost (higher entropy production). Conversely, if we want to minimise entropy production to maximise efficiency, this is expected to come at the cost of increasing the relative uncertainty. In Fig. 5(a) we plot the relative uncertainty as a function of phonon reorganisation energy for both the WCME and RCME. As the phonon coupling strength increases, the relative uncertainty decreases, corresponding to a reduction in the noise as compared to the increasing current (see Fig. 3(a)). In particular, in the regime of enhanced current above λβ R ∼ 1 the RCME predicts a significant suppression of υ. Here, power increases and fluctuations decrease, but the trade-off is that efficiency decreases as well. We see that the WCME overestimates υ in the same region, in some cases by over an order of magnitude, and it also fails to capture a very rapid increase in the relative uncertainty in the electron blockade regime for λβ R > 100. We now fix the reorganisation energy to λβ R = 3 and instead scan the voltage V across the relevant bias window for which the system acts as a thermoelectric, producing parametric plots of the relative uncertainty against power (Fig. 5(b)) and efficiency (Fig. 5(c)). The qualitative features of the relative uncertainty against power plot (Fig. 5(b)) are similar for both the RCME and WCME, though absolute values can vary substantially. For smaller values of the bias voltage, power can be increased at a relatively modest increase in υ, but once V approaches the stopping voltage, the relative uncertainty begins to diverge rapidly and the current becomes extremely noisy. In this regime the weak-coupling analysis predicts the entropy production rate to fall to zero and the efficiency to approach the Carnot bound, see Fig. 5(c). The RCME prediction is qualitatively different, however. As the stopping voltage is approached in the RCME case υ still diverges, but the efficiency now falls to zero, consistent with our earlier analysis in Fig. 4(a).

IV. THERMOELECTRIC REGIME II: PHONON RESOURCE
In the previous section we considered the left lead as the source of heat, with the phonons acting as a transport mediator as well as a heat sink. In this section we shall instead treat the phonons as a thermodynamic resource, allowing us to draw comparisons between the two regimes of operation. We thus set β ph < β el , where β el = β L = β R is the inverse temperature of the leads. We retain the original chemical potential gradient, µ R > µ L , so it is only by exploiting the temperature gradient between the phonons and the leads that we can generate a positive current flow.

A. Thermoelectric performance
Beginning with Fig. 6(a) we see that the variation of the current with phonon coupling strength is qualitatively similar to the case when the left lead was the resource [ Fig. 3(a)], though now the WCME overestimates the current in comparison to the more accurate RCME. As a result of the dual role played by the phonon bath, as both thermodynamic resource and transport mediator, the RCME predicts a maximised current flow that is much larger and occurs for a much smaller reorganisation energy than in Fig. 3(a), around λβ el = 1. However, the multi-phonon processes captured by the RC mapping do not enhance the current here in comparison to the WCME, due to competition between phonon-mediated excitation and de-excitation within the manifold of augmented system states caused by the phonon bath now being at the highest temperature.
Similar trends are also reflected in the current and power behaviour with changes in bias voltage shown in Figs. 6(b) and 6(c), respectively. Here we find that the stopping voltage is now reduced in the RCME treatment as compared to the WCME, and so the weak-coupling approximation overestimates the parameter range over which a thermoelectric current flows. In analogy with the left-lead resource regime, the stopping voltage is now given by V S = ∆(β el −β ph )/β el in the weak-coupling case. We noted previously that the RCME demonstrates an increased energy difference between the chemical potential of the left lead and the augmented system eigenvalues, thus increasing the stopping voltage when the left lead is the resource. In the phonon resource regime we find the reverse to be the case. It is now the energetic differences between the relevant tunnelling states that are important, and within the augmented system these are pushed closer together, resulting in a reduced stopping voltage as shown in Fig 6(d). Once again, the WCME incorrectly predicts no variation in stopping voltage with changing reorganisation energy.
Figs. 7(a) and (b) show the thermoelectric efficiency as a function of bias voltage and reorganisation energy, respectively. Again, the trend with varying bias voltage is similar to the left-lead resource regime, with the RCME predicting an increase in efficiency before a sudden drop to zero at the stopping voltage. The WCME generally underestimates the efficiency apart from in the regime beyond the point at which the RCME predicts the stopping voltage to occur. When fixing the bias voltage, as in Fig. 7(b), an interesting trend is seen whereby the efficiency increases with increasing phonon reorgan-isation energy, before decreasing sharply for large coupling strengths. Despite this, the efficiency at maximum power, shown in Fig. 7(c) decreases monotonically with increasing reorganisation energy. This is caused by a shift in the bias voltage that produces the maximum power output, which decreases for larger reorganisation energy, and so the efficiency does as well. In fact, considering the full bias range for which the system acts a a thermoelectric, as in Fig. 7(d), reveals that the maximum obtainable efficiency is greater for smaller reorganisation energies. The maximum obtainable power does increase for larger λ, but the behaviour is non-monotonic, such that for λβ el = 10 both the maximum power and efficiency that can be reached are lower than for λβ el = 1.
Turning to current fluctuations in the phonon-resource regime, we plot the relative uncertainty (υ) as a function of reorganisation energy in Fig. 8(a). The trend is similar to the case of a left-lead resource, with υ decreasing as the reorganisation energy increases and more current flows, until it reaches a minimum near the current maximum. For larger reorganisation energies the relative uncertainty increases rapidly primarily due to the suppression of current. As expected, the WCME only cap- , and the RCME once more predicts that the Carnot efficiency cannot be reached, even for diverging relative uncertainty and vanishing power as the stopping voltage is approached.

B. Breakdown of additivity
Finally, we shall explore the breakdown of environmental additivity within our thermoelectric model to highlight the importance of accurately capturing nonadditive effects for multiple environments beyond weakcoupling. This is the case even if only a single environ-ment is strongly coupled to the system, as we demonstrate through use of the aRCME. The aRCME is an additive simplification of the RCME wherein the electronphonon coupling is accurately modelled for strong reorganisation energies, but the lead couplings are treated phenomenologically, i.e. they are derived assuming no alterations to the system caused by the strong electronphonon coupling.
In Fig. 9 we consider the current flow (a) and relative uncertainty (b) predicted by the aRCME for varying phonon reorganisation energy, and compare them with the previous results for the WCME and the full (nonadditive) RCME in the phonon-resource regime. For the current flow we see that the aRCME completely misses the electron localisation predicted by the RCME at large reorganisation energies, and instead incorrectly predicts a continually increasing current with larger phonon coupling strength. This is also reflected in the relative uncertainty, which simply decreases with increasing reorganisation energy. Thus, it is not only the quantitative, but also the qualitative behaviour predicted by the aRCME that is incorrect. This occurs due to the inconsistent ap- plication of the RC mapping within the aRCME, whereby the leads count on the original system transport channels rather than those of the mapped (augmented) system. In contrast to the augmented system states, there is no suppression of vibrational overlaps between the original system states and the empty state, meaning that the aR- CME is then unable to capture the localisation regime. Furthermore, this implies that the power and heat flow definitions given by the aRCME are not thermodynamically consistent, as they are defined with respect to the original system states, yet the non-equilibrium steadystate is calculated with respect to the augmented system in the aRCME (albeit with phenomenological lead dissipators).
Though this failing is clear in the context of the RC treatment, more generally it highlights the importance of applying approximations consistently where multiple environments are concerned and at least one is to be treated beyond weak-coupling. A stark example of these effects is shown in Fig. 10, where the predictions for efficiency as a function of bias voltage match between the aRCME and WCME since they count on the same transport channels. Nevertheless, as the aRCME definitions also (incorrectly) give rise to an increased stopping voltage, it has the potential to breakdown completely, predicting unphysical thermoelectric efficiencies beyond the Carnot limit.

V. CONCLUSIONS
We have combined the RCME and counting statistics techniques to treat strong phonon-coupling in thermoelectric nanojunction models, avoiding the limitations of the standard Born-Markov WCME. Interestingly, we found discrepancies between the two approaches not only at very strong phonon coupling, as would be ex-pected, but also for relatively moderate phonon coupling strengths. For example, the WCME incorrectly predicts efficiencies and stopping voltages that are completely independent of the phonon reorganisation energy, highlighting the importance of an accurate representation of the electron-phonon coupling. We have also emphasised potential pitfalls of employing a phenomenological additive treatment of the influence of multiple baths when at least one couples strongly to the system. In the present context, this was shown to lead to an unphysical violation of the Carnot bound to the thermoelectric efficiency.
More generally, though there are common trends, the thermoelectric behaviour we have found appears also to be sensitive to the chosen parameter regimes, and whether strong coupling is beneficial or detrimental depends to a certain extent on which figure of merit is being optimised. Nevertheless, we do find that trade-offs seem to exist between power, efficiency, and fluctuations at strong coupling, and it would be interesting to explore whether this could be formalised in a more modelindependent manner. In future work, we also plan to explore the interplay between strong lead and phonon couplings by extracting RCs not only for the vibrational degrees of freedom but also for the fermionic environments [50,59].