#  In situ sub-50 attosecond active stabilization of the delay between
infrared and extreme ultraviolet light pulses

Martin Luttmann<EMAIL_ADDRESS>Université Paris-Saclay, CEA, CNRS,
LIDYL, 91191 Gif-sur-Yvette, France  David Bresteau  Université Paris-Saclay,
CEA, CNRS, LIDYL, 91191 Gif-sur-Yvette, France  Jean-François Hergott
Université Paris-Saclay, CEA, CNRS, LIDYL, 91191 Gif-sur-Yvette, France
Olivier Tcherbakoff  Université Paris-Saclay, CEA, CNRS, LIDYL, 91191 Gif-sur-
Yvette, France  Thierry Ruchon<EMAIL_ADDRESS>Université Paris-Saclay,
CEA, CNRS, LIDYL, 91191 Gif-sur-Yvette, France

###### Abstract

The blooming of attosecond science (1 as = $10^{-18}$ s) has raised the need
to exquisitely control the delay between two ultrashort light pulses, one of
them being intense and in the visible spectral range, while the second is weak
and in the Extreme Ultra-Violet spectral range. Here we introduce a robust
technique, named LIZARD (Laser-dressed IoniZation for the Adjustment of the
pump-pRobe Delay), allowing an active stabilization of this pump-probe delay.
The originality of the method lies in an error signal calculated from a two-
photon photoelectron signal obtained by photoionizing a gas target in an
electronic spectrometer with the two superimposed beams. The modulation of
sidebands in phase quadrature allows us to perform an in situ measurement of
the pump-probe phase, and to compensate for fluctuations with an uniform noise
sensitivity over a large range of delays. Despite an interferometer length of
several meters, we achieved a long term stability of 28 as RMS over hours.
This method could be applied to the stabilization of other types of two-color
interferometers, provided that one of the propagating beams is capable of
photoionizing a target.

††preprint: AIP/123-QED

##  I Introduction

Time resolved experiments targeting attosecond precision [1, 2] are commonly
based on pump-probe schemes, requiring a stability of the two optical paths
better than a few tens of nanometers. For instance, a stability of 3 nm
corresponds to a delay resolution of 10 as. These requirements would not be
too challenging to fulfill for standard optical interferometers. In attosecond
science experiments, the situation is much less favorable for several reasons.
In one of the arms of the interferometer, a laser delivers femtosecond light
pulses in the Visible-Infrared spectral range (Vis-IR), which are focused in a
gas target to undergo a highly non linear process, High Harmonic Generation
(HHG). They are thus upconverted to the Extreme Ultra-Violet (XUV) range [3,
4] and, upon spectral filtering, form attosecond light pulses in the temporal
domain [5]. The remaining IR is filtered out downstream the focus, at a
sufficient distance to avoid any damage on the filters and optics, usually
tens of centimeters to meters away [6]. In addition, in order to allow the XUV
beam propagation, the interferometers are under secondary vacuum. It is the
combination of XUV wavelength, highly non-linear stage, meter-size
interferometer, and vacuum environment that makes the nanometer-grade
stability a challenging requirement. On top of that, a strong demand is
emerging for performing experiments that require long accumulation times (over
tens of hours), for instance coincidence spectroscopy [7, 8, 9], or
photoemission from surfaces [10]. In order to meet those needs, two main
experimental strategies are combined: increasing the repetition rate of the
laser to reduce the acquisition time, and improving, passively or actively,
the mechanical stability of the interferometer.

An excellent control of the environment, including temperature, humidity and
vibrations, but also the mechanical decoupling of the optical elements from
the pumping system, is usually a prerequisite of this second strategy. In
addition, the length of the interferometer and the number of optical elements
should be reduced as much as possible. This trend culminated in some original
experiments [11, 12], in which the two beams were recombined _before_ the HHG
setup. Although allowing attosecond pump-probe experiments in fairly difficult
conditions, it reduces the flexibility of the beamline and leads to a
modulation of the attosecond pulses intensity, which may be difficult to cope
with. However those efforts are usually not sufficient. Even under such
controlled environment, the thermal expansion of the mirrors due to residual
absorption of laser light still causes a dramatic drift of the interferometer
arms lengths. This is especially true for high repetition rate lasers, asking
for active stabilization systems.

Active stabilization systems using feedback loops have been implemented in a
series of configurations, all requiring significant modifications of the main
setup. Indeed, since the intense IR laser beam generally has to be filtered
out from the generation arm after HHG (with e.g. metallic foils), it cannot be
recombined to form a standard interference pattern at the output of the
interferometer. Any copropagating beam of comparable wavelength cannot either
follow the same full path as the XUV beam. Therefore most existing active
stabilization techniques rely on an auxiliary continuous wave (cw) laser
propagating on a shifted path along the XUV and IR beams in the interferometer
[13, 14, 15, 16, 17, 18]. It yields interference fringes which are used to
perform feedback on the delay by displacing mirrors. Although Vaughan et al.
[19] developed a different technique allowing a coaxial propagation, with the
outer part of the cw beam passing around the filter, it does not result in a
simplification of the setup. Indeed, the two interfering cw beams do not exit
their interferometer coaxially and have to be realigned in a dedicated optical
setup, which is itself stabilized using a second auxiliary laser.

Pump-probe setups in collinear geometry do not use metallic filters in the
interferometer, since the annular generation and the inner dressing beams are
recombined before the HHG gas jet [5, 20, 11]. In this case the cw laser can
be directly superimposed to the main beams in the interferometer. However, in
the end, this cw laser still has to be redirected to the detectors for the
active stabilization, which requires at least one additional optic.

We can finally mention the active stabilization scheme proposed by Schlaepfer
et al. [21]. The presence of an Optical Parametric Amplifier (OPA) in one of
the arms avoids the need for an auxiliary laser. They instead collect the
outer part of the IR generation beam passing around the filter and recombine
it with the 1560 nm light exiting the OPA, in a BBO crystal optimized for
second harmonic generation. The interference of the generation IR and the OPA
light back-converted to 780 nm is used for the active stabilization. Yet, CEP
stability of the two pulses is needed for this technique to work.

All stabilization setups mentioned above require at least one additional
optic, and often many more; hence the measured phase might fluctuate
differently from the actual phase of interest between the IR and XUV beams.
For instance, when the additional laser is guided on specific optics attached
to the mounts of the main optics, the measured fluctuations do not include
those caused by a thermal expansion of the main mirrors. Such residual power
absorption of the femtosecond intense laser get all the more important as the
repetition rate of the laser increases. In addition to that, existing active
stabilization schemes increase the complexity of the optical setup.

Here we introduce a completely different approach, without any auxiliary beam
or dedicated optical system, based on an error signal that directly depends on
the XUV/IR delay in the active region of a gas phase electronic spectrometer.
Specifically, we form an optical image of the point where the experiment is
carried out, e.g. the sensitive region of a coincidence spectrometer. We
collect and condition the two-photon ionization signal delivered by this
auxiliary spectrometer, in order to obtain a measurement of the XUV/IR delay
which is then used in a feedback loop. This technique offers an in situ
measurement and servo-locking of the delay of interest. In the following, we
will refer to this new technique as the LIZARD (Laser-dressed IoniZation for
the Adjustment of the pump-pRobe Delay). This article is organized as follows.
In part II, theoretical elements are introduced. Part III presents the
experimental performances of the LIZARD, and the main conclusions and a few
outlooks are discussed in part IV.

##  II Theory

###  II.1 Pump-probe delay measurement

Active stabilization techniques rely on interferometric signals that depend on
the pump-probe delay. For instance, if we use a cw laser, a photodetector can
be placed at the output of the interferometer and will deliver a signal which
varies sinusoidally with the delay. A feedback loop can be performed so that
this output signal is locked at zero [16, 19]. However, for applications
requiring delay scans, it is of prime importance to be able to continuously
lock the pump-probe phase at arbitrary set-points, with a uniform stability. A
solution consists in using two signals modulated with the delay, but not in
phase with one another [22, 23, 24]. Calling $\phi(t)$ the phase difference
between the pump and the probe beams, that we want to measure and stabilize,
such signals write

$\displaystyle M_{1}(t)$ $\displaystyle=O_{1}+A_{1}\sin(\phi(t)+\Delta)$ (1)
$\displaystyle M_{2}(t)$ $\displaystyle=O_{2}+A_{2}\cos{\phi(t)}.$

In the following, we assume that $M_{1}$ and $M_{2}$ are dimensionless
quantities. A fast preliminary calibration scan over the delay range gives a
measurement of the parameters $(O_{1},O_{2},A_{1},A_{2},\Delta)$, that we
assume constant over the stabilization time. If $\Delta=0$, $M_{1}$ and
$M_{2}$ are in phase quadrature; hence for any value of $\phi(t)$, either
$M_{1}$ or $M_{2}$ will have a significant slope, ensuring good dynamics of
the delay measurement. In the general case, measurements of $M_{1}$ and
$M_{2}$ allow the reconstruction of the phase

$\phi(t)=\mathcal{U}\Bigg{(}\tan^{-1}\bigg{(}\frac{1}{\cos\Delta}\times\frac{A_{2}\cdot(M_{1}(t)-O_{1})}{A_{1}\cdot(M_{2}(t)-O_{2})}-\tan\Delta\bigg{)}\Bigg{)},$
(2)

with $\mathcal{U}$ a phase unwrapping operator dealing with phase jumps, which
works as long as the phase does not vary by more than $\pm\pi$ during the
sampling time. More complex solutions for the phase reconstruction exist, such
as using more than two modulated signals [14].

Experimentally, due to detection noise, the values of $M_{1}$ and $M_{2}$ are
known within some uncertainties $\delta M_{1}$ and $\delta M_{2}$. If $\Delta$
is not negligible, this leads to an uncertainty $\delta\phi$ on the phase
measurement that strongly depends on the set-point phase. Indeed, $\delta\phi$
can be expressed as (details in Appendix A)

