# A statistical perspective on microsolvation

Mohammad Rahbar  Christopher J. Stein<EMAIL_ADDRESS>Theoretische
Physik and CENIDE, Universität Duisburg-Essen, D-47048 Duisburg, Germany

###### Abstract

The lack of a procedure to determine equilibrium thermodynamic properties of a
small system interacting with a bath is frequently seen as a weakness of
conventional statistical mechanics. A typical example for such a small system
is a solute surrounded by an explicit solvation shell. One way to approach
this problem is to enclose the small system of interest in a large bath of
explicit solvent molecules, considerably larger than the system itself. The
explicit inclusion of the solvent degrees of freedom is obviously limited by
the available computational resources. A potential remedy to this problem is a
microsolvation approach where only a few explicit solvent molecules are
considered and surrounded by an implicit solvent bath. Still, the sampling of
the solvent degrees of freedom is challenging with conventional grand
canonical Monte Carlo methods, since no single chemical potential for the
solvent molecules can be defined in the realm of small-system thermodynamics.

In this work, a statistical thermodynamic model based on the grand canonical
ensemble is proposed that avoids the conventional system size limitations and
accurately characterizes the properties of the system of interest subject to
the thermodynamic constraints of the bath. We extend an existing
microsolvation approach to a generalized multi-bath "micro-statistical" model
and show that the previously derived approaches result as a limit of our
model. The framework described here is universal and we validate our method
numerically for a Lennard–Jones model fluid.

###### keywords:

American Chemical Society, LaTeX

Current affiliation: Technical University of Munich; TUM School of Natural
Sciences, Department of Chemistry, Lichtenbergstr. 4, D-85748 Garching,
Germany Current affiliation: Technical University of Munich; TUM School of
Natural Sciences, Department of Chemistry, Lichtenbergstr. 4, D-85748
Garching, Germany IR,NMR,UV

##  1 Introduction

When we take the limit of $N\longrightarrow\infty$ particles, thermodynamics
emerges from statistical mechanics. We can provide a microscopic justification
for thermodynamics when we evaluate the relative fluctuations in thermodynamic
quantities such as energy and show that they vanish as $1/\sqrt{N}$.1 When
fluctuations are ignored, thermodynamics is the statistical mechanics of near-
equilibrium systems. However, this interpretation is no longer valid when a
small system is examined. In other words, the lack of a procedure to determine
equilibrium thermodynamic variables of a small system interacting with a bath
is frequently seen as a weakness of conventional statistical mechanics.2, 3,
4, 5 One way to deal with this problem is to directly simulate a small system
enclosed in a bath in the thermodynamic limit, which typically requires the
explicit inclusion of a large portion of the bath compared to the size of the
small system. Such an approach is often constrained in practice by the
available computational resources 6. This is exacerbated by the generally high
computational cost of accurate ab initio calculations. For purely short-ranged
interactions an often used conventional technique is to adopt periodic
boundary conditions and empirically assess the system size dependence of the
computed thermodynamical properties by carrying out simulations with
successively larger baths. Typical problems in computational chemistry of
small systems interacting with a large environment, such as molecule solvation
or the study of nanoconfined species, however, can be substantially affected
by long-range electrostatic interactions7, 8, 9. As a result, implicit solvent
models that ignore the explicit molecular structure of the solvent bath in
favor of an approximate description of its effect became increasingly popular
over the last three decades.10, 11, 12, 13, 14, 15, 16 Among the drawbacks of
such implicit models is the fact that the accuracy of implicit solvent models
relies strongly on the geometric details of the solute-solvent interface,
which is often described by simple models rather than physical properties of
the solute-solvent interactions 17, 18, 19. Unlike implicit models, explicit
models offer a physical and spatially precise description of the solvent.
However, as stated above, due to the sampling problem, many of the explicit
models are computationally expensive and therefore may fail to reproduce
experimental results. A potential remedy are hybrid methods that combine
implicit and explicit schemes to reduce the computational cost while keeping
the solvent’s spatial resolution in the solvation sell. Despite the hybrid
approaches’ apparent effectiveness, there are unanswered questions concerning
the construction of microsolvation models20, 21, 22, 23 such as the subtleties
of the coupling to the bulk solvent, and the extension of these models to
solvent mixtures.

Understanding solvation phenomena hence requires quantitative knowledge of
occupancy statistics of solvent molecules around the solutes. Recently, Dixit
et al. 24 have shown that a single-bath grand canonical simulation does not
allow for an accurate description when the size of the solvation shell is
comparable to the size of the particle for a simple hard-sphere model system.
In this article, we propose and derive a statistical thermodynamic model based
on the grand canonical ensemble that avoids the conventional system size
limitations of typical simulations with periodic boundary conditions by
eliminating the need to explicitly include large parts of the thermodynamical
bath in the simulations. We extend an existing microsolvation approach to a
generalized multi-bath "micro-statistical" model to obtain the properties of
small systems and show that the previously derived approaches result as a
limit of our model. We then validate our method numerically for a
Lennard–Jones model fluid. From a statistical point of view, we propose a
model that allows us to sample the properties of a small core system
interacting with its bulk environment, the bath, when the core size is on the
order of the average inter-particle distance. This is achieved from a
simulation of the core alone.

##  2 Theory

Let us consider a system like the one shown on the left panel of Figure 1.

Figure 1: A) The composite system core+shell with the volume
$V=V_{\text{c}}+V_{\text{s}}$, in which the volumes of the core and shell
subsystems are $V_{\text{c}}$ and $V_{\text{s}}$, respectively, is in contact
with a bath, at the chemical potential $\mu$ and the reciprocal thermodynamic
temperature $\beta=1/k_{\text{B}}T$, where $k_{\text{B}}$ is the Boltzmann
constant. $V_{\text{s}}$ is very large compared to $V_{\text{c}}$, and
throughout this article, we consider $V\simeq V_{\text{s}}$ and
$V>>V_{\text{c}}$. B) A schematic representation of mapping the problem of
finding the microstate probability distribution for a small system to a system
coupled to theoretically (practically) infinitely (finitely) many artificial
thermodynamic baths with the corresponding artificial thermodynamic variables,
$\mu_{\text{c}}^{j}$ and $\beta_{\text{c}}^{j}$ (see Eq. (10)).

In the grand canonical ensemble, the composite system core+shell is in contact
with a bath and therefore at constant chemical potential $\mu$ and constant
temperature $T$. The microstates of the combined system are characterized by
the set $\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}$, in which $\vec{\mathbf{q}}$
and $\vec{\mathbf{p}}$ are position and momenta of all particles,
respectively. The probability distribution of the various microstates of this
composite system core+shell can be written as 25

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)\quad,$
(1)

with the volumes of the core and shell subsystem, $V_{\text{c}}$ and
$V_{\text{s}}$, respectively and the reciprocal thermodynamic temperature
$\beta=1/k_{\text{B}}T$ where $k_{\text{B}}$ is the Boltzmann constant. The
quantity
$p(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V)$
can be interpreted as the joint probability distribution of microscopic
degrees of freedom constrained by the macroscopic thermodynamic variables. So
we can write,

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)$
$\displaystyle=\frac{\exp\left(\beta\mu
N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}})-\beta\mathcal{H}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}{\mathcal{Q}\left(\mu,\beta,V\right)}$
(2)

with the grand partition function $\mathcal{Q}$, the total number of particles
$N$ and the Hamiltonian of combined system $\mathcal{H}$. The grand partition
function is defined as

$\displaystyle{\mathcal{Q}\left(\mu,\beta,V\right)}=$
$\displaystyle\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}}{\exp\left(\beta\mu
N\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)-\beta\mathcal{H}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}.$
(3)

Here and in the remainder of the article we write an integration over
microstates as an infinite sum. The Hamiltonian of the combined core+shell
system can be partitioned as follows

$\displaystyle\mathcal{H}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)$
$\displaystyle\equiv\mathcal{H}_{\mathrm{s}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}}\right)+\mathcal{H}_{\mathrm{c}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)+\mathcal{H}_{{\mathrm{int}}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right),$
(4)

where $\mathcal{H}_{\mathrm{c}(\mathrm{s})}$ is the Hamiltonian of the
core(shell) system, respectively and $\mathcal{H}_{\mathrm{int}}$ is the
interaction between the subsystems. Our actual quantity of interest, the
marginal distribution,25
$p(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V)$,
can be calculated by integrating the joint probability distribution,
$p(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V)$,
over all shell degrees of freedom
$\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}}$

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)=\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}}}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid,\mu,\beta,V\right)$
$\displaystyle=\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}}}\frac{\exp\left(\beta\mu
N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}})-\beta\mathcal{H}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}{\mathcal{Q}\left(\mu,\beta,V\right)}$
$\displaystyle={\exp\left(\beta\mu
N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}})-\beta\mu\mathcal{H}_{\mathrm{c}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}$
$\displaystyle\times\underbrace{\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}}}\frac{\exp\left(\beta\mu
N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}})-\beta\mathcal{H}_{\mathrm{s}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}}\right)-\beta\mathcal{H}_{{\mathrm{int}}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{s}}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}{{\mathcal{Q}\left(\mu,\beta,V\right)}}}_{\mathcal{W}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)}$
$\displaystyle={\exp\left(\beta\mu
N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}})-\beta\mathcal{H}_{\mathrm{c}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)\mathcal{W}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)}\,.$
(5)

In this expression for the marginal distribution, the interaction between core
and shell is exclusively described by
$\mathcal{W}(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V)$.
This interaction term is highly dependent on the characteristics of the
core+shell interaction Hamiltonian, $\mathcal{H}_{\text{int}}$, which in
general is difficult to approximate. As a result, while Eq. (2) is strictly
correct, it has no practical significance. An efficient method to calculate
the distribution of the microstates of the core system conditional to the
thermodynamic constraints of the bath — which is not in direct contact with
the core part of the system — without the need of an explicit calculation of
the interaction Hamiltonian $\mathcal{H}_{\text{int}}$ is therefore highly
desired.

The need to explicitly calculate the interaction is rooted in a microstate
interpretation of the joint probability distribution. One way to circumvent
this issue is to introduce $j$ macrostate degrees of freedom and corresponding
artificial thermodynamic variables, $\mu_{\text{c}}^{j}$ and
$\beta_{\text{c}}^{j}$, and link these infinite macrostate degrees of freedom
to the core system microstates. This link is provided by the law of total
probability, which can also be stated for conditional probabilities

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)=\sum_{j}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V,\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\mid\mu,\beta,V\right)\,.$
(6)

Now, using the chain rule for conditional probability (see the Appendix), we
have

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V,\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=\frac{p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}},\mu,\beta,V,\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}{p\left(\mu,\beta,V,\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}$
$\displaystyle=\frac{p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu,\beta,V\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)}{p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\mid\mu,\beta,V\right)p\left(\mu,\beta,V\right)}\,.$
(7)

For a sufficiently large shell compared to a small core, we can conclude that
the newly introduced artificial thermodynamic constraints on the core do not
affect the thermodynamic properties of the shell (which is still in
equilibrium with the bath). 24 As a result, the artificial thermodynamic
variables of the core system and the shell are almost independent. With this
assumption we arrive at

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\mid\mu,\beta,V\right)\approx
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$ (8)

and

$\displaystyle
p\left(\mu,\beta,V\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}},\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\approx
p\left(\mu,\beta,V\right)\,.$ (9)

Inserting the results of Eqns. (2)-(9) into Eq. (6), we find

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)$
$\displaystyle\approx\sum_{j}\frac{p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu,\beta,V\right)}{p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu,\beta,V\right)}p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle\approx\sum_{j}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\,.$
(10)

This constitutes an approximation of the core microstate distribution under
the thermodynamical constraints of a bath which is not in direct contact with
this core. The first term in the sum,
$p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$,
is the probability for a given core microstate under the constraint of the
artificial thermodynamic variables assigned to the core, which is the standard
grand canonical probability distribution for a chosen set of
$\mu^{j}_{\text{c}},\beta^{j}_{\text{c}}$ and $V_{\text{c}}$

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=\frac{\exp\left(\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}})-\beta^{j}_{\mathrm{c}}\mathcal{H}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}{\mathcal{Q}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\text{c}}\right)}\,.$
(11)

We can interpret Eq. (10) as a weighted average of the probability
distributions of infinitely many artificial thermodynamic baths with the
individual weights given by
$p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$,
Figure 1 (right panel). In the limit of the core volume approaching the shell
volume, $\lim{V_{\text{c}}}\rightarrow V_{\text{s}}$, we expect that we can
describe the probability distribution by a single bath with conventional grand
canonical distribution function. This seemingly obvious statement reveals an
important property of the probability distribution
$p(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})$. Taking this
limit on both sides of Eq. (10) we find

