Liouville-type equations for the n-particle distribution functions of an open system

TL;DR

This paper develops a mathematical framework for open systems exchanging particles and momentum with a reservoir, deriving Liouville-type equations for n-particle distributions that incorporate reservoir effects and particle exchanges.

Contribution

It introduces a set of coupled Liouville-type equations for open systems' n-particle distributions derived from Hamiltonian dynamics, extending previous models.

Findings

Derived hierarchy of equations for open system dynamics.

Accounted for particle exchange and external momentum forcing.

Compared with Bergmann-Lebowitz model for open systems.

Abstract

In this work we derive a mathematical model for an open system that exchanges particles and momentum with a reservoir from their joint Hamiltonian dynamics. The complexity of this many-particle problem is addressed by introducing a countable set of n-particle phase space distribution functions just for the open subsystem, while accounting for the reservoir only in terms of statistical expectations. From the Liouville equation for the full system we derive a set of coupled Liouville-type equations for the n-particle distributions by marginalization with respect to reservoir states. The resulting equation hierarchy describes the external momentum forcing of the open system by the reservoir across its boundaries, and it covers the effects of particle exchanges, which induce probability transfers between the n- and (n+1)-particle distributions. Similarities and differences with the Bergmann-Lebowitz model of open systems (P.G.Bergmann, J.L. Lebowitz, Phys.Rev., 99:578--587 (1955)) are discussed in the context of the implementation of these guiding principles in a computational scheme for molecular simulations.

Equations

$$\frac{\partial f(t,\mathbf{q},\mathbf{p})}{\partial t}+\sum _{i=1}^N\left(v_i\cdot {\nabla }_{q_i}+F_i\cdot {\nabla }_{p_i}\right)f(t,\mathbf{q},\mathbf{p})=0.$$

$$H_n=\sum _{i=1}^n\frac{{\vec{p}}_i^2}{2M}+\sum _{i=1}^n\sum _{j\ne i}^n\frac{1}{2}V({\vec{q}}_j-{\vec{q}}_i)({\vec{q}}_i,{\vec{q}}_j\in \Omega )$$

$$\begin{array}[]{rcrcl}f_{n}&:&\mathbb{R}^{+}\times S^{n}&\to&\mathbb{R}\\ &&(t,\boldsymbol{X}^{n})&\mapsto&f_{n}(t,\boldsymbol{X}^{n})\qquad\text{for}\qquad(n=0,...,N)\,,\end{array}$$

$${\mathbf{X}}^n=[X_1,...,X_n]\text{and}{X}_i=({\vec{q}}_i,{\vec{p}}_i)\in S=\Omega \times {\mathbb{R}}^3,$$

$$f_n(t,{\mathbf{X}}^n)=\left(\genfrac{}{}{0ex}{}{N}{n}\right)\underset{(S^c{)}^{N-n}}{\int }{F}_N(t,{\mathbf{X}}^n,{\mathbf{\Xi }}_n^N)d{\mathbf{\Xi }}_n^N\text{where}{S}^c={\Omega }_c\times {\mathbb{R}}^3$$

$$\boldsymbol{\Xi}_{n}^{N}\equiv[\Xi_{n+1},.....{{\color[rgb]{0,0,0}\Xi_{N}}}]\qquad\text{where}\qquad\Xi_{i}=(\vec{q}_{i},\vec{p}_{i})\in S^{c}\,.$$

$$f_n^{\ast }(t,{\mathbf{X}}^n)={\left(\genfrac{}{}{0ex}{}{N}{n}\right)}^{-1}{f}_n(t,{\mathbf{X}}^n)$$

$$\sum _{n=0}^N\underset{{\Omega }^n}{\int }\underset{({\mathbb{R}}^3{)}^n}{\int }{f}_n(t,(\mathbf{q},\mathbf{p}))d\mathbf{p}d\mathbf{q}=1.$$

$$\begin{array}[]{rrcl}F_{N}:&\displaystyle\mathbb{R}^{+}\times(U\times\mathbb{R}^{3})^{N}&\to&\displaystyle\mathbb{R}\\ &\displaystyle(t,\boldsymbol{X}^{N})&\mapsto&\displaystyle F_{N}(t,\boldsymbol{X}^{N})\end{array}\,,$$

$$\frac{d{F}_N}{dt}=0$$

$$\frac{\partial F_{N}}{\partial t}={\color[rgb]{0,0,0}-}\sum\limits_{i=1}^{N}\left[\nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}F_{N}\right)+\nabla_{\vec{p}_{i}}\cdot\left(-\nabla_{q_{i}}V_{\text{tot}}(\boldsymbol{q}^{N})F_{N}\right)\right]\,,$$

$$V_{\text{tot}}({\mathbf{q}}^N)=\sum _{i<j}^{N}V({\vec{q}}_i-{\vec{q}}_j)$$

$$\underset{B_r(0)}{\int }{\nabla }_{{\vec{p}}_k}\cdot \left(-{\nabla }_{q_k}{V}_{\text{tot}}({\mathbf{q}}^N)F_N\right)d^3{p}_k=-\underset{\partial B_r(0)}{\int }\vec{n}\cdot \left({\nabla }_{q_k}{V}_{\text{tot}}({\mathbf{q}}^N)\right)F_{N}d{\sigma }_{p_k}\to 0(r\to \infty )$$

$$\int\limits_{(S^{c})^{N-n}}\nabla_{\vec{p}_{i}}\cdot\Bigl{(}\nabla_{q_{i}}V(\vec{q}_{i}-\vec{q}_{j})\,F_{N}(t,\boldsymbol{X}^{n},\boldsymbol{\Xi}_{n}^{N})\Bigr{)}d\Xi_{n}^{N}=\nabla_{\vec{p}_{i}}\cdot\Bigl{(}\nabla_{q_{i}}V(\vec{q}_{i}-\vec{q}_{j})\,f^{*}_{n}(t,\boldsymbol{X}^{n})\Bigr{)}\,,$$

$$\begin{array}[]{rl}\displaystyle\int\limits_{S^{c}}\int\limits_{(S^{c})^{N-n-1}}&\displaystyle\nabla_{\vec{p}_{i}}\cdot\Bigl{(}-\nabla_{q_{i}}V(\vec{q}_{i}-\vec{q}_{j})\,F_{N}(t,\boldsymbol{X}^{n-1,},X_{i},(\vec{q}_{j},\vec{p}_{j}),\boldsymbol{\Xi}_{n+1}^{N})\Bigr{)}{{\color[rgb]{0,0,0}d\boldsymbol{\Xi}_{n+1}^{N}}}\ dp_{j}dq_{j}\,.\end{array}$$

