# NeuroTac: A Neuromorphic Optical Tactile Sensor
applied to Texture Recognition

Benjamin Ward-Cherrier, Member, IEEE, Nicholas Pestell, Student Member, IEEE,
Nathan F. Lepora, Member, IEEE BWC was supported by a University of Bristol
Vice-Chancellor’s fellowship, NP was supported by an EPSRC DTP studentship and
NL was supported in part by a Leverhulme Trust Research Leadership Award on ’A
biomimetic forebrain for robot touch’ (RL-2016-039).Authors are with the
Department of Engineering Mathematics, University of Bristol and Bristol
Robotics Laboratory, University of Bristol, UK.
Email: {b.ward-cherrier, n.pestell<EMAIL_ADDRESS>
###### Abstract

Developing artificial tactile sensing capabilities that rival human touch is a
long-term goal in robotics and prosthetics. Gradually more elaborate
biomimetic tactile sensors are being developed and applied to grasping and
manipulation tasks to help achieve this goal. Here we present the neuroTac, a
novel neuromorphic optical tactile sensor. The neuroTac combines the
biomimetic hardware design from the TacTip sensor which mimicks the layered
papillae structure of human glabrous skin, with an event-based camera
(DAVIS240, iniVation) and algorithms which transduce contact information in
the form of spike trains. The performance of the sensor is evaluated on a
texture classification task, with four spike coding methods being implemented
and compared: Intensive, Spatial, Temporal and Spatiotemporal. We found
timing-based coding methods performed with the highest accuracy over both
artificial and natural textures. The spike-based output of the neuroTac could
enable the development of biomimetic tactile perception algorithms in robotics
as well as non-invasive and invasive haptic feedback methods in prosthetics.

##  I INTRODUCTION

The long-term scientific goal of human-like artificial touch is being worked
towards with the creation of gradually more elaborate biomimetic tactile
sensors. The development of sensors which mimic aspects of biological touch
could lead to safer robots and more intuitive and versatile prosthetic
devices. An example of such sensors are neuromorphic tactile sensors, which
aim to replicate the spike-based representation of information found in the
nervous system.

The principal objective of neuromorphic engineering is to develop technologies
which can exploit efficient representations of information and operate in a
rapid, power-efficient manner. Using neuromorphic technologies with event-
based outputs also holds potential for integration with the human nervous
system, as has been demonstrated with artificial retinas [1]. Optical tactile
sensors have also recently demonstrated progress on a number of tactile tasks
[2, 3], and have the advantage of capitalizing on advances in image
recognition and machine learning techniques. These technologies will be
combined here to develop a tactile sensor for robotic manipulation and
prosthetics applications.

The aims of this paper are as follows:

  * •

Develop a novel neuromorphic optical tactile sensor to enable the development
and investigation of biomimetic spike-based information processing methods.

  * •

Validate this sensor on a texture classification task.

  * •

Investigate 4 spike coding methods (intensive, spatial, temporal and
spatiotemporal) and their effect on texture recognition performance.

Figure 1: Transduction, encoding and decoding mechanisms for the neuroTac
sensor. The sensor mimics biological processes by accumulating pixel events
(Potentials) from an event-based camera and combining them into taxel events
(Spikes).

The neuroTac sensor, described here, follows the tradition of neuromorphic
technologies [4] in seeking to produce and decode spike-based information
(Fig. 1). The spike-based sensor output is coded using 4 different methods:
intensive (overall number of spikes), spatial (number of spikes per taxel),
temporal (number of spikes per time window) and spatiotemporal (Van Rossum
metric [5]). Although neuromorphic devices may present certain advantages over
their non-neuromorphic counterparts (speed and energy efficiency), our
principal objective here is linked to biomimetism. The TacTip sensor emulates
the internal structure of human skin [6], and the neuroTac adds to that
biomimetic morphology by producing a neuromorphic, spike-based output. We aim
to develop the neuroTac and associated spike-based information processing
methods to investigate the advantages this approach might provide biological
organisms.

The performance of the sensor and its associated coding methods is validated
on texture recognition tasks, in which artificial and natural textures are
identified using a K-nearest-neighbour (KNN) algorithm. The sensor
successfully performs texture recognition, with temporal methods producing the
highest classification accuracy. This underlines the importance of spike
timing in texture recognition tasks.

##  II BACKGROUND AND RELATED WORK

One of the strongest motivations for studying the human sense of touch is its
important role in our interactions and in the manipulation of our environment.
Artificial tactile sensing thus often mimics particular features of biological
touch through biomimetic hardware or perception algorithms [7]. An emerging
area of biomimetic tactile sensing aims to replicate biological spike-based
signalling through event-based systems and study the encoding of tactile
information [8]. This is generally referred to as neuromorphic sensing [4].

Neuromorphic sensing arose from seminal work on silicon neurons [9] and in the
area of vision it has led to the successful integration of artificial retinas
with the human nervous system [1]. Evidence suggests that a neuromorphic
tactile sensor could similarly be used to restore a natural sense of touch to
amputees [10, 11].

Large-scale event-based tactile sensors have been developed for use as robotic
skins [12], with a focus on the efficient processing of large quantities of
asynchronous data. Bartolozzi et. al. also developed an architecture for use
with distributed off-the shelf tactile sensors, transforming raw data into
event-based signals [13]. Another example of a large-scale system for the
investigation of spiking outputs in neuromorphic touch is that developed by
Lee et. al. [14]. These systems represent crucial technological developments
in the area of event-based systems, and will be essential in the creation of
rapidly reacting fully tactile robot systems.

Here, our focus is not on large-scale event-based robotic skins or platforms
for investigating spike-based encoding but rather on an artificial fingertip
with high spatial resolution for fine in-hand manipulation tasks. Oddo et. al.
developed a tactile sensor more in line with this design objective [15]. Their
neuromorphic tactile fingertip comprises 16 taxels combining raw outputs from
4 microelectromechanical system (MEMS) sensors, fed to an Izhikevich model to
produce spikes. This system simulates the outputs from biological SA1
afferents, and has been proven capable of accurately distinguishing natural
textures [16]. A crucial difference with the sensor presented here is that the
neuroTac emulates the fast-adapting responses of FA1 afferents to dynamic
contact (rather than the slowly adapting SA1 afferent responses). Recently, a
sensor has been developed using similar hardware, with an event-based camera
capturing the deformations of a deformable membrane [17]. The system was
calibrated for use as a force sensor and to measure material hardness. The
neuroTac sensor design aligns more closely both with traditional tactile
sensor designs comprising taxels (internal markers which transduce contact),
and biological fingertip structures (with taxels functioning as artificial
mechanoreceptors).

We validate the neuroTac on a texture classification task, in which 3d-printed
and natural textures are classified by sliding the sensor horizontally across
each stimulus. Texture identification is a common task in tactile sensing due
to its important implications in object recognition and manipulation. Studies
on the tactile sensing of textures generally involve naturally occurring
textures [18, 16, 19, 20]. Here, we initially wish to investigate the sensor’s
response to highly structured textures. To achieve this, we utilize purposely
designed 3d-printed stimuli with regular grids of cylindrical bumps as has
been done in past studies of biological touch [21]. Following this, we
evaluate the sensor’s classification performance on a set of 20 natural
textures.

##  III METHODS

###  III-A NeuroTac design and operation

Figure 2: The neuroTac sensor. The tip contains internal pins treated as
mechanoreceptors, which produce pixel events at the event-based camera. These
are pooled and converted to taxel events (akin to biological spikes) upstream.

##### Sensor design

The NeuroTac is based on the TacTip sensor [6], a 3d-printed optical tactile
sensor with a compliant dome-shaped outer membrane comprising biomimetic
internal markers which emulate the internal structure of human fingertips
[22]. To convert the TacTip into a neuromorphic sensor, we replace its camera
module (ELP, USBFHD01M-L21) with an event-based camera (iniVation, DAVIS240)
which processes only dynamic events within the image frame. The DAVIS240 is an
evolution of the DVS128 [23] and comprises 240x180 pixels, which process
changes in brightness through independent electronic circuits to produce
events in the address-event representation (AER) format [24]. These events are
combined by taxel and transmitted by the sensor, analogously to biological
spike trains (see Section III-A for more details on the sensor operation).

Like the TacTip, the NeuroTac is made up of 3 main hardware elements (Fig. 2):

  * •

Tip: This is a compliant, 3d-printed modular part whose outer membrane comes
into contact with the environment. Its internal surface comprises white-tipped
pins which are displaced during contact. It is filled with silicone gel
(Techsil, RTV27905) and covered with an acrylic lens to protect electronic
components.

  * •

LED ring: This illuminates the tip’s internal surface.

  * •

Camera: Housed within the main part of the sensor (the 3d-printed Body), this
camera is event-based (iniVation, DAVIS240) and produces AER outputs in
response to movements of the sensor’s internal markers.

Figure 3: Sensor operation. Pixel events produced by the camera (iniVation,
DAVIS240) are initially filtered, then pooled into a single taxel event.
Finally, the position of taxels is updated (see Section III-A). Figure 4:
Experimental setup: The NeuroTac is attached to a 6-dof industrial robot arm
(ABB, IRB120) which is used to slide the sensor horizontally across the
3d-printed textures.

##### Sensor operation

As the NeuroTac’s compliant membrane deforms through contact, its 49 white-
tipped internal pins deflect, and their movement triggers events in the camera
(iniVation, Davis240) in the Address-Event Representation (AER) format. These
events are produced by thresholding brightness changes at each photodiode in
parallel [23, 24], leading to fast data transmission and high temporal
precision. We designate these events ’pixel events’ as they are created at the
pixel level. Data transduction by the neuroTac then occurs through the 3
following steps (illustrated in Fig. 3):

  * •

Noise filtering: each of the sensor’s 49 pins represents a taxel, and has a
receptive field assigned to it (6 px diameter). Pixel events that occur
outside a taxel’s receptive field, or which do not have another pixel event
occur within a given spatiotemporal window (neighbouring pixels, 5 ms) are
filtered out.

  * •

Pooling: pixel events are pooled over a short duration (20 ms) and combined
into a ’taxel event’ based on the receptive field they are located within.
Each taxel event is an array comprising 3 numbers: the number of pixel events
it contains, their average location and their average timing. Note that each
taxel event is interpreted as a spike associated with an artificial
mechanoreceptor.

  * •

Position update: Receptive fields are re-centred around each taxel by shifting
towards detected pixel events, to account for each pin’s movement across the
image.

###  III-B Experimental setup

The first experiment with artificial textures involves sliding the NeuroTac
across 11 3d-printed textures. Artificial textures consist of rectangular
grids of cylindrical bumps (1 mm height) with equal spacing and diameter. The
texture grid size varies from 0 mm (smooth surface) to 5 mm in steps of 0.5 mm
(Fig. 4).

Figure 5: Natural textures used for classification. 20 textures were used,
with a full list presented in Table II.

Data is collected by mounting the NeuroTac on a 6-dof robotic arm (ABB,
IRB120) and sliding it horizontally across the textures (Fig. 4). The robot
slides the NeuroTac across each texture at a speed of 15 mm/s, over a distance
of 60 mm, comprising one data sample. We collect 100 samples for each texture,
to obtain a dataset of 100 (number of runs) $\times$ 11 (number of textures)
$\times$ 49 (number of taxels) spike trains.

In the second experiment on natural texture classification, the data
collection procedure is repeated for 50 runs over 20 natural textures(see
Table II).

###  III-C Spike train encoding methods

The multi-taxel spike trains obtained from the sensor are denoted as a matrix
of spike times $t_{n}^{i}$, where $n=1,2...,N$ (N denotes the number of
taxels, $N=49$ in this case) and $i=1...I_{n}$ ($I_{n}$ denotes the number of
spikes for taxel $n$). A sample is denoted as the multi-taxel spike train
resulting from the sensor sliding 60 mm horizontally across one texture.
Following, we describe the 4 coding methods applied to transform the spike
trains $t_{n}^{i}$ to the encoded representation $R$.

####  III-C1 Intensive coding

This encoding consists of the average spike count per taxel for a given data
sample. We name it intensive, analogously to work on biological touch [25],
since it can be interpreted as an overall intensity of the sample signal with
an absence of any spatial or temporal resolution. The encoding produces a
single average spike count per sample, which is used for texture
classification.

$R=\frac{1}{N}\sum_{n=1}^{N}\sum_{i=1}^{I_{n}}t_{n}^{i}$ (1)

where $N$ is the total number of taxels and $I_{n}$ denotes the number of
spikes for taxel n.

####  III-C2 Spatial coding

Here we consider the topology of the sensor, and sum the spike counts
separately for each of the sensor’s 49 taxels. The resulting array of 49 spike
rates is used as an encoded representation of the texture

$R_{n}=\sum_{i=1}^{I_{n}}t_{n}^{i}$ (2)

####  III-C3 Temporal coding

Temporal coding considers a rolling window of width $\Delta t$ over each data
sample, which is rolled forward over the time domain in timesteps of 1 ms.
Within the window $\Delta t$, the average spikes per taxel are recorded before
proceeding to the next timestep.

$R_{n}(t)=\frac{1}{N}\sum_{t=t}^{t+\Delta t}t_{n}^{i}(t)$ (3)

The $\Delta t$ parameter in this encoding method affects the classification
accuracy, therefore we optimize it through brute-force optimization within a
limited range of 1-200 ms (see Section IV-A).

####  III-C4 Spatiotemporal coding

Spatiotemporal coding uses the spatial and temporal features of spike trains
produced by the sensor. A multi-neuron Van Rossum distance [26] is used as a
metric in the texture classification, to ensure that both spatial and temporal
dimensions of the data are considered.

The spatiotemporal encoded representation can be considered the convolution of
the multi-taxel spike train with an exponential kernel used in the Van Rossum
distance calculation (see Section III-D).

$R_{n}^{i}(t)=t_{n}^{i}h(t-t_{i})$ (4)

where

$h=\begin{cases}0,&t\leq 0\\\ \frac{1}{\tau}e^{-t/\tau},&t\geq 0\end{cases}$
(5)

$\tau$ is a time constant parameter which can be optimized to create the most
accurate clustering of texture representations.

###  III-D Classification

Texture classification is performed through a KNN algorithm (k=4), which
assigns a class (texture grid size 0-5 mm in 0.5 mm steps here) to a test
sample based on its 4 closest training samples. For the intensive, spatial and
temporal coding methods, a standard Euclidean distance is used to calculate
distances between samples.

For the spatiotemporal coding, we use a multi-neuron Van Rossum distance [26]
which involves the convolution of spike trains with an exponential kernel
(Section III-C4) before applying a distance metric.

In the original Van Rossum metric [5], which applies to single neuron spike
trains, the distance metric between two spike trains is simply:

$d_{2}(t_{1},t_{2};h)=\sqrt{\int dt(f_{1}-f_{2})^{2}}$ (6)

where $f_{1}$ and $f_{2}$ are the convolved functions.

The extension of the Van Rossum metric to multi-neuron cases, as described by
Houghton and Sen [26], introduces an additional parameter $\theta$ which
represents the correlation between neurons. The parameter’s effect is best
described by taking its cosine, which varies from 0 to 1 with:

  * •

$cos\theta=0$. Labelled line code: each neuron is considered independent and
their distances are summed.

  * •

$cos\theta=1$. Summed population code: spike trains are superimposed before
calculating the distance metric.

  * •

$1>cos\theta>0$: an intermediate level of inter-neuron correlation between
these two extremes.

The $\theta$ angle is an optimizable parameter, and thus we perform a two
parameter Bayesian optimization of $\tau$ (exponential kernel time constant)
and $\theta$ for spatiotemporal encoding as described in section IV-B.

##  IV RESULTS

###  IV-A Inspection of data - Artificial textures

Figure 6: Examples of the spike trains produced by the NeuroTac when sliding
across three different 3d-printed textures. From left to right, textures are 0
mm grid (smooth), 2.5 mm grid and 5 mm grid, with the corresponding spike
trains, spatial and temporal distributions displayed below them.

Data is gathered by sliding the sensor horizontally across the set of 11
artificial textures (see Section III-B). The output of the NeuroTac consists
of 49 spike trains which represent the events being produced at each taxel.
The spike trains will vary in overall intensity for different textures, as
well as having distinct spatial and temporal signatures.

Examples of spike trains obtained for 3 textures of grid sizes 0 mm (smooth),
2.5 mm and 5 mm are displayed here (Fig. 6, second row). It is visually
noticeable that these 3 multi-neuron spike trainsdiffer, with the number of
spikes produced increasing with texture coarseness. As expected, applying
intensive coding to these samples gives readily distinguishable average spike
counts over all taxels of $4.63$ spikes/taxel (0 mm grid), $9.84$ spikes/taxel
(2.5 mm grid) and $34.73$ spikes/taxel (5 mm grid).

Spatial coding (spike frequency per taxel) reveals an additional topological
structure to the data (Fig. 6, third row), wherein certain taxels seem to fire
at higher rates than others (taxels 10, 17 and 27 for the 2.5 mm texture for
instance). This could be due to their location on the sensor being directly in
the path of the texture’s raised bumps, leading to more events being produced
for those taxels.

We also illustrate temporal coding with a rolling window of size $\Delta
t=159\,ms$ (Fig. 6, bottom row), in which the increased response with texture
grid size is visibly apparent. There is also a sharp increase in activity at
the beginning of each sample which corresponds to the start of the sensor’s
horizontal sliding motion, and is likely linked to the coefficient of static
friction between the sensor and texture.

###  IV-B Texture classification - Artificial textures

Figure 7: Confusion matrices for artificial texture classification with each
of the 4 encoding methods: Intensive (Top Left), Spatial (Top Right), Temporal
(Bottom Left) and Spatiotemporal (Bottom Right).

For each encoding method (see Section III-C), the resulting data is classified
with a KNN algorithm to identify the 11 3d-printed textures.

Note that temporal and spatiotemporal encodings are both parameterised (see
Section III-C3, III-C4), and their parameters ($\Delta t$ for temporal
encoding, $cos\theta$ and $\tau$ for spatiotemporal encoding) are optimised to
maximize classification accuracy. Temporal encoding contains a single
parameter $\Delta t$, allowing us to perform a brute-force optimization over a
range of values from 1 to 200 ms. We set an upper threshold on the $\Delta t$
parameter both to accelerate the optimization procedure and ensure the method
is distinct from intensive coding (which corresponds to $\Delta t=T$, where T
is the full duration of the sample). This optimization process give us a value
of $\Delta t=159\,ms$.

In the case of the spatiotemporal encoding, we perform a Bayesian optimization
over the two parameters $cos\theta$ and $\tau$. The domain space we apply
optimization over is $cos\theta=0-1$ and $\tau=10-100ms$. We place an upper
limit on the $\tau$ parameter for faster convergence and to ensure
spatiotemporal coding is distinct from spatial coding. The parameters converge
over 1000 optimization epochs to values of $cos\theta=0.4$ and $\tau=76ms$.

Coding method | Performance - Artificial textures (%)
---|---
Intensive | 78.5 $\pm$ 41.1
Spatial | 95.5 $\pm$ 20.6
Temporal | 98.1 $\pm$ 13.7
Spatiotemporal | 98.3 $\pm$ 13
TABLE I: Comparison of different encoding methods using leave-one-out cross-
validation for texture classification.

First, we run a simple 80/20 train-test split on the gathered dataset, and
attempt to classify the textures using a KNN classifier (k=4). The results are
presented in the form of confusion matrices (Fig. 7), and provide insight into
each coding method’s performance. Visually, we can observe that the intensive
method is the least effective, particularly for the smoother textures (bump
diameter 0-2 mm). The spatial coding method seems more accurate, though there
is still some slight inaccuracy in the classification of smoother textures
(bump diameter 0-1.5 mm). Temporal and spatiotemporal coding appears to
provide near-perfect classification accuracy, indicating that sliding the
sensor over the textures produces discriminable time-dependent features.

To obtain a more complete picture of the performance of each coding method, we
perform a leave-one-out cross-validation over the gathered dataset. Results
are presented in table I and mirror the results displayed in the train-test
split confusion matrices. Intensive coding performs worst, with spatial
information improving performance significantly. However both temporal and
spatiotemporal coding have the highest accuracy, indicating the sensor
produces features in the temporal domain which are likely most representative
of the classified textures.

###  IV-C Texture classification - Natural textures

Here we seek to test the sensor’s performance more thoroughly on texture
recognition by replacing the 3d-printed textures with natural textures, listed
below in Table II.

We again run a KNN classifier (k=4) over 20 runs with an 80/20 train-test
split (giving 4 test runs for each texture). The results are presented in the
form of confusion matrices for each coding method (Fig. 8).

Figure 8: Confusion matrices for natural texture classification with each of
the 4 encoding methods: Intensive (Top Left), Spatial (Top Right), Temporal
(Bottom Left) and Spatiotemporal (Bottom Right).

Classification accuracy appears generally strong, with Intensive coding the
weakest method, and Temporal and Spatiotemporal coding performing the best, as
was the case for artificial textures. This is confirmed through a leave-one-
out cross-validation, with results displayed in Table III. One significant
difference between the natural and artificial textures is that natural
textures appear more irregular in their misclassifications. The artificial
textures lie along a regularly designed gradient of coarseness, and thus
classification errors occur most often between neighbouring classes. In the
case of natural textures, there is no such regularity or order and
classification errors appear to be spread more widely across the confusion
matrices (Fig. 8).

It is also interesting to note that classification errors vary between coding
methods. For instance for Intensive and Spatial coding, fake fur and felt
(textures 6 and 17) are confused, whereas Temporal and Spatiotemporal coding
can distinguish them. Equally, liquid satin (texture 7) is misclassified as
foam (texture 1) only when Spatiotemporal coding is applied (Fig. 8). The fact
that these misclassifications vary between coding methods likely stems from
the fact that each method produces a distinct set of features.

##  V Discussion

The neuroTac sensor developed here is a neuromorphic optical tactile sensor,
with a data output consisting of multi-taxel spike trains. The spike trains
were encoded using four different bio-inspired encoding mechanisms: Intensive,
Spatial, Temporal and Spatiotemporal. We validated the sensor through a
texture classification task, in which 11 3d-printed textures (grid sizes 0-5
mm in steps of 0.5 mm) and 20 natural textures were discriminated using a KNN
classification algorithm.

We chose texture discrimination as a validation task as it has been a key
experimental procedure in psychophysical and neuroscientific studies of touch.
Here we will attempt to contrast our results with existing theories in these
fields, and identify limitations as well as gaps for further investigation.

We found that applying spatial coding to the neuroTac data produced high
classification accuracy for coarser textures (2.5-5 mm grid size), but
performed less well for smoother textures (0-2.5 mm grid size). This result
seems to concur with Katz’s duplex theory of human tactile perception, with
spatial resolution being important for rougher textures, and vibration cues
taking over as textures get smoother [27]. The stronger performance of the
temporal and spatiotemporal coding methods for smooth textures could be linked
to their detection of high frequency vibrational cues.

Texture | Texture | Texture | Texture
---|---|---|---
number | name | number | name
1 | Foam | 11 | Flarefree net
2 | Plywood | 12 | Embroidered cotton
3 | MDF | 13 | Shirt canvas
4 | Acrylic | 14 | Fairydust organza
5 | Wool | 15 | Sequins allover
6 | Fake fur | 16 | Metallic mesh
7 | Liquid satin | 17 | Felt
8 | Shimmer organza | 18 | Needlecord
9 | Tweed | 19 | Fleece
10 | Lace | 20 | Microdot foil
TABLE II: Natural textures and their corresponding number. Coding method | Performance - natural textures (%)
---|---
Intensive | 69.8
Spatial | 85.5
Temporal | 93.0
Spatiotemporal | 92.8
TABLE III: Comparison of the 4 different coding methods using leave-one-out
cross-validation for texture classification of the 20 natural textures.

It has also been suggested in a recent study that the specific timing of
spikes carries significance for texture discrimination, as it could encode
spatial frequency features of the textures being contacted [28]. The
spatiotemporal coding method described here could capture these precise spike
timings, with a resolution dependent on the time constant parameter $\tau$.
Here $\tau$ was optimized to be 76 ms, which could indicate an approximate
timescale for the frequency features of the textures considered.

Human touch and proprioception provide much richer and more complex data than
that explored here, and most theories suggest that human texture recognition
relies on a complex coding strategy involving inputs from several
mechanoreceptors (Pacinian corpuscles, Merkel cells) [28]. Further studies
with the neuroTac could include methods for simulating input from these
biological afferents, for instance by feeding taxel positions as an input to
spiking neuron models. This would open the path to more complex coding
algorithms to be applied to the sensor’s neuromorphic data which could improve
generalizability and robustness.

##  VI Conclusion

We presented a neuromorphic optical tactile sensor and demonstrated its
performance on a texture classification task. Four bio-inspired spike encoding
mechanisms were investigated, which suggested information about texture
coarseness is encoded in the timing of spikes. The neuroTac’s fast spike-based
output could lead to a step forward in the areas of robotic manipulation and
prosthetics.

## References

  * [1] K. Zaghloul and K. Boahen.  A silicon retina that reproduces signals in the optic nerve.  Journal of neural engineering, 3(4):257, 2006.
  * [2] R. Li and E. Adelson.  Sensing and recognizing surface textures using a gelsight sensor.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1241–1247, 2013.
  * [3] Jasper Wollaston James, Nicholas Pestell, and Nathan F Lepora.  Slip detection with a biomimetic tactile sensor.  IEEE Robotics and Automation Letters, 3(4):3340–3346, 2018.
  * [4] S. Liu and T. Delbruck.  Neuromorphic sensory systems.  Current opinion in neurobiology, 20(3):288–295, 2010.
  * [5] M. van Rossum.  A novel spike distance.  Neural computation, 13(4):751–763, 2001.
  * [6] B. Ward-Cherrier, N. Pestell, L. Cramphorn, B. Winstone, M. Giannaccini, J. Rossiter, and N. Lepora.  The tactip family: Soft optical tactile sensors with 3d-printed biomimetic morphologies.  Soft robotics, 5(2):216–227, 2018.
  * [7] R. Dahiya, G. Metta, M. Valle, and G. Sandini.  Tactile sensing—from humans to humanoids.  Robotics, IEEE Trans. on, 26(1):1–20, 2010.
  * [8] R. Johansson and R. Flanagan.  Coding and use of tactile signals from the fingertips in object manipulation tasks.  Nature Reviews Neuroscience, 10(5):345–359, 2009.
  * [9] G. Indiveri, B. Linares-Barranco, T. Hamilton, A. Van Schaik, R. Etienne-Cummings, T. Delbruck, S. Liu, P. Dudek, P. Häfliger, S. Renaud, et al.  Neuromorphic silicon neuron circuits.  Frontiers in neuroscience, 5:73, 2011.
  * [10] S. Raspopovic, M. Capogrosso, F. Petrini, M. Bonizzato, J. Rigosa, G. Di Pino, J. Carpaneto, M. Controzzi, T. Boretius, E. Fernandez, et al.  Restoring natural sensory feedback in real-time bidirectional hand prostheses.  Science translational medicine, 6(222):222ra19–222ra19, 2014.
  * [11] L. Osborn, A. Dragomir, J. Betthauser, C. Hunt, H. Nguyen, R. Kaliki, and N. Thakor.  Prosthesis with neuromorphic multilayered e-dermis perceives touch and pain.  Science Robotics, 3(19):eaat3818, 2018.
  * [12] F. Bergner, P. Mittendorfer, E. Dean-Leon, and G. Cheng.  Event-based signaling for reducing required data rates and processing power in a large-scale artificial robotic skin.  In 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 2124–2129. IEEE, 2015.
  * [13] C. Bartolozzi, P. Ros, F. Diotalevi, N. Jamali, L. Natale, M. Crepaldi, and D. Demarchi.  Event-driven encoding of off-the-shelf tactile sensors for compression and latency optimisation for robotic skin.  In 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 166–173. IEEE, 2017.
  * [14] W. Lee, S. Kukreja, and N. Thakor.  A kilohertz kilotaxel tactile sensor array for investigating spatiotemporal features in neuromorphic touch.  In 2015 IEEE Biomedical Circuits and Systems Conference (BioCAS), pages 1–4. IEEE, 2015.
  * [15] C. Oddo, L. Beccai, M. Felder, F. Giovacchini, and M. Carrozza.  Artificial roughness encoding with a bio-inspired mems-based tactile sensor array.  Sensors, 9(5):3161–3183, 2009.
  * [16] U. Rongala, A. Mazzoni, and C. Oddo.  Neuromorphic artificial touch for categorization of naturalistic textures.  IEEE transactions on neural networks and learning systems, 28(4):819–829, 2015.
  * [17] Fariborz Baghaei Naeini, Aamna Alali, Raghad Al-Husari, Amin Rigi, Mohammad K AlSharman, Dimitrios Makris, and Yahya Zweiri.  A novel dynamic-vision-based approach for tactile sensing applications.  IEEE Transactions on Instrumentation and Measurement, 2019.
  * [18] Nawid Jamali and Claude Sammut.  Majority voting: Material classification by tactile sensing using surface texture.  IEEE Transactions on Robotics, 27(3):508–521, 2011.
  * [19] Alison I Weber, Hannes P Saal, Justin D Lieber, Ju-Wen Cheng, Louise R Manfredi, John F Dammann, and Sliman J Bensmaia.  Spatial and temporal codes mediate the tactile perception of natural textures.  Proceedings of the National Academy of Sciences, 110(42):17107–17112, 2013.
  * [20] Danfei Xu, Gerald E Loeb, and Jeremy A Fishel.  Tactile identification of objects using bayesian exploration.  In 2013 IEEE International Conference on Robotics and Automation, pages 3056–3061. IEEE, 2013.
  * [21] Mandayam A Srinivasan, JM Whitehouse, and Robert H LaMotte.  Tactile detection of slip: surface microgeometry and peripheral neural codes.  Journal of neurophysiology, 63(6):1323–1332, 1990.
  * [22] C. Chorley, C. Melhuish, T. Pipe, and J. Rossiter.  Development of a tactile sensor based on biologically inspired edge encoding.  In Advanced Robotics, ICAR Int. Conf. on, pages 1–6. IEEE, 2009\.
  * [23] P. Lichtsteiner, C. Posch, and T. Delbruck.  A 128x128 120 db 15 mu s latency asynchronous temporal contrast vision sensor.  IEEE journal of solid-state circuits, 43(2):566–576, 2008.
  * [24] C. Brandli, R. Berner, M. Yang, S. Liu, and T. Delbruck.  A 240$\times$ 180 130 db 3 $\mu$s latency global shutter spatiotemporal vision sensor.  IEEE Journal of Solid-State Circuits, 49(10):2333–2341, 2014.
  * [25] Tetsu Miyaoka, Tadaaki Mano, and Masahiro Ohka.  Mechanisms of fine-surface-texture discrimination in human tactile sensation.  The journal of the acoustical society of America, 105(4):2485–2492, 1999.
  * [26] C. Houghton and K. Sen.  A new multineuron spike train metric.  Neural computation, 20(6):1495–1511, 2008.
  * [27] M. Hollins and R. Risner.  Evidence for the duplex theory of tactile texture perception.  Perception & psychophysics, 62(4):695–705, 2000.
  * [28] E. Mackevicius, M. Best, H. Saal, and S. Bensmaia.  Millisecond precision spike timing shapes tactile perception.  Journal of Neuroscience, 32(44):15309–15317, 2012.


# Duality Symmetry, Two Entropy Functions, and an Eigenvalue Problem in Gibbs’
Theory

Jeffrey Commons  Department of Physics, University of Washington, Seattle, WA
98195-1560, USA  Ying-Jen Yang  Department of Applied Mathematics, University
of Washington, Seattle, WA 98195-3925, USA  Hong Qian  Department of Applied
Mathematics, University of Washington, Seattle, WA 98195-3925, USA

###### Abstract

We generalize the convex duality symmetry in Gibbs’ statistical ensemble
formulation, between Massieu’s free entropy $\Phi_{V,N}(\beta)$ and the Gibbs
entropy $\varphi_{V,N}(u)$ as a function of mean internal energy $u$. The
duality tells us that Gibbs thermodynamic entropy is to the law of large
numbers (LLN) for arithmetic sample means what Shannon’s information entropy
is to the LLN for empirical counting frequencies. Following the same logic, we
identify $u$ as the conjugate variable to counting frequency, a Hamilton-
Jacobi equation for Shannon entropy as an equation of state, and suggest an
eigenvalue problem for modeling statistical frequencies of correlated data.

###### pacs:

In recent years, through activities in stochastic thermodynamics and applying
the mathematical theory of large deviations to the asymptotic behavior of the
LLN, “the essence of a thermodynamic description is not found in its
connection to conservation laws, microscopic reversibility, or the equilibrium
state relation they entail, despite the central role those play” in the
current teaching; see esmith  and the references cited within. In addition to
the derivation of fluctuation relations in connection to various
nonequilibrium entropy productions, more recent results on the unification of
stochastic chemical kinetics and Gibbsian thermochemistry ge-qian-16 , as well
as the probabilistic elucidation of T. L. Hill’s thermodynamics of small
systems (1963) hill-book ; lu-qian-20 , have all further substantiated the
above claim.

The Large Deviations Theory (LDT) from probability oono ; dembo  provides the
mathematical underpinning for understanding statistical observables, of which
the arithmetic mean values and empirical relative frequencies (proportions)
are the two most important types. They have asymptotic behaviors that are
independent of the details of a specific probabilistic model vulpiani-book ;
in fact we shall show that the Gibbsian statistical thermodynamics can be
framed, independent of mechanics, as the asymptotic statistical laws for
counting frequencies. The mathematics reveals an important duality structure
between two functions, a large deviation rate function (LDRF) and its cumulant
generating function (CGF). In statistical mechanics, they correspond to Gibbs’
entropy as a function of all extensive variables and the Massieu-Guggenheim
entropy as a function of all intensive, conjugate variables lu-qian-20 .

LDT offers a clear contradistinction between the two best-known entropies,
that of Gibbs and that of Shannon shannon-book . Although writings on the
subject are already vast, here we give a more coherent presentation showing
that the Gibbs entropy, defined via temperature derivative of canonical
partition function as in standard textbooks, is to the LLN for arithmetic mean
values what the Shannon entropy is to the LLN for empirical counting
frequencies. They are known as Crámer’s theorem and the Sanov theorem in
mathematics respectively dembo .

With this newfound perspective and by explicating the duality for the Shannon-
Sanov relative entropy hobson ; shore-johnson , we show that the conjugate
variables to “relative frequencies” $\\{\nu_{i}\\}$ are a set of “energy
parameters” in $k_{B}T$ units, $\\{u_{j}\\}$. The duality turns out to be
precisely Gibbs’ canonical theory. We further consider repeated measurements
on a time-correlated Markov process and find a novel eigenvalue problem for
$\\{\sqrt{\nu_{i}}\\}$.

Duality symmetry between $\psi_{V,N}(\beta)$ and $\varphi_{V,N}(u)$ — It is
important to recall that as thermodynamic potential functions, the Gibbs
entropy is a function of an extensive variable, the internal energy, and
Helmholtz free energy or Massieu free entropy are functions of the intensive
conjugate variable $T$. There is a convex duality between these two functions
rockafellar .

Consider a mechanical system consisting of $N$ particles, ${\bf
x}=(x_{1},\cdots,x_{N})$, in a container with volume $V$ under temperature
$T=(k_{B}\beta)^{-1}$, with a potential energy function $U_{N}({\bf x})$. The
logarithm of the canonical partition function

$\Phi_{V,N}(\beta)=k_{B}\log\int_{V^{N}}e^{-\beta U_{N}({\bf x})}{\rm d}{\bf
x}$ (1)

has been known as Massieu’s free entropy. See lu-qian-20  for a very recent
discussion. Then the equilibrium entropy of the statistical mechanical system
according to Gibbs’ theory is

$\displaystyle S(\beta)$ $\displaystyle=$
$\displaystyle\frac{\partial[\beta^{-1}\Phi_{V,N}(\beta)]}{\partial\beta^{-1}}=-\beta\frac{\partial\Phi_{V,N}}{\partial\beta}+\Phi_{V,N}(\beta),$
(2)

which can be expressed in terms of Legendre-Fenchel transform (LFT) from
convex analysis:

$\displaystyle\varphi_{V,N}(u)$ $\displaystyle\equiv$
$\displaystyle-S\big{(}\beta(u)\big{)}=-\inf_{\beta}\Big{\\{}\beta
u+\Phi_{V,N}(\beta)\Big{\\}}$ (3a) $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{l}\displaystyle-\beta
u-\Phi_{V,N}(\beta)\\\\[6.0pt] \displaystyle
u(\beta)=-\partial\Phi_{V,N}(\beta)/\partial\beta,\end{array}\right.$ (3d)

in which function $u(\beta)$ is monotonic, therefore $\beta(u)$ is well
defined. Since $\Phi_{V,N}(\beta)$ is a convex function of $\beta$, Eq. 3 has
a dual expression:

$\Phi_{V,N}(\beta)=-\inf_{u}\Big{\\{}\beta u+\varphi_{V,N}(u)\Big{\\}},$ (4)

and a corresponding pair of differential equations showing $\varphi$ and
$\Phi$ as thermodynamic potential functions for their respective independent
variables, with corresponding “forces” $\beta$ and $u$:

$\displaystyle{\rm d}\varphi_{V,N}(u)$ $\displaystyle=$
$\displaystyle\left(\frac{{\rm d}\varphi_{V,N}(u)}{{\rm d}u}\right){\rm
d}u=-\beta(u){\rm d}u,$ (5a) $\displaystyle{\rm d}\Phi_{V,N}(\beta)$
$\displaystyle=$ $\displaystyle\left(\frac{{\rm d}\Phi_{V,N}(\beta)}{{\rm
d}\beta}\right){\rm d}\beta=-u(\beta){\rm d}\beta.$ (5b)

The duality in Eq. 5 can be further extended to intensive conjugate variables
$p$ and $\mu$ which correspond to extensive variables $V$ and $N$. This leads
to a duality symmetry between the fundamental thermodynamic relation and the
Hill-Gibbs-Duhem (HGD) equation in nanothermodynamics lu-qian-20 ,
respectively. In the classical macroscopic thermodynamic limit, the HGD
equation becomes the Gibbs-Duhem equation, with the duality symmetry broken:
When $\varphi_{V,N}(u)$ is an Eulerian homogeneous degree-$1$ function of
$u,V,N$, its complete LFT $\psi(\beta,p,\mu)$ becomes singular except on a
hyper-surface that yields the equation of state (EoS):

$\mu=\frac{1}{N}\sup_{V}\Big{\\{}pV-\beta^{-1}\Phi_{V,N}(\beta)\Big{\\}}.$ (6)

The $\mu$ is a function of $p$ and $\beta$; the $\psi(\beta,p,\mu)$ vanishes.
In the familiar standard thermodynamics notations, ${\Phi=-F/T=-(U-TS)/T}$
where ${1/T=\partial S/\partial U}$, ${\mu=G/N=(pV+F)/N}$ where $p=-\partial
F/\partial V$. They correspond to Eqs. 4 and 6, with the infimum and supremum
carried out by differentiating convex functions.

The Gibbs and the Shannon-Sanov entropies, currently their relation — For a
given probability distribution over a set of $n$ discrete events with
probability mass function $p_{i}$, $1\leq i\leq n$, Shannon’s entropy for
empirical frequencies $\\{\nu_{i}\\}$ is a special case of the more general
relative entropy hobson ; shore-johnson

$H\big{(}\\{\nu_{i}\\}\big{)}=\sum_{i=1}^{n}\nu_{i}\log\frac{\nu_{i}}{p_{i}}\
\text{ or }\int_{\mathbb{R}}f(x)\log\frac{f(x)}{p(x)}{\rm d}x.$ (7)

Shannon’s original formula is taken relative to a uniform probability:
$-H(\nu)=-\sum_{i=1}^{n}\nu_{i}\log\nu_{i}-\log n$. There is an extra, non-
consequential constant $-\log n$. The first term is interpreted as the
“missing information” w.r.t. the maximum information $-\log n$.

According to standard textbooks on equilibrium statistical mechanics
kersonhuang , the Gibbs entropy is defined by starting from a given density of
states in terms of energy $\Omega(u)$, $u\in\mathbb{R}$,

$F(\beta)=-\beta^{-1}\log\int_{\mathbb{R}}\Omega(u)e^{-\beta u}{\rm d}u,\
S=-\frac{\partial F(\beta)}{\partial(1/\beta)},$ (8)

where $S$ is the Gibbs entropy. $F$ is Helmholtz free energy and
$\beta=(k_{B}T)^{-1}$ is inverse temperature. A canonical distribution is

$f(E)=\Omega(u)e^{-\beta E}\left(\int_{\mathbb{R}}\Omega(E)e^{-\beta E}{\rm
d}E\right)^{-1},$

from which one has the well-known equation

$S=-\frac{\partial
F(\beta)}{\partial(1/\beta)}=-\int_{\mathbb{R}}f(E)\log\left(\frac{f(E)}{\Omega(E)}\right){\rm
d}E.$ (9)

If one identifies $\Omega(E)$ with the uniform density $1/N$ and integration
w.r.t. $E$ with summation over $n$; then the Gibbs entropy and the Shannon
entropy are considered to be the same.

Beyond the current understanding — One has found an agreement between the
mathematical expressions for the two entropies. However, a fundamental
question remains unanswered: Why should the mathematical function in (7) be
minimized? Though from thermodynamics originally, this principle must have a
deeper root in statistics; the answer has to be found in the mathematical
logic of probability.

There is a new principle that calls for a revisit of the two entropies: The
entropy function arises as the asymptotic exponent for a law of large numbers.
There are many different laws of large numbers, thus there are many different
forms of entropy function, c.f. Eq. 21 below. Gibbs’ and Shannon’s are merely
two special ones for the two most widely studied LLN. We shall consider two
i.i.d. sampling problems: (i) A discrete random event with probability mass
function $\pi_{k}$, $1\leq k\leq n$, and (ii) a continuous real-valued random
variable $X$ with probability density function $\Omega(x)$.

Shannon’s relative entropy as a LDRF — For the discrete problem, let us assume
that within total $M$ i.i.d. samples, there are $m_{k}$ number of observations
of the $k^{th}$ event: ${m_{1}+\cdots+m_{n}=}M$. It is clear that if one
carries out another $M$ observations, one would obtain another set of
${\\{m_{k};1\leq k\leq n\\}}$. So the $m_{k}$’s are themselves random
variables with a joint probability distribution

$\frac{M!}{m_{1}!\cdots m_{n}!}\pi_{1}^{m_{1}}\cdots\pi_{n}^{m_{n}}.$ (10)

As $M\to\infty$, the relative frequency $\nu_{k}:=\frac{m_{k}}{M}\to\pi_{k}$:
The proportion of observing the $k^{th}$ event approaches to its probability.
This is known as Borel’s LLN.

The Large Deviation Principle says that the distribution in (10), as
$M\to\infty$, has an asymptotic expression
$e^{-M\varphi(\nu_{1},\cdots,\nu_{n})}$, where
$\varphi(\nu_{1},\cdots,\nu_{n})$ is exactly the $H(\\{\nu_{i}\\})$ in Eq. 7,
with $p_{k}=\pi_{k}$.

The Gibbs entropy is also a LDRF — Now for the continuous problem, let
$X_{1},X_{2},\cdots,X_{M}$ be $M$ i.i.d. samples of a random variable $X$,
with probability density function $\Omega(x)$. Then the LLN for the mean value
says that

$\overline{X}_{M}:=\frac{X_{1}+\cdots+X_{M}}{M}\to\mathbb{E}[X]=\int_{\mathbb{R}}x\Omega(x){\rm
d}x,$

as $M\to\infty$. It is clear that if one obtains another $M$ samples, the
$\overline{X}_{M}$ will be different. So $\overline{X}_{M}$ itself is a random
variable with a probability density function $f_{\overline{X}_{M}}(x)$.

Again, the LDT says that as $M\to\infty$, $f_{\overline{X}_{M}}(x)$ has an
asymptotic expression $e^{M\psi(x)}$. To compute the exponent $\psi(x)$, we
employ the method of CGF,

$\Lambda(\beta)=\log\int_{\mathbb{R}}\Omega(x)e^{-\beta x}{\rm d}x,$ (11)

and note that the CGF of $\overline{X}_{M}$ is
$\Lambda_{\overline{X}_{M}}(\beta)=M\Lambda(\beta/M)$. Then

$\Lambda(\beta)=\frac{1}{M}\Lambda_{\overline{X}_{M}}\big{(}M\beta\big{)}=\sup_{x}\big{\\{}\psi(x)-\beta
x\big{\\}},$ (12)

which implies

$\psi(x)=\inf_{\beta}\big{\\{}\beta x+\Lambda(\beta)\big{\\}}.$ (13)

Now to carry out the optimization problem in (13), let us use calculus to find
the optimal $\beta^{*}$: $x+\Lambda^{\prime}(\beta)=0$. Therefore,

$\left\\{\begin{array}[]{ccl}\psi&=&-\beta\Lambda^{\prime}(\beta)+\Lambda(\beta)=\displaystyle\frac{\partial}{\partial(1/\beta)}\Big{[}\beta^{-1}\Lambda(\beta)\Big{]},\\\
x&=&\displaystyle-\frac{\partial}{\partial\beta}\Lambda(\beta).\end{array}\right.$
(14)

Eq. (14) gives the function $\psi(x)$ in a parametric form. We see that if we
identify $-\beta^{-1}\Lambda(\beta)$ as $F(\beta)$, then the $\psi$ in (14) is
the Gibbs entropy, and $x$ is the internal energy.

Therefore, as the asymptotic exponent for the probability of large deviations,
the Shannon entropy is to Borel’s LLN for empirical counting frequencies what
the Gibbs entropy is to the standard LLN for arithmetic mean values.

Gibbs-Kirkwood potential energy and the dual to the Shannon-Sanov entropy —
Through Legendre transforms, the theory of classical thermodynamics introduces
intensive conjugate variables to extensive quantities, such as $\beta$ to
internal energy $U$, $p$ to volume $V$, and $\mu$ to particle number $N$. The
same conjugate variables enter Gibbs’ theory of statistical mechanics via
Laplace transforms of various partition functions corresponding to extensive
variables. The result from the previous sections shows clearly that a duality
structure between the LDRF and its CGF plays a central role in a new
understanding of Gibbs’ and Shannon’s entropies.

By the same logic the following question immediately arises: What is the
corresponding CGF for Shannon’s relative entropy, i.e., Sanov’s LDRF in (7)?
We carry out its LFT using $u_{i}$ to denote the conjugate variable of
$\nu_{i}$ under the normalization constraint $\sum_{i}^{n}\nu_{i}=1$:

$\displaystyle\Xi\big{(}\\{u_{i}\\}\big{)}$ $\displaystyle=$
$\displaystyle\inf_{\\{\nu_{i}\\}}\left\\{\sum_{i=1}^{n}u_{i}\nu_{i}+H\big{(}\\{\nu_{i}\\}\big{)}\right\\}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\nu_{i}^{*}\log\left(\frac{\nu_{i}^{*}}{p_{i}e^{-u_{i}}}\right)=\lambda=-\log\sum_{i=1}^{n}p_{i}e^{-u_{i}},$

in which $\lambda$ appears as the Lagrange multiplier for the constrained
optimization ${\nu_{i}^{*}=p_{i}e^{-u_{i}+\lambda}}$. For fixed $u_{i}$, the
$\nu_{i}$ are determined from ${u_{i}=-\partial
H(\\{\nu_{i}\\})/\partial\nu_{i}+\lambda}$. The additive constant $\lambda$ to
the $u_{i}$’s represents a gauge freedom, a feature that is absent in the
traditional Legendre transform. Eq. Duality Symmetry, Two Entropy Functions,
and an Eigenvalue Problem in Gibbs’ Theory has the exact form of Helmholtz
free energy from standard statistical mechanics, when a set of “energy”
parameters in $k_{B}T$ units, $\\{u_{i}\\}$ is given:

$\Xi\big{(}\\{u_{i}\\}\big{)}=\underbrace{\sum_{i=1}^{n}\nu_{i}^{*}u_{i}}_{\text{mean
internal
energy}}-\underbrace{\left(-\sum_{i=1}^{n}\nu^{*}_{i}\log\frac{\nu^{*}_{i}}{p_{i}}\right)}_{\text{entropy}}.$
(16)

More precisely, the mean internal energy term in (16) should be understood as
the mean value of Kirkwood’s potential of mean force $u_{i}$ kirkwood . The
conjugate variables $u_{i}$ therefore are interpreted as a “potential
function” that represents the statistical deviations of the empirical
frequencies $\\{\nu_{i}\\}$ away from the true probability $\\{p_{i}\\}$ in
the large-deviation scenario.

We shall call the set of $\\{u_{i}\\}$ introduced in (Duality Symmetry, Two
Entropy Functions, and an Eigenvalue Problem in Gibbs’ Theory), the conjugate
variables to the counting frequencies $\\{\nu_{i}\\}$, Gibbs-Kirkwood
potential energies; they comes with a “natural” additive constant $\lambda$
due to the normalization of frequencies $\\{\nu_{i}\\}$. One then recovers the
duality:

$\displaystyle{\rm d}H\big{(}\\{\nu_{i}\\}\big{)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\left(\frac{\partial
H}{\partial\nu_{i}}\right){\rm
d}\nu_{i}=-\sum_{i=1}^{n}u_{i}\big{(}\\{\nu_{i}\\}\big{)}{\rm d}\nu_{i},\hskip
19.91684pt$ (17a) $\displaystyle{\rm d}\Xi\big{(}\\{u_{i}\\}\big{)}$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{n}\left(\frac{\partial\Xi}{\partial
u_{i}}\right){\rm d}u_{i}=\sum_{i=1}^{n}\nu_{i}\big{(}\\{u_{i}\\}\big{)}{\rm
d}u_{i},$ (17b)

where $u_{i}(\\{\nu_{i}\\})=$ $-\log(\nu_{i}/p_{i})+\lambda$.

The Massieu free entropy as a function of $\beta$ is the LFT of the Gibbs
entropy as a function of energy; the Helmholtz free energy as a function of
$\\{u_{i}\\}$ is the LFT of Shannon-Sanov entropy as a function of
frequencies.

Entropy analysis as in thermodynamics — To compare and contrast with
statistical thermodynamics, we consider a proper entropy function of non-
normalized absolute frequencies $\\{m_{i}\\}$,
$S(\\{m_{i}\\})=-MH(\\{m_{i}/M\\})$, which is an Eulerian homogeneous function
of degree $1$: ${S\big{(}\\{\alpha m_{i}\\}\big{)}=\alpha S(\\{m_{i}\\})}$.
Its Legendre transform yields a singular function:

$u_{n}=\frac{1}{m_{n}}\underbrace{\left(S\big{(}\\{m_{i}\\}\big{)}-\sum_{j=1}^{n-1}u_{j}m_{j}\right)}_{\text{
a function of $u_{1},\cdots,u_{n-1}$, and $m_{n}$}},$ (18)

which resembles the Gibbs-Duhem equation, with both integral and differential
forms:

$\displaystyle S\big{(}\\{m_{i}\\}\big{)}=\sum_{i=1}^{n}m_{i}u_{i},$ (19a)
$\displaystyle m_{i}\big{(}\\{u_{k}\\}\big{)}=\frac{p_{i}e^{-u_{i}}}{Z}\text{
where }Z=\frac{\sum_{i=1}^{n}p_{i}e^{-u_{i}}}{\sum_{j=1}^{n}m_{j}},\hskip
28.45274pt$ (19b)
$\displaystyle\sum_{i=1}^{n}m_{i}\big{(}\\{u_{k}\\}\big{)}{\rm
d}u_{i}=\sum_{i=1}^{n}{\rm d}m_{i}\big{(}\\{u_{k}\\}\big{)}=0.$ (19c)

We have just generalized the concepts of Gibbs function and chemical
potentials in thermochemistry to statistical absolute counting frequencies, in
terms of an entropy function $S(\\{m_{i}\\})$ and the Gibbs-Kirkwood potential
functions as the conjugate variables $\\{u_{i}\\}$. Fixing the gauge
$\lambda=0$ gives the EoS, a stationary Hamilton-Jacobi equation of the
entropy function $S$:

$\sum_{i=1}^{n}p_{i}e^{-u_{i}}-1=0,\ u_{i}=\frac{\partial
S(\\{m_{i}\\})}{\partial m_{i}}.$ (20)

Markov chain large deviations and conjugate process — Beyond the i.i.d.
empirical counting, most repeated measurements of a stationary stochastic
signal have a temporal correlation which can be represented by a Markov
process $X_{t}$, $t=0,1,\cdots$, with transition probability matrix
$\\{p_{ij}\\}$ between consecutive measurements. The LDRF of the counting
frequencies along a path is dembo :

$I(\\{\nu_{i}\\})=-\inf_{\text{all
}\mu_{i}>0}\sum_{i=1}^{n}\nu_{i}\log\left(\frac{1}{\mu_{i}}\sum_{j=1}^{n}\mu_{j}p_{ji}\right).$
(21)

The infimum in (21) can be solved from a system of nonlinear equations for the
$\mu_{i}$’s:

$\displaystyle\sum_{i=1}^{n}\nu_{i}\left(\frac{\mu_{k}p_{ki}}{\sum_{j=1}^{n}\mu_{j}p_{ji}}\right)=\nu_{k},\
1\leq k\leq n.$ (22a) $\displaystyle
I\big{(}\\{\nu_{i}\\})=\sum_{i=1}^{n}\nu_{i}\log\left(\frac{\nu_{i}}{v_{i}}\right),\
\nu_{i}=v_{i}\sum_{j=1}^{n}\frac{\mu_{j}p_{ji}}{\mu_{i}}.$ (22b)

Eq. 22 has two very different interpretations: In data statistical terms, the
observed relative frequency $\\{\nu_{i}\\}$ is the stationary distribution of
a Markov chain with transition probability $\hat{p}_{ij}$:

$\hat{p}_{ij}=\frac{\mu_{j}p_{ji}}{\sum_{k=1}^{n}\mu_{k}p_{ki}}=\frac{\mu_{j}\Pr\\{X_{t+1}=i|X_{t}=j\\}}{\Pr\\{X_{t+1}=i\\}},$
(23)

which is the Bayesian posterior distribution for $X_{t}$ given $X_{t+1}=i$,
under the original Markov process with
${\Pr\\{X_{t}=j,X_{t+1}=i\\}=\mu_{j}p_{ji}}$. (22a) finds the hidden
$\\{\mu_{i}\\}$ from the observed $\\{\nu_{i}\\}$. Eq. 22b also has a
thermodynamic interpretation: $\log(\nu_{i}/v_{i})$ in (22b) is a Gibbs-
Kirkwood potential representing the deviation of observed statistical
frequency $\\{\nu_{i}\\}$ from the stationary probability $\\{\pi_{i}\\}$. It
is zero when $\nu_{i}=\pi_{i}$ $\forall i$, which yields $\mu_{k}=\pi_{k}$ and
$\nu_{i}=v_{i}$ from (22). In terms of $\nu$’s and $v$’s, $\hat{p}_{ij}=$
$(\mu_{j}p_{ji}/\mu_{i})v_{i}/\nu_{i}$.

From a mathematical modeling standpoint, finding $\\{p_{jk}\\}$ for an
underlying Markov system and finding a conjugate process $\\{\hat{p}_{jk}\\}$
for the empirical frequency of measurements on the stochastic system are
equivalent. Any given set of observed relative frequencies $\\{\nu_{i}\\}$ is
represented as the stationary distribution of the conjugate Markov process
w.r.t. a hidden measure $\\{\mu_{i}\\}$ that serves as a tilting barato , and
represents the measurement on probability.

For a continuous-time Markov chain with $p_{ij}=(e^{{\bf Q}\Delta t})_{ij}$
where ${\bf Q}=\\{q_{ij}\\}$ is a transition rate matrix, The conjugate
process w.r.t $\\{\mu_{i}\\}$ has transition probability rate matrix

$\hat{q}_{ij}=\frac{\mu_{j}q_{ji}}{\mu_{i}}\text{ if }i\neq j,\
\hat{q}_{ii}=-\sum_{j\neq i}^{n}\hat{q}_{ij},$ (24)

and $\\{\nu_{i}\\}$ is its stationary measure. Actually
${\hat{q}_{ij}=\mu_{j}q_{ji}/\mu_{i}-r_{i}\delta_{ij}}$ where
$r_{i}=\nu_{i}/v_{i}$ has another interpretation: the instantaneous rate of
change in the logarithmic probability density of state $i$ following the
Markov process: ${\rm d}\log\mu_{i}/{\rm
d}t=\sum_{k=1}^{n}\mu_{k}q_{ki}/\mu_{i}=r_{i}$. When $\forall i$
$\mu_{i}=\pi_{i}$, $r_{i}=0$. The larger a $|r_{i}|$ value for a state $i$,
the faster the state probability relaxes to $\pi_{i}$. The Gibbs-Kirkwood
potential for a Markov chain also acquires a kinetic meaning.

Markov chain with detailed balance — The infimum in Eq. 21 can be solved if
$\\{q_{ij}\\}$ has detailed balance, ${\pi_{i}q_{ij}=\pi_{j}q_{ji}}$:

$\sum_{i,j=1}^{n}\frac{\mu_{j}q_{ji}\nu_{i}}{\mu_{i}}\geq\sum_{i=1}^{n}q_{ii}\nu_{i}+2\sum_{i>j}^{n}\big{(}q_{ji}\nu_{i}q_{ij}\nu_{j}\big{)}^{\frac{1}{2}},$
(25)

where the equality holds true when
$\mu_{j}q_{ji}\nu_{i}/\mu_{i}=\mu_{i}q_{ij}\nu_{j}/\mu_{j}$. That is if
$(\mu_{j}/\mu_{i})^{2}$ $=\nu_{j}\pi_{j}/(\nu_{i}\pi_{i})\;\forall i,j$.
Therefore

$\frac{I(\\{\nu_{i}\\})}{\Delta t}=\sum_{i,j=1}^{n}\sqrt{\nu_{i}}\
\hat{\Pi}_{ij}\sqrt{\nu_{j}},$ (26)

in which matrix $\hat{\mbox{\boldmath$\Pi$}}$ is semi-positive definite with
elements $\hat{\Pi}_{ij}=-\sqrt{\pi}_{i}q_{ij}/\sqrt{\pi_{j}}=\hat{\Pi}_{ji}$.
As a LDRF, $I(\\{\nu_{i}\\})\geq 0$ and the equality holds true when
$\nu_{i}=\pi_{i}$. This corresponds to matrix $\\{q_{ij}\\}$ having an
eigenvalue $0$ with eigenvector $(1,\cdots,1)$. The Legendre-Fenchel transform
of (26) yields an eigenvalue problem that solves the $\\{\nu_{i}\\}$ for given
$u_{i}=\partial I/\partial\nu_{i}$, $(1\leq i\leq n)$:

$\sum_{k=1}^{n}\big{[}\
\hat{\Pi}_{ik}-(u_{i}+\lambda)\delta_{ik}\big{]}\sqrt{\nu_{k}}=0,$ (27)

in which $\lambda$ again represents the gauge freedom that arises from the
constraint among $\nu$’s. The $\lambda$ is the largest eigenvalue of matrix
$\big{[}\hat{\mbox{\boldmath$\Pi$}}-\text{diag}(u_{1},\cdots,u_{n})\big{]}$,
with the corresponding eigenvector being $\\{\sqrt{\nu_{i}}\\}$.

Eq. 26 offers a novel narrative for the statistical behavior of correlated
empirical counting frequencies from a time-reversible Markov process: The
empirical relative frequencies are split into $\sqrt{\nu_{i}}$ and
$\sqrt{\nu_{j}}$, they are captured by a symmetric matrix
$\sqrt{\nu_{i}\nu_{j}}$ whose trace is $1$. This resembles von Neumann’s
density operator. (27) further suggests a novel, symmetric eigenvalue problem
for $\sqrt{\nu_{i}}$. The physical significance of this set of results remains
to be explored.

The nature of the Gibbsian theory — We identify Gibbs’ theories of statistical
mechanics and chemical thermodynamics, as a whole, as a very different program
from Boltzmann’s original objective of developing a mechanical theory of heat,
e.g., thermodynamic principles as mathematical consequences of classical
mechanics. Gibbs’ theory, with certain justifications, puts the notion of
temperature which represents stochastic mechanical motions as an a priori
concept and canonical ensemble theory as a first principle. This theory
reproduces all the known relationships in classical equilibrium thermodynamics
by taking the thermodynamic limit; one of the key results in this regard is
the Gibbs-Duhem equation. However, the new theory itself is applicable to
“small systems”. The entirety of statistical chemistry and soft matter
physics, including theories of non-ideal solutions, polymers,
polyelectrolytes, proteins, and the great successes these applications have
brought about, is testimony to the validity of the theory, which is far beyond
Boltzmann’s original vision, if one can surmise. Boltzmann’s theory has a
goal; Gibbs’ is open ended.

The macroscopic limit gives rise to theories of emergent phenomena, such as
phase transitions. This insight originated in Yang and Lee’s work in 1952
yanglee , who departed from then traditional approach of McMillan-Mayer
cluster expansion. By 1972, this had become one of the fundamental tenets for
emergence pwanderson ; vulpiani-book .

The current work finds the Gibbsian theory a new foundation and its
limitation. It is a mathematical consequence within the theory of probability,
not due to the large size thermodynamic limit, but rather the limit law of
large numbers of measurements whose empirical frequency converges to
probability. In other words, the “law” articulated in the Gibbsian theory is
due to universal statistical behavior of frequency as an empirical measurement
for probability, a mathematical rather than physical notion.

Maximum entropy principle re-interpreted — Maximum entropy principle (MEP)
jaynes  and its dynamic counterpart maximum calibre (MaxCal) dill  have a
great number of followers in engineering and data science. The theory of large
deviations offers the very statistical origin for such a principle: The key is
to distinguish the abstract mathematical concept of “probability meansure”
associated with $(\Omega,\mathcal{F},\mathbb{P})$ and the statistical counting
frequency that is widely used as the surrogate for the former. The latter can
be experimentally measured. Very recently, K. A. Dill has called for a
distinction between entropy as a scalar quantity of a system and entropy as a
“function of the relative
frequencies”.111https://math.ucr.edu/home/baez/SMB2021/SMB2021_dill.pdf

Gauge freedom, normalization and free energy of ground state. The Legendre-
Fenchel dual to the Shannon-Sanov entropy of counting frequency
$\\{\nu_{i}\\}$ naturally gives rise to the notion of Gibbs-Kirkwood potential
function $\\{u_{i}\\}$ together with a gauge freedom $\lambda$ which is
mathematically related to the normalization of the $\nu$’s. It defines through
a gauge fixing the “ground” state in which its Helmholtz free energy $\lambda$
is zero. It is important to recognize that the set of $\\{\nu_{i}\\}$ is a
result of measurements approaching infinity in number: Given a set of
$\\{\nu_{i}\\}$ or alternatively given a set of $\\{u_{i}\\}$, one defines a
statistical or “thermodynamic” system with permanence. In classical physics,
such a system is defined mechanically through its energy function. In the
present work, it is replaced by the underlying probability and the LD entropy
function.

Bottom up, top down, and bottom down — Gibbs’ theory, as understood from the
present work, provides a surprisingly penetrating analytic device for carrying
out the vision of R. B. Laughlin laughlin : In the past, the theory at the
“bottom” has always been mechanics, which provides the notion of mechanical
energy and mechanical work, from which thermodynamics arises. But as many
practitioners of biophysical chemistry have recognized schellman , the
relationship between the macromolecular Boltzmann weights (or Kirkwood’s
potential of mean force) one has at hand and the constitutive mechanical
Hamiltonian of point masses is strenuous even at its best. The present work,
in agreement with esmith , again suggests the traditional view might not be
the necessary path to thermodynamic behavior. Statistical law could provide an
alternative origin for the very notion of internal energy, via the Gibbs-
Kirkwood potential function. Further investigations into the counting
frequency of correlated events might reveal more unexpected behavior hoffmann
.

Acknowledgement — We thank R. H. Austin, J. C. Baez, K. A. Dill, C. Jia, and
H. Zhao for helpful discussions.

## References

  * (1) Smith, E. (2020) Intrinsic and extrinsic thermodynamics for stochastic population processes with multi-level large-deviation structure. Entropy, 22, 1137.
  * (2) Ge. H. and Qian, H. (2016) Mesoscopic kinetic basis of macroscopic chemical thermodynamics: A mathematical theory. Phys. Rev. E 94, 052150.
  * (3) Hill, T. L. (1963) Thermodynamics of Small Systems, Dover, New York.
  * (4) Lu, Z. and Qian, H. (2020) Emergence and breaking of duality symmetry in thermodynamic behavior: repeated measurements and macroscopic limit. arXiv:2009.12644.
  * (5) Oono, Y. (1989) Large deviations and statistical physics. Progr. Theor. Phys. Supp. 99, 165–205.
  * (6) Dembo, A. and Zeitouni, O. (1998) Large Deviations Techniques and Applications. 2nd ed., Springer, New York.
  * (7) Chibbaro, S., Rondoni, L. and Vulpiani, A. (2014) Reductionism, Emergence and Levels of Reality, Springer, New York.
  * (8) Shannon, C. E. and Weaver, W. (1949) The Mathematical Theory of Communication, Univ. Ill. Press, Champaign, IL.
  * (9) Hobson, A. (1969) A new theorem of information theory. J. Stat. Phys. 1, 383–391.
  * (10) Shore, J. E. and Johnson, R. W. (1980) Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy. IEEE Trans. Inf. Th. 26, 26–37.
  * (11) Rockafellar, R. T. (1996) Convex Analysis, Princeton Univ. Press, NJ.
  * (12) Huang, K. (1963) Statistical Mechanics, John Wiley & Sons, New York.
  * (13) Kirkwood, J. G. (1935) Statistical mechanics of fluid mixtures. J. Chem. Phys. 3, 300–313.
  * (14) Maes, C. and Netočný, K. (2007) Minimum entropy production principle from a dynamic fluctuation law. J. Math. Phys. 48, 053306.
  * (15) Barato, A. C. and Chetrite, R. (2015) A formal view on level 2.5 large deviations and fluctuation relations. J. Stat. Phys. 160, 1154–1172.
  * (16) Yang, C. N. and Lee, T. D. (1952) Statistical theory of equations of state and phase transitions. I. Theory of condensation. Phys. Rev. 87, 404–409.
  * (17) Anderson, P. W. (1972) More is different. Science 177, 393–396.
  * (18) Jaynes, E. T. (2003) Probability Theory: The Logic of Science, Cambridge Univ. Press, U.K.
  * (19) Ghosh, K., Dixit, P. D., Agozzino, L. and Dill, K. A. (2020) The maximum caliber variational principle for nonequilibria. Annu. Rev. Phys. Chem. 71, 21–238.
  * (20) Laughlin, R. B. (2006) A Different Universe: Reinventing Physics from the Bottom Down, Basic Books, New York.
  * (21) Schellman, J. A. (1997) Thermodynamics, molecules and the Gibbs conference. Biophys. Chem. 64, 7–13.
  * (22) Hoffmann, B. (2011) The Strange Story of the Quantum, Dover, New York.