Theraa TORK

#  $\psi(2S)$ production as a function of charged-particle multiplicity in pp
collisions at $\sqrt{s}$= 13 TeV and p–Pb collisions at $\snn$ = 8.16 TeV with
ALICE at the LHC

on behalf of the ALICE collaboration

###### Abstract

The self-normalised $\psip$ yields as a function of charged-particle
multiplicity in pp collisions at $\s$ = 13 TeV and p–Pb collisions at $\snn$ =
8.16 TeV, using the ALICE detector at the CERN LHC, are reported in these
proceedings. This measurement represents one of main tools to understand the
role of multiparton interactions and their influence on charmonium production
in small systems, i.e pp and pA collisions. The ratio of $\psip$-over-J/$\psi$
as a function of charged-particle multiplicity is also presented for both pp
and p–Pb collisions. The comparison to theoretical calculations is reported.

##  1 Introduction

Charmonium production is a complex mechanism that involves both perturbative
and non-perturbative QCD aspects. Due to their heavy masses, the charm quarks
are produced during the hard process of the collision, while the formation of
the bound state is a soft QCD process. Moreover, charmonia are useful tools to
investigate the medium formed in heavy-ion collisions (A-A), where a
deconfined state of quarks and gluons, known as quark-gluon plasma (QGP), is
expected. The study of charmonium production in smaller systems, i.e., pp and
pA collisions, is essential as a reference system where no QGP effects are
expected. However, several intriguing QGP-like behaviors were observed in
small systems, such as the non-zero elliptic flow of charged-particles through
their long-range angular correlation measurements in high multiplicity events
[1]. Multiparton interactions (MPI) are one of the main scenarios proposed to
explain these findings. In the MPI context, several parton-parton interactions
occur in the same hadron collision. The multiplicity dependence of charmonium
production is an indirect probe for MPI. It provides information on the
interplay between the soft and hard QCD process, as heavy quark production is
a hard process while charged particles are produced via soft QCD. In these
proceedings, the multiplicity dependence of the self-normalised $\psip$ yields
in pp and p–Pb collisions at $\s$ = 13 TeV and $\snn$ = 8.16 TeV,
respectively, are discussed. The $\psip$-over-J/$\psi$ ratio as a function of
charged-particle multiplicity is also shown.

##  2 Experimental setup

A detailed description of the ALICE detector and its performance can be found
in [2]. The analyzed data are recorded using the muon spectrometer, located in
the forward arm of the ALICE detector and covering the rapidity window of 2.5
< $\eta\rm{{}_{lab}}$ < 4.0. The muon spectrometer consists of the following
subdetectors: (i) a set of absorbers that are essential to reduce the
background contamination coming from the decay of hadrons and beam-gas
interaction, (ii) a dipole magnet that deflects particle trajectories to
measure their momentum and electric charge, (iii) five stations of muon
tracking chambers to reconstruct the particle trajectory and two stations of
muon trigger chambers which are essential for the muon identification. Other
detectors are also used in this analysis, such as the silicon pixel detector
(SPD) and the V0 detector. The SPD (|$\eta$| < 2) is important for the
reconstruction of the charged-particle tracks and the primary vertex. The V0
detector, which consists of two sets of scintillator arrays located at
backward -3.7 < $\eta$ <-1.7 and forward 2.8 < $\eta$ < 5.1 rapidity, serves
as a minimum bias trigger in this analysis.

##  3 Results

The multiplicity dependence of the self-normalised $\psip$ yield in pp
collisions at $\s$ = 13 TeV is reported in [3]. In this analysis, the $\psip$
mesons are reconstructed from their decay channel into dimuons in the rapidity
window 2.5 < $y\rm{{}_{lab}}$ < 4.0, while the charged-particle multiplicity
is measured at midrapidity |$\eta\rm{{}_{lab}}$| < 1\. The $\psip$ yields and
the charged-particle multiplicity are normalised to their respective average
values evaluated in the integrated multiplicity bin. As shown in Fig.1(a), the
$\psip$ normalised yields increase with charged-particle multiplicity within
uncertainties. The resulting trend is compatible with a linear increase
depicted with a dashed line in the figure. The result is compared to PYTHIA
8.2 calculations [4], with and without color reconnection (CR). PYTHIA 8.2
event generator describes the production of the charm quark during the initial
hard processes and with MPI. The number of charged particles produced is
proportional to the number of events containing MPI. In PYTHIA 8.2, the color
reconnection scenario represents the fusion of the final state partons
produced via MPI. In Fig.1(a) PYTHIA 8.2 calculations, both with and without
color reconnection, show similar behavior for the $\psip$ production as a
function of multiplicity. The trend of the calculation is compatible with
measurements up to 5 times the average multiplicity.

(a) $\psi(2S)$ self-normalised yield.

(b) $\psi(2S)$-to-J/$\psi$ self-normalised yield ratio.

Figure 1:  (a) self-normalised $\psi(2S)$ yields, measured at forward
rapidity, as a function of the midrapidity charged-particle multiplicity, in
pp collisions at $\s$ = 13 TeV [3]. (b) $\psi(2S)$-to-J/$\psi$ self-normalised
yield ratios, measured at forward rapidity, as a function of the midrapidity
charged-particle multiplicity [3]. The results are compared to PYTHIA 8.2
calculations and comovers model [4, 5].

The $\psi(2S)$-over-J/$\psi$ self-normalised yield ratio as a function of the
midrapidity charged-particle multiplicity is presented in Fig.1(b). The double
ratio measurement allows one to disentangle any possible differences in the
multiplicity dependence between the charmonium excited and ground states. The
resulting trend is compatible with a flat behavior, within uncertainties. The
result is compared to PYTHIA 8.2 event generator calculations which predict a
similar behavior for both particles within uncertainties. The result is also
compared to the comovers model calculation [5]. In the comovers model,
charmonia interact with the final state particles that are comoving with them.
The dissociation probability depends on the size of charmonioum state and the
density of comoving particles. The comovers model predicts a stronger
suppression for the $\psi(2S)$ with respect to J/$\psi$, as the former has a
larger size than the J/$\psi$ [6]. The multiplicity dependence of the J/$\psi$
production in pp collisions at $\s$ = 13 TeV can be found in [7].

(a) $\psi(2S)$ self-normalised yield in pp and p–Pb collisions.

(b) $\psi(2S)$ self-normalised yield Pb-going direction.

Figure 2: (a) $\psip$ self-normalised yields as a function of the charged-
particle multiplicity in pp collisions at $\s$ = 13 TeV and in p–Pb collisions
at $\snn$ = 8.16 TeV, for both p- and Pb-going direction. (b) $\psip$ self-
normalised yields as a function of the charged-particle multiplicity in p–Pb
collisions, for the Pb-going direction, at $\snn$ = 8.16 TeV [3]. The result
is compared to the Percolation + comover + EPS09 model [5, 6, 8].

From the theoretical point of view, the production of hadrons in pA collisions
is expected to be affected by several nuclear effects, so called cold nuclear
matter effects. For example, the gluon shadowing [9] is responsible for the
suppression of charmonia in p–Pb collisions with respect to pp collisions.
This is due to the modification of the nuclear parton distribution functions
in the nucleus compared to free nucleons. The percolation model [6] represents
the initial state partons as strings that have a specific size. These strings
can interact and overlap at sufficiently high energies, reducing the total
effective number of strings and suppressing the charmonioum production. The
effects above are expected to affect equally both J/$\psi$ and $\psi(2S)$
production, as they are sensitive to the initial state partons. The
multiplicity dependence of the $\psip$ production in p–Pb collisions at $\snn$
= 8.16 TeV is presented in [3]. The dataset used in this analysis was recorded
with two different beam configurations. The first one is obtained when the
proton beam is directed from the ALICE interaction point (IP) toward the muon
spectrometer, the so-called p-going direction. In this case, the charmonium
yields are reconstructed in the center of mass rapidity window 2.03 <
$y\rm{{}_{cms}}$ < 3.53. The second configuration is obtained by inverting the
direction of the colliding beams, i.e. the direction of the Pb beam is from IP
toward the muon spectrometer (Pb-going direction). The charmonium yields in
this case are reconstructed in the backward center of mass rapidity -4.46 <
$y\rm{{}_{cms}}$ < -2.96. Fig.2(a) shows the multiplicity dependence of
$\psip$ yields in p–Pb collisions for both p-going and Pb-going directions
compared to the result from pp collisions. The $\psip$ yield shows an
increasing trend with increasing charged-particle multiplicity for both
configurations. The trend of the result is compatible with the measurements
from pp collisions within uncertainties. Fig.2(b) presents the $\psip$ self-
normalised yields as a function of charged-particle multiplicity compared to
the percolation coupled to comover model and using the EPS09 nuclear parton
distribution function (nPDF). The Percolation + comover + EPS09 model is
compatible with the result within uncertainties. The self-normalised
$\psi(2S)$-over-J/$\psi$ yield ratio, shown in [3], is compatible with a flat
distribution within the large experimental uncertainties for both backward and
forward rapidity in p–Pb collisions. The comovers model suggests a stronger
suppression for the $\psi(2S)$ with respect to the J/$\psi$, especially at
high multiplicities (see [10] for the J/$\psi$ multiplicity dependence in p–Pb
collisions at $\snn$ = 8.16 TeV). A stronger suppression is expected at
backward rapidity as compared to forward rapidity due to the larger density of
soft particles, which increases the probability of interaction with charmonia
and, therefore, the dissociation probability. Within the large experimental
uncertainties, the comover predictions are compatible with the data (although
there is no evidence for a decreasing trend with multiplicity of the
$\psi(2S)$-over-J/$\psi$ ratio).

##  4 Conclusion

The $\psip$ self-normalised yields as a function of the charged-particle
multiplicity in pp and p–Pb collisions at $\s$ = 13 TeV and $\snn$ = 8.16 TeV,
respectively, are discussed. The results show an increase in the self-
normalised yields as a function of multiplicity for both data sets. This
increase is compatible with a linear trend within uncertainties. The trend of
the result is described by several calculations, which model initial state
effects (including MPI) and/or final state effects. The multiplicity
dependence of $\psi(2S)$-over-J/$\psi$ self-normalised yield ratio is also
reported. The trend of the double ratio is flat within uncertainties,
suggesting a similar behavior for the ground and excited charmonium states as
a function of multiplicity. The LHC Run 3 will provide a data set with higher
statistics allowing for significantly more precise measurements with an
increased multiplicity reach, therefore adding more constraints on the MPI
modeling.

