11institutetext:  1 Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir
km 4, Torrejón de Ardoz, E-28850, Madrid, Spain
2 Cosmic Dawn Center (DAWN), Denmark
3 DTU-Space, Technical University of Denmark, Elektrovej 327, DK-2800 Kgs.
Lyngby, Denmark
4 Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
5 Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800,
9700AV Groningen, The Netherlands
6 Department of Astronomy, Stockholm University, Oscar Klein Centre, AlbaNova
University Centre, 106 91 Stockholm, Sweden
7 Dept. of Physics and Astronomy, University College London, Gower Street,
London WC1E 6BT, United Kingdom
8 UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill,
Edinburgh EH9 3HJ, UK
9 Centro de Astrobiología (CAB), CSIC-INTA, Camino Bajo del Castillo s/n,
E-28692 Villanueva de la Cañada, Madrid, Spain
10 I.Physikalisches Institut der Universität zu Köln, Zülpicher Str. 77, 50937
Köln, Germany
11 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200
Copenhagen, Denmark
12 Aix Marseille Université, CNRS, LAM (Laboratoire d’Astrophysique de
Marseille) UMR 7326, 13388, Marseille, France
13 Telespazio UK for the European Space Agency, ESAC, Camino Bajo del Castillo
s/n, 28692 Villanueva de la Cañada, Spain
14 School of Physics & Astronomy, Space Park Leicester, University of
Leicester, 92 Corporation Road, Leicester LE4 5SP, UK
15 Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The
Netherlands
16 Centre for Extragalactic Astronomy, Durham University, South Road, Durham
DH1 3LE, UK
17 University of Vienna, Department of Astrophysics Türkenschanzs- trasse 17,
1180 Vienna, Austria
18 Institute of Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-
Str 27, 8093 Zurich, Switzerland
11email<EMAIL_ADDRESS>
# MIDIS. Near-infrared rest-frame morphology of massive galaxies
at $3<z<5.5$ in the Hubble eXtreme Deep Field

L. Costantin 11 S. Gillman 2233 L. A. Boogaard 44 P. G. Pérez-González 11 E.
Iani 55 P. Rinaldi 55 J. Melinder 66 A. Crespo Gómez 11 L. Colina 11 T. R.
Greve 223377 G. Östlin 66 G. Wright 88 A. Alonso-Herrero 99 J. Álvarez-Márquez
11 M. Annunziatella 11 A. Bik 66 K. I. Caputi 2255 D. Dicken 88 A. Eckart 1010
J. Hjorth 1111 O. Ilbert 1212 I. Jermann 2233 A. Labiano 1313 D. Langeroodi
1111 F. Peißker 1010 J. P. Pye 1414 T. V. Tikkanen 1414 P. P. van der Werf
1515 F. Walter 44 M. Ward 1616 M. Güdel 17171818 T. K. Henning 44

###### Abstract

Context. Thanks to decades of observations using the Hubble Space Telescope
(HST), the structure of galaxies at redshift $z>2$ has been widely studied in
the rest-frame ultraviolet regime, which traces recent star formation from
young stellar populations. But, we still have little information about the
spatial distribution of the older, more evolved, stellar populations,
constrained by the rest-frame infrared portion of galaxies’ spectral energy
distribution.

Aims. We present the morphological characterization of a sample of 21 massive
galaxies ($\log(M_{\star}/M_{\odot})>9.5$) at redshift $3<z<5.5$. These
galaxies are observed as part of the Guaranteed Time Observations program
MIDIS with the Mid-Infrared Instrument (MIRI) onboard the James Webb Space
Telescope (JWST). The deep MIRI 5.6 $\mu$m imaging (28.64 mag $5\sigma$
depth), allows us to characterize for the first time the rest-frame near-
infrared structure of galaxies beyond cosmic noon, at higher redshifts than
possible with NIRCam, tracing their older stellar populations.

Methods. We derive the galaxies’ non-parametric morphology, focusing on the
Gini, $M_{20}$, concentration, and deviation statistics. Furthermore, we model
the galaxies’ light distribution with a single Sérsic component and derive
their parametric morphology.

Results. We find that at $z>3$ massive galaxies show a smooth distribution of
their rest-infrared light, strongly supporting the increasing number of
regular disk galaxies already in place at early epochs. On the contrary, the
ultraviolet structure obtained from HST observations (WFC3/F160W) is generally
more irregular, catching the most recent episodes of star formation.
Importantly, we find a segregation of morphologies across cosmic time, having
massive galaxies at redshift $z>4$ later-type morphologies compared to $z\sim
3$ galaxies. These findings suggest a transition phase in galaxy assembly and
central mass build up already taking place at $z\sim 3-4$.

Conclusions. MIRI provides unique information about the structure of the
mature stellar population of high-redshift galaxies, allowing us to unveil
that massive galaxies beyond cosmic noon are prevalently compact disk galaxies
with smooth mass distribution.

###### Key Words.:

Galaxies: evolution, Galaxies: formation, Galaxies: high-redshift, Galaxies:
structure

##  1 Introduction

Galaxy morphology is a key proxy of galaxy diversity, since it provides a
first glimpse on the physical processes involved in galaxy evolution. Indeed,
the structural evolution of galaxies across cosmic time is shown to be
strongly related to stellar mass and star formation history (Kauffmann et al.,
2003; Baldry et al., 2004), merger history (Lotz et al., 2008), and
environment (Blanton & Moustakas, 2009; Kormendy et al., 2010).

In the last two decades, our knowledge of galaxy structure beyond the local
Universe was based on studies making use of the Hubble Space Telescope (HST),
which has given us access to the rest-frame ultraviolet to optical morphology
of galaxies up to redshift $z\sim 2-3$. At these wavelengths, galaxies at
$z>1$ appear more irregular in their light distribution than local galaxies
(Conselice et al., 2000, 2008), and such peculiar systems dominate the galaxy
population beyond $z\sim 2.5$ (e.g., Buitrago et al., 2013; Huertas-Company et
al., 2015).

Based on almost 30 years of HST observational campaigns, multiple results
pointed to the conclusion that the Hubble sequence is established around at
$z\sim 1-2$ (Mortlock et al., 2013; Huertas-Company et al., 2016). But,
properly quantifying the morphological transformation of galaxies using HST
datasets suffers some limitations. A combination of insufficient spatial
resolution and limited red wavelength coverage makes it difficult to properly
resolve the first complex structures assemblying in the first $2-3$ Gyr after
the Big Bang (i.e., $z>2-3$). Thus, this observational limitation leaves open
the question whether the Hubble sequence was already in place at earlier
cosmic times then previously thought. Indeed, spectroscopic observations
already suggest an epoch of early disk assembly (e.g., Wisnioski et al., 2015;
Simons et al., 2017; Rizzo et al., 2020), with the first complex structures
starting to build-up at and beyond cosmic noon (e.g., Tacchella et al., 2015;
Costantin et al., 2021, 2022; Jin et al., 2024; Jegatheesan et al., 2024).

Now, for the first time, the James Webb Space Telescope (JWST; Gardner et al.,
2023) opens the possibility to characterize the detailed structure of the bulk
of the stellar population of highest-redshift galaxies ($z>3$) with
unprecedented level of detail (e.g., Ferreira et al., 2022, 2023; Kartaltepe
et al., 2023; Costantin et al., 2023a; Treu et al., 2023; Jacobs et al.,
2023). The prevailing conclusion of these initial studies is that the fraction
of galaxies with disk-like morphologies is higher than that inferred with HST,
although their exact nature still needs to be investigated (e.g., Vega-Ferrero
et al., 2024; Pandya et al., 2024).

In this context, the Mid-Infrared Instrument (MIRI; Rieke et al., 2015;
Bouchet et al., 2015; Wright et al., 2023; Dicken et al., 2024) onboard JWST,
probing the observed near-to-mid infrared regime (4.9 to 27.9 $\mu$m),
provides a huge jump in sensitivity compared to previous observatories at
these wavelengths ($\sim$10 times deeper than IRAC/Spitzer; Fazio et al.,
2004). MIRI allows to target the rest-frame near-infrared morphology of
galaxies at $z>3$ (up to $z\sim 5.5$), resolving structures at a few kpc scale
($FWHM\sim 0.2$ arcsec at 5.6 $\mu$m, corresponding to $\sim 1.6$ kpc at
$z=3$). Furthermore, MIRI sensitivity could also allow us to probe the near-
infrared morphology of galaxies at $5<z<10$ (i.e., F770W and F1000W bands;
$FWHM\sim 0.3$ arcsec at 10 $\mu$m, corresponding to $\sim 1.4$ kpc at
$z=10$).

At high redshift, the MIRI Deep Imaging Survey (MIDIS; Ostlin et al.  _in
prep._) is the best available dataset to address key open questions about the
clumpy/irregular distribution of the bulk of the stellar mass, as traced by
the older stellar populations (see e.g., Boogaard et al., 2023, Gillman et al.
in prep.). It consists of $\sim 41.34$ hours of net exposure time for deep
imaging of the Hubble eXtreme Deep Field (XDF; Illingworth et al., 2013) at
5.6 $\mu$m, with parallel observations of the surrounding area (see e.g.,
Pérez-González et al., 2023) using the Near Infrared Camera (NIRCam) and the
Near Infrared Imager and Slitless Spectrograph (NIRISS). Furthermore, 8.5
hours of net exposure time were dedicated to deep imaging in the F1000W band
(Pérez-González et al., 2024b; Iani et al., 2024b).

In this work, we study the near-infrared rest-frame morphologies of massive
galaxies at redshift $z>3$ in the XDF using MIDIS observations. We derive
their non-parametric and parametric morphology, characterizing the structure
of their old stellar population and looking at the morphological
transformation of the bulk of their stellar mass across cosmic time. The paper
is organized as follows. In Sect. 2 we describe the selection of the sample of
galaxies. In Sect. 3 we derive the non-parametric and parametric morphology,
discussing the implications in the context of galaxy evolution. Finally, in
Sect. 4 we summarize our results and provide our conclusions.

Throughout this work we assume a Planck Collaboration et al. (2020) cosmology
with $H_{0}=67.4$ km s-1 Mpc-1, $\Omega_{m}=0.315$, and
$\Omega_{\Lambda}=0.685$. We quote magnitudes in the AB system (Oke & Gunn,
1983).

Figure 1: Near-infrared rest-frame morphology (MIRI F560W band) of the 21
galaxies studied in this work. The cutouts are $35\times 35$ pixel2, which
corresponds to $2.1\times 2.1$ arcsec2. We report the angular resolution as
the FWHM of the PSF.  Table 1: Sample galaxies and morphological statistics.