$\displaystyle\delta\phi=$
$\displaystyle\sqrt{\bigg{(}\frac{\partial\phi}{\partial M_{1}}\delta
M_{1}\bigg{)}^{2}+\bigg{(}\frac{\partial\phi}{\partial M_{2}}\delta
M_{2}\bigg{)}^{2}}$ $\displaystyle=$
$\displaystyle\sqrt{\bigg{(}\frac{\cos\phi}{A_{1}\cos\Delta}\delta
M_{1}\bigg{)}^{2}+\bigg{(}\frac{\sin(\phi+\Delta)}{A_{2}\cos\Delta}\delta
M_{2}\bigg{)}^{2}}.$ (3)

Figure 1: (color online) Uncertainty on the reconstructed pump-probe phase due
to the noise in the signals $M_{1}$ and $M_{2}$, for different values of
$\Delta$. Both graphs were plotted for $\delta M_{1}=\delta M_{2}=\delta M$
and $A_{1}=2$. a) Case where $A_{1}=A_{2}$ ($M_{1}$ and $M_{2}$ have the same
signal to noise ratio). b) An example with $A_{2}=A_{1}/2$ (the signal to
noise ratio is lower for $M_{2}$ than for $M_{1}$).

The behavior of this function is shown in Fig. 1 in some specific cases. If
the two signals have equal signal to noise ratio (Fig. 1.a), $\delta\phi$ is
uniform over the whole pump-probe phase range when $\Delta=0$. This does not
hold if $M_{1}$ and $M_{2}$ have different signal to noise ratios (Fig. 1.b).
However, in both cases, the average value of $\delta\phi$ is minimized when
$\Delta=0$. The phase quadrature is therefore an optimal situation for a
precise active stabilization.

The originality of the LIZARD lies in the origins of the phase quadrature
signals, which are extracted from the photoelectron spectrum obtained in a
laser-dressed photoionization experiment.

###  II.2 Principle of laser-dressed ionization by a HHG frequency comb

The setup is similar to generic pump-probe experiments in attosecond physics
(see Fig. 4). The IR laser (of angular frequency $\omega_{L}$) is separated
into two beams using a beam splitter. Most of the intensity is used for HHG,
but a small portion, called "dressing", propagates unchanged. The two beams
are then recombined and focused in a Time Of Flight Magnetic Bottle Electronic
Spectrometer (TOF-MBES). With multicycle driving laser pulses, the
photoelectron energy spectrum exhibits peaks located at the harmonics
energies. The field in the gas jet being a superposition of XUV and IR, two-
photon transitions are also allowed. They are represented for different
harmonics in Fig. 2. Two transitions – absorption of a photon from harmonic
$H_{q}$ and an IR photon; absorption of a photon from harmonic $H_{q+2}$ and
stimulated emission of an IR photon – yield the same photoelectron energy (red
dashed lines in Fig. 2), creating sidebands between the harmonic peaks in the
photoelectron spectrum.

Figure 2: (color online) Sketch of the two-photon transitions involved in
laser-dressed ionization. XUV and IR photons are represented with purple and
red arrows respectively. The photoelectron spectrum (right part) shows the
main harmonics peaks, spaced by twice the energy of the driving laser. Weaker
sidebands are located between those peaks, at even multiples of the driving
laser energy.

These two quantum paths interfere, which leads to a modulation of the
sidebands intensity with the time delay $\tau$ between the XUV and IR fields,
at twice the laser frequency [25]

$I_{q}(\tau)=O_{q}+A_{q}\sin{(2\omega_{L}\tau+\phi_{\text{at}}+\phi_{q+1}-\phi_{q-1})},$
(4)

with $I_{q}$ the intensity of sideband $SB_{q}$ (at energy $q\omega_{L}$),
$\phi_{\text{at}}$ an usually negligible intrinsic atomic phase, and
$\phi_{q}$ the spectral phase of the XUV field at energy $q\omega_{L}$. The
term $\phi_{q+1}-\phi_{q-1}$ in Eq. (4) is therefore the increment of the XUV
spectral phase between the two successive harmonics surrounding sideband
$SB_{q}$. $O_{q}$ and $A_{q}$ are the offset and amplitude of $SB_{q}$, which
depend on the intensity of harmonics $H_{q-1}$ and $H_{q+1}$, the intensity of
the dressing beam, and detection settings.

The modulation of the sidebands with the pump-probe delay $\tau$, coupled to
the fact that they are not in phase with one another, make them excellent
candidates for an active stabilization of the delay. Quite remarkably, the
physics of HHG, which gives rise to an ’attosecond chirp’[26], will often
provide sidebands sufficiently separated in phase so that at least two of them
can be found close to phase quadrature. Indeed, the three steps model of HHG –
in which an outer electron successively undergoes tunnel ionization,
acceleration in the continuum and recombination with the parent ion – predicts
that the recombination time of the electron will vary almost linearly with the
photon energy of the emitted harmonic (Fig. 3). The recombination time
difference between the highest and the lowest part of the emitted spectrum
(represented with two triangles in Fig. 3) is approximately $0.3\,T_{L}$, with
$T_{L}$ the IR optical cycle, and does not depend on the generating laser
intensity. The corresponding dephasing between the highest and the lowest
sidebands is thus given by $\Delta\Phi\approx 2\omega_{L}\times
0.3\,T_{L}=1.2\,\pi$. As it is greater than $\pi/2$, whatever the generating
medium and the laser intensity, one should always be able to find two
sidebands in the photoelectron spectrum whose phase difference is close enough
to $\pi/2$ to serve as the $M_{1}$ and $M_{2}$ signals in the LIZARD scheme.
As a downside, this also means that attosecond pulses which are perfectly
post-compressed (by dedicated devices) would not allow a straightforward
implementation of the LIZARD, since they do not possess chirp.

It may be noted that the high harmonic peaks can be used as well: they are
also modulated with the delay, with a phase just opposite to that of the
neighbouring sidebands. This is why, in the following, ’delay modulated
signals’ will indifferently refer to sidebands or harmonics.

Figure 3: (color online) Electron recombination time calculated within the
strong field approximation [27], as a function of the harmonic order. We
considered an Argon atom submitted to a 800 nm wavelength. We disregarded any
envelope effect of the driving pulse. The short and long trajectories are
represented with continuous and dotted lines respectively. The simulation was
performed for different values of the IR laser electric field. The two
triangles on the vertical axis indicate roughly the range of the XUV emission
time in the case of short trajectories. Figure 4: (color online) Experimental
setup. BS: beam splitter, PM: plane mirror, TM: toroidal mirror. The
interferometer arms are approximately 2.5 m long. The Cold Target Recoil Ion
Momentum Spectrometer (COLTRIMS) experiment in the end station is a possible
application of the LIZARD; it is not part of the stabilization system itself.

##  III Experimental implementation and results

To test the LIZARD scheme, experiments were conducted on the user-oriented
FAB10 beamline at the ATTOLab facility. The beamline is based on a Ti:Sapphire
chirped pulse amplification laser [28] delivering 23 fs pulses carrying 2 mJ,
with a 10 kHz repetition rate. Focused onto a continuous Argon gas jet, it
drives a HHG source emitting attosecond pulses trains (APT). The beamline is
optimized for experiments requiring fairly high recurrence and long
acquisition times, for instance coincidence measurement in gas phase [9] or
photoemission in solids [29]. RABBIT [5, 30] spectroscopy is permanently
operational on FAB10, which allows the characterization of the XUV light and
the experiments in the end station to be performed simultaneously. For that
purpose, the XUV beam is first focused in the TOF-MBES, then re-focused by a
gold mirror in the user’s experimental chamber. The parity of the number of
mirrors and focusing optics is identical in the two arms of the interferometer

To make sure that we are sensitive to delay fluctuations only, we have to
reduce as much as possible the noise on the signals $M_{1}$ and $M_{2}$. It
may be due to fluctuations of the dressing intensity and the XUV harmonic
intensity. The latter will cause changes in the total number of emitted
photoelectrons per second. To limit this effect, a beam pointing stabilization
of the IR the laser is performed. Moreover, the TOF-MBES signal will be
averaged over a large number of laser shots (which is possible thanks to the
high repetition rate of the beamline), and normalized. The resulting value at
energy $q\omega_{L}$ is denoted $\hat{I}_{q}$.

The time delay $\tau$ between the IR dressing and the XUV probe is controlled
by two mirrors mounted on a translation stage (model SmarAct SLC-24 150-P-R-S,
based on a stick-slip technology combining long travel range (10 cm) and piezo
accuracy (1 nm)) before the HHG chamber (Fig. 4). The electronic signal from
the magnetic bottle’s microchannel plate (MCP) detector is acquired with an
oscilloscope (model Lecroy Waverunner 6 Zi).

Before starting the stabilization, a fast delay scan is performed in order to
record the oscillations of the two selected sidebands. The parameters $O_{1}$,
$O_{2}$, $A_{1}$, $A_{2}$ and $\Delta$ of Eq. (1) are extracted from
sinusoidal fits of the data. During the stabilization process, the feedback
loop performed at each step consists in reconstructing the pump-probe phase
difference $\phi(t)$ using Eq. (2), calculating the error value
$\Delta\phi(t)=\phi(t)-\phi_{0}$ (with $\phi_{0}$ the desired set-point phase)
and performing feedback using a PID (proportional-integral-derivative) filter.
Those operations are performed by a GUI program based on PyMoDAQ [31], an open
source Python framework for data acquisition programs. Since we mostly need to
compensate for slow thermal delay drifts, a 10 Hz correction frequency is high
enough. It allows us to average the photoelectron signal over 1000 laser
shots, ensuring a reasonable level of noise.

###  III.1 Stabilization at a fixed set-point

Figure 5: (color online) Comparison between open loop (blue) and closed loop
(orange) operation. a) Long term measurement of the phase error
$\Delta\phi(t)=\phi(t)-\phi_{0}$. The inset shows the evolution of the delay
stage position when the setpoint optical paths difference is changed by 100
nm, in the case of a non-optimized (dashed line) and optimized (continuous
line) PID. As expected, the kick yields an overall 50 nm displacement of the
stage (because of the mirrors corner cube configuration). For an optimized
PID, we consider that the new setpoint value is reached in approximately 0.5
s. b) Phase error histogram. The full width at half maximum for the closed
loop error is indicated with a red line. Assuming a gaussian distribution, it
should be approximately $2\sqrt{\ln 2}\approx 1.67$ times the RMS error (224
mrad). c) Power spectral density of the phase error $\delta\phi(t)$. The low
frequency noise, corresponding to long term drifts, is efficiently suppressed
by the LIZARD.

