# Order and density fluctuations near the boundary in sheared dense
suspensions

Joia M. Miller  Daniel L. Blair  Jeffrey S. Urbach<EMAIL_ADDRESS>Department of Physics and Institute for Soft Matter Synthesis and Metrology,
Georgetown University, Washington, DC 20057.

(November 2, 2022)

###### Abstract

We introduce a novel approach to reveal ordering fluctuations in sheared dense
suspensions, using line scanning in a combined rheometer and laser scanning
confocal microscope. We validate the technique with a moderately dense
suspension, observing modest shear-induced ordering and a nearly linear flow
profile. At high concentration ($\phi=0.55$) and applied stress just below
shear thickening, we report ordering fluctuations with high temporal
resolution, and directly measure a decrease in order with distance from the
suspension’s bottom boundary as well as a direct correlation between order and
particle concentration. Higher applied stress produces shear thickening with
large fluctuations in boundary stress which we find are accompanied by
dramatic fluctuations in suspension flow speeds. The peak flow rates are
independent of distance from the suspension boundary, indicating that they
likely arise from transient jamming that creates solid-like aggregates of
particles moving together, but only briefly because the high speed
fluctuations are interspersed with regions flowing much more slowly,
suggesting that shear thickening suspensions possess complex internal
structural dynamics, even in relatively simple geometries.

##  I Introduction

Flowing dense suspensions of colloidal particles appear in a wide range of
important industrial processes, and the presence of flow can dramatically
modify the suspension microstructure, which in turn impacts flow properties
(reviewed in Vermant and Solomon (2005); Morris (2009); Lettinga (2016)). In
particular, the presence of simple shear modifies the well-understood phase
behavior of dense suspensions of nearly monodisperse colloidal particles in
equilibrium. Modest shear promotes layering which enhances crystallization,
while higher shear rates often disrupt crystalline order Chen et al. (1992);
Vermant and Solomon (2005); Holmqvist et al. (2005); Wu et al. (2009); Derks
et al. (2009); Richard and Speck (2015); Lettinga (2016). In addition, the
bulk response of the suspension can be very sensitive to the material close to
the confining boundaries Shereda et al. (2010); Cheng et al. (2012); Xu et al.
(2013); Mackay et al. (2014). Planar boundaries promote the formation of
layers of particles parallel to the boundaries Shereda et al. (2010); Villada-
Balbuena et al. (2022), enabling more efficient shearing as layers slide past
each other with reduced close particle interactions, and a resulting reduction
in the local viscosity Kulkarni and Morris (2009); Pieper and Schmid (2016);
Küçüksönmez and Servantie (2020). The layering also enhances ordering within
the layer, and in many circumstances the sliding layers show a high degree of
hexagonal ordering, similar to what is seen in two dimensional colloidal
crystals Myung et al. (2013); Mackay et al. (2014). In many dense suspensions,
the shear thinning arising from increased layering is followed by dramatic
shear thickening when the applied stress exceeds a critical value, with the
increase in viscosity attributable to a transition from primarily hydrodynamic
particle interactions at low stress to frictional interactions at high stress
(reviewed in Morris (2020)). The onset of frictional interactions will likely
disrupt layering Lee et al. (2018), a scenario confirmed by recent computer
simulations Goyal et al. (2022). A variety of evidence suggests that strong
shear thickening is accompanied by complex spatiotemporal dynamics, including
fluctuations in flow speed, local stress, and particle concentration Boersma
et al. (1991); Lootens et al. (2003); Nakanishi et al. (2012); Guy et al.
(2015); Nagahiro and Nakanishi (2016); Hermes et al. (2016); Rathee et al.
(2017); Saint-Michel et al. (2018); Ovarlez et al. (2020); Gauthier et al.
(2021); Maharjan et al. (2021); Rathee et al. (2020a, b, c, 2022). These
results highlight the need for methods that can probe the dynamics of dense
colloidal suspension with high spatial and temporal resolution. Here we
introduce a new approach to assessing dynamic ordering using a combined
rheometer and laser scanning confocal microscope Dutta et al. (2013), where
the laser scanning is limited to one spatial direction. In this case the scan
is perpendicular to the flow direction, and the suspension flow is primarily
responsible for the temporal evolution of the recorded intensity. We validate
the technique using a moderately dense suspension (volume fraction
$\phi=0.52$), where only modest shear-induced ordering is observed. At higher
concentration ($\phi=0.55$) and applied stress just below the onset of shear
thickening, we measure ordering fluctuations with high temporal resolution,
and directly measure the decrease in ordering with distance from the bottom
boundary of the suspension. We also observe a direct correlation between
ordering and particle concentration, with local regions of high order
corresponding to high particle concentrations, consistent with recent
observations from computer simulations Goyal et al. (2022). Higher applied
stress produces shear thickening, and we find that the large fluctuations in
boundary stress that underlie thickening are accompanied by dramatic
fluctuations in suspension flow rates. The peak flow rates are independent of
distance from the suspension boundary, indicating that they arise from
transient jamming that creates solid like aggregates of particles moving
together. Such aggregates must be short-lived because the high speed
fluctuations are interspersed with regions flowing much more slowly,
suggesting that shear thickening suspensions possess very complex internal
structural dynamics, even in relatively simple geometries.

##  II Approach:

For particulate suspensions, particle velocities are often extracted from
image sequences by particle tracking or correlation analysis (PIV) Besseling
et al. (2009). Particle tracking typically requires that the time between
images is small compared to the time for a particle to move an interparticle
separation (so that each particle can be unambiguously identified in
successive images). Thus for a frame rate $f$ and an interparticle spacing
$\delta$, the flow speed must satisfy $v_{max}^{imaging}<<f\delta$. The
typical small $\delta$ in dense suspensions requires either a vanishingly
small $v_{max}$ or extremely high frame rate $f$ for successful tracking. PIV
can, in principle, handle larger displacements between images, but only if the
particle configuration maintains a consistent distinctive pattern so that
there is a clear unique peak in the spatial cross-correlation between images,
which is typically not the case for the dense suspensions studied here.

We use an alternative approach here, inspired by fluid dynamics and machine
vision applications where the motion of material past point or line detectors
is used to infer speeds and structure Taylor (1938); He et al. (2017). We take
advantage of the rapid scan rate of a laser scanning confocal microscope,
where the point of focus is rapidly scanned back and forth in one direction,
which we will call the $x$ axis, with the normal rastering perpendicular to
the scan line in the focal plane disabled (Fig. 1A). We illustrate the
principle with a simple example: A fluorescent particle transiting the line
will produce a local intensity extremum, from which its $x$ position can be
determined. The time it takes the particle to transit the line can be measured
by the length of the ridge produced when a position-time intensity surface is
generated (Fig. 1B). If the center of the particle is in the plane of focus
and neglecting the effects of the finite optical resolution, the magnitude of
the component of the particle velocity perpendicular to the line can be
determined as $v=D/T$, where $D$ is the particle diameter and $T$ the transit
time. The requirement that the transit time is long compared to the time
between scans means that the flow speed must be $v_{max}^{line}<<sD$, where
$s$ is the scan rate and $D$ the particle diameter. Because the scan rate for
laser scanning confocal microscopy is 2-3 orders of magnitude larger than the
frame rate (the time to generate a 2D image is approximately the scan duration
times the number of lines in an image), and the interparticle spacing in dense
suspensions is comparable to the particle diameter, $\delta\sim D$,
$v_{max}^{line}>>v_{max}^{imaging}$. Figure 1C shows a representative
position-time intensity surface generate by our system (Leica SP5, scan rate
8000 Hz), where the passage of individual 1 $\mu$m diameter particles with
transit times $\sim 5$ ms indicating speeds of $\sim 200\mu m/s$ can be
clearly identified despite the high particle density $\phi\approx 0.55$. In
practice individual particles cannot be robustly identified at from the data,
and the position of the particles relative to the plane of focus is variable,
so we instead extract speeds from characteristic lengths extracted from
correlation analyses, as described below and detailed in Methods. The
technique works equally well for non-fluorescent particles in a fluorescent
suspending fluid, which is the approach we use in what follows.

Figure 1:  A) Schematic of linescan imaging of a sheared dense suspension. The
illumination laser rapidly scans back and forth, focused along the green line.
Particles transiting the line are illuminated (orange). The imaging region for
the usual 2D laser scan is indicated by the dashed line, and the arrows
represent the average velocity field of the sheared suspension. The depth of
focus (effective width of the image plane in $z$) depends on the optical
setup, and is $\sim 1\mu m$ here) B) Graphical representation of the linescan
data, where the continuous colored lines represent a small portion of each
individual scan, separated by a short interval. A fluorescent particle
transiting the scan region produces a local intensity maximum (marked with a
black line). C) Rendering of linescan data showing the passage of many
particles. The transit time can be determined by the length of the ridge
(black line). (Note that in what follows the solvent is fluorescent but the
particles are not, so the transit of particles produces local intensity
minima.)

The situation sketched above describes motion advection purely in the
direction of the shear flow, but the linescan measurement will also be
affected by other components of the flow velocity. Because of the geometry,
the speed in the flow direction will normally be of order $\dot{\gamma}z$,
where $\dot{\gamma}$ is the shear rate and $z$ is the distance in the gradient
direction from the bottom surface to the focus of the scan line (Fig. 1A). In
the cases considered here, the transverse velocity components are expected to
be much smaller, but it is important to keep this consideration in mind when
interpreting the linescan data.

A related approach has been used to measure flow profiles in capillary and
microfluidic systems (reviewed in Koynov and Butt (2012); Dong and Ren (2014))
where fluorescence correlation spectroscopy (FCS) is employed to extract
speeds using fluorescent single molecules or nanoparticles Kunst et al.
(2002). The exact shape of the intensity autocorrelation function depends on
the relative importance of diffusion and advection, but when diffusion is
negligible compared to advection (as is the case in the situations considered
here), the intensity autocorrelation function is a Gaussian with a width
determined by the flow speed and the size of the illumination volume Gösch et
al. (2000). By varying the scan orientation relative to the flow direction,
the full (average) three dimensional velocity vector can be determined Pan et
al. (2007, 2009). In those applications, however, the volume fraction of the
fluorescent species is always low, so they are reasonably assumed to be
distributed randomly. In the results presented below, spatial ordering
dramatically affects the spatiotemporal intensity fluctuations, complicating
the determination of the flow speed but revealing novel behavior in dense
colloidal shear flow.

Validation: To verify the accuracy of the linescan approach under these
conditions, we performed tests on a moderately dense suspension, volume
fraction $\phi=0.52$, where the viscosity is only weakly non-Newtonian, under
constant shear rate conditions. Figure 2 shows three representative
kymographs, each generated from 512 scans (65 ms total time), showing the
range of speeds from the slowest ($\dot{\gamma}=10s^{-1}$; $z=3\mu m$) to the
fastest ($\dot{\gamma}=30s^{-1}$; $z=20\mu m$). The wide range of transit
times is evident, as is the appearance of regions of hexagonal ordering, which
will be considered in the next section. Rather than identifying individual
particles, we can quantify the length and time scales of the intensity
variations using correlation analyses, with faster flow producing more rapid
fluctuations along the time axis. In particular, as detailed in Methods, we
calculate the speed of the flow by the rate of decay of the autocorrelation in
the time direction (vertical in Fig. 2). This produces a very precise
measurement of the fluctuation timescale, but conversion to speed requires
knowing the length scale associated with the decay in the density
autocorrelation. At the volume fractions considered here ($\phi\sim 0.5$) that
is of order the particle diameter, but the exact value is sensitive to the
local structure of the suspension. Thus, as discussed in the Methods section,
there is some uncertainty in the absolute value of the measured speeds, and
likely some correlation between the local order and the calculated speed. In
future work we plan to more fully characterize these effects with further
experimental tests and comparisons with computer simulations.