## References

  * [1] ALICE Collaboration, B. B. Abelev et al., “Long-range angular correlations of $\rm\pi$, K and p in p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV”, Phys. Lett. B 726 (2013) 164–177, arXiv:1307.3237 [nucl-ex].
  * [2] ALICE Collaboration, B. B. Abelev et al., “Performance of the ALICE Experiment at the CERN LHC”, Int. J. Mod. Phys. A 29 (2014) 1430044, arXiv:1402.4476 [nucl-ex].
  * [3] ALICE Collaboration, “Measurement of $\psi$(2S) production as a function of charged-particle pseudorapidity density in pp collisions at $\sqrt{s}$ = 13 TeV and p-Pb collisions at $\sqrt{s_{\rm NN}}$ = 8.16 TeV with ALICE at the LHC”, arXiv:2204.10253 [nucl-ex].
  * [4] T. Sjostrand, “An introduction to PYTHIA 8.2”, Comput. Phys. Commun. 191 (2015) 159–177, arXiv:1410.3012 [hep-ph].
  * [5] E. G. Ferreiro, “Excited charmonium suppression in proton–nucleus collisions as a consequence of comovers”, Phys. Lett. B 749 (2015) 98–103, arXiv:1411.0549 [hep-ph].
  * [6] E. G. Ferreiro and C. Pajares, “High multiplicity $pp$ events and $J/\psi$ production at LHC”, Phys. Rev. C 86 (2012) 034903, arXiv:1203.5936 [hep-ph].
  * [7] ALICE Collaboration, S. Acharya et al., “Forward rapidity J/$\psi$ production as a function of charged-particle multiplicity in pp collisions at $\sqrt{s}$ = 5.02 and 13 TeV”, JHEP 06 (2022) 015, arXiv:2112.09433 [nucl-ex].
  * [8] H. J. Drescher, “Parton based Gribov-Regge theory”, Phys. Rept. 350 (2001) 93–289, arXiv:hep-ph/0007198.
  * [9] R. Vogt, “Shadowing and absorption effects on $J/\psi$ production in dA collisions”, Phys. Rev. C 71 (2005) 054902, arXiv:hep-ph/0411378.
  * [10] ALICE Collaboration, S. Acharya et al., “J/$\psi$ production as a function of charged-particle multiplicity in p-Pb collisions at $\sqrt{\textit{s}_{\rm NN}}~{}=~{}8.16$ TeV”, JHEP 09 (2020) 162, arXiv:2004.12673 [nucl-ex].


# Direct molecular gas dynamics simulations of re-entry vehicles via the
Boltzmann equation

Tarik Dzanic111Postdoctoral fellow Center for Applied Scientific Computing,
Lawrence Livermore National Lab, Livermore, CA 94550, USA  Luigi
Martinelli222Professor; Associate Fellow, AIAA Department of Mechanical and
Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA

###### Abstract

This work explores the feasibility of performing three-dimensional molecular
gas dynamics simulations of hypersonic flows such as re-entry vehicles through
directly solving the six-dimensional nonlinear Boltzmann equation closed with
the BGK (Bhatnagar–Gross–Krook) collision model. Through the combination of
high-order unstructured spatial discretizations and conservative discrete
velocity models as well as their efficient implementation on large-scale GPU
computing architectures, we demonstrate the ability to simulate unsteady and
non-equilibrium three-dimensional high-speed flows at a feasible computational
cost through a unified numerical framework. We present the results of high-
order simulations of the Apollo capsule at realistic re-entry conditions from
the AS-202 mission flight path, including the steady non-equilibrium flow in
the high-altitude regime at a Mach number of 22.7 and a Reynolds number of
43,000 as well as the unsteady turbulent flow in the low-altitude regime at a
Mach number of 8 and a Reynolds number of 550,000. The results show the
validity of the approach over the entire range of a typical re-entry
trajectory from the rarefied to the continuum limit, the ability to directly
resolve strong shocks profiles without numerical shock capturing techniques,
and the ability of resolving small-scale unsteady flow structures in the
inertial range.

##  1 Nomenclature

$\alpha$ = | angle of attack
---|---
$C$ = | particle collision operator
$d$ = | spatial dimensionality
$D$ = | capsule diameter
$\Delta x$ = | mesh spacing
$\delta$ = | number of internal degrees of freedom
$E$ = | macroscopic total energy
$f$ = | particle distribution function
$\mathbf{f}$ = | discrete particle distribution function
$F$ = | reduced velocity particle distribution function
$F^{I}$ = | common interface flux
$\mathbf{F}$ = | flux vector
$g$ = | equilibrium particle distribution function
$g_{\mathbf{u}}$ = | reduced velocity equilibrium particle distribution function
$g_{\zeta}$ = | reduced internal energy equilibrium particle distribution function
$\mathbf{g}$ = | discrete equilibrium particle distribution function
$G$ = | reduced internal energy particle distribution function
$\mathbf{h}$ = | flux reconstruction correction function
$\mathbf{I}$ = | identity matrix
$K$ = | Kelvin
$Kn$ = | Knudsen number
$Kn_{\text{GLL}}$ = | gradient-length local Knudsen number
$\lambda$ = | particle mean free path
$\Lambda$ = | Lambda function
$m$ = | velocity dimensionality
$M$ = | Mach number
$\mu$ = | dynamic viscosity
$N_{e}$ = | number of mesh elements
$N_{f}$ = | number of flux nodes
$N_{s}$ = | number of solution nodes
$N_{u}$ = | number of velocity nodes
$\Omega^{\mathbf{u}}$ = | velocity domain
$\Omega^{\mathbf{x}}$ = | spatial domain
$\Omega^{\mathbf{x}}_{k}$ = | element in mesh
$P$ = | macroscopic pressure
$Pr$ = | Prandtl number
$\Pi$ = | pressure tensor
$\mathbf{q}$ = | macroscopic primitive variables
$\mathbf{Q}$ = | macroscopic conserved variables
$\mathbf{Q}^{\prime}$ = | modified macroscopic conserved variables
$Re$ = | Reynolds number
$\rho$ = | macroscopic density
$\rho\mathbf{U}$ = | macroscopic momentum
$t$ = | physical time
$\tau$ = | collision time scale
$\theta$ = | macroscopic scaled temperature
$\theta_{\text{ref}}$ = | reference macroscopic scaled temperature
$\theta_{S}$ = | Sutherland scaled temperature
$\mathcal{T}$ = | elliptic temperature tensor
$\mathbf{u}$ = | velocity variable
$u_{\max}$ = | maximum discrete particle velocity
$\mathbf{x}$ = | spatial variable
$\mathbf{x}^{s}$ = | spatial solution nodes
$\mathbf{x}^{f}$ = | spatial flux nodes
$\zeta$ = | internal energy variable

##  2 Introduction

The comprehensive investigation of aerothermodynamic behavior of spacecraft
during atmospheric re-entry remains a critical component of the engineering
development cycle, providing crucial information for structural and material
design as well as guiding the development of thermal protection systems for
the vehicles. However, accurate simulations of flows in high-altitude,
hypersonic conditions can be an exceptionally difficult undertaking, primarily
driven by the complexity of the multi-scale, non-equilibrium nature of the
flow in these conditions. Moreover, prediction of laminar-turbulent transition
in hypersonic boundary layers is paramount as it determines the heat transfer
and friction drag, and the prediction of the dynamics of vehicles at
hypersonic speeds may be substantially hampered by uncertainties in the
modeling of turbulence and turbulent transition, which is still poorly
understood [1].

Standard computational fluid dynamics approaches, which rely on solving the
Navier–Stokes equations governing the dynamics of fluid flows, inherently
assume that the fluid in question can be treated as continuum, such that the
flow can be adequately described through its macroscopic state (e.g., density,
momentum, energy). Under re-entry conditions, this assumption can become
severely ill-posed as the flow experiences strong thermodynamic non-
equilibrium, and its behavior can drastically vary from the predictions given
by governing equations that do rely on the continuum assumption. In these
cases, to provide an accurate description of the underlying flow physics, it
is necessary to revert to more general governing equations that are derived
from the kinetic theory of gases which can account for the non-equilibrium
nature of the flow.

The governing equations of molecular gas dynamics, which underpin the
macroscopic behavior of the flow, provide a means of accurately describing
non-equilibrium flows while seamlessly recovering the hydrodynamic equations
in the continuum limit. As such, these models are the only ones truly capable
of describing the flow physics across the entire flight path, from the high-
altitude, rarefied regime to the low-altitude, continuum regime. One such
approach is provided by approximation of the Boltzmann equation, which gives a
statistical description of particle transport and collision at the mesoscale
level. However, the complexity of the Boltzmann equation, primarily driven by
the computational cost associated with its high-dimensional nature and the
difficulty of approximating the collision process, severely limits its
applicability. While more complex flow problems have been attempted with
probabilistic methods such as direct simulation Monte Carlo approaches, the
direct deterministic approximation of the Boltzmann equation has been widely
considered computationally intractable for practical high-speed applications,
and their use has generally been restricted to simpler, two-dimensional flows.

The aim of this work is to demonstrate the feasibility of performing three-
dimensional molecular gas dynamics simulations of high-speed flows such as re-
entry vehicles through directly solving the Boltzmann equation. This
capability stems from recent developments by the authors for approximating the
Boltzmann equation closed with the Bhatnagar–Gross–Krook (BGK) model which
combines highly-efficient, high-order unstructured discontinuous spectral
element method with a conservative discrete velocity model, the combination of
which results in a low-memory, compute-heavy numerical algorithm which is
well-suited for large-scale GPU computing architectures. As such, the authors
have presented results of complex three-dimensional turbulent flows and high-
speed flows simulated without any numerical shock capturing approaches by
directly approximating the Boltzmann equation, both of which are entirely
feasible on modern computing hardware [2, 3]. These results pave the way for
investigating non-equilibrium flow around spacecraft in re-entry conditions
through deterministic approximations of molecular gas dynamics.

In this work, we present the results of high-order three-dimensional numerical
experiments of the flow around an AS-202 Apollo capsule at realistic re-entry
conditions obtained through direct approximation of the Boltzmann–BGK
equation. The first simulation is performed in the high-altitude region at a
Mach number of 22.7 and a Reynolds number of 43,000, which results in strong
non-equilibrium flow that is not accurately resolved by continuum fluid
dynamics approximations [4]. The second simulation is performed in the low-
altitude region at a Mach number of 8 and a Reynolds number of 550,000, which
results in complex unsteady flow separation and turbulent wake flow. While
three-dimensional re-entry flows have been attempted using kinetic/hybrid
schemes in works such as that of Yuan and Zhong [5], Zhang et al. [6], and
Chen et al. [7], these approaches have only been used for steady flows and
with low-order numerical schemes. To the authors’ knowledge, this work is the
first to showcase the ability to directly solve the Boltzmann equation for
unsteady hypersonic turbulent flows, achieving so with a high-order
unstructured numerical scheme. The results presented in this work, in
conjunction with previous work by the authors, highlight some unique
properties of the proposed model, such as: i) validity over the entire range
of a typical re-entry trajectory, from the rarefied to the continuum limit;
ii) ability of directly resolving strong shocks profiles without numerical
shock capturing techniques; and iii) ability of resolving small-scale unsteady
flow structures in the inertial range. The overarching goal of this work is to
showcase a novel approach for tackling complex high-speed flow problems which
is more suited for handling the high-temperature aerothermodynamic effects
that are encountered in these applications and may lead to enhancing our
fundamental understanding of flow phenomena such as transition to turbulence
through a radically different perspective of flow dynamics encoded in the
high-dimensional solution of the Boltzmann equation.

##  3 Methodology

###  3.1 Governing equations

The Boltzmann equation, which gives a deterministic description of molecular
gas dynamics, can be given as

$\partial_{t}f(\mathbf{x},\mathbf{u},\zeta,t)+\mathbf{u}{\cdot}\nabla
f=\mathcal{C}(f,f^{\prime}),$ (1)

