# arfpy: A python package for density estimation and generative modeling with
adversarial random forests

Kristin Blesch∗1,2 Marvin N. Wright1,2,3

(∗ corresponding author<EMAIL_ADDRESS>
1Leibniz Institute for Prevention Research and Epidemiology – BIPS, Germany;
2Faculty of Mathematics and Computer Science, University of Bremen, Germany;
3Department of Public Health, University of Copenhagen, Denmark)

## Abstract

This paper introduces arfpy, a python implementation of Adversarial Random
Forests (ARF) (Watson et al., 2023), which is a lightweight procedure for
synthesizing new data that resembles some given data. The software arfpy
equips practitioners with straightforward functionalities for both density
estimation and generative modeling. The method is particularly useful for
tabular data and its competitive performance is demonstrated in previous
literature. As a major advantage over the mostly deep learning based
alternatives, arfpy combines the method’s reduced requirements in tuning
efforts and computational resources with a user-friendly python interface.
This supplies audiences across scientific fields with software to generate
data effortlessly.
https://github.com/bips-hb/arfpy

## Keywords

Generative Modeling; Density Estimation; Machine Learning

## Introduction

Generative modeling is a challenging task in machine learning that aims to
synthesize new data which is similar to a set of given data. State of the art
are computationally intense and tuning-heavy algorithms such as generative
adversarial networks (GANs)[1, 2], variational autoencoders [3], normalizing
flows [4], diffusion models [5] or transformer-based models [6]. A much more
lightweight procedure is to use an Adversarial Random Forest (ARF) [7]. ARFs
achieve competitive performance in generative modeling with much faster
runtime [7], yet they do not require the practitioner to have extensive
knowledge of generative modeling.
Further, ARFs are especially useful for data that comes in a table format,
i.e., tabular data. That is because ARFs are based on random forests [8] which
leverage the advantages that tree-based methods have over neural networks on
tabular data [9] for generative modeling. Further, as part of the procedure,
ARFs give access to the estimated joint density, which is useful for several
other fields of research, e.g., unsupervised machine learning. For the task of
density estimation, ARFs have been demonstrated to yield remarkable results as
well [7]. In brief, ARFs are a promising methodological contribution to the
field of generative modeling and density estimation, providing a ready-made
solution to generate data for practitioners across fields.
ARFs have already gained some attention in the scientific community [10],
however, the paper by [7] provides the audience with a R software package
only. The machine learning and generative modeling community, however, is
mostly using python as a programming language and to reach a broad audience
more generally, a fast and user-friendly implementation of ARFs in python is
highly desirable. We aim to fill this gap with the presented python
implementation of ARFs.
arfpy is inspired by the R implementation called arf [11], but transfers the
algorithmic structure to match the class-based structure of python code and
takes advantage of computationally efficient python functions. Similar to the
R implementation, separate functions for first fitting the density (FORDE
algorithm [7]) and then generating new data samples (FORGE algorithm [7])
exist. However, in arfpy, the functions are called for an initialized ARF
class object, which is showcased in the usage example below.
Crucially, for practitioners working with python as programming language, the
direct python implementation is more robust and convenient to users than
calling fragile wrappers like rpy2 [12] that aim to make R code running in
python. The benefits of a direct python implementation of ARFs for the
generative modeling community have already been recognized by now. For
example, arfpy is integrated in the data synthesizing framework synthcity
[13].

## Implementation and architecture

### Module Design

The general workflow of generating data with arfpy is (1) to initialize an
object of class arf with real data, (2) estimate the density and (3) sample
new data. This procedure is visualized in Figure 1.
The architecture of arfpy reflects this workflow and we have class arf
building the backbone of the procedure. An instance of class arf takes the
real data set as input and trains an ARF, i.e., learns the actual data’s
structure. To this object, functions to estimate the density (FORDE algorithm
[7], function forde()) and generate data (FORGE algorithm [7], function
forge()) can be applied. This architecture allows users to learn the structure
of the real data once (when initializing the arf class object) and then
flexibly adapt density estimation, e.g., using different parameters, or
repeatedly sampling new data without having to refit the model.

Figure 1: Workflow of arfpy’s core functionalities.

### Methodology Overview

For interested readers, we want to briefly describe the methodological
foundations of ARFs, but refer to [7] for further details. From a given real
data set, first, naive synthetic data is generated (initial generation step)
by sampling from the marginal distributions of the features. Then, a random
forest [8] is fit to distinguish this synthetic from the real data (initial
discrimination step). This procedure, also known as fitting an unsupervised
random forest [14], guides the random forest to learn the dependency structure
in the data. Using this forest, we can sample observations from the leaves of
the trees to generate updated synthetic data (generation step). Subsequently,
a new random forest is fit to differentiate between synthetic and real data
(discrimination step). Drawing on the adversarial idea of GANs, this iterative
procedure of data generation and discrimination will be repeated until the
discriminator cannot distinguish between generated and real data anymore. At
this stage, the accuracy of the forest will be $\leq 0.5$ and the forest is
assumed to have converged, which implies mutually independent features in the
terminal nodes. This drastically simplifies density estimation and generative
modeling, as it allows us to formulate the univariate density for each feature
separately with data in the leaves of the fitted ARF (FORDE algorithm) and
then combine them to the joint density, instead of having to model
multivariate densities. For data generation, we can use this trait to sample a
new observation by drawing a leaf from the forest of the last iteration step
and use the data distributions with parameters estimated from that leaf to
sample each feature separately (FORGE algorithm).

### Example Usage

Let us illustrate the usage of arfpy with a visually intuitive example: We
create data using make_moons from sklearn.datasets, which results in data
along two continuous axes that looks like two moons from different categories.
Statistically speaking, this is a tabular dataset, consisting of both
continuous and categorical features that exhibit a dependency structure. For a
more intuitive understanding of the data, see Figure 2, Panel A. The task of
arfpy is to learn the structure of this given (real) data and generate new
data instances that appear similar.

To initialize the workflow, we need to run relevant imports, including the
import of class arf from the arfpy module, and create the real dataset. The
arf class takes a pandas DataFrame as input, so the real data is pre-processed
to match this requirement. This incorporates setting unique column names
(’dim_1’, ’dim_2’,’label’) and ensuring that ’label’ is stored in the correct
data type ’category’.

⬇

1import pandas as pd

2from sklearn.datasets import make_moons

3from arfpy import arf

4

5moons_X, moons_y = make_moons(n_samples = 3000, noise=0.1)

6df = pd.DataFrame({"dim_1" : moons_X[:,0], "dim_2" : moons_X[:,1], "label" :
moons_y})

7df[’label’] = df[’label’].astype(’category’)

8

9df.head()

10

11#> dim_1 dim_2 label

12#> 1.782717 0.099124 1

13#> 1.087497 0.298744 0

14#> -0.576695 0.801675 0

15#> 0.623931 -0.506896 1

With the real dataset preprocessed as needed, we can proceed with training the
ARF to learn the data’s structure. Creating an object of class arf will
trigger ARF model fitting using the data provided.

⬇

1my_arf = arf.arf(x = df)

2

3#> Initial accuracy is 0.82

4#> Iteration number 1 reached accuracy of 0.36

Because we have used the parameter default verbose = True, the training of
my_arf prints out some interesting information: The initial accuracy, which
corresponds to the accuracy of the random forest in distinguishing real data
from naive synthetic data, is $0.82$. This implies that the random forest is
doing very well in distinguishing real from naive synthetic data and
therefore, we can assume the model to have learned relevant dependencies that
allow the model to make this distinction. Using this forest to sample updated
synthetic data, and fitting a new random forest to distinguish this data from
real data leads to an accuracy of only $0.36$. This accuracy is below the
default threshold of $0.5$, so loosely speaking, the synthetic data generated
with the forest cannot be accurately distinguished from real data, i.e., the
generated data looks like real data, which is the goal the algorithm was
aiming for. In other words, the relevant dependency structures of the real
data have been learned by the forest in the first round of iteration already,
so the algorithm has converged and no further iterations need to be conducted.

After the ARF has converged, we can proceed to estimating the joint density.
Recap that in a converged ARF, the features are mutually independent in the
leaves, which simplifies the challenging multivariate density estimation task
into many simple univariate density estimation tasks. The joint density is
then a factorization of the individual density estimates across leaves in the
ARF. We can call function forde() on the my_arf object to estimate the density
and store the returned dictionary to explore the parameters. The FORDE
dictionary contains the estimated parameters for continuous (key ’cnt’) and
categorical features (key ’cat’). As mentioned in the above paragraph, the
parameters are estimated using the data points in the forest’s leaves, so we
will get estimates for each leaf individually. The parameters for the
categorical features simply correspond to the empirical frequency of
categories in the leaves, so for a more complex example, we can take a look at
the continuous feature’s parameter estimates in FORDE[’cnt’]. We have used the
default distribution (truncated normal distribution) to model continuous
features, so the output will reflect estimates for the mean and standard
deviation for each feature (’dim_1’, ’dim_2’) in each leaf, which is uniquely
identified by ’tree’ and ’nodeid’:

⬇

1FORDE = my_arf.forde()

2

3FORDE[’cnt’].iloc[:,:5].head()

4

5#> tree nodeid variable mean sd

6#> 0 3 dim_1 0.961437 0.214925

7#> 0 3 dim_2 -0.671571 0.028193

8#> 0 11 dim_1 1.040565 0.185581

9#> 0 11 dim_2 -0.621924 0.003328

With the parameters estimated, we can move on to the final step of the
generative modeling task and sample new data instances with the function
forge().
For each instance to be generated, the function randomly samples a leaf from
the forest with weighted probability according to the coverage of real data in
the leaves of the ARF and then uses the parameters estimated through forde()
to sample a new data instance.

⬇

1df_syn = my_arf.forge(n = 1000)

2

3df_syn.head()

4

5#> dim_1 dim_2 label

6#> -0.018004 0.283963 1

7#> 1.734200 -0.085115 1

8#> -0.009840 1.046872 0

9#> 0.868400 -0.352692 1

Calling forge() completes the task of generating synthetic data that mimics
real data. From the generated data table itself, the similarity is hard to
grasp, but we can visually inspect the quality of the synthetic data in Figure
2.

Figure 2: Comparison of real and synthesized data.

## Quality control

The software has been tested through unit tests, which includes testing of
relevant functionalities with various input data sets. The workflow of running
these tests is automated on GitHub actions, but can be run locally and with
customized data sets using the instructions provided in the software
repository. Further, the repository allows to publicly raise questions or
report bugs and gives clear guidelines on how to contribute to the open source
software project are lined out.

## (2) Availability

## Operating system

Platform Independent

## Programming language

Python $\geq$ 3.8

## Additional system requirements

No specific requirements

## Dependencies

numpy $\geq$ 1.20.3, pandas $\geq$ 1.5, scikit-learn $\geq$ 0.24, scipy $\geq$
1.4

## List of contributors

Blesch, Kristina,b;
Wright, Marvin N.a,b,c;
(a) Leibniz Institute for Prevention Research and Epidemiology – BIPS,
Germany;
(b) Faculty of Mathematics and Computer Science, University of Bremen,
Germany;
(c) Department of Public Health, University of Copenhagen, Denmark

## Software location:

Archive

Name:

arfpy

Persistent identifier:

https://pypi.org/project/arfpy/

Licence:

MIT

Publisher:

Kristin Blesch

Version published:

v0.1.1

Date published:

22/09/2023

Code repository

Name:

arfpy

Persistent identifier:

https://github.com/bips-hb/arfpy

Licence:

MIT

Date published:

06/09/2023

## Language

English

## (3) Reuse potential

ARFs have been introduced with a solid theoretical background, yet do not have
to compromise on a complex algorithmic structure and instead are a low-key
algorithm that does not require extensive hyperparameter tuning [7]. In
contrast to the typically deep learning based alternatives, ARF does not
require background knowledge of generative modeling, intense tuning efforts or
large computational resources. Given the theoretical foundation and
straightforward implementation with arfpy, the methodology is attractive for
both scholars conducting rather theoretical research in statistics, e.g.,
density estimation, as well as practitioners from other fields that need to
generate new data samples.
Typical use cases of such synthesized data samples are, for example, the
imputation of missing values, data augmentation or the conduct of analyses
that respect data protection rules. With the specialty of ARFs being
particularly suitable for tabular data, including a natural incorporation of
both continuous and categorical features, the straightforward python
implementation of ARFs provides a convenient algorithm to a broad audience
from different fields.
With the python programming language being widespread, arfpy can smoothly
integrate in the code of various applications. Further, usability is enhanced
by the intuitive documentation provided at https://bips-hb.github.io/arfpy/,
making arfpy an easily accessible tool to generate data.
In sum, arfpy introduces density estimation and generative modeling with ARFs
to python, which enables practitioners from a wide variety of fields to
generate fast and reliable synthetic data and density estimates with python as
a programming language.