Figure 2:  Representative space-time kymographs generated from linescans
(horizontal) stacked vertically taken in a sheared dense suspension of silica
particles of diameter $1\mu m$, $\phi=0.52$, at different shear rates. The
spheres are not fluorescent, and a fluorescent dye is added to the solvent. A)
$\dot{\gamma}=10s^{-1}$, $z=3\mu m$, B) $\dot{\gamma}=20s^{-1}$, $z=6\mu m$,
C) $\dot{\gamma}=30s^{-1}$, $z=20\mu m$. Each image consists of 512 scans, for
an elapsed time of 65 ms.

We acquired a series of 900 kymographs at heights above the bottom boundary of
$z={3,6,10,20}\mu\rm{m}$, for constant applied shear rates
of$\dot{\gamma}={10,20,30}\rm{s}^{-1}$. The resulting flow profiles are
presented in Fig. 3A, with error bars given by the standard deviation of each
series. Figure 3B shows the same data, with the speed scaled by the speed of
the top plate, $v_{p}=\dot{\gamma}d$, where $d\approx 150\mu m$ is the width
of the rheometer gap at the imaging position. (For the cone-plate geometry
used here, the gap is a continuous function of radius, but changes relatively
little over the $\approx 65\mu m$ scan range.) The data collapse, with no free
parameters, to a unique, nearly linear profile that demonstrates the wide
range of speeds over which the technique provides precise speed measurements.
The deviation from an affine profile close to the wall likely arises from a
combination of wall slip, layering, and ordering (because the ordering is
strongest close to the wall and, as indicated above, spatial order will have
an effect on the calculated speed). Finally, we note that these measurements
are taken at relatively high flow rates, so Brownian Motion or particle
sedimentation are not significant on the timescale of the decay of the
intensity autocorrelation function (see Discussion).

Figure 3: A) Flow profiles for constant shear rate generated from linescan
data at different heights ($\phi=0.52$), at different constant shear rates. B)
Particle speeds scaled by the speed of the top plate of the rheometer,
$\dot{\gamma}d$, as a function of height above the bottom surface, scaled by
the rheometer gap. Error bars are standard deviations (N=900) . Dashed line
corresponds to purely affine flow with no wall slip.

##  III Results

Here we report the results of a dense suspension of 1 $\mu m$ silica spheres
at a volume fraction $\phi\approx 0.55$. The flow curve from bulk rheology for
this sample is shown in Fig. 4, and exhibits behavior similar to previous
reports of shear thickening in similar systems Morris (2020), and in
particular matches that of our previous results measuring local boundary
stress fluctuations of the same system Rathee et al. (2017). The viscosity
shows substantial shear thinning (viscosity decreasing as applied shear stress
is increased) until a critical stress of $\sigma_{c}\approx 20$ Pa, at which
point strong shear thickening is observed. The transition occurs at a shear
rate of $\dot{\gamma}_{c}\approx 10s^{-1}$, which corresponds to a Peclet
number $Pe=\dot{\gamma}r^{2}/D_{o}\approx 500$, where $D_{o}=k_{B}T/6\pi\eta
r$ is the diffusion coefficient for an isolated sphere of radius $r$ in a
solvent of viscosity $\eta$ Vermant and Solomon (2005); Morris (2009);
Lettinga (2016) (approximately 80 $mPa\cdot s$ for the 80:20 glycerol:water
mixture used here. Note that the high concentrations used here substantially
modify particle diffusion). Figure 4B shows shear rate vs. time for constant
applied stresses of 20, 50, and 100 Pa, corresponding to the conditions for
the kymographs reported below. The increasing shear rate fluctuations with
increasing stress is characteristic of a suspension showing strong shear
thickening Morris (2020).

Figure 4: A) Viscosity ($\CIRCLE$) and shear rate ($\blacksquare$ ) as a
function of applied shear stress for the suspension used in this study,
composed of silica particles of diameter $1.0\mu m$, $\phi\approx 0.55$. B)
Representative plots of shear rate vs. time for a constant applied stress of
20 (bottom), 50 (middle) and 100 (top) Pa. Specifically, the time series shown
correspond to the kymographs reported at $z=3\mu m$.

###  III.1 Order and density fluctuations before the onset of shear thickening

As summarized in Section I, the shear thinning observed for applied stress
$\sigma<\sigma_{c}$ is generally attributed to shear induced ordering, where
organization of the monodisperse spheres into layers, and hexagonal ordering
with the layers, reduces the viscosity by reducing collisions as the layers
slide past one another. Figure 5 shows representative kymographs at $\sigma$ =
20 Pa, at three different heights. At $z=3\mu m$ above the bottom surface, the
kymograph (top) reveals a very regular hexagonal lattice, elongated along the
time axis because the particle speed is relatively small. Defects in the
hexagonal order tend to be extended in the flow (time) directions, e.g. at the
top right of the panel. The left side of the kymograph shows a region of
relatively low order, which is uncommon at this height, as quantified below.
At $z=10\mu m$ (middle), regions of hexagonal order show up as small islands
in a sea of disorder, with a tendency to be be longer in the time direction
than in the spatial (vorticity) direction. Particle extent and separations in
the time direction are proportional to the flow speed, so comparing dimensions
requires rescaling by the flow speed (Fig. S1), or, more simply, counting the
number of particle diameters in the different directions. The interpretation
of the kymographs as direct images of 2D particle arrangements presumes that
the rate of particle rearrangement is relatively slow compared to the rate at
which structures traverse the scan region. The persistence of the hexagonal
order suggests that this is a reasonable approximation.

In addition to revealing the existence and shape of small islands of order,
the kymographs show that the ordered regions are darker than regions of
disorder. This can be seen clearly in the islands present at $z=10\mu m$,
where the overall intensity is lower in the islands than around them (Fig. 5,
middle), but is also evident at $z=3\mu m$ (top), where the small regions of
disorder are brighter than their surroundings. Because the fluorescent dye is
in the solvent, brighter regions correspond to a lower concentration of
colloidal particles. Below we quantify the correlation between order and low
intensity, and show that is substantial at all heights at 20 Pa applied
stress. This observation is consistent with simulations of dense suspensions
of monodisperse particles, where a correlation between local ordering and high
particle concentration was observed Goyal et al. (2022).

Figure 5:  Representative space-time kymographs generated from linescans
(vertical) stacked horizontally at 20 Pa constant applied stress at heights of
3$\mu m$ (top), 10 $\mu m$ (middle), and 20 $\mu m$ above the bottom surface.
Each image consists of 4096 scans, for an elapsed time of 520 ms. (Note that
the orientation is different than Fig. 2, with time on the horizontal axis.)

Quantifying the degree of hexagonal order in the kymographs is complicated by
several factors. The images generated by individual particles depend on the
height of the particle relative to imaging position. As discussed above, the
particles tend to order in layers, and the patterns will be different when the
center of the particle layer is at the linescan height compared to when the
space between the layers aligns with the scan. A further complication arises
from the fact that the 1 $\mu m$ diameter of the particles is only a little
larger than the resolution of the optical microscope, so particle images are
not well separated. Finally, the distortion of the image in the flow direction
depends on the speed of the flow, which fluctuates in time. Despite these
complications we have found that we can extract a robust quantitative measure
of the flow speed and order in each kymograph by calculating the two-
dimensional intensity autocorrelation function, $g(\vec{r})$, for the entire
kymograph, rescaling distances in the time direction by the flow speed, and
then calculating the degree of hexagonal order present at $|\vec{r}|=r_{o}$,
where $r_{o}$ is the first peak in the angle averaged correlation function
$g(r)$. This measure, described in more detail in Section V, produces a
complex scalar that is analogous to the quantity $C_{6}$ used in Morris
(2009). We use that name here, but note that because of the complications
described above, even a fully ordered kymograph will produce a $C_{6}$ with a
magnitude that is considerably less than unity (in what follows, $C_{6}$
refers to the magnitude of the complex order parameter).

Figure 6:  Section of time series of speed, hexagonal order, and intensity
generated from kymographs at applied stress 20 Pa and height 10 $\mu m$. The
first 8 data points correspond to the image shown in Fig. 5 (middle).

Figure 6 shows a portion of the timeseries of speed, $C_{6}$, and average
intensity generated from kymographs at $z=10\mu m$. The first eight points in
the time series correspond to the image shown in Fig. 5 (middle). The decrease
in order in the image (moving right to left) shows as a decrease in $C_{6}$
from $\approx 0.5$ to $<0.3$. Over the same region, the intensity increases
modestly but significantly, as evident both in the image and the graph. The
full timeseries shows that there is a reasonably consistent anticorrelation
between $C_{6}$ and intensity, consistent with the qualitative impression from
small scale variations in individual images, but that there are also
significant fluctuations that appear uncorrelated. The correlations can be
quantified by the Pearson correlation coefficient, and we find
C($C_{6}$,intensity) = $\\{-0.19,-0.37,-0.34,-0.58\\}$ at heights
$z=\\{3,6,10,20\\}\mu m$. The effects are particularly significant at the
intermediate heights, where the fluctuations in $C_{6}$ are largest. This
negative correlation with intensity indicates that the ordered regions have a
higher particle concentration, presumably as a consequence of more efficient
packing.

A modest positive correlation can be observed between order and speed, which
is consistent with an ordering-induced increase in local velocity close to the
wall, but is complicated by the fact that the correlation analysis used to
determine the flow speed likely has a weak dependence on order (see Methods).
Finally, we note that the orientation of the hexagonal order is such that
intensity peaks and valleys always align in the flow direction, corresponding
to a phase angle for the complex order parameter $\approx 0$ (data not shown).
By contrast, the order seen at $\phi=0.52$, while considerably weaker (Fig.
2), has a phase angle $\approx\pi/6$ (Figure S3).

Figure 7: Histograms of speed and order generated from kymographs at applied
stress 20 Pa (Each height: 3900 kymographs, 65ms each).

The change in flow speed and order with distance from the bottom surface can
be seen most clearly from histograms, as shown in Fig. 7. The speeds increase
with height, as expected, with a significant increase in the spread of speeds
that is roughly proportional to the speed increase. Specifically, the ratio of
the standard deviation to the mean speed is $\\{0.14,0.19,0.21,0.17\\}$ for
heights $z=\\{3,6,10,20\\}\mu m$, respectively. The slightly higher fractional
spread at 6 and 10 $\mu m$ likely arises from the coupling between order and
speed mentioned above. The change in order with height is much more dramatic,
showing uniformly high order $3\mu m$ above the bottom surface, occasional
disordered regions at $6\mu m$, roughly equal order and disorder at $10\mu m$,
and mostly disordered at $20\mu m$. This behavior is summarized in Fig. 8,
which shows the average and standard deviations for speed and order as a
function of height. Figure 8 also includes the data from 50 Pa applied stress,
which shows very similar trends. The increase in curvature of the flow profile
close to the boundary indicates a lower local viscosity. Boundary-induced
ordering and an associated viscosity decrease in dense suspensions has been
seen previously, but the high speed imaging approach used here enables us to
quantify spatio-temporal fluctuations in previously inaccessible regimes.

Figure 8: Profiles of average velocity and order for 20 and 50 Pa applied
stress. Error bars are standard deviations (N=3900).

