# Early dynamical evolution of rotating star clusters in a tidal field

Maria A. Tiongco,1,2 Enrico Vesperini,2 and Anna Lisa Varri3,4
1JILA and Department of Astrophysical and Planetary Sciences, University of
Colorado, Boulder, CO 80309, USA
2Department of Astronomy, Indiana University, Bloomington, IN 47405, USA
3Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
4School of Mathematics, University of Edinburgh, Kings Buildings, Edinburgh
EH9 3JZ, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In order to explore how the early internal rotational properties of star
clusters are affected by the external potential of their host galaxies, we
have run a suite of $N$-body simulations following the early dynamical
evolution and violent relaxation of rotating star clusters embedded in a tidal
field. Our study focuses on models for which the cluster’s rotation axis has a
generic orientation relative to the torque of the tidal field. The interaction
between the violent relaxation process, angular momentum of the cluster, and
the external torque creates a complex kinematic structure within the cluster,
most prominently a radial variation in the position of the rotation axis along
both the polar and azimuthal directions. We also examine the cluster’s
velocity dispersion anisotropy and show that the projected anisotropy may be
affected by the variation of the rotation axis directions within the cluster;
the combination of projection effects and the complex kinematical features may
result in the measurement of tangential anisotropy in the cluster’s inner
regions. We also characterize the structural properties of our clusters as a
function of their initial rotation and virial ratio and find that clusters may
develop a triaxial morphology and a radial variation of the minor axis not
necessarily aligned with the rotation axis. Finally, we examine the long-term
evolution of these complex kinematic features.

###### keywords:

galaxies: star clusters: general – methods: numerical

††pubyear: 2015††pagerange: Early dynamical evolution of rotating star
clusters in a tidal field–Early dynamical evolution of rotating star clusters
in a tidal field

##  1 Introduction

Bulk internal rotation of star clusters is an important component of the
dynamical picture of these systems. On the theoretical side, several
pioneering studies of the evolution of rotating star clusters have
demonstrated how two-body relaxation over time redistributes angular momentum
from the inner regions of the cluster outward, and causes stars to escape and
carry away angular momentum, lowering the rotation in the cluster (see, e.g.,
Einsel & Spurzem 1999; Kim et al. 2002; Kim et al. 2004, 2008; Ernst et al.
2007; Hong et al. 2013; Tiongco et al. 2017). Most of the numerical studies
that included the effects of a tidal field explored the geometrically simplest
case where the cluster’s angular momentum vector is parallel or anti-parallel
to the torque of the host galaxy (e.g., see Boily et al. 1999), but additional
kinematic complexities have been shown to arise when more general
configurations are considered (Tiongco et al. 2018).

A recent dramatic increase in the number and quality of observational
kinematic studies has led to major progress towards a complete kinematic
portrait of star clusters. Such studies draw from proper-motion data from Gaia
(see, e.g, Bianchini et al. 2018; Sollima et al. 2019; Vasiliev 2019), HST
(Bellini et al. 2017), and line-of-sight velocities from spectroscopic studies
using Integral Field and/or multi-fiber spectrographs (see, e.g., Fabricius et
al. 2014; Kamann et al. 2018; Boberg et al. 2017; Ferraro et al. 2018; Lanzoni
et al. 2018a, b; Wan et al. 2021). In general agreement with the theoretical
models predicting a decrease of the strength of rotation over time,
observational studies are starting to show that globular clusters with longer
relaxation times (i.e., dynamically younger clusters) tend to have stronger
rotation (Kamann et al. 2018; Bianchini et al. 2018). With more detailed
kinematic studies combining proper motion data and spectroscopic data to get
three-dimensional velocities, we are closer to being able to build a complete
picture of the phase space of star clusters (see, e.g., Bellini et al. 2017;
Vitral & Mamon 2021; Dalessandro et al. 2021a).

The dynamical origin of the rotation observed in star clusters is still a
matter of investigation; a number theoretical studies have shown that angular
momentum inherited from the molecular cloud out of which the cluster formed or
acquired from external gas accreted during the cluster formation, dynamical
interactions among the various subclumps in the primordial cluster, and the
effects of the host galaxy tidal field may all play an important role for the
origin of clusters’ internal rotation (Vesperini et al. 2014; Mapelli 2017;
Tiongco et al. 2017; Lahén et al. 2020; Ballone et al. 2021; Chen et al.
2021).

In this study, we broaden the classic studies of violent relaxation of
isolated systems (e.g., see Lynden-Bell 1967; van Albada 1982; Aarseth et al.
1988 and also Gott 1973; Hohl & Zang 1979; Aguilar & Merritt 1990 for rotating
systems), to consider such a process in the case of rotating clusters in the
external tidal field of their host galaxy and explore the kinematic evolution
in the general case in which the cluster’s internal angular momentum vector is
not aligned with the orbital angular momentum vector. One of the issues we
will explore concerns the origin of the radial variation of the position angle
of rotation axis observed in some globular clusters (Gebhardt et al. 2000;
Bianchini et al. 2013; Boberg et al. 2017; Usher et al. 2021). In Tiongco et
al. (2018), we followed the long-term evolution of rotating star clusters in
an external tidal field and rotation axes initially misaligned with their
orbital angular velocity vectors. Because rotation of the cluster due to tidal
synchronization occurs parallel to orbital angular velocity vector, the
cluster may develop a radial variation in the direction of the rotation axis,
due to the combination of cluster’s intrinsic rotation and rotation that
develops from evolution towards tidal locking. In this study we explore
whether a radial variation in the rotation axis’ orientation can arise also
during the violent relaxation phase.

Many new and upcoming studies are also focusing on the kinematics of young
star clusters, including cluster members and stars within the cluster’s
vicinity that may be former members of the cluster (see, e.g., Kuhn et al.
2019; Armstrong et al. 2020; Lim et al. 2020; Jerabkova et al. 2021;
Dalessandro et al. 2021b). Because these star clusters are young, the escaped
stars may carry some memory of the cluster’s early dynamical history.
Motivated by these observational studies, we will also focus our attention on
the differences in the structural and kinematical properties between the
cluster’s inner, outer, and extratidal regions and on how these complexities
may appear in observations along generic lines of sight.

This paper is structured as follows: in Section 2, we discuss the setup of our
$N$-body simulations, in Section 3, we present the results of our simulations,
and in Section 4, we conclude with a summary of our study.

##  2 Method and Initial Conditions

We have run $N$-body simulations of star clusters on circular orbits about a
point-mass galaxy, using the GPU-accelerated $N$-body code, nbody6 (Aarseth
2003; Nitadori & Aarseth 2012). The equations of motion are solved in a co-
rotating reference frame centred on the cluster (Heggie & Hut 2003), rotating
with angular speed $\Omega$, equal to that of the angular speed of the
cluster’s orbit. The $x$-$y$ plane of coordinate system is coincident with the
cluster’s orbital plane, where the positive $x$-axis of the rotating reference
frame always points away from the galactic centre, and the positive $y$-axis
points in the direction of the cluster’s motion. The cluster’s orbital angular
velocity vector is parallel to the positive $z$-axis.

Similar to Tiongco et al. (2018, 2019), the analysis of our results is
executed in a non-rotating reference frame centred on the cluster. We will
refer to the angular velocity vector of the cluster, $\boldsymbol{\omega}$,
and we also define this vector’s direction in space using the angles
$\phi_{\omega}$ and $\theta_{\omega}$: $\phi_{\omega}$ is the angle between
the projection of $\boldsymbol{\omega}$ in the $x$-$y$ plane (the orbital
plane) and the positive $x$-axis, measured in a counter-clockwise direction
from the positive $x$-axis, and $\theta_{\omega}$ is the angle between
$\boldsymbol{\omega}$ and the $z$-axis, measured starting from the positive
$z$-axis (the axis pointing outwards from the orbital plane).

