2022The Authors. This work is an open access article under the CC BY-SA 4.0
license FOAM-extend 4.1 https://bitbucket.org/C_Habes/lcs4foam Preprint
submitted to OpenFOAM Journal

# lcs4Foam – An OpenFOAM Function Object to
Compute Lagrangian Coherent Structures

C. Habes1,∗ A. Kameke2,∗ M. E. Fadeli1 H. Marschall1 1Mathematical Modeling
and Analysis, Technical University of Darmstadt, 64287 Darmstadt, Germany
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>2Department of Mechanical Engineering and
Production Management, Hamburg University of Applied Sciences, 20099 Hamburg,
Germany<EMAIL_ADDRESS>
###### Abstract

To facilitate the understanding and to quantitatively assess the material
transport in fluids, a modern characterisation method has emerged in fluid
dynamics in the last decades footed in dynamical systems theory. It allows to
examine the most influential material lines which are called Lagrangian
Coherent Structures (LCS) and order the material transport into dynamically
distinct regions at large scales which resist diffusion or mixing. LCS reveal
the robust skeleton of material surfaces and are essential to assess material
transport in time-dependent flows quantitatively. Candidates of LCS can be
estimated and visualised from finite-time stretching and folding fields by
calculating the Finite-Time Lyapunov Exponents (FTLE).

In this contribution, we provide an OpenFOAM function object to compute FTLE
during CFD simulation. This enables the OpenFOAM community to assess the
geometry of the material transport in any flow quantitatively on-the-fly using
principally any OpenFOAM flow solver.

##  1 Introduction

Material transport and mixing in fluids is enhanced by advection. This
advection is usually described mathematically in an Eulerian view by a time-
dependent velocity field $\mathbf{u}(\mathbf{u},\,t)$. With this Eulerian
description, numerous important fluid mechanical characteristics can be
derived and assessed. For instance, a higher Reynolds number (higher
velocities) typically will go along with better overall mixing. However, such
intuition might be misleading as has been shown for example in [1] studying a
rising bubble. Here, a coherent structure has been found to arise for
intermediate Reynolds numbers, which causes material to move together and
locally hinders mixing and increases residence times in the vicinity of the
bubble rear. The example shows: a closer look at the coherent structures is
necessary to evaluate the details of the material transport in the specific
flow situation. Lagrangian Coherent Structures (LCS) are often observable in
fluid flows due to the shape that passive tracers take on, e.g. plankton in
the ocean [2] or dissolved oxygen in the wake behind a rising bubble.

That the classical Eulerian view on advection is not optimal for addressing
these issues was first noted in oceanography and atmospheric science [3, 4].
The transport analysis was therefore started from its roots, the Lagrangian
view, where the observer travels on the fluid parcels rather than watching
them move by (Eulerian frame). The Lagrangian analysis thus considers the
trajectories of individual fluid parcels and allows to draw conclusions on the
transport from their evaluation. Nowadays, computational and theoretical
advances allow for the calculation and analysis of the time-dependent
dynamical system that governs material transport. The underlying ideas for
Lagrangian analysis stem from dynamical systems theory. In time-independent
incompressible velocity fields, the dynamical system is the velocity field
itself and the streamlines of the velocity field coincide with the
trajectories of the fluid parcels. As such, trivially, structures in the
velocity fields represent governing structures for the material that is
transported (as long as molecular diffusion is comparably low and negligible)
[5]. In this setting, unstable and stable manifolds divide the flow into
different subdomains that move coherently (together) [6]. For time-dependent
flows however, the instantaneous streamlines and the trajectories of the fluid
parcels do not coincide. It is thus a misleading habit to draw any conclusion
about the material transport from the streamlines or any other material lines
of the mean velocity field of a fluid flow. The resulting transport structures
might have no relevance for the real dynamical system at all.

To obtain the lines that govern material transport in time-dependent flows the
Lagrangian Coherent Structures are calculated from the trajectories of
particles evaluated in the time-dependent velocity field
$\mathbf{u}=\mathbf{u}(\mathbf{x},\,t)$. LCS are those material lines and
surfaces that separate regions of particles with very different fates or
history for the time interval under consideration. Several different
approaches to evaluate LCS have been developed during the last years [7].

With this contribution we introduce an OpenFOAM function object that
calculates the three dimensional Finite Time Lyapunov Exponents (FTLE) on-the-
fly based on the general purpose numerical library _libcfd2lcs_ [8] with the
main computational details explained in [9]. The ridges in the FTLE-field are
then candidates for LCS and can be assumed to coincide with LCS if some
further conditions are met [10]. However, as also pointed out in [8], these
additional conditions are hard to evaluate in 3D and thus the FTLE-field will
be viewed as an approximate representation of the 3D LCS. The details about
the calculation of the FTLE-field and the underlying mathematical foundation
are set out in Section 2.

##  2 Theoretical background of LCS calculations

From time-resolved CFD simulations, the time-dependent velocity field
$\mathbf{u}(\mathbf{x},t)$ is known in space and time. From this information
the fluid parcel or passive particle trajectories

$\displaystyle\mathbf{x}(\mathbf{x_{0}},t)=\mathbf{x}_{{0}}+\int_{t_{0}}^{t}\mathbf{u}(\mathbf{x}(\tau),\tau)d\tau$
(1)

can be calculated, where $\mathbf{x}_{0}$ is the starting point of a
trajectory in 3D space at a starting time $t_{0}$. Note, that each trajectory
is now labelled by its start location in space and time. If a set of initially
close passive particles is released at the same time the distances between
them change over time due to the fluid motion. Passive particles initially
forming a tiny sphere will undergo a linear deformation towards an ellipse for
short times as would occur in a solid body under stress before it breaks.
Certainly, in a fluid, the deformation will progress, and non-linear higher-
order terms will play a role in causing stretching and folding which is
crucial for mixing. However, as a first approximation and for short times
these higher-order terms are neglected for the analysis of the deformation. If
we consider infinitesimal spheres of initially close particles around all mesh
cell centres of our simulation starting at the same initial time $t_{0}$, we
obtain a set of different ellipsoids. All these ellipsoids have differently
stretched and contracted principal axes which point in different directions at
a slightly later time $t_{1}$. The principal axes of each ellipse denominate
the final directions of maximal stretching (major axis) and maximal
contraction (minor axis) of the initially spherical particle blob. The
stretching factor $S$ is the length of the major axis of the final ellipse
divided by the initial radius of the sphere. If this stretching factor at each
initial grid point is plotted, a 3D stretching field results revealing the
regions at which stretching and thus particle separation for the time interval
of interest $[t_{0},t_{1}]$ is largest due to the local flow conditions.
Normally, the scaled logarithm of this stretching factor, defined by

$\displaystyle\sigma_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)=\frac{1}{\left|t_{1}-t_{0}\right|}\log(S)\;,$
(2)

is plotted. This scaled logarithmic stretching factor is called the Finite-
Time Lyapunov Exponent [11]. Connected areas or lines of large FTLE values
characterise the fluid transport as these denote the areas or lines along
which deformation and thus particle separation is largest. All these
geometrical considerations have their mathematical counterparts. The
stretching factor as described is the square root of the maximal eigenvalue of
the right Cauchy-Green deformation tensor $\mathbf{C}_{t_{0}}^{t_{1}}$. This
tensor can be calculated for every mesh cell as envisioned above for the
ellipsoid. As its name reveals it includes all the information about the
deformation of the fluid masses at this point for the short time interval
$t_{1}-t_{0}$, and notably it is an objective tensor such that high stretching
values and candidates for LCS derived from it will persist regardless of the
motion of the observer (invariant to a time-dependent translation and rotation
of the coordinate system of the observer) [12, 11].

The governing ordinary differential equation (ODE) for the evolution of a
fluid parcel or a passive particle reads

$\displaystyle\dot{\mathbf{x}}=\mathbf{u}(\mathbf{x}(t),t)\;.$ (3)

Therefore, the infinitesimal separation
$\mathbf{\gamma}=\mathbf{x}-\mathbf{x}^{*}$ of the passive particle, imagined
in the centre of a infinitesimal sphere, to a particle on the surface of the
sphere will be governed by the ODE

$\displaystyle\mathbf{\delta\dot{\mathbf{x}}}=\nabla\mathbf{u}\mathbf{\gamma}\;.$
(4)

The solution of this ODE is an exponential function, which explains why the
FTLE is defined as the logarithm of the stretching factor.

To analyse the stretching during short but finite time intervals, particles
distributed on a mesh are advected with the flow from an initial time $t_{0}$
over the time interval $T=\left|t_{1}-t_{0}\right|$ to $t_{1}$. From the
integral version of the governing ODE (Eq. 3) we obtain the definition of the
flow map, $\Phi_{t_{0}}^{t_{1}}$, which maps all the particles from their
initial positions onto their final positions at time $t_{1}$, viz.

$\Phi_{t_{0}}^{t_{1}}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n};\quad\mathbf{x}_{0}\mapsto\mathbf{x}_{0}+\int_{t_{0}}^{t_{1}}\mathbf{u}(\mathbf{x}(\tau),\tau)d\tau\;.$
(5)

To obtain the separation of two initially close particles after this time
interval a Taylor series

$\displaystyle\delta\mathbf{x}(t_{1})=\Phi_{t_{0}}^{t_{1}}(\mathbf{x}_{0})-\Phi_{t_{0}}^{t_{1}}(\mathbf{x_{0}+\delta
x(t_{0})})=\mathbf{D}\boldsymbol{\Phi}_{t_{0}}^{t_{1}}(\mathbf{x}_{0},t_{0})\delta
x(t_{0})+\mathcal{O}(|\delta x(t_{0})^{2}|)$ (6)

around the initial position can be employed. Where
$\mathbf{D}\boldsymbol{\Phi}_{t_{0}}^{t_{1}}(\mathbf{x}_{0},t_{0})$ is the
gradient (Jacobian) of the flow map with regard to the initial separation and
is also the normalised fundamental matrix solution of the equation of
variations above (Eq. 4) [12]. Therefore, the particle separation at time
$t_{1}$ can be written as