where $\mathbf{x}$ represents a $d$-dimensional physical space, $\mathbf{u}$
represents an $m$-dimensional physical space, $\zeta$ represents a scalar
internal energy, $f(\mathbf{x},\mathbf{u},\zeta,t)$ is a scalar particle
distribution function, and $\mathcal{C}(f,f^{\prime})$ is some collision
operator that accounts for intermolecular interactions [8]. The distribution
functions represents a probability measure for a particle existing at a given
location $\mathbf{x}$ traveling at a given velocity $\mathbf{u}$ with an
internal energy $\zeta$. The moments of this distribution function recover the
macroscopic state of the system, i.e.,

$\mathbf{Q}(\mathbf{x},t)=\left[\rho,\rho\mathbf{U},E\right]^{T}=\int_{\mathbb{R}^{m}}\int_{0}^{\infty}f(\mathbf{x},\mathbf{u},\zeta,t)\
\boldsymbol{\psi}(\mathbf{u},\zeta)\ \mathrm{d}\zeta\ \mathrm{d}\mathbf{u}$
(2)

where $\rho$ is the density, $\rho\mathbf{U}$ is the momentum vector, $E$ is
the total energy, and
$\boldsymbol{\psi}(\mathbf{u},\zeta)\coloneqq[1,\mathbf{u},(\mathbf{u}\cdot\mathbf{u})/2+\zeta]^{T}$
is the vector of collision invariants. From this, a vector of primitive
variables can be defined as $\mathbf{q}=[\rho,\mathbf{U},P]$, where
$\mathbf{U}=\rho\mathbf{U}/\rho$ is the macroscopic velocity,
$P=(\gamma-1)(E-\rho\mathbf{U}\cdot\mathbf{U}/2)$ is the pressure, and
$\gamma$ is the specific heat ratio. We use the notation $\mathbf{u}$ to
denote the microscopic velocities and $\mathbf{U}$ to denote the macroscopic
velocities.

The collision operator is approximated using the Bhatnagar–Gross–Krook (BGK)
model [9], given as

$\mathcal{C}(f,f^{\prime})\approx\frac{g(\mathbf{x},\mathbf{u},\zeta,t)-f(\mathbf{x},\mathbf{u},\zeta,t)}{\tau},$
(3)

where $g(\mathbf{x},\mathbf{u},\zeta,t)$ is the equilibrium distribution
function (EDF) and $\tau$ is the collision time scale. To recover a non-unit
Prandtl number and account for the effects of internal degrees of freedom, the
EDF is computed through the ellipsoidal BGK (ES-BGK) approach [10] in
conjunction with the internal energy model of Baranger et al. [11]. The EDF
can be represented as a product of the monatomic ES-BGK EDF $g_{\mathbf{u}}$
and the internal energy EDF $g_{\zeta}$ as

$g(\mathbf{x},\mathbf{u},\zeta,t)=g_{\mathbf{u}}(\mathbf{x},\mathbf{u},t)\times
g_{\zeta}(\mathbf{x},\zeta,t).$ (4)

The ES-BGK EDF can be computed as

$g_{\mathbf{u}}(\mathbf{x},\mathbf{u},t)=\frac{\rho(\mathbf{x},t)}{\sqrt{\text{det}(2\pi\mathcal{T})}}\exp\left[-\frac{1}{2}\left[\mathbf{u}-\mathbf{U}(\mathbf{x},t)\right]^{T}\mathcal{T}^{-1}\left[\mathbf{u}-\mathbf{U}(\mathbf{x},t)\right]\right],$
(5)

where the temperature tensor $\mathcal{T}$ is computed as a combination of the
temperature $\theta=P/\rho$ and the density-normalized pressure tensor $\Pi$,
i.e.,

$\mathcal{T}=\frac{1}{Pr}\theta\mathbf{I}+(1-\frac{1}{Pr})\frac{\Pi}{\rho},\quad\text{where}\quad\Pi_{ij}=\int_{\mathbb{R}^{m}}f(u_{i}-U_{i})(u_{j}-U_{j})\
\mathrm{d}\mathbf{u}.$ (6)

The Prandtl number is set as $Pr=0.71$. From this model, the collision time
scale can be related to the dynamic viscosity $\mu$ as

$\tau=\frac{1}{Pr}\frac{\mu}{P}.$ (7)

In this work, the Sutherland viscosity model [12] is used, given by the
relation

$\mu=\mu_{\text{ref}}\left(\frac{\theta}{\theta_{\text{ref}}}\right)^{\frac{3}{2}}\frac{\theta_{\text{ref}}+\theta_{S}}{\theta+\theta_{S}},$
(8)

where $(\cdot)_{\text{ref}}$ denotes the reference quantities and $\theta_{S}$
denotes the (scaled) Sutherland temperature.

The internal energy EDF is computed as

$g_{\zeta}(\mathbf{x},\zeta,t)=\Lambda(\delta)\left(\frac{\zeta}{\theta(\mathbf{x},t)}\right)^{\frac{\delta}{2}-1}\frac{1}{\theta(\mathbf{x},t)}\exp\left(-\frac{\zeta}{\theta(\mathbf{x},t)}\right),$
(9)

where $\delta\geq 0$ is the number of internal degrees of freedom and
$\Lambda(\delta)=1/\Gamma(\delta/2)$ is a normalization factor to ensure that
$g$ recovers the same macroscopic state as $f$. The number of internal degrees
of freedom are not necessarily fixed nor required to be an integer, such that
spatially varying and temperature-dependent internal degrees of freedom are
possible. However, the addition of an internal energy model further increases
the dimensionality of the Boltzmann equation, which would significantly
increase the associated computational cost. Since it is typically not
necessary to know the actual distribution within the internal energy domain,
only its effect on the macroscopic quantities, it is beneficial to utilize a
reduced distribution technique for the internal energy domain [11]. In this
approach, the distribution function $f(\mathbf{x},\mathbf{u},\zeta,t)$ (and
its evolution) is integrated with respect to the internal energy variable
$\zeta$ to yield a pair of distribution functions $[F,G]^{T}$ as

$\begin{bmatrix}F(\mathbf{x},\mathbf{u},t)\\\
G(\mathbf{x},\mathbf{u},t)\end{bmatrix}=\int_{0}^{\infty}\begin{bmatrix}1\\\
\zeta\end{bmatrix}f(\mathbf{x},\mathbf{u},\zeta,t)\ \mathrm{d}\zeta.$ (10)

Critically, these reduced distribution functions are only defined over the
spatial and velocity domains, not the internal energy domain, such that the
computational cost and memory requirements are only twice the cost of the
monatomic case. The evolution of these reduced distribution function can be
computed as

$\partial_{t}\begin{bmatrix}F(\mathbf{x},\mathbf{u},t)\\\
G(\mathbf{x},\mathbf{u},t)\end{bmatrix}+\mathbf{u}\cdot\nabla\begin{bmatrix}F(\mathbf{x},\mathbf{u},t)\\\
G(\mathbf{x},\mathbf{u},t)\end{bmatrix}=\frac{1}{\tau}\begin{bmatrix}\hphantom{\frac{\delta\theta}{2}}g_{\mathbf{u}}(\mathbf{x},\mathbf{u},t)-F(\mathbf{x},\mathbf{u},t)\\\
\frac{\delta\theta}{2}g_{\mathbf{u}}(\mathbf{x},\mathbf{u},t)-G(\mathbf{x},\mathbf{u},t)\end{bmatrix},$
(11)

where it can be seen that the equilibrium distribution function of the two
components varies only by a factor of $\frac{\delta\theta}{2}$. From this
formulation, the macroscopic quantities can be computed as

$\begin{bmatrix}\rho\\\ \rho\mathbf{U}\\\
E\end{bmatrix}=\int_{\mathbb{R}^{m}}\begin{bmatrix}1\\\
\mathbf{u}\cdot\mathbf{u}\\\
(\mathbf{u}\cdot\mathbf{u})/2\end{bmatrix}F+\begin{bmatrix}0\\\ 0\\\
1\end{bmatrix}G\ \mathrm{d}\mathbf{u}.$ (12)

Similarly, the pressure tensor $\Pi$ can be computed as

$\Pi_{ij}=\int_{\mathbb{R}^{m}}F(u_{i}-U_{i})(u_{j}-U_{j})\
\mathrm{d}\mathbf{u}.$ (13)

###  3.2 Numerical method

The underlying numerical scheme in this work is the method of Dzanic et al.
[2]. This approach combines a high-order spatial discretization with a
discretely-conservative velocity model, the combination of which results in a
low-memory, compute-heavy numerical algorithm which is well-suited for large-
scale GPU computing architectures. A brief overview of this scheme is
presented here, with a schematic of the approach shown in Fig. 1.

$\Omega^{\mathbf{x}}$$\Omega^{\mathbf{x}}_{k}$$\Omega^{\mathbf{u}}$$u$$v$

Figure 1:  Schematic of a two-dimensional phase space discretization using an
unstructured spatial domain $\Omega^{\mathbf{x}}$ with $\mathbb{P}_{2}$
elements and a velocity domain $\Omega^{\mathbf{u}}$ with $N_{u}=16^{2}$.
Circles denote the spatial solution nodes (red), interface flux nodes (blue),
and velocity space nodes (orange), respectively.

The spatial domain is discretized using a high-order flux reconstruction
scheme [13], a generalization of the nodal discontinuous Galerkin method [14].
When combined with a nodal velocity discretization, the coupled spatial and
velocity domains reduce to a set of decoupled advection equations for each
nodal velocity point $\mathbf{u}_{0}$ with a nonlinear source term $S$, i.e.,

$\partial_{t}f(\mathbf{x},\zeta,t)+\boldsymbol{\nabla}\mathbf{F}(f)f=S,\quad\mathbf{F}(f)=\mathbf{u}_{0}f.$
(14)

For this scheme, the spatial domain $\Omega^{\mathbf{x}}$ is partitioned into
$N_{e}$ elements $\Omega^{\mathbf{x}}_{k}$ such that
$\Omega^{\mathbf{x}}=\bigcup_{N_{e}}\Omega^{\mathbf{x}}_{k}$ and
$\Omega^{\mathbf{x}}_{i}\cap\Omega^{\mathbf{x}}_{j}=\emptyset$ for $i\neq j$,
shown on the left-hand side of Fig. 1. The distribution function within each
element is approximated by an order $\mathbb{P}_{p}$ polynomial as

$f(\mathbf{x})=\sum_{i=1}^{N_{s}}f(\mathbf{x}^{s}_{i}){\phi}_{i}(\mathbf{x}),$
(15)

where $\mathbf{x}^{s}_{i}\in\Omega^{\mathbf{x}}_{k}\ \forall\
i\in\\{1,...,N_{s}\\}$ is a set of $N_{s}$ solution nodes and
${\phi}_{i}(\mathbf{x})$ is a set of nodal basis functions. The flux is then
computed as

$\mathbf{F}(\mathbf{x})=\mathbf{u}_{0}f(\mathbf{x})+\sum_{i=1}^{N_{f}}\left[F^{I}_{i}-f(\mathbf{x}^{f}_{i})\mathbf{u}_{0}\cdot\mathbf{n}_{i}\right]\mathbf{h}_{i}(\mathbf{x}).$
(16)

