# Capillary nanowaves on surfactant-laden liquid films with surface viscosity
and elasticity

Yixin Zhang<EMAIL_ADDRESS>Physics of Fluids Group, Max Planck Center
Twente for Complex Fluid Dynamics and J. M. Burgers Centre for Fluid Dynamics,
University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands  Zijing
Ding  School of Energy Science and Engineering, Harbin Institute of
Technology, Harbin, China

###### Abstract

Thermal motions of molecules can generate nanowaves on the free surface of a
liquid film. As nanofilms are susceptible to the contamination of surfactants,
this work investigates the effects of surfactants on dynamics of nanowaves on
a bounded film with a finite depth, using both molecular dynamics simulations
and analytical theories. In molecular simulations, a bead-spring model is
adopted to simulate surfactants, where beads are connected by the finite
extensive nonlinear elastic potentials. Fourier transforms of the film surface
profiles $h(x,t)$ extracted from molecular simulations are performed to obtain
the static spectrum $|h_{q}|_{\mathrm{rms}}$ and temporal correlations of
surface modes $<h_{q}(0)h_{q}^{*}(t)>$. It is shown that the spectral
amplitude is increased for the contaminated liquid surface compared to the
clean surface because surfactants can decrease surface tension. A higher
concentration of surfactants on the surface not only decreases the surface
tension but also causes elastic energy to the free surface, as the scaling of
spectral amplitude with wavenumbers changes from $|h_{q}|_{\mathrm{rms}}\sim
q^{-1}$ to $|h_{q}|_{\mathrm{rms}}\sim q^{-2}$ for modes with large
wavenumbers. Regarding the temporal correlations of surface modes, it is
observed that the presence of surfactants leads to a slower decay, which,
however, cannot be predicted by only considering the decreased surface
tension. Based on the Boussinesq–Scriven model for surface viscosity, a linear
stability analysis of Stokes flow for films with arbitrary depth is conducted
and the obtained dispersion relation considering surface viscosity can justify
the simulation results.

##  I introduction

Surfactants are widely used in a number of industrial applications such as
foams and emulsions [1], detergents, inks [2], and oil recovery [3]. The
presence of surfactants can significantly alter the behaviors of liquid
surfaces. For example, surfactants in water may lead to a Marangoni stress on
the bubble surface that can slow the rising of bubbles [4]. Thus, it is of
great interest and importance to study the effects of surfactants on various
kinds of interfacial fluid dynamics elaborately.

Thermal capillary waves (TCW) on a free liquid surface are waves spontaneously
excited by thermal motions of molecules. These waves have been used in
experiments to measure fluid properties such as surface tension [5] and
viscoelasticity [6]. These waves are also critical to the instability of
nanofilms [7, 8, 9, 10], coalescence of nanodroplets [11] and bubbles, and
rupture of nanojets [12]. The roughness on a liquid surface created by TCW is
usually on the scale of nanometres, but micrometer roughness can also be
obtained and thus observed optically using ultra-law surface tension mixtures
($\gamma\sim 10^{-6}$ N/m) [13], since thermal roughness is usually
proportional to thermal length $\sqrt{k_{B}T/\gamma}$ ($k_{B}$, $T$, and
$\gamma$ are the Boltzmann constant, temperature, and surface tension
respectively). It is worth mentioning that the amplitude of TCW will diverge
with the increasing system size if only considering surface tension on the
free surface. The introduction of gravity or a certain binding potential,
however, can remedy this problem [14].

At thermal equilibrium, the amplitude of TCW, namely the static spectrum, can
be described by the famous capillary wave theory [15, 16]. Recently, the
capillary wave theory has been extended using a Langevin equation to describe
the transient dynamics of non-equilibrium TCW and their approach to thermal
equilibrium [17]. When in thermal equilibrium, it is well known that the
temporal correlations of overdamped surface modes of TCW show an exponential
decay [18, 19], with a decay rate given by the dispersion relation of the
system, which provides the basis for measuring the system’s property using TCW
[5, 6]. This technique has the advantage of being noninvasive since no
external forces are applied. As the surface of a nanometric liquid film can
often be contaminated by surface-active agents such as surfactants [20, 21],
it is crucial to understand how surfactants modify the properties of TCW.

Perhaps the most known effect of surfactants is their ability to decrease
surface tension. Thus, adding surfactants to a liquid surface may enhance the
amplitude of thermal waves and roughness. Apart from decreasing surface
tension, surfactants may form a monolayer on the liquid surface and lead to
the elastic properties of the liquid surface [22, 23]. For example, the
stability of emulsions, found in mixtures of water, oil, and amphiphilic
surfactants, may be influenced by the bending rigidity ($\kappa$) of the
surfactant monolayers [24, 25]. The bending rigidity of emulsions is on the
order of $\kappa\sim k_{B}T$ and may overwhelm the effects of surface tension
[25]. It is thus interesting to see how the bending rigidity affects thermal
capillary waves. In fact, a recent experiment used the static spectra of TCW
in an oil-water-surfactant system to measure the bending rigidity of the
water-oil interface of droplets [26].

The adsorption of surfactants on a fluid surface may also increase that
surface’s resistance to deformations and lead to a new surface property called
surface viscosity [27]. The importance of surface viscosity on free surface
flows has been generally acknowledged now. For example, Joye et al. [28]
experimentally showed that high surface viscosity can prevent asymmetric
drainage of thin films in foams. Ponce-Torres et al. [29] found surface
viscosity is responsible for the accumulation of surfactants in the satellite
droplets formed by the breakup of a large surfactant-laden droplet. Surface
viscosity can also change the instability and pinch-off profile of viscous
threads [30, 31]. For thin liquid films destabilized by the disjoining
pressure, surface viscosity is able to decrease the growth of perturbations
[32, 33]. Shen et al. [34] theoretically showed surface viscosity contributes
an overall damping effect to the amplitude of the capillary waves on films
with infinite depth. However, there are very few studies about the effects of
surface viscosity on TCW. Given the importance of surface viscosity on free
surface flows, seeking accurate methods to measure surface viscosity is
necessary. Though much progress has been made, it is still debatable how to
obtain reliable values of surface viscosity, since Marangoni effects and
surface viscosity can often coexist during measurements [35]. Experimental
measurements of surface viscosity have reported values that are orders of
magnitude apart and are controversial [36].

As experiments at the nanoscale are difficult to perform due to the small
temporal and spatial scales, molecular dynamics simulations (MD) become a
valuable tool to investigate nanoflows and can act as virtual experiments to
validate newly developed theories. In terms of nanoflows involving free
surfaces, there are increasing molecular studies for nanodroplets, nanobubbles
[37], nanofilms [8] and nanojets [12]. MD studies of TCW have been carried out
for a variety of problems including TCW for free liquid films [38], TCW for
films on no-slip [39] and anisotropic-slip [19] substrates, to name a few. The
introduction of surfactants to TCW in MD has been implemented in several
studies [25, 40] but previous works only considered the effects of surfactants
on surface tension and bending rigidity of liquid surface by examining the
static spectrum. The dynamics of TCW such as the relaxation of TCW
correlations associated with surface viscosity and elasticity have not been
studied.

In this work, MD simulations are used to study the effects of surfactants on
the overdamped thermal capillary waves for liquid films bounded on substrates.
A bead-spring model is adopted to simulate surfactants, where beads are
connected by the finite extensive nonlinear elastic potentials. The surface
concentrations of surfactants are varied to see how they change the behaviors
of the liquid interface. We obtain surface modes of surface waves from MD
simulations and calculate their static spectra $|h_{q}|_{\mathrm{rms}}$ and
temporal correlations $<h_{q}(0)h_{q}^{*}(t)>$. The static spectrum is used to
infer bending rigidity and surface tension. The Boussinesq–Scriven model for
the surface viscosity is adopted and we perform a linear stability analysis of
a liquid film with arbitrary depth in Stokes limits and obtain the
corresponding dispersion relation for the first time, which is validated
against the MD temporal correlations.

This paper is organized as follows. In Sec. II, we formulate the problem that
we are going to solve and derive the dispersion relation for films with any
depth and surface viscosity. In Sec. III, the MD model of nanoscale liquid
films with surfactants is introduced. Sec. IV, shows the comparison between MD
and analytical theories. We conclude our findings and outline future
directions of research in Sec. V.

##  II Mathematical modeling

As shown in Fig. 1, we consider a Newtonian liquid film on a solid surface.
The free surface of the liquid film is covered with insoluble surfactants.
Initially, the film has a size $\left(L_{x},L_{y},h_{0}\right)$ and is quasi-
two dimensional (2D) by letting $L_{x}\gg L_{y}$. Due to thermal motions of
molecules, the film surface is fluctuating weakly around $h_{0}$ and the
instantaneous surface profile is given by $h(x,t)$.

Figure 1:  Sketch of a (quasi-two dimensional) surfactant-laden liquid film on
a plate. The $h_{0}$ is the initial film thickness and $h=h(x,t)$ is the film
thickness under spontaneous perturbations due to thermal fluctuations. The
film has a small depth $L_{y}$ in the $y$ direction (into the page).

###  II.1 Static spectrum and temporal correlations for surfactant-laden
liquid films

The static spectrum of capillary waves on a clean film can be determined by
the equipartition theorem; this forms the basis of classical capillary wave
theory [16]. For a film contaminated with surfactants, the static spectrum has
to be modified to consider the contributions of the elastic energy arising
from the formed monolayer of surfactants. In terms of a quasi-2D film under
small perturbations in Fig. 1, the extra free energy $f$ due to the change of
surface area is [25, 24]

$f\approx\frac{L_{y}}{2}\gamma\int{{{\left(\frac{\partial h}{\partial
x}\right)}^{2}}}dx+\frac{L_{y}}{2}\kappa\int{{{\left(\frac{\partial^{2}h}{\partial
x^{2}}\right)}^{2}}}dx.$ (1)