## Acknowledgements

We thank David S. Watson and Jan Kapar for their contributions to establishing
the theoretical groundwork of adversarial random forests.

## Funding statement

This work was supported by the German Research Foundation (DFG), Emmy Noether
Grant 437611051.

## Competing interests

The authors declare that they have no competing interests.

## References

  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  _Advances in Neural Information Processing Systems_ , 27, 2014.
  * Xu et al. [2019] Lei Xu, Maria Skoularidou, Alfredo Cuesta-Infante, and Kalyan Veeramachaneni.  Modeling tabular data using conditional GAN.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Kingma and Welling [2014] DP Kingma and M Welling.  Auto-encoding variational bayes.  _International Conference on Learning Representations_ , 2014.
  * Rezende and Mohamed [2015] Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  _Proceedings of the 32nd International Conference on Machine Learning_ , 37:1530–1538, 07–09 Jul 2015.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _Advances in Neural Information Processing Systems_ , 33:6840–6851, 2020.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Watson et al. [2023] David S Watson, Kristin Blesch, Jan Kapar, and Marvin N Wright.  Adversarial random forests for density estimation and generative modeling.  _International Conference on Artificial Intelligence and Statistics. PMLR_ , 206:5357–5375, 2023.
  * Breiman [2001] Leo Breiman.  Random forests.  _Machine learning_ , 45:5–32, 2001.  10.1023/A:1010933404324.
  * Grinsztajn et al. [2022] Léo Grinsztajn, Edouard Oyallon, and Gaël Varoquaux.  Why do tree-based models still outperform deep learning on typical tabular data?  _Advances in Neural Information Processing Systems_ , 35:507–520, 2022.
  * Nock and Guillame-Bert [2023] Richard Nock and Mathieu Guillame-Bert.  Generative forests.  _arXiv Preprint 2308.03648_ , 2023.  10.48550/arXiv.2308.03648.
  * Wright and Watson [2023] Marvin N. Wright and David S. Watson.  _arf: Adversarial Random Forests_ , 2023.  URL https://CRAN.R-project.org/package=arf.  R package version 0.1.3.
  * Laurent Gautier et al. [2023] Laurent Gautier et al.  rpy2: Python-r bridge.  _GitHub repository_ , 2023.  URL https://github.com/rpy2/rpy2.
  * van der Schaar Lab [2023] van der Schaar Lab.  synthcity: A library for generating and evaluating synthetic tabular data.  _GitHub repository_ , 2023.  URL https://github.com/vanderschaarlab/synthcity.
  * Shi and Horvath [2006] Tao Shi and Steve Horvath.  Unsupervised learning with random forest predictors.  _Journal of Computational and Graphical Statistics_ , 15(1):118–138, 2006.  10.1198/106186006X94072.


# Nanothermodynamic description and molecular simulation of a single-phase
fluid in a slit pore

Olav Galteland  PoreLab, Department of Chemistry, Norwegian University of
Science and Technology<EMAIL_ADDRESS>Dick Bedeaux  PoreLab,
Department of Chemistry, Norwegian University of Science and Technology  Signe
Kjelstrup  PoreLab, Department of Chemistry, Norwegian University of Science
and Technology

###### Abstract

We describe the thermodynamic state of a highly confined single-phase and
single-component fluid in a slit pore using Hill’s thermodynamics of small
systems. This theory was more recently named nanothermodynamics. We start by
constructing an ensemble of slit pores for controlled temperature, volume,
surface area, and chemical potential. We present the integral and differential
properties according to Hill, and use them to define the disjoining pressure.
We identify all thermodynamic pressures by their mechanical counterparts in a
consistent manner, and investigate the identification by molecular dynamics
simulations. We define and compute the disjoining pressure, and show that it
contains the standard definition. We compute the entropy and energy densities,
and find in agreement with the literature, that the forces at the wall are of
an energetic, not entropic nature. The subdivision potential is zero for this
slit pore with large walls, but unequal to zero for related sets of control
variables. We show how Hill’s method can be used to find new Maxwell relations
of a confined fluid, in addition to a scaling relation, which applies when the
walls are separated far enough. By this expansion of nanothermodynamics, we
set the stage for further developments of the thermodynamics of confined
fluids, a field that is central in nanotechnology.

##  1 Introduction

The thermodynamic state of a fluid in confinement is important for the
understanding of adsorption to walls, chemical reactions, film formation and
transport in porous media [1, 2, 3, 4, 5, 6]. The molecular structuring at the
walls and the forces between particles and walls are central. The change in
thermodynamic properties upon confinement is substantial. This has been known
for a long time [7, 1]. Derjaguin considered the measurable force that
attracts or repels two walls that are close together, and defined from this
the disjoining pressure [7]. The disjoining pressure has also been called the
solvation pressure [1]. When the walls are far apart the disjoining pressure
vanishes. It is not well known, however, how size and shape as variables
affect the disjoining pressure or other properties of the confined fluid.

Confinement is considered to be important, for instance in the context of CO2
separation and sequestration by metal-organic frameworks [2] or for adsorption
in zeolites [3]. The disjoining pressure is of interest when studying
aggregation of colloidal particles, suspended or adsorbed [4, 5, 6]. It is
likely to be important also for film flow on the macroscale [8].

More knowledge of confined fluids on the nanoscale is therefore needed. It
may, for instance, help us solve the well-known up-scaling problem in porous
media science [9]. The central problem is to understand how to integrate
properties on the pore scale to the macroscale where Darcy’s law applies. In
order to account for shape and size effects it was recently proposed to use
the four Minkowski functionals [10, 11, 12]. This simplifies the description
of a representative elementary volume (REV). Another procedure using the grand
potential for average variables in the REV was also proposed [13, 14].

In this work we want to further examine this procedure [13, 14], by looking
for a way to describe the confined fluid in a pore. We are looking for a way
to deal with size- and shape-dependent variables in a systematic and general
manner. Two thermodynamic approaches are common. The approach following Gibbs
is most popular [15, 16, 17, 18, 19]. But the method of Hill may offer an
attractive alternative [20, 21, 22, 23, 24, 25], partly because it may provide
an independent check on Gibbs procedure, but also because general geometric
scaling relations are obtainable from Hill’s method. We shall see that this is
also the case in the present study.

We will pursue the method of Hill. This starts with the observation that a
small system has a surface energy comparable to its bulk energy. A consequence
is that the properties are not Euler homogeneous. Hill proposed to deal with
such systems in an original manner [26, 27]. His idea was to introduce an
ensemble of replicas of the small systems, on which standard thermodynamics
could be applied. Hill’s method has not gained much attention, in spite of a
renewed effort to spur interests [28, 29].

The long-range aim of this work is to contribute to the effort of finding
variables that characterize the confined fluid, for instance in a REV. The
grand potential offers one option to describe the pressure and other variables
[13, 14]. We will pursue this route and study a single-phase and single-
component fluid in a slit pore using Hill’s method. The so-called integral and
differential pressures introduced by Hill are central. Hill did not consider
the disjoining pressure, however it is a small system property and we believe
it is a concept that could benefit also from insights of nanothermodynamics.
The purpose of this paper is to clarify the use of Hill’s nanothermodynamics,
by applying the method to a fluid in a slit pore with walls of large surface
areas. This is a well studied case in the literature [1, 30], and is well
suited to bring out the properties of the method and its results. Grand
canonical Monte Carlo and molecular dynamics simulation techniques are well
suited to investigate thermodynamic relations. We will simulate a single-phase
fluid in a slit pore in the grand canonical ensemble. A particular advantage
of this technique is that the simulations offer a mechanical picture of the
system.

Solid-fluid and fluid-fluid interactions are considered here, but not solid-
solid interactions. The solid-solid interactions will have a large effect on
the thermodynamic state of the system at very small slit pore heights. We do
not consider quantum effects that follow from system smallness in this work.
We will consider slit pores of height $h/\lambda_{B}=6$ to $71$, where
$\lambda_{B}$ is the thermal de Broglie wavelength.

Section 2.1 introduces the reader to Hill’s nanothermodynamics. We apply this
theory in section 2.2 to define size- and shape-dependent properties, and new
Maxwell relations that follow from these. A definition of the disjoining
pressure follows naturally in section 2.3. In order to be able to verify
relations with molecular simulations as a tool, we need to identify the
integral pressures and surface tension in terms the mechanical pressure tensor
components. This is done in section 2.4. In section 3 we describe the
molecular simulations.

We proceed in section 4 to investigate relations in the theory, and illustrate
them with numerical results. We compute the local mechanical and thermodynamic
variables according to Hill; i.e. the integral and differential pressures, and
the integral and differential surface tensions. The grand potential, or the
replica energy, is equal to minus the integral pressure times the volume. The
set of thermodynamic variables of the nanothermodynamic framework, in terms of
mechanical properties, is found to be self-consistent. We offer concluding
remarks in section 5.

##  2 Theory

###  2.1 Hill’s nanothermodynamics

Consider an ensemble of $\mathcal{N}$ slit pores, where each slit pore is
filled with a single-phase and single-component fluid. The slit pores do not
interact with each other. The $j^{\text{th}}$ slit pore has two parallel plane
walls of area $\Omega_{j}$, separated by a distance $h_{j}$. The ensemble of
slit pores has the total internal energy $U_{t}$, total entropy $S_{t}$, total
volume $V_{t}=\sum_{j=1}^{\mathcal{N}}h_{j}\Omega_{j}$, total surface area
$2\Omega_{t}=2\sum_{j=1}^{\mathcal{N}}\Omega_{j}$, and total number of
particles $N$. The factor two in the total surface area arises because there
are two fluid-solid surfaces of equal area per slit pore. By construction the
ensemble variables $U_{t},S_{t},V_{t},\Omega_{t},N$, and $\mathcal{N}$ are
Euler homogeneous functions of the first order in the number $\mathcal{N}$ of
slit pores. The total differential of the total internal energy is

$\text{d}U_{t}=T\text{d}S_{t}-p_{\perp}\text{d}V_{t}+2\gamma\text{d}\Omega_{t}+\mu\text{d}N_{t}+\varepsilon\text{d}\mathcal{N}.$
(1)

This type of equation for the total internal energy we call the Hill-Gibbs
equation [25]. The last term was added by Hill. The partial derivatives of the
total internal energy define the temperature $T$, the normal pressure
$p_{\perp}$, the surface tension $\gamma$ and the chemical potential $\mu$

$\begin{split}\left(\frac{\partial U_{t}}{\partial
S_{t}}\right)_{V_{t},\Omega_{t},N_{t},\mathcal{N}}&=T,\qquad\left(\frac{\partial
U_{t}}{\partial
V_{t}}\right)_{S_{t},\Omega_{t},N_{t},\mathcal{N}}=-p_{\perp},\\\
\left(\frac{\partial
U_{t}}{\partial\Omega_{t}}\right)_{S_{t},V_{t},N_{t},\mathcal{N}}&=2\gamma,\qquad\left(\frac{\partial
U_{t}}{\partial
N_{t}}\right)_{S_{t},V_{t},\Omega_{t},\mathcal{N}}=\mu.\end{split}$ (2)

The control variables in subscripts are kept constant while taking the
derivatives. The volume derivative is taken while keeping the total surface
area constant, which implies that the volume is changed by changing the
distance between the surfaces $h_{j}\equiv V_{j}/\Omega_{j}$. The surface
derivative is taken while keeping the volume constant, which implies that both
the pore heights and surface areas are changed in such a way that the change
in the total volume is zero.

The new thermodynamic variable $\varepsilon$ is the subdivision potential. It
is defined by

$\left(\frac{\partial
U_{t}}{\partial\mathcal{N}}\right)_{S_{t},V_{t},\Omega_{t},N_{t}}=\varepsilon.$
(3)

The subdivision potential is defined here as the increase in the total
internal energy as the number of slit pores $\mathcal{N}$ increases while
keeping $S_{t},V_{t},\Omega_{t}$, and $N_{t}$ constant. This definition is
different from the definition in a previous article by us [23], where only the
entropy, volume and number of particles were kept constant and not the surface
area. This led to a different expression for the subdivision potential.

