††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Flexibility in noisy cell-to-cell information dynamics

Ismail Qunbar  Racah Institute of Physics, The Hebrew University of Jerusalem,
Jerusalem 9190401, Israel  Michael Vennettilli  AMOLF, Science Park 104, 1098
XG Amsterdam, The Netherlands  Amir Erez<EMAIL_ADDRESS>Racah
Institute of Physics, The Hebrew University of Jerusalem, Jerusalem 9190401,
Israel

###### Abstract

Exchange of molecules allows cells to exchange information. How robust is the
information to changes in cell parameters? We use a mapping between the
stochastic dynamics of two cells sharing a stimulatory molecule, and
parameters akin to an extension of Landau’s equilibrium phase transition
theory. We show that different single-cell dynamics lead to the same dynamical
response—a flexibility that cells can use. The companion equilibrium Landau
model behaves similarly, thereby describing the dynamics of information in a
broad class of models with coupled order parameters.

In the choreography of life, cells respond to environmental stimuli by
detecting chemical signals and secreting molecules. This process, fundamental
to cellular communication, has attracted much attention, and recent advances
in our understanding of information transfer in biological systems include:
quantifying the energy cost of information transmission in coupled receptors
[1]; interpreting positional information and spatial coupling [2, 3];
establishing bounds on the energy cost of transmitting information in various
settings [4]; outlining essential trade-offs between cost and predictive power
[5]. These studies underscore the complexity and efficiency of cellular
communication using molecular exchange and highlight the importance of
understanding how cells manage and optimize information exchange. However,
despite the dynamic nature of cellular environments, a scaling theory for
cell-to-cell information dynamics under changing conditions had yet to be
developed.

Our previous work described a minimal model of feedback in cellular sense-and-
secrete dynamics [6] and investigated how pairs of such cells share
information [7]. For example, upon antigen stimulation, T cells exhibit a
bimodal distribution of doubly phosphorylated ERK (ppERK), a critical protein
that initiates cell proliferation and determines the immune response [8, 9,
10]. We were able to show that a broad class of models, including Schögl’s
second model, can be used to explain and extract key features from single cell
measurements of ppERK [6, 11]. Here, we extend this line of inquiry to explore
how time-varying individual cellular properties impact overall information
exchange. How does information in a minimal sense-and-secrete two-cell system
respond to gradually changing cellular parameters? We demonstrate that
inherent flexibility allows one cell to match the dynamics of the other, in a
choreographed dance, so that in a broad class of models, cell-to-cell
information is maintained and key observables scale universally.

Our minimal model of cellular sense-and-secrete dynamics maps a class of well-
mixed stochastic biochemical feedback models, in steady state, to parameters
analogous to those in Landau’s thermal equilibrium phase transition theory.
Though applicable to other stochastic models, we focused specifically on
mapping the dynamics of Schlögl’s second model [12, 13] as if it were a Landau
theory. Instead of stochastic reaction rates, one may describe the dynamics
with an effective reduced temperature, $\theta$, magnetic field, $h$, and a
magnetization-like order parameter, $m$ [6, 7]. These parameters can be
extracted from biological data without fitting or knowledge of the underlying
molecular details [6]. The stochastic dynamics are never in thermal
equilibrium, and the noise is demographic in nature. For this reason we will
refer to the Schögl model at steady state and the Landau model at thermal
equilibrium as ‘companion’ systems. The critical transition from a finite to a
zero ‘magnetization’ in the Landau theory is equivalent to the bifurcation
point of the stochastic dynamics [14]—between a bimodal state, having both
high and low molecule counts as stable points the dynamics fluctuate about—to
a unimodal state with an intermediate molecule count (_cf_. SI Appendix Sec.
A). Near the transition point, the correlation time of the system diverges,
exhibiting critical slowing down as expected from the Landau theory [11].

A diverging timescale and critical slowing down do not matter much when
considering systems at equilibrium or steady state. In contrast, in a
dynamically changing setting, collective properties are expected to be
influenced by the slowing down of the collective dynamics. The two-cell system
(cells: $X$ and $Y$) has four main parameters: two ‘fields’ ($h_{X},h_{Y}$)
that control overall molecule count bias in each cell and two ‘reduced
temperatures’ ($\theta_{X},\theta_{Y}$) that act as a bifurcation parameter.
The steady state can be described by two collective coordinates: $H,T$. These
collective coordinates dictate both the steady state mutual information
between the cells and the autocorrelation time of the molecule count, which
are closely related [7, 15]. In this two-cell collective state, each
individual cell’s reaction rates can be set away from its critical point,
while the two-cell collective system remains at criticality. Because the
collective state is critical, both the mutual information and the correlation
time are maximized; to gain information at steady state one must ‘pay’ with an
increased correlation time [7]. Importantly, the collective coordinates $H,T$
are steady-state, linearized solutions that need not apply to the system’s
dynamics. Here, we consider the univesal properties of the two-cell response
to time-varying conditions.

Can the system’s dynamical changes be described using the collective
coordinates, disregarding each cell’s specific dynamics? In this manuscript,
we demonstrate that they can be. Therefore, diverse single-cell parameter
trajectories can lead to identical systemic responses, allowing substantial
flexibility in single-cell properties without impacting the two-cell
information. This flexibility in single-cell behavior enhances robustness and
evolutionary adaptability, facilitating the exploration of advantageous
strategies. Furthermore, while our focus is on cell-to-cell communication,
manifested through the stochastic non-equilibrium dynamics of the Schögl
model, our study of the companion Landau model reveals a similar degeneracy in
the analogous equilibrium system. This extension of equilibrium Landau theory
has been used in the context of phase transitions in minerals and surface and
hydration forces [16, 17, 18]. Therefore, this study describes a general class
of models with coupled order parameters, both equilibrium and non-equilibrium,
that bifurcate according to the Ising mean-field universality class.

The rest of this manuscript is organized as follows: after a brief recap of
the model, we compare the steady state of the stochastic dynamics of the
extended Schlögl model to its equilibrium Landau companion. We then
dynamically ramp both cells’ parameters across the bifurcation or critical
point, showing the lag in the dynamical response of the ramped system and its
scaling. Finally, we show how different single-cell trajectories, sharing the
same collective coordinates, lead to identical responses in the information
exchange between the two cells.

Figure 1: Schematic of the two-cell model and its steady state. (a)
Representation of the two-cell version of the stochastic dynamics, which
extend Schlögl’s second model [12], similarly to our previous manuscript [7].
(b) Varying the Monte-Carlo temperature $k_{B}T_{MC}$ in the Landau model
affects some observables, e.g., the variance of $m_{X}$. The Schlögl model
(black dots) exhibits a dependence on $H$ which cannot be captured by tuning
the Landau model’s $k_{B}T_{MC}$, reflecting its non-equilibrium nature. (c)
The mutual information in the Landau model depends only weakly on
$k_{B}T_{MC}$, and closely matches the Gaussian analytic result (black dots),
though at $T=H=0$ the mutual information in the Gaussian case diverges. (d)
Comparison of the mutual information of the two-cell system as a function of
the collective field, $H$. Blue—the non-equilibrium Schögl steady state;
Orange—the extended Landau model at thermodynamic equilibrium. Dashed
black—the analytic results from the Gaussian approximation. Here, both cells
are identical: $h_{x}=h_{y}=h$ and $\theta_{x}=\theta_{y}=0$, with
$n_{c}=1000$.

_Model and mapping to Landau’s equilibrium phase transition theory_ —Within
each cell, biochemical reactions in a complex signaling cascade result in the
net production and degradation of a molecular species of interest. As
illustrated in Fig. 1a, in the first (second) cell, species $X$ ($Y$) can be
produced spontaneously from bath species at rate $k^{+}_{1}$ ($q^{+}_{1}$),
and can be produced nonlinearly at rate $k^{+}_{2}$ ($q^{+}_{2}$) via a
trimolecular reaction involving two existing $X$ ($Y$) species and a bath
species. Species $X$ ($Y$) can be degraded linearly with molecule number at a
rate $k^{-}_{1}$ ($q^{-}_{1}$), or in a reaction involving three existing $X$
($Y$) molecules at rate $k^{-}_{2}$ ($q^{-}_{2}$). In addition to the internal
reactions, $X$ ($Y$) can be exchanged from the neighboring cell at rate
$\gamma_{xy}$ ($\gamma_{yx}$). Mechanistically, this can be through a gap
junction or through diffusion [7].

Reiterating our previous work [6, 11], we use a mapping from Schlögl to
Landau-like parameters. Without exchange, ($\gamma=0$), the deterministic
dynamics corresponding to the reactions in the left cell in Fig. 1a are
$dx/dt=k_{1}^{+}-k_{1}^{-}x+k_{2}^{+}x^{2}-k_{2}^{-}x^{3}$, where we have
neglected the small shifts of $-1$ and $-2$ for large $x$. Defining the order
parameter $m=(x-n_{c})/n_{c}$, we choose $n_{c}$ to eliminate the term
quadratic in $m$, leading to the Landau form [6],

$\frac{dm}{d\tau}=h-\theta m-\frac{m^{3}}{3},$ (1)

where we have defined $n_{c}=k_{2}^{+}/3k_{2}^{-}$,
$\tau=(k_{2}^{+})^{2}t/3k_{2}^{-}$,
$\theta=3k_{1}^{-}k_{2}^{-}/(k_{2}^{+})^{2}-1$, and
$h=9k_{1}^{+}(k_{2}^{-})^{2}/(k_{2}^{+})^{3}-3k_{1}^{-}k_{2}^{-}/(k_{2}^{+})^{2}+2/3$.
The number of molecules in the system is controlled by $n_{c}$, which controls
all scaling properties of the single-cell system, acting as a _finite system
size_ of the equivalent critical Ising system [6, 11].