$\displaystyle\underbrace{\lim_{V_{\mathrm{c}}\rightarrow
V}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)}_{p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V}\mid\mu,\beta,V\right)}$
$\displaystyle=\sum_{j}\lim_{V_{\mathrm{c}}\rightarrow
V}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\lim_{V_{{\mathrm{c}}}\rightarrow
V}p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
(12)
$\displaystyle=\iint_{\mu_{\mathrm{c}}\beta_{\mathrm{c}}}\text{d}\mu_{\mathrm{c}}\text{d}\beta_{\mathrm{c}}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V}\mid\mu_{\mathrm{c}},\beta_{\mathrm{c}},V\right)p\left(\mu_{\mathrm{c}},\beta_{\mathrm{c}},V\right)\,.$
(13)

We see that $p(\mu_{\mathrm{c}},\beta_{\mathrm{c}},V)$ filters the value of
the integrand
$p(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V}\mid\mu_{\mathrm{c}},\beta_{\mathrm{c}},V)$,
at the point of its occurrence $(\mu,\beta)$ — a property fulfilled by the
Dirac delta function

$\displaystyle
p\left(\mu_{\mathrm{c}},\beta_{\mathrm{c}},V\right)\approx\delta^{2}(\mu_{\mathrm{c}}-\mu,\beta_{\mathrm{c}}-\beta)\,.$
(14)

At this point we mapped the problem of finding the microstate probability
distribution of a core system coupled indirectly to a thermodynamic bath via a
shell that is in equilibrium with this bath, to a core system coupled to
infinitely many artificial thermodynamic baths (see Figure 1). We are left
with the difficulty of calculating the coupling strengths (or weights),
$p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$, of
the individual baths. However, in light of the thermodynamic limit, we require
the weight distribution to collapse to a Dirac delta function for large core
volumes. We arrived at this interpretation without any reference to
characteristic properties of the system. The above-mentioned behaviour at the
thermodynamic limit was predicted in previous works 26, 24, 27 but has not
been rigorously proven so far.

###  2.1 Identifying a functional form for the weight distribution

The function $p(\mu^{j}_{\text{c}},\beta^{j}_{\text{c}},V_{\text{c}})$ assigns
a weight to a bath with thermodynamic variables $\mu^{j}_{\text{c}}$ and
$\beta^{j}_{\text{c}}$ coupled to the core of size $V_{\text{c}}$.
Unfortunately, the functional form of the weight function is unknown a priori.
We now attempt to rewrite this unknown function in a more familiar form.
Therefore, we start by applying the fundamental law of probability, meaning
that for a given set of artificial thermodynamic variables
$(\mu^{j}_{\text{c}},\beta^{j}_{\text{c}},V_{\text{c}})$, the sum of all the
probabilities
$p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
on all possible degrees of freedom
${\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\text{c}}}}$ must be equal to
one,

$\displaystyle\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=1\,.$
(15)

Multiplying by $p(\mu^{j}_{\text{c}},\beta^{j}_{\text{c}},V_{\text{c}})$ on
both sides of Eq. (15),

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
(16)

and using Eq. (11), we can write

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle={\mathcal{Q}^{-1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\underbrace{p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\sum_{\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}}{\exp\left(\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}})-\beta^{j}_{\mathrm{c}}\mathcal{H}_{\mathrm{c}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)\right)}}_{X\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}$
(17)

By now, we can interpret
$p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$ as
a weighted distribution with a known part (the inverse of the grand canonical
partition function) and a yet to be determined part
$X\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=\mathcal{Q}^{-1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)X\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\,.$
(18)

For notational simplicity, we replace
$(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})$ with $y$

$\displaystyle
p\left(y\right)=\mathcal{Q}^{-1}\left(y\right)X\left(y\right)\,.$ (19)

We know that $\mathcal{Q}^{-1}(y)$ has an exponential form with known argument
$g(y)$. We can write the argument as the sum of two parts $g_{1}(y)$ and
$g_{2}(y)$, respectively. It is trivial that after choosing one part (for
which we have infinitely many choices), the other part is uniquely determined.
Now we can write

$\displaystyle
p\left(y\right)=\exp\left({g_{1}\left(y\right)}+g_{2}\left(y\right)\right)X\left(y\right)=\exp\left({g_{1}\left(y\right)}+g_{2}\left(y\right)+\ln{X\left(y\right)}\right)\,.$
(20)

The idea is that among all of the infinite number of options for picking
$g_{2}(y)$, we choose $X(y)$ for our ansatz, i.e. $X(y)=g_{2}(y)$, leading to

$\displaystyle
p\left(y\right)=\exp\left({\overbrace{g_{1}\left(y\right)}^{\text{determined}}}+\overbrace{g_{2}\left(y\right)}^{\text{choose}}+\ln{\overbrace{g_{2}\left(y\right)}^{\text{ansatz}}}\right)\,.$
(21)

As is evident from Eq. (21), $p(y)$ is fully defined. With the original
variables we find

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle=\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\right)g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle={\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)+\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\right)}\,.$
(22)

Now the proposed functional form is proportional to the inverse of the
partition function, which has a physical interpretation. In summary, the
functional form of the grand partition function has been used as a guide for
constructing the weighted distribution. The remaining question is which part
of the inverse partition function should be chosen as $g_{2}(y)$ and hence
define $X(y)$. To identify a potential expression, we first need to inspect
the grand partition function more closely. This is the topic of the next
section.

###  2.2 Grand canonical partition function for an interacting system of
particles

The grand canonical partition function
${\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}$
can be written as a sum over canonical partition functions
$\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})$

$\displaystyle\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})=\sum_{N^{{}^{\prime}}=0}^{\infty}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\sum_{\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid
N^{{}^{\prime}},V_{\mathrm{c}}\right)}\mathrm{e}^{-\beta^{j}_{\mathrm{c}}\mathcal{H}_{N^{{}^{\prime}}}\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\right)}=\sum_{N^{{}^{\prime}}=0}^{\infty}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\,.$
(23)

It is obvious that if we evaluate the canonical partition function
$\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})$, then
finding the grand canonical partition functions is straightforward. As
discussed previously, we attempt to understand the role of more complicated
interactions between particles such as the Lennard–Jones model and investigate
how to treat them in our proposed framework. We consider a one-component
system of $N^{{}^{\prime}}$ spherical particles with mass $m,$ positions
$\vec{\mathbf{q}}$ and momenta $\vec{\mathbf{p}}$. Particles $i$ and $j$
interact via a pair-wise potential
$\phi_{ij}\left(\mathbf{r}_{i},\mathbf{r}_{j}\right)$. The Hamiltonian of this
system of interacting particles can be written as,

$\displaystyle\mathcal{H}_{N^{{}^{\prime}}}=\sum_{i=1}^{N^{{}^{\prime}}}\frac{\vec{{p}_{i}}^{2}}{2m}+\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\,,$
(24)

in which $\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N}\right)$ is total internal
potential function which can be written as

$\displaystyle\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)=\frac{1}{2}\sum_{i\neq
j}^{N^{{}^{\prime}}}\phi_{ij}\left({\vec{q}}_{i},{\vec{q}}_{j}\right)\,.$ (25)

We can now write the canonical partition function of the system of interacting
particles 28. (see the Appendix for details).

$\displaystyle\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})={\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}{\left\langle\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]\right\rangle_{\textbf{Non-
Int}}}$ (26)

Taking the natural logarithm on both sides of Eq. (26),

$\displaystyle\ln\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\ln{\left\langle\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]\right\rangle_{\textbf{Non-
Int}}}$ (27)

and using Eq. (108), we arrive at

$\displaystyle\ln\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\ln{\left\langle\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]\right\rangle_{\textbf{Non-
Int}}}$
$\displaystyle=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{{\mathrm{c}}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\ln{\left\langle\sum_{m=0}^{\infty}\frac{(-\beta^{j}_{\mathrm{c}})^{m}}{m!}\phi^{m}\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right\rangle_{\textbf{Non-
Int}}}$
$\displaystyle=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\ln\left[{\sum_{m=0}^{\infty}\frac{(-\beta^{j}_{\mathrm{c}})^{m}}{m!}\left\langle\phi^{m}\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right\rangle_{\textbf{Non-
Int}}}\right]$
$\displaystyle=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\ln\left[\exp\left[{\sum_{m=1}^{\infty}\frac{(-\beta^{j}_{\mathrm{c}})^{m}}{m!}\left\langle\phi^{m}\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right\rangle_{\textbf{Non-
Int, c}}}\right]\right]$
$\displaystyle=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+{\sum_{m=1}^{\infty}\frac{(-\beta^{j}_{\mathrm{c}})^{m}}{m!}\left\langle\phi^{m}\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right\rangle_{\textbf{Non-
Int,c}}}\,.$ (28)

Considering explicitly only the first and second term of the expansion we
obtain

$\displaystyle\ln\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+{\frac{(-\beta^{j}_{\mathrm{c}})}{1!}\left\langle\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right\rangle_{\textbf{Non-
Int,c}}}$
$\displaystyle\frac{(-\beta^{j}_{\mathrm{c}})^{2}}{2!}\left\langle\phi^{2}\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right\rangle_{\textbf{Non-
Int,c}}+\text{Higher order terms}\,.$ (29)

We can substitute the internal potential, Eq. (25)

$\displaystyle\ln\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\frac{N^{{}^{\prime}}(N^{{}^{\prime}}-1)}{2V_{\mathrm{c}}}\int\mathrm{~{}d}^{3}\vec{r}\
\phi_{ij}({\vec{r}})$
$\displaystyle+\frac{N^{{}^{\prime}}(N^{{}^{\prime}}-1)}{2}\left[\int\frac{\mathrm{d}^{3}\vec{r}}{V_{\mathrm{c}}}\phi_{ij}^{2}(\vec{r})-\left(\int\frac{\mathrm{d}^{3}\vec{r}}{V_{\mathrm{c}}}\phi_{ij}(\vec{r})\right)^{2}\right]+\text{Higher
order terms}$ (30)

Using mathematical induction 28 for higher order terms, we can write,

$\displaystyle\ln\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})$
$\displaystyle=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\sum_{m=1}^{\infty}\frac{(-\beta)^{m}}{m!}\frac{N^{{}^{\prime}}(N^{{}^{\prime}}-1)}{2}\int\frac{\mathrm{d}^{3}\vec{r}}{V_{\mathrm{c}}}\phi_{ij}^{{m}}(\vec{r})+\mathcal{O}\left(\frac{{N^{{}^{\prime}}}^{3}}{V_{\mathrm{c}}^{2}}\right)$
$\displaystyle=\ln\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]+\frac{N^{{}^{\prime}}(N^{{}^{\prime}}-1)}{2V_{\mathrm{c}}}\int\mathrm{d}^{3}\vec{r}[\exp(-\beta\phi_{ij}(\vec{r}))-1]+\mathcal{O}\left(\frac{{N^{{}^{\prime}}}^{3}}{V_{\mathrm{c}}^{2}}\right)$
(31)

Terminating the expansion at order ${{N^{{}^{\prime}}}^{2}}/{V_{\mathrm{c}}}$
and applying the exponential function on both sides of Eq. (2.2), we can write

$\displaystyle\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\approx\left[{\frac{V_{\mathrm{c}}^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\right]\exp\left[\frac{N^{{}^{\prime}}(N^{{}^{\prime}}-1)}{2V_{\mathrm{c}}}\int\mathrm{d}^{3}\vec{r}[\exp(-\beta^{j}_{\mathrm{c}}\phi_{ij}(\vec{r}))-1]\right]\,.$
(32)

With Eq. (23), we finally obtain for the grand canonical partition function

$\displaystyle{\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}=\sum_{N^{{}^{\prime}}=0}^{\infty}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\approx\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\mathcal{C}_{2}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}$
(33)

in which

$\displaystyle\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\frac{1}{N^{{}^{\prime}}!}\left[{{V_{\mathrm{c}}}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3}}\right]^{N^{{}^{\prime}}}$
(34)

and

$\displaystyle\mathcal{C}_{2}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\exp\left[\frac{N^{{}^{\prime}}N^{{}^{\prime}}}{2\,V_{\mathrm{c}}}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})-\frac{N^{{}^{\prime}}}{2\,V_{\mathrm{c}}}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})\right]\,,$
(35)

where

$\displaystyle\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})=\int\mathrm{d}^{3}\vec{r}[\exp(-\beta^{j}_{\mathrm{c}}\phi_{ij}(\vec{r}))-1]\,.$
(36)

As we will see in the following section (precisely in Eq. (3)), we will
interpret the infinite sum over $N^{{}^{\prime}}$ as the series expansion of
the exponential function. To do this, we must eliminate the dependence on
$(N^{{}^{\prime}})^{2}$ in $\mathcal{C}_{2}$. We achieve this by replacing
$\mathcal{C}_{2}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})$ with
$\mathcal{C}_{2}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})$
defined as