Note that the effect of disjoining pressure is ignored as we are going to
study a stable film where the effect of disjoining pressure is weak. As you
will see, in our MD simulations, the thickness of the film is much larger than
the cut-off distance so that the effect of disjoining pressure is safely
neglected. Let us define the Fourier transform of $\delta h=h(x)-h_{0}$ as
$h_{q}=\int{\delta h{{e}^{-iqx}}dx}$. From Parseval’s theorem, we can express
Eq. (1) in terms of Fourier modes:

$f=\frac{1}{2}\frac{{{L}_{y}}}{{{L}_{x}}}\gamma\sum{{{q}^{2}}{{\left|h_{q}\right|}^{2}}}+\frac{1}{2}\frac{{{L}_{y}}}{{{L}_{x}}}\kappa\sum{{{q}^{4}}{{\left|h_{q}\right|}^{2}}}.$
(2)

As each summand appears quadratically, it has the same energy
$\frac{1}{2}k_{B}T$, from the equipartition theorem, so that

$\frac{1}{2}{{k}_{B}}T=\frac{{{L}_{y}}}{{{L}_{x}}}\left(\frac{1}{2}\gamma{{q}^{2}}+\frac{1}{2}\kappa{{q}^{4}}\right){{\left|h_{q}\right|}^{2}}.$
(3)

Thus, the static spectrum for a surfactant-laden film is derived:

$S_{s}=\sqrt{\left\langle\left|h_{q}\right|^{2}\right\rangle}=\sqrt{\frac{{{L}_{x}}}{{{L}_{y}}}\frac{{{k}_{B}}T}{\gamma{{q}^{2}}+\kappa
q^{4}}}.$ (4)

Here $\left\langle...\right\rangle$ represents an ensemble average.

In thermal equilibrium, the temporal correlations of overdamped capillary
waves usually decay exponentially to zero [5, 39]

$C_{h_{q}h_{q}^{*}}=\frac{\left\langle
h_{q}\left(q,t\right){{h_{q}}^{*}}\left(q,t^{\prime}\right)\right\rangle}{S_{s}^{2}}={{e}^{\Omega\left(q\right)\left|t-t^{\prime}\right|}},$
(5)

where the asterisk denotes conjugate values. The decay rate $\Omega$
($\Omega<0$) is given by the dispersion relation of the system (temporal
growth rate of a surface mode), which is derived by a linear stability
analysis (LSA) in the next subsection.

###  II.2 LSA of Stokes flow for films with arbitrary depth and surface
viscosity

The liquid is assumed to be incompressible

$\nabla\cdot\mathbf{u}=0,$ (6)

where $\mathbf{u}=(u_{x},u_{z})$, and $u_{x},u_{z}$ are the velocities in $x$
and $z$ directions, respectively. The momentum equation with the assumption of
Stokes flow (the Reynolds number is small), is as follows

$\mu\nabla^{2}\mathbf{u}=\nabla p,$ (7)

where $\mu$ and $p$ are the viscosity and pressure of the liquid. Note that
transient inertia may come into play depending on the frequency of the
perturbations, which leads to a critical wavenumber below which the
correlations are underdamped instead of being overdamped [38]. However, as one
can see from our following molecular simulations, all measured correlations
are overdamped.

For the dynamic boundary condition at the liquid-vapour interface $z=h(x,t)$,
we have the leading-order Boussinesq–Scriven model for one-dimensional surface
[34]

$\mathbf{n}\cdot\bm{\tau}=\nabla_{s}\tilde{\gamma}-\tilde{\gamma}\left(\nabla_{s}\cdot\mathbf{n}\right)\mathbf{n}+\mathbf{f},$
(8)

where $\bm{\tau}$ is the hydrodynamic stress tensor,
$\tau_{ij}=-p\delta_{ij}+\mu\left(\partial u_{i}/\partial x_{j}+\partial
u_{j}/\partial x_{i}\right)$. The $\tilde{\gamma}$ is the surface tension
modified by the surface viscosity,
$\tilde{\gamma}=\gamma+\left(\mu_{d}+\mu_{s}\right)\nabla_{s}\cdot\mathbf{u}$.
$\nabla_{s}$ is the surface gradient. The dilatational ($\mu_{d}$) and shear
($\mu_{s}$) components of surface viscosity occur in additive pairs in a 2D
system and will hereafter be written in terms of a single parameter,
$\eta=\mu_{d}+\mu_{s}$ in the rest of this paper [32, 33]. The
$\mathbf{f}=\kappa\frac{\partial^{4}h}{\partial x^{4}}\mathbf{n}$ is the
elastic pressure [41]. Finally, $\mathbf{n}$ is the outward normal to the free
surface:

$\mathbf{n}=\frac{\left(-\partial h/\partial
x,1\right)}{\sqrt{1+{{\left(\partial h/\partial x\right)}^{2}}}}.$ (9)

Under the assumption of small perturbations (${\partial h}/{\partial x}\ll 1$)
the dynamic boundary condition is reduced to (in the normal direction):

$-p+2\mu\frac{\partial u_{z}}{\partial z}=\gamma\frac{\partial^{2}h}{\partial
x^{2}}+\kappa\frac{\partial^{4}h}{\partial x^{4}},$ (10)

and in the tangential directions to the surface

$\mu\left(\frac{\partial u_{x}}{\partial z}+\frac{\partial u_{z}}{\partial
x}\right)=\eta\frac{\partial^{2}u_{x}}{\partial x^{2}},$ (11)

where we have assumed surface tension (concentration of surfactants) is
uniform along the surface so that the ${\partial\gamma}/{\partial x}$ term is
negligible. This is because, at the nanoscale, surface diffusion can counter
advective effects and result in a uniform surfactant distribution [42]. This
is also directly found in our MD simulations shown in the appendix.

The kinematic condition at the free surface is given by

$u_{z}=\frac{\partial h}{\partial t}+u_{x}\frac{\partial h}{\partial x}.$ (12)

At $z=0$, the no penetration condition and Navier’s slip condition are,

$\displaystyle u_{z}$ $\displaystyle=0,$ (13) $\displaystyle u_{x}$
$\displaystyle=\ell\frac{\partial u_{x}}{\partial z}.$ (14)

Equations (6-14) are linearised using normal modes

$\left(u_{x},u_{z},p,h-h_{0}\right)=\left(\tilde{u}_{x},\tilde{u}_{z},\tilde{p},\tilde{h}\right){{e}^{\Omega
t+i{{q}}x}}.$ (15)

With those, the expression for the single variable $\tilde{u}_{z}$ can be
obtained from Eq. (6) and Eq. (7) as:

$\frac{{{d}^{4}}\tilde{u}_{z}}{d{{z}^{4}}}-2{{q}^{2}}\frac{{{d}^{2}}\tilde{u}_{z}}{d{{z}^{2}}}+{{q}^{4}}\tilde{u}_{z}=0.$
(16)

The general solution for $\tilde{u}_{z}$ is thus [18]

$\tilde{u}_{z}={{C}_{1}}\cosh\left(qz\right)+{{C}_{2}}\sinh\left(qz\right)+{{C}_{3}}qz\cosh\left(qz\right)+{{C}_{4}}qz\sinh\left(qz\right),$
(17)

and one can also obtain the solution for $\tilde{u}_{x}$, and $\tilde{p}$:

$\displaystyle\tilde{u}_{x}=\frac{i}{q}\frac{d\tilde{u_{z}}}{dz},$ (18)
$\displaystyle\tilde{p}=2\mu
q\left[{{C}_{3}}\cosh\left(qz\right)+{{C}_{4}}\sinh\left(qz\right)\right].$
(19)

Here, $C_{1}-C_{4}$ are four coefficients to be determined by the boundary
conditions Eqs. (10-14), whose linearised forms are

$-\tilde{p}+2\mu\frac{d\tilde{u}_{z}}{dz}=-\gamma q^{2}\tilde{h}-\kappa
q^{4}\tilde{h},\quad\mbox{at\ }\quad z=h$ (20a)
$\mu\left(\frac{d\tilde{u}_{x}}{dz}+iq_{x}\tilde{u}_{z}\right)=-\eta
q^{2}\tilde{u}_{x},\quad\mbox{at\ }\quad z=h$ (20b)
$\Omega=\frac{\tilde{u}_{z}}{\tilde{h}},\quad\mbox{at\ }\quad z=h$ (20c)
$\tilde{u}_{z}=0,\quad\mbox{at\ }\quad z=0$ (20d)
$\tilde{u}_{x}=\ell\frac{d\tilde{u}_{x}}{dz}.\quad\mbox{at\ }\quad z=0$ (20e)

Substituting Eqs. (17-19) into Eqs. (20) leads to the required dispersion
relation

$\displaystyle\Omega=-\frac{\gamma{{q}^{2}}+\kappa{{q}^{4}}}{\mu
q}\frac{B_{1}}{B_{2}},$ $\displaystyle\begin{aligned}
B_{1}=&\mu\left[\sinh\left(2qh_{0}\right)-2qh_{0}+4ql{{\sinh}^{2}}\left(qh_{0}\right)\right]\\\
&+\left\\{{{\sinh}^{2}}\left(qh_{0}\right)-{{q}^{2}}{h_{0}^{2}}+ql\left[\sinh\left(2qh_{0}\right)-2qh_{0}\right]\right\\}{{\eta}}q,\\\
B_{2}=&\mu\left\\{4{{q}^{2}}{h_{0}^{2}}+4{{\cosh}^{2}}\left(qh_{0}\right)+4ql\left[2qh_{0}+\sinh\left(2qh_{0}\right)\right]\right\\}\\\
&+\left[2qh_{0}+\sinh\left(2qh_{0}\right)+4ql{{\cosh}^{2}}\left(qh_{0}\right)\right]{{\eta}}q.\end{aligned}$