$\displaystyle\left\|\delta\mathbf{x}\left(t_{1}\right)\right\|=\sqrt{\left\langle\delta\mathbf{x}(t_{0}),\left[\mathbf{D}\boldsymbol{\Phi}_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)\right]^{*}\left[\mathbf{D}\boldsymbol{\Phi}_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)\right]\delta\mathbf{x}(t_{0})\right\rangle}.$
(7)

The right Cauchy-Green deformation tensor is then defined as

$\displaystyle\mathbf{C}_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)=\left[\mathbf{D}\Phi_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)\right]^{*}\left[\mathbf{D}\Phi_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)\right]\;.$
(8)

In this way the Finite-Time Lyapunov Exponent $\sigma_{t_{0}}^{t_{1}}$ for the
time interval $t_{0}$ to $t_{1}$ can now be defined on the basis of this
tensor in a more thorough, mathematical way. Therefore, it is now defined by

$\displaystyle\sigma_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)=\frac{1}{\left|t_{1}-t_{0}\right|}\log\sqrt{\lambda_{\operatorname{max}}\left(\mathbf{C}_{t_{0}}^{t_{1}}\left(\mathbf{x}_{0},t_{0}\right)\right)}\;.$
(9)

Here $\lambda_{\operatorname{max}}$ is the maximum eigenvalue of the right
Cauchy-Green deformation tensor and can be calculated using standard solvers.
In the picture of the small ellipsoid, the square root of the eigenvalue is
just the above stretching rate $S$.

##  3 Computational details

The computation of flow maps within libcfd2lcs is described thoroughly in [9].
The following section presents a brief overview of how the computation is done
in practice and which different timescales play a role in the calculations.
Hereafter, we describe the structure and functionality of the newly developed
function object. We will focus on how the function object acts as an interface
between OpenFOAM and libcfd2lcs, how parallelisation is ensured and what has
to be considered for the output of the generated data.

###  3.1 Numerical flow map computation in libcfd2lcs

libcfd2lcs is able to calculate both forward-time and backward-time FTLE
fields. However, it uses two very different approaches for calculating the
respective flow maps. The general approach used for the computation of the
forward time flow-map $\Phi_{t_{0}}^{t_{0}+T}$ and the resulting forward-time
FTLE field is very straightforward. A set of tracer particles is initialised
on a grid with spacing $\Delta x_{lcs}$ by setting each initial tracer
coordinate to the cell centre coordinate of a corresponding mesh cell. Then
the flow map at each cell centre is computed by passively advecting these
tracers with the flow, which mathematically corresponds to an integration of
equation

$\displaystyle\frac{d\mathbf{x}}{dt}=\mathbf{u}(\mathbf{x},t)$ (10)

over the time interval $T$. Numerically this integration is done by utilising
Runge-Kutta methods, with step size $\Delta t_{lcs}$. The time and space
dependent velocity field $\mathbf{u}(\mathbf{x},t)$ results from the specific
fluid simulation under consideration and is passed to libcfd2lcs after each
simulation time step $\Delta t_{sim}$ (see Section 3.2). In order to save the
flow map $\Phi_{t_{0}}^{t_{0}+T}$, the location of each particle after the
integration is stored at its initial position.

As the evaluation of FTLE fields, indicating LCS candidates, is mainly
relevant for time-dependent flows, it is often important to animate their
evolution. At first glance, this would mean that a sequence of large particle
sets would have to be integrated, requiring a great amount of computation.
This problem is solved using a method developed by Brunton and Rowley [13].
With this method a flow map of the interval $T$ can be constructed from a
sequence of $k$ flow maps over a smaller interval $h$, where $T=kh$. Following
the notation of [13] this can be expressed as

$\Phi_{t_{0}}^{t_{0}+T}=\Phi_{t_{0}+(k-1)h}^{t_{0}+kh}\circ\cdots\circ\Phi_{t_{0}+h}^{t_{0}+2h}\circ\Phi_{t_{0}}^{t_{0}+h}\;.$
(11)

In practical terms, this means that the particle grid is reinitialised for
every new time interval $h$ after which they are advected again with the flow.
Then the sub-step flow map is stored and the complete flow map is constructed
when all needed sub-step flow maps are available. It is important to note that
since a discrete particle grid is used for the sub-step flow map computation,
interpolation of the sub-step flow maps is needed in order to match the
trajectories at different timelevels when reconstructing the flow map
$\Phi_{t_{0}}^{t_{0}+T}$ (see [13] for more details).

A different approach for constructing the backward-time flow maps is used.
This is due to the fact that using the Lagrangian approach would require to
store all computed velocity fields in the subset interval $h$ before the
integration of the tracers from $t_{0}+h$ to $t_{0}$ could be done backward in
time. Although this already includes Brunton’s and Rowley’s method for the
flow map construction, the Lagrangian approach would be ”cumbersome and
resource intensive” [9, p. 4]. Therefore, libcfd2lcs uses an Eulerian approach
for the flow map computation proposed by Leung [14]. In contrast to the
forward-time flow map, the backward-time flow map $\Phi_{t_{0}+T}^{t_{0}}$
describes for each grid point where a particle, that is at that point at time
$t_{0}+T$, originally was at time $t_{0}$. With Leung’s Eulerian approach this
backward-time flow map at time $t_{0}+T$ is computed by initialising a vector
field $\mathbf{\Psi}(\mathbf{x},t_{0})$ on a grid with the cell centre
coordinates at time $t_{0}$. The advection of this so called ”takeoff
coordinate field” in an Eulerian reference frame is then described by the
level set equation

$\frac{\partial\mathbf{\Psi}(\mathbf{x},t)}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{\Psi}(\mathbf{x},t)=\mathbf{0}\;.$ (12)

Solving this equation over the time Interval $[t_{0},t_{0}+T]$ in forward time
gives $\mathbf{\Psi}(\mathbf{x},t_{0}+T)$, which represents the takeoff
coordinates of a Lagrangian particle at $t_{0}$ reaching $\mathbf{x}$ at time
$t_{0}+T$. Thus, the backward-time flow map $\Phi_{t_{0}+T}^{t_{0}}$ is
equivalent to $\mathbf{\Psi}(\mathbf{x},t_{0}+T)$. libcfd2lcs solves equation
(12) by using a semi-Lagrangian advection approach with the time step size

$\Delta t_{lcs}=\frac{c_{cfl}\,\Delta x_{lcs}}{\mathbf{u}(\mathbf{x},t)}$ (13)

of this procedure being restricted by the CFL condition $c_{cfl}<1$ (see [9]
and [14] for more details). Furthermore, Brunton’s and Rowley’s flow map
construction method is also applied to the backward-time flow maps computed
with the Eulerian method. Hence, the takeoff coordinate field is reinitialised
after every sub-step time interval $h$ and the backward-time flow map

$\Phi_{t_{0}+T}^{t_{0}}=\Phi_{t_{0}+h}^{t_{0}}\circ\Phi_{t_{0}+2h}^{t_{0}+h}\circ\cdots\circ\Phi_{t_{0}+kh}^{t_{0}+(k-1)h}$
(14)

is constructed form $k$ sub-step backward-time flow maps.

Figure 1: Different timescales relevant for the LCS evaluation

Since a lot of different timescales are relevant in the practical FTLE field
computation described above, we try to differentiate and order them in the
following, before describing the structure and functionality of the newly
developed function object in the next section. The basis of the on-the-fly LCS
evaluation is a parallel running simulation that provides the velocity fields.
Here, three intervals are of interest (see Fig. 1): the overall simulation
time that spans from the simulation start time $t_{sim\\_start}$ to the
simulation end time $t_{sim\\_end}$, the time step size of the simulation
$\Delta t_{sim}$ and the write time interval of the simulation results $\Delta
t_{sim\\_write}$. The computed velocity fields represent a fluid flow for
which a reference timescale $\Delta t_{ref}$ can be identified. This reference
timescale characterises the dominant hydrodynamic timescale of the flow and is
typically larger than the simulation time step size. In order to save
computing resources the LCS evaluation of the simulated flow does not
necessarily have to start and end at the same time as the simulation.
Therefore, a separate start and end time for the LCS evaluation denoted as
$t_{lcs\\_start}$ and $t_{lcs\\_end}$ can be defined (see Fig. 1). During the
LCS evaluation, a series of FTLE fields are computed. These FTLE fields are
calculated from time $T$ flow maps, which themselves are calculated as
described earlier in this section. This means storing and constructing the
time $T$ flow maps from multiple sub-step flow maps after each LCS sub-step
integration interval $h$. Calculating the sub-step flow maps in turn requires
to numerically solve the equations (10) or (12) using the finite time step
$\Delta t_{lcs}$. While $\Delta t_{lcs}$ is set automatically according to
equation (13) and a specified CFL number, $T$ and $h$ have to be defined by
the user. In order to detect all LCS candidates, $T$ is usually chosen to be
larger than $\Delta t_{ref}$ of the investigated flow [9]. With the aim of
animating the evolution of the FTLE field, $h$ is typically set significantly
smaller than $\Delta t_{ref}$ while being in the order of magnitude of $\Delta
t_{sim\\_write}$.

###  3.2 Structure and functionality of the function object

In general, function objects can be used to generate additional data at
runtime of the simulation. In doing so, function objects can access data
generated by the flow solver at runtime, which offers a great advantage over
classical post-processing since it can only utilise the stored fields or
logged information. The newly developed function object incorporates the
functionalities of libcfd2lcs into OpenFOAM at runtime while acting as an
interface between both. This is achieved by processing the data generated by
OpenFOAM and the subsequent exchange of this data via the libcfd2lcs API (see
[8] for a detailed description of the libcfd2lcs API).

The calculation of the flow maps, the calculation of the resulting FTLE fields
and the subsequent saving of these fields is completely handled by libcfd2lcs.
The basic task of the function object is to pass the cell centre position
vectors of the computational grid as well as the velocity field calculated by
OpenFOAM to libcfd2lcs. Due to the very strict data structure requirements of
libcfd2lcs this is not a trivial task. libcfd2lcs can only use static
rectlinear grids for the calculation of forward-time and backward-time flow
maps and therefore needs the velocity fields on these grids. This means that
the mesh and velocity data has to be globally organised in an $(i,\,j,\,k)$
structured format [8, p. 6]. Since the LCS evaluation should also be available
for simulations on moving grids with general topology and adaptive grid
refinement, the function object offers several possibilities to deal with this
problem.