$\displaystyle\mathcal{C}_{2}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}}):=\exp\left[\frac{{\left\langle
N_{\mathrm{c}}\right\rangle}N^{{}^{\prime}}}{2\,V_{\mathrm{c}}}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})-\frac{{\left\langle
N_{\mathrm{c}}\right\rangle}}{2\,V_{\mathrm{c}}}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})\right]$
(37)

so that we can define ${\mathcal{Q}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}$
as

$\displaystyle{\mathcal{Q}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}=\sum_{N^{{}^{\prime}}=0}^{\infty}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{Z}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})$
$\displaystyle=\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\mathcal{C}_{2}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}\,.$
(38)

We know that
${\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}$ is
an infinite series of strictly non-negative terms, but it is also the grand
canonical partition function that equals the phase space volume of a physical
system and hence it converges to a finite number. From calculus, we know that
the comparison test is a way of deducing the convergence of an infinite
series. Based on the comparison test, since
${\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}$
converges and the following inequality holds for all sufficiently large
$N^{{}^{\prime}}$,

$\displaystyle
0<{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\mathcal{C}_{2}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}<{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\mathcal{C}_{2}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})},$
(39)

then the infinite series ${\mathcal{Q}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}$
also converges and we can write

$\displaystyle{\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}=\xi\,{\mathcal{Q}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})},\quad\xi>0$
(40)

where $\xi$ is a scaling factor. By doing so,

$\displaystyle{\mathcal{Q}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}=\xi\,{\mathcal{Q}^{\left\langle
N_{\mathrm{c}}\right\rangle}(\beta^{j}_{\mathrm{c}},\mu^{j}_{\mathrm{c}},V_{\mathrm{c}})}$
$\displaystyle\propto\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\exp\left[\frac{\rho
N^{{}^{\prime}}}{2}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})-\frac{\rho}{2}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})\right]}\,,$
(41)

where $\rho$ is the number density of particles inside the shell

$\displaystyle\rho=\frac{{\left\langle
N_{\mathrm{c}}\right\rangle}}{V_{\mathrm{c}}}\,.$ (42)

Obviously $\mathcal{Q}$ and $\mathcal{Q}^{{\langle N_{\text{C}}\rangle}}$ are
related via a simple scaling factor that becomes irrelevant once the weighted
distribution is normalized at fixed $\beta$.

##  3 Derivation of a functional form for
$p(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})$

As discussed before, to fully determine the weight function according to Eq.
(2.1), we need to choose
$g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$.
By inspection of Eq. (2.2), we suggest the following ansatz,
$g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
can be written as

$\displaystyle
g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=\frac{\rho}{2}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})$
(43)

With this choice, $g_{2}$ can now be interpreted as a dimensionless variable
with the same order of magnitude as the average number of particles in the
core. Following Eq. (2.1), the weighted distribution function can be written
as

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle={\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)+\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\right)}$
$\displaystyle={\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),\frac{\rho}{2}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})+\ln({\frac{\rho}{2}\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})})\right)}\,.$
(44)

Here, we want to remind the reader that after choosing $g_{2}$, $g_{1}$ is
uniquely determined. To make our mathematical derivation easier to follow, we
continue with the notation $g_{2}$ that is now defined, and insert Eq. (43) in
the final step. Ideally, we would simply strive for an analytical expression
of Eq. (3). Unfortunately, this is not possible as we show in the Appendix and
we need another approximation to continue. Based on the fact that $g_{2}$ is
of the same order of magnitude as the average number of particles
${\left\langle N_{\mathrm{c}}\right\rangle}$ in the core and our starting
condition that the core is small with a volume comparable to that of the
particles, it is reasonable to conclude that the observed average number of
particles in the core is a small number. For $\lim{\left\langle
N_{\mathrm{c}}\right\rangle}\rightarrow 1$, we see that
$\displaystyle\lim_{{\left\langle N_{\mathrm{c}}\right\rangle}\rightarrow
1}{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\approx
1$. Now, approximating
$g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\approx\ln\left(g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\right)+1$
at the point
$(g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right))=1$,
we can write

$g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)+\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\approx
2\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+1\,.$

To extend the range of validity for the large average number of particles in
the core, we adjust the approximation by introducing two parameters

$\displaystyle
g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)+\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\approx\lambda_{1}\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+\lambda_{2}\,.$
(45)

As can be seen, one would have $\lambda_{1}=2$ and $\lambda_{2}=1$ in the
vicinity of the average number of particles being equal to one. We emphasize
that this parameterization is neither rigorous nor unique but a reasonable
choice with the correct low particle number limit that allows us to move on.
With this, Eq. (3) can be reformulated to give

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)={\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)+\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\right)}$
$\displaystyle\approx{\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),\lambda_{1}\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+\lambda_{2}\right)}$
$\displaystyle\propto\left[\sum_{N^{{}^{\prime}}=0}^{\infty}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\times\exp\left[\left(\lambda_{1}\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+\lambda_{2}\right)\,N^{{}^{\prime}}-\left(\lambda_{1}\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+\lambda_{2}\right)\right]\right]^{-1}\,.$
(46)

The artificial chemical potential $\mu^{j}_{\mathrm{c}}$ is defined in terms
of the artificial entropy $S^{j}_{\mathrm{c}}$ as

$\displaystyle\mu^{j}_{\mathrm{c}}\equiv-T^{j}_{\mathrm{c}}\left(\frac{\partial
S^{j}_{\mathrm{c}}}{\partial N^{{}^{\prime}}}\right)_{\phi,V_{\mathrm{c}}}\,,$
(47)

with volume $V_{\mathrm{c}}$, internal energy $\phi$, and the number of
particles $N^{{}^{\prime}}$. The Sackur–Tetrode equation 29, 30 expresses the
entropy $S$ of a monoatomic ideal gas in terms of its thermodynamic states

$\displaystyle
S^{j}_{\mathrm{c}}=N^{{}^{\prime}}k_{B}\left[\ln\left(\frac{V_{\mathrm{c}}}{N^{{}^{\prime}}}\left(\frac{4\pi
m\phi}{3N^{{}^{\prime}}h^{2}}\right)^{3/2}\right)+\frac{5}{2}\right]\,.$ (48)

Inserting Eq. (48) into Eq. (47) allows us to obtain the chemical potential
for a monoatomic ideal gas

$\displaystyle\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}=-\ln\left[\frac{V_{\mathrm{c}}}{N^{{}^{\prime}}}\left(\frac{2\pi
mk_{\text{B}}T}{h^{2}}\right)^{3/2}\right]=\ln{\rho{{\Lambda^{j}_{\mathrm{c}}}^{3}}}\,,$
(49)

where ${{\Lambda^{j}_{\mathrm{c}}}^{3}}=h/\sqrt{2\pi mk_{\mathrm{B}}T}$ is the
thermal de Broglie wavelength and $h$ is Planck’s constant. For an interacting
system of particles with no internal structure, i.e. no intramolecular degrees
of freedom, we can write31

$\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}=\ln{\rho\Lambda^{3}}+\beta^{j}_{\mathrm{c}}\mu^{j,{\mathrm{ex}}}_{\mathrm{c}}\,,$

with the excess chemical potential $\mu^{j,{\mathrm{ex}}}_{\mathrm{c}}$, which
accounts for intermolecular interactions between molecules and is given by the
potential distribution theorem 32, 33. Equivalently, the number density can be
written as

$\displaystyle\rho=\gamma\frac{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}{{{\Lambda^{j}_{\mathrm{c}}}^{3}}},\quad\gamma=\exp\,(-{\beta^{j}_{\mathrm{c}}\mu^{j,{\mathrm{ex}}}_{\mathrm{c}}})>0\,.$
(50)

Now we can write

$\displaystyle
g_{2}(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})=\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}},\quad
k=\frac{\mathcal{C}_{3}(\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})}{{2\Lambda^{j}_{\mathrm{c}}}^{3}}$
(51)

and

$\displaystyle\lambda_{1}\ln{g_{2}(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}})}+\lambda_{2}=\lambda_{1}\,\ln\left[{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\right]+\lambda_{2}=\lambda_{1}\,{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}+\lambda_{1}\,\ln{\gamma\,k}+\lambda_{2}$
(52)

With this we can finally write

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle\propto\left[\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\exp\left[(\lambda_{1}\,{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}+\lambda_{1}\,\ln{\gamma\,k}+\lambda_{2})\,N^{{}^{\prime}}-(\lambda_{1}\,{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}+\lambda_{1}\,\ln{\gamma\,k}+\lambda_{2})\right]}\right]^{-1}$
$\displaystyle\propto\left[\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\left[\,\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{N^{{}^{\prime}}}\left[\mathrm{e}^{\lambda_{2}}\,(\gamma\,k)^{\lambda_{1}}\right]^{N^{{}^{\prime}}}\left[\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{-1}\left[\mathrm{e}^{\lambda_{2}}\,(\gamma\,k)^{\lambda_{1}}\right]^{-1}}\right]^{-1}$
$\displaystyle\propto\left[\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\frac{1}{N^{{}^{\prime}}!}\left[{{V_{\mathrm{c}}}\left[\frac{\sqrt{2mk_{\text{B}}T^{j}_{\mathrm{c}}}}{h}\right]^{3}}\right]^{N^{{}^{\prime}}}\left[\,\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{N^{{}^{\prime}}}\left[\mathrm{e}^{\lambda_{2}}\,(\gamma\,k)^{\lambda_{1}}\right]^{N^{{}^{\prime}}}\left[\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{-1}\left[\mathrm{e}^{\lambda_{2}}\,(\gamma\,k)^{\lambda_{1}}\right]^{-1}}\right]^{-1}$
$\displaystyle\propto\left[\left[\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{-1}\left[\mathrm{e}^{\lambda_{2}}\,(\gamma\,k)^{\lambda_{1}}\right]^{-1}\sum_{N^{{}^{\prime}}=0}^{\infty}\frac{1}{N^{{}^{\prime}}!}\left[\underbrace{{{V_{\mathrm{c}}}\left[\frac{\sqrt{2mk_{\text{B}}T^{j}_{\mathrm{c}}}}{h}\right]^{3}}\,(\gamma\,k)^{\lambda_{1}}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\lambda_{2}}}_{\alpha}\right]^{N^{{}^{\prime}}}\right]^{-1}\,.$
(53)

Interpreting the infinite sum over $N^{{}^{\prime}}$ as the power series of
the exponential function,

$\displaystyle\exp{\alpha}:=\sum_{N^{{}^{\prime}}=0}^{\infty}\frac{\alpha^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}$
(54)

we arrive at

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\propto\,(\gamma\,k)^{\lambda_{1}}\,\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\frac{(\gamma\,k)^{\lambda_{1}}\,V_{\mathrm{c}}}{{\Lambda^{j}_{\mathrm{c}}}^{3}}}\mathrm{e}^{(\lambda_{1}+1)\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}+\lambda_{2}}}\,,$
(55)

and with

$\displaystyle\quad\alpha_{1}=\,(\gamma\,k)^{\lambda_{1}},\quad\alpha_{2}=\frac{(\gamma\,k)^{\lambda_{1}}\,V_{\mathrm{c}}}{{\Lambda^{j}_{\mathrm{c}}}^{3}},\quad\gamma>0\,$
(56)

a simplified notation gives

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\propto\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}+\lambda_{2}}}\,.$
(57)

Since $\alpha_{1}$ and $\alpha_{2}$ are explicit functions of
$\beta^{j}_{\mathrm{c}}$ (via $\gamma$), for the sake of simplicity, we
consider all artificial baths to have a fixed temperature equal to the thermal
bath surrounding the composite system, i.e. $\beta^{j}_{\mathrm{c}}=\beta$,
for all $j$ and only the artificial chemical potentials are allowed to
fluctuate. It is important to note that Eq. (57) is valid for any set of
artificial baths with any arbitrary combination of artificial thermodynamic
variables. We choose to fix the temperature for simplicity. In this case, Eq.
(10) must be adjusted to give

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)$
$\displaystyle\approx\sum^{j}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{\prime,j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu^{\prime,j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle=\int_{\mu^{\prime}_{\mathrm{c}}}\text{d}\mu^{\prime}_{\mathrm{c}}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V}\mid\mu^{\prime}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu^{\prime}_{\mathrm{c}},V_{\mathrm{c}}\right)\,.$
(58)