$$\begin{array}[]{rl}\displaystyle\int\limits_{S^{c}}\int\limits_{(S^{c})^{N-n-1}}&\displaystyle\nabla_{\vec{p}_{i}}\cdot\Bigl{(}-\nabla_{q_{i}}V(\vec{q}_{i}-\vec{q}_{j})\,F_{N}(t,\boldsymbol{X}^{n-1},X_{i},(\vec{q}_{j},\vec{p}_{j}),\boldsymbol{\Xi}_{n+1}^{N})\Bigr{)}d\boldsymbol{\Xi}_{n+1}^{N}\ dp_{j}dq_{j}\,.\\ &\displaystyle=\nabla_{\vec{p}_{i}}\cdot\Bigl{(}\vec{F}_{\text{av}}(\vec{q}_{i})f^{*}_{n}(t,\boldsymbol{X}^{n-1},X_{i})\Bigr{)}\,,\end{array}$$

$${\vec{F}}_{\text{av}}({\vec{q}}_i)=-\underset{S^c}{\int }{\nabla }_{{\vec{q}}_i}V({\vec{q}}_i-{\vec{q}}_j)f_2^{\circ }(X_j|X_i)d{X}_j$$

$$\nabla_{\vec{p}_{i}}\cdot\Bigl{(}\vec{F}_{\text{av}}(\vec{q}_{i})f_{n}(t,\boldsymbol{X}^{n-1},X_{i})\Bigr{)}$$

$${{N}\choose{n}}\int\limits_{(S^{c})^{N-n}}\nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}F_{N}(t,\boldsymbol{X}^{n},\boldsymbol{\Xi}_{n}^{N}\right)\,d\boldsymbol{\Xi}_{n}^{N}=\nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}f_{n}\right)$$

$$\begin{array}[]{rl}&\displaystyle{{N}\choose{n}}\sum\limits_{i={{\color[rgb]{0,0,0}n+1}}}^{N}\int\limits_{S^{c}}\displaystyle\int\limits_{(S^{c})^{N-n-1}}\nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}F_{N}(t,\boldsymbol{X}^{n},\boldsymbol{\Xi}_{n+1}^{N},(\vec{q}_{i},\vec{p}_{i})\right)\,d\boldsymbol{\Xi}_{n+1}^{N}\ d\Xi_{i}\\[20.0pt] &={{N}\choose{n}}(N-n)\ \displaystyle\int\limits_{S^{c}}\nabla_{\vec{q}_{i}}\cdot\biggl{\{}\int\limits_{(S^{c})^{N-n-1}}\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}F_{N}(t,\boldsymbol{X}^{n},\boldsymbol{\Xi}_{n+1}^{N},(\vec{q}_{i},\vec{p}_{i})\right)\,d\boldsymbol{\Xi}_{n+1}^{N}\biggr{\}}d\Xi_{i}\\ &=-\displaystyle(n+1)\int\limits_{\partial\Omega}\int\limits_{\mathbb{R}^{3}}\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}\right)\ {\widehat{f}}_{n+1}(t,\boldsymbol{X}^{n},(\vec{q}_{i},\vec{p}_{i})\,d^{3}p_{i}\,d\sigma_{i}\,.\end{array}$$

$${\widehat{f}}_{n+1}=\left\{\begin{array}[]{l@{\quad}l}f_{n+1}&\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}>0\right)\\ f_{n}f_{1}^{\circ}&\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}<0\right)\end{array}\right.\,.$$

$${\widehat{f}}_{n+1}=\left\{\begin{array}[]{l@{\quad}l}f_{n+1}&\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}>0\right)\\ f_{n+1}^{\text{GC}}&\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}<0\right)\end{array}\right.\,.$$

$$\frac{\partial f_{n}}{\partial t}+\sum\limits_{i=1}^{n}\left(\nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}f_{n}\right)+\nabla_{\vec{p}_{i}}\cdot\left(\vec{F}_{i}f_{n}\right)\rule{0.0pt}{12.0pt}\right)=\Psi_{n}+\Phi_{n}^{n+1}\,,$$

$$\begin{array}[]{rcl}\displaystyle\vec{F}_{i}&=&\displaystyle-\sum\limits_{j=1;j\not=i}^{n}\nabla_{\vec{q}_{i}}V(\vec{q}_{i}-\vec{q}_{j})\end{array}$$

$$\begin{array}[]{rcl}\Psi_{n}&=&\displaystyle-\sum\limits_{i=1}^{n}\nabla_{\vec{p}_{i}}\cdot\Bigl{(}\vec{F}_{\text{av}}(\vec{q}_{i})f_{n}(t,\boldsymbol{X}^{i-1},X_{i},\boldsymbol{X}_{i}^{n-i})\Bigr{)}\\ \Phi_{n}^{n+1}&=&(n+1)\displaystyle\int\limits_{\partial\Omega}\int\limits_{({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n})>0}\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}\right)\ \left({f}_{n+1}\left(t,\boldsymbol{X}^{n},(\vec{q}_{i},\vec{p}_{i})\right)-{f}_{n}\left(t,\boldsymbol{X}^{n}\right){f}_{1}^{\circ}\left(\vec{q}_{i},-\vec{p}_{i}\right)\rule{0.0pt}{12.0pt}\right)\,d^{3}p_{i}\ d\sigma_{i}\end{array}$$

$${\vec{F}}_{\text{av}}({\vec{q}}_i)=-\underset{S^c}{\int }{\nabla }_{{\vec{q}}_i}V({\vec{q}}_i-{\vec{q}}_j)f_2^{\circ }(X_j|X_i)d{X}_j.$$

$$Q_N={\int }_U\frac{1}{h^{3N}N!}{e}^{-\beta H_N}{d}^N\vec{q}{d}^N\vec{p}$$

$$F_N({\mathbf{X}}^N)=\frac{\frac{1}{h^{3N}N!}{e}^{-\beta H_N}}{Q_N}$$