###  III.2 Fluctuations associated with shear thickening

At higher applied stresses, the suspension shows substantial shear thickening
and the nature of the fluctuations change dramatically. We have previously
shown that shear thickening is associated with large fluctuations in stresses
at the boundary of the sheared suspension Rathee et al. (2017, 2020a, 2020c,
2020b, 2022), and specifically for the particles used here shear thickening is
associated with a proliferation of the regions of high stress that propagate
in the flow direction with approximately half the speed of the top plate
($v_{p}=\dot{\gamma}d$) Rathee et al. (2017). Here local stress measurements
at the bottom boundary are limited to a small region (62x62 $\mu m^{2}$) but
still show clustered spikes of high stress that propagate in the flow
direction, plotted in figure 9A. The intermittent spikes are similar in
magnitude and duration to those previously observed, but the clustering was
not reported previously. The connection with those and other results is
discussed below, but for what follows the important observation is that the
boundary stress shows regularly spaced clusters of intermittent high stress
spikes. Figure 9B shows an expanded view or one cluster of events, clearly
showing the presence of large discrete spikes in the boundary stress with
durations $\approx 100-500$ ms (consistent with all time series reported, the
data points are separated by 65 ms).

Figure 9: A) Boundary stress (component in the flow direction) for 100 Pa
applied stress. B) Expanded view of one cluster of events from data show in
(A).

Figure 10 shows the velocity time series for 4 different heights at 100 Pa
applied stress. The clustering of intermittent events is consistent in every
data set. (Each time series is a different measurement run, and is different
from the run that produced the boundary stress shown in Fig. 9.) Thus we can
be confident that the intermittent spikes in speed represent the flow
fluctuations that are associated with the intermittent high boundary stress
spikes. Interestingly, the spacing between the clusters, $\sim 65s$, is
roughly 1/2 of the rotation period for the rheometer tool, suggesting that the
clusters reflect stable features propagating in the flow direction with speed
$v_{p}/2$.

Figure 10: Full timeseries of speeds for 100 Pa applied stress at heights of
3, 6, 10, and 20 $\mu m$.

Figure 11 shows the time series for a cluster of spikes at each observed
height, with the time axis shifted so that in each set the spikes start at the
same time. The kymograph at the top of the figure, showing a single speed
spike at $z=3\mu m$, indicates that before the spike, the particles exhibit a
high degree of hexagonal ordering, consistent with the behavior observed at
this height at lower applied stress. The transition to high speed flow is
evidenced by a remarkably sharp boundary that extends primarily along the
vorticity direction, but with some meandering in the flow direction,
consistent with the fluctuations in boundary stress reported previously Rathee
et al. (2017). A complete absence of ordering and a significant increase in
average intensity marks the period of high speed flow, suggesting a region of
reduced particle concentration, although there are correlations between speed
and intensity introduced by the instrumentation, and further calibration will
be required to separate measurement effects from concentration when the speed
variation is large. The transition back to lower speed flow is also quite
rapid, but the hexagonal ordering takes longer to recover. These features are
consistent across the high speed events, as shown by the time series of speed,
$c_{6}$, and intensity at $3\mu m$ (Figure S2).

Figure 11: Timeseries of speed for 100 Pa applied stress at heights of 3, 6,
10, and 20 $\mu m$ during one burst of high speed events. The image shows a
kymograph from a single event. The initial times are shifted so the bursts
start at approximately the same time.

Also evident in Fig. 10 is that the speed differential between the spikes and
the background is smaller at larger $z$, indicating that the events represent
a more dramatic speed-up close to the bottom boundary. Finally, the clusters
appear to be more spread out at larger $z$. Since the periodicity is
independent of $z$, the propagation speed of the clusters should be
independent of $z$, so the spreading of the clusters with height suggests that
the spatial extent of the region of high stress increases with depth.

We can use the clustering of high speed events in the velocity time series to
quantify the region of “normal” flow, away from the events (the specific
segmentation is shown in Figure S3). Figure 12A shows the flow profile
obtained during those time periods, and the overall shape closely matches that
seen at lower stresses (Fig. 8). The agreement is made more evident by scaling
the measured speeds by the average shear rate as reported by the rheometer
(Fig. 12B), where all three profiles overlap, with no free parameters.
Something very different happens during the high speed events. We find
empirically that selecting the highest 60 speeds from the 3900 measurement
points provides a reasonable measure of the peak speeds during the high speed
events (Figure S4), and we find that the average of those top speeds is nearly
height independent (Fig. 12). Scaling the peak speeds by average the speed of
the top plate shows that the height-independent speed is on the order of
$v_{p}$, consistent with a solid jammed aggregate moving with the speed of the
top plate, perhaps with some slip. (Note however, as discussed below, that the
high speed events are quite likely associated with substantial non-affine
flows, and thus it is possible that flows in the gradient direction contribute
to the decay of the correlation and thus to the calculated speed.)  Imaging of
the bottom layer of particles confirms that the the high speed fluctuations
are associated with large slip of that layer. Those apparent aggregates are
intertwined with regions of the suspension where the particles are flowing
much slower (Fig. 11), which is inconsistent with a single solid aggregate.
These results provide the first direct measurements of the structure of the
velocity field inside of high stress fluctuations in shear thickening
suspensions.

Figure 12: A) Average speed at different heights from portions of the time
series away from the high speed events, and average speed for the highest
speeds in each time series (N=60). B) Data in (A) scaled by the speed of the
top plate, with the scaled profiles from 20 and 50 Pa for comparison (Fig. 8A)
.

##  IV Discussion

Order and concentration fluctuations before the onset of shear thickening: The
dynamics revealed by high speed linescan imaging presented above are
consistent with previous observations of boundary induced ordering in sheared
dense monodisperse suspensions, but provide a more detailed picture of the
spatiotemporal dynamics and demonstrate the presence of a strong connection
between ordering and concentration fluctuations. Planar boundaries enhance
layering in colloidal suspensions, and the layering can reduces viscous
dissipation, producing shear banding Mackay et al. (2014); Küçüksönmez and
Servantie (2020); Goyal et al. (2022). Hexagonal ordering of monodisperse
particles within layers of sheared suspension has also been seen in a variety
of circumstances, starting with the seminal work of Ackerson and coworkers.
Ackerson (1990); Chen et al. (1992). Particularly relevant to this work, an
imaging study by Wu et al. observed fluctuation in hexagonal ordering during
the process of shear-induced melting in a confined suspension arising from the
nucleation of localized domains that temporarily lost and regained their
ordered structure Wu et al. (2009). Similar behavior has been reproduced in
computer simulations Mackay et al. (2014); Myung et al. (2013). Recent
simulation results reported that defects in crystalline order are associated
with local decrease in particle concentration Goyal et al. (2022). It is
perhaps not surprising that the interplay between shear flow, crystal
nucleation, and shear-induced crystal breakup produces a complex phase diagram
with complex dynamics Lettinga (2016). Here we show that for a monodisperse
suspension at high packing fraction and high Peclet number ($Pe\sim 500$), the
transition from mostly ordered to mostly disordered occurs $\sim 10$ particle
diameters from the boundary (Fig. 8), and is associated with rapid temporal
fluctuations in order, speed, and concentration (Fig. 6). The concentration
fluctuations are particularly significant, because the local increase in
particle concentration seen in ordered domains requires an compensating
decrease in solvent concentration, which implies the existence of local fluid
migration. These results indicate that an accurate quantitative model for
ordering fluctuations in sheared monodisperse suspensions at high Peclet
number must include the flow and pore pressure fields responsible for relative
flow between the particulate and fluid phases Meer (2016); O’Neill et al.
(2019); Maharjan et al. (2021); Rathee et al. (2022).

Fluctuations associated with high local stresses during shear thickening: The
high speed imaging approach employed here has revealed that localized high
boundary stresses are accompanied by large rapid speed increases near the
boundary, as well as a loss of order (Fig. 11). Recent simulation results
showed that the transition from hydrodynamic to frictional interparticle
interactions that is believed to underlie strong shear thickening Morris
(2020) is associated with a disruption of layering and ordering Goyal et al.
(2022). It seems likely that we are observing a similar transition, from
relatively low stress, lubricated particle interactions producing layered,
ordered low viscosity flow, to high stress, frictional particle interactions
producing high viscosity disordered flow in the regions that produce high
boundary stress. This transition appears to be remarkably sharp, with a
boundary that is only a few particles wide (Fig. 11, top). Also remarkable is
the observation that the particle speeds during these events are independent
of depth (Fig. 12). This suggests that the particles are moving together, as a
jammed aggregate, similar to the model proposed to explain propagating high
normal stresses observed in cornstarch suspensions Gauthier et al. (2021). It
is important to note, however, that the linescan measurements performed here
do not distinguish between flow components perpendicular to the scan direction
(i.e. in the flow or gradient directions). As discussed above, under usual
flow conditions the component of the velocity in the flow direction, of order
$\dot{\gamma}z$, will be large compared to any non-affine flows, but the
transient jamming observed during shear thickening is clearly unusual, so more
complicated large-scale non-affine flows are likely. Further measurements, for
example addinge measurements with linescans oriented in the flow direction,
are needed to get a clearer picture of the flow profile during the high stress
events.

More generally, measurements of nearly affine flow away from the localized
high stresses, both in this system and in cornstarch Rathee et al. (2022),
combined with the rapid fluctuation in velocity measured during the clusters
of events (Fig. 11), suggests a complicated flow field. Furthermore, although
instrumental effect preclude a direct measurement of particle concentration
during the high speed events, the longstanding connection between frictional
interactions and dilatancy Morris (2020), evidenced by our observation of
fluid migration associated with high stress fluctuations in cornstarch
suspensions Rathee et al. (2022), suggests that relative flow between the
particulate and fluid phases likely plays an important role in the observed
speed fluctuations.

These observations contribute to a growing body of evidence indicating that
the shear thickening transitions typically involved complex spatiotemporal
dynamics with structures propagating in the flow direction, including dilatant
fronts Nakanishi et al. (2012); Nagahiro and Nakanishi (2016), local
deformations of the air-sample interface at the edge of the rheometer tool
Hermes et al. (2016); Maharjan et al. (2021), local normal stresses in sheared
cornstarch Gauthier et al. (2021), concentration fluctuations appearing as
periodic waves moving in the direction of flow Ovarlez et al. (2020), and high
shear stress at the suspension boundary Rathee et al. (2017, 2020b, 2020a,
2022). The specifics of the dynamics vary considerably, presumably indicating
a sensitivity to the details of the suspension and the measurement geometry.
This sensitivity is perhaps not surprising, given that shear thickening
involves instabilities that can produce discontinuities in material
parameters. The location of those discontinuities in a uniform extended system
represent a broken symmetry, and thus in any physical realization will be very
sensitive to the parameter variations (e.g. of shear rate, rheometer gap,
distance from suspension boundary) that are present in all experimental
systems.

Here we have presented an initial application of a powerful new approach that
reveals spatiotemporal dynamics of sheared dense suspensions with high spatial
and temporal resolution. With further testing and validation the approach has
the potential to provide accurate and precise measurements of local speed,
structure, and particle concentration, and should provide a new avenue to
answer open questions about dense suspensions under flow.

##  V Methods

