# Self-referencing cellular automata: A model of the evolution of information
control in biological systems

Theodore P. Pavlic1, Alyssa M. Adams1, Paul C. W. Davies1 Sara Imari Walker1

1Arizona State University, Tempe, AZ 85287
<EMAIL_ADDRESS>
###### Abstract

Cellular automata have been useful artificial models for exploring how
relatively simple rules combined with spatial memory can give rise to complex
emergent patterns. Moreover, studying the dynamics of how rules emerge under
artificial selection for function has recently become a powerful tool for
understanding how evolution can innovate within its genetic rule space.
However, conventional cellular automata lack the kind of state feedback that
is surely present in natural evolving systems. Each new generation of a
population leaves an indelible mark on its environment and thus affects the
selective pressures that shape future generations of that population. To model
this phenomenon, we have augmented traditional cellular automata with
state/dependent feedback. Rather than generating automata executions from an
initial condition and a static rule, we introduce mappings which generate
iteration rules from the cellular automaton itself. We show that these new
automata contain disconnected regions which locally act like conventional
automata, thus encapsulating multiple functions into one structure.
Consequently, we have provided a new model for processes like cell
differentiation. Finally, by studying the size of these regions, we provide
additional evidence that the dynamics of self/reference may be critical to
understanding the evolution of natural language. In particular, the rules of
elementary cellular automata appear to be distributed in the same way as words
in the corpus of a natural language.

## Introduction

Cellular automata (CA) are model complex systems that combine spatial memory
with relatively simple update rules to produce rich dynamic patterns
(Wolfram,, 2002). In this regard, they can be viewed as models for life.
“Rules” in DNA encode policies that iteratively react to the environment by
modifying it. Thus, over evolutionary time scales, natural selection can
explore the nucleic/acid rule space and amplify those rules which provide
useful functions. With this narrative in mind, researchers have used _in
silico_ artificial selection to explore the CA rule space and amplify rules
with certain computational abilities (Das et al.,, 1994; Mitchell et al.,,
1994; Breukelaar and Bäck,, 2005; Hordijk,, 2013). Moreover, there has been
much interest in understanding the demographic dynamics of these CA rule
populations over their evolutionary history (Mitchell et al.,, 1994; Hordijk,,
2013). That is, artificial selection of these evolving CA’s has itself become
a model for the innovation intrinsic to natural selection.

One major difference between evolving cellular automata (EvCA) and evolving
natural organisms is the lack of feedback in the fitness channel of the
former. In EvCA, each new generation faces the same selective pressures as
prior generations. However, with new generations of natural organisms, there
is feedback between the current demographics and the selective pressures
shaping future demographics. Goldenfeld and Woese, (2011) point out that these
self/referential dynamics are a unique characteristic of life – making life
distinctly different from any other physical system. Two of us (SIW and PCWD)
have proposed that self/referential dynamics are one of the hallmarks of life,
emerging with its origin (Walker and Davies,, 2013). Where EvCA’s will only
innovate in the presence of external “abiotic” pressures, natural organisms
put pressure on themselves to re/organize even without an external fitness
driver. Similarly, at everyday and ontogenetic time scales, conventional CA’s
do not easily embed the regulatory mechanisms that permeate throughout life. A
single genome gives rise to a wide variety of differentiated cell phenotypes
which, locally, appear to follow a consistent set of operational rules but
globally appear to have no shared program. By modeling feedback explicitly,
these phenomena can be explained using gene regulatory network (GRN)
frameworks (Schlitt and Brazma,, 2007; Davidson,, 2010), where the expression
level of one gene promotes or inhibits the expression level of another and
thus “latches” cells into different types. Thus, by making feedback between
state and dynamics explicit in CA’s, it may be possible to enrich their
ability to model how organisms emerge, evolve, develop, and react to both
their external environment and their internal state.

Here, we introduce self/referencing cellular/automaton framework we call
PICARD where _P ICARD Implements CA Rules Differently_ (PICARD). Like a
traditional one/dimensional CA, PICARD executions move from one iteration to
another by some rule. However, whereas traditional CA’s require the rule to be
static and externally specified, PICARD infers the iteration rule from the
current state of the CA itself. As we will show, executions from multiple
static CA’s can be embedded within a single PICARD – a PICARD can be identical
to one static CA from certain initial conditions and another static CA from
other initial conditions. Thus, a PICARD can combine the computational
abilities of different CA’s within one entity. Moreover, whereas CA’s differ
by their rule, PICARD’s differ by their state-to-rule mapping. Because there
are many more state-to-rule mappings than there are CA rules, the PICARD
parameter space is potentially much richer for later evolutionary
investigations.

To some extent, PICARD is a simple attempt to add dynamical feedback that is
missing in traditional evolutionary cellular automata. However, because
iterations are generated by rules that are encoded in previous iterations,
PICARD feedback is the kind of self/reference that is thought to be a
characteristic feature of life (Hofstadter,, 1979; Kataoka and Kaneko, 2000a,
; Kataoka and Kaneko, 2000b, ; Goldenfeld and Woese,, 2011; Walker and
Davies,, 2013). Our CA approach shares many similarities with self/referencing
functional/dynamics developed by Kataoka and Kaneko, 2000a to model the
evolution of rules (Kataoka and Kaneko, 2000a, ; Kataoka and Kaneko, 2000b, ).
Attempting to avoid the biochemical complexities of the evolution of nucleic
acids, they turn their focus on the evolution of natural language.
Furthermore, they draw connections between attractors in their
coupled/logistic/map landscape and words that accumulate in language. Although
their framework is very different than the automata we study here, we too have
results that appear to be strongly connected to the evolution of natural
language. Thus, augmenting cellular automata with self/reference widens our
ability to model to evolution of language in unanticipated ways.

## Cellular Automata with PICARD Mappings

Although PICARD implements CA rules differently, once each rule is defined, a
PICARD iteration is identical to an iteration of a conventional elementary CA.
As shown in Fig. 1(a), a traditional elementary CA generates each row based on
the pattern in the preceding row.

(a) Single iteration of rule 182

(b) Execution of rule 182

Figure 1: Traditional cellular automaton. In 1(a), a single iteration of a
traditional CA is shown in the top two rows, where row 0 is taken as an
initial condition on which a static CA rule operates to produce row 1. The
static CA rule is summarized in the eight groups of four cells in the bottom
of the figure. In each group, the top three cells match adjacent cells in row
0, and the bottom one cell represents the cell generated beneath the row-0
cells in row 1. One such grouping in the two rows is outlined in a broken line
along with the corresponding group from the elementary CA rule below it. This
CA rule can be equivalently summarized as an eight-bit binary string
0b10110110, a two/digit hexadecimal 0xB6, or a decimal 182. In 1(b), an
execution of rule 182 is shown that resembles a Sierpinski triangle.

The iteration rule is a lookup table that maps each triplet of bits in the
preceding row to a single bit in the following row. Thus, with the right
initial conditions, some rules can produce intricate patterns over many
generations, as shown in Fig. 1(b).