The subdivision potential is the work done on the system when adding a new
slit pore while keeping the other control variables constant. The subdivision
potential may be positive or negative, depending on whether work is needed or
gained by adding new slit pore replicas.

We are aiming to describe, see section 2.3, the disjoining pressure of an open
slit pore when the volume, surface area, temperature and chemical potential
are control variables. This set of variables is useful for describing
experiments and simulations. We use the average volume per slit pore
$V=V_{t}/\mathcal{N}$ and the average surface area per slit pore
$2\Omega=2\Omega_{t}/\mathcal{N}$, rather than the total volume $V_{t}$ and
the total surface area $2\Omega_{t}$. In order to obtain an appropriate Hill-
Gibbs equation for this case, we substitute the total volume with the average
volume $V_{t}=V\mathcal{N}$ and total surface area with the average surface
area $2\Omega_{t}=2\Omega\mathcal{N}$. The total differentials of the volume
and surface area are

$\text{d}(V\mathcal{N})=\mathcal{N}\text{d}V+V\text{d}\mathcal{N}\quad\text{and}\quad\text{d}(\Omega\mathcal{N)}=\mathcal{N}\text{d}\Omega+\Omega\text{d}\mathcal{N}.$
(4)

By introducing this into the Hill-Gibbs equation, see equation 1, we obtain

$\text{d}U_{t}=T\text{d}S_{t}-p_{\perp}\mathcal{N}\text{d}V+2\gamma\mathcal{N}\text{d}\Omega+\mu\text{d}N_{t}+(\varepsilon-
p_{\perp}V+2\gamma\Omega)\text{d}\mathcal{N}.$ (5)

The parenthesis define the replica energy

$X(T,V,\Omega,\mu)\equiv\varepsilon-p_{\perp}V+2\gamma\Omega.$ (6)

The subdivision potential $\varepsilon$, normal pressure $p_{\perp}$, and
surface tension $\gamma$ depend on the control variable set $T,V,\Omega$, and
$\mu$. The replica energy density will be used to define the disjoining
pressure.

The grand partition function covers all microstates available to the slit
pore. In order to calculate this partition function one chooses a volume $V$
and a surface area $\Omega$. Both can be varied independently. For large $h$
and $\Omega$ only a change in the volume $V=h\Omega$, and not in $h$ and
$\Omega$ separately, matters. For a small volume it is necessary to use the
volume $V$ and the surface area $\Omega$ as independent variables. This is
because the same volume change due to a change of the pore height $h=V/\Omega$
or due to a change of the surface area $\Omega$ produces different changes in
the partition function, and therefore in the thermodynamic variables.

The replica energy was identified by Hill as the grand potential, here equal
to minus the integral pressure times the volume [26],

$X(T,V,\Omega,\mu)=-\hat{p}V=-\hat{p}_{\perp}V+2\hat{\gamma}\Omega.$ (7)

In the last equality, we have chosen to identify the integral pressure as the
integral normal pressure times the volume minus the integral surface tension
times the surface area. We will refer to $\hat{p}_{\perp}$, $\hat{\gamma}$ and
$\hat{p}$ as the integral normal pressure, the integral surface tension and
the integral pressure, respectively, while $p_{\perp}$ and $\gamma$ are the
differential normal pressure and the differential surface tension,
respectively. The names integral and differential pressure where coined by
Hill to reflect that the differential pressure involves the differential of
the integral pressure. We have chosen a control variable set with volume and
surface area, such that we do not have a differential pressure but a
differential normal pressure $p_{\perp}$ and differential surface tension
$\gamma$ in its place.

From equations 6 and 7 it follows that the subdivision potential is

$\varepsilon=\left(p_{\perp}-\hat{p}_{\perp}\right)V-2\left(\gamma-\hat{\gamma}\right)\Omega.$
(8)

The subdivision potential $\varepsilon$ indicate that the integral normal
pressure and the integral surface tension may be different from the
corresponding differential variables.

###  2.2 Maxwell relations for a slit pore

Using that the total internal energy is Euler homogeneous of the first order
in the number of slit pores, see equation 5, it follows that the total
internal energy is

$U_{t}=TS_{t}+\mu N_{t}+X\mathcal{N}.$ (9)

We introduce the average internal energy, entropy and number of particles per
slit pore

$U_{t}=U\mathcal{N},\qquad S_{t}=S\mathcal{N},\quad\text{and}\quad
N_{t}=N\mathcal{N}.$ (10)

By introducing the average properties into the total internal energy in
equation 9 we obtain the internal energy per slit pore

$U=TS+\mu N+(\varepsilon-p_{\perp}V+2\gamma\Omega)=TS+\mu N+X.$ (11)

Substituting the averages properties per slit pore in equation 10 into the
corresponding Hill-Gibbs equation, see equation 5, and using the internal
energy in equation 11, we obtain the total differential of the internal energy

$\text{d}U=T\text{d}S-p_{\perp}\text{d}V+2\gamma\text{d}\Omega+\mu\text{d}N.$
(12)

By differentiating the internal energy in equation 11 and using the total
differential of the internal energy in equation 12, we obtain the total
differential of the replica energy

$\text{d}X=-\text{d}(\hat{p}_{\perp}V-2\hat{\gamma}\Omega)=-\text{d}(\hat{p}V)=-S\text{d}T-p_{\perp}\text{d}V+2\gamma\text{d}\Omega-N\text{d}\mu.$
(13)

This equation was termed the Hill-Gibbs-Duhem equation [25], because it
reduces to the Gibbs-Duhem equation for a large system. It follows that the
partial derivatives of the grand potential is

$\begin{split}\left(\frac{\partial\left(\hat{p}V\right)}{\partial
T}\right)_{V,\Omega,\mu}&=-\left(\frac{\partial X}{\partial
T}\right)_{V,\Omega,\mu}=S,\\\
\left(\frac{\partial\left(\hat{p}V\right)}{\partial
V}\right)_{T,\Omega,\mu}&=-\left(\frac{\partial X}{\partial
V}\right)_{T,\Omega,\mu}=p_{\perp},\\\
\left(\frac{\partial\left(\hat{p}V\right)}{\partial\Omega}\right)_{T,V,\mu}&=-\left(\frac{\partial
X}{\partial\Omega}\right)_{T,V,\mu}=-2\gamma,\\\
\left(\frac{\partial\left(\hat{p}V\right)}{\partial\mu}\right)_{T,V,\Omega}&=-\left(\frac{\partial
X}{\partial\mu}\right)_{T,V,\Omega}=N.\end{split}$ (14)

Rather than the names replica energy or grand potential, we name from now
$X=-\hat{p}V$ by minus the integral pressure times the volume. In section 2.4
we will define the integral pressure in terms of the average tangential
mechanical pressure which can be calculated from molecular simulations. In the
simulations we consider surface areas $\Omega$ much larger than the diameter
of the fluid particles and heights $h=V/\Omega$ comparable to the diameter of
the fluid particles. This implies that $\hat{p}$, $\hat{p}_{\perp}$,
$\hat{\gamma}$, $u\equiv U/V$, $s\equiv S/V$, $p_{\perp}$, $\gamma$, and
$\rho\equiv N/V$ do not depend on $\Omega$, but the variables will depend on
the height $h$ and therefore on the volume $V$.

The volume and surface derivative can be rewritten in terms of derivatives of
the slit pore height and surface area. These derivatives are needed to
calculate the differential normal pressure and differential surface tension.
The differential normal pressure is

$p_{\perp}=\left(\frac{\partial(\hat{p}V)}{\partial
V}\right)_{T,\Omega,\mu}=\hat{p}+V\left(\frac{\partial\hat{p}}{\partial
V}\right)_{T,\Omega,\mu}=\hat{p}+h\left(\frac{\partial\hat{p}}{\partial
h}\right)_{T,\Omega,\mu}.$ (15)

The differential surface tension is

$\gamma=-\frac{V}{2}\left(\frac{\partial\hat{p}}{\partial\Omega}\right)_{T,V,\mu}=\frac{h^{2}}{2}\left(\frac{\partial\hat{p}}{\partial
h}\right)_{T,\Omega,\mu}-\frac{V}{2}\left(\frac{\partial\hat{p}}{\partial\Omega}\right)_{T,h,\mu}.$
(16)

Combining the equations for the differential normal pressure and surface
tension it follows that the integral pressure is

$\hat{p}=p_{\perp}-\frac{2}{h}\gamma-\Omega\left(\frac{\partial\hat{p}}{\partial\Omega}\right)_{T,\mu,h}.$
(17)

The last term is the equal to minus the subdivision potential divided by
volume, see equations 7 and 8. The subdivision potential is

$\varepsilon=\Omega
V\left(\frac{\partial\hat{p}}{\partial\Omega}\right)_{T,\mu,h}.$ (18)

This implies that in the grand canonical ensemble with $T,V,\Omega,\mu$ as
control variables, the thermodynamic description of the slit pore with a small
height is the same as for the slit pore with a large height. However, it
changes when the integral pressure depends on the surface area. This may not
be the case for any other set of control variables. In general the properties
of a small system (confined fluid) depend on the set of control variables.

In this work we deal with large surface areas, such that the integral pressure
does not depend on it. As a consequence the subdivision potential is zero. We
will find that the integral and differential normal pressure are equal. Using
that $\varepsilon=0$ it follows that the integral and differential surface
tensions are also equal for a large surface area, see equation 8,

$\hat{p}_{\perp}=p_{\perp}\quad\text{and}\quad\hat{\gamma}=\gamma.$ (19)

We shall investigate the first identity for large surface areas using
molecular simulations.

From equation 13 we obtain the following Maxwell relations of the differential
surface tension,

$\begin{split}\left(\frac{\partial\gamma}{\partial
T}\right)_{V,\Omega,\mu}&=-\frac{1}{2}\left(\frac{\partial
S}{\partial\Omega}\right)_{T,\mu,V}=\frac{h^{2}}{2}\left(\frac{\partial
s}{\partial h}\right)_{T,\mu,\Omega},\\\
\left(\frac{\partial\gamma}{\partial\mu}\right)_{T,V,\Omega}&=-\frac{1}{2}\left(\frac{\partial
N}{\partial\Omega}\right)_{T,\mu,V}=\frac{h^{2}}{2}\left(\frac{\partial\rho}{\partial
h}\right)_{T,\mu,\Omega},\\\ \left(\frac{\partial\gamma}{\partial
V}\right)_{T,\Omega,\mu}&=-\frac{1}{2}\left(\frac{\partial
p_{\perp}}{\partial\Omega}\right)_{T,\mu,V}=\frac{h}{2\Omega}\left(\frac{\partial
p_{\perp}}{\partial h}\right)_{T,\mu,\Omega}.\end{split}$ (20)

In the second identity we have used that the entropy density, fluid number
density and differential normal pressure are independent of the surface area
for large surface areas. Other Maxwell relations are possible, see the Hill-
Gibbs-Duhem equation 13. The last equality can be written as

$\left(\frac{\partial\gamma}{\partial
h}\right)_{T,\mu}=\frac{h}{2}\left(\frac{\partial p_{\perp}}{\partial
h}\right)_{T,\mu,\Omega}.$ (21)

The important implication of this expression is that the constant nature of
one variable implies the constant nature of the other for large pore heights.
With the mechanical description of the integral properties given in section
2.4 all the above thermodynamic relations can be tested.

###  2.3 The disjoining pressure

We define the excess replica energy as the replica energy of a slit pore of
height $h$ minus the replica energy of the slit pore where the slit pore
height approaches infinity,

$X^{\text{ex}}\ =X-\lim_{h\rightarrow\infty}X.$ (22)

We denote the thermodynamic variables where the slit pore height approaches
infinity with $\infty$ in superscript,

$\lim_{h\rightarrow\infty}X=X^{\infty}.$ (23)

The excess replica energy can be written in terms of the excess integral
pressure $\hat{p}^{\text{ex}}$, excess integral normal pressure
$\hat{p}_{\perp}^{\text{ex}}$, and excess integral surface tensions
$\hat{\gamma}^{\text{ex}}$

$X^{\text{ex}}\ =-\hat{p}^{\text{ex}}\
V=-\hat{p}_{\perp}^{\text{ex}}V+2\hat{\gamma}^{\text{ex}}\Omega.$ (24)

We define now the disjoining pressure as the excess normal pressure,

