S

Yang W. Phillips Song J. Harder Cha<EMAIL_ADDRESS>Email<EMAIL_ADDRESS>
[a]Department of Engineering Science, University of Oxford, Parks Road,
Oxford, OX1 3PJ, UK [b]Advanced Photon Source, Argonne National Laboratory,
Argonne, IL 60439, USA

# Annealing of focused ion beam damage in gold microcrystals: An in situ Bragg
coherent X-ray diffraction imaging study

David  Nicholas  Kay  Ross  Wonsuk

###### Abstract

Focused ion beam (FIB) techniques are commonly used to machine, analyse and
image materials at the micro- and nanoscale. However, FIB modifies the
integrity of the sample by creating defects that cause lattice distortions.
Methods have been developed to reduce FIB-induced strain, however these
protocols need to be evaluated for their effectiveness. Here we use non-
destructive Bragg coherent X-ray diffraction imaging to study the in situ
annealing of FIB-milled gold microcrystals. We simultaneously measure two non-
collinear reflections for two different crystals during a single annealing
cycle, demonstrating the ability to reliably track the location of multiple
Bragg peaks during thermal annealing. The thermal lattice expansion of each
crystal is used to calculate the local temperature. This is compared to
thermocouple readings, which are shown to be substantially affected by thermal
resistance. To evaluate the annealing process, we analyse each reflection by
considering facet area evolution, cross-correlation maps of displacement field
and binarised morphology, and average strain plots. The crystal’s strain and
morphology evolve with increasing temperature, which is likely to be caused by
the diffusion of gallium in gold below $\sim$280°C and the self-diffusion of
gold above $\sim$280°C. The majority of FIB-induced strains are removed by
380-410°C, depending on which reflection is being considered. Our observations
highlight the importance of measuring multiple reflections to unambiguously
interpret material behaviour.

##  1 Introduction

Focused ion beam (FIB) techniques utilise a beam of energetic ions, often
gallium (Ga) extracted from liquid metal, to remove material at the micro- and
nanoscale [Rajput2015]. The ability to control sample features at these
previously inaccessible length scales has furthered the understanding of
material properties, e.g., by enabling the exploration of size-dependent
material behaviour. Most notably, the “smaller-is-stronger” effect has been
extensively studied using nanoindented FIB-machined specimens [Uchic2004,
Greer2011]. The deformation behaviour of these samples varied significantly,
depending on several competing mechanisms. One mechanism, exhibited by
polycrystalline copper films [Kiener2007] and Mo alloy single crystals
[Bei2007], shows that FIB damage can lead to a dense population of defects,
making the escape of dislocations at the surface more difficult, which causes
an apparent hardening. On the contrary, softening has been reported in FIB-
fabricated Mo alloy micropillars with a low initial defect density [Shim2009].
Here FIB-induced dislocation loops act as sources for glide dislocations,
leading to a reduction in yield stress [Shim2009]. Additionally, FIB can cause
other material changes such as amorphisation [Basnar2003], local
recrystallisation [Rubanov2004], twin nucleation [Lou2018], and intermetallic
phase formation [Babu2016], to name a few.

Experiments on gold (Au) microcrystals have shown that FIB milling can cause
lattice strains that extend more than 100 nm below the surface of the material
[Hofmann2017a, Hofmann2018a]. This can be a major concern as samples are no
longer in their original state and these FIB-induced loads will act in
addition to externally applied stresses. Notwithstanding, FIB remains a
regularly used tool for 3D material characterisation [Nan2019],
micromechanical test sample fabrication [Brousseau2010], and top-down
preparation of transmission electron microscopy (TEM) samples [Rubanov2004,
Rajput2015]. The severity of FIB-induced sample modification depends on
specific milling conditions [Hutsch2014]. Methods introduced to mitigate this
damage must be versatile and effective to allow FIB milling to reach its full
potential as a micromachining tool.

Numerous approaches have been suggested for removing FIB damage: use of a
plasma FIB [Ernst2017], post-milling low energy (2 keV) ion polishing
[Giannuzzi2005], sacrificial coatings [Rubanov2001], chemical etching
[Roediger2011], and post-manufacture heat treatment [Kiener2012, Lee2016].
Here we investigate the effects of annealing, a process that increases defect
mobility and thereby provides a pathway for the removal of FIB-induced defects
and, more importantly, the resulting strains. To evaluate the effectiveness of
annealing for FIB-induced strain removal, we require a technique that allows
nanoscale spatial resolution and high strain resolution. Previous in situ TEM
studies have been performed, showing that temperatures of half the material’s
melting point were sufficient to remove FIB-induced dislocation loops in Al
micropillars [Lee2016] and reduce defect density in FIB-made Cu samples
[Kiener2012]. However, TEM is limited to thin samples, where the close
proximity of free surfaces can dominate behaviour. Hence, we require a
technique capable of capturing high strain resolution information in 3D
micron-sized volumes.

Bragg coherent X-ray diffraction imaging (BCDI) allows for 3D-resolved,
nanoscale strain measurements and can provide 3D spatial resolution of less
than ten nanometers [Cherukara2018] as well as strain resolution on the order
of $\sim$2$\times 10^{-4}$ [Carnis2019, Hofmann2020]. BCDI involves fully
illuminating a crystalline sample inside the coherent volume of an X-ray beam,
which is approximately $1\mathrm{\mu m}\times 1\mathrm{\mu m}\times
1\mathrm{\mu m}$ at third generation synchrotron sources [Clark2012,
Hofmann2017b]. Once the Bragg condition is met for a specific $hkl$
reflection, the diffraction pattern is collected on a pixelated area detector
positioned perpendicular to the outgoing wave vector in the Fraunhofer
diffraction regime. By rotating the sample about a rocking axis, a 3D coherent
X-ray diffraction pattern (CXDP) is collected as the detector sequentially
intersects different parts of the chosen 3D Bragg peak in reciprocal space
[Robinson2001]. If the CXDP is oversampled by at least twice the Nyquist
frequency [Miao2000b], at least 4 pixels per fringe period, iterative phase
retrieval algorithms that apply constraints in real and reciprocal space can
be used to recover the phase [Fienup1982, Robinson2001]. The amplitude and
phase in reciprocal space are related to the real-space object via an inverse
Fourier transform [Miao2000b, Robinson2001, Clark2012] followed by a space
transformation from detector conjugated space to orthogonal lab or sample
space [Yang2019, Li2019b, Maddali2020]. The resulting amplitude,
$\mathbf{\rho(r)}$, where $\mathbf{r}$ is the position vector, is proportional
to the effective electron density of the crystalline volume associated with
the particular crystal reflection. The phase, $\mathbf{\psi(r)}$, corresponds
to the projection of the lattice displacement field, $\mathbf{u(r)}$, onto the
Bragg vector, $\mathbf{Q_{\mathit{hkl}}}$, of the $hkl$ crystal reflection
under consideration:

$\mathbf{\psi_{\mathit{hkl}}(r)}=\mathbf{Q_{\mathit{hkl}}}\cdot\mathbf{u(r)}$
(1)

Previously, BCDI has been used to investigate annealing of 3C-SiC
nanoparticles [Hruszkewycz2018], diamond nanoparticles [Hruszkewycz2017b],
magnetite crystals [Yuan2019] and metallic glass [Chen2020]. Furthermore, BCDI
has been extensively used to study Au subject to iron diffusion
[Estandarte2018], copper diffusion [Xiong2014a] and catalytic oxidation
[Suzana2019]. However, these experiments only considered one reflection, and
thus only the strain along the Bragg vector could be resolved. Our previous
BCDI studies found that FIB-milled Au crystals had large phase features
corresponding to lattice dilation and contraction [Hofmann2017b, Hofmann2017a,
Hofmann2018a]. These features were concentrated at the top and on the side of
the milled crystals–on the surfaces exposed to the FIB beam. Defects such as
stair-rod dislocations and their corresponding Burgers vectors were identified
in the FIB-affected sample volume [Hofmann2017a, Hofmann2018a]. This contrasts
the comparatively flat, homogeneous phase observed in the as-grown crystals
[Harder2013, Hofmann2017a, Hofmann2018a]. The study [Hofmann2018a] further
characterised the removal of substantial FIB damage by using a 5 keV polishing
step. The ability to capture these observations throughout an entire 3D volume
with high strain sensitivity and spatial resolution makes BCDI attractive for
the study of FIB-damaged materials.

Here we investigate the in situ annealing of FIB-milled Au microcrystals using
BCDI to examine how the 3D strain and crystal morphology evolve as a function
of temperature. We demonstrate the successful measurement of two unique
reflections for two separate FIB-milled Au crystals during an in situ
annealing cycle from 23°C (room/ambient temperature) to 587°C. After
recovering the amplitude and phase of the crystals, we identify temperatures
where structural changes occur through observations of crystal facet areas,
Pearson’s normalised cross-correlation coefficient, $r$, of
$\mathbf{u_{\textit{hkl}}(r)}$ and binary morphology masks, and average strain
in different crystal regions of interest (ROI), which are defined according to
their FIB exposure. The use of two reflections for each sample highlights
subtle differences in the lattice relaxation, which is observed to be
dependent on the crystallographic direction and the geometry of the sample. We
also consider the diffusion of Ga in Au and the self-diffusion of Au to
determine which mechanism is probably responsible for specific changes. With
four reflections, we faithfully identify temperatures where transitions in
crystal shape, $r$, and average strain occur to evaluate the effectiveness of
annealing in eliminating FIB-induced material damage.

##  2 Experimental methodology

###  2.1 Sample preparation

Au microcrystals were prepared on a [001]-oriented silicon substrate with a
100 nm thick thermally grown oxide layer. First, the wafer was spin-coated
with ZEP resist and patterns of 2 $\mu\mathrm{m}$ wide lines were produced on
the substrate using electron-beam lithography. Following the removal of the
resist, the wafer was coated with a 3 nm thick Cr adhesion layer followed by
40 nm of Au. A pattern of Au lines remained on the substrate after a “lift-
off” procedure was used to remove the unpatterned areas. Next, the sample was
annealed at 1000°C in air for 10 h, dewetting the lines to form microcrystals
ranging from 200 nm to 1 $\mathrm{\mu m}$ in diameter. Scanning electron
microscopy (SEM) was used to select suitable crystals for this study.

Using a Zeiss AURIGA FIB-SEM, candidate crystals were subjected to FIB milling
with Ga ion energy of 30 keV and beam current of 20 pA, as this closely mimics
the final milling cut during the manufacture of micromechanics test specimens
[Kiener2008, Li2015, Gong2015]. Approximately one third of the crystal was
removed using an incremental trench-milling mode. SEM imaging rather than low
dose FIB imaging was used to position the FIB milling scan, as our previous
results indicate that even a single FIB image can induce large lattice strains
in the Au microcrystals [Hofmann2017a]. Fig. 1(a) shows high-resolution SEM
micrographs of the chosen FIB-milled microcrystals, crystal A and B, in
excellent agreement with their respective morphologies recovered from BCDI at
room temperature in Fig. 1(b-c). Crystal A and B were $\sim$10 $\mathrm{\mu
m}$ apart on the substrate.

Figure 1: Comparison between SEM micrographs of crystals A and B and sample
morphology recovered from BCDI at room temperature (23°C). (a) SEM micrographs
of the FIB milled crystals. (b) Superimposed reconstructions of the two
measured linearly independent reflections for each crystal. Opaque
morphologies are rendered to allow examination of their agreement. The near
vertical line in crystal B corresponds to a twin domain boundary on the bottom
face. (c) Average morphology of the two reflections for each crystal. The
amplitude threshold for (b) and (c) is 0.30. The coordinate axes are
positioned in sample space and plotted with a length of 500 nm.

###  2.2 Bragg coherent X-ray diffraction

