# Confronting the thermodynamics knowledge gap: A short course on
computational thermodynamics in Julia

Luc T. Paoli  Department of Chemical Engineering, Imperial College London,
London SW7 2AZ, United Kingdom  Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, MA 02139, United States  Pavan K. Inguva
Department of Chemical Engineering, Imperial College London, London SW7 2AZ,
United Kingdom  Massachusetts Institute of Technology, 77 Massachusetts
Avenue, Cambridge, MA 02139, United States  Andrew J. Haslam*, Department of
Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom
Pierre J. Walker*, Department of Chemical Engineering, Imperial College
London, London SW7 2AZ, United Kingdom  Division of Chemistry and Chemical
Engineering, California Institute of Technology, Pasedena, CA 91125, United
States

###### Abstract

Computational elements in thermodynamics have become increasingly important in
contemporary chemical-engineering research and practice. However, traditional
thermodynamics instruction provides little exposure to computational
thermodynamics, leaving students ill-equipped to engage with the state-of-the-
art deployed in industry and academia. The recent rise of easy-to-use open-
source thermodynamic codes presents an opportunity for educators to help
bridge this gap. In this work, we present a short course that was developed
and rolled-out using the Clapeyron.jl package, the material of which is all
openly available on GitHub. The course can serve as a foundation for others to
similarly integrate computational material in thermodynamics education. The
course is structured into three sections. Section one serves as a refresher
and covers core material in numerical methods and thermodynamics. Section two
introduces a range of thermodynamic models such as activity-coefficient models
and cubic equations of state, outlining their implementation. In section three
the focus is moved to deployment, guiding students on how to implement
computational-thermodynamics methods covering volume solvers, saturation
solvers, chemical-stability analysis and flash problems. In a pilot study
conducted with both undergraduate and graduate students, participants found
the material engaging and highly relevant to their chemical-engineering
education.

11footnotetext: Corresponding authors. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

The application of thermodynamics in industry and research is becoming
increasingly diverse and sophisticated with advances in computational
capabilities, available resources, theoretical understanding, and use-cases
[1, 2, 3, 4]. The expanding role that thermodynamics plays in the progress of
science and engineering is not surprising considering its centrality in
describing physical systems. Nevertheless, it places a challenge on educators
to provide suitably advanced thermodynamics instruction to ensure students are
well-equipped to meet their professional needs after graduation. Accordingly,
the continued development of education opportunities and resources for
thermodynamics remains as important as ever.

One of the most common tasks encountered is the use of various thermodynamic
models and methods to analyze systems of interest, to generate estimates for
various properties (e.g., phase diagrams or activity coefficients), or as part
of a process / physics-based simulation at a larger scale. In that regard,
there are several currently available software packages of varying degrees of
functionality that have been used in professional and/or teaching settings
e.g., ASPEN [5], PYroMAT [6, 7], CoolProp [8], XSEOS [9], teqp [10] , and
custom code [11]. A key component that is implicit in the use of these models
and software packages is the computational aspect of thermodynamics, yet this
is often neglected in thermodynamics courses, even at a graduate level. Often,
one of the major selling points of many software packages is the abstraction
of the numerical implementation of the model and algorithms, which can be
complex, thus enabling users to rapidly employ advanced thermodynamic models
in their workflow (e.g., see [12, 10]).

In conceptualizing an advanced teaching course, it is imperative to not treat
the various models and packages as “black boxes” which “magically” generate
the desired answer [13, 14], as the objective of such advanced courses is to
help students develop deeper insights and become informed practitioners.
Therefore, course content and packages employed should not only allow students
to engage meaningfully with the theory and use of various thermodynamic models
(including state-of-the-art models where practicable), they should also
provide the opportunity to explore various back-end details such as the
numerics and algorithms needed to compute the desired result. Other
considerations that should be taken into account in relation to software for
use in educational settings include availability and ease-of-use [14, 15, 16].

In this work, we outline the development and roll-out of a short course for
introducing computational thermodynamics to advanced undergraduate / graduate
students. The course is intended to be a stand-alone course that can be
delivered on an ad-hoc basis, though we intend to integrate components of the
course into a graduate thermodynamics module. The course was designed to be
modular so that self-directed learners and/or educators can adapt specific
components for their own needs or develop an expanded course using the
material presented as a foundation. All course material and code are currently
publicly available in a GitHub repository and additional contributions are
planned to expand the material offered.

The remainder of this article is laid out as follows: the context in which the
course is / can be taught is provided in section 2. After a short discussion
on the computational environment (section 3), the structure of the course is
provided in section 4. The results and feedback obtained from a pilot run of
the course is presented and discussed in section 5. Finally, the overall
findings of this article are summarised in section 6.

##  2 Course Context

The short course was developed to serve a dual role of introducing advanced
undergraduate students and graduate students to advanced thermodynamic
modelling tools and real-world applications of the statistical thermodynamic
theories. Students taking this course should already be familiar with the
fundamentals of thermodynamics (for example, the laws of thermodynamics,
Maxwell relations, and phase equilibria) and, perhaps, some thermodynamic
models (such as the ideal gas, van der Waals equation of state, and ideal-
solution models), as well as having taken some form of introductory coding
class. The latter is not necessarily a strong requirement as the course
material is delivered in such a way that the students only do minimal coding.
Within the context of the Imperial College undergraduate degree program, by
the end of the second year, students will have taken two courses on
thermodynamics, covering all of the topics mentioned previously, and one
introductory course on MATLAB coding, which we deem more than sufficient for a
student intending to take part in this course.

For the benefit of those students who are not familiar with all of the above
pre-requisites or require a refresher for any of the material, the first
section of the course provides sufficient introductory content. Conversely, to
suit participants who are already familiar with large sections of this course,
or instructors who wish to use only specific sections in their own class, the
course is designed to be modular, allowing the user to engage with the
material flexibly.

##  3 Computational Environment

The main software package used in this course is the recently developed open-
source fluid thermodynamics toolkit Claperyon.jl [12] which is implemented in
the Julia language. The release of Claperyon.jl has enabled users to employ a
large variety of thermodynamic models (ranging from activity-coefficient
models and standard cubics to the state-of-the-art SAFT equations) in an easy-
to-use and extensible manner. The package has several built-in capabilities
for thermophysical-property estimation and is currently under active further
development to incorporate additional capabilities to make it the go-to
thermodynamics toolkit. Claperyon.jl has garnered significant interest from a
diverse group of users with over 135 stars on GitHub.

Although the Julia language is less mature than alternatives such as Python
and MATLAB, it has several features and packages that make it appealing for
scientific computing while also being well-suited in an educational setting
[17, 18]. In terms of ease-of-use, Julia has a similar syntax to MATLAB and
Python which helps reduce challenges students may have with reading and
writing code [19]. Julia also has native package-management tools that can
negate the need for setting up environments, reducing time spent introducing
students to the coding environment.

Figure 1: Example taken from the Pluto notebook used in the course.
Explanatory text is here highlighted in yellow; example code is highlighted in
green. To reveal the hidden code used to generate the embedded chart, students
would simply select the symbol highlighted in red. The chart represents the
isobaric heat capacity of carbon dioxide obtained using different equations of
state.[20, 21]. Students are invited to vary which species and equation of
state is used to examine how the ideal contribution plays a role in
determining real properties.

One particular package we would like to highlight is Pluto.jl [22], which
provides reactive notebooks. For the presentation of a course developed
previously on introducing partial-differential-equation solvers[15], a
combination of ready-made Python scripts for code and accompanying pdf files
was used for the course material. From student feedback, two main challenges
were noted: some students experienced difficulties in setting up the Python
environment needed to run the scripts; and the lack of active and structured
learning activities hampered student engagement with the material. Pluto.jl
can elegantly resolve both these issues. For this course, we decided to
integrate both the lecture notes and codes for a given section within a single
file. Illustrated visually in figure 1, within Pluto notebooks, we are able to
provide both the code needed to generate the results as well as the text
explaining these results. As seen in figure 1, the explanatory text along with
the main function call and plotting code are embedded seamlessly. Pluto has
the advantage of avoiding clutter in the notebook by allowing some sections of
the code (e.g., plotting) to be hidden; hidden text can be revealed if desired
by selecting the designated symbol (highlighted in red in figure 1).

Although similar functionalities are available in more commonly used
interactive environments such as Jupyter notebooks, which have already found
widespread use in education[23, 24, 25], the reactive nature of Pluto has one
major additional benefit to users and educators. In the example shown in
figure 1, students are invited to modify the provided code and re-run the
cell. In doing so, the entire notebook ‘reacts’ to update every other cell
block, including the hidden plotting code. As a result, the figure will
automatically update based on the changes they have made. Interactive
notebooks by virtue of their sequential nature do not allow for such dynamic
behaviour [26]. We believe reactive notebooks allow students to engage more
intently with the lecture material and, thereby, foster independent
learning[27]. Further details on setting up the computational environment for
this course are provided in the online repository.

##  4 Course Structure

Progression was used as the guiding principle when designing the course, to
help accelerate students’ learning in a short period of time [28]. The course
is divided into three sections, each of which is delivered in a three-hour
lecture (including two ten-minute breaks). An outline of the course structure
is presented in figure 2. The first section serves as an introduction and
provides a refresher on the fundamentals of thermodynamics; important
numerical methods (such as Newton’s algorithm and automatic differentiation)
are outlined. In the second section five classes of equations of state are
examined, providing students with a high-level overview of each class, as well
as their respective advantages and disadvantages. The final section builds on
material from the previous two sections and examines various computational
thermodynamic problems, such as rooting-finding algorithms, stability analysis
and flash algorithms. Details on theory and numerical methods for each problem
are also covered.

Figure 2: Course structure summary with bullets providing a brief description
of each subsection. Dotted purple arrows denote concepts carried over between
subsections of the course.

###  4.1 Section 1: Fundamentals

The first section of the course is intended to serve as a review to ensure
participants are familiar with the pre-requisites of the course. Two main
topics are covered: thermodynamics; and numerical methods. In the
thermodynamics review, we revisit fundamental concepts such as thermodynamic
state variables (temperature, pressure, volume, number of moles etc.) and
state functions. We emphasise the point that from a state function expressed
in terms of three state variables, one can obtain a complete thermodynamic
description of the system using derivatives (e.g., Maxwell relations).

Having set the context for how thermodynamic properties can be generated, the
core numerical methods that will be utilised are then briefly covered. We
first introduce automatic differentiation, which can be used to obtain
numerically exact derivatives of arbitrary functions. Its use in many
computational applications, including machine learning, is highlighted to
students. Subsequently, root-finding algorithms (mainly Newton–Raphson and
fixed-point iterations) and optimisation are discussed. At this stage,
students are expected only to be aware of the existence of these methods, how
they are implemented (conceptually), and their advantages and limitations.
Students should recognise that errors arising during the implementation of
thermodynamic models and methods could originate from the numerical methods
used and not just from a limitation in their thermodynamics knowledge.

With the exception of a short exercise on how to use automatic differentiation
in Julia, this section of the course was delivered through PowerPoint slides
only. The remainder of the allocated lecture time was spent ensuring all
students had working notebooks, to allow them to actively take part in the
remaining two sections.

###  4.2 Section 2: Implementing thermodynamic models

In this section fives classes of thermodynamic models are discussed; these are
intended to cover the range of models that are typically encountered in
academia or industry. The order in which they are taught was chosen so as to
start with more-familiar approaches (ideal gas and ideal solution) and
gradually approach more-advanced models such as SAFT and empirical equations
of state.

####  4.2.1 Ideal-gas equation

We felt that students would be most comfortable starting with something with
which they should all be familiar; the obvious choice is therefore the ideal-
gas equation,