ID | RA | DEC | $z$ | Stellar Mass |  mag5.6μm | Gini | $M_{20}$ | $C$ | $R_{\rm e}$ | $n$
---|---|---|---|---|---|---|---|---|---|---
| [degree] | [degree] |  |  [$\log(M_{\star}/M_{\odot}$)] | [mag] |  |  |  | [kpc] |
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
13994 | 53.1772202 | $-27.7873143$ | $3.05$ | $9.86$ | $23.65\pm 0.01$ | $0.49_{-0.01}^{+0.01}$ | $-1.67_{-0.08}^{+0.02}$ | $2.54_{-0.01}^{+0.01}$ | $0.86_{-0.03}^{+0.03}$ | $1.9_{-0.2}^{+0.2}$
13396 | 53.1492069 | $-27.7914675$ | $3.08$ | $9.63$ | $23.99\pm 0.01$ | $0.40_{-0.01}^{+0.01}$ | $-1.68_{-0.03}^{+0.03}$ | $2.47_{-0.03}^{+0.03}$ | $2.6_{-0.2}^{+0.2}$ | $2.2_{-0.2}^{+0.3}$
13073 | 53.1525720 | $-27.7938433$ | $3.12$ | $10.06$ | $23.76\pm 0.01$ | $0.50_{-0.01}^{+0.01}$ | $-1.67_{-0.02}^{+0.02}$ | $2.55_{-0.01}^{+0.01}$ | $0.77_{-0.01}^{+0.01}$ | $2.0_{-0.1}^{+0.1}$
14587 | 53.1784856 | $-27.7840349$ | $3.21$ | $10.11$ | $23.01\pm 0.01$ | $0.56_{-0.01}^{+0.01}$ | $-1.64_{-0.12}^{+0.03}$ | $2.77_{-0.02}^{+0.02}$ | $0.51_{-0.01}^{+0.01}$ | $1.34_{-0.04}^{+0.05}$
15458 | 53.1745867 | $-27.7778661$ | $3.31$ | $9.61$ | $25.10\pm 0.01$ | $0.53_{-0.01}^{+0.01}$ | $-1.70_{-0.02}^{+0.02}$ | $2.54_{-0.02}^{+0.02}$ | $0.47^{\star}$ | $0.59_{-0.06}^{+0.07}$
13846 | 53.1696650 | $-27.7880738$ | $3.34$ | $9.62$ | $24.83\pm 0.01$ | $0.43_{-0.01}^{+0.01}$ | $-1.56_{-0.04}^{+0.03}$ | $2.28_{-0.02}^{+0.02}$ | $0.98_{-0.04}^{+0.04}$ | $1.2_{-0.2}^{+0.2}$
17345 | 53.1770660 | $-27.7643540$ | $3.38$ | $9.80$ | $23.91\pm 0.01$ | $0.46_{-0.01}^{+0.01}$ | $-1.45_{-0.05}^{+0.04}$ | $2.72_{-0.37}^{+0.04}$ | $1.8_{-0.1}^{+0.2}$ | $4.5_{-0.7}^{+1.0}$
16368 | 53.1880571 | $-27.7715280$ | $3.41$ | $9.74$ | $24.42\pm 0.01$ | $0.48_{-0.01}^{+0.01}$ | $-1.67_{-0.03}^{+0.04}$ | $2.45_{-0.02}^{+0.02}$ | $0.62_{-0.05}^{+0.04}$ | $1.7_{-0.3}^{+0.3}$
15718 | 53.1595397 | $-27.7767456$ | $3.43$ | $9.64$ | $25.99\pm 0.02$ | $0.42_{-0.01}^{+0.01}$ | $-1.64_{-0.03}^{+0.04}$ | $2.60_{-0.02}^{+0.02}$ | $2.4_{-0.5}^{+0.8}$ | $9.0_{-1.2}^{+0.7}$
12448 | 53.1679012 | $-27.7979524$ | $3.45$ | $9.72$ | $24.32\pm 0.01$ | $0.50_{-0.01}^{+0.01}$ | $-1.63_{-0.02}^{+0.02}$ | $2.53_{-0.02}^{+0.01}$ | $0.46^{\star}$ | $0.63_{-0.08}^{+0.12}$
14463 | 53.1481266 | $-27.7844048$ | $3.58$ | $10.12$ | $23.42\pm 0.01$ | $0.48_{-0.01}^{+0.01}$ | $-1.65_{-0.02}^{+0.02}$ | $2.50_{-0.01}^{+0.01}$ | $0.88_{-0.04}^{+0.03}$ | $3.0_{-0.3}^{+0.4}$
16207 | 53.1741046 | $-27.7729638$ | $3.67$ | $9.74$ | $24.79\pm 0.01$ | $0.48_{-0.01}^{+0.01}$ | $-1.67_{-0.02}^{+0.02}$ | $2.44_{-0.01}^{+0.01}$ | $0.47^{\star}$ | $2.5_{-0.3}^{+0.3}$
16296 | 53.1564707 | $-27.7721989$ | $3.68$ | $9.58$ | $25.04\pm 0.01$ | $0.46_{-0.01}^{+0.01}$ | $-1.66_{-0.03}^{+0.05}$ | $2.56_{-0.02}^{+0.02}$ | $0.46^{\star}$ | $9_{-2}^{+1}$
16506 | 53.1617927 | $-27.7706417$ | $3.84$ | $9.90$ | $24.73\pm 0.01$ | $0.49_{-0.01}^{+0.01}$ | $-1.67_{-0.02}^{+0.08}$ | $2.52_{-0.01}^{+0.01}$ | $0.81_{-0.03}^{+0.02}$ | $1.3_{-0.2}^{+0.2}$
12384 | 53.1512255 | $-27.7981941$ | $3.90$ | $9.63$ | $24.15\pm 0.03$ | $0.39_{-0.01}^{+0.02}$ | $-1.35_{-0.09}^{+0.08}$ | $2.54_{-0.07}^{+0.07}$ | $2.6_{-0.2}^{+0.2}$ | $1.1_{-0.3}^{+0.4}$
15320 | 53.1987279 | $-27.7789023$ | $4.47$ | $9.73$ | $24.48\pm 0.02$ | $0.40_{-0.01}^{+0.01}$ | $-1.61_{-0.04}^{+0.07}$ | $2.32_{-0.03}^{+0.03}$ | $1.16_{-0.05}^{+0.06}$ | $1.2_{-0.2}^{+0.3}$
12781 | 53.1721209 | $-27.7950917$ | $4.83$ | $9.64$ | $25.62\pm 0.08$ | $0.31_{-0.02}^{+0.03}$ | $-1.4_{-0.1}^{+0.1}$ | $2.05_{-0.07}^{+0.07}$ | $0.7_{-0.2}^{+0.3}$ | $4_{-2}^{+2}$
14441 | 53.1926641 | $-27.7841464$ | $4.86$ | $9.89$ | $25.10\pm 0.04$ | $0.39_{-0.01}^{+0.01}$ | $-1.44_{-0.07}^{+0.09}$ | $2.25_{-0.04}^{+0.04}$ | $0.7_{-0.2}^{+0.2}$ | $4_{-1}^{+2}$
15552 | 53.1888042 | $-27.7770858$ | $4.92$ | $9.57$ | $24.49\pm 0.06$ | $0.36_{-0.04}^{+0.07}$ | $-0.80_{-0.07}^{+0.06}$ | $1.64_{-0.08}^{+0.15}$ | $1.7_{-0.1}^{+0.2}$ | $0.57_{-0.05}^{+0.11}$
17476 | 53.1702323 | $-27.7628467$ | $5.27$ | $9.88$ | $24.58\pm 0.04$ | $0.34_{-0.01}^{+0.01}$ | $-1.40_{-0.06}^{+0.07}$ | $2.04_{-0.03}^{+0.03}$ | $0.72_{-0.10}^{+0.08}$ | $0.9_{-0.3}^{+0.6}$
17541 | 53.1708041 | $-27.7622288$ | $5.28$ | $9.67$ | $24.69\pm 0.05$ | $0.38_{-0.02}^{+0.02}$ | $-1.49_{-0.12}^{+0.08}$ | $2.38_{-0.05}^{+0.04}$ | $1.1_{-0.3}^{+0.5}$ | $5_{-2}^{+3}$
111Columns: (1) CANDELS ID. (2) Right ascension. (3) Declination. (4)
Photometric redshift derived in Hainline et al. (2023). (5) Stellar mass
derived in Pérez-González et al. (2008) and Barro et al. (2019). (6) Magnitude
in the F560W band. (7)-(9) Gini, $M_{20}$, and $C$ statistics. (10) Parametric
effective radius. (11) Parametric Sérsic index. Galaxies marked by the ⋆
symbol have unreliable size or Sérsic index, on the edge of the parameter
space. Errors are reported as the 16th–84th percentile interval.

##  2 Data and Sample

We select galaxies from the MIRI $5.6\leavevmode\nobreak\ \mu$m imaging of the
JWST Guaranteed Time Observations (GTO) program MIDIS (PID: 1283; Ostlin et
al.  _in prep._), which represents the deepest image of the Universe at these
wavelengths (28.64 mag $5\sigma$ depth, calculated within circular apertures
with diameter 0.45 arcsec and corrected for drizzling correlation; Rinaldi et
al., 2023a, b; Iani et al., 2024a; Boogaard et al., 2023). The data were
calibrated using the procedure described in Ostlin et al.  _in prep._ , using
a modified version of the official JWST pipeline version 1.12.3 (pmap 1137)
and applying a background homogenization algorithm (including 1/f-noise
removal) prior to obtaining the final mosaic drizzled at a pixel scale of 0.06
arcsec (see also Pérez-González et al., 2024a, c, for more details).

Due to the complex observational strategy of MIDIS, and the lack of bright and
isolated stars in the field, we use a varying Point Spread Function (PSF)
model built using empirical PSFs at different positions on the MIRI detector.
In particular, the empirical PSFs are created by stacking stars in images from
the commissioning program ID 1473 and two calibration programs ID 1040 and ID
1024, which are then processed according to the MIDIS observational strategy.
All the details about these empirical PSFs are enumerated in Boogaard et al.
(2023).

Galaxies are initially selected for having redshift $z>3$, using the
photometric estimations of the JWST Advanced Deep Extragalactic Survey (JADES;
Eisenstein et al., 2023; Hainline et al., 2023; Williams et al., 2023) and
stellar masses $\log(M_{\star}/M_{\odot})>9.5$, according to estimations
detailed in Pérez-González et al. (2008) and Barro et al. (2019). The total
number of such galaxies in the MIDIS field is 28. Visually inspecting the
images of the initial sample, we discard four faint galaxies that are not
detected in the MIRI imaging and one galaxy that is at the edge of the MIRI
pointing. After this selection, the sample consists of 23 galaxies. We finally
discarded two galaxies at redshift $z>5.5$, probing for the first time the the
rest-frame regime $\gtrsim 0.9$ $\mu$m for 21 galaxies up to $z=5.5$, which is
unique of the MIDIS dataset and where there is essentially no variations of
morphological non-parametric indicators (e.g., Ren et al., 2024).

The main properties of the sample galaxies are enumerated in Table 1. In
particular, photometry is performed for all sources as described in Gillman et
al.  _in prep_. Briefly, galaxies are detected and modeled in F560W MIRI band,
and then forced photometry in HST, NIRCam, and MIRI bands is performed with
The Farmer (Weaver et al., 2022), allowing the flux to vary, whilst keeping
the structural parameters fixed. It is finally worth noticing that galaxies
15458 and 12448 are Ly$\alpha$ emitters identified in Iani et al. (2024a),
galaxy 14443 is presented in the ALMA-selected sample analyzed in Boogaard et
al. (2023), while galaxy 14587 is in the sample of X-ray Active Galactic
Nuclei detailed in Gillman et al.  _in prep._

##  3 Method and Results

In this work we present the non-parametric and parametric morphology for a
sample of 21 massive galaxies at $z>3$ using MIRI imaging at 5.6 $\mu$m. Given
the redshift range of our sample ($3<z\lesssim 5.5$), we are probing the rest-
frame near-infrared structure of galaxies ($\lambda_{\rm rest}\sim 0.9-1.4$
$\mu$m), which best trace more evolved stellar population that only MIRI can
probe up to this high redshift. To complement our analysis, we compare the
rest-frame near-infrared structure of these galaxies with their rest-frame
ultra-violet spatial distribution using HST/F160W imaging from CANDELS. Apart
from the different stellar populations probed by these two bands, our choice
is also justified by the similar width of the PSF of these two datasets,
minimizing effects related to spatial resolutions.

###  3.1 Non-parametric morphology

