# Probing interfacial effects with thermocapillary flows

M. Maza-Cuello  CNRS, Sciences et Ingénierie de la Matière Molle, ESPCI Paris,
PSL Research University, Sorbonne Université, 75005 Paris, France  C. Frétigny
CNRS, Sciences et Ingénierie de la Matière Molle, ESPCI Paris, PSL Research
University, Sorbonne Université, 75005 Paris, France  L. Talini  CNRS, Surface
du Verre et Interfaces, Saint-Gobain, 93300 Aubervilliers, France

###### Abstract

We report on the thinning of supported liquid films driven by
thermocapillarity. The liquids are oil films, of initial thicknesses of a few
tens of microns. A local and moderate heating of the glass substrate on which
they are spread on induces a thermocapillary flow, which allows the formation
of ultra-thin films. We show that, within a given time range, the thinning
dynamics of submicron thick films is governed by the thermocapillary stress.
It results in a simple dependency of the thickness profiles with time, which
is evidenced by the collapse of the data onto a master curve. The master curve
only depends on the liquid properties and on the thermal gradient, and allows
a measurement of the latter. As the films further thin down to thicknesses
within the range of molecular interactions, a deviation from the master curve
appears. Although all investigated oils are supposedly fully wetting glass,
the nature of the deviation differs between alkanes and silicone oils. We
attribute the observed behaviors to differing signs of disjoining pressures,
and this picture is confirmed by numerical resolution of the thin-film
equation. We suggest thermocapillary flows can be used to finely probe
molecular interactions.

##  I Introduction

The surface tension of a pure liquid decreases with increasing temperature and
a temperature gradient at its free surface results in a thermocapillary
stress, inducing a flow away from the hottest region. The migration of liquids
induced by a temperature gradient has been known and exploited for a long time
[1]. It is used for instance for the patterning of surfaces, in the context of
micromachining with molten metals [2, 3]. In this process, a strong laser beam
is shined on the surface of a body of metal, producing the melting of an upper
layer which is driven mainly by the surface tension gradient. In recent years,
a similar principle has been applied in the micropatterning of polymer melts
for lithography techniques [4]. In contrast, in other applications
thermocapillary flows may be undesired. For instance, they are known to cause
failure of lubrication in tribo-systems in which thermal gradients
spontaneously appear; therefore, the design of devices must be thought in
consequence [5]. From a more fundamental point of view, the behavior of a
supported liquid film submitted to thermal gradients has been studied in
different configurations. Theoretical works have focused on the dimpling of
the free surface of the film that results from a thermocapillary stress, and
to which both capillarity and gravity oppose [1, 6, 7, 8]. Experimentally,
controlled thermal gradients at the free surface can be created by locally
heating the liquid by radiative transfer [9], but more usually it is the
substrate they lie on that is heated. Focusing a laser at an absorbing
substrate provides an easy way to locally increase its temperature and create
a temperature gradient at a liquid interface by heat transfer [10, 11, 12,
13]. Whatever the effect inducing the thermocapillary stress, the latter
allows the formation and control of thin liquid films, in a region whose size
is given by the extent of the thermal gradient [6]. Within this region, in
which the films can reach submicron thicknesses even with a moderate
temperature rise, intermolecular forces have non negligible effects, which in
particular condition the rupture of the film. Thermocapillary-induced
dewetting of liquids on solid substrates has thus been investigated by
monitoring the thickness at the central hottest point [11, 13]. However, up to
date, poor attention has been paid to the dynamics of the whole thin film
region, in which intermolecular forces are large. Here, we show that
interfacial effects at very small scales can be probed through accurate
measurements of thickness variations in this region.

The paper is organized as follows. We describe the experimental setup and
methods in section II. A review of the hypotheses that are made leading to a
well-known thin film equation describing the evolution of a thin liquid film
under the action of thermocapillarity are discussed in section III. We further
present the experimental results in section IV. Finally, section V provides a
conclusion and outlook for future work.

##  II Experimental set-up and methods

###  II.1 Experimental set-up

The liquid films are spread by spin-coating on substrates that are rectangular
glass slides of thickness $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Prior
to deposition, the top side of the substrate is washed by spin-coating
isopropanol and acetone over it several times, letting it dry finally for a
couple of minutes before spreading the liquid. An electrical circuit is
printed on the bottom side of the slides, with a pattern consisting in two
squares connected by a line, as shown in Fig. 1a). The squares are connected
to a dc voltage generator and the resulting electrical current is measured.
The central line, of width $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$, acts
as a heat dissipator, generating a thermal field. We denote $O$ the centre of
the line and $Ox$ and $Oy$ the axes respectively perpendicular and parallel to
the line. We make the approximation that the induced thermal field does not
depend on $y$, which is justified in the vicinity of the centre of the line.
We define the thermal excess $\theta(x)=T(x)-T_{0}$ where $T_{0}$ is the room
temperature and $T(x)$ the thermal field. The dissipated power $P$ along the
line is a control parameter that sets the maximal temperature rise,
$\theta_{max}=\theta(0)$. In our experiments $P\sim$1\text{\,}\mathrm{mW}$$
and $\theta_{max}\sim$1\text{\,}\mathrm{K}$$. Infra-red measurements have
evidenced that the horizontal extent of the thermal gradient $w$ is a few
millimeters. The thermal gradient drives a flow in the liquid films away from
the central line, inducing a thinning of the films along that line. A typical
surface of a film once a very thin film has been formed at the central region
is shown in Fig. 1b). The profile was measured via an implementation of the
Free-surface synthetic Schlieren method [14], which we will not detail here
since we focus on measurements at smaller thicknesses, which are described in
the following. The region we focus on is that at the centre (where the heating
line lies), which has been emptied of liquid so that the thickness is below
the micrometre scale in an area spanning several millimetres. We emphasise
that, as expected, the measured profile only very weakly depends on $y$, which
justifies the 1D profile we will consider in the following. Experiments were
performed with different initial film thicknesses, $h_{0}$ and electrical
powers $P$. The initial thickness $h_{0}$ is measured by optical profilometry
just after spreading the films.

Figure 1: a) Top view of the glass substrate. On its bottom side, a
microprinted circuit connected to a voltage generator creates the thermal
gradient. b) Typical thickness surface $h(t,\mathbf{x})$ at an intermediate
stage of the dynamics, reconstructed surface via Free-surface synthetic
Schlieren method [14]. The thermocapillary effect promotes the apparition of a
very thin ($h<$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) layer that extends
over several millimetres.

###  II.2 Investigated liquids

The investigated liquids are silicone oil (ABCR AB112154), $n$-hexadecane
(Sigma-Aldrich) and a ramified isomer of hexadecane:
2,2,4,4,6,8,8-heptamethylnonane (HMN, Sigma-Aldrich). Their useful physical
properties are listed in table 1; in particular, they have very close values
of coefficient of surface tension variation,
$\gamma_{\theta}=-\partial_{T}\gamma|_{T=T_{0}}$. In contrast, we emphasise
that their melting temperatures are very different.

Since the temperature increase imposed to produce the thermocapillary effect
remains modest, the variations of viscosity and density are small. An
estimation from the thermal dependencies of these properties for the three
liquids investigated [15, 16, 17] gives a variation on viscosity and density
of around 2 % and 0.1 %, respectively, for a thermal rise of
$\theta_{max}=$1\text{\,}\mathrm{K}$$. Similarly, the variation of refractive
index of alkanes [18] and silicone oil [19] caused by this temperature change
is at most $0.01$ %. Therefore, the variations of viscosity, density and
refractive index produced by heating the liquid will be neglected in the
following.

|  $\eta$ ($\mathrm{mPa}\text{\,}\mathrm{s}$) |  $\rho$ ($\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$) |  $\gamma_{\theta}$ ($\mathrm{N}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$) |  $T_{\text{melt}}$ ($\mathrm{\SIUnitSymbolCelsius}$) | ref.
---|---|---|---|---|---
Silicone oil | $19$ | $950$ | $6.5\times 10^{-5}$ | $<-40$ |  [20],[21]
$n$-C16 | $2.2$ | $733$ | $8.7\times 10^{-5}$ | $18$ |  [22],[16],[23]
HMN | $3.2$ | $781$ | $8.7\times 10^{-5}$ | $-4$ |  [22],[17],[24]
Table 1: Investigated liquids with their main properties: viscosity, $\eta$,
density, $\rho$, thermal coefficient of the surface tension,
$\gamma_{\theta}=-\partial_{T}\gamma|_{T=T_{0}}$, and melting temperature,
$T_{\text{melt}}$.

###  II.3 Thickness measurements

The thickness of the film at the central region, where it will be below the
micrometre scale, is measured with 3-colour interferometry [25, 26]. The
measurement setup is sketched in Fig. 2a). The cell is shined from below using
a RGB-LED of peak wavelengths $\lambda_{R}=$635\text{\,}\mathrm{nm}$$,
$\lambda_{G}=$518\text{\,}\mathrm{nm}$$, and
$\lambda_{B}=$463\text{\,}\mathrm{nm}$$. Images are taken from above with a
CMOS monochrome camera. Images were taken typically each
$10\text{\,}\mathrm{min}$ for silicone oil during a couple of days, and each
$10\text{\,}\mathrm{s}$ during some hours for alkanes.

At any given measuring time, three images, one for each colour is taken. This
leads to interference patterns such as the one shown in Fig. 2b). From such an
image, an interference profile $I_{C}(x)$ is obtained at each time and for
each colour $C=R,\,G,\,B$ by averaging the intensity profiles along the
$y$-axis within a region spaning around $1\text{\,}\mathrm{mm}$, which
corresponds to the rectangle in Fig. 2b). To obtain the thickness $h(x)$ from
$I_{C}(x)$, one considers the theoretical linearised expression of
transmission interference through a thin film of thickness $h$,

$I_{C}^{th}(h)\propto 1-A\cos\left(\dfrac{4\pi nh}{\lambda_{C}}\right)$ (1)

where $n$ is the refractive index of the liquid and $A$ is a Fresnel
coefficient depending on the refractive indexes of the media considered, in
our experiments $A\simeq 0.012$. In the computations, the interference
profiles $I_{C}(h)$ are convolved with a Gaussian function mimicking the
broad-band emission of the LED source for each channel.

Figure 2: a) Sketch of the measurement setup. A RGB Panel illuminates from
below the experimental cell, consisting of a glass substrate over which a
liquid film is deposited via spin-coating. A monochrome camera on top takes
pictures at specified intervals. b) Typical interferometric profile of a
silicone oil film obtained after an initial transient period. Here
$\lambda_{R}=$635\text{\,}\mathrm{nm}$$ (red) and the height different between
intensity maxima is $\Delta h=$113\text{\,}\mathrm{nm}$$. The black rectangle
indicates the region of interest for measuring the film dynamics, the
interference intensity $I_{R}(x)$ is obtained by performing an average along
the vertical direction (and similarly for green, $I_{G}$ and blue, $I_{B}$).
c) Grey dots: interferometric data in intensity space $(I_{R},I_{G},I_{B})$
for the film at a given time. Black line: theoretical intensity curve, with
the thickness $h$ as parameter. The coloured points and lines are projections
into the orthogonal planes for better visualization. Fitting the data points
to the theoretical curve gives an absolute measurement of the thickness $h$ at
the corresponding $x$ coordinate. The squares corresponds to
$h=$80\text{\,}\mathrm{nm}$$, the stars to $h=$300\text{\,}\mathrm{nm}$$.