Where PICARD differs from a traditional CA is that no static rule is
specified. Instead, a map is provided from each row to the rule that will
operate on it. We are primarily interested in CA’s with more than eight cells;
consequently, this mapping is a coarse graining of the system – some rules
will necessarily correspond to multiple different row patterns. With this
multiple realizability of rules in mind, each PICARD row and rule can be
viewed as a _microstate_ -_macrostate_ pair. Consequently, in the following,
we will use the terms _rule_ and _macrostate_ interchangeably.

For example, Fig. 2 shows two arbitrary examples of microstate/to/rule
mappings.

(a) Sum of ones

(b) Density of ones

Figure 2: Example aggregation PICARD mappings. The CA rule induced by each
mapping is shown between rows 0 and 1, where the triplet above each position
represents the pattern of bits in row 0 that will result in the boxed bit in
row 1. Positions in the rule that result in a changed bit from row 0 to row 1
have been highlighted with broken lines; the resulting toggled bits in row 1
have also been highlighted. Each mapping represents some aggregate property of
the first row. In 2(a), the rule represents the total number of 1’s in row 0,
modulo 256. Because there are only 24 cells in the row, this PICARD mapping
can only induce 25 different elementary CA rules. In the case of a row with
more than 255 cells, every elementary rule is possible. In 2(b), the rule
represents the density of 1’s in row 0 scaled by 255.

Each of these two mappings extract an aggregate property of the corresponding
row. In Fig. 2(a), the focal macrostate is the number of 1’s in the row. In
Fig. 2(b), the focal macrostate is the density of 1’s in the row. For the
analogous case of interacting particles, macroscopic aggregates like
temperature and pressure are not typically viewed as having causal influence
on the microscopic states of the system. That is, macroscopic states and
microscopic states exist in two different levels of description. However, if a
collection of particles is put into contact with a heat bath characterized by
only its temperature, the macroscopic properties of the gas appear to drive
the evolution of the microscopic system. Likewise, these PICARD aggregation
mappings are simultaneously descriptive coarse grainings as well as
prescriptive rules governing the dynamics of the cells. Goldenfeld and Woese,
(2011) argue that self/reference appears in biological systems but apparently
not in physical systems because of similar reasoning – biological phenomena
are emergent and may require self/reference for analysis without appealing to
underlying interacting physical microstates. Self/reference therefore may be
intimately related to framework suggested to characterize emergence, such as
top-down causation (Davies,, 2012), and may play an important role in the
emergence of new organizational levels through major evolutionary transitions
(Walker et al.,, 2012).

Aggregation mappings may have intuitive connections to statistical mechanics,
but PICARD mappings may be built in entirely different ways. The block
mappings in Fig. 3 are formed by clustering bits of the CA microstate into
eight groups and matching each group to a predicate function that determines
the bit in the rule which has the same relative position as the 3-bit group.

(a) Block ones majority

(b) Block ones odd parity

(c) Block transition odd parity

Figure 3: Example block PICARD mappings. Rules are shown just as in Fig. 2. However, the 24 bits in row 0 have been broken into 8 3-bit clusters shown with separating fences. The eight bits of the rule are then induced from a simple 3-bit predicate from each group to the corresponding rule bit. In 3(a), each bit in the rule is a 1 if and only if the corresponding 3-bit cluster in row 0 has a majority of 1’s. In 3(b), the rule bit is set if there are an odd number of 1’s in the cluster. In 3(c), the rule bit is set if there are an odd number of zero/to/one or one/to/zero transitions in the 3-bit cluster (see Table 1). Microstate Block | Number of Transitions | Rule Bit
---|---|---
000 | 0 | 0
001 | 1 | 1
010 | 2 | 0
011 | 1 | 1
100 | 1 | 1
101 | 2 | 0
110 | 1 | 1
111 | 0 | 0
Table 1: Lookup table for transition odd parity. Reading from left to right,
the number of zero/to/one and one/to/zero transitions is counted. The odd
parity of this count is equivalent to the exclusive or of the outer two bits
in the microstate block. Figure 4: Block transition mapping with permutation.
This PICARD mapping is equivalent to the mapping from Fig. 3(c) composed with
a permutation that switches middle two bits of each nybble.

For example, in Fig. 3(a), the 24 bits of the microstate are broken into eight
mutually exclusive groups of three bits. Each 3-bit group is replaced by a
single 1 if it has a majority of 1’s and a 0 otherwise. In Fig. 3(b), groups
are mapped into 1 if they have an odd number of 1’s. Alternatively, in Fig.
3(c), the actual binary identity of the bits in each group is ignored and
instead the number of transitions is used (see Table 1). That is, from left to
right within each 3-bit group, a cell will transition from one level to
another 0, 1, or 2 times. The corresponding rule bit will be a zero unless the
3-bit group only contains a single transition. This predicate is equivalent to
the exclusive or (XOR) of the outer two bits in the microstate group. As these
examples show, PICARD executions starting from the same microstate can be very
different based on the chosen mapping from microstate to rule.

While the predicates used within block mappings may be relatively simple, the
selection of microstate and macrostate bits for each block mapping adds an
additional layer of complexity. For example, a single PICARD block mapping may
use different predicates for each microstate block. Additionally, microstate
blocks may not be mutually exclusive – they may overlap or leave some
microstate bits uncovered. Even when a a single predicate is used over a set
of mutually exclusive, equally sized microstate groups that cover all
microstate bits, the ordering of the microstate groups need not match the
ordering of the corresponding rule bits. In Fig. 4, a mapping of this sort is
shown. That is, the mapping from Fig. 3(c) has been composed with a
permutation. For the remainder of this paper, we will use this mapping to
demonstrate the richness of an individual PICARD.

## A PICARD Case Study

In the following, we consider the mapping described above in Fig. 4. Because
the rule induced by each row can be viewed as a macrostate of the system, it
is useful to consider both the CA microstate dynamics and the rule macrostate
dynamics. In Fig. 5, executions of the PICARD mapping are shown starting from
three different initial conditions. The executions of the CA are shown in the
left column, and the corresponding executions through the macrostate rule
space is shown in the right column.

(a) Execution leading to mutual oscillation

(b) Execution leading to mutual fixed points

(c) Microstate oscillations under invariant rule

Figure 5: Example PICARD executions. Cellular automata on the left are
adjacent to their corresponding rule histories on the right. Executions in
5(a) and 5(b) show CA and rule histories that have equilibria of the same
type. In 5(c), an execution is shown where oscillations in microstate do not
correspond to oscillations in rule. Figure 6: Macrostate representation over
macroexecutions. Each macrostate may be invariant across a range of
disconnected macroexecutions. Every macroexecution terminates in a fixed point
or a periodic cycle of microstates that each correspond to the same
macrostate. Each terminating fixed point or cycle may be reached from a number
of other transient microstates which themselves correspond to that same
macrostate. A complete macroexecution includes the terminating cycle and all
leading transient microstates which share the same macrostate. Using the
PICARD mapping from Fig. 4, the entire microstate space has been explored and
every complete macroexecution has been found. Displayed here as open circles
are the $\log_{2}$/scaled counts of macroexecutions within each macrostate,
omitting macrostates with no macroexecutions. The asterisks represent the
$\log_{2}$-scaled count of macroexecutions after all fixed/point
macroexecutions are removed. The horizontal axis shows the rank of the
corresponding macrostate when ordered by the total number of complete
macroexecutions within it. A broken piecewise/linear line of best fit is shown
with an upper slope near $-1$. Figure 7: Zipf plot of macrostate distribution
over macroexecutions. Displayed are the data from Fig. 6; however, the
horizontal axis is also $\log_{2}$ scaled. The broken piecewise/linear line of
best fit shows that the distribution is approximately Zipfian.

