Further author information:
E. D. Miller: E-mail<EMAIL_ADDRESS>
# Understanding the effects of charge diffusion in next-generation soft X-ray
imagers

Eric D. Miller  Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA
02139, USA  Gregory Y. Prigozhin  Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, 77 Massachusetts Ave,
Cambridge, MA 02139, USA  Beverly J. LaMarr  Kavli Institute for Astrophysics
and Space Research, Massachusetts Institute of Technology, 77 Massachusetts
Ave, Cambridge, MA 02139, USA  Marshall W. Bautz  Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology, 77
Massachusetts Ave, Cambridge, MA 02139, USA  Richard F. Foster  Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA  Catherine E. Grant
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA  Craig S. Lage
Department of Physics, UC Davis, One Shields Avenue, Davis, CA 95616, USA
Christopher Leitz  Lincoln Laboratory, Massachusetts Institute of Technology,
244 Wood St., Lexington, MA 02421, USA  Andrew Malonis  Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology, 77
Massachusetts Ave, Cambridge, MA 02139, USA

###### Abstract

To take advantage of high-resolution optics sensitive to a broad energy range,
future X-ray imaging instruments will require thick detectors with small
pixels. This pixel aspect ratio affects spectral response in the soft X-ray
band, vital for many science goals, as charge produced by the photon
interaction near the entrance window diffuses across multiple pixels by the
time it is collected, and is potentially lost below the imposed noise
threshold. In an effort to understand these subtle but significant effects and
inform the design and requirements of future detectors, we present simulations
of charge diffusion using a variety of detector characteristics and
operational settings, assessing spectral response at a range of X-ray
energies. We validate the simulations by comparing the performance to that of
real CCD detectors tested in the lab and deployed in space, spanning a range
of thickness, pixel size, and other characteristics. The simulations show that
while larger pixels, higher bias voltage, and optimal backside passivation
improve performance, reducing the readout noise has a dominant effect in all
cases. We finally show how high-pixel-aspect-ratio devices present challenges
for measuring the backside passivation performance due to the magnitude of
other processes that degrade spectral response, and present a method for
utilizing the simulations to qualitatively assess this performance. Since
compelling science requirements often compete technically with each other
(high spatial resolution, soft X-ray response, hard X-ray response), these
results can be used to find the proper balance for a future high-spatial-
resolution X-ray instrument.

###### keywords:

X-ray detectors, CCDs, charge diffusion, detector response

##  1 Introduction

Future silicon-based soft X-ray imaging instruments will require small pixels
to take advantage of advanced high-resolution optics, while also greatly
benefiting from thick detector bulk to maximize the sensitive band-pass. For
example, the wide-field imaging detectors baselined for both the AXIS Probe-
class mission concept[1] and the Lynx Great Observatory concept[2] require
$\sim$16-$\mu$m pixels to properly sample the $1^{\prime\prime}$ PSF goal, yet
they need to be backside-illuminated devices fully depleted to $\sim$100
$\mu$m to provide sufficient quantum efficiency over the full 0.2–12 keV band.
This pixel aspect ratio of $\sim$6:1 is unprecedented among similar X-ray
instruments, which are closer to 2:1 for Chandra ACIS-S3 and Suzaku XIS1.
“Taller” pixels allow for more lateral diffusion of the X-ray-induced
photoelectrons as they drift toward the collection gates under influence of
the bias electric field, as demonstrated schematically in Figure 1. This
diffusion is greatest for soft photons which interact with the silicon bulk
near the entrance window, and since these photons also produce fewer
photoelectrons than their harder counterparts, charge can be spread out across
many pixels and potentially lost below the imposed noise thresholds, having
substantial effects on the overall soft X-ray response. Here “response” refers
to the combined effects of our ability (or lack thereof) to reconstruct the
energy of individual detected photons, and mainly includes the spectral
resolution FWHM, knowledge of the energy scale, and reduced quantum efficiency
due to events lost below the noise threshold.