Our simulations follow the evolution of star clusters initially with some
rotation and also undergoing a phase of violent relaxation. The violent
relaxation phase is set up by setting the system’s initial virial ratio,
$Q=T/|V|$ (where $T$ and $V$ are the kinetic energy and potential energy of
the system, respectively) to a value less than 0.5, so that the system is not
in virial equilibrium. We also decompose $Q$ into its rotational and random
components, $Q_{\rm rot}$ and $Q_{\rm ran}$, where the kinetic energy of each
component is from bulk rotation of the system or random motions of the
particles, respectively (we adopt the nomenclature of the components from Gott
1973). The clusters are initially spherical with uniform density, with random
velocities drawn from a Gaussian distribution, and rotation following a solid-
body rotational profile. The initial size of our system is defined by $r_{\rm
L}$, the limiting radius, and we set the ratio of $r_{\rm L}$ to the Jacobi
radius, $r_{\rm J}$, to 0.5. Notable to our study, we also set the initial
orientation of the rotation axis (also defined as the direction of
$\boldsymbol{\omega}$) relative to the tidal field as a free parameter.

Table 1 summarizes the initial parameters for our simulations. We note that
$Q_{\rm rot}$ is given in the non-rotating reference frame. Also, we do not
aim to characterize in detail the evolution of systems that begin with
$\theta_{\rm\omega,i}$ = 0 or 180, i.e., parallel and anti-parallel with the
orbital angular velocity vector, as the evolution of these configurations is
straightforward and has been studied by the references mentioned in the
introduction.

The simulations are run with 131 072 equal-mass stars, for approximately 18
free-fall times:

$t_{\rm ff}=\sqrt{\frac{3\pi}{32G\rho}},$ (1)

(see, e.g., Binney & Tremaine 2008), except for simulation ‘Cold th45’ which
we have run until approximately 9 initial half-mass relaxation times:

$t_{\rm rh,i}=\frac{0.138N^{1/2}r_{\rm h}^{3/2}}{\langle
m\rangle^{1/2}G^{1/2}\log(0.11N)},$ (2)

(see, e.g., Heggie & Hut 2003), i.e., until core-collapse of the system, for a
preliminary study illustrating the long-term evolution of these clusters. For
part of our analysis we will also include stars outside of the Jacobi radius,
$r_{\rm J}=(GM/3\Omega^{2})^{1/3}$. Stars are removed from the simulation when
they move beyond a radius equal to two times $r_{\rm J}$.

Table 1: Information about the $N$-body simulations presented in this study. The symbol $\theta_{\rm\omega,i}$ denotes the initial value of the angle between the direction of the global angular velocity vector and the $z$-axis, measured starting from the $z$-axis (which points outwards from the orbital plane). For further details, see Section 2. Model ID | $Q_{\rm rot}$ | $Q_{\rm ran}$ |  $\theta_{\rm\omega,i}$ (degrees)
---|---|---|---
Even th45 | 0.1 | 0.1 | 45
Even th90 | " | " | 90
Even th135 | " | " | 135
Cold th45 | 0.1 | 0.05 | 45
Cold th90 | " | " | 90
Cold th135 | " | " | 135
Slow th45 | 0.05 | 0.1 | 45
Slow th90 | " | " | 90
Slow th135 | " | " | 135

##  3 Results

###  3.1 General properties and visualization of the violent relaxation
process

Figure 1: Final snapshot of the Even th45, Cold th45, and Slow th45
simulations (each simulation is featured in a set of 6 panels spanning two
columns). Shown in columns 1, 3, and 5 are the projected number density
(normalized to the densest region) along 3 different projections, and in
columns 2, 4, and 6, the line-of-sight (los) velocities in the corresponding
projections to the projected number density images, with the contours of the
projected number density overlaid; the levels of each contour (outermost to
innermost) are 10, 25, 125, and 500 stars per pixel. The los velocities are
normalized to $\Omega r_{\rm J,i}$, the speed of co-rotation at the initial
Jacobi radius. The vertical and horizontal axes are normalized to the 3d half-
mass radius. The red dotted line points from the centre of the cluster to the
galactic centre. White pixels in the bottom panels correspond to values either
outside the range of the colourbar, or an absence of stars in that pixel. Only
stars inside the Jacobi radius are represented in these images. Figure 2:
Final snapshot of the Cold th90 simulation, similar to Fig. 1, but zoomed out
to emphasize the difference in the kinematics between the inner and outer and
extratidal regions, e.g., see the bottom middle panel. The red circle
represents the Jacobi radius of the cluster, while the blue lines in the
leftmost panels represent a more accurate shape of the Jacobi surface in that
projection; the axis ratios of the Jacobi surface are approximately 2/3 for
the middle/major axis and 0.5 for the minor/major axis (see, e.g., Heisler &
Tremaine 1986). Figure 3: Time evolution of $\phi_{\omega}$ (the $\phi$
component of the angular velocity vector, given in degrees) as a function of
radius normalized to the half-mass radius, $r_{\rm h}$, of the cluster at the
time specified, shown for all of our models. The small circles connected with
a dashed line are $\phi_{\omega}$ calculated from stars outside of the Jacobi
radius. Times in the legend are given in free-fall times, $t_{\rm ff}$.
Typical error bands are shown for the final snapshot and are derived using
bootstrap sampling of the snapshot.  Figure 4: Time evolution of
$\theta_{\omega}$ (the $\theta$ component of the angular velocity vector,
given in degrees) as a function of radius normalized to the half-mass radius,
$r_{\rm h}$, of the cluster at the time specified, shown for all of our
models. The small circles connected with a dashed line are $\theta_{\omega}$
calculated from stars outside of the Jacobi radius. Times in the legend are
given in free-fall times, $t_{\rm ff}$. Typical error bands are shown for the
final snapshot and are derived using bootstrap sampling of the snapshot.
Figure 5: Time evolution of angular velocity vector magnitude multiplied by
the radius it was measured at to give a rotational velocity at that radius
(radius is normalized to the half-mass radius, $r_{\rm h}$). Figs. 3 and 4
show the direction of the angular velocity vector, $\boldsymbol{\omega}$, and
the colors of the lines correspond to the same times in these mentioned
figures. The rotational velocities are also normalized to the velocity of the
rotating coordinate system at the initial Jacobi radius. Typical error bands
are shown for the final snapshot and are derived using bootstrap sampling of
the snapshot. Figure 6: Averaged radial profiles (where radius is normalized
to the half-mass radius, $r_{\rm h}$) of the axial ratios of all of models.
The solid lines show the median profiles of the snapshots between
17.1–17.9$t_{\rm ff}$ and the dashed lines show the maximum and minimum
profiles over this timespan. Figure 7: For model ‘Cold th90’, a comparison of
the profiles of $\phi_{\omega}$ and $\theta_{\omega}$ with $\phi$ and $\theta$
of minor axis ($\phi_{\rm c}$ and $\theta_{\rm c}$, all in degrees) at the end
of the simulation. We find this misalignment between the rotation axis and the
minor axis appears in our other models and throughout the timescale of the
simulations. Figure 8: Velocity dispersion anisotropy radial profiles of model
‘Slow th90’ using various binning prescriptions and lines of sight (los) at
the end of the simulation. The top left panel shows the profile of the
traditional anisotropy parameter, $\beta$ (see Eq. 4), calculated in spherical
shells. In the remaining panels, the solid line represents $\beta$ calculated
in cylindrical shells aligning them parallel to the los stated in the legend
of the panel. The bottom-middle and bottom-right panel use a los parallel to
the angular velocity vector $\boldsymbol{\omega}$ determined from stars within
1.5$r_{\rm h}$ and greater than 1.5$r_{\rm h}$, respectively. The dashed lines
represent the projected anisotropy parameter $1-\sigma_{\rm T}/\sigma_{\rm
R}$, where the velocity dispersions are calculated from velocities projected
onto the plane perpendicular to the los. $R$ and $R_{\rm h}$ refer to the
cylindrical radius and projected half-mass radius. This model was chosen to
show the similarities between the top left panel (‘los:x’) and the bottom
middle panel (‘los:r<1.5rh’) as this model was initiated with its
$\boldsymbol{\omega}$ parallel to the $x$-axis.

Figure 9: continued

Velocity dispersion anisotropy at different lines of sight, different binning
prescriptions, and different anisotropy prescriptions (see previous set of
panels) for model ‘Slow 135.’ This model was chosen because its initial
$\boldsymbol{\omega}$ is not aligned with the $x$-, $y$-, or $z$-axes.