All experimental suspensions were composed of $0.9\mu m$ silica beads (Bang’s
Lab) in an index-matched (80/20 v/v) glycerol/water mixture. For imaging
purposes, fluorescein sodium salt was added to the suspension so that the
unlabeled spheres could be imaged as dark spots in a fluorescent background.
Rheological measurements were performed on a stress-controlled rheometer
(Anton Paar MCR 301) mounted on an inverted confocal (Leica SP5) microscope
Dutta et al. (2013) using a cone-plate geometry with a diameter 25 mm. The
linescan data was acquired with a 60X objective at a radius of 2/3 of the
plate diameter, where the rheometer gap is $d\approx 145\mu m$.
Linescan analysis: A typical data set is composed of a 1-2 million scans, with
the scan direction aligned along the vorticity axis (perpendicular to the flow
and the gradient). For visualization, we typically divide the set into a
series of 2D arrays, with one dimension (horizontal) given by the number of
pixels per scan (here 1024), and the other dimension (vertical) given by the
chosen number of scans per image, in this case 512 scans. Fig. 14A shows an
example image generated from $1\mu$m diameter non-fluorescent spheres in a
fluorescent background, 10 $\mu m$ above the bottom of the sheared suspension,
similar to those shown in Fig.5. Individual particles show up as dark ovals.
In principle, the vertical axis of each oval could be used to measure the
speed of each particle. In practice, in many places identifying individual
particles is challenging, and we have found correlation analyses more
reliable. (This arises in part because we have no control of the position of
particles relative to the plane of focus, and the image will include
contributions from particles that are $\approx 1\mu$m above or below the focal
plane). Note also that the orientation of the ovals provides a measure for the
component of the flow in the vorticity direction. Here that component is
always small, and is not included in our analysis. An alternative approach,
not employed here, would be to orient the scan line along the flow direction.
This would enable very precise measurements of flow speed, as particles move
along the line, but would not reveal spatial structure in the flow-vorticity
plane.

One approach to extract a characteristic transit time from the linescan
intensity data, $I(x,t)$, is to calculate the autocorrelation in the time
direction, $g(\Delta t)=<\delta I(x,t)\delta I(x,t+\Delta t)>_{x,t}/<\delta
I(x,t)^{2}>_{x,t}$, where $\delta I(x,t)=I(x,t)-<I(x,t)>_{x,t}$. The range of
$x$ and $t$ to include in the averages can be varied depending on the spatial
and temporal resolution required by the effective particle size. Slower shear
rates require averaging over larger time windows so as to capture at least one
full particle transit. For the conditions in this study we found that 512
scans provides enough data for robust correlation analysis while still
allowing adequate temporal resolution. Figure 13 shows $g(\Delta t)$ for the
three images displayed in Fig. 2.

Figure 13:  Normalized intensity autocorrelation function $g(\Delta t)$ as
defined in Methods for the three kymographs shown in Fig. 2,
$\dot{\gamma}=10s^{-1}$ , $z=3\mu m$ ($\CIRCLE$), $\dot{\gamma}=20s^{-1}$ ,
$z=6\mu m$ ($\blacksquare$ ), and $\dot{\gamma}=30s^{-1}$, $z=20\mu m$
($\blacklozenge$). The width of the central peak and the position of the first
minimum are proportional to the flow speed, for fixed positional
configurations, but both measures are have a some sensitivity to the presence
of positional order.

The location of the first minimum in the autocorrelation, $t_{min}$ can be
precisely identified algorithmically, and provides a measure that is
proportional to the speed of the flow. We expect that
$v_{flow}=l_{flow}/2t_{min}$, where $l_{flow}$ is approximately equal to the
spacing between particles in the flow direction. More precisely, $l_{flow}$
should be equal to the first minimum in the average instatntaneous density
autocorrelation calculated along the flow axis. Unfortunately that quantity is
unknown. Assuming it does not change with time, $v_{flow}\propto 1/t_{min}$.
Similarly, the minimum of the autocorrelation in the space direction provides
a direct measure of the spacing in the $x$ (vorticity) direction, specifically
$l_{vorticity}=2x_{min}.$ We find that this measure does not depend directly
on the flow speed, but is quite sensitive to local hexagonal ordering. In part
because of this sensitivity, and in part because we are interested in
quantitative measures of hexagonal order, we have instead employed a slightly
more complex analysis approach based on the 2D autorcorrelation function, but
we have found that the speed variations are essentially indistinguishable from
those calculated with the 1D correlation analysis.

Using a 2D correlation analysis, we can quantify the hexagonal order that is
evident in Fig. 14A. It is important to remember that the image is not a
snapshot of the two-dimensional particle arrangement, as the top of the image
is generated at an earlier time than the bottom. However, if the structure of
the suspension evolves slowly compared to the time to generate an image, it
will in fact be an accurate representation of the spatial arrangement of the
particles, with the vertical separations expanded (or compressed) by a factor
$\alpha=(s\cdot d)/v_{flow}$, where $s$ is the scan rate, $d$ is the actual
spatial separation, and $v$ is the speed of the flow. Figure 14B shows a
contour plot of the two-dimensional spatial autocorrelation, $G(\vec{r})$, of
the image shown in A. The sixfold symmetry is evident, as is the expansion in
the vertical direction. Note that $G(\vec{r})$ is large along the horizontal
axis ($y=0)$. This is likely a consequence of fluctuations in the laser
intensity on a sub-millisecond timescale. These fluctuations can be filtered
out, but for our purposes it is sufficient to simply exclude the horizontal
axis from our analysis of $G(\vec{r})$. Fitting the central peak of
$G(\vec{r})$ (highlighted region in Fig. 14B) to a 2D Gaussian provides an
accurate measure of the width in the flow direction, $\Delta y$ (number of
scan lines) which in inversely proportional to the speed,
$v_{flow}=\delta\cdot s/\Delta y$. Determining the scale factor $\delta$, the
width of the central peak in $G(\vec{r})$ in physical units, is a significant
source of uncertainty in the measurement. The particle diameter is known
precisely, but $\delta$ also depends on the particle separation, and is also
somewhat sensitive to the degree of order. We have estimated $\delta$ from
looking at $G(x)$ (the correlation along the scan direction) in flows with no
detectable order, where we find that the peak width is 0.2 times the distance
to the first peak in $G(x)$. Assuming the first peak in G(x) is at a location
$d_{peak}=d_{particle}/\phi^{-1/3}\approx 1.2*d_{particle}$, we find
$\delta\approx 1.2d_{particle}/5$. This conversion enables us to provide speed
measurements in physical units, but further tests will be required to validate
the accuracy of the results. For the measurements described here, the $\sim
10\%$ uncertainty in the exact value of the proportionality constant does not
impact the conclusions, but systematic effects of ordering on $\Delta y$ may
impact some results, e.g. the apparent correlation between order and speed in
the measurements at 20 Pa.

Using the measured width $\Delta y$, we can rescale the flow direction to
produce a symmetric $G(\vec{r})$. Figure 14C shows a contour plot of the
rescaled $G(\vec{r})$, and thus represents our measure of the average spatial
structure of the imaged region. In the results presented below, we are
primarily interested in measuring the degree of hexagonal order in the
suspension. We can quantify this by measuring the intensity as a function of
angle at a distance $r_{NN}$ from the origin, where $r_{NN}$ is the position
of the first maximum in $G(r)$, approximately the average nearest neighbor
separation. The points included in the calculation of are highlighted in Fig.
14, and the resulting data $G(r_{NN},\theta)$ is shown in Figure 14D. A
complex scaler $C_{6}$ representing the hexagonal order can be calculated
according to

$C_{6}=\frac{\int_{0}^{2\pi}G(r_{NN},\theta)\exp^{-i6\theta}d\theta}{\int_{0}^{2\pi}G(r_{NN},\theta)d\theta}.$
(1)

Figure 14:  A) Example kymograph showing showing moderate hexagonal ordering.
B) Two-dimensional autocorrelation of (A), with the central region used to
determine the peak widths highlighted. C) 2D autocorrelation with vertical
distances rescaled so that the central peak is circularly symmetric. The ring
represents points whose distance from the origin is equal to the expected
average inter-particle spacing. D) Intensity as a function of angle for the
ring of points identified in (C), resulting in a $C_{6}$ with magnitude 0.36
and angle 0.03 radians.