In general, the lower resolution macrostate executions are each visually
similar to their higher resolution microstate. In Fig. 5(a), the microstate
reaches a steady/state oscillation. Given that the rule that governs the
microstate evolution is derived from the microstate itself, it is not
surprising that the macrostate also reaches a steady/state oscillation.
Similarly, in Fig. 5(b), the microstate and macrostate both reach a fixed
point. However, while the CA does not reach its fixed point until generation
16, its macrostate becomes fixed in generation 14. Thus, at generation 14,
this PICARD degenerates into an elementary CA under rule 44 (i.e.,
0b00101100). If all generations before 14 were removed from the history of the
CA, it would be indistinguishable from an elementary CA. Moreover, as shown in
Fig. 5(c), this phenomenon is not restricted to only fixed points. By
generation 8, the rule trajectory becomes fixed on rule 5 (i.e., 0b00000101);
however, the CA trajectory reaches a steady/state oscillation where some cells
are fixed and others cycle. Again, the tail of this CA execution is
indistinguishable from one that would be generated under elementary CA rule 5.
Consequently, this single PICARD mapping is able to represent executions from
multiple elementary CA rules simultaneously while also potentially producing
novel executions.

When a PICARD mapping generates an execution that is in a region of macrostate
invariance, we call that execution a _macroexecution._ Once a macroexecution
has been reached, state feedback is not required to determine the future
trajectory. Thus, macroexecutions move through regions of _locally elementary_
CA space. For the mapping explored in this section, Fig. 6 shows the relative
“size” of each of these locally elementary regions. In particular, for each of
the 256 elementary CA rules, all of the disconnected macroexecutions for which
the rule is invariant were counted. The rules with non/zero counts were ranked
and plotted as the open circles in the semilog frequency distribution in Fig.
6. The asterisks on the plot represent the number of oscillatory
macroexecutions – that is, the asterisks represent the count when the
fixed/point macroexecutions are removed. Thus, the highest ranking elementary
CA rule includes all fixed/point macroexecutions, but the second ranking
elementary CA rule includes no fixed/point macroexecutions.

As shown by the piecewise linear broken line in Fig. 6, the highest ranking
elementary CA rules have counts which peak at $2^{16}$ and are halved with
each increase in rank over a significant range of rules. However, this trend
is not long lasting and may simply be an artifact of the particular PICARD
mapping we have chosen for this case study. In Fig. 7, the same distribution
is shown on a log–log plot and shows evidence that the greater trend is
Zipfian (Zipf,, 1949). Words in a language tend to be Zipfian distributed in
their use – the relative frequency of words is inversely proportional to their
rank. Following an analogy with language, elementary CA rules act like words
which are used in PICARD macroexecutions. Using a functional/dynamics
framework as a model for the evolution of natural language, Kataoka and
Kaneko, 2000b  show how self/referencing function dynamics can act like a
filter that produces a corpus of words that are each fixed points of the
dynamical model. The invariant macrostate CA rules that we describe are
qualitatively similar to this idea, and their representation is consistent
with measured distributions of words in actual natural language.

The data from Fig. 6 have been reproduced in Table 2 for the sixteen
highest/ranking macrostates.

Rule | Macroexecutions |  Non/trivial Macroexecutions
---|---|---
204 (0xCC) | 65536 | 0
51 (0x33) | 32768 | 32768
205 (0xCD) | 10196 | 0
76 (0x4C) | 9360 | 0
236 (0xEC) | 3420 | 0
200 (0xC8) | 2166 | 0
4 (0x04) | 1792 | 0
12 (0x0C) | 1072 | 0
132 (0x84) | 1072 | 0
68 (0x44) | 1040 | 0
140 (0x8C) | 640 | 0
196 (0xC4) | 640 | 0
108 (0x6C) | 603 | 219
201 (0xC9) | 603 | 219
232 (0xE8) | 384 | 0
77 (0x4D) | 320 | 0

Table 2: Sixteen highest ranking macrostates.

.

The two highest/ranking macrostates, rule 204 and rule 51, each have a number
of macroexecutions that is a power of two, $2^{16}$ and $2^{15}$. In the case
of rule 204, its macroexecutions exhaust all of the $2^{16}$ microstates that
map to it. Thus, rule 204 has no transients – it is a region filled entirely
with fixed points. Likewise, because there are no fixed points in rule 51 and
yet $2^{15}$ macroexecutions, each macroexecution must be a period/2 cycle. So
the rule-51 region also has no transient microstates, but it is filled
entirely with these cycles. Moving farther down the ranks, more microstates
become available for forming macroexecutions with longer/period oscillations
within them. For example, Fig. 8 shows a period/6 oscillation embedded within
a rule-109 macroexecution.

Figure 8: Period/6 macroexecution with short transient. A steady/state
rule-109 (0x6D) macroexecution is shown that has a period of six generations
preceded by its longest possible transient. Rule 109 has 31 macroexecutions,
which are all non/trivial.

Because each of the $2^{8}$ macrostates represents a finite number of
microstates ($2^{16}$), the upper bound on the number of distinct oscillations
which are macrostate invariant must decrease with the average length of those
oscillations. Likewise, rule 109 has only 31 macroexecutions, which are all
non/trivial and thus all end on cycles of period 2 or greater. Although this
mapping favors short/period oscillations, other mappings may allow for longer
cycles. For example, if the aggregate sum mapping in Fig. 2(a) is composed
with a function that maps a unity sum to rule 4, then PICARD will shift each
sum-1 row to the right by one position on each iteration. These shifts will
not alter the sum of the row, and thus this shifting will continue
indefinitely. Consequently, in this hypothetical case, the rule-4 region will
contain macroexecutions with periods as long as the length of each microstate
row. However, by the discussion above, rule 4 is unlikely to contain a high
number of macroexecutions due to the number of microstates consumed by these
long/period oscillations.

As shown above in Fig. 5(a), PICARD oscillations need not be contained within
a macroexecution – microstates can oscillate across rule boundaries.
Additionally, the transient components of PICARD executions can have rich
structure. In Fig. 9, there is a long sequence of transient execution leading
up to an eventual rule-0 fixed point.

Figure 9: Fixed point with long structured transient. A rule-0 fixed point is
shown that is preceded by a very long transient that appears to have
significant visual structure.