Figure 10: Time evolution of the radial profile of the projected velocity
anisotropy using a line-of-sight parallel to the $\boldsymbol{\omega}$ of
stars within 1.5 $r_{\rm h}$ to remove as much of the inner tangential bias as
possible. $R$ and $R_{\rm h}$ refer to the projected radius and projected
half-mass radius. Figure 11: Time evolution of the rotational velocity profile
about the $z$-axis in all systems. $R$ and $R_{\rm h}$ refer to the projected
radius and projected half-mass radius, respectively. The rotational velocities
are also normalized to the velocity of the rotating coordinate system at the
initial Jacobi radius. Figure 12: Left panel: time evolution of the Cartesian
components of $\boldsymbol{\omega}$ calculated on the whole cluster in model
‘Cold th45’ over a timespan of several initial half-mass relaxation times,
$t_{\rm rh,i}$. Middle and right panels: time evolution of the $\phi$ and
$\theta$ components of global $\boldsymbol{\omega}$. Figure 13: Left panel:
time evolution of the radial profiles of $\phi_{\omega}$ over a timespan of
about 1$t_{\rm rh,i}$ (the legend of each panel is time normalized to $t_{\rm
rh,i}$). Right panel: time evolution of the radial profiles of
$\theta_{\omega}$ over a timespan of about 1$t_{\rm rh,i}$, but also includes
profiles at 2, 4 and 8$t_{\rm rh,i}$ to show that the radial variation remains
over several relaxation times.

We begin with an overview and visualization of the violent relaxation phase of
our systems. Our systems are all initially cold homogeneous spheres of stars
spinning in solid-body rotation. They then collapse from lack of pressure
support, reaching a maximum density, shortly after 1$t_{\rm ff}$. The clusters
then expand and eventually settle into an equilibrium configuration. This
process causes many stars to escape the cluster, and leaves a configuration
characterized by a inner dense core and sparse outer halo; many of the halo
stars are energetically unbound but still within the Jacobi radius of the
cluster. By conservation of angular momentum, when the cluster collapses, it
rotates faster; when the cluster expands, the rotation speed of the stars that
move to the outer regions slows down, creating a differential rotation within
the cluster (see next section and also Vesperini et al. 2014; Tiongco et al.
2017; Livernois et al. 2021). Each clusters’ final shape (specifically near
the core to around the half-mass radius, $r_{\rm h}$, excluding the sparse
halo) at the end of the violent relaxation process is triaxial (see also Sylos
Labini et al. 2015); the elongation of the major axis is primarily influenced
by the initial amount of rotation, with faster rotating initial conditions
producing more elongated systems.

For the set of models starting with $\theta_{\rm\omega,i}$ equal to 45
degrees, Fig. 1 shows the kinematics and structure of the clusters after they
have settled into an equilibrium configuration at the end of the simulation.
Specifically, this figure shows three projected number density images along
the $x$-, $y$-, and $z$-axes, and also three images of the line-of-sight
(hereafter, los) velocities in those same projections with isodensity contours
overlaid.

Fig. 2 shows these same panels for simulation Cold th90, but zoomed out to
emphasize the kinematical differences between the inner and the outer regions;
these differences appear in all of the simulations111Animations showing the
evolution of these projected number density and line-of-sight velocity images
(zoomed in and out) are available online at https://mtiongco.github.io/vrrot-
anims.

The images selected here for illustration, along with the animations, offer a
representation of the complexity that emerges from our generic initial
conditions, particularly if a cluster is ‘observed’ along an arbitrary
direction of the line-of-sight and evolutionary stage, and we do not assume
knowledge of the initial orientation of the rotation axis. Especially in the
cases where clusters start with their $\theta_{\rm\omega,i}$ equal to 45 or
135 degrees, the cluster’s triaxiality and the rotation axis not being aligned
with any of the $x$, $y$, or $z$ directions, can lead to projected
configurations characterised by a variety of aspect ratios and apparent axes
about which the cluster rotates. The simulations also show that the los
velocities are not all aligned along one rotation axis, with a larger
deviation from the initial rotation axis direction in the outer regions, where
the effects of the tidal field are stronger. Though this may not be obvious in
all of the los velocity maps, we will further illustrate this point using
another diagnostic in the next section.

###  3.2 Three-dimensional Properties

The complexity of the kinematics that develops during the violent relaxation
phase in our clusters is reflected in the evolution of the radial profiles of
$\phi_{\omega}$ and $\theta_{\omega}$, as shown in Figs 3 and 4. In Fig. 5, we
also quantify the amount of rotation in each system at different times. In all
these figures, at each time the cluster is divided into spherical shells in
which the angular velocity vector, $\boldsymbol{\omega}$ is calculated by
solving the equation $\textbf{{L}}=\textbf{{I}}\boldsymbol{\omega}$, where L
is the angular momentum vector and I is the moment of inertia tensor. In Fig.
5, we illustrate, as a function of the spherical radius, the rotational
velocity which results from the product of the magnitude of the angular
velocity vector by the radius of the spherical shell. These profiles are
intended to represent the ordered motions averaged under the assumption of
spherical symmetry.

These figures show that each system has differential rotation (Fig. 5), but
between each spherical shell the direction of the rotation axis varies in both
$\phi_{\omega}$ and $\theta_{\omega}$, and these radial variations also change
throughout the simulation. The radial variation in $\phi_{\omega}$
distinguishes these results from the simulations in Tiongco et al. (2018,
which did not include an analysis of the violent relaxation phase), where no
radial variation in $\phi_{\omega}$ was observed. Such radial variation in
$\phi_{\omega}$ introduces a twist of the alignment of the rotation axis
throughout the spatial extension of the system. It can also be observed in
these figures that the evolution in the radial variations is most significant
in the outermost regions of the cluster, sweeping over a large range of
angles, while in the inner regions, the rotation axis directions stay closer
to their original values.

At the end of the violent relaxation process, $\boldsymbol{\omega}$ in the
outermost regions of each cluster points closer to the direction of the
orbital angular velocity vector than it did initially. Such a feature is more
evident in clusters that begin with $\theta_{\rm\omega,i}=45$, while in
clusters that begin with $\theta_{\rm\omega,i}=135$ the outermost regions (up
to the Jacobi radius) tend to point into the orbital plane
($\theta_{\omega}\approx 90$) by the end of the simulation.

Next, we proceed in the presentation of our results with a discussion of the
structure of the systems’ equilibrium properties emerging at the end of the
violent relaxation process. After the system relaxes into an approximately
steady state, where structural properties such as the half-mass radius,
$r_{\rm h}$, no longer change significantly on the free-fall timescale, we
calculate the ratios of the principal semi-axes $a$, $b$, and $c$ (in order
from longest to shortest) in ellipsoidal shells at different radii throughout
the cluster, where the method to calculate the shape of each shell is detailed
in Zemp et al. (2011, method S1). Briefly, the principal semi-axes are found
by diagonalizing the shape tensor,

$\textbf{{S}}_{ij}=\frac{\sum_{k}m_{k}(\textbf{{r}}_{k})_{i}(\textbf{{r}}_{k})_{j}}{\sum_{k}m_{k}},$
(3)

of a shell of particles, and the square roots of the eigenvalues are
proportional to the principal semi-axes. This process is iterated until the
ratios of the principal semi-axes converge to values that do not change by
more than 0.1 percent between iterations.

In Fig. 6, we show the axial ratios as a function of radius throughout each
cluster, averaged over several profiles between 17.1 and 17.9 $t_{\rm ff}$. In
general, the clusters have a triaxial shape, with differences depending on the
initial virial ratios. The ‘Cold’ series of clusters have a deeper collapse,
and thus faster rotation and flatter shape. These clusters are also flatter in
the center and rounder in the outer regions, where the transition region is
about 1–3 $r_{\rm h}$. The ‘Even’ series of clusters follows a similar trend
but are overall less flat than the ‘Cold’ series. Finally, the models in the
the ‘Slow’ series are those with morphology closest to spherical; this is a
consequence of the fact that these models are those with the lowest initial
rotational energy among those considered in this study.