$$f_n(t,{\mathbf{X}}^n)=\frac{N!}{n!(N-n)!}\frac{\frac{1}{h^{3N}N!}{e}^{-\beta H_n}\underset{(S^c{)}^{N-n}}{\int }{e}^{-\beta H_{N-n}+V_{cx}({\vec{q}}_{\Omega }-{\vec{q}}_{{\Omega }^c})}{d}^{N-n}\Xi }{Q_N}$$

$$f_n(t,{\mathbf{X}}^n)=\frac{1}{Q_N}\left[\frac{1}{n!(N-n)!}\frac{1}{h^{3n}{h}^{3(N-n)}}{e}^{-\beta H_n}\underset{(S^c{)}^{N-n}}{\int }{e}^{-\beta H_{N-n}}{d}^{N-n}\Xi \right]$$

Full text

In this work we derive a mathematical model for an open system that exchanges particles and momentum with a reservoir from their joint Hamiltonian dynamics. The complexity of this many-particle problem is addressed by introducing a countable set of $n$-particle phase space distribution functions just for the open subsystem, while accounting for the reservoir only in terms of statistical expectations. From the Liouville equation for the full system we derive a set of coupled Liouville-type equations for the $n$-particle distributions by marginalization with respect to reservoir states. The resulting equation hierarchy describes the external momentum forcing of the open system by the reservoir across its boundaries, and it covers the effects of particle exchanges, which induce probability transfers between the $n$- and $(n+1)$-particle distributions. Similarities and differences with the Bergmann-Lebowitz model of open systems (P.G. Bergmann, J.L. Lebowitz, Phys. Rev., 99:578–587 (1955)) are discussed in the context of the implementation of these guiding principles in a computational scheme for molecular simulations.

I Introduction

Open systems of microscopic particles that exchange energy and matter with a large environment are ubiquitous in the physical world. Strictly speaking, closed systems that exchange neither particles nor momentum or energy with their environment don’t even exist. Examples of open systems span from naturally-occurring systems to artificially engineered systems of modern and future technologyQuian (2007); Bordin et al. (2012); Abramo (2003); Politzer (1996); Herzog and Olshanii (1997); Delle Site (2018a). The description of the complexity of such systems often requires computer simulations for the prediction of their physical properties. Yet, the computational study of open systems comes with a higher degree of complexity than that of closed ones: Energy/momentum and particle exchanges with the environment render obsolete the established computational protocols, in particular those of Molecular Dynamics (MD), which rely on just these conservation laws by design Frenkel and Smit (2002); Tuckerman (2010).

Meanwhile, the latest generation of MD algorithms based on the concept of adaptive molecular resolution provide a generic root model for computational approaches to open systems Praprotnik, Delle Site, and Kremer (2008); Ensing et al. (2007); Heyden and Truhlar (2008); Wagoner and Pande (2013); Mones et al. (2015); Delle Site and Praprotnik (2017); Delle Site (2018b); Ciccotti and Delle Site (2019): Different regions of space are treated at different molecular resolution while allowing molecules/particles to freely move in space and change resolution accordingly. As a consequence, one has a prototype scheme for simulating a high resolution region coupled to a simplified environment (reservoir) with the exchange of energy and matter.

Applications of such schemes to real world systems led to accurate numerical results for a large class of physio- and bio-chemical systems Nielsen, Moore, and Ensing (2010); Delgado-Buscalioni, Sablić, and Praprotnik (2015); Sablić, Praprotnik, and Delgado-Buscalioni (2016); Fiorentini et al. (2017); Krekeler and Delle Site (2017); Netz, Potestio, and Kremer (2016); Shadrack Jabes et al. (2018); Shadrack Jabes, Klein, and Delle Site (2018); Shadrack Jabes and Delle Site (2018); Zavadlav et al. (2018); Wagoner and Pande (2018). In this perspective a mathematical formalization of the model is highly desirable as it helps to avoid situations in which largely empirical definitions of calculated quantities lead to artificial results or to their misleading interpretation. For example, a key question posed in recent years regards the proper definition of time correlation functions in the high resolution region for those degrees of freedom that are not present in the simplified reservoir Agarwal et al. (2015); Agarwal and Delle Site (2015); Delle Site (2016); Agarwal and Delle Site (2016); Agarwal, Clementi, and Delle Site (2017). The definition of a proper time evolution equation for the phase space probability densities, i.e., the analogue of a Liouvillian, is mandatory in this case. In other words, a mathematical model that formalizes the algorithms and makes their implied probability density evolution explicit became necessary.

In this perspective, the Adaptive Resolution Simulation scheme (AdResS) Praprotnik, Delle Site, and Kremer (2005) has been mapped onto a formal model of open systems known as the Bergmann-Lebowitz model Lebowitz and Bergmann (1957); Bergmann and Lebowitz (1955). This model, whose details are reported later on in this paper, defines a Liouville-type equation for the probability distribution of the open system via a stochastic coupling with the reservoir based on the idea of impulsive system-reservoir interactions. The mapping between AdResS and the Bergmann-Lebowitz model is unfortunately only qualitative Agarwal et al. (2015), but it provides an attractive reference for the structure of the algorithm and led to favorable comparisons of time correlations between AdResS and fully atomistic simulations Ciccotti and Delle Site (2019). The key characteristic of AdResS and of similar schemes is an exchange of particles and momentum/energy with some reservoir according to classical Newtonian dynamics. In contrast, the Bergmann-Lebowitz model features a stochastic coupling to describe this exchange. The aim of the present work is to derive a mathematical model that is similar in spirit to the Bergmann-Lebowitz model for open systems but based on the first principles of Newtonian mechanics, and that therefore accounts for interactions with the reservoir as they occur in the computational model.

We propose to describe the evolution of the probability distribution of the open system by introducing a total of $N$ $n$-particle phase space distribution functions $f_{n}$, where $N$ is the number of particles in the entire System and $n\in\{0,...,N\}$ denotes one of the possible numbers of particles that reside in the considered open subsystem at any instance of time. A hierarchy of Liouville-type equations, one for each of the $f_{n}$, is derived directly from the full Liouville equation for the probability distribution $F_{N}$ of the entire $N$-particle system as a whole. The derivation proceeds by integrating out, i.e., marginalizing w.r.t., the degrees of freedom of the $N-n$ particles that reside outside the open system subregion.