Here, $\mathbf{x}^{f}_{i}\in\partial\Omega^{\mathbf{x}}_{k}\ \forall\
i\in\\{1,...,N_{f}\\}$ is a set of $N_{f}$ interface flux nodes,
$\mathbf{n}_{i}$ is their associated outward-facing normal vector, $F^{I}_{i}$
is the upwind numerical flux at the interface, and $\mathbf{h}_{i}$ is the
correction function used to recover the nodal discontinuous Galerkin approach
[13, 14]. The upwind numerical flux is simply computed as

$F_{i}^{I}=\begin{cases}u_{n}f_{i}^{-},\quad\quad\mathrm{if}\ u_{n}>0,\\\
u_{n}f_{i}^{+},\quad\quad\mathrm{else},\end{cases}$ (17)

where $u_{n}=\mathbf{u}_{0}\cdot\mathbf{n}_{i}$ and the superscripts $-$ and
$+$ denote the interior value (from the element of interest) and the exterior
value (from the face-adjacent element) of the solution at the interface,
respectively.

The velocity domain is approximated by a set of discrete velocity nodes
distributed across a Cartesian grid in velocity space, shown on the right-hand
side of Fig. 1. The extent of the velocity domain is chosen such as to
adequately encompass the majority of the distribution function, typically
taken as some factor of the maximal thermal velocity and macroscopic velocity.
The velocity domain is then taken as $[-u_{\max},u_{\max}]^{m}$, where
$u_{\max}$ is the maximum velocity extent. The governing equations are then
solved for across the discrete velocity nodes, each corresponding to a
spatially decoupled advection equation as in Eq. 14. These spatially decoupled
systems are linked together by the collision operator which acts as a
nonlinear source term across the velocity space.

However, the BGK (and ES-BGK) collision operators require integrating the
particle distribution function across the velocity space to compute the
moments. If this integration is not exact for the moments corresponding to the
macroscopic state, the equilibrium distribution function $g$ may not possess
the same macroscopic state as $f$. As a result, the scheme would not be
conservative with respect to the macroscopic state, and these conservation
errors can ultimately result in numerical instabilities [2]. A typical remedy
for this problem is to increase the resolution in velocity space to reduce the
integration errors or choose the velocity basis and/or nodal quadrature points
judiciously such as to integrate particular moments of the distribution
function exactly (e.g., moment-based approximations/Gauss–Hermite quadrature).
However, this either results in a prohibitively large computational cost or a
velocity discretization which may be ill-suited, particularly so for high Mach
number flows [15, 16].

Instead, we utilize the discrete velocity model of Mieussens [17] which
ensures that the scheme remains discretely conservative regardless of the
resolution and choice of velocity nodal distribution, which can drastically
reducing the computational cost and resolution requirements of standard nodal
velocity approximations and offer a more general and robust approach than
moment/Gauss–Hermite based methods. In brief, this approach requires
iteratively solving for a modified equilibrium distribution function that
possesses the same macroscopic moments as $f$ under discrete integration. As a
result, this can transform a memory-heavy velocity space approximation into a
lower memory, compute-heavy task which is optimal for modern GPU
architectures.

The goal of the discrete velocity model is to find this modified EDF
$\mathbf{g}^{\prime}$ which satisfies the _discrete compatibility condition_ ,
i.e.,

$\mathbf{M}\cdot\left[\mathbf{f}\otimes\boldsymbol{\psi}\right]=\mathbf{M}\cdot\left[\mathbf{g}^{\prime}\otimes\boldsymbol{\psi}\right].$
(18)

where

$\mathbf{M}\cdot\mathbf{f}\approx\int_{\mathbb{R}^{m}}f(\mathbf{u})\
\mathrm{d}\mathbf{u}$ (19)

is a discrete integration operator with strictly-positive entries. Due to the
use of a uniform Cartesian distribution in the velocity domain, we utilize the
trapezoidal rule which offers spectral convergence for smooth, compactly-
supported, periodic functions on uniform grids [18]. It was shown in Mieussens
[17] that the EDF which satisfies both the discrete compatibility condition
and the discrete form of the H-theorem is equivalent to the analytic EDF
formed around a modified macroscopic state $\mathbf{Q}^{\prime}$ (i.e.,
$\mathbf{g}(\mathbf{Q}^{\prime})$) which converges to the true macroscopic
state in the limit of infinite resolution. As a result, finding this modified
EDF typically can be framed simply as the solution of a $d+2$ dimensional
nonlinear system, irrespective of the dimensionality of the velocity space –
note that in some scenarios, the ES-BGK model may require a nonlinear solve of
a marginally higher dimensionality if the velocity space does not possess
certain symmetry properties [19]. This low-dimensional, compute-heavy
optimization process can be rapidly performed on GPUs using root-finding
approaches such as Newton’s method.

###  3.3 Geometry & flow conditions

The numerical experiments focus on the three-dimensional simulation of the
AS-202 Apollo capsule at realistic operating conditions. A schematic of the
cross-sectional profile of the capsule geometry is shown in Fig. 2, normalized
to a unit non-dimensional diameter corresponding to a physical diameter of
$3.912$ m. The capsule geometry consists of a blunt spherical forebody
followed by a $33^{\circ}$ conical afterbody. We use the notation that
fore/aft is with respect to the flow-oriented reference frame, not with
respect to the reference frame of the astronauts which are positioned in the
opposite direction.

$1$$362/391$$469/391$19.6/391 R+2.3/391 R+++$33^{\circ}$

Figure 2: Cross-section schematic of the AS-202 Apollo capsule geometry
normalized to unit diameter.

The operating conditions are taken from the AS-202 experimental flight data of
Hillje [20] as reported by Wright et al. [4]. The flight data is reported over
a range of flight times corresponding to $4350-5200$ seconds after launch,
during which the capsule travels from an altitude of $122$ to $8$ kilometers.
During this time period, the operating Mach number steadily decreases from
$M=31$, while the Reynolds number (based on the capsule diameter) increases to
a maximum of $Re=10^{7}$. The first simulation operating conditions were
chosen to match the experimental flight data at 4700 seconds, which
corresponds to an altitude of $77.2$ kilometers, a Mach number of $M=22.7$, a
Reynolds number of $Re=43,000$, a freestream temperature of $203$ K, an angle
of attack of $\alpha=18.5^{\circ}$, and a zero side-slip angle. These
conditions were chosen as the high Mach number and moderate Reynolds number
result in strong non-equilibrium flow behavior, such that continuum fluid
dynamics approximations yield erroneous predictions [4]. The second simulation
operating conditions were chosen similarly to the experimental flight data at
5000 seconds, which corresponds to an altitude of $44.8$ kilometers and a Mach
number of $M=8$. The Reynolds number was set to $Re=550,000$ to yield unsteady
turbulent flow in the continuum regime, although this was marginally less than
the true experimental Reynolds number which was on the order of $Re=10^{6}$.
As the purpose of this second experiment was primarily to show the validity of
the approach for unsteady flows across the continuum regime instead of exactly
matching the experimental conditions, the freestream temperature, angle of
attack, and side-slip angle were kept the same as with the first simulation.

We note here that the numerical experiments in this work are primarily meant
to represent a preliminary proof of capability of directly solving the
Boltzmann equation for such flows. As such, certain aspects of the flow as
neglected, although these are strictly modeling assumptions and do not impact
the computational feasibility of the proposed approach. In particular, we
neglect the effects of multi-species reactions/chemical non-equilibrium and
vibrational degrees of freedom, both of which were considered in the numerical
experiments of Wright et al. [4], as well as ionization, shock-layer
radiation, and material ablation. Additionally, the wall is simply assumed to
be isothermal with a temperature equal to the freestream, whereas the work of
Wright et al. [4] assumes radiative equilibrium with surface catalysis effects
which is more representative of the experimental flow conditions. As such, the
results in this work are primarily used for qualitative comparison of the flow
fields to the work of Wright et al. [4]. A direct comparison of more complex
flow effects such as near-wall flow behavior and an evaluation of the accuracy
of predicting quantitative metrics such as wall heat transfer effects are
contingent on future work in employing these additional modeling assumptions
in the proposed approach.

###  3.4 Implementation & problem setup

Due to the zero side-slip angle, a half-domain was considered with symmetry
enforced along the pitch plane ($z=0$). An initial highly-resolved
unstructured mesh was generated upon which the Navier–Stokes equations were
solved to give an initial indication on the flow field in the problem. A
coarser shock-fitted mesh was then generated for the Boltzmann–BGK solver, the
results of which were used to guide further mesh improvements. After several
iterations of solving the flow field via the Boltzmann–BGK approach and
adapting the preliminary meshes, the mesh shown in Fig. 3 was generated and
used for the resulting numerical experiments. We note here that the mesh was
fitted primarily for the high-altitude/Mach simulation, and the same mesh was
used for both simulations.

Figure 3:  Visualization of the high-order finite element mesh on capsule
surface and symmetry plane.

The mesh consisted of a unstructured triangular surface region with a near-
wall prismatic boundary layer region. The domain was primarily filled with
tetrahedral elements except near the shock region where shock-aligned
prismatic elements were used. The mesh consisted of $14,033$ primastic
elements and $25,680$ tetrahedral elements with a minimum mesh spacing of
$\Delta x=4\cdot{10^{-3}}$. This resolution level was chosen to mimic the mesh
used in Wright et al. [4], with approximately 1.3 million spatial degrees of
freedom at the operating order for the high-altitude case and 2.4 million
spatial degrees of freedom for the low-altitude case, which, when combined
with the high-order spatial accuracy of the proposed approach, is expected to
yield more performant resolving power. Note that the mesh appears
significantly coarser than a finite volume mesh of equivalent resolution due
to the multiple degrees of freedom per element in high-order finite element
methods (e.g., a single prismatic element using a third-order polynomial
approximation contains 40 degrees of freedom), such that the effective
resolution is on the order of $\Delta x/(p+1)$. To account for the coarseness
of the mesh in high-order finite element methods, we utilize quadratic
elements in the near-wall region to better account for the curvature of the
geometry.

The operating approximation order was set to $\mathbb{P}_{3}$ for the high-
altitude case and $\mathbb{P}_{4}$ for the low-altitude case, yielding a
spatially fourth-order and fifth-order accurate solution, respectively, and
the solution/flux points were placed at the $\alpha$-optimized quadrature
points [14] for simplex elements and Gauss–Lobatto quadrature points for
tensor-product elements. For the high-altitude/Mach simulation, the velocity
space resolution was set to $N_{u}=24^{3}$, resulting in a total of 30 billion
degrees of freedom. The velocity domain extent was set to
$\Omega^{\mathbf{u}}=[-0.6,1.4]\times[-0.8,1.2]\times[-1,1]$. We note here
that results obtained with lower velocity space resolution $N_{u}=20^{3}$
showed marked similarities in the macroscopic flow fields. Therefore, it may
be surmised that the flow was effectively converged in the velocity domain due
to the spectral convergence properties of the approximation, although
quantitative convergence analysis was not performed. For the low-altitude/Mach
simulation, the velocity space resolution was set to $N_{u}=16^{3}$ due to the
lower resolution requirements at lower Mach numbers, and the velocity domain
extent was set identically to the high-altitude simulation. Due to the reduced
velocity space resolution in the low-altitude case, the operating order was
increased from $\mathbb{P}_{3}$ to $\mathbb{P}_{4}$, such that the total
degrees of freedom were approximately the same.