Observational studies of nuclear star clusters and galaxies have shown
misalignments of structural principal axes and kinematic rotation axes (Franx
& Illingworth 1988; De Rijcke et al. 2004; Geha et al. 2005; Thomas et al.
2006; Feldmeier et al. 2014; Tsatsi et al. 2017; del Pino et al. 2021). The
method used to find the shape of ellipsoidal shells also gives the 3D
orientation of the principal axes via the eigenvectors. At many times
throughout each simulation, we find that the orientation of the semi-minor
axis, $c$, differs from the orientation of rotation axis,
$\boldsymbol{\omega}$. Fig. 7 shows an example of this for one of the models.

In the final part of this section, we discuss the measurement of velocity
dispersion anisotropy in the system. Velocity anistotropy is often measured by
the parameter,

$\beta=1-\frac{\sigma_{\theta}^{2}+\sigma_{\phi}^{2}}{2\sigma_{\rm r}^{2}},$
(4)

(see, e.g., Binney & Tremaine 2008), and the radial profile of $\beta$ is
typically calculated in spherical shells. We show the results of using this
method in Fig. 8 at the end of simulation ‘Slow th90’ in the top left panel.
For the adopted definition of $\beta$, $\beta=0$ corresponds to an isotropic
velocity distribution, while negative (positive) values of $\beta$ correspond
to a tangentially (radially) anisotropic velocity distribution. The profile of
$\beta$ in the mid-to-outer regions is consistent with previous numerical
studies of violent relaxation of star clusters in an external tidal field
(e.g., see Boily et al. 1999; Vesperini et al. 2014) in that radial anisotropy
has a maximum in between 2–4$r_{\rm h}$, and then, going further outwards, the
radial anisotropy decreases and the velocity distribution becomes almost
isotropic. The shape of the profile in the mid-to-outer regions can be
explained by the expansion of the cluster during the rebounce phase of violent
relaxation which creates radial orbits in the outer regions of the cluster,
while presence of the tidal field makes the outermost regions of the cluster
more isotropic due to the preferential escape of stars on radial orbits (see,
e.g., Vesperini et al. 2014; Tiongco et al. 2016).

In our models, we find that the inner regions are characterized by tangential
anisotropy (for similar considerations in isolated models, see also Trenti et
al. 2005). However, this may be in part due to the inappropriate geometry
associated with the choice of spherical binning for rotating systems: a
spherical bin would include stars that have most of their velocity along the
azimuthal direction, but also stars that do not, creating an artificially
higher tangential velocity dispersion. We attempt to reconcile this by
calculating $\beta$ in cylindrical shells aligning them along a few different
directions. This particular model shown in Fig. 8, ‘Slow th90,’ was chosen
because its initial $\boldsymbol{\omega}$ points along the positive $x$-axis.
Since the most ideal subtraction of the mean velocity occurs when the
cylindrical bins are aligned with the rotation axis of the cluster, the top
right panel of Fig. 8 shows that using a los along the $x$-axis in the
velocity anisotropy calculation has a reduced inner tangential anisotropy
compared to setting the los along the $y$\- or the $z$\- axes. In the last two
panels of the second row of Fig. 8, we use lines of sight parallel to
$\boldsymbol{\omega}$ of stars within 1.5$r_{\rm h}$ and outside of 1.5$r_{\rm
h}$ (the dividing point was determined from slightly further out from where
the rotational radial profile approximately peaks). We find that using this
binning and los reduces the tangential velocity anisotropy in $\beta$ in the
innermost regions, although finding the most effective alignment of the
cylindrical bin is complicated by the radial variation of the alignment of the
rotation axis.

We also show a similar analysis in a second set of panels in Fig. 8 for
simulation ‘Slow th135’ where the initial $\boldsymbol{\omega}$ is not aligned
along any of the $x$-, $y$-, or $z$-axes. The largest reduction of the
tangential anisotropy in $\beta$ is obtained by adopting a los aligned with
the rotation axis of the cluster.

###  3.3 Projected Properties

In order to connect our results to observational studies, we now present the
results of our analysis of the cluster kinematic properties determined using
projected quantities. We first focus our attention on the velocity dispersion
anisotropy: Fig. 8 shows (dashed lines) the projected anisotropy parameter,
$1-\sigma_{\rm T}/\sigma_{\rm R}$; $\sigma_{\rm T}$ and $\sigma_{\rm R}$ are,
respectively, the tangential and radial velocity dispersion measured on the
projection plane perpendiclar to the line of sight. This parameter is measured
in the same cylindrical bins used for $\beta$ in the panels after the top left
one, and is effectively a different measure since it disregards the
contribution of line-of-sight velocities (i.e., this measure is similar to
that obtained with the velocities from proper motions). We find that using the
projected anisotropy parameter further reduces the innermost tangential
anisotropy bias, eliminating it in the cases where the los is aligned with the
rotation axis of the cluster.

Next, in Fig. 10, we show all of the models together plotting the evolution of
the projected anisotropy parameter, $1-\sigma_{\rm T}/\sigma_{\rm R}$. In Fig.
10, the line of sight is aligned with the rotation axis of the inner regions
(corresponding to ‘los:r<1.5rh’ in Fig. 8). The anisotropy profile shows that
the system is more radially anisotropic, particularly in the outer regions
during the violent relaxation process, and this reflects the expansion of the
cluster after the initial collapse phase. As more stars escape from the
cluster when the system is settling, the system loses some radial anisotropy,
but the evolution slows down at later times in the simulation.

Finally, in Fig. 11, we plot the evolution of the projected rotation curves
for all models from the same line of sight, chosen to be the $z$-axis. These
radial profiles are calculated from tangential velocities in the plane
perpendicular to the chosen line of sight (equivalent to rotation curves
calculated from proper motions in observations), making this figure different
from Fig. 5 where $\omega r$ was calculated from 3d velocities, and the
rotation axis direction was not constant with radius. Fig. 11 shows that a
variety of rotation curve shapes can be obtained looking at the cluster along
different lines of sight: the ‘th45’ and the ’th135’ series differ in the
direction of rotation but both show a typical rotation curve characterized by
rotational velocities increasing in the inner regions and decreasing in the
outer regions; the ‘th90’ series, on the other hand, shows very small
velocities due to line of sight being almost completely misaligned from the
rotation axis. These results highlight the importance of obtaining the full
phase space information to accurately measure, for example, the angular
momentum in a system.

###  3.4 Long-term Evolution

In this section, we present the results of a simulation following the long-
term evolution driven by two-body relaxation of our model, ‘Cold th45’. The
goal here is to explore how the various complexities emerging at the end of
the violent relaxation phase are affected by the cluster’s long-term evolution
and compare these results with those found in Tiongco et al. (2018).

One of the major results of Tiongco et al. (2018) was that the rotation axis
of the cluster precesses due to the torque of the host galaxy. In Fig. 12, we
show the evolution of the Cartesian components of the global
$\boldsymbol{\omega}$ for ‘Cold th45’, where the precession about the $z$-axis
is indicated by the oscillation of $x$\- and $y$\- components, and along with
the almost linear change with time in $\phi_{\omega}$. Also in line with our
previous results is the gradual loss of angular velocity (rotation) over
several relaxation times, as indicated by the magnitude of $\omega_{\rm z}$
and amplitude of the oscillations in $\omega_{\rm x}$ and $\omega_{\rm y}$
decreasing. The third panel of Fig. 12 shows the evolution of
$\theta_{\omega}$ decreasing a few degrees due to the loss of intrinsic
angular momentum and the system moving towards a tidally locked configuration
in which the rotation of the cluster would be predominantly about the
$z$-axis.

We also examine the long-term evolution of the profiles $\phi_{\omega}$ and
$\theta_{\omega}$; as discussed in the previous sections, at the end of the
violent relaxation phase the cluster is characterized by a radial variation in
$\phi_{\omega}$ and $\theta_{\omega}$. The left panel of Fig. 13 shows that
the radial variation in $\phi_{\omega}$ is gradually erased over the two-body
relaxation timescale (the profiles for $\phi_{\omega}$ move down along the
plot due to the precession of the cluster’s rotation axis as mentioned
previously). On the other hand, the other panel shows the evolution of
$\theta_{\omega}$, where the radial variation remains over several two-body
relaxation timescales.