In order to preserve the probability distribution as a dimensionless function,
we should define $\mu^{\prime}_{\mathrm{c}}:=\beta\mu_{\mathrm{c}}$. So we can
equivalently write,

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu,\beta,V\right)$
$\displaystyle\approx\sum^{j}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\beta\mu^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\beta\mu^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle=\beta\int_{\mu_{\mathrm{c}}}\text{d}\mu_{\mathrm{c}}p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V}\mid\beta\mu_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\beta\mu_{\mathrm{c}},V_{\mathrm{c}}\right)$
(59)

in which

$\displaystyle
p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\beta\mu^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=p\left(\\{\vec{\mathbf{q}},\vec{\mathbf{p}}\\}_{V_{\mathrm{c}}}\mid\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)$
(60)

and

$\displaystyle
p\left(\beta\mu^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)\,.$
(61)

For the remaining derivations in the fixed temperature scenario we use the
following integration convention,

$\displaystyle\sum^{j}\Longrightarrow\iint_{\mu_{\mathrm{c}}\beta_{\mathrm{c}}}\text{d}\mu_{\mathrm{c}}\text{d}\beta_{\mathrm{c}}$
(62)

which is adjusted to

$\displaystyle\sum^{j}\Longrightarrow\beta\int_{\mu_{\mathrm{c}}}\text{d}\mu_{\mathrm{c}}\,.$
(63)

The weighted distribution function can be written as

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}}}{\sum_{j}\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}}}\,,$
(64)

and the normalization factor of $p(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}})$
can be written as

$\displaystyle{\sum_{j}\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}}}=\beta\int_{-\infty}^{\infty}\text{d}\mu_{\mathrm{c}}\,\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta\mu_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu_{\mathrm{c}}+\lambda_{2}}}\,.$
(65)

We impose a change of variables
$(\lambda_{1}+1)\beta\mu_{\mathrm{c}}+\lambda_{2}=-u$, such that
$\text{d}\mu_{\mathrm{c}}=-{\text{d}u}/{(\lambda_{1}+1)\beta}$

$\displaystyle\beta\int_{-\infty}^{\infty}\text{d}\mu_{\mathrm{c}}\,\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta\mu_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu_{\mathrm{c}}+\lambda_{2}}}=\frac{\beta}{(\lambda_{1}+1)\,\beta}\int_{-\infty}^{\infty}\text{d}u\,\alpha_{1}\,\mathrm{e}^{-\frac{\lambda_{1}\,u}{\lambda_{1}+1}+\frac{\lambda_{2}}{\lambda_{1}+1}}\mathrm{e}^{-\alpha_{2}\mathrm{e}^{-u}}$
$\displaystyle=\frac{\beta}{(\lambda_{1}+1)\,\beta}\int_{-\infty}^{\infty}\text{d}u\,\alpha_{1}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{-\frac{\lambda_{1}}{\lambda_{1}+1}}\alpha_{2}^{\frac{\lambda_{1}}{\lambda_{1}+1}}\,\mathrm{e}^{-\frac{\lambda_{1}\,u}{\lambda_{1}+1}}\mathrm{e}^{-\alpha_{2}\mathrm{e}^{-u}}$
$\displaystyle=\frac{\beta}{(\lambda_{1}+1)\beta}\int_{-\infty}^{\infty}\text{d}u\,\alpha_{1}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{-\frac{\lambda_{1}}{\lambda_{1}+1}}(\alpha_{2}\,\mathrm{e}^{-{u}})^{\frac{\lambda_{1}}{\lambda_{1}+1}}\mathrm{e}^{-\alpha_{2}\mathrm{e}^{-u}}\,.$
(66)

Another change of variables $\alpha_{2}\mathrm{e}^{-u}=x$, such that
$\text{d}u=-{\text{d}x}/{\alpha_{2}x}$ gives

$\displaystyle\frac{\beta}{(\lambda_{1}+1)\beta}\int_{-\infty}^{\infty}\text{d}u\,\alpha_{1}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{-\frac{\lambda_{1}}{\lambda_{1}+1}}(\alpha_{2}\,\mathrm{e}^{-{u}})^{\frac{\lambda_{1}}{\lambda_{1}+1}}\mathrm{e}^{-\alpha_{2}\mathrm{e}^{-u}}$
$\displaystyle=\frac{\beta}{(\lambda_{1}+1)\beta}\alpha_{1}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\int_{0}^{\infty}\text{d}x\,x^{\frac{\lambda_{1}}{\lambda_{1}+1}-1}\mathrm{e}^{-x}=\frac{\beta}{(\lambda_{1}+1)\beta}\alpha_{1}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})\,,$
(67)

and we finally arrive at

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{\alpha_{1}\,\mathrm{e}^{\lambda_{1}\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}}}{\frac{\beta}{(\lambda_{1}+1)\beta}\alpha_{1}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})}=\frac{{(\lambda_{1}+1)\beta}\,\mathrm{e}^{\lambda_{1}\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}}}{{\beta}\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})}\,.$
(68)

And we finally obtain a valid analytic expression for the weighted
distribution function of the artificial thermodynamic baths with three
parameters $\lambda_{1},\lambda_{2}$ and $\alpha_{2}$

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{{(\lambda_{1}+1)}\,\mathrm{e}^{\lambda_{1}\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu^{j}_{\mathrm{c}}+\lambda_{2}}}}{\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})}\,$
(69)

where

$\displaystyle\alpha_{2}=\frac{(\gamma\,k)^{\lambda_{1}}\,V_{\mathrm{c}}}{{\Lambda}^{3}},\quad\gamma>0,\quad
k=\frac{\mathcal{C}_{3}(\beta,V_{\mathrm{c}})}{2{\Lambda}^{3}},\quad\lambda_{1}>0,\quad\Lambda=\frac{h}{\sqrt{2\pi
mk_{\mathrm{B}}T}}$
$\displaystyle\mathcal{C}_{3}(\beta,V_{\mathrm{c}})=\int\mathrm{d}^{3}\vec{r}[\exp(-\beta\phi_{ij}(\vec{r}))-1]\,.$

The purpose of this study is to capture the properties of an interacting
system in an explicitly simulated core region conditional to the thermodynamic
constraints of the bath. The computational low-cost strategy was proposed in
Eq. (10). Now, by starting from $p(N_{\mathrm{c}}\mid\mu,\beta,V)$ and
following the same steps as in Eqns. (6-10) the counterpart of Eq. (10) can be
written as

$\displaystyle
p\left(N_{\mathrm{c}}\mid\mu,\beta,V\right)\approx\sum_{j}p\left(N_{\mathrm{c}}\mid\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\,.$
(70)

In the following we will establish a workflow that allows us to accurately
capture the probability distribution of the number of particles inside the
core $p(N_{\mathrm{c}}\mid\mu,\beta,V)$ from grand canonical simulations. For
sake of simplicity, we again only allow artificial chemical potentials
$\mu^{j}_{\mathrm{c}}$ to change while we keep $\beta^{j}_{\mathrm{c}}$
constant and equal to $\beta$. Eq. (70) can be expanded by substituting Eq.
(69)

$\displaystyle
p\left(N_{\mathrm{c}}\mid\mu,\beta,V\right)\approx\beta\int_{-\infty}^{\infty}\,\text{d}\,\mu_{\mathrm{c}}\,p\left(N_{\mathrm{c}}\mid\mu_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)\,\frac{{(\lambda_{1}+1)}\,\mathrm{e}^{\lambda_{1}\,\beta\mu_{\mathrm{c}}+\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\beta\mu_{\mathrm{c}}+\lambda_{2}}}}{\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})}\,.$
(71)

We use a simplified parametric notation for the right hand side of Eq. (71) as
follows

$\displaystyle p\left(N_{\mathrm{c}}\mid\mu,\beta,V\right)\approx
p\left(N_{\mathrm{c}}\mid\beta,V_{\mathrm{c}},\alpha_{2},\lambda_{1},\lambda_{2}\right)\,.$
(72)

As we discussed in the introduction, in order to find
$p(N_{\mathrm{c}}\mid\mu,\beta,V)$, a trivial albeit computationally expensive
way is to simulate a system where a large fraction of the bath is included in
addition to the small core. We can treat the combined system as a micro-
canonical, canonical or grand canonical system under an appropriate
thermodynamic constraint. Then, by using Monte Carlo simulations, we can
sample the probability distribution of the number of particles in the core.
Alternatively, with our newly proposed framework that culminates in Eq. (69),
we have mapped the problem — provided the parameters $\alpha_{2}$,
$\lambda_{1}$ and $\lambda_{2}$ have been determined appropriately — to a
computationally cheaper method.

We will show in the following that we can determine a set of parameters that
leads to accurate results for a hard-sphere model system and a Lennard–Jones
fluid. To optimize these parameters, we opted to perform canonical reference
calculations on a very large system. Obviously, this route only leads to a
computationally cheaper method if the parameters are transferable for similar
situations (e.g. core volume, interaction strengths), such that the expensive
reference simulations are rarely required. We will come back to this point in
the discussion. For now, we describe how we optimized the parameters based on
reference data and demonstrate that our model is indeed capable to reproduce
the reference data from core-only simulations.

To generate our reference data, we run a canonical Monte Carlo simulation for
the combined system, i.e. bath+shell+core and then sample the probability
distribution of the number of particles in the core, termed
${p_{\text{ref}}(N_{\mathrm{c}}\mid\mu,\beta,V)}$. We then have to fit the
parametric model
$p(N_{\mathrm{c}}\mid\beta,V_{\mathrm{c}},\alpha_{2},\lambda_{1},\lambda_{2})$
to match ${p_{\text{ref}}(N_{\mathrm{c}}\mid\mu,\beta,V)}$. For a given
reference distribution
$({p_{\text{ref}}({N^{1}_{\mathrm{c}}})},\ldots,{p_{\text{ref}}({N^{k}_{\mathrm{c}}})})$,
where $({N^{1}_{\mathrm{c}}},\ldots,{N^{k}_{\mathrm{c}}})$ is a set of $k$
different numbers of particles inside the core, a class of candidates for the
parametric model can be defined as

$\displaystyle{p}({N^{i}_{\mathrm{c}}}\mid\beta,V_{\mathrm{c}},\theta)={p}\left({N^{i}_{\mathrm{c}}}\mid\beta,V_{\mathrm{c}},\alpha_{2},\lambda_{1},\lambda_{2}\right)\quad\forall
i=1,\ldots,k,$
$\displaystyle\theta=(\alpha_{2},\lambda_{1},\lambda_{2})\in\Theta,\quad\text{s.t.}\quad\Theta=\\{(\alpha_{2},\lambda_{1},\lambda_{2})\mid\alpha_{2}>0,\,\lambda_{1}>0,\,\lambda_{2}\in\mathbb{R}\\}\,,$
(73)

where the structure of the function
${p}({N^{i}_{\mathrm{c}}}\mid\beta,V_{\mathrm{c}},\theta)$ is completely
determined by Eq. (71) and the $3$-dimensional parameter vector $\theta$. A
reasonable and quite intuitive way to fit the parametric model to the
reference data is to use a quadratic approximation criterion. The average
squared approximation error w.r.t. the observed reference

$\displaystyle\delta(\theta):=\delta(\alpha_{2},\lambda_{1},\lambda_{2})$
$\displaystyle=\sum_{i=1}^{k}\left[{\log
p_{\text{ref}}({N^{i}_{\mathrm{c}}}\mid\mu,\beta,V)}-\log{p}({N^{i}_{\mathrm{c}}}\mid\beta,V_{\mathrm{c}},\theta)\right]^{2}$
(74)

will be used to quantify the difference between the proposed model and
canonical simulation results.

Note that in the simulation section, we denote $\delta$ as the error between
our proposed model, Eq. (71) – with optimized parameters
$\alpha_{2},\lambda_{1},\lambda_{2}$ – and the canonical reference
calculations, whereas we show the discrepancy between the canonical
simulations and the single-bath grand canonical simulations with $\Delta$ that
is defined as

$\displaystyle\Delta:=\sum_{i=1}^{k}\left[{\log
p_{\text{ref}}({N^{i}_{\mathrm{c}}}\mid\mu,\beta,V)}-\log{p}_{\text{GC}}({N^{i}_{\mathrm{c}}}\mid\mu^{*}_{\mathrm{c}},\beta,V_{\mathrm{c}})\right]^{2}\,,$
(75)

where
${p}_{\text{GC}}({N^{i}_{\mathrm{c}}}\mid\mu^{*}_{\mathrm{c}},\beta,V_{\mathrm{c}})$
is the single-bath grand canonical probability distribution with the chemical
potential $\mu^{*}_{\mathrm{c}}$ that is defined as