We measure the structure of the sample galaxies observed in the MIRI/F560W
filter. Firstly, we create segmentation maps of each galaxy with SEP (Barbary,
2018), a Python library implementing SExtractor (Bertin & Arnouts, 1996).
Then, we derive non-parametric morphological diagnostics using statmorph
(Rodriguez-Gomez et al., 2019). Following Crespo Gómez et al. (2024), we
quantify the uncertainties associated to each parameter by performing 500
Monte Carlo realizations for each galaxy. We perturb each image pixel-by-pixel
with gaussian noise ($\sigma=$ root mean square of the background level) and
create new segmentation maps at each iteration. The 16th-84th percentile
interval of each parameter is used to quantify their uncertainty.

In the following, we focus on the Gini and $M_{20}$ statistics (see Table 1
and Figs. 2-4). The Gini coefficient quantifies the light distribution of the
galaxy, ranging from $G=0$ for homogeneous brightness distribution to $G=1$
when the entire flux is concentrated in a single pixel (Abraham et al., 2003;
Lotz et al., 2004). The $M_{20}$ statistic is the normalized second-order
moment of the brightest 20% of the galaxy’s flux and it provides valuable
information about the spatial distribution of any sub-structure, such as
bright nuclei, bars, or spiral arms (Lotz et al., 2004). Thus, the
Gini-$M_{20}$ diagram has been largely employed to separate early, late-type,
and merging galaxies from low to high redshift (Lotz et al., 2004, 2008; Rose
et al., 2023; Crespo Gómez et al., 2024).

Figure 2: Gini-$M_{20}$ diagram. Orange stars correspond to non-parametric
morphology measured on MIRI F560W images. The green points and shaded region
correspond to non-parametric morphology measured on WFC3 F160W images. Green
and light red contours represent the $2\sigma$ levels at (observed) 1.6 $\mu$m
and 5.6 $\mu$m, respectively. ETGs, LTGs, and mergers are separated according
to Eqs. (4) in Lotz et al. (2008).  Figure 3: Gini-$M_{20}$ diagram based on
MIRI F560W morphology, separating galaxies at $3<z<4$ (red shaded region) from
those at $4<z<5.5$ (blue shaded region). The inset panel shows the redshift
and absolute magnitude distributions of the sample, with galaxies at $3<z<4$
as red dots and galaxies at $4<z<5.5$ as blue squares.

Furthermore, to better characterize the sample galaxies, using statmorph we
derive the concentration– asymmetry–smoothness (CAS; Conselice, 2003) and
multimode–intensity–deviation statistics (MID; Freeman et al., 2013), which
allow us to better discriminate the true nature of galaxies in the
Gini-$M_{20}$ diagram. In particular, we complement the information inferred
from the Gini-$M_{20}$ diagram with that provided by the concentration and
deviation statistics (see Table 1 and Fig. 4, top panels). The concentration
index is defined as the ratio of the radii which enclose 80% and 20% of the
light of the galaxy (Conselice, 2003; Lotz et al., 2004). As a consequence,
elliptical galaxies are the most concentrated systems, and the concentration
decreases for later Hubble types (Bershady et al., 2000). The deviation index
(Freeman et al., 2013) quantifies the variation of irregular or peculiar
morphologies from elliptical symmetry. In this case, it is expected that
elliptical or regular disk galaxies with no substructures (e.g., clumps) show
values of deviation clustering near zero.

It is worth noting that while both Gini and $M_{20}$ usually correlate with
concentration, they differ in important respects. On one side, the Gini
coefficient is independent of the large-scale spatial distribution of light,
and high values of Gini may also arise if the bright structures are not
located in the center of a galaxy (unlike $C$). On the other side, $M_{20}$
(scaling as $r^{2}$) is more heavily weighted than $C$ by the spatial
distribution of bright regions. Moreover, $M_{20}$ is more sensitive than $C$
to merger signatures, since it is not measured within circular or elliptical
apertures with a fixed center (see Lotz et al., 2004, for further details).

Figure 4: Gini-$M_{20}$ diagram, color-coded according to concentration (top
left panel), deviation (top right panel), effective radius (bottom left
panel), and Sérsic index (bottom right panel). Galaxies are separated in two
redshift bins: $3<z<4$ (circles) and $4<z<5.5$ (squares).

####  3.1.1 Gini-$M_{20}$

In Fig. 2, we show the distribution of rest-frame near-infrared morphologies
in the (Gini, $M_{20}$) diagram (orange stars and $2\sigma$ contour). The
first result is that all galaxies in our sample occupy the region of late-type
systems, even if some (i.e., 14587 and 15458) are quite close to the edge of
early-type objects. The most peculiar galaxy in the sample is 15552, which has
Gini = 0.36, $M_{20}=-0.80$ and lies on the bottom right side of the diagram.
Looking at its visual morphology (Fig. 1), it is the only galaxy showing an
irregular surface-brightness distribution (see also Fig. 4).

As a comparison, we derive the rest-frame ultraviolet morphologies of the
sample galaxies using HST/WFC3 imaging from CANDELS (green points, shaded
region, and $2\sigma$ contour). Already at $z>3$, galaxies appear more
irregular at shorter wavelengths, pointing to the presence of an underlying
already mature population with a smooth stellar mass distribution already in
place when the Universe was $\lesssim 2$ Gyr old. Despite the low statistics,
we find that in the rest-frame near-infrared regime galaxies have Gini${}_{\rm
NIR}=0.46_{-0.08}^{+0.03}$ and $M_{20,\rm NIR}=-1.64_{-0.04}^{+0.22}$, while
in the ultraviolet regime they have Gini${}_{\rm UV}=0.42_{-0.03}^{+0.10}$ and
$M_{20,\rm UV}=-1.4_{-0.2}^{+0.2}$. Thus, galaxies move from later-type
morphologies to earlier-type ones if observed at shorter or longer
wavelengths, respectively. This trend is consistent with previous results
based on ultraviolet and optical differences of typical star-forming galaxies
at $z\sim 3-4$ reported in Conselice et al. (2008) and Wuyts et al. (2012),
and could be explained as evidence for disk assembly through the inward
migration of clumps and gas accretion.

In Fig. 3, we show the (Gini, $M_{20}$) diagram, separating galaxies at
$3<z<4$ (red shaded region) from galaxies at $4<z<5.5$ (blue shaded region).
The former show median Gini${}_{\rm NIR}=0.48_{-0.06}^{+0.02}$ and $M_{20,\rm
NIR}=-1.66_{-0.02}^{+0.08}$, while the latter have Gini${}_{\rm
NIR}=0.37_{-0.04}^{+0.02}$ and $M_{20,\rm NIR}=-1.43_{-0.09}^{+0.15}$. We see
that galaxies in the higher redshift bin appear more irregular than galaxies
at lower redshift, suggesting a segregation of morphologies from irregular to
smoother light (and stellar mass) profiles in massive galaxies.

####  3.1.2 Gini-$M_{20}$: concentration and deviation

To better quantify the structure of the sample galaxies, we derive their
concentration (Fig. 4, top left panel) and deviation statistic (Fig. 4, top
right panel).

The concentration statistic reflects the trend previously highlighted, with
lower redshift galaxies being more concentrated (median
$C=2.54_{-0.09}^{+0.05}$) than higher redshift ones (median
$C=2.11_{-0.07}^{+0.21}$). As expected, the galaxy’s concentration increases
moving from LTGs to ETGs along empirical trend Gini $\propto-M_{20}$, which is
also used to separate normal from merging systems (Lotz et al., 2008). On the
other hand, the deviation statistic allows us to identify a peculiar galaxy,
that is, 15552. This galaxy is classified as a late-type galaxy, but its
morphology is highly irregular (see also Fig. 1).

###  3.2 Parametric morphology

We performed the two-dimensional photometric decomposition of the MIRI images,
modeling the surface brightness of each galaxy with a Sérsic (1968) law. We
make use of anduryl (Marrero de la Rosa et al.  _in prep._) a computational
tool designed for 2D photometric fitting of galaxies, leveraging Bayesian
inference as its foundational framework. Presently, the software accommodates
fits based on either a Sérsic model or a combination of Sérsic and exponential
models. At the heart of anduryl lies a robust parameter space exploration,
executed through Nested Sampling, a technique for approximating the posterior
probability integral implemented through
Nestle222http://kylebarbary.com/nestle/index.html. This method enables the
inference of the posterior probability distribution, thereby facilitating the
extraction of marginal posterior distributions for each parameter. From these
distributions, the mean values emerge as the most plausible estimates for each
parameter.

Looking at the size (half-light radius $R_{\rm e}$) and Sérsic index of the
sample galaxies, we find no clear trend (Fig. 4, bottom panels). Lower
redshift galaxies have median $R_{\rm e}=0.8_{-0.3}^{+1.4}$ kpc and Sérsic
index $n=1.8_{-0.7}^{+2.3}$, while higher redshift galaxies have median
$R_{\rm e}=0.9_{-0.2}^{+0.3}$ kpc and Sérsic index $n=2_{-2}^{+2}$. Overall,
the median values for the entire sample are $R_{\rm e}=0.8_{-0.0.3}^{+1.0}$
kpc and $n=1.8_{-0.9}^{+2.6}$, which could be considered typical of late-type
morphologies, but with a wide scatter to high values.

###  3.3 Comparison with previous works at shorter wavelengths

It is worth recalling that our work is pioneering in exploring the rest-frame
near-infrared morphology of high-redshift galaxies at wavelengths beyond 5
$\mu$m, thanks to the sensitivity and spatial resolution of MIRI at these
wavelengths. However, different studies already analyzed the evolution of the
morphological and structural properties of galaxies up to $z\sim 8-9$ at
shorter wavelength (up to $\sim 4\leavevmode\nobreak\ \mu$m) using NIRCam.
These studies are mainly based on the visual classification of galaxies or the
value of their parametric morphology (i.e., Sérsic index).

Kartaltepe et al. (2023) visually classified and measured the morphological
properties of 850 galaxies in CEERS (Cosmic Evolution Early Release Science;
Finkelstein et al., 2022; Bagley et al., 2023) using NIRCam/F277W for systems
at $3<z<4$, NIRCam/F356W for systems at $4<z<5$, and NIRCam/F444W for systems
at $z>4.5$. Our sample galaxies have smaller sizes and higher Sérsic indeces
than the pure disk class in Kartaltepe et al. (2023), which have median
$R_{\rm e}=1.19_{-0.48}^{+0.61}$ kpc and $n=1.16^{+0.88}_{-0.40}$. But, our
galaxies have larger sizes and smaller Sérsic indeces than the pure spheroid
class in Kartaltepe et al. (2023), which have median $R_{\rm
e}=0.54_{-0.32}^{+0.42}$ kpc and $n=2.46^{+2.32}_{-1.41}$. Similarly, Ferreira
et al. (2023) reported that disks and peculiar galaxies at $3<z<6$ exhibit
similar Sérsic profiles with $n\sim 1$, although spheroids show $n\sim 1.6$.
Interestingly, the population of compact disk-like galaxies is missing in
cosmological simulations (i.e.; TNG50 Pillepich et al., 2018; Nelson et al.,
2018; Costantin et al., 2023b), which tend to present larger sizes (median
$R_{\rm e,F356W}=2.06^{+0.89}_{-0.62}$) and lower values of Sérsic index
(median $n_{\rm F356W}=0.89^{+0.37}_{-0.26}$). It is worth noticing that in
Kartaltepe et al. (2023) the separation of disk- or spheroidal-dominated
galaxies based on the Sérsic index (i.e., $n=2$) identify 71% of the visually
identified disks and only 45% of the visually identified spheroids, implying
that the interpretation of such high-redshift galaxies based on one-component
models should be further investigated (see also Costantin et al., 2018; Gao et
al., 2020, for local galaxies).