##  4 Conclusions

We have explored via $N$-body simulations the early dynamical evolution of
rotating star clusters undergoing violent relaxation under the effects of the
tidal field of their host galaxy. By allowing the clusters to initially have a
rotation axis orientation not aligned with their orbital angular velocity
vector, we have found that the equilibrium configurations reached at the end
of the violent relaxation phase are characterized by a variety of complex
structural and kinematic features.

We have explored nine different initial conditions varying in initial virial
ratio in random and rotational kinetic energy: the ‘Even’ series having the
same virial ratio in random and rotational kinetic energy, the ‘Cold’ series
having a lower random virial ratio, and the ‘Slow’ series having a lower
rotational virial ratio. We have also explored the evolution of systems with
different values of the initial angle between the cluster’s rotation axis and
the cluster’s orbital angular velocity vector (45 degrees, 90 degrees, and 135
degrees).

The violent relaxation process leaves the cluster with a differential rotation
curve characterized by a rotational velocity rising in the inner regions,
peaking in the intermediate regions, and then falling in the outer parts of
the cluster.

Most notably, we find that the clusters end up with a radial variation in
rotation axis direction along both the polar angle $\theta$ and the azimuthal
angle $\phi$ directions. During the violent relaxation process, the rotation
axis directions oscillate along the $\theta$ and $\phi$ directions, with a
stronger amplitude in the outermost regions. The radial variation in the
rotation axis directions is found also in the final equilibrium configuration
at the end of the violent relaxation process, with the outermost regions
pointing closer to the direction of the orbital angular velocity vector.

Our clusters develop a triaxial morphology with the models in the ‘Cold’
series being those showing the strongest flattening; this is a consequence of
the fact that the models of this series are those that undergo the deepest
collapse and emerge from the violent relaxation with the most rapid rotation.
We find a radial variation in the morphological minor axis position angle, but
this position angle does not necessarily align with the rotation axis.

We also explored how the complex features of these clusters appear in
projection as they do in observational studies. From different perspectives, a
radial velocity field image combined with photometric isophotes may look like
the cluster is rotating about its major, minor, or another axis altogether.
The measurement of velocity dispersion anisotropy may have a tangential bias
near the centre of the cluster due to a suboptimal orientation of the line of
sight (i.e., a los that is not close to the rotation axis of cluster) that can
lead to anisotropy measurements contaminated by ordered rotational motion. The
tangential bias may be reduced when the projected anisotropy parameter is
measured on the plane normal to the mean position of the rotation axis;
however, the tangential bias may not be completely eliminated due to the
rotation axis positional variation throughout the cluster.

Finally, we have followed the evolution of one of our clusters over a
timescale of several half-mass relaxation times, and found that the radial
variation of the rotation axis position along the azimuthal direction is
erased over a relaxation time. Thus, observing a radial variation of the
position angle of the rotation axis along the azimuthal direction of a star
cluster may be a kinematical signature of a dynamically young cluster. On the
other hand, the radial variation of the rotation axis position along the polar
$\theta$ direction is not erased during the cluster’s long-term evolution.
This is a consequence of the evolution towards a tidally locked rotation that
drives the evolution of the cluster towards rotation about the $z$-axis (i.e.
parallel to the orbital angular velocity vector of the cluster).

Current morphological and kinematical profiles of Galactic globular clusters
do not yet have the ability to fully identify the subtle radial variations we
have highlighted here, mostly due to the astrometric limitations in the
reconstruction of the 3D velocity space of objects characterised by such
crowded fields. Yet, we wish to note that a detailed study of the Sagittarius
dwarf galaxy by del Pino et al. (2021) has recently brought to light the first
evidence of precession and nutation motions imparted by the Milky Way tidal
field onto a satellite stellar system. We are therefore hopeful that
forthcoming Gaia data releases will enable this kind of refined investigations
also for local globular clusters. More generally, the complex structural and
kinematic features revealed by our study provide new elements for the
interpretation of observations of star clusters and can shed light on their
formation and dynamical history. In future studies, we will extend the
investigation presented in this paper to consider more realistic galactic
potentials and the effects of a time-dependent external potential for clusters
on non-circular orbits.

## Acknowledgements

MT thanks the Eccentric Dynamics research group at CU Boulder for feedback on
the results and presentation. This research was supported in part by Lilly
Endowment, Inc., through its support for the Indiana University Pervasive
Technology Institute, and in part by the Indiana METACyt Initiative. The
Indiana METACyt Initiative at IU is also supported in part by Lilly Endowment,
Inc. ALV acknowledges support from a UKRI Future Leaders Fellowship
(MR/S018859/1).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations. Cambridge University Press
  * Aarseth et al. (1988) Aarseth S. J., Lin D. N. C., Papaloizou J. C. B., 1988, ApJ, 324, 288
  * Aguilar & Merritt (1990) Aguilar L. A., Merritt D., 1990, ApJ, 354, 33
  * Armstrong et al. (2020) Armstrong J. J., Wright N. J., Jeffries R. D., Jackson R. J., 2020, MNRAS, 494, 4794
  * Ballone et al. (2021) Ballone A., Torniamenti S., Mapelli M., Di Carlo U. N., Spera M., Rastello S., Gaspari N., Iorio G., 2021, MNRAS, 501, 2920
  * Bellini et al. (2017) Bellini A., Bianchini P., Varri A. L., Anderson J., Piotto G., van der Marel R. P., Vesperini E., Watkins L. L., 2017, ApJ, 844, 167
  * Bianchini et al. (2013) Bianchini P., Varri A. L., Bertin G., Zocchi A., 2013, ApJ, 772, 67
  * Bianchini et al. (2018) Bianchini P., van der Marel R. P., del Pino A., Watkins L. L., Bellini A., Fardal M. A., Libralato M., Sills A., 2018, MNRAS, 481, 2125
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Boberg et al. (2017) Boberg O. M., Vesperini E., Friel E. D., Tiongco M. A., Varri A. L., 2017, ApJ, 841, 114
  * Boily et al. (1999) Boily C. M., Clarke C. J., Murray S. D., 1999, MNRAS, 302, 399
  * Chen et al. (2021) Chen Y., Li H., Vogelsberger M., 2021, MNRAS, 502, 6157
  * Dalessandro et al. (2021a) Dalessandro E., Raso S., Kamann S., Bellazzini M., Vesperini E., Bellini A., Beccari G., 2021a, MNRAS, 506, 813
  * Dalessandro et al. (2021b) Dalessandro E., et al., 2021b, ApJ, 909, 90
  * De Rijcke et al. (2004) De Rijcke S., Dejonghe H., Zeilinger W. W., Hau G. K. T., 2004, A&A, 426, 53
  * Einsel & Spurzem (1999) Einsel C., Spurzem R., 1999, MNRAS, 302, 81
  * Ernst et al. (2007) Ernst A., Glaschke P., Fiestas J., Just A., Spurzem R., 2007, MNRAS, 377, 465
  * Fabricius et al. (2014) Fabricius M. H., et al., 2014, ApJ, 787, L26
  * Feldmeier et al. (2014) Feldmeier A., et al., 2014, A&A, 570, A2
  * Ferraro et al. (2018) Ferraro F. R., et al., 2018, ApJ, 860, 50
  * Franx & Illingworth (1988) Franx M., Illingworth G. D., 1988, ApJ, 327, L55
  * Gebhardt et al. (2000) Gebhardt K., Pryor C., O’Connell R. D., Williams T. B., Hesser J. E., 2000, AJ, 119, 1268
  * Geha et al. (2005) Geha M., Guhathakurta P., van der Marel R. P., 2005, AJ, 129, 2617
  * Gott (1973) Gott Richard J. I., 1973, ApJ, 186, 481
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational Million-Body Problem. Cambridge University Press
  * Heisler & Tremaine (1986) Heisler J., Tremaine S., 1986, Icarus, 65, 13
  * Hohl & Zang (1979) Hohl F., Zang T. A., 1979, AJ, 84, 585
  * Hong et al. (2013) Hong J., Kim E., Lee H. M., Spurzem R., 2013, MNRAS, 430, 2960
  * Jerabkova et al. (2021) Jerabkova T., Boffin H. M. J., Beccari G., de Marchi G., de Bruijne J. H. J., Prusti T., 2021, A&A, 647, A137
  * Kamann et al. (2018) Kamann S., et al., 2018, MNRAS, 473, 5591
  * Kim et al. (2002) Kim E., Einsel C., Lee H. M., Spurzem R., Lee M. G., 2002, MNRAS, 334, 310
  * Kim et al. (2004) Kim E., Lee H. M., Spurzem R., 2004, MNRAS, 351, 220
  * Kim et al. (2008) Kim E., Yoon I., Lee H. M., Spurzem R., 2008, MNRAS, 383, 2
  * Kuhn et al. (2019) Kuhn M. A., Hillenbrand L. A., Sills A., Feigelson E. D., Getman K. V., 2019, ApJ, 870, 32
  * Lahén et al. (2020) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2020, ApJ, 904, 71
  * Lanzoni et al. (2018a) Lanzoni B., et al., 2018a, ApJ, 861, 16
  * Lanzoni et al. (2018b) Lanzoni B., et al., 2018b, ApJ, 865, 11
  * Lim et al. (2020) Lim B., Hong J., Yun H.-S., Hwang N., Kim J. S., Lee J.-E., Park B.-G., Park S., 2020, ApJ, 899, 121
  * Livernois et al. (2021) Livernois A., Vesperini E., Tiongco M., Varri A. L., Dalessandro E., 2021, MNRAS,
  * Lynden-Bell (1967) Lynden-Bell D., 1967, MNRAS, 136, 101
  * Mapelli (2017) Mapelli M., 2017, MNRAS, 467, 3255
  * Nitadori & Aarseth (2012) Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545
  * Sollima et al. (2019) Sollima A., Baumgardt H., Hilker M., 2019, MNRAS, 485, 1460
  * Sylos Labini et al. (2015) Sylos Labini F., Benhaiem D., Joyce M., 2015, MNRAS, 449, 4458
  * Thomas et al. (2006) Thomas D., Brimioulle F., Bender R., Hopp U., Greggio L., Maraston C., Saglia R. P., 2006, A&A, 445, L19
  * Tiongco et al. (2016) Tiongco M. A., Vesperini E., Varri A. L., 2016, MNRAS, 455, 3693
  * Tiongco et al. (2017) Tiongco M. A., Vesperini E., Varri A. L., 2017, MNRAS, 469, 683
  * Tiongco et al. (2018) Tiongco M. A., Vesperini E., Varri A. L., 2018, MNRAS, 475, L86
  * Tiongco et al. (2019) Tiongco M. A., Vesperini E., Varri A. L., 2019, MNRAS, 487, 5535
  * Trenti et al. (2005) Trenti M., Bertin G., van Albada T. S., 2005, A&A, 433, 57
  * Tsatsi et al. (2017) Tsatsi A., Mastrobuono-Battisti A., van de Ven G., Perets H. B., Bianchini P., Neumayer N., 2017, MNRAS, 464, 3720
  * Usher et al. (2021) Usher C., Kamann S., Gieles M., Hénault-Brunet V., Dalessandro E., Balbinot E., Sollima A., 2021, MNRAS,
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 489, 623
  * Vesperini et al. (2014) Vesperini E., Varri A. L., McMillan S. L. W., Zepf S. E., 2014, MNRAS, 443, L79
  * Vitral & Mamon (2021) Vitral E., Mamon G. A., 2021, A&A, 646, A63
  * Wan et al. (2021) Wan Z., et al., 2021, MNRAS, 502, 4513
  * Zemp et al. (2011) Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2011, ApJS, 197, 30
  * del Pino et al. (2021) del Pino A., Fardal M. A., van der Marel R. P., Łokas E. L., Mateu C., Sohn S. T., 2021, ApJ, 908, 244
  * van Albada (1982) van Albada T. S., 1982, MNRAS, 201, 939