The initial conditions of the flow were normalized such as to yield a unit
density and velocity,

$\displaystyle\rho$ $\displaystyle=1,$ (20a) $\displaystyle u$
$\displaystyle=\cos(\alpha),$ (20b) $\displaystyle v$
$\displaystyle=\sin(\alpha),$ (20c) $\displaystyle w$ $\displaystyle=0,$ (20d)
$\displaystyle P$ $\displaystyle=P_{\text{ref}},$ (20e)

with the pressure and reference viscosity set such as to yield the appropriate
freestream Mach and Reynolds numbers, i.e.,

$\displaystyle P_{\text{ref}}=\frac{1}{\gamma M^{2}},$ (21a)
$\displaystyle\mu_{\text{ref}}=\frac{1}{Re}.$ (21b)

The reference scaled temperature was then set to
$\theta_{\text{ref}}=P_{\text{ref}}/\rho_{\text{ref}}=P_{\text{ref}}$ and the
Sutherland scaled temperature was set to
$\theta_{S}=\theta_{\text{ref}}\frac{110.4}{203}$ which corresponds to a
physical Sutherland temperature of $110.4$ K and a physical freestream
temperature of $203$ K. Diffuse wall boundary conditions were enforced on the
capsule surface, with the wall temperature set equal to the
freestream/reference temperature. At the symmetry plane, symmetry was enforced
using specular (reflecting) wall boundary conditions. At the farfield,
Dirichlet boundary conditions were used with the boundary state set equal to
the initial conditions. All boundary conditions were enforced weakly through
the standard approach in discontinuous spectral element methods where an
exterior state is given and common interface flux is computed via the Riemann
solver. As such, upwinding was performed naturally at the particle level.

To reduce the unnecessary computational cost in resolving the initial
transients in the solution, a Navier–Stokes flow field was first computed on
the mesh and advanced until a quasi-steady state was achieved using the HLL
Riemann solver [21] and the entropy filtering approach [22] for shock
capturing. The Boltzmann solution was then initialized from this flow field by
setting the initial distribution function state to the EDF corresponding to
the macroscopic state of the Navier–Stokes flow field. The operating order was
then gradually increased and the flow as computed by the Boltzmann–BGK
approach was allowed to develop until the operating order was reached. For the
high-altitude case, the flow was then advanced to a steady state, after which
post-processing of the results was performed. For the low-altitude case, the
primary purpose was to showcase the ability to compute unsteady flows. As
such, we simply present results in terms of an instantaneous flow snapshot
after the initial flow transients have been convected away.

Computations were performed using the PyFR codebase [23] as the underlying
solver. For the given problem and number of degrees of freedom, the memory
requirements could be achieved with 16 NVIDIA H100 GPUs. The simulations in
this work were however performed with 80 NVIDIA V100 GPUs. The total compute
time for the high-altitude steady case was approximately 2 days, corresponding
to roughly 3800 GPU hours. We note here that a sizeable portion of this
compute time is spent advancing the solution to a steady state using an
explicit time integration method. As such, the approach is better suited
towards computing unsteady flow features [3], and for the low-altitude
unsteady case, the compute time was approximately 11 hours (900 GPU hours) per
characteristic flow time which corresponds to one flow over diameter.

##  4 Results

###  4.1 High-altitude conditions

####  4.1.1 Macroscopic flow cross-sections

The results for the steady, high-altitude, high Mach case were first analyzed
with respect to the macroscopic flow state on the symmetry plane cross-section
($z=0$), with the contours of density, pressure, velocity magnitude, and
temperature shown in Fig. 4 normalized by the freestream values. The
quintessential bow shock was evident in the flow field with large jumps in the
density, pressure, and temperature across the shock. Due to the ability of the
Boltzmann equation to directly resolve shock structures without any ad hoc
numerical stabilization approach, the solution in the vicinity of the shock
was well-behaved even with the use of a high-order numerical scheme in the
presence of strong shocks. Density ratios of $\mathcal{O}(10)$ and pressure
jumps of $\mathcal{O}(10^{2})$ were observed across the shock, which align
with the predictions of the normal shock equations. The highest density and
pressure ratios were observed at the capsule surface on the windward side,
with density and pressure ratios of over 200 and 800, respectively.

A strong separation region was observed aft of the capsule surface on the
leeward side, with low velocity magnitude and high density and pressure. This
separation region largely corresponded to the location of the recompression
shocks in the wake, with strong density and pressure gradients aligning with
the edge of the separation region. These strong gradients were also well-
resolved by the numerical scheme.

(a) Density

(b) Pressure

(c) Velocity

(d) Temperature

Figure 4: Contours of macroscopic normalized density (top left), pressure (top
right), velocity magnitude (bottom left), and temperature (bottom right) on
the pitch plane ($z=0$) and the capsule surface for the high-altitude flow
conditions.

####  4.1.2 Surface flow behavior

The behavior of the flow at the surface was then visualized through surface
velocity streamlines (i.e., oil flow). The contours of surface slip velocity
as well as the oil flow streamlines are shown in Fig. 5. Due to the local non-
equilibrium nature of the flow, non-negligible levels of surface slip velocity
were observed, with the highest magnitude seen in the shoulder region. In this
region, the slip velocity was on the order of $5\%$ of the freestream
velocity. Another region of high slip velocity was observed on the windward
side of the afterbody, although the magnitude was not to the extent of the
slip velocity at the shoulder.

The surface streamlines show a pronounced separation region on the leeward
side of the afterbody. The separation region was largely convex and circular
with a size on the order of half of the diameter of the capsule. This surface
flow behavior was notably different than the predictions of Wright et al. [4],
although these discrepancies are expected due to the different governing
equations and wall boundary conditions (e.g., the present work does not assume
no slip walls) as well as the differences in the physics being modeled.
However, the location and size of the separation region was somewhat similar
between the present work and the work of Wright et al. [4].

Figure 5: Surface oil flow streamlines overlaid on contours of surface slip
velocity magnitude for the high-altitude flow conditions.

####  4.1.3 Wake structure

To observe the flow structure in the wake region, the velocity field on the
pitch plane was first visualized through surface streamlines, shown in Fig. 6.
The recompression shocks are evident in the streamlines, with two distinct
regions of high velocity gradients aft of the capsule surface. Furthermore, a
complex wake structure can be observed starting with the separation of the
flow at the shoulder. This then develops into a small recirculation region on
the leeward side of the capsule afterbody. The interaction between this
recirculating region and the recompression shocks causes a complex rollup in
the wake with a distinct critical point.

Figure 6: Surface streamlines of velocity in the wake region on the pitch
plane ($z=0$) for the high-altitude flow conditions.

The three-dimensional behavior of the flow was then visualized through flow
streamlines and velocity cross-sections, shown in Fig. 7. The streamlines and
velocity cross-sections evidently show the separation region aft of the
capsule surface on the leeward side. Furthermore, the rollup of the
recirculation region is also shown clearly in the streamlines. The sonic
surface, the isosurface where $M=1$, is shown in blue. This sonic surface
encapsulates a majority of the wake region and shows similarity to the results
of Wright et al. [4].

Figure 7: Three-dimensional visualization of velocity streamlines and cross-
sections (placed 1, 2, and 3 capsule diameters downstream, respectively)
colored by velocity magnitude for the high-altitude flow conditions. Sonic
surface shown by light blue isosurface.

####  4.1.4 Shock structure & non-equilibrium flow behavior

Finally, the shock structure and effects of flow non-equilibrium were
analyzed. To visualize shocks in the flow, a Schlieren-type representation of
the density gradient norm is shown in Fig. 8, visaulized on the pitch plane
and capsule surface. The high density gradients in the bow shock are clearly
visible as well as in the recompression region. Furthermore, it can be seen
that even with a high-order scheme on an unstructured mesh, the numerical
solution is still very well behaved around strong discontinuities, showcasing
the ability of the approach in directly resolving shock structures without
shock capturing schemes.

(a) Schlieren

(b) Local Knudsen number

Figure 8: Schlieren-type visualization of the density gradient norm (left) and
contours of the density gradient-length local Knudsen number (right) on the
pitch plane ($z=0$) and the capsule surface for the high-altitude flow
conditions.

The behavior of the flow across the bow shock is shown more explicitly in the
flow profiles in front of the capsule nose in Fig. 9. Prior to the compression
in the flow near the capsule surface, the behavior of the shock is similar to
a static one-dimensional normal shock. It can be seen that the large jumps in
density and pressure are well-resolved by the numerical scheme without any
spurious oscillations. Furthermore, from the validation work performed in
Dzanic et al. [2], it is expected that the scheme predicts the physically
correct shock thickness and behavior. As can be seen in Fig. 9, the shock
thickness is of $\mathcal{O}(10^{-2}D)$, which is consistent with an upstream
mean free path of $\mathcal{O}(10^{-3}D)$ and falls well within the resolving
capability of the given mesh and approximation order.

$-0.05$$-0.10$$-0.15$$0$$2$$4$$6$$8$$10$$x/D$ $\rho/\rho_{\infty}$

(a) Density

$-0.05$$-0.10$$-0.15$$0$$0.2$$0.4$$0.6$$0.8$$1$$x/D$ $\|U\|/\|U\|_{\infty}$

(b) Velocity

$-0.05$$-0.10$$-0.15$$0$$200$$400$$600$$800$$x/D$ $P/P_{\infty}$

(c) Pressure

Figure 9:  Profiles of normalized density (left), velocity magnitude (middle),
and pressure (right) across the bow shock in front of the capsule nose ($y=0$,
$z=0$) for the high-altitude flow conditions.

To quantify the degree of non-equilibrium in the flow, the density gradient-
length local Knudsen number was computed, given as

$Kn_{\text{GLL}}=\lambda\frac{\|\nabla\rho\|_{2}}{\rho},$ (22)

where $\lambda$ is the local mean free path defined as

$\lambda=\frac{\mu(\theta)}{\rho}\sqrt{\frac{\pi}{2\theta}.}$ (23)

From Boyd et al. [24], it is expected that continuum breakdown begins when
$Kn_{\text{GLL}}>0.05$. The contours of the local Knudsen number are also
shown on the pitch plane and capsule surface in Fig. 8. It can be seen that
this breakdown condition is exceeded in a large portion of the flow. In
particular, the bow shock shows a very large local Knudsen number, with a
maximum of nearly $0.5$. Furthermore, the capsule surface on the leeward side
near the shoulder also shows flow regions which exceed this condition as well
as the recompression shock region. These results are in line with the
observations of Wright et al. [4] which may explain the discrepancies in the
results obtained by Wright et al. [4] using the continuum fluid dynamics
equations in comparison to the experimental data.

###  4.2 Low-altitude conditions

####  4.2.1 Macroscopic flow cross-sections