The result is a coupled system of Liouville-type equations for the entire hierarchy of distributions $(f_{n})_{n=0}^{N}$. The coupling arises naturally in the form of external forcing across the boundaries of the open system, and it is due to the exchange of particles and momentum/energy with the outside reservoir. A solution of this coupled system describes the time evolution of all $n$-particle probability distributions of the open system simultaneously. The probabilistic exchange of particles with the reservoir is accounted for by the probabilities of transitions from $n$- to $n+1$-particle states across the hierarchy of equations. Thus, at each time $t$ the set of distributions $(f_{n})_{n=0}^{N}$ characterizes the state of the open system in the sense of statistical mechanics.

If one is not interested in the detailed particle-based description of the reservoir, then the key theoretical step is the realization that in absence of detailed information regarding the time evolution of the reservoir state, probability distributions for the state of the open system can only be described by statistical averages over the reservoir’s microscopic degrees of freedom. Adopting the simplest possible approach, the derivation in this paper considers the expectations of the $n$-particle open system distributions with respect to the statistics of the reservoir.

The key point in connecting this model with the AdResS computational approach is the introduction of a boundary/surface layer around the open system whose particles may enter the system or exit fully into the reservoir, and which are close enough to the open system boundary that their pair potential momentum exchange with particles inside the system is non-negligible during a time step. This boundary layer, which is the central signature of the AdResS scheme, is responsible for modelling the statistics of the outside reservoir which enters the theory only through its single- and two-particle distributions.

We follow a procedure similar to the so-called Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy scheme (BBGKY model) Bogoliubov (1946), where the Liouville equation in a system of $N$ particles is written for a subset $s$ of particles whose probability distribution function is the $N-s$-marginalized probability distribution function of the total system. However, our equations are substantially different in form and meaning from those of the BBGKY model, due to the partitioning of position space in our model into the domain of the open system and its complement.

The paper is organized as follows: First we define the terms of the problem, that is we discuss the Liouville equation for a hierarchy of $N$-particle distribution functions specifically aimed at the description of an open system. Next we show how this hierarchy is constructed and present the explicit expression of the resulting Liouville-type equations. A verification for its physical consistency is provided by discussing the situation of equilibrium and the resulting Grand Canonical distribution. To stress the difference between the proposed approach and fully stochastic system-reservoir coupling models, we discuss similarities and differences with the Bergmann-Lebowitz model Lebowitz and Bergmann (1957); Bergmann and Lebowitz (1955). Finally, the model is put in relation with computational approaches aiming at simulating open systems, in particular we will discuss the correspondence between our model and the principles of the AdResS algorithm.

II Liouville-type equations for the $N$-particle distribution functions of an open system

II.1 Liouville-type equations

In general, a Liouville equation is the mathematical formalization of the concept that the statistical weight of an individual member of an ensemble of realizations of a dynamical system does not change in time. For a closed system, this is equivalent to stating that the probability density of initial states of the system is advected along trajectories of the dynamics, and, for an $N$-particle system, it is expressed by the phase space density transport equation:

[TABLE]

with $f(t,{\boldsymbol{q}},{\boldsymbol{p}})$ the distribution function in the $6N$-dimensional $({\boldsymbol{q}},{\boldsymbol{p}})$ phase space system, $v_{i}$ the velocity and $F_{i}$ the force acting upon the $i$-th particle, and analogously $(\vec{q}_{i},\vec{p}_{i})$ the $6$-dimensional position and momentum phase space point of the of the $i$-th particle. In this work we shall go beyond the standard application of Eq. 1 and explore its consequences for the statistics of a finite size open subsystem which, in the course of time, may host in principle any number of particles from [math] to $N$ as part of the dynamics. We will therefore have to address the time evolution of an entire hierarchy of $n$-particle distributions $f_{n}(t,\boldsymbol{X}^{n})$  $(n=0,...,N)$.

Notice that this hierarchy will differ from the well-known BBGKY hierarchy for $n$-particle marginal distributions of the original system in that we restrict to $n$ particle distributions conditioned on their and only their residence in the open system subdomain $\Omega$.

II.2 General set-up for a hierarchy of phase space density functions

As underlined before, the interest of this work is in the mathematical modelling of an open system of particles subject to Hamiltonian dynamics internally and subject to external forcing from its boundaries. The Hamiltonian for a situation that has $n$ particles in the system’s spatial domain $\Omega$ reads:

[TABLE]

where $V(\vec{q})$ is a two-particle potential, $X_{i}\equiv(\vec{q}_{i},\vec{p}_{i})\in S=\Omega\times\mathbb{R}^{3}$ are the $i$th particle’s physical and momentum space coordinates, respectively, and $M$ is the mass of an individual particle.

By saying that the system is “open” we mean that it can exchange particles with a, typically large, surrounding reservoir. Furthermore, we aim to account for the possibility that particles which reside inside the system’s boundary can interact with reservoir particles outside via inter-particle forces. These interactions are responsible for an exchange of kinetic energy/momentum between the system and its surroundings.

By definition of an open system, details of the “state of the world” at any time are known only for the particles inside the system. Properties of the “rest of the world” have to be assumed or prescribed by the modeller, and they depend on which kind of embedding she or he is thinking of. In the present work we assume that the outside world is described by suitable statistical information. That being the case, the open system description to be developed can only be statistical with respect to the outside world, the simplest case being a description in terms of expectations.

In this sense, the present suggestion of a “Liouville-type” system of equations for open systems differs substantially from the classical Liouville equation for a closed particle system with deterministic dynamics: the evolution of phase space densities is, in this latter case, fully deterministic, while here we have in principle even a statistical ensemble of hierarchies of open system phase space densities. Restricting to expectations relative to the outside world, we reduce this complexity to the task of handling the evolution of just a finite, or in the worst case countable, set of $n$-particle phase space distribution functions. This sequence of probability distributions is similar to the BBGKY hierarchy of $n$-particle probability densities with the key distinction that we consider $n$-particle densities conditional upon $n$ and only $n$ particles residing inside the open system’s domain.

A mathematical characterization of the system is given through an extended Liouville equation for the probability density:

[TABLE]

with

[TABLE]

of finding, at time $t$, any $n$ particles in the system with locations $(X_{i})_{i=1}^{n}$ in the particles’ phase space. For a “universe” with a total of $N$ particles in a domain $U$ that contains $\Omega$ and that has an $N$-particle density $F_{N}(t,\boldsymbol{X}^{N})$ for identifiable particles, we let

[TABLE]

