Thermodynamic bounds on current fluctuations Von der Fakultät Mathematik und Physik der Universität Stuttgart zur Erlangung der Würde eines Doktors der Naturwissenschaften (Dr. rer. nat.) genehmigte Abhandlung vorgelegt von Patrick Pietzonka aus Stuttgart Hauptberichter: Prof. Dr. Udo Seifert Erster Mitberichter: Prof. Dr. Siegfried Dietrich Zweiter Mitberichter: Prof. Dr. Andreas Engel Vorsitzender: Prof. Dr. Tilman Pfau Tag der Einreichung: 26. März 2018 Tag der mündlichen Prüfung: 11. Juni 2018 II. Institut für Theoretische Physik der Universität Stuttgart 2018 Contents Publications 7 Zusammenfassung 9 Summary 15 1 Introduction 19 2 Stochastic thermodynamics on discrete state spaces 23 2.1 Markovian dynamics in closed systems . . . . . . . . . . . . . . . . . . 23 2.2 Markovian dynamics in open systems . . . . . . . . . . . . . . . . . . 26 2.3 Energetics of open systems . . . . . . . . . . . . . . . . . . . . . . . . 29 2.4 Entropy production . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Time extensivity and time reversal . . . . . . . . . . . . . . . . . . . . 33 2.6 Fluctuation relations . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.7 Unicyclic molecular motor as paradigmatic example . . . . . . . . . . . 38 2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3 Fluctuations in stationary Markov processes 43 3.1 Cumulant generating function for time-additive observables . . . . . . . 43 3.2 Calculation of low-order cumulants . . . . . . . . . . . . . . . . . . . . 49 3.3 Large deviation functions . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4 Level 2.5 large deviation theory . . . . . . . . . . . . . . . . . . . . . . 58 3.5 Continuous state spaces . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4 Dissipation-dependent bound on current fluctuations 67 4.1 General setup for bounds on current fluctuations . . . . . . . . . . . . . 67 4.2 Linear response regime . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.3 Beyond linear response: Unicyclic case . . . . . . . . . . . . . . . . . 73 4.4 Beyond linear response: multicyclic case . . . . . . . . . . . . . . . . . 75 4.5 The thermodynamic uncertainty relation . . . . . . . . . . . . . . . . . 77 4.6 Dissipation-dependent bound in driven diffusive one dimensional systems 78 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 3 Contents 5 Application: Bounds on the efficiency of molecular motors 81 5.1 Thermodynamically consistent motor model . . . . . . . . . . . . . . . 81 5.2 The bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.3 Variants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.4 Numerical case study . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6 Affinity- and topology-dependent bound on current fluctuations 89 6.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 6.2 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.3 Numerical illustration . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.4 Uniform cycle decomposition . . . . . . . . . . . . . . . . . . . . . . . 95 6.5 Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.6 Vanishing currents and equilibrium . . . . . . . . . . . . . . . . . . . . 101 6.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Appendices to chapter 6 103 6.A Proof of the upper bound on the generating function for unicyclic networks103 6.B Proof of the concavity of ψ(ζ, s) in s . . . . . . . . . . . . . . . . . . . 105 6.C Proof of the monotonic decrease of ψ(ζ, s)/s in s . . . . . . . . . . . . 106 7 Activity-dependent and asymptotic bound 107 7.1 Activity-dependent bound . . . . . . . . . . . . . . . . . . . . . . . . . 107 7.2 Asymptotic bound for unicyclic networks . . . . . . . . . . . . . . . . 110 7.3 Asymptotic bound for multicyclic networks . . . . . . . . . . . . . . . 112 7.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 8 Finite-time generalizations 117 8.1 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.2 Experimental illustration . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.3 Bound on the generating function . . . . . . . . . . . . . . . . . . . . . 121 8.4 Short-time and linear response limit . . . . . . . . . . . . . . . . . . . 123 8.5 Illustration for unicyclic networks . . . . . . . . . . . . . . . . . . . . 124 8.6 Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 8.7 Uncertainty relation for transient non-equilibrium states . . . . . . . . . 129 8.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 9 Concluding perspective 133 References 137 Danksagungen 151 4 List of Figures 1.1 Uncertainty associated with a molecular motor walking along a filament 21 2.1 Illustration of the network representation of an open system . . . . . . . 28 2.2 Models for unicyclic motors . . . . . . . . . . . . . . . . . . . . . . . 39 3.1 Generic spectrum of the tilted operator . . . . . . . . . . . . . . . . . . 46 3.2 Large deviation function for the coin-toss experiment . . . . . . . . . . 53 3.3 Large deviation function for the asymmetric random walk . . . . . . . . 55 4.1 Unicyclic and multicyclic networks . . . . . . . . . . . . . . . . . . . . 68 4.2 Unicyclic illustration of the dissipation-dependent bound. . . . . . . . . 74 4.3 Multicyclic illustration of the dissipation-dependent bound. . . . . . . . 76 4.4 Schematic illustrations of lattice transport models . . . . . . . . . . . . 79 4.5 Dissipation dependent bound for lattice transport models . . . . . . . . 80 5.1 Molecular motor walking along a periodic track . . . . . . . . . . . . . 82 5.2 Efficiency bound applied to experimental data for kinesin . . . . . . . . 85 5.3 Numerical illustration for a motor-cargo complex . . . . . . . . . . . . 87 6.1 Bounds on the generating function for a unicyclic network with high affinity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2 Bounds on the generating function for unicyclic network with various affinities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.3 Affinity- and topology-dependent bound for a multicyclic network . . . 94 6.4 Illustration of the uniform cycle decomposition . . . . . . . . . . . . . 96 6.5 Illustration of the existence of uniform cycles . . . . . . . . . . . . . . 97 6.6 Properties of the function ψ(ζ, s) . . . . . . . . . . . . . . . . . . . . . 99 6.A.1Characteristic polynomial for unicyclic networks. . . . . . . . . . . . . 104 7.1 Example for the activity-dependent bound in a unicyclic network . . . . 109 7.1 Illustration of the asymptotic bound for a multicyclic network . . . . . . 112 7.1 Summary of bounds for a unicyclic network . . . . . . . . . . . . . . . 115 8.1 Experimental setup and data for the dichotomously switching optical trap 119 8.2 Work distributions for the dichotomously switching trap . . . . . . . . . 120 5 List of Figures 8.1 Illustration for the finite-time bound on the generating function for a unicyclic network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 8.2 Illustration for the finite-time bound on the generating function for a multicyclic network . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 6 Publications Parts of this thesis have been published in: • P. Pietzonka, A. C. Barato, and U. Seifert, “Universal bounds on current fluctuations”, Phys. Rev. E 93, 052145 (2016). © 2016 American Physical Society. Chapter 4 and parts of Chapters 6 and 7 are based on this publication. • P. Pietzonka, A. C. Barato, and U. Seifert, “Affinity- and topology-dependent bound on current fluctuations”, J. Phys. A: Math. Theor. 49, 34LT01 (2016). © 2016 IOP Publishing Ltd. Chapter 6 is based on this publication. • P. Pietzonka, A. C. Barato, and U. Seifert, “Universal bound on the efficiency of molecular motors”, J. Stat. Mech.: Theor. Exp. 124004 (2016). © 2016 IOP Publishing Ltd and SISSA Medialab srl. Chapter 5 is based on this publication. • P. Pietzonka, F. Ritort, and U. Seifert, “Finite-time generalization of the thermodynamic uncertainty relation”, Phys. Rev. E 96, 012101 (2017). © 2017 American Physical Society. Chapter 8 is based on this publication. Further publications I have co-authored: • P. Pietzonka, E. Zimmermann, and U. Seifert, “Fine-structured large deviations and the fluctuation theorem: Molecular motors and beyond”, EPL 107, 20002 (2014). • P. Pietzonka, K. Kleinbeck, and U. Seifert, “Extreme fluctuations of active Brownian motion”, New J. Phys. 18, 052001 (2016). 7 Publications • P. Pietzonka andU. Seifert, “Entropy production of active particles and for particles in active baths”, J. Phys. A: Math. Theor. 51, 01LT01 (2018). • M. Uhl, P. Pietzonka, and U. Seifert, “Fluctuations of apparent entropy production in networkswith hidden slowdegrees of freedom”, J. Stat. Mech.: Theor. Exp. 023203 (2018). • L. P. Fischer, P. Pietzonka, and U. Seifert, “Large deviation function for a driven underdamped particle in a periodic poten- tial”, Phys. Rev. E 97, 022143 (2018). • P. Pietzonka and U. Seifert, “Universal trade-off between power, efficiency, and constancy in steady-state heat engines”, Phys. Rev. Lett. 120, 190602 (2018). 8 Zusammenfassung Prozesse in lebendigen Systemen und nützlichen künstlichen Maschinen laufen unter Nichtgleichgewichtsbedingungen ab. Das bedeutet, dass diese Systeme sich im Kontakt mit mehreren Reservoiren befinden, die nicht in wechselseitigem thermodynamischem Gleichgewicht sind. Die Reservoire stellen Ressourcen wie Nahrung oder Treibstoff bereit oder fungieren als thermische Umgebung die Wärme absorbiert. Unter der Vor- aussetzung, dass das betrachtete System hinreichend klein gegenüber den Reservoiren ist kann man davon ausgehen, dass der Zustand der Reservoire sich auf der relevanten Zeit- skala nur unwesentlich ändert. Wenn außerdem das System nicht von außen manipuliert wird stellt sich eine zeitunabhängige Dynamik ein, die als stationärer Nichtgleichge- wichtszustand bezeichnet wird. Durch zufällige thermische Einflüsse aus der Umgebung ist der künftige Zustand des Systems unvorhersehbar, sodass dessen Dynamik als stochastischer Prozess modelliert werden kann. Besonderes Augenmerk liegt in dieser Dissertation auf den Strömen, die zwischen den Reservoiren und dem System im stationären Nichtgleichgewicht fließen. Dazu gehören z.B. der Verbrauch oder die Produktion einer chemischen Spezies oder die mit demAnheben eines Gewichtes verknüpfte mechanische Arbeit. Von besonderer ther- modynamischer Bedeutung ist der Strom der mit der gesamten Entropieproduktion des Systems verknüpft ist und damit dessen Nichtgleichgewichtscharakter quantifiziert. Im deutlichem Gegensatz zu Systemen im thermodynamischen Gleichgewicht sind Nicht- gleichgewichtssysteme in der Lage, Ströme aufrechtzuerhalten, die im Mittel von Null verschieden sind. Insbesondere die Rate der Entropieproduktion ist dabei, wie vom zwei- ten Hauptsatz der Thermodynamik gefordert, immer größer als Null. Jedoch unterliegen auch diese Ströme dem thermischen Einfluss der Umgebung, sodass ihre zeitliche Ent- wicklung von Fluktuationen überlagert wird. Das bedeutet, dass auf kurzen Zeitskalen Ströme von ihren Mittelwerten abweichen können und sogar die Entropieproduktion negativ werden kann. Das Ziel dieser Arbeit ist, eine umfassende Charakterisierung der Statistik von Strom- fluktuationen zu entwickeln. Die exakte Berechnung dieser Statistik ist zwar möglich, hängt aber ab von allen mikroskopischen Details des Systems und den von den Reservoi- ren aufgebrachten treibenden Kräften. Da derart detaillierte Informationen in der Praxis nicht verfügbar und auch nicht von Interesse sind, setzen wir uns zum Ziel, Schran- ken an die Statistik von Stromfluktuationen herzuleiten, die von wenigen, möglichst aussagekräftigen thermodynamischen Eigenschaften des Systems abhängen. Ausgangspunkt für diese Arbeit is eine bemerkenswerte Ungleichung, die als “thermo- 9 Zusammenfassung dynamische Unschärferelation” bekannt ist [A.C. Barato und U. Seifert, Phys. Rev. Lett. 114, 158101 (2015)]. In diese Relation gehen zwei wichtige Kenngrößen ein. Die eine ist eine statistische Kenngröße, nämlich die Unschärfe des Stromes, welche sich aus der Amplitude von Fluktuationen im Vergleich zum Mittelwert ergibt. Zweitens betrach- ten wir die mittlere Entropieproduktionsrate als thermodynamische Kenngröße. Das Produkt dieser zwei Größen muss immer größer als zwei sein, was eine fundamentale Abwägung zwischen Präzision und den thermodynamischen Kosten für einen Nicht- gleichgewichtsprozess zum Ausdruck bringt. Diese Abwägung gilt für jeden Strom und für die große Klasse an Systemen die sich durch Markovprozesse beschreiben lassen. Wir stellen dieses Ergebnis in einen breiteren mathematischen Kontext, indem wir es als Konsequenz aus einer ebenso universellen Schranke an das gesamte Spektrum der Stromfluktuationen herleiten. Dazu verwenden wir als mathematisches Werkzeug die sogenannte “large deviation theory”, welche die Statistik großer Abweichungen vomErwartungswert einer geeigneten Zufallsgröße beschreibt. Diese Herangehensweise erlaubt es, etliche Verfeinerungen und Verallgemeinerungen dieser Schranke herzuleiten und liefert neue, ergänzende Schranken an Stromfluktuationen. Kapitel 2: Stochastische Thermodynamik in diskreten Zustandsräumen. In die- sem Kapitel erläutern wir die Formulierung einer thermodynamisch konsistenten sto- chastischen Dynamik in einem diskreten Zustandsraum. Dessen Herleitung beginnt mit der Identifikation von Mesozuständen in einem geschlossenen Gesamtsystem, welches die Reservoire beinhaltet. Innerhalb dieser Mesozustände erreicht das System ein loka- les Gleichgewicht, was dazu führt, dass die Übergänge zwischen den Mesozuständen einemMarkovprozess entsprechen. Da das Gesamtsystem langfristig ein thermodynami- sches Gleichgewicht erreichen muss, können wir detailliertes Gleichgewicht (“detailed balance”) für die Übergangsraten zwischen den Mesozuständen voraussetzen. Bei der Einschränkung auf ein kleines offenes System im Kontakt mit Reservoiren, kann diese Bedingung als ein lokales detailliertes Gleichgewicht ausgedrückt werden. Maßgeb- lich für diese Relation sind dann die Entropien der Mesozustände des kleinen Systems (bzw. deren freie Energien in einem isothermen Kontext) und die mit einem Übergang einhergehenden Änderungen in den Reservoiren. Die resultierende stochastische Dyna- mik beschreibt den stationären Nichtgleichgewichtszustand, welchen das kleine System erreicht. Den Ideen der stochastischen Thermodynamik folgend definieren wir dann thermodynamische Observablen wie Wärme, Arbeit uns insbesondere die Entropiepro- duktion entlang einzelner Trajektorien des Systems. Diese Observablen erweisen sich als Zeit-additiv und als antisymmetrisch unter Zeitumkehr, was sie als Ströme im oben ein- geführten Sinn identifiziert. Einige dieser Observablen erfüllen Fluktuationsrelationen, die bislang bekanntesten universellen Eigenschaften der Statistik von Stromfluktuatio- nen. Zum Abschluss werden die in diesem Kapitel eingeführten Konzepte anhand einer Auswahl einfacher, unizyklischer Modelle für molekulare Motoren illustriert. 10 Zusammenfassung Kapitel 3: Fluktuationen in stationären Markovprozessen. In diesem Kapitel erläu- tern wir die mathematischen Werkzeuge die nötig sind um Fluktuationen in stationären Markovprozessen zu beschreiben.Wir beginnenmit der allgemeinenDefinition von Zeit- additiven Observablen und definieren die kumulantenerzeugende Funktion, welche die Fluktuationen solcher Observablen charakterisiert. Im Grenzfall großer Zeitintervalle wird die Skalierung der kumulantenerzeugenden Funktion durch den größten Eigenwert einer modifizierten Ratenmatrix festgelegt. Kumulanten niedriger Ordnung der Fluk- tuationen, insbesondere deren Mittelwert und die Varianz, können jedoch auch ohne explizite Kenntnis dieses Eigenwerts berechnet werden. Hierfür stellen wir zwei ver- schiedene Methoden vor. Die “large deviation theory”, welche sowohl typische als auch extreme Fluktuationen von Strömen charakterisiert, stellt eine weitere Möglichkeit zur mathematischen Be- schreibung dar. Wir stellen kurz die Grundkonzepte dieser Theorie vor, deren zentrales Element die sogenannte Ratenfunktion ist. Sie erfasst den exponentiellen Zerfall der Wahrscheinlichkeit untypischer Fluktuationen und steht in engem Zusammenhang mit der kumulantenerzeugenden Funktion. Als wesentliche Methode für später Beweise von Schranken an Stromfluktuationen, stellen wir die “level 2.5 large deviation theory” vor, die die gemeinsamen Fluktuationen empirischer Ströme und Dichten in Markovnetzwer- ken beschreibt. Zum Abschluss zeigen wir durch einen geeigneten Kontinuumslimes, dass der zunächst für diskrete Zustände entwickelte Formalismus sich auch auf über- dämpfte Langevindynamik in einem kontinuierlichen Zustandsraum anwenden lässt. Kapitel 4: Dissipationsabhängige Schranke an Stromfluktuationen. In diesem Ka- pitel leiten wir eine dissipationsabhängige Schranke an Stromfluktuationen her, welche unser am vielseitigsten anwendbares Ergebnis darstellt. Es handelt sich dabei um eine parabolische obere Schranke an die Ratenfunktion die einen beliebigen fluktuierenden Strom beschreibt bzw. um eine parabolische untere Schranke an die entsprechende ku- mulantenerzeugende Funktion. Entscheidend ist dabei, dass die Entropieproduktionsrate des Gesamtsystems als einziger Parameter in diese Schranke eingeht. Nahe dem ther- modynamischen Gleichgewicht kann diese Schranke aus der positiven Definitheit der Onsagermatrix abgeleitet werden. Fernab vom Gleichgewicht kann die Schranke in- tuitiv anhand von unizyklischen Netzwerken erklärt werden, den multizyklischen Fall illustrieren wir numerisch. Bei der Auswertung für typische Fluktuationen geht unsere allgemeine Schranke über in die thermodynamische Unschärferelation, deren wichtigs- ten Implikationen wir kurz diskutieren. Abschließend illustrieren wir die Gültigkeit der dissipationsabhängigen Schranke anhand analytischer Lösungen für eindimensionale getriebene diffusive Transportmodelle. Kapitel 5: Anwendung: Schranken an die Effizienz molekularer Motoren. Dieses Kapitel ist als Einschub gedacht, in dem eine wichtige Implikation der thermodyna- 11 Zusammenfassung mischen Unschärferelation für die Effizienz molekulare Motoren hergeleitet wird. In dieser Anwendung ist der relevante Strom die von einem Motor zurückgelegt Strecke, die man typischerweise experimentell beobachten kann. Die Messung der Fluktuatio- nen in diesem Strom ermöglicht eine Abschätzung der Dissipationsrate, die wiederum dazu verwendet werden kann eine obere Schranke an die thermodynamische Effizienz des Motors anzugeben. Diese Schranke gilt für beliebige thermodynamisch konsistente Modelle für molekulare Motoren und benötigt ausschließlich experimentell leicht zu- gängliche Messgrößen. Eine ähnliche Schranke erhalten wir für die Stokes-Effizienz, die einen Motor ohne externe mechanische Last charakterisiert. Kapitel 6: Affinitäts- und topologieabhängige Schranke anStromfluktuationen. Mit- hilfe der durch die level 2.5 large deviation theory bereitgestellten Werkzeuge leiten wir eine Verfeinerung der dissipationsabhängigen Schranke her. Das Aufstellen dieser Schranke benötigt zusätzlich Kenntnis über die Affinitäten die das System antreiben und über die Topologie des zugrunde liegenden Markov-Netzwerkes. Entscheidend ist dabei der Zyklus mit der kleinsten Affinität pro Anzahl der Knoten. Ausgedrückt als unte- re Abschätzung für die kumulantenerzeugende Funktion hat diese Schranke die Form eines Kosinus hyperbolicus, ähnlich der kumulantenerzeugenden Funktion des asym- metrischen Random Walks. Wir illustrieren diese Schranke numerisch. Als wesentliche Voraussetzung für den Beweis entwickeln wir eine Zerlegung der Netzwerkströme in “gleichförmige” Zyklen, die so definiert sind, dass der mit einem einzelnen Zyklus ver- knüpfte Strom entlang jeder Kante parallel zum gesamten Strom ist. Die Affinitäts- und topologieabhängige Schranke kann auch auf Systeme im thermodynamischen Gleich- gewicht angewandt werden, wobei die Anzahl der Zustände im gesamten Netzwerk die Entropieproduktion als entscheidende Kenngröße ersetzt. Kapitel 7: Aktivitätsabhängige und asymptotische Schranke. In diesem Kapitel führen wir weitere Schranken an Stromfluktuationen ein, die von anderen Kenngrö- ßen des Systems als dessen Entropieproduktionsrate abhängen. Zunächst leiten wir eine Schranke her die nur von der Aktivität des Systems abhängt, d.h. von der Anzahl der Übergänge pro Zeiteinheit. Ausgedrückt als Schranke an die kumulantenerzeugende Funktion hat diese Schranke die Form einer Exponentialfunktion, die die tatsächliche Funktion an einer Seite berührt. Mithilfe des Fluktuationstheorems für Ströme kann diese Schranke symmetrisiert werden, sodass sie für beide Seiten der kumulantener- zeugenden Funktion relevant wird. Diese Schranke ist fernab vom thermodynamischen Gleichgewicht oftmals stärker als die dissipationsabhängige Schranke. Ihr Beweis ba- siert auf einem einfachen linearen Ansatz für den empirischen Fluss in der level 2.5 Ratenfunktion. Das zweite Hauptergebnis in diesem Kapitel ist eine Schranke die eher schwach ist für typische Fluktuationen, die aber asymptotisch exakt wird für extreme Fluktuationen 12 Zusammenfassung und damit in einem Bereich in dem die anderen Schranken typischerweise eher schwach sind. Diese Schranke benötigt im Wesentlichen Information über die Topologie des dem Modell zugrunde liegenden Netzwerks, in dem derjenige Zyklus identifiziert wird der für eine bestimmte extreme Fluktuation am relevantesten ist. Zusammen mit den anderen Schranken erhält man so eine umfassende Charakterisierung des Spektrums der Stromfluktuationen. Kapitel 8: Verallgemeinerung für endliche Zeitintervalle. Alle so weit diskutierten Schranken beziehen sich auf Stromfluktuationen die im Grenzfall langer Zeitintervalle auftreten. Die stochastische Thermodynamik hat sich jedoch besonders dann als nütz- lich herausgestellt, wenn relevante Ergebnisse bereits auf endlichen Zeitskalen gewonnen werden können. In diesem Kapitel stellen wir eine Verallgemeinerung der dissipations- abhängigen Schranke für endliche Zeitintervalle vor, die wiederum auch die thermodyna- mische Unschärferelation für endliche Zeiten verallgemeinert. Diese Verallgemeinerung illustrieren wir anhand experimenteller Daten für die verrichtete Arbeit an einemKolloid das sich in einer stochastisch zwischen zwei Positionen springenden optischen Falle be- findet. Überraschenderweise kann dieMessung für endliche Zeitintervalle sogar stärkere Schranken an die Entropieproduktionsrate liefern als der entsprechende Grenzfall langer Zeiten. Wir besprechen eine Beweismethode für diese Verallgemeinerung der thermodyna- mischen Unschärferelation, die auf der level 2.5 large deviation theory für ein Ensemble von unabhängigen Kopien des Systems basiert. Dieser Beweis gilt zunächst für Anfangs- bedingungen die der stationären Verteilung entstammen. Wir entwickeln eine weitere Verallgemeinerung, die sich auch auf transiente Prozesse in endlichen Zeitintervallen mit beliebigen Anfangsverteilungen anwenden lässt. 13 Summary Living systems, as well as useful artificial machines, operate under non-equilibrium conditions. This means that they are in contact with several reservoirs that are not in mutual thermodynamic equilibrium. These reservoirs provide resources such as food or fuel, or act as a thermal environment absorbing heat. Provided that the system under consideration is sufficiently small compared to the reservoirs, the state of the reservoirs will change negligibly on relevant time scales. If additionally the system is not manipulated externally, its dynamics becomes time-invariant, which is called a non-equilibrium steady state (NESS). Due to the thermal influence from its environment, the state of the system becomes erratic, which allows us to model its dynamics as a stochastic process, for which one can define several thermodynamic observables. Of particular interest throughout this Thesis are the input- and output-currents associated with a NESS. Examples for such currents include the consumption or production of a specific chemical species or thework associated with lifting a weight. A current of particular thermodynamic importance is the production of entropy in the total system, which quantifies its non-equilibrium character. In stark contrast to equilibrium systems, non-equilibrium systems are capable of maintaining non-zero average currents. In particular, the rate of entropy production is, due to the second law of thermodynamics, always greater than zero on average. However, again due to the thermal influence from the environment, the temporal evolution of these currents is superimposed by fluctuations. This means, that on short time scales, currents can deviate from the average intensity, and the entropy production can even become negative. The main objective of the work documented in this Thesis is to provide a comprehen- sive characterization of the statistics of current fluctuations. While an exact calculation of these statistical properties is possible, the results typically depend on all microscopic details of the system and on the driving forces associated with the reservoirs. Since such detailed information is practically neither available nor relevant, we focus on the derivation of bounds on the statistics of current fluctuations, which ideally depend on only a few thermodynamic properties of the system. Starting point for our work is a prominent inequality known as the “thermodynamic uncertainty relation” [A.C. Barato and U. Seifert, Phys. Rev. Lett. 114, 158101 (2015)]. It considers the uncertainty of a current, comparing the amplitude of its fluctuations to its mean, as a statisticalmeasure and on the other hand the average rate of entropy production as a thermodynamic measure. The product of these two key quantities must always be 15 Summary greater than two, expressing a trade-off between precision and the thermodynamic cost for a non-equilibrium process. It holds for any current and for the huge class of systems that can be described in terms of Markov processes. We put this relation in a wider mathematical context, employing large deviation theory to derive it as a result of an equally general bound on the whole spectrum of current fluctuations. Our formalism allows for several refinements and generalizations of that bound and yields complementary, novel bounds on current fluctuations. Chapter 2: Stochastic thermodynamics on discrete state spaces. We review a for- malism for thermodynamically consistent stochastic dynamics on a discrete set of states. Its derivation starts with an identification of mesostates in the total, closed system in- cluding the reservoirs. Within these mesostates, the system equilibrates locally, such that the transitions between mesostates become Markovian. Since the total system must ulti- mately reach thermodynamic equilibrium, we can stipulate a detailed balance condition for the transition rates between mesostates. Focusing on a small system in contact with reservoirs, this detailed balance condition gets transformed to a local detailed balance condition involving the entropies of the mesostates of the small systems (or their free energies in an isothermal setting) and the changes in the reservoirs upon a transition. The resulting stochastic dynamics describes the NESS associated with the small system. Fol- lowing the ideas of stochastic thermodynamics, we define thermodynamic observables, such as heat, work and in particular the entropy production along individual trajectories of the system. These observables turn out to be time-additive and antisymmetric under time reversal, which qualifies them as currents in the above defined sense. Fluctuation relations, as the most prominent constraints on current fluctuations so far, are briefly reviewed. Finally, we illustrate the concepts introduced in this Chapter for a selection of simple, unicyclic models for molecular motors. Chapter 3: Fluctuations in stationaryMarkov processes. In thisChapter, we present the mathematical toolbox that is necessary for the description of fluctuations in station- ary Markov processes. We start with a general definition of time-additive observables and define the cumulant generating function that characterizes the fluctuations of such an observable. In the long-time limit, the scaling of the cumulant generating function is given by the largest eigenvalue of a modified transition matrix. Low order cumulants of the fluctuations, in particular their mean and variance, can be calculated without explicit knowledge of this eigenvalue. We present two different methods for carrying out this calculation. Large deviation theory characterizes both typical and extreme current fluctuations, which is complementary to the description via the scaled cumulant generating function. We briefly review the basic concepts of this theory and work out the connection to the algebraic approach using generating functions. A key element of this theory is the 16 Summary so-called large deviation function, or rate function, which captures the exponential decay of the probability of unlikely fluctuations. As the central method for proving bounds on current fluctuations later on, we re-derive the level 2.5 large deviation function, which describes the joint fluctuations of empirical currents and densities in a Markovian net- work. Finally, we show through a suitable continuum limit that the formalism developed so far for discrete state-spaces applies also to overdamped Langevin dynamics on a continuous state space. Chapter 4: Dissipation-dependent bound on current fluctuations. In this Chapter we derive a dissipation-dependent bound on current fluctuations as our most universally applicable result. It can be formulated as a quadratic upper bound on the large deviation function or, equivalently, as a quadratic lower bound on the generating function associ- ated with any current. Crucially, the only parameter entering in this bound is the average rate of entropy production in the total system. In linear response, this bound can be derived from the positive definiteness of the Onsager matrix. Beyond linear response, we give an intuitive explanation for the bound for unicyclic networks and illustrate the general, multicyclic case numerically. When evaluated for typical current fluctuations, our general bound implies the thermodynamic uncertainty relation, whose major impli- cations we briefly review. Finally, we illustrate the validity of the dissipation-dependent bound for analytic solutions for one-dimensional driven diffusive transport-models. Chapter 5: Application: Bounds on the efficiency ofmolecularmotors. ThisChap- ter is intended as an intermezzo, deriving an important implication of the thermodynamic uncertainty relation for the efficiency of molecular motors. The current of interest in this example is the displacement of a motor, which can typically be observed experi- mentally. The measurements of its fluctuations allows for an estimate of the dissipation rate, which can, in turn, be used to derive an upper bound on the thermodynamic effi- ciency of the motor. This bound holds for any thermodynamically consistent model for a molecular motor and uses only experimentally accessible quantities. A similar bound can be derived for the Stokes efficiency that characterizes a motor in the absence of a mechanical load. We illustrate the bound numerically for a motor-cargo complex and for experimental data for kinesin. Chapter 6: Affinity- and topology-dependent bound on current fluctuations. We employ the toolbox provided by level 2.5 large deviation theory to derive a refinement of the dissipation-dependent bound. It requires additional knowledge of the driving affinities and the network topology of the Markov model underlying the system of interest. In particular, the bound is characterized by the cycle with the minimal affinity per length. Formulated as a bound on the generating function, this bound has the shape of a hyperbolic cosine, akin to the generating function for the asymmetric random walk, 17 Summary which we illustrate numerically. As an essential ingredient to the proof of the bound, we derive a decomposition of the network currents into what we call uniform cycles. Along every edge of such a cycle, the cycle current is parallel to the total current. The affinity- and topology-dependent bound can also be applied to current fluctuations in equilibrium systems, where the total number of states enters instead of the entropy production rate. Chapter 7: Activity-dependent and asymptotic bound. In this Chapter, we introduce further bounds on current fluctuations, which depend on other characteristic features of the system than its rate of entropy production. First, we derive a bound that depends only on the activity of the system, i.e., the total number of transitions per unit time. Expressed as a bound on the generating function, this bound has the shape of an exponential function that touches the actual generating function at one flank. Using the fluctuation theorem for currents, this bound can be symmetrized so that it becomes relevant for both flanks of the generating function. This bound is often stronger than the dissipation-dependent one in systems that are driven far away from equilibrium. The proof of the activity-dependent bound is based on a simple, linear ansatz for the empirical flow entering the level 2.5. large deviation function. The second main result of this Chapter is a bound that is rather loose for the typical current fluctuations, but which becomes asymptotically tight for extreme fluctuations, where the other bounds are typically weak. This bound essentially requires information about the topology of the underlying network of states, from which one identifies the cycle that is most relevant for producing a specific extreme fluctuation. Together with the other bounds one obtains a comprehensive characterization of the spectrum of current fluctuations. Chapter 8: Finite-time generalizations. The bounds that have been discussed so far characterize the fluctuations that occur in the limit of long time intervals. However, stochastic thermodynamics usually provesmost useful in cases where relevant results can be obtained on finite time-scales. In this Chapter, we present a finite time-generalization of the dissipation-dependent bound, which likewise generalizes the thermodynamic uncertainty relation to finite times. This generalization is illustrated using experimental data for the work performed on a colloid in a stochastically switching optical trap. Surprisingly, a measurement on finite time intervals can yield stronger bounds on the rate of entropy production than the corresponding long-time limit. Using level 2.5 large-deviation theory for an ensemble of independent copies of the system, we review a rigorous proof for the finite-time uncertainty relation in its original formulation applying to initial conditions sampled from the steady-state distribution. Moreover, this method allows us to derive a further generalization to transient finite- time processes with initial conditions sampled from an arbitrary initial distribution. 18 1 Introduction The statistical mechanics for systems in thermodynamic equilibrium, as developed in the 19th century mainly by Gibbs, Maxwell and Boltzmann, is arguably one of the most successful concepts in physics. It provides universally consistent formal definitions for such intuitive quantities as temperature or pressure, which are based on the microscopic mechanics of a system. In doing so, one needs to make only minimal assumptions about the microscopic mechanics, namely that it conserves energy and that it leads to every microscopic state of the system being equally likely. Thus, statistical mechanics can be applied to virtually every closed many-particle system, including large systems comprising a heat reservoir and an open system of interest. A distinctive feature of a system in thermodynamic equilibrium is its time-reversal symmetry: Any realization of a dynamic process the system might undergo must occur equally likely in a time-reversed fashion. However, this property is in conflict to what one would expect from a useful machine or a living organism. To name a few examples, an engine burns fuel to perform mechanical work, a fridge uses electricity to cool down its content, and a living organism processes food in order to maintain a multitude of life- supporting cycles. For all of these examples, the reverse process would be impossible or at least undesirable. Indeed, all of these systems operate under non-equilibrium conditions. This means that the total system, including the thermal environment and any fuel-, electricity-, or food-reservoirs, has not yet reached thermodynamic equilibrium. In this Thesis, we focus on autonomous systems that show persistent non-equilibrium behavior. Such a persistence requires the system proper to be small compared to the reservoirs, such that the supply of resources does not vary on time-scales relevant for the observation. If this condition is met, the system will reach what is called a non- equilibrium steady state (NESS). How can one quantify the non-equilibrium character of such a NESS? While the macroscopic state of an equilibrium system is fully characterized by a few thermody- namic observables such as temperature, entropy or internal energy, the situation is far less clear in non-equilibrium. Two important characteristics will be recurringly con- sidered in this Thesis. One of these characteristics is a thermodynamic one, namely the production of entropy. It describes the effect of the consumption of resources and of dissipation into the environment transducing the total system from a prepared initial state with low entropy to the final equilibrium state with high entropy. In the long run, the entropy saturates towards its final equilibrium value, but on the time scale relevant for the small system, one observes a constant rate of increase of the entropy. Due to 19 Introduction the second law of thermodynamics, this rate must always be positive. In general, one expects that a high rate of entropy production corresponds to a strongly driven system operating far away from equilibrium. The second characteristic is a statistical one and is less concerned with the driving of the system than with its output. It is based on observables, called “currents”, that are antisymmetric under time reversal. This means, that when a movie of the system is played backwards, the value of such an observable reverses its sign. For the example of a molecular motor walking along a filament, the velocity would qualify as a current observable, which can be integrated to yield the traveled distance. In equilibrium, the time-reversal symmetry stipulates that any current must be zero on average. Nonetheless, the current exhibits unbiased fluctuations. In the example, a non-driven molecular motor would perform a randomwalk taking equally likely forward- and backward-steps, so that it remains on average at the same position. Hence, we can identify a non-vanishing average current as a further signature of a non-equilibrium system. Accordingly, a bias in the motion of a molecular motor is a clear indication of the presence of a driving mechanism. However, the average value of the current has an arbitrary dimension and is therefore not suitable for a universal quantification of the non-equilibrium character. We will therefore compare the average current to the amplitude of its fluctuations, thus defining a dimensionless measure of precision or, as its inverse, the uncertainty. For driven systems very close to equilibrium any current must be rather small on average, such that the current fluctuations, which are present already in equilibrium, prevail. This leads to low precision or, equivalently, high uncertainty. In contrast, strongly driven systems can be expected to produce substantial currents whose averages exceed the respective fluctuations, yielding high precision. These considerations suggest that precision in a current generally requires a system to be driven sufficiently far from equilibrium, which should be reflected in its entropy production as well. Indeed, in a seminal paper [1], Barato and Seifert have reported an inequality that es- tablishes a universal relation between the uncertainty in any current of a non-equilibrium system and the total entropy production. This relation, which has been dubbed the ther- modynamic uncertainty relation, can be stated as ε2C ≥ 2. The squared uncertainty ε2 is defined via the variance and mean of a time-integrated current, and C is the on average produced entropy in units of the Boltzmann constant, characterizing the thermodynamic “cost” for the driving. The thermodynamic uncertainty relation can be interpreted at least in two ways. The first one is from the perspective of a system that needs to achieve a certain degree of precision. For example, the distance traveled by the molecular motor shown in Fig. 1.1 in a fixed time interval fluctuates around an average value. Smaller fluctuations require a higher consumption of fuel that comes with an increased rate of dissipation. As another 20 Introduction START +- heat ε Figure 1.1: A molecular motor walking along a filament. The chemical driving due to the hydrolysis of ATP leads to a bias in the direction of motion and ultimately to dissipation of heat. The motor is repeatedly let run from a fixed starting position and for a fixed time interval. In each run it reaches a slightly different final position, as indicated by the overlaid snapshots and characterized by the uncertainty ε. example, circadian clocks of organisms have evolutionarily adapted to anticipate the end of nighttime even under deprivation of any external stimuli [2]. According to the uncertainty relation, the organism has to expend the more energy, the more precise this prediction is required to be. Second, the thermodynamic uncertainty relation is also promising from an experimen- tal point of view. With the advent of new experimental techniques for the manipulation and observation of small non-equilibrium systems, measuring the fluctuations of an accessible degree of freedom has become a fairly simple procedure [3]. However, measuring the rate of dissipation that characterizes such a system remains a major chal- lenge [4, 5], especially when chemical reactions are involved. The uncertainty relation yields lower bounds on the rate of entropy production based on the visible fluctuations of any current observable, without requiring any knowledge of the microscopic details of the system. Such bounds will potentially prove useful for testing hypotheses one may have for the yet unknown functioning of, e.g., biomolecular systems. In the work documented in this Thesis, we put the thermodynamic uncertainty relation on a broad theoretical basis, integrating it into the framework of stochastic thermody- namics [6]. In doing so, we derive various generalizations, refinements and implications of the uncertainty relation. Moreover, we explore which other properties of a system in a steady state provide bounds on its current fluctuations. As key result, we introduce a universal, dissipation-dependent bound on the whole spectrum of current fluctuations [7]. While the mean and variance of a fluctuating current, forming the uncertainty ε, characterize only typical fluctuations, we will use a description that applies also to untypical, extreme fluctuations. It is based on the so-called large deviation function that captures the exponential contributions to the distribution of current fluctuations. This function turns out to be bounded by a quadratic function, in which the total dissipation enters as the only parameter. This bound implies 21 Introduction the thermodynamic uncertainty relation as the special case of typical fluctuations. This augmented level of description in terms of large deviation theory has indeed proven crucial for a general proof of the uncertainty relation [8]. The dissipation-dependent bound is particularly strong for weak driving of the system close to equilibrium. Far away from equilibrium the bound is still valid, but presents a rather loose estimate of the actual current fluctuations. Quite in general, more details about a system will lead to stronger bounds on its fluctuations. For instance, we present a major refinement of the dissipation-dependent bound that depends additionally on the topology of the Markovian network that models the system and on the strengths of the individual driving forces. Two other bounds that are unrelated to dissipation complement the dissipation-dependent bound. One of them depends on the overall rate of transitions in the system and is often stronger far from equilibrium, the other one depends essentially on the topology of the network and becomes strong for the asymptotics of extremely untypical fluctuations. The initial formulation of the thermodynamic uncertainty relation considers exclu- sively current fluctuations that occur on large time scales. On these time scales cor- relations vanish, such that the central limit theorem requires the distribution of typical current fluctuations to be Gaussian. In contrast, current fluctuations on finite time scales are distributed in a rather different, non-Gaussian fashion. Despite these differences, the uncertainty relation along with the bounds on the fluctuation spectra turn out to be valid on finite time scales. For experimental applications, this insight not only reduces the amount of required data, it can also lead to stronger bounds on the dissipation rate. All these bounds taken together provide a comprehensive picture of how only a few crucial thermodynamic properties coin the extend and shape of current fluctuations of a system in a NESS. 22 2 Stochastic thermodynamics on discrete state spaces The past few decades have seen a rise in new experimental techniques to observe and manipulate micro- and mesoscopically small systems like colloids, biomolecules or nanoelectronic devices [3]. For such systems, the ubiquitous thermal noise plays a non-negligible role. Stochastic thermodynamics provides a theoretical framework for the description of such systems, that merges concepts from classical thermodynamics, non-equilibrium statistical mechanics and stochastic processes [6, 9–11]. In this Chapter, we review the concepts of stochastic thermodynamics on the basis of a universally applicable description of thermal systems using Markovian processes on discrete state spaces. 2.1 Markovian dynamics in closed systems On the most detailed level of description, we consider a closed system with microstates ξ that undergo Hamiltonian dynamics with some Hamiltonian function H (ξ) and a total energy Etot. Here, the microstate ξ represents a high dimensional vector that contains the configuration and momenta of all N particles within the system. For example, for a biomolecular system, this would comprise all atoms of the biomolecules, of the solvent, and of dissolved reactants. In equilibrium, the microcanonical distributions for these microstates is peq(ξ) = δ(H (ξ) − Etot)/(h3NΩ(Etot)), (2.1) with the normalization constant Ω(Etot) ≡ ∫ dξ h3N δ(H (ξ) − Etot) ≡ eS/δE, (2.2) where the integration is carried out over the whole phase space. The quantities h and δE are arbitrarily chosen constants with dimensions of an action and an energy, respectively, which allows for a dimensionless definition of themicrocanonical entropy S. Throughout this Thesis, we set Boltzmann’s constant kB = 1, such that the conventional definition of entropy becomes equivalent to S. 23 2 Stochastic thermodynamics on discrete state spaces Since a full solution of the microscopic equations of motion for ξ is neither feasible nor desirable, we aim for a coarse-grained description of the system and its dynamics by grouping the microstates into mesostates x. The mapping has to be such that every microstates belongs to exactly one mesostate, moreover, we group microstates with the same spacial configuration of particles but different momenta into the same mesostate. This way, in equilibrium there is no net flow between two mesostates x and y, i.e., p[ξ (0) ∈ x, ξ (t) ∈ y] = p[ξ (0) ∈ y, ξ (t) ∈ x] (2.3) for arbitrary times t, as can be checked via time reversal andLiouville’s theorem. Notably, Eq. (2.3) already implies micro-reversibility, stating that whenever a transition from x to y is possible, also the reverse transition from y to x must be possible. The probability to observe the microstate x in equilibrium is given by peqx = ∫ ξ∈x dξ peq(ξ) = ∫ ξ∈x dξ h3N δ(H (ξ) − Etot)/Ω(Etot) ≡ eSx−S, (2.4) which defines the constrained entropies Sx . The sequence of mesostates visited by a trajectory ξ (t) defines a coarse-grained trajectory x(t). Due to the large number of microstates underlying one mesostate the dynamics of x(t) becomes effectively stochastic even if, for a closed system, the dynamics of ξ (t) is a fully deterministic solution of Hamilton’s equations. In order to describe the dynamics of the coarse grained trajectory using stochastic methods, a further assumption is necessary. Since the knowledge of the sequence of visited mesostates allows one in principle to narrow down the set of possible microstates ξ (t), the probability of transitions from one mesostate x to another one y depends a priori on the whole history of previously visited mesostates. However, we will assume that within a mesostate the distribution of microstates and states of the heat bath relax into a local equilibrium distribution on a time scale that is much smaller than the typical time scale for transitions between mesostates. Such an assumption is justified if the mesostates are chosen such that they are separated by sufficiently high potential or entropic barriers. As a consequence, we canmodel the dynamics of x(t) as a “memoryless” stochastic process, i.e., conditional probabilities for trajectories with x(ti) = xi at times t1 > t2 > . . . satisfy the Markov property [12] p(x1, t1 |x2, t2; x3, t3; . . .) = p(x1, t1 |x2, t2) = p(x1, t1 − t2 |x2, 0), (2.5) where the second equality reflects the time-homogeneity of the process, which is ulti- mately a consequence of the explicit time-independence of the Hamiltonian H (ξ). The propagator p(y, t |x, 0) thus fully characterizes the dynamics of the mesostates, however, this characterization is not yet minimal. Since the propagator for a large time interval can be calculated from the propagators for smaller time intervals via the consistency 24 2.1 Markovian dynamics in closed systems condition, known as Chapman-Kolmogorov equation, p(y, t |x, 0) = ∑ x′ p(x, t |x′, t′)p(x′, t′|y, 0) (2.6) with arbitrary 0 < t′ < t, we can focus on the transition rates kxy that are defined from the propagator in the small t limit as kxy ≡ lim t→0 1 t p(y, t |x, 0). (2.7) Throughout this Thesis, we denote in the index of transition rates first the old state and then the new state. With this definition, Eq. (2.6) can be written in a differential form, leading to the so-called Master equation ∂tpx (t) = ∑ y [py (t)kyx − px (t)kxy] (2.8) for the temporal evolution of a probability distribution px (t) on the set of mesostates {x}. Assuming that the full system is ergodic, the dynamics (2.8) transforms for large t any initial distribution into the equilibrium distribution peqx in Eq. (2.4). Its stationarity requires 0 = ∑ y [peqy kyx − peqx kxy]. (2.9) Considering Eq. (2.3) for small times t, we see that every term of the sum in Eq. (2.9) vanishes individually, which can be expressed as the detailed balance condition peqy kyx = p eq x kxy . (2.10) Inserting the microcanonical distribution (2.4) for peqx leads to kxy kyx = eSy−Sx . (2.11) Thus, the mere existence of an equilibrium distribution, as postulated by statistical mechanics, yields a thermodynamic constraint on the transition rates kxy, that is itself independent of the actual probability distribution. The main concern of this Thesis is on systems that are not in an equilibrium state. For a closed system as described above, a non-equilibrium state is any distribution px that does not satisfy Eq. (2.10), leading to non-zero probability currents jxy ≡ pxkxy − pykyx, (2.12) 25 2 Stochastic thermodynamics on discrete state spaces while the transition rates still satisfy the constraint (2.11). Such a distribution may arise in various ways. First, the distribution may represent the knowledge of an external observer who has performed a measurement on the system in equilibrium. Without further measurements ormanipulations of the system, the distribution ofmesostates conditioned on the outcome of the measurement is governed by the master equation (2.8) and relaxes back into the equilibrium distribution. Second, the system may be prepared in a non-equilibrium state. For example, such a preparation can be achieved by letting the system equilibrate for some Hamiltonian function which is then suddenly changed. After this quench, the system is in a transient non-equilibrium state that relaxes towards the Boltzmann distribution corresponding to the new Hamiltonian. Third, the Hamiltonian function H (ξ, λ) with a control parameter λ can be subject to a time-dependent protocol λ(t) that is performed by an external agent. Provided that the decomposition into time-scale separated mesostates is still independent of λ, one obtains time-dependent entropies Sx (λ) and S(λ) in Eq. (2.4) and subsequently time-dependent transition rates in Eq. (2.11) and the master equation (2.8). If the external protocol modulates λ(t) in a time-periodic fashion, the dynamics of the master equation will lead to an equally time-periodic probability distribution px (t), as realized in so-called stochastic pumps, see, e.g., [13, 14]. If the control parameter λ(t) itself undergoes an athermal, time-homogeneous Markov process, as realized experimentally for example in [15, 16], the dynamics in the joint state space of x and λ becomes a non-equilibrium stationary state (NESS), which has a unique stationary distribution with non-vanishing probability currents. 2.2 Markovian dynamics in open systems Without external manipulation, closed systems (or systems in contact with a single heat bath) can only be in transient non-equilibrium states that are bound to relax into equilibrium after some equilibration time. In contrast, open systems are often observed to be in a steady state and yet maintain steady fluxes, which is only possible out of equilibrium. Examples include an electric circuit with a bulb that steadily produces light or biological systems like chloroplasts producing organic molecules. Common to these systems is the contact to various reservoirs providing mechanical work, heat, chemicals, or electrons. Following the ideas put forward in Refs. [17–20], we deduce the (thermo-)dynamics of such an open system by embedding it along with all the reservoirs in a large closed system. The reservoirs (labeled by r) are fully characterized by their entropies Sr , other extensive quantities Xr , and an equation of state that gives the energy Er = Er (Sr,Xr ). Associated with these extensive observables are the temperature Tr ≡ β−1r ≡ ∂Er/∂Sr 26 2.2 Markovian dynamics in open systems and the generalized force Kr ≡ ∂Er/∂Xr . For instance, for a heat and particle reservoir one would have the number of particles as extensive quantity Xr = N and the chemical potential as intensive quantity Kr = µ. A pure work reservoir could be a solid lifted to height Xr = h against the gravitational force mg. Such pure work reservoirs with a single degree of freedom are formally assigned zero entropy and temperature. Likewise, pure heat reservoirs are assigned zero Xr and Kr . For reservoirs that provide more than one type of resource, the quantities Xr and Kr are understood as vectors with implied dot product multiplication. The reservoirs are weakly coupled to the open system, with which they exchange heat, work, particles or other resources. Thus, in the long run, the full system will relax into an equilibrium state characterized by all temperatures Tr and all forces Kr becoming pairwise equal. However, we assume that the reservoirs are sufficiently large, so that that there is a clear time scale separation between this slow relaxation process and the time scale relevant for the open system. In particular, the intensive quantities Tr and Kr can be assumed to be constant on the time scale of interest. Under these conditions, the open system can effectively be described by a dynamics that leads locally to a NESS. Mesostates of the full system can now be identified as x = (i, {Er }, {Xr }), where i denotes a mesostate of the open system alone. The constrained entropy (2.4) of such a state is then composed of an intrinsic entropy Si and the entropy of the reservoirs, Sx = Si + ∑ r Sr (Er,Xr ). (2.13) We assume Markovian transitions with rate kxy between two states x = (i, {Er }, {Xr }) → y = ( j, {Er + ∆Er }, {Xr + ∆Xr }), (2.14) which requires a sufficiently fast local equilibration of the degrees of freedom of both the open system and the reservoirs. The detailed balance condition (2.11) now reads with the entropies (2.13) ln ( kxy/kyx ) = Sj − Si + ∑ r (βr ∆Er − βrKr ∆Xr ). (2.15) Since the reservoirs are assumed to be large, the rates kxy will be independent of the absolute values of Er and Xr during the relevant observation times. This fact allows us to reduce the state space necessary for the description from the mesostates {x} of the full system to only the mesostates {i} of the open system, while still keeping track of the changes occurring in the reservoirs. In this reduced description, the transition rates kxy are denoted as kai j , where the upper index a ≡ ({∆Ear }, {∆Xar }) (2.16) 27 2 Stochastic thermodynamics on discrete state spaces Figure 2.1: Illustration of the network representation of an open system. If multiple transi- tions between two states or within a single state are possible (colored arrows), the corresponding transition rates are labeled with upper indices. denotes the changes in the reservoirs associated with a particular transition. This way, we arrive at a network description of the open system, in which each vertex represents a state i and each directed edge represents a non-zero transition rate kai j , as illustrated in Fig. 2.1. As a consequence of micro-reversibility, each transition with rate kai j comes in pair with a reversed transition with rate k a¯ji, where a¯ ≡ ({−∆Ear }, {−∆Xar }). For this pair of rates, the constraint (2.15) appears as the so-called local detailed balance condition ln ( kai j/k a¯ ji ) = Sj − Si + ∑ r (βr ∆Ear − βrKr ∆Xar ). (2.17) If the changes in the reservoirs are a unique function of the states i and j, we will drop the upper index in the transition rates for notational convenience. A further constraint on the rates in Eq. (2.17) is imposed by the conservation of the energy in the total system, leading for every transition to E j − Ei + ∑ r ∆Ear = 0, (2.18) where Ei is the intrinsic energy associated with state i of the open system. We can now formulate, analogously to Eq. (2.8), the master equation for the distribu- tion pi (t) for the mesostates of the open system only as ∂tpi (t) = ∑ j,a [p j (t)kaji − pi (t)kai j] ≡ ∑ i,a jai j (t), (2.19) which can be regarded as a conservation law for the probability currents jai j (t). Using the Perron-Frobenius theorem [21], one can prove that in the long-time limit, pi (t) tends 28 2.3 Energetics of open systems to a stationary distribution psi . It can be obtained as the unique solution of the linear system of equations 0 = ∑ j,a [psj k a ji − psi kai j] ≡ ∑ i,a ja,si j , (2.20) which can be solved either directly using linear algebra or using a diagram method for cycle fluxes [22]. Characteristically for a NESS, the stationary probability currents ja,si j are typically non-vanishing. What vanishes instead is the sum of all net influx currents jai j into a state j, which can be regarded as an instance of Kirchhoff’s law. For isothermal systems with all reservoir temperatures βr ≡ β, the constraint (2.18) allows one to restate the local detailed balance condition (2.17) as ln ( kai j/k a¯ ji ) = −β(Fj − Fi + ∑ r Kr ∆Xar ), (2.21) where we have defined the intrinsic free energy Fi ≡ Ei − Si/β. In the presence of only a heat reservoir, the stationary distribution solving Eq. (2.20) is given by the Boltzmann distribution peqi ≡ exp (−βFi)/Z (β), (2.22) which is normalized by the canonical partition function Z (β). For this equilibrium distribution, all probability currents vanish. Thus, in order to generate a genuine NESS, at least one further reservoir that is not in equilibrium with the existing one is required. 2.3 Energetics of open systems The basic concept of stochastic thermodynamics is to define thermodynamic quanti- ties like heat, work, entropy, or internal energy not only macroscopically, but also as trajectory-dependent and hence fluctuating quantities [6, 9]. The internal energy of the open system can straightforwardly be defined as a trajectory- dependent observable by evaluating the intrinsic energy Ei along the trajectory. Upon an individual transition i → j of type a, the change of this energy can be expressed with the conservation law (2.18) as ∆Ei j ≡ E j − Ei = − ∑ r ∆Ear = − ∑ r (Tr ∆Sar + Kr ∆Xar ). (2.23) In between the transitions, the internal energy and the reservoir observables stay constant, since we assume that there is no time-dependent external manipulation of the system or the reservoirs. 29 2 Stochastic thermodynamics on discrete state spaces The identification of work and heat in stochastic thermodynamics is highly context- dependent. In systems that are externally manipulated one typically defines a “protocol work” along a microscopic trajectory ξ (t) as the change of the Hamiltonian H (ξ, λ) of the total system for a time-dependent protocol λ(t). For instance, this definition of work is the one occurring in the Jarzynski relation [23], and it can be identified along trajectories of mesostates i(t) through appropriate coarse-graining techniques [24]. In contact with a heat reservoir, heat is identified through the first law of thermodynamics as the difference between the protocol work and the change of internal energy. Since we focus in this Thesis on autonomously operating systems, we refer for a discussion of this definition of work and heat to the literature [6, 9]. In the context of open systems that are time-independently driven by several reser- voirs, we consider both work and heat as energies transferred to or from the reservoirs. Tentatively stating the first law for work and heat associated with a transition i → j of type a as ∆Ei j = −W a −Qa = ∑ r (−W ar −Qar ) (2.24) (with the sign convention of work performed by the system and heat dissipated into the reservoirs), leads to the question of an appropriate splitting of the terms on the right- hand side of Eq. (2.23) to the amounts of work and heat exchanged with the respective reservoirs. For the case of a pure heat reservoir (which does not have a variableXr), it is obvious to attribute the term Qar ≡ Tr ∆Sar fully to the heat. Conversely, for a pure work reservoir consisting of a single mechanical degree of freedom X one readily identifies the mechanical workW ar ≡ Kr ∆Xar . One might be tempted to use the same identification for other types of reservoirs, e.g., for particle reservoirs that unavoidably act as heat reservoirs as well. However, this identification can be misleading. For instance, consider a reservoir providing particles that can have two different states A and B that affect neither the internal energies of the particles nor their interactions. Now suppose that the open system catalyzes a reaction A ⇌ B without changing its own state. Such a reaction step alters the configurational entropy of the reservoir but leaves its energy unaffected. Since this process obviously involves neither heat nor mechanical work, the Tr ∆Sar cannot be the heat. In order to arrive at a consistent and unambiguous definition of heat, we consider each such mixed reservoir as consisting of a finite container of the resource Xr , which is still large enough for the thermodynamic limit to be applicable, and which is embedded in a larger, purely thermal heat bath at equal temperature Tr . Denoting the small part of the reservoir as r′ and the large part as r′′, the total change of energy in the combined reservoir r upon a transition of type a becomes ∆Ear = ∆Ear ′ + ∆Ear ′′ = Tr ∆Sar ′ + Kr ∆Xar + Tr ∆Sar ′′ . (2.25) 30 2.4 Entropy production We now define the heat associated with this transition and this reservoir as the energy Qar ≡ ∆Ear ′′ = Tr ∆Sar ′′ (2.26) that is dissipated into the heat reservoir r′′. Using the free energy Fr ′ of the reservoir r′, we find for its change of entropy in the isothermal environment thermostatted by the reservoir r′′ ∂Sr ′ ∂Xr Tr = − ∂ 2Fr ′ ∂Tr∂Xr = −∂Kr ∂Tr . (2.27) Therefore, the heatQar can be identified without explicit reference to the formal splitting of the reservoir as Qar = ∆Ear + ( Tr ∂Kr ∂Tr − Kr ) ∆Xar . (2.28) For the example considered above with the reservoir containing two species with particle numbersX = (NA, NB), the chemical potentials derived from the configurational entropy are K = (µA, µB) = (T ln NA,T ln NB). A transition does not change the reservoir energy and the term in parentheses in Eq. (2.28) vanishes, therefore the associated heat Qa vanishes as expected. Mechanical work is a well defined concept only for athermal systems with one or a few controllable mechanical degrees of freedom. Using the then obvious definitions W ar = Kr ∆Xar andW a = ∑ r W ar , the first law of thermodynamics for a transition i → j of type a can be stated generally as ∆Ei j + ∑ r ′ ∆Ear ′ = W a −Qa . (2.29) Here, the left-hand side of this equation then subsumes the internal energies of the open system and of all finite parts of combined reservoirs. Occasionally, we will also consider non-mechanical types of work, in particular a “chemical work” that is associated with the change of free energy ∆µ in a particle reservoir. It should be noted, however, that this type of work involves also a change of entropy and should therefore not be interpreted in the sense of the first law (2.29). 2.4 Entropy production For the change of entropy along a trajectory, one can identify three contributions. First, there is the intrinsic entropy Si(t) of the current state i, as introduced in Eq. (2.13). As the intrinsic energy, this entropy is constant while the system remains in a mesostate 31 2 Stochastic thermodynamics on discrete state spaces and changes discontinuously with increment Sj − Si upon a transition i → j. A sec- ond contribution is the stochastic entropy − ln[pi (t)], where pi (t) is a solution of the Master equation (2.19) for a given initial distribution evaluated for the current state i [25]. This entropy is defined such that its ensemble average matches the Gibbs entropy −∑i pi (t) ln pi (t) of the distribution pi. The intrinsic and the stochastic contributions to the entropy can be subsumed as the system entropy Ssys(t) ≡ Si − ln[pi (t)] (2.30) while the system is in state i. The third contribution is the entropy production S(t) in the reservoirs. It changes upon a transition i → j of type a by ∆Sa = ∑ r ∆Sar = ∑ r (βr ∆Ear − βrKr ∆Xar ) = ln ( kai j/k a¯ ji ) − (Sj − Si), (2.31) where we have used the local detailed balance condition (2.17). Provided that there are only pure work and pure heat reservoirs, this change of entropy is connected to the heat Qa = ∑ r ∆Sar /βr exchanged with the reservoirs. The total entropy production follows as the sum Stot(t) ≡ Ssys(t) + S(t) (2.32) of the system entropy and the entropy of the reservoirs. Its change upon a transition i → j of type a at time t is ∆Satot = ∆Sa + ln pi (t) p j (t) = ln pi (t)kai j p j (t)k a¯ji . (2.33) If the probability distribution pi (t) is not yet the stationary distribution, Stot(t) changes also continuously in between transitions due to the contribution from the stochastic entropy. Importantly, the total entropy production can be identified solely through the transition rates and the corresponding solution of the Master equation and is thus inde- pendent of the identification of heat and work. Eq. (2.33) could therefore be regarded as a premise that defines entropy production for any Markovian jump process. This “math- ematical” entropy production becomes the physical one through the thermodynamic consistency imposed by local detailed balance, provided that the Markovian model fol- lows the principles detailed in Secs. 2.1 and 2.2. In particular, at times in between transitions, the reservoirs must not exchange (or “leak”) energy or other resources, ex- cept for the small fluctuations with zero mean of the local equilibrium associated with a mesostate. 32 2.5 Time extensivity and time reversal 2.5 Time extensivity and time reversal Building on the concepts developed above for open systems in contact with reservoirs, we consider a trajectory as a sequence of transitions between mesostates. At times {tn}, the trajectory jumps from state i−n to state i+n with the concomitant changes an = ({∆Eanr }, {∆Xanr }) in the reservoirs. In order to keep the notation simple, we assume that these changes in a transition i → j are a unique functions ∆Ei jr , ∆Xi jr , allowing us to replace the upper index a by i j and a¯ by ji. Under this assumption, trajectories can simply be referred to as functions i(τ), where τ parameterizes the time. The generalization to trajectories in networks with several transition pathways between two states labeled additionally by the sequence of indices an is obvious. In total, we identify two types of observables in stochastic thermodynamics. Observ- ables like the internal energy, intrinsic entropy or stochastic entropy are state functions depending only on the mesostate and possibly explicitly on time. Their change along a stochastic trajectory i(τ) with 0 ≤ τ ≤ t can therefore be stated using only boundary terms, e.g., the change of system entropy is given by ∆Ssys[i(τ)] = Si(t) − Si(0) − ln pi(t) (t)pi(0) (0) . (2.34) In the steady state, these observables fluctuate around an average value and do not scale with time. In contrast, observables like heat, work, reservoir entropy and the exchange of energy or other resources with an individual reservoir are typically time-extensive observables (see also [26]). They depend not only on the boundary but on the whole course of a trajectory, as they are accumulated with every transition. The total (or “integrated”) work along a trajectory is given by the sum over all transitions W [i(τ)] = ∑ n W i + ni − n (2.35) and likewise for the other time-extensive observables defined through their respective changes upon a transition. The average rate of change of such an observable in the steady state is given by⟨ W˙ (t) ⟩ = ∑ i j psi ki jW i j = ∑ i> j jsi jW i j, (2.36) where psi ki j gives the average number of transitions from i to j per unit time. Using the antisymmetry of W i j , this average can also be represented through the stationary currents jsi j and a sum over unique pairs of states (indicated by i > j, such that pairs are not counted twice). Throughout this Thesis, we use the notation ⟨. . .⟩ for averages over the ensemble of trajectories under consideration, in this case the ensemble of possible trajectories in the steady state. 33 2 Stochastic thermodynamics on discrete state spaces The total entropy production (2.32) along a trajectory is given by Stot[i(τ)] = − ln pi(t) (t)pi(0) (0) + ∑ n ln ki+ni−n ki−ni+n . (2.37) Since the contributions from the intrinsic and stochastic entropy are not time-extensive, the rate of entropy production in the steady state is the same for the total entropy and the reservoir entropy, leading to [27] σ = ⟨S˙(t)⟩ = ⟨S˙tot(t)⟩ = ∑ i> j jsi j ln ki j k ji = ∑ i> j jsi j ln psi ki j psj k ji ≥ 0. (2.38) The non-negativity of σ follows from the non-negativity of the individual terms in the sum of the last equality with the stationary current jsi j = p s i ki j − psj k ji. Using Eq. (2.31), the rate of entropy production can be written as a sum over contributions from the individual reservoirs σ = ∑ r βr (⟨E˙r⟩ − Kr⟨X˙r⟩) (2.39) with the rates of energy transfer ⟨E˙r⟩ and the turnover rates of the resources ⟨X˙r⟩. Crucially, the observables we have discussed so far are odd under time-reversal, i.e., the sign of an increment is flipped for the opposite transition, e.g., for the increment of the work, W i j = −W ji must hold. Expressed in a time-integrated fashion, the time antisymmetry reads W [i(τ)] = −W [i˜(τ)], (2.40) where we denote the time-reversed counterpart of a trajectory of length t as i˜(τ) ≡ i(t − τ). (2.41) Analogous relations hold for the heat, the integrated changes in the reservoirs and essentially all intrinsic observables of the state of the open system. The only exception is the stochastic entropy, which has a non-symmetric explicit time dependence when the system is not in the steady state. Three other types of time-extensive observables are conceivable and have indeed been considered is various contexts. First, one can define time-symmetric observables like the “dynamical activity”, “traffic” or “frenecy” that count a number of transitions irrespective of the direction, see, e.g., [28, 29]. A second time-symmetric observable is the time spent in a state or set of states [30, 31]. Third, a “directed current” or “flow” counting transitions between states only in a specific direction is neither odd nor even under time reversal. Despite their technical relevance for Markov processes in general, such types of observables play a subordinate role in stochastic thermodynamics, since they have no physical counterpart in classical thermodynamics. 34 2.6 Fluctuation relations 2.6 Fluctuation relations The central role that entropy production plays throughout stochastic thermodynamics is largely due to its quantification of the concept of irreversibility in a series of prominent results that are known as fluctuation relations [6]. A straightforward derivation of these results for systems in a steady state starts with the consideration of path weights. The statistical weight of a trajectory i(τ) underMarkovian dynamicswith transition rates ki j andwith the initial state from the stationary distribution psi is given by p[i(τ)] = psi(0) *, N∏ n=1 exp [ −(τn − τn−1)ri−n ] ki−ni+n+- exp [ −(t − τN )ri+N ] , (2.42) for a trajectory of length t performing jumps at times {τ1, . . . , τ2}. We use the same notation for the visited states as in Sec. 2.5, set τ0 ≡ 0, and define the exit rate ri ≡ ∑ j ki j . (2.43) This exit rate describes in Eq. (2.42) the exponential decay of the probability to stay in an individual state. The full path weight is made up of the probability of the initial state psi(0) drawn from the steady state and the product of the probabilities for the times spent in states and the transition rates along the path. For the time reversed trajectory i˜(τ) = i(t − τ), the time spent in the individual states is the same as for i(τ). This contribution therefore cancels in the ratio of weights of the forward and time-reversed path, leading to ln p[i(τ)] p[i˜(τ)] = ln ⎡⎢⎢⎢⎢⎣ psi(0) psi(t) N∏ n=1 ki−ni+n ki+ni−n ⎤⎥⎥⎥⎥⎦ = ∆Stot[i(τ)], (2.44) where we identify the total entropy production of the forward trajectory from Eq. (2.37) (for the case of a unique type a of a transition considered here). Hence, the total entropy production of a trajectory quantifies its degree of irreversibility by comparing its likelihood to that of the time-reversed trajectory. With this property, the integral fluctuation theorem for the total entropy production [25] is readily derived as ⟨ e−∆Stot ⟩ = ∑ i(τ) p[i(τ)]e−∆Stot[i(τ)] = ∑ i(τ) p[i(τ)] p[i˜(τ)] p[i(τ)] = ∑ i˜(τ) p[i˜(τ)] = 1. (2.45) Here, we have represented the ensemble average using a path integral, denoted as a sum over trajectories i(τ) weighted with the path weight (2.42). Since forward and 35 2 Stochastic thermodynamics on discrete state spaces time-reversed trajectories must have the same integration measure in the path integral, we can switch to an integration over time-reversed trajectories. The normalization of the path weights finally allows for a derivation of Eq. (2.45) without need for an explicit specification of the measure underlying the path integrals. While we have focused on the case of systems in a steady state, the integral fluctuation theorem (2.45) and variants thereof (e.g., the Jarzynski relation [23]) can be derived along similar lines for time- dependently driven systems and relaxations from non-stationary distributions [6, 25]. Since the exponential is a convex function, one can apply Jensen’s inequality to Eq. (2.45), which yields the positivity of entropy production ⟨∆Stot⟩ = σt ≥ 0 (2.46) as required by the second law of thermodynamics (see also Eq. (2.38)). In this sense, the fluctuation theorem (2.45) can be regarded as a refinement of the second law, replacing the inequality by an equality. Other fluctuation relations relevant for systems in the steady can be derived from a “master fluctuation theorem”, as has been discussed in a more general setting in Ref. [6]. We consider an arbitrary time-antisymmetric and possibly vectorial functional X[i(τ)] = −X[i˜(τ)] of the trajectory and an arbitrary function Γ (X ) thereof. For the average of this function we then derive the relation ⟨Γ (X )⟩ = ∑ i(τ) p[i(τ)] Γ (X[i(τ)]) = ∑ i(τ) p[i˜(τ)] e−∆Stot[i˜(τ)] Γ (−X[i˜(τ)]) = ⟨ e−∆StotΓ (−X ) ⟩ . (2.47) Setting X = ∆Stot and Γ (X ) = δ(X + s) yields the so-called detailed fluctuation theorem p(−∆Stot) p(∆Stot) = e−∆Stot, (2.48) which expresses that the probability to observe a negative entropy production during the time t is exponentially less likely than to observe the corresponding positive entropy production. Similarly, one derives for the joint fluctuations of the entropy production and of a vector of time-antisymmetric observables X the relation p(−∆Stot,−X ) p(∆Stot, X ) = e−∆Stot (2.49) by setting X = (∆Stot, X ) and using a delta-like function Γ (X ). The fluctuation relations we have discussed so far all hold for arbitrary time intervals t. In contrast, fluctuation relations that do not directly involve the total entropy production 36 2.6 Fluctuation relations often hold only asymptotically in the limit of large time intervals. Consider an observable of the type X = ∆(Stot(τ) + ϕi(τ)), which is the change of the total entropy up to a state function ϕi (the notation using ∆ implies the difference of the argument evaluated at τ = t and τ = 0). The general idea is that on long time scales the values of X and ∆Stot can scale with t, whereas the contribution from ϕi remains of order O(1). The most prominent example of such an observable is the entropy production ∆S in the reservoirs (2.33), for which we set ϕi = ln psi − Si. With Γ (X ) = δ(X + x) we then obtain the relation [6] p(−X ) p(X ) = e−X ⟨eϕi (t )−ϕi (0) |X⟩, (2.50) where ⟨. . . |X⟩ denotes an average conditioned on the observation of X . Since this average does not scale with t, the asymptotic fluctuation relation ln p(−X ) p(X ) ≈ −X (2.51) holds in linear order for large t. For X = ∆S, this relation is known as Gallavotti- Cohen symmetry, which was found numerically [32] and has been proven for chaotic dynamics [33] and stochastic dynamics [34, 35]. It should be noted, however, that for values of X which do not scale with t, the relation (2.51) is generally not exact even in the limit t → ∞. For instance, for one-dimensional periodic systems, the sub-exponential contributions lead to a periodic modulation or “fine-structure” in the fluctuation relation [36]. Fluctuation relations of the type (2.51) may also hold for observables that are unrelated to the entropy production [37]. A similar asymptotic fluctuation relation is the fluctuation theorem for currents [38– 41], which follows for a vector of currents X that determine the total entropy production as a linear form ∆Stot =A · X + ∆ϕi (2.52) with a vector of affinities A up to a state function ϕi. Such a linear form can be identified in Eq. (2.31) with the currents X = ({∆Er }, {∆Xr }) denoting the exchange with the reservoirs along a trajectory and the affinitiesA = ({βr }, {−βrKr }). Along the same lines as for Eq. (2.51), one obtains the fluctuation relation ln p(−X ) p(X ) ≈ −A · X (2.53) in linear order for large t. 37 2 Stochastic thermodynamics on discrete state spaces 2.7 Unicyclic molecular motor as paradigmatic example In order to illustrate the concepts of open systems and stochastic thermodynamics intro- duced above, we discuss simple models of a molecular motor that is fueled by chemical energy and that produces mechanical work. These models are based on considerations of thermodynamic consistency, which are also the basis of more elaborate models mo- tivated by experimental observations, e.g., widely used models for kinesin [42] or the F1-ATPase [43, 44]. Unicyclic models, such as the ones presented here, have recently been used to infer design principles for molecular machines that optimize the output under some natural constraints [45]. We consider an isothermal setup thermostatted by a heat reservoir at temperature T . Embedded in this reservoir is a particle reservoir for substrate and product molecules. To be specific, we take ATP as substrate, which is hydrolyzed to ADP and inorganic phosphate (Pi). The particle reservoir is a dilute solution with concentrations cT, cD, cP of ATP, ADP, and phosphate, respectively. The corresponding chemical potentials are then given by µρ = µ 0 ρ + T ln ( cρ/c0 ) (2.54) with ρ ∈ {T,D, P} and where µ0ρ denotes the chemical potential of species ρ at a concentration of reference c0. In addition, the enzyme is coupled to an athermal work reservoir containing a single degree of freedom X that can be displaced against a mechanical force f . In the formalism introduced in Sec. 2.2, we identify the resources (Xr ) = (NT, ND, NP, X ) with the particle numbers Nρ in the reservoir and the generalized forces (Kr ) = (µT, µD, µP, f ). As the first version of a model for the enzyme proper, we consider a set of three mesostates i ∈ {0,T,D}, comprising the bare enzyme (0), a state with ATP bound to the enzyme (T), and a state with ADP bound to the enzyme (D), as shown in Fig. 2.2a. Each of these states has an intrinsic energy Ei, entropy Si and free energy Fi = Ei −TSi. We allow for six different transitions with rates ki j that are fully identified by the initial state i and the final state j and that form a cycle in the network of states. For simplicity, only the transitions between D and 0 come with a transfer of work to or from the work reservoir. The transition rates for the transition between 0 and T must satisfy the local detailed balance condition (2.21) ln(k0T/kT0) = −(FT − F0 − µT)/T, (2.55) which involves the free energy change of the enzyme and of the particle reservoir providing one ATP molecule, i.e., ∆N0TT = −1. In the step 0 → T, the energy ∆E0T = E0 − ET is transferred to the joint heat and particle reservoirs. Subtraction of the energy that is transferred to the particle reservoir leads according to Eq. (2.28) and with Eq. (2.54) to the heat Q0T = E0 − ET + µ0T that is dissipated into the isothermal 38 2.7 Unicyclic molecular motor as paradigmatic example (a) 0 TD (b) 0 D (c) 0 Figure 2.2:Models for unicyclic motors. (a) Three-state model, (b) two-state model, (c) one- state model. Transitions are indicated by black arrows, along which the particle and work exchanges are illustrated. ATP is colored red, ADP and Pi are colored pink, and work transfered to the work reservoir is indicated as a lifted weight. environment. The transition from T to D involves the release of a phosphate molecule. The local detailed balance condition (2.21) therefore reads ln(kTD/kDT) = −(FD − FT + µP)/T (2.56) and we identify the heat QTD = ET − ED − µ0P. Finally, upon a transition between D and 0, an ADP molecule is released into the reservoir and there is an exchange of work displacing the mechanical degree of freedom by ∆XD0 ≡ d against the force f . The local detailed balance condition (2.21) for this transition is ln(kD0/k0D) = −(F0 − FD + µD + f d)/T (2.57) and work and heat follow as WT0 = f d and QT0 = ED − E0 − µ0P + f d, respectively. Thus, for given chemical potentials and intrinsic free energies, we arrive at a model for which the six transition rates are reduced through the constraints (2.55), (2.56), (2.57) to a set of three model parameters. The steady state of this model system is determined through the stationary master equation (2.20), leading to a stationary cycle current Js ≡ js0T = jsTD = jsD0 (2.58) 39 2 Stochastic thermodynamics on discrete state spaces that flows, according to Kirchhoff’s law in Eq. (2.20), through all links. Due to the coupling of the transitions to the exchange with the reservoirs, this current is connected to the chemical turnover rates ⟨N˙ρ⟩ and the mean velocity v ≡ ⟨X˙⟩ of the mechanical degree of freedom as Js = v/d = −⟨N˙T⟩ = ⟨N˙D⟩ = ⟨N˙P⟩. (2.59) The entropy production rate follows from Eq. (2.38) as σ = Js ( ln k0T kT0 + ln kTD kDT + ln kD0 k0D ) = Js(∆µ − f d)/T ≡ JsA, (2.60) where we define ∆µ ≡ µT − µD − µP and identify the cycle affinity A as the amount of entropy that is produced with every completion of the cycle (see Chap. 6 for a more general treatment of cycle affinities and fluxes). The sign of A therefore determines the direction of the cycle current Js: For ∆µ > f d, ATP is hydrolyzed and work is performed by the system, whereas for the case ∆µ < f d, work is performed on the system resulting in the synthesis of ATP. In the special case A = 0, i.e., ∆µ = f d, there is a balance between the mechanical and the chemical driving, leading to vanishing current Js and entropy production σ. Even for non-equilibrium concentrations of the chemicals in the reservoir, i.e. ∆µ , 0, the total system is then in equilibrium. The mean work current or mechanical output power of the motor is given by P ≡ ⟨W˙ ⟩ = Js f d = f v (2.61) and the rate of heat dissipation reads ⟨Q˙⟩ = Js(µ0T − µ0D − µ0P − f d), (2.62) which is different from the entropy production rate σ comprising also the entropy production in the particle reservoir. The rate of change of the free energy of the solution is described by the chemical work ⟨W˙ chem⟩ = Js∆µ. (2.63) The three-state model is not yet the simplest model for a unicyclic molecular motor. If one of the three states is only short-lived, say, the state T with bound ATP, it is possible to use a description employing only the remaining long-lived mesostates 0 and D, as shown in Fig. 2.2b. The microstates belonging to the state T are then considered as part of the barrier between 0 and D and as such it is irrelevant to which of the two mesostates they are formally attributed. For this model, transitions are no longer uniquely identified by the states of the open system, it is therefore necessary to distinguish them by the 40 2.7 Unicyclic molecular motor as paradigmatic example concomitant changes in the reservoir. The transition 0 → D via the former mesostate T comes with the change a = (∆Ea = E0 − ED, ∆NaT = −1, ∆NaD = 0, ∆NaP = 1, ∆X a = 0) (2.64) and the direct transition from D to 0 is associated with b = (∆Eb = ET − E0 − f d, NbT = 0, ∆NbD = 1, ∆NbP = 0, ∆X b = d), (2.65) using the notation from Eq. (2.16) with the corresponding reverse transitions a¯ and b¯. In total, this model has four transition rates ka0D, k a¯ D0, k b D0, k b¯ 0D, which satisfy the local detailed balance conditions ln ( ka0D/k a¯ D0 ) = −(F0−FT− µT+ µP)/T, ln ( kbD0/k b¯ 0D ) = −(FT−F0+ µD+ f d)/T, (2.66) leaving two free model parameters. A connection between the rates of the original three-state model and the simplified two-state model can be made through appropriate methods of coarse-graining, see, e.g., [44, 46, 47]. The stationary distribution does not depend on the multiplicity of the transition pathways and is therefore determined by the combined rates k0D = ka0D + k b¯ 0D and kD0 = k b D0 + k a¯ D0, leading to ps0 = kD0 k0D + kD0 , psT = k0D k0D + kD0 . (2.67) The cycle current follows as Js = ja,s0D = j b,s D0 = kD0kaD0 − k0Dk a¯0D k0D + kD0 (2.68) and determines the average entropy, heat, and work currents in the same way as in Eqs. (2.60)-(2.63). As the simplest version for a unicyclic molecular motor, for which both mesostates D and T are shortlived, one can reduce the state space of the enzyme to only a single mesostate 0, see Fig. 2.2c. For this model, only one type of transition is possible, which involves the change c = (∆Ec = − f d, N cT = −1, ∆N cD = 1, ∆N cP = 1, ∆X c = d). (2.69) The forward and backward rates for this transition then satisfy the detailed balance relation ln ( kc00/k c¯ 00 ) = (∆µ − f d)/T = A. (2.70) Eqs. (2.60)-(2.63) hold as before with the cycle current Js = kc,s00 − k c¯,s00 . 41 2 Stochastic thermodynamics on discrete state spaces 2.8 Conclusion In this Chapter we have recalled the physical foundations for the stochastic description of systems out of equilibrium. Key concept is the detailed balance condition for transi- tions between mesostates of the total systems, which includes reservoirs. Focusing on transitions in the system proper leads to the local detailed balance relation (2.17). On this level of description, the reservoirs provide constant forces that drive the system into a non-equilibrium steady state. We have illustrated these concepts for minimal models of molecular motors, which are driven by chemical and mechanical forces. Stochastic thermodynamics allows one to identify thermodynamic quantities, such as heat, work, energy and, most notably, the entropy along individual fluctuating trajecto- ries. The key role that the entropy production plays is due to the fact that it quantifies the time-reversibility of trajectories. This property is expressed in terms of a set of fluc- tuation relations, which provide first constraints on the spectrum of current fluctuations. These constraints will turn out to be crucial for the design of new universal bounds later on. 42 3 Fluctuations in stationary Markov processes The objective of this Chapter is to introduce the general mathematical formalism that is used to quantify the fluctuations of the thermodynamic observables introduced in Chap. 2 for non-equilibrium steady states. Two main approaches are discussed: A direct computation of cumulants using algebraic methods [48, 49] and a description using large deviation theory [35, 50–53]. While the former is particularly useful for the exact computation of current fluctuations, the latter provides a powerful tool for proving universal bounds. 3.1 Cumulant generating function for time-additive observables We consider a Markovian network with a discrete and finite state-space {i}. Transitions occur at continuous times with transition rates ki j ≥ 0 (with kii ≡ 0 for all i), where i denotes the initial state and j the final state, giving rise to an ensemble of trajectories i(τ) with τ ≥ 0 on the state space. The time-evolution of the occupation probabilities pi (t) of the states i at some time τ = t is given by the Master equation (2.19), reading ∂tpi (t) = ∑ j Li jp j (t) (3.1) with the matrix Li j ≡ k ji − δi jri (3.2) and the exit rate ri ≡ ∑ j ki j as before in Eq. (2.43). The normalized stationary distribu- tion psi follows as the solution of 0 = ∑ j Li jpsj . (3.3) Along each trajectory i(τ), we define an observable X (t) as X (t) ≡ ∑ i aiTi (t) + ∑ i, j di jmi j (t), (3.4) 43 3 Fluctuations in stationary Markov processes where Ti (t) denotes the time the trajectory has spent in state i up to the time t and mi j (t) counts the number of transitions from i to j that have occurred up to the time t. The vector ai and the “distance-matrix” di j (with dii ≡ 0 for all i) have real entries that specify the observable. This type of observable is called “time-additive”, since for a concatenation of two trajectories (with matching final and initial state), X (t) is the sum of the values of X evaluated along the individual parts. In the mathematical literature, such processes are known as “Markov additive processes” [54–56]. The fluctuating currents associated with heat, work, entropy, etc., as introduced in Chap. 2 form an important subclass of time-additive observables, for which di j = −d ji and ai = 0 hold due to time-antisymmetry. Ultimately, we are interested in the statistics of the fluctuations of the observable X . Since X (t) itself is not Markovian, we start the calculation for the joint process of X (t) and i(t). The evolution of the joint probability p(i, X, t) of i and X in an infinitesimal time interval [t, t + dt] is p(i, X, t + dt) = ∑ j p(i, t + dt | j, t) p( j, X − δX ji, t), (3.5) with the propagator p(i, t + dt | j, t) = δi j + Li jdt for the dynamics of the states {i}. The change of the observable X during the time interval is δXii = aidt for the case that there is no transition and δXi j = di j + O (dt) for a transition from i to j. We thus obtain p(i, X, t + dt) = (1 − ridt)p(i, X − aidt, t) + ∑ i, j k jidt p( j, X − d ji, t) = p(i, X, t) + dt [ − (ri + ai∂X )p(i, X, t) + ∑ j k jie−d j i∂X p( j, X, t) ] , (3.6) where the operator exp(d∂X ) conveys a Taylor series shifting the argument X by d. With this, the partial differential equation governing the evolution of the joint probability becomes ∂tp(i, X, t) = ∑ j ( Li je−d j i∂X − aiδi j∂X ) p( j, X, t) ≡ ∑ j Li j (X ) p( j, X, t), (3.7) in which we identify the operator Li j (X ) acting on both j and X . Due to the time-independent and additive character of the Markov process for i and X , the propagator p(i, X ′′, t′′| j, X ′, t′) depends only on the differences X = X ′′ − X ′ and t = t′′ − t′ and can thus be written as p(i, X, t | j) = ⟨i |etL(X )δ(X ) | j⟩ , (3.8) using the matrix exponential and a bra-ket notation |i⟩ for the distribution corresponding to the system being localized in state i, with the usual rules for multiplication. The 44 3.1 Cumulant generating function for time-additive observables distribution of X for an ensemble of trajectories of length t sampled from the steady state follows by using the stationary distribution psj for the initial state and integrating out the initial and final state. This leads to p(X, t) = ∑ i, j p(i, X, t | j)psj = ⟨− etL(X )δ(X ) ps⟩ , (3.9) where ⟨−| ≡ ∑i ⟨i | is the vector containing 1 in every entry and |ps⟩ ≡ ∑ j psj | j⟩. An explicit calculation of p(X, t) using Eq. (3.9) is typically not possible. Nonetheless, this equation can be used to calculate the cumulants of the distribution of X . For this purpose, we define the moment generating function g(z, t) ≡ ⟨ ezX (t) ⟩ = ∫ dX p(X, t) ezX (3.10) with a typically real parameter z. This function yields the cumulants of the distribution of X by taking the logarithm and repeated derivation at z = 0, cn(t) ≡ ∂nz ln g(z, t) |z=0. (3.11) In particular, themean is given by c1 = ⟨X (t)⟩ and the variance by c2 = ⟨X (t)2⟩−⟨X (t)⟩2. The normalization of p(X, t) is expressed in c0 = ln g(0) = 0. Using Eq. (3.9) and integrating by parts, the cumulant generating function follows as g(z, t) = ⟨− etL(z) ps⟩ (3.12) with the so-called “tilted” matrix Li j (z) ≡ Li jed j i z + aiδi j z = k jied j i z + δi j (−ri + zai). (3.13) In contrast to L(X ), the entries of this matrix are not operators but simple algebraic functions. For small time intervals t, the matrix exponential (3.12) can be evaluated in a Taylor series. For arbitrary t, we instead rely on an expansion of the vectors |ps⟩ and ⟨−| in eigenvectors of L(z). Since L(z) is generally not Hermitian, one has to distinguish the right eigenvector |qν (z)⟩ and the left eigenvector ⟨q˜ν (z) | corresponding to an eigenvalue λν (z) of L(z). A generic spectrum λν (z) as a function of the parameter z is shown in Fig. 3.1. Only at singular branching points some eigenvalues are degenerate1, for all other values of z we can exploit the orthogonality of left and right eigenvectors corresponding to different eigenvalues, yielding etL(z) = ∑ ν etλν (z) qν (z)⟩ ⟨q˜ν (z) , (3.14) 1At the singular values of z where some eigenvalues are degenerate, one can use the Jordan canonical form of L(z) in order to evaluate the matrix exponential. Alternatively, one can evaluate the matrix exponential for non-singular z and perform a limit for singular z. 45 3 Fluctuations in stationary Markov processes ✲ ✲✁ ✲✂ ✵ ✂ ✁ ③ ✲ ✲✁ ✲✂ ✵ ✂ ✁ ❘ ✄ ☎ ✻ ✽ ✭ ✆ ✝ ✞ ✟ ✠ ✡ ☎ ✻ ✽ ✭ ✆ ✝ ✞ Figure 3.1: Generic spectrum of the matrix L(z) as function of the parameter z. Real parts of the eigenvalue are plotted as solid lines and the imaginary part as dashed lines. The system has four states, all transition rates have been drawn at randoma from a uniform distribution between 0 and 1, and the observable is the current along the link 1 → 2, specified through ai = 0, d12 = 1 = −d21. The Perron-Frobenius eigenvalue is shown in blue. aRate matrix L = *...., −2.2814 0.1560 0.6011 0.8324 0.9507 −1.0803 0.7081 0.2123 0.7320 0.0581 −2.2791 0.1818 0.5987 0.8662 0.9699 −1.2266 +////- . where we use the normalization ⟨q˜µ(z) |qν (z)⟩ = δµν and ⟨−|qν (z)⟩ = 1. For large time-intervals t, only the eigenvalue with the largest real part dominates the expansion (3.14). This eigenvalue has several special properties, which are due to the Perron-Frobenius theorem [21]. This theorem applies to matrices with entirely non- negative entries. Since the matrix L(z) contains negative entries only in its diagonal, it can be readily transformed into such a type of matrix by adding a sufficiently large multiple of the identity matrix. This transformed matrix L¯i j (z) = Li j (z) + ηδi j (3.15) with η > maxi[ri − zai] has a spectrum that is shifted by η to λν (z)+ η and has the same eigenvectors as L(z). The matrix L¯(z) is primitive, meaning that for a sufficiently large power L¯(z)n all entries [L¯(z)n]i j = *., n∏ k=2 ∑ ik +/- L¯ii2L¯i2i3 . . . L¯in j (3.16) 46 3.1 Cumulant generating function for time-additive observables become positive. This property becomes obvious by considering the graph spanned by the non-zero transition rates ki j . We can assume that this graph is connected, otherwise one would treat two or more unconnected sets of states as independent systems. We furthermore recall the micro-reversibility introduced in Sec. 2.2, which excludes the existence of “dead-ends”. As a consequence, any two states i, j are connected by a finite path associated with positive transition rates and thereby positive entries of L¯(z), leading to at least one positive term among the non-negative terms of Eq. (3.16). For such a primitive matrix L¯(z), the Perron-Frobenius theorem states that the eigen- value ϱwith the largest absolute value (the so-called spectral radius) is real and positive. We therefore identify λ(z) = ϱ − η as the eigenvalue with the largest real part of L(z), for which we drop the index ν. Moreover, this eigenvalue is non-degenerate and the corresponding right and left eigenvectors are positive for all entries qi (z) = ⟨i |q(z)⟩ and q˜i (z) = ⟨q˜(z) |i⟩. The eigenvectors to all other eigenvalues each have both positive and negative entries. For the special case z = 0, the matrix L(0) = L becomes the generator of the time evolution of themaster equation (3.1)with eigenvalue λ(0) = 0 and eigenvectors ⟨q˜(0) | = ⟨−| and |q(0)⟩ = |ps⟩. The expansion (3.14) for large t then reads exp(tL) ≈ |ps⟩ ⟨−| with exponentially decaying corrections and transforms any positive initial distribution into the stationary one in the long-time limit. For arbitrary z and large t, the matrix exponential (3.14) becomes dominated by etL(z) ∼ etλ(z) |q(z)⟩ ⟨q˜(z) | , (3.17) such that the generating function (3.10) reads ln g(z, t) ≈ tλ(z) + ln ⟨q˜(z) ps⟩ , (3.18) with exponentially decaying corrections. The leading order contribution to g(z, t) is called the “scaled cumulant generating function” λ(z) = lim t→∞ 1 t ln g(z, t) = lim t→∞ 1 t ln ⟨ eXz ⟩ , (3.19) which is equal to the Perron-Frobenius eigenvalue of L(z). Along with the scaling of g(z, t), all cumulants of X scale linearly for large t as well and can be rescaled as Cn ≡ lim t→∞ 1 t cn(t) = ∂nz λ(z) |z=0. (3.20) In particular, we have the mean rate of change Js ≡ C1 and the effective diffusion coefficient D ≡ lim t→∞ 1 2t ( ⟨X (t)2⟩ − ⟨X (t)⟩2 ) = C2/2. (3.21) 47 3 Fluctuations in stationary Markov processes When taking the limit t → ∞, contributions to the observable X (t) that do not scale with time become irrelevant for its fluctuations. In order to make this fact explicit, consider the time-additive observable X¯ (t) = X (t) + ϕi(t) with a state-function ϕi and accordingly the distance matrix d¯i j = di j + ϕ j − ϕi. The corresponding tilted matrix (3.13) can then be written as L¯(z) = V L(z) V−1 (3.22) with the diagonal matrix Vi j ≡ exp(ϕiz)δi j , which has the same spectrum as L(z). Hence, λ(z) and all cumulants Cn must be pairwise equal for X (t) and X¯ (t), thus forming a class of equivalent observables. For instance, such a relation holds for the total entropy production (2.33) and the entropy production in the reservoirs (2.31). However, in sub-leading order in t, the generating functions within such a class of observables differ due to the differences in the eigenvectors. The formalism we have discussed so far is also applicable to Markov processes with several transition pathways with transition rates labeled by kai j , as introduced in Sec. 2.2. The additive observable X (t) is then specified with a generalized distance matrix dai j . Along the same lines as above, the scaled cumulant generating function for the fluctuations of X (t) is derived as the Perron-Frobenius eigenvalue of the tilted matrix Li j (z) ≡ ∑ a kajie zdaji + δi j (−ri + aiz) (3.23) with the exit rates ri ≡ ∑a, j kai j . As a brief illustration, we revisit the simple motor models from Sec. 2.7 and calculate the scaled cumulant generating function for the displacement of the mechanical degree of freedom X (t). The one-state model is described by two rates kc00 ≡ k+ and k c¯00 ≡ k−, which satisfy the local detailed balance condition (2.66) reading ln ( k+/k− ) = A. The observable of interest is specified by dc00 = −d c¯00 = 1. This setup is equivalent to the asymmetric random walk X (t) along a uniform one-dimensional track with step size d = 1. The tilted matrix (3.23) has in this case only one entry which is thus equal to its eigenvalue λ(z) = L00(z) = k+ez + k−e−z − k+ − k− = 2 √ k+k− [ cosh ( z + A 2 ) − cosh (A 2 )] . (3.24) The rate of change and the diffusion coefficient follow as Js = λ′′(0) = k+ − k− and 2D = λ′(0) = k+ − k−, respectively. The two-state model with rates ka0D ≡ k+, k a¯D0 ≡ k−, kbD0 ≡ w+ and k a¯D0 ≡ k− corresponds to an asymmetric random walk on a one-dimensional lattice with two 48 3.2 Calculation of low-order cumulants alternating states with periodicity d = 1 and accordingly alternating rates. Here, the tilted matrix reads L(z) = ( −k+ − w− k− + w+ez k+ + w−e−z −k− − w+ ) (3.25) and has the Perron-Frobenius eigenvalue λ(z) = −Σ/2 + √ (Σ/2)2 + k+w+(ez − 1) + k−w−(e−dz − 1) (3.26) with Σ ≡ k+ + k− + w+ + w−. The rate of change becomes, in accordance with Eq. (2.68), Js = λ′(0) = (k+w+ − k−w−)/Σ and the diffusion coefficient is 2D = λ′′(0) = [k+w+ + k−w− − 2(Js)2]/Σ . 3.2 Calculation of low-order cumulants An analytic calculation of the scaled cumulant generating function as eigenvalue of the tilted matrix L(z) is possible for less than three states or particularly sparse transition matrices. In most other cases, one has to resort to a numerical computation of the eigenvalues in order to obtain the full generating function. Often one is interested only in low-order cumulants Cn of the fluctuations of an additive observable X (t), in particular the rate of change Js = C1 and the diffusion coefficient D = C2/2. However, a numerical calculation of the derivatives of the generating function is typically prone to discretization errors and round-off errors. Here, we discuss two methods that allow for an exact computation of the low-order cumulants, avoiding the direct evaluation of eigenvalues. Koza method. The method devised by Koza [48] starts by considering the character- istic polynomial χ(y, z) ≡ det(L(z) − y1) (3.27) of the tiltedmatrix, with the identitymatrix 1. This functionmust be zerowhen evaluated for the eigenvalue λ(z), i.e., χ(λ(z), z) = 0. (3.28) The characteristic polynomial is expanded in the power series χ(y, z) = ∞∑ n,m=0 bnmynzm (3.29) 49 3 Fluctuations in stationary Markov processes with coefficients bnm = n!m!∂ny ∂mz χ(y, z) |y,z=0. The power series of λ(z) is by definition λ(z) = ∞∑ k=1 Ck zk/k!. (3.30) Plugging these two series into Eq. (3.28) leads to 0 = ∞∑ n,m=0 bnmzm(C1z + C2z2/2 + . . .)n = b00 + (b01 + b10C1)z + (b02 + b10C2/2 + b20C21 + b11C1)z 2 + O(z3). (3.31) Collecting equal powers in z allows for a recursive determination of the cumulants Ck as a function of the coefficients bnm. In particular, one obtains for the two lowest cumulants Js = C1 = −b01/b10 (3.32) and D = C2/2 = −(b02 + b20C21 + b11C1)/b10. (3.33) The advantage of this method is that it yields analytic expressions for Js and D with- out needing the stationary distribution or other eigenvectors. However, with increasing number of states, these expressions soon become very lengthy. For a numerical evalu- ation of the cumulants, the Koza method is not suitable, since the determination of the coefficients bnm requires again error-prone numerical derivatives. Perturbative method. A second method for the calculation of low-order cumulants uses perturbation theory for a recursive determination of the coefficients of the power series for λ(z). This method has been used in applied contexts [57] and has been discussed in generality in Ref. [49]. We use the expansion L(z) = L + zL′ + z2L′′/2 + O(z3), (3.34) where L = L(0) is the Markov transition matrix (3.1) and the derivatives of L(z) at z = 0 are L′i j = k jid ji + aiδi j and L′′i j = k jid2ji, and so forth for higher derivatives. The eigenvalue is expanded as in Eq. (3.30), likewise we write for the right eigenvector |q(z)⟩ = ps⟩ + z q′⟩ + z2 q′′⟩ /2 + O(z3), (3.35) 50 3.3 Large deviation functions wherewe have used that |q(0)⟩ = |ps⟩. The normalization condition ⟨−|q(z)⟩ is respected by requiring ⟨−|q(n)⟩ = 0 for all n ≥ 1. Plugging these power series into the eigenvalue equation L(z) |q(z)⟩ = λ(z) |q(z)⟩ leads to a recursive scheme for the calculation of the coefficients Cn and |q(n)⟩. The zeroth order gives L |ps⟩ = 0, which is transformed into a solvable system of equations for |ps⟩ by adding the normalization condition ⟨−|ps⟩. In linear order, one gets L′ ps⟩ + L q′⟩ = C1 ps⟩ . (3.36) Multiplication from the left with ⟨−| makes the term with ⟨−| L = 0 vanish, leading to Js = C1 = ⟨− L′ ps⟩ =∑ i j d jik jipsj + ∑ i psi ai, (3.37) which we have already intuitively used to calculate rates of change in Sec. 2.5. Next, the vector |q′⟩ is determined through the linear set of equations L q′⟩ = (C1 − L′) ps⟩ . (3.38) The matrix L with eigenvalue 0 is not invertible. Nonetheless, Eq. (3.38) can be transformed into a solvable system of equations by replacing one line with the condition ⟨−|q′⟩ = 0. The second order of the eigenvalue equation gives 1 2 L′′ ps⟩ + L′ q′⟩ + 12L q′′⟩ = 12C2 ps⟩ + C1 q′⟩ . (3.39) Multiplication with ⟨−| yields the diffusion coefficient 2D = C2 = ⟨− L′′ ps⟩+ 2 ⟨− L′ q′⟩ =∑ i j [d2jik jip s j + 2d jik jiq ′ j]+ 2 ∑ i q′iai . (3.40) The term with |ps⟩ may be interpreted as the instantaneous quadratic change of X in the stationary state, whereas the term with |q′⟩ conveys the contribution from temporal correlations. Higher order cumulants and vectors |q(n)⟩ can subsequently be derived from the corresponding powers of z in the eigenvalue equation. The advantage of the perturbative method is that it can directly be implemented numerically, requiring only the solutions of linear systems of equations. The resulting cumulants are exact up to machine precision. 3.3 Large deviation functions The mathematical theory of large deviations provides a powerful toolbox to describe and compute the probability of atypical events. In this section, we give a brief introduction 51 3 Fluctuations in stationary Markov processes into this theory with a special focus on additive observables. More in-depth discussions of the applications of large deviation theory in physics can be found in the review articles [26, 50]. In general, large deviation theory is concerned with random variables A, whose probability distribution p(A, N ) depends on a scaling parameter N . Such a distribution is said to satisfy a large deviation principle if the limit I (A) ≡ − lim t→∞ 1 t ln p(A, N ) (3.41) exists2. The function I (A) is then called “large deviation function” or “rate function”. We use the short-hand notation p(A, N ) ∼ e−t I (A) . (3.42) The variable A may be a scalar, a vector, or a more complex object such as a field or a trajectory. Scaling parameters that are typically encountered in physics are the time (then called t), numbers of particles or the reciprocal of a noise intensity. Since the distribution must be normalized, the large deviation function must be non-negative for all values of A. Values of A for which I (A) = 0 are identified as “typical values”. Their probability scales sub-exponentially with N . For atypical values, I (A) > 0 gives the rate of the exponential decay of the probability. Values of A whose probability decays even faster or which are entirely impossible are formally assigned I (A) = ∞. A classic example for a random variable satisfying a large deviation principle is the ratio A = M/N of outcomes “head” in a coin-toss experiment that is repeated N times. If the probability for “head” in an individual experiment is p, the probability for having M such events is p(M, N ) = N! M!(N − M)!p M (1 − p)N−M . (3.43) Using Stirling’s approximation of the factorials, the logarithm of this distribution be- comes ln p(M, N ) = (M−N ) ln ( 1 − M N ) −M ln M N +M ln p+ (N−M) ln(1 − p)+O (ln N ) , (3.44) such that the large deviation function for p(A, N ) can be read off directly as I (A) = (1 − A) ln(1 − A) + A ln(A) − A ln p − (1 − A) ln(1 − p), (3.45) 2In order to define the large deviation principle also for distributions with singular points (e.g., sums of δ-functions), the distributions in Eq. (3.41) are understood to be integrated in an interval [A, A + δA] with small but non-zero δA. The limit δA → 0 is then performed after taking the limit N → ∞. In addition, this convention ensures that the logarithm has always dimensionless arguments. 52 3.3 Large deviation functions 0 0.2 0.4 0.6 0.8 1 A 0 0.2 0.4 0.6 0.8 1 I(A ) p=0.5 p=0.3 Figure 3.2: Large deviation function (3.45) (solid) for the coin-toss experiment with p = 0.5 (blue) and p = 0.3 (red). The dashed curves are the quadratic approximations for the typical values of A. with 0 ln 0 ≡ 0, and is defined for all possible values 0 ≤ A ≤ 1. As shown in Fig. 3.2, this function has the minimum I (A) = 0 at A = p = ⟨M⟩ /N , i.e., the average value is also the typical one. Around this typical value, the large deviation function can be approximated in quadratic order as I (A) ≈ I′′(p)(A − p)2/2 with I′′(p) = 1/[p(1 − p)]. This leads to the Gaussian approximation I (A) ≈ 1√ 2πp(1 − p)/N exp [ −N (A − p) 2 2p(1 − p) ] (3.46) with variance N/I′′(p) = Np(1−p). This distribution is also the result of the central limit theorem, treatingM as a sumof independent randomvariables that are 1 (with probability p) or 0. Thus, the large deviation principle may be regarded as a generalization of the central limit theorem for extreme events. We now focus on time-additive observables X (t) in aMarkov process, as introduced in Sec. 3.1. Scaling such a time-extensive observable by time t leads to the time-intensive observable J (t) ≡ X (t)/t, which expresses the average change of X along an individual trajectory of length t. The distribution of this observable satisfies a large deviation principle of the form p(X, t) ∼ p(J, t) ∼ e−t I (J), (3.47) 53 3 Fluctuations in stationary Markov processes where we note that the distributions of X and J, p(X, t) = p(J, t)/t, differ only on a sub-exponential scale. This large deviation function proves consistent with the scaling of the cumulant generating function (3.19), as can be seen by writing the exponential average of X as etλ(z) ∼ ⟨ eXz ⟩ = ∫ dJ p(J, t)et Jz ∼ ∫ dJ e−t I (J)+t Jz . (3.48) The integral on the right hand side can be evaluated using a saddle-point approximation of the integrand, which reduces in the leading exponential order to taking the maximum of the exponent, leading to λ(z) = max J [Jz − I (J)]. (3.49) For convex functions I (J), this relation is equivalent to the Legendre transform of I (J), which reads λ(z) = J (z)z − I (J (z)) (3.50) with J (z) determined by z = I′(J). If I (J) was non-convex, λ(z) would be the Legendre transform of the convex envelope of I (J) and would have non-analytic kinks at values of z that correspond to the linear slopes of the convex envelope. However, for Markovian systems with a finite state space, λ(z) follows as a non-degenerate root of the analytic characteristic polynomial (3.27) and must therefore be smooth for all z. Thus, the large deviation function must be convex. Nonetheless, singular points in generating functions are commonly observed in Markovian systems with infinite state-spaces and are studied in the context of dynamical phase transitions [26, 28, 58]. For convex I (J), the Legendre transform (3.49) presents an invertible mapping be- tween I (J) and the likewise convex λ(z). The former can therefore be gained from the scaled cumulant generating function by inverting the Legendre transform to I (J) = max z [Jz − λ(z)] = Jz(J) − λ(z(J)) (3.51) with z(J) determined by J = λ′(z), which is the inverse of the mapping J (z) above. This relation between scaled cumulant generating function and large deviation function, known as Gärtner-Ellis Theorem, allows one to calculate the large deviation function of an additive observable from the Perron-Frobenius eigenvalue of the tilted matrix (3.13) without need for an explicit knowledge of the probability distribution of the observable. The typical behavior of an additive observable is described by the mean rate of change Js = λ′(0), therefore z(Js) = 0 and I (Js) = 0, identifying Js as the unique minimum of 54 3.3 Large deviation functions -10 -5 0 0 5 (z ) 0 2 4 6 0 5 10 15 I(J ) Figure 3.3: Plot of the scaled cumulant generating function (3.24) (left) and the large deviation function (3.54) (right) for the asymmetric random walk with k+ = 1 and the affinity per stepA = 8. Selected pairs of corresponding points of z and J (z) are marked by symbols. the large deviation function. The curvature of the large deviation function at this point is I′′(Js) = 1 λ′′(0) = 1 2D , (3.52) such that the Gaussian approximation of the distribution p(X, t) becomes p(X, t) = 1√ 4πDt exp [ − (X − J st)2 4Dt ] , (3.53) which matches a normal biased diffusion in one dimension with drift velocity Js and diffusion coefficient D. For the example of the asymmetric random walk, the Legendre transform of the generating function (3.24) can be performed analytically [35]. The mapping between J and z reads here J = λ′(z) = J0 sinh(z − A/2) with J0 ≡ 2 √ k+k− = Js/ sinh(A/2). Inserting this relation in Eq. (3.51) leads to the large deviation function I (J) = J [ arsinh(J/J0) − arsinh(Js/J0)] −√(J0)2 + J2 +√(J0)2 + (Js)2. (3.54) Plots of this scaled cumulant generating function and large deviation function are shown in Fig. 3.3. As a typical feature for time-asymmetric observables, the large deviation function has a high curvature or “kink” around J = 0 [51, 59]. This value of J corresponds to the minimum of λ(z) at z = −A/2 with a particularly low curvature. 55 3 Fluctuations in stationary Markov processes If the distribution of the observable of interest obeys an asymptotic fluctuation relation of the type (2.51), the corresponding large deviation function satisfies I (−J) = I (J) + J . (3.55) In the Legendre transformed picture (3.49), this fluctuation relation expresses as the symmetry λ(z) = max J [zJ − I (J)] = max J [zJ + J − I (−J)] = λ(−1 − z) (3.56) of the scaled cumulant generating function. As an important example, consider the entropy production (2.31) by setting di j = ln ( ki j/k ji ) and ai = 0. Contributions from the intrinsic entropy, which is a state function that cannot scale with time, do not have to be taken into account (see also the discussion leading to Eq. (3.22)). Following the proof of the Gallavotti-Cohen symmetry provided in Ref. [35], the tilted matrix (3.13) for this observable can be written as Li j (z) = k ji ( k ji ki j ) z − riδi j = L ji (−1 − z). (3.57) Thus, not only the Perron-Frobenius eigenvalue, but also the whole spectrum of L(z) satisfies the symmetry (3.56) around the center at z = −1/2, for which L(z) becomes symmetric. The left and right eigenvectors exchange their roles upon inflection of z, qi (z) ∝ q˜i (−1 − z), (3.58) which holds due to the transposition of L(z) in Eq. (3.57). Focusing again on general additive observables X (t), we conclude the introduction to large deviations with a discussion of sub-exponential contributions to probability distributions [36,60]. Writing the distribution of J (t) as3 p(J, t) = p0(J, t) exp[−t I (J)], where p0(J, t) does not scale exponentially with t, the saddle-point approximation of the generating function (3.48) can be refined to eg(z,t) = ⟨ eXz ⟩ = ∫ dJ p0(J, t)et[Jz−I (J)] ≈ √ 2π t I′′(J∗) p0(J∗, t)et[J ∗z−I (J∗)] (3.59) with the saddle-point J∗ ≡ J (z) as above. Comparison to the scaling (3.18) of the cumulant generating function reveals p0(J, t) ≈ √ t I′′(J) 2π ⟨−|q(z)⟩ ⟨q˜(z) ps⟩ (3.60) 3In this context, it is again essential to evaluate distributions at first on a finite interval [J, J + δJ] and to take the limit δJ → 0 after any long-time limit t → ∞. Otherwise, the absolute value of X (t) = t J (t) may contain information about the initial or finite state of trajectories, leading to a fine-structure on the sub-exponential scale [36]. 56 3.3 Large deviation functions for any matching pair of J and z = z(J), with an error that decays exponentially. Using similar methods, it is possible to infer the probability of state i at an intermediate time τ conditioned on the observation of a large deviation J in a long time interval t, which we denote as p(i, τ |J, t). The corresponding conditioned exponential average reads analogously to (3.59) ⟨ eXzδi,i(τ) ⟩ = ∫ dJ p(i, τ |J, t)p(J, t)eJz ≈ √ 2π t I′′(J∗) p(i, τ |J∗, t)p0(J∗, t)et[J∗z−I (J∗)]. (3.61) On the other hand, using the methods for the propagator from Sec. 3.1, we calculate⟨ eXzδi,i(τ) ⟩ = ∫ dX1 ∫ dX2 ∑ i0,i1 e(X1+X2)zp(i1, X2, t − τ |i)p(i, X1, τ |i0)psi0 = ⟨−|e(t−τ)L(z) |i⟩ ⟨i eτL(z) ps⟩ ≈ eλ(z) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ⟨−|q(z)⟩ ⟨q˜(z) |i⟩ ⟨i |ps⟩ , τ = 0 ⟨−|q(z)⟩ ⟨q˜(z) |i⟩ ⟨i |q(z)⟩ ⟨q˜(z) |ps⟩ , 0 < τ/t < 1 ⟨−|i⟩ ⟨i |q(z)⟩ ⟨q˜(z) |ps⟩ , τ = t , (3.62) where we have used Eq. (3.17) for evaluating the matrix exponentials in the limit of long time intervals. Comparing this to Eqs. (3.61) and (3.60) yields the following picture for the ensemble of trajectories with J (t) = J for large t [26]: The initial state of these trajectories is distributed as p(i, 0|J, t) = ⟨−|q(z)⟩ ⟨q˜(z) |i⟩ ⟨i |p s⟩ ⟨−|q(z)⟩ ⟨q˜(z) |ps⟩ = q˜i (z)psi∑ i q˜i (z)psi , (3.63) with z = z(J) throughout, and the distribution of final states is p(i, t |J, t) = ⟨−|i⟩ ⟨i |q(z)⟩ ⟨q˜(z) |p s⟩ ⟨−|q(z)⟩ ⟨q˜(z) |ps⟩ = qi (z)∑ i qi (z) . (3.64) For most times τ in between that are sufficiently far form the boundaries, i.e., τ = O (t) and t − τ = O (t), the distribution of states is p(i, τ |J, t) = ⟨−|q(z)⟩ ⟨q˜(z) |i⟩ ⟨i |q(z)⟩ ⟨q˜(z) |p s⟩ ⟨−|q(z)⟩ ⟨q˜(z) |ps⟩ = q˜i (z)qi (z). (3.65) As we have seen in the derivation, the same distributions of i(τ) show up in the ensemble of trajectories that is re-weighted by the exponential exp(Xz) but which is otherwise unconstrained. This correspondence of the ensembles of trajectories in the long time limit is reminiscent of the equivalence of the micro-canonical and canonical ensemble in the thermodynamic limit [61]. 57 3 Fluctuations in stationary Markov processes 3.4 Level 2.5 large deviation theory So far, we have focused on the large deviation theory for an individual observable, which is called the “level 1” of large deviation theory. When dealing with the large deviations of the joint probabilities of several observables, a useful tool is the following contraction principle. Consider a joint probability p(J1, J2, t) of two observables J1 and J2 satisfying a large deviation principle of the type p(J1, J2, t) ∼ exp[−t I (J1, J2)] with a two-dimensional large deviation function I (J1, J2) (the generalization to more than two variables will be obvious). The distribution of a third, dependent observable defined as a function K = K (J1, J2) is then obtained for finite times through the integration p(K, t) = ∫ dJ1 ∫ dJ2 p(J1, J2, t) δ[K − K (J1, J2)], (3.66) which is typically hard to perform analytically. For large times t, only the values of J1 and J2 that lead to the smallest rate of decay in the integrand contribute significantly to the overall integral. Therefore, the large deviation function for K can be identified through a constrained minimization as I (K ) ≡ lim t→∞ 1 t p(K, t) = min J1,J2 |K (J1,J2)=K I (J1, J2), (3.67) which should be easier to perform than the full integration (3.66). This contraction principle allows for a complementary approach to the calculation of the large deviation function I (J) of a single variable. The basic idea is to write J as a function of several basic observables whose large deviation function has a closed analyt- ical form. The function I (J) then follows as the result of a constrained minimization of this “master”-large deviation function. Several levels of such large deviation functions are known for Markov processes. Higher levels correspond to higher-dimensional sets of basic observables, allowing for a larger class of observables that is covered through contraction. However, performing the contraction from higher levels generally becomes more tedious. A large deviation functional taking whole trajectories as arguments is commonly known as “level 3”. Less powerful is the “level 2” large deviation function, which takes as arguments the fractions of time spent in each of the states. Most relevant for the analysis of non-equilibrium steady states is the large deviation function of the intermediate “level 2.5”, which describes the joint fluctuations of the time spent in each state Ti (t) and the number of transitions mi j (t) between pairs of states [62, 63]. Scaling these time-additive observables by time leads to the empirical distribution pi ≡ Ti (t)/t (3.68) and the empirical flow µi j ≡ mi j (t)/t, (3.69) 58 3.4 Level 2.5 large deviation theory which fluctuate for an ensemble of trajectories of fixed length t around the stationary values psi and µ s i j ≡ psi ki j , respectively. The level 2.5 large deviation function can now be written as I ({pi}, {µi j }). A scaled general time-additive observable J = X (t)/t as defined in Eq. (3.4) is a linear combination J = J ({pi}, {µi j }) ≡ ∑i j di j µi j +∑i aipi of these basic observables, such that its (level 1) large deviation function follows through the contraction principle I (J) = min {pi },{µi j }|J ({pi },{µi j })=J I ({pi}, {µi j }). (3.70) For a derivation of the closed form of the level 2.5 large deviation function, we follow loosely the presentations in Refs. [29, 64, 65]. At first, we note that the empirical flows are not independent. Since the number of jumps into a state and out of it must be balanced,∑ j [m ji (t) − mi j (t)] ∈ {−1, 0, 1} (3.71) must hold for every trajectory and every state i. Accordingly, the empirical flows are constrained by the Kirchhoff-type law lim t→∞ ∑ j [µ ji − µi j] = 0. (3.72) For values of {µi j } that do not satisfy this constraint, we assign I ({pi}, {µi j }) = ∞. The importance of the empirical distribution and flow is founded in their determining the path weight (2.42) as p[i(t)] = psi(0) exp ⎡⎢⎢⎢⎢⎢⎣− ∑ i riTi + ∑ i j mi j ln ki j ⎤⎥⎥⎥⎥⎥⎦ = psi(0) exp ⎡⎢⎢⎢⎢⎢⎣−t *., ∑ i ripi − ∑ i j µi j ln ki j+/- ⎤⎥⎥⎥⎥⎥⎦ . (3.73) A direct marginalization of this path weight in order to obtain the joint distribution p({pi}, {µi j }, t) is rather tedious, since it requires to specify an integration measure in path space. Nonetheless, it is possible to compare these distributions and thus the level 2.5 large deviation function to those of a second, auxiliary process with rates {kˆi j }. Labeling quantities referring to this process by pˆ, Iˆ, etc., one obtains I ({pi}, {µi j })− Iˆ ({pi}, {µi j }) = − lim t→∞ 1 t ln p({pi}, {µi j }, t) pˆ({pi}, {µi j }, t) = ∑ i (ri−rˆi)pi− ∑ i j µi j ln ki j kˆi j . (3.74) 59 3 Fluctuations in stationary Markov processes If the set of rates {kˆi j } can be chosen such that {pi} = {pˆsi } and {µi j } = { µˆsi j } are the typical distribution and flow for the auxiliary process, respectively, one will have Iˆ ({pi}, {µi j }) = 0, such that Eq. (3.74) yields the large deviation function for the process of interest. Luckily, even though it is typically computationally hard to determine the stationary state and flow for a given set of rates, the reverse is performed straightforwardly by setting kˆi j = µˆsi j/pˆ s i . Plugging this into Eq. (3.74) finally yields the level 2.5 large deviation function I ({pi}, {µi j }) = ∑ i *.,ri − ∑ j µi j pi +/- pi− ∑ i j µi j ln piki j µi j = ∑ i ripi+ ∑ i j µi j ( ln µi j piki j − 1 ) (3.75) An alternative representation of this result is obtained by writing the empirical flow as µi j = ( ji j + ti j )/2 with the antisymmetric empirical current ji j ≡ µi j − µ ji = − j ji (3.76) and the empirical traffic ti j ≡ µi j + µ ji = t ji . (3.77) Eq. (3.75) reads with this substitution I ({pi}, { ji j }, {ti j }) = ∑ i ripi + ∑ i< j [ −ti j + ti j + ji j2 ln ti j + ji j 2piki j + ti j − ji j 2 ln ti j − ji j 2p j k ji ] . (3.78) Kirchhoff’s law (3.72) constrains the currents as∑ j ji j = 0 (3.79) for large times, whereas the traffic is a set of positive but otherwise unconstrained variables. In stochastic thermodynamics, we are mostly interested in the fluctuations of time- antisymmetric observables, which are formed as a linear combination of the currents { ji j }. It would therefore be desirable to eliminate the traffic {ti j } and the empirical distribution {pi} from Eq. (3.78) through contraction. The mutual independence of the traffic variables allows one to perform the minimization for {ti j } term by term, which leads to the traffic t∗i j ≡ √ 4piki jp j k ji + j2i j (3.80) 60 3.4 Level 2.5 large deviation theory that minimizes Eq. (3.78) to I ({pi}, { ji j }) = ∑ i< j ⎡⎢⎢⎢⎢⎢⎢⎣piki j + p j k ji − √ 4piki jp j k ji + j2i j + ji j ln ji j + √ 4piki jp j k ji + j2i j 2piki j ⎤⎥⎥⎥⎥⎥⎥⎦ . (3.81) Using the identity arsinh x = ln[x + √ x2 + 1] and defining api j ≡ 2 √ piki jp j k ji and jpi j ≡ piki j − p j k ji leads to a representation of Eq. (3.81) as I ({pi}, { ji j }) = ∑ i< j ⎡⎢⎢⎢⎢⎣ √ (api j ) 2 + ( jpi j ) 2 − √ (api j ) 2 + ( ji j )2 + ji j *,arsinh ji j api j − arsinh jpi j api j +- ⎤⎥⎥⎥⎥⎦ , (3.82) which makes its minimum at pi = psi and ji j = j s i j = p s i ki j − psj k ji evident. A similar expression applies to systems with several transition pathways between pairs of states, as introduced in Sec. 2.2. In that case, transitions are additionally labeled by the upper index a and the sum in Eq. (3.82) is understood to run over every pair of forward and backward transition once. For the asymmetric random walk with a single state 0, the empirical density must be p0 = 1 due to normalization and there is a single independent current J ≡ ja00. The sum in Eq. (3.82) then has only a single term, which is equal to the result from Eq. (3.54). The level 2 large deviation function I ({pi}) for the empirical distribution pi only follows by contraction of Eq. (3.82) over the empirical current { ji j }. In the general case of a driven system, this contraction cannot be performed analytically, since the constraint ∑ j ji j = 0 prohibits an individual minimization of the terms in Eq. (3.82). Nonetheless, for equilibrium systems with rates satisfying ki j/k ji = peqj /p eq i , the current ji j = 0 presents a unique minimum of I ({pi}, { ji j }) for fixed pi. Indeed, the linear term in the expansion I ({pi}, {ε ji j }) = I ({pi}, {0}) − ε2 ∑ i ln pi peqi ∑ j ji j + O ( ε2 ) (3.83) vanishes for all viable ji j and the convexity of the rate function ensures that this stationary point is a unique minimum. The corresponding level 2 large deviation function then follows from Eq. (3.81) as I ({pi}) = I ({pi}, {0}) = ∑ i< j (√ piki j − √ p j k ji )2 . (3.84) A further contraction of Eq. (3.82) to I ({ ji j }) is similarly not possible in a closed analytic form, because of the constraint ∑ i pi = 1. The main benefit of Eq. (3.82) for 61 3 Fluctuations in stationary Markov processes non-equilibrium systems is therefore its providing upper bounds on the large deviation function for a time antisymmetric observable. Such a bound follows by choosing an arbitrary, typically non-optimal, trial distribution {pti } in the contraction principle I ({ ji j }) = min{pi } I ({pi}, { ji j }) ≤ I ({p t i }, { ji j }). (3.85) Bounds of this type have been introduced in Ref. [8] and will turn out to be crucial for the proof of several thermodynamic bounds on current fluctuations. 3.5 Continuous state spaces So far, we have focused on Markovian processes on discrete state spaces. As motivated in Chapter 2, this description is appropriate for systems with well separated mesostates, such as enzymes or quantum dots. In this section, we briefly discuss the generalization of the present formalism to Markovian processes in continuous state spaces. Such a model applies to systems with a few continuous mesoscopic degrees of freedom, provided that their dynamics exhibit a clear time-scale separation from the faster dynamics of the microscopic degrees of freedom. The latter then reach local equilibria for any constellation of the mesoscopic degrees of freedom. A paradigmatic example for such a system is the Brownian motion of a colloidal particle suspended in a fluid that acts as a heat bath with temperature T . A coarse- grained, stochastic description of this system involves only the position of the colloid, whereas its velocity and the degrees of freedom of the fluid particles are assumed to be locally equilibrated. Here, we focus on the Brownian dynamics of a single, one-dimensional degree of freedom x in a homogenous environment. More general cases are discussed in the literature, e.g., in Refs. [6, 52, 53, 66]. We assume that there is a position-dependent potential V (x) and that a constant mechanical force f is driving the system. In order to make the connection to the discrete state spaces, we finely discretize the state space as a lattice with equidistant sites xi = ih with spacing h. We allow for transitions between neighboring sites xi and xi+1 with transition rates k+i for the forward transition and k − i+1 for the backward transition. The local detailed balance condition (2.21) k+i k−i+1 = e[ f h+V (xi )−V (xi+1)]/T (3.86) for these transitions can be satisfied by setting k+i = k0e θi [ f h+V (xi )−V (xi+1)]/T and k−i+1 = k0e −(1−θi )[ f h+V (xi )−V (xi+1)]/T (3.87) with a rate of reference k0 and arbitrary real numbers θi. Subsequently expanding in h leads to the current ji (t) ≡ pi (t)k+i −pi+1(t)k−i+1 = k0h T [pi (t)(g−∂xV (x))− (pi+1(t)−pi (t))/h]+O ( h2 ) , 62 3.5 Continuous state spaces (3.88) which is independent of the parameters θi (provided that they do not scale with h). We smoothen the discrete probability pi (t) and current ji (t) to the density fields p(x, t) and j (x, t), such that p(xi, t) = pi (t)/h and j (xi, t) = ji (t), respectively. The continuum limit exists if k0 scales inverse quadratically in h, leading to j (x, t) = µ( f − ∂xV (x))p(x, t) − D∂xp(x, t), (3.89) where we identify the mobility µ ≡ limh→0 k0h2/T and the diffusion constant D ≡ limh→0 k0h2. The Einstein relation D = µT emerges naturally from the local detailed balance condition (3.86). The conservation of probability, formerly expressed in the master equation (2.19), becomes in the continuum limit the Fokker-Planck equation ∂tp(x, t) = −∂x j (x, t) = −∂xµF (x)p(x) + D∂2x p(x, t), (3.90) with the total force F (x) ≡ f − ∂xV (x). For periodic boundary conditions, these dynamics lead to a steady state with a stationary probability distribution ps(x) and a stationary current js(x) ≡ js that is independent of x for the one-dimensional system. An equivalent description of the stochastic dynamics underlying the Fokker-Planck equation is given by the Langevin equation. Following standard methods [12,67], it can be read-off from Eq. (3.90) as x˙(t) = µF (x(t)) + ζ (t) (3.91) with a Gaussian white noise ζ (t) with zero mean and correlations ⟨ζ (t)ζ (t′)⟩ = 2Dδ(t − t′). A time-additive observable X (t) is defined in the discretized model by specifying the tensors ai and di j in Eq. (3.4). Preparing for the continuum limit, we write ai ≡ a(xi)h and di,i+1 ≡ d(xi)h with continuous functions a(x) and d(x). Since the total number of transitions, or “traffic”, diverges and is thus an ill-defined concept for trajectories on a continuous state-space, we restrict the distance matrix to be antisymmetric, i.e., di,i+1 = −di+1,i. After taking the limit h → 0, the dynamics of X (t) is given by the Langevin equation X˙ (t) = a(x(t)) + d(x(t)) x˙(t) = a(x(t)) + d(x(t))µF (x(t)) + d(x(t))ζ (t). (3.92) Due to the possible dependence of d(x) on x, this Langevin equation contains a multi- plicative noise term, which is to be interpreted in the Stratonovich sense. This way, the time-antisymmetric character of the d-dependent term is preserved. The mean rate of change of X (t) in the stationary state is given by ⟨X˙⟩ = ∫ dx [ps(x)a(x) + js(x)d(x)]. (3.93) 63 3 Fluctuations in stationary Markov processes For instance, the entropy production in the medium (assuming a constant intrinsic entropy) is described by ai = 0 and the generalized distance matrix di,i+1 = ln ( k+i /k − i+1 ) , see Eq. (2.31), which becomes d(x) = [ f − ∂xV (x)]/T in the continuum limit. For the total entropy production (2.32), the term −∂x ln ps(x) is added, leading with Eq. (3.89) to d(x) = js(x)/Dps(x) [25]. For both definitions, the mean entropy production rate in the steady state becomes σ = ∫ dx js(x)[ f − ∂xV (x)]/T = f js/T, (3.94) reflecting that the mechanical power delivered by the external force, on average f js, must ultimately be dissipated into the medium. Eqs. (3.91) and (3.92) present a system of two coupled stochastic differential equations with fully correlated noise terms. The corresponding Fokker-Planck equation for the joint probability of x and X is given by [12] ∂tp(x, X, t) = [−∂xµF (x) + a(x)∂X − d(x)µF (x)∂X + D(∂x + d(x)∂X )2]p(x, X, t) ≡ L(X )p(x, X, t), (3.95) which is analogous to Eq. (3.7) for discrete state spaces. In the picture of the generating function (3.12), the operator L(X ) acting on both x and X gets transformed to L(z) ≡ −∂xµF (x) − a(x)z + d(x)µF (x)z + D(∂x − d(x)z)2 (3.96) acting on x only. The scaled cumulant generating function λ(z), as defined in Eq. (3.19), the follows as the largest eigenvalue of L(z), or, equivalently and typically more conve- niently, of its adjoint L†(z) ≡ µF (x)∂x − a(x)z + d(x)µF (x)z + D(∂x + d(x)z)2. (3.97) The framework of level 2.5 large deviations can be used for the continuous state space as well. The fluctuating empirical density is here defined for an individual trajectory x(τ) as p(x) ≡ ∫ t 0 dτ δ(x − x(τ)), (3.98) and the empirical current density is j (x) ≡ ∫ t 0 dτ x˙(τ) δ(x − x(τ)), (3.99) againwith implied Stratonovich convention. In analogy toKirchhoff’s rule in the discrete case, the empirical current must be divergence-free in leading order in t, which amounts 64 3.6 Conclusion to being a scalar quantity in one dimension. The fluctuating observable X (t) can be expressed through these densities as X (t) = ∫ dx [p(x)a(x) + j (x)d(x)]. (3.100) For the steady-state average, where ⟨p(x)⟩ = ps(x) and ⟨ j (x)⟩ = js(x), one recovers Eq. (3.93). The probability for joint fluctuations of the fields p(x) and j (x) satisfies a large deviation principle. The corresponding level 2.5 large deviation function is a functional of these fields and reads [52, 68] I[p(x), j (x)] = 1 4T ∫ dx p(x) [ j (x) − jp(x) p(x) ]2 (3.101) where jp(x) ≡ µF (x)p(x) − D∂xp(x) (3.102) is the Fokker-Planck current (3.89) corresponding to the empirical density p(x). The same result for the level 2.5 large deviation function can be obtained as a quadratic approximation of Eq. (3.82) that becomes exact in the continuum limit h → 0, for which api j ≈ 2p(xi)hk0 ≫ ji j, jpi j holds. Thus, even thoughwe focus in the following on discrete state spaces, one can expect that all results that do not make explicit reference to the number of states or “traffic-like” quantities hold also for overdamped Brownianmotion in a continuous state space as a special, limiting case. The generalization to underdamped dynamics, however, is somewhat more subtle, requiring a distinction between reversible and irreversible currents [69]. 3.6 Conclusion In this Chapter we have introduced the mathematical concepts that will be used to describe the fluctuations of thermodynamic currents and, in general, of time-additive observables in stationary Markov processes on discrete state-spaces. We have discussed two different approaches to such a description. First, the scaled cumulant generating function can be calculated as the solution of an eigenvalue problem involving the tilted rate matrix. This approach is particularly useful when low order cumulants of the distribution of a fluctuating observable are to be calculated exactly. Second, large deviation theory provides a description of the exponential decay of the probability of untypical fluctuations. While the corresponding large deviation function of an individual current can generally not be calculated analytically, it is nonetheless possible to calculate the large deviation function associated with joint fluctuations of all possible currents and sojourn times in a formalism that has been dubbed level 2.5 large 65 3 Fluctuations in stationary Markov processes deviation theory. The large deviation function for an individual current can then be gained from the contraction principle. This approach turns out to be particularly useful for proving bounds on current fluctuations. Ultimately, both descriptions, though technically different, are interrelated through Legendre transformations and thus equivalent. Finally, we have shown that systems involving continuous degrees of freedom are covered by the present formalism through a suitable continuum limit. This insight allows us to focus largely on a discrete description throughout this Thesis, without loss of generality. 66 4 Dissipation-dependent bound on current fluctuations In this chapter, we introduce a bound on the fluctuations of current-like observables in non-equilibrium steady states. This bound holds universally for driven systems that are modeled in aMarkovian, thermodynamically consistent fashion, as introduced inChap. 2. Crucially, this bound depends on no other details of the system under consideration than its average rate of total entropy production. Technically, we employ large deviation theory, as reviewed in Sec. 3.3, to formulate this bound as a constraint on the large deviation function capturing the whole spectrum of fluctuations. Evaluated for the typical fluctuations, this bound implies the ther- modynamic uncertainty relation, expressing a universal trade-off between the entropy production, or thermodynamic cost, and the smallness of fluctuations, or precision, of a thermodynamic current [1]. Here, we mainly follow the presentation of the publication [7], where we have first described the dissipation-dependent bound as a conjecture based on numerical evidence. We start with a discussion for systems driven weakly within linear response. In this regime, the bound is typically strongest and can be derived straightforwardly, allowing the reader to build intuition. Beyond linear response, the bound is illustrated numerically for sets of randomly generated Markovian networks. A full mathematical proof of the dissipation-dependent bound has been accomplished by T. Gingrich et al. [8]. We will employ their techniques in Chap. 6 to derive an even stronger, affinity- and topology-dependent bound on current fluctuations. 4.1 General setup for bounds on current fluctuations We consider a Markovian network consisting of N discrete states {i} and allow for tran- sitions with rates ki j ≥ 0 from state i to j. All transitions are taken to be reversible, i.e., ki j > 0 implies k ji > 0. For notational convenience, we exclude the possibility of several transition pathways between any two states and transitions with i = j. The time-dependent probability distribution pi (t) of state i at time t evolves according to the master equation (3.1) with the transition matrixLi j from Eq. (3.2). External forces drive the system into a non-equilibrium steady state, with a stationary distribution psi . Ther- modynamic consistency for the system and the reservoirs providing the driving forces 67 4 Dissipation-dependent bound on current fluctuations ✭❛✮ Figure 4.1: (a) Unicyclic network with affinity A and five states. (b) Multicyclic network with two fundamental cycles, one with three states and affinityA1 and the other with four states and affinity A2. The red dashed lines indicate a cycle with affinity A1 +A2 and five states. (c) Multicyclic network with three fundamental cycles with three states each. The affinities of these cycles areA1,A2, andA3. The red dashed lines indicate a cycle with affinity A1 +A2 and four states. is implemented through the detailed balance relation (2.17) for the pairs of transitions along all links of the network. Following Schnakenberg, we identify a complete set of fundamental cycles {Cβ} within the network [27], as exemplified in Fig. 4.1. For each cycle, we define a time- antisymmetric observable X β (t) that counts cycle completions after time t (the so-called integrated current) and the corresponding average current Jsβ ≡ ⟨Jβ⟩ ≡ ⟨X β (t)⟩/t . (4.1) The average is independent of t for initial conditions drawn from the steady state distribu- tion. The set of cycle currents is complete in the sense that any other time-antisymmetric observable can be written as a linear combination of the X β (t) up to a state function that cannot scale with time. Unlike the currents through individual links of the network, which are constrained through Kirchhoff’s law (3.79), there are no constraints on the realizations of the X β (t). In order to formally define the cycle currents as time-additive observables of the type (3.4), we set ai = 0 and define generalized distances d βi j = −d βji that are constrained to add up to one for every closed loop that completes the cycle once in forward direction. The affinityA β is defined as the change of entropy in the reservoirs upon completion of an individual cycle Cβ. Using the local detailed balance condition (2.17), this affinity is identified as A β ≡ ∑ (i, j)∈Cβ ln ki j k ji = ∑ (i, j)∈Cβ ∑ r (∆Ei jr − Kr ∆Xi jr )/Tr, (4.2) where the first sum runs over the directed links of the cycle and the second one runs over the reservoirs. For isothermal systems with equal reservoir temperature Tr = T , the 68 4.1 General setup for bounds on current fluctuations affinity simplifies to A β = − ∑ (i, j)∈Cβ ∑ r Kr ∆Xi jr /T, (4.3) For instance, recalling unicyclic motor models from Sec. 2.7, one identifies a single cycle with an affinity given by the difference between chemical work and mechanical work per cycle. The cycle current then corresponds to the average motor speed in units of the step size. The fluctuating entropy production Stot(t), as defined in Eq. (2.32), is given by the time-additive observable with increments di j = si j ≡ ln psi ki j psj k ji . (4.4) The average entropy production (2.38) can be expressed with the cycle affinities and cycle currents as σ ≡ ⟨Stot(t)⟩ /t = ∑ β A β Jsβ . (4.5) Adopting a vector notation X (t) for the set of all integrated cycle currents X β (t) the scaled cumulant generating function (3.19) is defined as λ(z) ≡ lim t→∞ 1 t ln ⟨ exp[z · X (t)]⟩ , (4.6) where z is a real vector. As an abbreviation, we will refer to λ(z) simply as the “generating function”. As shown in Sec. 3.1, λ(z) is the largest eigenvalue of the tilted Markov generator (3.13), which reads here Li j (z) ≡ Li j exp ( z · d ji ) , (4.7) where d ji is a vector with components d βji. The probability distribution of the variable Jβ ≡ X β/t satisfies a large deviation principle (see Sec. 3.3) with a large deviation function I (J ). Since stationary cycle currents are typically non-zero in driven systems, we can also define the scaled variable ξ β ≡ X β/(t Jsβ) = Jβ/Jsβ, (4.8) with the corresponding large deviation function h(ξ ) ≡ I ({ξ β Jsβ}). The Legendre- Fenchel transform (3.51) relating the generating function and the large deviation function then reads h(ξ ) = max z ⎡⎢⎢⎢⎢⎢⎣ ∑ β zβ Jsβξ β − λ(z) ⎤⎥⎥⎥⎥⎥⎦ . (4.9) 69 4 Dissipation-dependent bound on current fluctuations The fluctuation theorem for currents (2.53) imposes on the large deviation function function the symmetry −h(ξ ) + h(−ξ ) = ∑ β ξ β JsβA β . (4.10) In terms of the generating function this symmetry reads [39] λ(z) = λ(−A − z). (4.11) An arbitrary fluctuating current is a linear combination X (t) ≡ ∑β cβX β (t) of the cycle currents with coefficients cβ. Its generating function (3.19) follows by evaluating the generating function (4.6) in the direction c, leading to λc (z) = lim t→∞ 1 t ln ⟨ exp[zX (t)] ⟩ = λ(zc). (4.12) In cases where the specific choice of c is not relevant, we will drop the index on the l.h.s. and distinguish generating functions with scalar and vectorial arguments. As a special case, the generating function for an individual cycle current follows as λα (z) ≡ λ(zeα), (4.13) where eα is the unit vector associated with the current in cycle Cα. Generally, the generating function for an individual current does not exhibit a symmetry of the form (4.11) as extensively discussed in [37] (see also [70, 71]). In contrast, the evaluation of λ(z) along the vectorA yields λs (z) ≡ λ(zA) = lim t→∞ 1 t ln ⟨ exp[zStot(t)] ⟩ (4.14) as the generating function of the entropy change. It is symmetric with respect to z = −1/2, which expresses the fluctuation theorem that holds for this observable. The large deviation functions associated with the probability distributions of these variables read hα (ξα) = max z [zJsαξα − λα (z)] (4.15) and, introducing the scaled entropy change s ≡ Stot(t)/(σt) in analogy to Eq. (4.8), hs (s) = max z [zσs − λs (z)], (4.16) respectively. In the following, we distinguish between unicyclic and multicyclic networks of states, illustrated in Fig. 4.1. For unicyclic networks, where there is only a single affinity 70 4.1 General setup for bounds on current fluctuations A ≡ Aα and a single fluctuating current X ≡ Xα, we no longer have to distinguish between the different types of generating functions and can simply write λ(z) ≡ λα (z) = λs (z/A). (4.17) Throughout this Thesis, we will be interested in functions b(z) that bound the gener- ating function λ(z) from below, i.e., b(z) ≤ λ(z) (4.18) for all z. As special cases, the relation (4.18) can be used to extract bounds for individual fluctuating currents, e.g., λα (z) ≥ bα (z) ≡ b(zeα), and the entropy change, λs (z) ≥ bs (z) ≡ b(zA). Such bounds immediately imply upper bounds on the corresponding large deviation functions hα (ξα) ≤ max z [zJsαξα − bα (z)]. (4.19) and hs (s) ≤ max z [zσs − bs (z)]. (4.20) For any generating function the coefficients of the Taylor expansion around z = 0 correspond to the cumulants (3.20). The Fano factor that quantifies the dispersion of the distribution is defined as F ≡ lim t→∞ ⟨ X (t)2 ⟩ − ⟨X (t)⟩2 ⟨X (t)⟩ = 2D Js = λ′′(0) λ′(0) , (4.21) where X (t) is an arbitrary current with average Js and diffusion coefficient D. We denote the Fano factor associated with an individual cycle current by Fα ≡ 2Dα/Jsα and the one associated with the entropy change in the medium by Fs ≡ 2Ds/σ. Since global lower bounds b(z) with b(0) = 0 must share a tangent with λ(z) at z = 0 while having a stronger curvature, every such bound implies with F ≥ b ′′(0) b′(0) (4.22) a bound on the Fano factor. Since λ(0) = 0 holds trivially for all networks, we usually require that our bounds are saturated for z = 0. Hence, if λ(z) is analytic, b(z) must have the same gradient as λ(z). 71 4 Dissipation-dependent bound on current fluctuations 4.2 Linear response regime In the limit of small affinitiesA β, the average current Jsα depends linearly on the affinities, Jsα = ∑ β LαβA β + O ( A2 ) (4.23) with the Onsager matrix Lαβ ≡ ∂J s α ∂A β A=0 = ∂ 2λ ∂zα∂A β (0, 0), (4.24) where we make the dependence of the generating function on the affinities explicit by writing λ(z,A). Following Ref. [39], we derive the fluctuation symmetry (4.11) for zα and A β, which leads to ∂2λ ∂zα∂A β (z,A) = ∂2λ ∂zα∂zβ (−A − z,A) − ∂ 2λ ∂zα∂A β (−A − z,A). (4.25) Letting z andA to zero yields Lαβ = 1 2 ∂2λ ∂zα∂zβ (0, 0), (4.26) which implies the Einstein relation Dα = Lαα (4.27) and the Onsager reciprocal relations Lαβ = L βα . (4.28) Moreover, since the generating function must be convex (see Sec. 3.3), we see that the Onsager matrix must be positive definite. In the region z ≲ O (A), the generating function λ(z) can be expanded as a quadratic form around its center of symmetry, which is, due to Eq. (4.11), located at z = −A/2. The requirement λ(0) = 0 and ∇λ(0) = Js fixes this expansion as λ(z) = ∑ β,γ (zβ +A β/2) L βγ (zγ +Aγ/2) − σ/4, (4.29) where the entropy production σ is given in Eq. (4.5), which becomes for linear response σ = ∑ βγ A βL βγAγ . (4.30) 72 4.3 Beyond linear response: Unicyclic case Evaluating this function for z = zeα yields as generating function related to the individual current λα (z) = zJsα + z 2Lαα . (4.31) The positive definiteness of the matrix G βγ ≡ (LααL βγ − LαβLαγ) [1, supplemental material] (with α fixed) yields∑ β,γ G βγA βAγ = Lαασ − (Jsα)2 ≥ 0. (4.32) Hence λα (z) is bounded from below by λα (z) ≥ zJsα (1 + zJsα/σ). (4.33) Using the Legendre transform (4.19), this bound can be transformed into a bound for the large deviation function hα (ξα) = Lαα 4(Jsα)2 (ξα − 1)2 ≤ σ4 (ξα − 1) 2. (4.34) Since the direction eα can be chosen arbitrarily, the bound (4.33) can be stated in a multidimensional formulation as λ(z) ≥ z · Js (1 + z · Js/σ). (4.35) Equality holds along the line z ∝ A, which corresponds to the generating function λs (z) = λ(Az) associated with entropy change. Within linear response, the large deviation function for the scaled entropy change s is thus given by hs (s) = σ 4 (s − 1)2. (4.36) Eq. (4.34) shows that the knowledge of the average entropy production is sufficient to bound the fluctuations of any individual current in the linear response regime. It should be noted, however, that Eq. (4.34) is still restricted to the fluctuations of the current Jα that are on the order of Jsα. Surprisingly, as we discuss next, this quadratic bound is also valid for both driving and fluctuations beyond the linear response regime. 4.3 Beyond linear response: Unicyclic case The quadratic shape of the generating function for z ≲ O(A) and of the large deviation function for ξ ≲ O(1) can be regarded as a signature of linear response. It arises only for nearly vanishing affinities or for freely diffusing particles, where the linearity between 73 4 Dissipation-dependent bound on current fluctuations ✵ ✵✳ ✶ ✲✶✳ ✲✶ ✲✵✳ ✵ ✵✳ ✕ ❆ ❘ ❲ ✭ ③ ✱ ✁ ✱ ✂ ✮ ✄❂☎ ☎✪ ✵ ✶ ✷ ✸ ✹ ✺ ✲✹ ✲✷ ✵ ✷ ✹ ✻ ❤ ❆ ❘ ❲ ✭ ✘ ✱ ✱ ✁ ✮ ❂ ✂ ✄✪ Figure 4.2: The generating function λ(z) and the large deviation function h(ξ) of the asymmetric random walk for selected affinitiesA (2, 5, 10, 50) and N = 1. Black arrows indicate the direction of increasingA. The quadratic bound for the generating function (4.42) and for the large deviation function (4.43) are shown as black dashed curves. affinity and current persists even for high affinities. Beyond this regime, one universally observes two characteristic changes in the large deviation function [51,59,72]. First, the tails for large values of |ξα | grow no longer quadratically but with a scaling somewhere between linear and quadratic. Second, there is a formation of a “kink” around the value ξα = 0. For finite numbers of states, the large deviation function is still analytic in this region, but it exhibits a significantly enhanced curvature. In the Legendre transformed picture of the generating function λα (z), these two effects show up as a faster than quadratic growth for large z and a pronounced plateau around the minimum of λα (z). This behavior of the generating function is best illustrated with an asymmetric random walk (ARW), as shown in Fig. 4.2. Consider a network consisting of a single cycle with N vertices and affinity A, as shown in Fig. 4.1a. The hopping rates in forward and backward directions k+ and k− are uniform with ln k+ k− = A/N . (4.37) The average current in this model is Js = (k+ − k−)/N and the entropy production is σ = JsA. As derived in Eq. (3.24) (with a different affinity per step), the generating function is given by λ(z) = k+ [ ez/N + e−(z+A)/N − 1 − e−A/N ] (4.38) = JsλARW(z,A, N ), where λARW(z,A, N ) ≡ cosh[(z +A/2)/N] − cosh[A/(2N )](1/N ) sinh[A/(2N )] . (4.39) 74 4.4 Beyond linear response: multicyclic case Similarly, the large deviation function corresponding to the generating function is given by Eq. (3.54), which reads for the present setup h(ξ) = Js hARW(ξ,A, N ), (4.40) where hARW(ξ,A, N ) ≡ N [ ξ arsinh(ξ/ξ0) − ξA/(2N ) − √ ξ20 + ξ 2 + √ ξ20 + 1 ] (4.41) and ξ0 ≡ 1/ sinh[A/(2N )]. As shown in Fig. 4.2, the generating function (4.38) is bounded from below by the quadratic function λ(z,A, N ) ≥ z Js (1 + z/A) = z Js (1 + zJs/σ), (4.42) and the large deviation function is bounded from above by the quadratic function h(ξ) ≤ JsA(ξ − 1)2/4 = σ(ξ − 1)2/4. (4.43) Thus, we see that the quadratic bounds (4.33) and (4.34) hold for the asymmetric random walk not only in the linear response regime but also for arbitrarily strong driving. The bound gets weaker when the affinity per stepA/N increases, i.e., for cycles with large affinity or few states. Likewise, the bound becomes weaker for non-uniform choices of the transition rates along the cycles, which we will illustrate in Sec. 6.3 in the context of an even stronger, affinity dependent bound. Taken together, these observations suggest that the quadratic bounds hold in fact for arbitrary unicyclic networks. 4.4 Beyond linear response: multicyclic case As has been conjectured in Ref. [7] and proven in Ref. [8], the quadratic, dissipation- dependent bound λ(z) ≥ z · Js (1 + z · Js/σ) (4.44) on the generating function (4.6) holds globally for all vectors z and for all types of Markovian networks. In terms of the individual current in a cycle Cα, this bound can be formulated as λα (z) ≥ zJsα (1 + zJsα/σ), (4.45) or as a bound on the large deviation function hα (ξα) ≤ σ(ξα − 1)2/4. (4.46) 75 4 Dissipation-dependent bound on current fluctuations -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -1.5 -1 -0.5 0 0.5 λ α (z )/ σ zJsα/σ Figure 4.3: Generating functions λα (z) for an individual current in fully connected networks with random transition rates. The red and green curves correspond to networks with N = 4 and N = 6 vertices, respectively. The dissipation-dependent bound (4.45) is shown as a black curve. The bound (4.44) becomes the same as (4.35) in the linear response regime. However, (4.44) is also valid beyond this regime where the currents Js are the actual average currents in the steady state, as determined from∇λ(0), which are different from the linear response currents (4.23). Likewise, this bound is no longer restricted to fluctuations on the order of a small stationary current Js. It can therefore be stated without scaling the fluctuating current as I (Jα) ≤ σ(Jα/Jsα − 1)2/4. (4.47) Analogous bounds of the form λ(z) ≥ zJs(1 + zJs/σ), λI (J) ≤ σ(J/Js − 1)2/4. (4.48) hold for all other currents that are linear combinations of the cycle currents with average Js ≡ c · Js, generating function (4.12) and the corresponding large deviation function I (J). In particular for the entropy production we have λs (z) ≥ zσ(1 + z). (4.49) A numerical illustration for this bound is provided in Fig. 4.3. For this illustration, we have generated a set of networks with N vertices and random transition rates ki j = exp [ 5(ϕi j + ϕ ji)/2 + 2θi j ] , (4.50) where ϕi j and θi j are independent Gaussian random numbers with zero mean and variance 1. As the affinity increases, generating functions globally deviate in a positive 76 4.5 The thermodynamic uncertainty relation direction from the parabolic shape. Only at the trivial points z = 0 and z = −A does the generating function in Eq. (4.44) acquire with zero the same value for all networks. For larger networks (N = 6) the left hand side of the plot becomes less populated, since the probability of the vectors eα andA being nearly parallel becomes smaller in higher dimensions. 4.5 The thermodynamic uncertainty relation The local evaluation (4.22) of the parabolic bound (4.48) for an individual current yields the relation Fσ/Js ≥ 2 (4.51) for the Fano factor associated with the current X (t). This result, dubbed the “thermody- namic uncertainty relation”, has first been conjectured in Ref. [1]. Hence, the quadratic bound on the large deviation function is a generalization of this relation. Applied to the current associated with the entropy production in the medium, the bound (4.49) leads to Fs ≥ 2. (4.52) In order to appreciate its implications, it is instructive to re-formulate the rela- tion (4.51). First, writing Dσ/(Js)2 ≥ 1, (4.53) the relation is expressed as a dissipation-dependent bound on the diffusion coefficient associated with X (t). From this relation, measurements of the dispersion and average of an individual current can provide a lower bound on the average entropy production σ ≥ 2(Js)2/D. This bound makes it possible to estimate the entropy production by measuring a single individual current, which is complementary to the methods for inferring the entropy production described inRefs. [73,74]. The bound (4.53) is saturated for the biased Brownian motion X (t) of a colloidal particle with mobility µ that is pulled in one dimension with a constant force f . In this case, we have the entropy production σ = f /T , the single current Js = µ f , and the diffusion coefficient given by the Einstein relation as D = µT . If, however, the particle is pulled along a periodically patterned substrate, there will be finite-time correlations in the current X˙ (t) that lead to deviations of the effective diffusion coefficient D from µT . These deviations must be such that the relation (4.53) is respected. Second, we can define the uncertainty of a process as ε2 ≡ ⟨ (X (t) − ⟨X (t)⟩)2 ⟩ ⟨X (t)⟩2 , (4.54) 77 4 Dissipation-dependent bound on current fluctuations which is a dimensionless quantity characterizing the amplitude of fluctuations compared to the mean. For large time intervals t fluctuations typically average out, such that ε2 decreases like t−1. On the other hand, the thermodynamic cost associated with maintaining the underlying stationary process is quantified by the produced entropy C ≡ σt. Thus, the thermodynamic uncertainty relation can be stated as [1] ε2C ≥ 2, (4.55) imposing a minimal energetic cost that must be paid for small uncertainty in the output of, e.g., an enzymatic reaction. In Chap. 8, we will show that this relation holds in fact for arbitrary time intervals t. For instance, measuring a certain time interval t with a biomolecular, “Brownian” clock in a thermal environment with a precision of ε = 0.01 costs at least 20 000 kBT [75]. 4.6 Dissipation-dependent bound in driven diffusive one dimensional systems The hydrodynamic fluctuation theory for driven diffusive systems in contact with two reservoirs by Bertini et al. [76–78] has been a major development in nonequilibrium statistical physics. This theory leads to a (typically hard) variational problem that, if solved, leads to the exact rate function of the current of particles or heat between reservoirs. The additivity principle derived in [79] is a more direct method that allows for the calculation of the generating function related to the current in driven diffusive systems. For example, this method has been used to calculate this function for the symmetric simple exclusion process (SSEP) [58, 79], the Kipnis-Marchioro-Pressuti (KMP) model [80, 81], and the weakly asymmetric simple exclusion process (WASEP) [82]. These results are valid in the limit of large system size, for which the number of states diverges, and have been verified numerically [81,82]. Even though the dissipation- dependent bound appears to be restricted to the case of a finite number of states, we show that the generating functions obtained from the additivity principle for these three models lies inside our parabolic bound. First we consider the WASEP and the SSEP, which is a particular case of the WASEP. These models are illustrated in Fig. 4.4 and their precise definition can be found in [58]. In the WASEP particles flow from the left reservoir with constant density ϱL to the right reservoir with density ϱR < ϱL. The current of particles in the system is proportional to the entropy production, and the affinity that drives the process out of equilibrium is given by [58, 82] AWASEP = − ln 1 − ϱL ϱL + ln 1 − ϱR ϱR + (L − 1) ln 1 − ν/L 1 + ν/L . (4.56) 78 4.6 Dissipation-dependent bound in driven diffusive one dimensional systems (a) (b) Figure 4.4: Schematic illustrations of the WASEP (a) and the KMP (b) models. For the WASEP, in the bulk the particles jump with rates p ≡ 1/2 + ν/(2L) to the right and q ≡ 1/2 − ν/(2L) to the left, where the SSEP corresponds to ν = 0. At the boundaries particles are exchanged with the reservoirs. The model also has the exclusion principle, i.e., the maximum number of particles in a site is one. For the KMP model energy flows from a hot reservoir at temperature TL to a cold reservoir at temperature TR . In the bulk a randomly chosen pair of sites exchange energy, which is a continuous variable, in such a way that the total energy is conserved. At the boundaries energy is exchanged with the reservoirs. The precise rules of these models can be found in [82] for the WASEP and [81] for the KMP model. The weak asymmetry of the bulk rates, which scales with 1/L, guarantees that in the thermodynamic limit L → ∞ the affinity is finite. In Fig. 4.5, we have calculated the generating function using the additivity principle for the SSEP, as explained in [58], and for the WASEP, as explained in [82]. In both cases the generating functions are inside the dissipation-dependent bound. TheKMPmodel is a driven diffusive system for the transport of energy from a reservoir at temperature TL to a reservoir at temperature TR < TL, as illustrated in Fig. 4.4. A key feature of the KMP model is that there is no dissipation in the bulk. The precise definition of the model can be found in [81]. The heat transfer from the left to the right reservoir is proportional to the entropy production with the affinity given by [81] AKMP = (T−1R − T−1L ). (4.57) The generating function for this model, which is obtained from the additivity principle as explained in [81], also satisfies the dissipation-dependent bound in Fig. 4.5 within the finite support −T−1R < z < T−1L of λ(z). As a consequence, the large deviation function satisfies the corresponding dissipation-dependent bound globally. These results demonstrate that our dissipation-dependent bound is even more univer- sal: it seems to be valid for these driven diffusive systems in the thermodynamic limit, for which the number of states diverges. Another interesting issue will be to explore whether the bound is still valid in the L → ∞ limit if the system undergoes a dynamical phase transition as the KMP model in a ring-like geometry [83]. 79 4 Dissipation-dependent bound on current fluctuations ✵ ✵✳ ✶ ✲✶✳ ✲✶ ✲✵✳ ✵ ✵✳ ✕ s ✭ ③ ✮ ❂ ✛ ✁✂❆ ♣❛r❛❜♦❧✐❝ ❜♦✉♥❞ ❙❙❊P ❲✄ ❙❊P ❑▼P Figure 4.5: Comparison of the dissipation-dependent bound (4.48) with generating functions for driven diffusive systems. For the SSEP the densities of the left and right reservoirs were chosen as ϱL = 0.99 and ϱR = 0.01. For theWASEP the parameters are ν = 10, ϱL = 4/7, and ϱR = 5/18, as in Ref. [82]. For the KMP model the parameters are TL = 2 and TR = 1, as in Ref. [81]. 4.7 Conclusion In this chapter, we have introduced the dissipation-dependent bound on current fluc- tuations as the most universal result covered in this Thesis. The total rate of entropy production, which characterizes the non-equilibrium character of a stationary driven system, is the only quantity that is required to formulate this bound. Using the lan- guage of large deviation theory, the bound is stated as a quadratic lower bound (4.44) on the generating function, or, in the Legendre transformed picture, as a quadratic upper bound (4.46) on the large deviation function. The bound is typically strongest for currents that are proportional to the entropy production. Here, the upper bound on the large deviation function states that spontaneous fluctuations of the current are at least as likely as predicted by a Gaussian distribution withmatching average andwhich displays the symmetry of the fluctuation relation (2.48). In the linear response regime, where the generating function and the large deviation function are themselves quadratic, the bound is saturated for currents proportional to the entropy production. With increasing driving forces, the generating function and the large deviation function develop characteristic deviations from this quadratic shape, which can be best studied for the example of the asymmetric random walk. These deviations typically lead to the dissipation-dependent bound becoming weaker for strong driving. A local evaluation of the dissipation-dependent bound for typical fluctuations yields the so-called thermodynamic uncertainty relation, which establishes a trade-off between precision in any fluctuating current and the total rate of entropy production. 80 5 Application: Bounds on the efficiency of molecular motors Molecular motors are small machines that can transform free energy liberated in a chemical reaction into mechanical work. Their thermodynamic efficiency η is defined as the ratio between the work exerted against an opposing mechanical force f (or torque for a rotary motor) and the free energy consumed in the chemical reaction driving the motor. This efficiency is universally constrained through the second law by one [84–100]. Molecular motors are isothermal machines for which η ≤ 1 replaces the Carnot expression applicable to heat engines [6]. In this Chapter, which is based on the publication [101], we discuss the consequences of the thermodynamic uncertainty relation for the efficiency of molecular motors. As a main result, we show that η is bounded by η ≤ 1 1 + vT/D f , (5.1) where v is the (mean) velocity of the motor, D its diffusion coefficient, and T is the temperature. The intriguing aspect of this bound arises from the fact that v, D, and f are experimentally accessible quantities. No knowledge of the underlying chemical reactions scheme is necessary for applying this bound. Moreover, it holds for a huge class of motor models, arguably essentially for all models that are thermodynamically consistent whether based on discrete states or on a continuous potential as often used in ratchet models [85, 86, 102, 103]. Likewise, it holds for complexes of motors pulling cooperatively a single cargo particle as, e.g., investigated in Refs. [98, 104–106]. Re- cently, the tightness of the thermodynamic uncertainty relation has been investigated for a number of motor models that match experimental data [107]. 5.1 Thermodynamically consistent motor model The molecular motor (complex) has an arbitrary number of internal states {i} that describe distinct conformations. Moreover, a possible change in conformation can be related to the binding of a solute molecule Aρ. Here, ρ label the species, like ATP whose hydrolization to ADP and Pi can drive the motor. The motor steps along a periodic track 81 5 Application: Bounds on the efficiency of molecular motors Figure 5.1: Schematic network of transitions for a molecular motor walking along a periodic track against an external force f . For each periodic interval the internal states i ∈ {1, 2, 3, 4} of the motor can be grouped into a “cell”. of periodicity d, which means that the set of states {i}, called a “cell”, is attached to a discrete linear sequence of spatial positions for the center of the motor, see Fig. 5.1 [36]. Transitions from state i to state j can occur within a cell, i.e., without net advancement of the motor or can be accompanied with a forward or a backward step. In the first case, we denote the rate by k0i j , in the two latter by k + i j and k − i j , respectively. In the spirit of general setup introduced in Sec. 2.2, the upper index distinguishes three sets of transition rates by the concomitant change in the work reservoir providing the mechanical force f . The changes in the chemical reservoirs providing the solute molecules are assumed to be uniquely determined by the initial states and final states i, j of a transition. As usual, the transitions are microscopically reversible, which means that whenever k0i j , 0, k 0 ji cannot vanish either. Likewise, k+i j , 0 implies k − ji , 0 and vice versa. However, there may be transitions within a cell, which do not occur with spatial motion, i.e., k+,−i j = 0 is allowed even if k0i j , 0 and vice versa. Thermodynamic consistency imposes local detailed balance (2.21) as constraints on both types of transitions. First, for transitions within one periodicity cell, k0i j k0ji = exp ⎡⎢⎢⎢⎢⎢⎣ ( Fi − Fj + ∑ ρ bρi j µρ ) /T ⎤⎥⎥⎥⎥⎥⎦ . (5.2) Here, Fi, j are the free energies of the two states, which may comprise contributions from the displacement against the external force for incomplete steps of the motor. If the transition from i to j requires binding a solute of species ρ (like, e.g., ATP), then bρi j = 1. Likewise, if this transition leads to a release of such a molecule, b ρ i j = −1. 82 5.1 Thermodynamically consistent motor model In both cases, the chemical potential µρ of this species enters the expression in the exponent providing a contribution to the total free energy involved in such a transition. If a transition additionally involves a step in the forward direction, against the applied force, then the ratio between such a step and the corresponding backward step becomes k+i j k−ji = exp ⎡⎢⎢⎢⎢⎢⎣ ( Fi − Fj + ∑ ρ bρi j µρ − f d ) /T ⎤⎥⎥⎥⎥⎥⎦ . (5.3) For a fixed applied external force f and externally maintained chemical potentials {µρ}, the motor reaches a steady state with velocity v = ∑ i j psi (k + i j − k−i j )d (5.4) with psi the steady state probability to find the motor in the internal state i. This velocity v = ⟨X (t)⟩ /t = ⟨n+(t) − n−(t)⟩ d/t (5.5) is given by the steady state average of the stochastic displacement X (t), where n+,−(t) are the number of forward and backward steps along the track in time t. In order to apply the thermodynamic uncertainty relation, besides the diffusion coefficient D for the displacement of the motor, as defined in Eq. (3.21), we need the rate of thermodynamic entropy production σ. Summing in Eq. (2.38) over all types of transition, one obtains σ = ∑ i j psi ( k0i j ln k0i j k0ji + k+i j ln k+i j k−ji + k−i j ln k−i j k+ji ) . (5.6) Using the detailed balance conditions (5.2,5.3), the steady state condition (2.20), reading here ∑ j psi (k 0 i j + k + i j + k − i j ) = ∑ j psj (k 0 ji + k − ji + k + ji), (5.7) and the symmetry bρi j = −bρji leads to the expression σ = ∑ i j psi (k 0 i j + k + i j + k − i j ) ∑ ρ bρi j µρ/T − v f /T . (5.8) The chemical work put into the motor arises from reactions of the type∑ ρ r ρν A ρ ⇌ ∑ ρ sρν A ρ, (5.9) 83 5 Application: Bounds on the efficiency of molecular motors where r ρν = 1 if species ρ is an educt of the forward reaction of type ν and s ρ ν = 1 if species ρ is a product of this reaction. In the simplest case of just ATP hydrolysis, as for F1-ATPase [84, 93, 95, 99] or kinesin [91, 108, 109], we have only one net reaction ATP ⇌ ADP + Pi. (5.10) In general, the rate of consumption of species ρ is given by Rρ = ∑ ν γν (r ρ ν − sρν ) (5.11) where γν is the effective net rate of the reaction of type ν with γν < 0 if the reaction goes backward on average. Likewise, Rρ < 0 if, on average, species ρ is rather produced than consumed. Since we assume that all reactions are catalyzed by the motor acting as an enzyme, thus requiring binding and release as introduced above, this net rate can also be written as Rρ = ∑ i j psi (k 0 i j + k + i j + k − i j )b ρ i j . (5.12) The rate of consumption of chemical (free) energy, i.e., the rate with which chemical work is put into the motor, becomes W˙ chem = ∑ ρ Rρµρ. (5.13) Inserting these expressions into the entropy production rate (5.8), we get σ = (W˙ chem − f v)/T = f v(1/η − 1)/T . (5.14) with the thermodynamic efficiency [85] η ≡ f v/W˙ chem = f v/( f v + Tσ). (5.15) 5.2 The bound Inserting the uncertainty relation (4.53) for the displacement current Js = v in (5.15), we obtain our main result stated in the introduction, copied here for convenience as η ≤ 1 1 + vT/D f . (5.16) Remarkably, this universal bound has been derived without any assumptions about the specific molecular mechanism driving the motor. It holds for a single motor as well 84 5.2 The bound 1 10 100 1.5 1.0 0.5 0.0 1000 0 1 2 3 4 5 6 7 1.5 1.0 0.5 0.0 Figure 5.2: Efficiency bound (5.16) for experimental data for kinesin from Ref. [111]. The randomness parameter r = 2D/vd was measured for various ATP-concentrations at constant force f = 3.59 pN (left panel) and for various forces f at constant ATP- concentration cATP = 2 mM (right panel). The shaded regions are labeled with the maximal efficiency η that applies to data points falling into the respective region. Reprinted from Ref. [65], Copyright (2018), with permission from Elsevier. as for complexes of several motors pulling a cargo to which an external force is applied. Hence, it is sufficient to measure the velocity and the diffusion constant of the motor to bound its thermodynamic efficiency. This bound can also be written in terms of the often used randomness parameter r ≡ 2D/vd [88,108,110–112]. With this quantity, the universal bound on efficiency (5.16) reads η ≤ r r + 2T/ f d . (5.17) In Fig. 5.2, we apply this bound to experimental data for kinesin obtained by Visscher et al. [111]. For strong forces and low fuel concentrations, the bound on the efficiency is rather close to unity, whereas for high fuel concentrations and weak forces the maximal efficiency is about 60% or even less. If only one type of chemical reaction, e.g, ATP hydrolysis is involved, the rate of chemical work becomes W˙ chem = RATP∆µ where ∆µ ≡ µATP − µADP − µPi . The rate of ATP-consumption can then be bounded as RATP ≥ ( f + vT D ) v ∆µ . (5.18) This relation makes it possible to infer the minimal rate of ATP consumption by mea- suring the first and second moments of the displacement of the motor. 85 5 Application: Bounds on the efficiency of molecular motors 5.3 Variants So far, we have focused on the thermodynamic efficiency η where the motor runs against an external load. Formotors pulling cargo through a viscous environment with no further external load, often the Stokes efficiency ηS ≡ γ (0)v2/W˙ chem (5.19) is used [103]. It compares the chemical work spent to the (mechanical) power required by a fictitious force to pull the cargo with a bare friction coefficient γ (0) at the same velocity through the viscous medium. Using σ = W˙ chem/T in the uncertainty relation (4.53), one easily gets ηS ≤ γ (0)D/T = D/D(0), (5.20) where the second equality follows from the Einstein relation between the bare friction coefficient and the bare diffusion coefficientD(0) = T/γ (0) of the cargo. Hence the Stokes efficiency is universally constrained by the ratio between the full diffusion coefficient of the motor cargo complex and the bare diffusion coefficient of the cargo. Again, this bound holds independently of all molecular details. So far, we have assumed that the motor can be described by a set of internal states replicated along a spatially periodic track using a discrete master equation dynamics. For motors modeled in a continuous potential possibly changing due to internal (chemi- cal) transitions [85] and for motors pulling big probe particles described by a Langevin equation [96, 113], these universal results remain valid. Formally, one has to discretize the spatial coordinate and introduce (biased) transition rates between neighboring spa- tial positions, as detailed in Sec. 3.5. After such a discretization, which can become arbitrarily fine, one is back at the model investigated above. For a motor that performs chemical work W˙ chemout > 0 by synthesizing molecules with high chemical potential driven by an applied external force f dr > 0 [114] similar results hold with velocity v > 0. The efficiency for this situation reads η ≡ W˙ chemout / f drv = ( f drv − Tσ)/ f drv ≥ 0. (5.21) The inequality (4.53) leads to η ≤ 1 − vT/ f drD. (5.22) Again, the efficiency of this molecular machine can be bounded without measuring the rate at which it produces molecules. 86 5.4 Numerical case study (a) A T P AD P + P i (b) 0 0.2 0.4 0.6 0.8 1 f d /∆µ 0 0.2 0.4 0.6 0.8 1 Th er m od yn am ic ef fic ie nc y η bound ( ∆µ=4 k B T) bound ( ∆µ=10 k B T) bound ( ∆µ=16 k B T) bound ( ∆µ=22 k B T) (c) 5 10 15 20 25 30 35 40 ∆µ [ k B T] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 St ok es e ffi cie nc y ηS bound Figure 5.3: (a) Illustration of the model for an ATP driven motor (blue) elastically coupled to a probe particle (red), to which the external force f is applied. (b) Comparison of the thermodynamic efficiency η as a function of the scaled force f d/∆µ (black line) and the bound (5.16) at various values of ∆µ (colored curves). (c) Comparison of the Stokes efficiency as a function of ∆µ (black) and the bound (5.20) (red). For large ∆µ, where the bound surpasses 1, it is continued as a dashed line. Model parameters were chosen as detailed in Ref. [96], different values of ∆µwere generated by varying the concentration of ATP while keeping the concentrations of ADP and Pi fixed. 5.4 Numerical case study We illustrate the bound for a simple model of the rotary motor F1-ATPase coupled elastically to a colloidal probe particle [96, 115]. As shown in Fig. 5.3a, after mapping the rotary motion to a linear one, the motor performs discrete steps of length d, with each forward step hydrolizing an ATP molecule and each backward step synthesizing an ATP molecule. The stepping rates depend on the chemical potential difference ∆µ of the ATP reaction and the elongation of the linker to the probe particle. The continuous dynamics of the probe particle with friction coefficient γ (0) is modeled as overdamped Brownian motion subject to the external force f and the potential force of the linker [96]. The average velocity v and the diffusion constant D for this model can be calculated numerically from the Master equation that is obtained by finely discretizing the possible 87 5 Application: Bounds on the efficiency of molecular motors elongations of the linker and truncating at large elongations. Due to the tight coupling between the chemical reaction and the motion of the motor, the rate of the chemical work is W˙ chem = ∆µv/d [96], so that the thermodynamic efficiency (5.15) becomes trivially η = f d/∆µ (5.23) while the Stokes efficiency (5.19) is ηS = γ (0)vd/∆µ. (5.24) In Fig. 5.3b, we compare the thermodynamic efficiency (5.23) to the bound (5.16) for forces ranging from 0 to ∆µ/d and selected values of ∆µ. The bound is trivially saturated for f = 0 and for the case f d = ∆µ, where the chemical and themechanical force balance each other, leading effectively to equilibrium conditions with vanishing velocity. Since the uncertainty relation becomes exact within linear response for unicyclic systems [1,8], the bound is also saturated for forces close to ∆µ/d. Similarly, for ∆µ → 0, the bound approaches η for all 0 < f < ∆µ/d. In Fig. 5.3c, the Stokes efficiency (5.24) is compared to the bound (5.20) for a range of ∆µ and constant f = 0. While this bound is also saturated in the linear response regime for ∆µ→ 0, it becomes rather loose for large values of ∆µ. It even surpasses the bound ηS < 1, proven in Ref. [103], when D becomes larger than the bare diffusion coefficient D(0) [116, 117]. Obviously, the bound (5.20) is not useful in this range of parameters. 5.5 Conclusion In this Chapter, we have shown that the thermodynamic uncertainty relation implies a universal bound on the efficiency ofmolecularmotors that depends only on the fluctuating displacement of the motor and an external load force. The sole knowledge of both the average and the dispersion of the experimentally accessible displacement yields a bound that is independent of the underlying mechano-chemical reaction scheme. This result applies to any nano or micro machine operating in an environment of fixed temperature. 88 6 Affinity- and topology-dependent bound on current fluctuations The main strength of the quadratic bound on current fluctuations covered in Chapter 6 lies in its dependence on solely the rate of entropy production – other details of the system under consideration do not have to be known in order to apply this result. However, as a consequence of this universality, the bound can be fairly loose when applied to specific systems, especially when they are driven far away from linear response. One may expect that with more knowledge on the details of a system, stronger bounds on current fluctuations can be devised. Indeed, as we will derive in this chapter, a refined variant of the bound exists, which requires additional information on the driving affinities and the topology of the Markovian network underlying the system. This information will typically be much better accessible than the full set of transition rates that would be required to calculate the spectrum of fluctuations exactly. We have conjectured this affinity- and topology-dependent bound along with the dissipation dependent bound in Ref. [7], where it is referred to as “hyperbolic cosine” bound according to its functional form as a bound on the generating function. We could publish a full proof of the bound in Ref. [118], where we build on the techniques developed for the dissipation-dependent bound in Ref. [8] and construct a suitable cycle- decomposition. For unknown affinities and network topology, the dissipation-dependent bound follows naturally as the weakest instance of the affinity- and topology-dependent bound. Thus, the proof presented in this chapter serves also as a proof of the former. As applied to typical fluctuations, the affinity- and topology-dependent bound implies a refinement of the thermodynamic uncertainty relation, which has first been conjectured in Ref. [119]. 6.1 Setup Just as the dissipation-dependent bound, the affinity- and topology-dependent bound is stated for a stationary Markov process on a finite network of N states, using the same notation as introduced in Sec. 4.1. We recall the definition of the stationary currents along links of the network, jsi j ≡ psi ki j − psj k ji = − jsji, (6.1) 89 6 Affinity- and topology-dependent bound on current fluctuations which are the expectation values of the fluctuating currents ji j , defined in Eq. (3.76) as the net number of transitions in a stochastic trajectory along the link (i, j) divided by the observation time t. An individual transition i → j contributes [6] si j = ln psi ki j psj k ji (6.2) to the entropy production in a NESS, leading to the average total entropy production σ = ∑ i< j jsi j si j . (6.3) We denote by ∑ i< j the sum over those links (i, j) in the network graph for which transitions are possible, i.e., for which ki j > 0 and k ji > 0. The decomposition of networks into cycles is an important concept providing a link between (bio-)physical properties of a system and its description as a stochastic process [22, 27, 120, 121]. After the intuitive introduction given in Sec. 4.1, we here formalize the concept of a cycle decomposition. A cycle Cα of length nα ≥ 3 (smaller cycles are not relevant for our discussion) is defined as a directed, self-avoiding, closed path [ℓ(1) → ℓ(2) → · · · → ℓ(nα) → ℓ(1)] along links of the network with kℓ(n)ℓ(n+1) > 0. We define the directed adjacency matrix χα = ( χαi j ) of such a cycle as χαi j ≡ nα∑ n=1 δℓ(n),iδℓ(n+1), j − δℓ(n), jδℓ(n+1),i, (6.4) which is +1 for links (i, j) where Cα passes in forward direction, −1 for the backward direction and zero otherwise. A set of cycles {Cα} is called complete if any set of currents { ji j } along the links of the network consistent with Kirchhoff’s law (3.79)∑ j ji j = 0 for all i (6.5) can be decomposed into a set of cycle currents {Jα} such that ji j = ∑ α Jα χαi j . (6.6) In particular, the stationary currents jsi j can be decomposed into stationary cycle currents Jsα. The affinity Aα of a cycle Cα is, as in Eq. (6.7), given by Aα = ∑ i< j χαi j ln ki j k ji = ∑ i< j χαi j si j = ∑ (i, j)∈Cα si j, (6.7) which can be used to write the average entropy production (6.3) as in Eq. (4.5) as σ = ∑ α JsαAα . (6.8) 90 6.2 Main result 6.2 Main result With the above definitions at hand, we can now state the central result from our publica- tion [118]. The generating function, defined in Eq. (4.6) as λ(z) = lim t→∞ 1 t ln ⟨ ez·Jt ⟩ (6.9) for the fluctuations of the vector J = (Jα) of cycle currents is bounded from below by λ(z) ≥ σ cosh[(z · J s/σ + 1/2)A∗/n∗] − cosh[A∗/(2n∗)] (A∗/n∗) sinh [A∗/(2n∗)] , (6.10) where A∗/n∗ is defined as the smallest positive value for the affinity per cycle length over all cycles in the network. In a scalar version, this bound implies for the generating function (4.12) of any generic current J, which can be written as an expansion in the cycle currents, the bound λ(z) ≥ σ cosh[(zJ s/σ + 1/2)A∗/n∗] − cosh[A∗/(2n∗)] (A∗/n∗) sinh [A∗/(2n∗)] , (6.11) where Js is the stationary value of J. In particular, for the entropy current ∑ α JαAα we have Js = σ, which further simplifies (6.11). Based on strong numerical evidence we have conjectured the relation (6.10) in Ref. [7]. The proof from Ref. [118] is given in Sec. 6.5 below. As shown in Appendix 6.C, the right hand side of Eq. (6.10) decreases monotonically for decreasing A∗ and all other quantities fixed. Taking the limit A∗ → 0 leads to the quadratic, dissipation-dependent bound (4.44), which Eq. (6.10) is therefore a refinement of. For the typical fluctuations, a local evaluation of the from (4.22) of the bound (6.11) leads to [119] F ≥ A ∗Js n∗σ coth (A∗ 2n∗ ) (6.12) as a bound on the Fano factor. This bound can be used to estimate the number of intermediate states in enzymatic schemes from measurements of the Fano factor in single molecule experiments, as discussed in Ref. [119]. Overdamped Brownian dynamics in a continuous state space is covered by our discrete formalism through the continuum limit described in Sec. 3.5. However, in this limit the affinity per cycle length vanishes with the discretization length h, such that the bound (6.10) degenerates to the dissipation-dependent bound. The latter can also be derived directly for overdamped dynamics from the level 2.5 large deviation function (3.101), as described in Refs. [64, 122]. 91 6 Affinity- and topology-dependent bound on current fluctuations -10 -5 0 5 10 15 -100 -80 -60 -40 -20 0 λ (z )/ Js z SD = 2 1 0.5 0.2 0.1 0.05 0.01 Figure 6.1: Generating function λ(z) scaled by the steady state current Js in unicyclic networks with N = 10 states and fixed affinityA = 100. The color-scale encodes the standard deviation (SD) used for sampling the transition rates (see Ref. [7, Appendix B]). Blue corresponds to a nearly uniform distribution of transition rates and yellow to a broad distribution of transition rates. The lower and upper bounds in Eq. (6.15) are shown as red dashed lines. 6.3 Numerical illustration Physical intuition for the result (6.11) can be gained from numerical case studies. To begin with, consider an arbitrary unicyclic model with N states and periodic boundary conditions. The transition rates are denoted by ki,i+1 = k+i and ki,i−1 = k − i , (6.13) where i = 1, 2, . . . , N . A fixed affinity A implies the constraint∏N i=1 k + i∏N i=1 k − i = eA (6.14) on the transition rates. Different choices of the transition rates that fulfill this restriction can lead to different generating functions, as shown in Fig. 6.1 for a randomly generated set of unicyclic networks. In particular, if the transition rates are uniform, i.e., k+i = k + and k−i = k −, as considered in Sec. 4.3, the generating function divided by the average cycle current λ(z)/Js becomes λARW(z,A, N ), which is given in Eq. (4.39). The opposite extreme choice for the transition rates is the case where the network behaves effectively like there was only one link between states (N = 1) and all the affinity is concentrated in this single link. In this case, λ(z)/Js becomes λARW(z,A, 1), which fulfills λARW(z,A, 1) ≥ λARW(z,A, N ). 92 6.3 Numerical illustration -1.5 -1 -0.5 0 0.5 -8 -7 -6 -5 -4 -3 -2 -1 0 1 λ (z )/ Js z A = 7 A = 5 A = 3 Figure 6.2: Generating function λ(z) scaled by the steady state current Js in randomly generated unicyclic networks with N = 10 for three families of affinities A. The black curves refer to the lower and upper bound from Eq. (6.15). These considerations suggest for unicyclic networks the bounds λARW(z,A, 1) ≥ λ(z)/Js ≥ λARW(z,A, N ). (6.15) The lower bound corresponds to Eq. (6.11), which becomes particularly simple in the unicyclic case, where A∗ and n∗ are uniquely determined by the single cycle affinity A and length N , and where there is only a single cycle current that relates to the entropy production rate as σ = AJs. The upper bound in Eq. (6.15) is specific for unicyclic networks and is proven in Appendix 6.A. In Fig. 6.2, we show the bounds in Eq. (6.15) for different values of the affinity. For smaller A the bounds are closer to each other. In the linear response regime, up to quadratic order in z, they become the same and equal to the dissipation-dependent bound in Eq. (4.42). The bounds in Eq. (6.15) lead to the bounds coth (A 2 ) ≥ F ≥ 1 N coth ( A 2N ) (6.16) on the Fano factor defined in Eq. (4.22). The lower bound is an affinity dependent bound on the Fano factor that has been obtained in [1]. In the formal limit A → ∞ it becomes F ≥ 1/N . This bound for formally divergent affinity is a key result in statistical kinetics [112] as it allows for an estimate on the number of states N from measurements of the Fano factor. The upper bound on F in Eq. (6.16) is a new result. It is a refinement of the known result F ≤ 1, which is also valid in the limit A → ∞ [112]. Next, we illustrate the bound (6.11) for multicyclic networks. In order to explain how to identify A∗/n∗ we consider the network of states in Fig. 6.3a. We arbitrarily choose 93 6 Affinity- and topology-dependent bound on current fluctuations (a) -0.2 0 0.2 -1 0 λ s( z) /σ z (b) -0.2 0 0.2 -1 0 λ s( z) /σ z (c) -0.2 0 0.2 -1 0 λ s( z) /σ z (d) -0.2 0 0.2 -1 0 λ α (z )/ σ zJsα/σ (e) Figure 6.3: Generating functions for entropy change (b-d) and randomly selected individual currents Jα (e) for the network shown in panel (a). The affinities are A = (1, 2, 3) in (b), A = (11, 12, 13) in (c) and (e), and A = (−11, 12, 13) in (d). The black curves correspond to the dissipation-dependent bound and the red dashed curves correspond to the affinity- and topology-dependent bound. The random transition rates were generated as explained in Ref. [7, Appendix B]. In panel (b)A∗/n∗ = 1/3, in panels (c) and (e)A∗/n∗ = 11/3, and in panel (d)A∗/n∗ = 1/4. For small values of A∗/n∗, as in panels (b) and (d), the parabolic and hyperbolic cosine bound are closer to each other. the cycles (1, 4, 2, 1) with affinity A1, (2, 4, 3, 2) with affinity A2, and (1, 3, 4, 1) with affinity A3 as the three fundamental cycles. Any other cycle in the network is just a composition of these fundamental cycles; for example, the cycle (1, 4, 3, 2, 1), which is marked with a red dotted line in Fig. 6.3a, with affinity A1 +A2 is the sum of the first and second fundamental cycles. If the affinities areA = (1, 2, 3) then the cycle with the smallest affinity per number of states is the fundamental cycle (1, 4, 2, 1). Therefore, in this case A∗ = 1 and n∗ = 3. If the affinities areA = (−11, 12, 13), where the negative affinity is with respect to counter-clockwise direction in Fig. 6.3a, the cycle with minimal affinity per number of states is (1, 4, 3, 2, 1). In this case A∗ = 1 and n∗ = 4. The bound in Eq. (6.11) for the current associated with the entropy production can be interpreted as follows. Given a network of states with fixed affinities, the transition rates that lead to the smallest possible generating function λs (z)/σ are those for which the cycle with smallestA/n dominates the network. This cycle dominates the network if the 94 6.4 Uniform cycle decomposition transition rates within the cycle are large, transition rates to leave the cycle are small, and transition rates to return to the cycle are large. With this choice for the transition rates, the multicyclic network is effectively a unicyclic network with affinity A∗ and number of states n∗, for which the lower bound in Eq. (6.15) holds. Any other choice of rates will just add cycles with a smaller affinity per number of states, which cannot decrease fluctuations. In Fig. 6.3 we illustrate the bound (6.11) along with the dissipation-dependent bound for sample realizations of a four-state network with various combinations of cycle affini- ties. If the cycle relevant for the bound has a rather small affinity per number of states, which is often the case in a large network of states, the bound (6.11) is only slightly stronger than the dissipation-dependent bound (4.45), as visible in Fig. 6.3b and 6.3d. An often tighter bound for this situation is derived in Sec. 7.1. Our numerics indicates that an affinity-dependent upper bound on the generating function in the multicyclic case does not exist. For fixed affinities, the generating function can become arbitrarily close to the trivial bound λs (z) < 0 for −1 < z < 0, visible in Fig. 6.3b-d. 6.4 Uniform cycle decomposition For multicyclic networks, the choice of a minimal complete set of cycles is not unique. We make use of this ambiguity to expand the stationary current jsi j in a specific set of cycles {Cβ}. This set is constructed in a way that every cycle Cβ contributing to the stationary current, i.e., for which Jsβ , 0, satisfies the following properties, which will turn out to be essential for our proof of the affinity dependent bound (6.10): (i) The links of every cycle are aligned with the stationary current, i.e., sgn χ βi j = sgn j s i j = sgn si j (6.17) for all i, j and β. We call such cycles uniformwith respect to jsi j . The equality with sgn si j follows directly from the definitions (6.1) and (6.2) and shows via Eq. (6.7) that uniform cycles have positive affinity. (ii) All stationary cycle currents are strictly positive, Jsβ > 0 for all β. It can be easily checked that an arbitrary cycle decomposition, for example the de- composition into fundamental cycles by Schnakenberg [27], does not necessarily meet these conditions. A construction of a cycle decomposition that satisfies at least condi- tion (ii) was introduced by J. MacQueen [123], its physical relevance and proof is also discussed in Ref. [124]. This decomposition, however, refers only to networks with irre- versible transitions. The network of our present setup with reversible transitions could be mapped on such a network by replacing every link by two irreversible links with 95 6 Affinity- and topology-dependent bound on current fluctuations 1 1 3 4 2 2 2 3 3 2 2 1 2 2 2 2 Figure 6.4: Example for a decomposition of a 5-state network into a complete set of three independent uniform cycles. Links are labeled by the strength of the currents j (β)i j , whose direction is marked by arrows. For each of the three iteration steps, a link with minimal current is labeled in red and the uniform cycle Cβ including this link is shown as a yellow polygon with red border. antiparallel direction, yet this procedure does not necessarily lead to a decomposition satisfying condition (i). Here, we present a variation of the algorithmofMacQueen for networkswith genuinely reversible transitions, which generates a set of cycles satisfying both conditions (i) and (ii), as exemplified in Fig. 6.4. • Initialization. Set j (1)i j ≡ jsi j . • Iteration over β ≥ 1. Locate a link (i∗, j∗) with minimal positive current j (β)i∗ j∗ and set Jsβ ≡ min i, j | j (β)i j >0 j (β)i j = j (β) i∗ j∗ . (6.18) Construct Cβ as a closed self avoiding path starting with the link (i∗, j∗) and passing only along links in the direction for which j (β)i j > 0. This path is not necessarily unique, one of several such paths can be chosen freely as Cβ. Using the directed adjacency matrix χ βi j corresponding to Cβ, as defined in Eq. (6.4), set j (β+1)i j = j (β) i j − Jsβ χ βi j . (6.19) • Terminate when j (β+1)i j = 0 for all i, j. By construction, the currents j (β)i j satisfy Kirchhoff’s law (6.5) for every β. As a consequence, it is always possible to construct an appropriate uniform cycle Cβ following the direction of j (β)i j : For every node where there is a way in, there must also be a way 96 6.5 Proof a) b) c) Figure 6.5: Illustration of a uniform cycle Cβ starting with the link (i∗, j∗). The direction of j (β)i j is indicated by the arrows, the (attempted) cycle is shown as solid blue lines. Neither dead ends (a) nor “dead-cycles” (b) are consistent with Kirchhoff’s law (6.5). out, which rules out dead ends (see Fig. 6.5a). Similarly, it is not possible to run into a cycle without exit that does not include the starting node i∗ (see Fig. 6.5b). On a finite set of states, every self avoiding path must reach the starting point i∗ again after at most N steps. Typically, there is more than one cycle starting with (i∗, j∗) meeting these conditions. Any of these cycles can be selected as Cβ. Since we use in every iteration the minimal current as new cycle current, the individual currents j (β+1)i j are either zero or have the same sign as jsi j . Thus, every cycle that is uniform with respect to the currents j (β)i j is also uniform with respect to j s i j , as required by condition (i). The matrices χ β i j are linearly independent, since for every β the link (i∗, j∗) is no longer contained in all the subsequent cycles. By construction, the algorithm leads to the decomposition jsi j = ∑ β Jsβ χ β i j . (6.20) Typically, the algorithm terminates when β reaches the number of fundamental cycles. It cannot terminate later since all cycles Cβ are linearly independent. It may happen that the algorithm terminates earlier, so that the set of cycles {Cβ} is not complete, i.e., it cannot be used to represent arbitrary currents different from the stationary current. In this case the set {Cβ} can be completed by adding further linear independent non-uniform cycles. The stationary current of these cycles is always zero. 6.5 Proof Equipped with the decomposition in uniform cycles, we can now prove the affinity dependent bound (6.10) on the large deviation function. As a starting point we use, building on Ref. [29] and analogously to Refs. [8, 64], the universally valid bound for the large deviation function of the distribution of the vector j = ( ji j ) of all fluctuating currents along links, I ( j) ≤ ∑ i< j ⎡⎢⎢⎢⎢⎣ √ (asi j ) 2 + ( jsi j ) 2 − √ (asi j ) 2 + ( ji j )2 + ji j *,arsinh ji j asi j − arsinh jsi j asi j +- ⎤⎥⎥⎥⎥⎦ , (6.21) 97 6 Affinity- and topology-dependent bound on current fluctuations with asi j ≡ 2 √ psi p s j ki j k ji. This relation holds for all values of j that are consistent with Kirchhoff’s law (6.5), otherwise I ( j) = ∞. Eq. (6.21) follows from the exact expression (3.82) for the level 2.5 large deviation function I ( j, p) for the joint distribution of the fluctuating current j and the fluctuating density p = (pi). Using the contraction principle (3.85) with p = ps as a trial function yields immediately Eq. (6.21) as bound. This particular choice of the trial function ensures that the bound is saturated for the most likely, stationary current and that one can ultimately identify physically meaningful steady-state averages such as the entropy production σ in the bound. At first, we assume for simplicity that the stationary current jsi j is non-zero along all links with ki j > 0. Then, we can rewrite the bound (6.21) as I ( j) ≤ ∑ i< j jsi jψ( ji j/ j s i j, si j ), (6.22) where the function ψ(ζ, s) is defined as ψ(ζ, s) ≡ ζ arsinh(ζ/b) − ζ s/2 − √ b2 + ζ2 + √ b2 + 1 (6.23) with b ≡ [sinh(s/2)]−1. (6.24) The functionψ(ζ, s) is the large deviation function (4.41) for the current in an asymmetric random walk with uniform transition rates, with affinity per step s, and with stationary current Js = 1 [35]. For this highly symmetric network, there is only one possible fluctuating current J and psi = p ∗ i (J) = 1/N holds independently of J, so that Eq. (6.22) becomes an equality. For positive s and fixed ζ , ψ(ζ, s) is concave in s, moreover, ψ(ζ, s)/s decreases monotonically with increasing s, as proven in Appendix 6.B and Appendix 6.C, respectively. These essential properties are demonstrated in Fig. 6.6. The Legendre transform of ψ(ζ, s) is λ(y, s) ≡ max ζ [yζ − ψ(ζ, s)] = [cosh(y + s/2) − cosh(s/2)]/ sinh(s/2). (6.25) Using the decomposition in the uniform cycles {Cβ} from Eq. (6.20), for which all si j along a cycle become positive, we can write Eq. (6.22) as I ( j) ≤ ∑ β Jsβ ∑ (i, j)∈Cβ ψ( ji j/ jsi j, si j ). (6.26) Next, we set ji j = ξ jsi j for all links (i, j) of the network. This choice simplifies the following steps and ensures that the fluctuating currents ji j are consistentwithKirchhoff’s 98 6.5 Proof 0 0.5 1 1.5 2 −2 −1 0 1 2 3 4 ψ (ζ ,s )/ s ζ s → 0 s → ∞ 0 5 10 15 20 0 5 10 15 ψ (ζ ,s ) s ζ = −1.6 ζ = −1 ζ = 0 ζ = 1 ζ = 5 Figure 6.6: Left: The function ψ(ζ, s)/s as a function of ζ for values of s ranging from 1 (violet) to 101.2 (red) with logarithmic spacing. For fixed ζ and s > 0, ψ(ζ, s)/s decreases monotonically in s, as required in Eq. (6.28). The limiting curves for lims→∞ ψ(ζ, s)/s = (|ζ | − ζ )/2 and for lims→0 ψ(ζ, s)/s = (ζ − 1)2/4 are shown in black. Right: The function ψ(ζ, s) as a function of s for values of ζ ranging from −2 (violet) to 6 (red). Selected values of ζ are shown as solid lines, in particular ψ(−1, s) = s and ψ(1, s) = 0. Obviously, ψ(ζ, s) is concave for s ≥ 0. law (6.5). Having ensured via the uniform cycle decomposition that all si j are positive, we can apply Jensen’s inequality to the second sum in Eq. (6.26), in which the summation index runs over the nβ links of the cycle Cβ. Since the Jsβ are positive, this procedure leads to I (ξ js) ≤ ∑ β Jsβnβ ∑ (i, j)∈Cβ 1 nβ ψ(ξ, si j ) ≤ ∑ β JsβA β (nβ/A β) ψ(ξ,A β/nβ), (6.27) where we have used Eq. (6.7). For a coarser bound, we identify the uniform cycle with minimal A β/nβ ≡ A∗/n∗, which is positive as explained below Eq. (6.17). Using the monotonic decrease of ψ(ζ, s)/s, the positivity of JsβA β, and Eq. (6.8), we obtain I (ξ js) ≤ σ (n∗/A∗) ψ(ξ,A∗/n∗). (6.28) Finally, we adopt a more elegant formulation of the large deviation function taking the cycle currents Jβ as argument instead of ( ji j ) = j, i.e., we write I(J ) ≡ I *., ∑ β Jβχ β +/- . (6.29) We thus avoid currents that are inconsistent with Kirchhoff’s law, for which the large deviation function would become infinite. The generating function corresponding to the 99 6 Affinity- and topology-dependent bound on current fluctuations large deviation function (6.29) turns out to be bounded by λ(z) = max J ⎡⎢⎢⎢⎢⎢⎣ ∑ β zβ Jβ − I( j) ⎤⎥⎥⎥⎥⎥⎦ ≥ max J=ξJs [ z · Jsξ − I (ξ js)] ≥ max ξ [ z · Jsξ − σ (n∗/A∗) ψ(ξ,A∗/n∗)] = σ A∗/n∗ λ(z · J s(A∗/n∗)/σ,A∗/n∗) = σ cosh[(z · Js/σ + 1/2)A∗/n∗] − cosh[A∗/(2n∗)] (A∗/n∗) sinh [A∗/(2n∗)] ≡ B(z · Js,A∗/n∗, σ), (6.30) where we have used Eq. (6.28) for the third line and Eq. (6.25) for the fifth line. Although the uniform cycle decomposition was essential in the steps leading to this result, its final statement is no longer dependent on the specific choice of independent cycles, as can be seen from the following argument. In a different complete set of cycles the cycle currents are represented as J˜α = ∑ β Cαβ Jβ with some transformation matrix Cαβ. The corresponding generating function then reads λ˜( z˜) = lim t→∞ 1 t ln ⟨ exp ⎡⎢⎢⎢⎢⎣ ∑ α z˜α J˜αt ⎤⎥⎥⎥⎥⎦ ⟩ = lim t→∞ 1 t ln ⟨ exp ⎡⎢⎢⎢⎢⎢⎣ ∑ α,β z˜αCαβ Jβt ⎤⎥⎥⎥⎥⎥⎦ ⟩ = λ ( CT z˜ ) (6.31) and satisfies the bound λ˜( z˜) ≥ B(CT z˜ · Js,A∗/n∗, σ) = B( z˜ · J˜s,A∗/n∗, σ), (6.32) which finally proves Eq. (6.10) for arbitrary cycle decompositions. It should be noted that, a priori, A∗/n∗ still refers to the minimal A β/nβ in a set of uniform cycles. However, in a possible application, where the network topology and the cycle affinities are known but individual transition rates are not known, the more general result (6.10), where A∗/n∗ refers to the minimum over all cycles, might lead to a weaker bound that is more useful. Since uniform cycles always have affinityA β > 0, the minimization can be restricted to cycles with positive affinity. Thus, in the physically important case of a network containing cycles with zero affinity, Eq. (6.10) provides a bound that is stronger than the dissipation-dependent bound that is obtained for A∗/n∗ → 0. This fact is in accordance with the insight from Ref. [125], stating that coarse graining such cycles with zero affinity has little effect on the fluctuations of the entropy production. 100 6.6 Vanishing currents and equilibrium 6.6 Vanishing currents and equilibrium So far, we have considered only the case of non-vanishing stationary currents jsi j along all links. If this is not the case for some links (although ki j > 0), we have to split the bound for the large deviation function in Eq. (6.21) in two parts as I ( j) ≤ ∑ i< j | jsi j,0 jsi jψ( ji j/ j s i j, si j )+ ∑ i< j | jsi j=0 ⎡⎢⎢⎢⎢⎣ ji j arsinh ji j asi j − √ (asi j ) 2 + j2i j + a s i j ⎤⎥⎥⎥⎥⎦ . (6.33) Evaluating this bound along currents j = ξ js, the second sum vanishes while the first sum leads along the same lines as before to the bound (6.30). However, in the equilibrium case, where all stationary currents vanish, this procedure leads merely to the trivial statement λ(z) ≥ 0. While typical fluctuations in equilibrium systems can be well described within linear response theory, rare fluctuations exhibit many features akin to systems far from equilibrium [126]. A non-trivial bound for these rare fluctuations in equilibrium systems, similar to the one derived in Ref. [7] for unicyclic networks, is obtained by letting the affinityAα of a single fundamental cycle go to zero while keeping the affinities of the other fundamental cycles fixed at zero. Then, none of the cycles can have an affinity smaller than Aα, so that A∗ = Aα. The length of the relevant cycle n∗ can be bounded by the total number of states N ≥ n∗. Due to the Einstein relation (4.27), the stationary current Jsα is given for small Aα by Jsα = DαAα + O(A2α), (6.34) where Dα is the diffusion coefficient associated with the current Jα. With the entropy production (6.8) reducing to σ = JsαAα = DαA2α, we thus obtain as bound on the generating function for the equilibrium fluctuations of Jα λα (z) ≥ limAα→0 B(zDαAα,Aα/N, DαA 2 α) = 2N2Dα (cosh(z/N ) − 1). (6.35) This bound is saturated for z ≪ 1, where λα (z) depends quadratically on z according to the Gaussian distribution of typical fluctuations within linear response. The probability of extreme fluctuations beyond linear response deviates the more from this Gaussian shape the smaller the number of states in the network is. 6.7 Conclusion In this Chapter, we have presented the bound (6.10) as a major refinement of the dissipation-dependent bound (4.44) on the spectrum of current fluctuations in a Markov network. This refined bound depends additionally on the minimumA∗/n∗ of the affinity per cycle length in any cycle of the network, which is accessible when the network 101 6 Affinity- and topology-dependent bound on current fluctuations topology and cycle affinities are known. Just as the dissipation-dependent bound, the refined bound is universally valid for arbitrary currents in arbitrary networks and entails a bound on the Fano factor in enzyme kinetics, previously conjectured in Ref. [119]. The bound (6.10) on the generating function of current fluctuations has the shape of a hyperbolic cosine, which is also the shape of the generating function of the current in an asymmetric random walk with uniform rates. Indeed, for this simple system the bound is saturated. Making the rates non-uniform and adding further cycles to the network typically drives the generating function away from this bound, as we have illustrated numerically. As a preliminary for the proof of the bound, we have shown that for every Markov network it is possible to construct a set of independent cycles such that all affinities and stationary cycle currents are positive and that all stationary currents along links are aligned with the direction of the cycles. We call this a decomposition in uniform cycles. Based on such a decomposition of the stationary currents, the bound (6.10) can be proven within the framework of level 2.5 large deviations. As a corollary, we have also proven a universal bound on equilibrium fluctuations, given in Eq. (6.35), that depends only on the diffusion coefficient and on the number of states. 102 Appendices to Chapter 6 6.A Proof of the upper bound on the generating function for unicyclic networks For unicyclic networks with N states, the tilted Markov generator (3.13) has the tridiag- onal shape L(z) = *............, −r1 k˜21 0 0 . . . k˜N1 k˜12 −r2 k˜32 0 . . . 0 0 k˜23 −r3 k˜43 . . . 0 0 0 k˜34 −r4 . . . ... ... ... . . . . . . . . . k˜N,N−1 k˜1N 0 0 . . . k˜N−1,N −rN +////////////- (6.36) with k˜i j ≡ ki jezdi j . The displacements di j must satisfy di j = −d ji. If the observable of interest is the number of turnovers, the displacements must add up to the cycle affinity d12 + d23 + · · · + dN1 = 1. It should be kept in mind that any statistical quantity in the long time limit (in particular the generating function and the rate function) must be independent of the specific choice of the individual di j . The characteristic polynomial associated with the matrix (6.36) reads χ(z, x) ≡ det(L(z) − x1N ) = ∑ π (−1)π N∏ i=1 (Liπ(i) (z) − xδiπ(i)), (6.37) where the sum runs over all permutations π of the indices i = 1, . . . , N and 1N is the N × N identity matrix. We identify 0 ≡ N and N + 1 ≡ 1 for the indices of matrix entries. There are two types of terms in (6.37) that contain a specific rate k˜i+1,i: the contribution from the next row i + 1 can either be k˜i,i+1 or k˜i+2,i+1. For the former type the z-dependence cancels out due to di+1,i = −di,i+1 and we end up with the constant factor k˜i+1,i k˜i,i+1 = ki+1,iki,i+1. Terms of the latter type must also contain k˜i,i−1 as the only possible contribution from the previous column i − 1. Iteratively, we see that there can be only one term of this type, namely the one that contains all forward transitions k˜12 k˜23 . . . k˜N1 = k12k23 . . . kN1ez ≡ Γ+ez . (6.38) 103 6 Affinity- and topology-dependent bound on current fluctuations ① P✭①✮ ② ✶ ① ✶ Figure 6.A.1: The polynomial P(x) for a generic unicyclic network with N = 6 states. The tangent at x = 0 and the values y = y1 and x = x1 are shown as dashed lines. An analogous argument can be set up for the lower off-diagonal of the matrix with the z-dependent term k˜21 k˜32 . . . k˜1N = k21k32 . . . k1Ne−z ≡ Γ−e−z = Γ+e−(z+A) . (6.39) All other terms in the determinant (6.37) do not depend on z and we can write χ(z, x) = (−1)N+1[Γ+ez + Γ−e−z − (Γ+ + Γ−) − P(x)] (6.40) with some polynomial P(x) that is independent of z and the specific choice of the di j . The alternating prefactor is due to the fact that the permutations associated with the terms (6.38) and (6.39) are either odd or even, depending on the number of states N . The generating function is thus given by λ(z) = P−1 ( Γ+ez + Γ−e−z − Γ+ − Γ−) = P−1 ( 2 √ Γ+Γ−[cosh(z +A/2) − cosh(A/2)] ) , (6.41) where the function P−1(y) returns the root of the polynomial P(x)− y that has the largest real part. Due to the Perron-Frobenius theorem, this root must be real for all arguments occurring in (6.41), i.e., for all y ≥ y1 ≡ 2 √ Γ+Γ−[1− cosh(A/2)]. The root associated with the minimal argument y1 is x1 ≡ P−1(y1) = minz λ(z) = λ(−A/2). Obviously, the polynomial P(x) (see Fig. 6.A.1) has the properties P(0) = χ(0, 0) = 0 and lim x→∞ P(x) = (−1) N lim x→∞ χ(z, x) = lim x→∞(−1) N det(−x1N ) = +∞. (6.42) 104 6.B Proof of the concavity of ψ(ζ, s) in s Since the matrix L(−A/2) can be brought to a symmetric form by choosing di j = ln ( ki j/k ji ) /A, the corresponding characteristic polynomial P(x) − y1 has only real roots xi with x1 denoting the largest one. The second derivative of P(x) is P′′(x) = d2 dx2 ⎡⎢⎢⎢⎢⎣y1 + N∏ i=1 (x − xi) ⎤⎥⎥⎥⎥⎦ = N∑ i=1 N∑ j=1 j,i N∏ ℓ=1 i,ℓ, j (x − xℓ). (6.43) For x > x1 this expression is positive so that P(x) is convex. As a consequence, the inverted function P−1(y) is concave for the relevant arguments y > y1. Hence it satisfies P−1(y) ≤ (P−1)′(0) y (6.44) with equality for y = 0. This relation leads to the upper bound λ(z) ≤ 2√Γ+Γ−(P−1)′(0) [cosh(z +A/2) − cosh(A/2)] (6.45) holding with equality for z = 0. From Eq. (6.41), we see that the prefactor in this bound is equal to the stationary current Js = λ′(z) = 2 √ Γ+Γ−(P−1)′(0) sinh(A/2), (6.46) which leads to the upper bound in Eq. (6.15). 6.B Proof of the concavity of ψ(ζ, s) in s The derivative of ψ(ζ, s) in Eq. (6.23) can be written as ∂sψ(ζ, s) = 1 2 √ b2 + ζ2 √ b2 + 1 − b 2 2 − ζ 2 . (6.47) As a function of b = [sinh(s/2)]−1 (which decreases monotonically in s), this expression increases monotonically for b > 0, since ∂b∂sψ(ζ, s) = b 2 *., √ b2 + ζ2 b2 + 1 + √ b2 + 1 b2 + ζ2 − 2+/- ≥ 0 (6.48) (note that x + 1/x ≥ 2 for x > 0). Therefore, ∂sψ decreases monotonically in s and ∂2s ψ(ζ, s) < 0 (6.49) holds for all ζ ∈ R and s > 0. 105 6 Affinity- and topology-dependent bound on current fluctuations 6.C Proof of the monotonic decrease of ψ(ζ, s)/s in s The monotonic decrease of ψ(ζ, s)/s in s for s > 0 is equivalent to the monotonic increase of its Legendre transform µ(z, s) ≡ max ζ [ζ z−ψ(ζ, s)/s] = λ(z f , s)/s = cosh[(z + 1/2)s] − cosh(s/2) s sinh(s/2) ≡ A(z, s) B(s) . (6.50) The derivative ∂sµ(z, s) is non-negative if 0 ≤ C(z) ≡∂sA(z, s) B(s) − A(s) ∂ f B(s) = [(z + 1/2) sinh((z + 1/2)s) − (1/2) sinh(s/2)] s sinh(s/2) − [cosh((z + 1/2)s) − cosh(s/2)] [sinh(s/2) + (s/2) cosh(s/2)] . (6.51) Equating ∂zC(z) to zero leads to 2(z + 1/2) tanh(s/2) = tanh[(z + 1/2)s], (6.52) which has the three solutions z = −1,−1/2, 0. The stationary points z = 0 and z = −1 are minima since ∂2zC(z = 0) = s2 2 (sinh s − s) > 0 (6.53) and C(z) is symmetric with respect to −1/2. Thus, the center of symmetry at z = −1/2 is a maximum of C(z) and the global minimum of C(z) is given by C(0) = C(−1) = 0. 106 7 Activity-dependent and asymptotic bound In this Chapter, we explore further universal bounds on fluctuations of time-additive observables that depend on few characteristic features of the system and its steady state. These bounds are complementary to the dissipation-dependent bound, as they can become strong in cases where that bound is rather weak, such as for strong driving far away from linear response and for extreme fluctuations far away from typical fluctuations. The first main result in this Chapter is a bound that depends on the overall activity of the system, i.e., the total number of transitions per unit time. As a bound on the generating function, this bound has the shape of an exponential function. Evaluated for the typical fluctuations, it implies a relation akin to the thermodynamic uncertainty relation, where the activity replaces the entropy production rate. The second main result is a bound which becomes strong for extreme fluctuations as described by large values of the parameter in the generating function. It captures the asymptotic behavior of the generating function and depends essentially on the topology of the Markov network. 7.1 Activity-dependent bound In Ref. [7], we have derived an activity-dependent, exponential bound on the generating function for currents in a stationary state. The proof therein builds on inequalities for the Perron-Frobenius eigenvalue of the tilted Markov generator [127]. A simpler proof using the formalism of level 2.5 large deviations has been devised by Garrahan [128] and will be followed here. We consider a Markovian network with a finite set of states and use the same notation as introduced in Sec. 4.1. The observable of interest, X (t), is specified according to Eq. (3.4) with a distance matrix di j and ai = 0. Unlike for the previous bounds, we do not require di j to be antisymmetric. Hence, the result to follow holds not only for current observables but also for counting observables [128]. The scaled fluctuating observable is given by J = X (t)/t = ∑ i j µi jdi j with the empirical flow defined in Eq. (3.69). Its expectation value in the stationary state is Js = ⟨X (t)⟩ /t = ∑ i j µsi jdi j = ∑ i j psi ki jdi j . (7.1) 107 7 Activity-dependent and asymptotic bound The large deviation function for J can be obtained through contraction from the level 2.5 large deviation function (3.75) as I (J) = min∑ i j µi jdi j=J I ({pi}, {µi j }). (7.2) The ansatz pi = psi and µi = ξµ s i with ξ ≥ 0 satisfies Kirchhoff’s law and leads to a scaled current J = ξ Js. The incomplete contraction (3.85) using this ansatz leads to the bound h(ξ) ≡ I (J = ξ Js) ≤ I ({psi }, {ξµsi j }) = ∑ i ripsi + ∑ i j psi ki jξ (ln ξ−1) = R[1−ξ−ξ ln ξ] (7.3) with the overall stationary activity R ≡ ∑ i ripsi = ∑ i j psi ki j (7.4) that gives the average number of transitions in the whole network per unit time. The Legendre transformation of this bound yields as bound on the generating function λ(z) ≥ R(ezJs/R − 1). (7.5) In particular, this bound holds for any current observable, such that the generating function (4.6) for the set of cycle currents can be written in a vectorial notation as λ(z) ≥ R(ez·Js/R − 1). (7.6) Due to the Gallavotti-Cohen symmetry (4.11) the activity-dependent bound can also be written as λ(z) ≥ R(e(−A−z)·Js/R − 1). (7.7) This bound is sharper than (7.6) for z · Js < −A · Js/2 = σ/2. Combining Eqs. (7.6) and (7.7) we obtain the activity-dependent bound λ(z) ≥ R [ e(|σ/2+z·J s |−σ/2)/R − 1 ] . (7.8) For an individual current Jα, the activity-dependent bound reads λα (z) ≥ R [ e( |σ/2+zJ s α |−σ/2)R − 1 ] . (7.9) The choice z = zA in Eq. (7.8) leads to λs (z) ≥ R [ e(|σ/2+zσ |−σ/2)/R − 1 ] (7.10) 108 7.1 Activity-dependent bound −10 0 10 20 −20 −15 −10 −5 0 5 0 10 20 30 40 50 −1 0 1 2 3 z λ(z) bound ξ h(ξ) bound Figure 7.1: Generating function (left) and rate function (right) for a five-state unicyclic network with rates ln k+i = (3, 3, 2, 3, 2) and ln k − i = (−1, 0,−1, 0, 0), affinityA = 15, current Js ≃ 2.25, and activity R ≃ 12.9. The functions are shown as solid lines and the activity-dependent bound (7.8) as dashed lines. Analytic continuations of the piecewise defined functions are shown as dotted curves. for the entropy change. In terms of the large deviation function for an individual current Jα, the application of the fluctuation refines Eq. (7.3) to hα (ξ) ≤ ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ R[1 + ξ − ξ ln |ξ |] − σξ, ξ ≤ −e−σ/(2R), R[1 − ξ + ξ ln |ξ |], ξ ≥ e−σ/(2R), R[1 − e−σ/(2R)] − σξ/2, otherwise, (7.11) where the intermediate, linear part stems from the construction of the convex envelope. An illustration of the activity-dependent bound (7.8) and (7.11), respectively, is pro- vided in Fig. 7.1. This bound is typically tighter than the dissipation-dependent bound for far from equilibrium conditions, i.e., for large affinity. For example, for a random walk on a unicyclic network with N sites, a uniform forward stepping rate k and a van- ishing backward stepping rate, which implies divergent affinity, the bound in Eq. (7.8) is saturated. Specifically, in this case the generating function is λ(z) = k ( ez/N − 1 ) , (7.12) the activity is R = k and the cycle current Js = k/N . For vanishing current at equilibrium, the activity-dependent bound reduces to the trivial statement λ(z) ≥ 0. 109 7 Activity-dependent and asymptotic bound Our numerics indicates that the affinity- and topology-dependent bound is always tighter than the activity-dependent bound in unicyclic networks. For multicyclic net- works the activity-dependent bound can be tighter. Furthermore, contrary to the affinity- and topology-dependent bound, the activity-dependent bound does not require knowl- edge of the topology of the network of states, only the average activity and, when using the variant (7.8), the average entropy production are required. Using (4.22) in the activity-dependent bound (7.8) for an individual current leads to F ≥ Js/R. (7.13) This relation provides a lower bound on the dispersion of an individual current, charac- terized by the Fano factor F, in terms of its average Js and the activity R. 7.2 Asymptotic bound for unicyclic networks The asymptotic bounds discussed in the following are exact results that become tighter than all previous bounds for large values of |z |. First we consider a unicyclic network with N states and affinity A. In this case, we can prove the following bound on the generating function: λ(z) ≥ JλARW(z,A, N ) + rARW − 1N N∑ i=1 ri, (7.14) where λARW(z,A, N ) is the generating function of the asymmetric random walk defined in Eq. (4.39), rARW ≡ k+ + k−, and k± ≡ *, N∏ i=1 k±i +- 1/N . (7.15) Our numerics indicate that with increasing |z | the difference between this bound and the actual generating function tends to zero. This fact is quite remarkable given the exponential growth of both functions. Unlike all other bounds presented so far, the bound (6.13) is not saturated at z = 0. Only for the case of uniform rates, i.e., k±i = k ±, the generating function (4.38) saturates the bound (7.14) globally. In order to prove this asymptotic bound, we consider the path weight (3.73) of a trajectory i(τ) in the unicyclic network, which is characterized by sojourn times Ti and the number of jumps m±i out of the state i in forward and backwards direction, respectively. The path weight then reads p[i(τ)] = psi(0) exp ⎡⎢⎢⎢⎢⎣− ∑ i (riTi + m+i ln k+i + m−i ln k−i ) ⎤⎥⎥⎥⎥⎦, (7.16) 110 7.2 Asymptotic bound for unicyclic networks where we use the same notation for the rates as in Eq. (6.13). The weight of the same paths in a process with modified transition rates kˆ±i and corresponding exit rates rˆi and stationary distribution pˆsi becomes pˆ[i(τ)] = pˆsi(0) exp ⎡⎢⎢⎢⎢⎣− ∑ i (rˆiTi + m+i ln kˆ+i + m−i ln kˆ−i ) ⎤⎥⎥⎥⎥⎦ . (7.17) For these modified transition rates we choose an asymmetric random walk with the uniform rates k±i = k ± from Eq. (7.15), which leads to rˆ ≡ k+ + k−. Ensemble averages using the path weight pˆ[i(τ)] are denoted as ⟨. . .⟩ARW. Since we are ultimately interested in the long-time limit, we are free to choose a fixed initial state i(0) = 1 for all trajectories, without loss of generality. The generating function for the integrated cycle current X (t) can be rewritten as λ(z) = lim t→∞ 1 t ln ⟨ ezX[i(τ)] ⟩ = lim t→∞ 1 t ln ∑ i(τ) ezX[i(τ)] p[i(τ)] pˆ[i(τ)] pˆ[i(τ)] = lim t→∞ 1 t ln ⟨ ezX[i(τ)] N∏ i=1 ( k+i k+ )m+i ( k−i k− )m−i e−(ri−rˆi )Ti ⟩ ARW , (7.18) where the sum in the first line denotes a path integral over all stochastic trajectories. The path dependent variables m±i count the jumps out of state i in forward or backward direction and Ti is the total sojourn time in state i. These variables are identically distributed in theARW-ensemble. The constant contribution from the initial probabilities ps1 and pˆ s 1 vanishes in the long-time limit. The probability pˆ(X ) is the probability that the fluctuating current is X in the ARW- ensemble. It is the sum of the weight of all trajectories for which the current is X . Using this pˆ(X ), Eq. (7.18) can be written as λ(z) = 1 t ln ∑ X pˆ(X )ezX ⟨ exp [ N∑ i=1 m+i ln ( k+i /k + ) + N∑ i=1 m−i ln ( k−i /k −) − N∑ i=1 (ri − rˆ)Ti ] X ⟩ ARW ≥ 1 t ln ∑ X pˆ(X )ezX exp [ N∑ i=1 ⟨ m+i |X ⟩ ARW ln ( k+i /k + ) + N∑ i=1 ⟨ m−i |X ⟩ ARW ln ( k−i /k −) − N∑ i=1 (ri − rˆ)t/N ] , (7.19) where the conditioned average in the first line represents a functional integration over all trajectories with fluctuating current equal to X and we used Jensen’s inequality from 111 7 Activity-dependent and asymptotic bound ③ ✷ ✶ ✲✁✂ ✲✄✂ ✲☎✂ ✲✆✂ ✂ ✆✂ ☎✂ ✄✂ ✁✂ ✲✁✂ ✲✄✂ ✲ ☎ ✂ ✲ ✆ ✂ ✂ ✆ ✂ ☎ ✂ ✄✂ ✁✂ ✂ ✂✵☎ ✂✵✁ ✂✵ ✝ ✂✵✞ ✆ Figure 7.1: Asymptotic bound for the house-shaped network with two fundamental cycles shown in Fig. 4.1b. The color code represents the ratio (7.24) between the generating function and the bound. Black dashed lines indicate the borders between sectors with constant relevant cycles Cˆ(z). For each sector, the relevant cycle Cˆ(z) is shown in white. The affinity of the three-cycle is A1 = 8 and the affinity of the four-cycle is A2 = 6. the first to the second line. Due to (7.15) the terms with the logarithms vanish, leading to the final result in Eq. (7.14). 7.3 Asymptotic bound for multicyclic networks In order to obtain an asymptotic bound that is also valid for multicyclic networks we define an arbitrary closed path C, which is a sequence of jumps that finishes at the state it started, as C ≡ [i(1) → i(2) → · · · → i(nC) → i(1)], (7.20) where nC is the length of the closed path. With this path we associate a geometric mean of the transition rates γC ≡ (ki(1),i(2)ki(2),i(3) . . . ki(nC ),i(1))1/nC (7.21) and integer winding numbers m βC that count how often the elementary cycle β is com- pleted within the path C. Applying a theorem valid for arbitrary non-negative matri- ces [21, Lemma 3.5.3] to the matrix Li j (z) + δi j maxℓ rℓ we obtain λ(z) + max ℓ rℓ ≥ f (z, C) ≡ γC exp *., 1 nC ∑ β m βC zβ +/- (7.22) 112 7.3 Asymptotic bound for multicyclic networks for any closed path C. The best bound on λ(z) in Eq. (7.22) is obtained by choosing an optimal path Cˆ(z), which in principle depends on z, that maximizes the r.h.s of Eq. (7.22) in the large z regime. First we consider this optimal path for the unicyclic network. In this case, the optimal path is a single cycle in the forward direction with mC = 1 if z > 0. If we consider a path C with two cycles, i.e., mC = 2, the bound remains the same as the number of states nC also doubles. If the closed path is not a direct cycle but contains, for example, one backward jump, then γC can become larger. However, such a backward jump also makes nC larger and hence, the exponent in Eq. (7.22) smaller. Since we are interested in the large z regime, this second effect should be dominant. Hence, for z > 0 the bound (7.22) leads to λ(z) ≥ k+ exp(z/N ) −max ℓ rℓ . (7.23) Even though this bound is different from (7.14), they both predict the same exponential growth, with the same prefactor, for large z > 0. The same reasoning is valid for z < 0, with the optimal path being a single cycle in the negative direction. For multicyclic networks we consider the house-shaped network with five states shown in Fig. 4.1b. This network consists of a cycle with three states and affinity A1 and a cycle with four states and affinity A2. We choose these cycles to be the fundamental cycles. This network also has a third cycle, which is the cycle with five states and affinity A1 +A2. Given a vector z = (z1, z2), the optimal path is the cycle that maximizes the r.h.s of Eq. (7.22). For large enough |z |, this optimal path depends only on the direction of the vector. Clearly a path that includes other cycles will lead to a weaker bound. A contour plot of the ratio f (z, Cˆ(z))/[λ(z) + maxℓ] for this house-shaped network is shown in Fig. 7.1. Remarkably, the r.h.s. of (7.22) captures the leading order of the asymptotics for large |z |, i.e., for large |z | f (z, Cˆ(z)) λ(z) + maxℓ rℓ → 1. (7.24) Only in the lines separating regions dominated by different cycles in Fig. 7.1 does this ratio tend to slightly lower values. Along this line the dominant cycle is degenerate. As we show below, relation (7.24) is valid for any multicyclic network. Hence, we conclude that the asymptotic bound predicts the exponential growth of the generating function, apart from exceptional regions in z where the optimal cycle is degenerate. In order to prove the relation (7.24), we write the characteristic polynomial (3.27) of the tilted Markov generator (4.7) as 0 = χ(z, λ(z)) = N∏ i=1 [−ri − λ(z)] + ∑ C (−1)CγnCC emC ·z ∏ j