##  III Molecular dynamics simulations

Molecular dynamics simulations are used to simulate surfactant-laden liquid
films on substrates. These simulations are performed using the open-source
code LAMMPS [43]. The molecular system contains fluid atoms (liquid and its
vapor), surfactant molecules, and solid atoms, as shown in Fig. 2(a). A
surfactant molecule $H_{m}T_{n}$ consists of $m+n$ atoms and is amphiphilic
with $m$ hydrophilic atoms (the head group) and $n$ hydrophobic atoms (the
tail group).

Figure 2:  Snapshots of MD simulations of a surfactant-laden liquid film on a
substrate. (a) Initial setting. The fluid atoms are colored in orange and the
solid atoms are navy blue. The surfactant is made of chained beads with a
hydrophobic tail (in bronze) and a hydrophilic head (in red). (b) Equilibrium
configurations of a surfactant-laden interface.

The non-bonded intermolecular potentials $U$ between $i$-type atoms and
$j$-type atoms are simulated with the standard Lennard-Jones (LJ) 12-6
potential:

$U({{r}_{ij}})=\begin{cases}4\varepsilon_{ij}\left[{{\left(\frac{\sigma_{ij}}{{{r}_{ij}}}\right)}^{12}}-{{\left(\frac{\sigma_{ij}}{{{r}_{ij}}}\right)}^{6}}\right]&\text{if}\,\,\,{{r}_{ij}}\leq{{r}_{c,ij}},\\\
0&\text{if}\,\,\,{{r}_{ij}}>{{r}_{c,ij}},\\\ \end{cases}$ (21)

where $r_{ij},\varepsilon_{ij},\sigma_{ij}$ and ${r}_{c,ij}$ are the pairwise
distance, energy parameter, length parameter, and cut-off distance
respectively. The complete lists of parameters among the liquid (l), solid
(s), head group of surfactants (h), and tail group of surfactants (t) are
given in Table 1.

Atom type | Atom type | $\varepsilon_{ij}/\varepsilon_{ll}$ | $\sigma_{ij}/\sigma_{ll}$ | $r_{c,ij}/\sigma_{ll}$
---|---|---|---|---
L | L | 1 | 1 | 5.5
L | S | 0.65 | 0.8 | 5.5
L | H | 0.80 | 1 | 5.5
L | T | 1 | 1 | $2^{1/6}$
H | H | 0.5 | 1 | 5.5
H | T | 1 | 1 | $2^{1/6}$
H | S | 1 | 1 | $2^{1/6}$
T | T | 0.35 | 1 | $2^{1/6}$
T | S | 1 | 1 | $2^{1/6}$
Table 1: Interaction parameters among the liquid (L), solid (S), head group of
surfactants (H), and tail group of surfactants (T).

The liquid is simulated as argon whose $\varepsilon_{ll}$, $\sigma_{ll}$, and
atomic mass are $1.67\times{10^{-21}}$ J, 0.34 nm, and $6.63\times{10^{-26}}$
kg, respectively. The temperature of this system is kept at $T=85$ K or
$T^{*}=0.7\varepsilon_{ll}/k_{B}$ (* henceforth denotes LJ units and $k_{B}$
is the Boltzmann constant) using the Nosé-Hoover thermostat. At this
temperature, the mass density of liquid argon is $1.40\times{10^{3}}$ kg/m3
and number density $n_{l}^{*}=0.83/\sigma_{ll}^{3}$. The cut-off distance for
liquid-liquid interactions, beyond which the intermolecular interactions are
omitted, is chosen as ${r_{c,ll}}^{*}=5.5\sigma_{ll}$. Under above parameters
and conditions, the surface tension of the clean liquid is
$\gamma_{0}=1.52\times{10^{-2}}$ N/m obtained by the mechanical route [44].
The dynamic viscosity of liquid is $\mu=2.87\times{10^{-4}}$ kg/(ms)
calculated by the Green-Kubo relation.

The substrate is platinum with a face-centred cubic (fcc) structure, and we
use its isotropic $\left\langle 100\right\rangle$ surface to contact the fluid
[19]. The platinum mass density is $21.45\times{10^{3}}$ kg/m3 with an atomic
mass of $3.24\times{10^{-25}}$ kg. The solid substrate is rigid in MD
simulations. The liquid-solid interactions are modelled by the same 12-6 LJ
potential with $\varepsilon_{ls}=C\varepsilon_{ll}$ and
$\sigma_{ls}=0.8\sigma_{ll}$. One may vary $C$ to obtain different amounts of
slip. Here we choose $C=0.65$ so that there is nearly no slip between liquid
and solid, since the effects of slip on capillary wave dynamics have been
examined elaborately in our previous works [17, 19].

To simulate surfactants, a coarse-grained model called the ‘bead-spring’ model
is adopted [45, 46]. The surfactant molecule $H_{m}T_{n}$ consists of $m+n$
atoms connected by the finite extensible nonlinear elastic (FENE) potential
[45]

$U=-0.5KR_{0}^{2}\mathrm{ln}\left[1-\left(\frac{r}{R_{0}}\right)^{2}\right].$
(22)

Here $K^{*}=12\varepsilon/\sigma^{2}$ and $R_{0}^{*}=1.4\sigma$ are used. The
tail group of a surfactant molecule is hydrophobic and this is achieved by
setting the $r^{*}_{c,lt}=2^{1/6}\sigma_{ll}$ so that interactions between
liquid and the tail group are purely repulsive. In our simulations, the type
of surfactant $H_{4}T_{4}$ is used.

The initial dimensions of the liquid film $(L_{x},L_{y},h_{0})$ in Fig. 1(a)
are chosen as $L_{x}=31.4$ nm, $L_{y}=3.14$ nm, and $h_{0}=3.14$ nm. The solid
substrate has the same lateral size as that of the liquid film and has a
thickness $h_{s}=0.78$ nm. As $L_{y}\ll L_{x}$, the system is quasi-2D.
Periodic boundary conditions are applied in both the $x$ and $y$ directions of
the system whilst vapour particles are reflected specularly in the $z$
direction at the top boundary of the system.

In one single simulation of the surfactant-laden film, the simulation with the
initial setting shown in Fig. 2(a) is run for 20 ns with a timestep 8.57 fs to
reach its equilibrium state, see Fig. 2(b). After the surface has reached
equilibrium, the simulation is run to output the positions of atoms every 2000
steps. The free surface position is defined as the usual equimolar surface and
the way to extract the surface profile $h(x,t)$ from positions of atoms in MD
simulations is detailed in our previous work [8]. After obtaining $h(x,t)$,
Fourier transforms are performed to obtain the amplitude of surface modes
$h(q,t)$.

##  IV Results and Discussions

In this section, we present and discuss the results of MD simulations and
their comparison to analytical solutions.

###  IV.1 Enhanced static spectrum, surface roughness and elasticity

The symbols in Fig. 3 represent the static spectra $|h_{q}|_{rms}$ of surface
waves obtained from simulations, which are the root mean square (rms) of
surface modes $h_{q}(q,t)$ (averaged over 20000 times).

Figure 3:  Static spectrum for capillary waves laden with different numbers
($N$) of surfactants.

One can see from MD results that the wave amplitude is enhanced by increasing
the number of surfactants on the liquid surface, since surfactants can
decrease surface tension. The values of surface tension for the film surface
contaminated with different numbers (concentrations) of surfactants can be
calculated from independent MD simulations through the mechanic route [44] and
are given in Fig. 4 (see the black squares). The linear relationship between
surface tension and surfactant concentrations is found only be valid for low
surfactant concentrations while surface tension decreases more rapidly for
high surfactant concentrations.

Figure 4:  Effects of the number of surfactants on surface tension and bending
rigidity. The values of surface tension are obtained by independent MD
simulations while the values of bending rigidity are obtained by fitting the
static spectra with predictions.

In Fig. 3, for a clean surface ($N=0$, green squares) and a surface with a low
concentration of surfactants (e.g., $N=54$, blue triangles), their spectra
show the conventional scaling $|h_{q}|_{rms}\sim q^{-1}$ for all wavenumbers,
which is the scaling when only considering surface tension on the surface.
Using the independently measured values of surface tension as shown in Fig. 4,
the MD spectra of $N=0$ and $N=54$ agree well with the theory, i.e., Eq. (4)
without the elastic term. However, when a surface is deposited with a high
concentration of surfactants (e.g., $N=250$, blue triangles), the scaling
changes from $-1$ to $-2$ for waves with large wavenumbers. The new scaling
can be justified by the elasticity of the surfactant monolayer formed on the
liquid surface, which gives an energy contribution $\sim\kappa
q^{4}|h_{q}|^{2}$ as described by Eq. (4). By fitting the MD results with Eq.
(4), the value of bending rigidity is obtained and its relation with
surfactant concentrations is given in Fig. 4, which shows that the bending
rigidity increases with the number of surfactants. We note that the
measurements of TCW spectra at large wavenumbers from MD simulations can
depend on how to define the free surface [38], which makes it more difficult
to measure the precise value of bending rigidity if it is weak.

The competition between elasticity and surface tension defines a length scale
${\lambda_{\kappa}}=2\pi\sqrt{\kappa/\gamma}$ and elasticity is dominant at
lengths smaller than $\lambda_{\kappa}$. This can be seen from the case of
$N=250$ in Fig. 3 where the scaling $-2$ is dominant at large wavenumbers.
Using the prediction for the strength of bending rigidity $\kappa\approx
k_{B}T$ based on the thermodynamic arguments [24], one finds
$\lambda_{\kappa}=2\pi\sqrt{k_{B}T/\gamma}$, which is on the order of thermal
length and indicates that elasticity plays more important roles at the
nanoscale.