Working on the intensity space $\mathbf{I}=(I_{R},I_{G},I_{B})$, as shown in
Fig. 2c), it is possible to use $h$ as a fitting parameter of the theoretical
curve for the data points, resulting in an absolute-thickness curve $h(x)$ at
the time of measurement. This procedure allows to obtain thicknesses below the
first-order fringe, i.e. smaller than $100\text{\,}\mathrm{nm}$, thus in the
range where intermolecular forces are to be expected. In our current
implementation, overall thickness profiles ranging from the tens of nanometres
to the micron scale can be measured with a nanometric accuracy. The floor
sensitivity of the measurement is limited by experimental noise, which mainly
originates from the camera sensor; as a result, thicknesses below $10$ nm
cannot be determined. Due to the opacity of the heating line, the film
thickness cannot be measured for $|x|\lesssim$0.1\text{\,}\mathrm{mm}$$.

###  II.4 Numerical integration

Numerical integration of the thin film equation discussed below is done using
the pdepe numerical solver of Matlab
[*[See][forthedocumentationpageofpdepe.RunwithMatlabR2022a.]matlab_pdepe].
Since the equation is fourth order on the thickness gradient, the second
derivative $u(t,x)\equiv\partial_{x}^{2}h$ is considered as an independent
field and a simple closure expression
$u-\partial_{x}\left\\{\partial_{x}h\right\\}=0$ is simultaneously solved. The
initial condition is a film of homogeneous thickness, $h(0,x)=h_{0}$,
$u(0,x)=0$. For ease of convergence, the $x$ spacing is logarithmic, allowing
to better describe the evolution close to the centre as well as capturing the
long range evolution. The numerical integration typically lasted some hours.

##  III Theoretical background

As schemed in figure 3, we consider a homogeneous layer of incompressible
liquid of constant viscosity $\eta$ and density $\rho$ with initial thickness
$h_{0}$, that is submitted to a symmetric temperature excess $\theta(x)$. The
excess temperature remains small: $\theta_{max}=\theta(0)\ll T_{0}$, with
$T_{0}$ the room temperature and all the properties of the liquid are assumed
to be uniform. The resulting dynamics is assumed to be invariant along the
$y$-axis. In addition, since we use poorly volatile fluids we consider that
evaporation is negligible.

###  III.1 Hypotheses and dimensionless numbers

The order of magnitude of the thermocapillary stress is
$\gamma_{\theta}\theta_{max}/w\sim$0.1\text{\,}\mathrm{Pa}$$. The resulting
velocity in a liquid film of thickness $10\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ and viscosity $10^{-2}$ $\mathrm{Pa}\text{\,}\mathrm{s}$ is
$v\sim$100\text{\,}\mathrm{\SIUnitSymbolMicro m}\text{\,}{\mathrm{s}}^{-1}$$.
Since the liquids considered in the experiments have a thermal diffusivity
$\kappa\sim 10^{-7}${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$$ [28, 29], the
thermal Péclet number $\text{Pe}=h_{0}v/\kappa\sim 10^{-2}$ is small enough to
consider thermal convection negligible. Therefore, the transport of heat from
the substrate to the liquid is purely diffusive. In addition, the thermal
field is established in a time $h_{0}^{2}/\kappa\sim 10^{-2}$ $\mathrm{s}$,
much faster than the timescale of hours involved in the dynamics. Hence, we
assume that the thermal field is established instantaneously and that it is
constant throughout the film’s thickness. The Reynolds number
$\text{Re}=h_{0}v/\eta\sim 10^{-7}$ is also very small.

In the literature, a dimensionless number has been defined to measure the
effect of the thermocapillary stress compared with that of hydrostatic
pressure, $R=3\Delta\gamma/2\rho gh_{0}^{2}$, where
$\Delta\gamma=\gamma_{\theta}\theta_{max}$ [6]. We typically work with $R\sim
10$, a larger thermocapillary effect than in Ref. [30] so that the thickness
at the centre will always reach the region where intermolecular effects are
relevant [6].

###  III.2 Thin-film equation

In experiments, the thermal field decays over a typical length
$w\sim$2\text{\,}\mathrm{mm}$$ much larger than the initial thickness,
$h_{0}\ll w$, allowing to consider a small gradient approximation. As a
result, it is possible to obtain a lubrication approximation of the Navier-
Stokes equations leading to a partial differential equation for the thickness
profile $h(t,x)$ belonging to the family of thin film equations [31].

Figure 3: Sketch of the film dynamics with main variables: a thin liquid film
of density $\rho$ and viscosity $\eta$, with initial thickness $h_{0}$, is
deposited on top of a glass substrate which is heated locally, $w$ being the
decay length of the thermal field.

For the system defined by the conditions above and schematized in Fig. 3, the
thin film equation reads [6]

$\partial_{t}h+\dfrac{1}{3\eta}\partial_{x}\left\\{h^{3}\partial_{x}\left[\gamma_{0}\partial_{x}^{2}h-\rho
gh+\Pi(h)\right]\right\\}-\dfrac{\gamma_{\theta}}{2\eta}\partial_{x}\left\\{h^{2}\partial_{x}\theta(x)\right\\}=0,$
(2)

where the pressure contributions considered on the second term of the l.h.s.
are the capillary pressure, the hydrostatic pressure and $\Pi(h)$ is a
disjoining pressure term accounting for the intermolecular interactions. The
third term represents the thermocapillary force driving film thinning.
Capillarity and gravity oppose film thinning whereas disjoining pressure may
oppose or amplify thinning, according to its sign. The relative importance of
the pressure terms varies as the film thins down and this variation has been
investigated theoretically [30]. If the initial thickness of the film is
larger than the range of molecular interactions, then in a first stage, the
film dimples and both gravity and capillarity contribute to slow down its
evolution. The central region of the film further empties and a central area
of very small thickness is formed. In this region of small thickness and small
curvature, both capillary and hydrostatic pressures are negligible. As a
result, film thinning results from an interplay of thermocapillary stress and
disjoining pressure. If the intermolecular forces are attractive then
spontaneous rupture of the film occurs. In the case of repulsive
intermolecular forces, a fully equilibrated stationary state can be reached
and has been observed in experiments [6, 32]; it corresponds to a balance of
the thermocapillary stress and disjoining pressure. In such a state, it would
be possible, in principle, for the disjoining pressure to be measured if the
thickness of the ultra-thin liquid film at the central region is measured and
the imposed thermal field is known. However, most experimental works have not
aimed at measuring disjoining pressure and rather focused on the question of
the film rupture [30, 11]. In addition, an accurate determination of the
thermal field is required to measure the disjoining pressure, which is a
challenging task. In the geometry we investigate, the thermal gradients extend
over millimetric lengths and we observe a very slow dynamics of the central
region, which thins down during a time scale of days. In practice, this
implies that the full equilibrium cannot be experimentally reached since other
effects, such as sample contamination or liquid evaporation, appear at these
very long time scales. Nevertheless, the thinning dynamics itself provides
information on the intermolecular forces. Therefore, rather than trying to
reach a stationary state, we take profit of the thinning dynamics to probe the
intermolecular interactions at stake.

###  III.3 Pure thermocapillary regime

Before focusing on the effects of disjoining pressure, we first examine the
preceding stage, in which the central area of small thickness has formed but
its thickness is large enough for disjoining pressure to be neglected, i.e.
larger than roughly $100$ nm. In this stage, the capillary and gravitational
pressures can as well be neglected in the description of the evolution,
consistently with former works [6, 30]. Film thinning is then entirely driven
by thermocapillarity, which no effect opposes. Equation (2) then becomes

$\partial_{t}h-\dfrac{\gamma_{\theta}}{2\eta}\partial_{x}\left\\{h^{2}\partial_{x}\theta(x)\right\\}=0.$
(3)

Equation (3) can be solved by separation of variables, giving

$h(t,x)=\dfrac{h_{0}}{1+(t/t^{*})}\dfrac{|\partial^{2}_{x}\theta(0)|}{2\sqrt{|\partial_{x}\theta(x)|}}\int_{0}^{x}\dfrac{ds}{\sqrt{|\partial_{s}\theta(s)|}},$
(4)

where $t^{*}=2\eta/(\gamma_{\theta}|\partial^{2}_{x}\theta(0)|h_{0})$. Since
$t^{*}$ is at most a few seconds and the timescale of the experiments is
hours, the limit $t\gg t^{*}$ is quickly reached. Equation (4) then simplifies
and the product $t\times h(t,x)$ can be expressed as

$t\times
h(t,x)=\dfrac{\eta}{\gamma_{\theta}}\dfrac{1}{\sqrt{|\partial_{x}\theta(x)|}}\int_{0}^{x}\dfrac{ds}{\sqrt{|\partial_{s}\theta(s)|}}.$
(5)

The product $t\times h(t,x)$ is therefore independent of time and solely a
function of the spatial coordinate, $x$. Since in addition it only depends on
the liquid properties and not on the initial film thickness, for a given
liquid, $t\times h(t,x)$ curves are expected to collapse onto a master curve
that is a function of $x$ only, whatever the initial conditions of the
experiment. Inversely, at any given location $x$, the thickness has a simple
dependency on time: $h(t,x)\propto t^{-1}$. This temporal scaling has been
reported in experiments performed by heating a substrate with a laser beam
[11]. The variations with time of the thickness of the oil film spread on the
substrate were found to follow $h(t,0)\propto t^{-1}$ for a wide range of
laser power and beam waists under half a micrometre. It confirms the relevancy
of a regime in which both capillarity and gravity are negligible. In the
former work, the role of capillarity was found to be apparent only for very
small beam waists at and below the micrometre scale, where the curvature of
the film is sufficiently large. We demonstrate in the following that, in our
experiments, capillarity and gravity are always negligible in the central
area.

##  IV Experimental results

###  IV.1 Master evolution at intermediate thicknesses

In figure 4, we show the result of a typical measurement of the evolution at
the central area, extending roughly a few millimeters away from the central
line, similar as the one in figure 1b). The evolution is displayed, between 1
and 6 hours after turning the thermal gradient on. Since the evolution is
symmetric with respect to the origin of the $x$-axis, the thickness of only
one side of the film is reported. The experiment uses a film of silicone oil
with initial thickness $h_{0}=$36(2)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$
and the dissipated power is $P=$0.81(1)\text{\,}\mathrm{mW}$$.

The film thickness increases very slowly from the centre to the edge: it
reaches a few micrometers only over a distance of a few millimeters. The
spatial variations of $t\times h(t,x)$ are also displayed in the figures:
remarkably, all curves nicely collapse onto a master curve. According to the
analysis of section III, the collapse indicates that the effects of both
capillary and hydrostatic pressures are negligible and that the dynamics is
driven by thermocapillarity only, following the simple time dependency given
by equation (5).

Figure 4: a) Experimental profiles of an experiment with silicone oil with
$h_{0}=$36(2)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$P=$0.81(1)\text{\,}\mathrm{mW}$$. The time between profiles is
$10\text{\,}\mathrm{min}$. The inset shows the same profiles in linear scale.
b) The products $t\times h(t,x)$ for the experimental profiles of panel a) at
corresponding times collapse onto a master curve, which indicates that the
local dynamics is solely described by the thermocapillary effect. The dashed
line is the average of the experimental curves that is used to determine the
thermal gradient.