L PREPRINT - accepted by IEEE International Symposium on Hardware Oriented
Security and Trust (HOST) 20223.5cm,1cm1.001.8

# pHGen: A pH-Based Key Generation Mechanism Using ISFETs

Elmira Moussavi1, Dominik Sisejkovic1, Fabian Brings2, Daniyar Kizatov1,
Animesh Singh2,
Xuan Thang Vu2, Rainer Leupers1, Sven Ingebrandt2, Vivek Pachauri2, and Farhad
Merchant1 {moussavi, sisejkovic, kizatov, leupers<EMAIL_ADDRESS>aachen.de
{brings, singh, vu, ingebrandt<EMAIL_ADDRESS>1Institute for
Communication Technologies and Embedded Systems, 2Institute of Materials in
Electrical Engineering 1  RWTH Aachen University, Germany

###### Abstract

Digital keys are a fundamental component of many hardware- and software-based
security mechanisms. However, digital keys are limited to binary values and
easily exploitable when stored in standard memories. In this paper, based on
emerging technologies, we introduce pHGen, a potential-of-hydrogen (pH)-based
key generation mechanism that leverages chemical reactions in the form of a
potential change in ion-sensitive field-effect transistors (ISFETs). The
threshold voltage of ISFETs is manipulated corresponding to a known pH buffer
solution (key) in which the transistors are immersed. To read the chemical
information effectively via ISFETs, we designed a readout circuit for stable
operation and detection of voltage thresholds. To demonstrate the
applicability of the proposed key generation, we utilize pHGen for logic
locking—a hardware integrity protection scheme. The proposed key-generation
method breaks the limits of binary values and provides the first steps toward
the utilization of multi-valued voltage thresholds of ISFETs controlled by
chemical information. The pHGen approach is expected to be a turning point for
using more sophisticated bio-based analog keys for securing next-generation
electronics.

###### Index Terms:

Hardware security, emerging technologies, logic locking, key generation, ion-
sensitive field-effect transistor.

##  I Introduction

Security has recently become particularly important for hardware designs [1].
Especially, due to the high cost of manufacturing and maintaining nanoscale
integrated circuit (IC) foundries, a significant portion of chips are designed
and fabricated offshore [2]. This globalization of the IC design flow often
leaks the vital information of the chips to untrustworthy parties and creates
opportunities for intellectual property (IP) theft, overproduction, reverse
engineering, counterfeiting, and insertion of hardware Trojans [3]. Techniques
to ensure the security of hardware integrity through IC design and fabrication
have been one of the emerging areas of information technology in recent years,
where hardware security approaches utilize conventional digital keys to offer
protection from an attack. However, there remain several ways to attack, such
as exploitation of the digital data (e.g., secret keys in cryptographic
applications) enabled by the highly vulnerable storage.

In this work, we present pHGen, a pH-based key generation mechanism that
utilizes emerging transistor technologies to generate secure keys from
(bio)chemical information. The bio-based key generation offers multiple
advantages over standard, binary keys, including a multi-valued
representation, enhanced security in terms of physical key inspection, and a
path towards personalized, biochemical digital keys. The main contributions of
this paper are as follows (Fig. 1):

  * •

The implementation of pHGen; a pH-based key generation mechanism that
manipulates the voltage thresholds of ISFETs immersed into specific pH buffer
solutions to extract analog chemical keys.

  * •

The design and development of a stable readout circuit to read the intrinsic
voltage threshold shift of ISFETs for integration with field-programmable gate
arrays (FPGAs).

  * •

The demonstration of pHGen for the generation of secure analog keys in logic-
locked circuits.

Figure 1: The flow of pHGen.

##  II pHGen Implementation

###  II-A ISFET- an Emerging Technology

To perform the transition from digital to analog keys, we utilize hydrogen ion
(H+) concentration or pH values as they reflect chemical information.
Miniaturized ion-selective sensors such as micro and nanoISFETs can be used as
rapid analytical devices to detect small variations of free hydrogen and
hydroxyl ions[4]. In this work, we use microISFETs as a test-board for well-
established, reproducible, clean-room fabrication at wafer-scale, ion-
recognition, and circuit-integration possibilities. Sensitive to surface-
charge, ISFETs are utilized as chemical and biosensors to analyze solution
samples in contact with its gate [5]. In addition, ISFETs offer possibilities
for scalable integration using an unmodified CMOS process towards industrial
application. ISFETs operate similar to MOSFETs, whose gate is extended with a
passivation layer in contact with an external reference electrode (Ag/AgCl)
via a sample solution (Fig. 2) also known as electrochemical gate
configurations. The output current $I_{DS}$ of an n-type transistor is
expressed as follows:

Figure 2: (a) ISFET as a modified MOSFET technology and (b) equivalent
electrical model.

$I_{DS}=\mu_{n}c_{ox}\frac{W}{L}(V_{GS}-V_{TH(ISFET)})V_{DS}-\frac{1}{2}V_{DS}^{2}.$
(1)

In comparison to MOSFET, $V_{TH(MOSFET)}$ is replaced by $V_{TH(ISFET)}$ in
eq. 1, which is influenced by surface charge changes as a function of the pH
value of the solution in contact with the ISFET gate. Thus, ISFETs measure the
solution pH in the form of threshold voltage change ($V_{TH}$) of a
corresponding floating-gate MOSFET device. The electrochemical interface
($V_{CHEMICAL}$) consists of a reference built-in potential ($V_{ref}$)
together with the double-layer capacitance known as Gouy–Chapman ($C_{Gouy}$)
and Helmholtz ($C_{Helm}$) plane[6]. $C_{Gouy}$ and $C_{Helm}$ define the pH
potential drop ($V_{pH}$) at the MOSFET gate and the $C_{pass}$ is the
passivation capacitance. The electrochemically induced voltage in the
$V_{TH(ISFET)}$ can be defined as follows:

$V_{TH(ISFET)}=V_{TH(MOSFET)}-V_{CHEMICAL},$ (2)

where $V_{CHEMICAL}$ is defined as:

$V_{CHEMICAL}=E_{i}+\frac{RT}{n_{i}F}ln(\alpha_{i}).$ (3)

In eq. 3, $E_{i}$ is the standard electrode potential, R: gas constant , T:
absolute temperature (K), F: Faraday constant, $n_{i}$: charge of ion $i$, and
$\alpha_{i}$: activity of ions. The only variable in eq. 3 is $\alpha_{i}$.
Thereby, to monitor the ion activity it is sufficient to keep $V_{DS}$ and
$I_{DS}$ constant. In this case, $V_{GS}$ is proportional to the natural
logarithm of ion activity. Hence, the ISFET needs to be connected to a readout
circuit to provide this condition.

###  II-B Readout Circuit Design

A constant voltage constant current (CVCC) circuit is implemented as a readout
circuit (ROC) to monitor the ion activity in response to an electrolyte
solution that is in contact with the ISFET gate. For robustness, the sensor
operates in the linear region, and the gate voltage is connected to the
ground. Thus, the source voltage serves as the output signal proportional to
the internal ISFET’s threshold voltage ($V_{TH(ISFET}$). Fig. 3 shows the
designed ROC. The voltage across ISFET’s drain and source ($V_{DS}$) is
constant and equal to $R1*I$. In addition, the source current is half to that
of the sink current, which according to Kirchhoff’s current law, keeps
$I_{DS}$ constant and is equivalent to $I$. A current source circuit is used
to generate the current sources. This circuit includes an instrumentation
amplifier (In-amp) and an operational amplifier. Therefore, the CVCC circuit
ensures a constant operating current over ISFET’s drain-source terminals
($I_{DS}=I=$ $50\text{\,}\mathrm{\SIUnitSymbolMicro A}$). Consequently, R1
constantly maintains $V_{DS}=R1*I=$ $10\text{\,}\mathrm{k\SIUnitSymbolOhm}$
$*$ $50\text{\,}\mathrm{\SIUnitSymbolMicro A}$ $=$ $0.5\text{\,}\mathrm{V}$.
In addition to current source circuits and the R1 resistance, the CVCC circuit
comprises one operational amplifier (U1) connected to ISFET’s drain providing
feedback to the source and ensuring the functionality of the CVCC circuit. A
buffer stage is employed to improve driving capability. Drawbacks related to
temperature dependency [7] can be addressed by using a temperature sensor to
compensate for the temperature effect.

Figure 3: Readout circuit designed for ISFET integration.

##  III Demonstration

(a) ISFET pH sensitivity with (Ta2O5) gate at the temperature of 30°C.

(b)  Reading the pH response and temperature information from the
(bio)chemical-system.

Figure 4: Results for (a) source voltage changes vs. pH solutions, and (b) gate source voltage change over the time at constant drain source voltage of 0.5V.  TABLE I:  MSFET 3330-2 pH sensor dimensions. | Width | Height | Length
---|---|---|---
Chip dimensions | $1.2\text{\,}\mathrm{mm}$ | $0.3\text{\,}\mathrm{mm}$ | $3\text{\,}\mathrm{mm}$
Packaged sensor | $5\text{\,}\mathrm{mm}$ | $1\text{\,}\mathrm{mm}$ | $50\text{\,}\mathrm{mm}$

Key-driven countermeasures, such as logic locking, have been proposed to
protect the integrity of hardware designs through the IC design and
fabrication flow [8, 9, 10]. Logic locking inserts a locking mechanism that
ensures a correct circuit behavior only upon the application of a correct
activation key. These inserted key-dependent gates are referred to as key
gates. As the key is loaded post-fabrication by the IP owner, the design
remains concealed while in the hands of untrusted third parties such as the
external design house and foundry. To demonstrate the applicability of pHGen,
we present a logic-locked circuit that uses the pH solution as the source of
the input key to drive an XOR key gate. The demonstration was implemented with
an MSFET 3330-2 ISFET, which is a commercial ultra-miniature pH ISFET sensor
from MICROSENS S.A [11]. It consists of silicon- and polysilicon-based
materials and a Ta2O5 gate as a pH-sensitive metal with 4′′ planar CMOS
process technology. The ISFET sensor dimensions are listed in Table I. The
voltage threshold shift per unit change in pH is defined as a pH sensor
sensitivity [12]. With this sensor, we achieved a sensitivity of the
conventional Nernst limit of $\Delta V=$
$59.2\text{\,}\mathrm{mV}$$/$$\mathrm{p}\mathrm{H}$ at a pH range of 4 to 9
and at the temperature of 30°C (Fig. 4(a)).

Figure 5: Locked circuit schematic and implementation on the FPGA board with
an FMC connector.

The output voltage (the source voltage of the ISFET ($V_{S}$)) changes
proportionally to the pH value (pH 4 to 9), as presented in Fig. 4(a). The pH
and temperature response of the ISFET over time are presented in Fig. 4(b) for
two pH (4 and 7) and temperature levels (30°C and 40°C). However, more
sensitivity can be attained by enlarging the sensing gate of an extended gate
ISFET [13] or double gate ISFET design (front and back oxide capacitance)
[14]. To measure the potential developed in the ion-sensitive gate, an ISFET
together with a reference electrode (Ag/AgCl) must be immersed into the
solution. For the solution, a phosphate buffer is prepared to adjust the pH
value well between the pH range of 5.8 to 8.0, the example buffer that we used
during the measurements is pH 5.8 with a voltage of $2.07\text{\,}\mathrm{V}$
(Fig. 4(a)). The reason to stay in the pH range of 5.8 to 8.0 is that many
macromolecular properties are pH-dependent and stable within this range [15],
thus enabling stable bio-based keys in the next stage. Table II summarizes the
specifications of the designed readout circuit for $\mathrm{p}\mathrm{H}$
$=5.8$. A microfluidic system can be used to control the flow of the pH
solution to drive more logic gates.

Figure 6: Demonstration of logic locking with pHGen.

To demonstrate the integration of pHGen with logic locking, a locked circuit
was implemented by inserting a key gate (KG) in a selected gate-level netlist,
as illustrated in Fig. 5. The pH-based key for the key gate is provided by a
pH buffer solution ($\mathrm{p}\mathrm{H}$ $=5.8$). The gate of the ISFET was
immersed into the buffer solution to read ion activity in the form of a
voltage shift. Thus, the CVCC circuit reads the amount of voltage drift and
buffers it to isolate the driver from the load. A voltage comparator
constantly compares the voltage with the known pH voltage value. Because the
IP owner designs the circuit and knows the correct voltage change from the pH-
based key. To have a functional circuit, the output voltage obtained from the
pH solution needs to be equal to the same designed voltage value. This voltage
sets the value of the input pin of the key gate (XOR KG) via an FPGA Mezzanine
Card (FMC) connector. The locked circuit was implemented on a ZedBoard
Zynq-7000 FPGA board connected to an FMC. Four switches (SW 0-3) are set as
input signals to two XOR gates and combined with a NOR gate to implement a
pHGen logic locked circuit. An LED is connected to the output pin of a NOR
gate to indicate the output signal. The FMC connector is used to bring
modularity to change the I/O configuration of the design without FPGA
intervention. Hence, a quad 2-input XOR gate IC is placed on the FMC
connector. One of the XOR gates is used to establish a connection between the
output voltage of the readout circuit and the input pin of the key gate. The
FPGA-based test setup is shown in Fig. 6. The pH-driven key gate is added to
the netlist, through the readout circuit, via the FMC connector. To
demonstrate the key generation, the ISFET together with the reference
electrode are immersed in the specific pH solution (secret key) to create a
particular voltage shift and activate the XOR gate with the dedicated key
value. The switches (SW 0-3) are considered as input signals and the results
for all input patterns are shown in Table III. The obtained results state that
the pH-key generation is feasible, CMOS-compatible, and can produce the same
behavior as a conventional binary key in the form of voltage shift. Due to the
robustness of the CVCC readout circuit, the ISFET operates in the linear
region at a fixed drain-source voltage, independent of temperature and the
power supply variation. Under these conditions, the pH solution’s ion activity
was monitored, and the output voltage was observed as a change in the ISFET
threshold voltage. With the presented use case, we demonstrated the
feasibility of using chemical information as a source to secure analog keys.

TABLE II:  Readout circuit characteristics at 24°C, $\mathrm{p}\mathrm{H}$: 5.8. | Typ
---|---
Supply voltage ($V_{CC}$) | $0.5\text{\,}\mathrm{V}$
Operating current ($I_{DS}$) | $50\text{\,}\mu\mathrm{A}$
Internal drain voltage ($V_{D}$) | $0.7\text{\,}\mathrm{V}$
Internal source voltage ($V_{S}$) | $0.2\text{\,}\mathrm{V}$
Output voltage ($V_{OUT})$ | $2.07\text{\,}\mathrm{V}$
CVCC circuit power consumption | $50\text{\,}\mu\mathrm{W}$
TABLE III:  Results of the demonstrated locked circuit. Input patterns |  Output ($\mathrm{p}\mathrm{H}$ $=5.8$, Voltage = $2.07\text{\,}\mathrm{V}$, key =1) |  Output ($\mathrm{p}\mathrm{H}$ $\neq 5.8$, Voltage = $0\text{\,}\mathrm{V}$, key $\neq$ 1)
---|---|---
0000 | 0 | 1
0001 | 1 | 0
0010 | 1 | 0
0011 | 0 | 1
0100 | 0 | 0
0101 | 0 | 0
0110 | 0 | 0
0111 | 0 | 0
1000 | 0 | 0
1001 | 0 | 0
1010 | 0 | 0
1011 | 0 | 0
1100 | 0 | 1
1101 | 1 | 0
1110 | 1 | 0
1111 | 0 | 1

##  IV Next steps

This novel approach of utilizing chemical information for securing ICs is
expected to enable the application of next-in-line DNA-based nanoISFETs [16]
[17] for personalized, secure systems. The work presented here establishes a
basis for exploring the design of such chemically inspired secure ICs. In the
next steps, we plan to extend the work to enhance the security aspects of
logic locking in the following directions: (1) the design of multi-valued
logic, (2) the implementation of massively parallel ISFET arrays to enable
more key gates, and (3) the usage of nanoscale ISFETs to deploy super
sensitive characteristics combined with high density and CMOS-compatible
integration. (4) Finally, a microfluidic system will be developed to
automatically control the flow of the solution.

##  V Conclusion

In this work, we presented pHGen, a novel analog key generation approach that
can be utilized as an alternative to conventional binary keys. Emerging
technologies, such as the ISFETs, ameliorate the change towards the use of
analog keys. The key is generated by shifting the intrinsic voltage threshold
of the ISFET, which is observed based on the changes in the pH-soluble ion
activities. To demonstrate the proposed mechanism, we integrated pHGen with
logic locking—a prominent technique to protect hardware designs against
manipulations and piracy. According to the results, pH-based analog keys have
the ability to substitute conventional digital keys. Hence, pHGen provides a
stepping stone to using a chemical-based key with ISFETs.

## Acknowledgement

This work was partially funded by Deutsche Forschungsgemeinschaft (DFG –
German Research Foundation) under the priority programme SPP 2253.

## References

  * [1] J. A. Roy _et al._ , “Ending piracy of integrated circuits,” _Computer_ , vol. 43, no. 10, pp. 30–38, 2010.
  * [2] S. Ravi _et al._ , “Security in embedded systems: Design challenges,” _ACM Transactions on Embedded Computing Systems (TECS)_ , vol. 3, no. 3, pp. 461–491, 2004.
  * [3] M. Rostami _et al._ , “A primer on hardware security: Models, methods, and metrics,” _Proceedings of the IEEE_ , vol. 102, no. 8, 2014.
  * [4] R. Caenn _et al._ , “Chapter 3 - evaluating drilling fluid performance,” in _Composition and Properties of Drilling and Completion Fluids_ , seventh edition ed., R. Caenn _et al._ , Eds. Boston: Gulf Professional Publishing, 2017, pp. 55–91.
  * [5] Y. Hu _et al._ , “A robust ISFET pH-measuring front-end for chemical reaction monitoring,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 8, no. 2, pp. 177–185, 2014.
  * [6] R. Van Hal _et al._ , “A novel description of ISFET sensitivity with the buffer capacity and double-layer capacitance as key parameters,” _Sensors and Actuators B: Chemical_ , vol. 24, no. 1-3, pp. 201–205, 1995\.
  * [7] W. Gui-Hua _et al._ , “ISFET temperature characteristics,” _Sensors and actuators_ , vol. 11, no. 3, pp. 221–237, 1987.
  * [8] M. Yasin _et al._ , _Trustworthy Hardware Design: Combinational Logic Locking Techniques_. Springer International Publishing, 2020.
  * [9] D. Sisejkovic _et al._ , “A secure hardware-software solution based on RISC-V, logic locking and microkernel,” in _SCOPES_. New York, NY, USA: ACM, 2020, p. 62–65.
  * [10] S. Rai _et al._ , “Vertical IP protection of the next-generation devices: Quo vadis?” in _DATE_ , 2021, pp. 1905–1914.
  * [11] “Microsens SA - pH-ISFET sensor,” accessed 06-Jan-2022. [Online]. Available: http://microsens.ch/products/ISFET.htm
  * [12] P. Bergveld, “Thirty years of isfetology: What happened in the past 30 years and what may happen in the next 30 years,” _Sensors and Actuators B: Chemical_ , vol. 88, no. 1, pp. 1–20, 2003.
  * [13] K. B. Parizi _et al._ , “Exceeding nernst limit (59mV/pH): CMOS-based pH sensor for autonomous applications,” in _2012 International Electron Devices Meeting_ , 2012, pp. 24.7.1–24.7.4.
  * [14] S. Sanjay _et al._ , “Super-nernstian ion sensitive field-effect transistor exploiting charge screening in wse2/mos2 heterostructure,” _npj 2D Materials and Applications_ , vol. 5, no. 1, pp. 1–8, 2021.
  * [15] K. Talley _et al._ , “On the pH-optimum of activity and stability of proteins,” _Proteins: Structure, Function, and Bioinformatics_ , vol. 78, no. 12, pp. 2699–2706, 2010.
  * [16] V. Pachauri _et al._ , “Biologically sensitive field-effect transistors: from isfets to nanofets,” _Essays in biochemistry_ , vol. 60, no. 1, 2016\.
  * [17] C. Toumazou _et al._ , “A new era of semiconductor genetics using ion-sensitive field-effect transistors: the gene-sensitive integrated cell,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 372, no. 2012, p. 20130112, 2014.