# An ancient massive quiescent galaxy found in a gas-rich z $\sim 3$ group

Boris S. Kalita CEA, Irfu, DAp, AIM, Universitè Paris-Saclay, Universitè de
Paris, CNRS, F-91191 Gif-sur-Yvette, France  Emanuele Daddi CEA, Irfu, DAp,
AIM, Universitè Paris-Saclay, Universitè de Paris, CNRS, F-91191 Gif-sur-
Yvette, France  Chiara D’Eugenio CEA, Irfu, DAp, AIM, Universitè Paris-Saclay,
Universitè de Paris, CNRS, F-91191 Gif-sur-Yvette, France  Francesco Valentino
Cosmic Dawn Center (DAWN), Copenhagen, Denmark  Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, DK-2200, Copenhagen, Denmark  R.
Michael Rich Department of Physics & Astronomy, Univ. of California Los
Angeles, PAB 430 Portola Plaza, Los Angeles, CA 90095-1547, USA  Carlos Gómez-
Guijarro CEA, Irfu, DAp, AIM, Universitè Paris-Saclay, Universitè de Paris,
CNRS, F-91191 Gif-sur-Yvette, France  Rosemary T. Coogan Max-Planck-Institut
für Extraterrestrische Physik (MPE), Giessenbachstr.1, 85748 Garching, Germany
Ivan Delvecchio INAF—Osservatorio Astronomico di Brera, via Brera 28, I-20121,
Milano, Italy  David Elbaz CEA, Irfu, DAp, AIM, Universitè Paris-Saclay,
Universitè de Paris, CNRS, F-91191 Gif-sur-Yvette, France  James D. Neill 1
California Institute of Technology, 1200 East California Boulevard, MC 278-17,
Pasadena, CA 91125, USA  Annagrazia Puglisi Centre for Extragalactic
Astronomy, Department of Physics, Durham University, Durham, UK  Veronica
Strazzullo University of Trieste, Piazzale Europa, 1, 34127 Trieste TS, Italy
INAF—Osservatorio Astronomico di Brera, via Brera 28, I-20121, Milano, Italy
INAF - Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131, Trieste,
Italy

(Received June 01, 2021; Revised July 13, 2021; Accepted July 22, 2021)

###### Abstract

Deep ALMA and HST observations reveal the presence of a quenched massive
galaxy within the $z=2.91$ galaxy group RO-1001. With a mass-weighted stellar
age of $1.6\pm 0.4\,$Gyr this galaxy is one of the oldest known at $z\sim 3$,
implying that most of its $10^{11}\rm\,M_{\odot}$ of stars were rapidly formed
at $z>6$–8. This is a unique example of the predominantly passive evolution of
a galaxy over at least $3<z<6$ following its high-redshift quenching and a
smoking-gun event pointing to the early imprint of an age-environment
relation. At the same time, being in a dense group environment with extensive
cold-gas reservoirs as betrayed by a giant Ly$\alpha$ halo, the existence of
this galaxy demonstrates that gas accretion shutdown is not necessary for
quenching and its maintenance.

Galaxy evolution; High-redshift galaxies; Star formation; Galaxy Groups

††journal: ApJL††software: SExtractor (Bertin & Arnouts, 1996), GALFIT (Peng
et al., 2002, 2010), grizli (v0.4.0; Brammer, 2018), CASA (McMullin et al.,
2007), HYPERZ (Bolzonella et al., 2000)

##  1 Introduction

Hydro-dynamical simulations and semi-analytical models have found it difficult
to reproduce massive ($\gtrsim 10^{11}\,\rm M_{*}$) quiescent galaxies (QGs)
beyond $\rm z\gtrsim 3$ (Steinhardt et al., 2016; Schreiber et al., 2018b;
Cecchi et al., 2019), while observations have been pushing the redshift
boundary by detecting such systems at $\rm z\sim 3-4$ (Gobat et al., 2012;
Glazebrook et al., 2017; Schreiber et al., 2018; D’Eugenio et al., 2020a;
Valentino et al., 2020; Forrest et al., 2020a, b; Saracco et al., 2020). Most
of the systems so far observed belong to a population of recently quenched
‘post-starburst’ (PSB) galaxies with their last star-formation episode
occurring over the last $<0.8\,$Gyr. However, we might still be missing a
fraction of older ($\gtrsim 1\,$Gyr) QGs (D’Eugenio et al., 2020b; Forrest et
al., 2020b). Moreover, there has not been any insight on the influence of the
environment over such galaxies, as most of those studied are field objects.

Studying a younger sub-sample of a larger population of high-z QGs may have
influenced our understanding of how quiescence occurs and persists. Currently,
for star-formation to be suppressed in the $z>2$ epoch which usually features
gas-rich star-forming galaxies, there are a variety of possible channels.
These include merger driven starbursts (Puglisi et al., 2021), AGN feedback
(Brennan et al., 2018), gas strangulation (Peng et al., 2015), halo quenching
(Feldmann & Mayer, 2015), and morphological quenching (Martig et al., 2009).
Hence, further probing the high-z QG population, including those in dense
environments, will provide opportunities to investigate the earliest mechanics
of quenching.

With this in mind, we report the detection of an extremely red and old
quiescent galaxy within a dense environment at $z\sim 3$. Throughout, we adopt
a concordance $\Lambda$CDM cosmology, characterized by $\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$, and $H_{0}=70$ km s${}^{-1}\rm Mpc^{-1}$. We use a
Chabrier initial mass function. Magnitudes and colors are on the AB scale.

##  2 Galaxy-D and RO-1001

Figure 1: The core of RO-1001 as seen is HST/WFC3 (F606W, F125W, F160W). The
green contours denote the Ly$\alpha$ nebula luminosity with the locations of
the associated filaments possibly powered by cold-flow accretion (Daddi et
al., 2021). Galaxy-D is marked with a red square.

Three massive (log(M⋆/M⊙)$\gtrsim$11) galaxies (A, B, and C in Fig. 1 and
Daddi et al., 2021) in the $z=2.91$ RO-1001 group have a combined star
formation rate (SFR) $\sim 1250\,\rm M_{\odot}\,yr^{-1}$. NOEMA observations
of the CO(3-2) transition showed that they lie at the same redshift of
$z=2.91$. Furthermore, the presence of extensive cold gas reservoirs is
suggested by Ly$\alpha$ observations that revealed extended filaments
converging onto a single massive halo with a combined luminosity of $1.3\pm
0.2\times 10^{44}$ erg s-1, most likely still affected by cold-accretion
persisting in this system with estimated dark matter halo mass of $4\times
10^{13}M_{\odot}$. These conclusions are based on the observed line profile,
the velocity field, and the energetics underlying the Lya emission (Daddi et
al., 2021). Moreover, also pointing to the availability of cold gas is the
redshift of the group, which places it in the epoch where cold flow accretion
is expected to dominate (Valentino et al., 2015; Overzier, 2016). Galaxy-D was
reported earlier as a possible passive member of RO-1001 (Fig. 1; Daddi et
al., 2021).

##  3 Observations

###  3.1 Optical and Near-IR imaging

RO-1001 was observed with HST/WFC3 imaging in 3 bands (F160W, F125W and F606W)
over a total of 11 orbits during Cycle 27 (Proposal ID: 15190, PI: E. Daddi).
The data reduction was executed using the pipeline grizli
111https://github.com/gbrammer/grizli. The $5\sigma$ point-source
sensitivities reached are 26.25 (F160W), 26.47 (F125W) and 26.39 (F606W) with
a pixel scale of $0.06^{\prime\prime}$ and a half-power beam-width of
$0.24^{\prime\prime}$ for F160W. Public COSMOS F814W imaging was also
incorporated into the analysis. Due to a lack of coverage in z and y-bands, we
include Subaru Hyper Suprime-Cam images were from the “COSMOS2015” database
(Laigle et al., 2016). Furthermore, we used the Ks band image from data
release 4 of Ultra-VISTA. Finally, IRAC $3.6\,\mu$m and $4.5\,\mu$m images
were taken from the Spitzer Matching Survey of the Ultra-VISTA Deep Stripes
(SMUVS), while those at $5.8\,\mu$m were taken from the COSMOS2015 database.

###  3.2 Sub-mm imaging

To have an estimate of the obscured star formation down to very faint levels
in RO-1001, and as a result on Galaxy-D, we also obtained sub-mm data. Atacama
Large Millimetre Array (ALMA) band 7 observations were taken in Cycle 7
(Project ID: 2019.1.00399.S, PI: R.M. Rich). The data reduction was mainly
carried out using the Common Astronomy Software Application (CASA). The final
mosaics have a maximum sensitivity of $\rm 28\,\mu Jy\,beam^{-1}$, with a
synthesised beam size of $0.49^{\prime\prime}\times 0.46^{\prime\prime}$.

##  4 Analysis and results

###  4.1 Photometry

Figure 2: Panels showing Galaxy-D (at the centre) along with its surrounding
foreground galaxies in 3 out of 5 wavelength bands for which GALFIT was used
to extract the photometry. As detailed in Sec. 4.1, the results of H-band
($\rm 1.6\,\mu m$) were used to deblend at the other bands ($\rm 2.2\,\mu
m,3.6\,\mu m,4.5\,\mu m,5.8\,\mu m$, the first two of which are shown in the
middle and lower rows respectively). More sources were fit and subtracted at
longer wavelengths due to the increase in PSF size. The residual images are
found to not have any fluctuations $>1\sigma$.

In the COSMOS2015 catalogue (Laigle et al., 2016) Galaxy-D is blended with a
galaxy located 1′′ to the North (Fig. 1). Hence, we make new flux measurements
for it in each of the wavelength windows we have access to. We first use
SExtractor for the 4 HST/WFC3 images since the the resolution is sufficient to
clearly separate galaxy-D from its neighboring foreground galaxy. We use F160W
as the detection image to carry out matched aperture photometry using
isophotal radii. To ensure equivalent resolutions, we also convolve the F606W,
F814W and F125W images with Gaussians such that their respective final point-
spread-functions (PSF) have the same size as that in F160W. From the catalogue
hence created, we find it to be undetected in F606W and F814W (which therefore
provide upper-limits) while F125W and F160W have clear detections with signal-
to-noise (S/N) of $\sim 25$ and $65$ respectively. The z and y-bands also
return non-detections, which are hence used to set flux upper-limits.

In order to ensure that all the flux from Galaxy-D is accounted for, we use
galaxy morphological model fitting with the software GALFIT. Since the
$1.6\,\mu$m (F160W) image has the highest S/N, we use it for this attempt. We
fit sérsic profiles to both Galaxy-D and the foreground galaxy, along with
Gaussian profiles to multiple surrounding objects (Fig. 2). We measure a total
AB-magnitude of 24.00 for Galaxy-D and get a sérsic index of $2.0\pm 0.2$. The
error on the latter value includes the results from a co-addition of the F125W
image to the F160W image to improve S/N as well as an independent analysis of
the F125W image, all of which are in agreement. The sérsic index suggests that
this galaxy might feature a disky component with a significant bulge
component, within the scatter expected from high-z QGs (Lustig et al., 2021).
The galaxy is also found to be extremely compact with an effective radius,
$\rm r_{e}=0.14\pm 0.03^{\prime\prime}$, which at z=2.9 is $1.1\pm 0.1\,$kpc.
We further use the flux difference between the GALFIT measurement and that
from SExtractor, to scale up the catalogue (isophotal) flux measurement at
$1.25\,\mu$m (F125W) as well as the upper-limits at $0.6\,\mu$m (F606W) and
$0.8\,\mu$m (F814W).

The morphological properties of Galaxy-D determined through the model fitting
at $1.6\,\mu$m are further used as priors in Ks band ($2.1\,\mu$m) as well as
the IRAC $3.6\,\mu$m, $4.5\,\mu$m and $5.8\,\mu$m. This is done in order to
deblend the galaxy from its foreground neighbors, which are also
simultaneously fit with models based on $1.6\,\mu$m priors (Fig. 2). Images
beyond this wavelength window do not show a detection, while the implied upper
limits are not constraining, and are hence not used.