$pV=nRT,$ (1)

where $p$ is pressure, $V$ is volume, $n$ is the number of moles, $R$ is the
universal gas constant and $T$ is temperature, we quickly highlight that the
equation alone does not constitute the full picture of an ideal gas. By
integrating (1) with respect to volume in order to obtain the Helmholtz free
energy, $A$,

$A=-\int pdV=-nRT\ln{V}+c(n,T),$ (2)

we notice that there are missing contributions corresponding to the
translational, rotational and vibrational modes of motion. These contributions
are vital when obtaining real properties such as heat capacities and speeds of
sound, particularly in the gas phase [29]. We illustrate that these missing
contributions can be obtained using ideal-heat-capacity correlations, with
which students may already be familiar from mass- and energy-balances. The
importance of these contributions is highlighted in the case of real gases
(shown in the embedded chart in figure 1). Although there are no exercises in
this section, students are invited to vary the species and conditions for
which this chart is generated in order to appreciate the importance of these
contributions.

Enabled by the simplicity of the ideal-gas model, we take the opportunity to
introduce students to other important concepts in computational
thermodynamics. For instance, while correlations provide simple models to
address typically complex properties, one does need to be careful to only use
these for interpolation as the accuracy of these models can degrade quickly
when extrapolating. We also introduce students to the concept of group-
contribution approaches, which can be very useful when trying to model species
for which no (molecular) parameters exist. While such approaches are more
common in relation to other models discussed later in the course, it is
helpful to first introduce these approaches in the context of a more-familiar
model.

####  4.2.2 Activity coefficient models

Although some students may already be familiar with the ideal-solution model,
we begin this section by describing the assumptions implicit in this approach,
and how they lead to the more-familiar Raoult’s law,

$p_{i}=x_{i}p_{\mathrm{sat},i}\,$ (3)

where $p_{i}$, $p_{\mathrm{sat},i}$ and $x_{i}$ are the partial pressure,
saturation pressure and molar composition of species $i$, respectively. It is
then easy to introduce the concept of activity coefficients ($\gamma_{i}$) as
deviations from the ideal-solution model, motivating the interest in activity-
coefficient models. We also highlight some of the limitations of activity-
coefficient models, primarily being that they can only be used to obtain
mixture properties.

Having now established the ideal-solution model, it is then simple to explain
the local-composition model, from which all the activity models we present in
this section are derived. We use the Wilson model[30] where the key
interaction parameter, $\Lambda_{ij}$, between two species $i$ and $j$ is
given by:

$\Lambda_{ij}=\frac{v_{j}}{v_{i}}\exp{-\Delta\lambda_{ij}/(RT)}\,$ (4)

where $v_{i}$ is the molar volume of species $i$ and $\Delta\lambda_{ij}$ is
the interaction energy between $i$ and $j$. We then highlight the two primary
effects captured by local-composition models: size effects ($v_{j}/v_{i}$) and
enthalpic effects ($\Delta\lambda_{ij}$). As a result, it is straightforward
to explain the assumptions that underly the other models (NRTL[31] and
UNIQUAC[32]) discussed in this section, allowing us to focus on the advantages
and limitations of each. One of the key messages we highlight is that, to
decide which of these approaches is ‘best’ for a given application, one needs
to compare the predictions to experimental data.

Finally, as we have already introduced the concept of group-contribution
approaches in relation to ideal-gas models, it is straightforward to now
introduce the student to perhaps the most-commonly used activity coefficient
model: UNIFAC[33]. We also highlight another issue students may encounter when
using thermodynamic models, particularly UNIFAC: differences in models used by
specific individuals but often referred to by the same name, exemplified here
by the two versions of UNIFAC (the ‘original’ and Dortmund versions), which
are both referred to by the same name in literature. We emphasise the
importance of verifying the source of parameters, as well as their predictive
accuracy.

####  4.2.3 Cubic equations of state

Figure 3: Example of exercise given in the cubic equation of state section.
Each student is tasked with completing the code at the commented sections. A
textbox below the code block will dynamically update based on the student’s
progress.

As the most-common class of thermodynamic models used in industry, we dedicate
a significant portion of section two to cubic equations of state. To give some
physical insight behind cubic equations of state, we begin with a
phenomenological discussion of the van der Waals equation[34], and remind
students how the parameters can be obtained from the critical pressure and
temperature of any species, allowing for a wide range of applicability. From
here, we highlight the limitations of the van der Waals equation and how the
subsequent cubic equations of state represent improvements. For conciseness,
we introduce only the Redlich–Kwong (RK)[35], Soave-Redlich–Kwong (SRK)[36],
and Peng–Robinson (PR)[37] equations but do point out that many others have
been developed. In the case of the latter two equations of state, we also
introduce the concept of an $\alpha$-function and its relationship to the
acentric factor ($\omega$). We then emphasise one of the key advantages of
cubic equations of state: they have a simple and universal form,

$p=\frac{nRT}{V-nb}-\frac{n^{2}a\alpha(T)}{(V+r_{1}nb)(V+r_{2}nb)}\,,$ (5)

where $\alpha(T)$ is the $\alpha$-function and, $a$ and $b$ are substance-
specific parameters. The coefficients $r_{1}$ and $r_{2}$ are equation of
state-specific. We then give students the opportunity to implement their own
generalised cubic equation of state. We provide them with a code template with
parts left incomplete (shown in figure 3). They are then asked to define the
coefficients $r_{1}$ and $r_{2}$ for the van der Waals, SRK and PR equations
of state and, write out the residual Helmholtz free energy corresponding to
equation 5, accounting for each equation of state. Aside from just gaining
experience in coding, there are implicit learning objectives to this exercise.
In implementing a generalised cubic equation of state, the students make use
of the multiple-dispatch feature in Julia and, in order to ensure that the
automatic derivatives of the Helmholtz free energy are correct, they will have
to explicitly define all the state variables (an easy mistake to make whenever
implementing any equation of state). A text box is added below the code block
which provides real-time feedback based on their current implementation and
provides hints based on what part of their code is incorrect. This feature
arises from the reactive nature of Pluto notebooks.

Subsequently, we illustrate one of the key benefits of cubic equations of
state: modularity. We highlight how one can improve the modelling capabilities
of cubic equations of state by replacing the $\alpha$-functions and adding a
volume-translation correction. The impact on real properties is also
illustrated with multiple examples.

Figure 4: Pressure–composition ($p,xy$) diagram of benzene+methanol at 433.15
K obtained using PSRK and a custom cubic equation of state.[38, 39] Each
student is invited to attempt to surpass the accuracy of PSRK by developing
their own cubic equation of state, varying the $\alpha$-function, mixing rule
and activity coefficient model.

Having covered how one can use cubic equations of state to model pure species,
we now introduce students to the concept of mixing rules to model mixtures of
species. We illustrate how these mixing rules map the pure-component
parameters to mixture parameters in order to represent the mixture. We cover
two classes of mixing rules. The first is the standard van der Waals one-fluid
mixing rule, which is perhaps the most common one used in industry. It is
emphasised that this mixing rule is generally inaccurate without introducing a
binary interaction parameter and that these typically need to be fitted by
adjustment using experimental data. The second class of mixing rule we
introduce is that of the EoS/$G^{E}$ mixing rules. The previously-established
effectiveness of activity models in modelling mixture phase equilibrium is
used to motivate the concept of mixing rules which ‘pair’ an activity-
coefficient model with a cubic equation of state (such as the Huron–Vidal[40]
and Wong–Sandler[41] mixing rules).

Finally, uniting all of the previous concepts, we introduce the final and
most-modern class of cubic equations of state, the predictive cubics[38, 42].
These equations of state combine everything the students will have encountered
up until now: $\alpha$-functions; volume-translation corrections; and mixing
rules. We highlight that these represent the gold-standard for cubic equation
of state modelling in industry. However, to emphasise the importance of
validating these models against experimental data, we challenge students to
try to develop their own cubic equation of state, using all the tools they
have learnt up until now, to improve upon an existing predictive cubic
equation of state. As shown in figure 4, obtaining a better model is not
necessarily challenging with the tools students have at their disposal. Note
that this figure and, indeed, all the subsequent charts presented in this and
the following section, is shown here exactly as seen by the students in the
Pluto notebook.

Figure 5: Vapour–liquid envelope of methanol obtained using different
equations of state[37, 43, 21, 44, 45], demonstrating the significant improve
of SAFT equations over cubics, but the minimal difference between individual
SAFT-type equations.

####  4.2.4 SAFT equations of state

The Statistical Associating Fluid Theory (SAFT)[46, 47] is not typically
covered at any level within chemical-engineering curricula, primarily due to
do the complexity of its derivation. However, in the past few years, SAFT
equations of state have been shown to describe experimental data more
accurately for a larger range of properties and species than cubic equations
of state. As part of this course, we introduce students to SAFT equations of
state by providing a very high-level description of how these equations
represent species. Starting from the van der Waals equation, we highlight the
improvements made with each term added within the SAFT equation. Before going
into the individual SAFT variants, we summarise the overall SAFT approach.
Some details are provided on the implementation of SAFT equations, primarily
the association term which requires special treatment in order to be used. We
believe that a simple, high-level understanding is sufficient for students to
confidently use SAFT equations in their own work; should they want to find out
more, we provide references and invite them to take an in-depth statistical-
thermodynamics course provided within the college.

The remainder of this section is focussed on introducing the more-common SAFT
equations of state, explaining their differences at a high level and giving
their respective advantages and disadvantages. We emphasise that, as with the
activity-coefficient models, to decide which approach is best, one needs to
compare predictions to experimental data. However, in the case of SAFT
equations, as shown in figure 5, a lot of the model predictions are of
comparable accuracy. As a final point, we add that there is another factor to
take into account when considering whether to use a SAFT equation: with their
increased complexity, the computational cost of these approaches becomes
significant, meaning that the students would need to weigh the added accuracy
of these models with this additional cost.

####  4.2.5 Empirical equations of state

The last class of equations of state we cover in this course is empirical
equations of state. Like the ideal-heat-capacity correlations introduced in
section 4.2.1, these approaches are purely intended for interpolations.
However, unlike the approaches discussed previously, while limited to
relatively few systems, these equations are extremely accurate for the range
of conditions in which they were regressed. As such, if students ever need to
model a substance for which there already is an empirical equation of state,
they might be better off using this, instead of one of the approaches
discussed hitherto. As there is a large range of these models, we only provide
two examples: IAPWS-95[48], for water, and GERG-2008[45], for liquified
natural-gas systems. We use these to illustrate the only real limitation of
these models: extrapolation beyond the conditions over which they were
parameterised. In short, this subsection is intended to make students aware
that these empirical equations of state exist and in what cases they might
wish to use them.

###  4.3 Section 3: Employing thermodynamic models

Having now developed a good understanding of the equations of state and
thermodynamic models, in this section students are introduced to the
computational methods used to obtain more-familiar and useful properties. One
of the changes made in content presentation and delivery in this section,
relative to section two, is that larger code blocks are left visible for the
students to examine, given that the focus is now on the numerical
implementation rather than just a high-level understanding.

####  4.3.1 Volume Solvers

(a) Example of a $(p,V)$ isotherm (pink dot-dashed line) with three possible
roots (purple markers). The pressure construction is denoted by the black
dashed line. The vapour–liquid envelope (blue continuous curve) and critical
point (green marker) are included for reference.

(b) Convergence of the liquid volume using Newton’s method from an initial
guess (green) to the true solution (purple). The volume at each iteration is
denoted by the orange markers. The continuous blue curve denotes objective
function being solved.

Figure 6: Key illustrations used when discussing volume solvers within the
Pluto notebooks.