Synchrotron X-ray diffraction measurements were carried out at the Advanced
Photon Source (APS), Argonne National Laboratory, USA. Before BCDI
measurements, micro-beam Laue diffraction at beamline 34-ID-E was used to
determine the lattice orientation of the crystals on the silicon wafer. A
procedure for the alignment of BCDI measurements, based on Laue diffraction
orientation measurements, is provided elsewhere [Hofmann2017b].

BCDI measurements at beamline 34-ID-C used an X-ray beam with 9 keV photon
energy $(\lambda=0.138\mathrm{\ nm})$ and a spot size of $760\mathrm{\
nm}\times 750\mathrm{\ nm}$ (full width at half maximum) using
Kirkpatrick–Baez mirrors [Kirkpatrick1948]. The Si substrate was clipped onto
a heater cell available at the beamline shown in Appendix A. A CXDP was
collected on a $256\times 256$ pixel quadrant of a Timepix detector (Amsterdam
Scientific Instruments) with a GaAs sensor and pixel size, $p$, of
$55\mathrm{\ \mu m}\times 55\mathrm{\ \mu m}$ positioned at a distance, $D$,
of $1.2$ m away from the sample to ensure that the diffraction pattern was
oversampled [Sayre1952]. The lower bound for $D$ was determined by
$D=\frac{2dp}{\lambda}$, where $d$ is the sample size. The $(002)$ and
$(1\bar{1}1)$ reflections were measured for crystal A and the $(200)$ and
$(11\bar{1})$ reflections for crystal B at every temperature step.

At ambient temperature, the $(1\bar{1}1)$ reflection for crystal A was rotated
through an angular range of $0.7\degree$ while the remaining reflections
occupied a range of $0.6\degree$. For each rocking curve, a slice through the
CXDP was captured every $0.004\degree$ with 0.2 s exposure time and 20
accumulations. A total of 175 increments were collected for the $(1\bar{1}1)$
reflection and 150 increments for the other reflections. Each CXDP measurement
was repeated multiple times to optimise the signal-to-noise ratio. For crystal
A, 8 scans were collected for the $(002)$ reflection and 14 for the
$(1\bar{1}1)$ reflection, while for crystal B, 9 scans were collected for the
$(200)$ reflection and 10 for the $(11\bar{1})$ reflection at room
temperature.

The angular range, number of accumulations, number of repeated scans and
exposure time were adjusted at higher temperatures to sufficiently sample the
temperature range as discussed in Section 2.2.1. From 90°C onward, an angular
range of $0.7\degree$ with $0.004\degree$ steps was used, i.e. 175 increments
were collected in the rocking scan for CXDPs above room temperature, and 10
accumulations were taken instead of 20. Only 3 repeated scans were collected
between 90°C and 209°C inclusive. From 213°C onward, a 0.1 s exposure time was
used and only 2 scans were measured, with the exception of the $(1\bar{1}1)$
reflection at 213°C where 3 scans were measured. After performing flat-field
and dead time corrections, repeated scans were aligned using a 3D version of
the approach proposed by [Guizar-Sicairos2008] to maximise their cross-
correlation and only scans with $r$ greater than 0.975 were averaged to return
the CXDP used for phasing. All scans were summed and used for phase retrieval.

####  2.2.1 In situ annealing

The temperature was increased using the step and hold temperature profile
shown in Fig. 2. Here we note that the calibrated temperature of the Au
microcrystals, $\mathrm{T_{calibrated}}$, is different from the temperature
recorded using the thermocouple, $\mathrm{T_{measured}}$. This discrepancy
arises from the thermal resistance between the thermocouple (attached to the
sample clip), and the microcrystals in the heating setup, shown in Appendix A.
To calculate $\mathrm{T_{calibrated}}$ from $\mathrm{T_{measured}}$, we use
the samples as local measurement devices. We consider the thermal expansion of
pure Au, $\epsilon_{\mathrm{thermal}}=\alpha_{\mathrm{Au}}\Delta T$, where
$\Delta T$ is the temperature change and $\alpha_{\mathrm{Au}}$ is the linear
thermal expansion coefficient [Howatson2009], to define
$\mathrm{T_{calibrated}}$. We assume that implanted Ga ions after FIB milling
have negligible effect on the Au thermal expansion coefficient, as the ions
only penetrate 20 nm below the surface [Hofmann2018a].

Figure 2: The annealing temperature profile as a function of time. Changes in
the temperature step size are marked by the red dashed lines. (a)
$\mathrm{T_{calibrated}}$ was raised from room temperature to 209°C using
$\sim$50°C steps with a ramp rate of 2°C/min. (b) From 209°C to 409°C,
$\sim$5°C steps were used with a ramp rate of 50°C/min. Finally, from 409°C to
587°C, $\sim$20°C steps were used with a ramp rate of 50°C/min. The ramp rate
reported here is with respect to $\mathrm{T_{measured}}$ and was lower in (a)
to ensure that the alignment macro could update the positions of the
reflections for larger temperature steps.

Since Au has fcc crystal structure, lattice thermal expansion acts uniformly
in all directions, causing the Bragg angle to decrease with increasing
temperature. To track the Bragg peak, an alignment macro was run during the
temperature ramp that ensured the diffractometer was positioned such that the
maximum intensity was always centred in the top left quadrant of the detector
immediately before the start of the rocking scans. This protocol ensured that
the strain measured by BCDI was relative to the zero-strain reference region
for each reflection as explained in Section 2.3, allowing a comparison between
temperatures. By tracking the change in detector position at every temperature
step, $i$, the lattice parameter, $a_{i}$, at each step can be calculated:

$a_{i}=\frac{\lambda}{2\sin(\theta_{i})}\sqrt{h^{2}+k^{2}+l^{2}}\mathrm{,\hskip
5.69054ptwhere\hskip
5.69054pt}\theta_{i}=\frac{1}{2}\cos^{-1}\left(\frac{\mathbf{S}_{i}\cdot\mathbf{S}_{\mathrm{0}}}{\mathbf{S}_{\mathrm{0}}\cdot\mathbf{S}_{\mathrm{0}}}\right)\\\
$ (2)

where $\theta_{i}$ is the Bragg angle, $\mathbf{S}_{i}$ is the outgoing
wavevector (the detector position) and $\mathbf{S}_{\mathrm{0}}$ is the
incident wavevector. Using $a_{i}$ and the literature value of
$\alpha_{\mathrm{Au}}$, we calculate $\mathrm{T_{calibrated}}$:

$T_{\mathrm{calibrated},i}=\frac{\epsilon_{\mathrm{thermal},i}}{\alpha_{\mathrm{Au}}}+23\degree\mathrm{C,\hskip
5.69054ptwhere\hskip
5.69054pt}\epsilon_{\mathrm{thermal},i}=\frac{a_{i}-a_{23\degree\mathrm{C}}}{a_{23\degree\mathrm{C}}}$
(3)

Fig. 3 shows the discrepancy between $\epsilon_{\mathrm{thermal}}$ plotted
against $\mathrm{T_{calibrated}}$ and $\mathrm{T_{measured}}$. The correction
for $\mathrm{T_{calibrated}}$ was determined by averaging over the four
reflections, with a mean standard deviation of $\sim$8°C, indicating good
agreement amongst the four reflections.

Figure 3: A comparison of $\epsilon_{\mathrm{thermal}}$ plotted against (a)
$\mathrm{T_{measured}}$ and (b) $\mathrm{T_{calibrated}}$ for each reflection.
The $\epsilon_{\mathrm{thermal}}$ was calculated using Eq. 2 and 3 by tracing
the angular position of the detector during annealing. In (a), the measured
$\epsilon_{\mathrm{thermal}}$ averaged over all reflections (red line) does
not agree with the predicted $\epsilon_{\mathrm{thermal}}$ (blue line) when
plotted using $\mathrm{T_{measured}}$. In (b), the measured
$\epsilon_{\mathrm{thermal}}$ plotted using $\mathrm{T_{calibrated}}$
(determined by $\alpha_{\mathrm{Au}}$) overlaps with the the predicted
$\epsilon_{\mathrm{thermal}}$ as intended.

In this paper all temperatures noted refer to $\mathrm{T_{calibrated}}$ rather
than $\mathrm{T_{measured}}$ unless specified. The temperature profile (shown
in Fig. 2) has three varying temperature steps so as to capture the entirety
of the annealing process whilst retaining a high sampling rate between
209-409°C. At low temperatures the migration energy was believed to be too
high for noticeable interstitial and vacancy movement, hence larger
temperature steps of $\sim$50°C were used. Below $\mathrm{T_{H\ddot{u}ttig}}$,
the Hüttig temperature, which marks the threshold for surface atom mobility
(127°C for Au) [Odarchenko2018], no changes are expected. Indeed in this
experiment, we focus on atomic mobility over tens of nanometers below the
surface, as this is the thickness of the layer in which in which Ga
implantation causes atomic displacement damage [Hofmann2018a]. Such solid-
state diffusion is expected to occur above $\mathrm{T_{Tammann}}$, the Tammann
temperature (327°C for Au) [Odarchenko2018], thus from 209°C $\sim$5°C steps
were taken to increase the temperature resolution of our data set in order to
capture structural transitions. It is important to note that the Hüttig and
Tammann temperatures are semi-empirical
($\mathrm{T_{H\ddot{u}ttig}=0.3T_{melting}}$) and
($\mathrm{T_{Tammann}=0.5T_{melting}}$) and depend on surface texture,
orientation, morphology and particle size [Moulijn2001] based on the Gibbs-
Thomson equation [Thomson1871]. After reaching 409°C, the temperature was
increased using $\sim$20°C steps. This was because much of the FIB damage was
expected to already have been removed, based on a study showing the near
complete recovery of nanoindented Au microcrystals annealed at $\sim$0.65
times the melting temperature (596°C for Au) [Kovalenko2017].

###  2.3 Phasing of the diffraction data

Before the real-space, complex-valued electron density of the crystals can be
reconstructed, the phase of the diffracted wavefield must be recovered using
an iterative phase retrieval algorithm [Robinson2001]. The use of established
phase revival algorithms enables the recovery of the object by iteratively
applying real and reciprocal space constraints to an initial guess
[Robinson2001, Williams2003, Marchesini2003, Clark2012]. Here, the starting
guess for each reflection is of particular interest, as seeding with a good
start guess generally results in a faster convergence. At each temperature
(except the initial and final temperatures), we could seed the reconstruction
with the reconstruction from the previous or succeeding temperature step (see
Appendix B). We instead seed each reconstruction with itself through multiple
rounds [Hofmann2020], with the initial round being seeded with a random guess
(see details in Appendix C).

Once the phases were retrieved, the strain fields, $\epsilon\mathbf{(r)}$,
projected along the Bragg vector were calculated using:

$\epsilon_{hkl}\mathbf{(r)}=\nabla\mathbf{\psi_{\mathit{hkl}}(r)}\cdot\frac{\mathbf{Q_{\mathit{hkl}}}}{|\mathbf{Q_{\mathit{hkl}}}|^{2}}$
(4)

These strain fields are relative to the zero-strain reference of the
reconstruction. The assumption of close to zero lattice strain beyond a
distance of $\sim$150 nm from the FIB-milled surface has previously been
validated [Hofmann2017a, Hofmann2018a]. Thus, we establish a $40\times
40\times 40$ voxel volume centred on the bulk ROI’s (Fig. 9) centre of mass, a
region relatively unaffected by FIB-induced strain, as the strain reference
for every reconstruction [Phillips2020]. The amplitude threshold is 0.30 for
all reconstructions. Furthermore, each final object is eroded by 1 voxel in an
attempt to remove a phasing surface artefact caused by steep phase gradients
across the air-sample interface convolved with the resolution function
[Carnis2019].