In steady state, $dm/d\tau=0$. We can thus interpret $m$ as an order parameter
for the single-cell system, $\theta\equiv(T-T_{c})/T_{c}$ as a ‘reduced
temperature’, and $h$ as a dimensionless field. Analogous to the Ising model,
when $h=0$ in the single-cell system, $\theta>0$ corresponds to a unimodal
steady-state distribution, and $\theta<0$ to a bimodal distribution.
Similarly, tuning $h$ biases the distribution to high or low molecule count.
Applying the same mapping to two coupled cells (with $k\to q$ for $Y$) results
in the Landau form,

$\displaystyle\frac{dm_{X}}{d\tau}$ $\displaystyle=$ $\displaystyle
h_{X}-\theta_{X}m_{X}-\frac{m_{X}^{3}}{3}+g_{XY}m_{Y}-g_{YX}m_{X}\,,$
$\displaystyle\frac{dm_{Y}}{d\tau}$ $\displaystyle=$ $\displaystyle
h_{Y}-\theta_{Y}m_{Y}-\frac{m_{Y}^{3}}{3}+g_{YX}m_{X}-g_{XY}m_{Y}\,.$ (2)

In this context, $g_{XY}=3\gamma_{xy}k_{2}^{-}/(k_{2}^{+})^{2}$ and
$g_{YX}=3\gamma_{yx}q_{2}^{-}/(q_{2}^{+})^{2}$ represent intercellular
exchange terms, set to unity throughout this manuscript.

Although Eq. Flexibility in noisy cell-to-cell information dynamics is a
reparameterization of the nonequilibrium stochastic dynamics of the extended
Schögl model (Fig. 1a), one may separately consider it as an extension of the
well-known Landau model (Eq. 1), computed at thermal equilibrium. This
extension of equilibrium Landau theory to two bilinearly-coupled order
parameters has been considered before: in the context of phase transitions in
minerals [16]; surface and hydration forces [17]; and more generally [18].
Here we focus on the biological system, with its non-equilibrium stochastic
dynamics, while in parallel demonstrating similar dynamical scaling results
for the companion Landau theory computed at thermal equilibrium.

Linearizing the deterministic steady state ($\frac{d}{d\tau}=0$) of the Landau
form (Eq. Flexibility in noisy cell-to-cell information dynamics) gives the
collective coordinates,

$\displaystyle T$ $\displaystyle=$
$\displaystyle\theta_{X}\theta_{Y}+g(\theta_{X}+\theta_{Y})\,,$ $\displaystyle
H$ $\displaystyle=$ $\displaystyle
g(h_{X}+h_{Y})+\left(h_{X}\theta_{Y}+h_{Y}\theta_{X}\right)/2\,.$ (3)

At the critical point, $T=H=0$, the two cells exhibit critical slowing down
and maximal mutual information between $X$ and $Y$, regardless of each cell’s
individual $h,\theta$ values [7].

_Results_ —We simulated the Schögl system’s stochastic dynamics using the
Gillespie algorithm [19, 20], and the equilibrium fluctuations of the
companion Landau model with the Metropolis–Hastings Monte Carlo algorithm. The
dynamics in the Landau case are taken as ‘model A’ in the Halperin-Hohenberg
classification [21, 22],

$\frac{dm_{X}}{d\tau}=-\Gamma\frac{\delta L}{\delta m_{X}}+\zeta(\tau)\,,$ (4)

and similarly for $m_{Y}$, with $\Gamma$ setting the relaxation timescale of
the system and
$L=-h_{X}m_{X}+\frac{1}{2}\theta_{X}m_{X}^{2}+\frac{1}{12}m_{X}^{4}-h_{Y}m_{Y}+\frac{1}{2}\theta_{Y}m_{Y}^{2}+\frac{1}{12}m_{Y}^{4}+\frac{1}{2}\,g\,(m_{X}-m_{Y})^{2}$.
The noise, $\zeta(t)$, is Gaussian white noise obeying $\langle\zeta\rangle=0$
and
$\langle\zeta(\tau)\zeta(\tau^{\prime})\rangle=\delta(\tau-\tau^{\prime})D$.
In the equilibrium statistical mechanics sense, to calculate the partition
function, or alternatively to compute expectation values of observables using
the Monte Carlo simulation accept/reject, we require a simulation temperature,
which we refer to as $k_{B}T_{MC}$. The long-time limit of the dynamics as a
Fokker-Planck equation approaches the equilibrium solution with
$P(m_{X},m_{Y})\propto\exp\left(-\frac{2\Gamma L}{D}\right)$ and accordingly a
fluctuation-dissipation relation, $D=2\Gamma k_{B}T_{MC}$ [22]. To set the
equilibrium temperature $k_{B}T_{MC}$ we resort to a heuristic argument: near
the bifurcation point of the dynamics, ($\theta=0,h=0$), we have
$\mbox{mean}[X]\sim\mbox{Var}[X]\sim n_{c}$ and therefore
$\mbox{Var}[m_{X}]=\mbox{Var}[X]/n_{c}^{2}\sim 1/n_{c}$. Away from $H=0$ but
at $|H|\ll g$, the Gaussian approximation gives $\mbox{Var}[m_{X}]\sim
k_{B}T_{MC}$ (_cf._ SI Appendix Eq. S9). Therefore, $k_{B}T_{c}\sim 1/n_{c}$.
However, tuning the precise value of $k_{B}T_{MC}$ in comparison to the
stochastic simulations of the Schögl model is pointless since the steady state
of the Schögl model does not comply with the equilibrium form,
$\exp\left(-\frac{2\Gamma L}{D}\right)$. Therefore, the probability
distribution of $m_{X}$, and an observable such as the variance,
$\mbox{Var}[m_{X}]$, depend on $k_{B}T_{MC}$ in the thermalized Landau model
(Fig. 1b, colored curves). No value of $k_{B}T_{MC}$ in the Landau simulation
will truly capture the nonequilibrium Schögl steady state (Fig. 1b, dashed
black). In this manuscript, setting the precise value of $k_{B}T_{MC}$ is not
necessary since we focus on the mutual information between the two cells, a
quantity that changes only very weakly with $k_{B}T_{MC}$ (Fig. 1c, colored
curves).