In the simplest case, where the simulation mesh is already a static rectlinear
mesh, the function object does not need to process the grid and velocity data,
but can directly transfer it to libcfd2lcs as basic C++ arrays. This is the
preferred method when the flow domain can be represented by a static
rectlinear mesh and e.g. immersed boundary methods are used. If a moving mesh,
a mesh of general topology or adaptive mesh refinement is used for the
simulation a different approach is needed in order to prepare the data for its
use in libcfd2lcs. Here, an additional static rectlinear mesh needs to be
constructed in the preprocessing step, which can be done e.g. by using the
blockMesh utility. This mesh has to contain the region for which the LCS
diagnostic should be performed, meaning that it can cover the whole simulation
domain as well as only a part of it. However, since libcfd2lcs also requires
boundary conditions for the FTLE field calculations, the boundary patches of
the additional LCS mesh must be set accordingly. The user can choose between
empty, symmetryPlane, wedge, cyclic and the generic patch patch types which
the function objects translates into the corresponding libcfd2lcs boundary
types. Then, during runtime, the velocity fields are mapped from the
simulation mesh of general topology to the static rectlinear LCS mesh, from
which the data can again be transferred to libcfd2lcs as basic C++ arrays.
Although this implies that interpolation errors are made during the mapping
process, the LCS evaluation is hardly affected by this. Haller showed in [15]
that LCS are very robust against errors in the velocity field. Also, the
additional computational overhead due to the mapping can be neglected compared
to the overhead caused by the flow map computations. The function object also
implements a third approach in which no additional LCS mesh is needed. This
approach utilises the ability to construct complex, moving mesh geometries out
of simple unconnected mesh regions in OpenFOAM with the oversetMesh approach.
Using this approach the function object can utilise any specified static
rectlinear mesh region of the oversetMesh for the LCS evaluation, meaning that
the background mesh as well as any other static rectlinear mesh region can be
used. In doing so, the function object extracts the mesh and velocity data
from the specified mesh region of the oversetMesh and passes it to libcfd2lcs
analogously to the previous approaches. Here the overset type patches are
generally passed on as inlet or outlet, as they are treated the same by
libcfd2lcs.

As libcfd2lcs also uses the domain decomposition approach and MPI for the
parallelisation of the computations, the integration within the
parallelisation of OpenFOAM is done in a straightforward manner. The local
subdomains of the rectlinear LCS mesh and its velocity data are passed to
libcfd2lcs together with an offset, which describes the position of the cell
data in the globally $(i,\,j,\,k)$ structured data array (see Fig. 2). For the
MPI communication, the same MPI communicator as used for OpenFOAM is shared
with libcfd2lcs. Therefore, the function object can be used for simulations
running in parallel or serial. However, if the approach involving an
additional LCS mesh is used, special attention is required for the domain
decomposition in the preprocessing step. Here the simulation mesh, as well as
the LCS mesh, must be cut along the same surfaces to make sure that the
mapping of the velocity fields from one mesh to the other works properly.

Figure 2: Subdomain offset of a $i,\,j,\,(k=0)$ structured rectlinear mesh

As already mentioned, the data output of the flow-map and FTLE field data is
completely handled by libcfd2lcs. This is due to the fact that the data output
interval defined by $h$ can differ from the solver write interval $\Delta
t_{sim\\_write}$ (see section 3.1). Therefore, the results generated by the
function object are not stored in corresponding time directories but in a
separate folder in the case directory called cfd2lcs_output. Additionally, a
directory named cfd2lcs_temp is created inside of which all the sub-step data
is stored. All data is stored in the Tecplot ASCII data file format (*.dat)
and therefore can be visualised in ParaView when opened with its internal
Tecplot reader or other common visualisation programs. In addition to this
data, the computational overhead generated by the use of the function object
with respect to the actual simulation is also output in the solver log file
after each simulation time step. This enables the user to examine the
computational costs of the LCS evaluation.

##  4 Examples of usage

In this section a few examples are presented which are designed to show the
functionality and capabilities of the function object. Therefore, example
cases are presented in which only a rectlinear simulation mesh, a separate
simulation and LCS mesh and a single oversetMesh are used.

###  4.1 Steady ABC flow

The Arnold-Beltrami-Childress (ABC) flow is an exact periodic solution of the
Euler equations and is often used in the literature to verify LCS calculation
methods. Therefore this case is also being reviewed here. The velocity field

$\displaystyle\mathbf{u}=\nabla\times[-\Psi\mathbf{k}+\nabla\times(\Phi\mathbf{k})]$
(15)

of the ABC flow can be described using 2 scalar potentials $\Psi$ and $\Phi$
[16] which themselves are defined as

$\displaystyle\begin{aligned} &\Psi=-[C\sin(y)+B\cos(x)]\\\
&\Phi=A[-x\cos(z)+y\sin(z)]-\Psi\;.\end{aligned}$ (16)

In (15), $\mathbf{k}$ can be any unit vector but is commonly chosen to be the
vertical unit vector. This leads to the three expressions of the velocity
components

$\displaystyle\begin{aligned} &u=A\sin(z)+C\cos(y)\\\ &v=B\sin(x)+A\cos(z)\\\
&w=C\sin(y)+B\cos(x)\;.\end{aligned}$ (17)

(a) Forward-time FTLE field for the whole computational domain

(b) Backward-time FTLE field for the whole computational domain

(c) Forward-time field for a cross-section at $z=2\pi$

(d) Backward-time field for a cross-section at $z=2\pi$

Figure 3: FTLE fields of a stationary ABC flow

The parameters $A$, $B$ and $C$ can be freely selected and influence the
properties of the ABC flow. In order to create comparability with literature
values, $A=0.5$, $B=0.8$, $C=0.8$ is chosen. In order to test the newly
developed function object on this flow configuration a dedicated ABC flow
OpenFOAM solver was written. This solver does not solve the Euler equations in
the usual sense, but sets the velocity components on a given computational
mesh according to (17). Due to the periodicity of the flow solution, the
dimensions of the computational mesh used in this case setup are specified as
$x,y,z\in[0,2\pi]$ with a mesh size of $100\times 100\times 100$. Since the
described mesh is rectlinear no additional LCS mesh is used. Again for reasons
of comparability, a LCS integration time of $T=10\mathrm{s}$ is selected for
the LCS evaluation. The results of the LCS evaluation, both in forward- and
backward-time, can be seen in Figure 3.

In these results the FTLE ridges, which indicate the LCS candidates in the ABC
flow, can be seen very clearly. Furthermore, the results agree very well with
the results from [16], both qualitatively and quantitatively, which suggests
that the new function object calculates the FTLE ridges reliably.

###  4.2 Time dependent double gyre

Another frequently used flow for the verification of LCS computing algorithms
is the time periodic Rayleigh-Bénard convection flow, or often called double
gyre, proposed by Solomon and Gollub [17]. The velocity field of this flow can
be describe by using a stream function $\psi$

$\displaystyle\begin{aligned} &u=-\frac{\partial\psi}{\partial y}\\\
&v=\frac{\partial\psi}{\partial x}\;.\end{aligned}$ (18)

Here $\psi$ is defined by

$\displaystyle\psi(x,y,t)=A\sin(\pi f(x,t))\sin(\pi y)$ (19)

with

$\displaystyle\begin{aligned} &f(x,t)=a(t)x^{2}+b(t)x\\\
&a(t)=\epsilon\sin(\omega t)\\\ &b(t)=1-2\epsilon\sin(\omega t)\end{aligned}$
(20)

This leads to the expressions for two-dimensional velocity components

$\displaystyle\begin{aligned} &u=-\pi A\sin(\pi f(x))\cos(\pi y)\\\ &v=\pi
A\cos(\pi f(x))\sin(\pi y)\frac{\mathrm{d}f}{\mathrm{~{}d}x}\;.\end{aligned}$
(21)

As the name double gyre suggests, this model defines the flow of two two-
dimensional gyres enclosed in a rectangle which expand and contract
periodically along the x-axis. Therefore, the periodic motion is controlled by
$\epsilon$ if $\epsilon\neq 0$. Then $\epsilon$ describes approximately how
far the line separating the gyres moves to the left or right from its centre
position [11]. Otherwise ($\epsilon=0$), no periodic motion is happening.
Furthermore, $A$ specifies the magnitude of the velocity vectors and
$\omega/2\pi$ determines the oscillation frequency of the gyres.

(a) Forward-time FTLE field

(b) Backward-time FTLE field

Figure 4: FTLE fields of a time-dependent double gyre at $t=15\mathrm{s}$.

Similar to the ABC flow example, a dedicated OpenFOAM solver was written for
this case, which sets the velocity field on a given computational mesh
according to (21). For comparability, a mesh with the same specifications as
in [18],[14] and [11] was used. It has the dimensions
${[0,2]\times[0,1]\times[0,0.1]}m$ and a resolution of $512\times 256\times 1$
cells. As this mesh is also static and rectlinear no additional LCS mesh was
used. For the mathematical model of the flow the parameter values are chosen
to be $\epsilon=0.1$, $A=0.1\mathrm{m}\,\mathrm{s}^{-1}$ and
$\omega=2\pi/10\mathrm{s}$. Since the oscillation frequency is known, the
hydrodynamic time scale can be easily determined by
$t_{ref}=2\pi/\omega=10\mathrm{s}$. As described in section 3.1, the LCS
integration time interval $T$ should be set larger than $t_{ref}$. Therefore,
it is set to $T=1.5\cdot t_{ref}=15\mathrm{s}$. Figure 4 shows the forward-
and backward-time FTLE fields at $t=15\mathrm{s}$ of the previously described
double gyre flow. Again, the results match very well with the results from
[18],[14] and [11]. This confirms that the function object is able to
calculate the correct FTLE fields from velocity fields generated by OpenFOAM.