The CA appears to have significant structure visually, but the rule evolution
is less predictable and shows clear path dependence. Thus, the feedback in
PICARD both constrains its executions and generates novel patterns that may be
fodder for further analysis.

## Summary and Future Work

By feeding a conventional one/dimensional elementary cellular automaton’s own
state back into the iteration rule that generates its next state, we have
developed a new modeling framework for investigating information control in
natural systems. Although open/loop cellular automata that evolve under
strictly static rules have been useful tools for studying evolution
artificially, these classic frameworks are not sufficiently rich to model the
dependence of fitness on current state. Moreover, by adapting the methods of
evolving cellular automata so that PICARD mappings can be the targets of
artificial selection, there is a potentially richer set of behaviors that can
be selected for.

We have shown how state feedback can generate self/reinforcing regions of
behavior. Thus, this framework provides a model of how functional diversity
can be embedded within a single automaton. This functional diversity is
similar to the diverse differentiation possibilities for cells in living
organisms. Moreover, as evidenced by the distribution of macrostates over
macroexecutions, this framework apparently shares similarities with the
processes responsible for the generation of natural language.

A possible criticism of PICARD is that it introduces non/local effects to
cellular automata. Traditional CA are built from the assumptions that cells
update based entirely on local information and information about neighbors. By
feeding the state of an entire row back into the CA rule, this locality
assumption is broken. A complete reductionist model of the dynamics of life
would necessarily have to incorporate dynamics of both the focal entity and
the environment around it – despite any significant differences in time
scales. By feeding the state of the focal entity back into the rules that
govern how the entity evolves, we avoid these complications and still allow
for entities to alter _and be altered by_ their environment. Thus, just as
coarse/grained descriptions have utility in modeling physical phenomena,
non/local effects in PICARD mappings can be viewed as useful coarse/grainings
of otherwise less tractable models.

As PICARD Implements CA Rules Differently, there are myriad directions for
future work – including directions which parallel investigations from
conventional cellular automata and directions which are specific to PICARD
self/reference. At a minimum, the distribution of macrostates over
macroexecutions needs to be investigated for a wide variety of PICARD mappings
to determine whether the Zipfian distribution is widespread. Although it is
attractive to view the PICARD automata in this paper as a mosaic of different
elementary CA that each govern small patches of local elementarity, it is
unlikely that such regions exist that contain executions with the relatively
open/ended complexity of open/loop elementary CA patterns. However, PICARD
mappings have the potential of generating new patterns and new functionality
that may be out of reach of conventional CA’s. Alternatively, PICARD mappings
may have the ability to increase the robustness of elementary CA function by
reducing the sensitivity to minute but functionally inconsequential changes in
initial condition. Extending PICARD mappings to be able to map from prior
microstates may allow for incorporating functionality normally associated with
higher/dimensional cellular automata. For example, Langton,’s loops (Langton,,
1984) are able to self replicate because growing patterns in space can turn
and interact with earlier patterns. A PICARD mapping can similarly connect
iterations across time and draw a connection between two/dimensional self
replication and one/dimensional pattern generation. In general, an important
future direction is to connect one/dimensional PICARD insights to higher
dimensional cellular automata.

## Acknowledgements

Kunihiko Kaneko and Larissa Albantakis provided useful feedback on earlier
versions of this work. Interaction with them was possible due to a workshop on
Information, Complexity, and Life organized by the BEYOND Center for
Fundamental Concepts in Science at Arizona State University. Andrea Richa also
provided helpful comments regarding the presentation of this material. This
project/publication was made possible through support of a grant from
Templeton World Charity Foundation. The opinions expressed in this publication
are those of the author(s) and do not necessarily reflect the views of
Templeton World Charity Foundation.

## References

  * Breukelaar and Bäck, (2005) Breukelaar, R. and Bäck, T. (2005).  Using a genetic algorithm to evolve behavior in multi dimensional cellular automata: emergence of behavior.  In Proceedings of the 7th Annual Conference on Genetic and Evolutionary Computation (GECCO ’05), pages 107–114, Washington, DC, USA.
  * Das et al., (1994) Das, R., Mitchell, M., and Crutchfield, J. P. (1994).  A genetic algorithm discovers particle computation in cellular automata.  In Proceedings of the Third Conference on Parallel Problem Solving from Nature, volume 866 of Lecture Notes in Computer Science, pages 344–353. Springer.
  * Davidson, (2010) Davidson, E. H. (2010).  The Regulatory Genome: Gene Regulatory Networks In Development And Evolution.  Academic Press.
  * Davies, (2012) Davies, P. C. W. (2012).  The epigenome and top-down causation.  Interface Focus., 2(1):42–48.
  * Goldenfeld and Woese, (2011) Goldenfeld, N. and Woese, C. (2011).  Life is physics: evolution as a collective phenomenon far from equilibrium.  Annu. Rev. Condens. Matter Phys., 2(1):375–399.
  * Hofstadter, (1979) Hofstadter, D. R. (1979).  Gödel, Escher, Bach: an Eternal Golden Braid.  Basic Books, Inc.
  * Hordijk, (2013) Hordijk, W. (2013).  The EvCA project: a brief history.  Complexity, 18(5):15–19.
  * (8) Kataoka, N. and Kaneko, K. (2000a).  Functional dynamics. I: articulation processes.  Phys. D: Nonlinear Phenom., 138(3–4):225–250.
  * (9) Kataoka, N. and Kaneko, K. (2000b).  Natural language from function dynamics.  Biosystems, 57(1):1–11.
  * Langton, (1984) Langton, C. G. (1984).  Self-reproduction in cellular automata.  Phys. D: Nonlinear Phenom., 10(1–2):135–144.
  * Mitchell et al., (1994) Mitchell, M., Crutchfield, J. P., and Hraber, P. T. (1994).  Evolving cellular automata to perform computations: mechanisms and impediments.  Phys. D: Nonlinear Phenom., 75(1–3):361–391.
  * Schlitt and Brazma, (2007) Schlitt, T. and Brazma, A. (2007).  Current approaches to gene regulatory network modelling.  BMC Bioinformatics, 8(6).
  * Walker et al., (2012) Walker, S. I., Cisneros, L., and Davies, P. C. W. (2012).  Evolutionary transitions and top-down causation.  In Proceedings of the Thirteenth International Conference on the Simulation and Synthesis of Living Systems (Artificial Life 13), volume 13, pages 283–290, East Lansing, Michigan.
  * Walker and Davies, (2013) Walker, S. I. and Davies, P. C. W. (2013).  The algorithmic origins of life.  J. R. Soc. Interface, 10(79).
  * Wolfram, (2002) Wolfram, S. (2002).  A New Kind of Science.  Wolfram Media.
  * Zipf, (1949) Zipf, G. K. (1949).  Human Behavior and the Principle of Least Effort.  Addison-Wesley Press, Oxford, England.


Also at ]Department of Physics & Department of Chemistry, Ulsan National
Institute of Science and Technology, Ulsan, 44919, Republic of Korea

date22042021

# Positional information as a universal predictor of freezing