Boundary stress microscopy. BSM measurements (Fig. 9 were performed as
described in Rathee et al. (2017, 2020b), but employing elastic films of
relatively high modulus ($G\sim 1MPa$) to minimize possible effects of
boundary compliance. Briefly, elastic films thickness 50 $\pm 3\mu$m were
deposited by spin coating PDMS (Sylgard 184; Dow Corning) and a curing agent
on 40 mm diameter glass cover slides (Fisher Sci) that were cleaned thoroughly
by plasma cleaning and rinsing with ethanol and deionized water Rathee et al.
(2017). After deposition of PDMS, the slides were cured at 85 oC for 2 hours.
After curing, the PDMS was functionalized with 3-aminopropyl triethoxysilane
(Fisher Sci) using vapor deposition for 40 min. Carboxylate-modified
fluorescent spherical beads of radius 0.5 $\mu$m with excitation/emission at
520/560 nm were attached to the PDMS surface. Before attaching the beads to
functionalized PDMS, the beads were suspended in a solution containing PBS
solution (Thermo-Fisher). The concentration of beads used was 0.006 $\%$
solids. The surface stresses at the interface were calculated using an
extended traction force technique and codes given in ref. Style et al. (2014).
Taking the component of the surface stress in the flow (velocity) direction,
we obtain the scalar field $\sigma_{BSM}(\vec{r},t)$, representing the
spatiotemporally varying surface stress. The data shown in Fig. 9 are a
spatial average of each field, $<\sigma_{BSM}(\vec{r},t)>_{t}$.

##  VI Acknowledgements

The authors thank Peter Olmsted and Emanuela Del Gado for helpful discussions.
This work was supported by the National Science Foundation (NSF) under Grant
No. DMR-1809890. J. S. U. is supported, in part, by the Georgetown
Interdisciplinary Chair in Science Fund.

## References

  * Vermant and Solomon (2005) J. Vermant and M. J. Solomon, Journal of Physics: Condensed Matter 17, R187 R216 (2005), ISSN 0953-8984.
  * Morris (2009) J. F. Morris, Rheologica Acta 48, 909 (2009), ISSN 0035-4511.
  * Lettinga (2016) M. P. Lettinga, pp. 81–110 (2016).
  * Chen et al. (1992) L. B. Chen, C. F. Zukoski, B. J. Ackerson, H. J. M. Hanley, G. C. Straty, J. Barker, and C. J. Glinka, Physical Review Letters 69, 688 (1992), ISSN 0031-9007.
  * Holmqvist et al. (2005) P. Holmqvist, M. P. Lettinga, J. Buitenhuis, and J. K. G. Dhont, Langmuir 21, 10976 (2005), ISSN 0743-7463.
  * Wu et al. (2009) Y. L. Wu, D. Derks, A. v. Blaaderen, and A. Imhof, Proceedings of the National Academy of Sciences of the United States of America 106, 10564 (2009), ISSN 0027-8424.
  * Derks et al. (2009) D. Derks, Y. L. Wu, A. v. Blaaderen, and A. Imhof, Soft Matter 5, 1060 (2009), ISSN 1744-683X.
  * Richard and Speck (2015) D. Richard and T. Speck, Scientific Reports 5, 14610 (2015), eprint 1509.08687.
  * Shereda et al. (2010) L. T. Shereda, R. G. Larson, and M. J. Solomon, Physical Review Letters 105, 228302 (2010), ISSN 0031-9007.
  * Cheng et al. (2012) X. Cheng, X. Xu, S. A. Rice, A. R. Dinner, and I. Cohen, Proceedings of the National Academy of Sciences 109, 63 (2012), ISSN 0027-8424.
  * Xu et al. (2013) X. Xu, S. A. Rice, and A. R. Dinner, The Journal of Physical Chemistry Letters 4, 3310 (2013), ISSN 1948-7185, eprint 1310.1331.
  * Mackay et al. (2014) F. E. Mackay, K. Pastor, M. Karttunen, and C. Denniston, Soft Matter 10, 8724 (2014), ISSN 1744-683X.
  * Villada-Balbuena et al. (2022) A. Villada-Balbuena, G. Jung, A. B. Zuccolotto-Bernez, T. Franosch, and S. U. Egelhaaf, Soft Matter 18, 4699 (2022), ISSN 1744-683X.
  * Kulkarni and Morris (2009) S. D. Kulkarni and J. F. Morris, Journal of Rheology 53, 417 (2009), ISSN 0148-6055.
  * Pieper and Schmid (2016) S. Pieper and H.-J. Schmid, Journal of Non-Newtonian Fluid Mechanics 234, 1 (2016), ISSN 0377-0257.
  * Küçüksönmez and Servantie (2020) E. Küçüksönmez and J. Servantie, Physical Review E 102, 012604 (2020), ISSN 2470-0045.
  * Myung et al. (2013) J. S. Myung, S. Song, and K. H. Ahn, Journal of Non-Newtonian Fluid Mechanics 199, 29 (2013), ISSN 0377-0257.
  * Morris (2020) J. F. Morris, Annual Review of Fluid Mechanics 52, 121 (2020), ISSN 0066-4189.
  * Lee et al. (2018) J. Lee, Z. Jiang, J. Wang, A. R. Sandy, S. Narayanan, and X.-M. Lin, Physical Review Letters 120 (2018), ISSN 0031-9007, not sure what this means.
  * Goyal et al. (2022) A. Goyal, E. D. Gado, S. Z. Jones, and N. S. Martys, Journal of Rheology 66, 1055 (2022), ISSN 0148-6055.
  * Boersma et al. (1991) W. H. Boersma, P. J. M. Baets, J. Laven, and H. N. Stein, Journal of Rheology 35, 1093 (1991), ISSN 0148-6055.
  * Lootens et al. (2003) D. Lootens, H. Damme, and P. Hébraud, Physical Review Letters 90, 178301 (2003), ISSN 1079-7114.
  * Nakanishi et al. (2012) H. Nakanishi, S.-i. Nagahiro, and N. Mitarai, Physical Review E 85, 011401 (2012), ISSN 1539-3755.
  * Guy et al. (2015) B. M. Guy, M. Hermes, and W. C. K. Poon, Physical Review Letters 115, 088304 (2015), ISSN 0031-9007.
  * Nagahiro and Nakanishi (2016) S.-I. Nagahiro and H. Nakanishi, Physical review. E 94, 062614 (2016), ISSN 2470-0053.
  * Hermes et al. (2016) M. Hermes, B. M. Guy, W. C. K. Poon, G. Poy, M. E. Cates, and M. Wyart, Journal of Rheology 60, 905 (2016), ISSN 0148-6055, eprint 1511.08011.
  * Rathee et al. (2017) V. Rathee, D. L. Blair, and J. S. Urbach, Proceedings of the National Academy of Sciences 114, 8740 (2017), ISSN 0027-8424.
  * Saint-Michel et al. (2018) B. Saint-Michel, T. Gibaud, and S. Manneville, Physical Review X 8, 031006 (2018), ISSN 2160-3308.
  * Ovarlez et al. (2020) G. Ovarlez, A. V. N. Le, W. J. Smit, A. Fall, R. Mari, G. Chatté, and A. Colin, Science Advances 6, eaay5589 (2020), ISSN 2375-2548.
  * Gauthier et al. (2021) A. Gauthier, M. Pruvost, O. Gamache, and A. Colin, Journal of Rheology 65, 583 (2021), ISSN 0148-6055.
  * Maharjan et al. (2021) R. Maharjan, E. O’Reilly, T. Postiglione, N. Klimenko, and E. Brown, Physical Review E 103, 012603 (2021), ISSN 2470-0045, eprint 2004.14316.
  * Rathee et al. (2020a) V. Rathee, D. L. Blair, and J. S. Urbach, Soft Matter 17, 1337 (2020a), ISSN 1744-683X.
  * Rathee et al. (2020b) V. Rathee, D. L. Blair, and J. S. Urbach, Journal of Rheology 64, 299 (2020b), ISSN 0148-6055.
  * Rathee et al. (2020c) V. Rathee, S. Arora, D. L. Blair, J. S. Urbach, A. K. Sood, and R. Ganapathy, Physical Review E 101, 040601 (2020c), ISSN 2470-0045, eprint 1906.06356.
  * Rathee et al. (2022) V. Rathee, J. Miller, D. L. Blair, and J. S. Urbach, Proceedings of the National Academy of Sciences 119, e2203795119 (2022), ISSN 0027-8424.
  * Dutta et al. (2013) S. K. Dutta, A. Mbi, R. C. Arevalo, and D. L. Blair, Rev Sci Instrum 84, 063702 (2013).
  * Besseling et al. (2009) R. Besseling, L. Isa, E. R. Weeks, and W. C. Poon, Advances in Colloid and Interface Science 146, 1 (2009), ISSN 0001-8686.
  * Taylor (1938) G. I. Taylor, Proceedings of the Royal Society of London. Series A - Mathematical and Physical Sciences 164, 476 (1938), ISSN 0080-4630.
  * He et al. (2017) G. He, G. Jin, and Y. Yang, Annual Review of Fluid Mechanics 49, 51 (2017), ISSN 0066-4189.
  * Koynov and Butt (2012) K. Koynov and H.-J. Butt, Current Opinion in Colloid & Interface Science 17, 377 (2012), ISSN 1359-0294.
  * Dong and Ren (2014) C. Dong and J. Ren, ELECTROPHORESIS 35, 2267 (2014), ISSN 0173-0835.
  * Kunst et al. (2002) B. H. Kunst, A. Schots, and A. J. W. G. Visser, Analytical Chemistry 74, 5350 (2002), ISSN 0003-2700.
  * Gösch et al. (2000) M. Gösch, H. Blom, J. Holm, T. Heino, and R. Rigler, Analytical Chemistry 72, 3260 (2000), ISSN 0003-2700.
  * Pan et al. (2007) X. Pan, H. Yu, X. Shi, V. Korzh, and T. Wohland, Journal of Biomedical Optics 12, 014034 (2007), ISSN 1083-3668.
  * Pan et al. (2009) X. Pan, X. Shi, V. Korzh, H. Yu, and T. Wohland, Journal of Biomedical Optics 14, 024049 (2009), ISSN 1083-3668.
  * Ackerson (1990) B. J. Ackerson, Journal of Rheology 34, 553 (1990), ISSN 0148-6055.
  * Meer (2016) D. v. d. Meer, Annual Review of Fluid Mechanics 49, 1 (2016), ISSN 0066-4189.
  * O’Neill et al. (2019) R. E. O’Neill, J. R. Royer, and W. C. K. Poon, Physical Review Letters 123, 128002 (2019), ISSN 0031-9007, eprint 1808.09950.
  * Style et al. (2014) R. W. Style, R. Boltyanskiy, G. K. German, C. Hyland, C. W. MacMinn, A. F. Mertz, L. A. Wilen, Y. Xu, and E. R. Dufresne, Soft Matter 10, 4047 (2014), ISSN 1744-683X.


11institutetext:  Instituto de Astrofísica, Universidad Andrés Bello,
Fernández Concha 700, Las Condes, Santiago, Chile
11email<EMAIL_ADDRESS>
#  Spatially resolved mid-infrared observations of the circumstellar
environment of the born-again object FG Sge ††thanks:  Based on observations
made with the Very Large Telescope Interferometer of the European Southern
Observatory. Program ID: 081.D-0244

K. Ohnaka 11 B. A. Jara Bravo 11

(Received / Accepted )

###### Abstract

Context.  FG Sge has evolved from the hot central star of the young planetary
nebula Hen 1-5 to a G–K supergiant in the last 100 years. It is one of the
three born-again objects that has been identified as of yet, and they are
considered to have undergone a thermal pulse in the post-asymptotic giant
branch evolution.

Aims.  We present mid-infrared spectro-interferometric observations of FG Sge
and probe its dusty environment.

Methods.  FG Sge was observed with MIDI at the Very Large Telescope
Interferometer at baselines of 43 and 46 m between 8 and 13 $\mu$m.

Results.  The circumstellar dust environment of FG Sge was spatially resolved,
and the Gaussian fit to the observed visibilities results in a full width at
half maximum of $\sim$10.5 mas. The observed mid-infrared visibilities and the
spectral energy distribution can be fairly reproduced by optically thick
($\tau_{V}$ $\approx 8$) spherical dust shell models consisting of amorphous
carbon with an inner radius $r_{\rm in}$ of $\sim$30 $R_{\star}$
(corresponding to a dust temperature of $1100\pm 100$ K). The dust shell is
characterized with a steep density profile proportional to $r^{-3.5\pm 0.5}$
from the inner radius $r_{\rm in}$ to $(5-10)\times\mbox{$r_{\rm in}$}$,
beyond which it changes to $r^{-2}$. The dust mass is estimated to be
$\sim\\!7\times 10^{-7}$ $M_{\sun}$, which translates into an average total
mass-loss rate of $\sim\\!9\times 10^{-6}$ $M_{\sun}$ yr-1 as of 2008 with a
gas-to-dust ratio of 200 being adopted. In addition, the 8–13 $\mu$m spectrum
obtained with MIDI with a field of view of 200 mas does not show a signature
of the polycyclic aromatic hydrocarbon (PAH) emission, which is in marked
contrast to the spectra taken with the Spitzer Space Telescope six and 20
months before the MIDI observations with wide slit widths of
3$\aas@@fstack{\prime\prime}$6–10″. This implies that the PAH emission
originates from an extended region of the optically thick dust envelope.

Conclusions.  The dust envelope of FG Sge is much more compact than that of
the other born-again stars’ Sakurai’s object and V605 Aql, which might reflect
the difference in the evolutionary status. The PAH emission from the extended
region of the optically thick dust envelope likely originates from the
material ejected before the central star became H-deficient, and it may be
excited by the UV radiation from the central star escaping through gaps among
dust clumps and/or the bipolar cavity of a disk-like structure.

###### Key Words.:

infrared: stars – techniques: interferometric – stars: circumstellar matter –
stars: carbon – stars: AGB and post-AGB – stars: individual: FG Sge

††offprints: K. Ohnaka

##  1 Introduction

FG Sge is the central star of the young planetary nebula Hen 1-5 and a very
rare object whose evolution can be observed in a lifetime of a human being. In
the last 100 years, it has rapidly evolved redward on the Hertzsprung-Russell
(H-R) diagram, from the hot central star of the planetary nebula (effective
temperature $T_{\rm eff}\approx 40\,000$ K) to a G supergiant ($T_{\rm
eff}\approx 5500$ K, Jeffery & Schönberner jeffery06 (2006)). The rapid
evolution of FG Sge is considered to be caused by a late thermal pulse (LTP)
that occurs after the asymptotic giant branch (AGB), when the star is evolving
blueward and almost horizontally on the H-R diagram. The nuclear energy
produced by the thermal pulse results in the expansion of the thin outer
layers on a timescale of 50 to a few 100 years, and the star evolves rapidly
redward on the H-R diagram, turning into a so-called born-again object (e.g.,
Langer et al. langer74 (1974); Iben iben76 (1976); Blöcker bloecker01 (2001)).
When the temperature of the expanding outer layers becomes low enough for
convection to penetrate to deep layers, the nuclear-processed material, such
as C, He, and s-process elements, is dredged up to the surface. Since the mass
of the H-rich outer layers is already small, this dredge-up process makes the
atmosphere H-deficient and C-rich. After the surface temperature reaches its
minimum of several 1000 K, the star evolves again blueward on a timescale of
$10^{3}$–$10^{4}$ years finally to a white dwarf.

The final flash episode is so brief that objects such as FG Sge are very rare.
Only three objects are known to date: FG Sge, Sakurai’s object (V4334 Sgr),
and V605 Aql. The evolution of FG Sge back to the AGB over $\sim$100 years is
consistent with the LTP scenario, while the return of Sakurai’s object and
V605 Aql to the AGB is even more rapid on the order of several years. They are
considered to have experienced the so-called very late thermal pulse (VLTP) –
a thermal pulse that occurs when the star is evolving downward to a white
dwarf on the H-R diagram (Clayton et al. clayton06 (2006)). In addition,
Reindl et al. (reindl17 (2017)) recently found that SAO244567, the central
star of the Stingray nebula, has been rapidly returning to the AGB since 2002,
and they conclude that it has undergone an LTP. The final flash episode is
considered to be a fairly common phenomenon, which occurs in 10–20% of stars
evolving off the AGB (Blöcker bloecker01 (2001); Herwig herwig05 (2005)).
Therefore, studies of final flash objects provide a unique opportunity to
obtain insight into the evolution of post-AGB stars.

Although FG Sge is discussed as an LTP object in the literature (e.g.,
Schönberner schoenberner08 (2008)), the evolution of its chemical composition
is not yet fully understood. The detailed spectroscopic analysis of Jeffery &
Schönberner (jeffery06 (2006)) shows that the atmosphere of FG Sge became
H-deficient some time between 1960 and 1995, and the detection of C2 (Iijima &
Strafella iijima93 (1993); Kipper & Kipper kipper93 (1993)) reveals that the
atmosphere is C-rich. However, this photospheric chemical composition of FG
Sge and its evolution on the H-R diagram cannot be consistently explained by
the present stellar evolutionary theory as discussed by Jeffery & Schönberner
(jeffery06 (2006)). The calculations of born-again giants (Blöcker 2001)
predict that the atmosphere becomes H-deficient only after the star begins to
evolve blueward again (i.e., when the temperature is increasing again).
However, FG Sge became H-deficient between 1960 and 1995, when the surface
temperature was still decreasing. In the case of a VLTP, hydrogen in the outer
layers becomes significantly depleted due to the hydrogen burning during the
helium flash (Herwig et al. herwig99 (1999)). This means that the star is
already H-deficient when it begins to evolve redward. However, FG Sge was
H-normal in 1960 when it was evolving redward. The alternative double-loop
scenarios proposed by Lawlor & MacDonald (lawlor03 (2003)) and Miller
Bertolami et al. (millerbertolami06 (2006)), in which FG Sge, Sakurai’s
object, and V605 Aql can be considered to represent different evolutionary
stages of the same evolutionary path, cannot entirely explain the observed
chemical composition of FG Sge either, as Jeffery & Schönberner (jeffery06
(2006)) have discussed.

Another remarkable change occurred in FG Sge in 1992, when its visual
brightness decreased suddenly by $\Delta V\approx 5$ magnitude and recovered
gradually over the course of $\sim$1 year. Since then, FG Sge has shown
sudden, unpredictable deep declines at an interval of 1–2 years; although,
such events have become so frequent that the central star has been seldom
visible since 2008 (Arkhipova et al. arkhipova22 (2022)). At the same time, it
started to show a noticeable infrared excess characterized by carbonaceous
dust with $\sim$1000 K (Gehrz et al. gehrz05 (2005)). The light curve of FG
Sge is similar to those of R CrB stars – a class of F–G supergiants in which
optically thick dust clouds are ejected in random directions – and an
occasional, deep decline is observed when a cloud happens to form in the line
of sight (Clayton clayton96 (1996)). While some authors propose that FG Sge is
undergoing such a random dust cloud formation (Gonzalez et al. gonzalez98
(1998)), Gehrz et al. (gehrz05 (2005)) have argued for a continuous outflow
with abrupt, episodic changes of the mass-loss rate. In either case, however,
it is unknown where and how dust can form around FG Sge, whose present
temperature (4500–5500 K) is deemed to be too high for dust to condensate.

To obtain a clearer picture of the dusty environment of FG Sge, high spatial
resolution mid-infrared observations are effective. In this paper, we present
the first spatially resolved mid-infrared (8–13 $\mu$m) observations of FG Sge
and radiative transfer modeling of its circumstellar envelope.

##  2 Observations

FG Sge was observed with the MIDI instrument (Leinert et al. leinert03 (2003))
at ESO’s Very Large Telescope Interferometer (VLTI) on 2008 June 22 (MJD54639)
using the unit telescope (UT) configuration UT2–UT3 with projected baseline
lengths of 43–46 m (Program ID: 081-D0224, P.I.: K. Ohnaka). A prism with a
spectral resolution of $\lambda/\Delta\lambda\simeq 30$ (at 10 $\mu$m) was
used to obtain spectrally dispersed fringes in the $N$ band between 8 and 13
$\mu$m. The data were taken in SCI-PHOT mode, in which photometric data of
each telescope are recorded simultaneously with interferometric data. A
detailed description of the observing procedure is given in Przygodda et al.
(przygodda03 (2003)), Leinert et al. (leinert04 (2004)), and Chesneau et al.
(chesneau05 (2005)). Table LABEL:table_obs summarizes the observations’ FG Sge
and the calibrator.

Table 1: MIDI observations of FG Sge and the calibrator: time of observation (coordinated universal time=UTC), projected baseline length $B_{\rm p}$, position angle of the projected baseline on the sky (P.A.), seeing in the visible, and the airmass at the times of the observations of FG Sge and the calibrator.  # | Object | $t_{\rm obs}$ | $B_{\rm p}$ | P.A. | Seeing | Airmass
---|---|---|---|---|---|---
|  | (UTC) | (m) | (°) | (″) |
2008 June 22 (MJD2454639) |
CAL1 | HD168723 | 07:27:42 | 46.0 | 46 | 1.0 | 1.35
SCI1 | FG Sge | 07:46:02 | 43.0 | 50 | 1.0 | 1.47
CAL2 | HD168723 | 08:36:15 | 43.3 | 43 | 1.0 | 1.84
SCI2 | FG Sge | 08:57:21 | 45.9 | 47 | 1.0 | 1.74
SCI3 | FG Sge | 09:06:47 | 46.1 | 46 | 1.1 | 1.82

We used the MIA+EWS package ver.2.0111Available at
http://www.strw.leidenuniv.nl/~nevec/MIDI to reduce the MIDI data (Leinert et
al. leinert04 (2004); Jaffe jaffe04 (2004)). We adopted a uniform-disk
diameter of 3.1 mas for the calibrator HD168723 (Bourges et al. bourges17
(2017)). It should be noted that there are significant differences in airmass
among the data sets. To avoid systematic errors, we calibrated the data sets
of FG Sge with those of the calibrator taken at a similar airmass: the data
set SCI1 was calibrated with CAL1, while the data sets SCI2 and SCI3 were
calibrated with CAL2.