Finally, our results are broadly compatible with the trend derived in Huertas-
Company et al. (2023) using $\sim 20,000$ galaxies in the CEERS survey and
Convolutional Neural Networks trained on HST/WFC3 labeled images and domain-
adapted to JWST/NIRCam. They found that at $3<z<6$ the disk-like systems
already dominate over the early-type population and bulge-dominated galaxies
also exist at these early epochs. This result holds when looking at different
cosmological fields (Lee et al., 2023), but with different fractions of
disturbed and regular galaxies among different studies. In Huertas-Company et
al. (2023), the fraction of galaxies with disturbed optical rest-frame
morphologies is predominant (80% disturbed, 20% regular) at $M_{\star}\sim
10^{9-10}\leavevmode\nobreak\ M_{\odot}$. On the other hand, in Ferreira et
al. (2023) disks and peculiar galaxies are present in similar abundance (40%),
while disks totally dominate the galaxy population in Lee et al. (2023). This
apparent discrepancy raises a note of caution in comparing morphological
studies using different classifications schemes.

##  4 Summary and Conclusions

In this work, we present the first detailed rest-frame near-infrared
morphological study of a sample of 21 massive galaxies at $3<z<5.5$, observed
as part of the MIRI Deep Imaging Survey in the XDF at 5.6 micron (F560W).

We employ non-parametric morphological diagnostics to classify the sample
galaxies, mainly focusing on the Gini and $M_{20}$, as well as concentration,
deviation, and parametric Sérsic modelling of the galaxies.

From this study we can draw two main conclusions. From one side, massive
galaxies at $z>3$ show more regular structures in the restframe near-infrared
regime (MIRI F560W) compared to a more peculiar morphology at shorter
wavelengths. We measure the morphology of galaxies using the F160W band
(WFC3/HST), probing the restframe ultraviolet regime of these galaxies,
finding that their structure is more irregular. Second, separating galaxies at
$4<z<5.5$ from those at $3<z<4$, at similar absolute magnitude we find a
segregation of galaxy morphologies across cosmic time in the Gini-$M_{20}$
diagram, from later to earlier types. This could be interpreted as a
transition phase in galaxy assembly, where the complex structures (i.e.,
bulges, bars) start to assemble first. Overall, the visual, non-parametric
(Gini-$M_{20}$, concentration, deviation), and parametric (Sérsic and $R_{\rm
e}$) analysis points to a compact population of disk-like galaxies, with a
mostly regular mass distribution.

Building on the results of this work, deep MIRI campaigns in the near future
will open the possibility to expand the morphological characterization of
galaxies in the early Universe to a larger sample and possibly to higher
redshift (e.g., up to $z=10$ using MIRI F770W and F1000W).

###### Acknowledgements.

The project that gave rise to these results received the support of a
fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code
is LCF/BQ/PR24/12050015. LC and PGPG acknowledge support from grant
PID2022-139567NB-I00 funded by Spanish Ministerio de Ciencia e Innovación
MCIN/AEI/10.13039/501100011033, FEDER _Una manera de hacer Europa_. LC, JA-M,
AC-G, and LC acknowledge support by grant PIB2021-127718NB-100 from the
Spanish Ministry of Science and Innovation/State Agency of Research
MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. SG
acknowledges financial support from the Villum Young Investigator grant 37440
and 13160 and the Cosmic Dawn Center (DAWN), funded by the Danish National
Research Foundation under grant DNRF140. LAB acknowledge support by ERC AdG
grant 740246 (Cosmic-Gas). EI acknowledges funding from the Netherlands
Research School for Astronomy (NOVA). JM and AB acknowledge support from the
Swedish National Space Administration (SNSA). AAH acknowledges support from
grant PID2021-124665NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by
ERDF A way of making Europe”. AE and KIC acknowledges funding from the
Netherlands Research School for Astronomy (NOVA) and from the Dutch Research
Council (NWO) through the award of the Vici Grant VI.C.212.036. IJ acknowledge
support from the Carlsberg Foundation (grant no CF20-0534) and the Cosmic Dawn
Center is funded by the Danish National Research Foundation under grant No.
140. FP acknowledge support through the German Space Agency DLR 50OS1501 and
DLR 50OS2001 from 2015 to 2023. JPP and TVT acknowledge financial support from
the UK Science and Technology Facilities Council, and the UK Space Agency.
This work was supported by research grants (VIL16599, VIL54489) from VILLUM
FONDEN. The work presented is the effort of the entire MIRI team and the
enthusiasm within the MIRI partnership is a significant factor in its success.
MIRI draws on the scientific and technical expertise of the following
organisations: Ames Research Center, USA; Airbus Defence and Space, UK; CEA-
Irfu, Saclay, France; Centre Spatial de Liége, Belgium; Consejo Superior de
Investigaciones Científicas, Spain; Carl Zeiss Optronics, Germany; Chalmers
University of Technology, Sweden; Danish Space Research Institute, Denmark;
Dublin Institute for Advanced Studies, Ireland; European Space Agency,
Netherlands; ETCA, Belgium; ETH Zurich, Switzerland; Goddard Space Flight
Center, USA; Institute d’Astrophysique Spatiale, France; Instituto Nacional de
Técnica Aeroespacial, Spain; Institute for Astronomy, Edinburgh, UK; Jet
Propulsion Laboratory, USA; Laboratoire d’Astrophysique de Marseille (LAM),
France; Leiden University, Netherlands; Lockheed Advanced Technology Center
(USA); NOVA Opt-IR group at Dwingeloo, Netherlands; Northrop Grumman, USA;
Max-Planck Institut für Astronomie (MPIA), Heidelberg, Germany; Laboratoire
d’Etudes Spatiales et d’Instrumentation en Astrophysique (LESIA), France; Paul
Scherrer Institut, Switzerland; Raytheon Vision Systems, USA; RUAG Aerospace,
Switzerland; Rutherford Appleton Laboratory (RAL Space), UK; Space Telescope
Science Institute, USA; Toegepast- Natuurwetenschappelijk Onderzoek (TNO-TPD),
Netherlands; UK Astronomy Technology Centre, UK; University College London,
UK; University of Amsterdam, Netherlands; University of Arizona, USA;
University of Cardiff, UK; University of Cologne, Germany; University of
Ghent; University of Groningen, Netherlands; University of Leicester, UK;
University of Leuven, Belgium; University of Stockholm, Sweden; Utah State
University, USA. A portion of this work was carried out at the Jet Propulsion
Laboratory, California Institute of Technology, under a contract with the
National Aeronautics and Space Administration. We would like to thank the
following National and International Funding Agencies for their support of the
MIRI development: NASA; ESA; Belgian Science Policy Office; Centre Nationale
D’Etudes Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur
Luft-und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y
Competitividad; Netherlands Research School for Astronomy (NOVA); Netherlands
Organisation for Scientific Research (NWO); Science and Technology Facilities
Council; Swiss Space Office; Swedish National Space Board; UK Space Agency.
For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to the Author Accepted Manuscript version arising
from this submission.

## References

  * Abraham et al. (2003) Abraham, R. G., van den Bergh, S., & Nair, P. 2003, ApJ, 588, 218
  * Bagley et al. (2023) Bagley, M. B., Finkelstein, S. L., Koekemoer, A. M., et al. 2023, ApJ, 946, L12
  * Baldry et al. (2004) Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681
  * Barbary (2018) Barbary, K. 2018, SEP: Source Extraction and Photometry, Astrophysics Source Code Library, record ascl:1811.004
  * Barro et al. (2019) Barro, G., Pérez-González, P. G., Cava, A., et al. 2019, ApJS, 243, 22
  * Bershady et al. (2000) Bershady, M. A., Jangren, A., & Conselice, C. J. 2000, AJ, 119, 2645
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Blanton & Moustakas (2009) Blanton, M. R. & Moustakas, J. 2009, ARA&A, 47, 159
  * Boogaard et al. (2023) Boogaard, L. A., Gillman, S., Melinder, J., et al. 2023, arXiv e-prints, arXiv:2308.16895
  * Bouchet et al. (2015) Bouchet, P., García-Marín, M., Lagage, P. O., et al. 2015, PASP, 127, 612
  * Buitrago et al. (2013) Buitrago, F., Trujillo, I., Conselice, C. J., & Häußler, B. 2013, MNRAS, 428, 1460
  * Conselice (2003) Conselice, C. J. 2003, ApJS, 147, 1
  * Conselice et al. (2000) Conselice, C. J., Bershady, M. A., & Jangren, A. 2000, ApJ, 529, 886
  * Conselice et al. (2008) Conselice, C. J., Rajgor, S., & Myers, R. 2008, MNRAS, 386, 909
  * Costantin et al. (2018) Costantin, L., Corsini, E. M., Méndez-Abreu, J., et al. 2018, MNRAS, 481, 3623
  * Costantin et al. (2023a) Costantin, L., Pérez-González, P. G., Guo, Y., et al. 2023a, Nature, 623, 499
  * Costantin et al. (2021) Costantin, L., Pérez-González, P. G., Méndez-Abreu, J., et al. 2021, ApJ, 913, 125
  * Costantin et al. (2022) Costantin, L., Pérez-González, P. G., Méndez-Abreu, J., et al. 2022, ApJ, 929, 121
  * Costantin et al. (2023b) Costantin, L., Pérez-González, P. G., Vega-Ferrero, J., et al. 2023b, ApJ, 946, 71
  * Crespo Gómez et al. (2024) Crespo Gómez, A., Colina, L., Álvarez-Márquez, J., et al. 2024, arXiv e-prints, arXiv:2402.18672
  * Dicken et al. (2024) Dicken, D., García Marín, M., Shivaei, I., et al. 2024, arXiv e-prints, arXiv:2403.16686
  * Eisenstein et al. (2023) Eisenstein, D. J., Willott, C., Alberts, S., et al. 2023, arXiv e-prints, arXiv:2306.02465
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
  * Ferreira et al. (2022) Ferreira, L., Adams, N., Conselice, C. J., et al. 2022, ApJ, 938, L2
  * Ferreira et al. (2023) Ferreira, L., Conselice, C. J., Sazonova, E., et al. 2023, ApJ, 955, 94
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2022, ApJ, 940, L55
  * Freeman et al. (2013) Freeman, P. E., Izbicki, R., Lee, A. B., et al. 2013, MNRAS, 434, 282
  * Gao et al. (2020) Gao, H., Ho, L. C., Barth, A. J., & Li, Z.-Y. 2020, ApJS, 247, 20
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001
  * Hainline et al. (2023) Hainline, K. N., Johnson, B. D., Robertson, B., et al. 2023, arXiv e-prints, arXiv:2306.02468
  * Huertas-Company et al. (2016) Huertas-Company, M., Bernardi, M., Pérez-González, P. G., et al. 2016, MNRAS, 462, 4495
  * Huertas-Company et al. (2023) Huertas-Company, M., Iyer, K. G., Angeloudi, E., et al. 2023, arXiv e-prints, arXiv:2305.02478
  * Huertas-Company et al. (2015) Huertas-Company, M., Pérez-González, P. G., Mei, S., et al. 2015, ApJ, 809, 95
  * Iani et al. (2024a) Iani, E., Caputi, K. I., Rinaldi, P., et al. 2024a, ApJ, 963, 97
  * Iani et al. (2024b) Iani, E., Rinaldi, P., Caputi, K. I., et al. 2024b, arXiv e-prints, arXiv:2406.18207
  * Illingworth et al. (2013) Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6
  * Jacobs et al. (2023) Jacobs, C., Glazebrook, K., Calabrò, A., et al. 2023, ApJ, 948, L13
  * Jegatheesan et al. (2024) Jegatheesan, K., Johnston, E. J., Häußler, B., & Nedkova, K. V. 2024, arXiv e-prints, arXiv:2402.00959
  * Jin et al. (2024) Jin, Y., Zhu, L., Zibetti, S., et al. 2024, A&A, 681, A95
  * Kartaltepe et al. (2023) Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54
  * Kormendy et al. (2010) Kormendy, J., Drory, N., Bender, R., & Cornell, M. E. 2010, ApJ, 723, 54
  * Lee et al. (2023) Lee, J. H., Park, C., Hwang, H. S., & Kwon, M. 2023, arXiv e-prints, arXiv:2312.04899
  * Lotz et al. (2008) Lotz, J. M., Davis, M., Faber, S. M., et al. 2008, ApJ, 672, 177
  * Lotz et al. (2004) Lotz, J. M., Primack, J., & Madau, P. 2004, AJ, 128, 163
  * Mortlock et al. (2013) Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2013, MNRAS, 433, 1185
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624
  * Oke & Gunn (1983) Oke, J. B. & Gunn, J. E. 1983, ApJ, 266, 713
  * Pandya et al. (2024) Pandya, V., Zhang, H., Huertas-Company, M., et al. 2024, ApJ, 963, 54
  * Pérez-González et al. (2024a) Pérez-González, P. G., Barro, G., Rieke, G. H., et al. 2024a, arXiv e-prints, arXiv:2401.08782
  * Pérez-González et al. (2023) Pérez-González, P. G., Costantin, L., Langeroodi, D., et al. 2023, ApJ, 951, L1
  * Pérez-González et al. (2008) Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234
  * Pérez-González et al. (2024b) Pérez-González, P. G., Rinaldi, P., Caputi, K. I., et al. 2024b, arXiv e-prints, arXiv:2402.16942
  * Pérez-González et al. (2024c) Pérez-González, P. G., Rinaldi, P., Caputi, K. I., et al. 2024c, arXiv e-prints, arXiv:2402.16942
  * Pillepich et al. (2018) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Ren et al. (2024) Ren, J., Liu, F. S., Li, N., et al. 2024, arXiv e-prints, arXiv:2404.16686
  * Rieke et al. (2015) Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584
  * Rinaldi et al. (2023a) Rinaldi, P., Caputi, K. I., Costantin, L., et al. 2023a, ApJ, 952, 143
  * Rinaldi et al. (2023b) Rinaldi, P., Caputi, K. I., Iani, E., et al. 2023b, arXiv e-prints, arXiv:2309.15671
  * Rizzo et al. (2020) Rizzo, F., Vegetti, S., Powell, D., et al. 2020, Nature, 584, 201
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez, V., Snyder, G. F., Lotz, J. M., et al. 2019, MNRAS, 483, 4140
  * Rose et al. (2023) Rose, C., Kartaltepe, J. S., Snyder, G. F., et al. 2023, ApJ, 942, 54
  * Sérsic (1968) Sérsic, J. L. 1968, Atlas de Galaxias Australes
  * Simons et al. (2017) Simons, R. C., Kassin, S. A., Weiner, B. J., et al. 2017, ApJ, 843, 46
  * Tacchella et al. (2015) Tacchella, S., Carollo, C. M., Renzini, A., et al. 2015, Science, 348, 314
  * Treu et al. (2023) Treu, T., Calabrò, A., Castellano, M., et al. 2023, ApJ, 942, L28
  * Vega-Ferrero et al. (2024) Vega-Ferrero, J., Huertas-Company, M., Costantin, L., et al. 2024, ApJ, 961, 51
  * Weaver et al. (2022) Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11
  * Williams et al. (2023) Williams, C. C., Tacchella, S., Maseda, M. V., et al. 2023, ApJS, 268, 64
  * Wisnioski et al. (2015) Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209
  * Wright et al. (2023) Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003
  * Wuyts et al. (2012) Wuyts, S., Förster Schreiber, N. M., Genzel, R., et al. 2012, ApJ, 753, 114