A typical long-term LIZARD measurement of the pump-probe phase difference, in
open loop, is shown in Fig. 5.a (blue curve). We observe a drift which is
remarkably linear, with a rate of $36\text{\,}\rad\mathrm{/}\mathrm{h}$ $\sim$
$7.6\text{\,}\mathrm{fs}\text{\,}{\mathrm{h}}^{-1}$ $\sim$
$2.3\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{h}}^{-1}$, that we
attribute to the thermal expansion of the mirrors. In other words, the
difference of optical paths varies by a sideband spatial period in 10 minutes.
This passive stability is too poor to maintain attosecond precision during
experiments lasting several hours, forcing the use of convoluted post-analysis
protocols [9].

The measurement in closed loop is shown as the orange curve in Fig. 5.a. The
two modulated signals chosen for the active stabilization are $SB_{14}$ and
$SB_{22}$. They are separated in phase by $0.46\,\pi\,$ rad, very close to
phase quadrature. We reached a phase stability of 134 mrad RMS, equivalent to
a 28 as time delay stability and a 9 nm optical path stability, which is
competing with state of the art performance [11, 21, 19, 13, 14, 15, 17]. Our
active stabilization system turns out to be robust and easy to operate. After
a proper alignment of the beamline, it only takes one minute to select two
energy bands, make a calibration scan and launch a long-term stabilization.

Figure 6: (color online) Comparison between a classical open loop RABBIT scan
(a,b,c) and a RABBIT scan stabilized by LIZARD (d,e,f). (a,d) Modulation of
the photoelectron spectrum with the delay $\tau$, after background subtraction
and normalization. In the stabilized scan (d), $\hat{I}_{\text{14}}$ and
$\hat{I}_{\text{22}}$ (white dotted bands) were used as the two modulated
signals $M_{1}$ and $M_{2}$ for the active stabilization. (b,e)
$\hat{I}_{\text{16}}$ (blue) and $\hat{I}_{\text{20}}$ (green), as a function
of the delay. The signals are integrated over the width of the bands with the
corresponding colors shown in a) and d). The black curves are sinusoidal fits.
(c,f) Squared modulus of the Fourier transforms of $\hat{I}_{\text{16}}(\tau)$
and $\hat{I}_{\text{20}}(\tau)$. The average value of each of the two
modulated signal was first subtracted in order to remove very low frequency
components. The vertical dashed line shows the theoretical $2\omega_{L}$
sidebands angular frequency. The raw Fourier transforms (dotted lines) were
smoothed using zero padding (solid lines).

###  III.2 Stabilized pump-probe delay scans

Fig. 5 presented the case of an interferometer locked at a constant value of
the pump-probe phase. The phase can however be locked at any arbitrary desired
set-point. In order to demonstrate the capabilities of the technique, we
performed open loop and closed loop RABBIT scans. In a classical open loop
RABBIT scan, the pump-probe optical path difference is typically scanned over
$4\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ with $40\,\mathrm{n}\mathrm{m}$
steps (Fig. 6.a). Such a scan ideally spans over exactly 10 sideband periods,
but is distorted by delay drifts (Fig. 6.b). As the delay drift in our
experiment is primarily linear, it results in a shift of the frequency of the
sidebands oscillation (Fig. 6.c). Using a least squared method, we fitted the
data with the sinusoidal function defined by Eq. (4). Importantly, the
sidebands frequency is not a parameter of the fit which adjusts only the
offset, amplitude and phase.

Fig. 6.d is a LIZARD-stabilized RABBIT scan, obtained with the following
procedure. The optical paths difference is stabilized and the set-point is
regularly adjusted by steps of 40 nm. The optimization of the PID ensures that
the new set-point value is reached in only a few tenth of seconds. Data
corresponding to the same set-point are averaged in post-analysis. In the
experiment presented here, the phase remained locked at each set-point during
50 photoelectron spectrum acquisitions. For the relevancy of the comparison,
both the non stabilized and the stabilized scan had the same duration,
approximately 5 minutes. With the active stabilization, the delay scan matches
the sinusoidal fit almost perfectly (Fig. 6.e). The sidebands Fourier
transform (Fig. 6.f) contains a thin peak nicely centered at $2\omega_{L}$.
The stabilized scan approaches the perfect sinusoidal oscillations of Eq. (4).
Another way to visualize the effect of the LIZARD stabilization is to look at
the correlation between two sidebands intensities (Fig. 7), similarly to
diagnosis performed with FEL [32] and the recently developed TURTLE [33].
Since the pump-probe delay is now precisely adjusted, the stabilization
results in a squeezing of the data points around the phase set-points.

Figure 7: (color online) Intensity correlation between four different
sidebands extracted from the spectrum of Fig. 6, and sideband $SB_{22}$. Data
corresponding to the open loop scan and to the closed loop scan are in blue
and orange respectively. Black squares indicate the average positions of the
sets of data corresponding to the same set-point phase. All axis have equal
unit lengths. The increasing tilt of the ellipses is due to the attochirp. a)
$SB_{20}$ versus $SB_{22}$. b) $SB_{18}$ versus $SB_{22}$. c) $SB_{16}$ versus
$SB_{22}$. d) $SB_{14}$ versus $SB_{22}$. Those two particular sidebands were
used for the LIZARD stabilization, they are in phase quadrature and yield a
non-tilted correlation ellipse. Figure 8: (color online) a) RABBIT scan in
argon stabilized by LIZARD that lasted approximately 2 hours. Between each 20
nm optical paths change, spectra were accumulated during 65 s, which
corresponds to $6.5\times 10^{5}$ laser shots. The yellow curve shows the
$2\omega$ amplitude. The very sinusoidal oscillations of two of the sidebands
(SB18 and SB20) are plot as an example in blue and green. (b-c-d-e) Spectral
resolution of some sidebands phases (Rainbow RABBIT technique [34]), in the
case of the 2 hours long stabilized scan (orange) and a fast (30 s)
calibration scan recorded just before (blue). f) RMS noise of the spectral
phase measurement, with respect to the number of laser shots over which the
spectra were averaged during each scan step.

One could argue that we should not be surprised to observe almost perfect
sinusoidal sidebands with 2$\omega_{L}$ angular frequency in the case of
closed loop RABBIT scans, since we used the RABBIT photoelectron spectrum
itself to perform the stabilization. An important point is that the active
stabilization is not using the whole photoelectron signal but only two
specific energy bands. $SB_{16}$ and $SB_{20}$, studied as an example in Fig.
6.d-f, are not the sidebands used for the LIZARD as the input signals $M_{1}$
and $M_{2}$. We can therefore consider that the information carried by the
remaining sidebands is preserved in a stabilized scan. The comparison between
Fig. 6.b and Fig. 6.e confirms that the sidebands sinusoidal fit, calculated
with an imposed angular frequency of $2\omega_{L}$, matches the experimental
data more precisely in the case of a stabilized RABBIT scan. Thus the
estimation of the phase difference $\Delta\Phi_{\text{16-20}}$ between
$SB_{16}$ and $SB_{20}$ (which corresponds to a photoemission time delay
between two harmonics) will be more accurate in the case of a stabilized scan.
To verify this, we use our sinusoidal fit to estimate the uncertainty on the
measurement of $\Delta\Phi_{\text{16-20}}$. The standard deviation associated
to the measured phase of a given sideband is obtained by considering the
variance of the sideband data points equal to the variance of the residual
(the distance between the data and the fit). Let us call
$\Delta\Phi_{\text{16-20}}^{\text{o-l}}$ and
$\Delta\Phi_{\text{16-20}}^{\text{c-l}}$ the phase difference between the two
sidebands extracted from the open loop and the closed loop RABBIT scan
respectively. We found, with a $2\sigma$ uncertainty,
$\Delta\Phi_{\text{16-20}}^{\text{o-l}}=$0.64\text{\,}$(\pm$0.55\text{\,}\mathrm{)}$\rad$
and
$\Delta\Phi_{\text{16-20}}^{\text{c-l}}=$0.67\text{\,}\mathrm{(}$\pm$0.06\text{\,}\mathrm{)}$\rad$.
Here, the active stabilization reduces the uncertainty on the measurement by a
factor of 8.

Finally, Fig. 8 illustrates the feasibility of delay scans lasting several
hours. It is readily visible from Fig. 8.b-e that accumulating long spectra
while stabilizing with LIZARD drastically increases the precision of spectral
phase retrieval, especially at low energies, where the $2\omega$ amplitude is
weak. To quantify this, we computed the RMS distance between the spectral
phase resulting from the averaging over intermediate accumulation times, and
the spectral phase obtained by averaging over 65 s, considered as a reference
(orange curves in Fig. 8.b-e). The results are displayed in Fig. 8.f. As
expected, although increasing the accumulation time has a negligible effect
where the signal to noise ratio is high ($SB_{16}$ and $SB_{18}$), the benefit
is dramatic for the weak signals of $SB_{12}$ and $SB_{14}$. Because it allows
such long averaging times, we can expect the LIZARD to become a precious tool
for the measurement of spectral phases when the signal to noise ratio is
particularly poor, due for instance to weak cross sections of the channels of
interest [35].

##  IV Conclusion and outlook