Following equation (5), the obtained master curve depends on the liquid
properties, which are known, and on the thermal gradient driving the dynamics.
The latter can hence be determined from the experimental master curve by
inverting the relation given by (5) to get

$\partial_{x}\theta(x)=-\dfrac{2\eta}{\gamma_{\theta}}\dfrac{1}{\left\langle
t\times h(x,t)\right\rangle^{2}_{exp}}\int_{0}^{x}\left\langle t\times
h(s,t)\right\rangle_{exp}ds.$ (6)

where $\left\langle t\times h(x,t)\right\rangle_{exp}$ is the experimentally
determined master curve from the collapse of profiles as in figure 4b), and is
only a function of $x$. Figure 5 shows the thermal gradient
$\partial_{x}\theta(x)$ obtained in this way. Since we do not have data up to
$x=0$, we extrapolate $\left\langle t\times h(x,t)\right\rangle_{exp}$ by a
constant for $x\leq$0.1\text{\,}\mathrm{mm}$$ in order to numerically compute
the integral in equation (6). This extrapolation satisfies the symmetry
requirement $\partial_{x}h(t,0)=0$ and we have checked that it results in a
negligible correction for larger values of $x$.

Figure 5: Derivative of the thermal field $\partial_{x}\theta(x)$ obtained
experimentally from the collapse of curves shown in panel b) of Fig. 4 using
equation (6). The dashed line represents a fit to the experimentally
determined gradient following the derivative of the thermal excess given by
equation (7) with the parameter values given in the text. Inset: integrated
thermal excess inferred from the experimentally determined gradient. The
dashed line represents equation (7).

A heuristic expression for $\partial_{x}\theta(x)$ can further be determined
from the temperature gradient found from the master curve. Once integrated, it
corresponds to a thermal field given by

$\theta(x)=\theta_{max}\left[\operatorname{sech}\left(\dfrac{x}{\delta}\right)\right]^{\delta/w}$
(7)

where $\delta$ and $w$ are two characteristic lengths. The length $\delta$
accounts for the finite length of the zone of maximum temperature. This length
results from both the finite size of the central line of the circuit
dissipating heat and the diffusion of heat through the thickness of the
substrate. Since the line width and substrate thickness are both of
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$, $\delta$ is expected to be close
to this value. The length $w$ is the decay length of the thermal field at
distances from the central line larger than $\delta$ and smaller than roughly
$1\text{\,}\mathrm{c}\mathrm{m}$. We emphasise that at a longer range, the
thermal excess is expected to go to zero asymptotically, with a characteristic
decay length one order of magnitude larger than $w$. However, since we are
interested in the central area of the film, for which submicron thicknesses
are reached, equation (7) provides a satisfactory description of the thermal
gradient in this area and further detail on the thermal field is not needed.
In the particular case of Fig. 6, a fit of equation (7) to the experimentally
determined gradient yields the following values:
$\theta_{max}=$0.37\text{\,}\mathrm{K}$$,
$\delta=$110\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$w=$2.3\text{\,}\mathrm{mm}$$. These values are consistent with independent
measurements of the thermal field at the bare substrate measured with an IR
camera. We emphasise that the modest value of $\theta_{max}$ confirms that the
liquids can be considered as having constant physical properties (i.e.
viscosity, optical index and density), although they are heated.

In summary, the master evolution followed by the film thickness within
intermediate times allows the determination of the thermal gradient. The
empirical expression that is obtained can be used in numerical resolutions to
test the consistency of the results, as described in the next section.

###  IV.2 Comparison of experimental data with numerical resolution

Since the thermal field is determined, equation (2) that includes capillary
and gravity terms can be solved, allowing the comparison of experimental data
with numerical results obtained without any assumption. We emphasise that, in
the following analysis, we do not consider the disjoining pressure term since
we focus here on large thickness values, for which intermolecular forces are
negligible. In Fig. 6, we can see that the numerical integration faithfully
represents the measured local dynamics, even for thicknesses well below the
micron scale. Therefore, the determination of the thermal gradient based on
the assumption of negligible effects of capillarity and gravity is fully
validated.

Figure 6: a) Same experimental thickness profiles as in figure 4. The solid
lines were obtained by numerical integration of equation (2), with the thermal
gradient of figure 5 and a zero disjoining pressure. b) Zoomed view of the
curves at small thicknesses showing the discrepancy between numerical and
experimental data at this small scale (for better visualization, only one out
of two curves is shown).

As expected, at thicknesses below $100$ nm, we observe that the numerical
profiles deviate from experimental data. In this thickness range, both
capillarity and gravity can still be neglected and the dynamics results from
thermocapillary and intermolecular effects. We have not tried to
quantitatively describe the disjoining pressure term in silicone oil. Instead,
we focus on a qualitative effect resulting from intermolecular forces that we
have evidenced.

###  IV.3 Transition to the ultra-thin film

We now examine the long term regime, when the thickness is on the order of
tens of nanometres, and when the effect of intermolecular forces becomes
important. In this regime, the reduced thermocapillary evolution (4) is not
valid any more, and the curves $t\times h(t,x)$ are expected to deviate from
the master curve.

Figure 7 shows the long-term regime of the evolution when the thickness enters
the range in which intermolecular forces are relevant, for the same experiment
as Fig. 6. The corresponding $t\times h(t,x)$ curves clearly lie above the
master curve, which demonstrates that the thinning induced by
thermocapillarity has been slowed down. This deceleration can be explained by
the presence of an effective repulsive intermolecular interaction that opposes
the thermocapillary effect. A repulsion between the interfaces of the silicone
oil film is expected since silicone oil totally wets the glass surface [33].
Therefore, the thinning slows down and the final stationary state would
correspond to an equilibrium of disjoining pressure and thermocapillary
stress.

Figure 7: a) and b): same plots as in Fig. 4 (dimmed curves), with the
addition of profiles at later times $t=13$, $23$, $33$,
$43\text{\,}\mathrm{h}$. The deviation of these last profiles from the master
curve in panel b) reveals the effect of intermolecular forces. c) Experimental
profiles $h(t,x)$, and d) curves $t\times h(t,x)$ for an experiment using
n-hexadecane ($n$-C16), with
$h_{0}=$13.0(2)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$P=$0.95(9)\text{\,}\mathrm{mW}$$. The time between profiles is
$1\text{\,}\mathrm{min}$. Observe that the time arrow on panel d) points in
the opposite direction as the corresponding arrow on panel b).

We have performed similar experiments with alkanes, and panels c) and d) of
figure 7 display the thickness profiles measured in a film of $n$-hexadecane.
Because of the smaller viscosity of $n$-hexadecane, the time scales of the
experiment are much shorter. The curves at earlier times are nevertheless
qualitatively similar to those of silicone oil and collapse as well onto a
master curve. In contrast, at later times, the experimental profiles obtained
with $n$-hexadecane reveal an acceleration of film thinning since they lie
below the master curve. The acceleration is associated with an unstable
behavior: the thickness profiles of $n$-hexadecane clearly collapse and form a
foot receding from the centre. We were however unable to evidence the presence
or absence of an ultra-thin film in the central region, since, if it exists,
its thickness lies below the sensitivity level of the measurement
($10\text{\,}\mathrm{nm}$). We can nevertheless state that, in the range of
investigated thicknesses, the disjoining pressure for $n$-hexadecane is
attractive since it accelerates thinning instead of opposing it, in contrast
to silicone oil.

This result is in contradiction with reports in the literature of hexadecane
fully wetting glass surfaces [34], hence suggesting that the interactions
between the interfaces of $n$-hexadecane films should be repulsive. The
attractive interactions we evidence rather evoke a partial or pseudo-partial
wetting of $n$-hexadecane on glass. In the case of partial wetting, the
spreading parameter is negative and the equilibrium situation is a drop
sitting on the substrate with a finite angle of the contact line; in the case
of pseudo-partial wetting, the spreading parameter is postive and the
equilibrium situation is also a drop with a finite angle but that sits on a
film of nanometric thickness. The latter situation results from non-monotonic
variations of the disjoining pressure with thickness [33]. Since glass is a
material of high surface energy, it is likely that the spreading parameter is
positive, favoring the hypothesis of pseudo-partial wetting.

We emphasise that we observe the same unexpected behavior with the ramified
alkane HMN as with the linear $n$-hexadecane, and hence that it does not
result from crystallisation effects. Actually, it is well known that surface
freezing occurs in alkanes at temperatures up to $3\text{\,}\mathrm{K}$ above
their solidification temperature [35], which affect the wetting behavior of
$n$-hexadecane films [36, 37]. However, HMN has a crystallisation temperature
which is well below room temperature, and remains liquid at room temperature
even when present on sub-nanometric layers [38]. As a result, surface freezing
cannot be invoked to explain the observed behaviors that are common to
$n$-hexadecane and HMN. We rather suggest that these alkanes only pseudo-
partially wet glass. In this case, the distinction from total wetting is
difficult in macroscopic experiments if the equilibrium angles are small. In
addition, times to attain an equilibrium state may reach the scale of the week
[39]. These difficulties are likely to explain why $n$-hexadecane has been
suggested to totally wet glass, and show that static experiments at the
macroscopic scale may not be sufficient to determine wetting behaviors
properly.

Additional insight is provided by numerical integration of the thin-film
equation with a given disjoining pressure. In the apolar liquids and the range
of film thicknesses we investigate, intermolecular forces are expected to be
non-retarded van der Waals interactions $\Pi(h)=-A_{H}/(6\pi h^{3})$, where
$A_{H}$ is the Hamaker constant [40]. A repulsive interaction then corresponds
to a negative Hamaker constant and an attraction to a positive Hamaker
constant. Using typical values for the thermal gradient, we have solved
equation (2) with disjoining pressure terms of respectively attractive and
repulsive nature. Because of numerical instabilities, full spatial variations
as in figure 7 cannot be obtained in the case of attractive intermolecular
forces. We nonetheless show in figure (8) the dynamics at the centre,
$h(t,0)$, for attractive ($A_{H}>0$) and repulsive ($A_{H}<0$) interactions,
as well as for zero disjoining pressure. Since the dynamics at the centre
cannot be measured in the experiment, the dynamics measured at a distance of
$0.2$ mm from the centre $h(t,x=$0.2\text{\,}\mathrm{mm}$)$ is displayed for
experiments with silicone oil and both $n$-hexadecane and HMN. In order to
compare the dynamics measured in different experimental conditions with
liquids that have different viscosities, we have represented the dimensionless
thickness a function of the ratio of the time and a thermocapillary timescale,
$t_{0}=2\eta w^{2}/(\Delta\gamma h_{0})$. As expected, the case $A_{H}=0$
results in a pure-thermocapillary evolution $h(t,0)\propto t^{-1}$. In
contrast, for non-zero interactions, a deviation from $t^{-1}$ variations
appears for thicknesses in the range of tens on nanometres. For $A_{H}<0$ the
thinning is slowed down, whereas for $A_{H}>0$ it is abruptly accelerated,
leading to film rupture.