It is at the softest energies where some of the most ground-breaking high-
energy astrophysics will be done, including probing the seeds of the first
supermassive black holes at high redshift, tracing the faintest hot gas in the
outskirts of galaxy clusters and the tenuous intergalactic medium, and
understanding the role of stellar and AGN feedback in shaping galaxies like
our own Milky Way[1, 2]iiiSee the Concept Study Reports for Lynx
(https://wwwastro.msfc.nasa.gov/lynx/docs/LynxConceptStudy.pdf) and AXIS
(https://smd-prod.s3.amazonaws.com/science-red/s3fs-
public/atoms/files/AXIS_Study_Rpt.pdf).. For these reasons, understanding the
processes governing the movement of charge in silicon detectors is an area of
active research.[3, 4, 5, 6, 7, 8] In this work, we extend our recent
investigation of charge diffusion[7] in real detectors under development for
strategic missions like AXIS and Lynx[9] by using simulations. These are
validated with real data from different sources, and we use the results to
explore the effects of different detector design parameters on the X-ray
response across a range of energies from 0.3 to 6 keV. We emphasize that these
results are applicable to any silicon-based X-ray imager, including CCDs[9],
CMOS devices[10], and DEPFETs[11].

Figure 1: The effects of thicker detectors with smaller pixels on X-ray event
energy reconstruction. The colored panels show the vertical electric field
strength in two simulated back-illuminated silicon detectors: a 45-$\mu$m
thick device with 24-$\mu$m pixels, similar to Chandra ACIS-S3 or Suzaku XIS1;
and a 100-$\mu$m thick device with 8-$\mu$m pixels, similar to what might be
flown on a future mission. A soft X-ray interacting near the entrance window
has a longer distance to drift in the thicker detector, resulting in more
charge diffusion over a larger number of smaller pixels, as shown in the top-
down view of pixel grids in the lower left.

##  2 Simulations and data analysis

###  2.1 Simulating the detector electric field and charge diffusion

Figure 2: An example run tracking the charge diffusion in Poisson CCD for a
single 5.9 keV Mn K$\alpha$ photon interacting $\sim$6 $\mu$m from the
entrance window in a device with $V_{\mathrm{sub}}=-50$ V. The left panel
shows tracks for individual photoelectrons as they drift and diffuse in the
device’s electric field. For clarity, only 100 of the $\sim$1600 electron
tracks are shown. Note the separation into 8-$\mu$m pixels near $Z=0$. The
right panel shows a top view of the final electron locations, with all
$\sim$1600 electrons represented by blue dots. The orange dot indicates the
projected $X,Y$ location of the initial photon interaction. With 8-$\mu$m
pixels, signal ends up in many pixels, with several containing only a handful
of electrons.

We performed simulations using Poisson CCD [12], a software package that,
given a set of structural and operational characteristics for the detector,
(1) solves Poisson’s equation for the electric field in three dimensions at
all points within the detector volume, and (2) tracks the drift and diffusion
of electrons placed within this volume until they reach the buried channel,
where they are collected in pixels defined by the channel stop and barrier
gate fields. The simulated detector is a backside-illuminated, 100-$\mu$m
thick CCD with 8-$\mu$m pixels, a three-phase gate structure with one
collecting gate held at +1.5 V and two barrier gates held at $-$1.5 V, and
implant structure and doping levels based on the MIT Lincoln Lab CCID-93
device described in our recent comparison of lab data and simulation[7]. The
simulation volume covered 9$\times$9 pixels, with non-linear electric field
grid spacing in the vertical direction allowing finer grid sampling of 30 nm
at the front and back sides of the device to properly capture the electric
field structure, and coarser grid sampling (up to 300 nm) throughout the bulk
of the device. The electric field was simulated at an operating temperature of
$-50^{\circ}$C at substrate bias voltages of $V_{\mathrm{sub}}=-30$, $-50$,
and $-100$ V, all sufficient to fully deplete the silicon bulk and eliminate
field-free regions. Examples of vertical electric field structure resulting
such simulations are shown in Figure 1.

Charge diffusion was simulated for each $V_{\mathrm{sub}}$ using 100,000
photons each at energies of 5.89 keV (Mn K), 1.25 keV (Mg K), 0.525 keV (O K),
and 0.277 keV (C K). This was done by introducing small clouds of electrons
into the simulated detector volume and allowing them to drift and diffuse
under the influence of the 3-D electric field until collected in pixels. The
number of electrons in each cloud was drawn from a Fano noise distribution
appropriate for the photon energy in Si (e.g., 1615$\pm$14 electrons for 5.9
keV), and the interaction depth was drawn from an exponential distribution
with appropriate energy-dependent attenuation length. Each interaction
location was generated from a uniform distribution across the central pixel of
the 9$\times$9-pixel simulation volume. An example diffusion run is shown in
Figure 2.

We simulated the quality of backside passivation by adjusting the
TopAbsorptionProb parameter, which is the probability that a diffused
photoelectron encountering the backside entrance window is absorbed rather
than reflected back into the bulk; this ranged from 0% (perfect passivation)
to 25%. The final pixelized 2-D spatial distribution of electrons was binned
2$\times$2 and 3$\times$3 to simulate 16-$\mu$m and 24-$\mu$m pixels,
respectively; although this is technically different than simulating drift and
diffusion in a detector with the correct pixel size, it greatly reduced the
required simulation time, and several test runs showed there was no noticeable
difference between the methods, as the majority of the diffusion occurs in the
bulk before the charge cloud feels the effects of the gate barrier and channel
stop electric fields (see Figures 1 and 2, left panel). Finally, Gaussian
readout noise ranging from 1 to 4 e- RMS was added to each pixel, and event
detection and characterization were performed in the same way as flight
software on Chandra ACIS and Suzaku XIS, as described in our previous work.[7]
For each single-photon frame, we used a 5 $\sigma$ threshold for event
candidate identification, and summed all pixels including neighbors above 4
$\sigma$ to calculate the total event energy. These thresholds are similar to
values tuned for previous instruments to eliminate noise events and inclusion
of noisy pixels in event pulse height summation. At a temperature of
$-50^{\circ}$C, we expect a conversion of 3.71 eV per e-[13], and so in the
highest noise case of 4 e- RMS, the event threshold corresponds to 74 eV and
the neighbor (or “split”) threshold to 59 eV. These are high enough to lose a
significant fraction of the signal in a low-energy event split over several
pixels. In addition, a 5 $\sigma$ event threshold has a 25% chance of
introducing a noise event in a 1-Mpix array; for the expected $\sim$20 fps
frame rate for a mission like AXIS, this will introduce many such events, and
thus an event threshold that is many factors of the readout noise and still
well below the science band of interest is desirable.

As also discussed in our previous work[7], the default charge diffusion
parameters in Poisson CCD produced much more diffusion in the simulations than
we observe in the lab data, with much larger event sizes. We used similar
DiffMultiplier parameter values that produced simulated event size
distributions closely matched to the real data, and also similar to values
supported by theoretical studies. The validation with real data presented in
our previous work and here in Section 2.2 support use of this diffusion
factor, and in fact more diffusion would magnify the importance of factors
such as readout noise, pixel size, and bias voltage in the resultant soft
response. The slight discrepancy between our tuned DiffMultiplier value and
the expected value will be explored in a future paper.

Figure 3: Validation of the charge diffusion simulations with real data. The
top panel compares lab data using a very similar but thinner detector to that
used in our other simulations here. As described in our previous work[7], we
fit a 2-D Gaussian to the pixel islands of the measured and simulated events
at a range of substrate bias voltage $V_{\mathrm{s}ub}$, and find that the
distribution of Gaussian $\sigma$ (a proxy for lateral diffusion magnitude) is
similar. The bottom panel shows the fraction of events spanning 1, 2, or 3–4
pixels in real Suzaku XIS1 flight data compared to simulations of that same
device. The fractions match remarkably well at these soft energies, suggesting
that the charge diffusion is accurately simulated.

###  2.2 Validation of simulations with real data

To validate the simulations, we simulated charge diffusion in two real X-ray
detectors for which we have substantial data, both BI CCDs fabricated by MIT
Lincoln Laboratory: the CCID-93 small-format, small-pixel (8-$\mu$m)
device[7]; and the CCID-41 device flown as Suzaku XIS1[14]. The CCID-93
simulations are similar to the unbinned 5.9-keV runs described in Section 2.1,
although with 50-$\mu$m thickness, and as in our previous work[7], we fit 2-D
Gaussians to pixelated event islands to measure the charge diffusion. As shown
in the top panel of Figure 3, the simulated sizes match the lab data very
well. The Suzaku data includes flight data from a number of bright, un-piled
sources, including the Perseus Cluster and the bright blazars PKS 2155 and 3C
273, the latter taken in 1/8 window mode. This data from these bright
continuum sources is dominated by X-ray photons across a wide energy band,
with very little contamination from particle events. We parametrize charge
diffusion using the ASCA grade specifying how many pixels in an event are
above the noise threshold. As shown in the bottom panel of Figure 3, the
simulations match the observed charge diffusion exceptionally well at soft
energies. While both real devices are $\sim$50 $\mu$m thick, about half that
of the notional simulated detector, the dependence of diffusion with drift
distance is well-defined, and this validation gives us confidence that our
simulation results are reliable.

##  3 Results

###  3.1 Spectral response

The effects of charge diffusion on spectral response are apparent in all of
our simulations. Figures 4 and 5 show histograms of reconstructed event energy
for the four monochromatic energies under study: C K (0.28 keV), O K (0.53
keV), Mg K (1.25 keV), and Mn K (5.9 keV). Colors show results for different
values of readout noise, ranging from 1 to 4 e- RMS. In each figure, an
“ideal” detector from the standpoint of maximizing charge collection is shown
in the top panels. A lower performance but still realistic detector is
represented in the lower panels. Pixel size, substrate bias voltage, and back
surface passivation quality contribute to the width and shape of the response.
These effects are greatest at low energies, where peak energy shifts and non-
Gaussian shapes appear even for low readout noise. Quantum efficiency (QE) is
affected as well, in the worst case resulting in no events being detected for
C K with 4 e- RMS noise (Figure 4, lower left).

Figure 4: Spectral response of the simulated detector to monochromatic photons
of (left) C K at 0.28 keV and (right) O K at 0.53 keV. Top panels show high-
quality charge collection with large pixels, high bias voltage, and perfect
backside passivation. Bottom panels show less-than-ideal charge collection,
with small pixels, low bias voltage, and imperfect passivation, diminished
performance that may be imposed by other instrument constraints. The event
energies are calculated by summing pixels above threshold, as described in the
text. Higher readout noise increases the spectral FWHM in all cases, and in
the less-than-ideal case causes the histogram to shift to much lower energies
and all events to be lost in the worst case at C K. Figure 5: Same as Figure
4, but for (left) Mg K at 1.25 keV and (right) Mn K at 5.9 keV. The effects of
the less-than-ideal detector and higher readout noise are moderate at these
energies, although non-Gaussian features are obvious in the higher noise
spectra.

Each spectral histogram in Figure 4 and 5 was fit with a simple Gaussian to
estimate the spectral resolution FWHM and the peak energy shift. These are
shown in Figure 6, with top and bottom panels showing different pixel sizes,
and substrate bias voltage increasing from right to left. In effect, the upper
left panel is the lower performance detector, and the bottom right panel is
the ideal. Colored bands show the range of FWHM and peak shift for a range of
“good” backside passivation. In all cases, readout noise of 2–3 e- or lower is
necessary to meet the spectral resolution requirements of a mission like AXIS,
shown by the gray band. Even at this level, there is little allowance for
additional FWHM degradation that may be caused by (for example) poorer
backside treatment or charge-transfer inefficiency (CTI) in CCDs[15].

Figure 6: Spectral FWHM and peak shift as a function of readout noise for our
simulated detector with different pixel size (top: 16 $\mu$m, bottom: 24
$\mu$m) and bias voltage ($V_{\mathrm{sub}}=-30$ to $-100$ from left to
right). Colored bands span a range of high passivation quality, 0–10%
probability of a photoelectron being absorbed if it diffuses to the entrance
surface. Theoretical limits are shown by dotted lines, and assume the charge
is collected in a single pixel with the specified readout noise. Nominal
requirements for a mission like AXIS are shown by the gray band. All notional
detectors require noise $<$3 e- to reach the required FWHM at 0.5 keV; for
smaller pixels and lower bias voltage, the noise limits are even more strict.
Poorer backside performance will degrade the response further. Performance of
the the Suzaku XIS1 CCD is shown as stars at the measured readout noise of 2.5
e-. These CCDs benefit from a small pixel aspect ratio and excellent backside-
passivation[16], both enabling significantly better charge collection than the
thick, small-pixel devices simulated here.

###  3.2 Assessing the back surface quality

To maximize QE at X-ray photon energies where the attenuation length in
silicon is very short ($<$0.5 $\mu$m below 0.5 keV), X-ray imagers must be
thinned, fully depleted, and illuminated on the backside. The thinning process
introduces open bonds on the surface that can absorb photoelectrons, and for
this and other reasons the surface must be passivated by introducing a thin
layer of dopant. The CCID-93 device currently under study uses molecular beam
epitaxy (MBE) to introduce this layer[17], and its effectiveness at
passivating the surface is of great interest as the science goals of future
strategic missions demand good soft X-ray QE and spectral resolution.

As we have shown, a large number of factors contribute to the overall detector
response, and these can be difficult to separate from each other in lab data.
The simulations we have presented here provide a method for doing just that,
since we can easily change the magnitude of these various factors, and since
they have been validated for a real-life instrument with exquisite soft
response (see Figure 3). To demonstrate this, we have used the 1 e- RMS
readout noise, 16-$\mu$m pixel, $V_{\mathrm{s}ub}=-50$ V simulations in an
attempt to quantify the fraction of photons that suffer charge loss at the
entrance surface. As described in Section 2.1, the backside quality is
captured in Poisson CCD in a parameter called TopAbsorptionProb, the
probability that a photoelectron encountering the entrance window will be
absorbed and thus lost from the charge packet created by the X-ray. We show in
Figure 7 monoenergetic spectra for these simulations at three energies
(columns) and for four values of the backside quality (rows). Perfect
passivation (TopAbsorptionProb = 0%) is shown at the top, and we fit each of
these spectra with a single Gaussian, which is a good fit. For the poorer
passivation simulations (TopAbsorptionProb $>$ 0%), we introduce two
Gaussians; the first, shown in blue, represents events that have no surface
losses, and in these fits the $\sigma$ and mean are fixed to the perfect
passivation Gaussian from the top panel. The second, orange Gaussian
represents the ensemble of events that have lost some fraction of charge to
the surface, and this is allowed to fit freely. The blue Gaussian amplitude is
allowed to vary to account for events that migrate to this ensemble.

Figure 7: Demonstration of a method to assess backside passivation quality by
separating it from other response-broadening effects. The black curve in each
panel shows the response to monochromatic X-rays from our simulations, with
different energies shown in columns and different backside quality shown in
rows. Here backside “quality” is parameterized by the TopAbsorptionProb
probability that a photoelectron encountering the entrance window will be
absorbed. We use the simulations with 1 e- RMS readout noise, 16-$\mu$m
pixels, and $V_{\mathrm{s}ub}=-50$ V. The top panels represent perfect
passivation, and we fit these histograms with a single Gaussian. The lower
panels represent increasingly poor passivation, which we model with two
Gaussians as explained in the text. As expected, the fraction of events
affected by surface losses (represented by the orange Gaussian) depends
strongly on energy.

At the lowest energy of 0.28 keV, where the attenuation length is only 0.11
$\mu$m, even a mild decrease in passivation quality results in a majority of
photons losing some charge to the surface: 64% for 5% absorption probability,
and 93% for 25%. At 0.53 keV and 0.49 $\mu$m attenuation length, the non-
Gaussianity of the response is easier to see, and about half of the photons
lose charge to the surface in the 25% absorption probably case. At 1.25 keV,
the attenuation length is $\sim$5 $\mu$m, and while a low-energy tail can be
seen, only about 6% of photons are contained in the orange Gaussian and thus
identified as losing charge to the surface. This fraction is the same for all
absorption probabilities, and may indicate something physical or a limitation
of this technique for photon ensembles which have a range of penetration
depths extending well into the silicon bulk. The results for 5.9 keV are not
shown, as with a 29-$\mu$m attenuation length there is little surface loss;
the vast majority of photons interact away from the surface and have their
charge quickly caught up in the bias-induced drift.

We have used the 1 $e^{-}$ readout noise simulations in this test case because
it is easiest to see the effects of backside quality when it dominates the
response broadening compared to other factors. In future work, we will expand
this method to more demanding cases with higher noise and other factors that
degrade resolution, such as CTI. We will also explore a more physical and
quantitative method of assessing the backside quality and other surface
effects than the simplified TopAbsorptionProb parameter. In particular, we
don’t expect real events with lost charge to follow a Gaussian distribution,
but rather have an extended non-Gaussian or even constant tail to lower
energies. The Gaussianity is likely an artifact of this simplified
probabilistic charge loss that ignores (for example) initial high-energy
photoelectrons completely leaving the silicon. We will also probe the
dependence at different energies and as a function of interaction depth of
individual photons rather than as an ensemble.

##  4 Summary

We have simulated charge diffusion in thick, small-pixel, back-illuminated
silicon X-ray imagers to understand the effects on soft X-ray response. The
simulations show that, while larger pixels, higher bias voltage, and optimal
backside passivation improve performance, reducing the readout noise has a
dominant effect in all cases. Since compelling science requirements often
compete technically with each other (high spatial resolution, soft X-ray
response, hard X-ray response), these results can be used to find the proper
balance for a future high-spatial-resolution X-ray instrument.

This variety of response-degrading effects in such detectors complicates
assessment of the backside passivation quality, vital for the good soft X-ray
response demanded by the science goals of future strategic missions. We have
presented a method to qualitatively assess the backside treatment, using our
validated simulations to estimate the relative size of the ensembles of
photons that lose charge to the surface and those that do not. We will explore
this more in future work, with the aim of developing a spectral resolution
budget that can be adapted to help formulate future mission requirements.

###### Acknowledgements.

We gratefully acknowledge support from NASA through the Strategic Astrophysics
Technology (SAT) program, grants 80NSSC18K0138 and 80NSSC19K0401 to MIT, and
from the Kavli Research Infrastructure Fund of the MIT Kavli Institute for
Astrophysics and Space Research.

## References

  * [1] Mushotzky, R., “AXIS: a probe class next generation high angular resolution x-ray imaging satellite,” in [Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray ], Proc. SPIE 10699 (August 2018).
  * [2] Gaskin, J. A., Allured, R., Bandler, S. R., Basso, S., Bautz, M. W., Baysinger, M. F., Biskach, M. P., Boswell, T. M., Capizzo, P. D., Chan, K.-W., Civitani, M. M., Cohen, L. M., Cotroneo, V., Davis, J. M., DeRoo, C. T., DiPirro, M. J., Dominguez, A., Fabisinski, L. L., Falcone, A. D., Figueroa-Feliciano, E., Garcia, J. C., Gelmis, K. E., Heilmann, R. K., Hopkins, R. C., Jackson, T., Kilaru, K., Kraft, R. P., Liu, T., McClelland, R. S., McEntaffer, R. L., McCarley, K. S., Mulqueen, J. A., Özel, F., Pareschi, G., Reid, P. B., Riveros, R. E., Rodriguez, M. A., Rowe, J. W., Saha, T. T., Schattenburg, M. L., Schnell, A. R., Schwartz, D. A., Solly, P. M., Suggs, R. M., Sutherlin, S. G., Swartz, D. A., Trolier-McKinstry, S., Tutt, J. H., Vikhlinin, A., Walker, J., Yoon, W., and Zhang, W. W., “Lynx Mission concept status,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10397, 103970S (Sept. 2017).
  * [3] Prigozhin, G., Butler, N. R., Kissel, S. E., and Ricker, G. R., “An experimental study of charge diffusion in the undepleted silicon of X-ray CCDs,” IEEE Transactions on Electron Devices 50, 246–253 (Jan. 2003).
  * [4] Miller, E. D., Foster, R., Lage, C., Prigozhin, G., Bautz, M., Grant, C., LaMarr, B., and Malonis, A., “The effects of charge diffusion on soft x-ray response for future high-resolution imagers,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10699, 106995R (July 2018).
  * [5] Haro, M. S., Fernandez Moroni, G., and Tiffenberg, J., “Studies on small charge packet transport in high-resistivity fully depleted ccds,” IEEE Transactions on Electron Devices 67(5), 1993–2000 (2020).
  * [6] Prigozhin, I., Dominici, S., and Bellotti, E., “Fbmc3d—a large-scale 3-d monte carlo simulation tool for modern electronic devices,” IEEE Transactions on Electron Devices 68(1), 279–287 (2021).
  * [7] LaMarr, B. J., Prigozhin, G. Y., Miller, E. D., Thayer, C., Bautz, M. W., Foster, R., Grant, C. E., Malonis, A., Burke, B. E., Cooper, M., Donlon, K., and Leitz, C., “Measurement and simulation of charge diffusion in a small-pixel charge-coupled device,” Journal of Astronomical Telescopes, Instruments, and Systems 8, 016004 (Jan. 2022).
  * [8] Avalos, N. E. and Sofo Haro, M., “G-CoReCCD: a GPU-based simulator of the charge transport in fully-depleted CCDs,” in [SPIE Future Sensing Technologies ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 11525, 115250I (Nov. 2020).
  * [9] Bautz, M. W., Foster, R., Grant, C. E., LaMarr, B., Malonis, A., Miller, E. D., Prigozhin, G., Burke, B., Cooper, M., Donlon, K., Lambert, R., Warner, K., Young, D., Chattopadhyay, T., Herrmann, S., Morris, R. G., Leitz, C., and Allen, S., “Performance of high frame-rate CCDs for future strategic missions,” in [Space Telescopes and Instrumentation 2022: Ultraviolet to Gamma Ray ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 12181, 12181–85 (2022).
  * [10] Chattopadhyay, T., Falcone, A. D., Burrows, D. N., Hull, S., Bray, E., Wages, M., McQuaide, M., Buntic, L., Crum, R., O’Dell, J., and Anderson, T., “X-ray hybrid CMOS detectors: recent development and characterization progress,” in [Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10699, 106992E (July 2018).
  * [11] Meidinger, N., Barbera, M., Emberger, V., Fürmetz, M., Manhart, M., Müller-Seidlitz, J., Nandra, K., Plattner, M., Rau, A., and Treberspurg, W., “The Wide Field Imager instrument for Athena,” Proc. SPIE 10397, 103970V (Aug. 2017).
  * [12] Lage, C., Bradshaw, A., Anthony Tyson, J., and LSST Dark Energy Science Collaboration, “Poisson_CCD: A dedicated simulator for modeling CCDs,” Journal of Applied Physics 130, 164502 (Oct. 2021).
  * [13] Groom, D. E., Bebek, C. J., Fabricius, M., Karcher, A., Kolbe, W. F., Roe, N. A., and Steckert, J., “Quantum efficiency characterization of LBNL CCD’s: Part I. The quantum efficiency machine,” in [Sensors, Cameras, and Systems for Scientific/Industrial Applications VII ], Blouke, M. M., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 6068, 133–143 (Feb. 2006).
  * [14] Koyama, K., Tsunemi, H., Dotani, T., Bautz, M. W., Hayashida, K., Tsuru, T. G., Matsumoto, H., Ogawara, Y., Ricker, G. R., Doty, J., Kissel, S. E., Foster, R., Nakajima, H., Yamaguchi, H., Mori, H., Sakano, M., Hamaguchi, K., Nishiuchi, M., Miyata, E., Torii, K., Namiki, M., Katsuda, S., Matsuura, D., Miyauchi, T., Anabuki, N., Tawa, N., Ozaki, M., Murakami, H., Maeda, Y., Ichikawa, Y., Prigozhin, G. Y., Boughan, E. A., Lamarr, B., Miller, E. D., Burke, B. E., Gregory, J. A., Pillsbury, A., Bamba, A., Hiraga, J. S., Senda, A., Katayama, H., Kitamoto, S., Tsujimoto, M., Kohmura, T., Tsuboi, Y., and Awaki, H., “X-Ray Imaging Spectrometer (XIS) on Board Suzaku,” PASJ 59, 23–33 (Jan. 2007).
  * [15] LaMarr, B., Bautz, M., Miller, E., Prigozhin, G., Foster, R., Cooper, M., and Donlon, K., “Charge transfer efficiency of CCDs with a single-polysilicon gate structure,” in [X-Ray, Optical, and Infrared Detectors for Astronomy X ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 12191, 12191–89 (2022).
  * [16] Bautz, M. W., Kissel, S. E., Prigozhin, G. Y., LaMarr, B., Burke, B. E., and Gregory, J. A., “Progress in x-ray CCD sensor performance for the Astro-E2 X-ray imaging spectrometer,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], A. D. Holland, ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 5501, 111–122 (Sept. 2004).
  * [17] Ryu, K. K., Leitz, C. W., Clark, H. R., Chen, X., Cooper, M. J., Zhu, M., Welander, P. B., Lambert, R. D., Bolkhovsky, V., Yost, D. R. W., Burke, B. E., Gregory, J. A., and Suntharalingam, V., “Oxide-bonded molecular-beam epitaxial backside passivation process for large-format CCDs,” in [Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray ], den Herder, J.-W. A., Nikzad, S., and Nakazawa, K., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10699, 106993P (July 2018).


# Dynamics of vortex defect formation in two dimensional Coulomb crystals

Michael Arnold<EMAIL_ADDRESS>Centre for Complex Systems,
Faculty of Engineering, The University of Sydney, Sydney, NSW 2006, Australia.
Ramil Nigmatullin<EMAIL_ADDRESS>Center for Engineered Quantum
Systems, Dept. of Physics & Astronomy, Macquarie University, 2109 NSW,
Australia.

###### Abstract

We study the non-equilibrium dynamics of two dimensional planar ion Coulomb
crystals undergoing a structural buckling transition to a three plane
configuration, driven by a reduction of the transverse confining frequency.
This phase transition can be theoretically modeled using a mapping to a two
dimensional Ginzburg-Landau theory with complex order parameter field. We
demonstrate that finite rate quenches result in creation of stable topological
vortices, which are localized point regions around which the phase of the
order parameter field winds a multiple of $2\pi$. The density of the defects
as a function of quench rate is investigated using molecular dynamics
simulations and its scaling is shown to be consistent with Kibble-Zurek theory
of defect formation. Following the quench, the annihilation of vortex and
anti-vortex pairs results in the relaxation of defect density that follows a
diffusive scaling with a logarithmic correction. This work highlights the
potential for investigating complex non-equilibrium statistical physics of
topological defects in an experimentally accessible ion trap setting.

##  I Introduction

Trapped ions is one of the most prominent quantum technologies. Doppler laser
cooling can bring ions to low temperatures at which the ions crystallize
forming regular structures, known as Coulomb crystals, whose shape is
determined by the trapping parameters. Linear chain crystals in Paul traps
have the simplest phonon spectrum and have been widely used in metrology [1],
quantum computing [2, 3] and quantum simulations [4]. In Penning traps, large
two dimensional planar crystals can be readily created and these structures
have also been used in metrology and quantum information processing
applications [5].

Beyond these simple Coulomb crystal geometries, there is great interest in
exploring more complex structural phases and the transitions between them [6,
7, 8, 9, 10]. The strong long-range Coulomb interactions between particles
leads to highly non-linear non-trivial dynamics, whose investigation is of
fundamental interest in the fields of nonlinear science, complex systems and
solid-state physics; it is useful as a platform for studying complex non-
linear and non-equilibrium dynamics in areas including the simulation of
Klein-Gordon fields on a lattice [11], Kibble-Zurek mechanism of defect
formation [12, 13, 14], dynamics of discrete solitons [15, 16, 17], dry
friction [18], energy transport [19, 20] and synchronization [21]. Coulomb
crystals with more complexity also provide new lattice geometries for quantum
simulations [22].

In this paper, we investigate numerically the non-equilibrium structural phase
transition from a quasi-two-dimensional 1-plane crystal to a 3-plane crystal.
We focus on the Kibble-Zurek mechanism of the formation of topological defects
and the subsequent coarsening dynamics of annihilation of defects and anti-
defect pairs. Previous studies of KZ mechanism in ion crystal systems have
focused on the linear to zigzag phase transition in a quasi-one-dimensional
system [12, 13]. This are symmetry breaking phase transitions and the
resulting defects are either kinks [12, 13] if the $Z_{2}$ symmetry is broken
or helical twists [14] if the $U(1)$ symmetry is broken. The 1-plane to
3-plane structural phase transition considered in this paper can be described
by an XY 6-clock model with an intermediate Kosterlitz-Thouless phase. The
defects are $U(1)$ point vortices whose physics is considerably richer [23].

Previously, simulations of finite quenches in a two dimensional XY spin model
have shown that the density of vortices is dictated by both the KZ mechanism
and coarsening dynamics of annihilation of vortex/anti-vortex pairs [24]. This
observation was corroborated in an experimental study of colloids undergoing a
phase transition via Kosterlitz– Thouless–Halperin–Nelson–Young (KTHNY)
melting scenario [25]. The molecular dynamics simulations of microscopic ion
crystal dynamics presented in this paper are consistent with the dynamics
observed in [24, 25]. An observed deviation from the theoretical coarsening
scaling may be attributed to the pinning effect of the emergent 6-clock
potential. Our work introduces a new platform for investigating the non-
equilibrium KT phase transition and more generally the collective dynamics of
interacting vortices.

The paper is organized as follows. Section II presents the microscopic model,
reviews its mapping to the Ginzburg-Landau field theory and the phase diagram
of a 1-plane to 3-plane structural transition. Section III uses molecular
dynamics simulations to demonstrate the existence of topological defects in
the 3-plane phase and presents an algorithm for determining their location.
Section IV focuses on molecular dynamics investigation of finite rate
quenches, where the KZ mechanism and coarsening dictate the evolution of the
average number of defects in the system.

##  II Ginzburg-Landau model of 1- to 3-plane transition

Ion traps confine repulsively interacting ions in space either by rapidly
varying oscillatory electric fields, as in the Paul traps, or a combination of
static electric and magnetic fields, as in the Penning traps [26]. The
dynamics of ion Coulomb crystal can be often approximated with the so-called
pondermotive approximation or pseudopotential theory, which replaces the time-
varying trapping fields experienced by particles by an effective time-
independent harmonic potential. The laser cooling reduces the temperature of
the ions such that they can form regular crystal-like configurations, whose
overall shape is determined by the trap parameters. We will consider a system
of $N$ ions confined to a periodic cell in the $x$-$y$ plane and by the
harmonic confinement in the $z$-direction. The potential energy is given by

$V=\frac{1}{2}m\omega_{z}\sum_{j}^{N}z_{j}^{2}+\mathcal{K}\sum_{i<j}\frac{1}{|\textbf{r}_{i}-\textbf{r}_{j}|}$
(1)

where $\textbf{r}_{j}=(x_{j},y_{j},z_{j})$ are the coordinates of the $j$th
ion, $\mathcal{K}\equiv q^{2}/4\pi\epsilon_{0}$, $q$ is the charge of the ion,
$\epsilon_{0}$ is the vacuum permittivity, $m$ is the mass of the ion and
$\omega_{z}$ is the trapping frequency in the $z$-direction. The periodic
boundary conditions results in a homogeneous spacing in the ion crystal. In a
real experimental system, the open boundary conditions and the harmonic
confinement in the $x$ and $y$ direction would result in inhomogeneous spacing
between the ions; the ions are closer together in the centre of the crystal
and further apart near the edges. The system with periodic boundary condition
can be viewed as an approximation to a central region of a large ion crystal,
where the spacing is approximately homogeneous and the boundary effects can be
neglected.

Above a certain critical value of $\omega_{z}=\omega_{z}^{(c)}$ the lowest
energy configuration is a planar triangular lattice. When the confining
frequency is reduced to below $\omega_{z}^{(c)}$, the 1-planar crystal
configuration undergoes a buckling structural transition into a 3-planes, all
of which in triangular lattice geometry but with double the lattice spacing
(see Figure 1). This buckling instability has been predicted in an early
theoretical work by Dubin [27] and has been observed experimentally in [28].
Recently, Podolsky _et. al._ [23] derived a Ginzburg-Landau (GL) field theory
for this transition thereby proving that it is in the universality class of a
two dimensional XY model [23].

The GL field theory is derived by Taylor expanding the non-linear Coulomb
interaction term in equation (1) in displacements around the equilibrium
lattice positions. In [23] it was shown that one must keep the terms up to the
sixth order in the expansion to correctly capture the critical properties of
the structural phase transition. The GL free energy density is give by

$\frac{f}{\mathcal{K}}=\frac{\gamma}{2}|\nabla\psi|^{2}+\epsilon|\psi|^{2}+u|\psi|^{4}+v|\psi|^{6}+\frac{w}{2}\left[\psi^{6}+(\psi^{*})^{6}\right]$
(2)

where
$\epsilon=\frac{1}{\sqrt{3}}\left(\frac{m\omega_{z}^{2}a^{2}}{2\mathcal{K}}-I_{2}\right)$,
$u=3/\sqrt{3}/4I_{4}$, $w=\frac{5}{8\sqrt{3}}I_{6}$,
$v=-\frac{25}{4\sqrt{3}}I_{6}$, $I_{2}=6.683$, $I_{4}=3.56$, $\gamma=0.223$,
$I_{6}=2.558$ and $a$ is the lattice spacing. The order parameter field at a
lattice point with coordinates $(x_{j},y_{j}$, $\psi(x_{j},y_{j})$, is an
implicit function of the transverse displacement

$z_{j}=\textrm{Re}\left[\psi e^{i\textbf{K}\cdot\textbf{r}_{j}}\right].$ (3)

Here K is the base vectors of the first Brillouin zone of the triangular
lattice given by $\textbf{K}=(4\pi/3,0)$,
$\textbf{r}_{i}=n_{1}\textbf{a}_{1}+n_{2}\textbf{a}_{2}$, with
$\textbf{r}_{1,2}=(\frac{1}{2},\pm\sqrt{3}/2)$. The order parameters $\psi$ is
complex and can be expressed as $\psi=|\psi|e^{i\theta}$.

Figure 1: Structural buckling transition between a 1-plane and 3-plane
configuration. (a) The transverse displacement of the ions in a crystal at
different values of $\omega_{z}$ in the vicinity of the critical
$\omega_{z}^{(c)}$. (b) The triangular lattice structure of the 3-plane phase.
Figure 2: (a) Single well potential of GL theory for $\epsilon>0$ (b) Mexican
hat potential for $\epsilon<0$.

The potential energy for the mean field configuration, which neglects the
spatial fluctuations in the order parameter,
$V(\psi)/\mathcal{K}=\epsilon|\psi|^{2}+u|\psi|^{4}+v|\psi|^{6}+\frac{w}{2}\left[\psi^{6}+(\psi^{*})^{6}\right]$,
is shown in Figure 2. For $\epsilon>0$ the order parameter is zero, $\psi=0$,
and the system is in the 1-plane phase. The potential is a single well and
since the order parameter has no preferred direction, the 1-plane phase is
disordered. For $\epsilon<0$ the order parameter is non-zero and the system is
in the 3-plane phase. The potential is a Mexican hat but with 6 equally spaced
wells which correspond to the local order of the 6 degenerate lattice
arrangements shown in Figure 2(b). This corresponds to the 6-clock phase which
has the discrete $Z_{2}\times Z_{3}$ symmetry. At higher energies the order
parameter can easily overcome the energy barrier between the neighboring
minima, and the symmetry changes to the broken $U(1)$ continuous symmetry. Two
dimensional systems with the broken $U(1)$ symmetry in the order parameters
support topological defect vortex configuration, which are localized regions
where the field winds around the Mexican hat potential. The presence of these
topological defects drastically alters the physics of the system leading to
the existence of the KT phase [29, 30]. Thus, near the critical point of the
1-plane to 3-plane structural phase transition the system can exist in 3
phases, disordered, KT and the 6-face clock phases, depending on the value of
$\epsilon$ and temperature $T$. The phase diagram was derived in [23] and is
sketched in Figure 3. The KT phase is characterized by a change in behavior in
correlation length. For $T>T_{KT}$, the system is disordered, the correlation
length decays exponentially and there is a finite density of unbound vortices.
For $T<T_{KT}$, there is a quasi-long range order with a power law decays of
the correlation length and vortices and anti-vortices form bound pairs. For
$T<T_{6}$, the system is in the 6-clock phase, which again exhibits a long-
range order with exponentially decaying correlation length.

Figure 3: Phase diagram of a structural 1-plane to 3-plane phase transition of
a Coulomb crystals. The control parameter space determined by the temperature
$T$ and coefficient $\epsilon$ which is a function of the ratio
$m\omega_{z}^{2}a^{2}/\mathcal{K}$. The line $T_{KT}$ indicates a transition
between 1-plane phase and a 3-plane phase in quasi-long-range ordered KT
phase. The line $T_{6}$ indicates a transition between the KT phase and the
long range ordered 3-plane phase with the $Z_{6}$ symmetry.  Figure 4: a)
Example of a configuration with topological defects in the 3-plane structural
phase represented as color plot (top) or arrow plot (bottom). The observed
defects are asters and vortices (topological charge, $s=+1$) or cross-hairs
and anti-vortices (topological charge, $s=-1$) (b) Instantaneous
configurations of a quasi-two-dimensional ion crystal with 9858
$\textrm{Be}^{+}$ ions at a temperature of $20\mu$K and $\gamma=8.96\times
10^{-7}\textrm{ps}^{-1}$ in the 3-plane structural phase at the end of a non-
equilibrium quench from a disordered single plane phase at different rate
$\tau_{Q}^{-1}$. The quench is implemented by reducing the transverse
confining frequency $\omega_{z}$. The density of topological defects increases
with increasing quench rate. Each color corresponds to a phase angle $\Theta$
of the local order parameter $\psi=|\psi|e^{i\Theta}$.

##  III Topological defects in the 3-plane phase

To verify the prediction of the existence of the topological defects, we have
performed molecular dynamics simulations of ion Coulomb crystals confined in a
box with periodic boundary conditions in the $x$ and $y$ direction. The
topological defects are produced by quenching a system from a 1-plane
disordered phase into the 3-plane phase. We use a Langevin thermostat to
simulate the interaction of the ions with the cooling laser beam, which
thermalized the ions. The equations of motion for the $j$th ion are given by

$\displaystyle m\partial_{tt}x_{j}$ $\displaystyle=$
$\displaystyle-m\gamma\partial_{t}x_{j}-\partial_{x_{j}}V_{c}+\theta_{xj}(t)$
(4) $\displaystyle m\partial_{tt}y_{j}$ $\displaystyle=$
$\displaystyle-m\gamma\partial_{t}y_{j}-\partial_{y_{j}}V_{c}+\theta_{yj}(t)$
(5) $\displaystyle m\partial_{tt}z_{j}$ $\displaystyle=$
$\displaystyle-m\omega(t)^{2}-m\gamma\partial_{t}z_{j}-\partial_{z_{j}}V_{c}+\theta_{zj}(t),$
(6)

where $m$ is the mass of the ion, $\omega(t)$ is the transverse confining
frequency, $V_{c}$ is the Coulomb interaction energy, $\gamma$ is the damping
coefficient. The force $(\theta_{xj},\theta_{yj},\theta_{zj})$ is is the
stochastic thermal force satisfying $\langle\theta_{\alpha,j}(t)\rangle=0$ and
$\langle\theta_{\alpha,j}(t)\theta_{\beta,k}(t^{\prime})\rangle=2m\gamma
k_{B}T\delta_{\alpha\beta}\delta_{jk}\delta(t-t^{\prime})$, where
$\langle...\rangle$ denotes ensemble averaging. The integration of the
equation of motion is performed using GPU accelerated OpenMM [31] framework,
and Ewald sums are used to approximate the Coulomb interactions in the $x$ and
$y$ directions.

To determine the location of defects in a given ion crystal configuration, one
must compute the local order parameter field $\psi$ using the individual ion
coordinates. Using equation (3) the order parameter in the 6-clock phase can
be written as

$\displaystyle z_{i}$ $\displaystyle=$ $\displaystyle\textrm{Re}[\Psi
e^{\textbf{K}\cdot\textbf{r}_{i}}]$ (7) $\displaystyle=$
$\displaystyle|\psi|\cos\left(\frac{\pi(2n_{i}+1)}{6}+\delta\Theta_{i}+\textbf{K}\cdot\textbf{r}_{i}\right),$
(8)

where $n_{i}\in\\{1,...,6\\}$ determines the clock state at the position of
the $i$th ion and $\delta\Theta_{i}$ is the fluctuation about this phase.
Denoting the splitting between the planes as $h\equiv\textrm{max}_{i}(|\langle
z_{i}\rangle|))$, where $z_{i}$ is the $z$-coordinate of an ion either in the
+ or - sublattice of the 3-plane structural phase, one finds that
$h=|\psi|\cos(\pi/6)\langle\cos(\delta\Theta)\rangle$ and equation (8) can be
written as

$z_{i}=\frac{h}{\cos(\pi/6)}\cos\left(\delta\Theta_{i}+\frac{\pi(2n_{i}+1)}{6}+\textbf{K}\cdot\textbf{r}_{i}\right).$
(9)

The values of clock state at each point, $n_{i}$, are determined by allocating
a phase value to an arbitrary chosen patch of three adjacent ions and then
assigning all other patches the best matching value relative to this chosen
reference. After assigning $n_{i}$, the correction term, $\delta\Theta_{i}$,
is obtained by solving numerically the non-linear equation (9) using gradient
descent algorithm.

Figure 4(a) shows a typical configuration with topological defects. The
locations of the defects are determined by finding localized regions on the
boundary of which the phase $\Theta$ winds an integer multiple of $2\pi$. The
topological charge $s$ of a defect is defined as a winding number along the
contour $C$ encircling the defect i.e.
$s=\frac{1}{2\pi}\int_{C}\partial_{l}\Theta dl$ where $l$ is the position
along the path of the chosen contour. In our simulations, we observe 4 types
of point defects: asters and vortices with charge +1, and cross-hairs and
anti-vortices with charge -1. The existence of such defects may be predicted
on general homotopy theoretic arguments [32]. Here, we have demonstrated that
in 3-plane Coulomb crystal all four types of defects are energetically stable.

The crystal configuration shown in the Figure 4(a) is in the 6-face clock
phase rather than the KT phase, as there are clear boundary lines between the
domains separating the regions with phase angles of the possible minimum
energy configuration. Increasing the temperature in the system blurs the
boundary lines until the system reaches a KT phase, where there is no
energetically preferred phase angles $\Theta$. In the KT phase there is no
long-range order in the system - the defects are mobile and oppositely charged
defects tend to annihilate. The long range order established by entering the
clock phase and the energy of the $Z_{6}$ domain walls, appears to have a
stabilizing effect on the point defects, reducing significantly their
mobility.

##  IV Kibble-Zurek mechanism of defect formation

Quenching the system from a disordered phase into the ordered phase at
different rates results in different defect densities as can be seen in Figure
4(b). The relationship between the number of defects and the quench rates
across symmetry breaking phase transitions was initially investigated by
Kibble in the context of cosmology [33] and Zurek in the context of condensed
matter [34], in what became known as Kibble-Zurek (KZ) theory. Experimentally,
the KZ mechanism for a 6-clock model has been investigated in the context of
ferroelectric materials [35, 36]. One should note, however, that in
ferroelectrics the transition happens in three dimensions where there is no KT
phenomena.

Lets review the KZ mechanism as applied to continuous second order phase
transitions. Consider approaching the critical point of a symmetry breaking
second order phase transition. The correlation length $\xi$, defined by
$\langle\psi(0,t)\psi(r,t)\rangle\sim e^{-r/\xi}$, diverges as a power law of
the control parameter $\xi=\xi_{0}/|\epsilon|^{\mu}$, where $\mu$ is the
critical exponent. The system also slows down on the approach to the critical
point i.e. the relaxation time, $\tau$, defined as
$\langle\psi(0)\psi(t)\rangle\sim e^{-t/\tau}$ diverges as
$\tau=\tau_{0}/|\epsilon|^{\mu}$. KZ mechanism proposes that the correlation
length freezes out when the relaxation time is equal to the time left until
the crossing of the point i.e. the freeze out time $\hat{t}$ is found by
solving $\tau(\hat{t})=\tau_{Q}$. This “cross-over” time $\hat{t}$ marks a
transition between the adiabatic dynamical regime, where the correlation
lengths adjust to its equilibrium value, and impulsive regime, where the
correlation length is fixed. For a linear quench $\epsilon=t/\tau_{Q}$, one
observes $\hat{t}=(\tau_{0}\tau^{\mu}_{Q})^{1/(1+\mu)}$ and the freeze-out
correlation length is
$\hat{\xi}=\xi(\hat{t})=\xi_{0}(\tau_{q}/\tau_{0})^{\nu/(1+\mu)}$. In the two
dimensional system, the number of defects $n$ is inversely proportional to the
square of the correlation length scale in the system i.e. $n\propto\xi^{-2}$.
Thus, the KZ prediction for the defect density in the end of the quench is
$n_{f}=\hat{\xi}^{-2}\sim(\tau_{q}/\tau_{0})^{-2\nu/(1+\mu)}$.

For KT phase transition, the same arguments applies except the correlation
length diverge exponentially at the critical point rather than algebraically
i.e. $\xi\sim\exp(a|\epsilon|^{-\mu}|)$ and
$\tau\sim\exp(b|\epsilon|^{-\nu}|)$, where $\epsilon$. In this case, the
equation for the freeze-out time, $\hat{t}=\tau_{Q}$ cannot be solved exactly,
but can still be determined numerically [24, 25]. An additional complication
is that the generated vortex defects can be highly mobile and the process of
annihilation of defects with opposite topological charges leads to the growth
of correlation length in the system. Thus, the number of defects in the
“impulsive” regime is not fixed but continuously decreases as the defects
annihilate. In the context of sudden quenches, this grows of correlation
length is known as coarsening.

In [37, 24], the authors propose that for a quench into the KT phase, one
should account for both KZ mechanism and coarsening expressing the time
evolution of correlation length as

Figure 5: Evolution of the average number of defects following quenches from
1-plane disordered phase into a 3-lane 6 face clock phase at finite rate
$1/\tau_{Q}$. The simulated system contained 9858 $\textrm{Be}^{+}$ ions at a
temperature of $20\mu$K and $\gamma=8.96\times 10^{-7}\textrm{ps}^{-1}$
confined in a periodic box in the $xy$ directions of size $1395.00\times
1376.98$ $\mu$m. The starting and ending frequencies were set to
$\omega_{z}^{(i)}=7.70\textrm{ Mhz}$ and $\omega_{z}^{(f)}=7.42\textrm{ MHz}$,
and the critical frequency is $\omega_{z}^{(i)}=7.60\textrm{ MHz}$. Figure 6:
Late time evolution of the average number of defects for three quenches at
different rates $1/\tau_{Q}$. The slopes of the dashed lines is -1, indicating
a $n\ln n\sim t^{-1}$ coarsening scaling.  Figure 7: Plot of the number of
defects at the end of the quench, $n_{f}$ as a function of quench rates.

$\displaystyle\xi(t)$ $\displaystyle=\begin{cases}\xi_{eq}(t),&\textrm{for
}t<\hat{t}\\\ \hat{\xi}+f(t),&\textrm{for }t\geq\hat{t}\end{cases},$ (10)

where $f(t)$ is the function representing the growth of the correlation length
due to coarsening. Equation (10) expresses the idea that before the KZ freeze-
out time, $\hat{t}$, the correlation length adopts its thermal equilibrium
value, and after crossing $\hat{t}$ the correlation length is growing via
coarsening. In a two dimensional system, which is quenched from disordered
into ordered phase in the presence of linear damping, one expects the
coarsening to proceed via diffusing law, where the correlation length grows as
the square root of time $f(t)\sim t^{1/2}$ [38]. Several studies noted that
the approach to diffusive law can be slow and that the coarsening is more
accurately described by including a logarithmic correction,
$f(t)\sim\left(t/\ln t\right)^{1/2}$ [39, 40].

In [24], a good agreement between the analytic expression (10) and the
numerical simulations of quenches in a two dimensional XY model is found. We
have carried out molecular dynamics simulations to verify whether the
dynamically crossing the 1-plane to 3-plane structural transition will result
in the same defect scaling behavior. We simulate the dynamics of $N=9858$
$\textrm{Be}^{+}$ ions at a temperature of $20\mu\textrm{K}$ and
$\gamma=8.96\times 10^{-7}\textrm{ps}^{-1}$ confined to a periodic box in an
$xy$ directions and a harmonic potential in the $z$ directions by numerically
integrating equation (4)-(6). The system is first thermalized at a confinement
frequency $\omega_{z}^{(i)}$ sufficiently far from the critical frequency
$\omega_{z}^{(c)}$ such that the correlation length is small i.e. of the order
of the lattice spacing. After that the confining frequency is decreased
linearly at a rate,
$\omega_{z}=\omega_{z}^{(i)}+t(\omega_{z}^{(f)}-\omega_{z}^{(i)})/\tau_{Q}$,
such that the system undergoes a transition between a 1-plane and 3-plane
structural phases at a rate $1/\tau_{Q}$. The defects are counted using the
method presented in Section III, and in order to obtain the ensemble averaged
defect number $\langle n(t)\rangle$ the simulations are carried out $\sim 140$
times for each $\tau_{Q}$.

Figure 5 shows the evolution of the defect density as function of time
following phase transition at different quench rates. Several dynamical
regimes can be seen in the figure. Initially, the system is in the 1-plane
phase far from phase transition point, the correlation length is small and the
density of defects is large. As $\omega_{z}$ approaches the critical
frequency, the defect density decreases adjusting to the new equilibrium
values. This regime crosses over to rapid relaxation, where the density of
vortices decreases significantly. According to the KZ theory, the crossover
occurs at $t>\hat{t}$. The relaxation continues at a slower pace consistent
with the power-law coarsening dynamics. Finally, the relaxation slows further
as the system enters the clock phase and the domains are stabilized.

In Figure 6 we zoom in into the late time evolution at three different quench
rates to verify whether the relaxation in the system is due to coarsening. The
growth of the correlation length due to coarsening in the KT phase is
predicted to follow $\xi\sim(t\ln t)^{1/2}$ and consequently the number of
defects follow a powerlaw with a logarithmic correction $n\ln n\sim t^{-1}$
[39]. In Figure 6, one can see that there is a region where the powerlaw with
the logarithmic correction is valid. However, this regime crosses over fairly
quickly to a regime of slower relaxation, which we believe is either due to
the stabilizing effect of entering into the clock phase or the pinning of the
vortices by the discreteness of the lattice.

Finally, Figure 7 shows the number of defects at the end of the quench at
$t=\tau_{Q}$ as a function of $\tau_{Q}$. We observe a powerlaw scaling with
an exponent of -0.80 i.e. $n(t=\tau_{Q})\propto\tau_{Q}^{-0.80}$. While an
accurate determination of the KZ scaling law involves numerical estimation of
$\hat{t}$, previous simulation results for quenches of the two dimensional XY
spin model found a powerlaw scaling of -0.72 at intermediate values of
$\tau_{Q}$. The exponent of -0.80 observed in our simulations is consistent
with these previous findings.

##  V Conclusion

In this paper we have studied the non-equilibrium dynamics of a planar Coulomb
crystals undergoing a structural transition from a 1-planar to 3-planar
configurations. The mapping to the Ginzburg-Landau theory reveals that this
phase transition corresponds to a transition from a disordered paramagnetic
phase to an ordered 6-clock phase with an intermediate KT phase. We used
molecular dynamics simulations to confirm that the KT and the 6-clock phase
support stable topological defect structures: vortices, anti-vortices, cross-
hairs and asters. The density of defects depends on the quench rate of the
structural transition as predicted by the Kibble-Zurek theory of defect
formation. We have verified that the defect scaling law is consistent with KZ
scaling previously observed in numerical simulation of two dimensional XY spin
model. Moreover, we have observed signatures of coarsening due to defect/anti-
defect annihilation, which follows a relaxation powerlaw with a logarithmic
correction.

Our work demonstrates that large planar ion Coulomb crystals can be used as
model for studying non-equilibrium statistical physics of vortices in an
effectively quasi-two-dimensional system. Such Coulomb crystals with open
boundary conditions can be realized in Penning traps. While the focus of this
paper was on the system with repulsive Coulomb interactions, similar
transition is expected for a lattice with dipolar interactions [23]. We hope
that our work will stimulate research towards experimental study of the
predicted topological defects and their rich collective dynamics.

## References

  * Keller  _et al._ [2019] J. Keller, D. Kalincev, T. Burgermeister, A. P. Kulosa, A. Didier, T. Nordmann, J. Kiethe, and T. Mehlstäubler, Probing time dilation in coulomb crystals in a high-precision ion trap, Phys. Rev. Applied 11, 011002 (2019).
  * Häffner  _et al._ [2008] H. Häffner, C. Roos, and R. Blatt, Quantum computing with trapped ions, Physics Reports 469, 155 (2008).
  * Bruzewicz  _et al._ [2019] C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage, Trapped-ion quantum computing: Progress and challenges, Applied Physics Reviews 6, 021314 (2019), https://doi.org/10.1063/1.5088164 .
  * Kim  _et al._ [2010] K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan, and C. Monroe, Quantum simulation of frustrated ising spins with trapped ions, Nature 465, 590 (2010).
  * Britton  _et al._ [2012] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of spins, Nature 484, 489 (2012).
  * Dubin and O’Neil [1999] D. H. E. Dubin and T. M. O’Neil, Trapped nonneutral plasmas, liquids, and crystals (the thermal equilibrium states), Rev. Mod. Phys. 71, 87 (1999).
  * Piacente  _et al._ [2010] G. Piacente, G. Q. Hai, and F. M. Peeters, Continuous structural transitions in quasi-one-dimensional classical wigner crystals, Phys. Rev. B 81, 024108 (2010).
  * Fishman  _et al._ [2008] S. Fishman, G. De Chiara, T. Calarco, and G. Morigi, Structural phase transitions in low-dimensional ion crystals, Phys. Rev. B 77, 064111 (2008).
  * Galván-Moya  _et al._ [2014] J. E. Galván-Moya, V. R. Misko, and F. M. Peeters, Generic ordering of structural transitions in quasi-one-dimensional wigner crystals, Phys. Rev. B 90, 094111 (2014).
  * Cornelissens  _et al._ [2000] Y. Cornelissens, B. Partoens, and F. Peeters, Transition from two-dimensional to three-dimensional classical artificial atoms, Physica E: Low-dimensional Systems and Nanostructures 8, 314 (2000).
  * Retzker  _et al._ [2008] A. Retzker, R. C. Thompson, D. M. Segal, and M. B. Plenio, Double well potentials and quantum phase transitions in ion traps, Phys. Rev. Lett. 101, 260504 (2008).
  * Pyka  _et al._ [2013] K. Pyka, J. Keller, H. L. Partner, R. Nigmatullin, T. Burgermeister, D. M. Meier, K. Kuhlmann, A. Retzker, M. B. Plenio, W. H. Zurek, A. del Campo, and T. E. Mehlstäubler, Topological defect formation and spontaneous symmetry breaking in ion coulomb crystals, Nature Communications 4, 2291 (2013).
  * Ulm  _et al._ [2013] S. Ulm, J. Roßnagel, G. Jacob, C. Degünther, S. T. Dawkins, U. G. Poschinger, R. Nigmatullin, A. Retzker, M. B. Plenio, F. Schmidt-Kaler, and K. Singer, Observation of the kibble-zurek scaling law for defect formation in ion crystals, Nature Communications 4, 2290 (2013).
  * Nigmatullin  _et al._ [2016] R. Nigmatullin, A. del Campo, G. De Chiara, G. Morigi, M. B. Plenio, and A. Retzker, Formation of helical ion chains, Phys. Rev. B 93, 014106 (2016).
  * Partner  _et al._ [2013] H. L. Partner, R. Nigmatullin, T. Burgermeister, K. Pyka, J. Keller, A. Retzker, M. B. Plenio, and T. E. Mehlstäubler, Dynamics of topological defects in ion coulomb crystals, New Journal of Physics 15, 103013 (2013).
  * Landa  _et al._ [2013] H. Landa, B. Reznik, J. Brox, M. Mielenz, and T. Schaetz, Structure, dynamics and bifurcations of discrete solitons in trapped ion crystals, New Journal of Physics 15, 093003 (2013).
  * Mielenz  _et al._ [2013] M. Mielenz, J. Brox, S. Kahra, G. Leschhorn, M. Albert, T. Schaetz, H. Landa, and B. Reznik, Trapping of topological-structural defects in coulomb crystals, Phys. Rev. Lett. 110, 133004 (2013).
  * Kiethe  _et al._ [2017] J. Kiethe, R. Nigmatullin, D. Kalincev, T. Schmirander, and T. E. Mehlstäubler, Probing nanofriction and aubry-type signatures in a finite self-organized system, Nature Communications 8, 15364 (2017).
  * Pruttivarasin  _et al._ [2011] T. Pruttivarasin, M. Ramm, I. Talukdar, A. Kreuter, and H. Häffner, Trapped ions in optical lattices for probing oscillator chain models, New Journal of Physics 13, 075012 (2011).
  * Timm  _et al._ [2020] L. Timm, H. Weimer, L. Santos, and T. E. Mehlstäubler, Energy localization in an atomic chain with a topological soliton, Phys. Rev. Research 2, 033198 (2020).
  * Lee and Cross [2011] T. E. Lee and M. C. Cross, Pattern formation with trapped ions, Phys. Rev. Lett. 106, 143001 (2011).
  * Bermudez  _et al._ [2011] A. Bermudez, J. Almeida, F. Schmidt-Kaler, A. Retzker, and M. B. Plenio, Frustrated quantum spin models with cold coulomb crystals, Phys. Rev. Lett. 107, 207209 (2011).
  * Podolsky  _et al._ [2016] D. Podolsky, E. Shimshoni, G. Morigi, and S. Fishman, Buckling transitions and clock order of two-dimensional coulomb crystals, Phys. Rev. X 6, 031025 (2016).
  * Asja and Leticia [2011] J. Asja and F. C. Leticia, Quench dynamics of the 2d xy model, J. Stat. Mech. (2011).
  * Deutschländer  _et al._ [2015] S. Deutschländer, P. Dillmann, G. Maret, and P. Keim, Kibble-zurek mechanism in colloidal monolayers, Proc Natl Acad Sci USA 112, 6925 (2015).
  * Gosh [1995] P. K. Gosh,  _Ion Traps_ (Claredon Press, Oxford, 1995).
  * Dubin [1993] D. H. E. Dubin, Theory of structural phase transitions in a trapped coulomb crystal, Phys. Rev. Lett. 71, 2753 (1993).
  * Mitchell  _et al._ [1998] T. B. Mitchell, J. J. Bollinger, D. H. E. Dubin, X.-P. Huang, W. M. Itano, and R. H. Baughman, Direct observations of structural phase transitions in planar crystallized ion plasmas, Science 282, 1290 (1998), https://science.sciencemag.org/content/282/5392/1290.full.pdf .
  * Berezinsky [1971] V. L. Berezinsky, Destruction of long range order in one-dimensional and two-dimensional systems having a continuous symmetry group. I. Classical systems, Sov. Phys. JETP 32, 493 (1971).
  * Kosterlitz and Thouless [1973] J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Eastman  _et al._ [2017] P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, R. P. Wiewiora, B. R. Brooks, and V. S. Pande, Openmm 7: Rapid development of high performance algorithms for molecular dynamics, PLOS Computational Biology 13, 1 (2017).
  * Mermin [1979] N. D. Mermin, The topological theory of defects in ordered media, Rev. Mod. Phys. 51, 591 (1979).
  * Kibble [1980] T. Kibble, Some implications of a cosmological phase transition, Physics Reports 67, 183 (1980).
  * Zurek [1985] W. H. Zurek, Cosmological experiments in superfluid helium?, Nature 317, 505 (1985).
  * Griffin  _et al._ [2012] S. M. Griffin, M. Lilienblum, K. T. Delaney, Y. Kumagai, M. Fiebig, and N. A. Spaldin, Scaling behavior and beyond equilibrium in the hexagonal manganites, Phys. Rev. X 2, 041022 (2012).
  * Lin  _et al._ [2014] S.-Z. Lin, X. Wang, Y. Kamiya, G.-W. Chern, F. Fan, D. Fan, B. Casas, Y. Liu, V. Kiryukhin, W. H. Zurek, C. D. Batista, and S.-W. Cheong, Topological defects as relics of emergent continuous symmetry and higgs condensation of disorder in ferroelectrics, Nature Physics 10, 970 (2014).
  * Biroli  _et al._ [2010] G. Biroli, L. F. Cugliandolo, and A. Sicilia, Kibble-zurek mechanism and infinitely slow annealing through critical points, Phys. Rev. E 81, 050101 (2010).
  * Mondello and Goldenfeld [1990] M. Mondello and N. Goldenfeld, Scaling and vortex dynamics after the quench of a system with a continuous symmetry, Phys. Rev. A 42, 5865 (1990).
  * Yurke  _et al._ [1993] B. Yurke, A. N. Pargellis, T. Kovacs, and D. A. Huse, Coarsening dynamics of the xy model, Phys. Rev. E 47, 1525 (1993).
  * Bray  _et al._ [2000] A. J. Bray, A. J. Briant, and D. K. Jervis, Breakdown of scaling in the nonequilibrium critical dynamics of the two-dimensional $\mathit{XY}$ model, Phys. Rev. Lett. 84, 1503 (2000).