$\Pi(h)\equiv\hat{p}_{\perp}^{\text{ex}}\equiv\hat{p}_{\perp}-p_{\perp}^{\infty}.$
(25)

Other possible names are minus the excess replica density, or the excess grand
potential density.

This definition of the disjoining pressure is different from the typical
definition found in the literature [7, 30, 31]. Typically the disjoining
pressure is defined to be equal to what we in this work call the excess
differential normal pressure, where the excess is relative to a bulk fluid.
For a large surface area our molecular dynamics simulations find that
$\hat{p}_{\perp}=p_{\perp}$. We will furthermore show that the integral normal
pressure, as the slit pore height approaches infinite separation, equals the
bulk fluid pressure $p_{\perp}^{\infty}=p^{b}$. This shows that our choice is
equivalent to the usual definition.

###  2.4 A mechanical description of the slit pore

The thermodynamic description of integral and differential pressures and
surface tensions has its mechanical equivalent description in terms of
components of the mechanical pressure tensor. A recent discussion on this
topic clarified the challenge of translating the mechanical pressure tensor
into a thermodynamic scalar variable in a meaningful manner [32, 33, 34]. We
will identify the integral normal pressure, integral surface tension, and
integral pressure in such a way that thermodynamic framework is self-
consistent. However, we do not claim that this is the only valid choices of
thermodynamic pressures and tensions in terms of the mechanical pressure
tensor.

The mechanical pressure tensor of a heterogeneous system is ambiguously
defined. This has been known for a long time, at least since the work by
Irving and Kirkwood in the 1940s. It was shown by Schofield and Henderson [35]
that the ambiguity is due to the arbitrary choice of the integration contour
$C_{ij}$, which is needed to calculate the configurational contribution to the
pressure tensor. The local mechanical pressure tensor is calculated in a
subvolume $V_{l}$ as a sum of the kinetic and the configurational
contributions,

$P_{\alpha\beta}(x)=P_{\alpha\beta}^{k}+P_{\alpha\beta}^{c}.$ (26)

Upper case $P$ is used to denote mechanical pressure tensor components in
order to distinguish them from the thermodynamic pressures, which are denoted
by lower case $p$. The kinetic contribution is the ideal contribution to the
mechanical pressure and is calculated as

$P_{\alpha\beta}^{k}=\frac{1}{V_{l}}\left\langle\sum_{i\in
V_{l}}m_{i}v_{i,\alpha}v_{i,\beta}\right\rangle.$ (27)

The sum with subscript $i\in V_{l}$ represents a sum over all particles in the
subvolume $V_{l}$. The particle mass is $m_{i}$ and $v_{i,\alpha}$ is the
velocity in the $\alpha$-direction. The solid walls do not have a velocity and
consequently do not contribute to the kinetic pressure. The brackets
$\langle\dots\rangle$ represent the grand canonical ensemble average. The
configurational contribution is the non-ideal contribution to the mechanical
pressure and is calculated as

$P_{\alpha\beta}^{c}=-\frac{1}{2V_{l}}\left\langle\sum_{i=1}^{N}\sum_{\begin{subarray}{c}j=1\\\
j\neq i\end{subarray}}^{N}f_{ij,\alpha}\int_{C_{ij}\in
V_{l}}\text{d}l_{\beta}\right\rangle.$ (28)

The sums represents a sum over all particle pairs. The $\alpha$-component of
the force vector acting on particle $i$ due to particle $j$ is
$f_{ij,\alpha}$. The fluid-fluid and fluid-solid interactions contribute to
the configurational pressure, the solid-solid interaction is zero and does not
contribute to the pressure. The line integral is the $\beta$-component of the
part of the contour $C_{ij}$ contained in the subvolume $V_{l}$.

The contour $C_{ij}$ is the source of the ambiguity of the mechanical pressure
tensor, it can be any continuous line from the centers of particles $i$ to
$j$. The Harasima [36] and the Irving-Kirkwood [37] contours are two common
choices for $C_{ij}$. The Harasima contour is defined as two continuous line
segments, a line from the center of particle $i$ parallel to the surface and a
line normal to the surface to the center of particle $j$. The Irving-Kirkwood
is the straight line from particle $i$ to $j$. For flat surfaces they are
equal. However for spherical surfaces the Harasima contour does not obey
momentum balance [38]. There are cases where the Harasima contour is more
useful than the Irving-Kirkwood contour [39]. In this work we will use the
equations by Ikeshoji et al. [40] with the Irving-Kirkwoood contour to
calculate the mechanical pressure tensor.

We will only consider surface areas $\Omega$ much larger in both directions
than the fluid particle diameter. This implies the mechanical pressure tensor
does not depend on $\Omega$. It does however depend on the height $h\equiv
V/\Omega$ and therefore on the volume $V$.

Using the translational symmetry of the slit pore in the $y$\- and
$z$-direction, where the $x$-direction is normal to the solid surface, the
equilibrium mechanical pressure tensor in the slit pore has the form

$\mathbf{P}(x;h)=P_{\perp}(h)\mathbf{e}_{x}\mathbf{e}_{x}+P_{\parallel}(x;h)(\mathbf{e}_{y}\mathbf{e}_{y}+\mathbf{e}_{z}\mathbf{e}_{z}).$
(29)

Where $\mathbf{e}_{x}$, $\mathbf{e}_{y}$ and $\mathbf{e}_{z}$ are the unit
vectors in the $x$-, $y$\- and $z$-directions. The normal pressure tensor
component is equal to the $xx$-component, and tangential pressure tensor
component is the average of the $yy$\- and $zz$-components,

$P_{\perp}(h)=P_{xx}\quad\text{and}\quad
P_{\parallel}(x,h)=\frac{1}{2}\left(P_{yy}+P_{zz}\right).$ (30)

Mechanical equilibrium requires that the tangential pressure is independent of
the $y$\- and $z$-coordinates, but depends on the $x$-coordinate. The normal
pressure is independent of all spatial coordinates.

We identify the thermodynamic integral normal pressure in terms of the volume
integral of the normal mechanical pressure divided by the volume. However,
since the normal mechanical pressure and the area are constant everywhere,
this simplifies to

$\hat{p}_{\perp}(h)\equiv\frac{1}{h}\int_{0}^{h}P_{\perp}(h)x=P_{\perp}(h).$
(31)

The integral normal pressure in a large pore is equal to the pressure in a
bulk phase in equilibrium with the pore. The difference in the integral normal
pressure in a liquid and vapor phase in a slit pore is described the Young-
Laplace equation, while the integral pressure is the same in both the liquid
and vapor phase [23]. An alternative route to the integral normal pressure is
via the local fluid number density profile [41]

$P_{\perp}(h)=\int_{0}^{h}f_{fs}(x)\rho(x,h)x,$ (32)

where $f_{fs}(x)$ is the fluid-solid force and $\rho(x,h)$ is the local fluid
number density. The fluid density and fluid-solid force are uniquely defined,
and do not have the inherent problem that the mechanical pressure tensor has.
We can use this equation to validate our method of calculating the integral
normal pressure.

We identify further the integral surface tension as the integral of the normal
minus tangential pressure tensor components,

$\hat{\gamma}(h)\equiv\frac{1}{2}\int_{0}^{h}\left(P_{\perp}(h)-P_{\parallel}(x;h)\right)\text{d}x,$
(33)

where the factor half is due to the fact that there are two fluid-solid
surfaces. It follows from equation 7 together with equations 31 and 33 that
the integral pressure is,

$\hat{p}(h)=\frac{1}{h}\int_{0}^{h}P_{\parallel}(x;h)\text{d}x.$ (34)

As shown by Harasima [36] and Schofield and Henderson [35] a sufficiently
large volume integral of the mechanical pressure tensor components does not
depend on the choice of the integral contour $C_{ij}$. We have identified the
integral normal pressure, integral surface tension and integral pressure in
terms of the mechanical pressure tensor. The local mechanical pressure tensor
has an inherent problem, specifically that the contour $C_{ij}$ can be any
continuous line from $i$ to $j$. However the thermodynamic variables do not
have this inherent problem. This is because we integrate the the local
mechanical pressure tensor across the whole volume $V$. This volume integral
includes all interactions.

The internal energy can be calculated as the sum of the kinetic and potential
energy,

$U=E_{k}+E_{p}.$ (35)

By dividing the internal energy in equation 11 by the volume we obtain the
internal energy density $u=U/V$. We also use the entropy density $s=S/V$, and
fluid number density $\rho=N/V$. By rearranging the equation we obtain the
entropy density as

$s=\frac{1}{T}\left(u+\hat{p}-\mu\rho\right),$ (36)

where the internal energy density, integral pressure and fluid number density
are known. The volume, chemical potential and temperature are imposed on the
system.

##  3 Simulation details

The thermodynamic state of slit pores of varying heights $h$ was investigated
by using grand-canonical Monte Carlo (GCMC) [42] in combination with molecular
dynamic (MD) simulations with the Nosé-Hoover thermostat [43]. This produced
the grand canonical ensemble, i.e. constant chemical potential, temperature,
volume and surface area. The GCMC method inserted and deleted fluid particles
to and from the simulation box from an imaginary fluid particle reservoir at
the same temperature and chemical potential. This controlled the chemical
potential of the fluid in the slit pore. The MD procedure updated the
positions and velocities of the fluid particles and controlled the temperature
with the Nosé-Hoover thermostat.

The simulations were carried out using LAMMPS [44]. The local mechanical
pressure tensor was calculated by post-processing the particle trajectories
with in-house software. The chemical potential and temperature were kept
constant at $\mu^{\ast}=1$ and $T^{\ast}=2$. The critical temperature of the
Lennard-Jones/spline fluid is $T_{c}^{\ast}=0.855$ [45]. All units in this
work are in reduced Lennard-Jones units, see table 1 for a definition.

Table 1: The reduced units are denoted with an asterisk in superscript, for example $T^{\ast}$. The variables are reduced using the molecular diameter $\sigma$, potential well depth $\epsilon$, fluid particle mass $m$ and Boltzmann constant $k_{\text{B}}$ Description | Definition
---|---
Energy | $E^{\ast}=E/\epsilon$
Entropy | $S^{\ast}=S/k_{\text{B}}$
Temperature | $T^{\ast}=Tk_{\text{B}}/\epsilon$
Distance | $x^{\ast}=x/\sigma$
Pressure | $p^{\ast}=p\sigma^{3}/\epsilon$
Chemical potential | $\mu^{\ast}=\mu/\epsilon$

The simulation box was a rectangular cuboid of side lengths
$L_{x},L_{y}=L_{z}$. The side lengths $L_{y}=L_{z}$ were chosen such that the
surface area was large. Large in this context indicates large enough for
$\hat{p}$, $\hat{p}_{\perp}$, $\hat{\gamma}$, $p_{\perp}$, and $\gamma$ to be
independent of the surface area $\Omega=L_{y}L_{z}$. The simulation box size
was decided such that the average number of fluid particles was approximately
$2\times 10^{4}$. The simulation box was periodic in the $y$\- and
$z$-directions, and non-periodic in the $x$-direction. This implies that the
particles did not interact across the simulation box boundary in the
$x$-direction.

The fluid-fluid and fluid-solid interaction was modeled with the Lennard-
Jones/spline potential [45]. The fluid-fluid and fluid-solid interactions were
equal. The potential energy of a fluid-fluid or fluid-solid pair separated by
a distance $r$ was

$u^{LJ/s}(r)=\begin{cases}4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]&\text{if
}r<r_{s}\\\ a(r-r_{c})^{2}+b(r-r_{c})^{3}&\text{if }r_{s}<r<r_{c}\\\
0&\text{else},\end{cases}$ (37)

where $r_{s}$, $a$, $b$ and $r_{c}$ were chosen such that the potential energy
and the force were continuous at the inflection point $r=r_{s}$ and the cut-
off $r=r_{c}$. The solid walls were placed at the simulation box boundaries
$x=-L_{x}/2$ and $x=L_{x}/2$. The distance between the fluid and solids were
$r=|x_{f}-L_{x}/2|$ and $r=|x_{f}+L_{x}/2|$, where $x_{f}$ is the $x$-position
of the fluid particle.

The dividing surfaces of the fluid-solid surfaces were chosen to be at
$x=-L_{x}/2$ and $x=L_{x}/2$. The slit pore height was consequently determined
to be $h=L_{x}$. Other choices of the dividing surface are possible, for
example the Gibbs dividing surface or the surface of tension. When
$L_{x}<2r_{c}$ the fluid particle can interact with both solid walls. See
figure 1 for an visualization of the simulation box for the case
$L_{x}=4\sigma$.