# Elasticity versus phase field driven motion in the phase field crystal model

Amit Acharya1 1Department of Civil and Environmental Engineering, and Center
for Nonlinear Analysis, Carnegie Mellon University, Pittsburgh, Pennsylvania
15213, USA.
2PoreLab, Njord Centre, Department of Physics, University of Oslo, P.O. Box
1048, 0316 Oslo, Norway.
3School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA.  Luiza Angheluta2 1Department of Civil and Environmental
Engineering, and Center for Nonlinear Analysis, Carnegie Mellon University,
Pittsburgh, Pennsylvania 15213, USA.
2PoreLab, Njord Centre, Department of Physics, University of Oslo, P.O. Box
1048, 0316 Oslo, Norway.
3School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA.  and Jorge Viñals3 1Department of Civil and Environmental
Engineering, and Center for Nonlinear Analysis, Carnegie Mellon University,
Pittsburgh, Pennsylvania 15213, USA.
2PoreLab, Njord Centre, Department of Physics, University of Oslo, P.O. Box
1048, 0316 Oslo, Norway.
3School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA.

###### Abstract

The inherent inconsistency in identifying the phase field in the phase field
crystal Theory with the material mass and, simultaneously, with material
distortion is discussed. In its current implementation, elastic relaxation in
the phase field crystal occurs on a diffusive time scale through a dissipative
permeation mode. The very same phase field distortion that is included in
solid elasticity drives diffusive motion, resulting in a non physical
relaxation of the phase field crystal. We present two alternative theories to
remedy this shortcoming. In the first case, it is assumed that the phase field
only determines the incompatible part of the elastic distortion, and therefore
one is free to specify an additional compatible distortion so as to satisfy
mechanical equilibrium at all times (in the quasi static limit). A numerical
solution of the new model for the case of a dislocation dipole shows that,
unlike the classical phase field crystal model, it can account for the known
law of relative motion of the two dislocations in the dipole. The physical
origin of the compatible strain in this new theory remains to be specified.
Therefore, a second theory is presented in which an explicit coupling between
independent distortion and phase field accounts for the time dependence of the
relaxation of fluctuations in both. Preliminary details of its implementation
are also given.

##  1 Introduction

The phase field crystal (PFC) model has become a successful and widely used
description of the mesoscale response of a nonequilibrium crystalline phase
[1]. Since, by construction, the need to resolve the time scale associated
with lattice vibration is eliminated, its computational implementation can
efficiently access long time phenomena that are difficult to describe by
molecular dynamics simulation [2, 3]. One of the strengths of the model is the
deceivingly simple description of defected solids. A longstanding debate, to
which we wish to contribute here, is the nature of the phase field itself, and
the proper long wavelength description of a defected elastic material.

A widely used interpretation is that the phase field represents the mass
density of the solid [4, 5]. In fact, phase field crystal model free energies
have been derived by using Density Functional Theory methods, with the
expectation of obtaining functionals of the mass density as the single,
relevant order parameter that describe the structure and elastic response of a
crystalline solid in equilibrium [6, 7]. In addition, the free energies
obtained have been shown to provide a reasonable description of the freezing
phase transition [8]. However, extensions so that the set of slow variables in
the Density Functional Theory also includes the momentum density have not been
considered to date (except for colloidal systems [9]). Therefore, the
constitutive laws relating momentum currents (stresses) and their conjugate
strains, as well as their coupling to the phase field, remain completely
phenomenological.

Most phase field crystal model implementations begin with the separation of
the mass density into two components: A slowly varying part in space $\rho$ –
which is then identified with the standard hydrodynamic mass density – and a
second component, $\psi$, which is obtained from the location of the Bragg
peaks of the crystalline structure under consideration [5], and that describes
crystalline order. The mass density is decomposed as,

$\rho(\mathbf{x})+\psi(\mathbf{x})=\rho(\mathbf{x})+\sum_{n}A_{n}e^{i\mathbf{q}^{(n)}\cdot\mathbf{x}}.$

where the $A_{n}=|A_{n}|e^{i\theta_{n}}$ are slowly varying complex
amplitudes. The phase $\theta_{n}$ represents a lattice distortion
$\theta_{n}=-i\mathbf{q}^{(n)}\cdot\mathbf{u}(\mathbf{x},t))$. The vectors
$\mathbf{q}^{(n)}$ are the critical eigenmodes of the model’s free energy at
the bifurcation point to the modulated phase, and $\mathbf{u}(\mathbf{x},t)$
is the displacement field [4, 5, 10]. Although this relation is only expected
to hold for weakly distorted lattices, it is commonly used to study defected
configurations as well. It is one of the goals of this paper to discuss the
limitations inherent in this mapping between the phase field $\psi$ and the
displacement field $\mathbf{u}$.

The phase field crystal is one among several theories that allow a short scale
regularization of defect core divergences inherent in classical elasticity,
while allowing for the treatment of defect topology and motion, grain boundary
energetics and motion, and explicit expressions for the associated mobilities.
For equilibrium studies, the only constitutive input required is the (defect
free) equilibrium free energy, a non convex functional of the phase field,
with a minimizer that corresponds to a spatially periodic configuration. For
nonequilibrium problems, current implementations of the phase field crystal
generally assume gradient flows driving the evolution of a configuration
towards a local minimizer of the free energy functional. Importantly for the
regularization, the phase field is smooth everywhere, including in defected
configurations, but its phase is singular at defect cores. The proper
treatment of this singularity is the main focus of this paper.

##  2 Thermodynamics and Linear Response

The defining feature of a solid material subject to a small distortion is its
elastic response under the application of an external body force
$\mathbf{f}_{ext}(\mathbf{x})$. Within linear elasticity, the energy of the
distortion is

${\cal F}=\frac{1}{2}\int
d^{d}x\;C_{ijkl}\bm{\epsilon}_{ij}\bm{\epsilon}_{kl}-\int
d^{d}x\;\mathbf{f}_{ext}\cdot\mathbf{u}\quad{\rm
with,}\quad\bm{\epsilon}_{ij}=\frac{1}{2}\left(U_{ij}+U_{ji}\right)$ (1)

with the distortion tensor $U_{ij}=\partial_{i}u_{j}$, and $\mathbf{u}$ the
displacement vector, with the usual definitions of strain $\epsilon_{ij}$ and
elastic constants $C_{ijkl}$. In contrast, subjecting a phase field crystal to
strain involves an indirect transformation, as it involves imposing an
externally fixed chemical potential $\mu_{ext}$, and requiring that the
distorted equilibrium configuration minimizes ${\cal F}={\cal F}_{sh}-\int
d^{d}x\;\mu_{ext}(\mathbf{x})\psi(\mathbf{x})$, where ${\cal F}_{sh}$ is
typically chosen to be of the Swift-Hohenberg type [11]. Alternatively, one
can directly impose the deformation on the order parameter $\psi$ directly. In
neither case is it possible to, say, specify body forces or tractions on the
boundary of the phase field crystal.

This is not a shortcoming in equilibrium. Let $\delta\psi$ be the local
variation of the phase field $\psi$ following a small distortion of the
equilibrium lattice. Since mass is conserved, this variation is
$\delta\psi=-\partial_{k}(\psi\delta u_{k})$, where $u_{k}$ is the $k$-th
component of the displacement vector. From this relation, the variation
$\delta{\cal F}/\delta u_{k}=\psi\partial_{k}(\delta{\cal F}/\delta\psi)$
follows as the change in free energy due to the distortion. The thermodynamic
stress $\mathbf{T}$ is defined from the same energy as
$\partial_{j}T_{ij}=-\delta{\cal F}/\delta u_{i}$, and therefore
$\partial_{j}T_{ij}=-\psi\partial_{i}(\delta{\cal F}/\delta\psi)$. This
relation is correct in equilibrium where both sides of the equation
independently vanish. However, this is not the case outside of equilibrium:
The right hand side of the equation is but one contribution to the reversible
stress, and the equality does not hold in general. Note that the equality
between what is an elastic body force $\partial_{j}T_{ij}$ and a
configurational force $-\psi\partial_{i}(\delta{\cal F}/\delta\psi)$ follows
directly from the assumption that lattice displacement and variation of $\psi$
are not independent, but rather related through
$\delta\psi=-\partial_{k}(\psi\delta u_{k})$.