We report the 3D spatial resolution to be $\sim$23 nm. This was determined by
differentiating two amplitude line profiles drawn across the air-sample
interface and fitting each of them with a Gaussian. The final spatial
resolution was the average of the full width at half maximum of the two
Gaussians for all reconstructions. We note that the resolution is somewhat
dependent on the measurement directions in reciprocal space and the geometry
of the sample. Appendix D shows a plot of the spatial resolution of each line
profile for each crystal reflection with respect to temperature.

##  3 Results and Discussion

###  3.1 Crystal strain and morphology

Fig. 4 shows a comparison of the strain found along each reflection,
$\mathbf{Q_{\mathit{hkl}}}$, of the FIB-milled crystals between the initial
and final temperature.At room temperature, the reconstructed crystal
morphologies are rough, especially on the surface of the milled face where
complex sputtering and redeposition processes compete, which results in a
pitted surface morphology. In Fig. 4 the four trench-like features visible on
the surface of the $(200)$ reflection for crystal B are near-surface
dislocations caused by the milling process. These manifest as pipes of missing
intensity in the amplitude of the reconstruction and are at the centre of a
phase vortex [Clark2015, Ulvestad2015a, Hofmann2017a, Dupraz2017,
Hofmann2020]. The phase vortices of the $(200)$ reflection are identified in
Appendix E. Dislocations are visible only when
$\mathbf{Q_{\mathit{hkl}}}\mathbf{\cdot b}\neq 0$ [Williams2009] and therefore
are not necessarily visible for all reflections. For crystal B, the
dislocations are visible only for the $(200)$ data. From two reflections alone
it is not possible to uniquely identify the Burgers vector of these
dislocations, although it is likely that it is of $(a/3)<110>$ type
[Hofmann2017a].

Figure 4: Rendering of the strain in crystal A and B for each reflection at
23°C and 587°C. In the top row of each quadrant the milled face is shown with
a slice perpendicular to the x-axis. The bottom row of each quadrant is the 2D
view of the corresponding slice. The small edge at the bottom of crystal B is
a twin domain boundary marking a region of missing intensity. The black arrows
represent the Bragg vector direction. The coordinate axes are positioned in
sample space and plotted with a length of 250 nm. Numbered white arrows
highlight trench-like structures that are reminiscent of dislocations. The
blue hashed region at 587°C on the top face of the $(200)$ reflection for
crystal B is a reconstruction artefact (see Section 2.3) and does not
physically represent a compressive layer.

The room temperature 2D strain slices for each reflection are comparable in
magnitude, with a maximum strain of $\sim$1.5$\times 10^{-3}$, to published
data for glancing incidence FIB damage in Au [Hofmann2018a]. The same
publication reports a thin compressive region followed by a tensile region
along the milled face, which is prominently observed here in the $(200)$
reflection for crystal B. A previous normal incidence FIB-milling study on Au
microcrystals at room temperature [Hofmann2017a] reported a compressive layer
at the top of the samples, caused by normal incident angle ions impacting the
top of the crystal from the tails of the FIB probe. This is seen in the
$(002)$ reflection for crystal A and both the $(200)$ and $(11\bar{1})$
reflections for crystal B.

After heating to 587°C, the strain and the crystal morphology have changed
dramatically (Fig. 4) and presents a relatively constant strain in the 2D
slices. Additionally, the surface morphology is observed to become smooth and
evolution of the facets is observed. Supplementary videos (SV) SV1 and SV2
show the evolution of the surface morphology and internal strain for each
crystal. In Fig. 4 the dislocations or pits on the milled face in the $(200)$
reflection disappear after roughly (1) 296°C, (2) 352°C, (3) 372°C and (4)
380°C, indicated by the elimination of phase vortices, as surface
restructuring is well underway. These temperatures are mostly consistent with
being above the Tammann temperature (327°C) for Au.

###  3.2 Facet area

With increasing temperature, both milled faces expanded outwards along the
z-direction and became rounder at the bottom, demonstrating the formation of
low-index facets, shown in Fig. 5, as a means to reduce the overall free
surface energy. Facet areas were computed by calculating the area of the
crystal whose isonormal is within a $19\degree$ angular threshold of a given
facet normal unit vector. In principle a maximum angular threshold of
$27.4\degree$ could be used, which is half the angle between {111} and {100}
planes. Our chosen tolerance is below this maximum value since not all facets
are fully developed, and the edges of the crystal are rounded and thus should
not be classified as part of a facet.

Figure 5: Normal vectors for each separately coloured facet are shown at
$23\mathrm{\degree C}$ and $587\mathrm{\degree C}$ for both crystals. Note
only vectors pointing out of the page in the top and bottom views are drawn
for clarity. The coordinate axes are positioned in sample space and plotted
with a length of 250 nm. Figure 6: Low-index facet areas of crystal A and B as
a function of temperature. See Fig. 5 for facet locations on the crystal.

Fig. 6 and supplementary video SV3 show facet area evolution as temperature
increases for both crystals. Fig. 4 shows that the facet at the bottom right
edge of crystal A that appears after 410°C (see supplementary video SV1) is
$(00\bar{1})$ and the facet at the bottom left edge of crystal B that appears
after 441°C (see supplementary video SV2) is $(1\bar{1}\bar{1})$. The
identification of these facets is based on Fig. 5 and supplementary video SV3.
The formation of low-index facets has been suggested to relieve strain
[Kim2018], and contribute to thermodynamic equilibrium [Kovalenko2015].
However, we do not fully reach the Winterbottom thermodynamic equilibrium
shape [Winterbottom1967] due to finite measurement time. For Au, the order of
free surface energies for low index facets is $\\{111\\}<\\{100\\}<\\{110\\}$
based on ab initio techniques that concur with the number of nearest-neighbour
broken bonds [Galanakis2002]. Thermodynamically, the formation of a
$(1\bar{1}\bar{1})$ facet is most favourable, shown in crystal B. However, the
production of the $(00\bar{1})$ facet in crystal A and the disappearance of
the $(00\bar{1})$ facet in crystal B after 527°C are not improbable as normal
movement of facets in dislocation‐free single crystals can be restricted by
the magnitude of the nucleation barriers if the crystal is larger than a few
nanometers [Mullins2000]. In crystal B, the $(00\bar{1})$ facet shrinks after
527°C but the $(1\bar{1}\bar{1})$ facet grows, suggesting that the latter is a
lower energy configuration.

###  3.3 Normalised displacement field cross-correlation matrix

To pinpoint transitions in strain, we consider a Pearson $r$ correlation
matrix of the displacement fields between different temperatures for each
reflection. When analysing single reflections, this method has been used
[Yau2017a, Ulvestad2015c] to infer structural transitions. The value of $r$
between each combination of temperatures was computed using Eq. 5
[Ulvestad2015c]:

$r(x,y)=\frac{\sum\limits_{n}(x_{n}-\bar{x})(y_{n}-\bar{y})}{\sqrt{\sum\limits_{n}(x_{n}-\bar{x})^{2}}\sqrt{\sum\limits_{n}(y_{n}-\bar{y})^{2}}}$
(5)

where $x_{n}$ is the value for a given voxel and temperature, $\bar{x}$ is the
mean of the entire array for that temperature, $y_{n}$ is the value for the
same voxel at another temperature and $\bar{y}$ is the mean of the entire
array at that temperature. By definition, $r=1$ for the diagonal because the
same array is perfectly correlated to itself. Before computing $r$, each
crystal was shifted such that its was in the middle of the $256\times
256\times 256$ voxel array, which was subsequently cropped to a $98\times
186\times 126$ voxel array, still fully encapsulating the crystal, to speed up
$r$ computation.

The $r$ matrix is shown in Fig. 7 for $\mathbf{u_{\textit{hkl}}(r)}$. Before
using Eq. 1 to calculate $\mathbf{u_{\textit{hkl}}(r)}$, phase ramps were
removed by re-centring the Fourier transform of the complex electron density
and a phase unwrapping algorithm [Cusack2002] was used to remove phase wraps.
Furthermore, we accounted for arbitrary phase offsets by establishing a
$40\times 40\times 40$ voxel volume centred about the bulk ROI’s centre of
mass as the zero-phase reference for every reconstruction, similar to the
procedure adopted for strain in Section 2.3.

Figure 7: Reflection-specific $r$ matrices between
$\mathbf{u_{\textit{hkl}}(r)}$ for different temperatures. The difference in
$r$ values corresponds to structural evolution. The $(002)$ reflection for
crystal A shows an abrupt $r$ change after 410°C, suggesting the abrupt
removal of most defects. This is less evident in the $(11\bar{1})$ reflection
for crystal B, which shows many subtle structural changes as a function of
temperature. The red dashed lines indicate where the temperature sampling
changed. Black arrows represent the most noticeable changes in displacement.

Structural transitions in the sample are identified by sharp changes in $r$
presented in Fig. 7. Surprisingly the two reflections from crystal A have
different temperatures where the most significant changes occur: abruptly
after 410°C in the $(002)$ reflection and less abruptly at 380°C in the
$(1\bar{1}1)$ reflection. These transitions are also reflected in the facet
area plot for crystal A in Fig. 6. The discrepancy between transition points
could be due to the $(00\bar{1})$ facet forming at 380°C but grows
significantly larger after 410°C.

Some of these abrupt changes in $r$ occur as a result of sudden relaxation
events associated with the removal of a defect, such as dislocations (3),
removed at 372°C, and (4), removed at 380°C, from Fig. 4 and supplementary
video SV2. In crystal A the $(002)$ reflection has reached a steady annealed
state after 410°C, but this annealed state is less defined in the other
reflections. Examining the same family of reflections can lead to different
conclusions, as crystal B shows less abrupt transitions within the $r$
matrices (especially the $(200)$ reflection) that suggest a more gradual
recovery of the lattice. In addition to these prominent transitions within the
$r$ matrix, less pronounced fluctuations are present and here are attributed
to the oscillation between intermediate annealing states.

The differences in observations between reflections from the same crystal
arises because not all defects are visible in every reflection and the
projection of associated displacement fields varies for each combination of
Burgers vector and scattering vector. This variance in material transformation
points highlights the importance of multi-reflection analysis, as it can be
misleading to derive conclusions from just a single $r$ matrix.

###  3.4 Normalised morphology cross-correlation matrix

To separate morphology-related changes from lattice displacement changes, the
$r$ correlation matrix of the morphology is shown in Fig. 8. The morphology of
a crystal was created by applying a 0.30 threshold to the reconstructed
amplitude.

Figure 8: Reflection-specific $r$ matrices for the morphology at different
temperatures. The difference in $r$ values corresponds to morphology
evolution. There is a more gradual $r$ transition as the temperature increases
compared to Fig. 7. The red dashed lines indicate where the temperature
sampling changed. Black arrows represent the most noticeable changes in
displacement as observed in Fig. 7.

There is good agreement between the morphology $r$ matrices and the facet area
plots. For instance, the disappearance of the $(00\bar{1})$ facet in crystal B
shown in Fig. 6 after 527°C can also be clearly identified in the
$(11\bar{1})$ reflection in Fig. 8. Notable changes in $r$ for
$\mathbf{u_{\textit{hkl}}(r)}$ in Fig. 7, marked by hollow black arrows in
Fig. 7 and Fig. 8 are different to changes in $r$ for the morphologies shown
in Fig. 8. There are less abrupt $r$ changes in Fig. 8, suggesting that
annealing causes smaller incremental changes in morphology, different from the
sudden changes in $\mathbf{u_{\textit{hkl}}(r)}$ (see Fig. 7). However, there
are changes in $\mathbf{u_{\textit{hkl}}(r)}$ that cannot be seen in the facet
area plots in Fig. 6 or the morphology matrices in Fig. 8, namely around
$\sim$310°C in the $(1\bar{1}1)$ reflection in Fig. 7. Changes in strain that
cannot be explained by changes in morphology can be examined directly in the
next section.

###  3.5 Average lattice strain in regions of interest