The results of the unsteady, low-altitude, low Mach simulation were then
analyzed with respect to the _instantaneous_ macroscopic flow state on the
symmetry plane cross-section ($z=0$). The contours of instantaneous density,
pressure, velocity magnitude, and temperature are shown in Fig. 10 normalized
by the freestream values. It can be seen that the structures in the flow are
significantly more complex than the high-altitude case, largely as a result of
the higher Reynolds number driving an unsteady turbulent wake flow. A similar
separation region could be seen aft of the capsule, but the recompression
shocks were noticeably sharper and showed more complex features which may also
be attributed to the higher Reynolds number and, by extension, lower Knudsen
number in the flow. The contours of density and pressure show a transitional
boundary layer aft of the shoulder which develop into unsteady flow features
in the wake. Similarly to the high-altitude case, the bow shock was evident
with large jumps in the macroscopic flow variables across the shock.
Furthermore, even when using a mesh not aligned with the shock, the scheme was
still able to robustly and accurately resolve the strong shock without any
explicit shock capturing approach. This showcases the robustness of the
approach, particularly for complex flows and flows which may span many flow
regimes as the method is not noticeably sensitive to mesh quality/design which
is typically a drawback for high-order numerical schemes.

(a) Density

(b) Pressure

(c) Velocity

(d) Temperature

Figure 10: Instantaneous contours of macroscopic normalized density (top
left), pressure (top right), velocity magnitude (bottom left), and temperature
(bottom right) on the pitch plane ($z=0$) and the capsule surface for the low-
altitude flow conditions.

####  4.2.2 Wake structure

The flow structure in the wake region was similarly analyzed by visualizing
the _instantaneous_ surface streamlines on the pitch plane, shown in Fig. 11.
It can be seen that the unsteady turbulent wake flow results in significantly
more complex behavior in the flow streamlines. A noticeably larger separation
region on the leeward side was observed than the high-altitude case, likely as
a result of the lower wall slip velocity in the near-continuum regime, but the
flow separated at a similar location aft of the shoulder. The streamlines show
very chaotic flow in this wake, with large flow gradients around the
recompression shocks and the distinct rollup of vortical flow features. On the
windward side of the capsule, the flow streamlines were primarily smooth, but
the near-wall region showed signs of transition to turbulence. These results
highlight the ability of the approach in resolving complex unsteady flow
features which high-order numerical schemes are particularly well-suited for.

Figure 11: Instantaneous surface streamlines of velocity in the wake region on
the pitch plane ($z=0$) for the low-altitude flow conditions.

##  5 Conclusion

We present the results of high-order numerical experiments of hypersonic re-
entry vehicles computed by performing three-dimensional molecular gas dynamics
simulations through directly solving the six-dimensional Boltzmann –BGK
equation. The underlying numerical approach combines a highly-efficient high-
order unstructured discontinuous spectral element method for the spatial
domain with a discretely-conservative discrete velocity model for the velocity
domain. This approach was used to simulate the AS-202 Apollo capsule at
realistic re-entry conditions in both the high-altitude regime at a Mach
number of 22.7 and a Reynolds number of 43,000 and the low-altitude regime at
a Mach number of 8 and a Reynolds number of 550,000. These results demonstrate
a proof of feasibility of the approach on modern computing architectures,
showcasing the ability to directly resolve strong shocks, non-equilibrium flow
phenomena, and unsteady turbulent flow features. Future work will focus on the
implementation of various models to better approximate the underlying physics
in high-speed re-entry flows.

## Acknowledgments

This work was partially performed under the auspices of the U.S. Department of
Energy by Lawrence Livermore National Laboratory under contract
DE–AC52–07NA27344. Release number LLNL–CONF–858000. The authors would like to
acknowledge the computational resources provided by the Princeton Institute
for Computational Science and Engineering and by Professor Freddie Witherden
at Texas A&M University.

## References

  * Schneider [2006] Schneider, S. P., “Laminar-Turbulent Transition on Reentry Capsules and Planetary Probes,” _Journal of Spacecraft and Rockets_ , Vol. 43, No. 6, 2006, pp. 1153–1173.  10.2514/1.22594.
  * Dzanic et al. [2023a] Dzanic, T., Witherden, F., and Martinelli, L., “A positivity-preserving and conservative high-order flux reconstruction method for the polyatomic Boltzmann–BGK equation,” _Journal of Computational Physics_ , Vol. 486, 2023a, p. 112146.  10.1016/j.jcp.2023.112146.
  * Dzanic et al. [2023b] Dzanic, T., Witherden, F. D., and Martinelli, L., “Validation of wall boundary conditions for simulating complex fluid flows via the Boltzmann–BGK equation: Momentum transport and skin friction,” _arXiv_ , 2023b.
  * Wright et al. [2006] Wright, M. J., Prabhu, D. K., and Martinez, E. R., “Analysis of Apollo Command Module Afterbody Heating Part I: AS-202,” _Journal of Thermophysics and Heat Transfer_ , Vol. 20, No. 1, 2006, pp. 16–30.  10.2514/1.15873.
  * Yuan and Zhong [2020] Yuan, R., and Zhong, C., “A conservative implicit scheme for steady state solutions of diatomic gas flow in all flow regimes,” _Computer Physics Communications_ , Vol. 247, 2020, p. 106972.  10.1016/j.cpc.2019.106972.
  * Zhang et al. [2023] Zhang, Y., Zeng, J., Yuan, R., Liu, W., Li, Q., and Wu, L., “Efficient parallel solver for high-speed rarefied gas flow using GSIS,” , 2023.  10.48550/ARXIV.2310.18916.
  * Chen et al. [2020] Chen, Y., Zhu, Y., and Xu, K., “A three-dimensional unified gas-kinetic wave-particle solver for flow computation in all regimes,” _Physics of Fluids_ , Vol. 32, No. 9, 2020.  10.1063/5.0021199.
  * Cercignani [1988] Cercignani, C., _The Boltzmann Equation and Its Applications_ , Springer New York, 1988.  10.1007/978-1-4612-1039-9.
  * Bhatnagar et al. [1954] Bhatnagar, P. L., Gross, E. P., and Krook, M., “A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems,” _Physical Review_ , Vol. 94, No. 3, 1954, pp. 511–525.  10.1103/physrev.94.511.
  * Holway [1966] Holway, L. H., “New Statistical Models for Kinetic Theory: Methods of Construction,” _Physics of Fluids_ , Vol. 9, No. 9, 1966, p. 1658.  10.1063/1.1761920.
  * Baranger et al. [2020] Baranger, C., Dauvois, Y., Marois, G., Mathé, J., Mathiaud, J., and Mieussens, L., “A BGK model for high temperature rarefied gas flows,” _European Journal of Mechanics - B/Fluids_ , Vol. 80, 2020, pp. 1–12.  10.1016/j.euromechflu.2019.11.006.
  * Sutherland [1893] Sutherland, W., “LII. The viscosity of gases and molecular force,” _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , Vol. 36, No. 223, 1893, p. 507–531.  10.1080/14786449308620508.
  * Huynh [2007] Huynh, H. T., “A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods,” _18th AIAA Computational Fluid Dynamics Conference_ , American Institute of Aeronautics and Astronautics, 2007.  10.2514/6.2007-4079.
  * Hesthaven and Warburton [2008] Hesthaven, J. S., and Warburton, T., _Nodal Discontinuous Galerkin Methods_ , Springer New York, 2008.  10.1007/978-0-387-72067-8.
  * Cai and Torrilhon [2019] Cai, Z., and Torrilhon, M., “On the Holway-Weiss debate: Convergence of the Grad-moment-expansion in kinetic gas theory,” _Physics of Fluids_ , Vol. 31, No. 12, 2019.  10.1063/1.5127114.
  * Loek Van Heyningen [2021] Loek Van Heyningen, R., “Discontinuous Galerkin solutions of the Boltzmann equation: spectral collocation and moment methods,” Master’s thesis, Massachusetts Institute of Technologyy, 2021.
  * Mieussens [2000a] Mieussens, L., “Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamic,” _Mathematical Models and Methods in Applied Sciences_ , Vol. 10, No. 08, 2000a, pp. 1121–1149.  10.1142/s0218202500000562.
  * Trefethen and Weideman [2014] Trefethen, L. N., and Weideman, J. A. C., “The Exponentially Convergent Trapezoidal Rule,” _SIAM Review_ , Vol. 56, No. 3, 2014, pp. 385–458.  10.1137/130932132.
  * Mieussens [2000b] Mieussens, L., “Discrete-Velocity Models and Numerical Schemes for the Boltzmann-BGK Equation in Plane and Axisymmetric Geometries,” _Journal of Computational Physics_ , Vol. 162, No. 2, 2000b, pp. 429–466.  10.1006/jcph.2000.6548.
  * Hillje [1967] Hillje, E. R., “Entry flight aerodynamics from Apollo mission AS-202,” Tech. Rep. 19670027745, Manned Spacecraft Center, NASA, Houston, Texas, October 1967.
  * Harten et al. [1983] Harten, A., Lax, P. D., and Leer, B. v., “On Upstream Differencing and Godunov-Type Schemes for Hyperbolic Conservation Laws,” _SIAM Review_ , Vol. 25, No. 1, 1983, p. 35–61.  10.1137/1025002.
  * Dzanic and Witherden [2022] Dzanic, T., and Witherden, F., “Positivity-preserving entropy-based adaptive filtering for discontinuous spectral element methods,” _Journal of Computational Physics_ , Vol. 468, 2022, p. 111501.  10.1016/j.jcp.2022.111501.
  * Witherden et al. [2014] Witherden, F., Farrington, A., and Vincent, P., “PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach,” _Computer Physics Communications_ , Vol. 185, No. 11, 2014, pp. 3028–3040.  10.1016/j.cpc.2014.07.011.
  * Boyd et al. [1995] Boyd, I. D., Chen, G., and Candler, G. V., “Predicting failure of the continuum fluid equations in transitional hypersonic flows,” _Physics of Fluids_ , Vol. 7, No. 1, 1995, p. 210–219.  10.1063/1.868720.


11institutetext: University Grenoble-Alps, INRIA, Grenoble, France 11email:
<EMAIL_ADDRESS>
# Can Authoritative Governments Abuse the Right to Access?

Cédric Lauradoux 11

###### Abstract

The right to access is a great tool provided by the GDPR to empower data
subjects with their data. However, it needs to be implemented properly
otherwise it could turn subject access requests against the subjects privacy.
Indeed, recent works have shown that it is possible to abuse the right to
access using impersonation attacks. We propose to extend those impersonation
attacks by considering that the adversary has an access to governmental
resources. In this case, the adversary can forge official documents or exploit
copy of them. Our attack affects more people than one may expect. To defeat
the attacks from this kind of adversary, several solutions are available like
multi-factors or proof of aliveness. Our attacks highlight the need for strong
procedures to authenticate subject access requests.

###### Keywords:

Subject Access Right Authentication Impersonation Forgery.

##  1 Introduction

In 2011, Max Schrems, an Austrian citizen, sent a request to the US company
Facebook111http://europe-v-facebook.org to access all the data held by the
company. He received a huge PDF file (1200 pages) describing the data Facebook
held on him. The analysis of the data shows that Facebook was not respecting
European data protection laws. In 2018, David Caroll submitted a subject
access request to UK company Cambridge Analytica. The company first refused to
comply with Caroll’s request, and it escalated until the UK data protection
authority [13] asked the company to comply with Caroll’s access request. In
the end, Cambridge Analytica’s actions during 2016’s US election were exposed.
In both cases, Facebook and Cambridge Analytica had to provide data to
individuals in order to respect their right to access. In Europe, this right
is defined in Article 15 of the GDPR. The right to access data is considered
as the root for all the other rights defined by the GDPR. Once, a data subject
has accessed the data collected by a data controller, he/she can exercise the
other rights (to object, to data portability, to be forgotten, etc.).