In nonequilibrium studies, this difficulty can be remedied, although somewhat
artificially, by assuming that the slow component $\rho$ is defined to be the
mass density, which is indeed conserved, but that the crystalline component
$\psi$ is not, regardless of the fact that both represent a mass density
(e.g., Ref. [5]). As a consequence, both fields are now subject to independent
variation. Under these conditions, the authors of Ref. [5] find that the mass
velocity $\mathbf{v}$ is different than the lattice velocity
$\partial_{t}\mathbf{u}$,

$\mathbf{v}=\partial_{t}\mathbf{u}+\frac{1}{2A_{0}^{2}}\mu_{\eta}\frac{\delta{\cal
F}}{\delta\mathbf{u}},$

where $A_{0}$ is the value of the equilibrium amplitude $A_{n}$ in the
hexagonal lattice studied in that reference, and $\mu_{\eta}$ is the mobility
associated with the gradient flow of $\psi$. The difference between the two
velocities is a dissipative term, and it corresponds to a permeation mode that
decouples lattice and mass motion.

Within this decomposition, and with explicit consideration of kinetic energy
terms $\frac{1}{2}\int\rho v^{2}d^{d}x$ and momentum conservation, the
following dispersion relation for the transverse displacement
$\mathbf{u}_{\perp}\cdot\mathbf{q}=0$ in the overdamped limit results,

$\omega_{\perp}(q)=-\frac{i}{2}B\mu_{\eta}q^{2}\pm\frac{q}{2}\sqrt{\frac{B}{\rho}\left(12A_{0}^{2}-B\mu_{\eta}^{2}\rho
q^{2}\right)}.$ (2)

where $B/2$ is the coefficient of the gradient kernel in the Swift-Hohenberg
free energy. This is an interesting result, as the decoupling allows for the
emergence of propagating modes at low $q$, not just the purely diffusive
response of the classical phase field crystal. This result needs to be
contrasted with the classical dispersion relation for transverse modes of an
isotropic crystal (with Newtonian viscosity) which is
$\omega=-i(\eta/2\rho)q\pm q\sqrt{\mu/\rho-\eta^{2}/(4\rho^{2})}$, with $\mu$
the shear modulus, and $\eta$ the shear viscosity. The speed of propagation of
transverse distortions is $c_{\perp}=\sqrt{\mu/\rho}$. The associated
transverse susceptibility is
$\chi(\mathbf{q},\omega)=\mathbf{u}_{\perp}(\mathbf{q},\omega)/\mathbf{f}_{\perp,ext}(\mathbf{q},\omega)=(1/\rho)\left(-\omega^{2}+(\mu-i\omega\eta)q^{2}/\rho\right)^{-1}$.
At zero frequency, the response is that of an elastic solid $1/(\mu q^{2})$.
For the phase field crystal the inverse susceptibility at zero frequency is
$\chi^{-1}(\omega=0)\propto(3BA_{0}^{2}/\rho)q^{2}+{\cal O}(q^{4})$ and hence
the phase has a solid response at low frequency with shear modulus
proportional to $3BA_{0}^{2}/\rho$. Two features are of interest. First, as
$q$ increases, the square root in Eq. (2) becomes imaginary, and the response
of the phase field crystal is purely diffusive. This is relevant for defected
solids because gradients of the phase field amplitude are large near cores.
Second, the imaginary part of $\omega_{\perp}$ is due entirely to permeation.
In a conventional elastic solid, the imaginary part is the shear kinematic
viscosity, which is essentially zero for practical purposes. Hence an elastic
solid remains distorted under tension, and it can store elastic energy. In the
phase field crystal, on the other hand, the ratio of the imaginary to the real
part of the susceptibility is $\omega(\rho\mu_{\eta}q)/(3\phi^{2})$. The
denominator is small (on the order of the dimensionless distance to the
ordering transition for the phase field crystal), and hence permeation becomes
negligible only for very small values of the mobility $\mu_{\eta}$, or as
$q\rightarrow 0$. Making the mobility very small is tantamount to eliminating
dissipative relaxation of $\psi$, the raison d’etre of the model in the first
place. Alternatively, elastic response is recovered only in the limit
$q\rightarrow 0$ which is not appropriate for the study of elastic solids.
Therefore, whereas an elastic solid reversibly stores elastic energy under
traction, the response of the phase field crystal is primarily diffusive.
Ultimately, the elastic response and the permeation mode in the phase field
crystal both result from the curvature of the planes of constant $\psi$, both
proportional to the transverse diffusion coefficient, which is proportional to
the same coefficient $B$. Although this difficulty has been recognized for a
long time, and a number of modified models have been introduced to allow for
relaxation of the phase field in a time scale faster than diffusion [4, 12,
11, 5, 13], it still remains an open question.

Finally, as concerns the dissipative mechanical response of crystalline
solids, it must be recognized that the phenomenon of plasticity is strongly
dissipative (with greater than $50\%$ of applied power being converted to
heat, more often close to $90\%$) and this dissipation arises primarily due to
energy taken up in incoherent atomic vibrations produced due to the motion of
crystal defects through the body.

##  3 Defected phase field crystal

The difficulties just described become more serious when the crystal is
allowed to have topological singularities. In this case, the permeation mode
just described acquires topological content. Let $\alpha_{ij}$ be the
dislocation density tensor defined so that its integral over a surface is the
sum of the Burgers vectors of the dislocation lines piercing the surface [14].
For an isolated dislocation line it can be written as
$\alpha_{ij}=t_{i}b_{j}\delta(\bm{\zeta})$, where $\mathbf{t}$ is the unit
tangent to the line, $\mathbf{b}$ is the Burgers vector of the line, and
$\bm{\zeta}$ is the two dimensional radial coordinate (with origin at the
line) on the plane normal to the line tangent. The tensor $\bm{\alpha}$
describes the incompatibility of the lattice distortion as

$\alpha_{ik}=-\epsilon_{ilm}\partial_{l}U_{mk},$ (3)

where $\bm{\epsilon}$ here is the antisymmetric tensor. An immediate
consequence of this result is that even though the phase field $\psi$ is
regular everywhere, since its phase $\theta_{n}$ is not, the displacement
field is not uniquely determined by a configuration of the phase field $\psi$,
as was given, for example, in the relation between the phase of $\psi$ and the
displacement in the previous section. Furthermore, if one were to assume that
$\mathbf{v}=\partial_{t}\mathbf{u}$, or that
$\partial_{i}v_{k}=\partial_{t}U_{ik}$, then simply by taking the time
derivative of Eq. (3), one has $\partial_{t}\bm{\alpha}=0$, the dislocation
line does not move. Dislocation motion requires the explicit uncoupling of
mass and lattice motion as in,

$\partial_{i}v_{k}=\partial_{t}U_{ik}-{\cal J}_{ik}.$ (4)

Substitution into the time derivative of Eq. (3) yields,

$\partial_{t}\alpha_{ik}+\epsilon_{ilm}\partial_{l}{\cal J}_{mk}=0,$ (5)

which is precisely the statement of conservation of topological charge under
the motion of the dislocation line. This result is standard in continuum
dislocation theory [15, 14] (for extensions to finite deformations, see [16],
[17, App. B]). The dislocation charge density current can be written as as
${\cal J}_{ik}=\epsilon_{irs}\alpha_{rk}(v_{s}+V_{s})$, where $\mathbf{V}$ is
the velocity of the dislocation line relative to the material velocity
$\mathbf{v}$ [18]. Unlike the purely dissipative character of the permeation
mode of Sec. 2, the dislocation density current has both reversible and
dissipative contributions.

In physical terms, it remains to be resolved whether the phase field crystal
model (an intrinsically dissipative theory) can approximate to some degree
both reversible topological charge density currents as well as the local
dissipative processes around dislocation cores during line motion, including
therefore Peierls stresses and barriers. The phase field crystal gives an
explicit expression for ${\cal J}_{mk}$ in terms of the time rate of change of
the phase field amplitude [19]. To the extent that in this latter treatment
the evolution of the phase field is purely dissipative, it follows that the
defect current generated by the motion of defects is a source of energy
dissipation.

###  3.1 Decoupling through a compatible distortion

One first possible approach at making the phase field $\psi$ independent from
distortion is given in Ref. [20]. Given that the phase field crystal can
describe topological singularities and their motion, but not elasticity, it
seems reasonable to use $\psi$ to determine the incompatible part of the
distortion only, and then introduce a new, independent, compatible distortion.
The former relaxes in plastic time scales as $\psi$ evolves diffusively,
whereas the second relaxes quickly and accounts for elastic response.

While the phase field is regular at defects, its phase is singular, fact that
can be used to define the Burgers vector density in two dimensions, or the
dislocation density tensor in three dimensions. In two dimensions, one has

$\oint d\theta_{n}=q_{j}^{(n)}b_{j}$

where the Burgers vector is defined as $\oint d\mathbf{u}=-\mathbf{b}$ (Fig.
1, see Ref. [10] for details).

Figure 1: Amplitudes and phases of the phase field $\psi$ for a dislocation
dipole in a two dimensional hexagonal lattice. Only two of the three wave
amplitudes of the order parameter go to zero at the defect core. The
corresponding phases exhibits a discontinuity along a line joining the cores.
The asymmetry in the amplitude contours reflects the different orientation of
the two base waves which are rotated by $120^{o}$ relative to each other.

The singular phases can be used to directly compute a configurational
distortion tensor, as will be described in Sec. 3.2 (see also Ref. [21]).

Since the phase field crystal satisfies its own evolution equation, it is not
necessary to compute the incompatible distortion explicitly. However, the
compatible part must be determined, and one can do so by using the requirement
that the crystal is in mechanical equilibrium at all times (which is
reasonable in the diffusive time scale of dislocation motion determined by the
evolution of $\psi$). In order to determine the compatible distortion, it is
possible, within linear elasticity, to use the incompatible part as the source
of stress. A configurational stress is first defined, which is regular even at
dislocation cores. Consider the change in phase field free energy $\Delta{\cal
F}={\cal F}[\psi(\mathbf{x}^{\prime})]-{\cal F}[\psi(\mathbf{x})]$ due to a
small affine distortion
$\mathbf{x}^{\prime}=\mathbf{x}+\mathbf{u}(\mathbf{x})$. Then a
configurational stress $T_{ij}^{\psi}=\langle\frac{\partial
f}{\partial(\partial_{i}u_{j})}\rangle$ is defined where $f$ is the free
energy density, and $\langle.\rangle$ stands for an average over the unit cell
of the phase field crystal in order to remove spatially oscillatory
components. As is well known, $\nabla\cdot T_{ij}^{\psi}\neq 0$ in general.
The method of Ref. [20] aims at introducing an independent (and compatible)
distortion $\mathbf{U}^{\delta}$, so that the combined corresponding stresses
satisfy $\nabla\cdot(\mathbf{T}^{\psi}+\mathbf{T}^{\delta})=0$. By
construction, this distortion is compatible, and can be integrated to find the
corresponding displacement $\mathbf{u}^{\delta}$ everywhere. The compatible
displacement is computed therefore from $\nabla\cdot({\cal
C}:\nabla\mathbf{u}^{\delta})=-\nabla\cdot\mathbf{T}^{\psi}$ where ${\cal C}$
is the tensor of elastic constants. In the case of a two dimensional hexagonal
lattice (elastically isotropic) considered in Ref. [20], the total stress
$\mathbf{T}=\mathbf{T}^{\psi}+\mathbf{T}^{\delta}$ was expressed in terms of
the Airy stress function $\chi$ as
$T_{ij}=\epsilon_{ik}\epsilon_{jl}\partial_{k}\partial_{l}\chi$. Then,