Our technique is particularly suited to attosecond beamlines that are already
equipped with an electron time of flight spectrometer. It is fundamentally
different from existing active stabilization schemes, which use the
interference of photons, and measure a pump-probe phase difference accumulated
during the propagation of a cw laser in the interferometer. The LIZARD
technique does not require any additional cw laser or even additional optical
elements but a refocusing optic that is common to both beams. Its principle
relies on the interference of photoelectron wavepackets, therefore providing
an in situ measurement of the XUV/IR phase difference in the active region of
the electronic spectrometer. The underlying physics of HHG guarantees that the
spectrum contains two sidebands in phase quadrature in the vast majority of
generating conditions. This allowed us to stabilize the delay with a 28 as RMS
accuracy during several hours - which is comparable to state of the art
stability performance - and to adjust the set-point delay over a few microns
range, thus enabling pump-probe experiments requiring long acquisition times.
We also showed that this technique can bring higher precision in the
measurement of photoemission time delays by the RABBIT technique. From a
practical perspective, the simplification of the optical system is
counterbalanced by the need of an electron spectrometer. This device is
however quite common on attosecond beamlines, and could even be used in a
simpler version collecting only two time of flight windows, one for each
modulated signal.

Further improvements of the method are intended. We could use the whole
photoelectron spectrum instead of only two sidebands, so that the pump-probe
delay could be deduced from a least squared optimization of all the sidebands
intensity. A statistical study showed that the accuracy of the delay
measurement increases with the number of sidebands considered. We could also
increase the correction frequency by using analog electronics instead of a
numerical code. However, since this would reduce the number of averaged shots
per delay measurements, it is not certain that it would improve the stability.
This novel method could possibly be extended to active stabilization schemes
with isolated attosecond pulses (IAP). Although they exhibit a continuous
frequency spectrum, the standard IAP temporal reconstruction techniques, FROG-
CRAB [36] and PROOF [37], are also based on laser-dressed photoionization, and
yield delay-modulated spectrograms. By integrating the spectrum over well
chosen energy bands, one might extract the required phase quadrature delay-
modulated signals, and use the LIZARD to actively stabilize pump-probe
experiments with IAP.

Finally, we believe that this approach could be applied to the stabilization
of two-color interferometers in other spectral domains, as long as one of the
beams is able to photoionize a target. Here we described in details the
situation where one of the two pulses is in the XUV spectral range, therefore
photoionization is possible thanks to the high energy of the photons involved.
However, even at lower photons energies (for instance, consider an infrared
pulse and its second harmonic), photoelectrons and oscillating sidebands could
be obtained by above-threshold ionization of gas targets [38, 39].

###### Acknowledgements.

The authors acknowledge Bertrand Carré, Carlo Spezzani, François Polack, David
Dennetiere, Ismaël Vadillo-Torre, Michel Bougeard and Julien Lenfant who
contributed to the design and construction of the FAB10 beamline; and Pascal
D’Oliveira and Fabien Lepetit for the design, implementation and operation of
the laser. The laser system and the experimental setup are supported by the
French “Investments for the Future” of the Agence Nationale pour la Recherche
(Contract No. 11-EQPX0005-ATTOLAB and ANR-14-CE32-0010 - Xstase), the
Scientific Cooperation Foundation of Paris-Saclay University through the
funding of the OPT2X research project (Lidex 2014), by the Île-de-France
region through the Pulse-X project and by the European Union’s Horizon 2020
Research and Innovation Programme No. EU-H2020-LASERLAB-EUROPE-654148.

##  Appendix A Derivation of Eq. (II.1)

The tangent of the reconstructed phase $\phi$ writes, according to Eq. (1)

$\tan\phi=\frac{1}{\cos\Delta}\times\bigg{(}\frac{M_{1}-O_{1}}{M_{2}-O_{2}}\times\frac{A_{2}}{A_{1}}\bigg{)}-\tan\Delta$
(5)

Let us differentiate this equation with respect to $M_{1}$

$\frac{\partial\tan\phi}{\partial
M_{1}}=\frac{1}{\cos\Delta}\bigg{(}\frac{A_{2}}{A_{1}(M_{2}-O_{2})}\bigg{)}=\frac{1}{A_{1}\cos\phi\cos\Delta}$
(6)

Similarly, we obtain for the derivative with respect to $M_{2}$

$\frac{\partial\tan\phi}{\partial
M_{2}}=\frac{-\sin(\phi+\Delta)}{A_{2}\cos^{2}\phi\cos\Delta}$ (7)

The derivatives of $\phi$ with respect to $M_{1}$ and $M_{2}$ are obtained
using $\text{d}\phi=\cos^{2}\phi\times\text{d}(\tan\phi)$, yielding the
following equations

$\displaystyle\frac{\partial\phi}{\partial M_{1}}$
$\displaystyle=\frac{\cos\phi}{A_{1}\cos\Delta}$ (8)
$\displaystyle\frac{\partial\phi}{\partial M_{2}}$
$\displaystyle=-\frac{\sin(\phi+\Delta)}{A_{2}\cos\Delta}$

By definition, the uncertainty on the phase measurement is

$\delta\phi=\sqrt{\bigg{(}\frac{\partial\phi}{\partial M_{1}}\delta
M_{1}\bigg{)}^{2}+\bigg{(}\frac{\partial\phi}{\partial M_{2}}\delta
M_{2}\bigg{)}^{2}}$ (9)

We obtain the last term of Eq. (II.1) by injecting Eq. (8) in Eq. (9).