The simplest and perhaps most-useful method with which the students should be
familiar with is that for obtaining the volume ($V$) of a system at a given
temperature ($T_{0}$), pressure ($p_{0}$) and composition ($\mathbf{n}_{0}$):

$p(V,T_{0},\mathbf{n}_{0})=p_{0}\,.$ (6)

In the first instance, for simplicity of illustration, we focus on using the
van der Waals equation of state, although we also highlight that any any cubic
equation of state could be used. The power of cubic equations of state becomes
immediately apparent as equation 6 can be re-arranged to solve for the roots
of a cubic polynomial:

$Z^{3}-\left(1+\frac{bp_{0}}{RT_{0}}\right)Z^{2}+\left(\frac{ap_{0}}{(RT_{0})^{2}}\right)Z-\frac{abp_{0}^{2}}{(RT_{0})^{3}}=0\,,$
(7)

where $Z=p_{0}V/(n_{0}RT_{0})$ is the compressibility factor. The above
equation can be solved analytically, obtaining all three roots. In the case
where there is only one real root, no further considerations need to be made.
However, if more than one real root is obtained, the most stable root is
selected as that which corresponds to the lowest Gibbs free energy. The
results from this procedure are demonstrated within a reactive $(p,V)$ diagram
of carbon dioxide; students are invited to vary the temperature of the plotted
isotherm (as shown in figure 6(a)) and observe the resulting changes in the
diagram. The metastable and unstable sections of the isotherm are included so
as to illustrate the pressure construction (dashed black line) where the fluid
will phase split.

We note that this problem only becomes more challenging when using equations
of state with more-complex functional forms. We consider the SAFT equations of
state, which must be solved numerically. However, instead of solving equation
6, we remind students of the definition of the isothermal compressibility,
$\beta_{T}$:

$\beta_{T}=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{T}\,;$ (8)

if we assume that $\beta_{T}$ is constant, and integrate between states 1 and
2, we obtain:

$\ln{V_{2}}=\ln{V_{1}}+\beta_{T}(p_{2}-p_{1})\,.$ (9)

If one were to set $p_{2}=p_{0}$, $\beta_{T}=\beta_{T}(V_{1})$ and
$p_{1}=p(V_{1})$, then $\ln{V_{2}}$ would be our next-best guess for the
volume, starting from some initial guess $V_{1}$, resulting in a recursive
relationship:

$\ln{V_{i+1}}=\ln{V_{i}}+\beta_{T}(V_{i})\cdot(p_{0}-p(V_{i}))\,.$ (10)

We point out that this iterative scheme is equivalent to solving equation 6
with respect to the logarithm of the volume. It is preferable to solve
equation 6 with this approach as it offers greater numerical stability,
particularly in the liquid phase. However, as with any iterative scheme, we
need good initial guesses for each phase. For the vapour phase, the natural
choice would simply be the volume of an ideal gas at the same conditions:

$V_{0}=\frac{n_{0}RT_{0}}{p_{0}}\,.$ (11)

However, the situation for the liquid phase is more complicated. Thankfully,
and the reason why we choose to use the SAFT equations as an example here, we
can use physical intuition to obtain an initial guess for the liquid volume.
We know that the minimum volume taken up by a species will be at least the
summed total volume of the molecules. In the case of SAFT equations, the total
volume taken up by one mole of a pure species is given by:

$V_{\mathrm{min.}}=\frac{\pi}{6}N_{\mathrm{A}}m\sigma^{3}\,,$ (12)

where $N_{\mathrm{A}}$ is Avogadro’s number and the parameter $\sigma$
represents the diameter of each of the $m$ spherical segments comprising the
model molecule. One can then set the initial guess of the volume phase to be
some (low) multiple of this lower bound (we choose a multiple of 1.25). We do
highlight that optimising this initial guess is of importance when
accelerating the convergence of the algorithm (as demonstrated in figure
6(b)). Naturally, one must solve for both the vapour and liquid roots to then
determine which of the two phases are stable.

We finish this exercise by pointing out that for more-complicated equations of
state, such as the empirical equations of state, there are sometimes more than
three roots at a given pressure and temperature. In such cases, it is
important to ensure that all roots have been found before deciding which is
most stable, as there is no guarantee that our initial guesses will lead to
the two most-stable solutions.

####  4.3.2 Saturation solvers for pure components

One of the most-important classes of property for chemical engineers to obtain
from equations of state is the saturation properties of pure components.
Solving for the saturation conditions at a given temperature is more
complicated than the volume solvers mentioned previously as it involves
solving for the volume ($v^{\mathrm{liq.}}$ and $v^{\mathrm{vap.}}$) of each
phase using the conditions for phase equilibrium between two phases:

$\displaystyle p^{\mathrm{liq.}}-p^{\mathrm{vap.}}$ $\displaystyle=0\,,$
$\displaystyle\mu^{\mathrm{liq.}}-\mu^{\mathrm{vap.}}$ $\displaystyle=0\,,$
(13)

where the superscripts liq. and vap. denote properties relating to the liquid
and vapour phases, respectively. Inherently, this problem amounts to simply
using the Newton–Raphson method to converge the volumes, allowing us to
illustrate important concepts when solving such problems. For example, we re-
iterate the importance of initial guesses, using the approaches we developed
in the volume solver, or using empirical correlations. In addition, we also
highlight the value of re-scaling the variables such that the numerical solver
is more stable (including solving for the logarithm of the volumes rather than
the absolute values, and normalising equation 13 by the characteristic length
scales of pressure and energy).

Students will now be able to trace the saturation curve. However, we point out
that, for some equations of state, such as SAFT-type equations, the end-point
of the saturation curve, i.e., the critical point, isn’t known a priori; this
must be solved for numerically. We remind students of the definition for the
critical point of a pure substance:

$\displaystyle\left(\frac{\partial p(v,T)}{\partial v}\right)_{T}$
$\displaystyle=0\,,$ $\displaystyle\left(\frac{\partial^{2}p(v,T)}{\partial
v^{2}}\right)_{T}$ $\displaystyle=0\,.$ (14)

Students can re-use the tricks we highlighted in relation to the saturation
solver to solve the above system of equations. As a result, the student can
now trace the entire vapour–liquid envelope.

####  4.3.3 Chemical Stability Analysis

While the previous exercises pertained to single-component systems, for most
applications, chemical engineers will be dealing with mixture systems. The
phase space of mixtures becomes more complex as we introduce wide regions
where the system may exist in two or more phases. As such, if we wish to study
the properties of a system at given set of conditions
($p_{0},T_{0},\mathbf{z}_{0}$), we must first answer the question: “Does a
phase split occur?”. This problem is non-trivial and involves carrying out
chemical-stability analysis. The standard approach was developed by
Michelsen[49] and involves minimising what is referred to as the Gibbs tangent
plane distance ($TPD$):

$TPD(\mathbf{x})=\sum_{i}x_{i}(\mu_{i}(\mathbf{x})-\mu_{i}(\mathbf{z}_{0}))\,,$
(15)

where $x_{i}$ is the mole fraction of species $i$ in a candidate phase. The
function is demonstrated visually in figure 7(a). If $TPD<0$, then we know
that a phase split will occur.

Once again, we take the students through the process of re-defining the
problem in a way that is more numerically stable. We also take the opportunity
to illustrate that, unless a global optimiser is used, there is a chance that
the chemical-stability analysis might fail. As shown in figure 7(b), if a
local optimiser is used, conditions near the black marker would not be seen as
unstable as, despite being locally stable, there exists a tangent plane with a
lower Gibbs free energy (denoted by the dashed line in the figure).

(a) Definition of the Gibbs tangent plane distance (black line). The
continuous blue curve denotes the Gibbs free energy of mixing and the dashed
orange line denotes the test tangent plane.

(b) Example conditions where chemical-stability analysis may fail. Local
optimisers would not recognise that initial conditions denoted by the black
dot are unstable. The dashed orange curve here denotes the co-existence
tangent plane.

Figure 7: Key illustrations used when discussing chemical stability analysis
within the Pluto notebooks.

####  4.3.4 Isothermal Flash Problem

Having now established when a phase split will occur, the next step is to
determine what the various phases will be. Balancing both simplicity and
rigour, we choose to focus solely on the isothermal flash problem, the
backbone of multiple vital processes in chemical engineering.

To build up the necessary tools to perform the flash calculations, we derive
the relationship between $K$-factors and fugacity coefficients from the
definition for chemical equilibrium. From here, we impose the mass-balance
equation. As a result, we obtain the venerable Rachford–Rice[50] equation:

$f(\beta)=\sum_{i}\frac{(K_{i}-1)z_{i}}{1+\beta(K_{i}-1)}=0\,,$ (16)

where $\beta$ is the phase fraction and $K_{i}$ is the $K$-factor for species
$i$. Solving equation 16 is at the heart of almost all flash algorithms. As
such, we task the students with solving the Rachford–Rice equation. In
introducing the problem, we highlight one of the benefits of this equation: it
is monotonically decreasing. The implementation of the solver is left to the
students to fill-in, within the sample code. The exercise is presented the
same way as the previous exercise shown in figure 3.

Students will thus be able to obtain the phase fraction, for a given set of
$K$-factors, representing the first iteration of the flash algorithm. For most
real systems, $K$-factors depend on the composition in each phase.
Accordingly, once we’ve solved for the phase fraction, we must update the
$K$-factors and verify if the composition of each phase has changed
significantly. If so, we must solve the Rachford–Rice equation again with the
updated $K$-factors, leading to an iterative scheme.

While the implementation we present uses a simple successive-substitution
iterative scheme, we highlight that more-advanced flash algorithms, such as
the Michelsen flash algorithm[49], use methods such as Newton–Raphson or
Anderson Acceleration to improve the convergence (shown visually in figure
8(b)). We highlight that, while these algorithms may require fewer iterations
to convergence, it is important to consider the computational cost of each
iteration as a Newton–Raphson step can be significantly more costly than an
Anderson Acceleration step.

By this point, students will have a strong foundation to flash algorithms,
allowing them to comfortably explore the literature on this topic at their own
leisure. We provide multiple references on more-advanced flash algorithms for
them to explore independently.

(a) Rachford–Rice equation as a function of phase fraction.

(b) Convergence of flash algorithm using different iterative methods.

Figure 8: Key illustrations used when discussing chemical stability analysis
within the Pluto notebooks.

##  5 Course pilot: results and feedback

A remote pilot run of the course took place in the summer of 2022 with 16
participants. As the course is intended for students who have already had at
least one introductory thermodynamics class, both undergraduate (beyond second
year) and graduate students were invited to participate. A pre-course and a
post-course survey were delivered. Both surveys were vetted and approved by
the Ethics Research Committee at Imperial College London and were anonymised
with the assistance of the StudentShapers programme[51]. Some questions
presented in the second survey were repeats of questions presented in the
first, so as to evaluate the change in student self-efficacy. The pre-course
survey contained questions designed to determine students’ general coding
proficiency prior to taking the course. The post-course survey contained more
questions related to the students’ experience and thoughts on the course.
Open-ended questions were also added to the post-course survey in order to
allow students to fully express any feedback they might have.

Figure 9: Summary of student response when asked to qualify their confidence
in programming using different languages. Figure 10: Student rating to aspects
of programming in the Julia language which present an activation barrier.
Symbols represent individual response and solid lines represent the mean
response.

###  5.1 Introduction and code familiarisation