$(1-\frac{\lambda}{2(\lambda+\mu)})\nabla^{4}\chi=\epsilon_{ik}\epsilon_{jl}\partial_{i}\partial_{j}T_{kl}^{\psi}-\frac{\lambda}{2(\lambda+\mu)}\nabla^{2}T_{kk}^{\psi},$
(6)

in terms of the Lame coefficients of the lattice. Given an instantaneous
configuration of $\psi$, the right hand side of Eq. (6) is known. Solving for
$\chi$, one then has

$T_{ij}^{\delta}=\epsilon_{ik}\epsilon_{jl}\partial_{k}\partial_{l}\chi-
T_{ij}^{\psi},$ (7)

and, from it, the corresponding distortion $\mathbf{U}^{\delta}$.

It is the case then that if a fixed configuration of the phase field is
transformed as
$\psi(\mathbf{x})\rightarrow\psi(\mathbf{x}+\mathbf{u}^{\delta})$, then the
new stress $\mathbf{T}^{\psi}$ is divergenceless. Furthermore, since the
distortion $\mathbf{u}^{\delta}$ is compatible, this transformation does not
change the topological content of the field $\psi$. The temporal evolution of
the phase field is assumed in Ref. [20] to include the dissipative evolution
given by the gradient flow induced by ${\cal F}$, plus and adiabatic
distortion by $\mathbf{u}^{\delta}$ at each time, so that the combined motion
takes place on the restricted manifold defined by
$\nabla\cdot\mathbf{T}^{\psi}=0$. Clearly, this evolution is arbitrary, and
unphysical. It fails to describe the physical mechanism by which the phase
field relaxes in a fast time scale to reach mechanical equilibrium for a
given, fixed, distribution of defects. In this model, the relaxation is
infinitely fast and accomplished by a simple displacement on $\psi$.

The resulting defect core motion observed deviates significantly from the
results given by the classical phase field crystal method (Fig. 2). Similar
conclusions were already reported in Ref. [22] that noted that different
methods of lattice distortion application led to different dislocation motion.
Two limiting cases were considered: a uniform, affine distortion of the
reference lattice (termed “rigid displacement”), and a “relaxational method”
in which a thin layer of boundary sites were displaced, and the rest of the
lattice allowed to relax under phase field diffusion. The motion of a
dislocation in the same two dimensional, hexagonal lattice studied here showed
a substantial dependence on the method of lattice distortion used.

Figure 2: Relative velocity of two edge dislocations of opposite Burgers
vector in a hexagonal lattice as a function of their dimensionless separation.
Left: evolution given by a purely dissipative phase field crystal. The
parameter $r$ in the box equals the coefficient of the quadratic term in the
phase field free energy. Right: evolution given by the model described in this
section [10]. Oscillations for large values of $r$ in both cases (far from the
bifurcation threshold to the hexagonal phase) correspond to transient pinning
effects described earlier [23].

###  3.2 The phase field as a configuration tensor

A different scheme in which to treat the distortions of the phase field
crystal and of the underlying material as independent has been given in Ref.
[21]. It combines the Field Dislocation Mechanics model (FDM, cf. [24, 25, 26,
27]) – demonstrated to represent features like large ‘stress-free’ shear bands
in the presence of significant material inertia, sub-inter-supersonic
dislocation motion in accord with molecular dynamic simulations in [28],
dislocation annihilation, dissociation, dipole and pile-up interactions and
equilibria, and Peierls stress (all without any parameter adjustment) – with
the phase field crystal methodology. The strength of the FDM methodology is in
providing a framework for calculating (nonlinear) elastic fields in bodies of
arbitrary shape and elastic anisotropy, and for evolution accounting for the
topological charge conservation of dislocations. Suitably interpreted, it
provides flexibility in the representation of lattice periodicity in a
continuum pde model. The work in [29] utilized Generalized Stacking fault
energies (GSFs) (determinable from interatomic potentials [30], as also
utilized in phase Field methods, e.g. [31]), a building block analogous to the
choice of a free energy functional from the Swift-Hohenberg family in the
phase field crystal. A goal of [21] is to use Swift-Hohenberg type models
instead of GSFs.

In Ref. [21], the phase field crystal is used to define a configurational
distortion tensor $\mathbf{P}$, independent of the material distortion, and
that is well defined in a defected medium. Given a general coordinate
transformation defined as the following transformation of the basis vectors,
$\overline{e}_{i}=P^{m}_{.\;i}e_{m}$, the induced transformation in reciprocal
space is

$q^{i}=(P^{-1\ T})^{i}_{.j}\ \overline{q}^{\ j}$

One can use this result to define a configurational distortion tensor from the
local phase gradients of the phase field (the distorted configuration), and a
the reciprocal basis of the equilibrium, undistorted configuration,

$\mathbf{q}_{0}^{i}=P^{-1\ T}\nabla\theta_{i}.$ (8)

The configurational distortion tensor $\mathbf{P}$ is a point wise functional
of the phase field $\psi$, which is defined so as to coincide with the inverse
elastic distortion of the medium $\mathbf{W}$ only in equilibrium. Away from
equilibrium, one allows relative fluctuations between the two so that the
elastic response of the medium is captured by $\mathbf{W}$ whereas the
diffusive relaxation of the phase field is described by $\mathbf{P}$.

In order to further motivate this definition, it is useful to consider the
small deformation limit. Note that the authors of Ref. [21] use a different
sign convention for the Burgers vector $\oint d\mathbf{u}=\mathbf{b}$, and
also the transpose of the distortion above ($U_{ij}=\partial_{j}u_{i}$ for a
compatible distortion). The inverse elastic distortion satisfies
$\mathbf{W}=\mathbf{I}-\mathbf{U}$, where $\mathbf{U}$ is the fundamental
measure of elastic distortion. Analogously, one defines for small deformations
$\mathbf{Q}=\mathbf{I}-\mathbf{P}$. Now if the phase of the order parameter is
tied to the distortion,
$\theta_{n}=\mathbf{q}^{(n)}\cdot\mathbf{r}-\mathbf{q}^{(n)}\cdot\mathbf{u}$,
and therefore
$\partial_{j}\theta_{n}=q^{(n)}_{j}-q^{(n)}_{i}\partial_{j}u_{i}=q^{(n)}_{j}-q^{(n)}_{i}U_{ij}=q^{(n)}_{i}(\delta_{ij}-U_{ij})$.
Inverting it to put in the form of Eq. (8) one finds
$P_{ij}=\delta_{ij}-U_{ij}$. Therefore $\mathbf{Q}=\mathbf{U}$ in the case in
which the phase of the phase field can be expressed as indicated. However,
they are treated in the theory as independent. The theory considers an elastic
energy functional of the distortion, a phase field free energy of the Swift-
Hohenberg type, and a penalty term which is a functional of relative
fluctuations between the elastic distortion $\mathbf{U}$ and the
configurational distortion $\mathbf{Q}$,

${\cal F}(\mathbf{U},\psi)={\cal F}_{el}(\mathbf{U})+C_{sh}{\cal
F}_{sh}(\psi)+\frac{C_{w}}{2}\int d^{d}x\;(\mathbf{U}-\mathbf{Q})_{sym}^{2}$
(9)

where $C_{sh}$ and $C_{w}$ are two coupling constants. In the governing
equations for the nonequilibrium evolution of the fields to be shown below,
the constant $C_{w}$ corresponds to an inverse time scale over which the
configurational distortion and the elastic distortion relax towards each
other.

An initially defected configuration will be described by an order parameter
field $\psi$ so that topological defects will be located in regions of non
zero curl of $\mathbf{Q}$, defined by the point wise oriented triad in
reciprocal space, generally not orthonormal, as given. If $C_{w}$ and the
scale of ${\cal F}_{sh}$, $C_{sh}$, are simultaneously large, the phase field
will relax quickly (and diffusively) to a local minimum of the constrained
free energy, relatively independently of the changes in elastic energy that
the relaxation might impose. This process will be accompanied by a relaxation
of the elastic distortion in phonon time scales, very quickly if the
quasistatic limit is invoked. Further evolution will be slow, driven by a
Peach-Kohler force resulting from both mechanical stress and phase field
configurational distortion (14).

Given the free energy in Eq. (9), a (mechanical) dissipation inequality is
introduced in Ref. [21] by requiring that the external power supplied to the
body through the applied tractions minus the rate of change of the free energy
and kinetic energy be non-negative, which reduces to the statement

$\int d^{d}x\mathbf{T}:\nabla\mathbf{v}\geq\frac{d{\cal F}}{dt},$ (10)

assuming balance of linear momentum is satisfied. Further, one also introduces
the fundamental kinematic relation of Field Dislocation Mechanics theory

$\dot{\bm{U}}=\nabla\mathbf{v}-\bm{\alpha}\times\mathbf{V},$ (11)

where $\mathbf{v}$ is the material velocity,
$\bm{\alpha}=\mbox{curl}\,\mathbf{U}$ and $\mathbf{V}$ is the dislocation
velocity (compare with Eq. (4), with a change of sign in $\bm{\alpha}$ due to
the sign convention in the definition of the Burgers vector). The dissipation
inequality and the kinematic constraint lead to constitutive equations for the
fields. The reversible part of the stress is given by

$\mathbf{T}=\partial{\cal
F}_{el}/\partial\bm{\epsilon}+C_{w}(\mathbf{U}-\mathbf{Q})_{sym}.$ (12)

The evolution equation for the phase field is given by

$\dot{\psi}=-L\left[\frac{\delta{\cal
F}_{sh}}{\delta\psi}+\frac{C_{w}}{2}\frac{\delta}{\delta\psi}\int
d^{d}x\;(\mathbf{U}-\mathbf{Q})_{sym}^{2}\right],$ (13)

where $\dot{\psi}$ is the material time derivative of $\psi$ and $L$ is a
kinetic mobility. The variational derivative in the right hand side of Eq.
(13) can be evaluated explicitly from the phase field as follows. We begin by
using the result [19]

$\sum_{n}q_{0\;i}^{(n)}q_{0\;j}^{(n)}=\frac{Nq_{0}^{2}}{d}\delta_{ij}$

where $N$ is a constant that depends on the lattice symmetry, $d$ is the space
dimension [19], and, as before, $q_{0}=1$. For example, $N=3$ for a two
dimensional hexagonal lattice, or $N=12$ for a three dimensional bcc lattice.
Therefore, the configurational distortion tensor is given by,

$P_{ij}=\frac{d}{N}\sum_{n}q_{0\;j}^{(n)}\partial_{i}\theta$

Since

$\nabla\theta_{n}=\frac{1}{|A_{n}|^{2}}{\rm Im}(A_{n}^{*}\nabla A_{n}),$

then ($Q=I-P$)

$Q_{ij}=-\frac{d}{N}\sum_{n}\frac{q_{0\;j}^{(n)}}{|A_{n}|^{2}}{\rm
Im}(A_{n}^{*}\nabla A_{n}).$

This is the relationship between the configurational distortion and the
amplitudes of the phase field. Finally, the required variational derivative
follows as,

$\delta{\cal I}=\delta\frac{1}{2}\int
d\mathbf{x}\;\left(U-Q\right)^{2}=-\sum_{n}\int
d\mathbf{x}\;\left[\frac{\delta I}{\delta A_{n}}\delta A_{n}+\frac{\delta
I}{\delta A_{n}^{*}}\delta A_{n}^{*}\right].$