$\mu^{*}_{\mathrm{c}}:=\underset{\mu^{j}\in[\mu_{\text{min}},\mu_{\text{max}}]}{\operatorname{argmin}}\overbrace{\bigg{|}\left\langle{N^{i}_{\mathrm{c}}}\,p_{\text{ref}}({N^{i}_{\mathrm{c}}}\mid\mu,\beta,V)\right\rangle-\left\langle{N^{i}_{\mathrm{c}}}\,{p}_{\text{GC}}({N^{i}_{\mathrm{c}}}\mid\mu^{j},\beta,V_{\mathrm{c}})\right\rangle\bigg{|}}^{W(\mu^{j})}\,,$
(76)

where ${\operatorname{argmin}}$ stands for the value of the chemical potential
in the domain of integration $[\mu_{\text{min}},\mu_{\text{max}}]$ that
produces the minimum value of $W(\mu^{j})$ over the the domain of integration.
Note that the bracket $\left\langle\dots\right\rangle$ indicates the average
number of the particles inside the core.

##  4 Simulations

Algorithm 1 Canonical Monte Carlo

$\bullet$ Set up the particles in some initial configuration $\mathbf{r}$.

$\bullet$ pick a maximum displacement vector $\left(d_{x},d_{y},d_{z}\right)$.

$\bullet$ Initialize the Monte Carlo step counter $n_{\text{MCS}}=0$ and
choose the maximum number of Monte Carlo steps $n_{\max}$.

while $n_{\mathrm{MCS}}<n_{\max}$ do

$\bullet$ From all $N$ particles, choose particle $i$ at random.

$\bullet$ Calculate the energy of the current system
$\phi_{\text{ini}}(\mathbf{r})$.

$\bullet$ Generate three uniform random numbers $u_{1},u_{2},u_{3}$ between
$-{1}$ and ${1}$.

$\bullet$ Through random displacement,
$\Delta\mathbf{r}=\left(u_{1}d_{x},u_{2}d_{y},u_{3}d_{z}\right)$,

shift atom $i$ to its new location
$\mathbf{r}_{i}=\mathbf{r}_{i}+\Delta\mathbf{r}$.

$\bullet$ Calculate the energy of the system after the movement
$U_{\text{fin}}(\mathbf{r})$.

$\bullet$ Find the difference in potential energy
$\Delta\phi=\phi_{\text{fin}}-\phi_{\text{ini}}$,

caused by the displacement of atom $i$,

then the Metropolis algorithm’s fundamental principle

is that a motion is accepted or rejected based on.

if $\Delta\phi\leq 0$ then

$\bullet$ Accept the state

else if $\Delta\phi>0$ then

$\bullet$ Generate a random number $k$ between 0 and 1.

if $k<e^{-\Delta\phi/k_{\text{B}}T}$ then

$\bullet$ Accept the state.

else if $k>e^{-\Delta\phi/k_{\text{B}}T}$ then

$\bullet$ The trial is rejected, the configuration that

existed before the rejected move is returned:
$\mathbf{r}_{i}=\mathbf{r}_{i}-\Delta\mathbf{r}$.

$\bullet$ Count the number of particles in the core for positions
$\mathbf{r}$, $N_{\text{c}}$.

$\bullet$ Set $n_{\mathrm{MCS}}=n_{\mathrm{MCS}}+1$

$\bullet$ Calculate probability distribution of the number of particles in the
core based on the histogram for all $N_{\text{c}}$.

Algorithm 2 Grand Canonical Monte Carlo

$\bullet$ Initialize the Monte Carlo step counter $n_{\text{MCS}}=0$ and
choose the maximum number of maximum Monte Carlo steps $n_{\max}$.

while $n_{\mathrm{MCS}}<n_{\max}$ do

$\bullet$ Select a uniform random number $\eta_{1}$ between 0 and 1 .

if $\eta_{1}<0.5$ (particle removal attempt) then

$\bullet$ Select a uniform random number $\eta_{2}$ between 0 and 1

$\bullet$ perform the Metropolis test given by  $\pi(N\rightarrow
N-1)=\operatorname{Min}\left(1,\frac{\Lambda^{3}N}{V}\exp(-\beta(\mu+\phi(N-1)-\phi(N)))\right)$

if $\eta_{2}<\pi(N\rightarrow N-1)$ then

$\bullet$ Remove the particle. Otherwise, keep the old configuration.

else if $\eta_{1}>0.5$ (particle insertion attempt) then

$\bullet$ Select a uniform random number $\eta_{2}$ between 0 and 1 .

$\bullet$ perform the Metropolis test given by  $\pi(N\rightarrow
N+1)=\operatorname{Min}\left(1,\frac{V}{\Lambda^{3}(N+1)}\exp(\beta(\mu-\phi(N+1)+\phi(N)))\right)$

if $\eta_{2}<\pi(N\rightarrow N+1)$ then

$\bullet$ Add the particle.

else if $\eta_{2}>\pi(N\rightarrow N+1)$ then

$\bullet$ keep the old configuration.

$\bullet$ Count the number of particles in the core.

$\bullet$ Set $n_{\mathrm{MCS}}=n_{\mathrm{MCS}}+1$

$\bullet$ Calculate probability distribution of the number of particles in the
core based on the histogram for all $N_{\text{c}}$.

Having derived an expression for the particle number distribution in the core,
$p(N_{\text{c}}|\mu,\beta,V)$ (see Eq. (71)), we will now demonstrate that for
a proper choice of parameters this expression is indeed flexible enough to
accurately convert the grand canonical simulation results to those obtained
from a canonical simulation with a large simulation box for an interacting
system. Only if this is achieved, we can obtain accurate results for the small
system of interest at the low computational cost of a grand canonical
simulation of the core system alone.

As a first check, we will test our approach for the case of a non-interacting
hard-sphere model. This system has been studied before by Dixit et al.24 and
we compare to their results in the following. We simulate the system at the
reduced density $8\rho r^{3}_{\text{p}}=\rho^{\prime}=0.9$ and a fixed radius
for the spherical core or solvation shell $R=2.2r_{\text{p}}$, where
$r_{\text{p}}$ is the radius of the hard-sphere particle. In all simulations
we fix one particle in the center of the core and count the number of
additional particles in the solvation shell or core. For the reference
canonical Monte-Carlo calculations we choose a cubic box with an edge length
of

$\displaystyle d=(8Nr^{3}_{\text{p}}/\rho^{\prime})^{1/3}\,,$ (77)

employing periodic boundary conditions and $N=100$.

We obtain reference canonical simulation results employing a modified version
of the Metropolis algorithm 1, where acceptance of a random Monte-Carlo step
is only possible if there is no overlap between the particles. A similar
modification is used for the grand canonical simulations following algorithm
2, since no inner energy $\phi$ is associated with the configurations of the
hard-sphere model.

Figure 2:  Particle number distributions for reduced density
$\rho^{\prime}=0.9$ of a three-dimensional hard sphere liquid. The blue dots
show the results for a reference canonical Monte Carlo simulation (MC), a
grand canonical Monte Carlo simulation for a single bath is shown in black
(GC) and results for the improved grand canonical description with coupling to
multiple baths is shown in green (Dixit, Ref. 24) and red (this work). The
parameters above the graph show the optimized parameters for the model
suggested in this work, the corresponding entropy is denoted by $S$ and
$\Delta$ denotes the discrepancy between the reference calculation and the
single-bath grand canonical simulation according to Eq. (75). The error
between our proposed model, Eq. (71) – with optimized parameters
$\alpha_{2},\lambda_{1},\lambda_{2}$ – and the reference calculations has been
denoted by $\delta$ according to Eq. (74) .

In Figure 2 we show the particle number distribution obtained from the
canonical reference calculations (blue dots), the standard grand canonical
simulation assuming a single bath (black dots), the corrected prediction
according to Dixit et al. (green dots) and the corrected prediction according
to Eq. (71) (red dots). We implemented a numerical integration scheme to
approximate the integral in Eq. (71). To do this, we had to transform the
domain of integration in $\mu_{\text{c}}$ from an indefinite to a definite
one, $\mu^{j}\in[\mu_{\text{min}},\mu_{\text{max}}]$. The integrand has to be
evaluated at a finite set of integration points i.e. chemical potentials and a
weighted sum of these values is employed to approximate the integral via the
trapezoidal rule. We carefully checked convergence with respect to the bounds
of the integral and the intervals applied for the numerical integration for
each individual calculation. Note that the integration points and weights rely
on the specific scheme employed and the accuracy demanded from the
calculation. As can be seen in Eq. (71), the integrand comprises two parts,
the first part is
$p\left(N_{\text{c}}\mid\mu^{j}_{\text{c}},\beta,V_{\text{c}}\right)$, which
is evaluated with Metropolis grand canonical Monte Carlo sampling at each
integration point, Algorithm 2, and the second part is the weight function
that can be analytically evaluated as a function of the parameters
$\theta=(\alpha_{2},\lambda_{1},\lambda_{2})$. Now, $\alpha_{2}$,
$\lambda_{1}$ and $\lambda_{2}$ can be determined such that they minimize the
error between $\log{p_{\text{ref}}({N^{i}_{\mathrm{c}}}\mid\mu,\beta,V)}$ and
$\log{p}({N^{i}_{\mathrm{c}}}\mid\beta,V_{\mathrm{c}},\alpha_{2},\lambda_{1},\lambda_{2})$,
Eq. (74). We employed a similar strategy to obtain optimal parameters for the
model of Dixit et al.24 (green dots).

While there is a discrepancy of $\Delta=21.17$, Eq. (75), between the single
bath grand canonical result and the reference simulations, both optimized
multi-bath descriptions (green and red dots) are perfectly capable to
reproduce the reference results, reducing the squared error to only
$\delta=2.30$. That our approach produces roughly the same result as the one
of Dixit et al. is to be expected since for the choice $\alpha_{2}=1.0$ the
functional form of our model equals that of Dixit et al.24 as we show in the
following.
Their weighted distribution function reads

$\displaystyle
P(\mu)=\frac{\mathrm{e}^{\left(-\lambda_{1}^{\prime}\,\mathrm{e}^{-\mu^{\prime}}-\lambda_{2}^{\prime}\mu^{\prime}\right)}}{\Gamma}\left(\lambda_{2}^{\prime}\right)\lambda_{1}^{{\prime}^{-\lambda_{2}^{\prime}}}=\frac{\mathrm{e}^{-\lambda_{1}^{\prime}\mathrm{e}^{-\mu^{\prime}}}\mathrm{e}^{-\lambda_{2}^{\prime}\mu^{\prime}}}{\Gamma\left(\lambda_{2}^{\prime}\right)\lambda_{1}^{{\prime}^{-\lambda_{2}^{\prime}}}}$
(78)

in which $\mu^{\prime}$ is the chemical potential in units of $\beta$. Writing
our proposed model, Eq. (69), in the units of Dixit et al.
($\mu^{\prime}=\beta\mu_{c}^{j}$) we obtain

$\displaystyle
p\left(\mu^{\prime},\beta,V_{\mathrm{c}}\right)=\frac{{(\lambda_{1}+1)}\,\mathrm{e}^{\lambda_{1}\,\mu^{\prime}}\mathrm{e}^{\lambda_{2}}\mathrm{e}^{{-\alpha_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\mu^{\prime}}\mathrm{e}^{\lambda_{2}}}}{\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}\alpha_{2}^{\frac{-2\,\lambda_{1}-1}{\lambda_{1}+1}}\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})}\,.$
(79)

Setting $\alpha_{2}=1.0$ considerably simplifies the expression to

$\displaystyle
p\left(\mu^{\prime},\beta,V_{\mathrm{c}}\right)=\frac{{(\lambda_{1}+1)}\mathrm{e}^{\lambda_{2}}\,\mathrm{e}^{-\mathrm{e}^{\lambda_{2}}\mathrm{e}^{(\lambda_{1}+1)\,\mu^{\prime}}}\,\mathrm{e}^{\lambda_{1}\,\mu^{\prime}}}{\Gamma(\frac{\lambda_{1}}{\lambda_{1}+1})\,\mathrm{e}^{\frac{\lambda_{2}}{\lambda_{1}+1}}}\,.$
(80)

With the following change of variables

$\lambda_{1}+1=A\quad\text{and}\quad\mathrm{e}^{\lambda_{2}}=B\,,$

we can write

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{A\,B\,\mathrm{e}^{-B\mathrm{e}^{A\,\mu^{\prime}}}\,\mathrm{e}^{(A-1)\,\mu^{\prime}}}{\Gamma(\frac{A-1}{A})\,B^{\frac{1}{A}}}\,$
(81)

and arrive at the same functional form as in Eq. (78). However, it is not
possible to map our set of parameters to those of Dixit et al. In the
following we will only present results obtained with our model and note here
that the model in Ref. 24 is capable to provide results of similar quality
also for an interacting system. In practice, however, we found it slightly
easier to optimize the parameters of our model due to well-defined boundaries
and a smaller range of values.