We will now consider the average strain within particular regions of the
crystal, illustrating key differences in the local lattice environment. The
average lattice strain for a reflection was computed by calculating the strain
for each voxel using Eq. 4, and taking the mean over all voxels with an
electron density magnitude greater than 0.30. Each crystal has been divided
into four ROIs, shown in Fig. 9, according to how the Ga-beam has interacted
with the sample. The ROIs are labelled glancing, normal, corner and bulk. The
glancing ROI is defined as 100 nm (20 voxels) from the surface of the milled
face where the Ga beam impinged on the sample at glancing incidence. Although
the FIB milling was intended to be normal to the substrate, the milled face is
$\sim$80° from the local sample surface normal. This is caused by the tail of
the Gaussian-like FIB probe profile not being able to fully reach the bottom
of the crystal. The normal ROI is 100 nm from the top of the crystal,
corresponding to the region where the FIB beam arrives approximately normal to
the sample surface. The corner ROI is the volume that joins the glancing and
normal regions. Previously the corner was seen to have the most defects
[Hofmann2018a]. The bulk ROI is the remaining part of the crystal. The average
lattice strain for the different ROIs at every temperature point is shown in
Fig. 9.

Figure 9: Average lattice strain in different ROIs as a function of
temperature. Vertical bars correspond to the standard deviation of the strain.
Over the annealing cycle, the average strain in each reflection gradually
decreases and stagnates after $\sim$410°C, as does the magnitude of the
standard deviation. The magnitude and direction of the average strain is
reflection-dependent, especially in the corner ROIs. The coordinate axes are
positioned in sample space and plotted with a length of 250 nm. The red dashed
lines indicate where the temperature sampling changed. The rendered crystals
show the locations of the ROIs.

In Fig. 9, the strain within the bulk ROI remained approximately zero for all
reflections and temperatures because no FIB-induced strain was present prior
to annealing. The room temperature data is similar to the values reported
previously [Hofmann2018a]. As the temperature increases, we see that FIB
damage is removed, and the average strain converges to zero. We observe that
the glancing and normal ROIs behave similarly regardless of reflection, with
the normal ROI having more compressive strain than the glancing ROI. This
suggests that the FIB tail imparts more damage to the top face of the crystal
than does glancing incidence milling to the milled face, in agreement with
previous work [Hofmann2018a, Hofmann2017a]. The most reflection-dependent ROI
is the corner, where both the average magnitude and the sign of strain vary.
For the corner ROI in crystal A, the $(002)$ reflection shows tensile strain
whereas the $(1\bar{1}1)$ reflection shows compressive strain of similar
average magnitude for the same crystal. The average strain magnitude can be
dependent on the reflection. This is observed in the corner ROI for crystal B
at room temperature, where the $(200)$ reflection has nearly twice the average
strain magnitude of the $(11\bar{1})$ reflection. These significant
differences would not be evident from single reflection BCDI experiments.

The average strain curves approach zero at $\sim$410°C in Fig. 9. This is
above the Tammann temperature (327°C) and matches the temperature where facet
change is observed for crystal A in Fig. 6 and where structural changes occur
for the $(002)$ reflection for crystal A in Fig. 7. Above $\sim$410°C there is
little strain evolution, indicating that most of the damage in the crystal has
been annealed out. Further evidence for this transition point can be seen in
the standard deviation magnitudes of Fig. 9. Below $\sim$410°C, the standard
deviation magnitudes for the glancing, normal and corner ROIs (plotted in
Appendix F) are greater than that of the bulk ROI. Above $\sim$410°C, the
standard deviations across all ROIs are similar, indicating that the average
strain across all ROIs has become equally homogeneous. The removal of strain
inside the crystals can also be seen by comparing the maximum intensity slices
of the CXDPs in Appendix G and noting the change in fringe sharpness before
and after annealing. The evolution of these CXDP slices can be seen in
supplementary videos SV4 and SV5.

We also note that the differences in the rate of average strain removal seem
to depend on the severity of the initial FIB damage. It appears that having a
higher initial average strain magnitude leads to a greater initial change in
average strain. This is most apparent when comparing the rate of average
strain change in the corner, glancing and normal ROIs for any reflection, as
well as solely the corner ROIs in crystal B.

###  3.6 Analysis technique comparisons

It is interesting to compare information extracted from the facet area plots,
$r$ matrices and average strain plots. While we see that the facet area plots
and displacement $r$ matrices can highlight distinct transition points, only
major changes can be confidently identified. In particular, crystal
transitions can be readily identified in the facet area plots and displacement
$r$ matrices but remain harder to resolve when observing average strain plots
and morphology $r$ matrices. The primary advantage of the average strain plots
over the consideration of $r$ matrices is the ability to identify the
magnitude and the direction of strains in various ROIs, which can confirm
whether the FIB-induced strain has indeed been removed.

The seemingly abrupt structural evolution in the facet area plots and
displacement $r$ matrices contrasts the more uniform rate of strain change
observed in the average strain plots and morphology $r$ matrices. The facet
area plots and displacement $r$ matrices portray annealing as a process
controlled by sudden substantial changes, implying a rapid evolution in
morphology and removal of strain after a certain temperature is reached, as
shown by the formation of facets in Fig 6 and the $r$ for
$\mathbf{u_{\textit{hkl}}(r)}$ in Fig 7. For the same reflection in the
morphology $r$ matrices in Fig 8 and in the average strain plots in Fig. 9,
annealing appears as a gradual process, by which the change in temperature
leads to a progressive change in lattice relaxation and facet edge softening.
The trench-like dislocations observed in Fig. 4 gradually decrease in size but
also fully disappear once certain temperatures are reached, seen in
supplementary video SV2, exhibiting both gradual and sudden changes in crystal
structure as a result of annealing.

The accuracy of these observations is limited by our spatial resolution and
sampling time. To better resolve defect dynamics, faster scans and higher
photon flux are required, which the Advanced Photon Source Upgrade Project
promises to deliver. Regardless, the conclusions drawn from each of the
analysis techniques are complementary and reveal different aspects of the
crystal behaviour, but point towards an annealing temperature of 380-410°C.
These methods are suitable for the analysis of individual reflections, but can
benefit from the incorporation of multiple reflections, enabling more
comprehensive comparisons to be made.

###  3.7 Diffusion considerations

Scan times on the order of minutes for BCDI rocking curves limit fast
structural changes from being resolved. However, we are able to capture
structural evolution on the timescale of the rocking scans, such as the
formation of facets, disappearance of dislocations and relative changes in
$\mathbf{u_{\textit{hkl}}(r)}$. These changes can be interpreted in terms of
the diffusion of Ga in Au, the self-diffusion of Au at the surface or within
the bulk of the crystal and the escape of dislocations to the free surface.

Fig. 10 shows the temperature and time required for diffusion of Ga over
lengths from 1 nm to 20 nm, as well as self-diffusion of Au over 1 nm to 100
nm. We use diffusion coefficients for these two systems as they provide an
indication of relative atomic mobility at various temperatures. The diffusion
coefficient of Ga in Au was obtained from [Gupta1987], while the self-
diffusion of Au was obtained from [Makin1957]. We determine the time required
for 1 nm of diffusion as a function of temperature, for each scenario, as
shown by the exponential curves in Fig. 10. The penetration depth of Ga under
the FIB conditions used is $\sim$20 nm [Hofmann2018a], thus we calculate the
time required for a diffusion length of 20 nm as the upper limit of Ga
diffusion. Diffusion lengths between 1 and 20 nm are shown as a grey region
between the 1 nm and 20 nm curves. Similarly, for the self-diffusion of Au,
diffusion lengths between 1 and 200 nm are relevant, as the greatest change in
facet area is about $4000\mathrm{\ nm^{2}}$ (Fig. 6). Superimposed on this
time-temperature plot is the temperature history of the present Au samples
(dotted blue line), indicating the time spent at each temperature during the
annealing process.

Figure 10: The time required for the bulk diffusion of Ga over lengths from 1
to 20 nm (light grey region) and Au diffusion over distances from 1 to 200 nm
(dark yellow region) throughout the annealing. Superimposed are temperature
points where dislocations in Fig. 4 disappeared, temperature points where
facets have formed in Fig. 6, and temperature regions where the
$\mathbf{u_{\textit{hkl}}(r)}$ $r$ values changed noticeably (see Fig. 7). The
dotted blue line shows the time-temperature trajectory of the Au sample in
this study, indicating the time spent at each annealing temperature. The red
dashed lines indicate where the temperature sampling changed.

To Fig. 10 we can further add the changes in $\mathbf{u_{\textit{hkl}}(r)}$
cross-correlation observed as a function of temperature (Fig. 7). From
$\sim$190°C until $\sim$410°C there is a gradual reduction in the cross-
correlation value. Considering the diffusion lengths of Ga and Au, this change
at lower temperatures up to $\sim$280°C is likely to be driven by bulk Ga
diffusion and Au surface diffusion. Above this temperature bulk Au self-
diffusion may be responsible. At $\sim$410°C Fig. 7 shows a rapid change in
cross correlation coefficient, suggesting major structural changes as
discussed above. The diffusion curves in Fig. 10 suggest that this is
associated with the rearrangement of Au rather than Ga atoms.

On Fig. 10 we have superimposed the temperatures at which removal of
dislocations is observed (black lines labelled 1 to 4, consistent with
dislocation numbers in Fig. 4), a well as the temperatures at which the
$(00\bar{1})$ and $(0\bar{1}\bar{1})$ facets are formed (green and purple
lines respectively, consistent with Fig. 5 and 6). Apart from the removal of
dislocation 1, which is at the boundary between Ga diffusion and Au diffusion
dominated areas, all other structural changes appear to be driven by Au
migration.

##  4 Conclusion

We have presented an in situ BCDI study of the annealing of FIB-induced damage
in Au microcrystals. Two separate crystals with two reflections each were
measured, confidently switching between reflections and samples, during an
annealing cycle. By tracking the change in detector position, we could use the
thermal strain of the microcrystals for temperature calibration. We have also
presented a new approach that monitors the reduction of free surface energy by
following morphology evolution through the formation of facets. The annealing
process was analysed using facet areas, cross-correlation matrices and average
strain plots, which provide complementary information. Surprisingly, two
different reflections from the same crystal can show different behaviour. This
highlights the need to measure multiple reflections in order to gain a
complete picture of sample evolution. Full annealing was achieved at a
temperature of 380-410°C. However, the complexity of the process becomes
apparent when trying to evaluate the rate of strain removal, as the facet area
plots and displacement $r$ matrices portray annealing as an abrupt process,
while the morphology $r$ matrices and average strain plots show much more
gradual changes. The gentle evolution of the crystal below $\sim$280°C is
likely the result of bulk Ga diffusion and Au surface diffusion whilst more
substantial changes (i.e. removal of dislocations and formation of facets)
above $\sim$280°C are likely dominated by the self-diffusion of Au. This
experiment demonstrates that annealing is a viable option for the removal of
FIB damage as evident by a decrease in strain heterogeneity inside the
crystal, removal of defects on the crystal surface and reduction of mean
strain in FIB-damaged areas. Annealing can thus be used as a method to remove
FIB-induced sample damage prior to further modification of a sample e.g.
nanoindentation, bending etc.

We anticipate that the techniques developed here could be used to pioneer
novel in situ BCDI measurements in complex environments for multiple samples
and reflections simultaneously. The processed diffraction patterns, final
reconstructions and data analysis scripts are publicly available at:
https://www.doi.org/10.5281/zenodo.4423216.

##  Appendix A 34-ID-C heater cell

Fig. 11(a) shows the heater within a furnace enclosure. Fig. 11(b) shows a
rendering of the enclosure cover, consisting of a stainless-steel cover plate
attached to a base using stainless-steel threaded rods and nuts. Only two rods
were used to minimise the obstruction of incident and reflected X-rays. The
circular hole in the plate was covered with a quartz window to allow
observation of the sample with a confocal microscope with a $5\times$
objective lens during the annealing experiment. The sides of the furnace
enclosure cover were covered with metalized Kapton (not shown) to allow for
transmission of X-rays and favorable thermal transport for cooling of the
film.