with,

$\frac{\delta I}{\delta A_{n}}=(U_{ij}-Q_{ij})\frac{\partial Q_{ij}}{\partial
A_{n}}-\partial_{k}\left[(U_{ij}-Q_{ij})\frac{\partial
Q_{ij}}{\partial(\partial_{k}A_{n})}\right].$

The set of evolution statements is completed by balance of linear momentum
that controls the evolution of the material velocity (of course coupled to
(11) and (13)):

$\rho\,\dot{\mathbf{v}}=\mbox{div}\,\mathbf{T}$

or, in the quasi-static case,

$\mbox{div}\,\dot{\mathbf{T}}=\mathbf{0};$

in the latter case, the equation for the velocity is akin to Stokes flow with
source terms. Equation (11) is the statement of the partition of the
(compatible) material velocity gradient into the elastic strain rate and the
plastic strain rate, the latter arising from the transport of dislocations; in
geometrically exact finite deformation analysis, this is a consequence of
Burgers vector conservation (with a minimalistic assumption of setting a free
gradient to zero). This framework allows large permanent deformations with
minimal elastic energy cost, encoding the plastic strain produced due to
dislocation motion, to be an outcome of the theory.

Finally, we mention that an expression for the Peach-Kohler force also follows
directly from the dissipation inequality [21],

$\mathbf{V}=\mathbf{M}:\mathbf{X}\left(\mathbf{T}^{T}\bm{\alpha}\right);\qquad
V_{i}=M_{ij}\epsilon_{jrk}T_{sr}\alpha_{sk},$ (14)

where $\mathbf{M}$ is a (positive semi-definite) mobility tensor. This
expression illustrates that, in this model, dissipative motion of a
dislocation follows from elastic stress on the dislocation (the stress depends
on the defect configuration).

##  4 Outlook

The model described in Sec. 3.2 has not yet been validated against the
benchmark configuration involving the motion of a dislocation dipole as in
Sec. 3.1. Whereas in the quasi-static limit of $C_{w}$ large we expect a
similar evolution, this latter model has additional solutions that have not
yet been investigated. In particular, there will be additional propagating
modes that will mix with diffusive modes in the vicinity of the dislocation
cores.

It is also apparent from the foregoing discussion, that a derivation of the
phase field crystal model from a density functional theory, including stresses
as a slow variable, could help resolve the issues of independence, nature of
the currents, and their coupling. For example, while permeation is a
dissipative mode, the stress driving it originates entirely from the
reversible contribution to the stress from the free energy. Dissipation
appears indirectly through the dissipative relaxation of the order parameter,
but not directly through a dissipative current. The latter could, in fact, be
of lower order in gradients. Along the same lines, the stress currents in Sec.
3.1 are only dissipative until the additional compatible mode is introduced.
That the reversible distortion needs to be introduced through a unphysical
transformation (the displacement of the field by the compatible displacement)
is difficult to justify when, in general, reversible currents should follow
from considerations of Galilean invariance within the model.

Finally, experimental advances driven by increased fluxes at Synchrotron
facilities already allow the combination of near field high energy diffraction
microscopy (nf-HEDM) with Bragg coherent diffractive imaging [32, 33, 34, 35].
This technique has the potential to extend the range of length scales probed
by 3D X-ray diffraction techniques into the range of a few tens of nanometers,
and hence to image individual lattice defects and their motion in small
polycrystalline samples. These experiments will provide key information to
help elucidate many of the issues regarding the phase field crystal theory
discussed in this paper, and to make the technique predictive in terms of
defect motion at the nanoscale.

## Acknowledgments

We are indebted to Audun Skaugen for the preparation of the figures, and to
him, Vidar Skogvoll, and Marco Salvalaglio for many interesting and
stimulating conversations. The research of JV has been supported by the
National Science Foundation under Grant No. DMR-1838977, and the work of AA by
Grant No. OIA-DMR 2021019.

## References

  * [1] K. Elder, M. Katakowski, M. Haataja, and M. Grant, “Modeling elasticity in crystal growth,” _Phys. Rev. Lett._ , vol. 88, no. 24, p. 245701, 2002.
  * [2] H. Emmerich, H. Löwen, R. Wittkowski, T. Gruhn, G. I. Tóth, G. Tegze, and L. Gránásy, “Phase-field-crystal models for condensed matter dynamics on atomic length and diffusive time scales: an overview,” _Adv. in Physics_ , vol. 61, no. 6, pp. 665–743, 2012.
  * [3] K. Momeni, Y. Ji, K. Zhang, J. A. Robinson, and L.-Q. Chen, “Multiscale framework for simulation-guided growth of 2d materials,” _npj 2D Materials and Applications_ , vol. 2, no. 1, pp. 1–7, 2018.
  * [4] P. Stefanovic, M. Haataja, and N. Provatas, “Phase-field crystals with elastic interactions,” _Phys. Rev. Lett._ , vol. 96, no. 22, p. 225504, 2006.
  * [5] V. Heinonen, C. V. Achim, J. M. Kosterlitz, S.-C. Ying, J. Lowengrub, and T. Ala-Nissila, “Consistent hydrodynamics for phase field crystals,” _Phys. Rev. Lett._ , vol. 116, p. 024303, Jan 2016.
  * [6] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant, “Phase-field crystal modeling and classical density functional theory of freezing,” _Phys. Rev. B_ , vol. 75, no. 6, p. 064107, 2007.
  * [7] Z.-F. Huang, K. Elder, and N. Provatas, “Phase-field-crystal dynamics for binary systems: Derivation from dynamical density functional theory, amplitude equation formalism, and applications to alloy heterostructures,” _Phys. Rev. E_ , vol. 82, no. 2, p. 021605, 2010.
  * [8] A. J. Archer, D. J. Ratliff, A. M. Rucklidge, and P. Subramanian, “Deriving phase field crystal theory from dynamical density functional theory: Consequences of the approximations,” _Phys. Rev. E_ , vol. 100, p. 022140, Aug 2019.
  * [9] A. J. Archer, “Dynamical density functional theory for molecular and colloidal fluids: A microscopic approach to fluid mechanics,” _J. Chem. Phys._ , vol. 130, no. 1, p. 014509, 2009.
  * [10] A. Skaugen, L. Angheluta, and J. Viñals, “Dislocation dynamics and crystal plasticity in the phase-field crystal model,” _Phys. Rev. B_ , vol. 97, p. 054113, Feb 2018.
  * [11] V. Heinonen, C. V. Achim, K. R. Elder, S. Buyukdagli, and T. Ala-Nissila, “Phase-field-crystal models and mechanical equilibrium,” _Phys. Rev. E_ , vol. 89, p. 032411, Mar 2014.
  * [12] S. Majaniemi and M. Grant, “Dissipative phenomena and acoustic phonons in isothermal crystals: A density-functional theory study,” _Phys. Rev. B_ , vol. 75, no. 5, p. 054301, 2007.
  * [13] W. Zhou, J. Wang, Z. Wang, and Z.-F. Huang, “Mechanical relaxation and fracture of phase field crystals,” _Phys. Rev. E_ , vol. 99, p. 013302, Jan 2019.
  * [14] A. Kosevich, “Crystal dislocations and the theory of elasticity,” in _Dislocations in Solids_ , F. Nabarro, Ed. New York: North-Holland, 1979, vol. 1, p. 33.
  * [15] T. Mura, “Continuous distribution of moving dislocations,” _Philosophical Magazine_ , vol. 8, no. 89, pp. 843–857, 1963.
  * [16] N. Fox, “A continuum theory of dislocations for single crystals,” _IMA Journal of Applied Mathematics_ , vol. 2, no. 4, pp. 285–298, 1966.
  * [17] A. Acharya, “Microcanonical entropy and mesoscale dislocation mechanics and plasticity,” _Journal of Elasticity_ , vol. 104, no. 1, pp. 23–44, 2011\.
  * [18] ——, “Jump condition for gnd evolution as a constraint on slip transmission at grain boundaries,” _Phil. Mag._ , vol. 87, no. 8-9, pp. 1349–1359, 2007\.
  * [19] V. Skogvoll, L. Angheluta, A. Skaugen, M. Salvalaglio, and J. Viñals, “A phase crystal theory of the kinematics of dislocation lines,” _J. Mech. Phys. Solids_ , vol. 166, p. 104932, 2022.
  * [20] A. Skaugen, L. Angheluta, and J. Viñals, “Separation of elastic and plastic timescales in a phase field crystal model,” _Phys. Rev. Lett._ , vol. 121, no. 25, p. 255501, 2018.
  * [21] A. Acharya and J. Viñals, “Field dislocation mechanics and phase field crystal models,” _Phys. Rev. B_ , vol. 102, no. 6, p. 064109, 2020.
  * [22] J. Berry, M. Grant, and K. R. Elder, “Diffusive atomistic dynamics of edge dislocations in two dimensions,” _Phys. Rev. E_ , vol. 73, p. 031609, Mar 2006. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.73.031609
  * [23] D. Boyer and J. Viñals, “Weakly nonlinear theory of grain boundary motion in patterns with crystalline symmetry,” _Phys. Rev. Lett._ , vol. 89, no. 5, p. 055501, 2002.
  * [24] A. Roy and A. Acharya, “Finite element approximation of field dislocation mechanics,” _Journal of the Mechanics and Physics of Solids_ , vol. 53, no. 1, pp. 143–170, 2005.
  * [25] X. Zhang, A. Acharya, N. J. Walkington, and J. Bielak, “A single theory for some quasi-static, supersonic, atomic, and tectonic scale applications of dislocations,” _Journal of the Mechanics and Physics of Solids_ , vol. 84, pp. 145–195, 2015.
  * [26] R. Arora, X. Zhang, and A. Acharya, “Finite element approximation of finite deformation dislocation mechanics,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 367, p. 113076, 2020.
  * [27] A. Acharya, “An action for nonlinear dislocation dynamics,” _Journal of the Mechanics and Physics of Solids_ , vol. 161, p. 104811, 2022.
  * [28] P. Gumbsch and H. Gao, “Dislocations faster than the speed of sound,” _Science_ , vol. 283, no. 5404, pp. 965–968, 1999.
  * [29] C. Zhang, X. Zhang, A. Acharya, D. Golovaty, and N. Walkington, “A non-traditional view on the modeling of nematic disclination dynamics,” _Quarterly of Applied Mathematics_ , vol. 75, p. 309, 2016.
  * [30] V. Vitek, “Intrinsic stacking faults in body-centred cubic crystals,” _Philosophical Magazine_ , vol. 18, no. 154, pp. 773–786, 1968.
  * [31] C. Shen and Y. Wang, “Phase field model of dislocation networks,” _Acta materialia_ , vol. 51, no. 9, pp. 2595–2610, 2003.
  * [32] R. Suter, “Multiscale measurements for materials modeling,” _Science_ , vol. 356, no. 6339, pp. 704–705, 2017. [Online]. Available: http://science.sciencemag.org/content/356/6339/704
  * [33] J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” _Nature_ , vol. 400, pp. 342–344, 1999.
  * [34] J. Miao, T. Ishikawa, I. K. Robinson, and M. M. Murnane, “Beyond crystallography: Diffractive imaging using coherent x-ray light sources,” _Science_ , vol. 348, no. 6234, pp. 530–535, 2015.
  * [35] A. Yau, W. Cha, M. W. Kanan, G. B. Stephenson, and A. Ulvestad, “Bragg coherent diffractive imaging of single-grain defect dynamics in polycrystalline films,” _Science_ , vol. 356, no. 6339, pp. 739–742, 2017\. [Online]. Available: http://science.sciencemag.org/content/356/6339/739