with $\Omega_{c}=U\backslash\Omega$ and

[TABLE]

The combinatorial prefactor in (5) results from the fact that the (classical) particles are indistinguishable, so that the probability density $F_{N}$ is symmetric w.r.t any permutation of the particles. For future reference, we let

[TABLE]

denote the full $N$-particle distribution marginalized with respect to its last $N-n$ coordinates, and without account of the indistinghuishability of the particles.

The distributions in (3) form a hierarchy of coupled $n$-particle density functions since, for any realization of the system, the total number of particles inside of its domain is generally time dependent. As a consequence, we need distributions for any particle number $n\in\{0,1,...,N\}$ for a complete statistical description of the system. The normalization condition for $F_{N}$, i.e., $\int_{S^{N}}F_{N}\,d\boldsymbol{X}^{N}=1$, implies the normalization condition for the $n$-particle hierarchy

[TABLE]

The normalization condition of (8) is routinely used in popular textbooks of statistical mechanics when partitioning a large system in (open) small subsystems, see e.g. Refs.37; 8. In section IV it will be explicitly verified that (8) is consistent with the derivation of the grand canonical distribution function. In physics the sum is usually extended to $\infty$; such an approach is based on qualitative arguments (thermodynamic limit), and it will be used here in section IV. In mathematics the extension of the sum to $\infty$ is not straightforward and it represents an interesting problem. However given the current focus, this discussion goes beyond the scope of this paper.

III Derivation of the Liouville-type equation hierarchy

Here we consider an open subsystem occupying the spacial domain $\Omega$ embedded in a much larger “Universe” $U\subset\mathbb{R}^{3}$ of $N$ particles. The starting point of the derivation is the Liouville equation for the $N$-particle distribution function. This is a standard Liouville equation for a closed system of $N$ particles and probability distribution $F_{N}(t,\boldsymbol{X}^{N})$, that is, the conservation law for the probability of the universe. In mathematical terms this means a definition of $F_{N}(t,\boldsymbol{X}^{N})$ as:

[TABLE]

and the corresponding equation of conservation:

[TABLE]

and as a consequence, in the form known as Liouville equation:

[TABLE]

where

[TABLE]

is the total potential energy of the system as a function of the positions of all of its $N$ particles. True to the definition of the open system $n$-particle distributions in (5), we marginalize this general evolution equation w.r.t. $N-n$ particles residing outside of the open system subdomain $\Omega$.

III.1 Marginalizing w.r.t. outside particle momentum spaces

Here we consider all possible configurations conditioned upon $n$ and only $n$ particles residing within the open subsystem and integrate the Liouville equation (11) over the $N-n$ degrees of freedom of the remaining particles outside the open subsystem. After this operation and by the definition of $f_{n}$ in (5), the left hand side yields the time derivative $\frac{\partial f_{n}}{\partial t}$ of the $n$ particle distribution. The terms on the right of the equation will be shown below to yield the expected Liouville-type transport terms for $f_{n}$ itself as well as coupling terms between $f_{n}$ and other functions $f_{k}$ for $k\not=n$.

Let $\Xi_{k}=(\vec{q}_{k},\vec{p}_{k})$ for some $k$ such that {{\color[rgb]{0,0,0}n+1}}<k\leq N and consider the marginalization of the term $\nabla_{\vec{p}_{k}}\cdot\left(\vec{F}_{k}\,F_{N}\right)$ with respect to $\vec{p}_{k}$. Integrating the term over any ball $B_{r}(0)$ of momentum space, we find

[TABLE]

In the limit of the radius of the ball tending to infinity the term will generally vanish because the probability $F_{N}$ will decay sufficiently rapidly for large momenta. For the equilibrium distribution this is obvious, since in this case $F_{N}$ factorizes into the momentum and position terms, and the momentum distribution is proportional to $\exp(-(\vec{p}_{k}/M)^{2}/2kT)$. All the terms in (13) with {{\color[rgb]{0,0,0}n+1}}\leq k\leq N follow this reasoning, so that they all vanish identically.

In contrast, for $1\leq i\leq n$, we have to carefully account for the dependencies of the total potential energy on the open system and universe coordinates. To this end, we assume here that the total potential is the sum of pair interactions as described in (12). Under this assumption, let us analyze the behavior of one of the terms under this double-sum when the momentum flux divergence in the Liouville equation is marginalized. If both, $i,j\leq n$, then

[TABLE]

provided the probability density is sufficiently smooth to allow for an exchange of integration and differentiation w.r.t. complementary coordinates. Thus, taking into account the particle exchange symmetry again (see (7)), we obtain the standard momentum flux divergence contribution that we would expect in an $n$-particle Liouville equation.

If, however, in one of the terms in (12) only $\vec{q}_{i}\in\Omega$ but $\vec{q}_{j}\in\Omega_{c}$, then the marginalization of the corresponding term in the Liouville equation reads

[TABLE]

Due to the dependence of the force on $\vec{q}_{j}$ and the fact that we are marginalizing over that variable, this term cannot be straightforwardly expressed in terms of the hierarchy of distributions $(f_{n})_{n=0}^{N}$. An explicit step of modelling is required for closure.

As this paper is to describe a general framework but is not focused on some specific system, let us suggest here one plausible closure model and leave the discussion of more sophisticated options for future work: Suppose that the pair interaction $V(r)$ is relatively short-range in comparison with $\text{diam}(\Omega)$, so that the pair interactions are felt only close to the open system’s boundary. Suppose further that the probability density of finding $n$ particles in states $\left(\boldsymbol{X}^{n-1},X_{i}\right)\in S^{n}$ and one other outer particle in $X_{j}$ is given by $f_{n}\left(\boldsymbol{X}^{n-1},X_{i}\right)f^{\circ}_{2}(X_{j}|X_{i})$, where $f^{\circ}_{2}(X_{\text{out}}|X_{\text{in}})$ is a known or modelled conditional distribution for joint appearances of an outer particle given the state of an inner one. The key modelling assumption here is that an inner particle $i$ whose position $q_{i}$ is sufficiently close to the boundary $\partial\Omega$ for it to feel the pair interaction from the outer particles, this conditional distribution is approximately independent of the states $\boldsymbol{X}^{n-1}$ of the other $n-1$ particles within the open system.

Given this modelling assumption, (15) results in a mean-field expression for the momentum flux divergence term,

[TABLE]