## References

  * Calegari  _et al._ [2016] F. Calegari, G. Sansone, S. Stagira, C. Vozzi, and M. Nisoli, Journal of Physics B Atomic Molecular Physics 49, 062001 (2016).
  * Ueda  _et al._ [2019] K. Ueda, E. Sokell, S. Schippers, F. Aumayr, H. Sadeghpour, J. Burgdörfer, C. Lemell, X.-M. Tong, T. Pfeifer, F. Calegari, A. Palacios, F. Martin, P. Corkum, G. Sansone, E. V. Gryzlova, A. N. Grum-Grzhimailo, M. N. Piancastelli, P. M. Weber, T. Steinle, K. Amini, J. Biegert, N. Berrah, E. Kukk, R. Santra, A. Müller, D. Dowek, R. R. Lucchese, C. W. McCurdy, P. Bolognesi, L. Avaldi, T. Jahnke, M. S. Schöffler, R. Dörner, Y. Mairesse, L. Nahon, O. Smirnova, T. Schlathölter, E. E. B. Campbell, J.-M. Rost, M. Meyer, and K. A. Tanaka, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 171001 (2019).
  * Ferray  _et al._ [1988] M. Ferray, A. L'Huillier, X. F. Li, L. A. Lompre, G. Mainfray, and C. Manus, Journal of Physics B: Atomic, Molecular and Optical Physics 21, L31 (1988).
  * McPherson  _et al._ [1987] A. McPherson, G. Gibson, H. Jara, U. Johann, T. S. Luk, I. A. McIntyre, K. Boyer, and C. K. Rhodes, J. Opt. Soc. Am. B 4, 595 (1987).
  * Paul  _et al._ [2001] P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Augé, P. Balcou, H. G. Muller, and P. Agostini, Science 292, 1689 (2001).
  * Ye  _et al._ [2020] P. Ye, T. Csizmadia, L. G. Oldal, H. N. Gopalakrishna, M. Füle, Z. Filus, B. Nagyillés, Z. Divéki, T. Grósz, M. Dumergue, P. Jójárt, I. Seres, Z. Bengery, V. Zuba, Z. Várallyay, B. Major, F. Frassetto, M. Devetta, G. D. Lucarelli, M. Lucchini, B. Moio, S. Stagira, C. Vozzi, L. Poletto, M. Nisoli, D. Charalambidis, S. Kahaly, A. Zaïr, and K. Varjú, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 154004 (2020).
  * Cirelli  _et al._ [2018] C. Cirelli, C. Marante, S. Heuser, C. L. M. Petersson, Á. J. Galán, L. Argenti, S. Zhong, D. Busto, M. Isinger, S. Nandi, S. Maclot, L. Rading, P. Johnsson, M. Gisselbrecht, M. Lucchini, L. Gallmann, J. M. Dahlström, E. Lindroth, A. L’Huillier, F. Martín, and U. Keller, Nature Communications 9, 955 (2018).
  * Gallmann  _et al._ [2017] L. Gallmann, I. Jordan, H. J. Wörner, L. Castiglioni, M. Hengsberger, J. Osterwalder, C. A. Arrell, M. Chergui, E. Liberatore, U. Rothlisberger, and U. Keller, Structural Dynamics 4, 061502 (2017).
  * Joseph  _et al._ [2020] J. Joseph, F. Holzmeier, D. Bresteau, C. Spezzani, T. Ruchon, J.-F. Hergott, O. Tcherbakoff, P. D’Oliveira, J. Houver, and D. Dowek, J. Phys. B: At. Mol. Opt. Phys. 53 (2020).
  * Dendzik  _et al._ [2020] M. Dendzik, R. P. Xian, E. Perfetto, D. Sangalli, D. Kutnyakhov, S. Dong, S. Beaulieu, T. Pincelli, F. Pressacco, D. Curcio, S. Y. Agustsson, M. Heber, J. Hauer, W. Wurth, G. Brenner, Y. Acremann, P. Hofmann, M. Wolf, A. Marini, G. Stefanucci, L. Rettig, and R. Ernstorfer, Phys. Rev. Lett. 125, 096401 (2020).
  * Weber  _et al._ [2015] S. J. Weber, B. Manschwetus, M. Billon, M. Böttcher, M. Bougeard, P. Breger, M. Géléoc, V. Gruson, A. Huetz, N. Lin, Y. J. Picard, T. Ruchon, P. Salières, and B. Carré, Review of Scientific Instruments 86, 033108 (2015).
  * Zaïr  _et al._ [2018] A. Zaïr, E. Mével, E. Cormier, and E. Constant, J. Opt. Soc. Am. B 35, A110 (2018).
  * Huppert  _et al._ [2015] M. Huppert, I. Jordan, and H. J. Wörner, Review of Scientific Instruments 86, 123106 (2015).
  * Böttcher  _et al._ [2008] F. Böttcher, B. Manschwetus, H. Rottke, N. Zhavoronkov, Z. Ansari, and W. Sandner, Appl. Phys. B , 287 (2008).
  * Sabbar  _et al._ [2014] M. Sabbar, S. Heuser, R. Boge, M. Lucchini, L. Gallmann, C. Cirelli, and U. Keller, Review of Scientific Instruments 85, 103113 (2014).
  * Fieß  _et al._ [2010] M. Fieß, M. Schultze, E. Goulielmakis, B. Dennhardt, J. Gagnon, M. Hofstetter, R. Kienberger, and F. Krausz, Review of Scientific Instruments 81, 093103 (2010).
  * Chini  _et al._ [2009] M. Chini, H. Mashiko, H. Wang, S. Chen, C. Yun, S. Scott, S. Gilbertson, and Z. Chang, Opt. Express 17, 21459 (2009).
  * Mashiko  _et al._ [2010] H. Mashiko, M. J. Bell, A. R. Beck, M. J. Abel, P. M. Nagel, C. P. Steiner, J. Robinson, D. M. Neumark, and S. R. Leone, Opt. Express 18, 25887 (2010).
  * Vaughan  _et al._ [2019] J. Vaughan, J. Bahder, B. Unzicker, D. Arthur, M. Tatum, T. Hart, G. Harrison, S. Burrows, P. Stringer, and G. M. Laurent, Opt. Express 27, 30989 (2019).
  * Mairesse  _et al._ [2003] Y. Mairesse, A. de Bohan, L. J. Frasinski, H. Merdji, L. C. Dinu, P. Monchicourt, P. Breger, M. Kovačev, R. Taïeb, B. Carré, H. G. Muller, P. Agostini, and P. Salières, Science 302, 1540 (2003).
  * Schlaepfer  _et al._ [2019] F. Schlaepfer, M. Volkov, N. Hartmann, A. Niedermayr, Z. Schumacher, L. Gallmann, and U. Keller, Opt. Express 27, 22385 (2019).
  * Krishnamachari  _et al._ [2006] V. V. Krishnamachari, E. R. Andresen, S. R. Keiding, and E. O. Potma, Opt. Express 14, 5210 (2006).
  * Jotzu  _et al._ [2012] G. Jotzu, T. J. Bartley, H. B. Coldenstrodt-Ronge, B. J. Smith, and I. A. Walmsley, Journal of Modern Optics 59, 42 (2012).
  * Freschi and Frejlich [1995] A. A. Freschi and J. Frejlich, Opt. Lett. 20, 635 (1995).
  * Véniard  _et al._ [1996] V. Véniard, R. Taïeb, and A. Maquet, Phys. Rev. A 54, 721 (1996).
  * Varjú  _et al._ [2005] K. Varjú, Y. Mairesse, B. Carré, M. B. Gaarde, P. Johnsson, S. Kazamias, R. López-Martens, J. Mauritsson, K. J. Schafer, P. Balcou, A. L’huillier, and P. Salières, Journal of Modern Optics 52, 379 (2005).
  * Lewenstein  _et al._ [1994] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, Phys. Rev. A 49, 2117 (1994).
  * Golinelli  _et al._ [2017] A. Golinelli, X. Chen, E. Gontier, B. Bussière, O. Tcherbakoff, M. Natile, P. d’Oliveira, P.-M. Paul, and J.-F. Hergott, Opt. Lett. 42, 2326 (2017).
  * Fanciulli  _et al._ [2020] M. Fanciulli, J. Schusser, M.-I. Lee, Z. E. Youbi, O. Heckmann, M. C. Richter, C. Cacho, C. Spezzani, D. Bresteau, J.-F. Hergott, P. D’Oliveira, O. Tcherbakoff, T. Ruchon, J. Minár, and K. Hricovini, Phys. Rev. Research 2, 013261 (2020).
  * Muller [2002] H. Muller, Appl Phys B 74 (2002).
  * [31] https://github.com/CEMES-CNRS/PyMoDAQ.
  * Maroju  _et al._ [2020] P. K. Maroju, C. Grazioli, M. Di Fraia, M. Moioli, D. Ertel, H. Ahmadi, O. Plekan, P. Finetti, E. Allaria, L. Giannessi, G. De Ninno, C. Spezzani, G. Penco, S. Spampinati, A. Demidovich, M. B. Danailov, R. Borghes, G. Kourousias, C. E. Sanches Dos Reis, F. Billé, A. A. Lutman, R. J. Squibb, R. Feifel, P. Carpeggiani, M. Reduzzi, T. Mazza, M. Meyer, S. Bengtsson, N. Ibrakovic, E. R. Simpson, J. Mauritsson, T. Csizmadia, M. Dumergue, S. Kühn, H. Nandiga Gopalakrishna, D. You, K. Ueda, M. Labeye, J. E. Bækhøj, K. J. Schafer, E. V. Gryzlova, A. N. Grum-Grzhimailo, K. C. Prince, C. Callegari, and G. Sansone, Nature 578, 386 (2020).
  * Krüger  _et al._ [2020] C. Krüger, J. Fuchs, L. Cattaneo, and U. Keller, Opt. Express 28, 12862 (2020).
  * Gruson  _et al._ [2016] V. Gruson, L. Barreau, Á. Jiménez-Galan, F. Risoud, J. Caillat, A. Maquet, B. Carré, F. Lepetit, J.-F. Hergott, T. Ruchon, L. Argenti, R. Taïeb, F. Martín, and P. Salières, Science 354, 734 (2016).
  * Klünder  _et al._ [2011] K. Klünder, J. M. Dahlström, M. Gisselbrecht, T. Fordell, M. Swoboda, D. Guénot, P. Johnsson, J. Caillat, J. Mauritsson, A. Maquet, R. Taïeb, and A. L’Huillier, Phys. Rev. Lett. 106, 143002 (2011).
  * Mairesse and Quéré [2005] Y. Mairesse and F. Quéré, Phys. Rev. A 71, 011401 (2005).
  * Chini  _et al._ [2010] M. Chini, S. Gilbertson, S. D. Khan, and Z. Chang, Opt. Express 18, 13006 (2010).
  * Agostini  _et al._ [1979] P. Agostini, F. Fabre, G. Mainfray, G. Petite, and N. K. Rahman, Phys. Rev. Lett. 42, 1127 (1979).
  * Beaulieu  _et al._ [2017] S. Beaulieu, A. Comby, A. Clergerie, J. Caillat, D. Descamps, N. Dudovich, B. Fabre, R. Géneaux, F. Légaré, S. Petit, B. Pons, G. Porat, T. Ruchon, R. Taïeb, V. Blanchet, and Y. Mairesse, Science 358, 1288 (2017).


# Simulating Crystallization in a Colloidal System Using State Predictive
Information Bottleneck based Enhanced Sampling

Vanessa J. Meraz  Institute for Physical Science and Technology, University of
Maryland, College Park 20742, USA.  Ziyue Zou  Department of Chemistry and
Biochemistry, University of Maryland, College Park 20742, USA.  Pratyush
Tiwary<EMAIL_ADDRESS><EMAIL_ADDRESS>Institute for Physical Science and
Technology, University of Maryland, College Park 20742, USA.  Department of
Chemistry and Biochemistry, University of Maryland, College Park 20742, USA.
University of Maryland Institute for Health Computing, Rockville 20852, USA.

(August 28, 2024)

###### Abstract

## Abstract

We investigate crystal nucleation in supersaturated colloid suspensions using
enhanced molecular dynamics simulations augmented with machine learning
techniques. The simulations reveal that crystallization in the model colloidal
system studied here, with particles interacting through a repulsive screened
Coulomb Yukawa potential, proceeds from vapor to dense liquid droplet to
crystalline phases across multiple high barriers. Employing a one-dimensional
reaction coordinate derived from the State Predictive Information Bottleneck
framework, our simulations capture back-and-forth phase transitions across
multiple barriers effectively in biased metadynamics simulations. We obtain
relative free energy differences between different phases and also quantify
the roles of different molecular level features in driving the phase changes.

##  I Introduction

Fundamental to physical chemistry, nucleation is critical in the study of
kinetics as the key initiation event of crystallization. Effectively
controlling the nucleation process of molecules or particles within a solvent
is of great importance amongst different scientific disciplines [1, 2, 3, 4].
In all cases, a cluster forms within the solvent and when it is sufficiently
large enough, the free energy barrier is overcome by the gain of emerging
crystalline order. While ubiquitous in nature and well observed
macroscopically, the mechanism behind this event is not so well understood.
Even at the start of the process, the initial crystallization greatly
influences the polymorph produced [5, 6]. Generally, the basis for studying
nucleation is built upon classical nucleation theory (CNT), a one-step
theoretical model, which considers only one variable necessary to capture the
complex dynamics resulting in phase transformations: the size of the largest
cluster, assumed to be spherical [7, 8].

Though simple and effective, CNT has been continuously improved upon, the most
prominent and extensively studied modification of which has been the two-step
model [7]. Through the simplistic view of CNT, there is a single free energy
barrier that needs to be overcome for crystallization to be energetically
favorable. The modified two-step mechanism accounts for metastable states
between the starting supersaturated solution and the emerging crystal. Hence,
there is more than one free energy barrier whose cost needs to be offset by
crystal growth. This idea of multistep crystallization is motivated by the
importance of density fluctuations when clusters are being formed. For
example, the colloidal system studied in this work starts as a vapor phase
which can then grow into a dense liquid droplet (DLD) if local density
fluctuations promote the clustering of particles which then may or may not
favor the ordering necessary for crystal emergence. While we can generalize
the steps for a system to nucleate, the path for crystallization itself is not
unique. As such, molecular dynamics (MD) simulations could be particularly
useful in directly observing crystallization [9]. However, nucleation is a
rare event whose timescale varies amongst systems (within seconds to minutes)
and is typically longer than what is accessible and achievable by MD.

In this work, we study crystallization in supersaturated colloid suspensions
[10]. The particles in our system interacted through a repulsive screened
Coulomb Yukawa potential [11], which are commonly used to describe charge-
stabilized colloidal suspensions [12]. Due to their larger size, colloids make
for a challenging model system in the study of nucleation. In some sense,
their slow movement allows for an easier look into dynamics that may not be
all apparent in atomistic systems. However the interaction between colloidal
particles themselves presents unique challenges. Specifically, in a stable
suspension, there is weak interaction between colloidal particles lending to
the difficulty of spontaneous aggregation into crystals.