###  4.2 SED fitting

Figure 3: (Top-left) The best-fit SED template using the available photometry
for Galaxy-D. The red bars show the observed photometry while the cyan points
mark the expected values from the SED. The blue dashed-lines mark the
locations of the Balmer and 4000 Å breaks. (Top-right) The probability
distribution function for the photometric redshift of Galaxy-D. The red dashed
line marks the spectroscopic redshift of RO-1001 while the blue shaded region
marks the $90\%$ confidence region of the photometric redshift for Galaxy-D.
(Bottom) The UVJ colors of high-z QGs ($\rm z>2.5$) adapted from D’Eugenio et
al. (2020a). The grey arrow shows the expected color evolution with increasing
age. Galaxy-D is marked with a star, clearly outside the PSB region specified
in Belli et al. (2019). The marker sizes are proportional to the stellar mass.

We fit the photometry derived in the previous section with BC03 stellar
population models (Bruzual & Charlot, 2003) to derive the properties of
Galaxy-D (Fig. 3, top-left). We use HYPERZ to determine the best fitting
spectral energy distribution (SED) templates through a $\chi^{2}$ minimisation
procedure.

We investigate a redshift range of $0-6$ with a step size of 0.01. For the
BC03 templates, a range of metallicities are allowed
($\rm\frac{1}{5}Z_{\odot}$, $\rm\frac{1}{2}Z_{\odot}$, $\rm Z_{\odot}$ and
$\rm\frac{5}{2}\,Z_{\odot}$). We also implement the Calzetti et al. (2000)
along with SMC and LMC laws for dust attenuation, with extinction values,
A${}_{\rm V}=0-7$ in steps of 0.025. Finally, a wide variety of star formation
histories (SFH) are used:

  1. 1.

Simple stellar population (SSP): This is the most basic model with all of the
stellar mass formed during an infinitesimal short burst of star formation at
$t=0$.

  2. 2.

Constant star formation models: Although we have enough evidence for the
quiescence of the galaxy (discussed in detail later in Sec. 5.1), we use this
model featuring no decline in star formation (it remains constant throughout)
to check if a high attenuation star-forming scenario is formally consistent
with the available photometry.

  3. 3.

Truncated models: These involve a constant SFR from $t=0$ up until a
$t_{stop}$, beyond which $\rm SFR=0\,M_{\odot}\,yr^{-1}$. The $t_{stop}$ is
varied from 0.1 Gyr to a generous upper limit of $5.0\,\rm Gyr$ in steps of
0.1 dex.

  4. 4.

Delayed $\tau$-models: Finally, we use an SFH that has been widely preferred
to model QGs. The exponentially declining SFH here is characterised as
$\propto(t/\tau^{2})\,e^{-t/\tau}$ with a peak of star formation at $t=\tau$.
The $\tau$ varies within [0.1, 5.0] Gyr with steps of 0.05 dex.

We investigate the reduced $\chi^{2}$ values from the fits using the various
BC03 models, over the above mentioned grid of AV. We get the best fit with the
delayed $\tau$-models (with a reduced $\chi^{2}\sim 0.4$), closely followed by
SSP and truncated models. In comparison, we find a complete lack of fit using
the constant star formation templates ($\Delta\chi^{2}>30$), with any level of
reddening, reinforcing our claim of Galaxy-D being quiescent.

###  4.3 Galaxy Parameters

The photometric redshift (photo-z) determined from the SED fitting is found to
be $2.9\pm 0.1$ (Fig. 3, top-right), with the limits giving the $1\sigma$
confidence interval ($\Delta\chi^{2}<1$). We also derive a stellar mass,
log($\rm M_{*}/M_{\odot})=11.0\pm 0.2$, in agreement with the previous results
in Daddi et al. (2021). Moreover, we find minimum levels of attenuation with
the fitting procedure returning a value of $\rm A_{V}\sim 0.1$. We further
confirm that our choice of gridding and ‘edge-effects’ do not influence our
results by investigating a variety of grids for $\rm A_{V}$ while ensuring
that the parameter space on both sides of the best fit value is well sampled.

Following the analysis in similar studies (Gobat et al., 2012; Schreiber et
al., 2018b; Valentino et al., 2020; D’Eugenio et al., 2020a), we define
$t_{50}$, which is the time elapsed since the epoch of “half-mass formation”
up until the time of observation. We determine this quantity to be $1.6\pm
0.4\,\rm Gyr$ at 90% confidence, after marginalising on all the unknown
parameters (SFH, metallicity, dust attenuation law and $A_{V}$, redshift).
Fixing metallicity to solar does not change our results. Studying each of the
values, we find that $t_{50}$ gets progressively lower as we go from sub-solar
to super-solar metallicities. At $z=2.9$, for $\rm\frac{1}{5}Z_{\odot}$, we
derive a best-fit $t_{50}\sim 2.2\,\rm Gyr$ which is almost equal to the age
of the universe at that redshift, while for the super-solar
$\rm\frac{5}{2}Z_{\odot}$, the $t_{50}$ is found to be $\sim 1.0\,\rm Gyr$.
Although these results incorporate all previously mentioned SFH models
together, using each model separately we find that the results from
delayed-$\tau$ and truncated models are in agreement, while those using SSP
(which can be considered as the light-weighted age) push the age upper-limit
by $\sim 0.2\,$Gyr. Finally, fixing the redshift to that of RO-1001 ($\rm
z=2.91$) does not change the constrain on age estimate.

We also investigate the robustness of our measurements with simulations. We
perturb the flux measurements with values drawn within Gaussians with sigma
equal to the $1\sigma$ uncertainty of the respective flux measurement to
generate 1000 sets of photometric datapoints. Besides a consistent average
best-fit age of $1.7\,$Gyr, we find the $1\sigma$ scatter to be $0.15\,$Gyr
(once accounting for the gridding effect). This uncertainty is slightly lower
than the $0.2\,$Gyr ($1\sigma$ or $68\%$ confidence interval; Fig. 4, left)
that we get from the $\Delta\chi^{2}$ analysis. We repeated the exercise
perturbing the best fit SED photometry instead of the observed photometry,
finding the same results. Furthermore, for both cases, we also check the
possible consequences of asymmetry in the range of acceptable ages in each fit
($\Delta\chi^{2}<1$) by taking an average of their values. This also returns
an age of $1.7\,$Gyr. Hence, we do not revise our more conservative results
from the latter technique quoted throughout this work.

Furthermore, in order to have a SFH model-independent assessment of Galaxy-D
within the framework of QGs at high-z, we also plot the (UVJ)rest colors (Fig.
3, bottom). These are determined using the F160W, Ks and IRAC $4.5\,\mu$m
photometry, that match nearly perfectly to the rest-frame U, V and J bands,
with small correction for the residual differences in effective wavelengths
adopted from the best fitting SED.

Figure 4: (Left) The variation of $\Delta\chi^{2}$ in the stellar population
age ($\rm t_{50}$) vs redshift space. The dot-dash lines demarcate the
lookback times for various redshifts along with the thickest of them showing
the age of the Universe, all of which are functions of redshift (Right) The
best fit delayed $\tau$-model SFH ($\tau=0.2\,$Gyr) as a function of age of
the universe as well as redshift. Also shown are the SFH models with
$\tau=0.1,0.3\,$Gyr which are within the 90% confidence interval. Finally, the
SFR $2\sigma$ upper-limit of $\rm 13\,M_{\odot}\,yr^{-1}$ derived from the
ALMA $870\,\mu$m non-detection is shown with the black downward arrow. Figure
5: The stellar population age ($\rm t_{50}$) vs redshift distribution for
log(M${}_{*}/\rm M_{\odot})>10.5$ QGs adapted from D’Eugenio et al. (2020b),
with the references therein. Although the definition of the ‘mass-weighted
ages’ for a subset of the works add a level of ambiguity, they are found to be
in agreement with our definition of $\rm t_{50}$ within a maximum of
$0.1\,$Gyr (although in case of Morishita et al. (2019), their ‘mass-weighted
age’ is equal to $\rm t_{50}$). We add Galaxy-D (light-blue star) with the
$90\%$ confidence intervals shown as the errorbars for both axes.

###  4.4 Uniqueness and RO-1001 membership

The photo-z measurement of $2.9\pm 0.1$ which is in remarkable agreement with
the z${}_{spec}=2.91$ of RO-1001, strongly suggests group membership. This is
further reinforced by the vanishingly low probability of $\sim 3\times
10^{-5}$ for a chance alignment of Galaxy-D within $10^{\prime\prime}$ from
the center of the RO-1001 group (as defined by the peak of the Ly$\alpha$
emission). This is estimated by the implied surface density of the 20 galaxies
with H-band magnitudes $\leq 24.0$, J-H color $\geq 1.1$ and H-Ks color $\geq
2.2$ over the whole COSMOS area (Laigle et al., 2016). The color-selection
thresholds determined by the photometry of Galaxy-D also incorporate the
uncertainties in the measurements. However, 16/20 of these galaxies are
consistent with being selected due to the $1\sigma$ scatter in their color
measurements (recall that the COSMOS catalog contains more than half million
sources). Hence, this number of objects found can be considered as a generous
upper-limit. The lack of deep HST data like those we have for Galaxy-D
primarily influences this conclusion. Out of the remaining 4/20, two lack
photometry in most of the COSMOS observation bands and therefore do not have a
photo-z reported in the catalogue. The final two galaxies (Laigle IDs: 749600
and 824982 222RA, dec: (149.47703, 2.43692) and (149.52635, 2.55155)) are the
only objects in the 2-deg-square area of COSMOS, which have been robustly
found to have similar observed characteristics within the framework of the
state-of-the-art multi-wavelength database available for the field. Also as a
consequence, these galaxies have similar redshift, mass and age estimates as
for Galaxy-D. However, due to a lack of sub-mm data (and therefore a
confirmation of quiescence) along with any environmental information, a
thorough study for these objects cannot be currently undertaken. However, this
lack of counterparts should not be regarded as a suggestion of similar objects
being absent. Rather, this could simply point to the inefficiency of current
surveys to detect such galaxies with a high completeness, making Galaxy-D
extremely rare within the current detection limits of surveys.

Hence, besides being very unique, it is extremely unlikely that this galaxy is
at a redshift different from that of the group RO-1001. Although, it cannot be
currently determined with certainty whether Galaxy-D is exactly within the
core of RO-1001 or part of the associated large-scale structure. In either
case Galaxy-D can still be considered to be within a dense environment.

##  5 Discussion

###  5.1 Quiescence and the last epoch of star-formation

Confirmation of quiescence at such high redshifts can be challenging
(Glazebrook et al., 2017; Schreiber et al., 2017; Simpson et al., 2017).
However for galaxy-D, the optical and NIR colors suggest a lack of star-
formation based on a complete disagreement with the star-forming SFH models
(Sec. 4.2). Taking into account all possible SFH models (within 90% confidence
interval), we determine a SFR upper-limit of $4\,\rm M_{\odot}/yr$, $\times
30$ below the MS.

The non-detection in our deep ALMA $870\,\mu$m data (Sec. 3.2) is used to
determine a $2\sigma$ upper-limit SFR (conservative; fully accounting for its
spatial extent of NIR emission of the galaxy based on its Sérsic profile) of
$13\,\rm M_{\odot}/yr$ assuming conservatively a main-sequence template
appropriate for $z=3$ (Béthermin et al., 2012). This is already a factor $\sim
10$ below the main sequence at $\rm z\sim 3$ (Santini et al., 2017). The limit
would be a further factor of 2 deeper if assuming instead a colder IR SED
following Gobat et al. (2018), i.e. at $\times 20$ below the MS. Hence, it can
be unambiguously concluded that Galaxy-D is truly quiescent, consistently from
both the optical and IR side.

