# High-Resolution Depth Maps Based on TOF-Stereo Fusion

Vineet Gandhi†, Jan Čech, and Radu Horaud
Inria Grenoble Rhône-Alpes, Montbonnot Saint-Martin, France
{vineet.gandhi, jan.cech<EMAIL_ADDRESS>
${{\dagger}}$ V. Gandhi acknowledges support from the Erasmus Mundus CIMET
master program.

###### Abstract

The combination of range sensors with color cameras can be very useful for
robot navigation, semantic perception, manipulation, and telepresence. Several
methods of combining range- and color-data have been investigated and
successfully used in various robotic applications. Most of these systems
suffer from the problems of noise in the range-data and resolution mismatch
between the range sensor and the color cameras, since the resolution of
current range sensors is much less than the resolution of color cameras. High-
resolution depth maps can be obtained using stereo matching, but this often
fails to construct accurate depth maps of weakly/repetitively textured scenes,
or if the scene exhibits complex self-occlusions. Range sensors provide coarse
depth information regardless of presence/absence of texture. The use of a
calibrated system, composed of a time-of-flight (TOF) camera and of a
stereoscopic camera pair, allows data fusion thus overcoming the weaknesses of
both individual sensors. We propose a novel TOF-stereo fusion method based on
an efficient seed-growing algorithm which uses the TOF data projected onto the
stereo image pair as an initial set of correspondences. These initial “seeds”
are then propagated based on a Bayesian model which combines an image
similarity score with rough depth priors computed from the low-resolution
range data. The overall result is a dense and accurate depth map at the
resolution of the color cameras at hand. We show that the proposed algorithm
outperforms 2D image-based stereo algorithms and that the results are of
higher resolution than off-the-shelf color-range sensors, e.g., Kinect.
Moreover, the algorithm potentially exhibits real-time performance on a single
CPU.

##  I Introduction

An advanced computer vision system should be able to provide both accurate
color and depth information for each pixel at high resolution. Such a system
can be very useful for automated vision problems especially in the context of
robotics, e.g., for building dense 3D maps of indoor environments.

The 3D structure of a scene can be reconstructed from two or more 2D views via
a _parallax_ between corresponding image points. However, it is difficult to
obtain accurate pixel-to-pixel matches for scenes of objects without textured
surfaces, with repetitive pattern, or in the presence of occlusions. The main
drawback is that stereo matching algorithms frequently fail to reconstruct
indoor scenes composed of untextured surfaces, e.g., walls, repetitive
patterns and surface discontinuities, which are typical in man-made
environments.

(a) A TOF-stereo setup

(b) The TOF image is shown in the upper-left corner of a color image.

(c) The proposed method delivers a high-resolution depth map.

Figure 1: (a) Two high-resolution color cameras (2.0MP at 30FPS) are combined
with a single low-resolution time-of-flight camera (0.03MP at 30FPS). (b) A
$144\times 177$ TOF image and a $1224\times 1624$ color image are shown at the
true scale. (c) The depth map obtained with our method. The technology used by
both these camera types allows simultaneous range and photometric data
acquisition with an extremely accurate temporal synchronization, which may not
be the case with other types of range cameras such as the current version of
Kinect.

‘

Alternatively, active-light range sensors, such as time-of-flight (TOF) or
structured-light cameras, can be used to directly measure the 3D structure of
a scene at video frame-rates. However, the spatial resolution of currently
available range sensors is lower than high-definition (HD) color cameras, the
luminance sensitivity is poorer and the depth range is limited. The range-
sensor data are often noisy and incomplete over extremely scattering parts of
the scene, e.g., non-Lambertian surfaces. Therefore it is not judicious to
rely solely on range-sensor estimates for obtaining 3D maps of complete
scenes. Nevertheless, range cameras provide good initial estimates
independently of whether the scene is textured or not, which is not the case
with stereo matching algorithms. These considerations show that it is useful
to combine the active-range and the passive-parallax approaches, in a mixed
system. Such a system can overcome the limitations of both the active- and
passive-range (stereo) approaches, when considered separately, and provides
accurate and fast 3D reconstruction of a scene at high resolution, e.g.
$1200\times 1600$ pixels at 30 frames/second, as in Fig. 1.

###  I-A Related Work

The combination of a depth sensor with a color camera has been exploited in
several robotic applications such as object recognition [14, 21, 2], person
awareness, gesture recognition [11], simultaneous localization and mapping
(SLAM) [3, 16], robotized plant-growth measurement [1], etc. These methods
have to deal with the difficulty of noise in depth measurement and the
inferior resolution of range data as compared to the color data. Also, most of
these systems are limited to RGB-D, i.e., a single color image combined with
the range data. Interestingly enough, the recently commercialized
Kinect111http://www.xbox.com/en-US/kinect camera falls in the RGB-D family of
sensors. We believe that extending this model to an RGB-D-RGB sensor is
extremely advantageous because it can incorporate stereoscopic matching and
hence better deal with the problems mentioned above.

Stereo matching has been one of the most studied paradigms in computer vision.
Several papers, e.g., [19, 20] present an overview of existing techniques and
highlight recent progress in stereo matching and stereo reconstruction.
Algorithms based on a greedy local search are typically fast but frequently
fail to reconstruct the poorly textured regions or ambiguous surfaces. Global
methods formulate the matching task as a single optimization problem which
leads to minimization of an Markov random field (MRF) energy function of the
image similarity likelihood and a prior on the surface smoothness. These
algorithms solve some of the aforementioned problems of local methods but are
very complex and computationally expensive since optimizing an MRF-based
energy function is an NP-hard problem in the general case.

A practical tradeoff between the local and the global methods in stereo is the
seed growing class of algorithms [5, 6, 4]. The correspondences are grown from
a small set of initial correspondence seeds. Interestingly, they are not
particularly sensitive to wrong input seeds. They are significantly faster
than the global approaches, but they have difficulties in presence of non
textured surfaces; Moreover, in these cases they yield depth maps which are
relatively sparse. Denser maps can be obtained by relaxing the matching
threshold but this leads to erroneous growth, which is a natural tradeoff
between the accuracy and density of the solution. Some form of regularization
is necessary in order to take full advantage of these methods.

Recently an external prior-based generative probabilistic model for stereo
matching was proposed in [13, 18] for reducing the matching ambiguities. The
prior used was based on surface-triangulation on initially-matched distinctive
interest points in the images. Again, in the absence of textured regions, such
support points are either not available or are not reliable enough and the
priors are erroneous. Consequently, the methods produce artifacts in cases the
priors win over the data and the solution is biased towards the incorrect
priors. This clearly shows the need for more accurate prior models. [22]
integrates a regularization term based on the depth values of initially
matched ground control points in a global energy minimization framework. The
ground control points are gathered using an accurate laser scanner. A laser
scanner is difficult to operate and cannot provide range information fast
enough such that it can be used in a practical robotic application.

A TOF camera is based on an active sensor principle222http://www.mesa-
imaging.ch that allows 3D data acquisition at video frame rates, e.g., 30FPS
as well as accurate synchronization with any number of color
cameras333http://www.4dviews.com. A modulated near infrared light from the
camera’s internal lighting source is reflected by objects in the scene and
travels back to the sensor, where its precise time of flight is measured
independently at each of the sensor’s pixel by calculating the phase delay
between the emitted and the detected wave. A complete depth map of the scene
can be obtained using this sensor at the cost of very low spatial resolution
and coarse depth accuracy.

The fusion of TOF data with stereo data has been recently studied. For
example, [8] obtained a higher quality depth map, by a probabilistic ad-hoc
fusion of TOF and stereo data. Work in [23] merges the depth probability
distribution function obtained from TOF and stereo. However both these methods
are meant for improvement over the initial data gathered with the TOF camera
and the final depth-map result is still limited to the resolution of the TOF
sensor. The method proposed in this paper increases the resolution from 0.03MP
to the full resolution of the color cameras being used, e.g., 2MP.

The problem of depth map upsampling has been previously addressed. In [7] a
noise-aware filter for adaptive multi-lateral upsampling of TOF depth maps is
presented. The work described in [14] extends the model of [9] and
demonstrates that the object detection accuracy can be significantly improved
by combining a state-of-art 2D object detector with 3D depth cues. The
approach deals with the problem of resolution mismatch between range- and
color-data using an MRF-based super-resolution technique in order to infer the
depth at every pixel. The proposed method is slow: It takes around 10 seconds
to produce a $320\times 240$ depth image. All of these methods are limited to
depth-map upsampling using only a single color image and do not exploit the
added advantage offered by stereo matching, which can highly enhance the depth
map both qualitatively and quantitatively. Recently, [12] proposed a method
which combines TOF estimates with stereo in a semiglobal matching framework.
However, at pixels where TOF disparity estimates are available, the image
similarity term is ignored. This make the method quite susceptible to errors
in regions where TOF estimates are not precise, especially in textured regions
where stereo itself is reliable.

###  I-B Contributions

In this paper we propose a novel Bayesian method for incorporating range data
within a robust seed-growing algorithm for stereoscopic matching [5]. A
calibrated system composed of an active range sensor and a stereoscopic color-
camera pair [15], e.g., Fig. 1, allows the range data to be projected onto
each one of the two images, thus providing an initial sparse set of point-to-
point correspondences (seeds) between the two images. This initial seed set is
used in conjunction with the seed-growing algorithm proposed in [5]. The
novelty is that the range data are used as the vertices of a mesh-based
surface representation which, in turn, is used as a prior to regularize the
image-based matching procedure. The Bayesian fusion, between the mesh-based
surface (initialized from the sparse range data) and the seed-growing stereo
matching algorithm itself, combines the merits of the two 3D sensing methods
and overcomes the limitations outlined above. The proposed fusion model can be
incorporated within virtually any stereo algorithm that is based on energy
minimization and which requires proper initialization. It is, however,
particularly efficient and accurate when used in combination with match-
propagation methods.

The remainder of this paper is structured as follows: Section II describes the
proposed range-stereo fusion algorithm. Experimental results on a real dataset
and evaluation of the method, are presented in section III. Finally, section
IV draws some conclusions.

##  II The proposed algorithm

As outlined above, the TOF camera provides a low-resolution depth map of a
scene. This map can be projected onto the left and right images associated
with the stereoscopic pair, using the projection matrices estimated by the
calibration method described in [15]. Projecting a single 3D point $(x,y,z)$
gathered by the TOF camera onto the rectified images provides us with a pair
of corresponding points $(u,v)$ and $(u^{\prime},v)$ in the respective images.
Each element $(u,u^{\prime},v)$ denotes a point in the disparity space444The
disparity space is a space of all potential correspondences [19].. Hence,
projecting all the points obtained with the TOF camera gives us a sparse set
of 2D point correspondences. This set is termed as the set of initial support
points or TOF seeds.

These initial support points are used in a variant of the seed-growing stereo
algorithm [5, 4] which further grows them into a denser and higher resolution
disparity map. The seed-growing stereo algorithms propagate the
correspondences by searching in the small neighborhood of given seed
correspondences. Hereby, only a small fraction of disparity space is visited,
which makes the algorithm extremely efficient from a computational point of
view. The limited neighborhood also gives a kind of implicit regularization,
nevertheless the solution can be arbitrarily complex, since multiple seeds are
provided.

The integration of range data within the seed-growing algorithm required two
major modifications: (1) The algorithm is using TOF seeds instead of the seeds
obtained by matching distinctive image features, such as interest points,
between the two images, and (2) the growing procedure is regularized using a
similarity statistic which takes into account the photometric consistency as
well as the depth likelihood based on disparity estimate by interpolating the
rough triangulated TOF surface. This can be viewed as a prior cast over the
disparity space.

Algorithm 1 Growing algorithm with sensor fusion

0: Rectified images ($I_{L},I_{R}$), initial correspondence seeds $\cal S$,
image similarity threshold $\tau$.

1: Compute the prior disparity map $D_{p}$ by interpolating seeds $\cal S$.

2: Compute $\operatorname{simil}(s|I_{L},I_{R},D_{p})$ for each every seed
$s\in\cal S$.

3: Initialize empty disparity map $D$ of size $I_{L}$ (and $D_{p}$).

4: repeat

5: Draw seed $s\in\cal S$ of the best
$\operatorname{simil}(s|I_{L},I_{R},D_{p})$ value.

6: for each of the four best neighbors $i\\!\in\\!\\{1,2,3,4\\}$
$q_{i}^{\ast}=(u,u^{\prime},v)$ = $\underset{q\in{\cal
N}_{i}(s)}{\operatorname{argmax}}$ $\operatorname{simil}(q|I_{L},I_{R},D_{p})$
do

7: $c:=\operatorname{simil}(q_{i}^{\ast}|I_{L},I_{R},D_{p})$

8: if $c\geq\tau$ and pixels not matched yet then

9: Update the seed queue ${\cal S}:={\cal S}\cup\\{q_{i}^{\ast}\\}$.

10: Update the output map ${D}(u,v)=u-u^{\prime}$.

11: end if

12: end for

13: until $\cal S$ is empty

14: return disparity map $D$.

The growing algorithm is summarized in Sec. II-A. The processing of the TOF
correspondence seeds is explained in Sec. II-B, and the sensor fusion based
similarity statistic is described in Sec. II-C.

###  II-A The Growing Procedure

The growing algorithm is sketched in pseudo-code as Alg. 1. The input is a
pair of rectified images $(I_{L},I_{R})$, a set of (refined) TOF seeds $\cal
S$, and a parameter $\tau$ which directly controls a trade-off between
accuracy and density of the matching. The output is a disparity map $D$ which
relates pixel correspondences between the input images.

First, the algorithm computes the prior disparity map $D_{p}$ by interpolating
TOF seeds. Map $D_{p}$ is of the same size as the input images and the output
disparity map, Step 1. Then, a similarity statistic
$\operatorname{simil}(s|I_{L},I_{R},D_{p})$ of the correspondence, which
measures both the photometric consistency of the possible correspondence and
consistency with the prior, is computed for all seeds
$s=(u,u^{\prime},v)\in\cal S$, Step 2. Recall that the seed $s$ stands for a
correspondence $(u,v)\leftrightarrow(u^{\prime},v)$ between pixels in the left
and the right images. For each seed, the algorithm searches other
correspondences in the surroundings of the seeds by maximizing the similarity
statistic. This is done in a 4-neighborhood $\\{{\cal N}_{1},{\cal
N}_{2},{\cal N}_{3}.{\cal N}_{4}\\}$ of the pixel correspondence, such that in
each respective direction (left, right, up, down) the algorithm searches the
disparity in a range $\pm 1$ pixel from the disparity of the seed, Step 6. If
the similarity statistic of a candidate exceeds threshold $\tau$, then a new
correspondence is found, Step 8. It becomes a new seed, and output disparity
map $D$ is updated. The process repeats until there are no more seeds to be
grown.

The algorithm is fairly insensitive to wrong initial seeds. Since the seeds
compete to be matched in the best first strategy, the wrong seeds typically
have low score $\operatorname{simil}(s)$ and therefore when they are drawn in
Step 5, the involved pixels are likely to have been matched already. For more
details on the growing algorithm, we refer to [4, 5].

###  II-B TOF Seeds and Their Refinement

The original version of the seed-growing stereo algorithm [4] uses an initial
set of seeds $\cal S$ obtained by detecting interest points in both images and
matching them. Here, we propose to use TOF seeds. As already outlined, these
seeds are obtained by projecting the low-resolution depth map associated with
the TOF camera onto the high-resolution images. Likewise in the case of
interest points, this yields a sparse set of seeds, e.g., approximately 25,000
seeds in the case of the TOF camera used in our experiments. Nevertheless, one
of the main advantages of the TOF seeds over the interest points is that they
are regularly distributed across the images regardless of the presence/absence
of texture. This is not the case with interest points whose distribution
strongly depends on texture as well as lighting conditions, etc. Obviously,
regularly distributed seeds will provide a better coverage of the observed
scene.

Figure 2: Projection of TOF sensor data on left and right images. The points
are color coded and the color represents disparity such that colder colors are
closer to the cameras. The images are not in the true scale. Notice wrong
correspondences on the computer screen due to low reflectance and artifacts
along occlusion boundaries.  Figure 3: The effect of occlusions. A background
(BG) point $P$ is seen in the left image (IL) and in the TOF image, while it
is occluded by a foreground object (FG) and hence not seen in the right image
(IR). In the process of reprojection of 3D TOF points, a wrong correspondence
($P_{IL}\leftrightarrow P^{\prime}_{IR}$) is produced.

However, TOF seeds are not always accurate. Indeed, when a 3D point is
projected onto the left- and the right-image, it does not always yield a valid
stereo match. There may be several sources of error which make the TOF seeds
less reliable than one would have expected, e.g., Fig. 2 and Fig. 3. In
detail:

  1. 1.

Imprecision due to the calibration data. The transformations allowing to
project the 3D TOF points onto the 2D images are obtained via a complex sensor
calibration process [15]. This introduces a localization error up to 2-3
pixels.

  2. 2.

Outliers due to the physical/geometric properties of the scene. Range sensors
are based on active light and on the assumption that the active beam of light
travels from the sensor and back to it. There are a number of situations where
the beam is lost, such as specular surfaces, absorbing surfaces (such as
fabric), scattering surfaces (such as hair), slanted surfaces, bright surfaces
(computer monitors), faraway surfaces (limited range), or when the beam
travels in an unpredictable way, such a multiple reflections.

  3. 3.

The TOF- and 2D cameras observe the scene from slightly different points of
view. Therefore, it may occur that a 3D point that is present in the TOF image
is only seen into the left or right image, e.g., Fig. 3.

Therefore, a fair percentage of the TOF seeds are outliers. Although the seed-
growing stereo matching algorithm is robust to the presence of outliers in the
initial set of seeds, as already explained in section II-A, we implemented a
straightforward refinement step in order to detect and eliminate these kind of
bad seed data, prior to applying Alg. 1. Firstly, the seeds that lie in low-
intensity (very dark) regions are discarded since TOF-based range data are not
reliable in these cases. Secondly, in order to handle the background-to-
foreground occlusion effect just outlined, we detect seeds which are not
uniformly distributed across image regions. Indeed, projected 3D points lying
on smooth fronto-parallel surfaces form a regular image pattern of seeds,
while projected 3D points that belong to a background surface and which
project onto a foreground image region do not form a regular pattern. See
occlusion boundary in Fig. 4(a).

Such seeds are detected by counting the occupancy of small $5\times 5$ pixel
windows around every seed point in both images. If there is more than one seed
point, the seeds belonging to the background are discarded. Refined set of
seeds is shown in Fig. 4(b). The refinement procedure typically filters 10-15%
of all seed points.

(a) Original set of seeds

(b) Refined set of seeds

Figure 4: An example of the effect of correcting the set of seeds on the basis
that they should be regularly distributed.

(a) Delaunay Triangulation on original seeds

(b) Delaunay Triangulation on refined seeds

(c) Prior obtained on original seeds

(d) Prior obtained on refined seeds

Figure 5: Triangulation and prior disparity map $D_{p}$. These are shown using
both raw seeds (a), (c) and refined seeds (b), (d). A positive impact of the
refinement procedure is clearly visible.

###  II-C Similarity Statistic Based on Sensor Fusion

The original seed-growing matching algorithm [4] uses Moravec’s normalized
cross correalation [17] (MNCC),

$\mbox{simil}(s)=\mbox{MNCC}(w_{L},w_{R})=\dfrac{2\mbox{cov}(w_{L},w_{R})}{\mbox{var}(w_{L})+\mbox{var}(w_{R})+\epsilon}$
(1)

as the similarity statistic to measure the photometric consistency of a
correspondence $s:(u,v)\leftrightarrow(u^{\prime},v)$. We denote by $w_{L}$
and $w_{R}$ the feature vectors which collect image intensities in small
windows of size $n\times n$ pixels centered at $(u,v)$ and $(u^{\prime}v)$ in
the left and right image respectively. The parameter $\epsilon$ prevents
instability of the statistic in cases of low intensity variance. This is set
as the machine floating point epsilon. The statistic has low response in
textureless regions and therefore the growing algorithm does not propagate the
correspondences across these regions. Since the TOF sensor can provide seeds
without the presence of any texture, we propose a novel similarity statistic,
$\operatorname{simil}(s|I_{L},I_{R},D_{p})$. This similarity measure uses a
different score for photometric consistency as well as an initial high-
resolution disparity map $D_{p}$, both incorporated into the Bayesian model
explained in detail below.

The initial disparity map $D_{p}$ is computed as follows. A 3D meshed surface
is built from a 2D triangulation applied to the TOF image. The disparity map
$D_{p}$ is obtained via interpolation from this surface such that it has the
same (high) resolution as of the left and right images. Fig. 5(a) and 5(b)
show the meshed surface projected onto the left high-resolution image and
built from the TOF data, before and after the seed refinement step, which
makes the $D_{p}$ map more precise.

Let us now consider the task of finding an optimal high-resolution disparity
map. For each correspondence $(u,v)\leftrightarrow(u^{\prime},v)$ and
associated disparity $d=u-u^{\prime}$ we seek an optimal disparity $d^{\ast}$
such that:

$d^{\ast}=\underset{d}{\operatorname{argmax}}~{}P(d|I_{L},I_{R},D_{p}).$ (2)

By applying the Bayes’ rule, neglecting constant terms, assuming that the
distribution $P(d)$ is uniform in a local neighborhood where it is sought
(Step. 6), and considering conditional independence
$P(I_{l},I_{r},D|d)=P(I_{L},I_{R}|d)P(D_{p}|d)$, we obtain:

$d^{\ast}=\underset{d}{\operatorname{argmax}}~{}P(I_{L},I_{R}|d)~{}P(D_{p}|d),$
(3)

where the first term is the image likelihood and the second term is the range-
sensor likelihood. We define the image and range-sensor likelihoods as:

$P(I_{L},I_{R}|d)\propto\mbox{EXPSSD}(w_{L},w_{R})=\\\
=\mbox{exp}\left(-\frac{\sum_{i=1}^{n\times
n}(w_{L}(i)-w_{R}(i))^{2}}{\sigma_{s}^{2}\sum_{i=1}^{n\times
n}(w_{L}(i)^{2}+w_{R}(i)^{2})}\right),$ (4)

and as

$P(D_{p}|d)\propto\mbox{exp}\left(-\frac{(d-d_{p})^{2}}{2\sigma_{p}^{2}}\right)$
(5)

respectively, where $\sigma_{s}$ are $\sigma_{p}$ two normalization
parameters. Therefore, the new similarity statistic becomes

$\operatorname{simil}(s|I_{L},I_{R},D_{p})=\mbox{EPC}(w_{L},w_{R},D_{p})=\\\
=\mbox{exp}\left(-\frac{\sum_{i=1}^{n\times
n}(w_{L}(i)-w_{R}(i))^{2}}{\sigma_{s}^{2}\sum_{i=1}^{n\times
n}(w_{L}(i)^{2}+w_{R}(i)^{2})}-\frac{(d-d_{p})^{2}}{2\sigma_{p}^{2}}\right).$
(6)

Notice that the proposed image likelihood has a high response for
correspondences associated with textureless regions. However, in such regions,
all possible matches have similar image likelihoods. The proposed range-sensor
likelihood regularizes the solution and forces it towards the one closest to
the prior disparity map $D_{p}$. A tradeoff between these two terms can be
obtained by tuning the parameters $\sigma_{s}$ and $\sigma_{p}$.

##  III Experiments

Our experimental setup comprises one Mesa Imaging
$\mbox{Swissranger}^{\mbox{\tiny TM}}$ SR4000 TOF camera and a pair of high-
resolution Point Grey555http://www.ptgrey.com/ color cameras, e.g., Fig. 1.
The two color cameras are mounted on a rail with a baseline of about 49 cm and
the TOF camera is approximately midway between them. All three optical axes
are approximately parallel. The resolution of the TOF image is of $144\times
177$ and the color cameras have a resolution of $1224\times 1624$. Recall that
Fig. 1(b) highlights the resolution differences between the TOF and color
images.This camera system was calibrated using the method described in [15].

In all our experiments, we set the parameters of the method as follows:
Windows of $5\times 5$ pixels were used for matching ($n=5$), matching
threshold in Alg. 1 to $\tau=0.5$, the balance between the photometric and
range sensor likelihoods was set to $\sigma_{s}^{2}=0.1$ and to
$\sigma_{p}^{2}=0.001$ in (6).

We show both qualitatively and quantitatively (using datasets with ground-
truth) benefits of the range sensor and an impact of particular variants of
the proposed fusion model integrated in the growing algorithm. Namely, we
compare results of (i) the original stereo algorithm [4] with MNCC correlation
and Harris seeds (MNCC-Harris), (ii) the same algorithm with TOF seeds (MNCC-
TOF), (iii) the algorithm which uses EXPSSD similarity statistic instead with
both Harris (EXPSSD-Harris) and TOF seeds (EXPSSD-TOF), and (iv) the full
sensor fusion model of the regularized growth (EPC). Finally small gaps of
unassigned disparity in the disparity maps were filled by a primitive
procedure which assigns median disparity in the $5\times 5$ window around the
gap (EPC - gaps filled). These small gaps usually occur in slanted surfaces,
since Alg. 1 in Step. 8 enforces one-to-one pixel matching. Nevertheless this
way, they can be filled easily if needed.

###  III-A Real-Data Experiments

We captured two real-world datasets using the camera setup described above,
SET-1 in Fig. 6 and SET-2 in 7. Notice that in both of these examples the
scene surfaces are weakly textured.

(a) Input data RGB-TOF-RGB (the true scale is shown in Fig. 1(b)).

(b) MNCC-Harris

(c) MNCC-TOF

(d) EXPSSD-Harris

(e) EXPSSD-TOF

(f) EPC (proposed)

(g) EPC (gaps filled)

Figure 6:  SET-1. (a) A triplet composed of a pair of color images (left and
right) and a TOF image (middle), results obtained (b) using the seed growing
stereo algorithm [4] combined Harris seeds and MNCC statistic, (c) using TOF
seeds and MNCC statistic , (d) using Harris seeds and EXPSSD statistic, (e)
using TOF seeds with EXPSSD statistics. Results obtained with proposed full
stereo-TOF fusion model using EPC similarity statistic (f) and full model EPC
after filling small gaps (g).

(a) Input data RGB-TOF-RGB (the true scale is shown in Fig. 1(b)).

(b) MNCC-Harris

(c) MNCC-TOF

(d) EXPSSD-Harris

(e) EXPSSD-TOF

(f) EPC (proposed)

(g) EPC (gaps filled)

Figure 7: SET-2. Please refer to the caption of Fig. 6 for explanation.

####  III-A1 Comparisons between disparity maps

Results as disparity maps are shown color-coded, such that warmer colors are
further away from the cameras and unmatched pixels are dark blue.

In Fig. 6(b), we can see that the original algorithm [4] has difficulties in
low textured regions which results in large unmatched regions due to MNCC
statistic (1), and it produces several mismatches over repetitive structures
on the background curtain, due to erroneous (mismatched) Harris seeds. In Fig.
6(c), we can see that after replacing the sparse erratic Harris seeds with
uniformly distributed mostly correct TOF seeds, results improve significantly.
There are no more mismatches on the background, but unmatched regions are
still large. In Fig. 6(d), the EXPSSD statistic (4) was used instead of MNCC
which causes similar mismatches as in Fig. 6(b), but unlike MNCC there are
matches in textureless regions, nevertheless mostly erratic. The reason is
that unlike MNCC statistic the EXPSSD statistic has high response in low
textured regions. However, since all disparity candidates have equal (high)
response inside such regions, the unregularized growth is random, and produces
mismatches. The situation does not improve much using the TOF seeds, as shown
in Fig. 6(e). Significantly better results are finally shown in Fig. 6(f)
which is using the full proposed fusion model EPC (6). The EPC statistic
compared to EXPSSD has the additional regularizing range sensor likelihood
term which guides the growth in ambiguous regions and attracts the solution
towards the rough estimate of the TOF camera. Results are further refined by
filling small gaps in Fig. 6(g). Similar observations can be made in Fig. 7.
The proposed model clearly outperforms the other discussed approaches.

####  III-A2 Comparisons between reconstructed surfaces

For the proper analysis of a stereo matching algorithm it is important to
inspect the reconstructed 3D surfaces. Indeed the visualization of the
disparity/depth maps can sometimes be misleading. Surface reconstruction
reveals fine details in the quality of the results. This is in order to
qualitatively show the gain of the high-resolution depth map produced by the
proposed algorithm with respect to low-resolution depth map of the TOF sensor.

(a) Dense surface reconstruction using the disparity map $D_{p}$ corresponding
to a 2D triangulation of the TOF data points. Zoomed sofa chair and zoomed
T-shirt from SET-2 in Fig. 7(a).

(b) Surface reconstruction using the proposed algorithm (EPC) shown on the
same zoomed areas as above, i.e., Fig. 7(g).

Figure 8: The reconstructed surfaces are shown as relighted 3D meshed for (a)
the prior disparity map $D_{p}$ (2D triangulation on projected and refined TOF
seeds), and (b) for the disparity map obtained using the proposed algorithm.
Notice the fine surface details which were recovered by the proposed method.

In order to provide a fair comparison, we show the reconstructed surfaces
associated with the dense disparity maps $D_{p}$ obtained after 2D
triangulation of the TOF data points, Fig. 8(a), as well as the reconstructed
surfaces associated with the disparity map obtained with the proposed method,
Fig. 8(b). Clearly, much more of the surface details are recovered by the
proposed method. Notice precise object boundaries and fine details, like a
pillow on the sofa chair and a collar of the T-shirt, which appear in Fig.
8(b). This qualitatively corroborates the precision of the proposed method
compared to the TOF data.

###  III-B Ground-Truth Evaluation

To quantitatively demonstrate the validity of the proposed algorithm, we
carried out an experiment on datasets with associated ground-truth results.
Similarly to [8] we use Middlebury dataset [19] and simulated the TOF camera
by sampling the ground-truth disparity map.

We used the Middlebury-2006
dataset666http://vision.middlebury.edu/stereo/data/scenes2006/. On purpose, we
selected three challenging scenes with weakly textured surfaces: Lampshade-1,
Monopoly, Plastic. The input images are of size $1330\times 1110$ pixels. We
took every 10th pixel in a regular grid to simulate the TOF camera. This gives
us about 14k of TOF points, which is roughly the same ratio to color images as
for the real sensors. We are aware that simulation TOF sensor this way is
naive, since we do not simulate any noise or artifacts, but we believe that
for validating the proposed method this is satisfactory.

Results are shown in Fig. 9. We show left input image, results of the same
algorithms as in the previous section with the real sensor, and the ground-
truth disparity map. For each disparity, we compute the percentage of
correctly matched pixels in non-occluded regions. This error statistic is
computed as number of pixels for which the estimated disparity differs from
the ground-truth disparity by less than one pixel divided by number of all
pixels in non-occluded regions. Notice that unmatched pixels are considered as
errors of the same kind as mismatches. This is in order to allow a strict but
fair comparison between algorithms which deliver solution of different
density. The quantitative evaluation confirms the observation from the real-
world setup. The proposed algorithm which uses the full sensor fusion model
significantly outperforms all other tested variants.

For the sake of completeness we also report error statistics for the prior
disparity map $D_{p}$ which is computed by interpolating TOF seeds, see Step 1
of Alg.1. This is $92.9\%,~{}92.1\%,~{}96.0\%$ for Lampshade-1, Monopoly,
Plastic scene respectively. These results are already quite good, which means
the interpolation we use to construct the prior disparity map is appropriate.
These scenes are mostly piecewise planar, which the interpolation captures
well. On the other hand, recall that in the real case, not all the seeds are
correct due to various artifacts of a range sensor. Nevertheless in all three
scenes, the proposed algorithm (EPC with gaps filled) was able to further
improve the precision up to $96.4\%,~{}95.3\%,~{}98.2\%$ for respective
scenes. This is again consistent with the experiments with the real TOF
sensor, where higher surface details were recovered, see Fig. 8.

Left image | MNCC-Harris | MNCC-TOF | EXPSSD-Harris | EXPSSD-TOF | EPC | EPC (gaps filled) | Ground-truth
---|---|---|---|---|---|---|---
|  |  |  |  |  |  |
Lampshade-1 | 61.5% | 64.3% | 44.9% | 49.5% | 88.8% | 96.4% | –
|  |  |  |  |  |  |
Monopoly | 51.2% | 53.4% | 29.4% | 32.1% | 85.2% | 95.3% | –
|  |  |  |  |  |  |
Plastic | 25.2% | 28.2% | 13.5% | 20.6% | 88.7% | 98.2% | –
Figure 9: Middlebury dataset. We show left images, results of the same
algorithms as in Fig. 6 and 7, and the ground-truth disparity maps. There are
error statistics (percentage of correctly matched pixels) below the disparity
maps. Observations from the real-world cases are confirmed quantitatively. The
algorithm using the full model (EPC) gives clearly the best results.

###  III-C Computational Costs

The original growing algorithm [4] has low computational complexity due to
intrinsic search space reduction. Assuming the input stereo images of size
$n\times n$ pixels, the algorithm has the complexity of ${\cal O}(n^{2})$,
while any exhaustive algorithm has the complexity at least ${\cal O}(n^{3})$,
[6]. Factor $n^{3}$ is the size of the search space where the correspondences
are sought, i.e. the disparity space. The growing algorithm does not compute
similarity statistics of all possible correspondences, but efficiently traces
out components of high similarity score around the seeds. This low complexity
is beneficial especially for high resolution imagery, which allows precise
surface reconstruction.

The proposed algorithm with all presented modifications does not represent any
significant extra cost. Triangulation of TOF seeds and the prior disparity map
computation is not really expensive as well as the computing the new EPC
statistic instead of MNCC.

For our experiments, we use an academic (combined Matlab and C) implementation
which takes about 5 seconds with 2 MP images. Nevertheless, recently [10]
presented an implementation of the growing algorithm [4] which runs in real-
time in normal CPU, without parallel hardware. This indicates that a real-time
implementation of the proposed algorithm would be feasible.

##  IV Conclusions

We have proposed a novel correspondence growing algorithm with fusion of a
range sensor and a pair of passive color cameras to obtain accurate and dense
3D reconstruction of a given scene. The proposed algorithm is robust and
performs well on both textured and texture less surfaces and ambiguous
repetitive patterns. The algorithm exploits the strengths of TOF sensor and
stereo matching between color cameras in a combined way to compensate for
their individual weaknesses. The algorithm has shown promising results on
difficult real world data as well as on challenging standard datasets which
quantitatively corroborates its favourable properties. Together with the
strong potential for real-time performance that we discussed, the algorithm
would be practically very useful in many computer vision and robotic
applications.

## References

  * [1] G. Alenyà, B. Dellen, and C. Torras, “3d modelling of leaves from color and tof data for robotized plant measuring,” in _ICRA_ , 2011, pp. 3408–3414.
  * [2] M. Attamimi, A. Mizutani, T. Nakamura, T. Nagai, K. Funakoshi, and M. Nakano, “Real-time 3D visual sensor for robust object recognition,” in _IROS_ , 2010, pp. 4560–4565.
  * [3] V. Castaneda, D. Mateus, and N. Navab, “SLAM combining tof and high-resolution cameras,” _IEEE Workshkop on Motion and Video Computing_ , 2011.
  * [4] J. Čech and R. Šára, “Efficient sampling of disparity space for fast and accurate matching,” in _In Proc. BenCOS Workshop CVPR_ , 2007.
  * [5] J. Cech, J. Matas, and M. Perdoch, “Efficient sequential correspondence selection by cosegmentation,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 32, no. 9, pp. 1568–1581, 2010.
  * [6] J. Cech, J. Sanchez-Riera, and R. Horaud, “Scene flow estimation by growing correspondence seeds,” in _CVPR_ , 2011, pp. 3129–3136.
  * [7] D. Chan, H. Buisman, C. Theobalt, and S. Thrun, “A noise-aware filter for real-time depth upsampling,” in _ECCV Workshop on Multi-camera and Multi-modal Sensor Fusion Algorithms and Applications_ , 2008.
  * [8] C. Dal Mutto, P. Zanuttigh, and G. M. Cortelazzo, “A probabilistic approach to tof and stereo data fusion,” in _3DPVT_ , May 2010.
  * [9] J. Diebel and S. Thrun, “An application of Markov random fields to range sensing,” in _NIPS_ , 2005.
  * [10] M. Dobias and R. Sára, “Real-time global prediction for temporally stable stereo,” in _ICCV Workshop on Live Dense Reconstruction with Moving Cameras_ , 2011, pp. 704–707.
  * [11] D. Droeschel, J. Stückler, D. Holz, and S. Behnke, “Towards joint attention for a domestic service robot - person awareness and gesture recognition using time-of-flight cameras,” in _ICRA_ , 2011, pp. 1205–1210.
  * [12] J. Fischer, G. Arbeiter, and A. Verl, “Combination of time-of-flight depth and stereo using semiglobal optimization,” in _ICRA_ , 2011, pp. 3548–3553.
  * [13] A. Geiger, M. Roser, and R. Urtasun, “Efficient large-scale stereo matching,” in _ACCV_ , 2010, pp. 25–38.
  * [14] S. Gould, P. Baumstarck, M. Quigley, A. Y. Ng, and D. Koller, “Integrating visual and range data for robotic object detection,” in _ECCV Workshop on Multi-camera and Multi-modal Sensor Fusion Algorithms and Applications (M2SFA2)_ , 2008.
  * [15] M. E. Hansard, R. Horaud, M. Amat, and S. Lee, “Projective alignment of range and parallax data,” in _CVPR_ , 2011, pp. 3089–3096.
  * [16] I. Jebari, S. Bazeille, E. Battesti, H. Tekaya, M. Klein, A. Tapus, D. Filliat, C. Meyer, S.-H. S.-H. Ieng, R. Benosman, _et al._ , “Multi-sensor semantic mapping and exploration of indoor environments,” in _IEEE Conference on Technologies for Practical Robot Applications_. IEEE, 2011, pp. 151–156.
  * [17] H. P. Moravec, “Toward automatic visual obstacle avoidance,” in _ICAI_ , 1977, pp. 584–94.
  * [18] R. A. Newcombe and A. J. Davison, “Live dense reconstruction with a single moving camera,” in _CVPR_ , 2010, pp. 1498–1505.
  * [19] D. Scharstein and R. Szeliski, “A taxonomy and evaluation of dense two-frame stereo correspondence algorithms,” in _IJCV_ , vol. 47, no. 1-3, 2002, pp. 7–42.
  * [20] S. Seitz, B. Burless, J. Diebel, D. Scharstein, and S. R., “A comparison and evaluation of multi-view stereo reconstruction algorithms,” in _CVPR_ , vol. 1, 2006, pp. 519–528.
  * [21] J. Stückler and S. Behnke, “Combining depth and color cues for scale- and viewpoint-invariant object segmentation and recognition using random forests,” in _IROS_ , 2010, pp. 4566–4571.
  * [22] L. Wang and R. Yang, “Global stereo matching leveraged by sparse ground control points,” in _CVPR_ , 2011, pp. 3033–3040.
  * [23] J. Zhu, L. Wang, R. Yang, and J. Davis, “Fusion of time-of-flight depth and stereo for high accuracy depth maps,” in _CVPR_ , 2008.


# Free energy calculations along entropic pathways: I. Homogeneous vapor-
liquid nucleation for atomic and molecular systems.

Caroline Desgranges and Jerome Delhommelle  Department of Chemistry,
University of North Dakota, Grand Forks ND 58202

###### Abstract

Using the entropy $S$ as a reaction coordinate, we determine the free energy
barrier associated with the formation of a liquid droplet from a
supersaturated vapor for atomic and molecular fluids. For this purpose, we
develop the $\mu VT-S$ simulation method that combines the advantages of the
grand-canonical ensemble, that allows for a direct evaluation of the entropy,
and of the umbrella sampling method, that is well suited to the study of an
activated process like nucleation. Applying this approach to an atomic system
such as $Ar$ allows us to test the method. The results show that the $\mu
VT-S$ method gives the correct dependence on supersaturation of the height of
the free energy barrier and of the size of the critical droplet, when compared
to predictions from classical nucleation theory and to previous simulation
results. In addition, it provides insight into the relation between entropy
and droplet formation throughout this process. An additional advantage of the
$\mu VT-S$ approach is its direct transferability to molecular systems, since
it uses the entropy of the system as the reaction coordinate. Applications of
the $\mu VT-S$ simulation method to $N_{2}$ and $CO_{2}$ are presented and
discussed in this work, showing the versatility of the $\mu VT-S$ approach.

##  I Introduction

In recent years, computer simulations have considerably advanced our knowledge
on phase transitions and molecular assembly processes by shedding light on the
microscopic mechanisms underlying these phenomena Yasuoka and Matsumoto
(1998); Oxtoby (1992); Shen and Debenedetti (1999); Weakliem and Reiss (1993);
Schenter et al. (1999); Zeng and Oxtoby (1991); Yi et al. (2002); Kinjo et al.
(1999); Toxvaerd (2001); Ford (1996); Talanquer and Oxtoby (1995); Reiss et
al. (1990); Kalikmanov and Van Dongen (1995); Horsch et al. (2008); Neimark
and Vishnyakov (2005); Oxtoby and Evans (1988); Lutsko (2008); Wang et al.
(2008); Ten Wolde et al. (1999); Gonzalez et al. (2015); Loeffler et al.
(2015); Sosso et al. (2016); Xu et al. (2015); Keasler and Siepmann (2015);
Wilhelmsen et al. (2015); van Meel et al. (2015); Hale (1986, 2005); Hale and
Thomason (2010); Yuhara et al. (2015); Lauricella et al. (2015); Singh and
Müller-Plathe (2014); Ni and Dijkstra (2013); Reinhardt and Doye (2014).
Simulations have allowed to identify key parameters for these activated
processes, which allow to measure and analyze the progress of the system
towards the formation of a new phase. These parameters, which are often
complex functions of e.g. the coordinates of the atoms or of some features of
their collective behavior, are known as order parameters or reaction
coordinates (RCs) Ten Wolde and Frenkel (1998); Chen et al. (2001); Oh and
Zeng (1999); Chen et al. (2002); Zhukhovitskii (1995); Nishi et al. (2015);
Lupi et al. (2016); Santiso and Trout (2015); Berryman et al. (2016);
Zimmermann et al. (2015); Lam et al. (2015); Kratzer and Arnold (2015);
Bolhuis and Dellago (2015). RCs are generally tailored to a type of system or
to a type of phenomenon under study. This is the case, for instance, for
vapor-liquid nucleation, a process of great interest for many applications
such as for atmospheric nucleation. Possible choices of RCs are e. g. the
radius of the incipient droplet or its size (in terms of numbers of atoms
within the droplet) Ten Wolde and Frenkel (1998); Tanaka et al. (2005); Kraska
(2006); Oh and Zeng (2000); Senger et al. (1999); Lau et al. (2015); Toxvaerd
(2016), and the analysis of the underlying pathway is generally made through
the determination of the free energy barrier of nucleation and of the critical
size of the nucleus beyond which spontaneous growth occurs.

The entropy of the system $S$ is an especially appealing choice for a RC since
it is, by definition, directly connected to the onset of order (i.e. $S$
decreases as order increases). It also has the additional advantage of being
very general and, as such, can be used for a wide range of systems and
phenomena, regardless of the geometry, structure and phases involved. In the
first part of this series, our goal is to establish the use of $S$ as a RC for
the homogeneous vapor-liquid nucleation of atomic and molecular systems. For
this purpose, we develop a simulation method that allows to vary the entropy
of the system and, as a result, to drive the formation of the liquid droplet.
We then calculate the free energy barrier associated with this process. The
simulations carried out in this work therefore allow us to identify the
entropic pathway underlying the nucleation process, which starts with the
metastable supersaturated vapor and ends at the top of the free energy
barrier, when a liquid droplet of a critical size has formed. By varying the
conditions of nucleation, we are also able to analyze the impact of the
supersaturation on the nucleation process and to rationalize the interplay
between the free energy of nucleation, the size of the liquid droplet and the
range of entropies spanned during the nucleation process. The validity of the
results is assessed through a comparison with the predictions from the
classical nucleation theory McGraw and Laaksonen (1996, 1997), which provides
a relation between the free energy of nucleation and the supersaturation, and
with prior simulation work Ten Wolde and Frenkel (1998); Koga and Zeng (1999).
Results are presented for both atomic ($Ar$) as well as molecular ($N_{2}$ and
$CO_{2}$) systems to demonstrate the versatility and the transferability of
the method to molecular systems.

The paper is organized as follows. In the next section, we present the
simulation method and discuss how the entropy of the system is used as a RC.
We also explain how we calculate the free energy profile for the the
nucleation process. We then present the molecular force fields used to model
the three systems studied in this work, before discussing the technical
details. In particular, we give an account of how we obtain and choose the
input parameters for the simulations to control the extent of supersaturation
during the nucleation process. Then, we discuss the results obtained here on
the three systems considered in this work, $Ar$, $N_{2}$ and $CO_{2}$ and
focus on quantifying the impact of supersaturation on the height of the free
energy barrier, as well as on the size of the liquid droplet and on the
entropy of the system at the top of the free energy barrier. We finally draw
the min conclusions from this work in the last section.

##  II Simulation Method

###  II.1 Sampling entropic pathways

The $\mu VT-S$ approach is based on the grand-canonical ensemble, where $\mu$
is the chemical potential, $V$ the volume, $T$ the temperature and $S$ the
entropy of the system. Having $\mu$ as an input parameter for the simulations
leads to a direct comparison with the theoretical predictions from the
classical nucleation theory (CNT) McGraw and Laaksonen (1996, 1997). CNT
analyzes the nucleation process, and the free energy barrier that the system
has to overcome during this process, as a competition between two
contributions, which cancel out when a droplet of a critical size has formed,
i.e. when the system has reached the top of the free energy barrier. The first
contribution, the surface term, is energetically unfavorable and is the
predominant term for small droplets. It results from the cost of creating the
interface between the liquid droplet and the vapor and can be evaluated as the
product of the area of the incipient spherical droplet by the surface tension.
The second contribution, the volume term, is energetically favorable and
becomes predominant for larger droplets. The volume term corresponds to the
energetic gain in ’converting’ the parent phase, the metastable supersaturated
vapor, into a droplet of the stable phase, the liquid phase. In the CNT
approach McGraw and Laaksonen (1996, 1997), the volume term, and hence the
free energy of nucleation, is proportional to the supersaturation
$\Delta\mu=\mu-\mu_{v}$, which is the difference between the chemical
potential of the liquid ($\mu$) and the chemical potential of the
supersaturated vapor $\mu_{v}$ at the same pressure as the liquid. Carrying
out simulations in the grand-canonical ensemble allows us to choose $\mu$ and
hence the supersaturation $\Delta\mu$ for the nucleation event.

Working in the grand-canonical ensemble also has additional advantages. Since
the number of atoms/molecules is allowed to fluctuate, this ensemble provides
a direct acces to thermodynamic quantities that are difficult to calculate in
other ensembles, as e.g. in the $NVT$ or $NPT$ ensemble. This is the case, for
instance, for the entropy $S$ of the system. In the grand-canonical ensemble,
in which the chemical potential $\mu$, the temperature $T$ and volume $V$ are
fixed during the simulations, it is straightforward to evaluate $S$ during the
simulations through

$S={\bar{U}-\mu\over T}$ (1)

where $\bar{U}$ is the internal energy per atom/molecule, obtained by
calculating the potential energy due to the interactions between
atoms/molecules and by adding the kinetic (ideal gas) contribution to the
internal energy of $k_{B}T/2$ per degree of freedom.

This means that, in the grand-canonical ensemble, $S$ can be used as the
reaction coordinate to measure the onset of order in the system and, in the
case of nucleation, its progress towards the formation of a liquid droplet of
a critical size. To achieve this, we carry out $\mu VT-S$ simulations by
taking advantage of the well-known umbrella sampling method Torrie and Valleau
(1977) to allow the system to follow the RC, in this case $S$. From a
practical standpoint, we add a bias potential energy, function of the entropy
of the system, $U_{bias}(S)$, to the total potential energy of the system. In
this work, we use the following harmonic function of $S$ as the bias potential
energy

$U_{bias}={1\over 2}k(S-S_{0})^{2}$ (2)

where $k$ is a spring constant and $S_{0}$ is the target value for the entropy
of the system.

$\mu VT-S$ simulations are carried out within the Monte Carlo (MC) framework,
with the usual Metropolis criteria used to accept/reject the different types
of MC moves attempted on the atoms/molecules of the system. More specifically,
since these simulations are rooted in the grand-canonical ensemble, the
conventional Metropolis criteria, as obtained e.g. by Allen and Tildesley
Allen and Tildesley (1987) in the grand-canonical ensemble, are used in $\mu
VT-S$ simulations. For instance, for the MC steps involving the translation of
single atom/molecule, as well as for the rotation of a molecule, we use the
following Metropolis criterion

$acc(o\to n)=min\left[1,{\exp(-\beta\Delta U}\right]$ (3)

with $\Delta U$ denoting the change in the total potential energy (interaction
energy+bias potential) corresponding to the move from the old ($o$)
configuration to the new ($n$) configuration. Similarly, for the MC moves
involving the insertion of a new atom/molecule, we use the following
criterion:

$acc(o\to n)=min\left[1,{\exp(-\beta\Delta
U)+\ln\left({N\over{zV}}\right)}\right]$ (4)

where $z$ is the activity $z=\exp(\beta\mu)/\Lambda^{3}$, where $\Lambda$ is
the de Broglie wavelength and $N$ is the current number of atoms/molecules in
the system. In the case of molecular fluids, one also needs to take into
account the terms related to the other degrees of freedom. For linear
molecules like $N_{2}$ and $CO_{2}$, we include for the rotation a factor of
$8\pi Ik_{B}T/2h^{2}$, where $I$ is the moment of inertia of the linear
molecule, $k_{B}$ is the Boltzmann constant and $h$ is Planck’s constant.

For all systems ($Ar$, $CO_{2}$ and $N_{2}$), we carry out $\mu VT-S$
simulations in cubic cells with an edge of $100$ Å and apply the usual $3D$
periodic boundary conditions. To simulate the entire nucleation process, we
carry out a series of umbrella sampling simulations with decreasing values for
the target entropy $S_{0}$. This allows us to sample the underlying entropic
pathway, which goes from the high entropy parent phase (supersaturated vapor)
to the low entropy phase (liquid). During the course of the simulation, we
collect histograms for the number of times each entropy interval is visited.
These histograms are then used to build the free energy profile of nucleation
using the techniques developed to analyze the results from umbrella sampling
simulations (more details can be found in previous work Torrie and Valleau
(1977); Allen and Tildesley (1987); Desgranges and Delhommelle (2007, 2009);
Desgranges and Delhommelle (2014a)). For each value of the target entropy
$S_{0}$, simulations are first run for $100\times 10^{6}$ $MC$ steps to allow
the system to relax. Then, a second run, the production run, is performed for
$200\times 10^{6}$ $MC$ steps during which the data are collected and the
calculation of average properties is carried out. The different types of $MC$
steps are attempted with the following rates. In the case of of $Ar$, $75$ %
of the $MC$ steps consist in the translation of a single atom, while the $25$
% remaining moves are equally split between the $MC$ steps corresponding to
the insertion and deletion of Ar atoms. In terms of computational efficiency,
working in the grand-canonical ensemble allows the number of atoms to vary,
which means that the simulations for the higher $S$ values will be performed
on systems with few atoms. This is advantageous when compared to methods based
on the NPT ensemble, which simulate systems with a constant number of atoms.
Furthermore, calculating the reaction coordinate $S$ through Eq. 1 after every
time step is very fast. The CPU time for a $\mu VT-S$ run is therefore very
close to a conventional grand-canonical run. For instance, a production run
carried out on systems containing an Ar droplet of a critical size (i.e. the
largest systems sizes during nucleation) takes 30 CPU hours (CPU time given
for a single Intel Xeon processor clocked at 3.2 GHz). For $CO_{2}$ and
$N_{2}$, we take into account the $MC$ steps for the rotation of a single
molecule ($37.5$ %), in addition to the translation ($37.5$ %) and the
insertion/deletion ($25$ %). Throughout the simulations, we check that we have
sufficiently large acceptance rates for the insertion/deletion steps to ensure
an efficient sampling. For all systems and all target entropies $S_{0}$, we
obtain acceptance rates greater than $40$ % for the insertion/deletion steps.

###  II.2 Models

We model Argon with a Lennard-Jones (LJ) potential

$\phi(r)=4\epsilon\left[{\left({\sigma\over r}\right)^{12}-\left({\sigma\over
r}\right)^{6}}\right]$ (5)

where $\epsilon=117.05~{}K$ is the depth of the potential well, $\sigma=3.4$ Å
is the atom diameter and $r$ is the distance between two interacting atoms.
The interactions between $Ar$ atoms are calculated for all distances up to
$13.6$ Å and neglected beyond that cutoff distance. We use this large cutoff
distance ($4\sigma$) and do not apply tail corrections, as in previous work on
vapor-liquid nucleation Ten Wolde and Frenkel (1998). We add that an
alternative approach, relying on the calculation of long-range corrections as
a function of the local density, has been developed in recent years for
inhomogeneous systems Ghoufi et al. (2016).

For $N_{2}$ and $CO_{2}$, we consider both molecules to be rigid and model the
interactions with a distribution of point charges to account for the molecular
quadrupole and with a distribution of LJ sites to model the dispersion-
repulson interactions. For $N_{2}$, we consider two LJ sites per $N_{2}$
molecule (one on each atom) and three point charges (one negative charge on
each atom and one positive charge at the center of the $N-N$ bond) to model
the quadrupole of the $N_{2}$ molecule. We use the following parameters
Delhommelle (2000) for the LJ sites, $\epsilon_{NN}/k_{B}=36$ K and
$\sigma_{NN}=3.30$ Å, and for the point charges, $q_{center}=0.966~{}e$ and
$q_{N}=-q_{center}/2$. The $N-N$ bondlength is set to $0.549$ Å, and use a
spherical cutoff set to $4\sigma_{NN}=13.2$ Å for the LJ part of the potential
and to $15$ Å for the quadrupole-quadrupole interactions.

For $CO_{2}$, we use the TraPPE force field Potoff et al. (1999) which is
based on a distribution of three LJ sites (located on the atoms of the
molecule) and of three atomic charges, with the following parameters:
$\epsilon_{CC}/k_{B}=27~{}K$,$\epsilon_{OO}/k_{B}=79~{}K$, $\sigma_{CC}=2.80$
Å and $\sigma_{OO}=3.05$ Å, $q_{c}=0.7~{}e$ and $q_{O}=-0.35~{}e$. The
parameters for the LJ interactions between unlike atoms are given by the
Lorentz-Berthelot mixing rules. The bondlength is set to $1.16$ Å as in the
original paper Potoff et al. (1999). We use a spherical cutoff set to
$4\sigma_{OO}=12.2$ Å for the LJ part of the potential and to $14.6$ Å for the
quadrupole-quadrupole interactions.

###  II.3 Setting up the simulations

####  II.3.1 Argon

$\mu VT-S$ simulations for $Ar$ are carried out at $T=128.76$ K. We present in
Table 1 the conditions for the three supersaturations studied in this work.
They are obtained through the Expanded Wang-Landau (EWL) simulation method we
recently developed Desgranges and Delhommelle (2012a, b, 2014b); Desgranges
and Delhommelle (2016) (other methods can also be used to determine this
information Gazenm$\ddot{\mathrm{u}}$ller and Camp (2007); Liu et al. (2011);
Potoff et al. (1999); Nezbeda and Kolafa (1991); Singh and Errington (2006);
Rai et al. (2007); Rane et al. (2013); Escobedo and de Pablo (1996); Shi and
Maginn (2007); Eslami and Müller-Plathe (2007); Vogt et al. (2001); Widom
(1963); Siepmann and Frenkel (1992)). EWL simulations provide the grand-
canonical partition function and the probability distribution for the number
of $Ar$ atoms in the system for any value of $\mu$. From there, the value of
the chemical potential at the vapor-liquid coexistence can be directly
obtained by finding the value of $\mu$ leading to equal probabilities for the
vapor and liquid phases (more details are provided in previous work Desgranges
and Delhommelle (2012a)). The EWL simulations also give access to all other
thermodynamic properties, including the pressure $P$ and entropies at
coexistence, $S_{l}$ and $S_{v}$, also given in Table 1. This provides an
estimate for the range of entropies that needs to be sampled during the $\mu
VT-S$ simulations to cover the entire vapor$\to$liquid transition. We also
give in Table 1 the supersaturation $\Delta\mu$ that we apply to obtain the
supersaturated vapors, which serve as parent phases for the nucleation
processes. Increasing the value of the chemical potential by $\Delta\mu$
($\Delta\mu>0$) pushes the system further into the domain of the phase diagram
where the liquid phase is stable. We add that $\mu_{v}$ that appears in the
definition of the supersaturation $\Delta\mu$ is very close to $\mu_{coex}$
for the low pressures studied here (of the order of $10^{-3}$ kJ/kg). Here, we
consider 3 different supersaturations and label each set of conditions as
system $1$, $2$ and $3$.

Table 1: $Ar$ at $128.76$ K: Input parameters for the three supersaturations studied in this work (systems 1, 2 and 3). | $\mu$ | $\Delta\mu$ | $P$ | $P/P_{coex}$ | $S_{l}$ | $S_{v}$
---|---|---|---|---|---|---
| $(kJ/kg)$ | $(kJ/kg)$ | $(bar)$ |  | $(kJ/kg/K)$ | $(kJ/kg/K)$
coex | -300.19 | 0.00 | 21.58 | 1.0 | 1.856 | 2.701
system 1 | -298.12 | 2.07 | 43.16 | 2.0 | 1.834 | -
system 2 | -297.72 | 2.47 | 47.48 | 2.2 | 1.831 | -
system 3 | -297.35 | 2.84 | 51.80 | 2.4 | 1.830 | -

The $\mu VT-S$ simulations are initiated from a very low density vapor,
obtained e.g. by placing an $Ar$ atom in the simulation cell. We then carry
out a series of umbrella sampling simulations for decreasing values of the
target entropy $S_{0,i}$, where $i$ is an index labeling the $i^{th}$ umbrella
sampling simulation. A reasonable estimate for the starting value for the
target entropy is provided by $S_{v}$, the value of the vapor at coexistence.
We show in Fig. 1 the histograms $p_{i}(S)$ for the number of times a given
entropy interval is visited during the $\mu VT-S$ simulations for system $2$.
For each umbrella sampling simulation, comparing the relative values of
$S_{max,i}$, the entropy for which the histogram $p_{i}(S)$ reaches its
maximum, to the target entropy $S_{0,i}$, imposed during the $i^{th}$ umbrella
sampling simulation, provides a great deal of insight into the nucleation
process. For instance, at the beginning of the nucleation process, the system
tries to overcome the free energy barrier of nucleation. This means that
increasing the size of the liquid droplet (achieved here by decreasing the
entropy of the system) is associated with the large cost, in terms of free
energy, due to the predominant surface term. As a result, the histogram
$p_{i}(S)$ for the entropy of the system lags behind the target value for the
entropy and that $S_{max,i}>S_{0,i}$. On the other hand, when the system is
past the top of the free energy barrier, the droplet has a larger probability
to keep growing and we have the reverse situation with $S_{max,i}<S_{0,i}$.
The same applies to very low density vapors, when the target value for the
entropy is greater than the entropy of the metastable supersaturated vapor.
Therefore, the entropies $S_{max,i}$ and $S_{0,i}$ only coincide when the
system reaches a point corresponding to an extremum of the free energy
profile, This extremum can be either a maximum (i.e. the top of the free
energy barrier, when a liquid droplet of a critical size has formed) or a
minimum (i.e. the metastable supersaturated vapor or the stable liquid).

Figure 1: $Ar$ at $T=128.76~{}K$. Entropy histograms $p_{i}(S)$ for successive
sampling windows for system $2$. Figure 2: $Ar$ at $T=128.76~{}K$. Potential
energy (top) and number of atoms (bottom) for system $2$ and target entropies
of $2.43$ kJ/kg/K (black) and $2.38$ kJ/kg/K (red).

Looking at the first umbrella sampling window, starting from the right of Fig.
1, we find that $S_{max,1}$ ($2.494$ kJ/kg/K) is lower than $S_{0,1}$ ($2.5$
kJ/kg/K). This shows that, at this stage, the entropy of the system is greater
than the entropy of the metastable supersaturated vapor, or, in other words,
that the density of the system has not yet reached the density of the
supersaturated vapor. Decreasing the target entropy allows us to identify the
first value of the entropy for which $S_{max,i}=S_{0,i}$, where we have the
metastable supersaturated vapor for an entropy of $2.49$ kJ/kg/K. Then, during
the next umbrella sampling simulations for decreasing $S_{0,i}$, the system
climbs up the free energy barrier of nucleation and we have
$S_{max,i}<S_{0,i}$, until we observe again the coincidence of $S_{max,i}$
with $S_{0,i}$ for an entropy of $2.3$ kJ/kg/K. This occurs when the system
reaches the top of the free energy barrier and a liquid droplet of a critical
size has formed. These results show that the $\mu VT-S$ method allows us to
explore the entropic pathway underlying the nucleation process and identify
the extrema of the free energy profile. We also show in Fig. 2 the impact of
changing the target entropy on the properties of the system. Comparing the
results for target entropies of $2.43$ kJ/kg/K and $2.38$ kJ/kg/K shows a
decrease in the potential energy of the system by close to $30$ % and an
increase in the number of $Ar$ atoms in the system by about $10$ %. Both
results can be accounted for by the increase in the density of the system,
triggered by the decrease in target entropy, which leads to an increased
number of attractive interactions between $Ar$ atoms. This shows that the
decrease in entropy achieved through the $\mu VT-S$ simulations results in a
transition towards the liquid phase. We leave the detailed discussion of the
energetics and characteristics of the liquid droplet to the ’Results’ section.

####  II.3.2 Nitrogen

The formation of a liquid droplet of $N_{2}$ from a supersaturated vapor is
studied at $T=100~{}K$. As for $Ar$, we use the EWL simulation method to
obtain the properties at coexistence (coex) and at three supersaturations
(systems 4, 5 and 6), which are increasingly deeper inside the liquid domain
of the phase diagram of $N_{2}$. We list in Table 2 the data for these points,
from which we simulate the nucleation process.

Table 2: $N_{2}$ at $100$ K: Conditions for vapor-liquid coexistence (coex) and for three supersaturations (systems 4, 5 and 6). | $\mu$ | $\Delta\mu$ | $P$ | $P/P_{coex}$ | $S_{l}$ | $S_{v}$
---|---|---|---|---|---|---
| $(kJ/kg)$ | $(kJ/kg)$ | $(bar)$ |  | $(kJ/kg/K)$ | $(kJ/kg/K)$
coex | -403.97 | 0.00 | 10.12 | 1.0 | 3.374 | 4.912
system 4 | -402.16 | 1.81 | 21.84 | 2.2 | 3.352 | -
system 5 | -401.86 | 2.11 | 23.85 | 2.4 | 3.349 | -
system 6 | -401.57 | 2.40 | 25.86 | 2.6 | 3.345 | -

####  II.3.3 Carbon dioxide

In line with the two previous systems, we start from the conditions for vapor-
liquid coexistence (coex) provided by the EWL simulations. We then define
three supersaturations, for increasing $\Delta\mu$, leading to systems 7, 8
and 9. The thermodynamic data at coexistence and for the three
supersaturations are given in Table 3.

Table 3: $CO_{2}$ at $260$ K: Thermodynamic data at the vapor-liquid coexistence (coex) and for 3 supersaturations (systems 7, 8 and 9). | $\mu$ | $\Delta\mu$ | $P$ | $P/P_{coex}$ | $S_{l}$ | $S_{v}$
---|---|---|---|---|---|---
| $(kJ/kg)$ | $(kJ/kg)$ | $(bar)$ |  | $(kJ/kg/K)$ | $(kJ/kg/K)$
coex | -894.03 | 0.00 | 30.22 | 1.0 | 2.962 | 3.946
system 7 | -890.84 | 3.19 | 60.28 | 2.0 | 2.944 | -
system 8 | -890.20 | 3.83 | 66.53 | 2.2 | 2.940 | -
system 9 | -889.61 | 4.42 | 72.31 | 2.4 | 2.937 | -

##  III Results and Discussion

###  III.1 Argon

Figure 3: $Ar$ at $T=128.76$ K: Free energy profiles of nucleation for system
1 (black), system 2 (red) and system 3 (green).

We start by discussing the results obtained for the free energy profiles of
nucleation for Argon at $T=128.76$ K. We show in Fig. 3 the free energy
barriers obtained from the $\mu VT-S$ simulations for systems 1, 2 and 3
against the entropy of the system, which serves as the reaction coordinate for
the nucleation process. Two main trends appear in this plot. First, the height
of the free energy barrier increases as the supersaturation $\Delta\mu$
decreases. When $\Delta\mu=2.84$ kJ/kg (system 3), we obtain a free energy
barrier of nucleation of $26\pm 2~{}k_{B}T$, while for $\Delta\mu=2.47$ kJ/kg
(system 2), we have a free energy barrier of $41\pm 3~{}k_{B}T$. Moreover, as
the supersaturation further decreases to $\Delta\mu=2.04$ kJ/kg (system 1),
the free energy barrier becomes $62\pm 5~{}k_{B}T$. This behavior is
consistent with the predictions from classical nucleation theory and with the
findings from previous simulation work. Changing the supersaturation
$\Delta\mu$ has a direct impact on the volume term, which is the gain, in free
energy, resulting from the formation of the thermodynamically stable liquid
phase from the metastable vapor phase. This volume term is proportional to
$\Delta\mu$ and it therefore follows that for a small $\Delta\mu$ (system 1),
the top of the free energy barrier is reached later, for a larger liquid
droplet, leading to a higher free energy barrier. The values obtained in this
work for the free energy barriers of nucleation are also in good agreement
with those predicted by the classical nucleation theory McGraw and Laaksonen
(1996, 1997) and with those found for the Lennard-Jones system Ten Wolde and
Frenkel (1998). For instance, CNT leads to the following free energy of
nucleation McGraw and Laaksonen (1996, 1997)
$W=16\pi\gamma^{3}/(3\rho_{l}^{2}\Delta\mu^{2})$, where $\gamma$ is the
surface tension for a flat interface and $\rho_{l}$ is the density of the
(bulk) liquid. Applying this to e.g. System 1 with $\gamma=3.535\times
10^{-10}kJ/cm^{2}$ (determined from the EWL simulations Desgranges and
Delhommelle (2012a)) and $\rho_{l}=1.0882~{}g/cm^{3}$ leads to a free energy
barrier of $70~{}k_{B}T$, in reasonably good agreement with the barrier found
here of $62\pm 5~{}k_{B}T$. Our results are also consistent with ten Wolde and
Frenkel Ten Wolde and Frenkel (1998) who reported free energy barriers of
nucleation $5-15~{}k_{B}T$ lower than the CNT predictions for the Lennard-
Jones system. The second main trend is the following. The entropic pathway
spanned during nucleation covers a range of entropy that depends on the
supersaturation. More specifically, for a large supersaturation, the entropy
range becomes narrower. This is due to the fact that for a low supersaturation
(i.e. not too deep into the liquid domain of the phase diagram), the
metastable supersaturated vapor is less dense. This means that the starting
point for the nucleation process is associated with a larger entropy.
Furthermore, at low supersaturation, the size of the critical nucleus is
larger. This implies that, when the supersaturation is low, the entropy of the
system at the top of the free energy barrier is smaller. Both factors account
for the extended range of entropy spanned during the nucleation process at low
supersaturation. We finally add that the range of entropy spanned during
nucleation is, for all supersaturations, well within the interval defined by
the entropies of the vapor ($2.701$ kJ/kg/K) and liquid ($1.856$ kJ/kg/K) at
coexistence. This is the expected behavior, since the metastable
supersaturated vapors are more dense, with a smaller entropy, than the vapor
at coexistence, and since systems containing a critical liquid droplet are
considerably less dense, with a larger entropy, than a uniform liquid phase
like the liquid at coexistence.

(a) (b)

Figure 4: $Ar$ at $T=128.76$ K: (a) Distribution for the number of neighbors
in the vapor (dashed line) and in the liquid (solid line) and (b)
Correspondence between the entropy and the number of liquid-like atoms along
the nucleation (system 2).

(a) (b) (c)

Figure 5: $Ar$ at $T=128.76$ K: Snapshots for system 1 during the nucleation
process for decreasing values of the entropy: (a) $S=2.43$ kJ/kg/K, (b)
$S=2.38$ kJ/kg/K, (c) $S=2.31$ kJ/kg/K

The plots in Fig. 3 also allow us to characterize the state of the system at
the top of the free energy barrier of nucleation in terms of the critical
entropy associated with the formation of a critical nucleus. This critical
entropy $S_{c}$ is shown to depend on supersaturation, with the largest value
for $S_{c}$ being obtained for the lowest supersaturation $S_{c}=2.21\pm 0.02$
kJ/kg (system 1) and the smallest value $S_{c}=2.35\pm 0.01$ kJ/kg for the
highest supersaturation (system 3). To analyze further the characteristics of
the liquid droplet, we examine the relation between the entropy and the size
of the liquid droplet. Starting from the approach developed by ten Wolde and
Frenkel Ten Wolde and Frenkel (1998), we first determine the distribution for
the number of neighbors in the vapor and in the liquid. Here, neighboring
particles are defined as particles that are less than $1.6\sigma$ apart. This
distance of $1.6\sigma$ corresponds to the first minimum of the radial pair
distribution function of the liquid. Fig 4(a) shows the distribution so
obtained for the vapor and for the liquid phases. The two distributions show
little overlap, with the vast majority of the liquid atoms having at least $6$
neighbors, while atoms in the vapor phase very rarely have more than $5$
neighbors. This allows us to define liquid-like atoms as having at least $6$
neighbors and vapor-like atoms having fewer than $6$ neighbors and to follow
the variations of the number of liquid-like atoms during nucleation. Applying
this analysis during the $\mu VT-S$ simulations leads us to draw a
correspondence between the value of the entropy and the number of liquid-like
atoms of the system as nucleation proceeds. We present in Fig. 4(b) the result
for this correspondence for system 2. The number of liquid-like atoms
increases steadily with entropy as the nucleation process advances, implying
that the size of the liquid droplet steadily increases as the entropy of the
system decreases. This can best be seen by examining snapshots of the
configurations of the system during the nucleation process. Fig. 5 shows the
formation and development of the liquid droplet along the entropic pathway. As
can be seen on the snapshots, the droplet size steadily increases with entropy
and reaches a critical size of $1715\pm 50$ atoms for system 1. The snapshots
for the other 2 supersaturations reveal a similar behavior, and critical sizes
of $980\pm 25$ and $649\pm 15$ atoms are found for the two other
supersaturations systems 2 and 3, respectively. We obtain the expected
correlation between the height of the free energy barrier of nucleation and
the size of the critical nucleus, as the size of the critical droplet
decreases with the free energy of nucleation as we go from system 1 to system
2 and finally to system 3.

###  III.2 Nitrogen

Figure 6: $N_{2}$ at $T=100$ K: Free energy profiles of nucleation along $S$
for system 4 (black), system 5 (red) and system 6 (green).

We now turn to the assessment of the versatility of the method and test it on
a molecular fluid. We apply the $\mu VT-S$ approach to study the nucleation of
a liquid droplet of $N_{2}$ at $100$ K. We follow the same procedure as for
$Ar$ and perform a series of umbrella sampling simulations spanning the
entropic pathway underlying the formation of the droplet. We present in Fig. 6
the free energy profiles of nucleation of $N_{2}$ at $T=100$ for three
supersaturations (systems 4, 5 and 6). The results show features that are
qualitatively similar to those found for $Ar$. In terms of entropy range, the
entropies sampled during the process lie within the wider interval defined by
the entropy of the liquid and of the vapor at coexistence (from
$S_{liq}=3.374$ kJ/kg/K to $S_{vap}=4.912$ kJ/kg/K). For instance, for the
lowest supersaturation (system 4), the nucleation starts from a supersaturated
vapor with $S=4.6$ kJ/kg/K (lower than the entropy of the coexisting vapor,
that has a lower density) and ends with configurations containing a critical
nucleus at $S=4.08$ kJ/kg/K (higher than the entropy of the liquid at
coexistence, which has a much a higher density). The same reasoning applied to
systems 5 and 6, for which the nucleation process is associated to an even
narrower entropy range due to the higher supersaturations involved for these
systems. In particular, the entropy $S_{c}$ for which a liquid droplet of a
critical size has formed is shown to increase as the supersaturation
decreases, starting from $S_{c}=4.08\pm 0.02$ kJ/kg/k (system 4) and
increasing to $S_{c}=4.16\pm 0.02$ kJ/kg/K (system 5) and finally to
$S_{c}=4.20\pm 0.02$ kJ/kg/K (system 6) for the highest supersaturation. The
free energy barriers of nucleation for $N_{2}$ are also shown to decrease
gradually as the supersaturation is increased. The free energy barrier of
nucleation reaches $73\pm 5~{}k_{B}T$ for system 4 (lowest supersaturation),
$60\pm 4~{}k_{B}T$ for system 5 and $52\pm 4~{}k_{B}T$ for system 6 (highest
supersaturation).

(a) (b)

Figure 7: $N_{2}$ at $T=100$ K: (a) Distribution for the number of neighbors
in the vapor (dashed line) and in the liquid (solid line) and (b)
Correspondence between the entropy and the number of liquid-like molecules
during the nucleation process (system 6).

(a) (b) (c)

Figure 8: $N_{2}$ at $T=100$ K: Snapshots of system 6 during the nucleation
process for decreasing values of the entropy: (a) $S=4.5$ kJ/kg/K, (b)
$S=4.36$ kJ/kg/K, (c) $S=4.2$ kJ/kg/K

We further refine our analysis by determining the size of the droplet during
the nucleation process. For this purpose, we adopt a similar procedure to that
followed for Argon. We compute the distributions for the number of neighbors
for $N_{2}$ molecules in the supersaturated vapor and in the liquid, using the
first minimum of the pair distribution function (based on the distance between
the centers-of-mass of 2 $N_{2}$ molecules) at $5.7$ Å as a cutoff distance
used to define neighbors. The results plotted in Fig. 7(a) show that these
distributions allow us to identify vapor-like molecules (with less than $6$
neighbors) and liquid-like molecules (with $6$ or more neighbors). The
variation of the number of liquid-like molecules in the system as a function
of entropy during the nucleation process is shown in Fig. 7(b). As with
$Ar$,the number of liquid-like molecules exhibits a steady increase with
entropy, corresponding to the increase in the size of the droplet as
nucleation proceeds. An examination of the snapshots (see Fig. 8) obtained
during the $\mu VT-S$ simulations sheds light on the formation of the liquid
droplet as the entropy of the system decreases. This allows us to identify the
size of the critical droplet, which contains $1006\pm 30$ $N_{2}$ molecules
for system 4, $735\pm 25$ molecules for system 5 and $638\pm 20$ molecules for
system 6. We observe the expected correlation between the size of the critical
droplet and the height of the free energy barrier of nucleation, with the
lowest supersaturation of system 4 leading to the largest free energy of
nucleation and to the biggest critical droplet, that are reached for the
lowest critical entropy. This set of results show that $S$ can be used as the
reaction coordinate for the nucleation process in a molecular fluid, and that
no adjustment or redefinition of the RC is necessary.

###  III.3 Carbon dioxide

Figure 9: $CO_{2}$ at $T=260$ K: Free energy profiles of nucleation along $S$
for system 7 (black), system 8 (red) and system 9 (green).

We now apply the $\mu VT-S$ method to a second molecular fluid and study the
nucleation of a liquid droplet in $CO_{2}$ at $260~{}K$. $CO_{2}$ is similar
to $N_{2}$ in many respects, having a similar shape and exhibiting similar
types of intermolecular interactions (with a quadrupole moment for $CO_{2}$
more than 6 times greater than for $N_{2}$). Carrying out the $\mu VT-S$
simulations for the formation of the liquid droplet at $260$ K allows us to
draw a comparison with the results on $N_{2}$. This is because the conditions
of temperature for $N_{2}$ and $CO_{2}$ can be considered as similar, if one
considers reduced temperatures, with respect to the critical temperature, for
each compounds ($0.85$ for $CO_{2}$ and $0.83$ for $N_{2}$). As for $Ar$ and
$N_{2}$, we perform a series of $\mu VT-S$ simulations to sample the entropic
pathway leading towards the formation of a liquid droplet of $CO_{2}$ of a
critical size for 3 supersaturations. We show in Fig. 9 the free energy
profiles of nucleation so obtained. Similarly to the nucleation of liquid
droplet in $N_{2}$, the profiles show that decreasing the entropy of the
system during the $\mu VT-S$ simulations allows the system to overcome the
free energy barrier and that $S$ can be reliably used to drive the nucleation
process in the molecular fluid of $CO_{2}$.

(a) (b)

Figure 10: $CO_{2}$ at $T=260$ K: (a) Distribution for the number of neighbors
in the vapor (dashed line) and in the liquid (solid line) and (b)
Correspondence between the entropy and the number of liquid-like molecules
during nucleation (system 7).

The range of entropies sampled along the entropic pathway is shown to depend
on the extent of supersaturation. As for $Ar$ and $N_{2}$, the entropy for the
starting point, the metastable supersaturated vapor, decreases at high
supersaturations as the supersaturated vapor becomes more dense. The entropy
for which the top of the free energy barrier is reached also decreases as the
supersaturation becomes high, as a result of the smaller size of the critical
droplet. Both factors account for the reduced range of entropies sampled for
system 9 (high supersaturation) when compared to system 7 (low
supersaturation). This leads to the following values for the entropies $S_{c}$
corresponding to systems containing a liquid droplet of a critical size, with
$S_{c}=3.55\pm 0.02$ kJ/kg/k (system 7) and decreasing to $S_{c}=3.59\pm 0.02$
kJ/kg/K (system 8) and finally to $S_{c}=3.61\pm 0.02$ kJ/kg/K (system 9) for
the highest supersaturation. This increase in $S_{c}$ at high supersaturations
is correlated with a decrease in the height of the free energy barrier of
nucleation which reaches $26\pm 3~{}k_{B}T$ for system 7 (lowest
supersaturation), $21\pm 4~{}k_{B}T$ for system 8 and $14\pm 2~{}k_{B}T$ for
system 9 (highest supersaturation). Since the reduced temperatures and the
supersaturations for the nucleations of $N_{2}$ and $CO_{2}$ are similar, we
attribute the much smaller free energy barriers of nucleation, when compared
to those found for $N_{2}$, to the stronger intermolecular interactions in
$CO_{2}$.

We now analyze the nucleation process in terms of the number of liquid-like
molecules in the system. Closely following the process used to define this
order parameter for the other systems, we determine the distributions for the
number of neighbors in the supersaturated vapor and in the liquid using a
cutoff distance of $5.85$ Å to characterize neighboring molecules. The
distributions plotted in Fig. 10(a) allow us to distinguish between vapor-like
molecules ($6$ neighbors or less) and liquid-like molecules ($7$ neighbors or
more). Using this criterion, we determine the number of liquid-like molecules
during the nucleation process as the entropy changes along the nucleation
pathway. The results are plotted in Fig. 10(b) for system 7 and show a smooth
increase in the number of liquid-like molecules as a function of the entropy
during the formation of the liquid droplet. The formation of the liquid
droplet can be best captured by looking at the snapshots of Fig. 11, which
show the system for ddecreasing values of the entropy during the nucleation
process. These snapshots show the increase in the size of the nucleus and
reveal that the critical droplet contains $395\pm 33$ molecules for system 7.
The critical sizes obtained for the other systems are the following: $290\pm
28$ molecules for system 8 and $242\pm 20$ molecules for system 9. As for $Ar$
and $N_{2}$, the critical sizes are correlated with the free energy barriers
of nucleation, with the highest barrier of system 7 corresponding to the
biggest critical droplet and the lowest supersaturation.

(a) (b) (c)

Figure 11: $CO_{2}$ at $T=260$ K: Snapshots of the system during the
nucleation process for decreasing values of the entropy: (a) $S=3.66$ kJ/kg/K,
(b) $S=3.6$ kJ/kg/K, (c) $S=3.54$ kJ/kg/K

##  IV Conclusion

In this work, we develop the $\mu VT-S$ simulation method to calculate the
free energy of nucleation of atomic and molecular fluids along entropic
pathways. We achieve this by working in the grand-canonical ensemble, which
allows for the direct evaluation of entropy during the simulations and for the
study of the nucleation process at a given supersaturation $\Delta\mu$,
thereby providing a direct connection with the predictions from the classical
nucleation theory. The simulation protocol for $\mu VT-S$ consists of
performing a series of umbrella sampling simulations that use $S$ as the
reaction coordinate for the nucleation process. The results presented here
show that the method provides a picture of the nucleation process that is
consistent with the classical nucleation theory and with that found in
previous work, including e.g. the dependence of the height of the free energy
barrier and of the size of the critical nucleus on supersaturation. They also
establish that the $\mu VT-S$ approach can be applied to atomic as well as
molecular systems without any modification. Furthermore, our simulations
reveal the correlation between the entropy of the system and the droplet size,
and show how supersaturation impacts the entropic pathway visited during the
nucleation process, with a narrowing of the range of entropy visited at high
supersaturation. Our findings also allow us to identify the critical value for
the entropy associated with the formation of a nucleus of a critical size.
Further work presented in the next two papers of the series will focus on the
development of the method for multi-component systems and for heterogeneous
nucleation processes.

Acknowledgements Partial funding for this research was provided by NSF through
CAREER award DMR-1052808.

## References

  * Yasuoka and Matsumoto (1998) K. Yasuoka and M. Matsumoto, J. Chem. Phys. 109, 8451 (1998).
  * Oxtoby (1992) D. W. Oxtoby, J. Phys. Condens. Matt. 4, 7627 (1992).
  * Shen and Debenedetti (1999) V. K. Shen and P. G. Debenedetti, J. Chem. Phys. 111, 3581 (1999).
  * Weakliem and Reiss (1993) C. L. Weakliem and H. Reiss, J. Chem. Phys. 99, 5374 (1993).
  * Schenter et al. (1999) G. K. Schenter, S. M. Kathmann, and B. C. Garrett, Phys. Rev. Lett. 82, 3484 (1999).
  * Zeng and Oxtoby (1991) X. C. Zeng and D. W. Oxtoby, J. Chem. Phys. 94, 4472 (1991).
  * Yi et al. (2002) P. Yi, D. Poulikakos, J. Walther, and G. Yadigaroglu, Int. J. Heat Mass Transfer 45, 2087 (2002).
  * Kinjo et al. (1999) T. Kinjo, K. Ohguchi, K. Yasuoka, and M. Matsumoto, Comput. Mat. Sci. 14, 138 (1999).
  * Toxvaerd (2001) S. Toxvaerd, J. Chem. Phys. 115, 8913 (2001).
  * Ford (1996) I. Ford, J. Chem. Phys. 105, 8324 (1996).
  * Talanquer and Oxtoby (1995) V. Talanquer and D. Oxtoby, J. Phys. Chem. 99, 2865 (1995).
  * Reiss et al. (1990) H. Reiss, A. Tabazadeh, and J. Talbot, J. Chem. Phys. 92, 1266 (1990).
  * Kalikmanov and Van Dongen (1995) V. Kalikmanov and M. Van Dongen, J. Chem. Phys. 103, 4250 (1995).
  * Horsch et al. (2008) M. Horsch, J. Vrabec, and H. Hasse, Phys. Rev. E 78, 011603 (2008).
  * Neimark and Vishnyakov (2005) A. V. Neimark and A. Vishnyakov, J. Chem. Phys. 122, 054707 (2005).
  * Oxtoby and Evans (1988) D. W. Oxtoby and R. Evans, J. Chem. Phys. 89, 7521 (1988).
  * Lutsko (2008) J. F. Lutsko, J. Chem. Phys. 129, 244501 (2008).
  * Wang et al. (2008) Z.-J. Wang, C. Valeriani, and D. Frenkel, J. Phys. Chem. B 113, 3776 (2008).
  * Ten Wolde et al. (1999) P. R. Ten Wolde, M. J. Ruiz-Montero, and D. Frenkel, J. Chem. Phys. 110, 1591 (1999).
  * Gonzalez et al. (2015) M. A. Gonzalez, J. L. Abascal, C. Valeriani, and F. Bresme, J. Chem. Phys. 142, 154903 (2015).
  * Loeffler et al. (2015) T. D. Loeffler, A. Sepehri, and B. Chen, J. Chem. Theory Comput. 11, 4023 (2015).
  * Sosso et al. (2016) G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, and A. Michaelides, Chem. Rev. 116, 7078 (2016).
  * Xu et al. (2015) W. Xu, Z. Lan, B. Peng, R. Wen, and X. Ma, J. Chem. Phys. 142, 054701 (2015).
  * Keasler and Siepmann (2015) S. J. Keasler and J. I. Siepmann, J. Chem. Phys. 143, 164516 (2015).
  * Wilhelmsen et al. (2015) Ø. Wilhelmsen, T. T. Trinh, S. Kjelstrup, and D. Bedeaux, J. Phys. Chem. C 119, 8160 (2015).
  * van Meel et al. (2015) J. van Meel, Y. Liu, and D. Frenkel, Mol. Phys. 113, 2742 (2015).
  * Hale (1986) B. N. Hale, Phys. Rev. A 33, 4156 (1986).
  * Hale (2005) B. N. Hale, J. Chem. Phys. 122, 204509 (2005).
  * Hale and Thomason (2010) B. N. Hale and M. Thomason, Phys. Rev. Lett. 105, 046101 (2010).
  * Yuhara et al. (2015) D. Yuhara, B. C. Barnes, D. Suh, B. C. Knott, G. T. Beckham, K. Yasuoka, D. T. Wu, and A. K. Sum, Faraday Discuss. 179, 463 (2015).
  * Lauricella et al. (2015) M. Lauricella, S. Meloni, S. Liang, N. J. English, P. G. Kusalik, and G. Ciccotti, J. Chem. Phys. 142, 244503 (2015).
  * Singh and Müller-Plathe (2014) J. K. Singh and F. Müller-Plathe, Appl. Phys. Lett. 104, 021603 (2014).
  * Ni and Dijkstra (2013) R. Ni and M. Dijkstra, Soft Matter 9, 365 (2013).
  * Reinhardt and Doye (2014) A. Reinhardt and J. P. Doye, J. Chem. Phys. 141, 084501 (2014).
  * Ten Wolde and Frenkel (1998) P. R. Ten Wolde and D. Frenkel, J. Chem. Phys. 109, 9901 (1998).
  * Chen et al. (2001) B. Chen, J. I. Siepmann, K. J. Oh, and M. L. Klein, J. Chem. Phys. 115, 10903 (2001).
  * Oh and Zeng (1999) K. Oh and X. C. Zeng, J. Chem. Phys. 110, 4471 (1999).
  * Chen et al. (2002) B. Chen, J. I. Siepmann, K. J. Oh, and M. L. Klein, J. Chem. Phys. 116, 4317 (2002).
  * Zhukhovitskii (1995) D. Zhukhovitskii, J. Chem. Phys. 103, 9401 (1995).
  * Nishi et al. (2015) K. Nishi, S. Inoue, and Y. Matsumura, Chem. Phys. Lett. 634, 194 (2015).
  * Lupi et al. (2016) L. Lupi, B. Peters, and V. Molinero, J. Chem. Phys. 145, 211910 (2016).
  * Santiso and Trout (2015) E. E. Santiso and B. L. Trout, J. Chem. Phys. 143, 174109 (2015).
  * Berryman et al. (2016) J. T. Berryman, M. Anwar, S. Dorosz, and T. Schilling, J. Chem. Phys. 145, 211901 (2016).
  * Zimmermann et al. (2015) N. E. Zimmermann, B. Vorselaars, D. Quigley, and B. Peters, J. Am. Chem. Soc. 137, 13352 (2015).
  * Lam et al. (2015) J. Lam, D. Amans, C. Dujardin, G. Ledoux, and A.-R. Allouche, J. Phys. Chem. A 119, 8944 (2015).
  * Kratzer and Arnold (2015) K. Kratzer and A. Arnold, Soft Matter 11, 2174 (2015).
  * Bolhuis and Dellago (2015) P. Bolhuis and C. Dellago, Eur. Phys. J. Special Topics 224, 2409 (2015).
  * Tanaka et al. (2005) K. K. Tanaka, K. Kawamura, H. Tanaka, and K. Nakazawa, J. Chem. Phys. 122, 184514 (2005).
  * Kraska (2006) T. Kraska, J. Chem. Phys. 124, 054507 (2006).
  * Oh and Zeng (2000) K. Oh and X. C. Zeng, J. Chem. Phys. 112, 294 (2000).
  * Senger et al. (1999) B. Senger, P. Schaaf, D. Corti, R. Bowles, D. Pointu, J.-C. Voegel, and H. Reiss, J. Chem. Phys. 110, 6438 (1999).
  * Lau et al. (2015) G. V. Lau, P. A. Hunt, E. A. Müller, G. Jackson, and I. J. Ford, J. Chem. Phys. 143, 244709 (2015).
  * Toxvaerd (2016) S. Toxvaerd, J. Chem. Phys. 144, 164502 (2016).
  * McGraw and Laaksonen (1996) R. McGraw and A. Laaksonen, Physical Rev. Lett. 76, 2754 (1996).
  * McGraw and Laaksonen (1997) R. McGraw and A. Laaksonen, J. Chem. Phys. 106, 5284 (1997).
  * Koga and Zeng (1999) K. Koga and X. C. Zeng, J. Chem. Phys. 110, 3466 (1999).
  * Torrie and Valleau (1977) G. M. Torrie and J. P. Valleau, J. Comput. Phys. 23, 187 (1977).
  * Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley, _Computer Simulation of Liquids_ (Clarendon Press, Oxford, 1987).
  * Desgranges and Delhommelle (2007) C. Desgranges and J. Delhommelle, J. Chem. Phys. 127, 144509 (2007).
  * Desgranges and Delhommelle (2009) C. Desgranges and J. Delhommelle, J. Phys. Chem. C 113, 3607 (2009).
  * Desgranges and Delhommelle (2014a) C. Desgranges and J. Delhommelle, J. Am. Chem. Soc. 136, 8145 (2014a).
  * Ghoufi et al. (2016) A. Ghoufi, P. Malfreyt, and D. J. Tildesley, Chem. Soc. Rev. 45, 1387 (2016).
  * Delhommelle (2000) J. Delhommelle, Ph.D. thesis, Universite Paris XI-Orsay (2000).
  * Potoff et al. (1999) J. J. Potoff, J. R. Errington, and A. Z. Panagiotopoulos, Mol. Phys. 97, 1073 (1999).
  * Desgranges and Delhommelle (2012a) C. Desgranges and J. Delhommelle, J. Chem. Phys. 136, 184107 (2012a).
  * Desgranges and Delhommelle (2012b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 136, 184108 (2012b).
  * Desgranges and Delhommelle (2014b) C. Desgranges and J. Delhommelle, J. Chem. Phys. 140, 104109 (2014b).
  * Desgranges and Delhommelle (2016) C. Desgranges and J. Delhommelle, J. Chem. Phys. 144, 124510 (2016).
  * Gazenm$\ddot{\mathrm{u}}$ller and Camp (2007) G. Gazenm$\ddot{\mathrm{u}}$ller and P. J. Camp, J. Chem. Phys. 127, 154504 (2007).
  * Liu et al. (2011) Y. Liu, A. Z. Panagiotopoulos, and P. G. Debenedetti, J. Phys. Chem. B 115, 6629 (2011).
  * Nezbeda and Kolafa (1991) I. Nezbeda and J. Kolafa, Molec. Simul. 5, 391 (1991).
  * Singh and Errington (2006) J. K. Singh and J. R. Errington, J. Phys. Chem. B 110, 1369 (2006).
  * Rai et al. (2007) N. Rai, J. I. Siepmann, N. E. Schultz, and R. B. Ross, J. Phys. Chem. C 111, 15634 (2007).
  * Rane et al. (2013) K. S. Rane, S. Murali, and J. R. Errington, J. Chem. Theory Comput. 9, 2552 (2013).
  * Escobedo and de Pablo (1996) F. Escobedo and J. J. de Pablo, J. Chem. Phys. 105, 4391 (1996).
  * Shi and Maginn (2007) W. Shi and E. J. Maginn, J. Chem. Theory Comput. 3, 1451 (2007).
  * Eslami and Müller-Plathe (2007) H. Eslami and F. Müller-Plathe, J. Comput. Chem. 28, 1763 (2007).
  * Vogt et al. (2001) P. S. Vogt, R. Liapine, B. Kirchner, A. J. Dyson, H. Huber, G. Marcelli, and R. J. Sadus, Phys. Chem. Chem. Phys. 3, 1297 (2001).
  * Widom (1963) B. Widom, J. Chem. Phys. 39, 2808 (1963).
  * Siepmann and Frenkel (1992) J. Siepmann and D. Frenkel, Mol. Phys. 75, 59 (1992).