We assessed the data quality and derived the calibrated visibilities as
described by Ohnaka et al. (ohnaka08 (2008)). Because it is difficult to
estimate the errors in the calibrated visibilities properly from two
calibrator data sets taken at very different airmasses, we assumed a total
relative error of 10% as adopted by Zhao-Geisler et al. (zhao-geisler11
(2011)) for SCI-PHOT mode data. Because the projected baseline length and
position angle of SCI2 and SCI3 are very close, we averaged the calibrated
visibilities of these two data sets. The visibilities obtained with MIA and
EWS agree well, and we only provide the results derived with EWS in the
discussion below.

The absolutely calibrated $N$-band spectra of FG Sge were obtained from the
same three SCI–CAL pairs taken at similar airmasses (SCI1–CAL1, SCI2–CAL2, and
SCI3–CAL2), using the method described in Ohnaka et al. (ohnaka07 (2007)). The
average and the standard deviation of the three spectra were taken as the
final spectrum and its error, respectively.

We also obtained the 19.5 $\mu$m flux of FG Sge with the mid-infrared imager
and spectrograph VLT/VISIR (Lagage et al. lagage04 (2004)). The VISIR
observations of FG Sge and a photometric standard star, HD173780, were carried
out on 2008 May 30, about three weeks before the MIDI observations, with the
Q3 filter (central wavelength = 19.5 $\mu$m, filter full width at half maximum
= 0.4 $\mu$m) and a pixel scale of 0$\aas@@fstack{\prime\prime}$075\. The data
with chopping and nodding were first reduced using the VISIR pipeline
ver.3.0.0 to obtain the images with the sky background almost subtracted. From
these images, we further subtracted the residual sky background measured in an
annular region with an inner and outer radius of
1$\aas@@fstack{\prime\prime}$5 and 2$\aas@@fstack{\prime\prime}$25,
respectively. The flux was obtained within an aperture radius of
1$\aas@@fstack{\prime\prime}$5 both for FG Sge and HD173780, and the data of
FG Sge were flux-calibrated with a Q3-band flux of 1.7276 Jy for
HD173780222https://www.eso.org/sci/facilities/paranal/instruments/visir/tools/
zerop_cohen_Jy.txt. The error in the calibrated flux of FG Sge was estimated
by the quadratic sum of the photon noise of the source and that of the
background. The observed 19.5 $\mu$m flux of FG Sge is $6.2\pm 0.49$ Jy.

Figure 1:  VLTI/MIDI observations of FG Sge. a: Observed visibilities of the
data set SCI1 (red dots) and the average of SCI2 and SCI3 (blue dots). The
error bars are only shown for the data set SCI1 for visual clarity. b: FWHM
obtained by a Gaussian fit to the observed visibilities.

Figure 2:  Observed mid-infrared spectrum of FG Sge. The solid red line with
the error bars shows the MIDI spectrum. Black dots correspond to the
Spitzer/IRS spectrum. The green diamond is the VISIR Q3-band (19.5 $\mu$m)
flux.

##  3 Results

Figure 1a shows the measured $N$-band visibilities, which slightly increase
from 0.8 at 8 $\mu$m to 0.9 at 13 $\mu$m. Figure 1b shows the full width at
half maximum (FWHM) obtained by the Gaussian fitting, which suggests that the
observed visibilities correspond to an approximately constant FWHM of
$\sim$10.5 mas. The observed visibilities suggest that the dust envelope of FG
Sge is compact compared to the other born-again objects’ Sakurai’s object and
V605 Aql. Chesneau et al. (chesneau09 (2009)) showed that the MIDI
visibilities of Sakurai’s object obtained at baseline lengths of 42–46 m – the
same as in our observations of FG Sge – range from 0.05 to 0.2, which is much
lower (thus much more extended) than the 0.8–0.9 obtained for FG Sge;
although, the distance to two objects is comparable (FG Sge: 3.0 kpc,
Sakurai’s object: 3.8 kpc, Evans et al. evans20 (2020)). Clayton et al.
(clayton13 (2013)) report that V605 Aql was so large that it was resolved out
in the MIDI observations with UT2 and UT3 (the same as used for FG Sge) in
spite of its distance of 4.6 kpc. This is primarily because the effective
temperature of FG Sge (4500–5500 K, see Sect. 4) is noticeably lower than
those of Sakurai’s object (12 000 K, van Hoof et al. vanhoof07 (2007)) and
V605 Aql (95 000 K, Clayton et al. clayton06 (2006)). In addition, the
luminosity of these latter objects ($\sim\\!\\!10^{4}$ $L_{\sun}$, Schönberner
schoenberner08 (2008); Chesneau et al. chesneau09 (2009)) is higher than that
of FG Sge ($4850^{+4600}_{-2690}$ $L_{\sun}$, Sect. 4), which can also make
them appear larger than FG Sge.