where

[TABLE]

is the mean field force exerted by the outer particles onto the $i$th inner particle. The derived expression is invariant w.r.t. particle permutations which accounts for a multiplicity of $N!$. At the same time, we keep $f_{n}(t,\boldsymbol{X}^{n})$ as a symmetric function w.r.t. particle permutations inside the open system, so that the statistical weight of a particular configuration is only $1/n!$. Furthermore, the marginalization integral is overcounting the statistical weight of the expression by a factor of $(N-n)!$ which corresponds to all permutations of the outer particles. Therefore, multiplication of (16) by ${{N}\choose{n}}$ yields a contribution

[TABLE]

to the evolution equation of $f_{n}$, i.e., the outer particles exert a meanfield force onto the particles within the open system.

III.2 Marginalizing w.r.t. outside particle position spaces

Next we consider the physical space transport terms \sum_{i=1}^{N}\nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}F_{N}\right) from (11) w.r.t. the outer particles. Consider first those terms in the sum for which $i\in 1,...,n$. For these,

[TABLE]

and we obtain the analogous transport term for the $n$-particle open system Liouville equation. For {{\color[rgb]{0,0,0}n+1}}\leq i\leq N, however, one of the integrals will be over $\Xi_{i}\in S^{c}$, and for the corresponding terms

[TABLE]

Leaving the exact definition of ${\widehat{f}}_{n+1}$ to be discussed in the next paragraph, we notice that the factor $(N-n)$ in the second line arises because all the terms from the first line for i\in{{\color[rgb]{0,0,0}n+1}},...,N are identical, and that the negative sign in the third line arises after we have applied Gauß’ theorem to the integral over $\Omega_{c}$, whereas we have used the outer normal $\vec{n}$ of $\Omega$ (not its complement) in the formula. Moreover, we have used that ${{N}\choose{n}}(N-n)=(n+1){{N}\choose{n+1}}$ to let ${\widehat{f}}_{n+1}$ acquire the same scaling with $N$ and $n$ as ${f}_{n+1}$.

Next we need to distinguish the values of the marginalized density ${\widehat{f}}_{n+1}$ that are to be used for the two possible signs of \left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}\right). To this end, we recall that the theory of characteristics states that trajectories of the particle system carry with them their initial values of the density (see, e.g., Godlewski and Raviart (1996)Godlewski and Raviart (1996), section V, pp. 417–460). Therefore, when particles exit the open system, that is, \left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}\right)>0, then ${\widehat{f}}_{n+1}=f_{n+1}$, whereas for \left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\cdot\vec{n}\right)<0 information is entering the open system from the outer universe, and one needs to decide upon a model to represent the statistical information transport in that case. To emphasize this all-important degree of freedom in open system modelling, we follow two options in the sequel.

A most straightforward approach assumes an equilibrium state of the universe that is statistically independent of the open system, and that amounts to letting ${\widehat{f}}_{n+1}=f_{n}f_{1}^{\circ}$, where $f_{1}^{\circ}$ is the single particle (equilibrium) density assumed for the reservoir. Thus,

[TABLE]

A second option is geared towards consistency with the thermodynamic and large system size limits for open systems. This amounts to imposing grand canonical distribution for state space trajectories that enter the open system from outside and reads

[TABLE]

See an explicit representation of the grand canonical distribution $f_{n+1}^{\text{GC}}$ in (33) below.

III.3 Liouville-type equation for the $\boldsymbol{n}$-state density

From the previous subsections we collect the evolution equation

[TABLE]

where

[TABLE]

is the total force onto the $i$th particle exerted by the remaining $n-1$ particles within the system through the potential interaction, and the coupling terms within the $n$-particle hierarchy of distributions are

[TABLE]

where $f_{1}^{\circ}$ is the single particle (equilibrium) density assumed for the reservoir and

[TABLE]

IV Grand Canonical Equilibrium Distribution function as verification of consistency of the model

In equilibrium, and under the hypothesis of short range interaction, that is $V(\vec{q}_{i}-\vec{q}_{j})$ decays fast as $|\vec{q}_{i}-\vec{q}_{j}|$ increases, the standard stationary Grand Canonical distribution is automatically obtained by the very definition of $f(t,\boldsymbol{X}^{n})$ as verified next.

Let us consider the total system (i.e., the universe) with $N$ particles interacting through the Hamiltonian: $H_{N}=\sum_{i=1}^{N}\frac{\vec{p}_{i}^{2}}{2M}+\sum_{i=1}^{N}\sum_{j\not=i}^{N}\frac{1}{2}V(\vec{q}_{j}-\vec{q}_{i})$, then its partition function reads:

[TABLE]

with $h$ the Planck constant as usual, and $\beta=\frac{1}{k_{B}T}$, $k_{B}$ the Boltzmann constant. The probability distribution for a given configuration $\boldsymbol{X}^{N}$ is then written as:

[TABLE]

as a consequence, following the definition of (5), in equilibrium one obtains the probability distribution of a subsystem of $n$ particles in a domain $\Omega$ as:

[TABLE]

where $V_{cx}(\vec{q}_{\Omega}-\vec{q}_{\Omega^{c}})=\sum\limits_{i=1}^{n}\sum\limits_{j=n+1}^{N}V(\vec{q}_{i}-\vec{q}_{j})$ is the interaction potential between the $n$ particle inside and the $N-n$ particles outside.

Under the hypotheses that (i) $N-n\gg n$, (ii) $\Omega^{c}\gg\Omega$, and (iii) $V_{cx}(\vec{q}_{\Omega}-\vec{q}_{\Omega^{c}})$ is negligible compared to $\sum\limits_{i<j}^{n}V(\vec{q}_{i}-\vec{q}_{j})$ for $\vec{q}_{i},\vec{q}_{j}\in\Omega$ – as justified, e.g., in case of short range potentials – the total Hamiltonian can be factorized as: $H_{n}+H_{N-n}$. Then the expression in (29) can be reduced to:

[TABLE]

that is:

[TABLE]

with $Q_{N-n}=\frac{1}{h^{3(N-n)}(N-n)!}\int\limits_{(S^{c})^{N-n}}e^{-\beta H_{N-n}}d^{N-n}\vec{q}d^{N-n}\vec{p}$.

Next, following well established textbooks of statistical mechanics (see e.g. Huang (1986)), one has:

[TABLE]

