Coordination and thermodynamic properties of aqueous protactinium(V) by first-
principle calculations†

Hanna Oher,a,b Jeremy Delafoulhouze,c Eric Renault,c Valérie Vallet,d and Rémi
Maurice∗a,b

Protactinium ($Z$ = 91) is a very rare actinide with peculiar physico-chemical
properties. Indeed, although one may naively think that it behaves similarly
to either thorium or uranium by its position in the periodic table, it may in
fact follow its own rules. Because of the quite small energy gap between its
valence shells (in particular the 5$f$ and 6$d$ ones) and also the strong
influence of relativistic effects on its properties, it is actually a
challenging element for theoretical chemists. In this article, we combine
experimental information, chemical arguments and standard first-principle
calculations, complemented by implicit and explicit solvation, to revisit the
stepwise complexation of aqueous protactinium(V) with sulfate and oxalate
dianionic ligands (SO4^2- and C2O4^2-, respectively). From a methodological
viewpoint, we notably conclude that it is necessary to at least saturate the
coordination sphere of protactinium(V) to reach converged equilibrium constant
values. Furthermore, in the case of single complexations (i.e. with one
sulfate or oxalate ligand bound in the bidentate fashion), we show that it is
necessary to maintain the coordination of one hydroxyl group, present in the
supposed [PaO(OH)]^2+ precursor, to obtain coherent complexation constants.
Therefore, we predict that this hydroxyl group is maintained in the formation
of 1:1 complexes while we confirm that it is withdrawn when coordinating three
sulfate or oxalate ligands. Finally, we stress that this work is a first step
toward the future use of theoretical predictions to elucidate the enigmatic
chemistry of protactinium in solution.

††footnotetext: a Subatech, UMR CNRS 6457, IN2P3/IMT Atlantique/Université de
Nantes, 4 rue A. Kastler, 44307 Nantes Cedex 3, France††footnotetext: b Univ
Rennes, CNRS, ISCR (Institut des Sciences Chimiques de Rennes) – UMR 6226,
F-35000 Rennes, France; E-mail<EMAIL_ADDRESS>c
Nantes Université, CNRS, CEISAM UMR 6230, F-44000 Nantes,
France††footnotetext: d Univ. Lille, CNRS, UMR 8523-PhLAM-Physique des Lasers,
Atomes et Molécules, F-59000 Lille, France††footnotetext: † Electronic
Supplementary Information (ESI) available: Preparatory study on molecular
geometries, which includes comparisons of methods, basis sets and codes.

##  1 Introduction

Protactinium ($Z=91$) is an enigmatic element for the chemists1. While it is
located in between thorium and uranium in the periodic table, it may display
unique properties that severely hinder experimental studies and
interpretations. For instance, the five valence electrons of its free atom may
spread in numerous ways on the valence 5$f$, 6$d$, 7$s$ and 7$p$ shells, which
generates a quantum chaos of states even well before ionization2, 3, 4. In
solution, protactinium chemistry is expected to be dominated by its +V
oxidation state (formally free of valence electron), with potential occurrence
under specific reducing conditions of the unstable +IV one (displaying one
unpaired 5$f$ electron in the ground state) 5. Therefore, the aqueous
chemistry of protactinium could well have been dominated by easily
identifiable and computable species.

Unfortunately, Pa(V) and Pa(IV) are prone to hydrolysis, polymerization and
precipitation and other issues such as a strong tendency to adsorb into glass,
etc. 6. Furthermore, the intrinsic radiations of its isotopes imply specific
safety regulations to be followed. Consequently, the speciation of Pa(V) and
Pa(IV) in aqueous solution is largely unknown, which further complicates the
performance of a theoretical study dedicated to the coordination and
thermodynamic properties of these; thus, one has to rely on ad hoc research
hypotheses.

In this work, we will focus on protactinium(V) since it is expected to have a
closed-shell electronic ground state in most of its associated chemical
species. Therefore, single-reference approaches such as Kohn-Sham density
functional theory (DFT) or second-order Møller-Plesset perturbation theory
(MP2) may be a priori safely used, contrary to Pa(IV) systems. Despite this
apparent simplicity in the electronic structure theory part, challenges
remain: solvation and/or hydrolysis issues must be tackled and relativistic
effects, especially the scalar ones, should be accounted for. In fact, since
little experimental data regarding speciation and complexation constants
available, the choice of the Pa(V) chemical complexes will be severely
restrained. Anyway, we will show that it is necessary to combine experimental
information, chemical arguments and outcomes of calculations to reach
conclusions.

The paper is organized as follows. First, the choice of the systems and our
working hypotheses will be explained. Then, computational details will be
given, prior to results and discussion, organized by increasing solvation
and/or chemical complexity.

##  2 Choice of the systems and working hypotheses

The chemistry of Pa(V) is characterized by a general lack of experimental
and/or computational data 7. In this work, we aim at showing that standard
quantum chemical calculations may be of help to predict the relative
stabilities of Pa(V) complexes, based on comparisons between theoretical and
experimental data. Therefore, it is necessary to make a quick survey of known
experimental data to justify our choice of the systems under study. In
particular, the available experimental data should ideally help us to
formulate relevant hypotheses concerning the nature of the complexes that are
involved in the equilibrium processes at play. After this, we will introduce
additional working hypotheses that are required to define a practical
computational strategy, systematically improvable in terms of solvation and/or
chemical models.

The identification of Pa(V) species is a difficult task. Several types of
complementary experiments may be performed, none of them being independently
conclusive. We may classify them in two main types, (i) experiments at tracer
scale and (ii) spectroscopy experiments (ponderable scale). In (i), we may
encounter methods to determine the molecular charge (ion chromatography and/or
migration experiments under an electric field such as capillary
electrophoresis) and also methods to determine equilibrium constants, based on
competition between two distinct phases (e.g. liquid-liquid extraction). In
(ii), we may find X-ray absorption data for instance, based on near-edge or
extended structures (XANES and EXAFS, respectively). By taking profit of the
X-ray absorption spectrum dependence on the molecular structure, especially in
the EXAFS regime, one may refine guessed structural models by a fitting
procedure, leading in the end to an experimental structure.

It is important to understand that equilibrium constants are determined by
fitting procedures, implying that an underlying thermodynamic model is
selected. Because the activity of water in water is set to 1, all the
reactions are defined up to a given number of water molecules. This is a
strong limitation of the liquid-liquid extraction approach. For instance,
Mendes et al. 8 have shown that Pa(V) can form a 1:1 complex with DTPA
(diethylenetriaminepentaacetic acid) and have determined the associated
formation constant for the following reaction:

$\ce{PaO(OH)^{2+}+DTPA^{5-}+3H+<=>[\beta_{app}]Pa^{V}-DTPA+xH2O}$ (1)

where Pa^V-DTPA is either PaO(H_2DTPA) or Pa(DTPA) and $x$ = 1 or $x$ = 2,
respectively. To make a long story short, the neutral charge was determined by
capillary electrophoresis, and the potentially formed complexes justified by
chemical arguments and quantum mechanical calculations. Clearly, the absence
of X-ray absorption data rules out this system for a first computational study
on equilibrium constants, since we are stuck with an ambiguity in the nature
of the formed complex.

EXAFS data of Pa(V) have been acquired in fluoric acid (HF), sulfuric acid
(H_2SO_4) and oxalic acid (H_2C_2O_4) concentrated media 9, 10, 11.
Concentrated media are interesting because we may expect to saturate the first
coordination sphere of the Pa(V) ion and end up with a simpler speciation than
at intermediate concentrations, i.e. we may observe single well-defined
chemical species. While in fluoric acid the mono-oxo bond of the [PaO(OH)]^2+
starting species is cleaved (up to eight fluorides being coordinated 11, 12),
it is observed in both the sulfuric and oxalic media. Therefore, we may expect
similar complexes to be formed. Since in the oxalic acid case, the analysis of
the data led to the firm identification of the [PaO(C_2O_4)_3]^3- complex, we
thus hypothesize the formation of the [PaO(SO_4)_3]^3- complex in sulfuric
acid medium. We recall here that our objective is to determine relative
complexation constants, meaning that it is important to compare systems of
similar natures to obtain accurate values by error cancellations, as was done
in the past for exploring the astatine ($Z$ = 85) basic chemistry 13, 14, 15,
16, 17.

Having suspected that complexes of similar natures may be formed at high
H_2SO_4 and H_2C_2O_4 concentrations, in particular of the [PaOX_3]^3- type, X
being either (SO_4) or (C_2O_4), we continue our survey with the corresponding
liquid-liquid extraction data7, 9, 10, 18. In these studies, the formations of
1:1, 1:2 and 1:3 complexes of Pa(V) with sulfate or oxalate ligands were
hypothesized and described by the following chemical equilibria:

$\ce{PaO(OH)^{2+}+iX^{2-}+H+<=>[\beta_{app,i}][PaO(X)_{i}]^{3-2i}+H2O}$ (2)

where β_app,i are the (global) apparent formation constants and $i$ varies
from 1 to 3. Again, since the formed species are defined up to a given number
of water molecules, only apparent constants have been determined so far. By
data analysis, the formation of the 1:1 to 1:3 complexes were confirmed in
both the cases and the apparent formation constants have been determined at
$25\text{\,}\mathrm{\SIUnitSymbolCelsius}$ by a fitting procedure 18, 10 (see
Table 1).

Table 1:  Experimental apparent formation constants of Pa(V) complexes determined at $25\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and deduced ligand-exchange experimental equilibrium constants. $25\text{\,}\mathrm{\SIUnitSymbolCelsius}$ |  SO4^2- 18 |  C2O4^2- 10 |  log $K_{exc,i}^{exp}$
---|---|---|---
log $\beta_{app,1}$ | $3.9(2)$ | $7.4(3)$ | $3.5(4)$
log $\beta_{app,2}$ | $7.0(2)$ | $14.5(8)$ | $7.5(9)$
log $\beta_{app,3}$ | $8.6(2)$ | $19.1(7)$ | $10.5(8)$

From the log $\beta_{app,i}$ values reported in Table 1, we observe a similar
trend in both the series, with the log $\beta_{app,2}$ values being nearly
twice the log $\beta_{app,1}$ ones, and the log $\beta_{app,3}$ ones being
nearly 2.5 times larger than the log $\beta_{app,1}$ ones. Because of this, we
hypothesize that for each $i$ value, the two products, i.e. [PaO(SO4)_i]^3-2i
and [PaO(C2O4)_i]^3-2i, have similar natures. Therefore, even if the nature of
the ending complexes is not (yet) fully known, we propose to focus on three
ligand-exchange reactions:

$\ce{[PaO(SO_{4})_{i}]^{3-2i}+iC_{2}O_{4}^{2-}<=>[K_{exc,i}][PaO(C_{2}O_{4})_{i}]^{3-2i}+iSO_{4}^{2-}}$
(3)

The determination of the experimental exchange constants is in this case
trivial as it derives from the difference of the individual apparent
complexation constants:

$\text{log}\leavevmode\nobreak\ K_{exc,i}^{exp}=\text{log}\leavevmode\nobreak\
\beta_{app,i}\leavevmode\nobreak\
[\ce{C2O4^{2-}}]-\text{log}\leavevmode\nobreak\
\beta_{app,i}\leavevmode\nobreak\ [\ce{SO4^{2-}}]$ (4)

It is important to recall that the mono-oxo bond has been observed in the
[PaO(C_2O_4)_3]^3-, i.e. for the complex of highest ligation with the ligand
that best bind to Pa(V). As a consequence, we assume that the mono-oxo bond is
maintained in all the complexes of interest. Thus, the main point to be
addresses relates to the potential departure of the hydroxyl group, present in
the starting [PaO(OH)]^2+ species, as mentioned earlier.

Our computational study aims at fairly reproducing the
$\text{log}\leavevmode\nobreak\ K_{exc,i}^{exp}$ values by using solvation and
chemical models of increasing complexities. To achieve this goal, we formulate
the following working hypotheses:

  1. 1.

In all the complexes, the ligands bind in a bidentate fashion (see Figure 1);

  2. 2.

We can compute the log of the exchange constants by means of a thermodynamic
cycle19, with computation of the log of the gas phase constants at our
reference level of theory and computation of the solvation energies of the
species with a polarizable continuum model (PCM);

  3. 3.

We assume that the errors done by our PCM approach concerning the solvation of
the free SO4^2- and C2O4^2- ligands are similar, meaning that explicit
solvation will not be applied on the free ligands (note that the opposite
would potentially turn out to be a dead end owing the complexity of the
generated microsolvated structures20, 21).