Figure 8: a) Comparison of the thickness evolution near the centre,
$h_{c}(t)=h(t,x=$0.2\text{\,}\mathrm{mm}$)$, for the experiment with silicone
oil and with $n$-hexadecane of Fig. 7. The time has been scaled by a
thermocapillary timescale $t_{0}=2\eta w^{2}/(\Delta\gamma h_{0})$ and the
thickness by $h_{0}$. Silicone oil:
$h_{0}=$36(2)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\theta_{max}=$0.37\text{\,}\mathrm{K}$$, $w=$2.3\text{\,}\mathrm{mm}$$
($t_{0}\simeq$211\text{\,}\mathrm{s}$$). n-hexadecane:
$h_{0}=$13.0(2)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\theta_{max}=$0.44\text{\,}\mathrm{K}$$, $w=$2.2\text{\,}\mathrm{mm}$$
($t_{0}\simeq$42\text{\,}\mathrm{s}$$). HMN:
$h_{0}=$19.5(6)\text{\,}\mathrm{\SIUnitSymbolMicro m}$$,
$\theta_{max}=$0.54\text{\,}\mathrm{K}$$, $w=$2.2\text{\,}\mathrm{mm}$$
($t_{0}\simeq$30\text{\,}\mathrm{s}$$). b) Same as in a) but with profiles at
the centre $h(t,0)$ from the numerical integration of equation (2) with
physical properties of $n$-hexadecane, when the intermolecular forces generate
an effective repulsion ($A_{H}<0$) vs. attraction ($A_{H}>0$) of the
interfaces. In the case where there is no disjoining pressure (e.g.
$A_{H}=0$), the evolution follows equation (4), i.e. $t\times h(t,0)\propto
const.$ for $t\gg t_{0}$.

Changing the sign of the disjoning pressure in the equation therefore
satisfactorily accounts for the different behaviors observed. This result
confirms that attractive intermolecular forces, at least in the thickness
range investigated, are at stake between the interfaces of films of the
considered alkanes. Beyond this qualitative finding, the very similar shape of
the experimental and computed curves suggest that the expression for
disjoining pressure we use is close to the real one. However, we have found
that the dynamics in silicone oil could not be quantitatively described by the
simple law we have used: clearly in figure (8), the experimental thinning is
faster than in the computation. Consistently, the full thickness profiles, as
the ones reported in figure 6 could not be described, even when changing the
value of the Hamaker constant. This discrepancy may result from the fact that
the used silicone oil is not monodisperse but rather constituted by a mixture
of macromolecules of different lengths. Mixtures effects have actually been
predicted to change the exponent of the thickness-dependency of molecular
interactions [41]. These complex effects are out of the scope of the present
work and we have not tried to further determine the disjoining pressure at
stake. However, since silicone oil is widely used as a model wetting liquid,
it could be an interesting aim for future research.

##  V Conclusion

In conclusion, we have demonstrated that the dynamics of thin films induced by
a thermal gradient can serve as a tool for investigating intermolecular
forces. We predict and experimentally evidence a regime in which
thermocapillarity fully drives thinning. The thermal gradient can thus be
determined from thickness profile measurements and, in turn, be used to
compute the expected dynamics when intermolecular forces are at stake. We
demonstrate that thickness measurements allow a clear discrimination between
the cases of attractive and repulsive interactions between the interfaces of
the liquid films. The experiments reveal an unexpected behavior of both linear
and ramified alkanes when film thicknesses reach a few tens of nanometres, and
thus show that accurate probes of liquid films at very small scales are still
needed.

###### Acknowledgements.

This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No 754387.

## References

  * Hershey [1939] A. V. Hershey, Ridges in a liquid surface due to the temperature dependence of surface tension, Phys. Rev. 56, 204 (1939).
  * Anthony and Cline [1977] T. R. Anthony and H. E. Cline, Surface rippling induced by surface‐tension gradients during laser surface melting and alloying, J. App. Phys. 48, 3888 (1977).
  * Willis and Xu [2000] D. A. Willis and X. Xu, Transport phenomena and droplet formation during pulsed laser interaction with thin films, J. Heat Transfer 122, 763 (2000).
  * Singer [2017] J. P. Singer, Thermocapillary approaches to the deliberate patterning of polymers, J. Polym. Sci. Part B: Polym. Phys. 55, 1649 (2017).
  * Dai  _et al._ [2021] Q. Dai, W. Huang, X. Wang, and M. Khonsari, Directional interfacial motion of liquids: Fundamentals, evaluations, and manipulation strategies, Tribology International 154, 106749 (2021).
  * Tan  _et al._ [1990] M. J. Tan, S. G. Bankoff, and S. H. Davis, Steady thermocapillary flows of thin liquid layers. I. Theory, Phys. Fluids A 2, 313 (1990).
  * Pimputkar and Ostrach [1980] S. M. Pimputkar and S. Ostrach, Transient thermocapillary flow in thin liquid layers, The Physics of Fluids 23, 1281 (1980).
  * Yeo  _et al._ [2003] L. Y. Yeo, R. V. Craster, and O. K. Matar, Marangoni instability of a thin liquid film resting on a locally heated horizontal wall, Phys. Rev. E 67, 056315 (2003).
  * Hitt and Smith [1993] D. L. Hitt and M. K. Smith, Radiation‐driven thermocapillary flows in optically thick liquid films, Phys. Fluids A 5, 2624 (1993).
  * Chraïbi and Delville [2012] H. Chraïbi and J.-P. Delville, Thermocapillary flows and interface deformations produced by localized laser heating in confined environment, Physics of Fluids 24, 032102 (2012).
  * Wedershoven  _et al._ [2014] H. M. J. M. Wedershoven, C. W. J. Berendsen, J. C. H. Zeegers, and A. A. Darhuber, Infrared laser induced rupture of thin liquid films on stationary substrates, App. Phys. Lett. 104, 054101 (2014).
  * Pottier  _et al._ [2015] B. Pottier, C. Fretigny, and L. Talini, Boundary condition in liquid thin films revealed through the thermal fluctuations of their free surfaces, Phys. Rev. Lett. 114, 227801 (2015).
  * Klyuev  _et al._ [2021] D. S. Klyuev, V. M. Fliagin, S. V. Semenov, and N. A. Ivanova, Thermocapillary deformation induced by laser heating of thin liquid layers: Physical and numerical experiments, International Journal of Heat and Mass Transfer 172, 121020 (2021).
  * Moisy  _et al._ [2009] F. Moisy, M. Rabaud, and K. Salsac, A synthetic schlieren method for the measurement of the topography of a liquid interface, Exp. Fluids 46, 1021 (2009).
  * Roberts  _et al._ [2017] C. C. Roberts, A. Graham, M. Nemer, L. M. Phinney, R. M. Garcia, M. M. Soehnel, and E. K. Stirrup, Physical properties of low-molecular weight polydimethylsiloxane fluids, Office of Scientific and Technical Information, U.S. Department of Energy 10.2172/1343365 (2017).
  * Dubey and Sharma [2008] G. P. Dubey and M. Sharma, Temperature and composition dependence of the densities, viscosities, and speeds of sound of binary liquid mixtures of 1-butanol with hexadecane and squalane, J. Chem. Eng. Data 53, 1032 (2008).
  * Luning Prak  _et al._ [2014] D. J. Luning Prak, S. M. Alexandre, J. S. Cowart, and P. C. Trulove, Density, viscosity, speed of sound, bulk modulus, surface tension, and flash point of binary mixtures of n -dodecane with 2,2,4,6,6-pentamethylheptane or 2,2,4,4,6,8,8-heptamethylnonane, J. Chem. Eng. Data 59, 1334 (2014).
  * Camin  _et al._ [1954] D. L. Camin, A. F. Forziati, and F. D. Rossini, Physical properties of n-hexadecane, n-decylcyclopentane, n-decylcyclohexane, 1-hexadecene and n-decylbenzene, J. Phys. Chem. 58, 440 (1954).
  * Van Raalte [1960] T. L. Van Raalte, Measurement of the refractive index of silicone oils for the UKAEA, J. Opt. Soc. Am. 50, 85_2 (1960).
  * [20] abcr GmbH Germany, AB112154 product page.
  * Ricci  _et al._ [1986] E. Ricci, R. Sangiorgi, and A. Passerone, Density and surface tension of dioctylphthalate, silicone oil and their solutions, Surface and Coatings Technology 28, 215 (1986).
  * Korosi and Kovats [1981] G. Korosi and E. S. Kovats, Density and surface tension of 83 organic liquids, J. Chem. Eng. Data 26, 323 (1981).
  * Ocko  _et al._ [1997] B. M. Ocko, X. Z. Wu, E. B. Sirota, S. K. Sinha, O. Gang, and M. Deutsch, Surface freezing in chain molecules: Normal alkanes, Phys. Rev. E 55, 3164 (1997).
  * [24] Chemical Book, 2,2,4,4,6,8,8-heptamethylnonane technical information.
  * Schilling  _et al._ [2004] J. Schilling, K. Sengupta, S. Goennenwein, A. R. Bausch, and E. Sackmann, Absolute interfacial distance measurements by dual-wavelength reflection interference contrast microscopy, Phys. Rev. E 69, 021901 (2004).
  * de Ruiter  _et al._ [2015] J. de Ruiter, F. Mugele, and D. van den Ende, Air cushioning in droplet impact. i. dynamics of thin films studied by dual wavelength reflection interference microscopy, Phys. Fluids 27, 012104 (2015).
  * [27] https://www.mathworks.com/help/matlab/ref/pdepe.html.
  * Bluestar Silicones [2021] Bluestar Silicones, Rhodorsil Oils 47 Technical information (2021).
  * Vélez  _et al._ [2015] C. Vélez, M. Khayet, and J. M. Ortiz de Zárate, Temperature-dependent thermal properties of solid/liquid phase change even-numbered n-alkanes: n-Hexadecane, n-octadecane and n-eicosane, Appl. Energy 143, 383 (2015).
  * Burelbach  _et al._ [1990] J. P. Burelbach, S. G. Bankoff, and S. H. Davis, Steady thermocapillary flows of thin liquid layers. II. Experiment, Phys. Fluids A 2, 322 (1990).
  * Oron  _et al._ [1997] A. Oron, S. H. Davis, and S. G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod. Phys. 69, 931 (1997).
  * Clavaud  _et al._ [2021] C. Clavaud, M. Maza-Cuello, C. Frétigny, L. Talini, and T. Bickel, Modification of the Fluctuation Dynamics of Ultrathin Wetting Films, Phys. Rev. Lett. 126, 228004 (2021).
  * De Gennes and Brochard-Wyart [2015] P.-G. De Gennes and F. Brochard-Wyart,  _Gouttes, bulles, perles et ondes_ (Belin Éducation, 2015).
  * Shiri  _et al._ [2021] S. Shiri, S. Sinha, D. A. Baumgartner, and N. J. Cira, Thermal marangoni flow impacts the shape of single component volatile droplets on thin, completely wetting substrates, Phys. Rev. Lett. 127, 024502 (2021).
  * Wu  _et al._ [1993] X. Z. Wu, E. B. Sirota, S. K. Sinha, B. M. Ocko, and M. Deutsch, Surface crystallization of liquid normal-alkanes, Phys. Rev. Lett. 70, 958 (1993).
  * Lazar  _et al._ [2005] P. Lazar, H. Schollmeyer, and H. Riegler, Spreading and two-dimensional mobility of long-chain alkanes at solid/gas interfaces, Phys. Rev. Lett. 94, 116101 (2005).
  * Lu  _et al._ [2010] L. Lu, K. J. Zander, and Y. Cai, Stability of the parallel layer during alkane spreading and the domain structures of the standing-up layer, Langmuir 26, 5624 (2010).
  * Gosvami  _et al._ [2008] N. N. Gosvami, S. K. Sinha, and S. J. O’Shea, Squeeze-out of branched alkanes on graphite, Phys. Rev. Lett. 100, 076101 (2008).
  * Pérez  _et al._ [2001] E. Pérez, E. Schäffer, and U. Steiner, Spreading dynamics of polydimethylsiloxane drops: Crossover from laplace to van der waals spreading, J. Colloid Interface Sci. 234, 178 (2001).
  * Churaev [2003] N. V. Churaev, Surface forces in wetting films, Colloid Journal 65, 263 (2003).
  * Derjaguin and Churaev [1978] B. Derjaguin and N. Churaev, On the question of determining the concept of disjoining pressure and its role in the equilibrium and flow of thin films, J. Colloid Interface Sci. 66, 389 (1978).