Figure 1: Visualization of the fluid particles in a slit pore of height
$L_{x}=4\sigma$, chemical potential $\mu^{\ast}=1$, and temperature
$T^{\ast}=2$. The fluid particles are rendered in red, and their diameter is
rendered at $\sigma$. The solid is not rendered. The solid lines illustrates
the edges of the simulation box. The simulation was rendered with OVITO [46].

The mechanical pressure tensor was calculated in thin rectangular cuboids,
called layers $l$, of side lengths $\Delta x,L_{y},L_{z}$. The thickness of
the layers was $\Delta x=0.005\sigma$ and the number of layers was
$n_{l}=h/\Delta x$. The diagonal components of the mechanical pressure tensor
was calculated using equations 29, 27 and 28.

The kinetic energy was calculated as the sum of the kinetic energy for each
fluid particle and the potential energy was calculated as the sum of the
potential energy of each fluid-fluid and fluid-solid pairs,

$E_{k}=\frac{1}{2}\sum_{i=1}^{N}m_{i}(\mathbf{v}_{i}\cdot\mathbf{v}_{i})\quad\text{and}\quad
E_{p}=\sum_{i=1}^{N}\sum_{j>i}^{N}u^{LJ/s}(r_{ij}).$ (38)

The sum of the kinetic and potential energies were used to calculate the
internal energy and entropy densities.

##  4 Results and discussion

The results are presented in figures 2 to 10 and discussed in that order
before general remarks are offered.

(a)

(b)

(c)

Figure 2: (a) Normal mechanical pressure, (b) tangential mechanical pressure,
and (c) fluid number density $\rho$ as a function of the $x$-direction for
slit pore heights $h=2.04\sigma$, $2.59\sigma$, and $8\sigma$.

The normal mechanical pressure $P_{\perp}$ is presented in figure 2(a). It
does not depend on the position $x$, but it depends strongly on the slit pore
height $h$. The figure shows a straight line of various lengths for each of
the three heights, which reflect the slit pore height $h$. The integral normal
pressure was identified as this component, $\hat{p}_{\perp}=P_{\perp}$. We see
that it is always constant, as demanded by equation 31. For slit pore heights
$h>7\sigma$ we find that the normal mechanical pressure is equal to the bulk
pressure. At $h=2.04\sigma$ the normal mechanical pressure is at a global
maximum and at $h=2.59\sigma$ it has a local minimum for the given temperature
and chemical potential. The normal mechanical pressures divided by the bulk
pressure for the two cases are $P_{\perp}/p^{\text{b}}=2.5135\pm 0.0008$ and
$P_{\perp}/p^{\text{b}}=0.5539\pm 0.0002$, respectively. The global minimum,
which is zero, is at $h<1.8\sigma$ when no fluid particles fit in the slit
pore. This is because solid-solid interactions and quantum effects are not
considered in this work.

The tangential mechanical pressure $P_{\parallel}$, illustrated in figure
2(b), depends in contrast on the position $x$ as well as on the slit pore
height $h$. The integral pressure is the average of $P_{\parallel}$, see
equation 34. The tangential mechanical pressure follows the trend of the fluid
number density, compare figures 2(b) and 2(c). For pore sizes $h>7\sigma$ the
tangential mechanical pressure is constant and equal to the bulk pressure
$p^{\text{b}}$ in the center of the pore. This indicates that the pore is
large enough to accommodate bulk liquid in the center. The fluid is highly
structured close to the fluid-solid surface [1]. As the slit pore height is
decreased, fluid structures on the two sides overlap. When regions of
structured fluids overlap, repulsive and attractive forces between the
surfaces appear, and the disjoining pressure becomes non-zero.

The fluid number density $\rho=N/V$ is presented in figure 3 as a function of
the slit pore height $h$. The bulk fluid number density $\rho^{\text{b}}$ in
equilibrium with the slit pore is shown as a dashed line. The fluid number
density converges to the bulk value as the slit pore height approaches
infinity. The volume $V$ depends on the choice of the fluid-solid dividing
surfaces. We have chosen the dividing surfaces to be at $x=-L_{x}/2$ and
$x=L_{x}/2$. The choice of the dividing surface determines how rapidly the
slit pore values converges to the bulk values. A dividing surface closer to
the fluid phase will reduce the volume and consequently the slit pore values
will converge faster to the bulk values. Other choices of the dividing surface
are possible.

Figure 3: Fluid number density $\rho=N/V$ as a function of slit pore height
$h$. The bulk fluid number density is shown as a dashed line.

The entropy density is presented in figure 4 as a function of the slit pore
height $h$. The entropy density is a monotonically increasing function of the
height $h$. This confirms the observation by Israelachvili [1] that the origin
of the oscillations of the disjoining pressure as a function of the height is
not entropic. As a further confirmation of this point, we find that the
internal energy density oscillates with a period equal to the particle
diameter, see figure 5. The oscillating forces or pressures are thus of
energetic origin. The bulk entropy and internal energy densities are shown as
dashed lines. The entropy and internal energy densities of the slit pore
converges to the bulk values as the slit pore height is increased.

Figure 4: Entropy density $s=S/V$ as a function of slit pore height $h$. See
equation 36. The dashed line shows the entropy density of the bulk
$s^{\text{b}}$. Figure 5: Internal energy density $u=U/V$ as a function of
slit pore height $h$. See equation 35. The dashed line shows the internal
energy density of the bulk $u^{\text{b}}$.

The integral pressure is equal to the volume average tangential mechanical
pressure $\hat{p}=h^{-1}\int_{0}^{h}P_{\parallel}\text{d}x$. It is of special
interest because it is equal to minus the grand potential divided by volume
$\hat{p}=-X/V$, or the replica energy density. The grand potential is the
starting point for the definition of the average thermodynamic properties of
the REV [13]. The integral pressure is presented in figure 6 as a function of
the slit pore height $h$. The integral pressure converges to the bulk pressure
$p^{\text{b}}$ as the slit pore height approaches infinity as expected. The
bulk pressure is shown as a dashed line.

In previous works [13, 14, 21] we argued that the gradient of the integral
pressure is the driving force for mass flux. In another work [23] we found the
integral pressure of a two-phase system in a slit pore to be equal in the
liquid and vapor in equilibrium. The identification of the integral pressure
in this work is consistent with this interpretation. The gradient of the
integral pressure is the driving force of the mass flux. For fluid flows
tangential to the slit pore surfaces it is the gradient in the tangential
mechanical pressure tensor component that gives the driving force when the
system is out of equilibrium. In this work we identify the integral pressure
as the average of the tangential mechanical pressure tensor components.

As stated in section 2.2, the integral and differential normal pressures and
integral and differential surface tensions are expected to be equal when the
surface area is large. If this is correct the integral pressure can be
computed as

$\hat{p}(h)=\frac{1}{h}\int_{h_{0}}^{h}\hat{p}_{\perp}h^{\prime}\quad=\int_{h_{0}}^{h}\frac{2\hat{\gamma}}{h^{\prime
2}}h^{\prime}.$ (39)

The lower integration limit is $h_{0}=1.8\sigma$, at which point the integral
pressure is in good approximation zero. The integral pressure computed from
equation 39 is shown in figure 6. The curves are identical. This implies that
the integral and differential normal pressures are equal. As we have already
shown that the subdivision potential is zero $\varepsilon=0$ for large surface
areas with this set of control variables, it follows that the integral and
differential surface tensions are also equal. We will from now on refer to the
integral normal pressure and integral surface tensions as the normal pressure
and surface tension.

Figure 6: Integral pressure $\hat{p}$ as a function of slit pore height $h$.
The integral pressure is computed as the average tangential mechanical
pressure, see equation 34, and from the normal pressure and surface tension,
see equation 39.

The normal pressure was identified as the normal mechanical pressure
$p_{\perp}=P_{\perp}$ in equation 31. It is presented in figure 7 as a
function of the slit pore height $h$. The normal pressure was also calculated
from the local fluid density and the fluid-solid force using equation 32. The
two methods of calculating the normal pressure agree, which indicates that we
have calculated the mechanical pressure tensor correctly.

The normal pressure oscillates with a period equal to the fluid particle
diameter at small heights $h$. The oscillations decay as the height $h$
increases. Such oscillations have been observed in experiments and are well
known, see for example Israelachvili [1]. The oscillations are caused by the
structuring of the fluid particles between the surfaces, and by the fact that
the fluid particles all have the same diameter. These oscillations are not
present in a mixture of molecules with different diameters. As the height is
increased above $h>7\sigma$ the oscillations vanish and the normal pressure is
constant and equal to the bulk pressure. The bulk pressure is shown as a
dashed line in the figure. At heights $h>7\sigma$ the fluid structuring near
the walls do not overlap. At lower densities smaller oscillations are expected
with faster decay. The normal pressure shows a similar trend to previous works
[41, 15, 16]. At very small heights the solid-solid interaction will dominate
and completely overshadow the fluid-fluid and fluid-solid interactions
presented here. We have not included any solid-solid interaction in this work,
and as a consequence the normal pressure approaches zero because there is no
room for any fluid particles to enter the slit pore.

Figure 7: Normal pressure $p_{\perp}=P_{\perp}$ as a function of slit pore
height $h$. It is computed as the normal mechanical pressure tensor component,
see equation 31, and as the integral of the fluid-solid force times the local
fluid density, see equation 32.

When $p_{\perp}=\hat{p}_{\perp}$ and $\varepsilon=0$, it follows from equation
8 that the integral and differential surface tensions are equal
$\gamma=\hat{\gamma}$. The surface tension as a function of the slit pore
height is presented in figure 8, see equation 33. The surface tension at
infinite separation $\gamma^{\infty}$ is computed as the average surface
tension of the slit pore, with height $h>10\sigma$, at which point the surface
tension is independent of $h$.

Figure 8: Surface tension $\gamma$ as a function of slit pore height $h$, see
equation 33. The dashed line shows the surface tension at infinite separation
$\gamma^{\infty}$.

The disjoining pressure $\Pi$ was computed from equation 25, and is shown as a
function of the slit pore height $h$ in figure 9. The disjoining pressure was
here defined to be equal to the integral normal pressure minus the normal
pressure at infinite separation. Because the integral and the differential
normal pressures are the same and because the normal pressure at infinite
separation is equal to the bulk pressure in this case, the definition contain
the commonly used definition of the disjoining pressure [7, 30]. The normal
pressure is constant for slit pore heights $h>7\sigma$. Consequently, the
normal pressure at infinite separation can be calculated as the normal
pressure when $h>7\sigma$. The normal pressure at infinite separation is equal
to the pressure in a bulk fluid with the same temperature and chemical
potential.

Figure 9: Disjoining pressure as a function the slit pore height $h$, see
equation 25.

For the present case, we can claim that our definition of the disjoining
pressure is equivalent to the common definition. Our definition is general, as
it also covers the cases where the integral and the differential normal
pressures are unequal. Examples where this is the case are given below.

Figure 10 shows how the normal pressure minus the integral pressure scales
with inverse slit pore height. It is equal to the scaling law presented in a
previous work [23]. The slope of this ideal curve is equal to two times the
surface tension. When the inverse slit pore height approaches zero, i.e. when
the walls are far apart, the normal pressure and integral pressures are equal
as predicted. For pores that are so small that no bulk fluid can form in the
center, the structuring at the walls starts to overlap, and a fluctuating
difference is seen in the difference of the normal and integral pressures. In
the region of the straight line, we have a scaling law, that relates states of
different heights. Heights smaller than approximately $h<5\sigma$ the scaling
law breaks down, the difference of the two pressures starts to oscillate.
There are positive and negative deviations from the law.

Figure 10: The scaling of normal pressure minus integral pressure as a
function of the inverse slit pore height $h$.

In our earlier work [23], we studied liquid-vapor coexistence in a slit pore.
In that work the surface area was not a control variable, and as a consequence
the subdivision potential was found to be equal to two times the surface
tension divided by the slit pore height, see figure 10. Another reasonable set
of control variables is the height $h$ instead of the volume $V$. For a
control variable set consisting of temperature, height, surface area and
chemical potential, the differential normal pressure is equal to

$p_{\perp}=\frac{1}{\Omega}\left(\frac{\partial(\hat{p}V)}{\partial
h}\right)_{T,\Omega,\mu}=\hat{p}+h\left(\frac{\partial\hat{p}}{\partial
h}\right)_{T,\Omega,\mu}.$ (40)