Figure 5:  The relation of surface roughness with (a) number of surfactants
and (b) thermal length.

The rms roughness of a liquid surface is calculated as (in discrete forms):

$W=\sqrt{\left\langle\frac{1}{{{N}_{b}}}\sum\limits_{i=1}^{i={{N}_{b}}}{{{\left({{h}_{i}}-\bar{h}\right)}^{2}}}\right\rangle}=\sqrt{\frac{1}{L_{x}^{2}}\sum\limits_{i=1}^{i={{N}_{b}}}{{{\left|{{h}_{q}}\right|}_{rms}^{2}}}},$
(23)

where $N_{b}$ is the number of bins used to extract the surface profile from
MD simulations. Note that this sum of modes can be usually integrated by
taking $q=2\pi/L_{x}$ to approach infinitely small values [17, 25]. However,
here our film simulated is short making the integration meaningless.

The obtained $W$ for different numbers of surfactants is shown in Fig. 5(a)
and it is found surface roughness is enhanced by adding surfactants. In Fig.
5(b), for $N<176$, the surface roughness is linearly proportional to thermal
length $\sqrt{k_{B}T/\gamma}$ as predicted [17]. However, due to the strong
elasticity of the liquid surface for $N>176$, this linear relationship breaks
down and elasticity has the potential to reduce surface roughness, in
competition with the decreased surface tension.

###  IV.2 Surface viscosity of surfactant-laden liquid surface

Figure 6:  Decay of temporal correlations of surface modes at two different
wavenumbers ($q=2\pi/L_{x}$ and $q=4\pi/L_{x}$) for case (a) $N=54$ and (b)
$N=176$.

Fig. 6 shows the temporal correlations $C_{h_{q}h_{q}^{*}}$ of surface modes
at two different wavenumbers ($q=2\pi/L_{x}$ and $q=4\pi/L_{x}$) for the case
$N=54$ (see Fig. 6(a)) and the case $N=176$ (see Fig. 6(b)). As expected, the
temporal correlations decay to zero with time and temporal correlations of
waves with a larger wavenumber decay much faster [19].

For the case $N=54$, its MD results (represented by symbols) can be predicted
nicely using Eq. (5) with $\gamma=0.9\gamma_{0}$ and $\kappa=0$ independently
obtained above (see Fig. 4), which suggests vanishing effects of surface
viscosity. However, when the number of surfactants is increased to the case
$N=176$, the relaxation of correlations is significantly slowed evidenced by
comparing black squares in Fig. 6(b) to the black squares in Fig. 6(a). Such a
slow relaxation can not be explained only considering the decrease of surface
tension, since the predictions (see the dash lines in Fig. 6(b)) using
$\gamma=0.57\gamma_{0}$ and $\kappa=0.4k_{B}T$ (obtained above) still
underestimate the MD results a lot.

We attribute this inconsistency to surface viscosity and the MD results can be
fitted by Eq. (5) and the newly derived dispersion relation Eq. (21) with a
value of surface viscosity $\eta=2.25\times 10^{-8}$ kg/s. The inferred
surface viscosity is a reasonable value compared to experimentally measured
values of sodium dodecyl sulphate (SDS) solution [36], which ranges from
$10^{-8}$ to $10^{-6}$ kg/s in different measurement techniques. By fitting
temporal correlations with Eq. (5) and Eq. (21) for other surfactant
concentrations, we obtain the relation between surface viscosity and the
number of surfactants shown in Fig. 7(a), which shows the same trend as
experimental results [33, 47]. Fig. 7(b) shows the reduction of the dispersion
relation (for $q=2\pi/L_{x}$) by increasing the number of surfactants at the
free surface. The reduction is due to both the decreased surface tension (blue
line) and the increased viscosity(black line).

Figure 7: (a) The effects of the varying number of surfactants on the values
of surface viscosity measured in MD. (b) The effects of the varying number of
surfactants on the dispersion relation of the wave with $q=2\pi/L_{x}$.

The ability of surface viscosity to slow relaxations of temporal correlations
has its maximum, as plotted in Fig. 8(a), where a further increase of surface
viscosity beyond $\eta=1\times 10^{-6}$ kg/s leads to very small changes of
the correlations.

Figure 8:  (a) Damping effects of surface viscosity on the relaxation of
correlations. The other parameters used are $\gamma=0.57\gamma_{0}$ and
$\kappa=0.4k_{B}T$. (b) Effects of surface viscosity on the dispersion
relation with $\gamma=\gamma_{0}$ and $\kappa=0$. The inset shows the
dispersion relation with a different number of surfactants.

Simplifications of the newly derived dispersion relation can be done for two
asymptotic limits $qh_{0}\ll 1$ and $qh_{0}\gg 1$ (Notably $qh_{0}=0.63$ in
our work for $q=2\pi/L_{x}$) as shown in Fig. 8(b). In the limits of
$qh_{0}\ll 1$ and no slip, the dispersion relation is simplified to

$\Omega=-\frac{\left(\gamma{{q}^{2}}+\kappa{{q}^{4}}\right){{q}^{2}}h_{0}^{3}}{12\mu}\frac{4\mu+{{h}_{0}}\eta{{q}^{2}}}{\mu+{{h}_{0}}\eta{{q}^{2}}},$
(24)

which is the dispersion relation for thin liquid films obtained earlier by
Edwards and Oron [32] though they did not have the elastic term in Eq. (25).
In the limits of zero surface viscosity, the constant factor in Eq. (25) is
$1/3$ (suggesting a no-shear boundary condition at the free surface) while it
is $1/12$ (suggesting a tangentially immobile condition at the free surface)
in the limits of infinite surface viscosity. Thus, surface viscosity has at
most a fourfold ability to decrease the relaxation of waves with $qh_{0}\ll
1$, as demonstrated by the dispersion relation in Fig. 8(b). In our
simulations, due to the reduced surface tension, the decrease of dispersion
relation by increasing the number of surfactants is larger than four times
(see Fig. 7(b) and the inset of Fig. 8(b) for the case $N=250$).

However, in the opposite limits $qh_{0}\gg 1$, Eq. (II.2) is simplified to

$\Omega=-\frac{\gamma+\kappa q^{2}}{2\mu}q,$ (25)

which is independent of surface viscosity. This is can be seen from Fig. 8(b)
where all solid lines with different surface viscosity overlap each other for
wavenumbers $qh_{0}\gg 1$. Notably, this counter-intuitive result is due to
the assumption of negligible inertia terms in momentum equations, which is
usually valid at the nanoscale. For capillary waves on macroscopic films with
infinite depth ($qh_{0}\gg 1$), Shen et al. [34] considered both the effects
of surface viscosity and inertia and derived a dispersion relation depending
on surface viscosity. Equation (26) in our work can be recovered by neglecting
the inertia terms of the dispersion relation in Ref. [34].

Finally, we remark that the bending rigidity has the potential to speed up the
relaxation of correlations but in our simulations its values are too weak to
have noticeable effects on the wavenumber $q=2\pi/L_{x}$ (see the red line
with $\kappa=0$ and the blue line $\kappa=1.4k_{B}T$ in Fig. 9). However much
larger bending rigidity can be obtained in experiments [48] so that the
effects of elasticity can be significant (e.g. the black line with $100k_{B}T$
in Fig. 9 ).

Figure 9:  The effects of elasticity on the relaxation of temporal
correlations. The parameters used are $\gamma=0.57\gamma_{0}$ and $\eta=0$.

##  V Conclusions

In this work, the effects of surfactants on the statics and dynamics of
thermal capillary waves are studied both theoretically and numerically.
Molecular simulations are used to simulate surfactant-laden free surface
flows. The obtained static spectra show surfactants can enhance the elasticity
of the liquid surface apart from reducing surface tension. The slower decay of
temporal correlations of surface modes is observed in simulations by adding
surfactants and can be predicted by a newly derived dispersion relation
considering surface viscosity. We believe the new theories developed here and
molecular simulations performed can motivate new experiments using thermal
capillary waves to measure surface viscosity.

It is worth mentioning that our current molecular system is ideal, in the
future, it is interesting to use more realistic molecular models or map our LJ
system to water or other simple globular solvents in colloidal science. We
also remark that the choice of the model for surfactants or parameters in the
FENE model adopted here is not unique. For example, Laradji and Mouritsen [25]
use a different model for the surfactant. The original FENE model [45] use the
same cutoff for all beads unlike us. In the future, a systematic test of
models and their parameters is needed to see how they change the elasticity
and surface viscosity.

Thermal fluctuations of fluid surface have been shown to increase the
instability of nanofilms and nanojets [8, 12]. As surfactants can decrease
surface tension and increase thermal roughness, the presence of surfactants
may enhance the instability and rupture of nanofilms. For nanobubbles, one may
intuitively think surfactants can stabilize surface nanobubbles by reducing
surface tension and the pressure-driven diffusion of gas from bubbles.
However, experiments have shown the opposite effect of surfactants [20, 49]
and this may be due to the instability of the bubble surface caused by
extremely small surface tension and large surface fluctuations.

##  VI Acknowledgments

Zhang wishes to thank for the discussions with Detlef Lohse, Yibo Chen and
Hongguang Zhang. Ding wishes to thank the financial support from National
Natural Science Foundation of China (Grant No. 12102109).

## appendix: Distribution of surfactants

Figure A1: (a) Distribution of the number of surfactants along the surface.
(b) Change of the number of surfactants with time in the middle region of the
surface.

Fig. A1(a) shows the distribution of surfactants along the surface at three
consecutive time while Fig. A1(b) shows the change of the number of
surfactants in the middle region of the surface. Though the number of
surfactants fluctuates with time due to the nature of molecules, it can be
seen that the distribution of surfactants is uniform along the surface.