We now challenge our approach with an interacting Lennard–Jones fluid. For
such a system, we can modulate the reduced interaction strength
$\epsilon^{\prime}$ (and the distance at the potential minimum $\sigma$) in
addition to the reduced density $\rho^{\prime}$. We again obtain the optimal
parameters for our model by minimizing the error function $\delta$ defined in
Eq. (74). To understand how different choices of the parameters affect the
model performance and in order to provide a reliable initial guess for
subsequently applied standard numerical optimization algorithms, we prepared
contour plots of the sum-of-squares error function $\delta$ for a 2D-scan of
the $\lambda_{1}$ and $\alpha_{2}$ parameters, keeping $\lambda_{2}$ fixed. We
varied $\lambda_{1}$ from 0 to 6 in intervals of 0.1, and $\alpha_{2}$ from 0
to 2 in intervals of 0.1. Contour plots where then created for values of
$-100\leq\lambda_{2}\leq 10$ in steps of 10.

Figure 3: The upper part of panel A shows the contour plot of the sum-of-
squares error function $\delta$ according to Eq. (74), for a 2D-scan of
$\lambda_{1}$ varied between 0 and 6 in intervals of 0.1 and $\alpha_{2}$
varied between 0 and 2 in intervals of 0.1, keeping $\lambda_{2}=-80$ fixed
for a constant well depth $\epsilon^{\prime}=1.8$ of a Lennard–Jones liquid
with a reduced density of $\rho^{\prime}=0.9$. In the upper part of panel B we
show the contour plot of the sum-of-squares error function $\delta$, for a
2D-scan of $\lambda_{1}$ varied between 0 and 2 in intervals of 0.1, and
$\lambda_{2}$ varied between -100 and 10 in intervals of 10, keeping
$\alpha_{2}=1$ fixed for the same Lennard–Jones liquid. The lower panels are
the corresponding particle number distributions. The blue dots show the
results for a reference canonical Monte Carlo simulation (MC), a grand
canonical Monte Carlo simulation for a single bath is shown in black (GC) and
results for the improved grand canonical description with coupling to multiple
baths are shown in red dots (panel A, three-parameter optimization) and the
red dashed line (panel B, two-parameter optimization, $\alpha_{2}=1.0$). The
parameters above each graph show the optimized parameters for the model
suggested in this work for the initial guesses
$(\alpha_{2}=1.51,\lambda_{1}=5.01,\lambda_{2}=-80)$ and
$(\alpha_{2}=1,\lambda_{1}=1.51,\lambda_{2}=-33)$. The corresponding entropy
is denoted by $S$ and $\delta$ denotes the error between our proposed model,
Eq. (71) – with optimized parameters $\alpha_{2},\lambda_{1},\lambda_{2}$ –
and the reference calculations according to Eq. (74).

The upper part of panel A in Figure 3 shows an exemplary result of such
calculations for $\lambda_{2}=-80.0$ for a reduced density of
$\rho^{\prime}=0.9$ and $\epsilon^{\prime}=1.8$. The parameter $\sigma$ was
fixed to a value of one for all simulations in the current work. Similarly,
the upper part of panel B in Figure 3 shows a contour plot of the sum-of-
squares error function $\delta$ for a 2D-scan of the $\lambda_{1}$ and
$\lambda_{2}$ parameters, keeping $\alpha_{2}=1$ fixed for the same set of
$\rho^{\prime}$ and $\epsilon^{\prime}$. Keeping $\alpha_{2}$ fixed at 1.0 is
motivated by the analogy to the approach of Dixit et al.24 and by the fact
that the optimized value for $\alpha_{2}$ is close to 1.0 for all systems
studied in this work (vide infra). We varied $\lambda_{1}$ from 0 to 2 in
intervals of 0.1, and $\lambda_{2}$ from -100 to 10 in intervals of 10. We
observe structured contours and substantially different local minima for both
scans (red dots). These local minima then served as initial guesses for the
numerical parameter optimization for which we chose a conjugate gradient
algorithm with numerical gradients. It is worth mentioning that for panels A
and B in Figure 3, we implemented the minimization algorithm for a three-
parameter optimization problem and a two-parameter (keeping $\alpha_{2}=1$,
fixed) optimization problem, respectively. The lower panels of Figure 3 show
the particle number distribution function obtained for the optimized set of
parameters. We see that we arrive at the approximately same set of
corresponding optimized parameters within numerical precision of the
optimization algorithm and hence quasi-identical particle number distribution
functions and residual errors $\delta$ (lower graphs of both panels in Figure
3) no matter if we keep $\alpha_{2}$ fixed at 1.0 or optimize it. We can
further characterize the particle number distributions with an entropy

$\displaystyle S=-\sum_{i=1}^{k}p\ln p$ (82)

and find approximately identical entropies for both simulations. If we
consider the reference distribution to be the accurate equilibrium
distribution, then we can interpret $\delta$ as the error between the
distribution of our proposed model and the equilibrium distribution. A smaller
$\delta$ can be interpreted as how much our model distribution converges
toward the equilibrium distribution. However, from a statistical mechanics
perspective, the entropy attributed to the equilibrium distribution should be
at its maximum and any slight difference from the equilibrium distribution
should lead to a decrease in entropy. One interesting observation is that at
panel A we have higher entropy compared to panel B, which is in good agreement
with observing a lower $\delta$ for panel A compared to panel B. Despite these
clearly rugged hypersurfaces it is obviously possible to optimize the three
parameters or two parameters (keeping $\alpha_{2}=1$) of our model with
standard optimization schemes and arrive at a global minimum if a suitable
initial guess is provided. As expected, the excellent agreement between the
reference simulations and our model can not be obtained if only a single bath
is assumed for the grand canonical.

We investigate the performance of these different models further in the
following for different reduced densities $\rho^{{}^{\prime}}$ and interaction
strengths $\epsilon^{\prime}$.

Figure 4: Particle number distributions for different reduced potential well
depths $\epsilon^{\prime}=0.2,1.0$ and $1.8$ (panels A, B and C, respectively)
and constant reduced density $\rho^{\prime}=0.9$ of a Lennard–Jones liquid.
The blue dots show the results of a reference canonical Monte Carlo simulation
(MC), a grand canonical Monte Carlo simulation for a single bath is shown in
black (GC) and results for the improved grand canonical description with
coupling to multiple baths are shown in red dots and red dashed line (this
work with three parameters optimization and this work with two parameters
optimization keeping $\alpha_{2}=1$ fixed, respectively). The optimized
parameters for the model suggested in this work (with three parameters
optimization) are shown above each panel along with the corresponding entropy
$S$ and $\Delta$ which denotes the discrepancy between the reference
calculation and the single-bath grand canonical simulation according to Eq.
(75).

Figure 4 shows the probability distribution of the number of particles in the
solvation shell for different reduced depths of the Lennard–Jones potential
well $\epsilon^{\prime}=(0.2,1.0,1.8)$, and a fixed reduced density
$\rho^{\prime}=0.9$. We observe that with increasing $\epsilon^{\prime}$, the
discrepancy $\Delta$ between the single-bath grand canonical and the canonical
Monte Carlo simulations becomes more pronounced. In all cases, however, we are
able to identify a set of parameters for which the multi-bath model we propose
here results in excellent agreement with the reference simulation.

Figure 5: Particle number distributions for constant reduced well depth
$\epsilon^{\prime}=1.0$ of a Lennard–Jones liquid for a reduced density of
$\rho^{\prime}=0.7$ (panel A) and $\rho^{\prime}=0.8$ (panel B). The blue dots
show the results for a reference canonical Monte Carlo simulation (MC), a
grand canonical Monte Carlo simulation for a single bath is shown in black
(GC) and results for the improved grand canonical description with coupling to
multiple baths are shown in red dots and red dashed line (this work with three
parameters optimization and this work with two parameters optimization keeping
$\alpha_{2}=1$ fixed, respectively). The optimized parameters for the model
suggested in this work are shown above the graphs, along with the
corresponding entropy $S$ and $\Delta$, which denotes the discrepancy between
the reference calculation and the single-bath grand canonical simulation
according to Eq. (75).

In Figure 5 we show the probability distribution of the number of particles in
the solvation shell, for different reduced densities
$\rho^{\prime}=(0.7,0.8)$, and a fixed reduced depth of the Lennard–Jones
potential well $\epsilon^{\prime}=1.0$. We observe that the smaller
$\rho^{\prime}$, the smaller the discrepancies between the single-bath grand
canonical and the canonical Monte Carlo reference simulation.

Figures 4 and 5 illustrate how the density and the potential well depth affect
the difference between the solvent-shell particle number distribution of the
single-bath grand canonical simulation and the canonical reference
simulations. The single-bath grand canonical description for a small solvation
shell is increasingly inaccurate for higher densities and more attractive
potentials suggesting the average pairwise internal interaction between two
particles to be the relevant property. For a system in thermal equilibrium
with a heat bath at a fixed temperature, we can calculate the classical
average inter-particle distance $l$ as

$\displaystyle l\propto\rho^{-\frac{1}{3}}\,.$ (83)

Since the magnitude of a pairwise force $F_{r}$ experienced by a particle can
be defined as the negative of the derivative of the internal potential energy
function with respect to the inter-particle distance $r$, the average pairwise
interaction $\left\langle F_{r}\right\rangle$ for a system at thermal
equilibrium can be calculated as

$\displaystyle\left\langle
F_{r}\right\rangle=F_{r}(l)=-\frac{\partial\phi_{\text{LJ}}}{\partial{r}}\bigg{|}_{r=l}=-\zeta\bigg{|}_{r=l}\,,$

where $\zeta$ is the slope of the internal potential energy curve.

Figure 6: The pairwise Lennard–Jones potential energy $\phi_{\text{LJ}}$ as a
function of the distance $r$ between a pair of particles for different reduced
potential well depth $\epsilon^{\prime}=0.2,1.0$ and $1.8$ (curves blue, black
and red, respectively) and constant $\sigma=1$ that is the distance at which
the particle-particle potential energy is zero. $l$ is the classical inter-
particle spacing that is controlled by the density. $\zeta_{i}^{j}$ is the
slope of the internal potential energy curve at a specific density $j$ and a
specific depth of potential well $i$, such that its negative value equals the
force on the particle.

Figure 6 shows the dependence of the average pairwise interaction in a
Lennard–Jones system with density and depth of potential well. For a fixed
density, increasing the depth of the potential well increases the magnitude of
the average pairwise interaction, and for a fixed depth of the potential well,
increasing the density increases the average pairwise interaction. In light of
the previous observations we find that weaker pairwise interactions lead to a
smaller discrepancy between the single-bath grand canonical simulation and
canonical reference calculation, making the former a more suitable
approximation in these cases.

##  5 Discussion and outlook

We present a model for the simulation of the thermodynamic properties of a
small explicit core region — or solvation shell — from grand canonical Monte
Carlo simulations by coupling to an ensemble of thermodynamic baths rather
than a single bath. To achieve this, we introduced a weight distribution
function embedded in the definition of the marginal probability distribution.
We then considered a pairwise interaction between the particles using
Lennard–Jones’ internal potential function and derived an analytical
functional form for the weight distribution. In an initial proof-of-principle
numerical experiment, we showed that our model is capable to correctly predict
the particle number distribution of a non-interacting system of hard-sphere
particles. Then, we challenged our model with an interacting Lennard–Jones
fluid for different interaction strengths $\epsilon^{\prime}$ and reduced
densities $\rho^{\prime}$. The simulation results demonstrate that our
proposed model allows us to simulate the particle number distribution of an
interacting system in the framework of computationally economic grand
canonical simulations.

We note that we discarded higher-order terms in Eq. (2.2) in order to arrive
at a suitable analytic form for the canonical partition function. Potentially,
including more terms in Eq. (2.2) will result in a mathematically well-behaved
model and make the optimization process faster. Such an analysis is left for
future work if the current model proves to be insufficient or not flexible
enough for more complicated interactions.

For the numerical integration of the newly derived expression for the particle
number distribution in the solvation shell, Eq. (71), we had to transform the
indefinite integral over the chemical potential into a definite one.
Currently, we determine both the limits of the integral and the number of
integration points based on a manual convergence check. In an ongoing project,
we attempt to automate this process and reduce the computational effort of
this step.

An obvious drawback of our current study is certainly the need to optimize the
three (or two) parameters of our model, which requires reference data that we
generate from a canonical Monte Carlo simulation for the entire system which
is usually rather computationally expensive. In upcoming research, we will
investigate the transferability of these optimized parameters, i.e. if a
parameter set optimized for a given reference system can be used for a
slightly modified system (e.g. conformational changes) if critical properties
like the size of the solvation shell remain constant. If so, a single
reference calculation would allow us to perform the calculation of
thermodynamical properties or even reactivity studies in this new grand
canonical framework. In addition, we will explore if a parameter optimization
based on the principle of maximum entropy is a valid alternative that
eradicates the need for reference data.