Figure 11: (a) The sample clipped in the heater cell as seen through the
heater enclosure window. Note the thermocouple is attached to the clip, which
contacts the silicon substrate. This results in thermal resistance between the
Au crystals and the thermocouple, causing the discrepancy between the measured
and actual temperature of the sample. (b) A rendering of the furnace enclosure
cover. The scale bar corresponds to 20 mm.

##  Appendix B Effect of seeding on the retrieved electron density

This section addresses how the seeding of the phase retrieval algorithm
influences the final object’s amplitude and strain. For a CXDP at an arbitrary
temperature, we consider the following options:

  1. 1.

Seeding with the reconstruction from the temperature step immediately below
the current temperature, except if the current temperature is 23°C, in which
case a random guess is used.

  2. 2.

Seeding with the reconstruction from the temperature step immediately above
the current temperature, except if the current temperature is 587°C, in which
case a random guess is used.

  3. 3.

Seeding with a random guess regardless of the current temperature.

A simpler phase retrieval algorithm compared to the one listed in Appendix C
was implemented to exaggerate the differences in the final reconstructions. A
pattern of 20 error reduction (ER) and 180 hybrid input-output (HIO)
iterations was repeated three times followed by 20 ER iterations, for a total
of 620 iterations. This algorithm was applied to each CXDP padded to a size of
$256\times 256\times 176$ voxels. After phase retrieval, each object was
mapped into orthogonal sample space and the strain was calculated using Eq. 4
and the procedure found at the end of Appendix C.

To compute the error between two reconstructions generated using two different
seeding approaches, we sum the difference between the two reconstructions for
each voxel, $n$, and normalise it over the union of the morphologies (created
using an amplitude threshold of 0.30) for that particular reflection and
temperature, This yields Eq. 6 and Eq. 7 for the strain error and amplitude
error respectively:

$\mathrm{E}_{\epsilon}=\frac{\sum\limits_{n}|\epsilon_{\mathrm{seed\
a}}-\epsilon_{\mathrm{seed\ b}}|}{\sum\limits_{n}\mathrm{Mask}}$ (6)

$\mathrm{E}_{\omega}=\frac{\sum\limits_{n}|\omega_{\mathrm{seed\
a}}-\omega_{\mathrm{seed\ b}}|}{\sum\limits_{n}\mathrm{Mask}}\\\ $ (7)

where seed a and b are any two of the three seeding scenarios listed
previously. Fig. 12 shows how the strain and amplitude for reconstructions
generated using each seeding variable compare.

Figure 12: The amplitude and strain error, computed using Eq. 6 and 7
respectively, between reconstructions obtained using different seeding
approaches. Each line represents the error between any two of the three
seeding variables. The red error line has the lowest error on average. The
large error spikes are the result of poorly converged solutions. The red
dashed lines indicate where the temperature sampling changed.

In Fig. 12 we observe that no two seeds generate reconstructions that are
exactly the same, which would otherwise be represented by a line with
negligible differences. Spikes on the plots point out that at least one
reconstruction did not converge as well as the others. At temperatures where
all error comparison lines overlap, the reconstructions are equally different
from each other. The average strain error is $\sim$2$\times 10^{-4}$, which
agrees well with the uncertainty seen in previous studies [Hofmann2017a,
Hofmann2018a, Hofmann2020]. The red lines in Fig. 12 represent the error
between between the two non-random seeds, and are in general lower than the
error between either seed and the random guess. This shows that phasing with a
seed generates reconstructions that converge more rapidly, compared to phasing
with a random guess. However we cannot conclude whether phasing with a
reconstruction from the next or previous temperature step yields a more
accurate result. Thus, we decided for each reconstruction to be seeded with
itself at the same temperature from the previous phasing round (except for the
first round, where a random guess was used for all temperatures), instead of
seeding with a reconstruction at a different temperature in the same round as
explained in Section 2.3.

##  Appendix C Phase retrieval

The reconstruction process for every CXDP was executed in four stages, similar
to the approach used by [Hofmann2020], using the output from the previous
stage to seed the next phasing stage (listed below). Note the $(002)$
reflection for crystal A and $(200)$ and $(11\bar{1})$ reflections for crystal
B at room temperature were padded to a size of $256\times 256\times 176$
voxels to make them consistent with the other CXDPs in the study.

  1. 1.

Each CXDP was seeded with a random guess. A guided phasing approach with 40
individuals and four generations was used [Chen2007]. The first and second
generations used low resolution data while the remaining generations used full
resolution data. For each generation, a block of 20 ER and 180 HIO iterations,
with $\beta=0.9$, was repeated three times. This was followed by 20 ER
iterations, which were averaged to return the final object. The best
reconstruction was determined using a sharpness criterion, as it is the most
appropriate metric for crystals containing defects [Clark2015, Ulvestad2017b].
The shrinkwrap algorithm [Marchesini2003] was used to periodically update the
real-space support.

  2. 2.

The reconstruction was seeded with the output from stage 1. This stage was
identical to the previous except now we account for the longitudinal and
transverse partial coherence [Clark2012] using the Richardson-Lucy
deconvolution algorithm.

  3. 3.

The reconstruction was seeded with the output from stage 2. A block of 20 ER
and 180 HIO iterations was repeated 15 times, followed by 1000 iterations of
ER. The final 50 ER iterations were averaged to produce the final output. The
same partial coherence settings in stage 2 were included, but guided phasing
was omitted.

  4. 4.

Stage 3 was repeated, using the output from stage 3 as the seed.

Next, phase offsets of 0, $-\frac{\pi}{2}$, and $\frac{\pi}{2}$ were applied
to the electron density for each reconstruction and the resulting phases
outside $-\pi$ to $\pi$ were returned to this range by adding or subtracting
$2\pi$. The resulting reconstructions for each of the three offsets were
transformed from detector conjugated space to orthogonal sample space with a
voxel size of $5\times 5\times 5\mathrm{\ nm^{3}}$. Phase ramps were
eliminated by centring the Fourier transform of the complex electron density.
This set of phase offset data was differentiated and the minimum gradient for
each voxel was used to calculate the strain, preserving its sign, to avoid
discontinuities which otherwise arise at dislocation cores [Hofmann2020].

##  Appendix D Spatial resolution

We note that the spatial resolution is dependent on the measurement directions
in reciprocal space and the geometry of the sample. The spatial resolution was
calculated by considering two line profiles that cross the air-sample
interface for each reflection: one normal to the top face and the other normal
to the milled face. The spatial resolution was determined by differentiating
the profiles and fitting the resulting peak with a Gaussian. The average
$2\sigma$ of the fitted Gaussian at each temperature step is reported as the
spatial resolution. A plot of the spatial resolution vs. temperature is shown
in Fig. 13.

Figure 13: Spatial resolutions with respect to temperature for each
reflection. The overall average for each crystal is the mean of the point
averages across all temperatures. The spatial resolution for this experiment,
23 nm, is the overall average for all reflections. The red dashed lines
indicate where the temperature sampling changed.

##  Appendix E Dislocations in crystal B, (200) reflection

The phase for crystal B’s $200$ reflection at room temperature is shown in
Fig. 14. Vortex-like phase features that resemble dislocations are
highlighted. These features are present in the trench-like regions shown in
Fig. 4.

Figure 14: A phase isosurface rendering of crystal B’s $200$ reflection at
room temperature, showing phase vortices reminiscent of dislocations marked by
black arrows. The coordinate axes are positioned in sample space and plotted
with a length of 200 nm. Here the amplitude threshold is 0.04 and no phase
unwrapping was performed.

##  Appendix F Standard deviation of average strain

The standard deviation of the strain for a particular ROI provides information
about strain heterogeneity within the ROI. This can be used to evaluate the
extent of annealing, whereby a lowered standard deviation corresponds to a
more relaxed crystal lattice. The strain standard deviation vs. temperature
for all ROIs is plotted in Fig. 15.

Figure 15: Standard deviation of strain for different ROIs at each temperature
point. As annealing proceeds, the standard deviation follows a decreasing
trend in all ROIs, which shows a reduction in strain heterogeneity across all
reflections. After $\sim$410°C, all standard deviation values are fairly
constant, suggesting that FIB-induced strain has been removed.

##  Appendix G CXDP comparison

Differences in the collected CXDPs at different temperatures show evidence of
structural change. The asymmetry in the Bragg peak is indicative of strain
within the crystal lattice. A comparison of the CXDPs at 23°C and 587°C is
shown in Fig. 16.

Figure 16: Central slice of the unprocessed CXDPs at the lowest and highest
temperatures for each reflection. The asymmetry in the CXDP slices before
annealing indicates a high degree of strain. This strain is alleviated after
heat treatment, which is demonstrated by the centrosymmetric fringes at 587°C.
The evolution of the raw CXDPs for each reflection during heat treatment can
be seen in supplementary videos SV4 and SV5.

##  Appendix H Supplementary video descriptions

Below is a list of brief descriptions for supplementary videos.

  * •

SV1 shows the evolution of the surface morphology and internal strain for
crystal A based on Fig. 4.

  * •

SV2 shows the evolution of the surface morphology and internal strain for
crystal B based on Fig. 4.

  * •

SV3 shows facet area evolution as temperature increases for both crystals
based on Fig. 5 and Fig. 6.

  * •

SV4 shows the evolution of the CXDP slices for crystal A based on Fig. 16.

  * •

SV5 shows the evolution of the CXDP slices for crystal B based on Fig. 16.

Acknowledgements The authors would like to thank Ian K. Robinson for providing
the patterned Au microcrystals on silicon substrates and Ruqing Xu and Wenjun
Liu at 34-ID-E for help with Laue diffraction measurements. The furnace
enclosure was originally designed by Evan R. Maxey and Wonsuk Cha and the
solid angle enhancement and confocal compatibility was updated by Nicholas W.
Phillips. D.Y, N.W.P and F.H. acknowledge funding from the European Research
Council under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 714697). K.S. acknowledges funding from the
General Sir John Monash Foundation. The authors acknowledge use of
characterisation facilities at the David Cockayne Centre for Electron
Microscopy, Department of Materials, University of Oxford and use of the
Advanced Research Computing (ARC) facility at the University of Oxford
[Richards2015]. X-ray diffraction experiments were performed at the Advanced
Photon Source, a US Department of Energy (DOE) Office of Science User Facility
operated for the DOE Office of Science by Argonne National Laboratory under
Contract No. DE-AC02-06CH11357.

iucr


# Investigations into the characteristics and influences of nonequilibrium
evolution

Xiaobing Li  Key Laboratory of Quark and Lepton Physics (MOE) and Institute of
Particle Physics,
Central China Normal University, Wuhan 430079, China  Mingmei Xu
<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MOE) and
Institute of Particle Physics,
Central China Normal University, Wuhan 430079, China  Yanhua Zhang  Department
of physics and electronic engineering, Yuncheng University, Yuncheng 044000,
China  Zhiming Li  Key Laboratory of Quark and Lepton Physics (MOE) and
Institute of Particle Physics,
Central China Normal University, Wuhan 430079, China  Yu Zhou  Department of
Mathematics, University of California, Los Angeles, California 90095, USA
Jinghua Fu  Key Laboratory of Quark and Lepton Physics (MOE) and Institute of
Particle Physics,
Central China Normal University, Wuhan 430079, China  Yuanfang Wu
<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MOE) and
Institute of Particle Physics,
Central China Normal University, Wuhan 430079, China

###### Abstract