## References

  * Sadoc and Rivier [2013] J.-F. Sadoc and N. Rivier,  _Foams and emulsions_ , Vol. 354 (Springer Science & Business Media, 2013).
  * Deng  _et al._ [2018] Y. Deng, X. Zheng, Y. Bai, Q. Wang, J. Zhao, and J. Huang, Surfactant-controlled ink drying enables high-speed deposition of perovskite films for efficient photovoltaic modules, Nat. Energy 3, 560 (2018).
  * Negin  _et al._ [2017] C. Negin, S. Ali, and Q. Xie, Most common surfactants employed in chemical enhanced oil recovery, Petroleum 3, 197 (2017).
  * Takagi and Matsumoto [2011] S. Takagi and Y. Matsumoto, Surfactant effects on bubble motion and bubbly flows, Annu. Rev. Fluid Mech. 43, 615 (2011).
  * Kim  _et al._ [2003] H. Kim, A. Rühm, L. B. Lurio, J. K. Basu, J. Lal, D. Lumma, S. G. J. Mochrie, and S. K. Sinha, Surface dynamics of polymer films, Phys. Rev. Lett. 90, 068302 (2003).
  * Jiang  _et al._ [2007] Z. Jiang, H. Kim, X. Jiao, H. Lee, Y.-J. Lee, Y. Byun, S. Song, D. Eom, C. Li, M. Rafailovich, _et al._ , Evidence for viscoelastic effects in surface capillary waves of molten polymer films, Phys. Rev. Lett. 98, 227801 (2007).
  * Grün  _et al._ [2006] G. Grün, K. Mecke, and M. Rauscher, Thin-film flow influenced by thermal noise, J. Stat. Phys. 122, 1261 (2006).
  * Zhang  _et al._ [2019] Y. Zhang, J. E. Sprittles, and D. A. Lockerby, Molecular simulation of thin liquid films: Thermal fluctuations and instability, Phys. Rev. E 100, 023108 (2019).
  * Vrij [1966] A. Vrij, Possible mechanism for the spontaneous rupture of thin, free liquid films, Discuss. Faraday Soc. 42, 23 (1966).
  * Vrij and Overbeek [1968] A. Vrij and J. T. G. Overbeek, Rupture of thin liquid films due to spontaneous fluctuations in thickness, J. Am. Chem. Soc. 90, 3074 (1968).
  * Perumanath  _et al._ [2019] S. Perumanath, M. K. Borg, M. V. Chubynsky, J. E. Sprittles, and J. M. Reese, Droplet coalescence is initiated by thermal motion, Phys. Rev. Lett. 122, 104501 (2019).
  * Zhao  _et al._ [2019] C. Zhao, J. E. Sprittles, and D. A. Lockerby, Revisiting the rayleigh–plateau instability for the nanoscale, J. Fluid Mech. 861 (2019).
  * Aarts  _et al._ [2004] D. G. Aarts, M. Schmidt, and H. N. Lekkerkerker, Direct visual observation of thermal capillary waves, Science 304, 847 (2004).
  * Rowlinson and Widom [2013] J. S. Rowlinson and B. Widom,  _Molecular theory of capillarity_ (Courier Corporation, 2013).
  * Buff  _et al._ [1965] F. Buff, R. Lovett, and F. Stillinger Jr, Interfacial density profile for fluids in the critical region, Phys. Rev. Lett. 15, 621 (1965).
  * MacDowell [2017] L. G. MacDowell, Capillary wave theory of adsorbed liquid films and the structure of the liquid-vapor interface, Phy. Rev. E 96, 022801 (2017).
  * Zhang  _et al._ [2021a] Y. Zhang, J. E. Sprittles, and D. A. Lockerby, Thermal capillary wave growth and surface roughening of nanoscale liquid films, J. Fluid Mech. 915 (2021a).
  * Henle and Levine [2007] M. L. Henle and A. J. Levine, Capillary wave dynamics on supported viscoelastic films: single and double layers, Phys. Rev. E 75, 021604 (2007).
  * Zhang  _et al._ [2021b] Y. Zhang, D. A. Lockerby, and J. E. Sprittles, Relaxation of thermal capillary waves for nanoscale liquid films on anisotropic-slip substrates, Langmuir 37, 8667 (2021b).
  * Lohse  _et al._ [2015] D. Lohse, X. Zhang, _et al._ , Surface nanobubbles and nanodroplets, Rev. Mod. Phys. 87, 981 (2015).
  * Maali  _et al._ [2017] A. Maali, R. Boisgard, H. Chraibi, Z. Zhang, H. Kellay, and A. Würger, Viscoelastic drag forces and crossover from no-slip to slip boundary conditions for flow near air-water interfaces, Phys. Rev. Lett. 118, 084501 (2017).
  * Langevin [2014] D. Langevin, Rheology of adsorbed surfactant monolayers at fluid surfaces, Annu. Rev. Fluid Mech. 46, 47 (2014).
  * Jiang and Ding [2022] W. Jiang and Z. Ding, Thin liquid films down a vertical microfiber: Effect of curvature elasticity, Phys. Rev. E 105, 035104 (2022).
  * Helfrich [1978] W. Helfrich, Steric interaction of fluid membranes in multilayer systems, Z. Naturforsch. A 33, 305 (1978).
  * Laradji and Mouritsen [2000] M. Laradji and O. G. Mouritsen, Elastic properties of surfactant monolayers at liquid–liquid interfaces: A molecular dynamics study, J. Chem. Phys 112, 8621 (2000).
  * Bolognesi  _et al._ [2016] G. Bolognesi, Y. Saito, A. I. Tyler, A. D. Ward, C. D. Bain, and O. Ces, Mechanical characterization of ultralow interfacial tension oil-in-water droplets by thermal capillary wave analysis in a microfluidic device, Langmuir 32, 3580 (2016).
  * Scriven and Sternling [1960] L. Scriven and C. Sternling, The marangoni effects, Nature 187, 186 (1960).
  * Joye  _et al._ [1994] J.-L. Joye, G. J. Hirasaki, and C. A. Miller, Asymmetric drainage in foam films, Langmuir 10, 3174 (1994).
  * Ponce-Torres  _et al._ [2017] A. Ponce-Torres, J. Montanero, M. Herrada, E. Vega, and J. Vega, Influence of the surface viscosity on the breakup of a surfactant-laden drop, Phys. Rev. Lett. 118, 024501 (2017).
  * Martínez-Calvo and Sevilla [2018] A. Martínez-Calvo and A. Sevilla, Temporal stability of free liquid threads with surface viscoelasticity, J. Fluid Mech. 846, 877 (2018).
  * Wee  _et al._ [2020] H. Wee, B. W. Wagoner, P. M. Kamat, and O. A. Basaran, Effects of surface viscosity on breakup of viscous threads, Phys. Rev. Lett. 124, 204501 (2020).
  * Edwards and Oron [1995] D. A. Edwards and A. Oron, Instability of a non-wetting film with interfacial viscous stress, J. Fluid Mech. 298, 287 (1995).
  * Choudhury  _et al._ [2020] A. Choudhury, V. K. Paidi, S. K. Kalpathy, and H. N. Dixit, Enhanced stability of free viscous films due to surface viscosity, Phys. Fluids 32, 082108 (2020).
  * Shen  _et al._ [2018] L. Shen, F. Denner, N. Morgan, B. van Wachem, and D. Dini, Capillary waves with surface viscosity, J. Fluid Mech. 847, 644 (2018).
  * Elfring  _et al._ [2016] G. J. Elfring, L. G. Leal, and T. M. Squires, Surface viscosity and marangoni stresses at surfactant laden interfaces, J. Fluid Mech. 792, 712 (2016).
  * Stevenson [2005] P. Stevenson, Remarks on the shear viscosity of surfaces stabilised with soluble surfactants, J. Colloid Interface Sci. 290, 603 (2005).
  * Dockar  _et al._ [2018] D. Dockar, M. K. Borg, and J. M. Reese, Mechanical stability of surface nanobubbles, Langmuir 35, 9325 (2018).
  * Delgado-Buscalioni  _et al._ [2008] R. Delgado-Buscalioni, E. Chacon, and P. Tarazona, Hydrodynamics of nanoscopic capillary waves, Phys. Rev. Lett. 101, 106102 (2008).
  * Willis and Freund [2010] A. Willis and J. Freund, Thermal capillary waves relaxing on atomically thin liquid films, Phys. Fluids 22, 022002 (2010).
  * Rekvig  _et al._ [2004] L. Rekvig, B. Hafskjold, and B. Smit, Chain length dependencies of the bending modulus of surfactant monolayers, Phys. Rev. Lett. 92, 116101 (2004).
  * Young and Stone [2017] Y.-N. Young and H. A. Stone, Long-wave dynamics of an elastic sheet lubricated by a thin liquid film on a wetting substrate, Phys. Rev. Fluids 2, 064001 (2017).
  * Timmermans and Lister [2002] M.-L. E. Timmermans and J. R. Lister, The effect of surfactant on the stability of a liquid thread, J. Fluid Mech. 459, 289 (2002).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117, 1 (1995).
  * Kirkwood and Buff [1949] J. G. Kirkwood and F. P. Buff, The statistical mechanical theory of surface tension, J. Chem. Phys. 17, 338 (1949).
  * Kremer and Grest [1990] K. Kremer and G. S. Grest, Dynamics of entangled linear polymer melts: A molecular-dynamics simulation, J. Chem. Phys. 92, 5057 (1990).
  * Zhang  _et al._ [2022] H. Zhang, S. Chen, Z. Guo, and X. Zhang, The fate of bulk nanobubbles under gas dissolution, Phys. Chem. Chem. Phys. 24, 9685 (2022).
  * Tambe and Sharma [1994] D. E. Tambe and M. M. Sharma, The effect of colloidal particles on fluid-fluid interfacial properties and emulsion stability, Adv. Colloid Interface Sci. 52, 1 (1994).
  * Dressaire  _et al._ [2008] E. Dressaire, R. Bee, D. C. Bell, A. Lips, and H. A. Stone, Interfacial polygonal nanopatterning of stable microbubbles, Science 320, 1198 (2008).
  * Zhang  _et al._ [2012] X. Zhang, M. H. Uddin, H. Yang, G. Toikka, W. Ducker, and N. Maeda, Effects of surfactants on the formation and the stability of interfacial nanobubbles, Langmuir 28, 10471 (2012).