The right to access needs to be implemented carefully to avoid data breaches.
Indeed, an adversary can send a subject access request while impersonating a
legitimate subject to a data controller. If the data controller does not
authenticate the access request properly, the adversary can obtain the data by
abusing the right to access. This threat has to be considered seriously as
data protections officers (DPOs) have already reported in [14] suspicious
access requests. Researchers [6, 8, 9, 14, 15, 16] have conducted studies to
test if the right to access was implemented properly. They have used different
strategies [8, 9, 15, 16] to fool DPOs and succeeded to steal personal data in
numerous cases. In response to those issues, the European Data Protection
Board (EDPB) has recently published some guidelines [12] on the right of
access to clarify its implementation and to improve the situation. Still many
things have to be done to fix how subject access requests are handled by data
controllers.

In this work, we focus on data controllers who authenticate a subject access
request using an official document like an ID card or a passport. This
authentication is often encountered in practice despite the fact that it can
be irrelevant as pointed out in [6, 12]. We alert that there are more
elaborated attacks possible if the adversary benefits from the support of a
government. In this case, the adversary can forge official document and
impersonate any data subjects to forge subject access request. It allows
authoritative or corrupted states to abuse the right to access to spy on
dissidents who use online services eligible to subject access request. This
kind of attacks needs to be avoided because it could clearly harm the
reputation of the data protection regulations like the GDPR.

Official document forgery shows the limit of document-based authentication.
Additional measures needs to be implemented by data controllers to prevent ID
forgery. Multiple-factor authentication (MFA) is an interesting solution to
thwart forgery attacks. It requires that the data subject and the data
controller agree on an additional secure communication (an electronic mail
address or a phone number).

The paper is organized as follows. Section 2 discusses how data controllers
authenticate subject access requests. We describe the benefit of implementing
the right to access using privacy dashboard or by DPOs. We review the existing
attacks against the right to access in Section 3. We introduce our new attack
scenarios in Section 4. Solutions are discussed in Section 5.

##  2 How are subject access requests authenticated?

The right to access one’s personal data (Article 15 of GDPR [10]) via Subject
Access Requests (SAR) is fundamental to ensure data protection. It makes the
processing of personal data by organizations transparent and accountable
towards data subjects, whereby these are made aware of, and can verify the
lawfulness of the processing of their data (Recital 63 of GDPR [10]).

This right to access by individuals complements the mandate of the data
protection authorities to monitor and enforce the application of the GDPR
(Article 57.1.a of GDPR [10]). The right to access also enables the exercise
of all the other rights (erasure, rectification, restriction, etc). When a
data controller receives an access request, it needs to perform two types of
verification. A verification of the _eligibility_ needs to be done. A data
controller needs to verify if the subject sending a request is concerned by a
data protection regulation which includes a right to access. A verification of
the _legitimacy_ needs also to be done by the data controller. Verifying the
legitimacy has two purposes: i) establish the subject’s identity, and ii)
check that the subject’s identity matches the requested data. Verifying both
_eligibility_ and _legitimacy_ can be obtained by executing an appropriate
authentication protocol between the data subject sending the access request
and the data controller.

If no verification is implemented, the data of a subject can be collected by
anyone impersonating the subject to the data controller. This is particularly
demonstrated by Bufalieri et al. [8]. They submitted 334 subject access
requests to different data controllers, and 58 organizations provided data
without any verification of the subject sending the request. It confirms that
the right to access without authentication is harmful to privacy. We discuss
next how data controllers are actually handling subject access requests.

###  2.1 Dashboard versus DPO

Two methods to handle subject access requests (see Fig. 1) have been
implemented by data controllers. In the first method, a privacy dashboard has
been implemented by the data controller. The subject authenticates to the
controller’s website and then can obtain his/her data or submit securely an
access request to the controller’s DPO. This solution has been chosen by
Google or Facebook. They have created dashboard on their websites to let data
subjects access their data. To access the dashboard, the data subject needs
first to authenticate to the data controller websites. An extra verification
may be required to obtain the data from the data controller as part as a two-
factor authentication. Privacy dashboards are very convenient for the data
subjects because the creation and the processing of their request is
automatic. Google Takeout (https://takeout.google.com) is a good example of
dashboard which implements two-factor authentication. More details on privacy
dashboards can be found in [18].

In the second method, all the subject access requests are handled by the
controller’s DPO. The DPO is in charge of verifying both eligibility and
legitimacy of the subject access requests. In this method, a first issue is to
setup a secure communication channel to let the subject submits his/her
request to the DPO. It can be done through a web page or by sending and
electronic mail222Other more classical methods like postal mail are also
available but they are not considered in this work.. In any case, the
communication channel between the subject and the DPO needs to be secure
(using end-to-end encryption for instance). Once a secure channel has been
setup, the DPO can interact with the subject to verify his/her request. It is
important to notice that third parties can help the subject to contact the
DPOs. A data protection authority can be contacted in the case for an indirect
access. It can be also a third party which helps subjects to create their
subject access requests like Tapmydata333https://tapmydata.com.

Figure 1: Current options for the implementation of the right to access. Data
---
Subject
DashboardDPO Personal
---
Data
Third
---
Party
Data controller

In Table 1, we compare the main characteristics of the different
implementations one can come across. The security of the authentication scheme
used in a dashboard is based on a secret shared (i.e. a password) between the
data controller and the data subject. It is human-to-computer security
mechanism based on the following assumptions: (i) the data subject does not
share the secret with other individuals, (ii) the secret is not re-used to
other data controllers, (iii) the secret is difficult to guess and (iv) the
authentication is performed over secure channel. This technology is tried-and-
tested even if it has some limitations [7].

Table 1: Comparison of the properties of the different implementations of the right to access. | Dashboard | DPO
---|---|---
Creation of SAR | Automatic | Manual
Processing of SAR | Automatic/Manual | Manual
Security model | Password (+MFA) | Choice of the DPO

The authentication scheme implemented by DPO is a human-to-human interaction.
It is based on the assumptions that: (i) the communication channels used
between the subject and the DPO are secure, (ii) the data subject is the only
individual able to prove that he/she has created his/her request and (iii) the
DPO is trained to verify requests. The verification is a choice made by the
DPO. This choice is critical to authenticate subject access requests. DPOs
often ask the subject sending a request to provide a copy of an official
identification document. This method of verification has been criticized in
[12].

##  3 Threat Model and Known Attacks

A general threat model for the right to access was proposed by Boniface et al.
[6]. This model considers both malicious subjects and malicious controllers,
whereby three different categories of threats were defined: i) denial of
access, ii) privacy invasion and iii) data breach due to impersonation. We
briefly review denial of access and privacy invasion and then focus on data
breach.

#### Denial of access.

This situation corresponds to cases in which data controllers decline to
provide data to data subjects; one alleged reason refers to scenarios wherein
they are not in a position to verify the identity of subjects.It also occurs
often with companies collecting and processing cookies [6, 11], but also with
those processing IP addresses [4]. Denial of access was also observed in [5,
6, 11, 19].

#### Privacy invasion.

This attack is defined in [6] as a case in which the data controller is
honest-but-curious but requests additional information to authenticate the
SAR. These additional requested information can increase the exposure of the
subject to the controller and potentially violates the _minimization and
proportionality principles_ prescribing that personal data shall be is
adequate, relevant and limited to what is necessary in relation to the
purposes for which they are processed (Article 5 (1)(c) of the GDPR). For
example, the controller should not be obliged to acquire additional
information in order to identify the data subject for the sole purpose of
being complaint (Recital 57 of the GDPR).

#### Data breach.

Data breach refers to the case of a data controller which provides the data of
a given subject A to another subject B. A data breach can be the result of an
error or it can be the result of an impersonation attack. Example of
impersonation attacks against the right to access are given in [8, 9, 15, 16].
These studies have demonstrated that it is possible to use the right of access
to steal data by impersonating a data subject. Recent studies [8, 9, 15, 16]
show that data controllers are not ready to process correctly subject access
requests because they still use weak authentication procedures vulnerable to
impersonation attacks. Table 2 summarizes the different methods used to
impersonate a subject and to abuse data controllers with illegitimate subject
access requests. It is important to notice that the attacks were successful
against data controllers which rely their DPOs to authenticate subject access
requests.

Authentication method | Attacks | Target
---|---|---
Copy of an ID | Social engineering [16] | DPO
Falsification [15] | DPO
Email confirmation | Email spoofing [8, 9, 15] | DPO
Bills | Falsification [16] | DPO
Personal question | Social engineering [16] | DPO
Table 2: Known impersonation attacks abusing the right to access and their
implementation.

##  4 Advanced Forgery Attacks

We demonstrate in this section that verification of subject access request
based on official identification document can not resist to an adversary which
is backup by a government. Indeed, we focus on DPOs who verify subject access
requests by asking a copy of official identification document. Despite being
criticized by EDPB in its guidelines [12], this method of verification is
still considered valid depending on the information known by the data
controller on the data subject. We assume that if the adversary is able to
provide an official identification document then he/she can fool DPOs and
therefore abuse the right to access. Other authors (like [14] for instance)
have made the same conclusion but they have not described scenarios and their
impact as it is done in our work. Our attack is foreseen in three steps (as
depicted in Figure 2).

  1. The malicious state (represented by Charlie) needs first to identify the websites and services used by the target (Alice) located in the EU. He obtains the contact information of all the data controllers collecting and processing the data of their targets (including Bob).

  2. Charlie produces official documents (ID card, passport…) which match the identity of Alice. Charlie creates documents for a perfect namesake with the same first, last name, birth date, gender, etc. These documents are transmitted through a proxy (Eve).

  3. Eve submits the subject access request to Bob and to any other controllers. Eve can provide valid government-issued documents whenever they are requested by data controllers.

  4. Eve (and later Charlie) obtains Alice’s data from Bob.

CharlieEveAliceBob① Identification ④ Alice’s
---
data
② Forged
---
Alice’s ID
Alice’s
---
ID
③ SAR using
---
Alice’s ID
④ Alice’s
---
data
Figure 2: Impersonation attack by an authoritative state (Charlie) against a
dissident (Alice) using a proxy (Eve) to contact a data controller (Bob).

The first step is not mandatory. Indeed, Charlie and Eve can adopt a guess-
and-determine strategy and submit subject access requests to many data
controllers without knowing whether they collect and process personal data
from Alice. If the data controller (Bob) answers, Charlie and Eve obtain
Alice’s data. However, this strategy is not particularly discreet and some
data controllers may alert Alice that they have received strange requests.

The authoritative state may not necessarily need to use Eve as a proxy. Eve
can be useful since she can be located in the same country as the targeted
dissident: the request submitted by Eve would mimic the request that the
target could have submitted.