The crux of the problem in understanding nucleation lies in the atomistic
detail necessary to glean information as to what promotes or inhibits it, a
formidable feat even for MD with state-of-the-art instrumentation. For this
reason, we turn to enhanced MD algorithms to reproduce the crystallization
process. Unsurprisingly, there have been several efforts in attempt to sample
a larger portion of configuration space within the allotted simulation time
[13]. One enhanced sampling approach involves learning an approximation to the
true, a priori unknown, reaction coordinate (RC) that can represent the
progress along a generic reaction pathway. Generally, such an RC is
constructed as a function of collective variables (CVs), or order parameters
(OPs) [14, 15, 16, 17], which are functions of particles’ positional
coordinates. For enhanced sampling of rare events, one then constructs biasing
potential or forces as a function of the RC, thereby speeding up slow
processes. An important caveat to note about the use of OPs used to build the
RC is their requirement to be continuous and differentiable.

In this work, we bias the approximate RC through the well-tempered
metadynamics method [18, 19, 20, 21, 22, 13]. The algorithm artificially fills
free energy wells with a bias potential as a function of the approximate RC
enabling an enhanced exploration of the landscape. We use the State Predictive
Information Bottleneck (SPIB) framework to construct the RC [23]. In the most
broad sense and in methods like SPIB, the RC is constructed as a linear or
nonlinear combination of OPs [16, 24, 25, 26, 27]. Using OPs that are
computationally cost-effective to calculate and differentiate, as decided by
their ranking detailed in Ref. [14], a linear RC is constructed to best
enhance crystallization. We find that such a one-dimensional RC is sufficient
to observe back and forth phase transitions in metadynamics simulations.
Through this, we observe different metastable states relevant to
crystallization in the colloidal system. Starting from the vapor phase, we
find that the system first goes to a dense liquid droplet phase, followed by a
crystal phase. We observe multiple back-and-forth transitions between these
phases. Additionally, oftentimes the system continues to a more dense crystal
phase, at which point it gets trapped there. This behavior reflects both the
strengths and limitations of our SPIB based RC. By reweighting the well-
tempered metadynamics simulations, we also calculate the relative free energy
differences between these different metastable states.

While machine learning based reaction coordinate methods[13] have been
recently applied to a large number of biophysical systems, their application
to nucleation in colloidal systems is rather limited. Our work shows that
these methods can also be directly applied to the study of nucleation in model
colloids. We thus believe the protocol demonstrated in this work will be
useful for gathering thermodynamic insights into generic colloidal systems.

##  II Methods

The development of enhanced sampling methods has been of great use in the
atomistic simulation of crystal nucleation which would otherwise be improbable
to observe even with the best available supercomputers. In this work, we are
interested in studying nucleation in a system of supersaturated colloid
suspensions using a Yukawa potential [28]. To simulate nucleation and obtain
accurate associated thermodynamic observables, here we use the enhanced
sampling method well-tempered metadynamics [20] biasing an approximate RC
learnt from an SPIB model.

###  1 Nucleation Order Parameters

SPIB constructs an approximate reaction coordinate (RC) as a function of input
features or order parameters (OPs). Through there isn’t a clear-cut definition
for which OPs to use for studying nucleation, there have been extensive
studies on the relative importance of different possible OPs [14, 17, 29]. In
this text, we focus on OPs studied in ref. [14], and in particular those of
which we find to be cost effective in terms of simulation time used to compute
these OPs and forces acting on them. They are respectively classified into two
categories: average properties of all particles (energy, coordination number,
Steinhardt bond orders) and total number of particles according to local
geometry (particles in a condensed or solid-like phase). They are listed in
Table 1 for reference and are defined as follows:

Table 1: Order parameters (OPs) input to SPIB Label | OP
---|---
cn.mean | Mean particle first-sphere coordination number
ncl | Number of particles in a condensed phase (CN $>$ 3)
ncs | Number of particles in a solid-like phase (CN $>$ 6)
q4.mean | Mean Q4 Steinhardt bond order
q6.mean | Mean Q6 Steinhardt bond order
ene | Potential energy

  1. 1.

Coordination Number: cn.mean

Eq.1 is used to compute the mean first-sphere coordination number, $cn.mean$,
where $r_{ij}$ is the distance between particles and $r_{max}$ is the distance
within the first-coordination sphere. It follows from the general definition
of a coordination number for a particle $i$ as given by Eq. 2 where $r_{o}$ is
a cut-off radial distance.

$cn.mean=\frac{1}{N}\sum_{i,j\neq
i}^{N}\frac{f(r_{ij})-f(r_{max})}{1-f(r_{max})}$ (1)

$f(r_{ij})=\sum_{j}\frac{1-(r_{ij}/r_{o})^{6}}{1-(r_{ij}/r_{o})^{12}}$ (2)

Here the summation variable $j$ enumerates all other particles. A general
descriptor of a particles’ local environment, the coordination number is a
simple number OP. It fails to capture the symmetry of a crystal structure and
thus is not able to distinguish emerging phases making it inadequate as a
stand-alone OP.

  2. 2.

Number of particles in a condensed or solid-like phase: ncl, ncs

Calculating the number of particles in a phase is similar to the description
of the coordination number $cn$. The $cn$ is calculated over all particles,
for which there is a distribution, and a filter, 1- f(r = x), is then applied
to classify the phase as condensed or solid-like. Bounds x=3 and x=7 were used
to find ncl and ncs respectively. The term f(x=r) is the same as seen in Eq.
2.

  3. 3.

Steinhardt bond order parameters: q4.mean, q6.mean

To quantify the symmetries in a given cluster made of different particles,
Steinhardt bond order parameters are widely used [30]. These OPs are functions
of spherical harmonics and describe the ordering of the local environment
surrounding a given particle. Here, we set the $l^{th}$ order of the spherical
harmonics $\mathbf{Y_{lm}}$ to 4 and 6 as they result in interesting
characteristics of crystal structures [31]. Accordingly, the $m^{th}$ order
has ranges {-4,…,+4} and {-6,…,+6} with respect to the $l^{th}$ order. This
results in 9 or 13 spherical harmonic functions to be evaluated as seen in Eq.
3, where $i,j$ are the same as in Eq. 1:

$q_{lm}=\frac{\sum_{j}r(r_{ij})\bf{Y}_{lm}(\bf{r}_{ij})}{\sum_{j}r(r_{ij})}$
(3)

  4. 4.

Potential energy: ene

Each time step of a simulation requires the calculation of potential energy as
a function of the force field chosen. As we intend to study phase
transformations, the potential energy of a system is a good indicator one has
taken place. As seen in Fig. 1, the increase in order of a system (from vapor-
to-liquid or liquid-to-crystal) results in a significant potential energy
difference. In the system studied in this work, a colloid/Yukawa [28]
potential was used.

###  2 State-Predictive Information Bottleneck (SPIB)

At the heart of understanding nucleation, and countless other phenomena in
physical chemistry, sits the need to project very high-dimensional data onto a
much lower dimensional manifold. To effectively address this task, it’s
crucial to have a concise and informed representation capable of
distinguishing between different states within a system, as well as capturing
a correct state-to-state transition probability density. As previously
discussed, a multitude of variables are necessary for accurately describing
transitions, hence a combination of OPs is employed to craft such a
representation. Within this framework, we employ the State Predictive
Information Bottleneck (SPIB), functioning akin to a variational autoencoder
(VAE).

SPIB serves to preserve pertinent dynamics between metastable states by
compressing a high-dimensional Molecular Dynamics (MD) trajectory into a low-
dimensional representation of latent space which can be shown to approximate
certain critical requirements of the true RC [32]. This representation
extracts the essential information from an MD trajectory, minimizing
redundancy while maximizing its predictive power regarding future states
following a specified time delay, denoted as $\Delta t$.

For the RC to be able to learn the number of metastable states and
differentiate between them, a variational mixture of posteriors prior
(VampPrior)[33] is used to construct a multimodal probability distribution
prior for SPIB. The prior is approximated as a weighted combination of, in our
case, OPs with initialized pseudo-inputs to be updated on-the-fly during
training. For a particular choice of time delay $\Delta t$, iterations of SPIB
will converge to an $N$ number of metastable states whose population
probability is reflected in the updated pseudo-inputs. Generally speaking, the
number of converged metastable states decreases as the time-delay $\Delta t$
is increased. See Ref. [23] for further description of SPIB.

The choice of encoder, linear or not, is left to the application. While a non-
linear encoder can be more analytically expressive, the need to interpret
results from biased metadynamics simulations restricts our use to a linear
encoder. Furthermore, for poor quality training data a non-linear encoder can
also be prone to overfitting and less capable of extrapolation outside
training domain, which is often critical and especially so in biased
simulations.

###  3 Enhanced Sampling with Metadynamics

Since nucleation is a rare event which is difficult to simulate directly in
unbiased MD, one needs to use enhanced sampling methods to simulate them. Here
we use well-tempered metadynamics. Metadynamics is a technique used to explore
a broader portion of the configuration space by applying a history-dependent
bias in the free energy landscape, described as a function of one to two
carefully chosen biasing variables. In the well-tempered variation of
metadynamics [20], the height of the Gaussian added at any given point in
biasing variable space is gradually reduced during the simulation. Gaussian
biasing potentials are periodically introduced as function of these biasing
variables, or the approximated RC, to discourage the system from getting stuck
in local free energy minima. The net product of a converged metadynamics
simulation is the free energy as function of the biased or any other variables
through a reweighting procedure [34]. While metadynamics is mathematically
proven to converge for any biasing variable [35], in practice the speed of
convergence is improved if the variables being biased capture all relevant
slow degrees of freedom, i.e. approximate the true reaction coordinate [22].
Here we achieve this through the use of SPIB.

Specifically, we use SPIB to generate a one-dimensional RC as a linear
combination of input OPs provided in Table 1 which requires initial
trajectories that have seen at least one phase transition. This can be
obtained by either biasing a trial RC or by launching a very large number of
unbiased runs and selecting one where nucleation happens. Here we do so by
running 1000 independent unbiased simulations, 21 of which show nucleation,
and hence are used to create this RC. Having sampled this space, the
metastable states of the system will be able to be discretized with some
degree of accuracy.