The mutual information shown in Fig. 1c is very well approximated by a
Gaussian analytical approximation (dashed black). The approximation linearizes
the dynamics around the deterministic fixed point: whether of the Landau model
or the chemical Langevin description of the Schlögl model. This yields a
multi-dimensional Ornstein-Uhlenbeck process [23]:
$d\vec{X}_{t}=\mathcal{J}\left(\vec{X}_{t}-\vec{\mu}\right)dt+\sigma
d\vec{W}_{t}\,,$ where $\vec{X}$ is the column vector of variables of
interest, magnetizations or molecule numbers, $\mathcal{J}$ is a negative-
definite matrix, $\vec{W}$ is a vector of independent Brownian motions, and
$\sigma$ is a matrix. The solution is a Gaussian process, where the mutual
information is completely determined by the covariance matrix. By computing
the pathwise solution (_cf._ SI Appendix Sec. B and [23]), one finds that the
steady state covariance matrix $\mathcal{C}$ satisfies the Lyapunov equation
$\mathcal{J}\mathcal{C}+\mathcal{C}\mathcal{J}^{T}=-\sigma\sigma^{T}$. This is
a linear system for the coefficients of $\mathcal{C}$, so all that is left is
to determine the matrices $\mathcal{J}$ and $\sigma$ for the two models. For
the Landau model with $0<|H|\ll g$, we find
$I_{\text{Landau}}=-\frac{1}{2}\log\left(1-\rho^{2}\right)$, with
$\rho^{2}=\frac{g^{2}}{(g+\overline{m}^{2})^{2}}$ and
$\overline{m}=\left(3H/2\right)^{1/3}$. Under the Gaussian approximation, the
covariance matrix is proportional to $k_{B}T_{MC}$, so the temperature cancels
out when computing the correlation coefficient. Therefore, to Gaussian order,
the mutual information is insensitive to the choice of $k_{B}T_{MC}$. Since
$I_{\text{Landau}}$ is an even function of $\overline{m}$, the information is
symmetric in $H$. Repeating the process for the Schlögl model, restricting to
$h_{X}=h_{Y}$, we find
$I_{\text{Schl\"{o}gl}}=-\frac{1}{2}\log\left(1-\rho^{2}\right)$, with
$\rho^{2}=\frac{4g^{2}\left(-\overline{m}^{2}+\overline{m}+2\right)^{2}}{\left(2g\left(2\overline{m}^{2}+\overline{m}+2\right)+\left(\overline{m}\left(\overline{m}+2\right)+4\right)\overline{m}^{2}\right)^{2}}$.
This Gaussian approximation for the Schlögl model is not even in
$\overline{m}$, and therefore asymmetric in $H$. As $H\rightarrow 0$, the
Gaussian information diverges, indicating that the cubic terms are needed for
stabilization and suggesting that the mutual information attains a maximum.
For the mathematical details, _cf._ SI Appendix Sec. B and [23].

Setting $\theta_{x}=\theta_{y}=T=0$ and $h_{x}=h_{y}$, and therefore
$H=h_{x}+h_{y}$, we calculated numerically the mutual information between $X$
and $Y$ at steady state for the Schögl model and at equilibrium for the Landau
case. As expected, the mutual information is maximized at the critical point,
$H=0$ (Fig. 1d). The Gaussian approximation captures well the numerical
simulations away from the critical point, including the assymetry in the
Schlögl case, reflecting the slightly higher information at positive $H$ due
to a higher molecule count. Interestingly, despite the typical number of
molecules, $n_{c}$, being the same in both systems, the mutual information for
the equilibrium Landau model is higher than the companion Schlögl model.

When a system is gradually driven through a critical point, critical slowing
down causes a lagged response to the driving. This phenomenon, known as the
Kibble-Zurek (KZ) effect in the statistical physics literature, results in a
lag that follows scaling rules governed by the critical exponents of the
transition point [24, 25, 26, 27, 28, 29, 30, 31, 32]. Does our two-cell
system exhibit KZ scaling when the collective coordinate $H$ is driven across
the transition point? While various driving protocols are possible in a
biological context, when crossing the critical point, terms beyond the
leading-order linear term do not asymptotically alter the critical scaling
[27]. This theoretical advantage allows us to focus on a simple protocol of
$H(t)$ changing linearly with time without losing biological realism—any
reasonable trajectory with $H(t)$ crossing the transition point at a constant
rate should scale the same.

_Hysteresis in the magnetization_ —Maintaining $T=0$, we first let the system
relax at $H=H_{i}=0.1$, then varied $H(t)$ from $H_{i}$ to $H_{f}=-0.1$ at a
rate determined by $\tau_{d}$, i.e., $H(t)=H_{i}+(H_{f}-H_{i})/\tau_{d}$.
Similarly, we initiated at $H_{i}=-0.1$ and reversed the trajectory.
Importantly, to control for the damping effect of finite-molecule number on
the divergence of the correlations, the system size, $n_{c}$, must be scaled
according to $n_{c}\sim\tau_{d}^{4/5}$ [27, 11]. Consistent with KZ theory,
driving across the transition induced hysteresis loops in both the Schlögl and
Landau cases, with loop characteristics dependent on the ramp time,
$\tau_{d}$, and the critical exponents of the transition point. In both cases
when scaled as per KZ predictions, $M\tau_{d}^{\beta/(\nu z+\beta\delta)}\sim
H\tau_{d}^{\beta\delta/(\nu z+\beta\delta)}$, the Ising mean-field values for
the exponents led to the collapse of the disparate hysteresis loops to a
single curve (Fig. 2).

Figure 2: Lagged response of the collective magnetization of the two-cell
system, $M=(m_{X}+m_{Y})/2$, to time-varying $H(t)$ as it is driven across the
transition point at $H=0$. The collective field, $H(t)$ is ramped according to
$H(t)=H_{i}+(H_{f}-H_{i})/\tau_{d}$. The arrows indicate the direction of the
ramp, starting at either $H_{i}=0.1$ or $-0.1$. (a) Schlögl dynamics for
$\tau_{d}=2000,1000,500,250$. To control for finite molecule-number effects,
we scale $n_{c}=4\,\tau_{d}^{4/5}$ [27, 11]. (b) The Schlögl results scaled
according to the Kibble-Zurek scaling prediction. (c) The extended Landau
system with $\tau_{d}=100000,50000,25000,12500$ and
$n_{c}=\frac{1}{5}\,\tau_{d}^{4/5}$. (d) Kibble-Zurek collapse of the Landau
system response. In all figures, as in the Ising mean-field universality
class, $\beta=\frac{1}{2},\nu z=1,\delta=3$.

_Flexibility in dynamics_ —There are many different trajectories in
$\\{h_{X},h_{Y},\theta_{X},\theta_{Y}\\}$ that have the same values for
collective $\\{H,T\\}$. We wanted to discern whether the dynamical response of
the system is sensitive only to $\\{H,T\\}$, even when their derivation
assumed steady state. To compare the dynamical response of the system to
different trajectories, we considered Eq. Flexibility in noisy cell-to-cell
information dynamics with $H(t)=H_{i}+(H_{f}-H_{i})\,t/\tau_{d}$ and $T=0$,
setting $H_{i}=0.1$ and $H_{f}=-0.1$, thereby crossing the two-cell transition
point $H=T=0$. To maintain $T=0$, we set
$\theta_{Y}=-\frac{\theta_{X}}{\theta_{X}+1}$. Therefore, in terms of $H(t)$,
the $\\{h_{X}(t),h_{Y}(t)\\}$ trajectories must obey
$h_{Y}(t)=\frac{1}{1-\theta_{Y}/2}\left[H(t)-h_{X}(t)\left(1-\frac{\theta_{X}/2}{\theta_{X}+1}\right)\right]$.
For the special case of $\theta_{X}=\theta_{Y}=0$ we retrieve
$h_{Y}(t)=H(t)-h_{X}(t)$, where we considered two situations:
$h_{Y}(t)=h_{X}(t)$ as in Fig. 2, and $h_{Y}(t)=h_{X}(t)-0.1$. Furthermore, we
considered a third trajectory where $\theta_{x}(t)=0.04-0.08(t/\tau_{d})$ and
$h_{x}(t)=0.1-0.2(t/\tau)^{2}$. We chose the quadratic dependence on time to
make explicit the freedom in choosing the trajectories of two of the four
control parameters while maintaining the same collective dynamics. Though many
other trajectories are possible, we trust that demonstrating that these three
examples give identical information is enough to establish the generic
flexibility in the cell-to-cell dynamics.

Contrary to the steady-state calculation (Fig. 1) where one tracks the system
over a sufficiently long time—here, to compute the information at each time-
step we averaged over an ensemble of independent stochastic trajectories, a
significant computational effort. We found that all three driving protocols,
which share the same $H(t),T(t)$, indeed have the same mutual information. As
expected, the mutual information lags in response to the changing conditions
(Fig. 3). As with the other results in this manuscript, the degeneracy between
the three protocols is true for both the Schögl dynamics and the equilibrium
simulations of the companion Landau model (using Metropolis Monte-Carlo).

Figure 3: Mutual information for the two-cell system with dynamically changing
parameters. (a) Three different driving protocols are shown, with the same
collective dynamics, $T=0$ and $H(t)=H_{i}+(H_{f}-H_{i})/\tau_{d}$, but
different realizations of
$\\{h_{x}(t),h_{y}(t),\theta_{x}(t),\theta_{y}(t)\\}$. Protocol 1:
$h_{x}(t)=h_{y}(t)$ and $\theta_{x}=\theta_{y}=0$; Protocol 2:
$h_{x}(t)=0.1+h_{y}(t)$ and $\theta_{x}=\theta_{y}=0$; Protocol 3:
$\theta_{x}(t)=0.04-0.08(t/\tau_{d})$ and $h_{x}(t)=0.1-0.2(t/\tau)^{2}$ (see
text). (b) Mutual information calculated from the stochastic dynamics of the
Schlögl model, demonstrating the same behavior of all three protocols (P1, P2,
P3), with $\tau_{d}=250$, showing lag and hysteresis when the transition point
is crossed. Dashes: the steady-state mutual information from Fig. 1d. The
black arrow indicates the direction of the parameter ramp. (c) Mutual
information of the companion Landau model simulated using Metropolis Monte
Carlo, with $\tau_{d}=6400$, showing convergence of the three protocols.
Dashes: the equilibrium mutual information. In all simulations, $n_{c}=1000$.

_Discussion_ —In this manuscript we considered a simple theoretical question
with a minimal model of cell-to-cell communication: two cells that exchange a
molecule and thereby share information, while each cell’s parameters change
with time. We sought to discern whether cell-to-cell information sharing could
be robust to varying individual cell parameters. We showed that indeed there
are different single-cell driving protocols that lead to the same system
response as long as the collective coordinates have the same dynamics. This
suggests a robustness in cell-to-cell communication and a flexibility in the
routes cells can take to achieve the same information-sharing outcomes.
Despite the minimal nature of the model here studied, our results generalize
to all models with either demographic or thermal noise that belong to the
Ising mean-field universality class.

We compared the cell-to-cell stochastic dynamics of the Schögl model with the
companion equilibrium system—an extension of Landau theory, that had been used
in the past to describe phase transitions in minerals [16]; surface and
hydration forces [17]; and even more generally [18]. To the best of our
knowledge, this is the first example in the literature of KZ scaling of
4-parameter nonequilibrium dynamics, here applied in the context of cell-to-
cell information. Inspired by the analogy between the systems, one may
consider reversing the sign of $g$, making the exchange molecule inhibitory,
rather than excitatory, to the other cell. In the Landau formulation, this
corresponds to anti-ferromagnetic interactions between the two subsystems.
However, in the Schlögl case, one must be careful since a negative $g$ would
imply that when there are zero copies of molecule $X$, having nonzero $Y$
could lead to a negative number of $X$ molecules—a contradiction. Therefore,
mechanistically, one would need to model an ‘interaction’ where $X$ encounters
$Y$.

The dependence of the two-cell system on its low-dimensional representation
could potentially allow cells to optimize other aspects of their function
without compromising information exchange, a possible mechanism for
evolutionary adaptation. Indeed, such degeneracy in single-cell configurations
obeying collective coordinates could extend beyond two-cell systems to entire
tissues or even larger biological systems.

_Data availability_ —All code and data used for this manuscript are freely
available to download in https://github.com/AmirErez/TwocellInformationPy.

_Acknowledgments_ —Numerical simulations were paid for by AE’s startup funds.
MV was supported by the Dutch Research Council, NWO Groot grant, number
2019.085.

### References

  * Ngampruetikorn  _et al._ [2020] V. Ngampruetikorn, D. J. Schwab, and G. J. Stephens, Energy consumption and cooperation for optimal sensing, Nature Communications 11, 975 (2020).
  * Sokolowski and Tkačik [2015] T. R. Sokolowski and G. Tkačik, Optimizing information flow in small genetic networks. IV. Spatial coupling, Physical Review E 91, 062710 (2015).
  * Tkačik and Gregor [2021] G. Tkačik and T. Gregor, The many bits of positional information, Development 148, dev176065 (2021).
  * Bryant and Machta [2023] S. J. Bryant and B. B. Machta, Physical Constraints in Intracellular Signaling: The Cost of Sending a Bit, Physical Review Letters 131, 068401 (2023).
  * Tjalma  _et al._ [2023] A. J. Tjalma, V. Galstyan, J. Goedhart, L. Slim, N. B. Becker, and P. R. ten Wolde, Trade-offs between cost and information in cellular prediction, Proceedings of the National Academy of Sciences 120, e2303078120 (2023).
  * Erez  _et al._ [2019] A. Erez, T. A. Byrd, R. M. Vogel, G. Altan-Bonnet, and A. Mugler, Universality of biochemical feedback and its application to immune cells, Physical Review E 99, 022422 (2019).
  * Erez  _et al._ [2020] A. Erez, T. A. Byrd, M. Vennettilli, and A. Mugler, Cell-to-Cell Information at a Feedback-Induced Bifurcation Point, Physical Review Letters 125, 048103 (2020).
  * Vogel  _et al._ [2016] R. M. Vogel, A. Erez, and G. Altan-Bonnet, Dichotomy of cellular inhibition by small-molecule inhibitors revealed by single-cell analysis, Nature Communications 7, 12428 (2016).
  * Erez  _et al._ [2018] A. Erez, R. Vogel, A. Mugler, A. Belmonte, and G. Altan-Bonnet, Modeling of cytometry data in logarithmic space: When is a bimodal distribution not bimodal?, Cytometry Part A 93, 611 (2018).
  * Altan-Bonnet and Germain [2005] G. Altan-Bonnet and R. N. Germain, Modeling T Cell Antigen Discrimination Based on Feedback Control of Digital ERK Responses, PLOS Biology 3, e356 (2005).
  * Byrd  _et al._ [2019] T. A. Byrd, A. Erez, R. M. Vogel, C. Peterson, M. Vennettilli, G. Altan-Bonnet, and A. Mugler, Critical slowing down in biochemical networks with feedback, Physical Review E 100, 022415 (2019).
  * Schlögl [1972] F. Schlögl, Chemical reaction models for non-equilibrium phase transitions, Zeitschrift für Physik 253, 147 (1972).
  * Grassberger [1982] P. Grassberger, On phase transitions in Schlögl’s second model, Zeitschrift für Physik B Condensed Matter 47, 365 (1982).
  * Bose and Ghosh [2019] I. Bose and S. Ghosh, Bifurcation and criticality, Journal of Statistical Mechanics: Theory and Experiment 2019, 043403 (2019).
  * Vennettilli  _et al._ [2020] M. Vennettilli, A. Erez, and A. Mugler, Multicellular sensing at a feedback-induced critical point, Physical Review E 102, 052411 (2020).
  * Salje [1992] E. K. H. Salje, Application of Landau theory for the analysis of phase transitions in minerals, Physics Reports 215, 49 (1992).
  * Kornyshev  _et al._ [1992] A. A. Kornyshev, D. A. Kossakowski, and S. Leikin, Surface phase transitions and hydration forces, The Journal of Chemical Physics 97, 6809 (1992).
  * Kornyshev  _et al._ [1995] A. A. Kornyshev, D. A. Kossakowski, and S. Leikin, Landau Theory of a System with Two Bilinearly Coupled Order Parameter sin External Field: Exact Mean Field Solution, Critical Properties and Isothermal Susceptibility, Zeitschrift für Naturforschung A 50, 789 (1995).
  * Gillespie [1976] D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics 22, 403 (1976).
  * Gillespie [2000] D. T. Gillespie, The chemical Langevin equation, The Journal of Chemical Physics 113, 297 (2000).
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Reviews of Modern Physics 49, 435 (1977).
  * Goldenfeld [1992] N. Goldenfeld,  _Lectures on Phase Transitions and the Renormalization Group_ , Vol. 5 (Addison-Wesley Reading, MA, 1992).
  * Klebaner [2012] F. C. Klebaner,  _Introduction to Stochastic Calculus with Applications_ (World Scientific Publishing Company, 2012).
  * Kibble [1976] T. W. B. Kibble, Topology of cosmic domains and strings, Journal of Physics A: Mathematical and General 9, 1387 (1976).
  * Zurek [1985] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985).
  * Biroli  _et al._ [2010] G. Biroli, L. F. Cugliandolo, and A. Sicilia, Kibble-Zurek mechanism and infinitely slow annealing through critical points, Physical Review E 81, 050101 (2010).
  * Chandran  _et al._ [2012] A. Chandran, A. Erez, S. S. Gubser, and S. L. Sondhi, Kibble-Zurek problem: Universality and the scaling limit, Physical Review B 86, 064304 (2012).
  * Deffner [2017] S. Deffner, Kibble-Zurek scaling of the irreversible entropy production, Physical Review E 96, 052125 (2017).
  * Fujita  _et al._ [2020] H. Fujita, M. Nishida, M. Nozaki, and Y. Sugimoto, Dynamics of logarithmic negativity and mutual information in smooth quenches, Progress of Theoretical and Experimental Physics 2020, 073B02 (2020).
  * Proverbio  _et al._ [2022] D. Proverbio, A. N. Montanari, A. Skupin, and J. Gonçalves, Buffering variability in cell regulation motifs close to criticality, Physical Review E 106, L032402 (2022).
  * Tarantelli and Vicari [2022] F. Tarantelli and E. Vicari, Out-of-equilibrium dynamics arising from slow round-trip variations of Hamiltonian parameters across quantum and classical critical points, Physical Review B 105, 235124 (2022).
  * Berger  _et al._ [2023] A. Berger, M. G. Harari, A. Gross, and A. Erez, Diversity loss in microbial ecosystems undergoing gradual environmental changes, bioRXiv , 2023.07.07.548183 (2023).
  * Shannon [1948] C. E. Shannon, A mathematical theory of communication, The Bell system technical journal 27, 379 (1948).

## SI Appendix

###  Appendix A A: Visualizing the bifurcation transition

To make explicit the shape of the probability distributions for the molecule
counts of the two cells ($X$,$Y$), we show in Fig. S1 the effect of the
bifurcation parameter, $\theta$, on the two-cell collective state, when $H=0$.
Note the flat-topped shape at the critical point [7].

Figure S1: Examples of the joint distribution $P(X,Y)$ shown as a heatmap,
calculated from Gillespie simulations of the Schögl model [7]. The bifurcation
point at $\theta=0$ separates a bimodal from a unimodal steady state. Color
map corresponds to $\log P$. In all simulations, $h_{X}=h_{Y}=0$ and
$n_{c}=1000$.

###  Appendix B B: Gaussian mutual information

Linearizing a stochastic differential equation around a deterministic fixed
point $\vec{\mu}$ yields an Ornstein-Uhlenbeck process [23]:

$d\vec{X}_{t}=\mathcal{J}\left(\vec{X}_{t}-\vec{\mu}\right)dt+\sigma
d\vec{W}_{t},$ (S1)

where the matrix $\mathcal{J}$ is negative definite and $\vec{X}$ and
$\vec{W}$ are column vectors. Pathwise solutions may be computed as

$\vec{X}_{t}=\vec{\mu}+e^{\mathcal{J}t}\left(\vec{X}_{0}-\vec{\mu}\right)+\int_{0}^{t}e^{\mathcal{J}(t-s)}\sigma
d\vec{W}_{s}.$ (S2)

Note that, for normally distributed initial conditions, this can be
interpreted as a sum of normally distributed variables, so solutions are
normal at all times.

The mutual information [33] between two variables $V_{1}$ and $V_{2}$ is

$I(V_{1},V_{2})=\mathbb{E}\left[\log\left(\frac{P(V_{1},V_{2})}{P(V_{1})P(V_{2})}\right)\right].$
(S3)

If $V_{1}$ and $V_{2}$ have a joint normal distribution, then this becomes

$I(V_{1},V_{2})=-\frac{1}{2}\log\left(1-\frac{\text{Cov}(V_{1},V_{2})^{2}}{\text{Var}(V_{1})\text{Var}(V_{2})}\right).$
(S4)

It suffices to compute the covariance matrix of our process $\vec{X}_{t}$ in
steady state. This can be done by subtracting the mean from Eqn. S2, right
multiplying by it by its transpose, taking expectations and using the Itô
isometry [23], taking $t\rightarrow\infty$, and finally using integration by
parts. The result is that the covariance matrix $\mathcal{C}$ obeys the
Lyapunov equation:

$\mathcal{J}\mathcal{C}+\mathcal{C}\mathcal{J}^{T}=-\sigma\sigma^{T}.$ (S5)

This is a linear system of equations in the coefficients of $\mathcal{C}$ that
may be solved algebraically.

The matrices $\mathcal{J}$ and $\sigma$ can be calculated at the level of the
deterministic dynamics (Eq. Flexibility in noisy cell-to-cell information
dynamics), so it is insensitive to the distinction between the Landau and
Schlögl models. We have done this in previous work [7], where, in the limit
$0<|h_{X}|,|h_{Y}|\ll g=1$, we found

$\overline{m}_{X},\overline{m}_{Y}\sim\left(\frac{3(h_{X}+h_{Y})}{2}\right)^{1/3}=\overline{m}.$
(S6)

We start with the Landau dynamics. These take the form

$\begin{gathered}dm_{X}=-\Gamma\partial_{m_{X}}Ld\tau+\sqrt{2\Gamma
k_{B}T_{MC}}\,dW^{(1)}_{\tau},\\\
dm_{Y}=-\Gamma\partial_{m_{Y}}Ld\tau+\sqrt{2\Gamma
k_{B}T_{MC}}\,dW^{(2)}_{\tau},\\\ \end{gathered}$ (S7)

where the Wiener processes used here are dimensionless (with variance $\tau$).
Linearizing the deterministic dynamics around
$(m_{X},m_{Y})=(\overline{m},\overline{m})$ and evaluating the noise matrix at
that point, we find that

$\mathcal{J}=-\Gamma\begin{bmatrix}g+\overline{m}^{2}&-g\\\
-g&g+\overline{m}^{2}\end{bmatrix},\quad\sigma=\sqrt{2\Gamma
k_{B}T_{MC}}\,\mathbb{I}_{2},$ (S8)

where $\mathbb{I}_{2}$ is the $2\times 2$ identity matrix. Substituting Eq. S8
into the Lyapunov equation (Eq. S5) and solving gives

$\begin{gathered}\text{Var}(m_{X})=\text{Var}(m_{Y})=\frac{k_{B}T_{MC}(g+\overline{m}^{2})}{2g\overline{m}^{2}+\overline{m}^{4}},\\\
\text{Cov}(m_{X},m_{Y})=\frac{gk_{B}T_{MC}}{2g\overline{m}^{2}+\overline{m}^{4}}.\end{gathered}$
(S9)

Therefore, the mutual information gives

$I_{\text{Landau}}=-\frac{1}{2}\log\left(1-\frac{g^{2}}{(g+\overline{m}^{2})^{2}}\right).$
(S10)

Because the covariance matrix was proportional to $T_{MC}$ under this
approximation, this Monte-Carlo completely cancels out from the mutual
information at the Gaussian limit. Additionally, note that this expression is
even in $\overline{m}$, and therefore in $H$.

Now we turn to the case of the Schlögl model, which has demographic noise. The
dynamics under the chemical Langevin approximation [20] can be written as

$\displaystyle
dx_{t}=\left([)k_{1}^{+}-k_{1}^{-}x_{t}+k_{2}^{+}x_{t}^{2}-k_{2}^{-}x_{t}^{3}+\gamma_{YX}y_{t}-\gamma_{XY}x_{t}\right)dt$
(S11)
$\displaystyle+\sqrt{k_{1}^{+}+k_{1}^{-}x_{t}+k_{2}^{+}x_{t}^{2}+k_{2}^{-}x_{t}^{3}}\,dW_{t}^{(1)}$
$\displaystyle-\sqrt{\gamma_{XY}x_{t}}\,dW_{t}^{(3)}+\sqrt{\gamma_{YX}y_{t}}\,dW_{t}^{(4)},$
$\displaystyle
dy_{t}=\left(q_{1}^{+}-q_{1}^{-}y_{t}+q_{2}^{+}y_{t}^{2}-q_{2}^{-}y_{t}^{3}-\gamma_{YX}y_{t}+\gamma_{XY}x_{t}\right)dt$
$\displaystyle+\sqrt{q_{1}^{+}+q_{1}^{-}y_{t}+q_{2}^{+}y_{t}^{2}+q_{2}^{-}y_{t}^{3}}\,dW_{t}^{(2)}$
$\displaystyle+\sqrt{\gamma_{XY}x_{t}}\,dW_{t}^{(3)}-\sqrt{\gamma_{YX}y_{t}}\,dW_{t}^{(4)}.$

In previous work, we inverted the relations between the Landau parameters and
the chemical reaction rates [7]. If we have $n_{c}\gg 1$,
$\theta_{X}=\theta_{Y}=0$, the inverse mapping becomes

$\begin{gathered}q_{1}^{-}=k_{1}^{-},\quad\gamma_{XY}=\gamma_{YX}=gk_{1}^{-},\\\
k_{1}^{+}=n_{c}k_{1}^{-}(h_{X}+1/3),\quad
q_{1}^{+}=n_{c}k_{1}^{-}(h_{Y}+1/3),\\\
k_{2}^{-}=q_{2}^{-}=k_{1}^{-}/(3n_{c}^{2}),\quad
k_{2}^{+}=q_{2}^{+}=k_{1}^{-}/n_{c}.\end{gathered}$ (S12)

We define the mean reactive/ non-exchange propensity for $X$:

$\displaystyle\overline{R}_{X}$
$\displaystyle=k_{1}^{+}+k_{1}^{-}\overline{x}+k_{2}^{+}\overline{x}^{2}+k_{2}^{-}\overline{x}^{3},$
(S13) $\displaystyle=n_{c}k_{1}^{-}\left[h_{X}+1/3+(\overline{m}+1)\right.$
$\displaystyle\qquad\qquad\qquad\qquad\left.+(\overline{m}+1)^{2}+\frac{(\overline{m}+1)^{3}}{3}\right]$
$\displaystyle=n_{c}k_{1}^{-}\left[h_{X}+8/3+4\overline{m}+2\overline{m}^{2}+\frac{\overline{m}^{3}}{3}\right].$

Analogously for $Y$, we have

$\overline{R}_{Y}=n_{c}k_{1}^{-}\left[h_{Y}+8/3+4\overline{m}+2\overline{m}^{2}+\frac{\overline{m}^{3}}{3}\right].$
(S14)

Finally, we define the mean diffusive propensity:

$\overline{D}=\gamma_{XY}\overline{x}=gk_{1}^{-}n_{c}(\overline{m}+1).$ (S15)

The noise matrix in the linearization is the matrix evaluated at
$(x,y)=(n_{c}(\overline{m}+1),n_{c}(\overline{m}+1))$. We find that

$\sigma=\begin{bmatrix}\sqrt{\overline{R}_{X}}&0&-\sqrt{\overline{D}}&\sqrt{\overline{D}}\\\
0&\sqrt{\overline{R}_{Y}}&\sqrt{\overline{D}}&-\sqrt{\overline{D}}\end{bmatrix}.$
(S16)

Evaluating the derivative of the deterministic part at the fixed point gives

$\mathcal{J}=-k_{1}^{-}\begin{bmatrix}g+\overline{m}^{2}&-g\\\
-g&g+\overline{m}^{2}\end{bmatrix},$ (S17)

which is proportional to the result from the Landau case. It is readily
apparent that $\sigma$ is not even in $\overline{m}$, so we expect asymmetry
in $H$ with demographic noise. The case with $h_{X}=h_{Y}$ allows some
simplification:

$\displaystyle I_{\text{Schl\"{o}gl}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\log\left(1-\rho^{2}\right)$ (S18)
$\displaystyle\rho^{2}$ $\displaystyle=$
$\displaystyle\frac{4g^{2}\left(-\overline{m}^{2}+\overline{m}+2\right)^{2}}{\left(2g\left(2\overline{m}^{2}+\overline{m}+2\right)+\left(\overline{m}\left(\overline{m}+2\right)+4\right)\overline{m}^{2}\right)^{2}}$

In the general case where $h_{X}\neq h_{Y}$, the expression is messy and is
studied numerically.


# Influence of substitution on the optical properties of functionalized
pentacene monomers and crystals: Experiment and theory

Y. Saeed1, K. Zhao1, N. Singh1,2, R. Li1, J. E. Anthony3, A. Amassian1, and U.
Schwingenschlögl1 1KAUST, Physical Science & Engineering division, Thuwal
23955-6900, Kingdom of Saudi Arabia  2Solar and Photovoltaic Energy Research
Center, KAUST, Thuwal 23955-6900, Kingdom of Saudi Arabia  3Department of
Chemistry, University of Kentucky, Lexington, Kentucky 40506-0055

###### Abstract

The influence of solubilizing substitutional groups on the electronic and
optical properties of functionalized pentacene molecules and crystals have
been investigated. Density functional theory is used to calculate the
electronic and optical properties of pentacene, TIBS-CF3-pentacene, and TIPS-
pentacene. The results are compared with experimental absorption spectra of
solutions and the complex dielectric function of thin films in the 1 eV to 3
eV energy range. In all cases, the band gaps of the isolated molecules are
found to be smaller than those of the crystals. The absorption spectra and
dielectric function are interpreted in terms of the transitions between the
highest occupied molecular orbitals and lowest unoccupied molecular orbitals.
The bands associated to C and Si atoms connecting the functional side group to
the pentacene in the (6,13) positions are found to be the main contributors to
the optical transitions. The calculated dielectric functions of thin films
agree with the experimental results. A redshift is observed in crystals as
compared to molecules in experiment and theory both, where the amplitude
depends on the packing structure.

##  I Introduction

The design of molecular semiconductors is increasingly important for the
development of organic electronics and organic photovoltaics (OPV) Zaumseil ;
WANG ; DODABALAPOR . Early on pentacene had proven to be one of the best
performing molecular semiconductors, as vacuum-deposited organic thin film
transistors have achieved a mobility as high as 6 cm2 V-1 s-1 Li ; Jurchescu ;
Kim . However, it did not lend itself well to solution processing, which is
believed to be key for low-cost manufacturing of organic semiconductors. In
recent years, chemical modification of the acene has made it possible to
overcome the low solubility and poor stability in solution, whilst maintaining
or enhancing the inter-molecular orbital overlap Anthony ; Subramanian .

Functional substitution of pentacene has been shown to induce favorable
crystal packing motifs for both electronic and OPV applications Subramanian ;
Shu . For example, pentacene without substitution shows a two-dimensional (2D)
herringbone packing motif (see Fig. 1a), while the popular compound
(6,13)-bis(tri-iso-propyl-silyl-ethynyl)-pentacene (TIPS-Pn) shows a brickwork
2D crystal packing (see Fig. 1c) Ostroverkhova . The latter is currently one
of the organic semiconductors exhibiting the highest field-effect mobility
Jackson , with recent carrier mobility reports exceeding 4 cm2 V-1 s-1 Bao-
Nature . When changing substitution from tri-iso-propyl-silyl-ethynyl to tri-
iso-butyl-silyl-ethynyl and introducing a tri-fluoro-methyl group on the acene
backbone to modify the energy levels, to get 2-tri-fluoro-methyl-(6,
13)-bis-(tri-iso-butyl-silyl-ethynyl)-pentacene (TIBS-CF3-Pn), the crystal
packing changes to one-dimensional (1D) sandwich herringbone (see Fig. 1b).
This molecule is found to perform as one of the best non-fullerene acceptor
molecules when mixed with P3HT donor polymer, yielding a power conversion
efficiency of 1.28% Shu .

The substitutional chemistry employed affects the electronic properties of the
monomer as well as of the solid state material itself. Meng _et al._ Meng
have demonstrated that adjusting the alkyl substitution to the four terminal
positions (2, 3, 9, and 10) of the pentacene chromophore shifts the energies
of both the highest occupied molecular orbital (HOMO) and lowest unoccupied
molecular orbital (LUMO) without significantly changing the gap between these
two. When substituting all hydrogen atoms of pentacene with fluorine atoms
some interesting changes of the extinction coefficient are found as the
optical band gap is redshifted Hinderhofer . Recently, Lim _et al._ lim  have
shown that the HOMO–LUMO energy levels can be tuned by varying the number of
nitrile groups in cyano-pentacene substitution. From the above reports, a
close relationship appears to exist between the substitution on the pentacene
chromophore and its electronic and optical properties. To tailor and improve
these properties, one should first understand the correlation between the
chemical modification (like silyl-ethynyl substitution as in the cases of
TIPS-Pn and TIBS-CF3-Pn) and the physical properties of the derivatives both
in monomer and in crystalline states.

The electronic and optical properties of pentacene in solution have been
previously calculated using first principles methods Tiago-DFT . The
calculated optical spectra of the vapor phase are found to be in agreement
with the measurements performed on the thin film phase of pentacene. Doi _et
al._ Doi-DFT  have calculated the electronic band structures for both the
single crystal and thin film polymorphs of pentacene and concluded that the
effective mass of the electrons or holes is larger in the single crystal. A
first principles simulation of the thin film phase of pentacene shows a
crucial dependence of the bandwidths of the HOMO and LUMO and of the band gap
on the molecular stacking angles Parisse . The electronic structures of
iodine- and rubidium-doped pentacene molecular crystals have also been
investigated by ab-initio calculations based on the ultrasoft pseudopotential
method, predicting a metallic behavior Shichibu-DFT . Recently, the structural
and electronic properties of pentacene multilayers on the Ag(111) surface have
been studied, revealing that pentacene has no electronic contribution at the
Fermi level Mete .

Figure 1: Molecular structure (top) and crystal packing (bottom) of (a)
pentacene, (b) TIBS-CF3-Pn, and (c) TIPS-Pn .

Despite several experimental and theoretical investigations, the influence of
the solubilizing chemical substitutions and the resulting changes of the
crystal packing on the electronic and optical properties of pentacene have not
been reported. In this work, we study and compare the theoretical (single
molecule and single crystal) and experimental (dissolved and thin film
polycrystal) optical properties of pentacene, TIPS-Pn, and TIBS-CF3-Pn. The
results are analyzed in light of the calculated density of states (DOS) to see
the influence of different alkyl-silyl groups on the electronic and optical
properties of both monomer and crystal of these materials useful to in
electronic and OPV applications.

##  II Experiments and Characterization

Pentacene, toluene (anhydrous 99.8%) and 1,3,5-trichlorobenzene (anhydrous
99%) were purchased from Sigma Aldrich and used without further purification.
TIPS-Pn and TIBS-CF3-Pn were synthesized Chaung-exp-Anthony  and purified by
multiple recrystallization from acetone (TIPS-Pn) or ethanol (TIBS-CF3-Pn).
Pentacene was dissolved in 1,3,5-trichlobenzene at 100∘C with a concentration
of 0.5 wt.% and stirred overnight in the dark. TIPS-Pn and TIBS-CF3-Pn were
dissolved in toluene at room temperature and stirred overnight in the dark.
The solutions were filled in a 1 mm thick quartz cuvette and loaded in a Cary
5000 (Varian) instrument to aquire UV-vis absorption spectra. The measurements
were performed over a spectral range from 300 nm to 2000 nm with a 2.0 nm slit
width. Single crystal Si(100) wafers with a thermal oxide layer of 100 nm
thickness were used as substrate for the thin film deposition. Prior to
deposition, the substrates were cleaned in amonium hydroxide (30% NH4OH),
hydrogen peroxide (30% H2O2) and Milli Q (1:1:5 ratio) for 15 min at 70∘C.
Thin films of TIPS-Pn and TIBS-CF3-Pn were spin cast at 1000 rpm for 30
seconds in a N2-filled glove-box and left to dry in inert atmosphere at room
temperature. The optical properties of the spin-coated films were measured
using variable angle spectroscopic ellipsometry (VASE) based on the M-2000XI
rotating compensator configuration (J. A. Woollam Co. Inc). VASE spectra
ranging from 0.734 eV to 5.895 eV were recorded at a 18∘ angle of incidence
with respect to the substrate normal from 45∘ to 80∘ with 2∘ increment. In the
paper, we focus on the spectral range from 1 eV to 3 eV. Optical analysis of
VASE data was performed using the EASETM and WVASE32 software packages (J. A.
Woollam Co. Inc). Optical modeling was performed assuming a homogeneous thin
film exhibiting uniaxial anisotropy. To describe the dielectric behavior, a
general oscillator approach consisting of Gaussian peaks in the imaginary part
of the dielectric function $\varepsilon_{2}(E)$ was applied (more detailed
information about the fitting procedure and the Gaussian parameters can be
found in the supporting information). All optical measurements were performed
at room temperature in ambient air.

##  III Simulations

Our calculations are based on density functional theory, using the full-
potential linearized augmented plane wave (FP-LAPW) approach as implemented in
the WIEN2k code wien2k . This approach describes the ground state of the
present compound with high accuracy udo1 . On the other hand, calculation of
optical spectra, in principle, involves excited states. Thus, additional
approximations have to be introduced, which, however, do not compromise the
following line of reasoning nsingh1 . Exchange and correlation effects are
treated within the local density approximation LDA . In the FP-LAPW method,
the unit cell is divided into two parts: non-overlapping atomic spheres
centered at the atomic sites and the interstitial region. The convergence
parameter $R_{mt}\cdot K_{max}$, where $K_{max}$ is the plane wave cut-off and
$R_{mt}$ is the smallest of the atomic sphere radii, controls the size of the
basis set. It is set to $R_{mt}\cdot K_{max}=5$ with $G_{max}=24$. A mesh of
48 uniformly distributed k-points in the irreducible wedge of the Brillouin
zone is used for calculating the electronic properties and a dense mesh of 112
k-points is used to calculate the optical properties. A total energy
convergence of at least $10^{-5}$ Ry is achieved. The experimental lattice
parameters of pentacene ($a=5.959$ Å, $b=7.596$ Å, $c=15.610$ Å,
$\alpha=81.25^{\circ}$, $\beta=86.56^{\circ}$, and $\gamma=89.90^{\circ}$),
TIPS-Pn ($a=7.565$ Å, $b=7.750$ Å, $c=16.835$ Å, $\alpha=89.15^{\circ}$,
$\beta=78.42^{\circ}$, and $\gamma=83.63^{\circ}$), and TIBS-CF3-Pn
($a=17.203$ Å, $b=16.552$ Å, $c=18.168$ Å, $\alpha=90^{\circ}$,
$\beta=113.34^{\circ}$, and $\gamma=90^{\circ}$) are used.

Figure 2: Comparison of the pentacene DOS with TIBS-CF3-Pn and TIPS-Pn in both
molecular and crystalline forms.

##  IV Results and discussion

In Fig. 1, we show the molecular and single crystal packing structures of
pentacene, TIPS-Pn, and TIBS-CF3-Pn. Pentacene (C22H22) and TIPS-Pn
(C44H54Si2) both exhibit triclinic ($P\bar{1}$) crystal symmetries, while
TIBS-CF3-Pn (C57H65F3Si2) has monoclinic ($P21/c$) symmetry.

In Fig. 2, we show the calculated projected DOS for pentacene, TIBS-CF3-Pn,
and TIPS-Pn for the molecule and crystal in the energy range $\pm 3.5$ eV. The
calculated band gaps of pentacene in the monomer and crystalline phases are
found to be 0.84 eV and 0.74 eV, respectively, in agreement with previous
calculations. The low band gap in single crystal pentacene may be due to the
increase in the bandwidth of the HOMO and the LUMO as compared to the monomer
(Fig. 2b). The band gaps in TIBS-CF3-Pn and TIPS-Pn are found to be 1 eV and
0.85 eV, respectively, in the monomer phases, and 0.80 eV and 0.45 eV,
respectively, in the single crystal phases. The band gap of the TIPS-Pn single
crystal is the lowest amongst these molecules owing to the largest LUMO
bandwidth amongst the materials investigated. The calculated band gaps are
lower than in experiments, due to the well known drawback of the local density
approximation. The DOS shows more localized peaks for the pentacene molecule
than its derivatives (Fig. 2). Due to the increase of the HOMO and LUMO
bandwidths from monomers to crystals, the bands overlap (below $E_{F}$) in
agreement with previous calculations for pentacene Doi-DFT . The HOMO and LUMO
consist mainly of bands belonging to the C and Si atoms which connect the side
group to the pentacene chromophore. The H bands appear 3.5 eV below $E_{F}$
for pentacene and its derivatives, while the F bands in TIBS-CF3-Pn lie
between $-6.0$ eV and $-7.5$ eV (not shown). This means that the electronic
states of the H and F atoms do not contribute to the optical transitions in
the visible energy range and do not participate significantly in the
conversion of sunlight into electricity via absorption.

Figure 3: Spectra of $\alpha$, $\varepsilon_{2}$, and $\varepsilon_{1}$ for
pentacene, TIBS-CF3-Pn, and TIPS-Pn in both molecule (monomer) and crystal
form (experiment and simulation).

To establish the effect of the various functional groups on the optical
properties of pentacene and its derivatives, we have calculated their
absorption spectra and dielectric function. The peaks in the optical spectra
are determined by the electric-dipole transitions between the HOMO and LUMO.
Since the local density approximation underestimates the band gap, the
calculated spectra are shifted by difference between experimental and
theoretical band gap to facilitate visual comparison with experimental
spectra. In Fig. 3, we present the absorption coefficient $\alpha(E)$ as well
as the photon energy-dependent complex dielectric function,
$\widetilde{\varepsilon}=\varepsilon_{1}(E)-i\varepsilon_{2}(E)$. An average
of the computed optical spectra along the three coordinate axes is taken and
compared with the average of the experimental optical spectra along the axes
parallel and perpendicular to the plane of substrate. The absorption
coefficients, $\alpha(E)$, of pentacene and its derivatives are calculated for
a single molecule and compared with experimental data of a dilute solution
(see Figs. 3(a,b,c)). In the case of the pentacene solution, the experimental
absorption spectrum shows several peaks at 2.13 eV, 2.30 eV, 2.47 eV, and 2.86
eV. TIBS-CF3-Pn and TIPS-Pn monomers exhibit peaks at 1.94 eV, 2.10 eV, 2.26
eV, 2.45 eV and 2.82 eV. The first absorption peak in pentacene (2.13 eV) is
more intense than the other three peaks, while the first two peaks in TIBS-
CF3-Pn and TIPS-Pn are more intense than the others. A red shift of 0.22 eV
between the first absorption peak of pentacene and of both TIBS-CF3-Pn and
TIPS-Pn, may be due to the effect of the substitutional groups in the latter
two derivatives.

The absorption spectrum of pentacene for the thin film exhibits four
absorption peaks at 1.82 eV, 2.13 eV, 2.30 eV, and 2.86 eV which are
redshifted as compare to the monomer. In the case of TIBS-CF3-Pn, the HOMO-
LUMO absorption bands are shifted to lower energies with respect to the
monomer spectrum, having peaks at 1.87 eV, 2.04 eV, 2.22 eV, and 2.81 eV. A
similar phenomenon is observed for TIPS-Pn thin films which show peaks at 1.78
eV, 1.92 eV, 2.08 eV, and 2.76 eV. The shifts of the HOMO-LUMO absorption band
are different in thin films of these materials with respect to the monomers,
namely 0.20 eV for pentacene, 0.07 eV for TIBS-CF3-Pn, and 0.16 eV for TIPS-
Pn. The different shifts may be attributed to the different crystalline
packing structures. The CF3 group does not have any contribution because the C
and F states (see the DOS) are well below the Fermi level. The calculated and
experimental absorption coefficients also show a redshift between the monomer
and the crystal. The calculations for monomers exhibit a single absorption
peak, while the experiment shows more than one peak which is consistent with
the DOS. The DOS of monomers shows the single sharp LUMO and HOMO bands (allow
only one transition peak) while that of for thin film have wider LUMO and HOMO
bands, which can have more transitions in absorption spectra. This may be due
to the complete isolation of the molecule in the calculation, which may not be
the case in solutions.

The calculated $\varepsilon_{2}(E)$ spectra along with their experimental
counterparts for thin film pentacene and its derivatives are presented in
Figs. 3(d,e,f). The experimental $\varepsilon_{2}(E)$ spectra of TIBS-CF3-Pn
and TIPS-Pn show two initial peaks at 1.85 eV and 1.91 eV, respectively,
reflecting the optical band gap. The optical band gap is redshifted by 0.06 eV
in TIPS-Pn as compared to TIBS-CF3-Pn. The third peak is situated at 2.26 eV
and 2.19 eV for TIBS-CF3-Pn and TIPS-Pn, respectively. Another significant
difference is the intensity of the third peak, which dominates in the
$\varepsilon_{2}(E)$ spectrum of TIBS-CF3-Pn while in TIPS-Pn the first peak
is most prominent. The $\varepsilon_{2}(E)$ spectrum changes dramatically by
introducing a Si-branch in TIBS-CF3-Pn and TIPS-Pn, which is due to the
modified crystal packing. The calculated $\varepsilon_{2}(E)$ spectra of TIBS-
CF3-Pn and TIPS-Pn are in qualitative agreement with our experiments.

The experimental $\varepsilon_{1}(E)$ spectra of TIBS-CF3-Pn and TIPS-Pn thin
films show the first transition peaks at energies of 1.80 eV and 1.86 eV,
respectively, while second and third peaks position remain at the same
energies in both crystals. This reflects that the alky-silyl length results in
changes of the energy state of the first transition peak in
$\varepsilon_{1}(E)$. The subsequent peaks at 2.18 eV might be associated with
a vibronic energy state between the Si HOMO and C LUMO. The calculated
$\varepsilon_{1}(E)$ spectra of crystals and monomers of pentacene, TIBS-
CF3-Pn, and TIPS-Pn show similar characteristics. The $\varepsilon_{1}(E)$
spectra of the TIBS-CF3-Pn and TIPS-Pn crystals demonstrate three peaks
similar to the experimental results. Overall, the calculated
$\varepsilon_{1}(E)$ spectra of single crystals of TIBS-CF3-Pn and TIPS-Pn are
in agreement with our experimental thin film results.

In conclusion, the effects of substitution on the electronic and optical
properties have been discussed based on experiments and theoretical results.
In the monomer state, the alkyl-silyl substitutions result in an energy shift
of 0.22 eV (experimental) in TIBS-CF3-Pn and TIPS-Pn as compared to pentacene.
In the crystal state, the alkyl-silyl substitution contributes to different
packing structures, which leads to a redshift by 0.09 eV in TIPS-Pn as
compared to TIBS-CF3-Pn. The HOMO-LUMO absorption band in thin films is
shifted towards lower energies as compared to the monomer, by 0.07 eV and 0.16
eV for TIBS-CF3-Pn and TIPS-Pn, respectively. Our first principles calculation
of optical spectra have been analyzed in terms of the calculated DOS. The
optical transitions originate primarily from C and Si bands. A redshift is
observed from monomer to crystal for all compounds, where the extent of
redshift depends on the packing structure. Overall, experiment and theory show
a reasonable agreement for the optical spectra.

## References

  * (1) Z. Zaumseil and H. Sirringhaus, Chem. Rev. 107, 1296 (2007).
  * (2) H. Wang and D. Yan, NPG Asia Mater. 2, 69 (2010).
  * (3) A. Dodabalapor, H. E. Katz, L. Torsi, and R. C. Haddon, Science 269, 1560 (1995).
  * (4) X. Li, B. K. C. Kjellander, J. E. Anthony, C. W. M. Bastiaansen, D. J. Broer, and G. H. Gelinck, Adv. Funct. Mater. 19, 3610 (2009).
  * (5) O. D. Jurchescu, S. Subramanian, R. J. Kline, S. D. Hudson, J. E. Anthony, T. N. Jackson, and D. J. Gundlach, Chem. Mater. 20, 6733 (2008).
  * (6) S. H. Kim, M. Jang , H. Yang, J. E. Anthony, and C. E. Park, Adv. Funct. Mater. 21, 2198 (2011).
  * (7) J. E. Anthony, J. Gierschner, C. A. Landis, S. R. Parkin, J. B. Sherman, and R. C. Bakus, Chem. Commun. 45, 4746 (2007).
  * (8) S. Subramanian, S. K. Park, S. R. Parkin, V. Podzorov, T. N. Jackson, and J. E. Anthony, J. Am. Chem. Soc. 123, 9482 (2001).
  * (9) Y. Shu, PhD thesis, University of Kentucky (2011).
  * (10) O. Ostroverkhova, D. G. Cooke, F. A. Hegmann, R. R. Tykwinski, S. R. Parkin, and J. E. Anthony, Appl. Phy. Lett. 89, 192113 (2006).
  * (11) S. K. Park, T. N. Jackson, J. E. Anthony, and D. A. Mourey, Appl. Phys. Lett. 91, 063514 (2007).
  * (12) G. Giri, E. Verploegen, S. C. B. Mannsfeld, S. A. Evrenk, D. H. Kim, S. Y. Lee, H. A. Becerril, A. A. Guzik, M. F. Toney, and Z. Bao, Nature 480, 504 (2011).
  * (13) H. Meng, M. Bendikov, G. Mitchell, R. Helgeson, F. Wudl, Z. Bao, T. Siegrist, C. Kloc, and C. H. Chen, Adv. Mater. 15, 1090 (2003).
  * (14) A. Hinderhofer, U. Heinemeyer, A. Gerlach, S. Kowarik, R. M. J. Jacobs, Y. Sakamoto, T. Suzuki, and F. Schreiber, J. Chem. Phys. 127, 194705 (2007).
  * (15) Y.-F. Lim, Y. Shu, S. R. Parkin, J. E. Anthony, and G. G. Malliaras, J. Mater. Chem. 19, 3049 (2009).
  * (16) M. L. Tiago, J. E. Northrup, and S. G. Louie, Phys. Rev. B 67, 115212 (2003).
  * (17) K. Doi, K. Yoshida, H. Nakano, A. Tachibana, T. Tanabe, Y. Kojima, and K. Okazaki, J. Appl. Phys. 98, 113709 (2005).
  * (18) P. Parisse, L. Ottaviano, B. Delley, and S. Picozzi, J. Phys.: Condens. Matter 19, 106209 (2007).
  * (19) Y. Shichibu and K. Watanabe, Jpn. J. Appl. Phys. 42, 5472 (2003).
  * (20) E. Mete, I. Demiroglu, M. F. Danisman, and S. Ellialtioglu, J. Phys. Chem. C 114, 2724 (2010).
  * (21) Y. S. Chung, N. Shin, J. Kang, Y. Jo, V. M. Prabhu, S. K. Satija, R. J. Kline, D. M. DeLongchamp, M. F. Toney, M. A. Loth, B. Purushothaman, J. E. Anthony, and D. Y. Yoon, J. Am. Chem. Soc., 133, 412 (2011).
  * (22) P. Blaha, K. Schwarz, G. Madsen, D. Kvasicka, and J. Luitz, WIEN2k, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties (TU Vienna, Vienna, 2001).
  * (23) U. Schwingenschlögl and C. Schuster, Phys. Rev. Lett. 99, 237206 (2007); EPL 79, 27003 (2007).
  * (24) N. Singh and U. Schwingenschlögl, Chem. Phys. Lett. 508, 29 (2011).
  * (25) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  * (26) R. D. McCullough, Adv. Mater. 10, 93 (1998)

##  V Supporting Information

The complex dielectric function
$\widetilde{\varepsilon}=\varepsilon_{1}-i\varepsilon_{2}$ is related to the
complex refraction index $\widetilde{n}=n-ik$ by the following equations:
$\varepsilon_{1}=n^{2}-k^{2}$ and $\varepsilon_{2}=2nk$. Here, $n$ and $k$ are
the refractive index and extinction coefficient, respectively. Kramers-Kronig
transformation was used during the model fitting as a constraint:

$\varepsilon_{1}(\omega)=1+\frac{2}{\pi}P\int_{0}^{\infty}\frac{\omega^{{}^{\prime}}\epsilon_{2}(\omega^{{}^{\prime}})}{\omega^{{}^{\prime}2}-\omega^{2}}d\omega^{{}^{\prime}}$
(1)

$\varepsilon_{2}(\omega)=-\frac{2\omega}{\pi}P\int_{0}^{\infty}\frac{\epsilon_{1}(\omega^{{}^{\prime}})-1}{\omega^{{}^{\prime}2}-\omega^{2}}d\omega^{{}^{\prime}}$
(2)

The mean square error was used to quantify the difference between experimental
and model-generated data:

$MSE=\sqrt{\frac{1}{3n-m}\sum_{i=1}^{n}\left[(N_{E_{i}}-N_{G_{i}})^{2}+(C_{E_{i}}-C_{G_{i}})^{2}+(S_{E_{i}}-S_{G_{i}})^{2}\right]}\times
1000$ (3)

where $n$ is the number of wavelengths, $m$ is the number of fit parameters,
and $N=\cos(2\Psi)$, $C=\sin(2\Psi)\cos(\Delta)$, $S=\sin(2\Psi)\sin(\Delta)$.
Where, $\Psi$ and $\Delta$ are the amplitude ratio and phase shift,
respectively.

The $MSE$ generated is 12.12 and 15.5 for TIPS-Pn and TIBS-CF3-Pn with all
angles variable from 45∘ to 80∘, with 2∘ increment, respectively.

Gaussian oscillators produce a Gaussian line shape in $\varepsilon_{2}$:

$\begin{split}&\varepsilon_{2}=\sum_{i}^{n}A_{n}\Bigg{(}\left[\Gamma(\dfrac{E-E_{n}}{\sigma_{n}})+\Gamma(\dfrac{E+E_{n}}{\sigma_{n}})\right]+\\\
&i\cdot\left(\exp\left[-(\dfrac{E-E_{n}}{\sigma_{n}})^{2}\right]+\exp\left[-(\dfrac{E+E_{n}}{\sigma_{n}})^{2}\right]\right)\Bigg{)}\end{split}$
(4)

where $\sigma_{n}=B_{n}/(2\sqrt{\ln(2)})$ and $n$ is the oscillator number,
$A_{n}=\varepsilon_{2}(E_{n})$ is the amplitude, $E_{n}$ (eV) is the center
energy and $B_{n}$ (eV) is the full width at half maximum of the peak. The
function $\Gamma$ is a convergence series that produces a Kramers-Kronig
consistent line shape for $\varepsilon_{1}$.

Table 1: Parameters of the modified Gaussian model obtained by fitting the imaginary part of dielectric function $\varepsilon_{2}(E)$ of TIBS-CF3-Pn. $\varepsilon_{2xx}(E)$=$\varepsilon_{2yy}(E)$ | $\varepsilon_{2zz}(E)$
---|---
$\varepsilon_{\infty}$=1.901$\pm$0.006 |  $\varepsilon_{\infty}$=2.052$\pm$0.015
UV pole amplitude=11.361$\pm$0.315 | UV pole amplitude=4.831$\pm$0.546
UV pole energy=6.883$\pm$0.024 | UV pole energy=6.663$\pm$0.073
A1=0.887$\pm$0.009 | B1=0.125$\pm$0.001 | E1=1.868$\pm$0.001 | A1=0.285$\pm$0.009 | B1=0.159$\pm$0.007 | E1=1.872$\pm$0.003
A2=0.395$\pm$0.002 | B2=0.115$\pm$0.001 | E2=2.040$\pm$0.000 | A2=0.255$\pm$0.021 | B2=0.127$\pm$0.013 | E2=2.050$\pm$0.008
A3=0.156$\pm$0.002 | B3=0.137$\pm$0.002 | E3=2.204$\pm$0.001 | A3=0.563$\pm$0.037 | B3=0.155$\pm$0.011 | E3=2.244$\pm$0.006
A4=0.123$\pm$0.002 | B4=0.206$\pm$0.005 | E4=2.368$\pm$0.004 | A4=3.225$\pm$0.028 | B4=0.670$\pm$0.003 | E4=3.873$\pm$0.008
A5=0.133$\pm$0.003 | B5=0.510$\pm$0.018 | E5=2.838$\pm$0.006 | A5=2.177$\pm$0.099 | B5=0.125$\pm$0.034 | E5=3.795$\pm$0.002
A6=0.110$\pm$0.015 | B6=0.119$\pm$0.001 | E6=3.304$\pm$0.006 | A6=0.252$\pm$0.034 | B6=0.337$\pm$0.017 | E6=4.569$\pm$0.007
A7=0.453$\pm$0.002 | B7=0.119$\pm$0.001 | E7=3.744$\pm$0.000 | A7=0.308$\pm$0.019 | B7=0.657$\pm$0.081 | E7=4.993$\pm$0.018
A8=0.629$\pm$0.001 | B8=0.989$\pm$0.002 | E8=4.187$\pm$0.002 | A8=0.605$\pm$0.005 | B8=0.975$\pm$0.069 | E8=5.847$\pm$0.014
A9=0.235$\pm$0.002 | B9=0.814$\pm$0.011 | E9=5.076$\pm$0.006 |  |  |
A10=0.082$\pm$0.006 | B10=0.145$\pm$0.013 | E10=5.416$\pm$0.006 |  |  |
A11=0.617$\pm$0.004 | B11=0.531$\pm$0.006 | E11=5.649$\pm$0.003 |  |  |
Table 2: Parameters of the modified Gaussian model obtained by fitting the imaginary part of dielectric function $\varepsilon_{2}(E)$ of TIPS-Pn. $\varepsilon_{2xx}(E)$=$\varepsilon_{2yy}(E)$ | $\varepsilon_{2zz}(E)$
---|---
$\varepsilon_{\infty}$=2.012$\pm$0.030 |  $\varepsilon_{\infty}$=2.251$\pm$0.023
UV pole amplitude=4.360$\pm$1.999 | UV pole amplitude=3.889$\pm$0.443
UV pole energy=7.139$\pm$0.139 | UV pole energy=6.190$\pm$0.033
A1=0.830$\pm$0.002 | B1=0.080$\pm$0.002 | E1=1.906$\pm$0.002 | A1=1.087$\pm$0.027 | B1=0.105$\pm$0.003 | E1=1.891$\pm$0.001
A2=0.353$\pm$0.004 | B2=0.080$\pm$0.001 | E2=2.069$\pm$0.000 | A2=0.666$\pm$0.028 | B2=0.111$\pm$0.006 | E2=2.053$\pm$0.006
A3=0.266$\pm$0.002 | B3=0.497$\pm$0.005 | E3=2.223$\pm$0.002 | A3=1.566$\pm$0.055 | B3=0.154$\pm$0.005 | E3=3.503$\pm$0.002
A4=0.508$\pm$0.005 | B4=0.198$\pm$0.003 | E4=4.139$\pm$0.007 | A4=2.709$\pm$0.203 | B4=0.132$\pm$0.011 | E4=3.886$\pm$0.005
A5=1.232$\pm$0.314 | B4=0.204$\pm$0.003 | E5=4.236$\pm$0.009 | A5=1.532$\pm$0.141 | B5=0.159$\pm$0.012 | E5=4.256$\pm$0.013
A6=2.963$\pm$0.028 | B6=0.467$\pm$0.001 | E6=5.666$\pm$0.006 | A6=1.344$\pm$0.047 | B6=0.234$\pm$0.004 | E6=4.132$\pm$0.016
A7=0.396$\pm$0.004 | B7=0.963$\pm$0.022 | E7=3.506$\pm$0.002 | A7=0.467$\pm$0.005 | B7=1.954$\pm$0.045 | E7=5.419$\pm$0.006
A8=0.939$\pm$0.009 | B8=0.586$\pm$0.007 | E8=3.868$\pm$0.003 | A8=0.400$\pm$0.045 | B8=0.160$\pm$0.018 | E8=5.737$\pm$0.008
|  |  | A9=0.835$\pm$0.020 | B9=0.453$\pm$0.141 | E9=4.088$\pm$0.008