The differential surface tensions is

$\gamma=-\frac{1}{2}\left(\frac{\partial(\hat{p}V)}{\partial\Omega}\right)_{T,\mu,h}=-\frac{\hat{p}h}{2}-\frac{V}{2}\left(\frac{\partial\hat{p}}{\partial\Omega}\right)_{T,\mu,h}.$
(41)

For this set of control variables the height $h$ is kept constant instead of
the volume $V$. The subdivision potential is accordingly

$\begin{split}\varepsilon&=-\hat{p}V+p_{\perp}V+2\gamma\Omega\\\
&=\hat{p}V+hV\left(\frac{\partial\hat{p}}{\partial
h}\right)_{T,\Omega,\mu}+\Omega
V\left(\frac{\partial\hat{p}}{\partial\Omega}\right)_{T,h,\mu}.\end{split}$
(42)

These relations also help us characterize the smallness of the slit pore with
the large walls. The subdivision potential, introduced by Hill as a measure of
smallness, deviates from zero in the last relation, also when the integral
pressure does not depend on the size of the area, with height $h$ and surface
area $\Omega$ as control variables. A dependency on area is relevant when
adsorption takes place on small spheres [24]. A slit pore with large walls,
may be expected to be small for small heights $h$, since the confined fluid is
not bulk-like. However, we have seen that when we use volume $V$ and surface
area $\Omega$ as control variables, the subdivision potential is zero for
large surface areas. A zero subdivision potential means that the system also
can be described perfectly using regular thermodynamics [18]. It is
nevertheless meaningful to define a non-zero integral pressure, because the
integral pressure enters the grand potential. It therefore determines the
thermodynamic properties of a REV. Away from equilibrium, it will create a
driving force.

The grand potential or minus the integral pressure times the volume are equal
to the replica energy. The replica energy is not zero in the present case.
Clearly, we have here an example where smallness is brought out in Hill’s
terms through the replica potential.

##  5 Conclusion

We have developed a nanothermodynamic description based on the ideas of Hill
to describe slit pores. As environmental control variables we chose the
chemical potential, temperature, volume and surface area. We have seen that
the outcome varies with the set chosen, but the procedure can be used for
complex geometries and different sets.

Following Hill, we introduced the subdivision potential. It is non-zero only
when the integral pressure depends on the surface area $\Omega$. For large
surface areas we find the subdivision potential is zero, and that
$p_{\perp}=\hat{p}_{\perp}$ and $\gamma=\hat{\gamma}$. In this sense
nanothermodynamics is equivalent to the usual thermodynamic description for
all heights. The replica energy and therefore the integral pressure are non-
zero. This allows us to identify a scaling law, which confirms earlier results
[23, 22]. By choosing height $h$ rather than volume $V$ among the control
variables, a non-zero subdivision potential appears.

We have identified the thermodynamic properties by their mechanical
counterparts in a consistent manner. The integral pressure, which is equal to
minus the grand potential divided by volume, can be understood as the average
tangential mechanical pressure. The normal pressure is the normal mechanical
pressure, and the surface tension is the integral of the normal minus the
tangential mechanical pressure. The entropy and internal energy densities vary
with the slit pore height, confirming the observations of Israelachvili [1].
The entropy density increased monotonically with increasing height, while the
energy density oscillates. This confirms that the disjoining pressure is not
of entropic origin, it is of energetic origin [1].

By these investigations of the nanothermodynamic theory and the subsequent
simulations, we hope to have expanded on the knowledge on Hill’s method,
making it more available for further studies, for instance of flow and
reactions in porous media.

##  6 Acknowledgments

Computer resources have been provided by Faculty of Natural Science at NTNU
and by the HPC resources at UIT and NTNU provided by, www.sigma2.no. We are
grateful to the Research Council of Norway for funding through its Center of
Excellence funding scheme, project number 262644, PoreLab.

## References

  * [1] J. N. Israelachvili, Intermolecular and surface forces.  Academic press, 1985.
  * [2] T. M. McDonald, J. A. Mason, X. Kong, E. D. Bloch, D. Gygi, A. Dani, V. Crocella, F. Giordanino, S. O. Odoh, W. S. Drisdell, et al., “Cooperative insertion of $\text{CO}_{2}$ in diamine-appended metal-organic frameworks,” Nature, vol. 519, no. 7543, pp. 303–308, 2015.
  * [3] T. Vlugt, R. Krishna, and B. Smit, “Molecular simulations of adsorption isotherms for linear and branched alkanes and their mixtures in silicalite,” J. Phys. Chem. B, vol. 103, no. 7, pp. 1102–1118, 1999.
  * [4] F. Bresme and M. Oettel, “Nanoparticles at fluid interfaces,” J. Phys. Condens. Matter, vol. 19, no. 41, p. 413101, 2007.
  * [5] F. Bresme, H. Lehle, and M. Oettel, “Solvent-mediated interactions between nanoparticles at fluid interfaces,” J. Chem. Phys., vol. 130, no. 21, p. 214711, 2009.
  * [6] O. Galteland, F. Bresme, and B. Hafskjold, “Solvent-mediated forces between ellipsoidal nanoparticles adsorbed at liquid–vapor interfaces,” Langmuir, vol. 0, no. 0, p. null, 0.  PMID: 33236631.
  * [7] B. Derjaguin, “Untersuchungen über die reibung und adhäsion, iv,” Kolloid Z., vol. 69, no. 2, pp. 155–164, 1934.
  * [8] M. Moura, E. G. Flekkøy, K. J. Måløy, G. Schäfer, and R. Toussaint, “Connectivity enhancement due to film flow in porous media,” Phys. Rev. Fluids, vol. 4, no. 9, p. 094102, 2019.
  * [9] D. Das and S. Hassanizadeh, Upscaling multiphase flow in porous media.  Springer, 2005.
  * [10] H. H. Khanamiri, C. F. Berg, P. A. Slotte, S. Schlüter, and O. Torsæter, “Description of free energy for immiscible two-fluid flow in porous media by integral geometry and thermodynamics,” Water Resour. Res., vol. 54, no. 11, pp. 9045–9059, 2018.
  * [11] R. T. Armstrong, J. E. McClure, V. Robins, Z. Liu, C. H. Arns, S. Schlüter, and S. Berg, “Porous media characterization using minkowski functionals: theories, applications and future directions,” Transp. Porous Med., vol. 130, no. 1, pp. 305–335, 2019.
  * [12] P. A. Slotte, C. F. Berg, and H. H. Khanamiri, “Predicting resistivity and permeability of porous media using minkowski functionals,” Transp. Porous Med., vol. 131, no. 2, pp. 705–722, 2020.
  * [13] S. Kjelstrup, D. Bedeaux, A. Hansen, B. Hafskjold, and O. Galteland, “Non-isothermal transport of multi-phase fluids in porous media. the entropy production,” Front. Phys., vol. 6, p. 126, 2018.
  * [14] S. Kjelstrup, D. Bedeaux, A. Hansen, B. Hafskjold, and O. Galteland, “Non-isothermal transport of multi-phase fluids in porous media. constitutive equations,” Front. Phys., vol. 6, p. 150, 2019.
  * [15] P. B. Balbuena, D. Berry, and K. E. Gubbins, “Solvation pressures for simple fluids in micropores,” J. Phys. Chem., vol. 97, no. 4, pp. 937–943, 1993\.
  * [16] K. E. Gubbins, Y. Long, and M. Śliwinska-Bartkowiak, “Thermodynamics of confined nano-phases,” J. Chem. Thermodyn., vol. 74, pp. 169–183, 2014\.
  * [17] D. Bedeaux and S. Kjelstrup, “Hill’s nano-thermodynamics is equivalent with gibbs’ thermodynamics for surfaces of constant curvatures,” Chem. Phys. Lett., vol. 707, pp. 40–43, 2018.
  * [18] M. A. Gjennestad and Ø. Wilhelmsen, “Thermodynamic stability of volatile droplets and thin films governed by the disjoining pressure in open and closed containers,” Langmuir, 2020.
  * [19] M. A. Gjennestad and Ø. Wilhelmsen, “Thermodynamic stability of droplets, bubbles and thick films in open and closed pores,” Fluid Phase Equilibr., vol. 505, p. 112351, 2020.
  * [20] B. A. Strøm, J.-M. Simon, S. K. Schnell, S. Kjelstrup, J. He, and D. Bedeaux, “Size and shape effects on the thermodynamic properties of nanoscale volumes of water,” Phys. Chem. Chem. Phys., vol. 19, no. 13, pp. 9016–9027, 2017.
  * [21] O. Galteland, D. Bedeaux, B. Hafskjold, and S. Kjelstrup, “Pressures inside a nano-porous medium. the case of a single phase fluid,” Front. Phys., vol. 7, p. 60, 2019.
  * [22] M. Erdős, O. Galteland, D. Bedeaux, S. Kjelstrup, O. A. Moultos, and T. J. Vlugt, “Gibbs ensemble monte carlo simulation of fluids in confinement: Relation between the differential and integral pressures,” Nanomaterials, vol. 10, no. 2, p. 293, 2020.
  * [23] M. T. Rauter, O. Galteland, M. Erdős, O. A. Moultos, T. J. Vlugt, S. K. Schnell, D. Bedeaux, and S. Kjelstrup, “Two-phase equilibrium conditions in nanopores,” Nanomaterials, vol. 10, no. 4, p. 608, 2020.
  * [24] B. A. Strøm, J. He, D. Bedeaux, and S. Kjelstrup, “When thermodynamic properties of adsorbed films depend on size: Fundamental theory and case study,” Nanomaterials, vol. 10, no. 9, p. 1691, 2020.
  * [25] D. Bedeaux, S. Kjelstrup, and S. K. Schnell, Nanothermodynamics. General theory.  NTNU, Trondheim, Norway, 2020.
  * [26] T. L. Hill, Thermodynamics of small systems, part 1.  Benjamin, 1963.
  * [27] T. L. Hill, Thermodynamics of small systems, part 2.  Benjamin, 1964.
  * [28] T. L. Hill and R. V. Chamberlin, “Extension of the thermodynamics of small systems to open metastable states: An example,” Proc. Natl. Acad. Sci. U.S.A, vol. 95, no. 22, pp. 12779–12782, 1998.
  * [29] T. L. Hill and R. V. Chamberlin, “Fluctuations in energy in completely open small systems,” Nano Lett., vol. 2, no. 6, pp. 609–613, 2002.
  * [30] J.-P. Hansen and I. R. McDonald, Theory of simple liquids.  Elsevier, 1990.
  * [31] C. Radke, “Film and membrane-model thermodynamics of free thin liquid films,” J. Colloid Interf. Sci., vol. 449, pp. 462–479, 2015.
  * [32] Y. Long, J. C. Palmer, B. Coasne, M. Śliwinska-Bartkowiak, and K. E. Gubbins, “Pressure enhancement in carbon nanopores: a major confinement effect,” Phys. Chem. Chem. Phys., vol. 13, no. 38, pp. 17163–17170, 2011\.
  * [33] D. van Dijk, “Comment on ”pressure enhancement in carbon nanopores: a major confinement effect” by y. long, j. c. palmer, b. coasne, m. sliwinska-bartkowiak and k. e. gubbins, phys. chem. chem. phys., 2011, 13, 17163,” Phys. Chem. Chem. Phys., vol. 22, no. 17, pp. 9824–9825, 2020\.
  * [34] Y. Long, J. C. Palmer, B. Coasne, K. Shi, M. Śliwińska-Bartkowiak, and K. E. Gubbins, “Reply to the ‘comment on “pressure enhancement in carbon nanopores: a major confinement effect”’by d. van dijk, phys. chem. chem. phys., 2020, 22,” Phys. Chem. Chem. Phys., vol. 22, no. 17, pp. 9826–9830, 2020.
  * [35] P. Schofield and J. R. Henderson, “Statistical mechanics of inhomogeneous fluids,” Proc. R. Soc. Lon. Ser.-A, vol. 379, no. 1776, pp. 231–246, 1982\.
  * [36] A. Harasima, “Molecular theory of surface tension,” Adv. Chem. Phys, vol. 1, pp. 203–237, 1958.
  * [37] J. Irving and J. G. Kirkwood, “The statistical mechanical theory of transport processes. iv. the equations of hydrodynamics,” J. Chem. Phys., vol. 18, no. 6, pp. 817–829, 1950.
  * [38] B. Hafskjold and T. Ikeshoji, “Microscopic pressure tensor for hard-sphere fluids,” Phys. Rev. E, vol. 66, no. 1, p. 011203, 2002.
  * [39] K. Shi, Y. Shen, E. E. Santiso, and K. E. Gubbins, “Microscopic pressure tensor in cylindrical geometry: Pressure of water in a carbon nanotube,” J. Chem. Theory Comput., vol. 16, no. 9, pp. 5548–5561, 2020.
  * [40] T. Ikeshoji, B. Hafskjold, and H. Furuholt, “Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and spherical layers,” Mol. Simulat., vol. 29, no. 2, pp. 101–109, 2003.
  * [41] R. Evans, U. M. B. Marconi, and P. Tarazona, “Fluids in narrow pores: Adsorption, capillary condensation, and critical points,” J. Chem. Phys., vol. 84, no. 4, pp. 2376–2399, 1986.
  * [42] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, vol. 1.  Elsevier, 2001.
  * [43] W. Shinoda, M. Shiga, and M. Mikami, “Rapid estimation of elastic constants by molecular dynamics simulation under constant stress,” Phys. Rev. B, vol. 69, no. 13, p. 134103, 2004.
  * [44] S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys., vol. 117, no. 1, pp. 1–19, 1995.
  * [45] B. Hafskjold, K. P. Travis, A. B. Hass, M. Hammer, A. Aasen, and Ø. Wilhelmsen, “Thermodynamic properties of the 3d lennard-jones/spline model,” Mol. Phys., vol. 117, no. 23-24, pp. 3754–3769, 2019.
  * [46] A. Stukowski, “Visualization and analysis of atomistic simulation data with ovito–the open visualization tool,” Model. Simul. Mater. Sc., vol. 18, no. 1, p. 015012, 2009.