Also, we recall the previously formulated hypotheses still apply, meaning that
we assume that the mono-oxo bond is present in all the complexes and also that
for each $i$ reaction, the corresponding sulfate and oxalate complexes have a
similar nature.

Finally, we define the log of the exchange constants as follows ($T$ is fixed
at 298.15 K and $P$ at 1 atm):

$\text{log}\leavevmode\nobreak\ K_{exc,i}=-\frac{\Delta
G_{exc,i}^{0}}{2.303RT}$ (5)

Fig. 1: Schematic representation of the geometrical models for the Pa(V)
complexes that are used to compute the $\text{log}\leavevmode\nobreak\
K_{exc,i}$ values.

##  3 Computation details

Geometries and electronic structures of Pa(V) are expected to be well
described by single-reference approaches, due to the closed-shell nature of
the valence-free reference Pa^5+ ion. We have thus directly restricted our
preparatory study to the MP2 and Kohn-Sham DFT methods, with inclusion of
scalar relativistic effects (via the Douglas-Kroll-Hess Hamiltonian22, 23, 24
or a suited relativistic pseudopotential), and in the latter case the use of a
hybrid exchange-correlation functional. We have actually compared results
obtained with the MP225, DFT/B3LYP26 and DFT/PBE027, 28 methods, with
different basis sets (with an all-electron basis 29 vs. bases 30, 31, 32
combined with a relativistic pseudopotential for Pa33), and different codes
(Gaussian34, ORCA 4.2.1 35 and OpenMolcas36). The results and a more complete
discussion can be found in Supporting Information. Our main findings are:

  1. 1.

Basis sets combined with pseudopotentials generally give geometries in good
agreement with the employed all-electron basis, unless a specific and
unrecommended implementation is used.

  2. 2.

The use of the def2-TZVP and aug-cc-pVTZ basis sets on the light atoms (i.e.
all atoms besides Pa) leads to a significant improvement over the previous
results that were based on the former 6–31+G* ones 10.

  3. 3.

All the tested approaches are in good agreement with the EXAFS structure, with
perhaps a slight advantage of the DFT/PBE0 approach, essentially concerning
the mono-oxo bond distance.

In the remainder of the article, the gas-phase geometries, harmonic
frequencies, and gas-phase free energies of the oxalate and sulfate complexes
with Pa(V) were obtained using the Gaussian 16.A.03 software 34. All the
reported structures correspond to minima of the energy (no negative
frequency). The Kohn-Sham equation was solved by using the hybrid PBE0
exchange-correlation functional 27, 28. In these calculations, the def2-TZVP
(second generation of triple-$\zeta$ plus polarization quality) basis set 32
has been used for all the light elements (H, C, O, S). For the protactinium
atom, the small-core scalar relativistic effective core potential ECP60MWB was
used 33 along with its associated basis set, formed by contracting a
14s13p10d8f set into a 10s9p5d4f one 31.

The aqueous solvation free energies were computed at the Hartree-Fock (HF)
level with the same basis sets and pseudopotentials as used in the DFT/PBE0
calculations. To mimic the long-range solvent effects, the conductor-like
polarizable continuum model (CPCM) 37, 38 was used as it is implemented in the
Gaussian 16 program package. In the absence of any alternative, the UFF radius
was for the Pa atom ($1.712\text{\,}\mathrm{\text{\r{A}}}$), while the rest of
the cavity was based on UAHF radii39 for the other atoms/groups
(O_oxo=$1.59\text{\,}\mathrm{\text{\r{A}}}$,
O_ligand=$1.35\text{\,}\mathrm{\text{\r{A}}}$,
C=$1.68\text{\,}\mathrm{\text{\r{A}}}$,
S=$1.98\text{\,}\mathrm{\text{\r{A}}}$,
OH-=$1.59\text{\,}\mathrm{\text{\r{A}}}$,
H2O=$1.68\text{\,}\mathrm{\text{\r{A}}}$). We recall here that the UAHF radii
depend on the charge and environment of the atom of interest. For the
bidentate SO4^2- and C2O4^2- ligands, we have considered ad hoc equal $-$1/2
charges on the oxygen atoms, hence the O_ligand radius value. Also, note that
no sphere is created around the H atoms (united-atom approach), which explains
the OH- and H2O notations. Since those radii were initially optimized with the
alpha scaling factor of 1.2, it was manually set to this value in our
calculations.

Choice was made to apply the UAHF model instead of the UAKS after initial
calculations concerning the geometry of the [PaO(C2O4)3]^3- complex and the
determination of the $\text{log}\leavevmode\nobreak\ K_{exc,3}$ value. In
fact, a positive solvation contribution to $\text{log}\leavevmode\nobreak\
K_{exc,3}$ was obtained with the UAKS model, resulting in an unphysical
negative $\text{log}\leavevmode\nobreak\ K_{exc,3}$ value. In a previous study
related to astatine chemistry15, it was also found that the UAHF model was
more accurate than the UAKS one to compute ligand-exchange reaction constants,
the UAKS one being more accurate only if applied in a single-point fashion at
the UAHF geometries. In view of finding a compromise between accuracy and
simplicity of application (no need for extra single-point calculations), we
have thus here retained the UAHF model for the production calculations.

##  4 Results and discussion

In this section, we will introduce models of increasing complexity. First, we
will look at the gas-phase ligand-exchange reaction constants, and complement
these by an implicit solvation model. Then, we will also include explicitly
water molecules in the quantum chemical framework, meaning that we will use a
combined implicit and explicit solvation model. Finally, we will discuss the
possibility for maintaining the hydroxyl group on both the left-hand and
right-hand sides of the ligand-exchange reactions. At each stage, we will
compare the computed $\text{log}\leavevmode\nobreak\ K_{exc,i}$ values to the
$\text{log}\leavevmode\nobreak\ K_{exc,i}^{exp}$ ones.

###  4.1 Gas phase (GP) and implicit solvation (PCM)