Based on the mass-weighted ($t_{50}$) as well as light-weighted age limits
(Sec. 4.3), the primary episode of star-formation had already elapsed by
$z\gtrsim 6$ and most likely by $z\sim 8$ (Fig. 4, left), right after which it
likely experienced a rapid decline in its SFR characterised by the best-fit
delayed $\tau$-model with $\tau=0.2\,$Gyr (Fig. 4, right). However, this
should be considered as an approximate upper limit, since $\tau=0.1\,$Gyr and
the SSP (which has $\tau$ equivalent to 0) also return acceptable fits (within
an overall $\Delta\chi^{2}<1$ or 68% confidence interval).

This form of rapid quenching has been extensively proposed to explain
populations of QGs at $\rm z\sim 3-4$ (e.g., Schreiber et al., 2018b; Forrest
et al., 2020a, b; Valentino et al., 2020; D’Eugenio et al., 2020a, b).
However, these are usually found to be PSBs, with $t_{50}<0.8\,$Gyr (Fig. 5).
Such studies usually leave out galaxies similar to Galaxy-D that are expected
to have an older ($>1\,$Gyr) stellar population, due to preferential selection
of brighter targets. Forrest et al. (2020b) shows that galaxies with such ages
would be considerably less bright (by $\sim 1-2$ magnitudes) in near-IR,
compared to PSBs. For example, in comparison to Galaxy-D, the Valentino et al.
(2020) and D’Eugenio et al. (2020b) samples in Ks band are on average brighter
by 1.7 and 2.4 magnitudes, respectively. This is also the case for a very
recent detection of a $\sim 0.6\,\rm Gyr$ old massive QG detected in a
protocluster environment at $\rm z=3.1$, with a Ks band magnitude 1.3 brighter
(Kubo et al., 2021). This difference is also found to manifest in the SFH-
independent UVJ diagram Fig. 3 (bottom), which clearly shows a marked
separation in color and hence an inferred age between the well studied PSB
population and Galaxy-D.

Hence, spectroscopic analysis of such very-old galaxies is much more
difficult, with requirements of extremely long integration times. Possibly as
a result, similarly old galaxies at $z\sim 3$ have not yet been
spectroscopically studied yet (Fig. 5). We do acknowledge the presence of a
population of galaxies at $z\sim 2-2.7$ in Fig. 5 with seemingly similar ages
from Morishita et al. (2019). But this comparison is affected by systematics,
as their age estimates from bursty star-formation histories would be lowered
by a factor of $\sim 1.5-2.0$ if smooth SFH models are employed as in our
study (see Fig. 12 in Morishita et al., 2019), further emphasizing the unicity
of Galaxy-D in the panorama of known old galaxies at $z\sim 3$.

Larger photometric studies, also with higher uncertainties regarding galaxies
being truly quiescent, primarily feature young stellar populations with an
extremely small fraction of their samples having ages $>1\,\rm Gyr$ (e.g.,
Straatman et al., 2014; Carnall et al., 2020). We do have examples of galaxies
at lower redshifts $z\sim 0.5$–2 with old enough stellar populations
suggesting that they should be already $>1$ Gyr old at $z>3$ (Tacchella et
al., 2021, and references therein). If those ages inferences are correct (it
is increasingly difficult at lower-z to accurately project star formation
histories at the earliest times), it is still possible that the progenitors of
the lower-z oldest galaxies were not yet hierarchically assembled into a
single galaxy by $z\sim 3$ (Mancini et al., 2019). It is hence unclear whether
old massive galaxies at $z\sim 3$ are rare due to them experiencing re-
accretion of gas and returning back to a star-forming population, to not
having yet assembled, or they have simply not yet been detected at these
redshifts. Future time-intensive studies would hence be required to complete
the mapping of high-z QGs in this regime. A crucial diagnostic for such
attempts would be the Balmer and $4000\,$Å breaks characterised by DB (Kriek
et al., 2006) and Dn4000 (Balogh et al., 1999) respectively. The transition of
the DB/Dn4000 ratio towards values $<1$ is a tell-tale sign of ages $\gtrsim
1.0\,$Gyr (D’Eugenio et al., 2020b). As an example, we calculate DB/Dn4000 =
0.94 from the best-fit SED template (Fig. 3, top-left), suggesting a small
value for DB. This is in stark contrast to the PSB population, which is
usually characterised by much stronger Balmer breaks and hence higher
DB/Dn4000 ratios (e.g., D’Eugenio et al., 2020b, find an average value of 1.53
for their sample). Future spectroscopy, that should be within reach of JWST,
will be needed to validate these inferences.

###  5.2 Tracing the evolution

The short characteristic time-scales (Fig. 4) and compact nature of Galaxy-D
($\rm r_{e}=1.1\pm 0.1\,kpc$ at $z=2.9$) can be explained by scenarios like
mergers (e.g., Puglisi et al., 2021) that drive gas into a compact central
region of the galaxy, leading to rapid consumption of gas through star
formation. This quickly exhausts the available gas, thereby quenching the
system and leaving behind a compact QG. Compact star-forming cores possibly
related to mergers are already observed in the three primary star-forming
galaxies in RO-1001 which are nearly identical in mass and size to galaxy-D
(Daddi et al., 2021, Kalita, B. S. et al. in prep.), but forming over 1000
$\rm M_{\odot}$ yr-1 in stars complessively.

Furthermore, with 50% stellar mass ($\sim 5\times 10^{10}\,\rm M_{\odot}$)
being formed at $\rm z\gtrsim 6$, we can determine whether its progenitor
could have existed based on the expected galaxy mass-function, mainly
sensitive to star-forming galaxies, at such redshifts. To make a comparison,
we estimate a log number density of $-7.2\pm 0.5\,$Mpc-3 by assuming a single
detection within $2.6<\rm z<3.2$ over the COSMOS $1.7\,\rm deg\times 1.7\,deg$
area. Then a number density of $>5\times 10^{10}\,\rm M_{\odot}$ possible
progenitors at $\rm z>5.5$ is derived from the results of Grazian et al.
(2015). We find the Galaxy-D-based density to be almost an order of magnitude
lower than the latter value determined from the galaxy mass-function.
Including the two candidate ‘ancient’ QGs discussed in Sec. 4.4 brings this
down to a factor $\sim 3$, with the possibility of further decrease if a few
of rest are found to also have similar properties. Nevertheless, it cannot yet
be ascertained whether more galaxies with formation epochs at $\rm z\gtrsim 6$
within the COSMOS field do actually exist at similar redshifts. We can
conclude however, that at least $10\%$ of the $z>6$ star-forming galaxies with
stellar mass $>5\times 10^{10}\,\rm M_{\odot}$ can be expected to have
undergone quenching by $z\sim 6$ to form old QGs like Galaxy-D by $\rm z\sim
3$.

###  5.3 Quiescence in a dense environment

Given the estimated $t_{50}$, not only we need a quenching mechanism that
could form it very early, but any galaxy rejuvenation has to be prevented
thereafter. This primarily involves curbing accretion of gas or preventing the
gas from forming stars. Galaxy-D belongs to RO-1001 or associated sub-halos at
$z=2.91$ (Sec. 4.4). This halo is expected to be fed by copious cold gas
accretion (Daddi et al., 2021) still at $z=2.91$, all the more in its past
evolutionary and assembly phases and up to the $z>6$ quenching of Galaxy-D
that inescapably happened in a progenitor halo in which accretion must have
been prominent. This shows that one may have quenching and keep fully quenched
galaxies (that remain so for over 1 Gyr at $z>3$–6) even with large amounts of
diffuse cold gas and in presence of cold accretion. This is evidence that
environmental quenching scenarios like cosmological starvation (Feldmann &
Mayer, 2015) and gas strangulation (Peng et al., 2015) do not play a necessary
role in quenching and/or maintaining the quiescence of high-z galaxies.

One of the alternative channels for inhibiting star formation would be AGN-
driven feedback (Sanders et al., 1988; Di Matteo et al., 2005; Hopkins et al.,
2006; McCarthy et al., 2011). Although Galaxy-D is not hosting a detectable
AGN (Daddi et al., 2021), multiple works have suggested that there are
enhanced AGN fractions in high-z massive QGs (e.g., Olsen et al., 2013; Aird
et al., 2019). Hence, it is highly likely that this galaxy would have at least
experienced episodes of AGN activity in the past, inevitable given its large
stellar mass. These AGN episodes might have or not curbed star formation as in
QSO-mode quenching. However, AGN radio-mode activity was most likely crucial
to allow it to evolve passively for $>1$ Gyr, as supported by the high radio
detection rate of distant quenched systems (Gobat et al., 2018; D’Eugenio et
al., 2020b). Furthermore, other scenarios like morphological quenching
(stellar bulges preventing the collapse of gas for star formation; Martig et
al., 2009) could also be at play to prevent any significant late-time star
formation.

But a question still remains: what sets Galaxy-D apart from the other
similarly massive but intensively star-forming galaxies in RO-1001? It appears
to be simply at a later evolutionary phase, while its vigorously star-forming
counterparts have grown their large stellar masses later and are still
undergoing rapid growth. They might have an eventual fate as that of Galaxy-D
(similar to the scenario outlined in Puglisi et al., 2021), and be hardly
distinguishable in the stellar population properties once reaching $z\sim 0$
more than 10 Gyr later.

While Galaxy-D has a mass-weighted age ($t_{50}$) of $1.6\pm 0.4\,$Gyr, most
other quiescent galaxies at similar redshifts or above are younger PSBs and
can predominantly be characterised as field galaxies (due to a lack of
association with known overdensities). Therefore, the age difference between
Galaxy-D and the currently known PSB population $\sim 1.0\,\rm Gyr$ is one of
the first known evidences of an age-environment relation at high-z.
Interestingly, similar difference of 1 Gyr between clusters and field have
also been found for local populations, where this relation has already been
established (e.g., Bernardi et al., 2006; Gallazzi et al., 2021). But unlike
at low-z, the age offset at $\rm z\gtrsim 3$ is a sizable fraction of the
stellar ages themselves, thereby manifesting as much more pronounced
observational differences as has already been discussed in Sec. 5.1.

In any case, Galaxy-D satisfying such a relation suggests a rapid hierarchical
buildup of mass in dense environments, earlier than in the field. This in turn
leads to the quenching scenarios discussed in the previous sections, mostly
driven by mass and internal processes, while still agreeing with our
conclusions of direct environment-dependent quenching scenarios not being
necessary.