In order to estimate qualitatively the influence of nonequilibrium evolution
in relativistic heavy ion collisions, we use the three dimensional Ising model
with Metropolis algorithm to study the evolution from nonequilibrium to
equilibrium on the phase boundary. The evolution of order parameter approaches
its equilibrium value exponentially, the same as that given by Langevin
equation. The average relaxation time is defined which is demonstrated to well
represent the relaxation time in dynamical equations. It is shown that the
average relaxation time at critical temperature diverges as the $z$th power of
system size. The third and the fourth cumulants of order parameter during the
nonequilibrium evolution could be either positive or negative, depending on
the observation time, consistent with dynamical models at $T>T_{\rm c}$. It is
found that the nonequilibrium evolution at $T>T_{\rm c}$ lasts very short, and
the influence is weaker than that at $T<T_{\rm c}$. Those qualitative features
are instructive to determine experimentally the critical point and the phase
boundary of QCD.

##  I Introduction

One goal of current relativistic heavy ion collision experiments is to map the
QCD phase diagram as a function of temperature $T$ and baryon chemical
potential $\mu_{\rm B}$. QCD phase transition is a nonperturbative process
whose analytical solutions are difficult. Lattice QCD is restricted to
vanishing or small baryon chemical potential, where the hadron phase smoothly
transitions to QGP phase lab1 . QCD-based models predict the phase transition
is first order at low $T$ and high $\mu_{\rm B}$ first . The first order phase
transition line ends at a critical point cp . The critical sensitive
observables are high cumulants of conserved charges high-cumu ; high-cumu2 .

Currently, at the RHIC beam energy scan, non-monotonous behavior of the fourth
cumulants (kurtosis) of net protons has been observed STAR . While, the third
cumulants (skewness) of net protons are always negative below the Poisson
baseline STAR , in contrast to the positive values obtained by the equilibrium
expectations Stephanov-1999 ; Stephanov-2011 ; Asakawa . The negative values
may be caused by nonequilibrium effects Mukherjee-2015 ; Nahrgang-2019 ;
Song-2019 ; Rajagopal-2019 .

Usually, equilibrium is assumed in lattice QCD and various QCD-based models
Mukherjee-2015 ; Nahrgang-2019 ; Berdnikov-2000 . However, in reality,
nonequilibrium is inevitable. The fireball formed in a heavy ion collision
only spends a limited time. The initial state of the fireball may reach, or
not reach thermal equilibrium PBM-nature ; equi-exp .

If the collision energy is high enough and the initial state reaches the
temperature needed for the deconfinement phase transition, the fireball
transforms to a deconfined phase consisting of quarks and gluons. In this
case, the QGP phase is formed and the system reaches equilibrium equi-exp . It
is also possible that local equilibrium is reached and QGP is formed in
droplets. Another possibility is that some collisions in heavy ion experiments
do not experience the phase transition, and the system does not reach
equilibrium. Therefore, equilibrium, local equilibrium and nonequilibrium are
three possible initial states.

In addition, at the critical region, even if the initial state reaches
equilibrium, the critical slowing down prevents the system from equilibrium
Berdnikov-2000 .

For the statistics of equilibrium state, phase transition is due to
spontaneous symmetry breaking. Systems sharing the same symmetry have the same
critical exponents, and belong to a universality class. Critical points of
both QCD and 3D Ising model belong to the Z(2) symmetry group univers-1 ; lab5
; univers-3 ; univers-4 ; Z2-1 . So in order to learn the critical phenomena
of QCD, people often refer to the 3D Ising model.

The singular part of the QCD free energy density can be obtained by means of
the scaling function of the Ising model Mukherjee-2015 ; scaling . The mapping
of the Ising phase diagram to that of the QCD is also constructed
Stephanov-2011 . The mapping makes the temperature axis of the Ising model
tangent to the first order phase transition line of the QCD, ensuring the
coincidence of critical points and approximate coincidence of first order
phase transition lines mapping2 . High cumulants of conserved charges in
relativistic heavy ion collisions are approximately corresponding to those of
magnetization of 3D Ising model.

As far as nonequilibrium evolution is concerned, there are no ready theories
at present. The common feature of nonequilibrium evolution is that relaxation
time at critical point diverges by the power of dynamic exponent. The dynamic
exponent is system and algorithm dependent.

Dynamical evolution equations, e.g. Langevin dynamics Song-2019 ;
Berdnikov-2000  and various relaxational models Mukherjee-2015 ; Nahrgang-2019
, are usually used to estimate the effects of nonequilibrium in the QCD
critical region. However, the solutions of dynamical evolution equations are
restricted to the region of crossover.

As we know, the numerical simulation of the Ising model with Metropolis
algorithm is suitable for studying nonequilibrium evolution lab20 lab26 . The
relaxation processes can be easily realized on the whole phase boundary. We
find that the nonequilibrium evolution of the order parameter in the Ising
model indeed approaches exponentially to its steady value, the same as that of
the Langevin equation Landau . The average relaxation time at critical
temperature diverges as the $z$th power of system size, the same as relaxation
time in dynamical equations lab33 . Moreover, the third and fourth cumulants
of order parameter on the crossover side oscillate around zero, and then
converge to their equilibrium values. The sign of the third cumulants can be
negative, consistent with those obtained from refs. Mukherjee-2015 ;
Nahrgang-2019 .

In this paper, we first demonstrate the characteristics of nonequilibrium
evolution and introduce a time scale of nonequilibrium evolution. Then the
influences of nonequilibrium evolution on observables are presented. Section
II gives a short introduction to Ising model and Metropolis algorithm. In
Section III, the dependences of average relaxation time on temperature, system
size and initial configuration are investigated. Characteristics of
nonequilibrium evolution are presented. Section IV presents the time evolution
of cumulants of order parameter. Influences of nonequilibrium on observables
are discussed. A summary is given in section V.

##  II Ising Model and Metropolis algorithm

The 3D Ising model considers a three dimensional cubic lattice composed of
$N=L^{3}$ sites, where $L$ is called the system size. Every site $i$ is
occupied by a spin, $s_{i}$. The spins can be in one of two states, either
spin-up, $s_{i}=+1$, or spin-down, $s_{i}=-1$. The state of the system can be
represented by a series of spins, i.e.

$\\{s_{1},s_{2},\cdots,s_{N}\\}.$ (1)

A shorthand notation $\\{s_{i}\\}$ is used in the following.

The spins at positions $i$ and $j$ interact with one another. For a pair of
parallel spins we assign an interaction energy of $-J$, while for a pair of
anti-parallel spins we assign an interaction energy of $+J$. Only interactions
with the nearest neighbors are considered.

The spins also interact with an external magnetic field $H$. The total energy
of a system of $N$ spins with constant nearest-neighbor interactions $J$
placed in a uniform external field $H$ is

$E_{\\{s_{i}\\}}=-J\sum_{\langle ij\rangle}s_{i}s_{j}-H\sum_{i=1}^{N}s_{i},$
(2)

where the notation $\langle ij\rangle$ restricts the sum to run over all the
nearest neighbor spins.

Then the partition function is

$Z(T,H)=\sum_{\\{s_{i}\\}}{\rm exp}(-E_{\\{s_{i}\\}}/k_{\rm B}T),$ (3)

where $k_{\rm B}$ is Boltzmann’s constant. The free energy is evaluated by

$F(T,H)=-k_{\rm B}T{\rm ln}Z.$ (4)

The average total magnetization is

$M=-(\frac{\partial F}{\partial H})_{T}=\langle\sum_{i=1}^{N}s_{i}\rangle,$
(5)

and the per-spin magnetization is

$m=\frac{1}{N}\sum_{i=1}^{N}s_{i}.$ (6)

We here focus on the behavior of the Ising model as a function of temperature
$T$, which defines an energy scale $k_{\rm B}T$. For temperatures $k_{\rm
B}T\ll J$, the spin-spin interactions are relatively strong, so that the spins
tend to align with one another, and $|m|$ is close to unity. This is an
ordered phase. For temperatures $k_{\rm B}T\gg J$, the spin-spin interactions
are relatively weak, so that the spins are effectively non-interacting and
point up and down randomly, and $|m|$ is close to zero. This is a disordered
phase. A phase transition from a high-temperature disordered phase to a low-
temperature ordered phase is anticipated which is continuous and called the
Curie point. The Curie point is the critical point on the plane of variables
$T$ and $H$. The temperature of the critical point is denoted by $T_{\rm c}$.

The value of $T_{\rm c}$ of the two dimensional Ising model was exactly
calculated by Onsager in 1944 lab28 . $T_{\rm c}$ of the 3D Ising model is
estimated by the finite size scaling theory and analysis of magnetization
distribution. $K_{\rm c}=J/k_{\rm B}T_{\rm c}=0.2216544$ is obtained lab31 ,
which agrees very well with the results obtained from the renormalization
group theory lab32 . Usually, $J$ and $k_{\rm B}$ are set to 1, so $T_{\rm
c}=4.51$, which is used in the following.

Metropolis algorithm was introduced by Nicolas Metropolis and his
collaborators in their paper in 1953 lab27 . The main steps of Metropolis
algorithm in the 3D Ising model are as follows.

  1. (1)

Set the temperature $T$ and the number of sites $N$.

  2. (2)

Generate the initial configuration.

Two kinds of initial configuration are usually used. One is random
configuration with all spins pointing randomly up or down, while the other one
is polarized configuration with all spins pointing in the same direction.

  3. (3)

Test a single spin whether flips or not.

Whether a spin is flipped depends on the acceptance probability
$A({\boldsymbol{u}}\rightarrow{\boldsymbol{v}})$, which is given by