# On the set of Gâteaux differentiability of the $L^{1}$ norm

Julio Delgado  Universidad del Valle
Departamento de Matematicas
Calle 13 100-00
Cali
Colombia<EMAIL_ADDRESS>and  Andrés F. Muñoz-Tello
Universidad Santiago de Cali
Calle 5 62-00
Cali
Colombia<EMAIL_ADDRESS>
###### Abstract.

Let $(\Omega,\mathcal{M},\mu)$ be a measure space. In this paper we establish
the set of Gâteaux differentiability for the usual norm of $L^{1}(\Omega,\mu)$
and the corresponding derivative formulae at each point in this set.

###### Key words and phrases:

Gâteaux differentiability, Set of differentiability

###### 1991 Mathematics Subject Classification:

Primary 58C20; Secondary 58C25

The first author was supported by Grant CI71234 Vic. Inv. Universidad del
Valle.

##  1\. Introduction

The investigation on the sets of differentiability of Lipschitz functions goes
back to the classical Rademacher’s Theorem [1] of 1919. This theorem states
that:

Every Lipschitz function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is
differentiable almost everywhere.

Regarding the converse of Rademacher’s Theorem it has been recently
established by Preiss and Speight (cf. [2]) that: if $n>1$ then there exists a
Lebesgue null set in $\mathbb{R}^{n}$ containing a point of differentiability
of each Lipschitz function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n-1}$. A
more detailed discussion around the Rademacher’s Theorem and some converse
versions can be found in [3]. An extension of Rademacher’s Theorem to
functions defined on separable Banach spaces and valued into Banach spaces was
established by Phelps in [4]:

######  Theorem 1.1 (Phelps).

Let $E$ be a separable real Banach space and $F$ a real Banach space with the
Radon-Nikodým property (RNP). If $T:G\rightarrow F$ is a locally Lipschitz
mapping defined on a nonempty open subset $G$ of $E$. Then $T$ is Gâteaux
differentiable in the complement of a Gaussian null subset of $G$.

Other extensions of Rademacher Theorem to Lipschitz maps between infinite-
dimensional Banach spaces have been investigated by several authors e.g.
Aronszajn [5], Bongiorno [6], Mankiewicz [7], Phelps [4], Preiss and Zajíček
[8], among others). The research on the differentiability on Banach spaces has
been a very active field in the last decades. We just quote some few recent
examples: In Potapov and Sukachev’s work [9], they answer to a question in the
theory of Schatten- Von Neumann ideals of whether their norms have same
differentiability properties as the norms of their commutative counterparts.
On the other hand Lindenstrauss et al. [10] prove that a real valued Lipschitz
function on an Asplund space has points of Fréchet differentiability. In in
[11] Goldys and Gozzi apply the Phelps Theorem in the proof of existence and
uniqueness of solutions of second order parabolic Hamilton-Jacobi-Bellman
equations in Hilbert spaces, obtaining optimal feedback for first order
stochastic PDEs that arise in economic theory, the theory of population
dynamics and financial models.

The study of the set of differentiability of Lipschitz functions on Banach
spaces has been carried out for instance in Deville et al. [12], Diestel [13],
Dore and Maleva [14] among others. For the $L^{p}(\Omega)$ space where
$\Omega$ is a bounded domain and $1<p<\infty$, it is well known (cf. [13])
that the norm $\|f\|_{L^{p}}=(\int_{\Omega}|f(w)|^{p}dw)^{1/p}$ is Fréchet
differentiable and moreover its derivative is given by

$\partial_{F}\|f\|_{L^{p}}=\frac{|f|^{p-1}\operatorname*{sign\,}(f)}{\|f\|_{L^{p}}^{p-1}},$
(1.1)

for any $f\neq 0$. The case of Fréchet differentiability for $p=1$ is more
pathological. Here we will explicitly find the set of Gâteaux
differentiability for the $L^{1}(\Omega,\mu)$ norm, where
$(\Omega,\cal{M},\mu)$ is a measure space. The characterization of this set is
given in Theorem 3.3 of Section 3.

##  2\. Preliminaries

In this section we briefly give a basic background in our context. We start
with the two most classical notions of differentiability on Banach spaces. The
most simple one is that of Gâteaux differentiability. Let $X$, $Y$ be Banach
spaces and $\Omega\subseteq X$ an open set. A function $f:\Omega\rightarrow Y$
is Gâteaux differentiable ($G$-differentiable) at a point $a\in\Omega$, if

$\lim\limits_{t\rightarrow 0}\frac{f(a+th)-f(a)}{t}=u\left(h\right)$

exists for all $h\in X$, and $u$ defines a bounded linear operator from $X$
into $Y$. The operator $u$ is called the Gâteaux derivative of $f$ at the
point $a$ and is denoted by $\partial_{G}f(a)$.

A function $f:\Omega\rightarrow Y$ is Fréchet differentiable
($F$-differentiable) at $a$, if there is a bounded linear operator $u$ such
that

$f\left(a+h\right)-f\left(a\right)-u\left(h\right)=r\left(h\right)$

where $\lim\limits_{h\rightarrow
0}\frac{\left\|r\left(h\right)\right\|_{Y}}{\left\|h\right\|_{X}}=0$. The
operator $u$ is called the Fréchet derivative of $f$ at the point $a$ and is
denoted by $\partial_{F}f(a)$.

It is not difficult to see that, if $f$ is $F$-differentiable at $a$ then $f$
is $G$-differentiable at $a$ and in that case
$\partial_{F}f(a)=\partial_{G}f(a)$. For a given function
$\varphi:X\rightarrow Y$, we say that $B$ is the set of $G$-differentiability
or $F$-differentiability of $\varphi$, if respectively for every element of
the set $B$ the function $\varphi$ is $G$-differentiable or
$F$-differentiable.

On the other hand, a set is called a $G_{\delta}$ set if it can be expressed
as a countable intersection of open sets. This type of set is fundamental in
the definition of the F-Asplund and $G$-Asplund spaces, since these
respectively are Banach spaces in which all convex and continuous function
defined in an open and convex subset is $F$-differentiable (resp.
$G$-differentiable) in a dense and $G_{\delta}$ set. The following theorem
shows important examples of $G$-Asplund spaces.

######  Theorem 2.1 (Mazur).

Every separable Banach space is a $G$-Asplund space.

Moreover, it is well known that the set of all $F$-Asplund spaces is strictly
contained in the set of all $G$-Asplund spaces. In fact, as proved in Deville
et al. [12], the usual norm of $\ell^{1}(\mathbb{R})$ is not
$F$-differentiable at any point of the space $\ell^{1}(\mathbb{R})$. Then
$\ell^{1}(\mathbb{R})$ is not a $F$-Asplund space; but $\ell^{1}(\mathbb{R})$
is a separable Banach space and therefore a $G$-Asplund space by Mazur’s
theorem.

As we have quoted before, the $L^{p}$ norm is $F$-differentiable at any point
in the case $1<p<\infty$ and its corresponding derivative is given by (1.1).
Moreover, since the dual of $L^{p}(\Omega)$ is $L^{q}(\Omega)$, for
$1<p<\infty$ and $\frac{1}{p}+\frac{1}{q}=1$, and they are separable Banach
spaces, by Theorem 2.12 in Phelps [15] the space $L^{p}(\Omega)$ is a
$F$-Asplund space for $1<p<\infty$.

######  Remark 2.2.

In the special case of the $L^{1}(\Omega,\cal{M},\mu)$ norm we observe the
following regarding the set of $G$-differentiability. It is well known
that(cf. [16]) $L^{1}(\Omega,\cal{M},\mu)$ is separable provided $\mu$ is
$\sigma$-finite and the $\sigma$-algebra $\cal{M}$ is countably generated.
Since $\mathbb{R}$ has the RNP, then the $L^{1}(\Omega,\mu)$ norm satisfies
the conditions of Phelps’s Theorem. Therefore, the Gaussian measure of the
complement of the $G$-differentiability set is zero and by Mazur’s Theorem
this $G$-differentiability set is dense and a $G_{\delta}$ set.

In the next section we will give an explicit description of the set of
$G$-differentiability for the usual norm of $L^{1}(\Omega,\cal{M},\mu)$
without assuming the condition of separability.

##  3\. The $G$-differentiability set for the $L^{1}(\Omega,\mu)$ norm

Given a measure space $(\Omega,\cal{M},\mu)$, we denote by
$\mathcal{L}^{1}(\Omega,\mu)$ the space of the real valued integrable
functions over $\Omega$. On $\mathcal{L}(\Omega,\mu)$ one defines the
pseudonorm $\varphi(f)=\int_{\Omega}|f(x)|d\mu(x)$. One also defines on
$\mathcal{L}(\mu)$ an equivalence relation $\mathcal{R}$ by,
$f\mathcal{R}g\Longleftrightarrow f=g\,\mu-a.e.$ The equivalence class of $f$
is denoted by $\tilde{f}.$ Thus, the quotient
$\mathcal{L}(\Omega)/\mathcal{R}$ defines the classical Banach space
$L^{1}(\Omega,\mu)$ endowed with the norm
$\varphi(\tilde{f})=\int_{\Omega}|f(x)|d\mu(x)$.

We also recall recall that the signature $\operatorname*{sign\,}(x)$ of a real
number $x$, is defined by $\operatorname*{sign\,}(x)=\frac{x}{|x|}$ if $x\neq
0$, and $\operatorname*{sign\,}(x)=0$ if $x=0$.

In order to prepare the proof of the main theorem we first establish a mild
lemma which will help to obtain a pointwise formulae for the Gâteaux
derivative.

######  Lemma 3.1.

Let $f,\,h:\Omega\rightarrow\mathbb{R}$ be functions and $x\in\Omega$ such
that $h(x)\neq 0$. Then

$\lim\limits_{t\rightarrow 0}\frac{|f(x)+th(x)|-|f(x)|}{t},$ (3.1)

exists if and only if $f(x)\neq 0$.
Moreover, if that is the case

$\lim\limits_{t\rightarrow
0}\frac{|f(x)+th(x)|-|f(x)|}{t}=\operatorname*{sign\,}(f(x))h(x).$ (3.2)

###### Proof.

If the limit (3.1) exists and $f(x)=0$, then the limit

$\lim\limits_{t\rightarrow 0}\frac{|th(x)|}{t}=\lim\limits_{t\rightarrow
0}\frac{|t||h(x)|}{t}$

does not exists, which is a contradiction, hence $f(x)\neq 0$.

We now assume $f(x)\neq 0$. Since $h(x)\neq 0$, we can choose
$\delta=\frac{|f(x)|}{2|h(x)|}$ and observe that for $|t|<\delta$ we have