One of the main factors motivating the decision to use the Julia language and
Pluto notebooks for this course is to reduce the friction students encounter
when coding and/or using new software, which had been an issue in prior work
[15]. When asked at the start of the course how confident students felt in
programming with the MATLAB, Python, and Julia languages, the majority were
either not at all or only slightly confident with using Julia, as shown in
figure 9. In comparison, students were mostly quite or extremely confident in
using MATLAB and/or Python. The relative lack of confidence in using Julia is
not surprising considering that it is a newer language and many students would
not have encountered its use. In anticipation of students’ lack of
familiarity, a detailed set of instructions for installing and using the code
environment was prepared and made available on the Github repository. We also
offered a one-hour diagnostic and troubleshooting session prior to the course
– which, surprisingly, no student attended. In contrast to our previous
experience with Python, students mostly did not find either setting up the
environment or learning the syntax of Julia to be a significant activation
barrier (see figure 9). The latter may not be surprising as the syntax used in
MATLAB (the language which most students were familiar) is quite similar to
that used in Julia. The former speaks to the simplicity with which the Julia
package manager, and Pluto, abstract away the complexities typically
experienced when setting up a coding environment, demonstrating the benefits
of using Julia in an educational setting. One student commented on the
similarity between MATLAB and Julia, as well as their experience with the
Pluto notebooks:

“I enjoyed the interactive notebooks which allowed for a better learning
experience. Reminded me of MATLAB on-ramp course, which was a good
introduction/helper course.”

The lack of documentation, which was not a significant a issue with Python,
has become more of a barrier. This is primarily due to how comparatively young
Julia is as a language. As a result, packages in Julia often have
comparatively less comprehensive documentation and support than their
Python/MATLAB equivalents. Another common source of help for coders are online
forums like StackOverflow [52, 53] and blog posts are also similarly less
developed for Julia. Nonetheless, we believe that this issue will improve over
time considering the rapid development and adoption of Julia by many groups
such as educators and researchers. Also highlighted in figure 10 is another
limitation of introducing students to coding, which has persisted from our
prior work, is the lack of aim / problems to solve. This is something we hope
to have resolved with this and future courses we develop.

###  5.2 Content and structure

Figure 11: Student responses to various aspects of the course delivery.

(a) Student responses to their self-perceived confidence to apply concepts of
computational thermodynamics in real life.

(b) Student response to their self-perceived confidence to master concepts of
computational thermodynamics.

Figure 12: ‘Before-and-after” student responses related to their self-
perceived confidence in aspects of computational thermodynamics before and
after the course.

Overall, the feedback for the course content and structure was quite positive.
This is summarised in figure 11. As we can see, the main area of improvement
is the instructions provided prior to the course, which came in the form of an
e-mail providing information and links to the online repository where students
could obtain more details on setting-up the coding environment. Had the course
been held in-person, perhaps it would have been more instructive to formally
set-up an in-person session for this. However, we were restricted by the
remote delivery of the course. Nevertheless, as demonstrated in the previous
section, this did not seem to be a significant issue as most students
successfully set-up the coding environment and became familiarised with the
syntax quickly, despite over half ($N=8$) not using code frequently.

As evident from figure 11, the arrangement of content and purpose of each
section was very well received, highlighting the value of using progression as
a guiding principle in the course development[28]. We are particularly happy
that most students agreed that the examples given were relevant to chemical-
engineering problems. To this end, one of the more-emblematic responses we
received from students after being asked which part of the course they enjoyed
most:

“The history of EoS as it really makes you understand how far we have come.”

This was indeed one of the main motivations for the development of this
course: bringing students up-to-speed with the modern class of equations of
state available today[4]. This is made even more apparent in figure 12 where
it is illustrated that, when students were asked whether or not the course
improved their confidence with regards to mastering and applying the concepts
of computational thermodynamics, there was a near-total increase. Based on one
student’s feedback:

“…the code template made it easy to actually solve the exercises despite
lacking the necessary knowledge to set up the problem in code.”,

It appears that the integration of code templates and explanatory text was
effective in helping students develop a good understanding of computational
thermodynamics, despite not necessarily having a strong computational
background. Interestingly, one of the pieces of feedback we received from
students was that we “Could have more activity sessions…”. The number of
activity sessions (reactive exercises) was restricted to just two since we
intended to develop an accelerated introduction to computational
thermodynamics. However, to compensate for this, we embedded many examples
within the Pluto notebooks, so that students might have the opportunity to
examine the code after the course, by which time they would have a stronger
computational-thermodynamics background. Nevertheless, clearly the students
felt that, when interacting with code, exercises were more effective as they
felt more comfortable when under the supervision of the teaching staff.

Another improvement that was suggested related to the pace of the course. With
three, three-hour long lectures (including two ten-minute breaks), one per
day, it was apparent each day that by the final hour, students had a harder
time concentrating. As such, in future implementations of the course, it would
be advised to extend the duration to five or six sessions in order to spread
out the material and give more opportunities for questions and interaction
with the Pluto notebooks.

##  6 Conclusions

We have developed and rolled-out an introductory course on computational
thermodynamics for both undergraduate and graduate students using the open-
source thermo-fluids package, Clapeyron.jl, and the reactive notebooks
provided by Pluto.jl. Through the course, students have been able to develop a
sufficiently high level of competency in the use of a state-of-the-art
thermodynamics software package and a familiarity with concepts and methods in
computational thermodynamics. This result is demonstrated by the overall
increase in students’ self-efficacy and confidence. The course thus helps
position students to become informed users and also obtain a sufficient
foundation to engage with the thermodynamics literature. Two key components
enable such rapid and effective learning: first, the use of progression as a
guiding philosophy in course design where content and exercises are structured
such that content and exercises progressively reach the required level of
complexity and sophistication; second, the use of the Julia language and
reactive notebooks.

In a previous course developed by the authors [15], we used template Python
scripts for code implementation. While the use of the Python language
presented no issues due to its familiar syntax, setting up the computational
environment was a “pain point” for students and the use of static scripts
limited the opportunities for student engagement during the course. In this
course, the use of the Julia language and Pluto notebooks successfully
addressed these two issues while also not presenting any significant
challenges from the perspective of introducing a new coding language. Setup
was much more straightforward, even with none of the participants attending
the pre-course trouble-shooting session. The use of reactive notebooks also
enabled integration of course content with code in an environment that
facilitates student engagement and exploration.

This course represents an important first step in addressing a pertinent gap
in advanced thermodynamics instruction. Providing students with the
opportunity to understand how thermodynamic theory is employed in modern
applications is important for their own practice and also to accelerate the
promulgation of developments in thermodynamics to the rest of the scientific
and engineering community. The modular nature of the course content enables
self-directed learners and other educators to pick and choose components that
are relevant to them and also to easily develop new material for their own
use/courses. However, due to the comparatively small number of students who
participated in the pilot study, further evaluation is important for
evaluating whether the proposed course structure and delivery format will be
effective for a larger cohort or if a similar quality of learning experience
can be achieved when incorporating parts of the course into a larger
thermodynamics module.

## Availability of Code

All the course notes and code can be found at the following repository:
https://github.com/ClapeyronThermo/introduction-to-computational-
thermodynamics

## Acknowledgments

The authors would like to thank the students who participated in this course
and provided their valuable feedback. This project was supported by funding
from StudentShapers (Imperial College London) to enable partnership with
students.