###  4.3 Flow around cylinder

As it has already been shown in the previous examples that the function object
can calculate the correct FTLE fields from velocity fields provided by
OpenFOAM, this example will focus on how to deal with non-rectlinear
simulation meshes. For this purpose, a standard flow problem is selected that
is very well suited for an LCS evaluation: the flow around an infinitely long
cylinder. The general case setup contains a fluid domain with size
${[-20,30]\times[-20,20]\times[-0.5,0.5]}\mathrm{m}$ that surrounds a cylinder
with diameter $D=2\mathrm{m}$ and its centre axis at $x=y=0\mathrm{m}$. The
free-stream velocity and the fluids kinematic viscosity are set to
$\mathbf{u}^{\sf T}=(1\;0\;0)\mathrm{m}\,\mathrm{s}^{-1}$ and
$\nu=0.01\mathrm{m}^{2}\mathrm{s}^{-1}$, respectively. This results in a
Reynolds number of $\operatorname{Re}=200$ which indicates that vortex
shedding behind the cylinder occurs in a barely laminar regime. If we also
assume a Strouhal number of $\operatorname{St}=0.2$ at
$\operatorname{Re}=200$, the hydrodynamic time scale of this flow is
$t_{ref}=D/(u\cdot\operatorname{St})=10\mathrm{s}$.

Because of the cylinder in the middle of the domain, a computational mesh
discretising this domain is no longer rectlinear. Therefore, we consider two
different procedures in the LCS evaluation, the first of which is carried out
in two different ways.

(a) Non rectlinear simulation mesh

(b) LCS mesh enclosing the whole domain

(c) Forward-time FTLE field

(d) Backward-time FTLE field

Figure 5: Simulation mesh, LCS mesh and FTLE fields of a flow around a
cylinder with $\operatorname{Re}=200$ at $t=118\mathrm{s}$.

Starting with the procedure where an additional rectlinear computational mesh
is used for the LCS evaluation, the flow domain is discretised with a
simulation mesh consisting of 9200 hexahedra (see upper left mesh in Fig. 5).
The flow solver that is used to simulate the previously described flow from
$t=0\mathrm{s}$ to $t=120\mathrm{s}$ is pimpleFoam with the initial conditions
being calculated by potentialFoam. The first additional LCS mesh that is used
within this procedure encloses the whole flow domain (see upper right mesh in
Fig. 5). In order to minimise the loss of information during the mapping of
the velocity fields between the two grids, the resolution of the LCS mesh is
chosen in a way that it corresponds approximately to the finest resolution in
the simulation mesh. This leads to a LCS mesh with $200\times 160\times 1$
hexahedra. The boundary patch types are set to patch for the left and right
patch (inlet,outlet), to symmetryPlane for the bottom and top patch and to
empty for the front and back patch. The LCS integration time $T$ is again
based on $t_{ref}$ and is set to $T=1.5\cdot t_{ref}=15\mathrm{s}$. For a good
animation of the dynamics of the FTLE fields $h$ is chosen to be
$h=T/10=1.5\mathrm{s}$. The results of the forward- and backward-time FTLE
fields can be seen in Fig. 5. They show how the vortices behind the cylinder
form large coherent structures, where the FTLE ridges of the backward-time
FTLE fields separate different fluid packages that do not mix in the vortex
street.

(a) Non rectlinear simulation mesh

(b) LCS mesh enclosing the region of interest

(c) Forward-time FTLE field

(d) Backward-time FTLE field

Figure 6: Simulation mesh, LCS mesh and FTLE fields of a flow around a
cylinder with $\operatorname{Re}=200$ at $t=118\mathrm{s}$.

Since the FTLE ridges only appear in a fraction of the overall domain and the
LCS evaluation is computation-wise a quite costly operation, a second LCS mesh
is prepared. This second LCS mesh is a lot smaller than the first one and
encloses only the fraction of the flow domain where the FTLE ridges are
expected to show up (see Fig. 6). The boundary patches on the smaller LCS mesh
and its spacial resolution are also set analogous to its bigger counterpart,
leading to a LCS mesh of size $[-13,27]\times[-7.5,7.5]\times[-0.5,0.5]$
containing $160\times 60\times 1$ hexahedra. Repeating the computations with
the use of the smaller LCS mesh gives the results which are displayed in Fig.
6 and are found to match with the results from the bigger LCS mesh. This shows
that the LCS evaluation, when done with a separate LCS mesh, can be used in a
very targeted way. The advantages this brings in terms of computational costs
are discussed after considering the second procedure for the LCS evaluation of
this flow problem.

(a) oversetMesh containing a rectlinear background mesh zone, a rectlinear
refinement zone which is also used for the LCS evaluation, and a cylindrical
mesh zone containing the cylinder.

(b) Forward-time FTLE field

(c) Backward-time FTLE field

Figure 7: oversetMesh and FTLE fields of a flow around a cylinder with
$\operatorname{Re}=200$ at $t=118\,\mathrm{s}$.

The second procedure, which can be used on problems where no single static
rectlinear mesh can be constructed, utilises OpenFOAM’s oversetMesh
functionalities. With regard to the flow problem considered here, an
oversetMesh is constructed with the same dimensions as the simulation mesh
used previously. It consists of three mesh zones, namely a rectlinear
background mesh zone that spans the whole fluid domain, another finer and
smaller mesh zone that is used for a finer resolution of the flow and a
cylindrical mesh that surrounds the cylinder (see Fig. 7). For comparability
reasons the finer rectlinear mesh zone has the same dimensions and resolution
as the smaller additional LCS mesh considered previously and is therefore
specified as the cell zone for the LCS evaluation. Also all other LCS
evaluation settings are adopted. The only difference to the previously
considered simulations is the used flow solver. Here the flow solver is
pimpleDyMOversetFoam due to the used oversetMesh. The resulting forward- and
backward-time FTLE fields of this simulation can be found in Fig. 7. They
match with the results from the previously considered procedure which shows
that both approaches can be used equally well. The only thing that stands out
are the high FTLE values along some boundaries in the studied solutions. These
occur because of the way libcfd2lcs handles its inlet and outlet boundary
conditions. It fixes out-flowing Lagrangian particles/takeoff coordinates on
”open” boundaries and cannot generate new in-flowing particles during the flow
map computation. Therefore, high FTLE values occur in the forward-time FTLE
fields at ”open” boundaries where inflow occurs, since there the most
”stretching” happens. Vice versa, high FTLE values occur in the backward-time
FTLE fields at ”open” boundaries where outflow occurs, since there the most
”folding” happens. These high values at ”open” boundaries are just artefacts
and have to be neglected. The reason they appear more in the oversetMesh
approach is that all overset type patches are passed to libcfd2lcs as ”open”
boundaries whereas the user can specify all patches problem dependent in the
additional LCS mesh approach.

Looking at the computation times of the flow calculations including the LCS
evaluation, it becomes evident that LCS evaluation is a very costly operation
(see Tab. 1). When using the ”large” additional LCS mesh the simulation takes
approximately 30 time longer than without the LCS evaluation. This can be
improved by using the smaller additional LCS mesh. Here the simulation takes 9
times longer than without the LCS evaluation. Since the costs for the LCS
evaluations are almost independent of the underlying simulation for a constant
grid size, this factor becomes smaller and smaller for more complex
simulations. This can also be seen from the fact that the factor is only 2.5
when the oversetMesh approach is used because the computations of the pressure
and velocity fields take longer on an oversetMesh. At this point, however, it
must be emphasised that the flow considered here is not a highly complex
problem, which can also be seen from the simulation time of 1.5 min on a
normal mesh and 8 min on an oversetMesh.

| ”large” add. LCS mesh | ”small” add. LCS mesh | oversetMesh
---|---|---|---
Total time | 44 | 13.5 | 20
Simulation time | 1.5 | 1.5 | 8
Total FTLE comp. time | 42.5 | 12 | 12
Fwd.-FTLE comp. time | 8.5 | 2 | 2
Bkwd.-FTLE comp. time | 21 | 6 | 6
I/O, flow map const., … | 13 | 4 | 4
Table 1: Computing times in minutes for the flow simulation including LCS
evaluation on 2 cores of a Intel Core i5 processor

##  5 Summary & Conclusion

We provide an OpenFOAM function object based on libcfd2lcs to compute Finite-
Time Lyapunov Exponent (FTLE) fields that indicate candidates of Lagrangian
Coherent Structures (LCS) and allow to visualise finite-time stretching and
folding fields. LCS reveal the robust skeleton of material surfaces and are
key to quantitatively assess material transport in time-dependent flows. This
enables the OpenFOAM community to assess the geometry of the material
transport in any flow quantitatively on-the-fly using principally any OpenFOAM
flow solver.

Focusing on the practical aspects, we only give a brief overview of the
mathematical foundation as well as how the computation is done in practice. We
describe the structure and functionality of the newly developed function
object. Further focus is laid on how the function object acts as an interface
between OpenFOAM and libcfd2lcs, how parallelisation is ensured and what has
to be considered for the output of the generated data.

From validation of the presented function object using simple benchmark
problems, a notable computational overhead has been recognised. However, if
LCS evaluations are used for much more complex problems as the ones used here,
the computational overhead significantly drops and the LCS evaluation no
longer accounts for the largest proportion of the computation time.
Nevertheless, the user should be aware that the calculation of FTLE fields is
expensive and should therefore think carefully about the size and position of
the LCS mesh. In addition, consideration should also be given to whether both
forward and backward-time FTLE calculations are required or if one of them is
sufficient.