$\operatorname*{sign\,}(f(x)+th(x))=\operatorname*{sign\,}(f(x)).$ (3.3)

Hence

$\displaystyle\frac{|f(x)+th(x)|-|f(x)|}{t}$ $\displaystyle=$
$\displaystyle\frac{\operatorname*{sign\,}(f(x)+th(x))(f(x)+th(x))-\operatorname*{sign\,}(f(x))f(x)}{t}$
$\displaystyle=$ $\displaystyle\frac{\operatorname*{sign\,}(f(x))th(x)}{t}$
$\displaystyle=$ $\displaystyle\operatorname*{sign\,}(f(x))h(x).$

Therefore

$\lim\limits_{t\rightarrow
0}\frac{|f(x)+th(x)|-|f(x)|}{t}=\operatorname*{sign\,}(f(x))h(x).$

∎

We now consider a measure space $(\Omega,\cal{M},\mu)$. As a consequence of
the lemma above we obtain the following corollary.

######  Corollary 3.2.

Let $f,h:\Omega\rightarrow\mathbb{R}$ be measurable functions such that
$h(x)\neq 0$ $\mu$-a.e. Then

$\lim\limits_{t\rightarrow 0}\frac{|f(x)+th(x)|-|f(x)|}{t}\mbox{ exists
}\mu-a.e$ (3.4)

if and only if $f(x)\neq 0$ $\mu$-a.e. Moreover, if that is the case the
formula (3.2) holds.

###### Proof.

We suppose that $h(x)\neq 0$ for all $x\in A$, with $\mu(A^{c})=0$ and first
assume that the limit (3.4) exists $\mu$-a.e. Hence, $f(x)\neq 0$ for every
$x\in A$ by Lemma 3.1, thus $f(x)\neq 0$ $\mu$-a.e.
On the other hand, if $f(x)\neq 0$ $\mu$-a.e. Then there exists a set
$D\in\cal{M}$ such that $\mu(D^{c})=0$ and $f(x)\neq 0$ for all $x\in D$. By
taking $x\in A\cap D$ we have that $f(x),h(x)\neq 0$. Then by Lemma 3.1 the
limit (3.4) exists for every $x\in A\cap D$. Therefore the limit (3.4) exists
$\mu$-a.e.
The last part of the statement clearly follows from the last part of Lemma 3.1
by applying the formula (3.2) in a set $\Gamma$ where we have $f(x),h(x)\neq
0$ and $\mu(\Gamma^{c})=0.$ ∎

We are now ready to present our main result. We consider a measure space
$(\Omega,\cal{M},\mu)$ and the following condition:

(A1) For every $I\in\cal{M}$ with $\mu(I)=\infty,$ there exists a set
$F\subset I$ with $F\in\cal{M}$ such that

$0<\mu(F)<\infty.$

Every $\sigma$-finite measure that is not the trivial null measure satisfies
the condition (A1). The counting measure on the power set of $\mathbb{R}$
satisfies the condition (A1) and is not $\sigma$-finite. On the other hand,
every non-atomic measure also satisfies the condition (A1).

The following theorem gives a precise description of the $G$-differentiability
set for the usual $L^{1}(\Omega,\mu)$ norm and the corresponding formulae for
the Gâteaux derivative.

Let $G$ be the subset of $L^{1}(\Omega,\mu)$ consisting of equivalence classes
$\tilde{f}$ such that $f\neq 0$ $\mu$-a.e., for some representative $f$ of the
class $\tilde{f}$. It is clear that if $\tilde{f}\in G$, any representative of
the class $\tilde{f}$ satisfies the condition defining $G$.

######  Theorem 3.3.

Let $(\Omega,\cal{M},\mu)$ be a measure space satisfying the condition (A1).
Let $B$ be the subset of $L^{1}(\Omega,\mu)$ above defined. Then the set of
$G$-differentiability of the norm $\varphi$ is $G$. Moreover, the Gâteaux
derivative of $\varphi$ for $\tilde{f}\in G$ and the direction $h\in
L^{1}(\Omega,\mu)$ is given by

$\partial_{G}\varphi(\tilde{f})(h)=\int_{\Omega}\operatorname*{sign\,}(f(x))h(x)d\mu(x),$
(3.5)

where $f$ is any representative of the class $\tilde{f}$.

###### Proof.

We assume $\tilde{f}\in G$ and that $f$ is a representative of the class
$\tilde{f}$. We start by observing that for any $h\in L^{1}(\Omega,\mu)$ the
right hand side of the formula (3.5) is independent of the representative in
the class $\tilde{f}$.

We are now going to prove the $G$-differentiability of $\varphi$ at
$\tilde{f}$. Let $h\in L^{1}(\Omega,\mu)\setminus\\{0\\}$, an application of
the Lebesgue dominated convergence Theorem gives us

$\displaystyle\lim\limits_{t\rightarrow 0}\frac{\varphi(f+th)-\varphi(f)}{t}$
$\displaystyle=$ $\displaystyle\lim\limits_{t\rightarrow
0}\frac{\int\limits_{\Omega}|f(x)+th(x)|d\mu(x)-\int\limits_{\Omega}|f(x)|d\mu(x)}{t}$
$\displaystyle=$ $\displaystyle\lim\limits_{t\rightarrow
0}\,\frac{\int\limits_{\Omega}\left(|f(x)+th(x)|-|f(x)|\right)d\mu(x)}{t}$
$\displaystyle=$ $\displaystyle\int\limits_{\Omega}\lim\limits_{t\rightarrow
0}\frac{\,\left(|f(x)+th(x)|-|f(x)|\right)}{t}d\mu(x).$

Since $h(x)\neq 0$ $\mu$-a.e. and $f(x)\neq 0$ $\mu$-a.e., by applying
Corollary 3.2 we have

$\lim\limits_{t\rightarrow
0}\frac{|f(x)+th(x)|-|f(x)|}{t}=\operatorname*{sign\,}(f(x))h(x)$

$\mu$-a.e., and

$\lim\limits_{t\rightarrow
0}\frac{\varphi(f+th)-\varphi(f)}{t}=\int_{\Omega}\operatorname*{sign\,}(f(x))h(x)d\mu(x).$

We now assume that $\tilde{f}\in G^{c}$ and we will show that $\varphi$ does
not have $G$-derivative in $\tilde{f}$. For the set
$Z_{1}=\\{x\in\Omega:f(x)=0\\}$ we have $\mu(Z_{1})>0.$ If $\mu(Z_{1})=\infty$
by condition (A1) we can choose a subset $Z$ of $Z_{1}$ such that
$0<\mu(Z)<\infty.$ We consider $h=1_{Z}$, thus $h\in
L^{1}(\Omega,\mu)\setminus 0$. Since $h=0$ on $Z^{c}$, for $t\in\mathbb{R}$ we
have $|f+th|-|f|=|f|-|f|=0$ on $Z^{c}$. Hence, for the calculation of the
$G$-derivative of $\varphi$ at $f$ in the direction $h$, is enough to consider
the integration on $Z$. For $t\neq 0$ we have

$\displaystyle\frac{\varphi(f+th)-\varphi(f)}{t}$ $\displaystyle=$
$\displaystyle\frac{\left(\int_{Z}|f(x)+th(x)|d\mu(x)-\int_{Z}|f(x)|d\mu(x)\right)}{t}$
$\displaystyle=$ $\displaystyle\frac{\int_{Z}|t|d\mu(x)}{t}$ $\displaystyle=$
$\displaystyle\frac{|t|}{t}\mu(Z)$ $\displaystyle=$
$\displaystyle\operatorname*{sign\,}(t)\mu(Z).$

Hence, we have $(\partial_{G}^{+}\varphi)f(h)=\mu(Z)$ and
$(\partial_{G}^{+}\varphi)f(h)=-\mu(Z)$. Therefore the $G$-derivative of
$\varphi$ at $f$ does not exists.

It is clear that the formula (3.5) defines a bounded linear operator from
$L^{1}(\Omega,\mu)$ into $\mathbb{R}$ for every $\tilde{f}$ in $G$, and this
concludes the proof of the theorem. ∎

######  Remark 3.4.

(1) In the case $\Omega=\mathbb{N}$ and the counting measure, Theorem 3.3
absorbs a more classical result on the Gâteaux differentiability of the
$\ell^{1}(\mathbb{N})$-norm, namely that the set of Gâteaux differentiability
of the $\ell^{1}(\mathbb{N})$-norm is the set $G\subset\ell^{1}(\mathbb{N})$
of sequences $(x_{n})_{n}$ such that $x_{n}\neq 0$ for all $n$, and the
Gâteaux derivative of $\varphi$ at $x$ and in the direction $h$ is given by

$\partial_{G}\varphi(x)(h)=\sum\limits_{n}\operatorname*{sign\,}(x_{n})h_{n}.$

(2) The condition (A1) allows $L^{1}$ to be be non-separable. Therefore,
Theorem 3.3 gives a description of the set of differentiability beyond the
scope of Phelps’s and Mazur’s theorem.

In the special separable case, as a consequence of Theorem 3.3 and Phelps
Theorem 1.1 we have:

######  Corollary 3.5.

Let $(\Omega,\cal{M},\mu)$ be a measure space satisfying the condition (A1)
and such that $L^{1}(\Omega,\cal{M},\mu)$ is separable. Then, the set of
Gâteaux differentiability of the $L^{1}$-norm is $G$, $G$ is a $G_{\delta}$
set and $G^{c}$ is a Gaussian null set.

As quoted before $L^{1}(\Omega,\cal{M},\mu)$ is separable provided $\mu$ is
$\sigma$-finite and the $\sigma$-algebra $\cal{M}$ is countably generated. If
the measure $\mu$ is not the trivial null measure then $\mu$ satisfies the
condition (A1).

## References

  * [1] H. Rademacher, Über partielle und totale differenzierbarkeit von Funktionen mehrerer Variabeln und über die Transformation der Doppelintegrale, Mathematische Annalen 79 (4) (1919) 340–359.
  * [2] D. Preiss, and G. Speight, Differentiability of Lipschitz functions in Lebesgue null sets, Inventiones Mathematicae 199 (2) (2014) 517–559.
  * [3] M. C. G. Alberti, D.Preiss, Differentiability of Lipschitz functions, structure of null sets, and other problems., Proceedings of the International Congress of Mathematicians, New Delhi. 2010 III 1379–1394.
  * [4] R. Phelps, Gaussian null sets and differentiability of Lipschitz map on Banach spaces., Pacific J. Math. 77 (2) (1978) 523–531.
  * [5] N. Aronszajn, Differentiability of lipschitzian mappings between Banach spaces, Studia Mathematica 57 (1976) 147–190.
  * [6] D. Bongiorno, Metric differentiability of Lipschitz maps, Journal of the Australian Mathematical Society 96 (1) (2014) 25–35.
  * [7] P. Mankiewicz, On the differentiability of Lipschitz mappings in Fréchet spaces, Studia Mathematica 45 (1) (1973) 15–29.
  * [8] D. Preiss, and L. Zajíček, Directional derivatives of Lipschitz functions, Israel Journal of Mathematics 125 (2001) 1–27.
  * [9] D. Potapov, and F. Sukochev, Fréchet differentiability of Sp norms, Advances in Mathematics 262 (2014) 436–475.
  * [10] J. Lindenstrauss, D. Preiss, and J. Tišer, Fréchet differentiability of Lipschitz functions via a variational principle, Journal of the European Mathematical Society 12 (2) (2010) 385–412.
  * [11] B. Goldys, and F. Gozzi, Second order parabolic Hamilton-Jacobi-Bellman equations in Hilbert spaces and stochastic control: $L^{2}\mu$ approach, Stochastic Processes and their Applications 116 (12) (2006) 1932–1963.
  * [12] R. Deville, G. Godefroy, and V. Zizler., Smoothness and Renormings in Banach Spaces, 1st Edition, Longman Scientific&Technical, Essex, England, 1993.
  * [13] J. Diestel, Geometry of Banach spaces : selected topics, Springer-Verlag, 1975\.
  * [14] M. Doré, and O. Maleva, A universal differentiability set in Banach spaces with separable dual, Journal of Functional Analysis 261 (6) (2011) 1674–1710.
  * [15] R. Phelps, Convex Functions, Monotone Operators and Differentiability, Vol. 1364 of Lecture Notes in Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, 1989.
  * [16] D. Cohn., Measure Theory, Birkhauser, Basel, 2013.