Tools like the application Tapmydata444https://tapmydata.com/ can simplify the
task of the adversary. Tapmydata can be used to submit automatically requests
to data controllers. It aims at simplifying the creation and the submission of
subject access requests. As such, an adversary can register with an email
address and confirm her identity by uploading a scanned of his/her forged
passport. The application knows how to contact 1400 data controllers and takes
care of all the administrative details of the request. When a data controller
replies to a request, the data is stored in a docker. In our context, it
amplifies our attack.

A similar scenario can occur if an administrative who works in a department
delivering official identification documents is corrupted. Figure 3
illustrates this case wherein Alice is a data subject using an online social
network maintained by Bob (and located in the EU). Eve wants to obtain
information on Alice. Eve knows that Charlie can be corrupted. Eve bribes
Charlie to obtain the genuine identification document bearing Alice’s
information. Eve submits a SAR to Bob impersonating Alice. Eve can provide all
the identification documents.

CharlieEveAliceBobBribe Duplicata of
---
Alice’s ID
Alice’s
---
ID
SAR using
---
Alice’s ID
Alice’s
---
data
Figure 3: Impersonation attack when an administrative Charlie is corrupted by
Eve.

###  4.1 Impact Assessment

One may consider that this attack is unlikely to occur and that the adversary
model is too strong. Governments have other legal means to access data from
their citizens. Our attack scenario makes sens if we consider the case of
dissidents and authoritative states. Let us assume that some dissidents are
using online services located in Europe. The authoritative state they are
against can use our forgery attack to abuse the right to access and spy on
them. This scenario would demonstrate that the GDPR can be abused for
surveillance purpose. Creating on-demand official identification documents and
sending subject access requests are rather inexpensive. If the impersonation
succeeds, there are two types of consequences:

  * •

Direct consequences: a state can destabilize, denigrate or blackmail
dissidents online.

  * •

Indirect consequences: after a first impersonation attack, a state can further
impersonate the target to exercise other rights (depending on its goal), for
example:

    * –

rectification (Article 16) to modify data or content published by the target.
It can harm the reputation of the target or it can be used to denigrate
his/her speech.

    * –

erasure (Article 17) as an effective censorship method by removing all the
data of the dissident from a website (his/her blog for instance). This is a
form of denial of service attack.

The scope of the GDPR increases the impact of our attack to more categories of
people. First, we have dissidents who are using European online service. The
authoritative state they are against can forge official identification
documents at will and submit subject access request to data controllers which
are in EU. Second, European citizen who have visited an authoritative state.
The authoritative state can keep a copy of their passport when they cross the
border during the customs check. The authoritative state can submit subject
access requests to data controller anywhere who collect and process the data
of the target. Third, Foreigners (Non European) who are using European online
service who have visited an authoritative state. It demonstrates that an
authoritative state can abuse the right to access and exploit the scope of the
GDPR to reach as many people as possible.

###  4.2 Ethical considerations

We have not implemented our attack due to ethical considerations. To his
knowledge, the author does not collaborate with authoritative states. We have
neither attempted to pay a bribe to an official. Even if our work is
theoretical, the attack described is serious and could affect many data
subjects.

##  5 Countermeasures

The verification made by the DPO to verify a request need to be relevant to
the controller’s knowledge of the data subject. Actually, document-based
authentication can establish the identity of the subject securely even against
governmental adversaries if three conditions are met:

  1. 1.

The data controller knows the identity of the subject and an identifying
element (picture of the face or another biometrics).

  2. 2.

The data controller and the subject sending the request are met physically.

  3. 3.

The DPO is trained to verify official documents (passport or ID card).

If often fails because at least one of these conditions is not satisfied. If
condition 3 is verified it can easily. The DPO needs to be trained to detect
fake identification documents and document tampering. If the other conditions
(1 and 2) are not verified, the patch is more costly.

#### Remote identity Proofing –

Let us assume that condition 2 cannot be satisfied, i.e. the DPO and the
subject cannot met physically like in a police control. The DPO needs to use a
remote identity proofing scheme as described in [2, 3]. The verification can
be done using a video call or pictures. It relies on a liveness check and
authenticity check (see [2] for more details). it is important that the data
controller knows both the identity and one biometrics of the subject before
the verification. Otherwise a governmental adversary can succeed a forgery
attack. If Bob knows only the identity of Alice, Charlie can forge a passport
in the name of Alice with the picture and the biometrics of Eve. Eve can then
impersonate Alice to Bob. Another approach for remote identity proofing is to
rely on secure public electronic identification [1]. Such systems are
considered for deployment in the European Union. They could be used
straightforward to authenticate subject access requests sent by European
residents. However, data controllers would need to find another method to
authenticate requests from data subjects which are outside the European Union.

#### Multi-factor Authentication –

MFA [17] is a solution to use if condition 1 cannot be satisfied. The DPO
needs to consider all the information known on the subject to authenticate the
subject’s request. It can be a confirmation code sent on an email address or
an SMS sent to a phone number for instance. It is important to notice that
using a factor unknown by the controller before the request is not relevant
and it cannot be trusted.

If the data controller needs strong authentication to interact with the
subject, then the authentication mechanism can be re-used to verify subject
access requests. Both data processing and exercise of the subjects right
aligned. An important question rises if the data controller does not need
strong authentication to process subjects data. Does the data controller need
to implement strong authentication to verify subject access requests? If we
consider the use of multi-factor authentication, a data controller can ask the
subject to provide a phone number when he/she registers for the sole purpose
of authenticating later subject access requests. Recital 64 of the GDPR states
that A controller should not retain personal data for the sole purpose of
being able to react to potential requests. The data controller needs to
implement the data minimization principle during the verification. There is a
conflict between the application of the minimization principle and the
security of subject access requests. Asking additional data to the subject
seems necessary to have a safe exercise of the right to access.

##  6 Conclusion

Governments have opportunities to abuse the right to access when a data
controller authenticates subject access requests using only official
identification documents. This form of authentication is particularly weak
against forgery attacks. Dissidents who use online services proposed by
European companies can be targeted by authoritative governments they oppose.
The scope of the GDPR makes our attack reach more people than expected.
Forgery attacks can be even obtained by corruption. If such attacks occur,
these attacks would not only harm data subjects but they would also harm the
trust on the personal data regulations like the GDPR.

Finding a solution is not simple because data controllers need to find the
good trade-off between security and privacy. Implementing strong
authentication using multiple factors mitigates the risks related to handling
subject access requests. But the data controller may need to obtain more
information like the subject phone number or an email address to implement
multi-factor authentication. Such information may not need necessary for the
controller’s processing. It is hard to satisfy both the data minimization
principle (Recital 64 of the GDPR) and a useful and secure right to access.
Data protection authorities and board need to clarify if minimization or right
to access prevails.

## References

  * [1] Privacy Features of European eID Card Specifications. Tech. rep., ENISA (January 2009)
  * [2] Remote ID Proofing. Tech. rep., ENISA (March 2021)
  * [3] Remote ID Proofing: Attacks & Countermeasures. Tech. rep., ENISA (January 2022)
  * [4] Adhatarao, S., Lauradoux, C., Santos, C.: IP-based Subject Access Requests Denied (2021)
  * [5] Ausloos, J., Dewitte, P.: Shattering one-way mirrors – data subject access rights in practice. International Data Privacy Law 8(1), 4–28 (2018)
  * [6] Boniface, C., Fouad, I., Bielova, N., Lauradoux, C., Santos, C.: Security Analysis of Subject Access Request Procedures How to authenticate data subjects safely when they request for their data. In: Annual Privacy Forum. Lecture Notes in Computer Science, vol. 11498, pp. 1–20. Springer, Rome, Italy (Jun 2019)
  * [7] Bonneau, J., Herley, C., van Oorschot, P.C., Stajano, F.: Passwords and the evolution of imperfect authentication. Commun. ACM 58(7), 78–87 (2015)
  * [8] Bufalieri, L., Morgia, M.L., Mei, A., Stefa, J.: GDPR: When the Right to Access Personal Data Becomes a Threat. In: 2020 IEEE International Conference on Web Services, ICWS 2020. pp. 75–83. IEEE, Beijing, China (October 2020)
  * [9] Cagnazzo, M., Holz, T., Pohlmann, N.: GDPiRated - Stealing Personal Information On- and Offline. In: Computer Security - ESORICS 2019 - 24th European Symposium on Research in Computer Security. Lecture Notes in Computer Science, vol. 11736, pp. 367–386. Springer, Luxembourg (September 2019)
  * [10] Council of European Union: Council regulation (EU) no 2016/679 (2016), https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:32016R0679
  * [11] Degeling, M., Utz, C., Lentzsch, C., Hosseini, H., Schaub, F., Holz, T.: We Value Your Privacy…Now Take Some Cookies: Measuring the GDPR’s Impact on Web Privacy. In: 26th Annual Network and Distributed System Security Symposium, NDSS 2019. p. https://arxiv.org/abs/1808.05096. The Internet Society, San Diego, California, USA (February 2019)
  * [12] European Data Protection Board: Guidelines 01/2022 on data subject rights - Right of access. Tech. rep. (January 2022), \url{https://edpb.europa.eu/our-work-tools/documents/public-consultations/2022/guidelines-012022-data-subject-rights-right_en}
  * [13] Information Commissioner’s Office: Enforcement Notice. Tech. Rep. 2258812, ICO (2018), https://ico.org.uk/media/action-weve-taken/enforcement-notices/2258812/en-scl-elections-20180504.pdf
  * [14] Martino, M.D., Meers, I., Quax, P., Andries, K., , Lamotte, W.: Revisiting Identification Issues in GDPR ‘Right Of Access’ Policies: A Technical and Longitudinal Analysis. In: Privacy Enhancing Technologies, PETS 2022. Lecture Notes in Computer Science, vol. 8555, p. To appear. Springer
  * [15] Martino, M.D., Robyns, P., Weyts, W., Quax, P., Lamotte, W., Andries, K.: Personal Information Leakage by Abusing the GDPR ’Right of Access’. In: Fourteenth Symposium on Usable Privacy and Security (SOUPS). p. 371–386. ACM, USENIX Association, Santa Clara, CA, USA (August 2019)
  * [16] Pavur, J.: GDPArrrrr: Using Privacy Laws to Steal Identities. In: Blackhat USA. p. https://arxiv.org/abs/1912.00731. Arxiv, Las Vegas, NV, USA (2019)
  * [17] van Tilborg, H.C.A., Jajodia, S. (eds.): Multifactor Authentication, pp. 808–808. Springer US (2011)
  * [18] Tolsdorf, J., Fischer, M., Iacono, L.L.: A Case Study on the Implementation of the Right of Access in Privacy Dashboards. In: Privacy Technologies and Policy - 9th Annual Privacy Forum, APF 2021. Lecture Notes in Computer Science, vol. 12703, pp. 23–46. Springer (June 2021)
  * [19] Urban, T., Tatang, D., Degeling, M., Holz, T., Pohlmann, N.: A Study on Subject Data Access in Online Advertising After the GDPR. In: Data Privacy Management, Cryptocurrencies and Blockchain Technology - ESORICS 2019 International Workshops, DPM 2019 and CBT 2019. Lecture Notes in Computer Science, vol. 11737, pp. 61–79. Springer, Luxembourg (September 2019)