## References

  * [1] Eric Hendriks, Georgios M. Kontogeorgis, Ralf Dohrn, Jean-Charles Hemptinne, Ioannis G. Economou, Ljudmila Fele Žilnik and Velisa Vesovic  “Industrial Requirements for Thermodynamics and Transport Properties”  In _Ind. Eng. Chem. Res._ 49.22, 2010, pp. 11131–11141  DOI: 10.1021/ie101231b
  * [2] Georgios M. Kontogeorgis, Ralf Dohrn, Ioannis G. Economou, Jean-Charles Hemptinne, Antoon Kate, Susanna Kuitunen, Miranda Mooijer, Ljudmila Fele Žilnik and Velisa Vesovic  “Industrial Requirements for Thermodynamic and Transport Properties: 2020”  In _Ind. Eng. Chem. Res._ 60.13, 2021, pp. 4987–5013  DOI: 10.1021/acs.iecr.0c05356
  * [3] Georgios M. Kontogeorgis, Antoon ten Kate, Martha Hajiw-Riberaud and Jean-Charles Hemptinne  “Conclusions from Round Table Discussion during IUT of ESAT 2021 electrolyte thermodynamics challenges - From industrial needs to academic research”  In _Fluid Phase Equilib._ 556, 2022, pp. 113399  DOI: 10.1016/j.fluid.2022.113399
  * [4] Jean-Charles Hemptinne, Georgios M. Kontogeorgis, Ralf Dohrn, Ioannis G. Economou, Antoon Kate, Susanna Kuitunen, Ljudmila Fele Žilnik, Maria Grazia De Angelis and Velisa Vesovic  “A View on the Future of Applied Thermodynamics”  In _Ind. Eng. Chem. Res._ 61.39, 2022, pp. 14664–14680  DOI: 10.1021/acs.iecr.2c01906
  * [5] Stanley I. Sandler  “Using Aspen plus in thermodynamics instruction: a step-by-step guide”  Hoboken, New Jersey: Wiley, 2015
  * [6] Christopher Martin, Joseph Ranalli and Jacob Moore  “PYroMat: A Python package for thermodynamic properties”  In _J. Open Source Softw._ 7.79, 2022, pp. 4757  DOI: 10.21105/joss.04757
  * [7] Christopher Martin, Joseph Ranalli and Jacob Moore  “Problem-based Learning Module for Teaching Thermodynamic Cycle Analysis using PYroMat”  In _2017 ASEE Annual Conference & Exposition Proceedings_ Columbus, Ohio: ASEE Conferences, 2017, pp. 28757  DOI: 10.18260/1-2–28757
  * [8] Ian H. Bell, Jorrit Wronski, Sylvain Quoilin and Vincent Lemort  “Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp”  In _Ind. Eng. Chem. Res._ 53.6, 2014, pp. 2498–2508  DOI: 10.1021/ie4033999
  * [9] M. Castier and Mohamed M. Amer  “XSEOS: An evolving tool for teaching chemical engineering thermodynamics”  In _Educ. Chem. Eng._ 6.2, 2011, pp. e62–e70  DOI: 10.1016/j.ece.2010.12.002
  * [10] Ian H. Bell, Ulrich K. Deiters and Allan M.. Leal  “Implementing an Equation of State without Derivatives: teqp”  In _Ind. Eng. Chem. Res._ 61.17, 2022, pp. 6010–6027  DOI: 10.1021/acs.iecr.2c00237
  * [11] Angel Martin, María Dolores Bermejo, Fidel A. Mato and María José Cocero  “Teaching advanced equations of state in applied thermodynamics courses using open source programs”  In _Educ. Chem. Eng._ 6.4, 2011, pp. e114–e121  DOI: 10.1016/j.ece.2011.08.003
  * [12] Pierre J. Walker, Hon-Wa Yew and Andrés Riedemann  “Clapeyron.jl: An Extensible, Open-Source Fluid Thermodynamics Toolkit”  In _Ind. Eng. Chem. Res._ 61.20, 2022, pp. 7130–7153  DOI: 10.1021/acs.iecr.2c00326
  * [13] Amir Karimi  “Use of Interactive Computer Software in Teaching Thermodynamics Fundamental Concepts”  In _Innovations in Engineering Education: Mechanical Engineering Education, Mechanical Engineering/Mechanical Engineering Technology Department Heads_ Orlando, Florida, USA: ASMEDC, 2005, pp. 235–242  DOI: 10.1115/IMECE2005-81943
  * [14] Christopher R. Martin, Jacob P. Moore and Joseph A. Ranalli  “Teaching the foundations of thermodynamics with PYro”  In _2016 IEEE Frontiers in Education Conference (FIE)_ Erie, PA, USA: IEEE, 2016, pp. 1–6  DOI: 10.1109/FIE.2016.7757589
  * [15] Pavan Inguva, Vijesh J. Bhute, Thomas N.H. Cheng and Pierre J. Walker  “Introducing students to research codes: A short course on solving partial differential equations in Python”  In _Educ. Chem. Eng._ 36, 2021, pp. 1–11  DOI: 10.1016/j.ece.2021.01.011
  * [16] Pavan Inguva, Parth Shah, Umang Shah and Clemens Brechtelsbauer  “How to Design Experiential Learning Resources for Independent Learning”  In _J. Chem. Educ._ 98.4, 2021, pp. 1182–1192  DOI: 10.1021/acs.jchemed.0c00990
  * [17] Suraj Pawar and Omer San  “CFD Julia: A Learning Module Structuring an Introductory Course on Computational Fluid Dynamics”  In _Fluids_ 4.3, 2019, pp. 159  DOI: 10.3390/fluids4030159
  * [18] Fani Boukouvala, Alexander Dowling, Jonathan Verrett, Zachary Ulissi and Victor Zavala  “Computational Notebooks in Chemical Engineering Curricula”  In _Chem. Eng. Educ._ 54.3, 2020, pp. 143–150
  * [19] Jeffrey M. Perkel  “Julia: come for the syntax, stay for the speed”  In _Nature_ 572.7767, 2019, pp. 141–142  DOI: 10.1038/d41586-019-02310-3
  * [20] Robert C. Reid, John M. Prausnitz and Bruce E. Poling  “The Properties of Gases and Liquids” ISSN: 0031-9228  New York: McGraw-Hill, 1987  DOI: 10.1063/1.3060771
  * [21] Joachim Gross and Gabriele Sadowski  “Perturbed-Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules”  In _Ind. Eng. Chem. Res._ 40.4, 2001, pp. 1244–1260  DOI: 10.1021/ie0003887
  * [22] Fons Van Der Plas, Michiel Dral, Paul Berg, Rik Huijzer, Mikołaj Bocheński, Alberto Mengali, Connor Burns, Benjamin Lungwitz, Hirumal Priyashan, Jerry Ling, Eric Zhang, Felipe S.. Schneider, Ian Weaver, Xiu-Zhe (Roger) Luo, Shuhei Kadowaki, Gabriel Wu, Jelmar Gerritsen, Rok Novosel, Supanat, Zachary Moon, Luis Müller, Tecosaur, Michael Abbott, Nicholas Bauer, Pbouffard, Satoshi Terasaki and Shashank Polasa  “Pluto.jl: v0.19.27”  Zenodo, 2023  DOI: 10.5281/ZENODO.4792401
  * [23] Jeremiah W. Johnson  “Benefits and Pitfalls of Jupyter Notebooks in the Classroom”  In _Proceedings of the 21st Annual Conference on Information Technology Education_ Virtual Event USA: ACM, 2020, pp. 32–37  DOI: 10.1145/3368308.3415397
  * [24] J.C. Domínguez, M.V. Alonso, E.J. González, M.I. Guijarro, R. Miranda, M. Oliet, V. Rigual, J.M. Toledo, M.M. Villar-Chavero and P. Yustos  “Teaching chemical engineering using Jupyter notebook: Problem generators and lecturing tools”  In _Educ. Chem. Eng._ 37, 2021, pp. 1–10  DOI: 10.1016/j.ece.2021.06.004
  * [25] J. Bascuñana, S. León, M. González-Miquel, E.J. González and J. Ramírez  “Impact of Jupyter Notebook as a tool to enhance the learning process in chemical engineering modules”  In _Educ. Chem. Eng._ 44, 2023, pp. 155–163  DOI: 10.1016/j.ece.2023.06.001
  * [26] Jeffrey M. Perkel  “Reactive, reproducible, collaborative: computational notebooks evolve”  In _Nature_ 593.7857, 2021, pp. 156–157  DOI: 10.1038/d41586-021-01174-w
  * [27] Nikolaos Kalogeropoulos, Pierre Walker, Colin Hale, Klaus Hellgardt, Andrew Macey, Umang V Shah and Marsha P Maraj  “Facilitating Independent Learning: Student Perspectives on the Value of Student-Led Maker Spaces in Engineering Education”  In _IJEE_ 36.4, 2020, pp. 14
  * [28] Mingrou Xie, Pavan Inguva, Wenqian Chen, Nicholaus Prasetya, Andrew Macey, Peter DiMaggio, Umang Shah and Clemens Brechtelsbauer  “Accelerating Students’ Learning of Chromatography with an Experiential Module on Process Development and Scaleup”  In _J. Chem. Educ._ 97.4, 2020, pp. 1001–1007  DOI: 10.1021/acs.jchemed.9b01076
  * [29] Pierre J. Walker and Andrew J. Haslam  “A New Predictive Group-Contribution Ideal-Heat-Capacity Model and Its Influence on Second-Derivative Properties Calculated Using a Free-Energy Equation of State”  In _J. Chem. Eng. Data_ 65.12, 2020, pp. 5809–5829  DOI: 10.1021/acs.jced.0c00723
  * [30] Grant M. Wilson  “Vapor-Liquid Equilibrium. XI. A New Expression for the Excess Free Energy of Mixing”  In _J. Am. Chem. Soc._ 86.2, 1964, pp. 127–130  DOI: 10.1021/ja01056a002
  * [31] Henri Renon and J.. Prausnitz  “Local compositions in thermodynamic excess functions for liquid mixtures”  In _AIChE J._ 14.1, 1968, pp. 135–144  DOI: 10.1002/aic.690140124
  * [32] Denis S. Abrams and John M. Prausnitz  “Statistical thermodynamics of liquid mixtures: a new expression for the excess gibbs energy of partly or completely miscible systems”  In _AIChE J._ 21.1, 1975, pp. 116–128  DOI: 10.1002/aic.690210115
  * [33] Aage Fredenslund, Jürgen Gmehling, Michael L. Michelsen, Peter Rasmussen and John M. Prausnitz  “Computerized Design of Multicomponent Distillation Columns Using the UNIFAC Group Contribution Method for Calculation of Activity Coefficients”  In _Ind. Eng. Chem. Process Des. Dev._ 16.4, 1977, pp. 450–462  DOI: 10.1021/i260064a004
  * [34] J.. Waals  “Molekulartheorie eines Körpers, der aus zwei verschiedenen Stoffen besteht”  In _Zeitschrift für Physikalische Chemie_ 5U.1, 1890, pp. 133–173  DOI: 10.1515/zpch-1890-0514
  * [35] Otto Redlich and J… Kwong  “On the Thermodynamics of Solutions. V. An Equation of State. Fugacities of Gaseous Solutions.”  In _Chem. Rev._ 44.1, 1949, pp. 233–244  DOI: 10.1021/cr60137a013
  * [36] Giorgio Soave  “Equilibrium constants from a modified Redlich-Kwong equation of state” Publisher: Pergamon Press  In _Chem. Eng. Sci._ 27.6, 1972, pp. 1197–1203  DOI: 10.1016/0009-2509(72)80096-4
  * [37] Ding Yu Peng and Donald B. Robinson  “A New Two-Constant Equation of State”  In _Ind. Eng. Chem. Fundam._ 15.1, 1976, pp. 59–64  DOI: 10.1021/i160057a011
  * [38] T. Holderbaum and J. Gmehling  “PSRK: A Group Contribution Equation of State Based on UNIFAC”  In _Fluid Phase Equilib._ 70.2-3, 1991, pp. 251–265  DOI: 10.1016/0378-3812(91)85038-V
  * [39] K.. Butcher and M.. Medani  “Thermodynamic properties of methanol-benzene mixtures at elevated temperatures”  In _J. Appl. Chem._ 18.4, 2007, pp. 100–107  DOI: 10.1002/jctb.5010180402
  * [40] Marie-José Huron and Jean Vidal  “New mixing rules in simple equations of state for representing vapour-liquid equilibria of strongly non-ideal mixtures”  In _Fluid Phase Equilib._ 3.4, 1979, pp. 255–271  DOI: 10.1016/0378-3812(79)80001-1
  * [41] David Shan Hill Wong and Stanley I. Sandler  “A theoretically correct mixing rule for cubic equations of state”  In _AIChE J._ 38.5, 1992, pp. 671–680  DOI: 10.1002/aic.690380505
  * [42] Jens Ahlers and Jürgen Gmehling  “Development of an universal group contribution equation of state. Prediction of liquid densities for pure compounds with a volume translated Peng-Robinson equation of state”  In _Fluid Phase Equilib._ 191.1-2, 2001, pp. 177–188  DOI: 10.1016/S0378-3812(01)00626-4
  * [43] Georgios M. Kontogeorgis, Epaminondas C. Voutsas, Iakovos V. Yakoumis and Dimitrios P. Tassios  “An equation of state for associating fluids”  In _Ind. Eng. Chem. Res._ 35.11, 1996, pp. 4310–4318  DOI: 10.1021/ie9600203
  * [44] Thomas Lafitte, Anastasia Apostolakou, Carlos Avendaño, Amparo Galindo, Claire S. Adjiman, Erich A. Müller and George Jackson  “Accurate statistical associating fluid theory for chain molecules formed from Mie segments”  In _J. Chem. Phys._ 139.15, 2013, pp. 154504  DOI: 10.1063/1.4819786
  * [45] O. Kunz and W. Wagner  “The GERG-2008 wide-range equation of state for natural gases and other mixtures: An expansion of GERG-2004”  In _J. Chem. Eng. Data_ 57.11, 2012, pp. 3032–3091  DOI: 10.1021/je300655b
  * [46] W.G. Chapman, K.E. Gubbins, G. Jackson and M. Radosz  “SAFT: Equation-of-state solution model for associating fluids”  In _Fluid Phase Equilib._ 52, 1989, pp. 31–38  DOI: 10.1016/0378-3812(89)80308-5
  * [47] Walter G. Chapman, Keith E. Gubbins, George Jackson and Maciej Radosz  “New reference equation of state for associating liquids”  In _Ind. Eng. Chem. Res._ 29.8, 1990, pp. 1709–1721  DOI: 10.1021/ie00104a021
  * [48] Wolfgang Wagner and Hans-Joachim Kretzschmar  “International Steam Tables: Properties of Water and Steam Based on the Industrial Formulation of IAPWS-IF97”  Springer, 2008
  * [49] Michael L. Michelsen  “The isothermal flash problem. Part II. Phase-split calculation”  In _Fluid Phase Equilib._ 9.1, 1982, pp. 21–40  DOI: 10.1016/0378-3812(82)85002-4
  * [50] H.H. Rachford and J.D. Rice  “Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium”  In _J. Pet. Technol._ 4.10, 1952, pp. 19–3  DOI: 10.2118/952327-G
  * [51] “StudentShapers - Imperial College London”, https://www.imperial.ac.uk/students/studentshapers/ URL: https://www.imperial.ac.uk/students/studentshapers/
  * [52] Trishala Bhasin, Adam Murray and Margaret-Anne Storey  “Student Experiences with GitHub and Stack Overflow: An Exploratory Study”  In _2021 IEEE/ACM 13th International Workshop on Cooperative and Human Aspects of Software Engineering (CHASE)_ Madrid, Spain: IEEE, 2021, pp. 81–90  DOI: 10.1109/CHASE52884.2021.00017
  * [53] Pierpaolo Dondio and Suha Shaheen  “Is StackOverflow an Effective Complement to Gaining Practical Knowledge Compared to Traditional Computer Science Learning?”  In _Proceedings of the 2019 11th International Conference on Education Technology and Computers_ Amsterdam Netherlands: ACM, 2019, pp. 132–138  DOI: 10.1145/3369255.3369258