with $A$ the Helmholtz free energy at given number of particles, volume and temperature. Since we assume $N\gg n$ and $\Omega_{c}\gg\Omega$, this yields: $A(N-n,U\backslash\Omega,T)-A(N,U,T)\approx n\mu+\Omega P$, with $\mu=\frac{\partial A}{\partial n}$ the chemical potential and $P=-\frac{\partial A}{\partial\Omega}$ the pressure. As a consequence one gets the standard Grand Canonical probability distribution function:

[TABLE]

and, as usual, the normalization condition:

[TABLE]

implies the thermodynamic definition of the grand potential $\Phi=P\Omega=k_{B}T\ln Z(\mu,\Omega,T)$, with $Z(\mu,\Omega,T)=\sum_{n=0}^{\infty}\int\limits_{\Omega^{n}}\int\limits_{(\mathbb{R}^{3})^{n}}\frac{1}{n!h^{3n}}e^{-\beta(H_{n}-\mu n)}\ d^{n}\vec{p}\ \ d^{n}\vec{q}$.

One can also verify that given the physical and mathematical hypothesis stated in the derivation above, the Grand Canonical probability distribution function of (33), satisfies (23). In fact in equilibrium, one would have:

[TABLE]

the left hand side equals zero, because \nabla_{\vec{q}_{i}}\cdot\left({{\color[rgb]{0,0,0}\vec{v}_{i}}}\frac{e^{-\beta(H_{n}-\mu n+P\Omega)}}{n!h^{3n}}\right)=\vec{F}_{i}\cdot{{\color[rgb]{0,0,0}\vec{v}_{i}}}f_{n}(\boldsymbol{X}^{n}) and \nabla_{\vec{p}_{i}}\cdot\left(\vec{F}_{i}\frac{e^{-\beta(H_{n}-\mu n+P\Omega)}}{n!h^{3n}}\right)=-\vec{F}_{i}\cdot{{\color[rgb]{0,0,0}\vec{v}_{i}}}f_{n}(\boldsymbol{X}^{n}). Next, the hypothesis that $V_{cx}(\vec{q}_{\Omega}-\vec{q}_{\Omega^{c}})$ is negligible, implies that $\vec{F}_{\text{av}}(\vec{q}_{i})$ is negligible, thus in equilibrium, $\Psi_{n}$ is a negligible “surface” effect. Finally, for $\Phi_{n}^{n+1}$, which describes the net flux of $n$-particle probability into/out of the open system, is trivially zero for grand canonical statistics of $f_{n+1}$ when we adopt the model from (22). For the other option in (21) we see that exact cancellation of this term would only occur if $f^{\text{GC}}_{n+1}-f_{n}^{\text{GC}}f_{1}^{\circ}=0$, which is generally not the case, since the Grand Canonical statistics does not factorize with respect to its position space dependence. Still, for large system size we can argue as for $\Psi_{n}$ and accept that the term is negligible to the extend that the surface-to-volume ratio of the open system is small.

The important aspect to be registered here is that in reality the “true” $f_{n}(t,\boldsymbol{X}^{n})$ that one obtains in equilibrium, includes explicitly the surface effects – as small as they might be – without any hypothesis on the behavior of the potential. The Grand Canonical distribution, in turn, is only a physical approximation that becomes rigorous in the thermodynamic limit for large $n$ and large system size. In this sense it will be interesting to see, in section V, the differences between the present open system model and the BL model, which assumes $V_{cx}(\vec{q}_{i}-\vec{q}_{j})=0$ when $\vec{q}_{i}\in\Omega$ and $\vec{q}_{j}\in\Omega_{c}$ and relies on the impulsive interaction hypothesis. Similarly, in section VI we discuss differences with computational algorithms in which $V_{cx}(\vec{q}_{\Omega}-\vec{q}_{\Omega^{c}})$ plays a relevant technical role.

V Similarities and differences between the Bergmann-Lebowitz and the present open system model

The key idea of the Bergmann-Lebowitz open system model (Lebowitz and Bergmann, 1957; Bergmann and Lebowitz, 1955) is to discretize the interaction of an open system with its environment. They allow for impulsive interactions of the system with statistically undisturbed reservoirs at discrete points in time, the interactions being induced by a suitable interaction kernel. In essence, the interaction with the reservoir leads to a discrete transition of the system from a state characterized by $n$ particles $\boldsymbol{X}^{n}$ to a state characterized by $m$ particles ($\boldsymbol{Y}^{m}$). As a hypothesis/model, the macroscopic thermodynamic variables of the reservoir are not influenced by the system and molecules entering into the system from the reservoir can only have, in thermal equilibrium, velocities consistent with the temperature of the reservoir.

These concepts are formalized via the the introduction of a stochastic kernel, $K_{nm}(\boldsymbol{Y}^{m},\boldsymbol{X}^{n})$, whose action allows for the transition from a microscopic state, $\boldsymbol{X}^{n}$, to another state, $\boldsymbol{Y}^{m}$ (and vice versa). In essence, such kernel corresponds to the probability per unit time that the system at $\boldsymbol{X}^{n}$ makes a transition to $\boldsymbol{Y}^{m}$, caused by the interaction between the system and the reservoir. The total system-reservoir interaction becomes: $\sum_{m=0}^{\infty}\int d\boldsymbol{Y}^{m}[K_{nm}(\boldsymbol{X}^{n},\boldsymbol{Y}^{m})f_{m}(\boldsymbol{Y}^{m},t)-K_{mn}(\boldsymbol{Y}^{m},\boldsymbol{X}^{n})f_{n}(\boldsymbol{X}^{n},t)]$. The time evolution of the probability, $f_{n}(\boldsymbol{X}^{n},t)$, is described by a Liouville equation where the Poisson bracket on the r.h.s. is augmented by a term, in integral form, describing the exchange with the reservoir:

[TABLE]

The condition of flux balance defines the equilibrium for this model through

[TABLE]

and in equilibrium, the stationary solution for $f_{n}(\boldsymbol{X}^{n})$ is the Grand Canonical probability density: $f_{n}(\boldsymbol{X}^{n})=\frac{1}{Q}e^{-\beta H_{n}(\boldsymbol{X}^{n})+{\beta}\mu n}$ where $\beta=1/k_{B}T$ and $\mu$ the chemical potential. The condition of flux balance is necessary and sufficient for obtaining the Grand Canonical probability density as stationary solution.