## References

  * [1] A. Kameke, S. Kastens, S. Rüttinger, S. Herres-Pawlis, and M. Schlüter, “How coherent structures dominate the residence time in a bubble wake: An experimental example,” _Chemical Engineering Science_ , vol. 207, pp. 317–326, Nov. 2019. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0009250919305366
  * [2] F. Huhn, A. von Kameke, V. Pérez-Muñuzuri, M. J. Olascoaga, and F. J. Beron-Vera, “The impact of advective transport by the South Indian Ocean Countercurrent on the Madagascar plankton bloom: Jet in Madagascar Plankton Bloom,” _Geophysical Research Letters_ , vol. 39, no. 6, pp. n/a–n/a, Mar. 2012. [Online]. Available: http://doi.wiley.com/10.1029/2012GL051246
  * [3] F. d’Ovidio, V. Fernández, E. Hernández-García, and C. López, “Mixing structures in the Mediterranean Sea from finite-size Lyapunov exponents: Mixing Structures in the Mediterranean Sea,” _Geophysical Research Letters_ , vol. 31, no. 17, pp. n/a–n/a, Sep. 2004. [Online]. Available: http://doi.wiley.com/10.1029/2004GL020328
  * [4] J. B. Weiss and A. Provenzale, Eds., _Transport and Mixing in Geophysical Flows: Creators of Modern Physics_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008. [Online]. Available: http://link.springer.com/10.1007/978-3-540-75215-8
  * [5] Z. Neufeld and E. Hernández-García, _Chemical and Biological Processes in Fluid Flows: A Dynamical Systems Approach_. Imperial College Press, Sep. 2009. [Online]. Available: https://www.worldscientific.com/worldscibooks/10.1142/p471
  * [6] S. Balasuriya, N. T. Ouellette, and I. I. Rypina, “Generalized Lagrangian coherent structures,” _Physica D: Nonlinear Phenomena_ , vol. 372, pp. 31–51, Jun. 2018. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0167278917302750
  * [7] A. Hadjighasem, M. Farazmand, D. Blazevski, G. Froyland, and G. Haller, “A critical comparison of Lagrangian methods for coherent structure detection,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 27, no. 5, p. 053104, May 2017. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.4982720
  * [8] Justin Finn, Romain Watteaux, and Andrew Lawrie, “libcfd2lcs: A general purpose library for computation of Lagrangian coherent structures during CFD simulation,” Feb. 2017. [Online]. Available: https://www.archer.ac.uk/training/virtual/2016-06-22-libcfd2lcs/ArcherVT.pdf
  * [9] J. Finn and S. V. Apte, “Integrated computation of finite-time Lyapunov exponent fields during direct numerical simulation of unsteady flows,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 23, no. 1, p. 013145, Mar. 2013. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.4795749
  * [10] M. Farazmand and G. Haller, “Computing Lagrangian coherent structures from their variational theory,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 22, no. 1, p. 013128, Mar. 2012. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.3690153
  * [11] S. C. Shadden, F. Lekien, and J. E. Marsden, “Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows,” _Physica D: Nonlinear Phenomena_ , vol. 212, no. 3-4, pp. 271–304, Dec. 2005. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0167278905004446
  * [12] G. Haller, “Lagrangian Coherent Structures,” _Annual Review of Fluid Mechanics_ , vol. 47, no. 1, pp. 137–162, Jan. 2015. [Online]. Available: https://www.annualreviews.org/doi/10.1146/annurev-fluid-010313-141322
  * [13] S. L. Brunton and C. W. Rowley, “Fast computation of finite-time Lyapunov exponent fields for unsteady flows,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 20, no. 1, p. 017503, Mar. 2010. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.3270044
  * [14] S. Leung, “An Eulerian approach for computing the finite time Lyapunov exponent,” _Journal of Computational Physics_ , vol. 230, no. 9, pp. 3500–3524, May 2011. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0021999111000799
  * [15] G. Haller, “Lagrangian coherent structures from approximate velocity data,” _Physics of Fluids_ , vol. 14, no. 6, pp. 1851–1861, Jun. 2002. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.1477449
  * [16] M. H. Sulman, H. S. Huntley, B. Lipphardt, and A. Kirwan, “Leaving flatland: Diagnostics for Lagrangian coherent structures in three-dimensional flows,” _Physica D: Nonlinear Phenomena_ , vol. 258, pp. 77–92, Sep. 2013. [Online]. Available: https://linkinghub.elsevier.com/retrieve/pii/S0167278913001450
  * [17] T. Solomon and J. Gollub, “Chaotic Transport in Time-Dependent Rayleigh-Bénard Convection,” _Physical review. A_ , vol. 38, pp. 6280–6286, Jan. 1989.
  * [18] D. Lipinski and K. Mohseni, “A ridge tracking algorithm and error estimate for efficient computation of Lagrangian coherent structures,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 20, no. 1, p. 017504, Mar. 2010. [Online]. Available: http://aip.scitation.org/doi/10.1063/1.3270049


A versatile and transportable endstation for controlled molecule experiments

Wuwei Jin
Hubertus Bromberger
Lanhai He
Melby Johny
Ivo S. VinklárekKarol DługołęckiAndrey SamartsevFrancesca Calegari
Sebastian Trippel
Jochen Küpper

We report on a new versatile transportable endstation for controlled molecule (eCOMO) experiments
providing a combination of molecular beam purification by electrostatic deflection and
simultaneous ion and electron detection using velocity-map imaging (VMI). The $b$-type
electrostatic deflector provides spatial dispersion of species based on their
effective-dipole-moment-to-mass ratio. This enables selective investigation of molecular
rotational quantum states, conformers, and molecular clusters. Furthermore, the double-sided VMI
spectrometer equipped with two high-temporal-resolution event-driven Timepix3 cameras provides
detection of all generated ions independently of their mass-over-charge ratio and electrons. To
demonstrate the potential of this novel apparatus, we present experimental results from our
investigation of carbonyl sulfide (OCS) after ionization. Specifically, we provide the
characterization of the molecular beam, electrostatic deflector, and electron- and ion-VMI
spectrometer. The eCOMO endstation delivers a platform for ultrafast dynamics studies using a
wide range of light sources from table-top lasers to free-electron-laser and
synchrotron-radiation facilities. This makes it suitable for research activities spanning from
atomic, molecular, and cluster physics, over energy science and chemistry, to structural biology.

§ INTRODUCTION

The advancement of synchrotron and x-ray free-electron laser (XFEL) radiation sources, such as those
at DESY (FLASH), the linac coherent light source (LCLS), the free-electron laser radiation for
multidisciplinary investigations (FERMI), the SPring-8 Angstrom compact free-electron laser (SACLA),
the x-ray free-electron laser at the Paul Scherrer Institute (SwissFEL), the European x-ray
free-electron laser (EuXFEL), and the Shanghai soft x-ray free-electron laser (SXFEL), substantially
pushes the boundaries of existing research areas and breaks the ground for completely new
fields <cit.>. These sources generate ultrashort
x-ray pulses spanning from tens of femtoseconds down to the attosecond
range <cit.>.
Furthermore, they provide a high peak brilliance, coherence and photon energies ranging from 5 to
250 keV <cit.>. Table-top laser systems offer unique advantages in
other areas: These sources include variable-wavelength optical-parametric-amplifier (OPA)-based
systems <cit.>, high power lasers with repetition rates up
to 1 MHz using optical parametric chirped-pulse amplification (OPCPA) <cit.>, and lasers with attosecond pulse durations <cit.>. Each photon source has unique properties, functions, and parameters that
make it suitable for different investigations, ranging from ultrafast studies of chemical reactions
over bio-molecular dynamics to materials and life
science <cit.>. This progress led to particular
needs on user-side applications, notably atomic-level-resolution molecular dynamics experiments
based on ultra-cold molecular beam methods <cit.>.

Molecular beams are of significant importance in physical chemistry and molecular physics as they
provide a unique opportunity to obtain fundamental insights into molecular photophysics and chemical
processes <cit.>. Whereas substantial research was conducted on isolated small molecules
and atoms, investigating more complex molecular systems with molecular-level insight is still
restricted by many experimental challenges <cit.>. One of these challenges is to generate well-defined ensembles of the reactant
particles. In this regard, many molecular beams usually contain mixtures of isolated molecules,
various molecular clusters, or multiple isomers, , of chiral molecules or large biomolecules. The
ability to select a particular species plays a crucial role in gaining a deeper understanding of its
importance in chemical behavior that also impacts a wide range of biological
processes <cit.>. To overcome this obstacle, the electrostatic deflector was developed to a
achieve spatial separation of different neutral species in ultra-cold ($T_{\text{rot}}<1$ K)
molecular beams, ultimately obtaining pure particle ensembles <cit.>.

A second challenge is recording the molecular dynamics using experimentally accessible observables.
The reaction microscope (REMI) technique, for example, provides complete, correlated, 3D velocity
information of all particles detected in coincidence <cit.>. However, it is limited by low count rates, which can lead to poor
statistics. High-count-rate approaches such as VMI <cit.>, provide a significant
improvement in that respect. A significant drawback of VMI in high count rate conditions however is
that one cannot rely on coincidence detection to extract, , information about different electrons
and ions coming from the same target molecule, although covariance mapping
methods <cit.> can partially overcome this problem. Lastly, for VMIs,
photoelectron photoion coincidence (PEPICO) spectroscopy that is synchronizing VMI- and
time-of-flight (TOF) information, is also challenging <cit.>. These issues so far
limit access to the full kinematics of single reactions and, therefore, systematic studies of the
underlying physical mechanisms.

This article provides a comprehensive overview of the newly established endstation for controlled
molecule experiments. It addresses the aforementioned challenges of the current experimental
research. eCOMO provides pure, cold, dense, and pulsed molecular beams of selected species using
Even-Lavie valves <cit.> in combination with the electrostatic
deflector <cit.>. The endstation is equipped with a double-sided VMI spectrometer
combined with high spatio-temporal-resolution Timepix3 cameras combined with the software package
PymePix <cit.>. This lets us overcome the
limitations of standard camera technology for advanced imaging experiments with high event
rates <cit.>. Moreover, this also meets the high-frequency data acquisition requirements of
various light sources with high repetition rates mentioned earlier <cit.>.

To demonstrate the capabilities of the novel system, we present the results of its commissioning
with an optical laser in our laboratory. We used the fragmentation of OCS after strong-field
ionization as a benchmark to characterize the molecular beam, the electrostatic deflector, and the
electron-and-ion-VMI spectrometer.