In Table 2, we report the parameters used for well-tempered metadynamics
simulations. These are (i) the initial height of the deposited Gaussian
$\omega$ in units of $k_{B}T^{*}$, (ii) bias factor $\gamma$ (unitless), (iii)
width of the Gaussian $\sigma$ in same units as the SPIB RC being biased; and
the reduced temperature $T^{*}$ defined as $k_{B}T$ / $\epsilon$ . The
Gaussian width $\sigma$ was kept at 5 times from the standard deviation of the
RC measured in a highly paced biased run to simulate an unbiased one which
lasted 3.25 x $10^{5}$ steps with timestep $\Delta$t = 0.005.

Table 2: Parameters for metadynamics (Definitions and units are provided in Sec. 3) $\omega$ | $\gamma$ | $\sigma$ | $T^{*}$
---|---|---|---
2.0 | 50 | 1.3781 | 2

##  III Results

###  1 Observations from initial Molecular Dynamics Simulations

The colloidal simulations performed in this work were done using LAMMPS [36]
patched with PLUMED [37, 38], following the detailed set-up reported in Ref.
[14].

We began by launching 1000 independent unbiased MD simulations starting from
the vapor phase to gather instances where crystallization occurred. As in Ref.
[14], this was characterized through the polyhedral template matching (PTM
[39]) tool in LAMMPS. Following that same protocol, we impose the condition
that the crystal structure formed during the simulation must consist of a
minimum of 188 particles. This condition adheres to the threshold value for
potential energy per particle E* = -10 $\epsilon$ seen in Fig.1. Of the
initial 1000 simulations, only 21 trajectories formed a crystalline structure.
This is similar to the crystallization rate around 1% observed in Ref. [14],
in tune with the rare event nature of colloid nucleation.

Furthermore, none of these select successful unbiased trajectories in our work
or in Ref. [14], once crystallized, ever revisited the starting vapor phase.
This illustrates that unbiased simulations such as these are not capable of
providing equilibrium free energy profile for the colloidal system.

Figure 1: Mean potential energy in reduced units, $E^{*}$, per particle as a
function of time measured in MD integration timesteps for all 21 crystallized
trajectories illustrated in different shades of purple.

Figure 2: (a) Ratio of hcp/fcc crystals plotted against $E^{*}(\epsilon)$
which is the average of potential energy measured over the last 5,000 frames.
(b) Number of particles not classified as fcc or hcp out of total 388
particles plotted against $E^{*}(\epsilon)$.

In Fig. 1, we plot the mean potential per particle as a function of MD
integration timesteps for the 21 unbiased MD trajectories that showed
crystallization. As was also reported in Ref. [14], there are clear features
representative of phase change. Of the 21 trajectories, 3 are of particular
interest due to their formation of a densely packed structure, with minimal
particles excluded. In Fig. 1, these three instances can be seen where the
mean potential energy goes below E* = -20 $\epsilon$. To understand this
phenomenon better, we plot the ratio of particles classified as hexagonal
close-packed (hcp) versus face-centered cubic (fcc) in Fig. 2(a). The 3 lowest
energies recorded in Fig. 1 correspond to the 3 structures seen as outliers in
the bottom left of Fig. 2(a). They have the same relative hcp/fcc ratio but
markedly lower potential energy. This is further seen in Fig. 2b where the
number of particles not classified as fcc or hcp are plotted against the
reduced mean potential energy and the three points of interest seemingly
behave as outliers.

To construct the RC for biasing in metadynamics, the order parameter outputs
from these crystallized trajectories are given as training input data into
SPIB. The learned RC is constructed as a one-dimensional linear encoder as a
function of OPs described in Sec. 1. For simplicity, we choose a time lag long
enough ($\Delta t$ = 10,000 MD integration timesteps) for states to be
discretized between crystal and non-crystal basins, as seen with distinct
colors in Fig. 3. The projection of the RC along different order parameters
seen in the figure also illustrates the rationale behind the previously
defined threshold for crystallization (around E* = -10 $\epsilon$).

Figure 3: Projection of SPIB RC onto different OPs. In (a) we plot the average
potential energy ene versus the Steinhardt bond order parameter q6. In (b) we
plot the average potential energy versus the coordination number cn. In both
(a) and (b), the different shades of blue indicate the two states SPIB
converged to given a time lag of $\Delta t$ = 10,000.

###  2 Observations from Metadynamics biasing SPIB

As constructed by the one-dimensional linear SPIB encoder, the RC is expressed
as a combination of the OP weights seen in Fig. 4. From the figure, we note
the high contribution of the OP capturing number of particles in a solid-like
phase (ncs) followed by the mean coordination number (cn), average potential
energy, and Steinhardt bond order parameter q6. Conversely, the Steinhardt
bond order parameter q4 contributes the least.

Figure 4: Order parameter (OP) weights from RC obtained through SPIB.

We performed 10 independent well-tempered metadynamics simulations biasing
along the SPIB based RC, which we label $z$. In contrast to the original
$\sim$2% rate of crystallization observed in the initial unbiased simulations,
all metadynamics runs indicated crystal formation as well as back-and-forth
movement between some phases. Fig. 5 shows the 4 different phases of the
system observed during our metadynamics simulations. The categorization of
these 4 phases as a function of the SPIB RC, $z$, is as follows: vapor, for
$z<$ 20; dense liquid droplet, for 20 $<z<$ 40; crystal, for 50 $<z<$ 70, and
a dense crystal phase for $z>90$.

In Fig. 6 we plot the coordination number, cn, as a function of simulation
time for a typical metadynamics trajectory (see SI for other nine metadynamics
trajectories showing qualitatively similar behavior). This figure shows a
clear back and forth movement between 3 of the 4 phases mentioned in the
previous paragraph. In other words, our trajectories move back-and-forth
between the vapor, dense liquid and crystal phases. However, for 8 of the 10
trajectories, the system reaches a dense crystal phase with $z>90$.

If this phase is reached, our simulations do not revert back. This reflects
that the RC learnt through SPIB captured the physics relevant to moving
between the first three phases including the crystal, but not for movement
back from the dense crystal. One possible explanation for this is the under-
representation of trajectories visiting the dense crystal phase in our
training data, where only 3 out of 21 trajectories visited this phase. In
Section 3 we elaborate qualitatively on the basis of the free energy profile.

Figure 5: Phases seen during MD simulation as rendered by Ovito. Within sub-
figures (c) and (d), particles in red are in the hcp phase and particles in
green are in the fcc phase. Figure 6: Coordination number cn as a function of
simulation time (measured in MD timesteps) for typical biased and one of the
successful unbiased trajectories, shown in colors blue and black respectively.

###  3 Free Energy Profile

####  1 Initial Molecular Dynamics Simulations

In Fig. 7(a) we provide the free energy approximated from all the initial
unbiased simulations, irrespective of whether they crystallized or not,
projecting onto the energy and q6 OPs. Also shown in Fig. 7(a) are three
metastable crystal basins likewise reported in Ref. [14]. At
$E^{*}<-20\epsilon$, we observe the dense crystal structure (Fig. 5(d)) in
three instances of the crystallized trajectories. However, given enough
simulation time, the system inevitably reaches this phase as shown by the
unbiased trajectory in Fig. 6. The transition from crystal to dense crystal
($E^{*}\sim-20\epsilon$) also accounts for a high energy barrier in the free
energy profile.

Figure 7: Free energy projection onto q6 and energy for (a) 1000 and (b) 21
unbiased simulations.

In Fig. 7(b) we provide the same approximate free energy obtained by
considering only the 21 trajectories that displayed crystallization. These are
the trajectories used as input to SPIB. Compared to Fig. 7(a), there is a
singular crystal basin which has similar depth to the liquid basin and the
barrier between them is relatively high although not as pronounced as in Fig.
7(a). Interestingly, both Figs. 7(a) and 7(b) have a similar barrier height
for the crystal-to-dense crystal transition while both vapor-to-liquid and
liquid-to-crystal barriers are lower when only considering crystallized
trajectories. This is quite likely due to the poor sampling of this region as
the likelihood of the system exploring this space within the simulation time
is more improbable than the system crystallizing ($\sim 2\%$ vs $0.3\%$). This
point is highlighted in our metadynamics results. We emphasize that Fig. 7(a)
and (b) are just approximations to the free energy for the sake of
visualization and developing intuition, as these simulations display
hysteresis and do not correspond to equilibrium. We provide the full converged
free energy in the next sub-section, obtained using SPIB-biased metadynamics.

####  2 Metadynamics biasing SPIB

We calculate a one dimensional free energy profile for the colloidal system
from each of the 10 metadynamics trajectories by reweighting [34]. Fig. 8
provides the one-dimensional averaged free energy (in black), the error bars
between them (in red), and the classification of phases as determined by OPs.
The boundaries of the phases themselves are loosely based on the one-
dimensional free energy projection of the OPs, provided in more detail in SI.
Understandably, there is a correlation between the OP weights in Fig. 4 and
the free energy projections which allow us to reliably define the phase along
SPIB coordinate $z$.

Figure 8: One-dimensional free energy profile from metadynamics in reduced
units as a function of the SPIB coordinate $z$. Figure 9: Free energy
projection onto Stenihardt bond order parameter q6 and potential energy from
metadynamics simulations.

In Fig. 9 we provide a two-dimensional projection of the free energy along
coordinates q6, Steinhardt bond order parameter, and ene, average potential
energy. This shows significant differences relative to the initial unbiased
projection in Fig. 7. Remarkably, there are two main wells: that of the
crystal and dense crystal phase. Just above in the island where
$E^{*}>-5\epsilon$, the vapor phase is seen. The barrier height difference
between the crystal and dense crystal phase in Fig. 9 is much more pronounced
in Fig. 8. This also explains why once the system crosses into the dense
crystal phase it is markedly more difficult for it to go back to the crystal
phase.