# Exact method for calculating the current fluctuations and nonlinear response
of Markovian dynamics

David Andrieux

###### Abstract

We show that the current fluctuations and nonlinear response of Markovian
dynamics can be obtained from a system of polynomial equations. This offers
new opportunities for analytical and numerical results. As an example, we
derive new expressions for the current nonlinear response coefficients.

##  I Context and objectives

Current fluctuations and their nonlinear response capture important dynamical
and thermodynamical properties of stochastic systems. Both fluctuations and
their response are closely intertwined, as revealed by the fluctuation theorem
for currents [1]. More recently, we showed that the nonlinear response takes a
simple, fully symmetric form when the system parameters are varied along
dynamical equivalence classes [2, 3, 4]

However, calculating current fluctutations and these dynamical equivalence
classes remain difficult. In the traditional approach, it requires solving an
eigenvalue and eigenvector problem for each value of the counting parameters.
While this can usually be done numerically efficiently, it is difficult to
gain additional insights using this approach.

In this paper we formulate the calculation of the current generating function
and the associated nonlinear response as a system of polynomial (cubic)
equations. This method is based on a decomposition of Markov dynamics into
cycles [5, 6, 7] combined with a partition into dynamical equivalence classes
[2, 3]. This approach offers an alternative way to numerically evaluate these
quantities, as well as opportunities to derive new analytical insights. As an
example, we use this formulation to derive new expressions for the nonlinear
response coefficients.

##  II Markov chains and decomposition in cycle matrices

We consider a Markov chain characterized by a transition matrix
$P=\left(P_{ij}\right)\in\mathbb{R}^{N\times N}$ on a finite state space. We
assume that the Markov chain is primitive, i.e. there exists an $n_{0}$ such
that $P^{n_{0}}$ has all positive entries. The chain $P$ thus admits a unique
stationary distribution $\pi$.

It will be convenient to refer to the associated matrix $F$ such that
$F_{ij}=\pi_{i}P_{ij}$. Each element $F_{ij}$ corresponds to the steady state
probability flux between state $i$ and $j$. $F$ satisfies

$\displaystyle\sum_{i}F_{ij}=\pi_{j}\quad\text{and}\quad\sum_{j}F_{ij}=\pi_{i}$
(1)

so that $\sum_{ij}F_{ij}=1$. Every matrix $F$ corresponds to a unique chain
$P$ and vice versa.

Cohen [5] and Alpern [6] demonstrated that the matrix $F$ can always be
decomposed as a combination of cycle matrices. If $[a_{1},\cdots,a_{m}]$ is a
sequence of distinct integers chosen from $1,\cdots,N$, then we define the
corresponding cycle matrix as the $N\times N$ matrix
$C^{[a_{1},\cdots,a_{m}]}$ given by
$c_{a_{1}a_{2}}=c_{a_{2}a_{3}}=\cdots=c_{a_{m}a_{1}}=1/m$ and $0$ otherwise.
We say that $\ell=m$ is the length of $C$. Then, for some vector
$(\lambda_{1},\cdots,\lambda_{l})$ such that $\sum_{e}\lambda_{e}=1$ and
$\lambda_{e}\geq 0$, and some cycle matrices $C^{e}$, we have

$\displaystyle F=\sum_{e}\lambda_{e}\,C^{e}\,.$ (2)

Note that the decomposition (2) is not unique. In this paper, we will use a
decomposition of the form (2) with a clear thermodynamic interpretation [8, 7,
9].

We represent the chain $P$ by a graph where each vertex corresponds to a state
and each transition $P_{ij}>0$ by an edge $e$. Then, for a graph with $N$
vertices and $E$ edges with $i\neq j$ and $D$ edges with $i=j$, there exists
$M=E-N+1$ independent thermodynamic currents $J_{\alpha}$, with each current
defined by a fundamental cycle $C^{\alpha}$ (see Ref. [8] for details). We
then have the following result:

Theorem: Denote by $e$ the set of cycles $[i,j]$ with $P_{ij}>0$ (and thus
$P_{ji}>0$ since P is primitive) and by $\alpha$ the fundamental cycles. Every
matrix $F$ is a convex combination of $D+E+M=D+2E-N+1\leq 2E-1$ cycle matrices

$\displaystyle
F=\sum_{e}\lambda_{e}\,C^{e}+\sum_{\alpha}\lambda_{\alpha}\,C^{\alpha}\,.$ (3)

with $\sum_{e}\lambda_{e}+\sum_{\alpha}\lambda_{\alpha}=1$ and
$(\lambda_{e},\lambda_{\alpha})\geq 0$.

Demonstration: Combine results from Refs. [8, 7, 9, 11].

The representation (3) separates the nonequilibrium components
$\lambda_{\alpha}=J_{\alpha}\ell_{\alpha}$ from the ’local’ fluxes
$\lambda_{e}$ that do not contribute to the net currents. In particular, $F$
is symmetric and thus corresponds to an equilibrium dynamics if and only if
$\lambda_{\alpha}=0$ for all $\alpha$.

We can also obtain the time reversal of a Markov chain from the decomposition
(3). The time reversal of a Markov chain is defined as

$\displaystyle P^{R}={\rm diag}(\pi)^{-1}\,P^{T}\,{\rm diag}(\pi)$ (4)

where $\pi$ is the stationary distribution of $P$. Using (3) and noting that
$(C^{e})^{T}=C^{e}$ and $(C^{\alpha})^{T}=C^{\alpha^{R}}$ where $\alpha^{R}$
is the reverse cycle $[a_{m},\cdots,a_{1}]$, we obtain that the associated
flux matrix

$\displaystyle
F^{R}=\sum_{e}\lambda_{e}C^{e}+\sum_{\alpha}\lambda_{\alpha}C^{\alpha^{R}}\,.$
(5)

The time reversed chain $P^{R}$ therefore has the same thermodynamic currents
as $P$ but in the reverse directions. We will use this observation to reduce
the size of computation of the current fluctuations in half.

Example: Transport on a disordered ring. Throughout the paper, we will
illustrate our results by looking at a disordered ring. We consider the Markov
chain $P$ representing a ring of length $\ell=N$ with periodic boundary
condition and no self-transition: $P_{ij}>0$ if $|i-j|=1$ or if
$|i-j|=\ell-1$, and $0$ otherwise. In this case $D=0$ and $E=N$ so that the
system possesses one independent current $J_{\alpha}=J$ associated with the
cycle $\alpha=[1,2,\ldots,\ell]$ if $J>0$ and $[\ell,\ldots,2,1]$ if $J<0$. If
$k$ denotes the cycles $[k,k+1]$ (with periodic boundary conditions), we have
that any such ring dynamics can be written as

$\displaystyle
F=\lambda_{[1,2,...,\ell]}C^{[1,2,...,\ell]}+\sum_{k}\lambda_{[k,k+1]}C^{[k,k+1]}$

where $\lambda_{[1,2,...,\ell]}=J\ell$ and

$\displaystyle
C^{[1,2,...,\ell]}=\frac{1}{\ell}\begin{pmatrix}0&1&0&\cdots&0\\\
0&0&1&\ddots&0\\\ \vdots&&\ddots&\ddots&1\\\ 1&0&\cdots&0&0\\\
\end{pmatrix}\,,\quad\quad C^{[k,k+1]}=\frac{1}{2}\begin{pmatrix}\ddots&&&&\\\
&&0&1&\\\ &&1&0&\\\ &&&&\ddots\\\ \end{pmatrix}$

where $(C^{[k,k+1]})_{k,k+1}=(C^{[k,k+1]})_{k+1,k}=1$ and $0$ otherwise. A
similar decomposition applies when $J<0$ with $C^{[1,2,,...,\ell]}$ replaced
by $C^{[\ell,...,2,1]}$, i.e. the case $J<0$ corresponds to the time reversed
chain.

The steady state probabilities are readily obtained through Eqs (1), giving
$\pi_{i}=\lambda_{[1,2,...,\ell]}/\ell+\lambda_{i,i+1}/2+\lambda_{i-1,i}/2.$

##  III Current fluctutations and generating function

The Markov chain $P$ generates random trajectories $i_{0}\rightarrow
i_{1}\rightarrow\ldots\rightarrow i_{n}$. The fluctuating currents are then
measured by

$\displaystyle G_{\alpha}(n)=\sum_{l=1}^{n}j_{\alpha}(l)\,,$ (6)

where $j_{\alpha}(l)=\pm 1$ if the transition $i_{l-1}\rightarrow i_{l}$
corresponds to the fundamental edge $\alpha$ in the positive (negative)
direction, and $0$ otherwise [8, 12].

The fluctuations of the currents $j_{\alpha}$ can be characterized by the
cumulant generating function

$\displaystyle
q(\boldsymbol{s})=\lim_{n\rightarrow\infty}\frac{1}{n}\ln\left\langle\exp\left({\sum_{\alpha}s_{\alpha}G_{\alpha}(n)}\right)\right\rangle\,.$
(7)

All the cumulants are obtained by successive derivations with respect to the
counting parameters $\boldsymbol{s}$.

The cumulant generating function (7) is given by [12, 4]

$\displaystyle q(\boldsymbol{s})=\ln\rho[P\circ Z(\boldsymbol{s})]\,,$ (8)