The zeroth-order description of solvation consists in fully neglecting any
solvent effect, implying that gas-phase reaction constants are actually first
sought. In the process, the geometries of the sulfate and oxalate complexes of
interest are determined, together with the structures of the free ligands. For
all the complexes, we find that binding in bidentate fashion is by far
preferable (which justifies our previous hypothesis number #1).

As already mentioned, the EXAFS structure of the [PaO(C2O4)3]^3- complex has
been reported by Mendes et al.10. We note that attempts to determine the
structure of the analogous [PaO(SO4)3]^3- complex have been reported9, 40, but
that information concerning the binding modes is still lacking, meaning that
even if the presence of the short mono-oxo bond was confirmed, it is hard to
know how accurate the other Pa–O_ligand bond distances are. Therefore, we only
retain the EXAFS structure of the [PaO(C2O4)3]^3- complex for comparison
purposes.

As shown in Supporting Information, the DFT/PBE0 approach leads to a good
agreement with the experimental structure, despite a moderate discrepancy on
the Pa–O_oxo bond distance ($1.86\text{\,}\mathrm{\text{\r{A}}}$
vs.$1.75\text{\,}\mathrm{\text{\r{A}}}$ in the experiment). At this stage, the
computational/experimental origin of this discrepancy is not known, and it may
relate to experimental and/or computational uncertainty. Addressing this issue
would require generating new EXAFS experiments and/or new fits of existing
data and also performing an extensive electronic structure theory benchmark,
which is out of the scope of the present paper. In fact, we remark that for
all the Pa–O_ligand bond distances of interest (i.e. the ones of the bound
atoms) as well as for the Pa–C bond distances, the mean absolute deviation
between the computed and experimental values is below
$0.02\text{\,}\mathrm{\text{\r{A}}}$. We thus assume that the retained
DFT/PBE0 structure are overall accurate enough for further proceeding with the
computation of ligand-exchange reaction constants. Moreover, we remark that a
Pa–O_oxo bond distance of $1.86\text{\,}\mathrm{\text{\r{A}}}$ seems
significantly longer than the one in bare PaO^3+
($1.68\text{\,}\mathrm{\text{\r{A}}}$) and in the PaO(OH)^2+
($1.71\text{\,}\mathrm{\text{\r{A}}}$) precursor. The three studied reactions
have been renamed as R1, R2 and R3, depending on the number of ligands
involved in the reaction. The gas-phase structures of the complexes involved
in these reactions are displayed in Figure S1. Results are reported in Table 2
(see more specifically the “GP” column).

Table 2:  Computed values of the logarithms of the ligand-exchange reaction constants as functions of the employed solvation model (see text). The experimental values from Table 1 are repeated for convenience. $\text{log}\leavevmode\nobreak\ K_{exc,i}$ | GP | PCM |  PCM + $n$ H2O | Expt.
---|---|---|---|---
|  |  | FCS | FCS+1 |
R1 | 6.0 | 9.4 | - | - | 3.5
R2 | 2.6 | 19.9 | 12.0 | 8.8 | 7.5
R3 | 2.8 | 16.1 | 16.1 | 14.6 | 10.5

For all the computed species, the obtained lowest-energy structure was
separated by at least $10\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$
with any other conformer or isomer (in terms of Gibbs free energy at
$298.15\text{\,}\mathrm{K}$). Therefore, we have used only one structure per
species to compute the reaction constants. As can be seen in Table 2, the gas-
phase ligand-exchange reaction constants do not follow at all the experimental
trend, going from $\text{log}\leavevmode\nobreak\ K_{exc,1}$ = 6.0 (R1) to
$\text{log}\leavevmode\nobreak\ K_{exc,3}$ = 2.8 (R3). This may be a net
indication that solvent effects may be at play, suggesting that important
solute-solvent interactions must somehow be accounted for 38.

As a first solvation model, we have included a correction in the free energy
of each species arising from a PCM model (see the Computational details
section). In computing those solvation-induced corrections to the free
energies, it is important to check that the geometry that is optimized with
application of the PCM is of the same nature as the gas-phase one.
Furthermore, in our calculations, we have employed the UAHF model, meaning
that both the gas phase and PCM geometries are obtained at the HF level of
theory. Therefore, we not only need to check that those two match well, we
also need to check that the gas phase DFT/PBE0 and HF geometries are also in
qualitative accord. In all the cases, the GP(DFT/PBE0), GP(HF) and PCM(HF)
structures qualitatively matched, meaning that consistent corrections to the
free energies of the species of interest were obtained. The deduced values of
the logarithms of the reaction constants are listed in Table 2 (see in
particular the “PCM” column).

Again, even if now $\text{log}\leavevmode\nobreak\ K_{exc,1}$ is the weakest
of the three, the computational trend differs drastically from the
experimental one. Since we observe important changes by introducing the PCM,
we may already conclude that solvation is crucial in this case, even if its
treatment remains to be more precise. Within the framework of a static quantum
chemical approach, we have two main degrees of freedom, we may (i) explicitly
treat some of the solvent (here water) molecules or (ii) change the nature of
the chemical species that are involved in the studied equilibria, which will
be done in the next subsections. Note that the addition of explicitly treated
water molecules in the first coordination sphere of the Pa(V) ions may also
lead to better geometrical models for the complexes of interest.

###  4.2 Implicit plus explicit solvation (PCM + $n$ H2O)

The combination of implicit and explicit solvation usually improves the
quality of thermodynamic results at the price of an increased complexity in
the conformational search of the solute geometries 38. A key point lies in the
choice of the number of water molecules. At each addition of one water
molecule, one may introduce key interactions such as a coordination or a
hydrogen bond. Since we compute ligand-exchange reaction constants, it is also
important to ensure that the added interactions on the left-hand and right-
hand sides of the reaction are of similar nature16. As a consequence, the
computed $\text{log}\leavevmode\nobreak\ K_{exc,i}$ values may not
monotonously evolve with the number of explicitly treated water molecules, and
it may not be worth reporting each and every step.

At this stage, it is important to come back to our computational hypothesis
#3, according to which error cancellations allow us to avoid performing a
tricky microsolvation study of the free oxalate and sulfate ligands (many
water molecules being required to fill their first solvation spheres20, 21).
In fact, two independent PCM studies were successful in reproducing an
experimental trend that involves the free oxalate or the sulfate ligand:
Aquino et al.41 used a
$-986.6\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ solvation energy for
the free oxalate ligand to study its successive complexation with the Al^3+
ion and Lee and McKee42 recommended a
$-1028.8\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ solvation energy for
the free sulfate ligand to study the dissolution of salts. We note that the
resulting difference in solvation energies,
$-42.2\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$, is pretty close to
the $-40.0\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ that we have
obtained with the UAHF model, which supports our hypothesis. As a result,
water molecules are only added to the protactinium complexes. We thus rewrite
our generic ligand-exchange reaction as:

$\begin{split}\ce{[PaO(SO_{4})_{i}(H2O)_{n}]^{3-2i}+iC_{2}O_{4}^{2-}<=>[$K_{exc,i}$]}\\\
\ce{[PaO(C_{2}O_{4})_{i}(H2O)_{n}]^{3-2i}+iSO_{4}^{2-}}\end{split}$ (6)

where n has to be carefully chosen, as mentioned above.

In practice, a stepwise water addition was performed, with a conformational
search based on multiple starting structures, having for instance ligands
bound in monodentate or bidentate fashions, and also water molecules
interacting with either the Pa(V) ion (coordination bonds) or the ligands
(hydrogen bonds). Out of this stepwise study, we concluded that we needed to
at least saturate the first coordination sphere of the Pa(V) ion to obtain
similar structures for each pair of sulfate and oxalate complexes. Therefore,
the intermediate steps with an incomplete first coordination sphere are not
reported here.

####  4.2.1 First coordination sphere (FCS)

The structures of the Pa(V) complexes with sulfate and oxalate ligands,
displaying saturated first coordination spheres (by the coordinated ligands,
in fact, always in bidentate modes, and by the added water molecules) are
shown in Figure 2. This saturation was confirmed by calculations with the
addition of one additional water molecule, which always ended up to structures
for which one of the explicitly treated water molecules was hydrogen-bonded to
one of the anionic ligands. Therefore, we conclude that the coordination
number (CN) of the Pa(V) ion can be either 8 ($i$ = 1, 2) or 7 ($i$ = 3). For
the [PaO(C2O4)3]^3- complex, our CN notably matches the EXAFS structure10.

(a) [PaO(C2O4)(H2O)5]^+

(b) [PaO(C2O4)2(H2O)3]^-

(c) [PaO(C2O4)3]^3-

(d) [PaO(SO4)(H2O)5]^+

(e) [PaO(SO4)2(H2O)3]^-

(f) [PaO(SO4)3] ^3-

Fig. 2: Representations of the GP(DFT/PBE0) structures of the Pa(V) complexes,
obtained with just-saturated first-coordination spheres (FCS). Color code:
dark blue (Pa), yellow (S), red (O), light blue (C) and white (H).

It is clear from Figure 2 that for each reaction, both the complexes display
similar interactions. Thus, we can readily proceed with the determination of
the reaction constants (see “PCM + $n$ H2O” / "FCS” column in Table 2). For
the sake of pedagogy, we will now discuss the results in reversed order, i.e.
from R3 till R1. Since there is no added water molecule in the [PaOX3]^3-
complexes, note that the value $\text{log}\leavevmode\nobreak\ K_{exc,3}$ of
16.1 is simply reported from the (previous) “PCM” column.

Things start to be less trivial with the [PaOX2]^- complexes. In both the
sulfate and oxalate complexes, three water molecules are introduced. As
before, we have checked that each series of three structures, namely
GP(DFT/PBE0), GP(HF) and PCM(HF) qualitatively match, meaning that these
correspond to a same isomer and conformer. Also, we have excluded the
possibility for the need of considering an ensemble of conformers, since no
other conformer was found at less than
$10\text{\,}\mathrm{kJ}\text{\,}{\mathrm{mol}}^{-1}$ above the retained one in
terms of Gibbs free energy (both in the gas phase and after introducing the
PCM correction). We now see a significant improvement of the computed
$\text{log}\leavevmode\nobreak\ K_{exc,2}$ value, since it is now for the
first time below the $\text{log}\leavevmode\nobreak\ K_{exc,3}$ one. In fact,
the $\text{log}\leavevmode\nobreak\ K_{exc,3}-\text{log}\leavevmode\nobreak\
K_{exc,2}$ value is already in good agreement with the experiment (4.1 vs.
3.0). Therefore, we conclude that the improvement of our computational
description goes together with a better accord with the experiment, which is
quite encouraging for what is to follow.

In the case of R1, we have not reported any value in Table 2. In fact,
depending on the starting conformer, we have obtained
$\text{log}\leavevmode\nobreak\ K_{exc,1}$ values ranging from 10.0 to -11.6,
and we have been unable to find couples of GP(HF) and PCM(HF) structures that
match. In this particular case, the static approach combined with a PCM
correction seems to act as a random number generator and it is not recommended
to report any data. We thus need to figure out how to obtain more stable value
in the R1 case, which may come with a change in the chemical process (more
water molecules, change in chemical composition) or even of paradigm (by going
beyond the static approach for instance). The latter possibility being out of
the scope of the present paper, we will continue by proceeding with the other
two possibilities.

Despite the relative favorable correspondence with the experimental data for
R2 and R3, it is still important to check that a potential change in the
process does not destroy this mighty agreement. In particular, it is important
to check that addition of one water molecule outside the first coordination
shell (denoted FCS+1 in the remainder of the article), which may physically
lead to a better stabilization of one of the dianionic ligands, maintains this
successful trend.

####  4.2.2 Beyond the first coordination sphere (FCS+1)

With the addition of one more water molecule in each complex, an interaction
with one of the ligands is actually further introduced. As before, we find a
pathological behaviour in the R1 case, and regular behaviour in the R2 and R3
cases. As can be seen in Table 2, the computed $\text{log}\leavevmode\nobreak\
K_{exc,2}$ and $\text{log}\leavevmode\nobreak\ K_{exc,3}$ values become
practically closer to the experimental ones while the
$\text{log}\leavevmode\nobreak\ K_{exc,3}-\text{log}\leavevmode\nobreak\
K_{exc,2}$ difference remains in fair agreement with experiment. We thus
conclude that the previous encouraging results obtained for R2 and R3 were not
fortuitous, since improvement of the description maintains them. One should
point out that using FCS+2 configurations did not bring any substantial
changes to the computed $\text{log}\leavevmode\nobreak\ K_{exc,n}$ values,
further confirming the protactinium(V) coordination number and the
thermodynamic stability of considered species. Moreover, since addition of one
or two water molecules does not solve the pathological situation encountered
in the R1 case, we pursue by increasing the chemical complexity. Since we have
seen that the FCS+1 results are overall better than the FCS ones, we stick to
the FCS+1 level in this complementary study.

###  4.3 Probing the relevance of the hydroxyl-group scenario

All the previous results were based on the hypothesis that the hydroxyl group,
present in the [PaO(OH)]^2+ starting species, is withdrawn by any
complexation. While this is quite clear to us in the R3 case (see the Choice
of the systems and working hypotheses section), for which we have obtained a
good agreement with the experiment, this hypothesis may be questioned in the
R1 and R2 cases, in particular the R1 one for which we were so far unable to
report any reliable value. If we make the opposite hypothesis, i.e. that the
hydroxyl group is maintained, the global formation reactions of the complexes
of interest now write:

$\ce{[PaO(OH)]^{2+}+iX^{2-}<=>[\beta_{i}^{\prime}][PaO(OH)X_{i}]^{2-2i}}$ (7)

where we recall $i$ is either 1 or 2. These global formation constants
translate into the following ligand-exchange reactions:

$\begin{split}\ce{[PaO(OH)(SO_{4})_{i}(H2O)_{n}]^{2-2i}+iC_{2}O_{4}^{2-}<=>[$K_{exc,i}^\prime$]}\\\
\ce{[PaO(OH)(C_{2}O_{4})_{i}(H2O)_{n}]^{2-2i}+iSO_{4}^{2-}}\end{split}$ (8)

To distinguish these reactions from the previous ones (referred to as R1 and
R2), they will be denoted as R1-OH and R2-OH, respectively, in the rest of the
article.

The four new FCS+1 GP(DFT/PBE0) structures are displayed in Figure 3. The
hydroxyl group is always found to be in trans position with respect to the
mono-oxo one, and the sulfate and oxalate ligands are still bound in the
bidentate fashion, as in all the previous structures. In those four new
structures, the CN of the Pa(V) ion is 7. Luckily, this reduces the
possibilities for arranging the water molecules within the FCS and also
further away from it. In practice, we now have the expected correspondence
between the GP(DFT/PBE0), GP(HF) and PCM(HF) structures even in the complexes
that are involved in R1–OH, meaning that we may now propose properly
determined log $K_{exc,i}$ values for all reactions. Also, we now have another
log $K_{exc,2}$ value, associated with the R2-OH equilibrium. Those two new
values, together with our previous best estimates (FCS+1) are reported in
Table 3.

(a) [PaO(OH)(C2O4)(H2O)3], H2O

(b) [PaO(OH)(C2O4)2(H2O)]^2-, H2O

(c) [PaO(OH)(SO4)(H2O)3], H2O

(d) [PaO(OH)(SO4)2(H2O)]^2-, H2O

Fig. 3: Representations of the GP(DFT/PBE0) structures of the complexes involved in the R1-OH and R2-OH reactions, obtained with one more water than just saturated first-coordination spheres (FCS+1). Color code: dark blue (Pa), yellow (S), red (O), light blue (C) and white (H). Table 3:  Computed values of the logarithms of the ligand-exchange reaction constants as functions of the chemical reaction (R/R-OH). Calculations are performed with one more water than just-saturated first-coordination spheres (FCS+1). Values are repeated from Tables 1 and 2 for convenience. log $K_{exc,i}/K_{exc,i}^{\prime}$ | Comput. | Expt.
---|---|---
R1/R1-OH | \- / 7.0 | 3.5
R2/R2-OH | 8.8 / 12.1 | 7.5
R3/R3-OH | 14.6 / - | 10.5

As can be seen in this Table, the values obtained for R1-OH and R2-OH are both
quite consistent with the experiment. Moreover, the
$\text{log}\leavevmode\nobreak\ K_{exc,2}-\text{log}\leavevmode\nobreak\
K_{exc,1}$ difference is also in good agreement with the experiment. Even if
this agreement, based on the analysis of apparent constants for the
experimental values, may not be a definitive proof for the occurrence of the
hydroxyl-group scenario in the R1 case (more experimental information may be
required), we conclude by predicting that this scenario must be at play since
our static approach is found to fail at reproducing the trends otherwise.

Concerning the R2 reaction, both scenarios lead to computed
$\text{log}\leavevmode\nobreak\ K_{exc,2}$ values in decent accord with the
experiment. Therefore, it is hard to know a priori if the hydroxyl group is
maintained or withdrawn in the 1:2 complexation. Moreover, being in fact an
intermediate situation sandwiched between the two clear R1 and R3 cases, it
could happen that the hydroxyl group is maintained in the sulfate complex and
not in the oxalate one. Unfortunately, this cannot be well probed with our
calculations since we need to have similar pictures in both the left-hand and
right-hand sides of the reactions to favor error cancellations, especially on
the solvation contributions.

##  5 Concluding remarks

In this work, we performed a step-by-step study to deduce ligand-exchange
reaction constants. The studied equilibria involved Pa(V) complexes with
sulfate and oxalate ligands, displaying the characteristic Pa mono-oxo bond,
and in some cases hydroxyl groups. Our computational approach is based on the
DFT/PBE0 level to determine the gas-phase complexation constants, and on the
UAHF model for the solvation corrections. In the absence of an alternative,
note that the UFF radius has been used for protactinium.

We have seen that solvation of Pa(V) complexes may be tricky, in the sense
that the implicit model alone does not suffice to reach good values, even if
relative constants are looked for. Instead, we show that we need to at least
saturate the coordination sphere of the Pa(V) ion by water addition. Also, in
the case of low-coordinated complexes, we reveal that the speciation may be
more complex than expected and that the relevant complexes may still display
one hydroxyl group.

The determination of the speciation of Pa(V) in non-complexing and in
complexing media is still in its infancy. It is extremely hard to firmly
identify new chemical species, and we believe that our work, which has led to
the validation of a computational protocol, can open the way for future
studies assisted by computations. In particular, we have shown that
calculations of ligand-exchange reaction constants can be quite informative
concerning trends (which ligands is a stronger complexant) and also regarding
the occurrence of given chemical groups and/or bonds in the formed complexes.

## Author Contributions

H.O. and R.M. defined the working hypotheses and the computational strategy.
H.O. was responsible for performing and/or checking all the calculations. J.D.
performed preliminary calculations concerning the hydroxyl-group scenario. The
first draft of the manuscript was written by H.O. and R.M. All the authors
contributed to the scientific development of the work and were given the
opportunity to amend the manuscript.

## Conflicts of interest

The authors declare no conflict of interest.

## Acknowledgements

This work was supported by the Région des Pays de la Loire (RCT-PaPAn project)
and by the French National Research Agency (ANR-21-CE29-0027, LABEX
CaPPA/ANR-11-LABX-0005-01 and I-SITE ULNE/ANR-16-IDEX-0004 ULNE). HPC
resources from CCIPL (Centre de Calcul Intensif des Pays de la Loire) were
used.

## Notes and references

  * Wilson 2012 R. Wilson, _Nat. Chem._ , 2012, 4, 586.
  * Viatkina _et al._ 2017 A. V. Viatkina, M. G. Kozlov and V. V. Flambaum, _Phys. Rev. A_ , 2017, 95, 022503.
  * Naubereit _et al._ 2018 P. Naubereit, T. Gottwald, D. Studer and K. Wendt, _Phys. Rev. A_ , 2018, 98, 022505.
  * Naubereit _et al._ 2018 P. Naubereit, D. Studer, A. V. Viatkina, A. Buchleitner, B. Dietz, V. V. Flambaum and K. Wendt, _Phys. Rev. A_ , 2018, 98, 022506.
  * Marquardt _et al._ 2004 C. M. Marquardt, P. J. Panak, C. Apostolidis, A. Morgenstern, C. Walther, R. Klenze and T. Fanghänel, _Radiochim. Acta_ , 2004, 92, 445–447.
  * Le Naour _et al._ 2022 C. Le Naour, M. Maloubier and J. Aupiais, _Radiochim. Acta_ , 2022, 110, 481–493.
  * Le Naour _et al._ 2019 C. Le Naour, J. Roques, C. Den Auwer, P. Moisy and J. Aupiais, _Radiochim. Acta_ , 2019, 107, 979–991.
  * Mendes _et al._ 2013 M. Mendes, S. Leguay, C. Le Naour, S. Hamadi, J. Roques, P. Moisy, D. Guillaumont, S. Topin, J. Aupiais, C. Den Auwer and C. Hennig, _Inorg. Chem._ , 2013, 52, 7497–7507.
  * Le Naour _et al._ 2005 C. Le Naour, D. Trubert, M. V. Di Giandomenico, C. Fillaux, C. Den Auwer, P. Moisy and C. Hennig, _Inorg. Chem._ , 2005, 44, 9542–9546.
  * Mendes _et al._ 2010 M. Mendes, S. Hamadi, C. Le Naour, J. Roques, A. Jeanson, C. Den Auwer, P. Moisy, S. Topin, J. Aupiais, C. Hennig and M. V. Di Giandomenico, _Inorg. Chem._ , 2010, 49, 9962–9971.
  * De Sio and Wilson 2014 S. M. De Sio and R. E. Wilson, _Inorg. Chem._ , 2014, 53, 12643–12649.
  * Wilson _et al._ 2016 R. E. Wilson, S. De Sio and V. Vallet, _Eur. J. Inorg. Chem._ , 2016, 2016, 5467–5476.
  * Champion _et al._ 2011 J. Champion, M. Seydou, A. Sabatié-Gogova, E. Renault, G. Montavon and N. Galland, _Phys. Chem. Chem. Phys._ , 2011, 13, 14984–14992.
  * Champion _et al._ 2013 J. Champion, A. Sabatié-Gogova, F. Bassal, T. Ayed, C. Alliot, N. Galland and G. Montavon, _J. Phys. Chem. A_ , 2013, 117, 1983–1990.
  * Sergentu _et al._ 2016 D.-C. Sergentu, G. David, G. Montavon, R. Maurice and N. Galland, _J. Comput. Chem._ , 2016, 37, 1345–1354.
  * Sergentu _et al._ 2016 D.-C. Sergentu, D. Teze, A. Sabatié-Gogova, C. Alliot, N. Guo, F. Bassal, I. Da Silva, D. Deniaud, R. Maurice, J. Champion, N. Galland and G. Montavon, _Chem. Eur. J._ , 2016, 22, 2964–2971.
  * Guo _et al._ 2016 N. Guo, D.-C. Sergentu, D. Teze, J. Champion, G. Montavon, N. Galland and R. Maurice, _Angew. Chem. Int. Ed._ , 2016, 128, 15595–15598.
  * Di Giandomenico and Le Naour 2009 M. Di Giandomenico and C. Le Naour, _Inorg. Chim. Acta_ , 2009, 362, 3253–3258.
  * Ho 2015 J. Ho, _Phys. Chem. Chem. Phys_ , 2015, 17, 2859–2868.
  * Gao and Liu 2004 B. Gao and Z.-f. Liu, _J. Chem. Phys._ , 2004, 121, 8299–8306.
  * Gao and Liu 2005 B. Gao and Z.-f. Liu, _J. Phys. Chem. A_ , 2005, 109, 9104–9111.
  * Douglas and Kroll 1974 M. Douglas and N. M. Kroll, _Ann. Phys._ , 1974, 82, 89–155.
  * Hess 1986 B. A. Hess, _Phys. Rev. A_ , 1986, 33, 3742–3748.
  * Jansen and Hess 1989 G. Jansen and B. A. Hess, _Phys. Rev. A_ , 1989, 39, 6016–6017.
  * Møller and Plesset 1934 C. Møller and M. S. Plesset, _Phys. Rev._ , 1934, 46, 618–622.
  * Becke 1998 A. D. Becke, _J. Chem. Phys._ , 1998, 109, 2092–2098.
  * Perdew _et al._ 1996 J. P. Perdew, M. Ernzerhof and K. Burke, _J. Chem. Phys._ , 1996, 105, 9982–9985.
  * Adamo and Barone 1999 C. Adamo and V. Barone, _J. Chem. Phys._ , 1999, 110, 6158–6170.
  * Pantazis and Neese 2011 D. A. Pantazis and F. Neese, _J. Chem. Theory Comput._ , 2011, 7, 677–684.
  * Kendall _et al._ 1992 R. A. Kendall, T. H. Dunning and R. J. Harrison, _J. Chem. Phys._ , 1992, 96, 6796–6806.
  * Cao and Dolg 2004 X. Cao and M. Dolg, _J. Mol. Struct. THEOCHEM_ , 2004, 673, 203–209.
  * Weigend and Ahlrichs 2005 F. Weigend and R. Ahlrichs, _Phys. Chem. Chem. Phys._ , 2005, 7, 3297–3305.
  * Cao _et al._ 2003 X. Cao, M. Dolg and H. Stoll, _J. Chem. Phys._ , 2003, 118, 487–496.
  * Frisch _et al._ 2016 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, _Gaussian 16 Revision A.03_ , 2016, Gaussian Inc. Wallingford CT.
  * Neese 2017 F. Neese, _WIREs Comput. Mol. Sci._ , 2017, 8, e1327.
  * Aquilante _et al._ 2020 F. Aquilante, J. Autschbach, A. Baiardi, S. Battaglia, V. A. Borin, L. F. Chibotaru, I. Conti, L. De Vico, M. Delcey, I. Fdez. Galván, N. Ferré, L. Freitag, M. Garavelli, X. Gong, S. Knecht, E. D. Larsson, R. Lindh, M. Lundberg, P. Å. Malmqvist, A. Nenov, J. Norell, M. Odelius, M. Olivucci, T. B. Pedersen, L. Pedraza-González, Q. M. Phung, K. Pierloot, M. Reiher, I. Schapiro, J. Segarra-Martí, F. Segatta, L. Seijo, S. Sen, D.-C. Sergentu, C. J. Stein, L. Ungur, M. Vacher, A. Valentini and V. Veryazov, _J. Chem. Phys._ , 2020, 152, 214117\.
  * Barone and Cossi 1998 V. Barone and M. Cossi, _J. Phys. Chem. A_ , 1998, 102, 1995–2001.
  * Cossi _et al._ 2003 M. Cossi, N. Rega, G. Scalmani and V. Barone, _J. Comput. Chem._ , 2003, 24, 669–681.
  * Barone _et al._ 1997 V. Barone, M. Cossi and J. Tomasi, _J. Chem. Phys._ , 1997, 107, 3210–3221.
  * Di Giandomenico _et al._ 2008 M. V. Di Giandomenico, M. Mendes, C. Le Naour, C. DenAuwer, P. Moisy and C. Hennig, _Behaviour of Protactinium (V) in complexing media_ , Conference: Atalante 2008: Nuclear fuel cycle for a sustainable future, Montpellier (France), 19-23 May 2008, 2008.
  * Aquino _et al._ 2000 A. J. A. Aquino, D. Tunega, G. Haberhauer, M. Gerzabek and H. Lischka, _Phys. Chem. Chem. Phys._ , 2000, 2, 2845–2850.
  * Lee and McKee 2011 T. B. Lee and M. L. McKee, _Inorg. Chem._ , 2011, 50, 11412–11422.

See pages 1-8 of ESI.pdf


# Using geometric algebra to represent curvature in shell theory with
applications to Starling resistors

Alastair L Gregory  Department of Engineering, University of Cambridge,
Trumpington Street, Cambridge, CB2 1PZ<EMAIL_ADDRESS>Anurag Agarwal
Department of Engineering, University of Cambridge, Trumpington Street,
Cambridge, CB2 1PZ  Joan Lasenby  Department of Engineering, University of
Cambridge, Trumpington Street, Cambridge, CB2 1PZ

(31 October 2017)

###### Abstract

We present a novel application of rotors in geometric algebra to represent the
change of curvature tensor, that is used in shell theory as part of the
constitutive law. We introduce a new decomposition of the change of curvature
tensor, which has explicit terms for changes of curvature due to initial
curvature combined with strain, and changes in rotation over the surface. We
use this decomposition to perform a scaling analysis of the relative
importance of bending and stretching in flexible tubes undergoing self excited
oscillations. These oscillations have relevance to the lung, in which it is
believed that they are responsible for wheezing. The new analysis is
necessitated by the fact that the working fluid is air, compared to water in
most previous work. We use stereographic imaging to empirically measure the
relative importance of bending and stretching energy in observed self excited
oscillations. This enables us to validate our scaling analysis. We show that
bending energy is dominated by stretching energy, and the scaling analysis
makes clear that this will remain true for tubes in the airways of the lung.

##  1 Introduction

Self excited oscillations of flexible tubes driven by fluid flow have been a
subject of interest for some time, and there is a considerable literature on
the subject, which is reviewed by [1, 2, 3, 4, 5, 6]. Experimental rigs
designed to study this phenomenon are often called Starling resistors. We are
interested in this phenomenon because of its possible relevance to wheezing in
the lung [7], which is one of the most commonly heard lung sounds used for
diagnosis [8, 9]. Previous work on Starling resistors has largely used water
as the working fluid. In the lung the working fluid is air. This means that
the density ratio between the working fluid and the tube material in the lung
is significantly different from almost all of previously completed work on
Starling resistors. Previous modelling work usually neglects wall inertia,
using instead a “tube law” [10], but there is strong evidence that wall
inertia is significant in the lung from [8, 11], who showed that when the
density of the fluid breathed in is changed, the frequency of the wheezes is
not affected significantly. It is clear therefore that the change in density
ratio results in a qualitatively different mechanism. For this reason we have
been conducting our own experiments, and creating models to understand the
onset of oscillations.

The flexible tube itself is generally modelled as an elastic shell.
Traditional shell theories [12, 13, 14, 15, 16, 17] are well developed but
difficult to implement. We have found that linearised shell theories do not
provide good predictions of the frequencies of oscillation, and we believe
that this is due in part to the fact that we have observed that oscillations
start from a collapsed or partially collapsed state. To use current
geometrically non-linear shell theories would require a numerical simulation
of a very complex fluid structure interaction problem, which would be of
similar value to experimental results (though arguably harder to implement),
and would provide the same problem of being difficult to physically interpret
due to the complex nature of the oscillations observed. Instead we would like
to gain a physical understanding of the important mechanisms behind the
oscillations, and this is difficult with shell theories based in differential
geometry [14, 15], in particular due to the lack of physical interpretations
of the change of curvature tensor in general situations. We recently
introduced geometric algebra to shell theory [18], which allowed us to express
the fundamental laws in a component-free form and clarify the role of angular
velocity and moments through the use of bivector representation. For an
introduction to the basics of geometric algebra see [19]. One of the most
powerful aspects of geometric algebra lies in the use of rotors to represent
rotations. In [20] these have been used to simplify Simo and Vu Quoc’s
numerical algorithm [21] for modelling the nonlinear behaviour of rods. In
projective and conformal geometry [19, §10] rotors have allowed geometric
primitives to be represented in a more simple and lucid manner, and in
relativity [22] and relativistic analogies [23] rotors can simplify
transformation between frames of reference. In this paper we make use of
rotors to better understand the change of curvature of a shell, which is of
prime importance to the constitutive law of the shell, but whose
representation has long caused controversy. In [12] at least 10 different
linearised shell theories are presented, and the differences are primarily
caused by disagreements over how to represent changes of curvature. [14] has
provided a tensor definition of the change of curvature that has become
accepted, however, the utility of this expression is limited by its
complexity. We have been able to simplify the representation of this tensor
using rotors, allowing a more lucid and physical interpretation of changes of
curvature. We take advantage of this to allow us to understand the importance
of the change of curvature in the context of our Starling resistor
experiments.

In order to compare results from the shell theory to our experimental results
we need to be able to calculate the kinematic parameters associated with the
deformation of the flexible tube. To enable this we use stereoscopic imaging,
which to our knowledge is the first time it has been used in the study of
Starling resistors. We take high speed video of the tube at the onset of
oscillation, and are able to track the motion of the surface, and consequently
compare the predictions of shell theory with empirical calculation.

##  2 Understanding Changes in Curvature

There is an energy associated with any deformation of a shell, and Koiter [17]
proposes the following form for this energy,

$\rho_{0}U=\frac{Eh}{2(1-\nu^{2})}\mathinner{\left((1-\nu)\operatorname{tr}(\mathsf{E}^{2})+\nu\operatorname{tr}(\mathsf{E})^{2}\right)}+\frac{Eh^{3}}{24(1-\nu^{2})}\mathinner{\left((1-\nu)\operatorname{tr}(\mathsf{H}^{2})+\nu\operatorname{tr}(\mathsf{H})^{2}\right)}.$
(1)

$U$ is the internal energy per unit mass of the shell, defined on the
reference configuration, $\rho_{0}$ is the time independent area density of
the shell on the reference configuration, $E$ is Young’s modulus, $\nu$ is
Poisson’s ratio, $h$ is the shell thickness (which in Koiter’s theory is
assumed constant), $\mathsf{E}$ is the two-dimensional Green-Lagrange strain
tensor defined on the reference configuration, $\mathsf{H}$ is the change of
curvature tensor, and $\operatorname{tr}$ is the trace operator. From (1) we
can derive the governing equations of the shell (for more details see [18]).
The first term on the right hand side of (1) represents the stretching energy,
and the second term represents the bending energy.

In general a shell is a body in which the thickness is smaller than the other
relevant defining length scales. The use of (1) for the energy of deformation
implies that we additionally assume that the midsurface of the body remains
the midsurface under deformation, a material line that is normal to the
midsurface remains normal to it under deformation, the shell thickness remains
constant with time, the first and second moments of density relative to the
midsurface are zero, strains within the shell are small, and so is the normal
stress (see [18] for further discussion).

Following the notation of [18] we take $B$ and $S$ to be the reference and
spatial configurations of the shell, and $X\in B$ and $x\in S$ to be locations
on these configurations. $\phi_{t}$ is the motion of the shell, meaning that
at time $t$, the point $X\in B$ is at the position $\phi_{t}(X)\in S$.
$\mathsf{G}$ and $\mathsf{g}$ are the identity functions on the reference and
spatial configurations. $Y$ and $y$ are vectors within the tangent spaces of
$B$ and $S$ respectively. $\\{X^{i}\\},i=1,2$ is a coordinate system on the
reference configuration $B$, which we can then use to define the frame on the
reference configuration $\\{E_{i}=\partial X/\partial X^{i}\\}$. $\\{E^{i}\\}$
is the reciprocal frame that satisfies $E^{i}\cdot E_{j}=\delta^{i}_{j}$.
$\\{x^{i}\\}$, $\\{e_{i}\\}$ and $\\{e^{i}\\}$ are the similarly defined
coordinate system, frame, and reciprocal frame on the spatial configuration.
The shell undergoing deformation is embedded within a flat three-dimensional
Euclidean space $\mathbb{E}^{3}$.

If $A$ and $B$ are general multivectors, then $AB$ is the geometric product
between them, $A\cdot B$ is the inner product, $A\wedge B$ is the outer
product, and $A\times B$ is the commutator product, defined by $A\times
B=\frac{1}{2}(AB-BA)$ (see [19, §4.1.3]). We also take
$\boldsymbol{\bm{\times}}$ (compared to $\times$) to be the cross product
between two vectors. If $I$ is the pseudoscalar of a three-dimensional space,
and $a$ and $b$ are vectors, then $a\boldsymbol{\bm{\times}}b=-Ia\wedge b$.

We are particularly interested in the change of curvature that is encoded in
$\mathsf{H}$. To understand this we must understand the curvature tensors on
the reference and spatial configurations $\mathsf{B}$ and $\mathsf{b}$. If
$E_{3}$ and $e_{3}$ are the normal vectors to the reference and spatial
configurations respectively, then $\mathsf{B}(Y)$ and $\mathsf{b}(y)$ are
given by,

$\mathsf{B}(Y)=-Y\cdot\partial E_{3},\quad\mathsf{b}(y)=-y\cdot\partial
e_{3},$ (2)

where $\partial$ is the intrinsic vector derivative to any surface. The
relationship between $\partial$ and the vector derivative of $\mathbb{E}^{3}$
is explained in [18]. On $B$ we can expand $\partial$ as
$\partial=E^{i}\partial/\partial X^{i}$ and on $S$ we can expand it as
$\partial=e^{i}\partial/\partial x^{i}$ [19, §6.5.1]. $\mathsf{B}$ and
$\mathsf{b}$ give non-zero results if the surface is not flat. The change of
curvature tensor $\mathsf{H}$ is given by,

$\mathsf{H}(Y)=\bar{\mathsf{F}}\mathsf{b}\mathsf{F}(Y)-\mathsf{B}(Y),$ (3)

where $\mathsf{F}$ is the deformation gradient, defined by
$\mathsf{F}(Y)=Y\cdot\partial\phi_{t}(X)$. $\mathsf{F}$ maps from the tangent
space of $B$ to the tangent space of $S$, providing information about the
local deformation of the surface. $\bar{\mathsf{F}}$ is the adjoint of
$\mathsf{F}$, i.e. $\mathsf{F}(Y)\cdot y=Y\cdot\bar{\mathsf{F}}(y)$. The
strain tensor $\mathsf{E}$, used in the constitutive law (1), is given by,

$\mathsf{E}(Y)=\frac{1}{2}\mathinner{\left(\bar{\mathsf{F}}\mathsf{F}(Y)-Y\right)}.$
(4)

This much is well known, though in other treatments coordinate dependent
definitions of $\mathsf{H}$ are used (e.g. in [14]). To make further progress
we will now use rotors to better understand what will produce changes in
$\mathsf{H}$.

To begin, we note that we can perform a polar decomposition on $\mathsf{F}$
such that $\mathsf{F}(Y)=\mathsf{R}\mathsf{U}(Y)$ where
$\bar{\mathsf{R}}=\mathsf{R}^{-1}$, $\det\mathsf{R}=1$, and
$\bar{\mathsf{U}}=\mathsf{U}$. $\mathsf{R}$ encodes rotation, and $\mathsf{U}$
encodes stretching. We can choose to consider a frame $\\{E_{i}\\}$ on the
reference configuration that is locally orthonormal111The coordinate system
$\\{X^{i}\\}$ must be chosen such that this is the case, and it may be
necessary to use several overlapping coordinate systems to achieve this.. In
this case $E_{3}=E_{1}\boldsymbol{\bm{\times}}E_{2}=-I_{3}(E_{1}\wedge
E_{2})=-I_{3}(E_{1}E_{2})$, where $I_{3}$ is the pseudoscalar of three-
dimensional Euclidean space $\mathbb{E}^{3}$.

To find $e_{3}$ we need two unit vectors in the tangent space of the spatial
configuration that are oriented in the same way as the pair
$\mathsf{F}(E_{1}),\mathsf{F}(E_{2})$, and are orthonormal. This pair of unit
vectors is given by $\mathsf{R}(E_{1}),\mathsf{R}(E_{2})$, and so $e_{3}$ is
given by
$e_{3}=\mathsf{R}(E_{1})\boldsymbol{\bm{\times}}\mathsf{R}(E_{2})=-I_{3}(\mathsf{R}(E_{1})\wedge\mathsf{R}(E_{2}))=-I_{3}(\mathsf{R}(E_{1})\mathsf{R}(E_{2}))$.
The fact that the function $\mathsf{R}$ is a rotation means that it has an
associated rotor $R$ such that $\mathsf{R}(Y)=RY\\!\tilde{R}$, where $R$ is an
even multivector that satisfies $\tilde{R}R=R\tilde{R}=1$, where $\tilde{R}$
is the reverse of $R$. This allows us to write $e_{3}$ as,

$e_{3}=-I_{3}\mathinner{\left((RE_{1}\tilde{R})(RE_{2}\tilde{R})\right)}=-I_{3}\mathinner{\left(R(E_{1}E_{2})\tilde{R}\right)}=RE_{3}\tilde{R},$
(5)

where we have used the fact that any rotor will commute with $I_{3}$. Thus we
have shown that the rotation associated with the deformation is also the
rotation between the normal vectors $E_{3}$ and $e_{3}$, which makes intuitive
sense. We have also extended the range and domain of $\mathsf{R}$ to
$\mathbb{E}^{3}$, while the range and domain of $\mathsf{U}$ is still
constrained to the tangent space of the reference configuration, and the range
and domain of $\mathsf{F}$ are constrained to the tangent spaces of the
spatial and reference configuration respectively.

Two results that we will find useful are,

$\displaystyle Y\cdot\partial\tilde{R}=-\tilde{R}(Y\cdot\partial R)\tilde{R},$
(6a) $\displaystyle\mathsf{F}(Y)\cdot\partial e_{3}=Y\cdot\partial e_{3}.$
(6b)

(6a) follows from $\tilde{R}R=1$. (6b) has implicit assumptions that require
explanation. The expression on the left of (6b) tells us how $e_{3}$ varies
over the spatial configuration in the direction defined by $\mathsf{F}(Y)$,
which lies in the tangent space of the spatial configuration. On the right of
(6b), $e_{3}=e_{3}(x)$ has been mapped to a vector field on the reference
configuration such that $e_{3}(X)=e_{3}(\phi_{t}(X))\;\forall X\in B$. This
allows the expression on the right to tell us how $e_{3}$ varies over the
reference configuration in the direction defined by $Y$, which is tangent to
the reference configuration. The equality of these expressions is a standard
result when mapping derivatives between manifolds, which can be proven by
considering derivatives with respect to convected coordinates that satisfy
$x^{i}(x)=X^{i}(\phi_{t}^{-1}(x))$. From this point we will assume that
$\\{x^{i}\\}$ are convected coordinates.

Using these results we can write $\bar{\mathsf{F}}\mathsf{b}\mathsf{F}(Y)$ as
$\bar{\mathsf{F}}\mathsf{b}\mathsf{F}(Y)=-\bar{\mathsf{F}}(Y\cdot\partial
e_{3})$, and the argument of $\bar{\mathsf{F}}$ can be expressed as,

$\displaystyle Y\cdot\partial e_{3}=Y\cdot\partial(RE_{3}\tilde{R})$
$\displaystyle=(Y\cdot\partial R)E_{3}\tilde{R}+R(Y\cdot\partial
E_{3})\tilde{R}+RE_{3}(-\tilde{R}(Y\cdot\partial R)\tilde{R})$ (7)
$\displaystyle=R(Y\cdot\partial E_{3})\tilde{R}+[(Y\cdot\partial
R)\tilde{R}]RE_{3}\tilde{R}-RE_{3}\tilde{R}[(Y\cdot\partial R)\tilde{R}]$
$\displaystyle=R(Y\cdot\partial E_{3})\tilde{R}+[(Y\cdot\partial
R)\tilde{R}]e_{3}-e_{3}[(Y\cdot\partial R)\tilde{R}]$
$\displaystyle=R(Y\cdot\partial E_{3})\tilde{R}+[2(Y\cdot\partial
R)\tilde{R}]\times e_{3}$ $\displaystyle=\mathsf{R}(Y\cdot\partial
E_{3})+[2(Y\cdot\partial R)\tilde{R}]\times e_{3},$

where $\times$ is the commutator product. Hence we can express
$\bar{\mathsf{F}}\mathsf{b}\mathsf{F}(Y)$ as,

$\displaystyle\bar{\mathsf{F}}\mathsf{b}\mathsf{F}(Y)$
$\displaystyle=-\bar{\mathsf{F}}\mathsf{R}(Y\cdot\partial
E_{3})-\bar{\mathsf{F}}\mathinner{\left([2(Y\cdot\partial R)\tilde{R}]\times
e_{3}\right)}$ (8) $\displaystyle=-\mathsf{U}(Y\cdot\partial
E_{3})-\bar{\mathsf{F}}\mathinner{\left([2(Y\cdot\partial R)\tilde{R}]\times
e_{3}\right)}$
$\displaystyle=\mathsf{U}\mathsf{B}(Y)+\bar{\mathsf{F}}\mathinner{\left(e_{3}\times[2(Y\cdot\partial
R)\tilde{R}]\right)},$

and finally we obtain an expression for $\mathsf{H}$,

$\displaystyle\mathsf{H}(Y)$
$\displaystyle=(\mathsf{U}-\mathsf{G})\mathsf{B}(Y)+\bar{\mathsf{F}}\mathinner{\left(e_{3}\times[2(Y\cdot\partial
R)\tilde{R}]\right)}$ (9)
$\displaystyle=(\mathsf{U}-\mathsf{G})\mathsf{B}(Y)+\bar{\mathsf{F}}\mathinner{\left([RE_{3}\tilde{R}]\times[2(Y\cdot\partial
R)\tilde{R}]\right)}.$

This shows that there are two contributions to $\mathsf{H}$. Firstly, if the
reference configuration is at all curved (i.e. $\mathsf{B}(Y)$ is non-zero),
then the strain of the shell, encoded in $\mathsf{U}-\mathsf{G}$, will result
in a change of curvature. The second contribution is due to variation of the
rotor $R$ over the shell. These two kinds of change of curvature are
illustrated well by an inflating sphere and deformation of a flat plate. As a
sphere is inflated to become a larger sphere, the normal vector is unchanged,
i.e. $e_{3}=E_{3}$, hence $R=1$ everywhere. This means that the second term in
our expression for $\mathsf{H}$ will be zero. However, the surface of the
sphere will stretch, meaning that $\mathsf{U}-\mathsf{G}$ will be non-zero. In
addition, $\mathsf{B}$ will be non-zero for a sphere, which tells us that the
first term in our expression for $\mathsf{H}$ will be non-zero. By contrast,
for a flat plate $\mathsf{B}(Y)$ will be zero, meaning that $R$ must vary over
the plate in order for there to be any change of curvature.

A variant on this expression for $\mathsf{H}$ can be obtained if we express
the rotor $R$ as $R=\exp(-A/2)=\exp(-\hat{A}\theta/2)$, where $A$ is a
bivector aligned with the plane of rotation whose magnitude is equal to the
angle of rotation $\theta$. $\hat{A}$ is a unit bivector ($\hat{A}^{2}=-1$).
Note that the direction of rotation is defined by the sign of $\theta$ and the
orientation of $\hat{A}$ together. We can set the convention that
$\theta\geqslant 0$, in which case $\theta$ and $\hat{A}$ are uniquely defined
if $A$ is known. Given this definition we can express $Y\cdot\partial R$ as,

$Y\cdot\partial R=-\frac{Y\cdot\partial A}{2}\exp(-A/2)=-\frac{Y\cdot\partial
A}{2}R.$ (10)

Using this we can express $\mathsf{H}$ as,

$\mathsf{H}(Y)=(\mathsf{U}-\mathsf{G})\mathsf{B}(Y)-\bar{\mathsf{F}}\mathinner{\left((RE_{3}\tilde{R})\times(Y\cdot\partial
A)\right)}.$ (11)

The term that $\bar{\mathsf{F}}$ operates on is the commutator product of a
vector and bivector, so we can replace the commutator product with a dot
product,

$\mathsf{H}(Y)=(\mathsf{U}-\mathsf{G})\mathsf{B}(Y)-\bar{\mathsf{F}}\mathinner{\left((RE_{3}\tilde{R})\cdot(Y\cdot\partial
A)\right)}.$ (12)

Taking the inner product of a bivector with $RE_{3}\tilde{R}=e_{3}$ means that
only vectors tangential to the spatial configuration are retained, which then
means that $\bar{\mathsf{F}}$ can operate and return vectors tangential to the
reference configuration. Hence our description confirms that the range and
domain of $\mathsf{H}$ are both the tangent space of the reference
configuration.

The explicit expression for the two possible contributions to change of
curvature, shown in (12), gives a new decomposition of the change of curvature
tensor which will be of use when we try to understand the importance of
bending in the Starling resistor.

##  3 Experiment Description

###  3.1 Experimental Setup

Figure 1 shows a schematic of the experimental setup used to investigate the
oscillations of flexible tubes. Air flows into the system through (1), then
through a rotameter (2) used to monitor flowrate. The noise that the rotameter
introduces into the flow, and any other noise, is isolated from the flexible
tube by the upstream settling chamber (3). Air flows into the upstream clean
flow tube (5) section via a shaped inlet (4) that reduces separation. A
contraction (6) leads to the flexible tube (7), before an expansion (6’) leads
to the downstream clean flow tube (5’) that exits into the downstream settling
chamber (3’). Suction is provided by a fan (8). The downstream settling
chamber (3’) isolates the flexible tube from the noise from this fan.
Experiments were performed in the Acoustics Laboratory in the Department of
Engineering at the University of Cambridge.

In the experiments relevant to this paper the suction at (8) is gradually
increased until the flexible tube just starts oscillating. With the tube
oscillating in this way high speed video is recorded from 2 FASTCAM-ultima APX
cameras (produced by Photron https://photron.com) with a frame rate of
$12\,500\text{\,}\mathrm{fps}$ and a resolution of $512\times 256$ pixels
(grayscale). Our experiments require us to focus on a small flexible tube at
reasonably close range. The Photron camera has an adaptor for Nikon lens’s. We
use a $50\text{\,}\mathrm{mm}$ lens combined with a $7\text{\,}\mathrm{mm}$
extension tube to allow us to focus on the tube and have it fill most of the
frame. An aperture of f/2.8 is used.

Figure 1: Schematic of Starling resistor experiment. 1: Flow inlet, 2:
Rotameter, 3/3’: Settling chambers, 4: Clean flow inlet, 5/5’: Clean flow
tubes, 6/6’: Contraction and expansion, 7: flexible tube, 8: tube to suction
fan. The downstream settling chamber is approximately
$4\text{\,}{\mathrm{m}}^{3}$ while the upstream settling chamber is
$0.03\text{\,}{\mathrm{m}}^{3}$.

The flexible tubes used are made out of rubber latex for which
$E=$1\text{\,}\mathrm{MPa}$$ and $\nu=$0.5$$. The tube diameter is
$6\text{\,}\mathrm{mm}$, the wall thickness is $0.3\text{\,}\mathrm{mm}$ and
the unstrained length is $19\text{\,}\mathrm{mm}$. The tubes are held in an
axially strained state, so the length of the tubes during the experiment is
$25\text{\,}\mathrm{mm}$.

###  3.2 Image Processing

The high speed cameras record at $12\,500\text{\,}\mathrm{fps}$, and are
triggered together, so that every $\Delta
t=$80\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ two images of the flexible tube
are taken. A schematic of the two cameras and the flexible tube is shown in
Figure 2. Dots are drawn on the flexible tube (shown in white in Figure 2),
which indicate a set of material points we would like to track over time in
three dimensions.

It is possible to find the characteristics222The characteristics of an
imaginary pinhole camera that we replace the real camera with to allow us to
use methods from projective geometry to perform the triangulation. of two
cameras such that if a point appears in simultaneous images from both cameras,
the point’s position in three-dimensional space can be triangulated [24, 25].
Calibration involves taking at least $3$, and in general between $10$ and
$20$, simultaneous images of a chequerboard pattern in various orientations.
From this the position of the two cameras relative to each other, and their
internal parameters, can be calculated. In this method cameras are modelled as
pinhole cameras, meaning that the focal length, pixel size and skew are the
important internal parameters. In addition it is possible to account for
radial distortion of the image by the camera lens, and tangential distortion,
which occurs when the image sensor is not perfectly perpendicular to the line
of sight of the camera. This calibration is performed using the computer
vision system toolbox of Matlab® [26]. The images produced by these pinhole
camera models is what is illustrated in Figure 2.

Figure 2: Schematic of the high speed camera setup.

To find the three-dimensional tracks of the material points, we must first
find the locations of the dots within each image. We refer to these coordinate
pairs as points. The dots on the tube surface are drawn in white, are spaced
by approximately $1\text{\,}\mathrm{mm}$, and have a diameter of approximately
$0.7\text{\,}\mathrm{mm}$. We take the material points to be the centres of
the dots, and we find them by first taking a two-dimensional convolution of
the image with a “mexican hat” function of the form,

$\frac{1}{\pi\sigma^{4}}\mathinner{\left(1-\frac{x^{2}+y^{2}}{2\sigma^{2}}\right)}e^{-\frac{x^{2}+y^{2}}{2\sigma^{2}}},$
(13)

where $\sigma$ is the expected radius of the dot in the image. The convolution
effectively smooths the image, removing any artefacts from drawing the dots,
leaving peaks in the centroid of each dot. These peaks are then used as the
locations of each point.

Hence at each time instance, we have two collections of points, representing
the material points as seen from each camera. If there are $n$ points on the
tube, and $m$ frames in our video, then in total we will have found the
locations of $2nm$ points. To make use of this data, we need to identify each
unique material point in each camera and over time.

To associate points across time for a single camera’s set of images, the
points at $t$ and $t+\Delta t$ are compared, and if two points are within a
certain distance of each other, then it is assumed that these represent the
same point. This works because the frame rate ($12\,500\text{\,}\mathrm{fps}$)
is much larger than the frequency of the observed vibrations
($\sim$500\text{\,}\mathrm{Hz}$$), so motions between frames are small.

A pair of corresponding points in the two camera images are illustrated in
Figure 2, but finding these pairings at each instant in time is more complex.
First we consider the line drawn from the focal point of camera 2 to dot $i$
in the image, which we will call a ray. Anything on this ray in three-
dimensional space will appear at the same highlighted location in camera 2\.
However, from camera 1, the ray will appear as a line. Therefore, if a point’s
location is known in one camera image, then it must lie on a specific line in
the other image. This line is known as an epipolar line [25]. Hence, for a
pair of points, one in each camera image, to correspond to the same material
point, they must each lie on the epipolar line of the other. However, because
of the specific arrangement of the cameras and dots, this does not usually
provide a unique set of pairs. The relative positioning of the cameras means
that the epipolar lines are all approximately horizontal, and the dots drawn
on the tube are arranged in horizontal rows, so multiple dots can be very
close to a given epipolar line. To overcome this, we specify $10$
corresponding point pairs between the two images ($20$ points in total), and
these $10$ pairs are then used to find the best fitting projective
transformation from camera 1 to camera 2. Applying this transformation to the
image from camera 1 places each point close to the corresponding point in the
image from camera 2, allowing all the remaining point pairs to be found. This
result is then checked for consistency with the epipolar line condition.

In addition to tracing material points over time as the self excited
oscillations occur, we must also find the locations of the material points on
the tube when it is unstrained. This is necessary for the calculation of the
kinematic variables used in expression for energy of deformation (1). To
achieve this a single image is taken from each camera when the tube is held in
its unstrained state. Pairing of points must then also be completed between
the two images of the tube in its unstrained state and images from the high
speed video of self excited oscillations. This pairing is done using the
methods described in the previous paragraph.

Once point pairs are known over time, the camera calibration can be used to
find three-dimensional point traces over time. The spatial resolution of this
trace is limited by the size of the pixels in the high speed video. This
results in point traces with distinct jumps in position. These jumps are by no
more than $0.1\text{\,}\mathrm{mm}$ in three-dimensional space, compared to
variations in position on the order of $2\text{\,}\mathrm{mm}$ over the course
of the self excited oscillations. For this reason we smooth the three-
dimensional point traces by fitting functions of the form,

$\sum_{i=1}^{8}A_{i}\sin(\omega_{i}t),$ (14)

to the three position components, where $A_{i}$ and $\omega_{i}$ are chosen to
fit the empirical data. These fits work well because the videos are of quasi-
steady behaviour at the onset of oscillation, and the observed motions are
close to sinusoidal. $8$ terms have been found to be sufficient to match the
experimental data. This smoothing is necessary in order for derivatives of the
point traces to give meaningful results.

###  3.3 Kinematic Calculations

In Figure 3 we show a single frame from the video in which the dots on the
surface of the tube have been automatically detected, identified with the
corresponding dots in the other video image, and identified with the
corresponding dots in a stereoscopic image of the unstrained tube (not shown).
With this information we can reconstruct the points in three dimensions and
fit a surface to them. The process can be summarised as follows,

  * •

Locate points in images of unstrained tube.

  * •

Track points in high speed video frames.

  * •

Associate points between all images (as described in §3.2) and triangulate to
have the position of each material point as a function of time, and its
position on the unstrained tube.

  * •

Assign a pair of coordinate values $\\{x^{i}\\}$ to each point for use as the
convected surface coordinates.

  * •

Fit a smoothing curve to every three-dimensional point as a function of time
as described in §3.2.

  * •

At a chosen time take the positions of all of the points and fit a polynomial
surface such that we have position as a function of the surface coordinates
$\\{x^{i}\\}$. Repeat this for the unstrained surface.

  * •

Take all possible first and second derivatives of the surface position with
respect to $\\{x^{i}\\}$. This is done analytically using the polynomial
surface. Repeat this for the unstrained surface.

The surface fit and its derivatives are used to calculate all the kinematic
properties of the undeformed and deformed surfaces.

In Figure 3 we show a typical image of the principal strains of the flexible
tube, i.e. the eigenvectors and values of $\mathsf{E}$. An eigenvalue of $1$
corresponds to no strain, a value less than $1$ corresponds to compression,
and a value greater than $1$ corresponds to tension. The eigenvectors give the
direction in which the strain is occurring. We can see that the first
principal strain is approximately aligned with the longitudinal direction, and
is tensile. This is due to the dominant pre-strain of the elastic tube. By
contrast, the second principal strain, which is primarily in the azimuthal
direction, is a mixture of compression and tension, and is generally closer to
$1$. In Figure 4 we show the principal curvatures of the deformed surface,
i.e. the eigenvectors and values of $\mathsf{b}$. For a cylinder, which the
tube is in its undeformed state, the principal curvatures would be $0$ in the
longitudinal direction and $1/a=$0.33\text{\,}{\mathrm{mm}}^{-1}$$ in the
azimuthal direction. We see that in the deformed tube the curvatures are still
closely aligned to the longitudinal and azimuthal directions, and can see that
the squashing of the tube results in a slightly negative longitudinal
curvature and a slight reduction in the azimuthal curvature towards the centre
of the tube. But these effects are fairly small.

Figure 3: A typical image of the principal strains of the flexible tubes,
shown along with the principal strain directions (illustrated with unit
vectors). On the left the original high speed camera images are shown with the
tracked surface points, and a view of the three-dimensional triangulation with
the surface fitted to them. Figure 4: A typical image of the principal
curvatures of the flexible tube, shown along with the principal curvature
directions (illustrated with unit vectors).

##  4 Calculation of Bending and Stretching Energies

We now aim to gain more understanding of how the tube deforms. To do this we
will use a mixture of empirical and analytical techniques. More specifically
we can use the high speed video reconstructions combined with the mathematical
framework for shells already developed in §2.

###  4.1 Scaling Analysis

We start by estimating the bending and stretching energies analytically, which
requires us to estimate the values of $\operatorname{tr}(\mathsf{E}^{2})$,
$\operatorname{tr}(\mathsf{E})^{2}$, $\operatorname{tr}(\mathsf{H}^{2})$, and
$\operatorname{tr}(\mathsf{H})^{2}$.

If $\alpha_{i}$ are the eigenvalues of the linear function $\mathsf{A}$, then
$\operatorname{tr}(\mathsf{A})^{2}=\alpha_{1}^{2}+2\alpha_{1}\alpha_{2}+\alpha_{2}^{2}$
and $\operatorname{tr}(\mathsf{A}^{2})=\alpha_{1}^{2}+\alpha_{2}^{2}$. We know
that the eigenvalues of $\mathsf{E}$ are $\frac{1}{2}(\lambda_{i}^{2}-1)$
($\lambda_{i}$ are the principal strains) and that $\lambda_{1}\sim\lambda$,
$\lambda_{2}\sim 1$, where $\lambda$ is the initial axial strain of the tube.
This allows us to get an order of magnitude estimate for
$\operatorname{tr}(\mathsf{E}^{2})$ and $\operatorname{tr}(\mathsf{E})^{2}$ of
$(\lambda^{2}-1)^{2}/4$. If we take $\lambda=1.3$333The strained length of the
tube is $l=$25\text{\,}\mathrm{mm}$$ and the unstrained length is
$l_{0}=$19\text{\,}\mathrm{mm}$$, giving $\lambda=l/l_{0}\approx 1.3$. then we
have
$\operatorname{tr}(\mathsf{E}^{2})\sim\operatorname{tr}(\mathsf{E})^{2}\sim$0.1$$.

To understand bending in the tube we consider the representation of
$\mathsf{H}$ derived in §2, and given in (12). We can write the action of
$\mathsf{U}(Y)$ as
$\mathsf{U}(Y)=\lambda_{1}(Y\cdot\hat{W}_{1})\hat{W}_{1}+\lambda_{2}(Y\cdot\hat{W}_{2})\hat{W}_{2}$,
where $\hat{W}_{i}$ are the unit eigenvectors of $\mathsf{U}$. Using the
approximation $\lambda_{1}\approx\lambda$ and $\lambda_{2}\approx 1$ this
becomes
$\mathsf{U}(Y)\approx\lambda(Y\cdot\hat{W}_{1})\hat{W}_{1}+(Y\cdot\hat{W}_{2})\hat{W}_{2}$.
We can write $\mathsf{B}(Y)$ as
$\mathsf{B}(Y)=C(Y\cdot\hat{E}_{2})\hat{E}_{2}$, where $\hat{E}_{i}$ are the
unit eigenvectors of $\mathsf{B}$, and $C$ is the principal curvature of the
undeformed tube in the azimuthal direction. Combining these we have,

$\displaystyle\mathsf{U}\mathsf{B}(Y)-\mathsf{B}(Y)$
$\displaystyle=\mathsf{U}(CY\cdot\hat{E}_{2}\hat{E}_{2})-CY\cdot\hat{E}_{2}\hat{E}_{2}$
(15) $\displaystyle\approx
CY\cdot\hat{E}_{2}(\lambda\hat{E}_{2}\cdot\hat{W}_{1}\hat{W}_{1}+\hat{E}_{2}\cdot\hat{W}_{2}\hat{W}_{2})-CY\cdot\hat{E}_{2}\hat{E}_{2}.$

We know from Figure 3 that $\hat{W}_{1}$ and $\mathsf{W}_{2}$ are
approximately aligned with the longitudinal and circumferential directions, so
we can write $\hat{E}_{2}\cdot\hat{W}_{1}\approx 0$ and
$\hat{E}_{2}\cdot\hat{W}_{2}\approx 1$. Using this we obtain,

$\mathsf{U}\mathsf{B}(Y)-\mathsf{B}(Y)\approx
CY\cdot\hat{E}_{2}\hat{W}_{2}-CY\cdot\hat{E}_{2}\hat{E}_{2}\approx 0.$ (16)

This is saying that because the directions of principal strain and principal
curvature are approximately perpendicular, the influence of strain on change
of curvature is removed.

Figure 5: A visualisation of the rotation $\mathsf{R}$. The axis of rotation
is shown along with the absolute value of the angle of rotation (in radians).

We are now in a position to consider the rotor $R$, since it is changes in
this multivector over the surface of the shell that are responsible for the
change of curvature. $R$ represents the rotation, and as is shown in §2 is
characterised by the bivector $A=\theta\hat{A}$, whose magnitude gives the
rotation angle in radians, and whose plane gives the plane of rotation. In
Figure 5 we visualise the angle and axis of rotation encoded in the rotation
tensor $\mathsf{R}$, corresponding to a typical deformation. The axes of
rotation shown in Figure 5 are primarily tangential to the surface, so the
bivector $A$ will be dominated by the components $e_{1}\wedge e_{3}$ and
$e_{2}\wedge e_{3}$, with little rotation in the $e_{1}\wedge e_{2}$ plane,
i.e. about the normal vector $e_{3}$. Hence, we can write $A$ as,

$A=\theta_{1}e^{1}\wedge e_{3}+\theta_{2}e^{2}\wedge
e_{3}=\theta_{i}e^{i}\wedge e_{3}.$ (17)

We have used the reciprocal frame $\\{e^{i}\\}$ instead of of $\\{e_{i}\\}$
because it will allow us to use the property
$\bar{\mathsf{F}}(e^{i})=E^{i}$444The coordinate system is convected, so
$e_{i}=\mathsf{F}(E_{i})$, and as always $e^{i}\cdot e_{j}=E^{i}\cdot
E_{j}=\delta^{i}_{j}$. Hence
$\delta^{i}_{j}=e^{i}\cdot\mathsf{F}(E_{j})=\bar{\mathsf{F}}(e^{i})\cdot
E_{j}$, from which it is clear that $\bar{\mathsf{F}}(e^{i})=E^{i}$.. We can
extend the frames $\\{e^{i}\\}$ and $\\{e_{i}\\}$ to span $\mathbb{E}^{3}$ by
using the normal vector $e_{3}$. Because $e_{3}$ is a unit vector and
perpendicular to all of $e_{i},e^{i}$, we can also write $e^{3}=e_{3}$, and we
have the frame $\\{e_{a}\\},a=1,2,3$ and $\\{e^{a}\\}$. Using this we define
the Christoffel coefficients $\gamma^{a}_{ib}=e^{a}\cdot\partial
e_{b}/\partial x^{i}$, $i=1,2;a,b=1,2,3$. These also satisfy
$e_{b}\cdot\partial e^{a}/\partial x^{i}=-\gamma^{a}_{ib}$.

Substituting $A$ into the second part of the change of curvature tensor given
in (12), using the fact that $e_{3}$ is normal to $e^{1}$ and $e^{2}$, and
$\gamma^{3}_{i3}=0$, we obtain,

$\displaystyle\bar{\mathsf{F}}\mathinner{\left(e_{3}\cdot\mathinner{\left(Y\cdot\partial
A\right)}\right)}$
$\displaystyle=Y^{i}\bar{\mathsf{F}}\mathinner{\left(-\partial_{i}(\theta_{j})e^{j}+\theta_{j}\gamma^{j}_{ik}e^{k}\right)}$
(18)
$\displaystyle=-Y^{i}\partial_{i}(\theta_{j})E^{j}+Y^{i}\theta_{j}\gamma^{j}_{ik}E^{k}.$

Therefore, given that the first part of $\mathsf{H}$ in (12) is zero,
$E_{i}\cdot\mathsf{H}(E_{j})=\mathsf{H}_{ij}$ is given by,

$\mathsf{H}_{ij}=\partial_{j}\theta_{i}-\theta_{k}\gamma^{k}_{ji}.$ (19)

We know that $\mathsf{H}$ is symmetric, so from this we see that our earlier
assumption on the form of $A$ must be joined by the condition
$\partial_{i}\theta_{j}=\partial_{j}\theta_{i}$ to produce consistent results.

The frame $\\{E_{i}\\}$ can be chosen by us to be orthonormal. More
specifically, we can align $E_{1}$ with the longitudinal direction and $E_{2}$
with the azimuthal direction on the unstrained cylindrical tube. From Figure 3
we can see that under a typical deformation these basis vectors remain close
to the axial and azimuthal directions. In Figure 6 we give a schematic
illustration of how $E_{i}$ maps to $e_{i}$. In Figure 6 we have also labelled
the values of $\theta_{i}$ where they are obvious. The regions where
$\mathinner{\\!\left\lvert\theta_{2}\right\rvert}\approx\pi/4$ can be seen
empirically in Figure 5. There are two unknown values of $\theta_{1}$ shown as
question marks. We can estimate the largest value of these rotations by
assuming a straight line from the clamped tube end and the centre of the tube
when the tube collapses completely at the centre. In this case
$\theta_{1}\sim\arctan(a/(l/2))\lambda$ where $a$ is the tube radius and $l$
is the tube length in its deformed state. The multiplication by $\lambda$ is
necessary because $\theta_{1}$ is the $e^{1}\wedge e_{3}$ component, and
$e^{1}$ is shortened by a factor of $\lambda$ compared to the unit vector
$E^{1}$. Up to angles of $30\text{\,}\mathrm{\SIUnitSymbolDegree}$,
$\tan\theta$ is within $10\text{\,}\mathrm{\char 37\relax}$ of $\theta$, so we
will take $\theta_{1}\sim\lambda a/(l/2)=a/(l_{0}/2)$ at the point in
question, where $l_{0}$ is the unstrained length of the tube. This, and the
values of $\theta_{i}$ shown in Figure 6, allows us to make the following
estimates,

$\displaystyle\partial_{1}\theta_{1}$
$\displaystyle\sim\frac{a}{l_{0}/2}\frac{1}{l_{0}/2}=\frac{4a}{l_{0}^{2}},$
(20) $\displaystyle\partial_{1}\theta_{2}=\partial_{2}\theta_{1}$
$\displaystyle\sim\frac{\pi/4}{l_{0}/2}=\frac{\pi}{2l_{0}},$
$\displaystyle\partial_{2}\theta_{2}$ $\displaystyle\sim\frac{\pi/4}{2\pi
a/8}=\frac{1}{a}.$

Because of replacement of $\arctan$ with the identity function, our estimate
for $\partial_{1}\theta_{1}$ will be an overestimate when the tube is very
short.

Figure 6: A schematic of the deformation of the flexible tube.

We can also estimate the values of the coefficients $\gamma^{i}_{jk}$ using
Figure 6,

$\displaystyle\partial_{1}e_{1}$
$\displaystyle\sim\frac{a}{(l_{0}/2)^{2}}e_{3},$ (21)
$\displaystyle\partial_{2}e_{1}=\partial_{1}e_{2}$ $\displaystyle\sim 0,$
$\displaystyle\partial_{2}e_{2}$ $\displaystyle\sim\frac{1}{a}e_{3}.$

From this we see that the changes in the basis vectors are primarily in the
$e_{3}$ direction, meaning that they do not contribute to $\mathsf{H}_{ij}$.

If we take $l_{0}$ to be much larger than $a$, then the dominant term in
$\mathsf{H}_{ij}$ will be $1/a$, but even if $l_{0}$ and $a$ are a similar
order of magnitude, all of the $\mathsf{H}_{ij}$ terms will be of the order
$1/a$. Hence, we expect $\operatorname{tr}(\mathsf{H}^{2})$ and
$\operatorname{tr}(\mathsf{H})^{2}$ will scale as
$(1/a)^{2}=(1/$3\text{\,}\mathrm{mm}$)^{2}=$0.1\text{\,}{\mathrm{mm}}^{-2}$$.

Given these scalings for $\operatorname{tr}(\mathsf{E}^{2})$,
$\operatorname{tr}(\mathsf{E})^{2}$, $\operatorname{tr}(\mathsf{H}^{2})$, and
$\operatorname{tr}(\mathsf{H})^{2}$, and the values
$E=$1\text{\,}\mathrm{MPa}$$, $\nu=0.5$, and $h=$0.3\text{\,}\mathrm{mm}$$, we
can obtain scalings for the bending and stretching energy given in (1),

stretching energy
$\displaystyle\sim$0.02\text{\,}\mathrm{N}\text{\,}{\mathrm{mm}}^{-1}$,$ (22)
bending energy
$\displaystyle\sim$1.5\text{\times}{10}^{-4}\text{\,}\mathrm{N}\text{\,}{\mathrm{mm}}^{-1}$,$

This indicates that given the kind of deformation we have observed in our
Starling resistors at onset, i.e. where the strain energy is dominated by the
effects of pre-strain, the axis aligned with the largest strain remains close
to perpendicular to the axis aligned with the largest curvature, rotations are
mostly about axes tangential to the shell, and changes in the rotation scale
with the change of rotation about the longitudinal axis in the azimuthal
direction, stretching energy will dominate bending energy. Moreover, this
result remains valid even when the tube length gets close to the tube
diameter. This is significant for our considerations of the lung, since the
length to diameter ratio of tubes in the lung typically varies from $1$ to $6$
[27].

###  4.2 Direct Calculation from Data

We can use the high speed video data to calculate
$\operatorname{tr}(\mathsf{E}^{2})$, $\operatorname{tr}(\mathsf{E})^{2}$,
$\operatorname{tr}(\mathsf{H}^{2})$, and $\operatorname{tr}(\mathsf{H})^{2}$
from (1). Typical plots are shown in Figure 7, from which we see that in the
units chosen these have similar orders of magnitude. We also see that the
scalings obtained in the previous section agree with these plots well,
providing support for the assumptions made. We can also calculate the bending
and stretching energy, and this is shown in Figure 8. This agrees very well
with the scaling values of the previous section, again supporting our
conclusions.

Figure 7: Values of the kinematic variables
$\operatorname{tr}(\mathsf{E}^{2})$, $\operatorname{tr}(\mathsf{E})^{2}$,
$\operatorname{tr}(\mathsf{H}^{2})$, and $\operatorname{tr}(\mathsf{H})^{2}$.
Figure 8: The stretching energy (top image) and bending energy (bottom image)
associated with the deformation of the flexible tube.

##  5 Conclusions

We have developed a new method of representing the change of curvature tensor
using rotors (see (12)), increasing our understanding of bending in shells. We
have used this representation to explain results from stereographic imaging of
Starling resistors that demonstrate that the bending energy in these
deformations is around 2 orders of magnitude lower than the stretching energy.
We have been able to show that this relies on the fact that the strain energy
is dominated by the effects of pre-strain, the axis aligned with the largest
strain remains close to perpendicular to the axis aligned with the largest
curvature, rotations are mostly about axes tangential to the shell, and
changes in the rotation scale with the change of rotation about the
longitudinal axis in the azimuthal direction. Further to this, our scaling
analysis remains valid even when the tube length gets close to the tube
diameter. This is of significance to our work in understanding wheezing, since
the length to diameter ratio of tubes in the lung typically varies from $1$ to
$6$. Hence we have provided a scaling analysis, confirmed by experiment, that
allows us to say that bending energy is dominated by stretching energy during
self excited oscillations in the airways of the lung. This should allow the
use of membrane theory to model the tube, which reduces the order of the
equations of motion from $4$ to $2$.

#### Ethics

This paper adheres to the Royal Society publishing and research ethics policy.

#### Data Accessibility

Data for this paper can be accessed at
https://doi.org/10.17863/CAM.10363.

#### Author Contributions

The work was completed as part of the PhD project of AG, with AA supervising,
and JL advising, in particular on image processing and the uses of geometric
algebra.

#### Competing Interests

The authors declare no competing interests.

#### Funding

The authors would like to acknowledge funding from the EPSRC, the IMechE
Postgraduate Research Scholarship, and Engineering for Clinical Practice
(http://divf.eng. cam.ac.uk/ecp/Main/EcpResearch).

#### Acknowledgements

Our lab technician John Hazelwood, who has been instrumental in the
experimental work completed, and Holger Babinsky, who provided advice on the
use of high speed cameras.

## References

  * Heil & Jensen [2003] Heil, M. & Jensen, O. E. 2003 Flows in deformable tubes and channels.  In _Flow past Highly Compliant Boundaries and in Collapsible Tubes_ (eds. P. W. Carpenter & T. J. Pedley), pp. 15–49. Kluwer Academic Publishers.
  * Bertram [2003] Bertram, C. D. 2003 Experimental Studies of Collapsible Tubes.  In _Flow past Highly Compliant Boundaries and in Collapsible Tubes_ (eds. P. W. Carpenter & T. J. Pedley), pp. 51–65. Kluwer Academic Publishers.
  * Grotberg & Jensen [2004] Grotberg, J. B. & Jensen, O. E. 2004 Biofluid Mechanics in Flexible Tubes.  _Annu. Rev. Fluid Mech._ , 36(1), 121–147.
  * Heil & Hazel [2011] Heil, M. & Hazel, A. L. 2011 Fluid-Structure Interaction in Internal Physiological Flows.  _Annu. Rev. Fluid Mech._ , 43(1), 141–162.
  * Kisilova _et al._ [2012] Kisilova, N., Hamadiche, M. & Gad-El-Hak, M. 2012 Mathematical models of biofluid flows in compliant ducts.  _Arch. Mech._ , 64(1), 65–94.
  * Jensen [2013] Jensen, O. E. 2013 Instabilities of flows through deformable tubes and channels.  In _Mechanics Down Under - Proceedings of the 22nd International Congress of Theoretical and Applied Mechanics, ICTAM 2008_ , pp. 101–116.
  * Grotberg & Davis [1980] Grotberg, J. B. & Davis, S. H. 1980 Fluid-dynamic flapping of a collapsible channel: sound generation and flow limitation.  _J. Biomechanics_ , 13(3), 219–230.
  * Forgacs [1978] Forgacs, P. 1978 _Lung Sounds_.  Baillière Tindall.
  * Bohadana _et al._ [2014] Bohadana, A. B., Izbicki, G. & Kraman, S. S. 2014 Fundamentals of lung auscultation.  _N Engl J Med_ , 370(8), 744–751.
  * Whittaker _et al._ [2010] Whittaker, R. J., Heil, M., Jensen, O. E. & Waters, S. L. 2010 A Rational Derivation of a Tube Law from Shell Theory.  _Q. J. Mech. Appl. Math._ , 63(4), 465–496.
  * Shabtai-Musih _et al._ [1992] Shabtai-Musih, Y., Grotberg, J. B. & Gavriely, N. 1992 Spectral content of forced expiratory wheezes during air, He, and SF6 breathing in normal humans.  _J. Appl. Physiol._ , 72(2), 629–635.
  * Leissa [1973] Leissa, A. W. 1973 Vibration of shells.  Tech. Rep. NASA SP-288, National Aeronautics and Space Administration, Washington D.C.
  * Naghdi [1972] Naghdi, P. M. 1972 The theory of shells and plates.  In _Handb. der Phys. (Encyclopedia Physics), Band (Volume) VIa/2, Festkorpermechanik II (Mechanics Solids II)_ (eds. S. Flugge & C. Truesdell), pp. 425–640. Springer-Verlag.
  * Ciarlet Jr [2005] Ciarlet Jr, P. G. 2005 An introduction to differential geometry with applications to elasticity.  _J. Elast._ , 78-79, 1–215.
  * Antman [2005] Antman, S. S. 2005 _Nonlinear Problems of Elasticity_.  Berlin: Springer-Verlag, 2nd edn.
  * Lacarbonara [2012] Lacarbonara, W. 2012 The Nonlinear Theory of Plates.  In _Nonlinear Structural Mechanics_ , pp. 497–592. Boston, MA: Springer US.
  * Koiter [1966] Koiter, W. T. 1966 On the nonlinear theory of thin elastic shells.  _Koninklijke Nederlandse Akademie van Wetenschappen_ , 69B, 1–54.
  * Gregory _et al._ [2017] Gregory, A. L., Lasenby, J. & Agarwal, A. 2017 The elastic theory of shells using geometric algebra.  _R. Soc. Open Sci._ , 4(3), 170 065.
  * Doran & Lasenby [2003] Doran, C. J. L. & Lasenby, A. N. 2003 _Geometric Algebra for Physicists_.  Cambridge: Cambridge University Press.
  * McRobie & Lasenby [1999] McRobie, F. A. & Lasenby, J. 1999 Simo-Vu Quoc rods using Clifford algebra.  _Int. J. Numer. Meth. Eng._ , 45, 377–398.
  * Simo & Vu-Quoc [1988] Simo, J. C. & Vu-Quoc, L. 1988 On the dynamics in space of rids undergoing large motions - a geometrically exact approach.  _Comput. Method. Appl. M._ , 66, 125–161.
  * Lasenby _et al._ [1998] Lasenby, A., Doran, C. & Gull, S. 1998 Gravity, gauge theories and geometric algebra.  _Phil. Trans. R. Soc. Lond. A_ , 356(1737), 487–582.
  * Gregory _et al._ [2015] Gregory, A. L., Sinayoko, S., Agarwal, A. & Lasenby, J. 2015 An acoustic space-time and the Lorentz transformation in aeroacoustics.  _Int. J. Aeroacoust._ , 14(7), 977–1003.
  * Dorst _et al._ [2009] Dorst, L., Fontijne, D. & Mann, S. 2009 _Geometric Algebra for Computer Science_.  An Object-Oriented Approach to Geometry. Morgan Kaufmann, 2nd edn.
  * Hartley [2004] Hartley, R. 2004 _Multiple View Geometry in Computer Vision_.  Cambridge University Press, 2nd edn.
  * MATLAB [2017] MATLAB 2017 _version 9.2.0 (R2017a)_.  Natick, Massachusetts: The Mathworks Inc.
  * Horsefield & Cumming [1968] Horsefield, K. & Cumming, G. 1968 Morphology of the bronchial tree in man.  _J. Appl. Physiol._ , 24(3), 373–383.