# Multi-scale Imaging of Nuclear Deformation at the Electron Ion Collider

Heikki Mäntysaari  Department of Physics, University of Jyväskylä, P.O. Box
35, 40014 University of Jyväskylä, Finland  Helsinki Institute of Physics,
P.O. Box 64, 00014 University of Helsinki, Finland  Björn Schenke  Physics
Department, Brookhaven National Laboratory, Upton, NY 11973, USA  Chun Shen
Department of Physics and Astronomy, Wayne State University, Detroit, Michigan
48201, USA  RIKEN BNL Research Center, Brookhaven National Laboratory, Upton,
NY 11973, USA  Wenbin Zhao  Department of Physics and Astronomy, Wayne State
University, Detroit, Michigan 48201, USA

###### Abstract

We show within the Color Glass Condensate framework that exclusive vector
meson production at high energy is sensitive to the geometric deformation of
the target nucleus at multiple length scales. Studying $e+$U collisions and
varying the deformation of the uranium target, we demonstrate that larger
deformations result in enhanced incoherent vector meson production cross
sections. Further, different multipole deformation parameters affect different
regions of transverse momentum transfer. Employing JIMWLK evolution to study
the Bjorken-$x$ dependence of our results, we find that the ratio of
incoherent to coherent cross sections decreases with decreasing $x$, largely
independently of the quadrupole deformation of the target. Comparing results
for the same process using ${\rm{{}^{20}Ne}}$ targets with ${\rm{{}^{16}O}}$
targets, we find that differences in deformation are clearly visible in the
incoherent cross section. These findings show that certain observables at the
Electron-Ion Collider are very sensitive to nuclear structure. Consequently,
deformations need to be taken into account when interpreting experimental
results. More importantly, this also means that $|t|$-differential diffractive
vector meson production could become a powerful tool, enabling the most direct
measurements of nuclear structure at different length scales, ranging from
nuclear deformation at low $|t|$ to nucleon- and subnucleon-size scales at
higher $|t|$.

1\. Introduction. Understanding the geometric structures of the proton and
nuclei, including their event-by-event fluctuations, is of fundamental
interest. Deep inelastic scattering (DIS) of leptons on hadrons is one of the
most important tools to probe the partonic structure of protons and nuclei. In
the 2030’s, the Electron-Ion Collider in the US will provide access to
nuclear-DIS for the first time in collider kinematics Abdul Khalek  _et al._
(2022); Aschenauer  _et al._ (2019). In addition to the EIC, there are also
other longer-term proposals for future nuclear-DIS facilities at CERN Agostini
_et al._ (2021) and in China Anderle  _et al._ (2021).

Exclusive vector meson (e.g. $\mathrm{J}/\psi$) production off nucleons and
nuclei is a particularly clean and powerful process to probe the nuclear high-
energy structure at small longitudinal momentum fraction $x$ for several
reasons. First, in order to produce only the vector meson and nothing else,
net color charge cannot be transferred to the target, and the final state is
unambiguously identified by a large rapidity gap. This also requires that at
least two gluons are exchanged, which renders the cross section approximately
proportional to the squared gluon distribution Ryskin (1993) at leading order
(see Ref. Eskola  _et al._ (2022) for a recent analysis at next-to-leading
order accuracy). Additionally, only in such exclusive scattering it is
possible to determine the total momentum transfer to the target hadron or
nucleus, which is the Fourier conjugate to the impact parameter and as such
provides access to the target geometry.

In the hot QCD community, understanding the spatial shape of colliding objects
in heavy-ion collisions is a necessary input when the deconfined Quark-Gluon
Plasma (QGP) state is probed. The initial nucleon and nuclear geometry
determine the spatial distribution of nuclear matter in the collision, which
in turn determines the initial pressure anisotropies that are transformed into
observable momentum space correlations when the space-time evolution of the
QGP is simulated. The hydrodynamic modeling of heavy ion collisions with an
accurate initial state description in collisions of several different ion
species has revealed how the detailed structural properties of the colliding
nuclei are visible in final state particle correlations Filip  _et al._
(2009); Masui  _et al._ (2009); Hirano  _et al._ (2013); Shen  _et al._
(2016); Schenke  _et al._ (2014, 2020); Xu  _et al._ (2021); Giacalone  _et
al._ (2021); Zhang and Jia (2022); Nijs and van der Schee (2021); Cheng  _et
al._ (2023); Ryssens  _et al._ (2023).

The study of collisions involving heavy ions at the Relativistic Heavy Ion
Collider (RHIC), the Large Hadron Collider (LHC), and the future Electron-Ion
Collider (EIC) offers a synergistic approach to understanding the properties
of high energy nuclei and hot and dense nuclear matter. The EIC will probe in
detail the small-$x$ structure of the nuclei, obtaining fundamentally
interesting information and constraining the initial state of heavy ion
collisions Mäntysaari and Schenke (2016a, b); Mäntysaari (2020). In turn, as
mentioned above, heavy ion collisions can themselves constrain the high-energy
structure of nuclei, including deformation and nucleon clustering. The
complementary information from these facilities is essential for gaining a
complete picture of the properties of the QGP, and the influence of nuclear
geometry in the initial state.

We employ the Color Glass Condensate (CGC) framework Gelis  _et al._ (2010)
supplemented with a model that describes the nuclear geometry including
nucleon substructure in terms of gluonic hot spots Mäntysaari (2020).
Recently, in Ref. Mäntysaari  _et al._ (2022a), we performed a statistically
rigorous Bayesian analysis to extract the posterior likelihood distribution
for the model parameters describing the event-by-event fluctuating proton
geometry from HERA $\mathrm{J}/\psi$ production data Alexa  _et al._ (2013).
In this work, we extend our studies to $e+A$ collisions and for the first time
explore the effects of nuclear deformation on the diffractive cross sections
for large nuclei. Specifically, we study the dependence of the $|t|$
differential diffractive $\mathrm{J}/\psi$ cross sections on the deformation
parameters of uranium. To examine how energy (or Bjorken-$x$) evolution
modifies the effects of deformation, we perform numerical simulations solving
the JIMWLK equations (see e.g. Mueller (2001)) to evolve the nuclear
configurations to smaller Bjorken-$x$.

We further study vector meson production in electron scattering off smaller
nuclei. Of particular interest are ${\rm{{}^{20}Ne}}$ and ${\rm{{}^{16}O}}$,
which have a similar mass number, but are expected to differ in shape.
Calculations in different models Ebran  _et al._ (2012); Zhou  _et al._
(2016); Marević  _et al._ (2018); Frosini  _et al._ (2022a); Bally  _et al._
(2023) obtain a characteristic bowling pin shape for ${\rm{{}^{20}Ne}}$,
essentially forming a ${\rm{{}^{16}O}}$ like structure with an additional
$\alpha$ cluster in its periphery. We provide model predictions for future EIC
measurements, in particular for the ratio of the $|t|$ differential
diffractive $\mathrm{J}/\psi$ production cross sections in the two systems.

2\. Vector meson production at high energy. The total diffractive cross
section in DIS gives insight into the total small-$x$ gluon densities of the
target nuclei. More differential observables such as the exclusive production
of a vector meson, $\gamma^{*}+p/A\to V+p/A$, as a function of (squared)
momentum transfer $-t$ can provide more detailed information on the target
structure. The coherent cross section, corresponding to the process where the
target remains in the same quantum state, can be obtained by averaging over
the target color charge configurations $\Omega$ at the amplitude level Good
and Walker (1960):

$\frac{\differential\sigma^{\gamma^{*}+A\to
V+A}}{\differential|t|}=\frac{1}{16\pi}\left|\left\langle\mathcal{A}\right\rangle_{\Omega}\right|^{2}.$
(1)

The incoherent vector meson production cross section, for which the final
state of the target is different from its initial state, is obtained by
subtracting the coherent contribution from the total diffractive vector meson
production cross section Miettinen and Pumplin (1978); Caldwell and Kowalski
(2010); Mäntysaari (2020). The incoherent cross section thus has the form of a
variance

$\frac{\differential\sigma^{\gamma^{*}+A\to
V+A^{*}}}{\differential|t|}=\frac{1}{16\pi}\left[\left\langle\left|\mathcal{A}\right|^{2}\right\rangle_{\Omega}\right.-\left.\left|\left\langle\mathcal{A}\right\rangle_{\Omega}\right|^{2}\right]\,.$
(2)

Here $\mathcal{A}$ is the scattering amplitude for diffractive vector meson
production, which at high energy describes the splitting of the virtual photon
into a quark anti-quark pair, the pair’s subsequent interaction with the
target, followed by the formation of the vector meson. It can be written as
Kowalski  _et al._ (2006); Hatta  _et al._ (2017) (see also Refs. Mäntysaari
and Penttala (2022, 2021) for recent developments towards NLO accuracy)

$\mathcal{A}=2i\int\differential[2]{{\mathbf{r}_{\perp}}}\differential[2]{{\mathbf{b}_{\perp}}}\frac{\differential{z}}{4\pi}e^{-i\left[{\mathbf{b}_{\perp}}-\left(\frac{1}{2}-z\right){\mathbf{r}_{\perp}}\right]\cdot{\boldsymbol{\Delta}_{\perp}}}\\\
\times[\Psi_{V}^{*}\Psi_{\gamma}](Q^{2},{\mathbf{r}_{\perp}},z)N_{\Omega}({\mathbf{r}_{\perp}},{\mathbf{b}_{\perp}},x_{\mathbb{P}}).$
(3)

Here ${\mathbf{r}_{\perp}}$ is the transverse size of the $q\bar{q}$ dipole,
${\mathbf{b}_{\perp}}$ is the impact parameter measured relative to the target
center, and $Q^{2}$ is the photon virtuality. The fraction of the large photon
plus-momentum carried by the quark is given by $z$, $x_{\mathbb{P}}$ is the
fraction of the target longitudinal momentum transferred to the meson in the
frame where the target has a large momentum, and
${\boldsymbol{\Delta}_{\perp}}$ is the transverse momentum transfer, with
$-t\approx{\boldsymbol{\Delta}}_{\perp}^{2}$. The $\gamma^{*}\to q\bar{q}$
splitting is described by the virtual photon light front wave function
$\Psi_{\gamma}$ Kovchegov and Levin (2022). The vector meson wave function
$\Psi_{V}$ is non-perturbative and needs to be modeled, introducing some
uncertainty. Here, we use the Boosted Gaussian parametrization from Kowalski
_et al._ (2006), where the model parameters are constrained by the decay width
data.