The key difference with the proposed model of this work is that we do not assume a discontinuous impulsive interactions of the system with the reservoir, but we allow for a dynamical continuous transition between different states described by the Hamiltonian of the “entire universe”. Proper conditions of interaction of the open system with the reservoir are derived analytically, and yield the terms, $\Psi_{n}+\Phi_{n}^{n+1}$, which link neighboring levels of the $n$-hierarchy of open system state space densities. These terms represent the flux of information in and out of the open system’s domain, $\Omega$. An interesting application of the developed physico-mathematical open system model lies in its explanatory power of interpreting computational schemes that simulate an open subsystem with detailed atomistic particle representations embedded in a reduced complexity environment and the dynamical exchange of energy and particles between them. In the next section we will discuss in detail the connection between the proposed model and a well established computational algorithm for simulation of open molecular systems.

VI Relation between the Adaptive Resolution Molecular Dynamics scheme and the proposed open system model

Recent developments in computational physics have contributed to the construction of algorithms for molecular dynamics simulation of systems that exchange matter with a reservoir. Fig. 1 shows the set up of the latest version of the Adaptive Resolution Simulation method (AdResS) together with the corresponding spatial partitioning of the simulation box Krekeler et al. (2018); Delle Site et al. (2019). The fully atomistic region is interfaced with a transition region which in turn is interfaced with a region of non-interacting particles (tracers). Upon entering the transition region a tracer particle acquires an atomistic identity and while crossing this region the molecule equilibrates with the full atomistic environment. In the other direction, a molecule in the transition region entering into the tracer region becomes a non-interacting particle and loses its atomistic character; the process of changing identity occurs in an abrupt manner.

The similarity between the theoretical model introduced in these notes and the computational implementation of AdResS can be brought forward noticing that the transition region of AdResS is essentially a region of the reservoir, where, according to the potential cutoff used, the terms $\Psi_{n}$ are defined and where the surface effects of $\Phi_{n}^{n+1}$ take place.

Moreover, in the transition region of AdResS an additional condition is imposed for the proper statistical exchange of particles between the system and reservoir. Such condition is derived from first principles of statistical mechanics with the aim of conserving the equilibrium statistical properties in AT region. This condition is imposed in the AdResS model through an additional thermodynamic force that acts on particles only in the transition region Fritsch et al. (2012); Wang et al. (2013). In conclusion, the AdResS scheme can be interpreted as a particle-based implementation of the open system model proposed in section II above. Obviously, the atomistic region represents the open system, whereas the transition ($\Delta$) region is to represent the outside world. For distance-truncated interaction potentials, the only part of the outside world which the system sees is, in fact, a finite thickness layer of particles that covers the cut-off length. Therefore, the thickness of the $\Delta$-region should be comparable to this cut-off length. Particles that reside beyond the $\Delta$-region, which undergo simplified coarse-grained dynamics (they are just passive tracers in the present notes), merely serve as a particle reservoir that is needed to make sure the statistical balancing mechanisms active in the $\Delta$-region always have a sufficient supply of particles. The main task of the $\Delta$-region in AdResS is to generate the desired outside world statistics. For a system as described above, whose dynamics involves only particle pair interactions, this means that the two-particle distribution $f^{\circ}_{2}$ is to be established. Under the AdResS philosophy, this distribution is to emerge from explicit particle dynamics in the $\Delta$-region which, due to the activity of a thermostat, can be expected to be ergodic, and which, due to the action of the thermodynamic force, is guaranteed to at least provide the correct mean outside world particle density to particles within the open system. Moreover, since in the $\Delta$-region the particles follow, except for these two effects, the same Hamiltonian dynamics as the particles inside, it is plausible that a particle ensemble with the correct mean density and temperature will also adjust in space correctly so as to reflect the radial distribution or particle pair position distributions. All the mechanisms of how the coupling between distributions of different particle numbers interact and what are the proper boundary conditions carried by particles entering the domain are taken care of by AdResS automatically. This is because AdResS is particle based and therefore implements the phase space transport of the distribution functions directly in a “Lagrangian fashion” by generating stochastic trajectories. In summary, a clear comparison between the Bergmann-Lebowitz model and our model of open system regarding the formalization of real molecular simulation techniques, would lead to the following conclusions: the Bergmann-Lebowitz model is optimal for those simulation techniques based on stochastic sampling of the phase space configuration, (i.e. Monte Carlo methods), instead our model is optimal for techniques based on molecular dynamics schemes where equations of motions and molecular trajectories are determined at each instant of the simulation. In fact in the first case the system-reservoir coupling kernel, expressing the transition probability in the phase space, can be straightforwardly modeled with a Monte Carlo procedure, while in the second case the system-reservoir coupling term can be expressed through two-particle forces and particle-particle correlation functions, automatically calculated in a molecular dynamics scheme.

VII Conclusions

We have proposed a physico-mathematical model for open systems that exchange particles and energy with the external world. The procedure of reducing a large “universe” system to a smaller (sub)system is similar to the marginalization procedure of the BBGKY hierarchy scheme that defines a hierarchy of probability functions of the system. However, our model differs substantially from the BBGKY model because in the BBGKY model the marginalization is done independently of the position of the particles, whereas in our model the marginalization is applied with the constraint that particles are located outside a predefined region. As a consequence, in our case terms involving directly the partitioning of space (i.e. surface/boundary integrals) appear in the equations while this is not the case of the BBGKY model.

The equations obtained model a system of particles that continuously, in a dynamical fashion, exchanges energy and particles with the reservoir. This characteristic makes our model differ from the well-established Bergmann-Lebowitz model of open systems, which is based instead on the hypothesis of discrete impulsive and stochastic interaction between the system and the reservoir. The motivation for the development of our model lies on the fact that modern molecular dynamics schemes for open systems are based on a dynamical and continuous exchange of information with the exterior and thus they can only be approximately described by the Bergmann-Lebowitz model.

Instead, as we discuss in the paper, our model provides a closer formalization of the computational algorithms with the specific example of the mapping of our theoretical model to the AdResS method.

Acknowledgements.

This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 “Scaling Cascade in Complex Systems,” Project Number 235221301, Project C01 “Adaptive coupling of scales in molecular dynamics and beyond to fluid dynamics.” Data availability Statement: The data that supports the findings of this study are openly available on arXiv, reference number: 1907.07557