§ APPARATUS DESCRIPTION

§.§ Vacuum system

The vacuum setup, shown schematically in the upper left corner of <ref>,
consists of four differentially pumped vacuum chambers housing the molecular beam source, the
electrostatic $b$-type deflector, the double-sided VMI spectrometer, and the molecular beam dump.

Sketch of the experimental setup consisting of a molecular beam valve, an electrostatic
$b$-type deflector, a double-sided VMI, and two Timepix3 cameras. The inset provides a
side view from a CAD model, where colored frames depict the different parts of the setup. All
distances are provided in cm. See text for details.

All chambers are mounted on rails to allow easy translation of the setup along the $z$-axis, which
is parallel to the molecular beam propagation direction. The source chamber is pumped with two
turbomolecular pumps (Pfeiffer Vacuum ATH2303 M), which results in a total pumping speed of
$\ordsim3600$ l/s for helium. The source chamber is equipped with two quick-access DN250 CF doors
(Pfeiffer Vacuum 420KTU250), enabling easy sample replacement in the Even-Lavie valve sample
container <cit.>. Furthermore, the source chamber is separated from the deflector
chamber with home-built automatic DN250 CF gate valve <cit.> to allow sample
changes without breaking the vacuum of the deflector chamber and the spectrometer. The deflector
chamber is also pumped with a turbomolecular pump with a pumping speed in the order of
$\ordsim1800$ l/s for helium (Pfeiffer Vacuum ATH2303 M). Next, the VMI chamber is pumped by two
turbomolecular pumps (Pfeiffer Vacuum PM P06 302 and PM P06 312) resulting in an overall pumping
speed of $\ordsim4000$ l/s for helium. Finally, the beam dump is pumped with a turbomolecular pump
(Pfeiffer Vacuum PM P04 024) resulting in a pumping speed of $\ordsim470$ l/s for helium. In
addition, a 250 mm diameter copper shield connected to a dewar for liquid nitrogen cooling is
located in the VMI chamber to increase the pumping speed for molecules with freezing points above
the liquid nitrogen temperature. The Dewar is designed such that liquid nitrogen must be refilled
every $\ordsim7$ hours to keep the shield cold.

The first pre-vacuum line that provides vacuum to the source chamber and to the second low pressure
pre-vacuum line, is pumped by an pre-vacuum pump (Pfeiffer Vacuum PR P00 009) with a nominal pumping
speed of $\ordsim 130~\mathrm{m^3/h}$. The second pre-vacuum line is pumped by a turbomolecular pump
(Pfeiffer Vacuum PM P03 900) with a pumping speed of 255 l/s for helium. The second pre-vacuum line
provides vacuum to the deflector-, VMI-, and beam dump chambers. This line is needed to compensate
for the low helium compression ratio given for, , the pumps of the VMI chamber by
$3\times10^{5}$. When operating the molecular beam valve at a repetition rate of 100 Hz, the typical
pressures in the source-, deflector-, VMI-, and beam dump chambers are $\ordsim 10^{-6}$, $\ordsim
10^{-8}$, $\ordsim 10^{-10}$, and $\ordsim 10^{-10}$ mbar, respectively. Under such conditions, the
typical pressure in the first and second pre-vacuum lines is $\ordsim 10^{-3}$ and $\ordsim
10^{-7}$ mbar, respectively.  Further information on the vacuum system including the layout of gas
lines, the interlock system with its logic to control the gate valves, and voltage power supplies,
is provided in the supplementary information. The primary purpose of the interlock system is to
ensure operation monitoring of the eCOMO vacuum system and to automatically protect the setup and
spectrometer in case of sudden emergency situations.

§.§ Molecular beam setup

The main part of <ref> shows the molecular beam setup. All molecular beam
components including the valve, skimmers, deflector, and the knife edge are mounted on motorized
translation stages in the $x$–$y$ plane. The individual positions can therefore be adjusted
perpendicular to the molecular beam propagation direction using linear feedthroughs and stepper
motors (Thorlabs, DRV014), which can be controlled remotely. The maximum travel range for the whole
molecular beam is $\ordsim11$ and $\ordsim12$ mm in $x$- and $y$ direction, respectively. The
implemented motors enable the alignment and vertical scanning of the molecular beam with respect to
any general photon source, enabling the investigation of different sections of the
quantum-state-dispersed molecular beam.

An Even-Lavie valve <cit.> is used to deliver the molecular beam by a supersonic
expansion of 10–100 bar helium with typically a few mbar of sample molecules into
vacuum <cit.>. The valve temperature can be adjusted between room temperature and
250 to control the vapor pressure of the sample under investigation and can be operated
with a repetition rate of up to 500 Hz. The valve is mounted on a semi-motorized DN63 CF 3D
manipulator (Caburn-MDC E-PSM-2504) with a travel range of 7 cm along the molecular beam propagation
direction. The first conical skimmer (Beam Dynamics, model 50.8, $\varnothing=3$ mm) is used for
skimming the molecular beam and to obtain differential pumping. The distance between the first
skimmer and the valve can be adjusted between 9 and 16 cm by moving the valve along the molecular
beam propagation direction. The first skimmer is mounted on a home-built flange with a built-in 2D
translation stage to adjust the skimmer position in the $x$- and $y$ directions. A second skimmer
(Beam Dynamics, model 40.5, $\varnothing=1.5$ mm) is placed 13.01 cm behind the first skimmer tip
for further collimation of the molecular beam. An electrostatic $b$-type
deflector <cit.>, which disperses the molecules in the molecular beam with
respect to their quantum states <cit.> is located 4.4 cm behind the
tip of the second skimmer. Both electrodes of the electrostatic deflector are connected through
high-voltage feedthroughs (Pfeiffer 420XST040-30-30-1) allowing voltages of up to $\pm30$ kV on each
electrode. The second skimmer and the electrostatic deflector are mounted together on an $x$–$y$
translation stage to adjust their positions. The electrostatically dispersed molecular beam is then
cut by a vertically adjustable knife edge placed 4.1 cm downstream of the deflector exit. This
allows for both, an improved sample separation and a higher column density of the molecular
beam <cit.>. The molecular beam is further skimmed by a third conical skimmer
(Beam Dynamics model 50.8, $\varnothing=1.5$ mm) placed 2.37 cm downstream of the knife edge,
providing differential pumping against the VMI chamber. This third skimmer is again mounted on a 2D
translation stage. The molecular beam enters the VMI chamber after the third skimmer. The distance
between the tip of the third skimmer and the interaction region is 27.8 cm. The state-selected
molecules are ionized by the photons used for the investigation inside the double-sided VMI
spectrometer <cit.>. The resulting electrons and ions are velocity-mapped onto a
120 mm (Photonis Model: APD 3 PS 120/32/25 19 I 60:1 NR 10”FM P47) and a 75 mm (Photonis Model: APD
3 PS 75/32/25/8 I 60:1 NR MGO 8”FM P47) diameter position-sensitive detectors, respectively. Each
detector consists of three z-stacked multi-channel plates (MCPs) and a fast phosphor screen
(P-47). Both detectors are captured by time and position-sensitive Timepix3 cameras built by
Amsterdam Scientific Instruments based on homebuilt
chips <cit.>. The molecular beam enters the beam dump
chamber through a $\varnothing=2$ cm aperture after passing the VMI. The beam dump chamber pump is
mounted with a 22 mm offset regarding the z-axis, allowing the molecular beam to directly hit onto
the rotors to prevent scattering of the direct beam back into the VMI chamber.

§.§ Double-sided VMI

The double-sided VMI setup specifically designed for eCOMO is illustrated in
<ref>a]fig:VMI-setup along with its relevant dimensions and the electrode labels.

(a) Schematic of the double-sided VMI setup cut along the $y$–$z$ plane containing
the symmetry axis. The left side represents the electron detector, while the right side
represents the ion detector. All distances are in mm. Typical voltage settings for each
electrode can be found in <ref>. The interaction region is
indicated by the red star. (b) Simulation for the large velocity range mode calculated
with SIMION. Equipotential lines are displayed in red. The trajectories of electrons are
depicted in green, while ion trajectories are represented in blue.

It allows for the simultaneous detection of electrons and ions and evolved from the initial concept
of velocity imaging in photoionization coincidence studies <cit.>. The VMI's
adaptable design permits both, vertical and horizontal installations within the VMI chamber,
facilitating flexible use in different laboratories and photon-science facilities. By default, the
eCOMO spectrometer is installed vertically in the VMI chamber. However, the spectrometer can be
installed alternatively horizontally along the molecular beam direction replacing the beam dump. In
this case, only the long side of the spectrometer can be used for either electrons or ions. The
VMI's point of origin is defined in the middle between the electrodes EE1 and IE1 as indicated by
the red star. This is the location where the molecular beam interacts with the photons. Electrodes
EE1 and IE1 are fabricated from oxygen-free copper (OF-Cu) with at least 5 of gold plating to
obtain a higher work function compared to stainless steel for experiments using high energetic
photons, , UV light. The remaining electrodes are made of stainless steel (1.4429 ESU), with
surfaces finely polished. Additionally, a rim structure is designed around the outer edge of each
electrode to reduce the influence of the VMI holder on the shape of the electric field in the center
of the spectrometer. All electrodes have an outer diameter of 140 mm and are separated by ceramic
spacers (Macor) as electrical insulators. A cylindrical shield surrounding the VMI, depicted in
green, with holes for the photons, molecular beam, potential diagnostic paddles for x-ray
experiments as well as for pumping, is composed of $\mathrm{\mu}$-metal to provide shielding from
stray fields and external magnetic fields. Double-layered meshes EM1–EM2 and IM1–IM2, with 80%
transmission for each layer, are positioned on both sides of the spectrometer to ensure homogeneous
drift tube fields close to the detector surfaces. The first mesh on both sides with respect to the
interaction region is conductively connected to the corresponding drift tube. The VMI is completely
mounted on the DN200 CF ion detector flange. Its design allows for straightforward installation or
replacement entirely from the ion detector side. The electrode connections extend through wires to
the DN200 CF flange, where they proceed to the exterior via electrical feedthroughs
(CF16-SHV20-SH-SE-CE-NI, Vacom) with a maximum rated voltage of 20 kV each.