Dependence on the small-$x$ structure of the target is included in the dipole
amplitude
$N_{\Omega}({\mathbf{r}_{\perp}},{\mathbf{b}_{\perp}},x_{\mathbb{P}})$, which,
for a given target color charge configuration $\Omega$, is
$N_{\Omega}({\mathbf{r}_{\perp}},{\mathbf{b}_{\perp}},x_{\mathbb{P}})=1-\frac{1}{{N_{\mathrm{c}}}}\tr\left[V\left({\mathbf{b}_{\perp}}+\frac{{\mathbf{r}_{\perp}}}{2}\right)V^{\dagger}\left({\mathbf{b}_{\perp}}-\frac{{\mathbf{r}_{\perp}}}{2}\right)\right].$
The $V({\mathbf{x}_{\perp}})$ represents a Wilson line, depending on $\Omega$
and $x_{\mathbb{P}}$, and describing the color rotation of a quark state when
propagating through the target field at transverse coordinate
${\mathbf{x}_{\perp}}$. The Wilson lines are obtained in the same way as in
the IP-Glasma initial state description Schenke  _et al._ (2012a) used e.g. in
Refs. Mäntysaari  _et al._ (2022b, 2021); Mäntysaari and Schenke (2020);
Mäntysaari  _et al._ (2019); Mäntysaari and Schenke (2018, 2016a, 2016b). They
are computed by first relating the average square color charge density to the
local saturation scale extracted from the IPSat dipole-proton amplitude
Schenke  _et al._ (2012b). Then, by solving the Yang-Mills equations for the
gluon fields, one obtains

$V({\mathbf{x}_{\perp}})=\mathrm{P}_{-}\left\\{\exp\left({-ig\int_{-\infty}^{\infty}\differential{z^{-}}\frac{\rho^{a}(x^{-},{\mathbf{x}_{\perp}})t^{a}}{\boldsymbol{\nabla}^{2}-m^{2}}}\right)\right\\}\,,$
(4)

where $\mathrm{P}_{-}$ represents path ordering in the $x^{-}$ direction and
$\rho^{a}$ is the color charge density. Here, we introduced the infrared
regulator $m$, which is needed to avoid the emergence of unphysical Coulomb
tails.

We note that in the incoherent cross section (2), the square of the impact
parameter dependent scattering amplitude is equivalent to the Fourier
transform of the two-point function of the nuclear thickness function, which
clarifies the sensitivity of this quantity to the target structure at
different length scales, depending on $|t|$ (see also Caldwell and Kowalski
(2010); Blaizot and Traini (2022)).

In this work we use sub-nucleonic fluctuations of the nucleon, introducing an
event-by-event fluctuating density by following Refs. Mäntysaari and Schenke
(2016a, b) and writing the density profile of nucleons
$T_{p}({\mathbf{b}_{\perp}})$ as

$T_{p}({\mathbf{b}_{\perp}})=\frac{1}{N_{q}}\sum_{i=1}^{N_{q}}p_{i}T_{q}({\mathbf{b}_{\perp}}-{\mathbf{b}_{\perp,i}}),$
(5)

where the single hot spot density distribution
$T_{q}({\mathbf{b}_{\perp}})=\frac{1}{2\pi
B_{q}}e^{-{\mathbf{b}}_{\perp}^{2}/(2B_{q})}\,$ and the coefficient $p_{i}$
allows for different normalizations for individual hot spots. It follows the
log-normal distribution with the width $\sigma$ controlling the magnitude of
the density fluctuations. Our prescription corresponds to having $N_{q}$ hot
spots with hot spot width $B_{q}$. The hot spot positions
${\mathbf{b}_{\perp,i}}$ are sampled from a two-dimensional Gaussian
distribution whose width is denoted by $B_{qc}$, and the center-of-mass is
shifted to the origin at the end. In this work, we use the Maximum a
Posteriori (MAP) parameter set from Bayesian analysis where the geometry
parameters at $x_{\mathbb{P}}\approx 0.0017$ are constrained by the exclusive
$\mathrm{J}/\psi$ production data from HERA Mäntysaari  _et al._ (2022a).

To model the geometric shape of large nuclei, we first sample nucleon
positions from a Woods-Saxon distribution

$\rho(r,\theta)=\frac{\rho_{0}}{1+\exp[(r-R^{\prime}(\theta))/a]}\,,$ (6)

with
$R^{\prime}(\theta)=R[1+\beta_{2}Y_{2}^{0}(\theta)+\beta_{3}Y_{3}^{0}(\theta)+\beta_{4}Y_{4}^{0}(\theta)]$,
and $\rho_{0}$ is the nuclear density at the center of the nucleus. Here $R$
is the radius parameter and $a$ the skin diffuseness, and $\theta$ is the
polar angle. A random rotation is applied after the sampling process. The
spherical harmonic functions $Y_{l}^{m}(\theta)$ and the parameters
$\beta_{i}$ account for the possible deformation from a spherical shape. The
default Woods-Saxon parameters for uranium are $\beta_{2}=0.28$,
$\beta_{3}=0$, $\beta_{4}=0.093$, $a=0.55$ fm, and $R=6.81$ fm Filip  _et al._
(2009); Masui  _et al._ (2009); Hirano  _et al._ (2013); Shen  _et al._
(2016); Schenke  _et al._ (2014, 2020). Following Moreland  _et al._ (2015);
Schenke  _et al._ (2020), we further impose a minimal distance of $d_{\rm
min}=0.9\,{\rm fm}$ between nucleons when sampling in three dimensions.111When
a nucleon is added and violates the minimum distance criterion with one or
more already sampled nucleons, we resample its azimuthal angle $\phi$ to keep
the distributions of radial distances and polar angles unchanged Moreland  _et
al._ (2015).

We also study smaller nuclei below. For the case of the nucleon density
distribution of ${{\rm{}^{20}Ne}}$, we use results from the ab initio
Projected Generator Coordinate Method (PGCM) Frosini  _et al._ (2022b, a, c);
Bally  _et al._ (2023). We also compare to the case of a spherical
${{\rm{}^{20}Ne}}$ nucleus described by a Woods-Saxon distribution with
parameters obtained in low energy electron-nucleus scattering De Vries  _et
al._ (1987). In this case the parameters are the radius $R=2.8$ fm, and skin
depth $a=0.57$ fm. For ${{\rm{}^{16}O}}$ we employ the nucleon density
distribution used in Loizides  _et al._ (2015), which is obtained from a
variational Monte-Carlo method (VMC) using the Argonne v18 (AV18) two-nucleon
potential+UIX interactions Carlson and Schiavilla (1998).

Figure 1: Coherent and incoherent $\mathrm{J}/\psi$ photoproduction cross
sections at $x_{p}=1.7\times 10^{-3}$ in e+U collisions for different
$\beta_{2}$ (a), $\beta_{3}$ (b) and $\beta_{4}$ (c) values. The bands show
the statistical uncertainty of the calculation.  Figure 2: Incoherent-to-
coherent $\mathrm{J}/\psi$ photoproduction cross section ratio as a function
of $x_{\mathbb{P}}$ in e+U collisions at different initial $\beta_{2}$ values.
Figure 3: Coherent and incoherent $\mathrm{J}/\psi$ photoproduction cross
sections in e+${}^{20}{\rm Ne}$ and e+ ${}^{16}{\rm O}$ collisions (top
panel), and the incoherent cross section ratios between ${}^{20}{\rm Ne}$ and
${}^{16}{\rm O}$ (bottom panel).

3\. Sensitivity of exclusive scattering to nuclear deformations. Eq. (1) shows
that the coherent cross section is sensitive to the average scattering
amplitude and as such probes the average structure of the target. The
incoherent cross section, Eq. (2), measures the scattering amplitude
fluctuations between the different possible color charge configurations.
Measuring the total momentum transfer ${\boldsymbol{\Delta}_{\perp}}$ allows
to constrain the geometry fluctuations in the target at length scales $\sim
1/\sqrt{-t}$. To determine how sensitive the future EIC measurements are to
the nuclear deformations in the currently unexplored small-$x$ region, we vary
one of the Woods-Saxon deformation parameters $\beta_{2}$, $\beta_{3}$ or
$\beta_{4}$ of uranium, while keeping the others set to their default values.

Figure 1 shows the effect of different levels of quadrupole ($\beta_{2}$),
octupole ($\beta_{3}$) and hexadecapole ($\beta_{4}$) deformation on the
coherent and incoherent
$\gamma^{*}+\mathrm{U}\to\mathrm{J}/\psi+\mathrm{U}^{(*)}$ cross sections for
a photon virtuality of $Q^{2}=0\,\mathrm{GeV}^{2}$. Increasing the degree of
deformation leads to increasing incoherent cross sections. For example, with a
realistic $\beta_{2}=0.28$ the incoherent cross section at $|t|\sim
0.01~{}{\rm GeV^{2}}$ is increased by a factor $\sim 5$ compared to the case
with an assumed spherical uranium target with $\beta_{2}=0.0$. This large
enhancement is a result of the random orientations of the deformed nucleus in
the laboratory frame leading to larger density fluctuations of the
configurations projected onto the transverse plane, which causes larger event-
by-event fluctuations in the scattering amplitude.

Importantly, the different types of deformations manifest in different $|t|$
regions and can as such be probed separately at the EIC. The $\beta_{2}$
modifies the incoherent cross section in the smallest $|t|\lesssim
0.015~{}\mathrm{GeV}^{2}$ region. On the other hand, $\beta_{3}$ is important
in the range $0.005~{}\mathrm{GeV}^{2}\lesssim|t|\lesssim
0.025~{}\mathrm{GeV}^{2}$, and $\beta_{4}$ for
$0.015~{}\mathrm{GeV}^{2}\lesssim|t|\lesssim 0.035~{}\mathrm{GeV}^{2}$. These
effects take place in different $|t|$ ranges because, for example, the
quadrupole deformation $\beta_{2}$ controls the geometric deformation of the
target at the longest length scale among these three deformation parameters,
which translates into the smallest $|t|$ region in momentum space.

The slope of the coherent cross section at low $|t|$ and the position of the
first diffractive minimum are not significantly modified by the deformations,
which means that the average size of the uranium nucleus is not modified.
However, the deformations do affect the coherent spectrum, especially at
higher $|t|$. This is due to the fact that (the projection of) the average
transverse density profile is different from that of the spherical nucleus
when we have non-zero deformation parameters $\beta_{i}$. Consequently, in
addition to the incoherent cross section, the coherent cross section in the
relatively large $|t|$ region (after the first diffractive minimum) can be
used to access deformations at the EIC.