_Acknowledgements:_ We would like to express our gratitude to Raphael Gobat
for his valuable inputs. We also thank Gabriel Brammer for assistance with the
use of grizli for the HST data reduction. R.M.R. and J.D.N. acknowledge
support from GO15910. V.S. acknowledges the support from the ERC-StG
ClustersXCosmo grant agreement 716762. F.V. acknowledges support from the
Carlsberg Foundation Research Grant CF18-0388 “Galaxies: Rise and Death” and
from the Cosmic Dawn Center of Excellence funded by the Danish National
Research Foundation under then Grant No. 140. A.P. gratefully acknowledges
financial support from STFC through grants ST/T000244/1 and ST/P000541/1.
B.S.K. would like to thank William G. Hartley (ETH Zurich, Switzerland), for
providing crucial software-related files that were needed for this work. This
paper makes use of the following ALMA data: 2019.1.00399. ALMA is a
partnership of ESO (representing its member states), NSF (USA) and NINS
(Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ.

## References

  * Aird et al. (2019) Aird, J., Coil, A. L., & Georgakakis, A. 2019, MNRAS, 484, 4360. doi:10.1093/mnras/stz125
  * Avni (1976) Avni, Y. 1976, ApJ, 210, 642. doi:10.1086/154870
  * Balogh et al. (1999) Balogh, M. L., Morris, S. L., Yee, H. K. C., et al. 1999, ApJ, 527, 54. doi:10.1086/308056
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
  * Belli et al. (2019) Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17. doi:10.3847/1538-4357/ab07af
  * Bernardi et al. (2006) Bernardi, M., Nichol, R. C., Sheth, R. K., et al. 2006, AJ, 131, 1288. doi:10.1086/499522
  * Béthermin et al. (2012) Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23. doi:10.1088/2041-8205/757/2/L23
  * Bolzonella et al. (2000) Bolzonella, M., Miralles, J.-M., & Pelló, R. 2000, A&A, 363, 476
  * Brammer (2018) Brammer, G. 2018, Zenodo
  * Brennan et al. (2018) Brennan, R., Choi, E., Somerville, R. S., et al. 2018, ApJ, 860, 14. doi:10.3847/1538-4357/aac2c4
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000. doi:10.1046/j.1365-8711.2003.06897.x
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682. doi:10.1086/308692
  * Carnall et al. (2020) Carnall, A. C., Walker, S., McLure, R. J., et al. 2020, MNRAS, 496, 695. doi:10.1093/mnras/staa1535
  * Cecchi et al. (2019) Cecchi, R., Bolzonella, M., Cimatti, A., et al. 2019, ApJ, 880, L14. doi:10.3847/2041-8213/ab2c80
  * Cimatti et al. (2006) Cimatti, A., Daddi, E., & Renzini, A. 2006, A&A, 453, L29. doi:10.1051/0004-6361:20065155
  * Daddi et al. (2021) Daddi, E., Valentino, F., Rich, R. M., et al. 2021, A&A, 649, A78. doi:10.1051/0004-6361/202038700
  * D’Eugenio et al. (2020a) D’Eugenio, C., Daddi, E., Gobat, R., et al. 2020, ApJ, 892, L2. doi:10.3847/2041-8213/ab7a96
  * D’Eugenio et al. (2020b) D’Eugenio, C., Daddi, E., Gobat, R., et al. 2020, arXiv:2012.02767
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604. doi:10.1038/nature03335
  * Feldmann & Mayer (2015) Feldmann, R. & Mayer, L. 2015, MNRAS, 446, 1939. doi:10.1093/mnras/stu2207
  * Forrest et al. (2020a) Forrest, B., Annunziatella, M., Wilson, G., et al. 2020, ApJ, 890, L1. doi:10.3847/2041-8213/ab5b9f
  * Forrest et al. (2020b) Forrest, B., Marsan, Z. C., Annunziatella, M., et al. 2020, ApJ, 903, 47. doi:10.3847/1538-4357/abb819
  * Gallazzi et al. (2021) Gallazzi, A. R., Pasquali, A., Zibetti, S., et al. 2021, MNRAS, 502, 4457. doi:10.1093/mnras/stab265
  * Glazebrook et al. (2017) Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Nature, 544, 71. doi:10.1038/nature21680
  * Gobat et al. (2012) Gobat, R., Strazzullo, V., Daddi, E., et al. 2012, ApJ, 759, L44. doi:10.1088/2041-8205/759/2/L44
  * Gobat et al. (2018) Gobat, R., Daddi, E., Magdis, G., et al. 2018, Nature Astronomy, 2, 239. doi:10.1038/s41550-017-0352-5
  * Grazian et al. (2015) Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96. doi:10.1051/0004-6361/201424750
  * Hopkins et al. (2006) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1. doi:10.1086/499298
  * Kriek et al. (2006) Kriek, M., van Dokkum, P. G., Franx, M., et al. 2006, ApJ, 645, 44. doi:10.1086/504103
  * Kubo et al. (2021) Kubo, M., Umehata, H., Matsuda, Y., et al. 2021, arXiv:2106.10798
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24. doi:10.3847/0067-0049/224/2/24
  * Lustig et al. (2021) Lustig, P., Strazzullo, V., D’Eugenio, C., et al. 2021, MNRAS, 501, 2659. doi:10.1093/mnras/staa3766
  * Mancini et al. (2019) Mancini, C., Daddi, E., Juneau, S., et al. 2019, MNRAS, 489, 1265. doi:10.1093/mnras/stz2130
  * Martig et al. (2009) Martig, M., Bournaud, F., Teyssier, R., et al. 2009, ApJ, 707, 250. doi:10.1088/0004-637X/707/1/250
  * McCarthy et al. (2011) McCarthy, I. G., Schaye, J., Bower, R. G., et al. 2011, MNRAS, 412, 1965. doi:10.1111/j.1365-2966.2010.18033.x
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, Astronomical Data Analysis Software and Systems XVI, 376, 127
  * Morishita et al. (2019) Morishita, T., Abramson, L. E., Treu, T., et al. 2019, ApJ, 877, 141. doi:10.3847/1538-4357/ab1d53
  * Olsen et al. (2013) Olsen, K. P., Rasmussen, J., Toft, S., et al. 2013, ApJ, 764, 4. doi:10.1088/0004-637X/764/1/4
  * Overzier (2016) Overzier, R. A. 2016, A&A Rev., 24, 14. doi:10.1007/s00159-016-0100-3
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., et al. 2002, AJ, 124, 266. doi:10.1086/340952
  * Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., et al. 2010, AJ, 139, 2097. doi:10.1088/0004-6256/139/6/2097
  * Peng et al. (2015) Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192. doi:10.1038/nature14439
  * Puglisi et al. (2021) Puglisi, A., Daddi, E., Brusa, M., et al. 2021, Nature Astronomy, 5, 319. doi:10.1038/s41550-020-01268-x
  * Sanders et al. (1988) Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ, 325, 74. doi:10.1086/165983
  * Santini et al. (2017) Santini, P., Fontana, A., Castellano, M., et al. 2017, ApJ, 847, 76. doi:10.3847/1538-4357/aa8874
  * Saracco et al. (2020) Saracco, P., Marchesini, D., La Barbera, F., et al. 2020, ApJ, 905, 40. doi:10.3847/1538-4357/abc7c4
  * Schreiber et al. (2017) Schreiber, C., Pannella, M., Leiton, R., et al. 2017, A&A, 599, A134. doi:10.1051/0004-6361/201629155
  * Schreiber et al. (2018) Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30. doi:10.1051/0004-6361/201731506
  * Schreiber et al. (2018b) Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018, A&A, 618, A85. doi:10.1051/0004-6361/201833070
  * Simpson et al. (2017) Simpson, J. M., Smail, I., Wang, W.-H., et al. 2017, ApJ, 844, L10. doi:10.3847/2041-8213/aa7cf2
  * Steinhardt et al. (2016) Steinhardt, C. L., Capak, P., Masters, D., et al. 2016, ApJ, 824, 21. doi:10.3847/0004-637X/824/1/21
  * Straatman et al. (2014) Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2014, ApJ, 783, L14. doi:10.1088/2041-8205/783/1/L14
  * Tacchella et al. (2021) Tacchella, S., Conroy, C., Faber, S., et al. 2021, Extragalactic Spectroscopic Surveys: Past, 59. doi:10.5281/zenodo.4721715
  * Valentino et al. (2015) Valentino, F., Daddi, E., Strazzullo, V., et al. 2015, ApJ, 801, 132. doi:10.1088/0004-637X/801/2/132
  * Valentino et al. (2020) Valentino, F., Tanaka, M., Davidzon, I., et al. 2020, ApJ, 889, 93. doi:10.3847/1538-4357/ab64dc


# CP2 Skyrmions and Skyrmion Crystals in Realistic Quantum Magnets

Hao Zhang  Department of Physics and Astronomy, The University of Tennessee,
Knoxville, Tennessee 37996, USA  Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  Zhentao Wang
Department of Physics and Astronomy, The University of Tennessee, Knoxville,
Tennessee 37996, USA  School of Physics and Astronomy, University of
Minnesota, Minneapolis, Minnesota 55455, USA  David Dahlbom  Department of
Physics and Astronomy, The University of Tennessee, Knoxville, Tennessee
37996, USA  Kipton Barros  Theoretical Division and CNLS, Los Alamos National
Laboratory, Los Alamos, New Mexico 87545, USA  Cristian D. Batista  Department
of Physics and Astronomy, The University of Tennessee, Knoxville, Tennessee
37996, USA  Quantum Condensed Matter Division and Shull-Wollan Center, Oak
Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA

(August 27, 2024)

###### Abstract

Magnetic skyrmions are nanoscale topological textures that have been recently
observed in different families of quantum magnets. These textures are known as
CP1 skyrmions because the target manifold of the magnetization field is the 2D
sphere $S^{2}\cong$ CP1. Here we report the emergence of magnetic CP2
skyrmions in realistic spin-$1$ models. Unlike CP1 skyrmions, CP2 skyrmions
can also arise as metastable textures of quantum paramagnets, opening a new
road to discover emergent topological solitons in non-magnetic materials. The
quantum phase diagram of the spin-$1$ models also includes magnetic field-
induced CP2 skyrmion crystals that can be detected with regular momentum-
(diffraction) and real-space (Lorentz transmission electron microscopy)
experimental techniques.

##  I Introduction

Lord Kelvin’s vision of the atom as a vortex in ether [1] inspired Skyrme [2,
3] to explain the origin of nucleons as emergent topologically non-trivial
configurations of a pion field described by a 3+1 dimensional O(4) non-linear
$\sigma$-model. In the modern language, these “skyrmions” are examples of
topological solitons and Skyrme’s model has become the prototype of a
classical theory that supports these solutions. Besides its important role in
high-energy physics and cosmology, Skyrme’s model also led to important
developments in other areas of physics. For instance, the baby Skyrme model
[4, 5, 6] (planar reduction of the non-linear $\sigma$-model), which is as an
extension of the Heisenberg model [4, 5, 7], has baby skyrmion solutions in
the presence of a chiral symmetry breaking Dzyaloshinskii-Moriya (DM)
interaction [8, 9, 10, 11].

Periodic arrays of magnetic skyrmions and single skyrmion metastable states
were originally observed in chiral magnets, such as MnSi, Fe1-xCoxSi, FeGe and
Cu2OSeO3 [12, 13, 14, 15, 16]. This discovery sparked the interest of the
community at large and spawned efforts in multiple directions. Identifying
realistic conditions for the emergence of novel magnetic skyrmions is one of
the main goals of modern condensed matter physics. Novel mechanisms are
usually accompanied by new properties. For instance, while skyrmions of chiral
magnets have a fixed vector chirality, this is still a degree of freedom in
centrosymmetric materials, such as BaFe1-x-0.05ScxMg0.05O19,
La2-2xSr1+2xMn2O7, Gd2PdSi3 and Gd3Ru4Al12 [17, 18, 19, 20, 21, 22, 23], where
skyrmions arise from frustration, i.e., from competing exchange or dipolar
interactions [24, 25, 26, 27, 28, 29, 30].

The target manifold of the above-mentioned planar baby skyrmions is
$S^{2}\cong$ CP1, i.e., the usual 2D sphere, corresponding to normalized
dipoles. More generally, one may consider the complex projective space CPN-1
that represents the normalized $N$-component complex vectors, up to an
irrelevant complex phase. The topologically distinct, smooth mappings from the
base manifold $S^{2}$ (2D sphere $\cong$ compactified plane) to the target
manifold CPN-1 can be labeled by the integers:
$\Pi_{2}(\mathrm{CP}^{N-1})=\mathbb{Z}$. This homotopy group suggests
generalizations of the planar Skyrme’s model to $N>2$, such as the CP2 non-
linear $\sigma$-model [31, 32, 33] and in the Faddeev-Skyrme type model [34,
35]. In recent work, Akagi et al. considered the SU(3) version of the
Heisenberg model with a DM interaction, whose continuum limit becomes a gauged
CP2 nonlinear $\sigma$-model with a background uniform gauge field [36]. An
attractive aspect of this model is that it admits analytical solution by the
application of techniques developed for the gauged non-linear $\sigma$-model.
However, it may be challenging to find materials described by this model
because SU(3) can only be an accidental symmetry of spin-spin interactions of
real quantum magnets.

The main purpose of this work is to demonstrate that exotic CP2 skyrmions
readily emerge in simple and realistic spin-$1$ ($N=3$) models and their
higher-spin generalizations. Remarkably, isolated CP2 skyrmions can either be
metastable states of a quantum paramagnet (QPM) or of a fully polarized (FP)
ferromagnet. Unlike the “usual” CP1 magnetic skyrmions, the dipolar field of
the metastable CP2 skyrmions of quantum paramagnets vanishes away from the
skyrmion core. Moreover, the application of an external magnetic field to the
QPM induces stable triangular crystals of CP2 skyrmions in the field interval
that separates the QPM from the FP state.

##  II Model

To illustrate the basic ideas we consider a minimal spin-$1$ model defined on
the triangular lattice (TL):

$\\!\\!\\!{\hat{\mathcal{H}}}=\sum_{\langle
i,j\rangle}J_{ij}\left({\hat{S}}^{x}_{i}{\hat{S}}^{x}_{j}+{\hat{S}}^{y}_{i}{\hat{S}}^{y}_{j}+\Delta{\hat{S}}^{z}_{i}{\hat{S}}^{z}_{j}\right)-h\sum_{i}{\hat{S}}^{z}_{i}+D\sum_{i}\left({\hat{S}}^{z}_{i}\right)^{2}.$
(1)

The first term includes an easy-axis ferromagnetic (FM) nearest-neighbor
exchange interaction $J_{1}<0$ and a second-nearest-neighbor antiferromagnetic
(AFM) exchange $J_{2}>0$. For simplicity, we assume that the exchange
anisotropy, defined by the dimensionless parameter $\Delta>1$, is the same for
both interactions. The second and third terms represent the Zeeman coupling to
an external field and an easy-plane single-ion anisotropy ($D>0$),
respectively. $\hat{\mathcal{H}}$ is invariant under the space-group of the TL
and the U(1) group of global spin rotations along the field-axis. We will
adopt $|J_{1}|$ as the unit of energy (i.e. $J_{1}=-1$).

The first step is to take the classical limit of $\hat{\mathcal{H}}$, where
multiple approaches are possible [37, 38]. The traditional classical limit is
based on SU(2) coherent states, which retains only the spin dipole expectation
value, and produces the Landau-Lifshitz spin dynamics. This approach is
adequate for modeling systems with weak single-ion anisotropy $D\ll|J_{1}|$.
To model systems in the regime $D\gtrsim|J_{1}|$, and correctly capture low-
energy excitations, it is important to retain more information about spin
fluctuations in the $(N=3)$-dimensional local Hilbert space. Specifically, our
classical limit will assume that the many-body quantum state is a direct
product of SU(3) coherent states [39, 37, 40, 41, 42, 43, 38]:

$|\bm{Z}\rangle=\otimes_{j}|\bm{Z}_{j}\rangle\;\;{\rm
with}\;\;|\bm{Z}_{j}\rangle=\sum_{a}Z^{a}_{j}\left|x^{a}\right\rangle_{j},$
(2)

where $\bm{Z}_{j}=\left(Z^{1}_{j},Z^{2}_{j},Z^{3}_{j}\right)^{\mathrm{T}}$ is
a complex vector of unit length and
$\\{\left|x^{1}\right\rangle_{j},\left|x^{2}\right\rangle_{j},\left|x^{3}\right\rangle_{j}\\}$
is an orthonormal basis for the local Hilbert state on site $j$.

Local physical operators are represented by Hermitian matrices that act on
SU(3) coherent states. The space of $3\times 3$ traceless, Hermitian matrices
comprises the fundamental representation of the $\mathfrak{su}(3)$ Lie
algebra. A basis $\hat{T}^{\mu}$ ($\mu=1,\dots,8$) for this space is
characterized by the commutation relations,

$\left[{\hat{T}}_{j}^{\eta},{\hat{T}}^{\mu}_{j}\right]=if_{\eta\mu\nu}{\hat{T}}_{j}^{\nu},$
(3)

where we are using the convention of summation over repeated Greek indices. We
may additionally impose an orthonormality condition

$\operatorname{Tr}\left({\hat{T}}_{j}^{\alpha}{\hat{T}}_{j}^{\beta}\right)=2\delta_{\alpha\beta}.$
(4)

It is conventional to define the structure constants as
$f_{\eta\mu\nu}=-\frac{i}{2}\operatorname{Tr}\left({\lambda}_{\eta}\left[{\lambda}_{\mu},{\lambda}_{\nu}\right]\right),$
where ${\lambda}_{\mu}$ are the Gell-Mann matrices.

The spin dipole operators
$\hat{\bm{S}}_{j}=(\hat{S}^{x}_{j},\hat{S}^{y}_{j},\hat{S}^{z}_{j})^{T}$
acting on site $j$ are generators for a spin-$1$ representation of SU(2). It
is possible to formulate generators of SU(3) as polynomials of these spin
operators,

$\displaystyle\left(\begin{array}[]{l}{\hat{T}}_{j}^{7}\\\
{\hat{T}}_{j}^{5}\\\
{\hat{T}}_{j}^{2}\end{array}\right)=-\left(\begin{array}[]{c}{\hat{S}}_{j}^{x}\\\
{\hat{S}}_{j}^{y}\\\
{\hat{S}}_{j}^{z}\end{array}\right),\left(\begin{array}[]{l}{\hat{T}}_{j}^{3}\\\
{\hat{T}}_{j}^{8}\\\ {\hat{T}}_{j}^{1}\\\ {\hat{T}}_{j}^{4}\\\
{\hat{T}}_{j}^{6}\end{array}\right)=\left(\begin{array}[]{c}-\left({\hat{S}}_{j}^{x}\right)^{2}+\left({\hat{S}}_{j}^{y}\right)^{2}\\\
\frac{1}{\sqrt{3}}\left[3\left({\hat{S}}_{j}^{z}\right)^{2}-{\hat{\bm{S}}}_{j}^{2}\right]\\\
{\hat{S}}_{j}^{x}{\hat{S}}_{j}^{y}+{\hat{S}}_{j}^{y}{\hat{S}}_{j}^{x}\\\
-{\hat{S}}_{j}^{z}{\hat{S}}_{j}^{x}-{\hat{S}}_{j}^{x}{\hat{S}}_{j}^{z}\\\
{\hat{S}}_{j}^{y}{\hat{S}}_{j}^{z}+{\hat{S}}_{j}^{z}{\hat{S}}_{j}^{y},\end{array}\right),$
(21)

where $T_{j}^{7,5,2}$ are the dipolar components of the spin-1 degree of
freedom, while the other five generators are the quadrupolar components. Here
we have adopted the notation and conventions of Ref. [36] to make closer
contact with the literature on high-energy physics 111Our definitions for
$\hat{S}^{x}$ and $\hat{S}^{z}$ differ from these two in Ref [36] by a minus
sign..

Let $|1\rangle_{j}$, $|0\rangle_{j}$, and $|\bar{1}\rangle_{j}$ denote the
normalized eigenstates of ${\hat{S}}^{z}_{j}$, with eigenvalues, 1, 0 and -1,
respectively. In the Cartesian basis,

$\left|x^{1}\right\rangle_{j}=\frac{i\left[|1\rangle_{j}-|\bar{1}\rangle_{j}\right]}{\sqrt{2}},\left|x^{2}\right\rangle_{j}=\frac{\left[|1\rangle_{j}+|\bar{1}\rangle_{j}\right]}{\sqrt{2}},\left|x^{3}\right\rangle_{j}=-i|0\rangle_{j},$
(22)

the SU(3) generators given in Eq. (21) are the Gell-Mann matrices:

$\langle
x^{a}_{j}|{\hat{T}}_{j}^{\mu}|x^{b}_{j}\rangle=\left({\lambda}_{\mu}\right)_{ab}\quad\mu=1,2,\cdots,8.$
(23)

The orbit of coherent states $|\bm{Z}_{j}\rangle$ is obtained by applying
SU(3) transformations to the highest weight state $|1\rangle_{j}$ [37]:
$|\bm{Z}_{j}\rangle={\hat{U}}_{j}|1\rangle_{j}$. Since the global phase is a
gauge degree of freedom, the orbit of physical SU(3) coherent states is
$S^{5}/S^{1}\cong\mathrm{CP}^{2}$. The “SU(3) classical limit” of the spin
Hamiltonian (1) is obtained by replacing the Hamiltonian operator
${\hat{\mathcal{H}}}$ with its expectation value

${\cal H}\equiv\langle\bm{Z}|{\hat{\mathcal{H}}}|\bm{Z}\rangle,$ (24)

after rewriting ${\hat{\mathcal{H}}}$ in terms of SU(3) spin components,

${\hat{\mathcal{H}}}=\sum_{\langle
i,j\rangle}{I}^{\mu}_{ij}{\hat{T}}^{\mu}_{i}{\hat{T}}^{\mu}_{j}-\sum_{i}B^{\mu}{\hat{T}}^{\mu}_{i},$
(25)

where
${I}^{\mu}_{ij}=J_{ij}(\Delta\delta_{\mu,2}+\delta_{\mu,5}+\delta_{\mu,7})$
and $B^{\mu}=(-h\delta_{\mu,2}-D\delta_{\mu,8}/\sqrt{3})$. Because of the
direct product form of Eq. (2), ${\cal H}$ can be expressed as a function of
the “color field”

$n^{\mu}_{j}\equiv\langle\bm{Z}_{j}|{\hat{T}}_{j}^{\mu}|\bm{Z}_{j}\rangle=\left({\lambda}_{\mu}\right)_{ab}\bar{Z}^{a}_{j}Z^{b}_{j},$
(26)

which satisfies the constraints

$n^{\mu}n^{\mu}=\frac{4}{3},\quad
n^{\mu}=\frac{3}{2}d_{\mu\nu\eta}n^{\nu}n^{\eta},$ (27)

where
$d_{\mu\nu\eta}=\frac{1}{4}\operatorname{Tr}\left(\lambda_{\mu}\left\\{\lambda_{\nu},\lambda_{\eta}\right\\}\right)$.
This in turn leads to the Casimir identity:
$d_{mpq}n^{m}n^{p}n^{q}=\frac{8}{9}.$ In terms of this color field, we can
express

${\mathcal{H}}=\sum_{\langle
i,j\rangle}{I}^{\mu}_{ij}{n}^{\mu}_{i}{n}^{\mu}_{j}-\sum_{i}B^{\mu}{n}^{\mu}_{i}.$
(28)

To avoid an explicit use of the structure constants ($f_{\eta\mu\nu}$), we
introduce the operator field
$\bm{\mathfrak{n}}_{j}=n^{\mu}_{j}{\lambda}_{\mu}$. Topological soliton
solutions of the color field become well-defined in the continuum (long
wavelength) limit, where the Hamiltonian can be approximated by

$\mathcal{H}\simeq\int\mathrm{d}r^{2}\left[-\frac{{\cal
I}^{\mu}_{1}}{2}\partial_{a}{n}^{\mu}\partial_{a}{n}^{\mu}+\frac{{\cal
I}^{\mu}_{2}}{2}\left(\nabla^{2}{n}^{\mu}\right)^{2}-\mathcal{B}^{\mu}n^{\mu}\right].$
(29)

Here $\partial_{a}$ denotes the spatial derivatives and there is an implicit
summation over the repeated $a$ index. The coupling constants can be expressed
in terms of the parameters of the lattice model (25):

$\displaystyle{\cal I}^{\mu}_{1}$ $\displaystyle=$
$\displaystyle\frac{3}{2}({I}^{\mu}_{1}+3{I}^{\mu}_{2}),\quad{\cal
I}^{\mu}_{2}=\frac{3}{32}(I^{\mu}_{1}+9I^{\mu}_{2}),$ $\displaystyle{\cal
B}^{\mu}$ $\displaystyle=$ $\displaystyle
B^{\mu}-3(\Delta-1)(J_{1}+J_{2})\delta_{\mu,8}.$ (30)

Eq. (29) corresponds to an anisotropic CP2 model. For skyrmion solutions the
base plane $\mathbb{R}^{2}$ can be compactified to $S^{2}$ because the color
field takes a constant value $n_{\infty}$ at spatial infinity. These spin
textures can then be characterized by the topological charge of the mapping
$\bm{\mathfrak{n}}:\mathbb{R}^{2}\sim S^{2}\mapsto CP^{2}$:

$C=-\frac{i}{32\pi}\int\mathrm{d}x\mathrm{d}y\varepsilon_{jk}\operatorname{Tr}\left(\bm{\mathfrak{n}}\left[\partial_{j}\bm{\mathfrak{n}},\partial_{k}\bm{\mathfrak{n}}\right]\right).$
(31)

For the lattice systems of interest, the CP2 skyrmion charge can be computed
after interpolating the color fields on nearest-neighbor sites
$\bm{\mathfrak{n}}_{j}$ and $\bm{\mathfrak{n}}_{k}$ along the CP2 geodesic:

$C=\sum_{\triangle_{jkl}}\rho_{jkl}=\frac{1}{2\pi}\sum_{\triangle_{jkl}}\left(\gamma_{jl}+\gamma_{lk}+\gamma_{kj}\right),$
(32)

where $\triangle_{jkl}$ denotes each oriented triangular plaquette of nearest-
neighbor sites $j\to k\to l$ and
$\gamma_{kj}=\arg\left[\left\langle\bm{Z}_{k}\mid\bm{Z}_{j}\right\rangle\right]$
is the Berry connection on the bond $j\to k$ [45].

Figure 1: $T=0$ phase diagram of the classical Hamiltonian ${\mathcal{H}}$ as
a function of the single-ion anisotropy $D$ and the external field $h$, for
$J_{2}/|J_{1}|=2/(1+\sqrt{5})$ and $\Delta=2.6$. The two insets show the
phases for small-$D$ and large-$D$, respectively. The solid (dashed) lines
indicate 1st- (2nd-) order phase transitions.

Figure 2: a, b Real space distribution of the dipolar sector of the CP2
skyrmion crystals SkX-I and SkX-II. The length of the arrow represents the
magnitude of the dipole moment of the color field
$|\langle\hat{\bm{S}}_{j}\rangle|=\sqrt{(n_{j}^{7})^{2}+(n_{j}^{5})^{2}+(n_{j}^{2})^{2}}$.
The color scale of the arrows indicates
$\langle\hat{S}^{z}_{j}\rangle=-n_{j}^{2}$. The insets display the static spin
structure factors
$\mathcal{S}^{\perp}(\bm{q})=n^{7}_{\bm{q}}n^{7}_{\bar{\bm{q}}}+n^{5}_{\bm{q}}n^{5}_{\bar{\bm{q}}}$
and $\mathcal{S}^{zz}(\bm{q})=n^{2}_{\bm{q}}n^{2}_{\bar{\bm{q}}}$, with
$\bm{n}_{\bm{q}}=\sum_{j}e^{{i\mkern
1.0mu}\bm{q}\cdot\bm{r_{j}}}\bm{n}_{j}/L$. The CP2 skyrmion density
distribution $\rho_{jkl}$ [see Eq. (32)] is indicated by the color of the
triangular plaquettes

.

##  III  Phase diagram

The $T=0$ phase diagram is obtained by numerically minimizing the classical
spin Hamiltonian $\mathcal{H}$ (28) in the $4L^{2}$-dimensional phase space of
a magnetic cell of $L\times L$ spins (see Methods). The shape and size of this
unit cell is dictated by the symmetry-related magnetic ordering wave vectors
$\bm{Q}_{\nu}$ ($\nu=1,2,3$) [see Figs. 2a and b], which are determined by
minimizing the exchange interaction in momentum space:
$J({\bm{q}})=\sum_{jl}J_{jl}e^{{i\mkern
1.0mu}\bm{q}\cdot(\bm{r}_{j}-\bm{r}_{l})}.$ The ratio between both exchange
interactions, $J_{2}/|J_{1}|=2/(1+\sqrt{5})$, is tuned to fix the magnitude of
the ordering wave vectors, $|\bm{Q}_{\nu}|=|\bm{b}_{1}|/5$ [27], corresponding
to a magnetic unit cell of linear size $L=5$. A we will see later, the
relevant qualitative aspects of the phase diagram do not depend on the
particular choice of the model. The three ordering wave vectors, which are
related by the $C_{6}$ symmetry of the TL, are parallel to the $\Gamma$-Mν
directions (denoted in Fig. 2).

The resulting phase diagram shown in Fig. 1 includes multiple magnetically
ordered phases between the FP phase and the QPM phase, where every spin is in
the $|0\rangle$ state. For $D\gg|J_{1}|$, these phases include two field-
induced CP2 skyrmion crystals (SkX-I and SkX-II), separated by two modulated
vertical spiral phases (MVS-I and MVS-II), where the polarization plane of the
spiral is parallel to the $c$-axis and the magnitude of the dipole moment is
continuously suppressed as the moment rotates from $\hat{\bm{z}}$ to
$-\hat{\bm{z}}$ directions. The spiral phases have the same symmetry and are
separated by a first order metamagnetic transition. As shown in Fig. 2a, the
CP2 skyrmions of the SkX-I crystal have dipole moments that evolve
continuously into the purely nematic state ($|0\rangle$) as they move away
from the core. Conversely, Fig. 2b shows that the spins in the SkX-II phase
have a strong quadrupolar character (the small dipolar moment is completely
suppressed in the large $D/|J_{1}|$ limit) at the skyrmion core, and evolve
continuously into the magnetic state $|1\rangle$ as they move away from the
core. The CP2 skyrmion density distribution $\rho_{jkl}$ is also indicated
with colored triangular plaquettes in Fig. 2a, b for SkX-I and SkX-II,
respectively. As shown in the inset of Fig. 1, the phase SkX-II extends down
to $D/|J_{1}|\simeq 5$, while the phase SkX-I disappears near $D/|J_{1}|\simeq
8$.

New competing orderings appear in intermediate $D/|J_{1}|$ region. In
particular, a significant fraction of the phase diagram is occupied by the so-
called canted spiral (CS) phase,

$|\bm{Z}_{j}\rangle=\cos\theta|0\rangle_{j}+e^{i\bm{Q}\cdot\bm{r}_{j}}\sin\theta\cos\phi|1\rangle_{j}+e^{-i\bm{Q}\cdot\bm{r}_{j}}\sin\theta\sin\phi|\bar{1}\rangle_{j},$
(33)

where $\theta$ and $\phi$ are variational parameters, and $\bm{Q}$ can take
any values among {$\bm{Q}_{1}$, $\bm{Q}_{2}$, $\bm{Q}_{3}$}. Upon increasing
$D$, the magnitude of the dipole moment of each spin,
$|\langle\hat{\bm{S}}_{j}\rangle|$, is continuously suppressed to zero at the
boundary,

$D_{c}=h\sqrt{1-\frac{4J^{2}(\bm{Q})}{h^{2}+4J^{2}(\bm{Q})}}-2J(\bm{Q})\left(1-\frac{2J(\bm{Q})}{\sqrt{h^{2}+4J^{2}(\bm{Q})}}\right),$
(34)

that signals the second order transition into the QPM phase. As shown in Fig.
1, several competing phases appear above the CS phase upon increasing $h$.
These phases include three triple-$\bm{Q}$ spiral orderings [3$\bm{Q}$ I-III]
with dominant weight in one of three $\bm{Q}$ transverse components and a
staggered distribution of the CP2 skyrmion density $\rho_{jkl}$ [see Eq. (32)]
and three different modulated double-$\bm{Q}$ orderings (MDQ I - III) and two
triple-$\bm{Q}$ orderings. All of these phases are described in detail in the
supplementary information. In the rest of the paper we will focus on the SkX
phases and the single-skyrmion metastable solutions that emerge in their
proximity.

Figure 3: Spectrum of the single-ion model
$\hat{\mathcal{H}}_{\text{SI}}=D(\hat{S}^{z})^{2}-h\hat{S}^{z}$. The shaded
region denotes the effective regime with a quasi-degenerate doublet:
$\\{|0\rangle,|1\rangle\\}$

.

##  IV Large-$D$ limit

The origin of the CP2 skyrmion crystals can be understood by analyzing the
small $|J_{ij}|/D$ regime, where $\hat{\mathcal{H}}$ can be reduced via first
order degenerate perturbation theory in $J_{ij}/D$ to an effective pseudo-
spin-$1/2$ low-energy Hamiltonian,

$\hat{\mathcal{H}}_{\text{eff}}=\sum_{\langle
i,j\rangle}{\tilde{J}}_{ij}(\hat{s}^{x}_{i}\hat{s}^{x}_{j}+\hat{s}^{y}_{i}\hat{s}^{y}_{j}+\tilde{\Delta}\hat{s}^{z}_{i}\hat{s}^{z}_{j})-{\tilde{h}}\sum_{i}\hat{s}^{z}_{i}.$
(35)

The pseudo-spin-$1/2$ operators are the projection of the original spin
operators into the low-energy subspace ${\cal S}_{0}$ generated by the quasi-
degenerate doublet $\\{|0\rangle_{j},|1\rangle_{j}\\}$ (see Fig. 3):

$\hat{s}^{z}_{j}={\cal P}_{0}\hat{S}^{z}_{j}{\cal
P}_{0}-\frac{1}{2},\quad\hat{s}^{\pm}_{j}=\frac{{\cal
P}_{0}\hat{S}^{\pm}_{j}{\cal P}_{0}}{\sqrt{2}},$ (36)

where ${\cal P}_{0}$ is the projection operator of the low-energy subspace.
Importantly, the first state of the doublet has a net quadrupolar moment but
no net dipole moment, $\langle 0|\hat{\bm{S}}_{j}|0\rangle_{j}=0$, while the
second state maximizes the dipole moment along the $\hat{\bm{z}}$-direction
$\langle 1|\hat{\bm{S}}_{j}|1\rangle_{j}=\hat{\bm{z}}$. This means that three
pseudo-spin operators generate an SU(2) subgroup of SU(3) different from the
SU(2) subgroup of spin rotations.

$\hat{\mathcal{H}}_{\text{eff}}$ represents an effective triangular easy-axis
XXZ model with effective exchange, anisotropy and field parameters
${\tilde{J}}_{ij}=2J_{ij}$, $\tilde{\Delta}=\frac{\Delta}{2}$ and
${\tilde{h}}=h-D-3\Delta(J_{1}+J_{2})$, respectively. This model is known to
exhibit a field-induced CP1 SkX phase [25, 27] that survives in the long-
wavelength limit [26]:

$\mathcal{H}_{\text{eff}}\simeq\int\mathrm{d}r^{2}\left[-\frac{{\cal
J}^{\eta}_{1}}{2}\partial_{a}\tilde{n}^{\eta}\cdot\partial_{a}\tilde{n}^{\eta}+\frac{{\cal
J}^{\eta}_{2}}{2}\left(\nabla^{2}\tilde{n}^{\eta}\right)^{2}-\tilde{\cal
B}\tilde{n}_{z}+\tilde{\cal D}\tilde{n}_{z}^{2}\right],$ (37)

where $\eta=x,y,z$ denotes the three components of the unit vector field
$\tilde{\bm{n}}$ ($|\tilde{\bm{n}}|=1$), and

$\displaystyle{\cal J}^{\eta}_{1}$ $\displaystyle=$
$\displaystyle\frac{3s^{2}}{2}(\tilde{J}_{1}+3\tilde{J}_{2})[1+(\tilde{\Delta}-1)\delta_{\eta
z}],$ $\displaystyle{\cal J}^{\eta}_{2}$ $\displaystyle=$
$\displaystyle\frac{3s^{2}}{32}(\tilde{J}_{1}+9\tilde{J}_{2})[1+(\tilde{\Delta}-1)\delta_{\eta
z}]$ $\displaystyle\tilde{\cal B}$ $\displaystyle=$ $\displaystyle
s{\tilde{h}},\quad\tilde{\cal
D}=3s^{2}(\tilde{\Delta}-1)(\tilde{J}_{1}+\tilde{J}_{2}),$ (38)

where $s=1/2$. Although the target manifold of this theory is CP1 (orbit of
SU(2) coherent states that belong ${\cal S}_{0}$), we must keep in mind that
$\hat{\mathcal{H}}_{\text{eff}}$ describes the large $D/|J_{1}|$ limit where
the CP2 skyrmions of the original spin-$1$ model become asymptotically close
to CP1 _pseudo-spin_ skyrmions. In other words, the SkXs include a finite
$|\bar{1}\rangle$ component for finite $D/|J_{1}|$, which increases upon
decreasing $D$. This component, which only appears in the low-energy model
when second order corrections in $J_{ij}$ are included, is responsible for the
metamagnetic transition between the MVS-I and MVS-II phases (the transition
disappears in the $D\to\infty$ limit).

Figure 4: Dipolar sector of CP2 skyrmions. The color scale indicates the value
of $n_{j}^{2}$ ($\langle\hat{S}_{j}^{z}\rangle$). a, b Skyrmion excitation on
top of a QPM background. c, d Skyrmion excitation on top of a fully polarized
background. $J_{2}/|J_{1}|=2/(1+\sqrt{5})$ and $\Delta=2.6$ in (a), (c), and
(d). $J_{2}/|J_{1}|=2/(3+\sqrt{5})$ and $\Delta=2.2$ in b. In these panels, a.
$D=17.1|J_{1}|$, $H=13|J_{1}|$. b. $D=18.3|J_{1}|$, $H=14|J_{1}|$. c.
$D=7|J_{1}|$, $H=5|J_{1}|$. d. $D=4|J_{1}|$, $H=2|J_{1}|$.

Since $\hat{\mathcal{H}}_{\text{eff}}(h)$ and
$\hat{\mathcal{H}}_{\text{eff}}(-h)$ are related by a pseudo-time-reversal
(PTR) transformation ($\hat{s}_{j}\to-\hat{s}_{j}$ on the lattice and
$\tilde{n}\to-\tilde{n}$ in the continuum) their corresponding ground states
are related by the same transformation. In particular, the ground state
$(\tilde{\bm{n}}=\hat{\bm{z}})$ that is obtained above the saturation field
($\tilde{\cal B}>\tilde{\cal B}_{\rm sat}$) corresponds to the FP state
($\langle\hat{\bm{S}}_{j}\rangle=\hat{\bm{z}}$) in the original spin-$1$
variables, while the ground state $(\tilde{\bm{n}}=-\hat{\bm{z}})$ below the
negative saturation field ($\tilde{\cal B}<-\tilde{\cal B}_{\rm sat}$)
corresponds to the QPM phase ($|\bm{Z}_{j}\rangle=|0\rangle_{j}$).
Correspondingly, the SkX induced by a positive $h$ has pseudo-spins polarized
along the quadrupolar direction ($|0\rangle$) near the core of the skyrmions
and parallel to the dipolar one ($|1\rangle$) at the midpoints between two
cores. This explains the origin of the SkX-II crystals depicted in Fig. 2 b.
The negative ${\cal B}$ counterpart of this phase, which is obtained by
applying the PTR transformation, corresponds to the SkX-I crystal shown in
Fig. 2a. In this case the skyrmion cores are magnetic, while the midpoints are
practically quadrupolar (they become purely quadrupolar in the large
$D/|J_{1}|$ limit). This simple reasoning explains the origin of the novel SkX
phases included in the $T=0$ phase diagram of $\mathcal{H}$ shown in Fig. 1.
The intermediate phase between the SkX-I and SkX-II ground state of
${\mathcal{H}}$ induced by positive and negative values of $h$ is a
single-$\bm{Q}$ spiral with a polarization plane parallel to the $c$-axis
known as vertical spiral (VS). This explains the origin of the MVS-I and MVS-
II phases in between the two SkX phases (the first order transition between
both phases disappears in the large-$D$ limit [25]).

Single-skyrmion solutions. Besides the SkX phases shown in Fig. 4, the
effective field theory (37) is known to support metastable CP1 single-skyrmion
solutions beyond the saturation fields
$|\tilde{\mathcal{B}}|>\tilde{\mathcal{B}}_{\rm sat}$. The pseudo-spin
variable is anti-parallel to the external field at the core and it gradually
rotates towards the direction parallel to the field upon moving away from the
center. Interestingly, this pseudo-spin texture translates into metastable
single-skyrmion solutions of the QPM phase that have a magnetic core and a
nematic periphery, as it is shown in Figs. 4a and b for different sets of
Hamiltonian parameters. The CP2 skyrmions are metastable solutions in the QPM
phase for $D\gtrsim 14$, implying that these exotic magnetic-nematic textures
should emerge in real magnets under quite general conditions.

Similarly, the metastable pseudo-spin single-skyrmion solutions of the FP
phase ($\tilde{\mathcal{B}}>\tilde{\mathcal{B}}_{\rm sat}$) lead to a spin
texture with a nematic (non-magnetic) core and a magnetic (FP) periphery, like
the one shown in Fig. 4c. Interestingly, this exotic CP2 skyrmion solution
remains metastable down to $D\simeq 4|J_{1}|$ and it coexists with regular
(CP1) metastable skyrmion solutions, like the one shown Fig. 4d, that emerge
below $D\simeq 4.25|J_{1}|$.

##  V  Discussion

In summary, we have demonstrated that CP2 skyrmion textures emerge in
realistic models of hexagonal magnets out of the combination of geometric
frustration with competing exchange and single-ion anisotropies. It is
important to note that the skyrmion crystals and metastable solutions reported
in this work survive in the long wavelength limit [26], implying that the
above described CP2 skyrmion phases should also exist in extensions of the
model to honeycomb and Kagome lattice geometries. The generic spin-$1$ model
considered in this work describes a series of triangular antiferromagnets in
the form of ABX3, BX2, and ABO2 [46, 47, 48], where A is a an alkali metal, B
is a transition metal, and X is a halogen atom. Several Ni-based compounds,
including NiGa2S4 [49], Ba3NiSb2O9 [50], Na2BaNi(PO4)2 [51], are also found to
be realizations of spin-$1$ models on TLs. Some of these compounds have been
already identified as candidates to host crystals of CP1 magnetic skyrmions
that are stabilized by the combination of frustration and exchange anisotropy
[52]. Others, such as FeI2 [53, 54], are described by the the Hamiltonian (1)
but the sign of the single-ion and exchange anisotropies is opposite to the
case of interest in this work. Related compounds, such as CsFeCl3, are known
to be quantum paramagnets described by the Hamiltonian (1) [55] with a
dominant easy-plane single-ion anisotropy $D/J_{1}\simeq 10$. An alternative
route to find realizations of our spin-$1$ Hamiltonian is to consider
hexagonal materials comprising $4f$ magnetic ions with a singlet single-ion
ground state and an excited Ising-like doublet (the effective easy-plane
single-ion anisotropy $D$ is equal to the singlet-doublet gap). Ultracold
atoms are also powerful platforms to realize spin-$1$ models with tunable
single-ion anisotropy [56].

Our results demonstrate that magnetic field induced CP2 skyrmion crystals
should emerge in the presence of a dominant single-ion easy-plane anisotropy
$D$ that is strong enough to stabilize a QPM at $T=0$. In terms of SU(3)
spins, the single-ion anisotropy acts as an external SU(3) field that couples
linearly to one component (${\hat{T}}_{j}^{8}$) of the quadrupolar moment.
Correspondingly, the QPM can be regarded as a uniform quadrupolar state
induced by a strong enough ${\hat{T}}_{j}^{8}$ component of the SU(3) field.
The field-induced quantum phase transition between this uniform quadrupolar
state and the skyrmion crystals is presaged by the emergence of metastable CP2
single-skyrmion solutions consisting of a magnetic skyrmion core that decays
continuously into a quadrupolar periphery. These novel skyrmions can be
induced by increasing the strength of the magnetic field in the neighborhood
of a given magnetic ion of a frustrated quantum paramagnet with competing
exchange and single-ion anisotropies.

##  VI Acknowledgments

We acknowledge useful discussions with Xiaojian Bai, Antia Botana, Ying Wai
Li, Shizeng Lin, Cole Miles, Martin Mourigal, Sakib Matin, Matthew Wilson, and
Shang-Shun Zhang. D. D., K. B. and C.D.B. acknowledge support from U.S.
Department of Energy, Office of Science, Office of Basic Energy Sciences,
under Award No. DE-SC0022311. The work by H.Z. was supported by the Graduate
Advancement, Training and Education (GATE) fellowship. Z.W. was supported by
the U.S. Department of Energy through the University of Minnesota Center for
Quantum Materials, under Award No. DE-SC-0016371.

##  VII Methods

The numerical minimization for the phase diagram Fig. 1 is done in a cell of
$10\times 10$ spins containing four magnetic unit cells ($L=5$). Two crucial
steps are useful to improve the efficiency of the local gradient-based
minimization algorithms [57]. In the first step, we set multiple random
initial conditions $|\bm{Z}\rangle$ ($\sim$ 50 for our case), where
$|\bm{Z}_{j}\rangle$ on _every_ site $j$ is uniformly sampled on the
CP${}^{2}\simeq S^{5}/S^{1}$ manifold. After running the minimization
algorithm, we keep the solution with the lowest energy for a given set of
Hamiltonian parameters. In the next step, half of the initial conditions are
randomly generated, while the other half correspond to the lowest-energy
solutions obtained in the first step within a predefined neighborhood of the
Hamiltonian parameters. This procedure is iterated until the phase diagram
converges.

## References

  * [1] Thomson, W.  Ii. on vortex atoms.  _London, Edinburgh, Dublin Philos. Mag. J. Sci._ 34, 15–24 (1867).  URL https://doi.org/10.1080/14786446708639836.  eprint https://doi.org/10.1080/14786446708639836.
  * [2] Skyrme, T. H. R.  A non-linear field theory.  _Proc. R. Soc. A_ 260, 127–138 (1961).  URL http://rspa.royalsocietypublishing.org/content/260/1300/127.abstract.
  * [3] Skyrme, T. H. R.  The origins of skymions.  _International Journal of Modern Physics A_ 3, 2745–2751 (1988).  URL http://inis.iaea.org/search/search.aspx?orig_q=RN:21029176.
  * [4] Bogolubskaya, A. & Bogolubsky, I.  Stationary topological solitons in the two-dimensional anisotropic heisenberg model with a skyrme term.  _Physics Letters A_ 136, 485–488 (1989).
  * [5] Bogolyubskaya, A. A. & Bogolyubsky, I. L.  ON STATIONARY TOPOLOGICAL SOLITONS IN TWO-DIMENSIONAL ANISOTROPIC HEISENBERG MODEL.  _Lett. Math. Phys._ 19, 171–177 (1990).
  * [6] Leese, R., Peyrard, M. & Zakrzewski, W.  Soliton scatterings in some relativistic models in (2+1) dimensions.  _Nonlinearity_ 3, 773–807 (1990).
  * [7] Polyakov, A. M. & Belavin, A. A.  Metastable States of Two-Dimensional Isotropic Ferromagnets.  _JETP Lett._ 22, 245–248 (1975).
  * [8] Dzyaloshinsky, I.  A thermodynamic theory of weak ferromagnetism of antiferromagnetics.  _J. Phys. Chem. Solids_ 4, 241 (1958).  URL http://www.sciencedirect.com/science/article/pii/0022369758900763.
  * [9] Moriya, T.  New mechanism of anisotropic superexchange interaction.  _Phys. Rev. Lett._ 4, 228–230 (1960).  URL http://link.aps.org/doi/10.1103/PhysRevLett.4.228.
  * [10] Moriya, T.  Anisotropic superexchange interaction and weak ferromagnetism.  _Phys. Rev._ 120, 91–98 (1960).
  * [11] Bogdanov, A. & Hubert, A.  Thermodynamically stable magnetic vortex states in magnetic crystals.  _Journal of Magnetism and Magnetic Materials_ 138, 255 – 269 (1994).  URL http://www.sciencedirect.com/science/article/pii/0304885394900469.
  * [12] Mühlbauer, S. _et al._ Skyrmion lattice in a chiral magnet.  _Science_ 323, 915 (2009).  URL http://www.sciencemag.org/content/323/5916/915.
  * [13] Yu, X. Z. _et al._ Real-space observation of a two-dimensional skyrmion crystal.  _Nature_ 465, 901 (2010).  URL http://www.nature.com/nature/journal/v465/n7300/full/nature09124.html.
  * [14] Yu, X. Z. _et al._ Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe.  _Nature Materials_ 10, 106 (2011).  URL http://www.nature.com/nmat/journal/v10/n2/full/nmat2916.html.
  * [15] Seki, S., Yu, X. Z., Ishiwata, S. & Tokura, Y.  Observation of skyrmions in a multiferroic material.  _Science_ 336, 198 (2012).  URL http://www.sciencemag.org/content/336/6078/198.
  * [16] Adams, T. _et al._ Long-wavelength helimagnetic order and skyrmion lattice phase in $\mathrm{Cu_{2}OSeO_{3}}$.  _Phys. Rev. Lett._ 108, 237204 (2012).  URL http://link.aps.org/doi/10.1103/PhysRevLett.108.237204.
  * [17] Yu, X. _et al._ Magnetic stripes and skyrmions with helicity reversals.  _Proceedings of the National Academy of Sciences_ 109, 8856–8860 (2012).  URL https://www.pnas.org/content/109/23/8856.
  * [18] Yu, X. Z. _et al._ Biskyrmion states and their current-driven motion in a layered manganite.  _Nature Communications_ 5, 3198 (2014).  URL https://doi.org/10.1038/ncomms4198.
  * [19] Mallik, R., Sampathkumaran, E. V., Paulose, P. L., Sugawara, H. & Sato, H.  Magnetic anomalies in gd2pdsi3.  _Pramana - J. Phys._ 51, 505 (1998).
  * [20] Saha, S. R. _et al._ Magnetic anisotropy, first-order-like metamagnetic transitions, and large negative magnetoresistance in single-crystal Gd2PdSi3.  _Phys. Rev. B_ 60, 12162–12165 (1999).  URL https://link.aps.org/doi/10.1103/PhysRevB.60.12162.
  * [21] Kurumaji, T. _et al._ Skyrmion lattice with a giant topological hall effect in a frustrated triangular-lattice magnet.  _Science_ 365, 914–918 (2019).  URL https://science.sciencemag.org/content/365/6456/914.
  * [22] Chandragiri, V., Iyer, K. K. & Sampathkumaran, E. V.  Magnetic behavior of Gd3Ru4Al12, a layered compound with distorted kagomé net.  _Journal of Physics: Condensed Matter_ 28, 286002 (2016).  URL https://doi.org/10.1088%2F0953-8984%2F28%2F28%2F286002.
  * [23] Hirschberger, M. _et al._ Skyrmion phase and competing magnetic orders on a breathing kagomé lattice.  _Nature Communications_ 10, 5831 (2019).  URL https://doi.org/10.1038/s41467-019-13675-4.
  * [24] Okubo, T., Chung, S. & Kawamura, H.  Multiple-$q$ states and the skyrmion lattice of the triangular-lattice heisenberg antiferromagnet under magnetic fields.  _Phys. Rev. Lett._ 108, 017206 (2012).  URL https://link.aps.org/doi/10.1103/PhysRevLett.108.017206.
  * [25] Leonov, A. O. & Mostovoy, M.  Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet.  _Nature Communications_ 6, 8275 (2015).  URL https://doi.org/10.1038/ncomms9275.
  * [26] Lin, S.-Z. & Hayami, S.  Ginzburg-landau theory for skyrmions in inversion-symmetric magnets with competing interactions.  _Phys. Rev. B_ 93, 064430 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.93.064430.
  * [27] Hayami, S., Lin, S.-Z. & Batista, C. D.  Bubble and skyrmion crystals in frustrated magnets with easy-axis anisotropy.  _Phys. Rev. B_ 93, 184413 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.93.184413.
  * [28] Batista, C. D., Lin, S.-Z., Hayami, S. & Kamiya, Y.  Frustration and chiral orderings in correlated electron systems.  _Reports on Progress in Physics_ 79, 084504 (2016).  URL https://doi.org/10.1088%2F0034-4885%2F79%2F8%2F084504.
  * [29] Wang, Z., Su, Y., Lin, S.-Z. & Batista, C. D.  Skyrmion crystal from RKKY interaction mediated by 2D electron gas.  _Phys. Rev. Lett._ 124, 207201 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevLett.124.207201.
  * [30] Hayami, S. & Motome, Y.  Topological spin crystals by itinerant frustration.  _J. Phys.: Condens. Matter_ 33, 443001 (2021).  URL https://doi.org/10.1088/1361-648x/ac1a30.
  * [31] Golo, V. L. & Perelomov, A. M.  Solution of the Duality Equations for the Two-Dimensional SU(N) Invariant Chiral Model.  _Phys. Lett. B_ 79, 112–113 (1978).
  * [32] D’Adda, A., Luscher, M. & Di Vecchia, P.  A 1/n Expandable Series of Nonlinear Sigma Models with Instantons.  _Nucl. Phys. B_ 146, 63–76 (1978).
  * [33] Din, A. M. & Zakrzewski, W.  General classical solutions in the cp${}^{n}-1$ model.  _Nuclear Physics_ 174, 397–406 (1980).
  * [34] Ferreira, L. A. & Klimas, P.  Exact vortex solutions in a cpnskyrme-faddeev type model.  _Journal of High Energy Physics_ 2010, 8 (2010).  URL https://doi.org/10.1007/JHEP10(2010)008.
  * [35] Amari, Y., Klimas, P., Sawado, N. & Tamaki, Y.  Potentials and the vortex solutions in the $CP^{N}$ Skyrme-Faddeev model.  _Phys. Rev. D_ 92, 045007 (2015).  eprint 1504.02848.
  * [36] Akagi, Y., Amari, Y., Sawado, N. & Shnir, Y.  Isolated skyrmions in the ${CP}^{2}$ nonlinear sigma model with a dzyaloshinskii-moriya type interaction.  _Phys. Rev. D_ 103, 065008 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevD.103.065008.
  * [37] Gnutzmann, S. & Kus, M.  Coherent states and the classical limit on irreducible representations.  _Journal of Physics A: Mathematical and General_ 31, 9871–9896 (1998).  URL https://doi.org/10.1088/0305-4470/31/49/011.
  * [38] Zhang, H. & Batista, C. D.  Classical spin dynamics based on SU(N) coherent states.  _Phys. Rev. B_ 104, 104409 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.104.104409.
  * [39] Papanicolaou, N.  Unusual phases in quantum spin-1 systems.  _Nuclear Physics B_ 305, 367–395 (1988).  URL https://www.sciencedirect.com/science/article/pii/0550321388900739.
  * [40] Batista, C. D. & Ortiz, G.  Algebraic approach to interacting quantum systems.  _Advances in Physics_ 53, 1–82 (2004).  URL https://doi.org/10.1080/00018730310001642086.  eprint https://doi.org/10.1080/00018730310001642086.
  * [41] Zapf, V. S. _et al._ Bose-Einstein condensation of $s=1$ nickel spin degrees of freedom in NiCl${}_{2}\mathrm{\text{$-$}}4\mathrm{SC}$(NH${}_{2}{)}_{2}$.  _Phys. Rev. Lett._ 96, 077204 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevLett.96.077204.
  * [42] Läuchli, A., Mila, F. & Penc, K.  Quadrupolar Phases of the ${{S}}=1$ Bilinear-Biquadratic Heisenberg Model on the Triangular Lattice.  _Phys. Rev. Lett._ 97, 087205 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevLett.97.087205.
  * [43] Muniz, R. A., Kato, Y. & Batista, C. D.  Generalized spin-wave theory: Application to the bilinear–biquadratic model.  _Prog. Theor. Exp. Phys._ 2014, 083I01 (2014).  URL https://doi.org/10.1093/ptep/ptu109.
  * [44] Our definitions for $\mathaccentV{hat}05E{S}^{x}$ and $\mathaccentV{hat}05E{S}^{z}$ differ from these two in Ref [36] by a minus sign.
  * [45] See Supplemental Material.
  * [46] Collins, M. F. & Petrenko, O. A.  Review/Synthèse: Triangular antiferromagnets.  _Can. J. Phys._ 75, 605–655 (1997).  URL https://cdnsciencepub.com/doi/10.1139/p97-007.
  * [47] McGuire, M. A.  Crystal and magnetic structures in layered, transition metal dihalides and trihalides.  _Crystals_ 7 (2017).  URL https://www.mdpi.com/2073-4352/7/5/121.
  * [48] Botana, A. S. & Norman, M. R.  Electronic structure and magnetism of transition metal dihalides: Bulk to monolayer.  _Phys. Rev. Materials_ 3, 044001 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevMaterials.3.044001.
  * [49] Nakatsuji, S. _et al._ Spin Disorder on a Triangular Lattice.  _Science_ 309, 1697–1700 (2005).  URL https://www.science.org/doi/full/10.1126/science.1114727.
  * [50] Cheng, J. G. _et al._ High-Pressure Sequence of Ba3NiSb2O9 Structural Phases: New ${S}=1$ Quantum Spin Liquids Based on Ni2+.  _Phys. Rev. Lett._ 107, 197204 (2011).  URL https://link.aps.org/doi/10.1103/PhysRevLett.107.197204.
  * [51] Li, N. _et al._ Quantum spin state transitions in the spin-1 equilateral triangular lattice antiferromagnet Na${}_{2}\mathrm{BaNi}{({\mathrm{PO}}_{4})}_{2}$.  _Phys. Rev. B_ 104, 104403 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevB.104.104403.
  * [52] Amoroso, D., Barone, P. & Picozzi, S.  Spontaneous skyrmionic lattice from anisotropic symmetric exchange in a ni-halide monolayer.  _Nature Communications_ 11, 5784 (2020).  URL https://doi.org/10.1038/s41467-020-19535-w.
  * [53] Bai, X. _et al._ Hybridized quadrupolar excitations in the spin-anisotropic frustrated magnet FeI2.  _Nature Physics_ 17, 467–472 (2021).  URL https://doi.org/10.1038/s41567-020-01110-1.
  * [54] Legros, A. _et al._ Observation of 4- and 6-magnon bound states in the spin-anisotropic frustrated antiferromagnet FeI2.  _Phys. Rev. Lett._ 127, 267201 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.127.267201.
  * [55] Kurita, N. & Tanaka, H.  Magnetic-field- and pressure-induced quantum phase transition in CsFeCl3 proved via magnetization measurements.  _Phys. Rev. B_ 94, 104409 (2016).  URL https://link.aps.org/doi/10.1103/PhysRevB.94.104409.
  * [56] Chung, W. C., de Hond, J., Xiang, J., Cruz-Colón, E. & Ketterle, W.  Tunable Single-Ion Anisotropy in Spin-1 Models Realized with Ultracold Atoms.  _Phys. Rev. Lett._ 126, 163203 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevLett.126.163203.
  * [57] Johnson, S. G.  The NLopt nonlinear-optimization package, http://github.com/stevengj/nlopt.