# Robust Coherent Control of Bimolecular Collisions beyond the Ultracold
Regime

Adrien Devolder1, Paul Brumer1, and Timur V. Tscherbul2 1Chemical Physics
Theory Group, Department of Chemistry, and Center for Quantum Information and
Quantum Control, University of Toronto, Toronto, Ontario, M5S 3H6, Canada
2Department of Physics, University of Nevada, Reno, NV, 89557, USA

###### Abstract

Quantum coherent control of bimolecular collisions beyond the ultracold regime
can face a major challenge due to the incoherent addition of different partial
wave contributions to the total scattering cross section. These contributions
become increasingly numerous as the collision energy increases, leading to a
loss of overall control. Here, we overcome this limitation by leveraging the
recently discovered Partial Wave Phase Locking (PWPL) effect, which
synchronizes the oscillations of all partial wave contributions. By using
rigorous quantum scattering calculations, we demonstrate that PWPL enables
coherent control of spin exchange in ion-atom collisions, far outside the
ultracold regime, even with as many as 5000 partial wave contributions. The
predicted extent of control is sufficient to be measurable in cold atom-ion
hybrid experiments.

Introduction– Two-body collisions and chemical reactions of atoms and
molecules are responsible for a wide range of phenomena in physics and
chemistry[1, 2, 3], such as energy transfer, thermalization, relaxation,
decoherence, and spectral line broadening to name a few. As these phenomena
determine the properties of dilute gases, they play a central role in atomic
and molecular spectroscopy [4], atmospheric science [5], astrochemistry [6],
and ultracold chemistry [7, 8]. For this reason, controlling the quantum
dynamics of two-body collisions has long been a major thrust of physics and
chemistry, and led to the development of several vast fields of research,
including coherent control [9], laser control of chemical reactions [10],
mode-selective chemistry [1] and stereochemistry [11, 12].

A key challenge in controlling binary collisions and chemical reactions lies
in the random nature of scattering events. Specifically, under ambient
conditions, the integral collision cross section (or reaction rate) is
determined by many partial wave contributions, which are essentially random
functions of $\ell$, the orbital angular momentum for the collision [13, 1].
Quantum control protocols target a given partial wave contribution and rely in
the phase of the underlying scattering amplitude (or S-matrix). Thus, the
optimal values of the control parameter (be it the value of an external field
in field-based control schemes, laser pulse parameters in optimal control, or
superposition parameters in coherent control) are necessarily
$\ell$-dependent, and thus cannot be optimized for all values of $\ell$
contributing to the integral scattering cross section. This fundamental issue,
which we will refer to as "partial wave scrambling" has generally prevented
the application of quantum control techniques to collisions [14] and the
observation of scattering resonances [15] in the multiple partial wave regime.

A common way to combat partial wave scrambling is to cool the colliding
species down to the ultracold regime, where collisions are dominated by a
single initial partial wave. This eliminates the scrambling in the incident
collision channel and reduces its severity in the outgoing channels, allowing
for a high degree of control. The current state-of-art control techniques rely
on cold molecules/atoms prepared in a well-defined internal state interacting
with electric [16, 17, 18], microwave [19, 20, 21, 22, 23], optical [24], or
magnetic fields [25, 26]. However, the use of external fields may not be
suitable for some applications, especially when the collision partners lack
electric and/or magnetic dipole moments [27].

Coherent control is an attractive method that does not rely on external
fields. This technique involves preparing superpositions of the internal
states of colliding particles to create interference effects that can be
manipulated by changing the relative phase between the states [9, 28, 29].
Complete coherent control is possible over ultracold resonant exchange
processes, such as spin, charge, or excitation exchange, where only a single
partial wave is involved in both the incident and final collision channels
[27]. These processes can be completely suppressed (or activated), via
destructive (constructive) interference. By contrast, coherent control in the
multiple partial wave regime can face a major challenge due to the partial
wave scrambling.

Here, we show that efficient coherent control in the multiple partial wave
regime can be achieved using the partial wave phase locking (PWPL) effect [30,
31, 32, 33], which manifests in a coherent addition of different partial wave
contribution. The physical origin of the PWPL effect can be attributed to the
short-range nature of the spin-exchange interaction and the small magnitude of
the centrifugal kinetic energy compared to the well depth of the interaction
potential [34], enhancing quantum interference in the multiple partial wave
regime. Here, we show this robust PWPL-assisted coherent control of spin
exchange in ion-atom collisions (Sr+-Rb). Cold atom-ion hybrid systems have
been realized experimentally as a promising platform for quantum science [32].
Hence our results can be readily verified in the laboratory. To our knowledge,
this is the first approach to control collisions beyond the ultracold regime
and can be applied to a wide range of quasiresonant process [33].

Initial superposition and coherent control of cross section– Consider a binary
collision $A+B$, where $A$ and $B$ denote atoms or molecules initially
prepared in a coherent superposition of internal angular momentum states :

$\ket{\psi_{A}}=N\left(\sqrt{\cos\eta}\ket{j_{A},m_{1A}}+\sqrt{\sin\eta}e^{i\frac{\beta}{2}}\ket{j_{A},m_{2A}}\right),$
(1)
$\ket{\psi_{B}}=N\left(\sqrt{\sin\eta}e^{i\frac{\beta}{2}}\ket{j_{B},m_{1B}}+\sqrt{\cos\eta}\ket{j_{B},m_{2B}}\right),$
(2)

where $N=\frac{1}{\sqrt{\sin\eta+\cos\eta}}$ is a normalization factor,
$\eta\in[0,\pi/2]$ and $\beta\in[0,2\pi]$ are the parameters that determine
the relative population and phase of the superposition, $j_{A}$ and $j_{B}$
are the internal angular momenta of the colliding partners $A$ and $B$, and
$m_{1A}$, $m_{2A}$, $m_{1B}$, and $m_{2B}$ are the corresponding projections
on the space-fixed quantization axis $Z$, subject to the constraint
$m_{1A}+m_{2B}=m_{1B}+m_{2A}$ imposed by rotational symmetry, which is
required to obtain interference [35]. The initial state for the collision is
given by the product $\ket{\psi_{A}}\ket{\psi_{B}}$:

$\begin{split}\ket{\Psi_{sup}}=N^{2}\Bigg{(}\cos\eta\ket{m_{1A};m_{2B}}+\sin\eta
e^{i\beta}\ket{m_{2A};m_{1B}}\\\
+\sqrt{\cos\eta\sin\eta}e^{i\frac{\beta}{2}}(\ket{m_{1A};m_{1B}}+\ket{m_{2A};m_{2B}})\Bigg{)},\end{split}$
(3)

where we have defined
$\ket{m_{1A};m_{2B}}\equiv\ket{j_{A},m_{1A}}\ket{j_{B},m_{2A}}$ etc, for
brevity. Due to rotational symmetry, the states $\ket{m_{1A};m_{2B}}$ and
$\ket{m_{2A};m_{1B}}$ interfere with each other, while the states
$\ket{m_{1A};m_{1B}}$ and $\ket{m_{2A};m_{2B}}$ instead give rise to satellite
terms[9, 27]. For this reason, the superposition,

$\ket{\Psi_{ent}}=\cos\eta\ket{m_{1A};m_{2B}}+\sin\eta
e^{i\beta}\ket{m_{2A};m_{1B}},$ (4)

provides better control, as previously demonstrated [36]. While this
superposition is harder to prepare experimentally than that in eq. (3) because
the colliding partners are entangled, it provides a useful reference point for
comparing different control schemes, as shown below.

The cross section from the non-entangled superposition (3) to a final state
$\ket{f}$ can be split into two parts: one related to the cross section from
the entangled superposition (4), and one from the satellite terms:

$\sigma_{sup\rightarrow f}(\eta,\beta)=N^{4}\left(\sigma_{ent\rightarrow
f}+\sigma_{sat\rightarrow f}\right),$ (5)

where

$\begin{split}\sigma_{ent\rightarrow
f}(\eta,\beta)=\frac{\pi}{k^{2}}\sum_{\ell,m_{\ell}}\sum_{\ell^{\prime},m^{\prime}_{\ell}}\Big{|}&\cos\eta\
T_{m_{1A}m_{2B}\ell m_{\ell}\rightarrow f\ell^{\prime}m^{\prime}_{\ell}}\\\
+&\sin\eta\ T_{m_{2A}m_{1B}\ell m_{\ell}\rightarrow
f\ell^{\prime}m^{\prime}_{\ell}}\Big{|}^{2},\end{split}$ (6)

and

$\sigma_{sat\rightarrow
f}(\eta)=\cos\eta\sin\eta\left(\sigma_{m_{1A};m_{1B}\rightarrow
f}+\sigma_{m_{2A};m_{2B}\rightarrow f}\right).$ (7)

Here, $\ell$ and $\ell^{\prime}$ are the initial and final orbital angular
momenta for the collision, while $m_{\ell}$ and $m^{\prime}_{\ell}$ are the
projections of $\ell$ and $\ell^{\prime}$ on a space-fixed quantization axis
$Z$, $T_{m_{1A}m_{2B}\ell m_{\ell}\rightarrow
f\ell^{\prime}m^{\prime}_{\ell}}$ and $T_{m_{2A}m_{1B}\ell m_{\ell}\rightarrow
f\ell^{\prime}m^{\prime}_{\ell}}$ are the $T$-matrix elements associated with
the initial states $\ket{m_{1A};m_{2B}}$ and $\ket{m_{2A};m_{1B}}$, and $k$ is
the relative wavevector.

To achieve good control of $\sigma_{sup\rightarrow f}$, it is necessary to
exert efficient control over $\sigma_{ent\rightarrow f}$. However, even with
good control over the latter, the inclusion of large satellite terms can
significantly reduce overall control. This illustrates two important
requirements for the efficient coherent control of $\sigma_{sup\rightarrow
f}$: (i) achieving the best possible control over $\sigma_{ent\rightarrow f}$
and (ii) minimizing the value of $\sigma_{sat\rightarrow f}$. To quantify
these requirements, we define two control indices. First, the extent of
control over the cross section from the entangled initial superposition
$\sigma_{ent\rightarrow f}$ can be defined as:

$R_{c,ent}=\frac{\left|\sigma_{int}\right|}{\sqrt{\sigma_{m_{1A};m_{2B}\rightarrow
f}\sigma_{m_{2A};m_{1B}\rightarrow f}}},$ (8)

where

$\sigma_{int}=\frac{\pi}{k^{2}}\sum_{\ell,m_{\ell}}\sum_{\ell^{\prime},m^{\prime}_{\ell}}T_{m_{1A}m_{2B}\ell
m_{\ell}\rightarrow f\ell^{\prime}m^{\prime}_{\ell}}T_{m_{2A}m_{1B}\ell
m_{\ell}\rightarrow f\ell^{\prime}m^{\prime}_{\ell}}^{*},$ (9)

is the interference contribution to the integral cross section, and
$\sigma_{m_{1A};m_{1B}\rightarrow f}$ and $\sigma_{m_{2A};m_{2B}\rightarrow
f}$ are the cross sections from the states $\ket{m_{1A};m_{1B}}$ and
$\ket{m_{2A};m_{2B}}$ respectively, to the final state $\ket{f}$. The
definition relies on the Schwartz inequality, so that $R_{c,ent}$ lies between
zero and one. Second, the effect of the satellite terms is quantified by:

$R_{sat}=\frac{min\left(\sigma_{m_{1A};m_{2B}\rightarrow
f},\sigma_{m_{2A};m_{1B}\rightarrow
f}\right)}{max\left(\sigma_{m_{1A};m_{1B}\rightarrow
f},\sigma_{m_{2A};m_{2B}\rightarrow f}\right)}.$ (10)

Here, $min\left(\sigma_{m_{1A};m_{2B}\rightarrow
f},\sigma_{m_{2A};m_{1B}\rightarrow f}\right)$ and
($max\left(\sigma_{m_{1A};m_{1B}\rightarrow
f},\sigma_{m_{2A};m_{2B}\rightarrow f}\right)$) is the smallest (largest) of
the two values. Small satellite terms correspond to $R_{sat}>>1$, a favorable
condition for coherent control.

Finally, we define a global control index as $R_{c,sup}=max_{\eta}(V)$, where
$V$ represents the visibility that measures the oscillation of
$\sigma_{sup\rightarrow f}(\eta,\beta)$ when only the relative phase $\beta$
is varied:

$V(\tilde{\eta})=\frac{\sigma_{sup\rightarrow
f}(\tilde{\eta},\beta_{max}^{\tilde{\eta}})-\sigma_{sup\rightarrow
f}(\tilde{\eta},\beta_{min}^{\tilde{\eta}})}{\sigma_{sup\rightarrow
f}(\tilde{\eta},\beta_{max}^{\tilde{\eta}})+\sigma_{sup\rightarrow
f}(\tilde{\eta},\beta_{min}^{\tilde{\eta}})},$ (11)

where $\beta_{min}^{\tilde{\eta}}$ ($\beta_{max}^{\tilde{\eta}})$ is the value
of $\beta$ for which the cross-section is minimal (maximal) when
$\eta=\tilde{\eta}$. Like $R_{c,ent}$, the values of $R_{c,sup}$ are bound
between 0 and 1.

Partial wave scrambling in coherent control- The control of the entangled
cross section $\sigma_{ent\rightarrow f}$ is limited by incoherent addition of
the initial and final partial waves $(\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell})$ in Eq. (6). For a given
superposition determined by the parameters ($\eta$ and $\beta$), some partial
wave contributions $(\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell})$ may
experience constructive interference while others may exhibit destructive
interference, suppressing the interference term in Eq. (9) [37]. This partial
wave scrambling issue becomes more significant as the collision energy
increases, along with the number of $(\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell})$ contributions, and can result in
complete loss of control.