Last but not least, it is remarkable that Eq. (71) provides a discrete-state
partition function to represent the distribution of the number of particles
inside the core. Quasi-chemical theories are based on this feature 31,
although a discrete-state partition function is the primary obstacle in the
construction of a quasi-chemical model. This conceptual resemblance inspired
us to think about how we can link our proposed model with quasi-chemical
theories. Since we incorporate the excess chemical potential inside the
definition of $\gamma$ whereas in quasi-chemical theories, the chemical
potential of the solvated species is usually expressed in terms of ideal and
excess contributions 6, we consider this to be a good starting point for
linking the two models.

Most importantly, in future work we will explore how appropriate our grand
canonical model describes more complicated interactions provided by advanced
force fields or the full quantum-chemical Hamiltonian. In an ongoing project,
we therefore calculate thermodynamic properties of solvated systems and
analyze the computational advantage of our multi-bath grand canonical approach
compared to a standard canonical simulations or other alternatives. In
addition to solvation studies and the development of novel solvation models,
we assume that our model will prove useful for the simulation of properties in
metal-organic frameworks and other nano-confined structures, where the rather
small system of interest (e.g. an active site) is in indirect contact with its
environment.

##  6 Appendix

###  6.1 Chain rule for random variables

For a set of random variables with indices $X_{1},\ldots,X_{n}$, we can use
the conditional probability definition to determine the value of each
component of the joint distribution

$\displaystyle\mathrm{p}\left(X_{n},\ldots,X_{1}\right)=\mathrm{p}\left(X_{n}\mid
X_{n-1},\ldots,X_{1}\right)\cdot\mathrm{p}\left(X_{n-1},\ldots,X_{1}\right)\,.$
(84)

The chain rule for random variables is obtained by applying the aforementioned
rule for each last term

$\displaystyle\mathrm{p}\left(\bigcap_{k=1}^{n}X_{k}\right)=\prod_{k=1}^{n}\mathrm{p}\left(X_{k}\mid\bigcap_{j=1}^{k-1}X_{j}\right)\,.$
(85)

With three variables the chain rule produces the following product of
conditional probabilities

$\displaystyle\mathrm{p}\left(X_{3},X_{2},X_{1}\right)=\mathrm{p}\left(X_{3}\mid
X_{2},X_{1}\right)\cdot\mathrm{p}\left(X_{2}\mid
X_{1}\right)\cdot\mathrm{p}\left(X_{1}\right)\,.$ (86)

###  6.2 Canonical partition function of the system of interacting particles

The canonical partition function for a system of interacting particles with a
Lennard–Jones internal interaction potential can be written as

$\displaystyle\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})=\frac{1}{N^{{}^{\prime}}!}\int\prod_{i=1}^{N^{{}^{\prime}}}\left(\frac{\mathrm{d}^{3}\vec{p}_{i}\mathrm{~{}d}^{3}\vec{q}_{i}}{h^{3}}\right)\exp\left[-\beta^{j}_{\mathrm{c}}\sum_{i}\frac{\vec{{p}_{i}}^{2}}{2m}\right]\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$
$\displaystyle=\frac{1}{N^{{}^{\prime}}!}\frac{1}{h^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{d}^{3}\vec{p}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\sum_{i}\frac{\vec{{p}_{i}}^{2}}{2m}\right]\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$
$\displaystyle=\underbrace{\frac{1}{N^{{}^{\prime}}!}\frac{1}{h^{3N^{{}^{\prime}}}}\left[\int{\mathrm{d}{p}_{i}}\exp\frac{{{p}_{i}}^{2}}{2mk_{B}T^{j}_{\mathrm{c}}}\right]^{3N^{{}^{\prime}}}}_{\frac{1}{N!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$
$\displaystyle={\frac{1}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$

We can now multiply Eq. (6.2) with the sum of the probability density function
evaluated over the domain of the entire phase space of a non-interacting
system of particles, which yields 1:

$\displaystyle\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})={\frac{1}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]\times
1$
$\displaystyle={\frac{1}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$
$\displaystyle\times\int\prod_{i=1}^{N^{{}^{\prime}}}\left(\frac{\mathrm{d}^{3}\vec{p}_{i}\mathrm{~{}d}^{3}\vec{q}_{i}}{h^{3}}\right)\frac{1}{\mathcal{Z_{\textbf{Non-
Int}}}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}\exp\left[-\beta^{j}_{\mathrm{c}}\sum_{i=1}^{N^{{}^{\prime}}}\frac{p_{i}^{2}}{2m}\right]$
(88)

In Eq. (6.2), the integrand can be separated into two distinct parts, one over
the momentum and one over the position. The position integrand, which
corrsponds to a non-interacting system, has an exponential term that expresses
the total kinetic energy and can be merged with the integrand of the total
internal potential function of the interacting system. The result is a double-
integral that can be considered as the expectation value of the internal
potential function multiplied by $\beta$. It is also worth mentioning that
this expectation value has been calculated based on the probability density
function of a non-interacting system of particles.

$\displaystyle\mathcal{Z}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})={\frac{1}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$
$\displaystyle\times\int\prod_{i=1}^{N^{{}^{\prime}}}\left(\frac{\mathrm{d}^{3}\vec{p}_{i}}{h^{3}}\right)\frac{1}{\mathcal{Z_{\textbf{Non-
Int}}}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}\exp\left[-\beta^{j}_{\mathrm{c}}\sum_{i=1}^{N^{{}^{\prime}}}\frac{p_{i}^{2}}{2m}\right]\int\prod_{i=1}^{N^{{}^{\prime}}}\mathrm{~{}d}^{3}\vec{q}_{i}$
$\displaystyle={\frac{1}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}\mathrm{~{}d}^{3}\vec{q}_{i}$
$\displaystyle\times\int\prod_{i=1}^{N^{{}^{\prime}}}{\mathrm{~{}d}^{3}\vec{q}_{i}}\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]$
$\displaystyle\times\int\prod_{i=1}^{N^{{}^{\prime}}}\left(\frac{\mathrm{d}^{3}\vec{p}_{i}}{h^{3}}\right)\frac{1}{\mathcal{Z_{\textbf{Non-
Int}}}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}\exp\left[-\beta^{j}_{\mathrm{c}}\sum_{i=1}^{N^{{}^{\prime}}}\frac{p_{i}^{2}}{2m}\right]$
$\displaystyle={\frac{1}{N^{{}^{\prime}}!}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3N^{{}^{\prime}}}}\int\prod_{i=1}^{N^{{}^{\prime}}}\mathrm{~{}d}^{3}\vec{q}_{i}$
$\displaystyle\times\underbrace{\int\prod_{i=1}^{N^{{}^{\prime}}}\left(\frac{{\mathrm{~{}d}^{3}\vec{q}_{i}}\mathrm{d}^{3}\vec{p}_{i}}{h^{3}}\right)\frac{1}{\mathcal{Z_{\textbf{Non-
Int}}}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})}\exp\left[-\beta^{j}_{\mathrm{c}}\sum_{i=1}^{N^{{}^{\prime}}}\frac{p_{i}^{2}}{2m}\right]\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]}_{\left\langle\exp\left[-\beta^{j}_{\mathrm{c}}\phi\left(\vec{q}_{1},\cdots,\vec{q}_{N^{{}^{\prime}}}\right)\right]\right\rangle_{\textbf{Non-
Int}}}$ (89)

###  6.3 Attempt to write the weight distribution function without
approximation (pseudo-Sophomore’s dream)

Let us assume that there is a possibility to derive an analytical form for
$p\left(\mu^{j}_{\mathrm{c}},\beta^{j},V_{\mathrm{c}}\right)$ without any
approximation. By doing so, we can write

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)={\mathcal{Q}^{-1}\left(g_{1}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right),g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)+\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}\right)}$
$\displaystyle=\left[\sum_{N^{{}^{\prime}}=0}^{\infty}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\right.$
$\displaystyle\left.\exp\left[\left(\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\right)\,N^{{}^{\prime}}-\left(\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\right)\right]\right]^{-1}\,,$
(90)

where

$\displaystyle\ln{g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)}+g_{2}\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)=\ln\left[{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\right]+\left[{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\right]={\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}+\ln{\gamma\,k}+{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\,.$
(91)

With this we can write

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)$
$\displaystyle\propto\left[\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\exp\left[({\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}+\ln{\gamma\,k}+{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}})\,N^{{}^{\prime}}-({\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}+\ln{\gamma\,k}+{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}})\right]}\right]^{-1}$
$\displaystyle\propto\left[\sum_{N^{{}^{\prime}}=0}^{\infty}{\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}N^{{}^{\prime}}}\mathcal{C}_{1}(\beta^{j}_{\mathrm{c}},N^{{}^{\prime}},V_{\mathrm{c}})\left[\,\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{N^{{}^{\prime}}}\left[\gamma\,k\right]^{N^{{}^{\prime}}}\left[\mathrm{e}^{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\right]^{N^{{}^{\prime}}}\left[\,\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{-1}\left[\gamma\,k\right]^{-1}\left[\mathrm{e}^{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\right]^{-1}}\right]^{-1}$
$\displaystyle\propto\left[\left[\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\right]^{-1}\left[\gamma\,k\right]^{-1}\left[\mathrm{e}^{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\right]^{-1}\sum_{N^{{}^{\prime}}=0}^{\infty}\frac{1}{N^{{}^{\prime}}!}\left[\underbrace{{{V_{\mathrm{c}}}\left[\frac{\sqrt{2mk_{B}T^{j}_{\mathrm{c}}}}{h}\right]^{3}}\,\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\,\gamma\,k\,\mathrm{e}^{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}}_{\alpha}\right]^{N^{{}^{\prime}}}\right]^{-1}$
(92)

Interpreting the infinite sum over $N^{{}^{\prime}}$as the power series of the
exponential function,

$\displaystyle\exp{\alpha}:=\sum_{N^{{}^{\prime}}=0}^{\infty}\frac{\alpha^{N^{{}^{\prime}}}}{N^{{}^{\prime}}!}\,,$
(93)

we arrive at

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta^{j}_{\mathrm{c}},V_{\mathrm{c}}\right)\propto\,\gamma\,k\,\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{{-\frac{\gamma\,k\,V_{\mathrm{c}}}{{\Lambda^{j}_{\mathrm{c}}}^{3}}}\mathrm{e}^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\gamma\,k\,{e^{\beta^{j}_{\mathrm{c}}\mu^{j}_{\mathrm{c}}}}}}$
(94)

we consider all artificial baths to have a fixed temperature equal to the
thermal bath surrounding the composite system, i.e.
$\beta^{j}_{\mathrm{c}}=\beta$, for all $j$ and only the artificial chemical
potentials are allowed to fluctuate.

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{\gamma\,k\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\gamma\,k\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{{-\frac{\gamma\,k\,V_{\mathrm{c}}}{{\Lambda^{j}_{\mathrm{c}}}^{3}}}\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\gamma\,k\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}}}{\sum_{j}\,\gamma\,k\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\gamma\,k\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{{-\frac{\gamma\,k\,V_{\mathrm{c}}}{{\Lambda^{j}_{\mathrm{c}}}^{3}}}\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\gamma\,k\,{e^{\beta^{\mu^{j}_{\mathrm{c}}}}}}}}$
(95)

a simplified notation with

$\displaystyle\quad\alpha_{1}=\,\gamma\,k,\quad\alpha_{2}=\frac{\alpha_{1}\,V_{\mathrm{c}}}{{\Lambda^{j}_{\mathrm{c}}}^{3}},\quad\gamma>0$
(96)

leads to

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{\alpha_{1}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}}}{\sum_{j}\,\alpha_{1}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}}}\,.$
(97)

The normalization factor of $p(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}})$ can
be written as

$\displaystyle{\sum_{j}\,\alpha_{1}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}}}=\beta\int_{-\infty}^{\infty}\text{d}\mu_{\mathrm{c}}\,\alpha_{1}\,\mathrm{e}^{\beta\mu_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}}\,.$
(98)

We now impose a change of variables
$\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}=u$, such that
$\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}\text{d}\mu_{\mathrm{c}}={\text{d}\,u}/{\beta\,u}$
and ${e^{\beta\mu_{\mathrm{c}}}}={\ln{u}}/{\alpha_{1}}$. We also note that for
a hard sphere system of particles we can easily prove that $k<0$ and on the
other hand we know $\gamma>0$ then $\alpha_{1}<0$.