Figure 2 shows the MIDI spectrum (red solid line with the error bars) and the
19.5 $\mu$m flux measured with VISIR (green filled diamond). To complement the
MIDI $N$-band spectrum and the VISIR $Q$-band photometric data, we downloaded
the mid-infrared spectrum taken with the InfraRed Spectrograph333 The IRS was
a collaborative venture between Cornell University and Ball Aerospace
Corporation funded by NASA through the Jet Propulsion Laboratory and Ames
Research Center. (IRS, Houck et al. houck04 (2004)) onboard the Spitzer Space
Telescope (Werner et al. werner04 (2004))444Downloaded from
https://sha.ipac.caltech.edu/applications/
Spitzer/SHA/. The data were taken on 2007 December 11 (MJD 54446) in the short
and long low mode with spectral resolutions of 64–128, covering from 5 to 40
$\mu$m (Program ID: 40061, P.I.: A. Evans).

The 19.5 $\mu$m flux measured with VISIR is in good agreement with the Spitzer
data, which were taken approximately six months before the MIDI observations.
The MIDI spectrum is also consistent with the Spitzer/IRS spectrum longward of
$\sim$10 $\mu$m. However, the Spitzer/IRS spectrum shows a bump at $\sim$8
$\mu$m, which is not seen in the MIDI spectrum. Evans et al. (evans15 (2015))
reported this feature and the one at $\sim$6.2 $\mu$m as the unidentified
infrared emission, which is commonly attributed to polycyclic aromatic
hydrocarbons or PAH (e.g., Allamandola et al. allamandola99 (1999)). Evans et
al. (evans15 (2015)) show that the Spitzer spectrum taken on 2006 October 24
(MJD54032) also exhibit the PAH features at 6.2 and 8 $\mu$m, which suggests
that the PAH emission was stably present over a little more than a year from
October 2006 to December 2007. Therefore, it is not very likely – if not
entirely excluded – that the PAH emission disappeared entirely within six
months between the Spitzer observation in December 2007 and our MIDI
observations in June 2008. The absence of the PAH emission in the MIDI
spectrum can be due to a smaller field of view of 200 mas of MIDI compared to
the slit widths of
3$\aas@@fstack{\prime\prime}$6–10$\aas@@fstack{\prime\prime}$7 of Spitzer/IRS,
which suggests that the PAH emission may originate from an extended region.
Given that the central star is now H-deficient, the recently rejected inner
circumstellar matter should also be H-deficient, from which PAH emission is
not expected. Therefore, the PAH emission likely originates from the older,
outer circumstellar material ejected before the star became H-deficient. This
is consistent with the above observations. However, the UV radiation required
to give rise to the PAH emission cannot penetrate the optically thick dust
envelope suggested by our modeling presented below. We return to this point in
Sect. 4.

##  4 Modeling of the dust envelope

We carried out radiative transfer modeling with DUSTY (Ivezić et al. ivezic99
(1999)) for the observed $N$-band visibilities and the spectral energy
distribution (SED). The circumstellar envelope of FG Sge is considered to be
complex, given the R CrB-like light curve that suggests the formation of
clumpy dust clouds or episodic mass loss. However, the present MIDI data are
insufficient for constraining such complex structures. Therefore, we assumed
spherical shell models to constrain the global properties of the dusty
environment of FG Sge.

In the present work, we adopted a distance of $3.0^{+1.2}_{-1.0}$ kpc derived
by Chornay et al. (chornay21 (2021))555The errors in the distance correspond
to the 84th and 16th percentile of the distance posterior. based on the Gaia
Early Data Release 3 (Gaia Collaboration gaia21 (2021)). The bolometric flux
of the central star was obtained from the SED obtained before the onset of the
mass ejection in 1992. Jurcsik & Montesinos (jurcsik99 (1999)) presented the
photometric data from the $V$ to $N$ band obtained from June to August in
1983, which they fit with a photospheric model with $T_{\rm eff}$ = 5500 K and
$\log\varg$ = 1.5. We obtained a bolometric flux of $1.72\times 10^{-11}$ W
m-2 by integrating this SED, which translates into a luminosity of
$4850^{+4600}_{-2690}$ $L_{\sun}$ with the distance of $3.0^{+1.2}_{-1.0}$
kpc. It should be noted that the large uncertainty in the distance (hence in
the luminosity) does not affect the fit to the SED and visibility, because the
luminosity of the model star increases as $\propto d^{2}$, but the flux
predicted at the distance of the Earth decreases as $\propto d^{-2}$,
resulting in the same predicted flux. Likewise, the inner radius of the dust
shell increases as $\propto\sqrt{L}\propto d$, but its angular size remains
unaffected because it follows $\propto d^{-1}$.

To obtain the SED observed as contemporaneously as possible with our MIDI
observations, we used the visible ($V$, $R$, $I$, and $R_{c}$ bands)
photometric data taken on MJD54642 (just three days after our MIDI
observations) by Arkhipova et al. (arkhipova09 (2009)) and the infrared
($JHKLM$) photometric data taken on MJD54659.5 (20 days after the MIDI
observations) by Shenavrin et al. (shenavrin11 (2011)). As presented above,
the MIDI spectrum and VISIR 19.5 $\mu$m flux are in agreement with the Spitzer
spectrum obtained in December 2007 except for the PAH emission. Therefore, we
used these data to cover the mid-infrared region of the SED. We adopted
$E(g-r)=0.23$ based on the 3D extinction map of Green et al. (green19
(2019))666http://argonaut.skymaps.info for the distance of 3.0 kpc and
converted it to $E(B-V)=0.202$ using the relations given in Schlafly et al.
(schlafly11 (2011)). This translates into $A_{V}=0.63$ with $R_{v}=3.1$
assumed. The observed SED was dereddened by applying the wavelength dependence
of the interstellar extinction from Cardelli et al. (cardelli89 (1989)).

The effective temperature of the central star is controversial. Jeffery &
Schönberner (jeffery06 (2006)) concluded that FG Sge had been maintaining
$T_{\rm eff}$$\approx$5500 K since 1980 (as of 1995) and, therefore, proposed
that this $T_{\rm eff}$ might be the lowest limit in its redward evolution. On
the other hand, Fadeyev (fadeyev19 (2019)) derived $T_{\rm eff}$ = 4445 K
based on theoretical modeling of FG Sge’s pulsation. We treated $T_{\rm eff}$
as a free parameter and calculated models with $T_{\rm eff}$ = 4500, 5000, and
5500 K. Combined with the luminosity of 4850 $L_{\sun}$, these values of
$T_{\rm eff}$ correspond to a stellar radius of 115, 93, and 77 $R_{\sun}$,
respectively.

Given the C-rich nature of the photosphere and the circumstellar envelope, we
assumed amorphous carbon for the dust grains around FG Sge. We adopted a grain
size of 0.01 $\mu$m as concluded by Arkhipova et al. (arkhipova22 (2022)) to
explain the SED before 2019. The density distribution of dust grains was
assumed to be $\propto r^{-p}$, where $p$ is a free parameter, and
characterized by the dust condensation temperature ($T_{c}$), which determines
the inner radius of the envelope ($r_{\rm in}$), and the optical depth in the
radial direction at 0.55 $\mu$m ($\tau_{V}$). The outer radius ($r_{\rm out}$)
of the circumstellar envelope was estimated in the following manner. Gonzalez
et al. (gonzalez98 (1998)) measured an outflow velocity of $\sim$200 km s-1
from the Na I D lines. Because the mass ejection started in August 1992, the
material should have traveled $1.0\times 10^{16}$ cm (670 AU) at the time of
our MIDI observations in June 2008. On the other hand, DUSTY models show that
the inner radius of the envelope ranges from 25 to 40 $R_{\star}$ for $T_{\rm
eff}$ = 4500 to 5500 K for a condensation temperature of $\sim$1000 K. Using
the above values of $R_{\star}$, the inner radius is estimated to be
$(1-3)\times 10^{14}$ cm. Therefore, the maximum distance that the ejected
material traveled ($1.0\times 10^{16}$ cm) corresponds to 33–100 $r_{\rm in}$.
We adopted a fixed value of 70 $r_{\rm in}$.

We computed models with $T_{\rm eff}$ = 4500, 5000, and 5500 K, $p$ = 2.0 …
4.0 ($\Delta p$ = 0.5), $T_{c}$ = 800 … 1600 K ($\Delta T_{c}=200$ K), and
$\tau_{V}$ = 5 … 30 ($\Delta\mbox{$\tau_{V}$}=5$). However, we found no model
that can reasonably reproduce the observed SED and $N$-band visibilities
simultaneously. For example, the models with the density profile $\propto
r^{-2}$, which corresponds to stationary mass loss with a constant outflow
velocity, can fit the observed SED very well, but they systematically predict
the $N$-band visibility to be much lower ($\la$0.5) than the observed data,
because the dust envelope appears to be much more extended than the observed
data suggest. The models with steeper density profiles can explain the
observed $N$-band visibilities. However, these models predict the flux
longward of $\sim$10 $\mu$m to be too low compared to the observed data,
because there is too little cold material in the outer region of the envelope.

Figure 3:  Comparison of a spherical shell model with the observed SED and
$N$-band visibilities of FG Sge. The model is characterized with $T_{\rm in}$
= 1100 K, a piecewise power-law density profile with $r^{-3.5}$
($\mbox{$r_{\rm in}$}\leq r\leq\mbox{$r_{\rm tr}$}$ with $r_{\rm in}$ = 33
$R_{\star}$ and $r_{\rm tr}$ = 5 $r_{\rm in}$) and $r^{-2}$ ($r>\mbox{$r_{\rm
tr}$}$), $\tau_{V}$ = 8, and $T_{\rm eff}$ = 5000 K. a: Observed SED. Orange
filled triangles correspond to visible photometric data from Arkhipova et al.
(arkhipova09 (2009)). Orange filled diamonds are the infrared photometric data
from Shenavrin et al. (shenavrin11 (2011)). The photometric data were
dereddened with $A_{V}$ = 0.63. The thick solid red line is the Spitzer/IRS
spectrum. The thin solid blue line corresponds to the MIDI spectrum. The blue
cross corresponds to VISIR photometry. The DUSTY model is shown with the solid
black line. b: Observed visibilities are shown in the same manner as in Fig.
1a. The solid red and blue lines represent the model visibilities predicted
for the 43.0 and 46.0 m baselines, respectively.