To quantify partial wave scrambling, we examine the distribution of the
superposition parameters ($\eta,\beta$) for which each partial wave
contribution reaches a minimum value, $\eta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}=\arctan\left(\frac{\left|T_{m_{1A}m_{2B}\ell
m_{\ell}\rightarrow
f\ell^{\prime}m^{\prime}_{\ell}}\right|}{\left|T_{m_{2A}m_{1B}\ell
m_{\ell}\rightarrow f\ell^{\prime}m^{\prime}_{\ell}}\right|}\right)$ and
$\beta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}=arg\left(T_{m_{1A}m_{2B}\ell
m_{\ell}\rightarrow
f\ell^{\prime},m^{\prime}_{\ell}}\right)-arg\left(T_{m_{2A}m_{1B}\ell
m_{\ell}\rightarrow f\ell^{\prime},m^{\prime}_{\ell}}\right)$ depending on the
ratio of the magnitudes, and on the difference of phases (arguments) of the
$T$-matrix elements, respectively. Note that the maximum parameters are simply
obtained from: $\eta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}+\eta_{max}^{\ell,m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}=\pi/2$
and $\beta_{max}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}-\beta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}=\pi$. The distribution of the optimal
parameters $\eta_{min}^{\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$ and
$\beta_{min}^{\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$ determines the
degree of partial wave scrambling. A random distribution leads to a rapid
decrease in control as the number $N_{\ell}$ of significant partial wave
contributions $(\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell})$ increases,
because $R_{c,ent}$ scales as $1/\sqrt{N_{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}}$ [37]. The solution of the partial
wave scrambling problem lies in the clustering of the optimal parameters
$\eta_{min}^{\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$ and
$\beta_{min}^{\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$. As shown below,
the clustering of $\beta_{min}^{\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$
is manifest in the PWPL mechanism, which dramatically reduces the partial wave
scrambling and thereby paves the way to coherent control in the multiple
partial wave regime.

Application: Coherent control of spin relaxation in Sr+-Rb collisions– The
PWPL phenomenon was first predicted to occur in spin relaxation in ion-atom
collisions, making hybrid ion-atom systems [32] ideal for investigating PWPL
assisted coherent control. Consider a rubidium atom (87Rb) prepared in a
superposition of hyperfine states $\ket{2,-1}_{B}$ and $\ket{2,0}_{B}$,
colliding with a trapped strontium ion 88Sr+ prepared in a superposition of
Zeeman states $\ket{1/2,-1/2}_{A}$ and $\ket{1/2,1/2}_{A}$. We present the
results for relaxation to the final states
$\ket{1/2,-1/2}_{A}\ket{1,0}_{B}\equiv\ket{\downarrow}$ and
$\ket{1/2,+1/2}_{A}\ket{1,0}_{B}\equiv\ket{\uparrow}$, as these states have
larger cross-sections than the final states $\ket{1/2,\pm
1/2}_{A}\ket{1,-1}_{B}$ and $\ket{1/2,\pm 1/2}_{A}\ket{1,1}_{B}$.

To motivate experimental studies we carried out rigorous coupled-channel (CC)
calculations of Sr+-Rb collisions using state of the art ab-initio interaction
potentials and second-order spin-orbit interactions, as described in [30]. To
ensure numerical convergence of the results for collision energies ranging
from 1 $\mu$K to 50 mK, we used extended CC basis sets including up to 80
partial waves.

Figure 1: (a) Control indices, $R_{c,ent}$ and $R_{c,sup}$, for the transition
from the entangled and non-entangled superpositions to the final states
$\ket{\downarrow}$ and $\ket{\uparrow}$ . (b) Satellite term index $R_{sat}$
for the final states $\ket{\downarrow}$ and $\ket{\uparrow}$

The control indices, $R_{c,ent}$ and $R_{c,sup}$, for Sr+-Rb calculated from
exact CC results, are shown in Fig. 1 (a). The high value of the entangled
control index, $R_{c,ent}$, demonstrates that efficient control is possible in
the multiple partial-wave regime. For the final state $\ket{\uparrow}$,
$R_{c,ent}$ is close to 1, indicating complete control, whereas for the final
state $\ket{\downarrow}$, $R_{c,ent}$ is around 0.5-0.6. Remarkably, the high
entangled control index remains mostly independent of collision energy,
indicating robust control over a wide energy range. Variation of the control
can be caused by the presence of resonances, for example at $E$= 100 $\mu$K,
400 $\mu$K, and 5 mK.

As expected, the control index for the non-entangled superposition,
$R_{c,sup}$, is affected by satellite terms, which can be quantified using the
parameters $R_{sat}$ (eq. 10) shown in Fig. 1 (b). For the final state
$\ket{\uparrow}$, the satellite terms are significant, resulting in a large
difference between $R_{c,ent}$ and $R_{c,sup}$, with $R_{c,sup}$ being around
0.2-0.3. In contrast, the satellite terms are small for the final state
$\ket{\downarrow}$, resulting in a small difference between $R_{c,ent}$ and
$R_{c,sup}$, with $R_{c,sup}\approx$ 0.4-0.5. The reason for the smaller
impact of the satellite terms on the final state $\ket{\downarrow}$ is that
for this state, the interfering transitions conserve the total internal
angular momentum projection, $m_{A}^{i}+m_{B}^{i}=m_{A}^{f}+m_{B}^{f}$,
whereas the transitions in the satellite terms do not. The situation is the
opposite for the final state $\ket{\uparrow}$. These results highlight a
complex trade-off between the partial wave scrambling and satellite terms.
Even though partial wave scrambling is less significant for the final state
$\ket{\uparrow}$, the overall control is better for the final state
$\ket{\downarrow}$ due to the smaller effect of the satellite terms.

Figure 2: Coherent control of scattering cross sections and branching ratios
for Sr+-Rb collisions with 5000 partial wave contributions $\left(\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}\right)$. Cross sections (a) from
$\ket{\Psi_{ent}}$ to $\ket{\downarrow}$, (b) from $\ket{\Psi_{sup}}$ to
$\ket{\downarrow}$, (c) from $\ket{\Psi_{ent}}$ to $\ket{\uparrow}$ and (d)
from $\ket{\Psi_{sup}}$ to $\ket{\uparrow}$. Branching ratios
$\sigma_{\downarrow}/\sigma_{\uparrow}$ (e) from $\ket{\Psi_{ent}}$ and (f)
from $\ket{\Psi_{sup}}$. The parameter $\eta$ is kept constant, $\eta=\pi/2$
for $\sigma_{\downarrow}$, while $\eta=21\pi/32$ for $\sigma_{\uparrow}$ and
the branching ratios.

To illustrate coherent control in the multiple partial wave regime, we
consider the experimentally realistic case of Sr+-Rb collisions at 50 mK [30].
At this collision energy, approximately 5000 $(\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell})$ partial wave states contribute to
the cross sections. Figures 2(a)-(d) illustrate coherent control of the cross
section by varying the phase angle $\beta$ of the initial superposition (4).
To obtain the best visibility $V$, the value of $\eta$ is fixed at
$\eta=\pi/2$ and $21\pi/32$ for the final states $\ket{\downarrow}$ and
$\ket{\uparrow}$, respectively. We observe a remarkable instance of complete
control of scattering from the entangled superposition to the final state
$\ket{\uparrow}$, showing near vanishing of the cross section due to
destructive interference (see Fig. 2 (c)). The minimum value of the cross-
section is 6.2 a.u., which is three orders of magnitude smaller than the
maximum value of 2517.3 a.u. This complete control is made possible by the
clustering of the optimal control parameters
$\eta_{min}^{\ell,m_{\ell},\ell^{\prime},m^{\prime}\ell}$ and
$\beta_{min}^{\ell,m_{\ell},\ell^{\prime},m^{\prime}\ell}$ [see Fig. 3(c) and
(d)], caused by PWPL. Note that if these parameters were randomly distributed,
the entangled control index $R_{c,ent}$ would be equal to $\sim
1/\sqrt{5000}=0.01$, significantly smaller than the calculated value. When the
initial superposition is non-entangled, complete destructive interference is
countered by the presence of satellite terms, resulting in a variation from
1926.6 to 3350.2 a.u, which remains significant and certainly large enough to
be detected in modern hybrid trapped ion-atom collision experiments [30, 31,
32, 33].

Figure 3: Phase variation of the 5000 partial wave contributions $\left(\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}\right)$ to the cross sections for the
final states (a) $\ket{\downarrow}$ and (b) $\ket{\uparrow}$. The value of
$\eta$ is fixed at $\eta=\pi/2$ and $21\pi/32$ for the final states
$\ket{\downarrow}$ and $\ket{\uparrow}$, respectively. All partial wave
contributions are plotted to illustrate the phase locking, with less
efficiency observed for the final state $\ket{\downarrow}$.

For the final state $\ket{\downarrow}$, the clustering of
$\beta_{min}^{\ell,m_{\ell},\ell^{\prime},m^{\prime}\ell}$ is less effective
than for $\ket{\uparrow}$, as shown in Fig.3 (a) and 4 (a). More detrimential
for the control, the distribution of the other optimal parameter,
$\eta_{min}^{\ell m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$ [see Fig. 4 (b)]
is broad, highlighting the importance of the distribution of $\eta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$. The PWPL effect locks the
difference of T-matrix phases but not the ratio of their magnitudes, and thus
does not guarantee the clustering of $\eta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$. While this limitation can prevent
PWPL from completely solving the problem of partial wave scrambling, we
observe that despite the relatively broad distribution of ratios, a good
degree of control is still achievable in the multiple partial wave regime,
with $\sigma_{sup\rightarrow f}(\eta,\beta)$ ranging from 2085.5 to 6131.7
a.u. The inclusion of the satellite terms has a slight impact on the control,
resulting in a variation from 1827.7 to 3850.7 a.u.

Figure 4: Distribution of the optimal parameters $\eta_{min}^{\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell}}$ and
$\beta_{min}^{\ell,m_{\ell},\ell^{\prime},m^{\prime}\ell}$ for the final
states ((a) and (b)) $\ket{\downarrow}$ and ((c) and (d)) $\ket{\uparrow}$.
Each point corresponds to a pair $(\ell
m_{\ell},\ell^{\prime}m^{\prime}_{\ell})$. The weight $w$ is calculated by the
sum of squared magnitudes of T-matrix elements: $w=\left|T_{m_{1A}m_{2B}\ell
m_{\ell}\rightarrow
f\ell^{\prime}m^{\prime}_{\ell}}\right|^{2}+\left|T_{m_{2A}m_{1B}\ell
m_{\ell}\rightarrow f\ell^{\prime}m^{\prime}_{\ell}}\right|^{2}$

As noted above, the control of the scattering to the final states
$\ket{\downarrow}$ and $\ket{\uparrow}$ is optimized in different regions of
parameter space allowing for efficient control of the corresponding branching
ratio $\sigma_{\downarrow}/\sigma_{\uparrow}$, as shown in Fig. 2(e)-(f) at
$\eta=21\pi/32$. With the entangled superposition Eq. (4), extremely robust
control is achieved, with $\sigma_{\downarrow}/\sigma_{\uparrow}$ varying by
three orders of magnitude (from 0.94941 to 960.66), thanks to the complete
control of $\sigma_{\uparrow}$. When using the non-entangled superposition
(3), the presence of satellite terms limits the extent of control over the
branching ratio to 0.66-2.26, which is large enough to be experimentally
measurable.

Conclusion- In summary, we have shown that partial wave phase-locking enables
coherent control in the multiple partial wave regime via a dramatic reduction
of partial wave scrambling. The clustering of the optimal superposition
parameters
$\beta_{\text{min}}^{\ell,m_{\ell},\ell^{\prime},m_{\ell}^{\prime}}$ enabled
by PWPL allows for the synchronized control of different partial wave
contributions to the total scattering cross section. In cases where the
distribution of the optimal control parameters
$\eta_{\text{min}}^{\ell,m_{\ell},\ell^{\prime},m_{\ell}^{\prime}}$ is broad,
such as for the final state $\ket{\downarrow}$, partial-wave scrambling is
only partially eliminated by PWPL. Even though the satellite terms reduce the
control with non-entangled superpositions, our rigorous CC calculations show
that coherent control over state-to-state integral cross sections and of the
branching ratios is significant and measurable. Therefore, collisions between
87Rb and 88Sr+, observed in a series of recent experiments [30, 31, 32],
appear to be ideal for the first experimental observation of coherent control
of two-body scattering outside of the ultracold domain. As the PWPL phenomenom
was shown to apply to any quasiresonant scattering process [33], a wide range
of these processes could soon become amenable to robust coherent control.
Furthermore, the extreme sensitivity of coherent control to PWPL effect
implies that the study of collisions of atoms and molecules prepared in
coherent superposition of internal states provides an ideal approach for
investigating the PWPL effect.

This work was supported by the U.S. Air Force Office of Scientific Research
(AFOSR) under Contract No.FA9550-22-1-0361.

## References

  * Levine [2005] R. D. Levine,  _Molecular Reaction Dynamics_ (Cambridge University Press, 2005).
  * May and Kuhn [2011] V. May and O. Kuhn,  _Charge and Energy Transfer Dynamics in Molecular Systems_ (Wiley-VCH, 2011).
  * Nitzan [2006] A. Nitzan,  _Chemical Dynamics in Condensed Phases_ (Oxford University Press, 2006).
  * Hartmann  _et al._ [2008] J.-M. Hartmann, C. Boulet, and D. Robert,  _Collisional Effects on Molecular Spectra_ (Elsevier, 2008).
  * Massey and Bates [1982] H. S. W. Massey and D. R. Bates,  _Applied Atomic Collision physics: Atmospheric Physics and Chemistry_ (Academic Press, 1982).
  * Smith [2011] I. W. M. Smith, Laboratory astrochemistry: Gas-phase processes, Annu. Rev. Astron. Astrophus. 49, 29 (2011).
  * Krems [2019] R. Krems,  _Molecules in Electromagnetic Fields_ (Wiley, 2019).
  * Dulieu and Osterwalder [2017] O. Dulieu and A. Osterwalder,  _Cold Chemistry: Molecular Scattering and Reactivity Near Absolute Zero_ (Royal Society of Chemistry, 2017).
  * Shapiro and Brumer [2012] M. Shapiro and P. Brumer,  _Quantum Control of Molecular Processes_ (Wiley-VCH, Weinheim, Germany, 2012).
  * Rice and Zhao [2000] S. Rice and M. Zhao,  _Optical Control of Molecular Dynamics_ (Wiley, 2000).
  * Alexander  _et al._ [2011] A. J. Alexander, M. Brouard, K. Kalogerakis, and J. P. Simons, Chemistry with a sense of direction—the stereodynamics of bimolecular reactions, Chem. Soc. Rev. 27, 405 (2011).
  * Orr-Erwing [1996] A. J. Orr-Erwing, Dynamical stereochemistry of bimolecular reactions, J. Chem. Soc.,Faraday Trans. 92, 881 (1996).
  * Mott and Massey [1965] N. F. Mott and H. S. W. Massey,  _The Theory of Atomic Collisions_ (Oxford University Press, 1965).
  * Krems [2008] R. Krems, Cold controlled chemistry, Phys. Chem. Chem. Phys. 10, 4079 (2008).
  * Skodje  _et al._ [2000] R. T. Skodje, D. Skouteris, D. Manolopoulos, S.-H. Lee, F. Dong, and K. Liu, Observation of a transition state resonance in the integral cross section of the F+HD reaction, J. Chem. Phys. 112, 4536 (2000).
  * Wang and Quéméner [2015] G. Wang and G. Quéméner, Tuning ultracold collisions of excited rotational dipolar molecules, New J. Phys. 17, 035015 (2015).
  * Li  _et al._ [2021] J.-R. Li, W. G. Tobias, K. Matsuda, C. Miller, G. Valtolina, L. D. Marco, R. R. W. Wang, L. Lassablière, G. Quéméner, J. L. Bohn, and J. Ye, Tuning of dipolar interactions and evaporative cooling in a three dimensional molecular quantum gas, Nat. Phys. 17, 1144 (2021).
  * Quéméner  _et al._ [2022] G. Quéméner, J. F. E. Croft, and J. L. Bohn, Electric field dependence of complex-dominated ultracold molecular collisions, Phys. Rev. A 105, 013310 (2022).
  * Karman and Hutson [2018] T. Karman and J. M. Hutson, Microwave shielding of ultracold polar molecules, Phys. Rev. Lett. 121, 163401 (2018).
  * Lassablière and Quéméner [2018] L. Lassablière and G. Quéméner, Controlling the scattering length of ultracold dipolar molecules, Phys. Rev. Lett. 121, 163402 (2018).
  * Anderegg  _et al._ [2021] L. Anderegg, S. Burchesky, S. S. Yu, T. Karman, E. Chae, K.-K. Ni, W. Ketterle, and J. M. Doyle, Observation of microwave shielding of ultracold molecules, Science 373, 779 (2021).
  * Schindewolf  _et al._ [2022] A. Schindewolf, R. Bause, X.-Y. Chen, M. Duda, T. Karman, I. Bloch, and X.-Y. Luo, Evaporation of microwave-shielded polar molecules to quantum degeneracy, Nature 607, 677 (2022).
  * Chen  _et al._ [2023] X.-Y. Chen, A. Schindewolf, S. Eppelt, R. Bause, M. Duda, S. Biswas, T. Karman, T. Hilker, I. Bloch, and X.-Y. Luo, Field-linked resonances of polar molecules, Nature 614, 59 (2023).
  * Xie  _et al._ [2020] M. L. T. Xie, R. Vexiau, A. Orban, O. Dulieu, and N. Bouloufa-Maafa, Optical shielding of destructive chemical reactions between ultracold ground-state NaRb molecules, Phys. Rev. Lett. 125, 153202 (2020).
  * Yang  _et al._ [2019] H. Yang, D.-C. Zhang, L. Liu, Y.-X.Liu, J. Nan, B. Zhao, and J.-W. Pan, Observation of magnetically tunable Feshbach resonances in ultracold 23Na40K + 40K collisions, Science 363, 261 (2019).
  * Tscherbul and Klos [2020] T. V. Tscherbul and J. Klos, Magnetic tuning of ultracold barrierless chemical reactions, Phys. Rev. Res 2, 013117 (2020).
  * Devolder  _et al._ [2021] A. Devolder, P. Brumer, and T. V. Tscherbul, Complete quantum coherent control of ultracold molecular collisions, Phys. Rev. Lett. 126, 153403 (2021).
  * Shapiro and Brumer [1996] M. Shapiro and P. Brumer, Coherent control of collisional events: Bimolecular reactive scattering, Phys. Rev. Lett. 77, 2574 (1996).
  * Gong  _et al._ [2003] J. B. Gong, M. Shapiro, and P. Brumer, Entanglement-assisted coherent control in nonreactive diatom-diatom scattering, J. Chem. Phys. 118, 2626 (2003).
  * Sikorsky  _et al._ [2018] T. Sikorsky, M. Morita, Z. Meir, A. A. Buchachenko, R. Ben-Shlomi, N. Akerman, E. Narevicius, T. V. Tscherbul, and R. Ozeri, Phase locking between different partial waves in atom-ion spin-exchange collisions, Phys. Rev. Lett. 121, 173402 (2018).
  * Kartz  _et al._ [2022] O. Kartz, M. Pinkas, N. Akerman, and R. Ozeri, Quantum logic detection of collisions between single atom–ion pairs, Nat. Phys. 18, 533 (2022).
  * Tomza  _et al._ [2019] M. Tomza, K. Jachymski, R. Gerritsma, A. Negretti, T. Calarco, Z. Idziaszek, and P. S. Julienne, Cold hybrid ion-atom systems, Rev. Mod. Phys. 91, 035001 (2019).
  * Coté and Simbotin [2018] R. Coté and I. Simbotin, Signature of the s-wave regime high above ultralow temperatures, Phys. Rev. Lett. 121, 173401 (2018).
  * Gao [2001] B. Gao, Angular-momentum-insensitive quantum-defect theory for diatomic systems, Phys. Rev. A 64, 010701(R) (2001).
  * Omiste  _et al._ [2018] J. J. Omiste, J. Floß, and P. Brumer, Coherent control of penning and associative ionization: Insights from symmetries, Phys. Rev. Lett. 121, 163405 (2018).
  * Devolder  _et al._ [2022] A. Devolder, T. V. Tscherbul, and P. Brumer, Coherent multichannel optical theorem: Quantum control of the total scattering cross section, Phys. Rev. A 105, 052808 (2022).
  * Devolder  _et al._ [2023] A. Devolder, T. V. Tscherbul, and P. Brumer, Coherent control of ultracold molecular collisions: The role of resonances, J. Phys. Chem. Lett. 14, 2171 (2023).