The spectrometer can operate in diverse photo-electron and photo-ion detection modes by adjusting
the voltages across the electrodes as detailed in <ref>.

Imaging mode     SMI     SVR–VMI     LVR–VMI

Electron mesh 2 (EM2)    0     -134     6864

Electron drift tube (EDT)    0     1956     9816

Electron extractor 2 (EE2)    0     490     673

Electron extractor 1 (EE1)    0     -277     -553

Ion extractor 1 (IE1)    0     -370     -593

Ion extractor 2 (IE2)    -420     -821     -4672

Ion extractor 3 (IE3)    -1600     -2373     -4717

Ion drift tube (IDT)    -1600     -2380     -4800

Ion mesh 2 (IM2)    -1600     -821     -4672

The voltage settings (in Volt) required for the corresponding electrodes, as displayed
in <ref>(a)]fig:VMI-setup to enable the operation of the double-sided VMI in the
specified modes. See text for details.

This table lists, as examples, the implementation of a spatial map imaging (SMI)
mode <cit.> – optimized exclusively for ions – and two VMI modes for different
maximum energies for the velocity components parallel to the detection plane $E_\text{p}$. The VMI
modes images both, electron and ion velocities. Each voltage setting is tailored for specific
applications. In the SMI mode, the projected position on the detector correlates with the ion's
initial position and is in the first order unaffected by its initial velocity, enabling precise
imaging of the ionizing region. This functionality is particularly advantageous for the alignment of
the experimental setup, ensuring accurate positioning of the light focus and the molecular beam.
Conversely, the VMI mode allows electrons or ions with identical velocity vectors to converge to a
single spot on the corresponding detector, in the first order irrespective of their initial
positions. Different VMI modes allow different maximum kinetic energies to be collected on the
detector, which can be adjusted depending on the experimental situation. The detection range for the
electron kinetic energy parallel to the detector surface extends up to 700 eV in the
large-velocity-range VMI (LVR–VMI) mode. Single-charged ions can be mapped here with a parallel
energy up to 12 eV. The range for the parallel electron kinetic energy extends up to 100 eV and for
single-charged ions up to 10 eV in the small-velocity-range VMI mode (SVR–VMI).

The combination of MCP and phosphor screen amplifies each individual ion and electron after the
particles hit the detector, by creating an avalanche of electrons that then scintillates on the
phosphor screen. The light-flashes are captured with a Timepix3 camera <cit.> operated and controlled by our open source library
PymePix <cit.>. PymePix is also utilized to extract the raw
physics events from the Timepix3 data stream, converting them into centroided events recognized as
single particle events <cit.>. The data-driven cameras, with a temporal
resolution of $\ordsim 1.6$ ns, enable us to operate in multi-mass detection mode and to directly
obtain the VMI images of all fragments by slicing the TOF coordinate in software during
analysis. The performance evaluation of the Timepix3 camera is beyond the scope of this article and
was extensively discussed in other publications <cit.>. The relative positions and relative TOFs measured by the
cameras can be converted to velocities in VMI conditions.

The energy resolution of the double sided VMI was determined through SIMION simulations. It is in
the first approximation linearly dependent on the source size so the volume in which the ions and
electrons are produced. In addition, it depends on the specific mode of spectrometer operation. In
our specific simulation we used a cylinder for the focal volume with the length given by the
diameter of the molecular beam $D=2$ mm, and the radius given by the standard deviation of the
transverse laser intensity distribution $\sigma_{r}=21$ resembling our experimental conditions.
Taking LVR–VMI as an example, we obtain $\Delta E_\text{p}/E_\text{p} = 0.004$ for the ion side, which
is to first order constant for all energies that can be mapped. This results in an energy resolution
of $\Delta E_\text{p}=0.02$ eV for, , a parallel ion energy of $E_\text{p}=5$ eV. For electrons,
using the same settings, we obtain $\Delta E_\text{p}/E_\text{p} = 0.003$. Both resolutions are in
good agreement with results obtained on similar systems <cit.>. Further settings to
meet special experimental requirements are possible, but must be again found by simulations.

§.§ Portability and Installation

The eCOMO setup is designed for easy and convenient transport between different light sources, which
is facilitated by four integrated wheels. Once eCOMO is positioned, four adjustable home-built
support feet connected to a screw jack with standing screw (Nozag NSE25), capable of both,
horizontal and vertical adjustments, allow precise 2D positioning with a travel range of
$\pm3$ cm. The overall height of the interaction center can be adjusted between 103 and 135 cm. The
source, deflector, and VMI chambers are mounted on profile guide rails attached to the stainless
steel platform by guide carriages, allowing further adjustments along the molecular beam
direction. The eCOMO setup incorporates all necessary equipment, as depicted schematically in the
upper left corner of <ref> in the green highlighted area, except for the
computer used to control eCOMO and acquire data <cit.>. Integrated within
the setup are all required power supplies with ports allowing quick installation, including three
power cables capable of handling the necessary electrical load, a pair of inlet and outlet pipes for
water circulation to cool the molecular pumps, an exhaust pipe for the vacuum system, and multiple
gas sources such as nitrogen, compressed air, and the molecular carrier gas. The total weight of
eCOMO including the pumps and electronics is $\ordsim1800$ kg. The overall dimensions of eCOMO are
90, 175, and 240 cm in width, height, and length, respectively.

§ CHARACTERIZATION OF THE SETUP

§.§ Molecular beam spatial profiles

We used OCS as a sample to characterize the molecular beam and VMI performance to commission
eCOMO <cit.>. OCS was ionized by pulses from a commercial
Ti:sapphire femtosecond laser system (Coherent Astrella) operated at a repetition rate of 1 kHz. The
laser beam was compressed to approximately 40 fs with a standard grating-based compression setup.
The laser beam was directed through a 75 cm focal-length lens, which resulted in a peak intensity of
$I_0\approx3\times10^{14}$ . The laser focus was spatially and temporally overlapped with the
molecular beam in the interaction center of the VMI spectrometer. The pulsed molecular beam was
produced by expanding a mixture of 100 bar helium and $\ordsim10$ mbar OCS. The gas mixture was
expanded into the source chamber through the Even–Lavie valve operated at a temperature of
50 and at a repetition rate of 250 Hz. Laser shots without the molecular beam were used
for background measurements. Optimized operation conditions of the valve allowed for the generation
of a dense ($>10^{7}~\invccm$) molecular beam in the interaction volume. The OCS sample was
displaced toward negative $y$ positions when a total voltage of 25 kV was applied across the
deflector, , de facto the populated quantum states of OCS exhibited a strong-field-seeking
character at the relevant electric field strengths in the deflector. The measured deflection profile
of OCS is shown in <ref> as red dots.

Measured molecular beam density profiles (points) and simulation (solid lines) for the
deflector switched off (blue) and at a total voltage of 25 kV (red) applied across the
electrodes. The inset shows the number of OCS^+ ions $\Sigma$ per shot (black squares)
as a function of the laser peak-intensity natural logarithm $\ln I_{0}$. The purple dashed
line represents the fit of the asymptotic slope. The data in the inset was measured at
$y=-0.5$ mm in the deflected molecular beam profile, marked by the black arrow.

Solid lines represent simulated OCS beam profiles, computed through Monte-Carlo sampling combined
with trajectory calculations that incorporate the geometrical constraints of the mechanical
apertures in the experimental setup. These simulated deflection profiles of OCS match the
experimental data, assuming an initial rotational temperature of $T_\text{rot}=0.5$ K for the
molecular beam entering the deflector. The results affirm that the eCOMO setup can generate cold and
stable molecular beams.

The OCS beam sample density was determined using a strong-field ionization
model <cit.>. The asymptotic slope of an integral ionization
signal $\Sigma(I_{0})$, with respect to natural logarithm of the peak intensity $\ln I_{0}$
approaches a limiting value as $I_{0}$ approaches infinity. This slope can be expressed as:

\begin{equation}
\lim_{ I_{0} \to \infty} \left(\frac{d\Sigma }{d\ln I_{0} } \right) = 2\pi \alpha \sigma _{r}^{2} D\rho
\end{equation}

where $\rho$ represents the sample density, $\sigma _{r}$ the standard deviation of the transverse
laser intensity distribution, $\alpha$ the instrument sensitivity, $D$ the length of the focal
volume in the molecular beam, and $\Sigma$ the total number of detected ions. The inset of
<ref> illustrates the parent-ion signal count $\Sigma$ per shot as a
function of the laser peak intensity $I_0$ on a semi-logarithmic scale. The data was measured at
$y=-0.5$ mm in the deflection profile as indicated by the black arrow. The asymptotic slope and
saturation onset were obtained by fitting a straight line through the highest peak-intensity points,
resulting in a saturation onset of $I_\text{sat} = (9.4 \pm 0.3_\text{stat} \pm 1.9_\text{syst})
\times 10^{13}~$. The subscripts indicate statistical and estimated systematic errors. This
is in good agreement with the previously reported saturation intensity of OCS measured under similar
experimental conditions given by $I_\text{sat} = 7.2 \times
10^{13}~$ <cit.>. To extract the sample density we determined the combined
sensitivity of the meshes and the MCP detector to be $\alpha = 35.9\% \pm 1.3\%$ as described in
<ref>. In addition we use $\sigma _{r}=21$ and $D=2$ mm for the laser intensity
standard deviation and the molecular beam width, respectively. The resulting sample density at
$y=-0.5$ mm is given by $\rho = (2.86 \pm 0.58_\text{stat} \pm 0.57_\text{syst}) \times
10^{6}~\invccm$. The maximum density is $\rho \approx 10^{7}~\invccm$ in the direct beam as shown in
<ref>, which is one order of magnitude lower than previously reported
on a similar experimental setup <cit.>, which is attributed to a lower
concentration of OCS in the helium mixture.

§.§ Time-of-flight mass spectrum (TOF-MS)