##  IV Conclusion

In this article, we have used machine learning augmented molecular dynamics
simulations to study the emergence of different ordered phases in a system of
supersaturated colloid suspensions. Specifically, we use the State Predictive
Information Bottleneck (SPIB)[23] for learning a one-dimensional reaction
coordinate (RC). We see that such a simplified RC is effective in capturing
phase transitions in well-tempered metadynamics simulations. This allows us to
observe various stable states relevant to the crystallization process in
colloidal systems which are identified as rare events under conventional
molecular dynamics simulations. Generally, our simulations display multiple
back-and-forth transitions between the vapor, liquid and crystal phases. From
the vapor state, we see the system initially shifting to a dense liquid
droplet stage before entering a crystalline phase. In addition, there are
seemingly direct to-and-from vapor-crystal transitions with only a brief, if
any, time spent in the dense liquid phase. Furthermore, it’s common for the
system to progress into an even denser crystal state where it remains after
it’s been reached. This reflects both the strengths and the limitations of our
approach. We believe the methodology demonstrated in this study will prove
valuable in gaining thermodynamic insights into phase transitions in various
colloidal systems with the introduction of linear or non-linear combinations
of conventional order parameters. To address a broader interest in the study
of nucleation problems which, in general, require a larger simulation cell,
our protocol can be readily applicable.

## Acknowledgments

This work was entirely supported by the U.S. Department of Energy under Award
Number DE-SC0021009. We thank UMD HPC’s Zaratan and NSF ACCESS (project
CHE180027P) for computational resources. Pratyush Tiwary is an investigator at
the University of Maryland-Institute for Health Computing, which is supported
by funding from Montgomery County, Maryland and The University of Maryland
Strategic Partnership: MPowering the State, a formal collaboration between the
University of Maryland, College Park and the University of Maryland,
Baltimore.

## References

  * Sosso et al. [2016] Sosso, G. C.; Li, T.; Donadio, D.; Tribello, G. A.; Michaelides, A. Microscopic Mechanism and Kinetics of Ice Formation at Complex Interfaces: Zooming in on Kaolinite. _The Journal of Physical Chemistry Letters_ 2016, _7_ , 2350–2355.
  * Fersht [1997] Fersht, A. R. Nucleation mechanisms in protein folding. _Current Opinion in Structural Biology_ 1997, _7_ , 3–9.
  * Aizenberg [2004] Aizenberg, J. Crystallization in Patterns: A Bio‐Inspired Approach. _Advanced Materials_ 2004, _16_ , 1295–1302.
  * Baumgartner et al. [2013] Baumgartner, J.; Dey, A.; Bomans, P. H. H.; Le Coadou, C.; Fratzl, P.; Sommerdijk, N. A. J. M.; Faivre, D. Nucleation and growth of magnetite from solution. _Nature Materials_ 2013, _12_ , 310–314.
  * Zou et al. [2021] Zou, Z.; Tsai, S.-T.; Tiwary, P. Toward Automated Sampling of Polymorph Nucleation and Free Energies with the SGOOP and Metadynamics. _The Journal of Physical Chemistry B_ 2021, _125_ , 13049–13056.
  * Zou et al. [2023] Zou, Z.; Beyerle, E. R.; Tsai, S.-T.; Tiwary, P. Driving and characterizing nucleation of urea and glycine polymorphs in water. _Proceedings of the National Academy of Sciences_ 2023, _120_.
  * Karthika et al. [2016] Karthika, S.; Radhakrishnan, T. K.; Kalaichelvi, P. A Review of Classical and Nonclassical Nucleation Theories. _Crystal Growth & Design_ 2016, _16_ , 6663–6681.
  * Sosso et al. [2016] Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.; Michaelides, A. Crystal Nucleation in Liquids: Open Questions and Future Challenges in Molecular Dynamics Simulations. _Chemical Reviews_ 2016, _116_ , 7078–7116.
  * Finney and Salvalaglio [2022] Finney, A. R.; Salvalaglio, M. Multiple pathways in NaCl homogeneous crystal nucleation. _Faraday Discussions_ 2022, _235_ , 56–80.
  * Cacciuto et al. [2004] Cacciuto, A.; Auer, S.; Frenkel, D. Onset of heterogeneous crystal nucleation in colloidal suspensions. _Nature_ 2004, _428_ , 404–406.
  * Everaers and Ejtehadi [2003] Everaers, R.; Ejtehadi, M. R. Interaction potentials for soft and hard ellipsoids. _Physical Review E_ 2003, _67_.
  * Kremer et al. [1986] Kremer, K.; Robbins, M. O.; Grest, G. S. Phase Diagram of Yukawa Systems: Model for Charge-Stabilized Colloids. _Phys. Rev. Lett._ 1986, _57_ , 2694–2697.
  * Mehdi et al. [2024] Mehdi, S.; Smith, Z.; Herron, L.; Zou, Z.; Tiwary, P. Enhanced Sampling with Machine Learning. _Annual Review of Physical Chemistry_ 2024, _75_ , null, PMID: 38382572.
  * Finney and Salvalaglio [2023] Finney, A. R.; Salvalaglio, M. A variational approach to assess reaction coordinates for two-step crystallization. _The Journal of Chemical Physics_ 2023, _158_ , 094503.
  * Neha et al. [2022] Neha; Tiwari, V.; Mondal, S.; Kumari, N.; Karmakar, T. Collective Variables for Crystallization Simulations from Early Developments to Recent Advances. _ACS Omega_ 2022, _8_ , 127–146.
  * Pérez-Hernández et al. [2013] Pérez-Hernández, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noé, F. Identification of slow molecular order parameters for Markov model construction. _The Journal of Chemical Physics_ 2013, _139_.
  * Santiso and Trout [2011] Santiso, E. E.; Trout, B. L. A general set of order parameters for molecular crystals. _The Journal of Chemical Physics_ 2011, _134_ , 064109\.
  * Laio and Parrinello [2002] Laio, A.; Parrinello, M. Escaping free-energy minima. _Proceedings of the National Academy of Sciences_ 2002, _99_ , 12562–12566.
  * Laio and Gervasio [2008] Laio, A.; Gervasio, F. L. Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. _Reports on Progress in Physics_ 2008, _71_ , 126601.
  * Barducci et al. [2008] Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. _Phys. Rev. Lett._ 2008, _100_ , 020603.
  * Salvalaglio et al. [2016] Salvalaglio, M.; Tiwary, P.; Maggioni, G. M.; Mazzotti, M.; Parrinello, M. Overcoming time scale and finite size limitations to compute nucleation rates from small scale well tempered metadynamics simulations. _The Journal of Chemical Physics_ 2016, _145_ , 211925.
  * Bussi and Laio [2020] Bussi, G.; Laio, A. Using metadynamics to explore complex free-energy landscapes. _Nat. Rev. Phys._ 2020, 1–13.
  * Wang and Tiwary [2021] Wang, D.; Tiwary, P. State predictive information bottleneck. _The Journal of Chemical Physics_ 2021, _154_.
  * Tiwary and Berne [2016] Tiwary, P.; Berne, B. J. Spectral gap optimization of order parameters for sampling complex molecular systems. _Proceedings of the National Academy of Sciences_ 2016, _113_ , 2839–2844.
  * Ma and Dinner [2005] Ma, A.; Dinner, A. R. Automatic method for identifying reaction coordinates in complex systems. _The Journal of Physical Chemistry B_ 2005, _109_ , 6769–6779.
  * Best and Hummer [2005] Best, R. B.; Hummer, G. Reaction coordinates and rates from transition paths. _Proc. Natl. Acad. Sci._ 2005, _102_ , 6732–6737.
  * Peters and Trout [2006] Peters, B.; Trout, B. L. Obtaining reaction coordinates by likelihood maximization. _The Journal of Chemical Physics_ 2006, _125_ , 054108.
  * Safran [2018] Safran, S. A. _Statistical Thermodynamics of Surfaces, Interfaces, and Membranes_ ; CRC Press, 2018.
  * Karmakar et al. [2021] Karmakar, T.; Invernizzi, M.; Rizzi, V.; Parrinello, M. Collective variables for the study of crystallisation. _Molecular Physics_ 2021, _0_ , e1893848.
  * Steinhardt et al. [1983] Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. Bond-orientational order in liquids and glasses. _Phys. Rev. B_ 1983, _28_ , 784–805.
  * Lechner and Dellago [2008] Lechner, W.; Dellago, C. Accurate determination of crystal structures based on averaged local bond order parameters. _The Journal of Chemical Physics_ 2008, _129_.
  * Wang et al. [2022] Wang, D.; Zhao, R.; Weeks, J. D.; Tiwary, P. Influence of Long-Range Forces on the Transition States and Dynamics of NaCl Ion-Pair Dissociation in Water. _The Journal of Physical Chemistry B_ 2022, _126_ , 545–551.
  * Tomczak and Welling [2017] Tomczak, J. M.; Welling, M. VAE with a VampPrior. _CoRR_ 2017, _abs/1705.07120_.
  * Tiwary and Parrinello [2014] Tiwary, P.; Parrinello, M. A Time-Independent Free Energy Estimator for Metadynamics. _The Journal of Physical Chemistry B_ 2014, _119_ , 736–742.
  * Dama et al. [2014] Dama, J. F.; Parrinello, M.; Voth, G. A. Well-Tempered Metadynamics Converges Asymptotically. _Physical Review Letters_ 2014, _112_.
  * Thompson et al. [2022] Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in ’t Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. _Computer Physics Communications_ 2022, _271_ , 108171.
  * Bonomi et al. [2019] Bonomi, M.; Bussi, G.; Camilloni, C. C. Promoting transparency and reproducibility in enhanced molecular simulations. _Nat. Methods._ 2019, _16_ , 670–673.
  * Tribello et al. [2014] Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Comp. Phys. Comm._ 2014, _185_ , 604–613.
  * Larsen et al. [2016] Larsen, P. M.; Schmidt, S.; Schiøtz, J. Robust structural identification via polyhedral template matching. _Modelling and Simulation in Materials Science and Engineering_ 2016, _24_ , 055007.