In order to justify the use of deformation parameters extracted in low energy
experiments to the situation of high energy scattering, and in particular to
see whether the nuclear deformations are washed out at small-$x$, we apply the
perturbative JIMWLK evolution equation Mueller (2001) as in Refs. Mäntysaari
and Schenke (2018); Mäntysaari  _et al._ (2022b) to describe the Bjorken-$x$
dependence of the uranium structure. The uranium configurations at the initial
$x_{\mathbb{P}}=1.7\cdot 10^{-3}$ are generated using the same three different
$\beta_{2}$ values as above, and default values are used for $\beta_{3}$ and
$\beta_{4}$. Figure 2 shows the incoherent-to-coherent cross section ratio as
a function of $x_{\mathbb{P}}$, with the total (in)coherent cross sections
integrated within $0.0<|t|<0.5~{}{\rm{GeV^{2}}}$. This ratio effectively
suppresses the uncertainties originating from the modeling of the
$\mathrm{J}/\psi$ wave function Mäntysaari and Schenke (2017). The decreasing
ratio towards small $x_{\mathbb{P}}$ implies that the incoherent cross section
grows more slowly than the coherent one with increasing energy, because the
event-by-event fluctuations are reduced by the evolution Mäntysaari and
Schenke (2018) (see also Cepila  _et al._ (2018); Mäntysaari  _et al._
(2022b)). The difference in the cross section ratio between the different
initial quadrupole deformations remains similar throughout the evolution. We
conclude that the fluctuations in the nuclear geometry originating from the
deformed structure are not washed out by the JIMWLK evolution, and as such we
expect the deformations previously inferred from low-energy experiments to
also be visible in high-energy electron-ion collisions at the EIC.

Let us next demonstrate the possibility to probe deformations in the high-
energy structure of light nuclei at the EIC, focusing on ${\rm{{}^{20}Ne}}$
and ${\rm{{}^{16}O}}$ (see also Mäntysaari and Schenke (2020) where deuteron
and helium were considered). Figure 3 shows both the coherent and incoherent
$\mathrm{J}/\psi$ photoproduction cross sections at $x_{\mathbb{P}}=1.7\cdot
10^{-3}$ off a ${\rm{{}^{20}Ne}}$ nucleus computed from the ab initio PGCM
method Frosini  _et al._ (2022b, a, c), which resembles the shape of a bowling
pin. They are compared to the case where we neglect all deformations and use a
Woods-Saxon distribution as well as to the case of a ${\rm{{}^{16}O}}$ target,
described as discussed above.

Consistent with the uranium case shown in Fig. 1, the bowling pin like shape
deformations enhance the incoherent cross section when the PGCM
${\rm{{}^{20}Ne}}$ is compared to the spherical one at small $|t|\lesssim
0.05~{}\mathrm{GeV}^{2}$. To make the enhancement more visible, the bottom
panel of Fig. 3 shows the ratio of the incoherent cross sections for both
${\rm{{}^{20}Ne}}$ cases to the ${\rm{{}^{16}O}}$ case. We find that in case
of the PGCM ${\rm{{}^{20}Ne}}$ target, there is an enhancement of up to a
factor 2.7, compared to 1.5 when comparing spherical neon to oxygen. The
reason for the enhancement of spherical neon compared to oxygen may simply be
the larger size of ${\rm{{}^{20}Ne}}$, which is also visible in the locations
of the first diffractive minima in the coherent spectra. The additional
enhancement for the case of the PGCM ${\rm{{}^{20}Ne}}$ nucleus is a result of
its bowling pin shape, which essentially results from adding a fifth alpha
cluster to the four clusters in ${\rm{{}^{16}O}}$.

The largest enhancement is in the regime where the coherent cross section
dominates, however, even for $|t|\approx 0.05$ GeV2 where the incoherent cross
section begins to dominate the incoherent cross section of ${\rm{{}^{20}Ne}}$
is almost a factor of 2 larger than that of ${\rm{{}^{16}O}}$. At high $|t|$,
the incoherent cross sections only differ by a factor $20/16$ originating from
the different nuclear mass numbers that affect the overall normalization. This
means that the short distance scale fluctuations are identical in oxygen and
neon, as it was assumed when constructing the model.

If such measurements can be performed with enough precision at the EIC, they
will provide the most direct access to the structure of nuclei over all
relevant length scales, from nucleus to subnucleon size scales.

4\. Summary.  We have demonstrated that the coherent and incoherent exclusive
vector meson production measurements in $e+A$ collisions are affected by the
deformed structure of light and heavy ions. In particular, we have shown that
quadrupole, octupole, and hexadecapole deformations can significantly modify
the incoherent $\mathrm{J}/\psi$ production cross section. Different
deformations, on different length scales, affect different regions of momentum
transfer, and as such can be explored separately at the EIC. We also used
numerical solutions to the JIMWLK equation to describe the evolution of the
uranium fluctuating structure with decreasing $x_{\mathbb{P}}$, and found that
the evolution suppresses the incoherent cross section at high energies, but
does not significantly reduce the effects of the initial deformations.

The comparison between the ${\rm{{}^{20}Ne}}$ nucleus obtained in modern PGCM
calculations and the spherical ${\rm{{}^{20}Ne}}$ or ${\rm{{}^{16}O}}$ shows
that for light ions a non-trivial shape also results in enhanced incoherent
cross sections at low $|t|$. The strongest effect was observed in a $|t|$
region where the coherent cross section dominates, but even in the range where
the incoherent cross section is dominant, a factor of 2 enhancement was
observed between PGCM ${\rm{{}^{20}Ne}}$ and ${\rm{{}^{16}O}}$. We conclude
that correctly predicting incoherent cross sections for the EIC, even within
factors of 3 or more, requires the consideration and possibly even precision
calculation of the nuclear deformation of the studied target. Even more
importantly, the predicted large sensitivity of the incoherent cross section
to deformations at low $|t|$, along with the previously observed sensitivity
to the structure at nucleon and subnucleon size scales at higher $|t|$,
implies that this observable carries an unmatched amount of information on
nuclear structure over the entire range of relevant size scales.

This information will be complementary to that obtained from low energy
nuclear structure experiments. It will further have direct applications in
heavy ion collisions at RHIC, LHC, and other future facilities.

Acknowledgments. We acknowledge Benjamin Bally, Thomas Duguet, Jean-Paul
Ebran, Mikael Frosini, Giuliano Giacalone, Govert Nijs, Tomás Rodriguez,
Vittorio Somà, Wilke van der Schee for kindly sharing unpublished results on
the ground-state density of ${\rm{{}^{20}Ne}}$. We also thank Giuliano
Giacalone for useful discussions. B.P.S. and C.S. are supported by the U.S.
Department of Energy, Office of Science, Office of Nuclear Physics, under DOE
Contract No. DE-SC0012704 and Award No. DE-SC0021969, respectively. C.S.
acknowledges a DOE Office of Science Early Career Award. This material is
based upon work supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, within the framework of the Saturated Glue (SURGE)
Topical Theory Collaboration. H.M. is supported by the Academy of Finland, the
Centre of Excellence in Quark Matter, and projects 338263 and 346567, and
under the European Union’s Horizon 2020 research and innovation programme by
the European Research Council (ERC, grant agreement No. ERC-2018-ADG-835105
YoctoLHC) and by the STRONG-2020 project (grant agreement No 824093). W.B.Z.
is supported by the National Science Foundation (NSF) under grant number
ACI-2004571 within the framework of the XSCAPE project of the JETSCAPE
collaboration. The content of this article does not reflect the official
opinion of the European Union and responsibility for the information and views
expressed therein lies entirely with the authors. This research was done using
resources provided by the Open Science Grid (OSG) Pordes  _et al._ (2007);
Sfiligoi  _et al._ (2009), which is supported by the National Science
Foundation award #2030508.