<ref> shows the time-of-flight mass spectrum (TOF-MS) resulting from strong-field
ionization of OCS and from the background gas in the VMI chamber measured with the Timepix3
camera <cit.>.

TOF-MS spectrum of ions originating from the molecular beam containing OCS (blue) and
from VMI-chamber background gases (red). The inset provides a magnified view of the specific
region where OCS fragment ions and more highly charged ions occur.

The figure displays two distinct regions: The first region depicts the OCS parent ion and its
isotopologues O$^{13}$CS and OC$^{34}$S, which correspond to 60, 61, and 62,
respectively. The respective signal strength is in good agreement with the corresponding natural
abundance of OCS's most abundant atmospheric isotopologues <cit.>. The
isotopologues with masses 63, and 64 are also present in the data, but are not visible in
the graph shown because of their natural abundance below 0.1%. The second feature is the presence
of fragments, as also seen in the zoomed-in subfigure of the graph. Two broad peaks correspond to
CO^+ with 28 and S^+ with 32 ions, which exhibit a broad ion velocity
distribution due to Coulomb explosion of OCS^2+. The broad velocity distribution is directly
reflected in the broad TOF distribution. The relatively strong background here originates from
impurities of O_2 and N_2 gases present in the VMI chamber and the molecular
beam. Additionally, two narrow peaks are observed in the inset, which correspond to doubly ionized
OCS and singly ionized H_2S, which is an impurity also present in the molecular beam. The
significant water ion peak can also be mostly attributed to poor background pressure conditions,
which at the time, were not fully optimized and given by $P_\text{VMI} = 4 \times 10^{-9}$ mbar when
the molecular beam was switched off. Additionally, there are small contributions of water ions that
arose from impurities in the molecular beam. Recent modifications of the setup showed considerable
improvement, achieving pressures as low as $P_\text{VMI} = 3 \times 10^{-11}$ mbar. Furthermore, the
purity of the gas line was also improved by installing a cooling thermostat (ECO RE 1050 S, LAUDA).

§.§ Performance of the double-sided VMI

In order to experimentally characterize the double-sided VMI, both electrons and ions from the
molecular beam containing OCS were examined after strong-field ionization. <ref>a]fig:VMI and
<ref>c]fig:VMI show the background subtracted cation-fragment velocity maps for 28 and
32, respectively. The laser-peak-intensity here was $I_0 = 2.5 \times 10^{14}$ .

(a) and (c) Projected ion velocity maps for 28 and for 32 fragments,
respectively. The laser polarization is indicated by the red arrow. (b) and (d) Projected
electron velocity maps for electrons correlated with OCS^+ and H_2O^+,
respectively. (e) and (f) Projected radial velocity distributions for ions and electrons,
respectively. (g) Ion-electron Timepix3-camera-trigger correlation. See text for more details.

The velocity maps show the commonly known features of OCS fragments after strong field
ionization <cit.>. The central circular parts of both
velocity maps are attributed to fragmentation to CO^+ and S^+ following single ionization.
The intense central spots are attributed to N$_2^+$ and O$_2^+$ in the 28 and 32 velocity
maps, respectively. These were impurities present in the molecular beam. The outer rings are
ascribed to the Coulomb explosion channels OCS + $nh\nu$ →OCS^+2 →CO^+ +
S^+. The reduced counts right above the central regions are attributed to strong signals from
background ions such as H_2O^+ arriving shortly before the selected ions in combination with
the dead time of the Timepix3 camera <cit.>. In addition, the detector could
be less sensitive there because all background ions hit this area. The projected radial velocity
distributions are shown in <ref>e]fig:VMI. The multiple channels corresponding to the VMIs are
visible again.

The ion information can be used to determine the electron velocity map corresponding to specific
ionization or fragmentation channels. To do so, both Timepix3 cameras have to be first synchronized
as they independently record electrons and ions. To correlate the ion and electron cameras we
simultaneously took ion and electron data at low laser intensities ($I_0 = 3 \times
10^{13}$ ) at an average event rate of 0.02 ions or electrons per shot. The number of ions
$N_\text{i}(N_\text{trg;i})$ or electrons $N_\text{e}(N_\text{trg;e})$ were assigned to the
camera-specific trigger numbers $N_\text{trg;i}$ and $N_\text{trg;e}$, respectively. The normalized
convolution $N_\text{i} \ast N_\text{e}$ of this signal is shown in <ref>g]fig:VMI as a
function of the relative trigger number $\Delta N_\text{trg} = N_\text{trg;i}-N_\text{trg;e}$. The
peak indicates the required relative trigger number shift for both cameras to be synchronized.

To demonstrate ion-electron correlations, we ionized OCS with a laser-peak-intensity of
$I_0 = 8 \times 10^{13}$ to have both, parent ions and fragments present. In this case, the
mean number of detected ions per shot was $\ordsim2$ so specific ions and electrons from single
laser shots can not be correlated directly. We can, however, obtain the statistically averaged
electron image using the covariance of the electron velocity maps with the boolean operation on a
specific ion detected <cit.>. <ref>b]fig:VMI displays the photoelectron
velocity distribution after the covariance analysis where an OCS^+ was measured in
coincidence. The figure reveals intricate angular and radial structures as previously
reported <cit.>. Specifically, the momentum distributions exhibit wide radial
patterns that repeat at integer multiples of the photon energy, which is commonly referred to as
above-threshold ionization (ATI). A detailed discussion about the structure of the OCS electron
momentum map is beyond the scope of this manuscript and further information can be found
elsewhere <cit.>. <ref>d]fig:VMI displays the photoelectron momentum
distribution accordingly when an H_2O^+ from ionization of the background gas was measured in
coincidence with an electron. A clear difference in the image structure is observed in comparison
with <ref>b]fig:VMI, which is attributed to the different ionization channels. The differences
become even more obvious in the electron's speed distribution projected onto the detector shown in
<ref>f]fig:VMI. In comparison with previous work <cit.>, our short data
acquisition time led to smaller experimental statistics, resulting in less detailed ATI structures.
Further discussion about the underlying reason with respect to the differences in the electron
velocity distributions is again beyond the scope of the manuscript. The clear differences between
the two electron velocity maps however clearly show that eCOMO allows to correlate electron and ion

The low-intensity correlation measurements used for <ref>g]fig:VMI enabled us furthermore to
obtain the detection efficiency for electrons and ions in our setup (see SI for details). From our
data we obtain a final total detection efficiency of $\alpha_\text{i}=35.9\% \pm 1.3\%$ and
$\alpha_\text{e}=35.4\% \pm 1.3\%$. Assuming a detection efficiency of the MCP given by
$\alpha_\text{MCP}=50$% we obtain a transparency for a single mesh $\alpha_\text{mesh}=83.7\% \pm
1.5\%$, with the estimated transmission of 80%. The deviation is attributed to the fact that the
effective transmission of the second grid is higher due to the shadow cast by the first grid.

§ CONCLUSIONS

In summary, eCOMO was developed, set up, and commissioned using experiments with an ultrashort-pulse
table-top laser system. The performance of the cold ($T_\text{rot}\le1$ K), pulsed, and dense
($\rho \approx 1 \times 10^{7}~\invccm$) molecular beam was characterized. SIMION simulations of the
double-sided VMI indicated that under experimental conditions, an energy resolution of
$\Delta E/E < 0.4 \%$ is achievable. We showed that the combination of the double-sided VMI and two
synchronized Timepix3 cameras at eCOMO works and that the individual ions can be separated in the
data analysis based on their mass over charge ratio. Furthermore, specific molecular dynamics can be
separated for systems with many fragmentation channels as shown here for OCS after ionization. Using
a covariance analysis, electron VMIs and corresponding projected electron velocity distributions
were assigned to specific mass over charge ratios. The absolute detection efficiency of the setup
was measured by coincidence methods to $\alpha_\text{i}=35.9\% \pm 1.3\%$ and
$\alpha_\text{e}=35.4\% \pm 1.3\%$ for ions and electrons, respectively. The use of the $b$-type
electrostatic deflector allows to separate specific rotational quantum states, conformers and
molecular cluster species with different effective-dipole-moment-to-mass
ratios <cit.>. This feature of the eCOMO setup has great
potential in various areas, such as the separation and purification of complex clusters for chemical
dynamics studies <cit.>.

The eCOMO setup exhibits versatile applications in the production of purified cold gas-phase
molecular beams and high-precision ion-electron coincidence and covariance measurements. It can be
readily adapted to a wide range of photon-source facilities, making it suitable for diverse research
activities spanning atomic, molecular, and cluster physics, as well as materials science, energy
science, chemistry, and biology.

§ ACKNOWLEDGMENTS

We thank Ilya Maksimov for his experimental contributions in an early phase of the project and
Daniel Gusa for technical support.

We acknowledge financial support by Deutsches Elektronen-Synchrotron DESY, a member of the Helmholtz
Association (HGF), also for the provision of experimental facilities and for the use of the Maxwell
computational resources operated at DESY. This research was supported by the Center for Molecular
Water Science (CMWS) in an Early Science Project. The research was further supported by the Cluster
of Excellence “Advanced Imaging of Matter” (AIM, EXC 2056, ID 390715994) of the Deutsche
Forschungsgemeinschaft (DFG), the Helmholtz Foundation through funds from the Helmholtz-Lund
International Graduate School (HELIOS, HIRS-0018), the German Federal Ministry of Education and
Research (BMBF) and the Swedish Research Council (VR 2021-05992) through the Röntgen-Ångström
cluster project “Ultrafast dynamics in intermolecular energy transfer: elementary processes in
aerosols and liquid chemistry” (UDIET, 05K22GUA), and the European Union's Horizon 2020 research
and innovation program under the Marie Skłodowska-Curie Grant Agreement “Molecular Electron
Dynamics investigated by Intense Fields and Attosecond Pulses” (MEDEA, 641789). We acknowledge
CAMP@FLASH at DESY for providing some material for the spectrometer.

§ DATA AVAILABILITY

The data recorded in the laser lab at DESY are available from the corresponding author upon
reasonable request.

§ CODE AVAILABILITY

The code for the data analysis is available from the corresponding author upon reasonable request.