It is possible that the density profile is steeper than the $r^{-2}$ in the
inner region of the dust envelope before the velocity reaches its terminal
velocity. We attempted to explained the observed data with models in which the
density profile falls as $r^{-p}$ from the inner radius $r_{\rm in}$ to a
certain radius $r_{\rm tr}$ and then changes to $r^{-2}$. Because the above
model grid shows that $\tau_{V}$ should be below $\sim$10 to explain the
visible part of the observed SED, we calculated models with the following
parameter ranges: $T_{\rm eff}$ (K) = 4500, 5000, and 5000; $\tau_{V}$ = 4.0 …
10.0 ($\Delta\mbox{$\tau_{V}$}\ =2.0$); $T_{\rm in}$ (K) = 800 … 1600
($\Delta\mbox{$T_{\rm in}$}$ = 100 K); $p$ = 2.5 … 4.5 ($\Delta p=0.5$); and
$r_{\rm tr}$/$r_{\rm in}$ = 2.5, 5, 10, 15, and 20. Figure 3 shows one of the
best-fit models, which reproduces the observed SED and $N$-band visibilities
fairly well. The model is characterized by $T_{\rm eff}$ = 5000 K, $\tau_{V}$
= 8, $T_{\rm in}$ = 1100 K ($r_{\rm in}$ = 33 $R_{\star}$), $p$ = 3.5, and
$r_{\rm tr}$ = 5 $r_{\rm in}$. The model still predicts the flux to be
somewhat lower than the observed data longward of 20 $\mu$m, and the model
visibilities are slightly lower than the MIDI data. Decreasing $p$ and/or
$r_{\rm tr}$ makes the fit to the SED better, but such models predict the
$N$-band visibilities to be even lower, because their density profiles become
similar to $r^{-2}$. Our model grid constrains the parameter ranges as
follows: $\tau_{V}$ = $8\pm 2$, $T_{\rm in}$ = $1100\pm 100$ K, $p=3.5\pm
0.5$, and $r_{\rm tr}$ = $10\pm 5$. We cannot constrain $T_{\rm eff}$ of the
central star within the adopted range between 4500 and 5500 K, because it is
heavily obscured.

Recently, Arkhipova et al. (arkhipova22 (2022)) modeled the SED observed from
0.4 to 5.0 $\mu$m in 2019 with a density profile of $\propto r^{-2}$ and
obtained an optical depth of 4.5 at 0.55 $\mu$m, which is significantly lower
than the $8\pm 2$ derived in the present work. This is because they modeled
the photometric data taken in 2019, when the dust clouds temporarily cleared.
The central star became visible with $V\approx 13$, which is $\sim$4 magnitude
fainter than the $V\approx 9$ before the start of the dust formation in 1992.
With $\mbox{$\tau_{V}$}=0.921\times\Delta V$ (if scattered light and dust
thermal emission are ignored), $\tau_{V}$ is estimated to be 3.7, which is in
broad agreement with 4.5 as derived from their modeling. As mentioned above,
the models with the density profile $r^{-2}$ can explain the observed SED, but
not the $N$-band visibilities. Given that Arkhipova et al. (arkhipova22
(2022)) only modeled the SED taken at a different epoch (May 2019) from ours,
we cannot draw the conclusion that the difference in the density profile
reflects the lack of spatially resolved data in the modeling or a temporal
variation in the density profile.

The remaining discrepancy between the best-fit models and the observed data is
likely due to the complex, clumpy nature of the circumstellar environment of
FG Sge as inferred from the irregular, R CrB-like declines in the light curve.
Furthermore, a disk-like or ring-like structure may also be present, as
observed in the other two born-again objects’ Sakurai’s object and V605 Aql
(Chesneau et al. chesneau09 (2009); Hinkle et al. hinkle14 (2014), hinkle20
(2020); Clayton et al. clayton13 (2013)). As described in Sect. 2, if the PAH
emission indeed originates from the outer region of the circumstellar
envelope, the UV radiation should escape the optically thick envelope. The
clumpy and/or disk-like structures can provide a natural explanation for the
UV radiation to escape and excite the carrier of the PAH emission even if it
is located farther away from the star.

##  5 Discussion and concluding remarks

The dust mass obtained from the model shown in Fig. 3 is $7.4\times 10^{-7}$
$M_{\sun}$, adopting a bulk density of 2.25 g cm-3 for amorphous carbon
(Gilman gilman74 (1974)). If we assume a gas-to-dust ratio of 200, the total
envelope mass is $1.5\times 10^{-4}$ $M_{\sun}$. Given that the mass ejection
started in 1992, the average mass-loss rate at the time of the MIDI
observations in 2008 is $9.3\times 10^{-6}$ $M_{\sun}$ yr-1. The dust mass of
FG Sge seems to be lower than $6\times 10^{-5}$ $M_{\sun}$ as found in
Sakurai’s object (Chesneau et al. chesneau09 (2009)) and $\sim\\!2\times
10^{-3}$ $M_{\sun}$ as derived for V605 Aql (Clayton et al. clayton13 (2013)).
The difference in the ejected mass may be related to the difference in the
nature between FG Sge and the other two final-flash objects – Sakurai’s object
and V605 Aql are VLTP objects, while FG Sge is likely to be an LTP object.
However, follow-up observations of FG Sge would be necessary to study the
long-term variation of the mass of the ejected dust.

As the light curve published by Arkhipova et al. (arkhipova22 (2022)) shows,
FG Sge is still experiencing the continuous ejection of dusty material. If we
adopt the outflow velocity of 200 km s-1 as inferred from the Na I D lines,
the material at the outer radius of the dust envelope in 2008 should have
traveled 560 au ($\sim$1200 $R_{\star}$) by now, which translates into an
angular displacement of 190 mas at the distance of 3.0 kpc. While the farthest
(thus the coldest) material at the outer boundary may not be detectable in the
10 $\mu$m region, temporal variations in the circumstellar envelope can be
expected. Furthermore, given the presence of disk-like structures in Sakurai’s
object and V605 Aql, it is of great interest to study whether FG Sge has a
disk or bipolar structures. To better constrain the structure and evolution of
the circumstellar environment of FG Sge, interferometric measurements at more
baselines and at multiple epochs are indispensable, which can be carried out
with the VLTI/MATISSE instrument (Lopez et al. lopez22 (2022)).

###### Acknowledgements.

We thank the VLTI team at ESO and the MIDI team for carrying out the
observations and making the data reduction software publicly available. K.O.
acknowledges the support of the Agencia Nacional de Investigación Científica y
Desarrollo (ANID) through the FONDECYT Regular grant 1210652. This research
made use of the SIMBAD database, operated at the CDS, Strasbourg, France.

## References

  * (1) Allamandola, L. J., Hudgins, D. M., Sandford, S. A. 1999, ApJ, 511, L115
  * (2) Arenou, F., Grenon, M., & Gómez, A. 1992, A&A, 258, 104
  * (3) Arkhipova, V. P., Esipov, V. F., Ikonnikova, N. P., Komissarova, G. V., & Shugarov, S. Yu. 2009, Astronomy Letters, 35, 534
  * (4) Arkhipova, V. P., Ikonnikova, N. P., Shenavrin, V. I., et al. 2022, Astronomy Letters, in press, https://arxiv.org/abs/2206.02493
  * (5) Blöcker, T. 2001, Ap&SS, 275, 1
  * (6) Bourges, L., Mella, G., Lafrasse, S., et al. 2017, JMMC Stellar Diameters Catalogue Version 2
  * (7) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * (8) Chesneau, O., Meilland, A., Rivinius, T., et al. 2005, A&A, 435, 275
  * (9) Chesneau, O., Clayton, G. C., Lykou, F., et al. 2009, A&A, 493, L17
  * (10) Clayton, G. C. 1996, PASP, 108, 225
  * (11) Clayton, G. C., Kerber, F., Pirzkal, N., et al. 2006, ApJ, 646, L69
  * (12) Clayton, G. C., Bond, H. E., Long, L. A., et al. 2013, ApJ, 771, 130
  * (13) Chornay, N., & Walton, N. A. 2021, A&A, 656, A110
  * (14) Evans, A., Gehrz, R. D., Helton, L. A., & Woodward, C. E. 2015, EAS Publications Series, Vol. 71–72, 281
  * (15) Evans, A., Gehrz, R. D., Woodward, C. E., et al. 2020, MNRAS, 493, 1277
  * (16) Fadeyev, Yu. A. 2019, Astronomy Letters, 45, 655
  * (17) Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1
  * (18) Gehrz, R. D., Woodward, C. E., Temim, T., Lyke, J. E., Mason, C. G. 2005, ApJ, 623, 1105
  * (19) Gilman, R. C. 1974, ApJS, 28, 397
  * (20) Gonzalez, G., Lambert, D. L., Wallerstein, G., et al. 1998, ApJS, 114, 133
  * (21) Green, G. M., Schlafly, E., Zucker, G., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * (22) Herwig, F. 2005, ARA&A, 43, 435
  * (23) Herwig, F., Blöcker, T., Langer, N., & Driebe, T. 1999, A&A, 349, L5
  * (24) Hinkle, K. H., & Joyce, R. R. 2014, ApJ, 785, 146
  * (25) Hinkle, K. H., Joyce. R. R., Matheson, T., Lacy, J. H., & Richter, M. J. 2020, ApJ, 904, 34
  * (26) Houck, J., Roellig, T., van Cleve, J., et al. 2004, ApJS, 154, 18
  * (27) Iben, I., Jr. 1976, ApJ, 208, 165
  * (28) Iijima, T., & Strafella, F. 1993, Inf. Bull. Variable Stars, 3959, 1
  * (29) Ivezić, Ž., Nenkova, M., Elitzur, M. 1999, User Manual for DUSTY, University of Kentucky Internal Report, accessible at http://www.pa.uky.edu/m̃oshe/dusty
  * (30) Jaffe, W. 2004, SPIE Proc., 5491, 715
  * (31) Jeffery, C. S., & Schönberner, D. 2006, A&A, 459, 885
  * (32) Jurcsik, J., & Montesinos, B. 1999, New Astronomy Reviews, 43, 415
  * (33) Kipper, T., & Kipper, M. 1993, A&A, 276, 389
  * (34) Lagage, P. O., Pel, J. W., Authier, M., et al. 2004, The Messenger, 117, 12
  * (35) Langer, G. E., Kraft, R. P., & Anderson, K. S. 1974, ApJ, 189, 509
  * (36) Lawlor, T. M., & MacDonald, J. 2003, ApJ, 583, 913
  * (37) Leinert, Ch., Graser, U., Przygodda, F. 2003, Ap&SS 286, 73
  * (38) Leinert, Ch., van Boekel, R., Waters, L. B. F. M., et al. 2004, A&A, 423, 537
  * (39) Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192
  * (40) Miller Bertolami, M. M., Althaus, L. G., Serenelli, A. M., & Panei, J. A. 2006, A&A, 449, 313
  * (41) Ohnaka, K., Driebe, T., Weigelt, G., & Wittkowski, M. 2007, A&A, 466, 1099
  * (42) Ohnaka, K., Driebe, T., Hofmann, K.-H., Weigelt, G., & Wittkowski, M. 2008, A&A, 484, 371
  * (43) Przygodda, F., Chesneau, O., Graser, U., Leinert, Ch., & Morel, S. 2003, Ap&SS, 286, 85
  * (44) Reindl, N., Rauch, T., Miller Bertolami, M. M., Todt, H., & Werner, K. 2017, MNRAS, 464, L51
  * (45) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * (46) Schönberner, D. 2008, ASP Conf. Ser. 391, 139
  * (47) Shenavrin, V. I., Taranova, O. G., & Nadzhip, A. E. 2011, Astronomy Reports, 55, 31
  * (48) van Hoof, P. A. M., Hajduk, M., Zijlstra, A. A., et al. 2007, A&A, 471, L9
  * (49) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1
  * (50) Zhao-Geisler, R., Quirrenbach, A., Köhler, R., Lopez, B., & Leinert, C. 2011, A&A, 530, A120