## References

  * Abdul Khalek  _et al._ (2022) R. Abdul Khalek _et al._ , “Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report,” Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Aschenauer  _et al._ (2019) E. C. Aschenauer, S. Fazio, J. H. Lee, H. Mäntysaari, B. S. Page, B. Schenke, T. Ullrich, R. Venugopalan, and P. Zurita, “The electron–ion collider: assessing the energy dependence of key measurements,” Rept. Prog. Phys. 82, 024301 (2019), arXiv:1708.01527 [nucl-ex] .
  * Agostini  _et al._ (2021) P. Agostini _et al._ (LHeC, FCC-he Study Group), “The Large Hadron–Electron Collider at the HL-LHC,” J. Phys. G 48, 110501 (2021), arXiv:2007.14491 [hep-ex] .
  * Anderle  _et al._ (2021) Daniele P. Anderle _et al._ , “Electron-ion collider in China,” Front. Phys. (Beijing) 16, 64701 (2021), arXiv:2102.09222 [nucl-ex] .
  * Ryskin (1993) M. G. Ryskin, “Diffractive ${J}/\psi$ electroproduction in LLA QCD,” Z. Phys. C 57, 89–92 (1993).
  * Eskola  _et al._ (2022) Kari J. Eskola, Christopher A. Flett, Vadim Guzey, Topi Löytäinen, and Hannu Paukkunen, “Exclusive J/$\psi$ photoproduction in ultraperipheral Pb+Pb collisions at the CERN Large Hadron Collider calculated at next-to-leading order perturbative QCD,” Phys. Rev. C 106, 035202 (2022), arXiv:2203.11613 [hep-ph] .
  * Filip  _et al._ (2009) Peter Filip, Richard Lednicky, Hiroshi Masui, and Nu Xu, “Initial eccentricity in deformed ${}^{197}\mathrm{Au}+^{197}\mathrm{Au}$ and ${}^{238}\mathrm{U}+^{238}\mathrm{U}$ collisions at $\sqrt{{s}_{\mathit{NN}}}=200$ GeV at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 80, 054903 (2009).
  * Masui  _et al._ (2009) Hiroshi Masui, Bedangadas Mohanty, and Nu Xu, “Predictions of elliptic flow and nuclear modification factor from 200 GeV $\mathrm{U}+\mathrm{U}$ collisions at RHIC,” Phys. Lett. B 679, 440–444 (2009), arXiv:0907.0202 [nucl-th] .
  * Hirano  _et al._ (2013) Tetsufumi Hirano, Pasi Huovinen, Koichi Murase, and Yasushi Nara, “Integrated Dynamical Approach to Relativistic Heavy Ion Collisions,” Prog. Part. Nucl. Phys. 70, 108–158 (2013), arXiv:1204.5814 [nucl-th] .
  * Shen  _et al._ (2016) Chun Shen, Zhi Qiu, Huichao Song, Jonah Bernhard, Steffen Bass, and Ulrich Heinz, “The iEBE-VISHNU code package for relativistic heavy-ion collisions,” Comput. Phys. Commun. 199, 61–85 (2016), arXiv:1409.8164 [nucl-th] .
  * Schenke  _et al._ (2014) Bjoern Schenke, Prithwish Tribedy, and Raju Venugopalan, “Initial-state geometry and fluctuations in Au + Au, Cu + Au, and U + U collisions at energies available at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 89, 064908 (2014), arXiv:1403.2232 [nucl-th] .
  * Schenke  _et al._ (2020) Bjoern Schenke, Chun Shen, and Prithwish Tribedy, “Running the gamut of high energy nuclear collisions,” Phys. Rev. C 102, 044905 (2020), arXiv:2005.14682 [nucl-th] .
  * Xu  _et al._ (2021) Hao-jie Xu, Wenbin Zhao, Hanlin Li, Ying Zhou, Lie-Wen Chen, and Fuqiang Wang, “Probing nuclear structure with mean transverse momentum in relativistic isobar collisions,” (2021), arXiv:2111.14812 [nucl-th] .
  * Giacalone  _et al._ (2021) Giuliano Giacalone, Jiangyong Jia, and Chunjian Zhang, “Impact of Nuclear Deformation on Relativistic Heavy-Ion Collisions: Assessing Consistency in Nuclear Physics across Energy Scales,” Phys. Rev. Lett. 127, 242301 (2021), arXiv:2105.01638 [nucl-th] .
  * Zhang and Jia (2022) Chunjian Zhang and Jiangyong Jia, “Evidence of Quadrupole and Octupole Deformations in ${}^{96}\mathrm{Zr}+^{96}\mathrm{Zr}$ and ${}^{96}\mathrm{Ru}+^{96}\mathrm{Ru}$ Collisions at Ultrarelativistic Energies,” Phys. Rev. Lett. 128, 022301 (2022), arXiv:2109.01631 [nucl-th] .
  * Nijs and van der Schee (2021) Govert Nijs and Wilke van der Schee, “Inferring nuclear structure from heavy isobar collisions using Trajectum,” (2021), arXiv:2112.13771 [nucl-th] .
  * Cheng  _et al._ (2023) Yi-Lin Cheng, Shuzhe Shi, Yu-Gang Ma, Horst Stöcker, and Kai Zhou, “How does Bayesian analysis infer the nucleon distributions in isobar collisions?” (2023), arXiv:2301.03910 [nucl-th] .
  * Ryssens  _et al._ (2023) Wouter Ryssens, Giuliano Giacalone, Björn Schenke, and Chun Shen, “Evidence of Hexadecapole Deformation in Uranium-238 at the Relativistic Heavy Ion Collider,” (2023), arXiv:2302.13617 [nucl-th] .
  * Mäntysaari and Schenke (2016a) Heikki Mäntysaari and Björn Schenke, “Revealing proton shape fluctuations with incoherent diffraction at high energy,” Phys. Rev. D 94, 034042 (2016a), arXiv:1607.01711 [hep-ph] .
  * Mäntysaari and Schenke (2016b) Heikki Mäntysaari and Björn Schenke, “Evidence of strong proton shape fluctuations from incoherent diffraction,” Phys. Rev. Lett. 117, 052301 (2016b), arXiv:1603.04349 [hep-ph] .
  * Mäntysaari (2020) Heikki Mäntysaari, “Review of proton and nuclear shape fluctuations at high energy,” Rept. Prog. Phys. 83, 082201 (2020), arXiv:2001.10705 [hep-ph] .
  * Gelis  _et al._ (2010) Francois Gelis, Edmond Iancu, Jamal Jalilian-Marian, and Raju Venugopalan, “The Color Glass Condensate,” Ann. Rev. Nucl. Part. Sci. 60, 463–489 (2010), arXiv:1002.0333 [hep-ph] .
  * Mäntysaari  _et al._ (2022a) Heikki Mäntysaari, Björn Schenke, Chun Shen, and Wenbin Zhao, “Bayesian inference of the fluctuating proton shape,” Phys. Lett. B 833, 137348 (2022a), arXiv:2202.01998 [hep-ph] .
  * Alexa  _et al._ (2013) C. Alexa _et al._ (H1), “Elastic and Proton-Dissociative Photoproduction of $\mathrm{J}/\psi$ Mesons at HERA,” Eur. Phys. J. C 73, 2466 (2013), arXiv:1304.5162 [hep-ex] .
  * Mueller (2001) Alfred H. Mueller, “A Simple derivation of the JIMWLK equation,” Phys. Lett. B 523, 243–248 (2001), arXiv:hep-ph/0110169 .
  * Ebran  _et al._ (2012) J. P Ebran, E. Khan, T. Niksic, and D. Vretenar, “How atomic nuclei cluster,” Nature 487, 341–344 (2012).
  * Zhou  _et al._ (2016) E. F. Zhou, J. M. Yao, Z. P. Li, J. Meng, and P. Ring, “Anatomy of molecular structures in 20Ne,” Phys. Lett. B 753, 227–231 (2016), arXiv:1510.05232 [nucl-th] .
  * Marević  _et al._ (2018) P. Marević, J. P. Ebran, E. Khan, T. Nikšić, and D. Vretenar, “Quadrupole and octupole collectivity and cluster structures in neon isotopes,” Phys. Rev. C 97, 024334 (2018), arXiv:1802.02873 [nucl-th] .
  * Frosini  _et al._ (2022a) Mikael Frosini, Thomas Duguet, Jean-Paul Ebran, Benjamin Bally, Tobias Mongelli, Tomás R. Rodríguez, Robert Roth, and Vittorio Somà, “Multi-reference many-body perturbation theory for nuclei: II. Ab initio study of neon isotopes via PGCM and IM-NCSM calculations,” Eur. Phys. J. A 58, 63 (2022a), arXiv:2111.00797 [nucl-th] .
  * Bally  _et al._ (2023) B. Bally _et al._ , “Deciphering small system collectivity with bowling-pin-shaped 20Ne isotopes,” in preparation (2023).
  * Good and Walker (1960) M. L. Good and W. D. Walker, “Diffraction disssociation of beam particles,” Phys. Rev. 120, 1857–1860 (1960).
  * Miettinen and Pumplin (1978) Hannu I. Miettinen and Jon Pumplin, “Diffraction Scattering and the Parton Structure of Hadrons,” Phys. Rev. D 18, 1696 (1978).
  * Caldwell and Kowalski (2010) A. Caldwell and H. Kowalski, “Investigating the gluonic structure of nuclei via $\mathrm{J}/\psi$ scattering,” Phys. Rev. C 81, 025203 (2010).
  * Kowalski  _et al._ (2006) H. Kowalski, L. Motyka, and G. Watt, “Exclusive diffractive processes at HERA within the dipole picture,” Phys. Rev. D 74, 074016 (2006), arXiv:hep-ph/0606272 .
  * Hatta  _et al._ (2017) Yoshitaka Hatta, Bo-Wen Xiao, and Feng Yuan, “Gluon Tomography from Deeply Virtual Compton Scattering at Small-$x$,” Phys. Rev. D 95, 114026 (2017), arXiv:1703.02085 [hep-ph] .
  * Mäntysaari and Penttala (2022) Heikki Mäntysaari and Jani Penttala, “Complete calculation of exclusive heavy vector meson production at next-to-leading order in the dipole picture,” JHEP 08, 247 (2022), arXiv:2204.14031 [hep-ph] .
  * Mäntysaari and Penttala (2021) Heikki Mäntysaari and Jani Penttala, “Exclusive heavy vector meson production at next-to-leading order in the dipole picture,” Phys. Lett. B 823, 136723 (2021), arXiv:2104.02349 [hep-ph] .
  * Kovchegov and Levin (2022) Yuri V. Kovchegov and Eugene Levin,  _Quantum Chromodynamics at High Energy_, Cambridge Monographs on Particle Physics, Nuclear Physics and Cosmology (33), Vol. 33 (Cambridge University Press, 2022).
  * Schenke  _et al._ (2012a) Bjoern Schenke, Prithwish Tribedy, and Raju Venugopalan, “Fluctuating Glasma initial conditions and flow in heavy ion collisions,” Phys. Rev. Lett. 108, 252301 (2012a), arXiv:1202.6646 [nucl-th] .
  * Mäntysaari  _et al._ (2022b) Heikki Mäntysaari, Farid Salazar, and Björn Schenke, “Nuclear geometry at high energy from exclusive vector meson production,” Phys. Rev. D 106, 074019 (2022b), arXiv:2207.03712 [hep-ph] .
  * Mäntysaari  _et al._ (2021) Heikki Mäntysaari, Kaushik Roy, Farid Salazar, and Björn Schenke, “Gluon imaging using azimuthal correlations in diffractive scattering at the Electron-Ion Collider,” Phys. Rev. D 103, 094026 (2021), arXiv:2011.02464 [hep-ph] .
  * Mäntysaari and Schenke (2020) Heikki Mäntysaari and Björn Schenke, “Accessing the gluonic structure of light nuclei at a future electron-ion collider,” Phys. Rev. C 101, 015203 (2020), arXiv:1910.03297 [hep-ph] .
  * Mäntysaari  _et al._ (2019) Heikki Mäntysaari, Niklas Mueller, and Björn Schenke, “Diffractive Dijet Production and Wigner Distributions from the Color Glass Condensate,” Phys. Rev. D 99, 074004 (2019), arXiv:1902.05087 [hep-ph] .
  * Mäntysaari and Schenke (2018) Heikki Mäntysaari and Björn Schenke, “Confronting impact parameter dependent JIMWLK evolution with HERA data,” Phys. Rev. D 98, 034013 (2018), arXiv:1806.06783 [hep-ph] .
  * Schenke  _et al._ (2012b) Bjoern Schenke, Prithwish Tribedy, and Raju Venugopalan, “Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions,” Phys. Rev. C 86, 034908 (2012b), arXiv:1206.6805 [hep-ph] .
  * Blaizot and Traini (2022) Jean-Paul Blaizot and Marco Claudio Traini, “Effect of dipole size fluctuations on diffractive photo-production of vector mesons,” (2022), arXiv:2209.15545 [hep-ph] .
  * Moreland  _et al._ (2015) J. Scott Moreland, Jonah E. Bernhard, and Steffen A. Bass, “Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions,” Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Frosini  _et al._ (2022b) Mikael Frosini, Thomas Duguet, Jean-Paul Ebran, and Vittorio Somà, “Multi-reference many-body perturbation theory for nuclei: I. Novel PGCM-PT formalism,” Eur. Phys. J. A 58, 62 (2022b), arXiv:2110.15737 [nucl-th] .
  * Frosini  _et al._ (2022c) Mikael Frosini, Thomas Duguet, Jean-Paul Ebran, Benjamin Bally, Heiko Hergert, Tomás R. Rodríguez, Robert Roth, Jiangming Yao, and Vittorio Somà, “Multi-reference many-body perturbation theory for nuclei: III. Ab initio calculations at second order in PGCM-PT,” Eur. Phys. J. A 58, 64 (2022c), arXiv:2111.01461 [nucl-th] .
  * De Vries  _et al._ (1987) H. De Vries, C. W. De Jager, and C. De Vries, “Nuclear charge and magnetization density distribution parameters from elastic electron scattering,” Atom. Data Nucl. Data Tabl. 36, 495–536 (1987).
  * Loizides  _et al._ (2015) C. Loizides, J. Nagle, and P. Steinberg, “Improved version of the PHOBOS Glauber Monte Carlo,” SoftwareX 1-2, 13–18 (2015), arXiv:1408.2549 [nucl-ex] .
  * Carlson and Schiavilla (1998) J. Carlson and R. Schiavilla, “Structure and Dynamics of Few Nucleon Systems,” Rev. Mod. Phys. 70, 743–842 (1998).
  * Mäntysaari and Schenke (2017) Heikki Mäntysaari and Björn Schenke, “Probing subnucleon scale fluctuations in ultraperipheral heavy ion collisions,” Phys. Lett. B 772, 832–838 (2017), arXiv:1703.09256 [hep-ph] .
  * Cepila  _et al._ (2018) Jan Cepila, Jesus Guillermo Contreras, and Michal Krelina, “Coherent and incoherent $\mathrm{J/}\psi$ photonuclear production in an energy-dependent hot-spot model,” Phys. Rev. C 97, 024901 (2018), arXiv:1711.01855 [hep-ph] .
  * Pordes  _et al._ (2007) Ruth Pordes _et al._ , “The Open Science Grid,” J. Phys. Conf. Ser. 78, 012057 (2007).
  * Sfiligoi  _et al._ (2009) Igor Sfiligoi, Daniel C. Bradley, Burt Holzman, Parag Mhashilkar, Sanjay Padhi, and Frank Wurthwrin, “The pilot way to Grid resources using glideinWMS,” WRI World Congress 2, 428–432 (2009).