$\displaystyle\beta\int_{-\infty}^{\infty}\text{d}\mu_{\mathrm{c}}\,\alpha_{1}\,\mathrm{e}^{\beta\mu_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}}=-\beta\int_{0}^{1}\,\frac{\text{d}u}{\beta}\,u^{-1}\,u^{1}\,\mathrm{e}^{-\frac{\alpha_{2}}{\alpha_{1}}\,u\,\ln{u}}$
$\displaystyle=-\frac{\beta}{\beta}\int_{0}^{1}\,{\text{d}u}\,\mathrm{e}^{-\frac{\alpha_{2}}{\alpha_{1}}\,u\,\ln{u}}=-\int_{0}^{1}\,{\text{d}u}\,\sum_{n=0}^{\infty}({-\frac{\alpha_{2}}{\alpha_{1}}})^{n}\frac{u^{n}(\ln{u})^{n}}{n!}$
$\displaystyle=-\sum_{n=0}^{\infty}\,({-\frac{\alpha_{2}}{\alpha_{1}}})^{n}\,\int_{0}^{1}{\text{d}u}\frac{u^{n}(\ln{u})^{n}}{n!}$
(99)

The above integral was historically computed by Johann Bernoulli 34,

$\displaystyle\int_{0}^{1}{\text{d}u}\frac{u^{n}(\ln{u})^{n}}{n!}=(-1)^{n}(n+1)^{-(n+1)}\,.$
(100)

With this we arrive at

$\displaystyle{\beta}\int_{-\infty}^{\infty}\text{d}\mu_{\mathrm{c}}\,\alpha_{1}\,\mathrm{e}^{\beta\mu_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu_{\mathrm{c}}}}}}=-\sum_{n=0}^{\infty}\left(\frac{V_{\mathrm{c}}}{{\Lambda_{\mathrm{c}}}^{3}}\right)^{n}(n+1)^{-(n+1)}\,.$
(101)

The infinite sum can only be computed when the inter-particle distance
approaches the radius of the solvation shell, i.e.
${V_{\mathrm{c}}}/{{\Lambda_{\mathrm{c}}}^{3}}=1$. Only under this condition
we can determine the normalization as

$\displaystyle{\sum_{j}\,\alpha_{1}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{2}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}}}=-\sum_{n=0}^{\infty}(n+1)^{-(n+1)}\approx\,-1.29\,,$
(102)

and write the weighted distribution function as

$\displaystyle
p\left(\mu^{j}_{\mathrm{c}},\beta,V_{\mathrm{c}}\right)=\frac{\alpha_{1}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\,\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}\mathrm{e}^{-\alpha_{1}\,\mathrm{e}^{\beta\mu^{j}_{\mathrm{c}}}\mathrm{e}^{\alpha_{1}\,{e^{\beta\mu^{j}_{\mathrm{c}}}}}}}{-1.29}\quad\text{s.t.}\quad\alpha_{1}=\,\gamma\,k<0,\quad\gamma>0\,.$
(103)

Since this condition for which we can solve the expression is unphysical, we
decided to approximate $g_{2}+\ln g_{2}$ by introducing the two parameters
$\lambda_{1}$ and $\lambda_{2}$ as shown in Eq. 45 and arrive at a model for
which we can formulate an analytic form without imposing unphysical
conditions.

###  6.4 Moments in terms of the cumulants

28 The expectation value of $\mathrm{e}^{-\mathrm{i}kx}$ can be written as,

$\displaystyle\left\langle\mathrm{e}^{-\mathrm{i}kx}\right\rangle=\int\mathrm{d}xf(x)\mathrm{e}^{-\mathrm{i}kx}=\mathcal{F}_{x}[f(x)](k)=F(k)\,,$
(104)

in which $f(x)$ is the density probability function. We can interpret $F(k)$
as the Fourier transform of the probability density function $f(x)$, called
the characteristic function. If we expand $\mathrm{e}^{-\mathrm{i}kx}$ based
its power series,

$\displaystyle\left\langle\mathrm{e}^{-\mathrm{i}kx}\right\rangle=\left\langle\sum_{m=0}^{\infty}\frac{(-\mathrm{i}k)^{m}}{m!}x^{m}\right\rangle=\sum_{m=0}^{\infty}\frac{(-\mathrm{i}k)^{m}}{m!}\left\langle
x^{m}\right\rangle\,.$ (105)

The logarithm of the characteristic function, $F(k)$, is the cumulant
generating function. Expanding this expression gives the cumulants of the
distribution defined by

$\displaystyle\ln
F(k)=\sum_{m=1}^{\infty}\frac{(-\mathrm{i}k)^{m}}{m!}\left\langle
x^{m}\right\rangle_{c}$ (106)

and hence

$\displaystyle
F(k)=\exp\left[\sum_{m=1}^{\infty}\frac{(-\mathrm{i}k)^{m}}{m!}\left\langle
x^{m}\right\rangle_{c}\right]\,.$ (107)

Combining the Eqns. (104), (105) and (107), we obtain

$\displaystyle\exp\left[\sum_{m=1}^{\infty}\frac{(-\mathrm{i}k)^{m}}{m!}\left\langle
x^{m}\right\rangle_{c}\right]=\sum_{m=0}^{\infty}\frac{(-\mathrm{i}k)^{m}}{m!}\left\langle
x^{m}\right\rangle\,.$ (108)

The authors thank the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany´s Excellence Strategy – EXC 2033 – 390677874 –
RESOLV for funding. C. J. S. further acknowledges funding by the Ministry of
Innovation, Science and Research of North Rhine-Westphalia (“NRW
Rückkehrerprogramm”). M. R. is grateful to Dipesh Somvanshi and Siyavash
Moradi for help with the implementation of the algorithms. The authors also
thank Prof. Martin Head-Gordon for introducing them to the work of Prof. Dilip
Asthagiri and hence initiating this study.

## References

  * Sethna 2021 Sethna, J. _Statistical mechanics: entropy, order parameters, and complexity_ ; Oxford University Press, USA, 2021; Vol. 14
  * Mandelbrot 1989 Mandelbrot, B. B. Temperature fluctuation: a well-defined and unavoidable notion. _Phys. Today_ 1989, _42_ , 71
  * Dunkel and Hilbert 2014 Dunkel, J.; Hilbert, S. Consistent thermostatistics forbids negative absolute temperatures. _Nat. Phys._ 2014, _10_ , 67–72
  * Mandelbrot 1962 Mandelbrot, B. The role of sufficiency and of estimation in thermodynamics. _Ann. Math. Stat._ 1962, 1021–1038
  * Dixit 2015 Dixit, P. D. Detecting temperature fluctuations at equilibrium. _PCCP_ 2015, _17_ , 13000–13005
  * Pratt and Rempe 1999 Pratt, L. R.; Rempe, S. B. Quasi-chemical theory and implicit solvent models for simulations. AIP Conf. Proc. 1999; pp 172–201
  * Juffer and Berendsen 1993 Juffer, A.; Berendsen, H. Dynamic surface boundary conditions: A simple boundary model for molecular dynamics simulations. _Mol. Phys._ 1993, _79_ , 623–644
  * Remsing et al. 2016 Remsing, R. C.; Liu, S.; Weeks, J. D. Long-ranged contributions to solvation free energies from theory and short-ranged models. _Proc. Natl. Acad. Sci. U.S.A._ 2016, _113_ , 2819–2826
  * Hegemann et al. 2017 Hegemann, D.; Hocquard, N.; Heuberger, M. Nanoconfined water can orient and cause long-range dipolar interactions with biomolecules. _Sci. Rep._ 2017, _7_ , 2045–2322
  * Belch and Berkowitz 1985 Belch, A. C.; Berkowitz, M. Molecular dynamics simulations of TIPS2 water restricted by a spherical hydrophobic boundary. _Chem. Phys. Lett._ 1985, _113_ , 278–282
  * Klamt et al. 1998 Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. _J. Phys. Chem. A_ 1998, _102_ , 5074–5085
  * Sato et al. 2000 Sato, H.; Kovalenko, A.; Hirata, F. Self-consistent field, ab initio molecular orbital and three-dimensional reference interaction site model study for solvation effect on carbon monoxide in aqueous solution. _J. Chem. Phys._ 2000, _112_ , 9463–9468
  * Mennucci 2010 Mennucci, B. Continuum Solvation Models: What Else Can We Learn from Them? _J, Phys. Chem. Lett._ 2010, _1_ , 1666–1674
  * Stein et al. 2019 Stein, C. J.; Herbert, J. M.; Head-Gordon, M. The Poisson–Boltzmann model for implicit solvation of electrolyte solutions: Quantum chemical implementation and assessment via Sechenov coefficients. _J. Chem. Phys._ 2019, _151_ , 224111
  * Herbert 2021 Herbert, J. M. Dielectric continuum methods for quantum chemistry. _WIREs Comput. Mol. Sci._ 2021, _11_ , e1519
  * Ringe et al. 2022 Ringe, S.; Hörmann, N. G.; Oberhofer, H.; Reuter, K. Implicit Solvation Methods for Catalysis at Electrified Interfaces. _Chem. Rev._ 2022, _122_ , 10777–10820
  * Thomas et al. 2013 Thomas, D. G.; Chun, J.; Chen, Z.; Wei, G.; Baker, N. A. Parameterization of a geometric flow implicit solvation model. _J. Comput. Chem._ 2013, _34_ , 687–695
  * Lange et al. 2020 Lange, A. W.; Herbert, J. M.; Albrecht, B. J.; You, Z.-Q. Intrinsically smooth discretisation of Connolly’s solvent-excluded molecular surface. _Mol. Phys._ 2020, _118_ , e1644384
  * Rana and Nguyen 2022 Rana, M. M.; Nguyen, D. D. EISA-Score: Element Interactive Surface Area Score for Protein–Ligand Binding Affinity Prediction. _J. Chem. Inf. Model.J_ 2022, _62_ , 4329–4341
  * Li and Hartke 2013 Li, Y.; Hartke, B. Assessing Solvation Effects on Chemical Reactions with Globally Optimized Solvent Clusters. _Chem. Phys. Chem._ 2013, _14_ , 2678–2686
  * Kildgaard et al. 2018 Kildgaard, J. V.; Mikkelsen, K. V.; Bilde, M.; Elm, J. Hydration of Atmospheric Molecular Clusters: A New Method for Systematic Configurational Sampling. _J. Phys. Chem. A_ 2018, _122_ , 5026–5036
  * Simm et al. 2020 Simm, G. N.; Türtscher, P. L.; Reiher, M. Systematic microsolvation approach with a cluster-continuum scheme and conformational sampling. _J. Comput. Chem._ 2020, _41_ , 1144–1155
  * Steiner et al. 2021 Steiner, M.; Holzknecht, T.; Schauperl, M.; Podewitz, M. Quantum Chemical Microsolvation by Automated Water Placement. _Molecules_ 2021, _26_
  * Dixit et al. 2017 Dixit, P. D.; Bansal, A.; Chapman, W. G.; Asthagiri, D. Mini-grand canonical ensemble: chemical potential in the solvation shell. _J. Chem. Phys._ 2017, _147_ , 164901
  * Dixit 2013 Dixit, P. D. A maximum entropy thermodynamics of small systems. _J. Chem. Phys._ 2013, _138_ , 184111
  * Bansal et al. 2017 Bansal, A.; Valiya Parambathu, A.; Asthagiri, D.; Cox, K. R.; Chapman, W. G. Thermodynamics of mixtures of patchy and spherical colloids of different sizes: A multi-body association theory with complete reference fluid information. _J. Chem. Phys._ 2017, _146_ , 164904
  * Bansal et al. 2016 Bansal, A.; Asthagiri, D.; Cox, K. R.; Chapman, W. G. Structure and thermodynamics of a mixture of patchy and spherical colloids: A multi-body association theory with complete reference fluid information. _J. Chem. Phys._ 2016, _145_ , 074904
  * Kardar 2007 Kardar, M. _Statistical physics of particles_ ; Cambridge University Press, 2007
  * Tetrode 1912 Tetrode, H. Die chemische Konstante der Gase und das elementare Wirkungsquantum. _Annalen der Physik_ 1912, _343_ , 434–442
  * Sackur 1911 Sackur, O. Die Anwendung der kinetischen Theorie der Gase auf chemische Probleme. _Annalen der Physik_ 1911, _341_ , 958–980
  * Beck et al. 2006 Beck, T. L.; Paulaitis, M. E.; Pratt, L. R. _The potential distribution theorem and models of molecular solutions_ ; Cambridge University Press, 2006
  * Widom 1963 Widom, B. Some topics in the theory of fluids. _J. Chem. Phys._ 1963, _39_ , 2808–2812
  * Rowlinson and Widom 1982 Rowlinson, J.; Widom, B. Molecular Theory of Capillarity (Oxford University Press). 1982
  * Dunham 2005 Dunham, W. Touring the calculus gallery. _Am. Math. Mon._ 2005, _112_ , 1–19