$A({\boldsymbol{u}}\rightarrow{\boldsymbol{v}})=\left\\{\begin{array}[]{ll}e^{-(E_{\boldsymbol{v}}-E_{\boldsymbol{u}})/k_{\rm
B}T}&\text{if $E_{\boldsymbol{v}}-E_{\boldsymbol{u}}>0$},\\\
1&\text{otherwise.}\end{array}\right.$ (7)

${\boldsymbol{u}}$ and ${\boldsymbol{v}}$ represent the state of the system
before and after flipping this spin.

If $A({\boldsymbol{u}}\rightarrow{\boldsymbol{v}})=1$, the spin is flipped.

If $A({\boldsymbol{u}}\rightarrow{\boldsymbol{v}})<1$, a random number $r$
($0<r<1$) is generated. If $A({\boldsymbol{u}}\rightarrow{\boldsymbol{v}})>r$,
the spin is flipped. Otherwise, the spin keeps its original state.

The testing of one single spin is called a Monte Carlo step.

  4. (4)

When $N$ Monte Carlo steps are completed, every spin in the lattice has been
tested for flipping and we say one sweep is completed. In this way, the
configuration of the system is updated once a sweep.

  5. (5)

After evolving enough sweeps, the magnetization approaches a steady value and
the system reaches equilibrium. Thermodynamic quantities are usually measured
in the equilibrium state.

  6. (6)

Change the temperature or the system size and repeat the above steps.

Metropolis algorithm flips one single spin at a step to get a new state. The
Wolff algorithm flips one cluster at a step. While different algorithms give
the same equilibrium state, the ways in which the system comes to equilibrium
are different. In this sense, one algorithm acts as one dynamic, and realizes
one relaxation process from nonequilibrium to equilibrium.

The dynamic exponent $z$ of the Metropolis algorithm and the Wolff algorithm
for 3D Ising model is 2.02 and 0.33, respectively lab17 . Due to the large
dynamic exponent, the Metropolis algorithm is more suitable for studying
relaxation processes lab20 lab26 .

By using Metropolis algorithm, we simulate the evolution of the 3D Ising model
from nonequilibrium to equilibrium at vanished external field and obtain
samples for the next calculations.

##  III Characteristics of nonequilibrium evolution

Starting from an initial configuration, the Ising system can evolve to an
equilibrium state of a given temperature spontaneously. The evolution before
reaching equilibrium is called the relaxation process, or nonequilibrium
evolution.

Starting from random configurations, we simulate 5000 evolution processes for
each temperature at a system size $L=60$. Fig. 1(a)-(d) show the time
evolution of $|m|$ (the absolute value of $m$ is used because the sign of the
magnetization is random at $H=0$) at four temperatures, i.e. $T/T_{c}=0.93$,
0.99, 1.00 and 1.03, respectively. The first two temperatures are
representative values on the first order phase transition line near the
critical point. $T/T_{c}=1.00$ represents the critical point. $T/T_{c}=1.03$
is a representative value on the crossover side. The horizontal axis is time
which is defined by the number of sweeps introduced in the fourth step of
Metropolis algorithm in section II. The red curve and the blue curve are the
results of two evolution processes randomly selected from the sample.

Figure 1: (a)-(d): The evolution of $|m|$ with time at $L=60$ for random
initial configuration at temperatures $T/T_{\rm c}=0.93$, 0.99, 1.00 and 1.03,
respectively. The blue curve and the red curve represent two evolution
processes randomly selected from the sample. (e)-(h): The distribution of
relaxation time $\tau_{\rm eq}$ for random initial configuration at four
temperatures. (i)-(l): The distribution of relaxation time $\tau_{\rm eq}$ for
polarized initial configuration at four temperatures.

In Fig. 1(a), $|m|$ shows an increasing trend at the beginning, and then gets
close to a steady value. After that, $|m|$ fluctuates slightly around the
steady value. The steady value represents an equilibrium state, and the
approaching to the steady value represents the relaxation to equilibrium.

We use $\mu$ to denote the steady value, i.e. the equilibrium expectation.
When the difference with the equilibrium expectation is considerably larger
than the root of the variance at equilibrium
$\sigma=\langle(x-\mu)^{2}\rangle^{1/2}$, the system is far from equilibrium.

The period of time from nonequilibrium to equilibrium is called the relaxation
time lab17 lab18 . In our simulation, relaxation time of the $i^{\rm th}$
evolution process $\tau_{\rm eq}^{i}$ is estimated by the time when the value
of $|m|$ enters the interval ($\mu-\sigma$, $\mu+\sigma$), i.e. the band of
thermal fluctuations around the equilibrium expectation.

Relaxation time of the $i^{\rm th}$ process defined above represents the
number of iterations needed to achieve equilibrium. One iteration is counted
after each of the spins is examined to flip or not by given dynamics. We count
the number of iterations, which is an integer. The number of iterations
presents the steps, i.e. time, that the system needs to achieve equilibrium.

In Fig. 1(a) the steady values of the red curve and the blue curve are the
same, but relaxation time of the red curve is much longer than the blue curve.
The difference of relaxation time between the red curve and the blue curve is
significant at low temperature as Fig. 1(a) shows and seems diminishing at
high temperature as Fig. 1(b)-(d) show. In order to show the difference of
relaxation time of different evolution processes, the distributions of
$\tau_{\rm eq}$ are plotted in Fig. 1(e)-(h). For the sake of comparison, the
horizontal ordinate of the four figures are set to the same.

At the temperature $T/T_{\rm c}=0.93$, the distribution of $\tau_{\rm eq}$ has
a long tail, as shown in Fig. 1(e). The long tail means there are a fraction
of evolution processes whose relaxation time is very long. When the
temperature gets closer to $T_{\rm c}$, the distribution gets narrow, as shown
in Fig. 1(f). At the critical temperature, the distribution gets wide again,
as Fig. 1(g) shows. The number of systems with a long relaxation time
increases at the critical temperature. On the crossover side, i.e. $T/T_{\rm
c}=1.03$, the width of the distribution is the smallest and the distribution
is concentrated at very short relaxation time, as Fig. 1(h) shows.

To quantify the relaxation time of different evolution processes, we define
the average relaxation time as

$\bar{\tau}_{\rm eq}=\frac{1}{n}\sum_{i=1}^{n}\tau_{\rm eq}^{i},$ (8)

where $n$ is the total number of evolution processes, $\tau_{\rm eq}^{i}$ is
relaxation time of the $i^{\rm th}$ evolution process. At $\bar{\tau}_{\rm
eq}$, not all systems are at equilibrium. There is still a proportion of
systems in nonequilibrium state. $\bar{\tau}_{\rm eq}$ should represent the
relaxation time $\tau_{\rm eq}^{\rm dyn}$ in dynamical equations.

Generally, the relaxation time depends on the mechanism of dynamic process,
the system size, temperature, initial configurations, and so on. In Fig. 2 we
systematically illustrate how the average relaxation time varies with
temperature, the system size and the initial configuration.

Figure 2: Average relaxation time $\bar{\tau}_{\rm eq}$ as a function of
temperature (a) at system sizes $L=50$ (blue circles) and 60 (red squares)
starting from random initial configuration; (b) for random initial
configurations (red squares) and polarized initial configurations (blue
circles) at a fixed system size $L=60$.  Figure 3: The finite size scaling of
average relaxation time at critical temperature. The straight line is a linear
fit.

As Fig. 2(a) shows, in the neighborhood of $T_{\rm c}$, the average relaxation
time has a peak which increases with the system size. Fig. 2(b) also shows
peaks around the critical temperature no matter what initial configuration is.

Due to critical slowing down, the relaxation time near $T_{\rm c}$ for an
infinite system size is expected to diverge as lab33

$\tau_{\rm eq}^{\rm dyn}\propto\xi^{z}\propto|T-T_{\rm c}|^{-z\nu},$ (9)

where $\nu$ is the critical exponent of correlation length. $z$ is the dynamic
exponent, which governs the dynamic universality class 1977 . For a finite
system size, $\xi\sim L$, Eq. (9) reads

$\tau_{\rm eq}^{\rm dyn}\propto L^{z}.$ (10)

This gives the power-law behavior of the relaxation time at $T_{\rm c}$.

In order to test the power-law behavior of $\bar{\tau}_{\rm eq}$, the log-log
plot of $\bar{\tau}_{\rm eq}$ versus system size is presented in Fig. 3. Its
linear region is well fitted by a straight line with slope equal to
2.06$\pm$0.03, consistent with refs. lab17 ; z-ref1 ; z-ref2 .
$\bar{\tau}_{\rm eq}$ of critical temperature indeed diverges as the z-th
power of system size, the same as $\tau_{\rm eq}^{\rm dyn}$.

We also examine the dependence of the standard error and the relative standard
error of $\overline{\tau}_{\rm eq}$ on system size. The standard error
increases with system size, and the relative standard error is almost
constant. This indicates a violation of self-averaging at the critical
temperature self-averaging . The violation of self-averaging appears to be a
common property at criticality.

In Fig. 2(a), for a fixed system size, $\overline{\tau}_{\rm eq}$ shows a
decreasing trend with increasing temperature on both sides of $T_{\rm c}$
except the neighborhood of $T_{\rm c}$. This trend originates from that the
acceptance probability $A({\boldsymbol{u}}\rightarrow{\boldsymbol{v}})$ in Eq.
(7) is an increasing function of temperature $T$. Higher temperature, higher
acceptance probability and then shorter relaxation time.

For $T>T_{\rm c}$, the average relaxation time is weakly dependent on the
system size and the initial configuration, as demonstrated in both Fig. 2(a)
and (b). No matter what the system size is and what the initial configuration
is, the average relaxation time is always short and has a trend towards zero.
This is due to larger acceptance probability for higher temperature. With
larger acceptance probability, it is easier for the system to change from one
configuration to another and hence the relaxation time is extremely short.

In contrary, for $T<T_{\rm c}$, the average relaxation time has strong
dependences on system size and initial configuration. First, Fig. 2(a)
demonstrates that the larger system size, the larger average relaxation time.
It is natural that larger system needs more time to get equilibrium. Second,
Fig. 2(b) demonstrates that the average relaxation time of random
configurations is much longer than polarized configurations. This is because
equilibrium state of $T<T_{\rm c}$ is close to the ordered phase. Polarized
configurations are ordered, and random configurations are disordered. So
polarized initial configuration evolves to an ordered state faster than that
of random initial configuration.

The effect of initial configuration is also demonstrated by the distribution
of $\tau_{\rm eq}$ in the second row (for random initial configuration) and
the third row (for polarized initial configuration) in Fig. 1, respectively.
Wide distributions which extend to large $\tau_{\rm eq}$ result in large
$\bar{\tau}_{\rm eq}$, while narrow distributions which are concentrated at
small $\tau_{\rm eq}$ result in small $\bar{\tau}_{\rm eq}$.

It is natural that the further the initial configuration deviates from the
equilibrium state, the longer the relaxation time is. Therefore, on the left
side of $T_{\rm c}$, random configuration has longer average relaxation time,
while on the right side of $T_{\rm c}$ polarized configuration has longer
average relaxation time, as demonstrated in Fig. 2(b).

In order to observe the maximum influences of nonequilibrium, initial
configurations that have longer nonequilibrium evolution are selected, i.e.
random initial configuration for $T\leqslant T_{\rm c}$ and polarized initial
configuration for $T>T_{\rm c}$ are used in the following.

##  IV Influences of nonequilibrium evolution

In this section the influences of nonequilibrium evolution on observables are
demonstrated in two ways, i.e. how observables vary with evolution time and
with temperature.

It is shown in last section that the distribution of relaxation time varies
with temperature. Even for a fixed temperature and a fixed system size,
relaxation time of two simulations is likely to be different, as the red curve
and the blue curve show in Fig. 1(a). If measurements are made at a time when
not all systems reach thermal equilibrium, there is a certain proportion of
nonequilibrium systems in the sample, very similar to the initial states of
relativistic heavy ion collisions.

High order cumulants of order parameter are sensitive observables in the
search of a critical point. Letting $X=|m|$, $\delta X=|m|-\langle|m|\rangle$,
the first four cumulants read

$\displaystyle C_{1}$ $\displaystyle=$ $\displaystyle\langle X\rangle,$ (11)
$\displaystyle C_{2}$ $\displaystyle=$ $\displaystyle\langle(\delta
X)^{2}\rangle,$ (12) $\displaystyle C_{3}$ $\displaystyle=$
$\displaystyle\langle(\delta X)^{3}\rangle,$ (13) $\displaystyle C_{4}$
$\displaystyle=$ $\displaystyle\langle(\delta X)^{4}\rangle-3\langle(\delta
X)^{2}\rangle^{2}.$ (14)

$C_{1}$ is the mean value of order parameter distribution. $C_{2}$ is the
variance of the distribution. $C_{3}$, $C_{4}$ are related to skewness and
kurtosis, which can quantify non-Gaussianity of the distribution.

Figure 4: The time evolution of cumulants of $|m|$ at a given system size
$L=60$ at $T/T_{\rm c}=0.99$ (the upper row) and 1.01 (the lower row). Random
initial configuration is used at $T/T_{\rm c}=0.99$ and polarized initial
configuration is used at $T/T_{\rm c}=1.01$.

Fig. 4 demonstrates how cumulants vary with the evolution time at two given
temperatures. The evolution starts at $t=0$. At $t=0$, $C_{n=1,2,3,4}$ in Fig.
4(a)-(d) are all zero due to random initial configuration, while $C_{1}$ is 1
and $C_{2,3,4}$ are zero due to polarized initial configuration as Fig.
4(e)-(h) show.

The time evolution of $C_{1}$ in Fig. 4 (a) and (e) demonstrates a similar
behavior with Langevin dynamics, i.e. exponentially approaching the
equilibrium value which is the solution of a linear differential equation
Landau . Therefore, the single-spin-flip mechanism in Metropolis algorithm
equivalently describes the relaxation in Langevin dynamics.

At the temperature on the first order phase transition line, i.e. $T/T_{\rm
c}=0.99$, $C_{1}$ varies monotonously with time until it approaches a steady
value, as Fig. 4(a) shows. The steady value is the equilibrium expectation. In
Fig. 4(b), $C_{2}$ increases first and decreases later, forming a peak during
the evolution. In Fig. 4(c) and (d), both $C_{3}$ and $C_{4}$ experience
oscillations before approaching a steady value. The oscillation results in
sign changes during the evolution. The sign of $C_{3}$ and $C_{4}$ can be
either positive or negative, depending on the evolution time.

At the temperature on the crossover side, i.e. $T/T_{\rm c}=1.01$, the trend
of the nonequilibrium evolution of cumulants is similar to that of $T/T_{\rm
c}=0.99$, as Fig. 4(e)-(h) show. $C_{1}$ varies monotonously with time.
$C_{2}$ experiences nonmonotonic changes before approaching a steady value.
$C_{3}$ also shows sign change during the evolution, similar to Fig. 4(c). The
sign is negative at first and then becomes positive, in contrary to Fig. 4(c).
$C_{4}$ in Fig. 4(h) also shows sign change, being negative first and then
positive, in similar pattern with Fig. 4(d). The sign of $C_{3}$ and $C_{4}$
can be either positive or negative, depending on the evolution time. The sign-
change behavior of $C_{3}$ and $C_{4}$ on the crossover side is consistent
with ref. Mukherjee-2015 ; Nahrgang-2019 .

Apart from the trend, the nonequilibrium evolution of cumulants on both sides
of $T_{\rm c}$ has big differences in two aspects. First, the time needed to
approach equilibrium expectations for $C_{3}$ is about 4000 at $T/T_{\rm
c}=0.99$ and is less than 2000 at $T/T_{\rm c}=1.01$. The difference is more
than twice. That is to say, it is more difficult for systems at $T<T_{\rm c}$
to achieve equilibrium in the same amount of time. Second, the magnitude of
the oscillation in $C_{3}$ and $C_{4}$ at $T<T_{\rm c}$ is much larger than
that at $T>T_{\rm c}$, differing about two orders of magnitude. It means that,
when systems on both sides of $T_{\rm c}$ suffer from nonequilibrium effect,
systems on the low temperature side deviate further from equilibrium
expectations than that on the high temperature side.

Figure 5: Nonequilibrium cumulants of $|m|$ at a given system size $L=60$
measured at $\bar{\tau}_{\rm eq}$ (the upper row) and $3\bar{\tau}_{\rm eq}$
(the lower row), respectively. Random initial configuration is used for
$T\leqslant T_{\rm c}$ and polarized initial configuration is used for
$T>T_{\rm c}$.

Since the sign of $C_{3}$ and $C_{4}$ during the nonequilibrium evolution
depends on the evolution time, two observation times are used as examples to
demonstrate in Fig. 5 how observables vary with temperature for a given
observation time.

In Fig. 5, cumulants in the upper row are measured at $\bar{\tau}_{\rm eq}$
and that in the lower row are measured at $3\bar{\tau}_{\rm eq}$. The
equilibrium cumulants denoted by red circles are obtained after 250,000 Monte
Carlo sweeps, which ensures that all the systems in the sample reach
equilibrium. The blue triangles represent the measured value of cumulants at a
given observation time during the nonequilibrium evolution. We call them
nonequilibrium cumulants.

From the upper row of Fig. 5 we can see that the blue triangles are nearly
coincident with the red circles at $T>T_{\rm c}$ for all $C_{n=1,2,3,4}$. Near
$T_{\rm c}$, the blue triangles start to deviate from the red circles. It is
more clearly seen in the insert map of each figure which focuses on the
neighbourhood of $T_{\rm c}$. The deviations become larger at $T<T_{\rm c}$.
It means that nonequilibrium cumulants observed at $\bar{\tau}_{\rm eq}$ are
not much different from equilibrium cumulants on the high temperature side.
The deviations get larger at low temperature. That is to say, low temperature
system is more largely affected by the nonequilibrium effect.

By comparing different orders of $C_{n}$, we find that the deviation of
nonequilibrium cumulants from equilibrium cumulants is modest at low order,
especially the first cumulant $C_{1}$. Cumulants of higher order have greater
deviations from equilibrium fluctuations.

As time goes by, more systems are expected to reach equilibrium, and the
impact of nonequilibrium should become smaller. The lower row of Fig. 5 which
is measured at a later time $3\bar{\tau}_{\rm eq}$ indeed shows the expected
trend, i.e. more points coincide, including all points at $T>T_{\rm c}$ and
part of points at $T<T_{\rm c}$. Even in the neighborhood of $T_{\rm c}$
almost no deviation is seen as shown in the insert map of each figure.
Deviations only exist for part of points at $T<T_{\rm c}$.

At an observation time of $3\bar{\tau}_{\rm eq}$, nearly all points of $C_{1}$
coincide. Deviations of $C_{2}$ are modest. Deviations increase drastically in
$C_{3}$ and $C_{4}$. This proves that relaxation time varies with observables.
Relaxation times for higher order cumulants are larger, consistent with
expectations Mukherjee-2015 . When low order cumulants have reached their
equilibrium expectations, high order cumulants may have not reached
equilibrium expectations yet.

Observation time of $3\bar{\tau}_{\rm eq}$ is enough for most temperatures.
For $T\lesssim 0.97~{}T_{\rm c}$ blue points still deviate from red points as
shown in Fig. 5(f)(g)(h). This may be caused by the long tail of the
$\tau_{\rm eq}$ distribution at low temperature (see Fig. 1(e)).

##  V Summary and discussion

In this paper we simulate the nonequilibrium evolution of the 3D Ising model
with Metropolis algorithm at zero external magnetic field. The simulation well
produces the dynamical features of nonequilibrium evolution on the phase
boundary.

The relaxation time of each evolution is defined as the number of iterations
needed to achieve equilibrium. It represents the steps that the system needs
to achieve equilibrium. First, the trend of the order parameter approaching a
steady value is consistent with Langevin dynamics. It demonstrates that the
number of iterations in the numerical simulation is equivalent to evolution
time in reality. Second, the average relaxation time at critical temperature
indeed diverges as the theoretical prediction Eq. (10). It indicates that
$\bar{\tau}_{\rm eq}$ exactly corresponds to relaxation time $\tau_{\rm
eq}^{\rm dyn}$ in dynamical equations.

It is found that the average relaxation time depends on the system size and
temperature. In average, the average relaxation time is short at $T>T_{\rm
c}$, almost independent of initial configurations. At $T<T_{\rm c}$, the
average relaxation time is much longer than that of $T>T_{\rm c}$ when the
initial configuration is far from the equilibrium state.

The time evolution of the first four cumulants of order parameter is
presented. For both $T<T_{\rm c}$ and $T>T_{\rm c}$, $C_{3}$ and $C_{4}$ show
oscillations and could be either positive or negative, depending on the
observation time, which is consistent with the results of dynamical models.

Since the relaxation time at $T<T_{\rm c}$ is much longer than that at
$T>T_{\rm c}$, it is easy to approach equilibrium at the crossover side. While
on the line of the first order phase transition, the system is more difficult
to achieve equilibrium and is more susceptible to nonequilibrium. Moreover,
the influence of nonequilibrium on observables on the crossover side is much
weaker than that on the line of the first order phase transition.

Those qualitative features mentioned above imply that the influence of
nonequilibrium evolution in QCD system is most probably negligible at
$T>T_{\rm c}$, i.e. on the crossover side, even if there is a fraction of
nonequilibrium initial states in relativistic heavy ion collisions. While, the
influence should be treated with more caution at $T<T_{\rm c}$, i.e. on the
line of the first order phase transition. Indeed, the sign of the third
cumulant of order parameter may be negative due to the evolution of
nonequilibrium, consistent with the STAR measurements.

By the way, for the nonequilibrium evolution from one equilibrium state to
another equilibrium state, such as a high temperature state is cooled to a low
temperature state, we will study in a coming paper.

## Acknowledgement

We thank Yi Yin, Huichao Song and Shanjin Wu for interesting discussions. This
work is supported in part by the NSFC of China under Grant No. U1732271. The
numerical simulations have been performed on the GPU cluster in the Nuclear
Science Computing Center at Central China Normal University (NSC3).

## References

  * (1) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675 (2006).
  * (2) M. Asakawa, K. Yazaki, Nucl. Phys. A 504, 668 (1989).
  * (3) For a review, see C. Blume, Cent. Eur. J. Phys. 10(6), 1245 (2012) (talk given in CPOD 2011 Wuhan).
  * (4) M.A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
  * (5) S. Ejiri, F. Karsch, K. Redlich, Phys. Lett. B 633, 275 (2006).
  * (6) J. Adam et al. (STAR Collaboration), Phys. Rev. Lett. 126, 092301 (2021).
  * (7) M. A. Stephanov, K. Rajagopal and E. Shuryak, Phys. Rev. D 60, 114028 (1999).
  * (8) M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
  * (9) M. Asakawa, S. Ejiri and M. Kitazawa, Phys. Rev. Lett. 103, 262301 (2009).
  * (10) S. Mukherjee, R. Venugopalan and Y. Yin, Phys. Rev. C 92, 034912 (2015).
  * (11) M. Nahrgang, M. Bluhm, T. Schafer and S. A. Bass, Phys. Rev. D 99, 116015 (2019).
  * (12) S. Wu, Z. Wu and H. Song, Phys. Rev. C 99, 064902 (2019).
  * (13) K. Rajagopal, G. Ridgway, R. Weller and Y. Yin, Phys. Rev. D 102, 094025 (2020).
  * (14) B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000).
  * (15) P. Braun-Munzinger, J. Stachel, Nature 448, 302 (2007).
  * (16) J. Adams et al. (STAR Collaboration), Nucl. Phys. A 757, 102 (2005).
  * (17) R.D. Pisarski, F. Wilczek, Phys. Rev. D 29, 338 (1984).
  * (18) M. A. Stephanov, K. Rajagopal and E. Shuryak, Phys. Rev. Lett. 81, 4816 (1998).
  * (19) M. Asakawa, J. Phys. G: Nucl. Part. Phys. 36 064042 (2009).
  * (20) P. de Forcrand, O. Philipsen, Phys. Rev. Lett. 105 152001 (2010).
  * (21) F. Karsch, E. Laermann and Ch. Schmidt, Phys. Lett. B 520, 41 (2001).
  * (22) B. Friman, F. Karsch, K. Redlich and V. Skokov, Eur. Phys. J. C 71, 1694 (2011).
  * (23) A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, N. Xu, Physics Reports 853, 1 (2020).
  * (24) Cheng-Wei Liu, A. Polkovnikov and A. W. Sandvik, Phys. Rev. B 89, 054307 (2014).
  * (25) M. Acharyya, Phys. Rev. E 56, 2407 (1997).
  * (26) Eq. (118.5) in Statistical Physics Part 1, Third edition, L. D. Landau and E. M. Lifshitz, Pergamon Press, 1980.
  * (27) D. T. Son and M. A. Stephanov, Phys. Rev. D 70, 056001 (2004).
  * (28) L. Onsager, Phys. Rev. 65, 117 (1944).
  * (29) A L Talapov,and H W J. Bl$\ddot{o}$te, J. Phys. A 29,5727(1996).
  * (30) G. S. Pawlcy, R. H. Swendsen,D. J. Wallace and K. G. Wilson, Phys. Rev. B 29, 4030 (1984).
  * (31) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * (32) M E J Newman and G T Barkema. Monte carlo methods in statistical physics[M]. Oxford University Press, 1999.
  * (33) Helmut G. K arXiv:0905.1629(2009).
  * (34) P.C. Hohenberg, B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
  * (35) S. Wansleben and D. P. Landau, J. Appl. Phys. 61, 3968 (1987).
  * (36) M. Hasenbusch, Phys. Rev. E 101, 022126 (2020).
  * (37) A. Malakis, N. G. Fytas, Phys. Rev. E 73, 016109 (2006).