Tamoghna Das<EMAIL_ADDRESS>Center for Soft and Living Matter,
Institute for Basic Science (IBS), Ulsan,
44919, Republic of Korea  Tsvi Tlusty<EMAIL_ADDRESS>[

(date)

###### Abstract

Variation of positional information, measured by the two-body excess entropy
$\mathsf{S}_{\mathrm{2}}$, is studied across the liquid-solid equilibrium
transition in a simple two-dimensional system. Analysis reveals a master
relation between $\mathsf{S}_{\mathrm{2}}$ and the freezing temperature
$T_{\mathrm{f}}$, from which a scaling law is extracted:
$-\mathsf{S}_{\mathrm{2}}\sim|T_{\mathrm{f}}-T|^{-1/3}$. Theoretical and
practical implications of the observed universality are discussed.

Excess entropy, Phase transition, Liquid State Theory

## Introduction

Locating thermodynamic transition lines is an essential part of material
characterisation. However, broad scanning of the phase diagram of a system for
any sharp changes in its thermodynamic response, primarily, heat capacity, is
an exhaustive computationally-heavy procedure. As freezing of a liquid is a
global structural transition, we propose here that the positional information
alone can be a quick and accurate indicator of this phase transition. Such
positional information of particulate systems is often directly accessible via
simple X-ray scattering measurements, rendering it a practical means to
efficiently chart the phase diagram.

Liquids lack practically any apriori positional information, as all particles
are arranged in a random configuration, with only occasional short range
correlations. Thus, setting the totally random configuration as a reference
state of zero configurational entropy, $\mathsf{S}_{\mathrm{c}}=0$, a liquid
would have a very small negative $\mathsf{S}_{\mathrm{c}}$. During freezing,
as spatial correlation increases, information is gained leading to further
decrease in $\mathsf{S}_{\mathrm{c}}$. Finally, after freezing, all particles
settle down in a specific crystalline symmetry described by a certain finite
amount of information, corresponding to a large negative value of
$\mathsf{S}_{\mathrm{c}}$.

In general, a disordered yet correlated arrangement of $N$ particles can be
envisioned as a collection of $n$-tuples where each tuple is a possible
representation of local $n$-body correlation occurring in the system with
frequency $\rho^{n}$. The total configurational entropy
$\mathsf{S}_{\mathrm{c}}$ can then be computed as an integrated contribution
of such correlations as [1, 2, 3]
$\mathsf{S}_{\mathrm{c}}/N=\mathsf{S}_{\mathrm{1}}+\sum_{n=2}^{N}\mathsf{S}_{\mathrm{n}}$,
where $\mathsf{S}_{\mathrm{1}}$ is the contribution of a system of non-
interacting particles, namely, a classical ideal gas (up to an additive
constant) and

$\mathsf{S}_{\mathrm{n}}=-\frac{1}{n!}~{}\rho^{n-1}\int
d\mathbf{r}^{2}g^{\mathrm{[n]}}\ln\left[\frac{g^{\mathrm{[n]}}}{\mathsf{H}}\right]~{}.$
(1)

Here, $g^{\mathrm{[n]}}$ is the $n$-body correlation and $\mathsf{H}$
represents the combination of all possible $(n-1)$-tuples.

For a simple liquid, large $n$-tuples can be reasonably approximated as the
product of smaller $n$-tuples due to the absence of any long-lived local
structures. [4] Under such approximation,

$\mathsf{H}\approx
g^{\mathrm{[3]}}(1,2,3)=g^{\mathrm{[2]}}(1,2)g^{\mathrm{[2]}}(2,3)g^{\mathrm{[2]}}(3,1)$

and $g^{\mathrm{[2]}}$ becomes the sike relevant positional correlation. The
configurational entropy $\mathsf{S}_{\mathrm{c}}$, expressed in terms of only
two-body correlation, now reads as: [5]

$\mathsf{S}_{\mathrm{2}}={\textstyle\frac{1}{2}}+{\textstyle\frac{1}{2}}\rho\int
d{\mathbf{r}}~{}\left(g^{\mathrm{[2]}}-1\right)-{\textstyle\frac{1}{2}}\rho\int
d{\mathbf{r}}~{}\left(g^{\mathrm{[2]}}\ln g^{\mathrm{[2]}}\right)~{}.$ (2)

The positional information $-\mathsf{S}_{\mathrm{2}}$ is the excess
information over the random configuration due to pairwise interaction in the
system. Note that the upper bound of $\mathsf{S}_{\mathrm{2}}$ is
${\textstyle\frac{1}{2}}$ and not zero.

$\mathsf{S}_{\mathrm{c}}$, or its approximant $\mathsf{S}_{\mathrm{2}}$,
appears to be a natural candidate to relate properties of a material to the
statistical nature of its underlying internal structure. Such measures
physically approximate the logarithm of the phase space volume available to
the systems. [6] As in a classical Enskog picture, one may relate this
positional information with the timescale of a particle’s motion within the
cage formed by its neighbours. Stemming from this idea, scaling laws [7, 8, 9]
relating $\mathsf{S}_{\mathrm{2}}$ and transport coefficients of materials
have been proposed, which have been observed to certain degrees in past and
recent experiments [10, 11, 12, 13] and simulation studies. [14, 15, 16, 17,
18, 19] In the context of glassy systems, $\mathsf{S}_{\mathrm{c}}$ is linked
the non-Arhenius relaxation processes, as formalised in the Adam-Gibbs theory.
Specifically, the relation between $\mathsf{S}_{\mathrm{2}}$ and glassy
relaxation times was examined. [20, 21] Recently, $\mathsf{S}_{\mathrm{2}}$
has been recast to account for the angular correlation [22] in particulate
systems, and a local measure [23] was formulated to study the microscopic
features of materials’ structure. $\mathsf{S}_{\mathrm{2}}$ has even been used
to identify structural motifs of proteins. [24] Still, one obvious case of
structural transformation, freezing, has only attracted rare attention [25,
26] in this regard.

In this study, we explore the possibility of using $\mathsf{S}_{\mathrm{2}}$
as an indicator of freezing transition in a simple liquid. We recognise that
the two integrands contributing to Eq.(2) originate from two independent
physical processes. [19] The first integrand originates from fluctuations in
the number of particles in the grand-canonical ensemble. This positive-
definite contribution, $\Delta\mathsf{S}_{\mathrm{fluc}}$, termed as number
fluctuation entropy, is related to the system’s isothermal compressibility
$\kappa_{\mathrm{T}}$,

$\Delta\mathsf{S}_{\mathrm{fluc}}\equiv{\textstyle\frac{1}{2}}+{\textstyle\frac{1}{2}}\rho\int
d{\mathbf{r}}~{}\left(g^{\mathrm{[2]}}-1\right)={\textstyle\frac{1}{2}}\rho\mathsf{k}_{\mathrm{B}}\kappa_{\mathrm{T}}~{},$
(3)

where $\mathsf{k}_{\mathrm{B}}$ is the Boltzmann constant. The second
integrand is a negative-definite quantity termed loss entropy as it captures
the entropy lost due to pair-wise correlation that is absent from the
reference ideal gas,

$\Delta\mathsf{S}_{\mathrm{loss}}\equiv-{\textstyle\frac{1}{2}}\rho\int
d{\mathbf{r}}~{}\left(g^{\mathrm{[2]}}\ln g^{\mathrm{[2]}}\right)~{}.$ (4)

These two entropy represent two opposing effects. On one hand, increasing
spatial correlation would shrink the phase space volume resulting in more
negative $\Delta\mathsf{S}_{\mathrm{loss}}$. On the other hand, this would
also make exchange of particles with the environment increasingly difficult,
i.e., $\Delta\mathsf{S}_{\mathrm{fluc}}$ becomes more positive. The two-body
excess entropy $\mathsf{S}_{\mathrm{2}}$ is governed by the interplay of these
two opposing contributions.

## Results

Figure 1:  (A) Pair correlation function, $g(r)$ for $\rho=0.9$ is plotted for
various temperatures $T$. At $T=1.0$, a bifurcation of the second peak
indicates the transition point between liquid ($T>1.0$) and solid ($T<1.0$)
states. Maximal correlation for each $T$ is marked by a $\star$. Inset:
position of this maximum, $r^{\star}$ as a function of $T$, where the
transition is marked by a dashed line $T=T^{\star}=1.0$. Two entropy terms,
(B) $\Delta\mathsf{S}_{\mathrm{loss}}$ and (C)
$\Delta\mathsf{S}_{\mathrm{fluc}}$, plotted against $T^{\star}/T$ show a data
collapse for systems in the liquid phase, $T>T^{\star}$, while the density
dependence is evident in the solid phase, $T<T^{\star}$.  Figure 2:  Master
relation between $\mathsf{S}_{\mathrm{2}}$ and the freezing temperature
$T_{\mathrm{f}}$. Two distinct branches of the collapsed data are fitted with
power laws $x^{-\nu}$ where $x=|\Delta T/T_{\mathrm{f}}|=|\tau-1|$. The solid
line, with slope $\nu\sim 1/3$, corresponds to the systems in the liquid phase
and the crystalline phase data are fitted with $\nu\sim-2$ is shown by (dashed
line).

To understand the temperature dependence of structural probe
$\mathsf{S}_{\mathrm{2}}$, we examined an ensemble of simulated equilibrium
configurations of a simple system of unit-mass point particles over a wide
range of densities and temperatures. The particles in this model system
interact pair-wise via a Lennard-Jones potential, [27]
$\phi(r)=4\epsilon\left[(\sigma/r)^{12}-(\sigma/r)^{6}\right]$ where
$r=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ is the distance between a pair $\\{i,j\\}$
of particles. The potential is commonly used model of simple non-polar fluids.
[28, 29, 30] The energy, length, and time scales are set by $\epsilon$,
$\sigma$ and $\tau=\mathchoice{{\hbox{$\displaystyle\sqrt[\
]{\sigma^{2}/\epsilon\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\textstyle\sqrt[\
]{\sigma^{2}/\epsilon\,}$}\lower 0.4pt\hbox{\vrule
height=7.5pt,depth=-6.00003pt}}}{{\hbox{$\scriptstyle\sqrt[\
]{\sigma^{2}/\epsilon\,}$}\lower 0.4pt\hbox{\vrule
height=5.25pt,depth=-4.20003pt}}}{{\hbox{$\scriptscriptstyle\sqrt[\
]{\sigma^{2}/\epsilon\,}$}\lower 0.4pt\hbox{\vrule
height=3.75pt,depth=-3.00002pt}}}\kern 1.00006pt$, respectively. We have
chosen to work in two-dimensions ($\mathsf{2D}$) for computational
convenience. Equilibrium configurations of $N=$60\,000$$ particles within a
periodically bounded box of size fixed by the density are simulated using
canonical molecular dynamics method. [31] The temperature $T$, measured in
units of $\epsilon$, is fixed by a Langevin thermostat throughout the
simulation (as implemented in LAMMPS [32]). We note that the choice of
thermostat does not affect the equilibrium configuration of the system.

Since $\mathsf{S}_{\mathrm{2}}$, as defined in Eq.(2), is an intensive
ensemble-invariant thermodynamic property, we may conveniently replace
$g^{\mathrm{[2]}}$ by its canonical counterpart, the pair-correlation
function, $g(r)=\rho^{-1}\left\langle\sum_{i}\delta(r-r_{i})\right\rangle$.
For a typical simple liquid, [33] $g(r)$ shows an oscillatory behaviour at
small $r$, comparable to a few particle diameter, pointing to short-range
ordering, which die out exponentially with increasing $r$, as $g(r)$ tends
asymptotically to unity. Upon lowering $T$ at a constant $\rho$, as
crystalline order emerges in the system, the second peak of $g(r)$ splits
symmetrically into two halves. [34] The geometric origin of this split is the
$\mathsf{2D}$ crystalline configuration: due to the hexagonal packing of
particles, one can always find two particles in the second-nearest neighbour
shell positioned at angles $2\pi/3$ and $\pi$ (with respect to a suitably
chosen axis).

As an example, we plot $g(r)$ of the model system (FIG.1A) for a set of $T$ at
$\rho=0.9$. A symmetric splitting of the second peak can be observed at
$T=1.0$. We also note that the position of the first peak, $r^{\star}$,
gradually shifts towards larger values as the system condenses upon lowering
$T$. (FIG.1A Inset) A rapid change in $r^{\star}$ is indicative of a phase
transition. The crossover temperature, marked by $T^{\star}$, coincides with
the split of the second peak. For $T<T^{\star}$, $g(r)$ is decorated with
peaks at specific positions dictated by the crystalline symmetry. As $T$
crosses $T^{\star}$, $g(r)$ exhibits more and more liquid-like features.
Following this observational cue, we consider $T^{\star}$ as a primary
estimate of the freezing temperature ($T_{\mathrm{f}}$) and identify it for
each of the densities. Given the $g(r)$ data, it is now straightforward to
compute $\Delta\mathsf{S}_{\mathrm{fluc}}$ and
$\Delta\mathsf{S}_{\mathrm{loss}}$ using Eqs. (3) and (4). In the liquid
state, $T>T^{\star}$, the density dependence of both entropy terms can be
mostly eliminated by re-scaling the temperature as $T^{\star}/T$ (FIG.1B &
1C). However, the density dependence remains evident in the solids,
$T<T^{\star}$.

Figure 3:  Density dependence of the freezing temperature, $T_{\mathrm{f}}$
(blue circles). The dashed line is the theoretical prediction for high-density
Lennard-Jones system in $\mathsf{2D}$ [35]. Horizontal lines represent the
only known data for the coexistence region in this system [36].

Strikingly, the $\mathsf{S}_{\mathrm{2}}$-data for all densities collapse to a
single master curve when the values of $T^{\star}$ are slightly adjusted to
$T_{\mathrm{f}}$. The maximal adjustment, less than $5\%$ of $T^{\star}$, was
required for the densest configuration ($\rho=1.0$). In FIG.2, we plot
$-\mathsf{S}_{\mathrm{2}}$ shifted by $0.5$ as a function of the reduced
temperature difference, $x=|\Delta T/T_{\mathrm{f}}|$ and find two distinct
branches for the liquid phase (values $<0.5$) and crystalline phase (values
$>0.5$). Both branches are fitted with power laws $x^{-\nu}$. A dashed line
with $\nu=-2.0$ fits the crystalline branch. For the liquid branch, we find
$\nu=0.326\pm 0.003$ as shown by a solid line.

The thermal dependence of $\mathsf{S}_{\mathrm{2}}$ closely resembles the
critical behavior of heat capacity of a three-state Potts model with
$\nu=1/3$. But unlike heat capacity, which is a thermodynamic response of the
system, $\mathsf{S}_{\mathrm{2}}$ is an intensive geometric quantity. In
recent years, increasing evidences suggests a physical view of liquids as a
dynamic multi-state systems, [37, 38, 39] in contrast to the mean field
description based on homogeneous local density. Specifically, a three-state
geometric picture [40] appears to capture the local dynamical aspects of
liquids better than the current predictions. While the present finding
strengthens this view, it also questions the traditional liquid-state theory
and calls for further detailed investigations.

The values of the freezing temperature $T_{\mathrm{f}}$, as plotted in FIG.3,
are found to reside within the known boundaries of the coexistence region
determined by Barker et al. [36] (horizontal lines). Interestingly, these
values closely follow a newly-derived equation of state of $\mathsf{2D}$
Lennard-Jones systems. [35] The dashed line in FIG.3 marks the temperature
above which the anharmonic effects at a specific $\rho$ become prominent, and
the freezing temperature is expected to be below this curve. Fulfilling this
expectation, $T_{\mathrm{f}}$ appears as an overall fair estimate of freezing
temperature of the present model system.

## Conclusion

The phenomenological scaling relation,
$-\mathsf{S}_{\mathrm{2}}\sim|T_{\mathrm{f}}-T|^{-1/3}$, proposes
$\mathsf{S}_{\mathrm{2}}$ as a reliable universal predictor of the freezing
transition, at least in the case of simple liquids. Higher-order correlations
would be necessary to construct similar scaling relations for more complex
network forming systems. [41, 42, 43, 44] We expect that in systems such as
carbon-based polymeric mixtures, nano-particle composites, aggregate forming
colloidal and other bio-inspired materials, $\mathsf{S}_{\mathrm{2}}$-based
relations would provide a primary guidance for transition lines.

The results of this study come with two implications. First, from a
theoretical point of view, it suggests the Corresponding States Principle [45,
46, 47, 48] — so far used to advocate the universality in liquids — might now
be extended to the crystalline states of matter. Thus, a unified description
of the states of matter might be drawn based on the system’s information
content. Recalling the numerous efforts to relate $\mathsf{S}_{\mathrm{2}}$
with the dynamical response in condensed materials, [49] one might also hope
to gain insight regarding the links between the dynamics and thermodynamics of
matter in terms of the information exchange within the system. Second, from a
practical perspective, the scaling relations allow to predict the transition
temperature from simple X-ray scattering measurements of positional
correlation in the system at any density and temperature, at least for a
simple liquid. Since identification of a phase transitions often requires an
exhaustive scan of the parameter space, this scaling relation would
effectively reduce that huge search space into a single measurement.

## References

  * Green [1952] H. S. Green,  _The Molecular Theory of Fluids_ (North Holland, Amsterdam, 1952).
  * Nettleton and Green [1958] R. E. Nettleton and M. S. Green, “Expression in terms of molecular distribution functions for the entropy density in an infinite system,” J. Chem. Phys. 29, 1365 (1958).
  * Wallace [1989] D. C. Wallace, “Statistical theory for the entropy of a liquid,” Phys. Rev. A 39, 4843 (1989).
  * Kirkwood [1934] J. G. Kirkwood, “Quantum statistics of almost classical assemblies,” Phys. Rev. 45, 116 (1934).
  * Baranyai and Evans [1989] A. Baranyai and D. J. Evans, “Direct entropy calculation from computer simulation of liquids,” Phys. Rev. A 40, 3817 (1989).
  * Chapman and Cowling [1990] S. Chapman and T. G. Cowling,  _The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_ (Cambridge University Press, 1990).
  * Rosenfeld [1977] Y. Rosenfeld, “Relation between the transport coefficients and the internalentropy of simple systems,” Phys. Rev. A 15, 2545 (1977).
  * Rosenfeld [1999] Y. Rosenfeld, “A quasi-universal scaling law for atomic transport in simple fluids,” J. Phys.: Condens. Matter 11, 5415 (1999).
  * Dzugutov [1996] M. Dzugutov, “A universal scaling law for atomic diffusion in condensed matter,” Nature 381, 137 (1996).
  * Vaz  _et al._ [2012] R. V. Vaz, A. L. Magalhaes, D. L. Fernandes, and C. M. Silva, “Vaz rv, magalhães al, fernandes dl, silva cm. universal correlation of self-diffusion coefficients of model and real fluids based on residual entropy scaling law,” Chem. Eng. Sci. 79, 153 (2012).
  * Ma  _et al._ [2013] X. Ma, W. Chen, Z. Wang, Y. Peng, Y. Han, and P. Tong, “Test of the universal scaling law of diffusion in colloidal monolayers,” Phys. Rev. Lett. 110, 078302 (2013).
  * Hansen  _et al._ [2018] H. W. Hansen, A. Sanz, K. Adrjanowicz, B. Frick, and K. Niss, “Evidence of a one-dimensional thermodynamic phase diagram for simple glass-formers,” Nat. commun. 9, 1 (2018).
  * Galloway  _et al._ [2020] K. L. Galloway, X. Ma, N. C. Keim, D. J. Jerolmack, A. G. Yodh, and P. E. Arratia, “Scaling of relaxation and excess entropy in plastically deformed amorphous solids,” Proc Natl Acad Sci USA 117, 11887 (2020).
  * Hoyt, Asta, and Sadigh [2000] J. J. Hoyt, M. Asta, and B. Sadigh, “Test of the universal scaling law for the diffusion coefficient in liquid metals,” Phys. Rev. Lett. 85, 594 (2000).
  * Chakraborty and Chakravarty [2006] S. N. Chakraborty and C. Chakravarty, “Diffusivity, excess entropy, and the potential-energy landscape of monatomic liquids,” J. Chem. Phys. 124, 014507. (2006).
  * Mittal, Errington, and Truskett [2007] J. Mittal, J. R. Errington, and T. M. Truskett, “Relationships between self-diffusivity, pack-ing fraction, and excess entropy in simple bulk and confined fluids,” J Phys Chem B 111, 10054 ((2007)).
  * Agarwal  _et al._ [2011] M. Agarwal, M. Singh, R. Sharma, M. P. Alam, and C. Chakravarty, “Relationship between structure, entropy, and diffusivity in water and water-like liquids,” J. Phys. Chem. B 114, 6995 (2011).
  * Joy [2017] A. Joy, “Universal scaling of pair-excess entropy and diffusion in yukawa liquids,” Phys. Plasmas 24, 010702 (2017).
  * Gao and Widom [2018] M. C. Gao and M. Widom, “Information entropy of liquid metals,” J. Phys. Chem. B 122, 3550 (2018).
  * Gallo and Rovere [2015] P. Gallo and M. Rovere, “Relation between the two-body entropy and the relaxation time in supercooled water,” Phys. Rev. E 91, 012107 (2015).
  * Banerjee  _et al._ [2017] A. Banerjee, M. K. Nandi, S. Sastry, and S. M. Bhattacharya, “Determination of onset temperature from the entropy for fragile to strong liquids,” J. Chem. Phys. 147, 024504 (2017).
  * Ingebrigtsen and Tanaka [2018] T. S. Ingebrigtsen and H. Tanaka, “Structural predictor for nonlinear sheared dynamics in simple glass-forming liquids,” Proc Natl Acad Sci USA 115, 87 (2018).
  * Piaggi and Parinello [2017] P. M. Piaggi and M. Parinello, “Entropy based fingerprint for local crystalline order,” J. Chem. Phys. 142, 114112 (2017).
  * Young  _et al._ [2007] T. Young, R. Abel, B. Kim, B. J. Berne, and R. A. Friesner, “Motifs for molecular recognition exploiting hydrophobic enclosure in protein–ligand binding,” Proc Natl Acad Sci USA 104, 808 (2007).
  * Giaquinta, Giunta, and Giarritta [1992] P. V. Giaquinta, G. Giunta, and S. P. Giarritta, “Entropy and the freezing of simple liquids,” Phys. Rev. A 45, R6966 (1992).
  * Rosenfeld [2000] Y. Rosenfeld, “Excess-entropy and freezing-temperature scalings for transport coefficients: Self-diffusion in yukawa systems,” Phys. Rev. E 62, 7524 (2000).
  * Lennard-Jones [1924] J. E. Lennard-Jones, “On the determination of molecular fields. ii. from the equation of state of a gas,” Proc. R. Soc. Lond. A 106, 463 (1924).
  * Vliegenthart and Lekkerkerker [2000] G. A. Vliegenthart and H. N. W. Lekkerkerker, “Predicting the gas–liquid critical point from the second virial coefficient predicting the gas–liquid critical point from the second virial coefficient,” J. Chem. Phys. 112, 5364 (2000).
  * Dunikov, Malyshenko, and Zhakhovski [2001] D. O. Dunikov, S. P. Malyshenko, and V. V. Zhakhovski, “Corresponding states law and molecular dynamics simulations of the lennard-jones fluid,” J. Chem. Phys. 115, 6623 (2001).
  * Charpentier and Jakse [2005] I. Charpentier and N. Jakse, “Phase diagram of complex fluids using an efficient integral equation method,” J. Chem. Phys. 123, 204910 (2005).
  * Frenkel and Smit [1996] D. Frenkel and B. Smit,  _Understanding Molecular Simulation_ (Academic Press, 1996).
  * Plimpton [1995] S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comp. Phys. 117, 1 (1995).
  * Hansen and McDonald [2006] J.-P. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Academic Press, 2006).
  * Truskett  _et al._ [1998] T. M. Truskett, S. Torquato, S. Sastry, P. G. Debenedetti, and F. H. Stillinger, “Structural precursor to freezing in the hard-disk and hard-sphere systems,” Phys. Rev. E 58, 3083 (1998).
  * Shi  _et al._ [2018] K. Shi, K. Gu, Y. Shen, D. Srivastava, E. E. Santiso, and K. E. Gubbins, “High-density equation of state for a two-dimensional lennard-jones solid,” J. Chem. Phys. 148, 174505 (2018).
  * Barker, Henderson, and Abraham [1981] J. A. Barker, D. Henderson, and F. F. Abraham, “Phase diagram of the two-dimensional lennard-jones system; evidence for first-order transitions,” Physica A 106, 226 (1981).
  * Hoang, Teboul, and Odagaki [2015] V. V. Hoang, V. Teboul, and T. Odagaki, “New scenario of dynamical heterogeneity in supercooled liquidand glassy states of 2d monatomic system,” J. Phys. Chem. B 119 (2015).
  * Das and Douglas [2018] T. Das and J. F. Douglas, “Quantifying structural dynamic heterogeneity in a dense two-dimensional equilibrium liquid,” J. Chem. Phys. 149, 144504 (2018).
  * Krebs  _et al._ [2018] Z. Krebs, A. B. Roitman, L. M. Nowack, C. Liepold, B. Lin, and S. A. Rice, “Transient structured fluctuations in a two-dimensional system with multiple ordered phases,” J. Chem. Phys. 149, 034503 (2018).
  * Das and Douglas [2019] T. Das and J. F. Douglas, “Three-state heterogeneity in a model two-dimensional equilibrium liquid,” J. Mol. Liq. 293, 111466 (2019).
  * Baranyai and Evans [1990] A. Baranyai and D. J. Evans, “Three-particle contribution to the configurational entropy of simple fluids,” Phys. Rev. A 42, 849 (1990).
  * Molinero and Moore [2009] V. Molinero and E. B. Moore, “Water modeled as an intermediate element between carbon and silicon,” J. Phys. Chem. B 113, 4008 (2009).
  * Cao  _et al._ [2011] Q.-L. Cao, X.-S. Kong, Y. D. Li, X. Wu, and C. S. Liu, “Revisiting scaling laws for the diffusion coefficients in simple melts based on the structural deviation from hard-sphere-like case,” Physica B 406, 3114 (2011).
  * Dhabal  _et al._ [2015] D. Dhabal, A. H. Nguyen, M. Singh, P. Khatua, V. Molinero, S. Bandyopadhyay, and C. Chakravarty, “Excess entropy and crystallization in stillinger-weber and lennard-jones fluids,” J. Chem. Phys. 143, 164512 (2015).
  * Pitzer [1939] K. S. Pitzer, “Corresponding states for perfect liquids,” J. Chem. Phys. 7, 583 (1939).
  * Pitzer [1955] K. S. Pitzer, “The volumetric and thermodynamic properties of fluids. i. theoretical basis and virial coefficients,” J. Am. Chem. Soc. 77, 3427 (1955).
  * Pitzer  _et al._ [1955] K. S. Pitzer, D. Z. Lippmann, J. R. F. Curl, C. M. Huggins, and D. E. Petersen, “The volumetric and thermodynamic properties of fluids. ii. compressibility factor, vapor pressure and entropy of vaporization,” J. Am. Chem. Soc. 77, 3433 (1955).
  * Noro and Frenkel [2000] M. G. Noro and D. Frenkel, “Extended corresponding-states behavior for particles with variable range attractions,” J. Chem. Phys. 113, 2941 (2000).
  * Dyre [2018] J. C. Dyre, “Perspective: Excess-entropy scaling,” J. Chem. Phys. 149, 210901 (2018).