where $\rho(Y)$ denotes the spectral radius (i.e. the modulus of the largest
eigenvalue) of the operator $Y$, $\circ$ is the Hadamar product of two
operators, and

$\displaystyle
Z_{ij}(\boldsymbol{s})\equiv\begin{cases}\exp\left(+s_{\alpha}\right)&\text{if
the transition $i\rightarrow j$ corresponds to the edge $\alpha$ in the
positive direction,}\\\ \exp\left(-s_{\alpha}\right)&\text{if the transition
$i\rightarrow j$ corresponds to the edge $\alpha$ in the negative
direction,}\\\ 1&\text{otherwise}.\end{cases}$

The operator $P\circ Z$ is non-stochastic when $\boldsymbol{s}\neq 0$.
Therefore, its spectral radius, and thus the generating function, often cannot
be resolved analytically. Obtaining the generating function then requires
calculating the Perron eigenvalue of the operator $P\circ Z$. Note that, as
$|\boldsymbol{s}|$ increases, some elements of $P\circ Z$ become exponentially
large while other become exponentially small, making the numeral evaluation of
$q$ difficult.

The Perron eigenvectors $x_{\boldsymbol{s}}$ define a path in the space of
Markov dynamics through

$\displaystyle P^{*}(\boldsymbol{s})=\frac{1}{\rho(\boldsymbol{s})}{\rm
diag}(x_{\boldsymbol{s}})^{-1}\,[P\circ Z(\boldsymbol{s})]\,{\rm
diag}(x_{\boldsymbol{s}})\,.$ (9)

These dynamics play a special role in the nonlinear response theory [2, 3]. In
particular, the counting parameters $\boldsymbol{s}$ correspond to the
affinities of these dynamics: $\boldsymbol{A}[P^{*}]=2\boldsymbol{s}$. We will
return to these observations in the next sections.

Example (continued). For the disordered ring, the single independent current
can be measured along any edge $i\rightarrow i+1$. Choosing the fundamental
transition as $\alpha=\ell\rightarrow 1$ the current generating function is
given by the largest eigenvalues of

$\displaystyle P\circ Z(s)=\begin{pmatrix}0&P_{1,2}&0&\cdots&P_{1,\ell}\,{\rm
e}^{-s}\\\ P_{21}&0&P_{23}&\ddots&0\\\
\vdots&&\ddots&\ddots&P_{\ell-1,\ell}\\\ P_{\ell,1}\,{\rm
e}^{+s}&0&\cdots&P_{\ell,\ell-1}&0\\\ \end{pmatrix}\,.$

The affinity of the resulting dynamics (9) takes the value $A=2s$.

##  IV Dynamical equivalence classes and alternative method to calculate the
current fluctuations

In this section we derive an alternative method to calculate the current
generating function (7). To this end, we build on the fact that the dynamics
(9) are related through the relation

$\displaystyle
P^{*}\circ(P^{*})^{T}(\boldsymbol{s})=\gamma(\boldsymbol{s},\boldsymbol{s}^{\prime})\,P^{*}\circ(P^{*})^{T}(\boldsymbol{s}^{\prime})$
(10)

where $\gamma=\rho^{2}(\boldsymbol{s}^{\prime})/\rho^{2}(\boldsymbol{s})$.
This relation defines equivalence classes in the space of Markov dynamics [2].
In Ref. [2] we also proved that the factor $\gamma$ is related to the current
generating function by $q=-(1/2)\log\gamma+C$.

We now consider the equivalence class $[P]=\\{P^{*}(\boldsymbol{s})\\}$ and
how to generate it using the decomposition (3).

Theorem: Let $E=P\circ P^{T}$. The equivalence class $[P]$ is parametrized by
the independent currents $\lambda_{\alpha}=J_{\alpha}\ell_{\alpha}$. In
particular the associated dynamics $F(\lambda_{\alpha})\in[P]$ satisfy

$\displaystyle
F_{ij}F_{ji}=\gamma\,E_{ij}\,\Big{(}\sum_{j}F_{ij}\Big{)}\Big{(}\sum_{i}F_{ij}\Big{)}$
(11)

for some $\gamma$.

Demonstration: The equivalence class (10) is defined by
$[P]=\\{T|\,T_{ij}T_{ji}=\gamma E_{ij}\\}$ for some $\gamma$. Inserting the
definition $F_{ij}=\pi_{i}T_{ij}$ and using that $\sum_{j}F_{ij}=\pi_{i}$ we
obtain equation (11).

This result can now be used to obtain an alternative method to calculate the
current fluctuations. For simplicity we start from an equilibrium dynamics
$\bar{P}$, the method can be directly extended to cover the case where the
starting dynamics is out of equilibrium. In addition, we showed that the
generating function of an arbitrary nonequilibrium dynamics can be obtained as
the translation of an equilibrium generating function [3].

Alternative method for the calculation of current fluctuations

Let $\bar{P}$ be an equilibrium dynamics,
$\bar{F}=\sum_{e}\bar{\lambda}_{e}C^{e}$ its associated flux matrix, and
$E=\bar{P}\circ\bar{P}^{T}$. We will vary the thermodynamic currents
$\lambda_{\alpha}=J_{\alpha}\ell_{\alpha}\in[0,1]$. For each value of these
currents, we consider the $D+E+1$ variables
$(\boldsymbol{\lambda}_{e},\gamma)$.

According to (11), these $D+E+1$ variables satisfy the system of $D+E+1$
polynomial (cubic) equations

$\displaystyle\begin{cases}[\lambda_{e}/2+\sum_{\alpha}\max(\epsilon_{\alpha}(e),0)\lambda_{\alpha}/\ell_{\alpha}][\lambda_{e}/2-\sum_{\alpha}\min(\epsilon_{\alpha}(e),0)\lambda_{\alpha}/\ell_{\alpha}]=\gamma\,E_{ij}\pi(i)\pi(j)\\\
\sum_{e}\lambda_{e}+\sum_{\alpha}\lambda_{\alpha}=1\,,\end{cases}$ (12)

where for notation simplicity we introduced
$\pi(i)=\sum_{j}F_{ij}=\sum_{j}[\sum_{e}\lambda_{e}(C^{e})_{ij}+\sum_{\alpha}\lambda_{\alpha}(C^{\alpha})_{ij}]$.
Here the $\max/\min$ functions select the positive and negative currents along
edge $e$, respectively. While the equations (12) look complicated, in practice
they simplify considerably for a given network topology. Also, they can easily
be implemented for numerical evaluation [10].

The generating function for the dynamic
$F(\boldsymbol{\lambda})=F(\boldsymbol{\lambda}_{e},\boldsymbol{\lambda}_{\alpha})$
is then given by $q=-(1/2)\ln\gamma(\boldsymbol{\lambda})$. We can then solve
this system for multiple values of the currents
$\lambda_{\alpha}=J_{\alpha}\ell_{\alpha}$ to obtain the full generating
function. In particular, at equilibrium
$\boldsymbol{\lambda}_{\alpha}=\boldsymbol{J}_{\alpha}=0$ and the solution
takes the form $(\bar{\lambda}_{e},1)$.

To obtain the parametrization of the generating function in terms of the
counting parameters $\boldsymbol{s}$, we simply use the relation
$2\boldsymbol{s}=\boldsymbol{A}(\boldsymbol{J})$ to obtain
$q=-(1/2)\ln\gamma[\boldsymbol{A}(\boldsymbol{J})/2]$.

In other words, we can solve the Perron eigenvalue and eigenvector problem by
solving this polynomial system of $D+E+1$ equations. One advantage of this
procedure is that all variables are bounded, $0\leq\lambda_{e}\leq 1$ and
$0\leq\gamma\leq 1$ (Fig. 1). In this way the system of equation is well-
behaved even for large $\boldsymbol{s}$, which corresponds to
$\gamma\rightarrow 0$ and $\lambda_{k}\rightarrow 0$.

Example (continued). For the disordered ring, the system of $\ell+1$ equations
(12) takes the form

$\displaystyle\begin{cases}\lambda_{k}(\lambda_{k}+2\lambda_{\alpha}/\ell)=\gamma\,E_{k}(\lambda_{k}+\lambda_{k-1}+2\lambda_{\alpha}/\ell)(\lambda_{k}+\lambda_{k+1}+2\lambda_{\alpha}/\ell)\quad{\rm
for}\quad k=1,...,\ell\\\ \sum_{k}\lambda_{k}+\lambda_{\alpha}=1\end{cases}$
(13)

where we used that
$\pi_{k}=\lambda_{k}/2+\lambda_{k-1}/2+\lambda_{\alpha}/\ell$. The current
generating function calculated based on this system of $\ell+1$ variables is
illustrated in Fig. (1).

Figure 1: Transport on a disordered ring. (a) Numerical evaluation of the
system (13) as a function of $\lambda_{\alpha}=J\ell$. At equilibrium
$\gamma=1$ and decreases to $\gamma=0$ at maximal (irreversible) current
$\lambda_{\alpha}=1$. (b) Corresponding current generating function,
calculated as $q=-(1/2)\ln\gamma$ with the method (13) (solid line) and the
eigenvalue problem (8) (triangles). By the fluctuation theorem the generating
function is symmetric (not shown); the dynamics for negative values of the
current are obtained by the time-reversed dynamics (5). The parameters take
the value $N=\ell=5$, and the equilibrium dynamics is defined by
$\lambda_{[1,2]}=0.35,\lambda_{[2,3]}=0.04,\lambda_{[3,4]}=0.12,\lambda_{[4,5]}=0.1,\lambda_{[5,1]}=0.39,\gamma=1$.

##  V Expressions for the nonlinear response coefficients

The current response theory expands the currents as functions of the
affinities around equilibrium ($\boldsymbol{A}=0)$:

$\displaystyle
J_{\alpha}=\sum_{l=1}^{\infty}\frac{1}{l!}L_{\alpha,\beta_{1}...\beta_{l}}\,A_{\beta_{1}}\cdot\cdot\cdot\
A_{\beta_{l}}\,,$ (14)

where we sum over repeated indices. The linear response $L_{\alpha,\beta}$ can
be expressed as equilibrium correlations [1, 13]. Further away from
equilibrium however, the response coefficients depend on how the system
parameters are varied are not thus not uniquely defined [13]. This can be
easily understood since there are $M=N-E+1$ independent affinities while there
are $D+E+M>M$ independent transition probabilities that define a Markov chain
$P$ with the same network topology. Therefore, there are multiple dynamics
that can achieve a set of affinities $\boldsymbol{A}$ (Figure 2).

Figure 2:  Nonlinear response curves for different variations of parameters as
a function of the current. (a) Intrinsic dynamics (9) giving rise to the
reponse curve (16). (b) Linear decrease
$\lambda_{e}(\lambda_{\alpha})=\bar{\lambda}_{e}(1-\lambda_{\alpha})$. (c)
Quadratic interpolation between $\bar{\lambda}_{e}$ and
$\lambda_{e}(\lambda_{\alpha}^{\prime})=(1-\lambda_{\alpha}^{\prime})/N$ for
$\lambda_{\alpha}\leq\lambda_{\alpha}^{\prime}=0.1$ and
$\lambda_{e}(\lambda_{\alpha}^{\prime})=(1-\lambda_{\alpha})/N$ for
$\lambda_{\alpha}\geq\lambda_{\alpha}^{\prime}$. Note that for some parameters
variations the affinity diverges for current values $\lambda_{\alpha}<1$. For
example, for the path $\lambda_{e}(\lambda_{\alpha})=\bar{\lambda}_{e}$ for
all $e$ except
$\lambda_{1}(\lambda_{\alpha})=\bar{\lambda}_{1}-\lambda_{\alpha}$ the
dynamics become irreversible at $\lambda_{\alpha}=\bar{\lambda}_{1}$ and the
affinity diverges at that point (not shown). The equilibrium parameters
$\bar{\lambda}_{e}$ take the same values as in Figure (1).

By varying the affinities along the equivalence class (10) the nonlinear
response coefficients take the simple form [4]

$\displaystyle L_{\alpha,\beta_{1}...\beta_{l}}$ $\displaystyle=$
$\displaystyle 0\quad\quad\quad\quad\quad\quad\quad\quad\text{if {\it l} is
even}$ (15) $\displaystyle L_{\alpha,\beta_{1}...\beta_{l}}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{2}\right)^{l}K_{\alpha\beta_{1}...\beta_{l}}(\boldsymbol{0})\quad\text{if
{\it l} is odd.}$ (16)

This shows that all response coefficients $L_{\alpha,\beta_{1}...\beta_{l}}$
are fully symmetric in ($\alpha,\beta_{1},...,\beta_{l}$). In particular, we
recover the Onsager symmetry $L_{\alpha,\beta}=L_{\beta,\alpha}$ and the
corresponding Green-Kubo formula,
$L_{\alpha,\beta}=(1/2)K_{\alpha\beta}(\boldsymbol{0}).$ In addition, all
response coefficients are entirely expressed in terms of equilibrium
correlations.

In view of these properties, we propose that this form constitutes the
’intrinsic’ nonlinear reponse of a system. By construction, the system (12)
generates the equivalence class (9). We can thus use it to calculate the
response theory (15)-(16).

In addition, we now show how the decomposition (3) provides analytical
expressions for the response coefficients. It will be easier to consider the
series

$\displaystyle
A_{\alpha}=\sum_{l=1}^{\infty}\frac{1}{l!}B_{\alpha,\beta_{1}...\beta_{l}}\,J_{\beta_{1}}\cdot\cdot\cdot\
J_{\beta_{l}}\,,$ (17)

which expresses the affinities in terms of the currents rather than the
currents as function of the affinities as in (14). We can then reverse the
series expansion to express the coefficients
$L_{\alpha,\beta_{1}...\beta_{l}}$ in terms of the coefficients
$B_{\alpha,\beta_{1}...\beta_{l}}$ [14].

We start by noting that the affinities $A_{\alpha}$ can be written as

$\displaystyle A_{\alpha}=\sum_{e}A_{e}\epsilon_{\alpha}(e)$ (18)

where $\epsilon_{\alpha}(e)$ takes the value $\pm 1$ if the edge $e$ belongs
to $C^{\alpha}$ in the positive (negative) direction and $0$ otherwise [8].
The local affinities $A_{e}$ are in turn expressed in terms of the positive
and negative fluxes along the edge:

$\displaystyle
A_{e}(\boldsymbol{J})=\ln\left(\frac{\lambda_{e}/2+\sum_{\beta}\max[\epsilon_{\beta}(e),0]\lambda_{\beta}/\ell_{\beta}}{\lambda_{e}/2-\sum_{\beta}\min[\epsilon_{\beta}(e),0]\lambda_{\beta}/\ell_{\beta}}\right)\,.$
(19)

Importantly, the $\lambda_{e}$s are function of the currents, i.e.
$\lambda_{e}=\lambda_{e}(\boldsymbol{J})$, since
$\sum_{e}\lambda_{e}+\sum_{\alpha}\lambda_{\alpha}=1$.

We can then insert (19) into (18) and expand in power series of
$J_{\alpha}=\lambda_{\alpha}/\ell_{\alpha}$ to express the affinties in terms
of the currents. After some algebra we obtain the linear response coefficients

$\displaystyle
B_{\alpha,\beta}=\sum_{e}\left(\frac{2}{\bar{\lambda}_{e}}\right)\epsilon_{\alpha}(e)\epsilon_{\beta}(e)\,.$
(20)

These coefficients depend on the equilibrium fluxes $\bar{\lambda}_{e}/2$
only, and are thus uniquely defined. This formula was previously obtained by
Hill [13] and Schnakenberg [8].

The second-order coefficients read

$\displaystyle
B_{\alpha,\beta_{1}\beta_{2}}=-\sum_{e}\left(\frac{2}{\bar{\lambda}_{e}}\right)^{2}\epsilon_{\alpha}(e)\left(\epsilon_{\beta_{1}}(e)\left(\frac{\lambda^{\prime}_{e,\beta_{2}}}{2}\right)-\frac{1}{2}u[\epsilon_{\beta_{1}}(e),\epsilon_{\beta_{2}}(e)]\right)\,.$
(21)

where $\lambda^{\prime}_{\beta}=d\lambda/dJ_{\beta}$ and
$u(a,b)=\max(a,0)\max(b,0)+\min(a,0)\min(b,0)$. To the best of our knowledge,
formula (21) has not been derived before. We can derive similar expressions
for the higher order coefficients using the same procedure, which will then
depend on higher order variations $d^{n}\lambda/dJ^{n}$.

These results illustrates that the nonlinear response coefficients are not
uniquely defined. Indeed, we see that the response coefficients depend on the
path taken in the space of stochastic dynamics (here through
$\lambda^{\prime}$) and which can vary arbitrarily with the current $J$ (while
respecting the constraint
$\sum_{e}\lambda_{e}+\sum_{\alpha}\lambda_{\alpha}=1$). For example, if we
follow the ’intrinsic’ nonlinear response path defined by (9) the second order
coefficient vanishes, $L_{\alpha,\beta_{1}\beta_{2}}=0$, which translates into
$B_{\alpha,\beta_{1}\beta_{2}}=0$.

To recap, these results formulate the calculation of the current fluctuations,
equivalence classes, and nonlinear response using a polynomial system of
equations with bounded variables. This formulation is based on combining (1) a
cycle decomposition of stochastic matrices, (2) a representation of these
matrices in terms of thermodynamic currents, and (3) a partition of
nonequilibrium reponse based on dynamical equivalence classes. This opens the
way for new analyses of the behavior of nonequilibrium systems, here starting
with new expressions for the nonlinear response coefficients.

The decomposition (3) could also be combined with other methods to calculate
the generating function. For example, it could be used to solve the
minimization problem introduced by Shieh, which expresses the calculation of
the spectral radius as the minimization of the Kublack-Leibler distance
between two distributions [15, 16].

Disclaimer. This paper is not intended for journal publication.

## References

  * [1] D. Andrieux and P. Gaspard, J. Chem. Phys. 48, 571 (2004).
  * [2] D. Andrieux, arXiv 1208.5699 (2012).
  * [3] D. Andrieux, arXiv 1212.1807 (2012).
  * [4] D. Andrieux, arXiv:2205.10784 (2022).
  * [5] J. E. Cohen, Ann. Probability 9, 899 (1981).
  * [6] S. Alpern, The Annals of Probability 11, 789 (1983).
  * [7] S. Kalpazidou, Ann. Prob. 22, 703 (1994).
  * [8] J. Schnakenberg, Rev. Mod. Phys, 48, 571 (1976).
  * [9] S. Kalpazidou, Cycle Representations of Markov Processes (Springer, 2006).
  * [10] By construction the system admits a unique solution, although there doesn’t seem to be a simple proof.
  * [11] D. Andrieux, arxiv:2303.05930 (2023).
  * [12] D. Andrieux and P. Gaspard, J. Stat. Phys. 127, 107 (2007).
  * [13] T. L. Hill, Free Energy Transduction and Biochemical Cycle Kinetics (Dover, 2005).
  * [14] See for example https://mathworld.wolfram.com/SeriesReversion.html (retrieved March 2023). Other methods exist to calculate the reverse coefficients.
  * [15] Shi-Feng Shieh, Entropy 13, 2036 (2011).
  * [16] D. Andrieux, arXiv:1201.1232 (2012).