11institutetext: I. Physikalisches Institut, Universität zu Köln, Zülpicher
Str. 77, 50937 Köln, Germany
11email<EMAIL_ADDRESS>22institutetext: Dept. de Astrofísica
Molecular, Instituto de Física Fundamental (IFF-CSIC), Serrano 121, 28006
Madrid, Spain 33institutetext: Centro de Desarrollos Tecnológicos,
Observatorio de Yebes (IGN), 19141 Yebes, Guadalajara, Spain 44institutetext:
Observatorio Astronómico Nacional (OAN, IGN), Madrid, Spain 55institutetext:
Institut des Sciences Moléculaires (ISM), CNRS, Univ. Bordeaux, 351 cours de
la Libération, 33400, Talence, France 66institutetext: Laboratoire
d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-
Hilaire, 33615 Pessac, France 77institutetext: Instituto de Estructura de la
Materia, IEM-CSIC, Serrano 123, 28006 Madrid, Spain

# Discovery of H2CCCH+ in TMC-1 ††thanks: Based on observations carried out
with the Yebes 40m telescope (projects 19A003, 20A014, 20D15, and 21A011) and
the Institut de Radioastronomie Millimétrique (IRAM) 30m telescope. The 40m
radiotelescope at Yebes Observatory is operated by the Spanish Geographic
Institute (IGN, Ministerio de Transportes, Movilidad y Agenda Urbana). IRAM is
supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).

W. G. D. P. Silva 11 J. Cernicharo 22 S. Schlemmer 11 N. Marcelino 3344 J.-C.
Loison 55 M. Agúndez 22 D. Gupta 11 V. Wakelam 66 S. Thorwirth 11 C. Cabezas
22 B. Tercero 3344 J. L. Doménech 77 R. Fuentetaja 22 W.-J. Kim 11 P. de
Vicente 33 O. Asvany 11

(Received X, 2023; accepted X, 2023)

Based on a novel laboratory method, 14 mm-wave lines of the molecular ion
H2CCCH+ have been measured in high resolution, and the spectroscopic constants
of this asymmetric rotor determined with high accuracy. Using the Yebes 40 m
and IRAM 30 m radio telescopes, we detect four lines of H2CCCH+ towards the
cold dense core TMC-1. With a dipole moment of about 0.55 Debye obtained from
high-level ab initio calculations, we derive a column density of
5.4$\pm$1$\times$1011 cm-2 and 1.6$\pm$0.5$\times$1011 cm-2 for the ortho and
para species, respectively, and an abundance ratio N(H2CCC)/N(H2CCCH+)=
2.8$\pm$0.7. The chemistry of H2CCCH+ is modelled using the most recent
chemical network for the reactions involving the formation of H2CCCH+. We find
a reasonable agreement between model predictions and observations, and new
insights into the chemistry of C3 bearing species in TMC-1 are obtained.

###### Key Words.:

astrochemistry – ISM: molecules – ISM: individual objects: TMC-1 – line:
identification – molecular data – methods: laboratory: molecular

##  1 Introduction

Molecular ions are important intermediates in the chemistry of the
interstellar medium (ISM). These charged species can rapidly react with
neutral partners or recombine with electrons to form other ionic and neutral
molecules under astrophysical conditions (Larsson et al., 2012; Agúndez &
Wakelam, 2013). Despite their fundamental role in astrochemistry, many ions
remain elusive mainly due to their highly reactive character and lack of
accurate laboratory data to support astronomical detections (McGuire et al.,
2020). Consequently, many ions in astrochemical models and theories await
confirmation through spectroscopic detection in the ISM.

One example is the formation of the hydrocarbons C3H2 and C3H, whose both
cyclic (c-C3H2 and c-C3H) and linear (l-C3H2 and l-C3H, i.e. H2CCC and HCCC,
respectively) isomers were detected in the ISM (Thaddeus et al., 1985a, b;
Yamamoto et al., 1987; Cernicharo et al., 1991), and even deuterated versions
were observed (Bell et al., 1986; Spezzano et al., 2013, 2016; Agúndez et al.,
2019). The synthesis of the cyclic and linear forms of C3H2 and C3H is thought
to occur via the dissociative recombination of the respective isomers of C3H3+
with electrons, c-C3H3+ and H2CCCH+ (Maluendes et al., 1993). In turn, both
isomers of C3H3+ would be produced through the radiative association of C3H+
and H2 (Savić & Gerlich, 2005). The proof of this chemical pathway for the
cyclic variants is difficult because c-C3H3+ is a symmetric molecule and can
only be detected based on its rovibrational fingerprints in the infrared (Zhao
et al., 2014), which might be feasible with the James Webb Space Telescope
(JWST) in the near future. While the singly-deuterated version c-C3H2D+ could
be probed by radio astronomy, it has a predicted low dipole moment and low
column densities (Gupta et al., 2023). This leaves only H2CCCH+ as a good
candidate for radio astronomical searches.

In this Letter, based on a novel experimental method, we report first
laboratory mm-wave data of H2CCCH+ and its radio-astronomical detection
towards the cold dark core TMC-1. We derive its column density towards TMC-1
and discuss these results in the context of state-of-the-art chemical models.

##  2 Laboratory Work

H2CCCH+ is a closed-shell, planar and near-prolate asymmetric top molecular
ion (see sketch in Fig. 1). H2CCCH+ ions were generated in the Cologne
laboratory in a storage ion source via electron impact ionization ($E_{e}$
$\approx$30 eV) of the precursor gas allene (C3H4). By applying a novel trap-
based technique called leak-out-spectroscopy (LOS, Schmid et al., 2022) in the
cryogenic ion trap machine COLTRAP (Asvany et al., 2010, 2014), the
vibrational bands $\nu_{1}$ and $\nu_{3}$+$\nu_{5}$ were measured in the range
3180 - 3240 cm-1 in high resolution. The vibrational measurements, whose
details will be described in a forthcoming publication, enabled the ground
state spectroscopic parameters of H2CCCH+ to be determined. Subsequently, pure
rotational lines were detected using a vibrational-rotational double resonance
(DR) method. Such methods have been reviewed by Asvany & Schlemmer (2021), and
the particular scheme involving LOS has only recently been demonstrated by
Asvany et al. (2023). An example measurement for H2CCCH+ is shown in Fig. 1.

DR spectra were recorded in multiple individual measurements in which the mm-
wave frequency (blue arrow in Fig. 1) was stepped in an up-and-down manner
several times. Selected rovibrational lines from the $\nu_{1}$ or the
$\nu_{3}+\nu_{5}$ combination band were used for the IR excitation (red arrow
in Fig. 1). The frequency steps of the mm-wave radiation were kept constant in
individual experiments, and varied between 3 and 50 kHz (the larger steps
typically used for searching new lines). The spectroscopic data were
normalized employing a frequency switching procedure, i.e., by dividing the
H2CCCH+ counts monitored while scanning the spectral range of interest by the
counts at an off-resonant mm-wave reference frequency. Therefore, the baseline
in Fig. 1 is close to unity. The on-resonance signal enhancement is on the
order of 15 %. Transition frequencies were determined by adjusting the
parameters of an appropriate line shape function (typically a Gaussian) to the
experimental spectrum in a least-squares procedure. In total, 14 rotational
lines were detected in the laboratory and are summarised in Table 1. The
frequencies and their uncertainties in the table result from the weighted
average of several (up to eleven) independent line-center determinations for
each transition.

Figure 1:  Pure rotational transition ($J^{\prime}_{Ka^{\prime}Kc^{\prime}}$
$\leftarrow$ $J^{\prime\prime}_{Ka^{\prime\prime}Kc^{\prime\prime}}$ =
$8_{18}$ $\leftarrow$ $7_{17}$) of H2CCCH+ recorded using the double resonance
spectroscopic scheme. For this measurement, the IR laser frequency (red arrow)
was kept fixed on resonance with the $9_{19}$ $\leftarrow$ $8_{18}$
rovibrational transition within the $\nu_{1}$ band.  Table 1:  Ground state
rotational transition frequencies of H2CCCH+, determined by laboratory and
astronomical detections. The uncertainty in the last digits of each line is
provided in parentheses.

$J^{\prime}_{Ka^{\prime}Kc^{\prime}}$ $\leftarrow$ | $J^{\prime\prime}_{Ka^{\prime\prime}Kc^{\prime\prime}}$ | frequency / MHz | obs-calc / kHz | source
---|---|---|---|---
$2_{12}$ $\leftarrow$ | $1_{11}$ | 37702 | .627(10) | 0 | .3 |  astro a
$2_{02}$ $\leftarrow$ | $1_{01}$ | 38037 | .044(10) | 7 | .0 |  astro a
$2_{11}$ $\leftarrow$ | $1_{10}$ | 38368 | .581(20) | 33 | .0 |  astro a
$5_{15}$ $\leftarrow$ | $4_{14}$ | 94254 | .055(4) | 7 | .7 | lab
$5_{05}$ $\leftarrow$ | $4_{04}$ | 95086 | .000(10) |  $-$12 | .8 | lab
$5_{14}$ $\leftarrow$ | $4_{13}$ | 95918 | .753(3) | 7 | .9 | lab
$6_{16}$ $\leftarrow$ | $5_{15}$ | 113103 | .272(4) |  $-$2 | .6 | lab
$6_{06}$ $\leftarrow$ | $5_{05}$ | 114099 | .086(6) | 6 | .0 | lab
$6_{15}$ $\leftarrow$ | $5_{14}$ | 115100 | .840(3) |  $-$2 | .7 | lab
$7_{17}$ $\leftarrow$ | $6_{16}$ | 131951 | .635(3) |  $-$5 | .4 | lab
$7_{07}$ $\leftarrow$ | $6_{06}$ | 133109 | .888(6) |  $-$3 | .9 | lab
$7_{16}$ $\leftarrow$ | $6_{15}$ | 134282 | .033(4) |  $-$5 | .1 | lab
$8_{18}$ $\leftarrow$ | $7_{17}$ | 150799 | .008(4) | 5 | .3 | lab
$8_{08}$ $\leftarrow$ | $7_{07}$ | 152118 | .074(5) | 0 | .8 | lab
$8_{17}$ $\leftarrow$ | $7_{16}$ | 153462 | .176(3) |  $-$3 | .7 | lab
$10_{110}$ $\leftarrow$ | $9_{19}$ | 188490 | .152(12) | 1 | .6 | lab
$10_{19}$ $\leftarrow$ | $9_{18}$ | 191818 | .719(10) | 27 | .3 | lab
111
$a$$a$footnotetext: Observed frequency assuming a Local Standard of Rest
velocity (vLSR) of 5.83 km s-1 (Cernicharo et al., 2020a).

The fit of the assigned lines (lab + astro) was carried out using Watson’s
$S$-reduced Hamiltonian in the _I r_ representation as implemented in
Western’s PGOPHER program (Western, 2017). The resulting spectroscopic
parameters are given in Table 2. As H2CCCH+ is a near-prolate asymmetric top
($\kappa=-0.9976$) with an $a$-type spectrum, the $A_{0}$ rotational constant
is not well constrained from experiment. Overall, the experimental rotational
and centrifugal distortion constants compare favorably with those calculated
by Huang et al. (2011) and very well with the best estimate values obtained in
this study (given in the last three columns of Table 2). The obtained obs-calc
values for the rotational lines are given in Table 1. In total, the weighted
rms of the fit is on the order of 1.4, indicating somewhat optimistic
uncertainties of our measurements.

Table 2:  Spectroscopic parameters of H2CCCH+ (in MHz) obtained by fitting
Watson’s S-reduced Hamiltonian to the rotational transitions from Table 1.
Western’s PGOPHER program (Western, 2017) has been used for the fit. The
experimental uncertainties are reported in parentheses.

| Exp. | Calculated
---|---|---
Parameter | This study | Huang et al. (2011) |  This studya |  This study, BEb
$A_{0}$ | 281856 | .(247) | 281911 | .9 | 283318 | .43 | 282182 | .6
$B_{0}$ | 9675 | .841(1) | 9580 | .2 | 9670 | .18 | 9675 | .97
$C_{0}$ | 9342 | .877(1) | 9251 | .9 | 9337 | .82 | 9343 | .14
$D_{J}\times 10^{3}$ | 2 | .942(4) | 3 | . | 2 | .68 | 3 | .06
$D_{JK}$ | 0 | .4388(9) | 0 | .479 | 0 | .479 | 0 | .433
$D_{K}$ | 20 |  .678c | 20 | .862 | 20 | .678 | 20 | .678
$d_{1}\times 10^{3}$ |  $-$0 | .121(3) | 0 | . |  $-$0 | .097 |  $-$0 | .103
$d_{2}\times 10^{3}$ |  $-$0 |  .051c | 0 | . |  $-$0 | .037 |  $-$0 | .051
$\mu_{A}$ / Debye | $\cdots$ | 0 |  .524 d | 0 | .55 | $\cdots$
222
$a$$a$footnotetext: From equilibrium structure and force field calculated at
the ae-CCSD(T)/cc-pwCVQZ and fc-CCSD(T)/ANO1 levels, respectively.
$b$$b$footnotetext: Best estimate (BE) values, obtained through scaling of
calculated values via isoelectronic propadienylidene, H2CCC, see text for
details.
$c$$c$footnotetext: Fixed to best estimate value. $d$$d$footnotetext: From
Huang & Lee (2011)

Figure 2: Lines of H2CCCH+ observed in the 31–115 GHz domain towards TMC-1.
The abscissa corresponds to the rest frequency in MHz. The ordinate is the
antenna temperature corrected for atmospheric and telescope losses in mK. The
spectral resolution is 38 kHz below 50 GHz and 48 kHz above. Line parameters
are given in Table 4. The red lines show the computed synthetic spectra (see
text).

##  3 Quantum Chemical Calculations

The H2CCCH+ molecular ion has been a subject of several quantum-chemical
investigations in the past (e.g., Botschwina et al., 1993, 2011; Huang et al.,
2011; Marimuthu et al., 2020, and references therein). In the present study,
complementary high-level calculations were performed at the CCSD(T) level of
theory (Raghavachari et al., 1989) together with correlation consistent
(augmented) polarized weighted core-valence basis sets (Kendall et al., 1992;
Peterson & Dunning, 2002) as well as atomic natural orbital basis sets (Almlöf
& Taylor, 1987). All calculations were performed using the CFOUR program suite
(Matthews et al., 2020; Harding et al., 2008). Equilibrium rotational
constants were calculated at the all-electron (ae-)CCSD(T)/cc-pwCVQZ level of
theory that is known to yield molecular equilibrium structural parameters of
very high quality for molecules comprising first- and second-row elements
(e.g., Coriani et al., 2005).

Zero-point vibrational contributions
$\frac{1}{2}\sum_{i}\alpha_{i}^{A,B,C,calc}$ ($=\Delta A_{0},\Delta
B_{0},\Delta C_{0}$) to the equilibrium rotational constants and centrifugal
distortion parameters were calculated at the frozen core (fc-)CCSD(T)/ANO1
level. Best estimate (BE, Tables 2 and 3) rotational and centrifugal
distortion constants were finally obtained through empirical scaling of the
calculated rotational parameters using factors (i.e., the ratios
$X_{exp}/X_{calc}$ of a given parameter) derived from isoelectronic
propadienylidene, H2CCC, the pure rotational spectrum of which is known well
from previous study (Vrtilek et al., 1990). The technique of empirical scaling
using structurally closely related (isoelectronic) species known from
experiment may provide rotational parameters at a predictive power greatly
exceeding that of high-level calculations alone (see, e.g., Thorwirth et al.,
2005, 2020; Martinez, Jr. et al., 2013), and has been used recently to
identify species not studied in the laboratory using radio astronomy (see,
e.g., Cernicharo et al., 2020b).

The dipole moment of H2CCCH+ is not very large. At the ae-CCSD(T)/aug-cc-
pwCVQZ level of theory the (center-of-mass frame) equilibrium value $\mu_{e}$
amounts to 0.53 D in good agreement with earlier estimates. Zero-point
vibrational effects (fc-CCSD(T)/ANO1) have an almost negligible influence
resulting in $\mu_{0}=0.55$ D.

##  4 Observations

New receivers, built within the Nanocosmos333ERC grant
ERC-2013-Syg-610256-NANOCOSMOS.
https://nanocosmos.iff.csic.es/ project and installed at the Yebes 40 m radio
telescope, were used for the observations of TMC-1 ($\alpha_{J2000}=4^{\rm
h}41^{\rm m}41.9^{\rm s}$ and
$\delta_{J2000}=+25^{\circ}41^{\prime}27.0^{\prime\prime}$). The observations
of TMC-1 belong to the on-going QUIJOTE444Q-band Ultrasensitive Inspection
Journey to the Obscure TMC-1 Environment line survey (Cernicharo et al.,
2021a, 2023). A detailed description of the telescope, receivers, and backends
is given by Tercero et al. (2021). Briefly, the receiver consists of two cold
high electron mobility transistor amplifiers covering the 31.0-50.3 GHz band
with horizontal and vertical polarizations. The backends are $2\times 8\times
2.5$ GHz fast Fourier transform spectrometers with a spectral resolution of
38.15 kHz providing the coverage of the whole Q-band in both polarisations.

The observations, carried out during different observing runs, are performed
using the frequency-switching mode with a frequency throw of 10 MHz in the
very first observing runs, during November 2019 and February 2020, 8 MHz
during the observations of January-November 2021, and alternating these
frequency throws in the last observing runs between October 2021 and February
2023. The total on-source telescope time is 850 hours in each polarization
(385 and 465 hours for the 8 MHz and 10 MHz frequency throws, respectively).
The sensitivity of the QUIJOTE line survey varies between 0.17 and 0.25 mK in
the 31-50.3 GHz domain. The intensity scale used in this work, antenna
temperature ($T_{A}^{*}$), was calibrated using two absorbers at different
temperatures and the atmospheric transmission model ATM (Cernicharo, 1985;
Pardo et al., 2001). Calibration uncertainties have been adopted to be 10 %.
The beam efficiency of the Yebes 40 m telescope in the Q-band is given as a
function of frequency by $B_{\rm eff}$= 0.797 exp[$-$($\nu$(GHz)/71.1)2]. The
forward telescope efficiency is 0.97. The telescope beam size varies from
56.7′′ at 31 GHz to 35.6′′ at 49.5 GHz.

The data of TMC-1 taken with the IRAM 30 m telescope consist of a 3 mm line
survey obtained with the old ABCD receivers connected to an autocorrelator
that provided a spectral resolution of 40 kHz (Marcelino et al., 2007;
Cernicharo et al., 2012). Some additional high sensitivity frequency windows
observed in 2021 used the new 3 mm EMIR dual polarization receiver connected
to four fast Fourier transform spectrometers providing a spectral resolution
of 49 kHz (Agúndez et al., 2022; Cabezas et al., 2022a). All the observations
were performed using the frequency switching method. The final 3 mm line
survey has a sensitivity of 2-10 mK. However, at some selected frequencies the
sensitivity is as low as 0.6 mK.

##  5 Detection of H2CCCH+ in TMC-1

We searched for one para (202-101) and two ortho (212-111, 211-110) lines of
H2CCCH+ within the QUIJOTE line survey. The three lines are clearly detected
and are shown in Fig. 2. In the data at 3 mm we covered the frequencies of
three para and four ortho lines, with upper levels $J$= 4, 5, 6 and energies
below 30 K. Only one line, the $J=5_{14}-4_{13}$, falls in one of the high
sensitivity windows ($\sigma$= 0.6 mK) of our line survey, and it is also
clearly detected (see Fig. 2). Two other lines are within frequency ranges
with $\sigma$ below 2 mK, and are marginally detected. The derived line
parameters for all searched transitions of H2CCCH+ are given in Table 4. We
checked that the detected lines cannot be assigned to lines of other species
or isotopologues by exploring the spectral catalogues MADEX (Cernicharo,
2012), CDMS (Müller et al., 2005) and JPL (Pickett et al., 1998).

To estimate the column density of H2CCCH+, we considered the ortho and para
levels as belonging to two different species. For the dipole moment we use the
value of 0.55 D calculated in this work. No collisional rates are available
for this molecule. However, Khalifa et al. (2019) have computed the
collisional rates between He and H2CCC. This cumulenic species is
isoelectronic with our molecule and has a very similar structure. We therefore
adopted these rates, correcting for the abundance of He with respect to H2, to
estimate the excitation temperatures of the observed transitions of H2CCCH+.
Assuming a volume density of (1-3)$\times$104 cm-3 (Fossé et al., 2001; Lique
et al., 2006; Pratap et al., 1997), we derive excitation temperatures close to
10 K for the $J$= 2-1 lines, and $\sim$8-10 K for the lines in the 3 mm
domain, the largest value corresponding to $n$(H2)= 3$\times$104 cm-3. These
excitation temperatures are considerably larger than those obtained for H2CCC,
which are between 4 and 5 K, due to the larger dipole moment of this species
(4.1 versus 0.55 D). From the adopted rotational temperatures of 9 K, and
assuming a source of uniform brightness temperature with a radius of 40′′
(Fossé et al., 2001), we derive a column density for ortho-H2CCCH+ of
(5.4$\pm$1)$\times$1011 cm-2. For the para species the estimated column
density is (1.6$\pm$0.5)$\times$1011 cm-2, a value that is consistent with the
expected ortho to para ratio of 3/1. The computed synthetic spectra show an
excellent agreement with the observed line intensities (see Fig. 2). Hence,
the total column density of H2CCCH+ in TMC-1 is (7.0$\pm$1.5)$\times$1011
cm-2.

It is interesting to compare the abundance of H2CCCH+ to that of H2CCC. For
the latter species we detected, with an excellent signal to noise ratio, all
its ortho and para lines in the frequency range of our line surveys. The
derived line parameters are summarized in Table 5 and the lines are shown in
Fig. 4. The decline of the line intensity between the $J_{u}$= 2 and $J_{u}$=
5 lines is obvious in Fig. 4, which indicates that the lines are not
thermalized to the kinetic temperature of the cloud for this species. Using
the collisional rates of Khalifa et al. (2019), and adopting the same
assumptions on the source size than for H2CCCH+, we derive a column density
for the ortho and para species of (1.5$\pm$0.1)$\times$1012 and
(0.45$\pm$0.05)$\times$1012 cm-2, respectively. The best fit is obtained for a
density of $n$(H2)=8$\times$103 cm-3. The total column density of H2CCC is
(1.95$\pm$0.15)$\times$1012 cm-2 and the ortho to para ratio for this species
is 3.3$\pm$0.6. The H2CCC/H2CCCH+ abundance ratio is 2.8$\pm$0.7 which is on
the order of that found for C3O/HC3O+ (Cernicharo et al., 2020b), but much
smaller than the abundance ratio found in TMC-1 for other neutral species and
their protonated forms (Marcelino et al., 2020; Cernicharo et al., 2021b, c;
Cabezas et al., 2022b; Agúndez et al., 2022).

Figure 3: Continuous lines: Abundances of C3, C3H+, H2CCC and H2CCCH+ as a
function of time predicted by our model. The horizontal colored areas
represent the observations in TMC-1 assuming an uncertainty of 3 on the
measurements (Cernicharo et al. 2022 for C3H+ and this work for H2CCCH+).
Dashed lines: the same with a threefold decrease in H2 \+ C3H+ and a fivefold
increase in H2CCCH+ \+ e- rate coefficients. The vertical grey area represents
the values given by the most probable chemical age for TMC-1 given by the
better agreement between calculations and observations for 67 key species.

##  6 Discussion

To describe the chemistry of H2CCCH+, we used the Nautilus code (Ruaud et al.,
2016), a 3-phase (gas, dust grain ice surface, and dust grain ice mantle)
time-dependent chemical model with a chemical network for C3H${}_{x}^{+}$
species very similar to the one presented in Loison et al. (2017). To describe
the physical conditions in TMC-1, we use an homogeneous cloud with a density
equal to 2.5$\times$104 cm-3, a temperature equal to 10 K for both the gas and
the dust, a visual extinction of 30 mag and a cosmic-ray ionization rate of
1.3$\times$10-17 s-1. All elements are assumed to be initially in atomic form,
except for hydrogen, which is entirely molecular (Hincelin et al., 2011). The
calculated abundances relative to H2 for H2CCC and H2CCCH+, and also for C3
and C3H+, which are strongly linked to H2CCCH+, are shown in Fig. 3.

As can be seen in Fig. 3, the H2CCC and H2CCCH+ abundances observed are
relatively well reproduced by the model for a relatively early molecular cloud
age around 2$\times$105 years with, however, a smaller H2CCC/H2CCCH+ ratio
than the one observed. Looking in more detail at the chemistry of H2CCCH+ and
H2CCC (see Fig. 3 of Loison et al. 2017), it appears that H2CCC is a product
of H2CCCH+ (and also of c-C3H${}_{3}^{+}$) but that the flow of protonation of
H2CCC toward H2CCCH+ is a very minor pathway for the formation of H2CCCH+
which is almost essentially produced by the reaction C3H+ \+ H2. This inverted
link between H2CCC and H2CCCH+ explains the unusually high MH+/M ratio,
compared to those cases in which the protonated form comes from the
protonation of the neutral form (Agúndez et al., 2022).

Considering the link between C3H+, H2CCCH+ and C3, it is interesting to see if
the observations of H2CCCH+ (this work) combined to the observation of C3H+
(Cernicharo et al., 2022) allow us to estimate the abundance of C3, which in
our dense cloud model is the second carbon reservoir, accounting for up to 15%
of carbon. If C3 does not react with atomic oxygen, as calculated by Woon &
Herbst (1996), the protonation reactions of C3 producing C3H+ are by far the
main reactions of destruction of C3 and of production of C3H+. As these
protonation reactions control the destruction of C3, the uncertainties on the
rates affect the abundance of C3, but does not change the flux of these
reactions. The underestimation of C3H+ in the model is therefore not related
to these protonation rates, but to the rate of the reaction C3H+ \+ H2, which
controls the destruction of C3H+ (Maluendes et al., 1993; Savić & Gerlich,
2005). A decrease in the rate of the reaction C3H+ \+ H2 at 10 K by a factor
of 3 allows us to reproduce the abundance of C3H+ observed by Cernicharo et
al. (2022), as shown in Fig. 3. This change in the rate coefficient does not
affect the flux of the C3H+ \+ H2 reaction and therefore does not affect the
abundance of H2CCCH+, as long as the branching ratios to H2CCCH+ and
c-C3H${}_{3}^{+}$ are not varied. Indeed, since H2CCCH+ is mainly destroyed by
the reaction with electrons, its abundance is controlled by the rate of this
reaction, which is known only over the temperature range 172-489 K (McLain et
al., 2005) with a temperature dependency inconsistent with theory. An increase
in this rate at 10 K by a factor of 5, which is not impossible given the
uncertainties, allows us to reproduce the observation for H2CCCH+ with a ratio
between H2CCC/H2CCCH+ very close to the observed one. Considering the
uncertainties of the different chemical reactions linking C3 to C3H+ and
H2CCCH+, the observations of C3H+ (Cernicharo et al., 2022) and H2CCCH+ (this
work) validate the chemical scheme controlling the formation of cyclic and
linear C3H and C3H2 and the large gas-phase abundance of C3 in TMC-1, around
10-6 relative to H2.

It would also be very interesting to know the abundance of the more stable
cyclic isomer c-C3H${}_{3}^{+}$, which in the chemical model is predicted to
be slightly more abundant (a factor of three) than H2CCCH+. This species has
no dipole moment, but its deuterated version, c-C3H2D+, has a low dipole
moment of 0.225 D (Huang & Lee, 2011) and its rotational spectrum has been
recently measured in the 90-230 GHz frequency range in Cologne (Gupta et al.,
2023). We searched for c-C3H2D+ in the QUIJOTE line survey but at the current
level of sensitivity, this species is not detected and we derive an upper
limit to its column density of 4$\times$1012 cm-3. Assuming that the
c-C3H${}_{3}^{+}$/c-C3H2D+ ratio is 10, as found for the analog case of CH3CCH
with three equivalent H nuclei (Cabezas et al., 2021), the column density of
c-C3H${}_{3}^{+}$ is $<$4$\times$1013 cm-2, which is not very meaningful.

###### Acknowledgements.

The work has been supported by an ERC advanced grant (MissIons: 101020583),
the Deutsche Forschungsgemeinschaft (DFG) via SFB 956 (project ID 184018867),
sub-project B2, and the Gerätezentrum ”Cologne Center for Terahertz
Spectroscopy” (DFG SCHL 341/15-1). W.G.D.P.S. thanks the Alexander von
Humboldt foundation for funding through a postdoctoral fellowship. We also
acknowledge the support from the MICINN projects PID2020-113084GB-I00 and
PID2019-106110GB-I00, the CSIC project ILINK+ LINKA20353, the ERC grant
ERC-2013-Syg610256-NANOCOSMOS, and the Regional Computing Center of the
Universität zu Köln (RRZK) for providing computing time on the DFG-funded high
performance computing system CHEOPS.

## References

  * Agúndez et al. (2022) Agúndez, M., Cabezas, C., Marcelino, N., et al. 2022, A.&A., 659, L9
  * Agúndez et al. (2019) Agúndez, M., Marcelino, N., Cernicharo, J., Roueff, E., & Tafalla, M. 2019, A&A, 625, A147
  * Agúndez & Wakelam (2013) Agúndez, M. & Wakelam, V. 2013, Chemical Reviews, 113, 8710
  * Almlöf & Taylor (1987) Almlöf, J. & Taylor, P. R. 1987, J. Chem. Phys., 86, 4070
  * Asvany et al. (2010) Asvany, O., Bielau, F., Moratschke, D., Krause, J., & Schlemmer, S. 2010, Rev. Sci. Instr., 81, 076102
  * Asvany et al. (2014) Asvany, O., Brünken, S., Kluge, L., & Schlemmer, S. 2014, Appl. Phys. B, 114, 203
  * Asvany & Schlemmer (2021) Asvany, O. & Schlemmer, S. 2021, Phys. Chem. Chem. Phys., 23, 26602
  * Asvany et al. (2023) Asvany, O., Thorwirth, S., Schmid, P. C., Salomon, T., & Schlemmer, S. 2023, Phys. Chem. Chem. Phys., submitted
  * Bell et al. (1986) Bell, M. B., Feldman, P. A., Matthews, H. E., & Avery, L. W. 1986, ApJ, 311, L89
  * Botschwina et al. (1993) Botschwina, P., Horn, M., Flügge, J., & Seeger, S. 1993, J. Chem. Soc. Faraday Trans., 89, 2219
  * Botschwina et al. (2011) Botschwina, P., Oswald, R., & Rauhut, G. 2011, PCCP, 13, 7921
  * Cabezas et al. (2022a) Cabezas, C., Agúndez, M., Marcelino, N., et al. 2022a, A.&A., 657, L4
  * Cabezas et al. (2022b) Cabezas, C., Agúndez, M., Marcelino, N., et al. 2022b, A.&A., 659, L8
  * Cabezas et al. (2021) Cabezas, C., Endo, Y., Roueff, E., et al. 2021, A&A, 646, L1
  * Cernicharo (1985) Cernicharo, J. 1985, Iram internal report
  * Cernicharo (2012) Cernicharo, J. 2012, ECLA2011; Proc. of the European Conference on Laboratory Astrophysics, EAS Publications Series, 2012, Ed.: C. Stehl, C. Joblin, & L. d’Hendecourt (Cambridge: Cambridge Univ. Press), 251, https://nanocosmos.iff.csic.es/?page$\\_$id=1619
  * Cernicharo et al. (2021a) Cernicharo, J., Agúndez, M., Kaiser, R., et al. 2021a, A.&A., 652, L9
  * Cernicharo et al. (2022) Cernicharo, J., Agúndez, M., Cabezas, C., et al. 2022, A.&A., 657, L16
  * Cernicharo et al. (2021b) Cernicharo, J., Cabezas, C., Endo, Y., et al. 2021b, A.&A., 646, L3
  * Cernicharo et al. (2021c) Cernicharo, J., Cabezas, C., S., B., et al. 2021c, A.&A., 646, L7
  * Cernicharo et al. (1991) Cernicharo, J., Gottlieb, C. A., Guelin, M., et al. 1991, ApJ, 368, L39
  * Cernicharo et al. (2020a) Cernicharo, J., Marcelino, N., Agúndez, M., et al. 2020a, A.&A., 642, L8
  * Cernicharo et al. (2020b) Cernicharo, J., Marcelino, N., Agúndez, M., et al. 2020b, A.&A., 642, L17
  * Cernicharo et al. (2012) Cernicharo, J., Marcelino, N., Roueff, E., et al. 2012, ApJ, 759, L43
  * Cernicharo et al. (2023) Cernicharo, J., Pardo, J., Cabezas, C., et al. 2023, A.&A., 670, L19
  * Coriani et al. (2005) Coriani, S., Marcheson, D., Gauss, J., et al. 2005, J. Chem. Phys., 123, 184107
  * Fossé et al. (2001) Fossé, D., Cernicharo, J., Gerin, M., & Cox, P. 2001, ApJ, 552, 168
  * Gupta et al. (2023) Gupta, D., Dias de Paiva Silva, W. G., Doménech, J. L., et al. 2023, Faraday Discuss., http://dx.doi.org/10.1039/D3FD00068K
  * Harding et al. (2008) Harding, M. E., Metzroth, T., Gauss, J., & Auer, A. A. 2008, J. Chem. Theory Comput., 4, 64
  * Hincelin et al. (2011) Hincelin, U., Wakelam, V., Hersant, F., et al. 2011, A&A, 530, 61
  * Huang & Lee (2011) Huang, X. & Lee, T. J. 2011, Astrophys. J., 736, 33
  * Huang et al. (2011) Huang, X., Taylor, P. R., & Lee, T. J. 2011, J. Phys. Chem. A, 115, 5005
  * Kendall et al. (1992) Kendall, R. A., Dunning, T. H., & Harrison, R. J. 1992, J. Chem. Phys, 96, 6796
  * Khalifa et al. (2019) Khalifa, M. B., Sahnoun, E., Wiesenfeld, L., et al. 2019, PCCP, 21, 1443
  * Larsson et al. (2012) Larsson, M., Geppert, W. D., & Nyman, G. 2012, Rep. Prog. Phys., 75, 066901
  * Lique et al. (2006) Lique, F. ., Cernicharo, J., & Cox, P. 2006, ApJ, 653, 1342
  * Loison et al. (2017) Loison, J.-C., Agúndez, M., Wakelam, V., et al. 2017, MNRAS, 470, 4075
  * Maluendes et al. (1993) Maluendes, S. A., McLean, A. D., & Herbst, E. 1993, ApJ, 417, 181
  * Marcelino et al. (2020) Marcelino, N., Agúndez, M., Tercero, B., et al. 2020, A.&A., 643, L6
  * Marcelino et al. (2007) Marcelino, N., Cernicharo, J., Agúndez, M., et al. 2007, ApJ, 665, L127
  * Marimuthu et al. (2020) Marimuthu, A. N., Sundelin, D., Thorwirth, S., et al. 2020, J. Mol. Spectrosc., 374, 111377
  * Martinez, Jr. et al. (2013) Martinez, Jr., O., Lattanzi, V., Thorwirth, S., & McCarthy, M. C. 2013, J. Chem. Phys., 138, 094316
  * Matthews et al. (2020) Matthews, D. A., Cheng, L., Harding, M. E., et al. 2020, J. Chem. Phys., 152, 214108
  * McGuire et al. (2020) McGuire, B. A., Asvany, O., Brünken, S., & Schlemmer, S. 2020, Nature Reviews Physics, 2, 402
  * McLain et al. (2005) McLain, J. L., Poterya, V., Molek, C. D., et al. 2005, J. Phys. Chem. A, 109, 5119
  * Müller et al. (2005) Müller, H., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struc., 742, 215
  * Pardo et al. (2001) Pardo, J., Cernicharo, J., & Serabyn, E. 2001, IEEE Trans. Antennas and Propagation, 49, 12
  * Peterson & Dunning (2002) Peterson, K. A. & Dunning, T. H. 2002, J. Chem. Phys., 117, 10548
  * Pickett et al. (1998) Pickett, H., Poynter, R., Cohen, E., et al. 1998, J. Quant. Spectrosc. Radiat. Transfer, 60, 883
  * Pratap et al. (1997) Pratap, P., Dickens, J., Snell, R., et al. 1997, ApJ, 486, 862
  * Raghavachari et al. (1989) Raghavachari, K., Trucks, G. W., Pople, J. A., & Head-Gordon, M. 1989, Chem. Phys. Lett., 157, 479
  * Ruaud et al. (2016) Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756
  * Savić & Gerlich (2005) Savić, I. & Gerlich, D. 2005, PCCP, 7, 1026
  * Schmid et al. (2022) Schmid, P. C., Asvany, O., Salomon, T., Thorwirth, S., & Schlemmer, S. 2022, J. Phys. Chem. A, 126, 8111, pMID: 36278898
  * Spezzano et al. (2013) Spezzano, S., Brünken, S., Schilke, P., et al. 2013, ApJ, 769, L19
  * Spezzano et al. (2016) Spezzano, S., Gupta, H., Brünken, S., et al. 2016, A&A, 586, A110
  * Tercero et al. (2021) Tercero, F., López-Pérez, J. A., Gallego, J. D., et al. 2021, A.&A., A37, 552
  * Thaddeus et al. (1985a) Thaddeus, P., Gottlieb, C. A., Hjalmarson, A., et al. 1985a, ApJ, 294, L49
  * Thaddeus et al. (1985b) Thaddeus, P., Vrtilek, J. M., & Gottlieb, C. A. 1985b, ApJ, 299, L63
  * Thorwirth et al. (2020) Thorwirth, S., Harding, M. E., Asvany, O., et al. 2020, Mol. Phys., 118, e1776409
  * Thorwirth et al. (2005) Thorwirth, S., McCarthy, M. C., Dudek, J. B., & Thaddeus, P. 2005, J. Chem. Phys., 122, 184308
  * Vrtilek et al. (1990) Vrtilek, J. M., Gottlieb, C. A., Gottlieb, E. W., Killian, T. C., & Thaddeus, P. 1990, Astrophys. J., 364, L53
  * Western (2017) Western, C. M. 2017, J. Quant. Spectros. Rad. Transfer, 186, 221
  * Woon & Herbst (1996) Woon, D. E. & Herbst, E. 1996, ApJ, 465, 795
  * Yamamoto et al. (1987) Yamamoto, S., Saito, S., Ohishi, M., et al. 1987, ApJ, 322, L55
  * Zhao et al. (2014) Zhao, D., Doney, K. D., & Linnartz, H. 2014, Astrophys. J. Lett., 791, L28

##  Appendix A Structural calculations, internal coordinates

Bond lengths are given in Å, angles in degrees.

###  A.1 H2CCCH+

    H2C3H+, CCSD(T)/cc-pwCVQZ
    H
    C 1 r1
    X 2 rd 1 a90
    C 2 r2 3 a90 1 d180
    X 4 rd 2 a90 3 d0
    C 4 r3 5 a90 2 d180
    H 6 r4 4 a1 5 d0
    H 6 r4 4 a1 5 d180
    r1   =        1.073118130747218
    rd   =        1.000000818629780
    a90  =       90.000000000000000
    r2   =        1.228490857165352
    d180 =      180.000000000000000
    d0   =        0.000000000000000
    r3   =        1.346565030012042
    r4   =        1.086025065701484
    a1   =      120.382476682815607

    H2C3H+, CCSD(T)/ANO1
    H
    C 1 r1
    X 2 rd 1 a90
    C 2 r2 3 a90 1 d180
    X 4 rd 2 a90 3 d0
    C 4 r3 5 a90 2 d180
    H 6 r4 4 a1 5 d0
    H 6 r4 4 a1 5 d180
    r1   =        1.074705266639338
    rd   =        1.000000204657382
    a90  =       90.000000000000000
    r2   =        1.234261558515133
    d180 =      180.000000000000000
    d0   =        0.000000000000000
    r3   =        1.351397370372584
    r4   =        1.088292080166718
    a1   =      120.373581553220575

###  A.2 H2CCC

    H2CCC, CCSD(T)/cc-pwCVQZ
    C
    C 1 r1
    X 2 rd 1 a90
    C 2 r2 3 a90 1 d180
    H 4 r3 2 a1 3 d0
    H 4 r3 2 a1 5 d180
    r1   =        1.287645153115967
    rd   =        1.000000000000000
    a90  =       90.000000000000000
    r2   =        1.327897213723068
    d180 =      180.000000000000000
    r3   =        1.083584615827809
    a1   =      121.271597504509145
    d0   =        0.000000000000000

    H2CCC, CCSD(T)/ANO1
    C
    C 1 r1
    X 2 rd 1 a90
    C 2 r2 3 a90 1 d180
    H 4 r3 2 a1 3 d0
    H 4 r3 2 a1 5 d180
    r1   =        1.294632885150227
    rd   =        1.000000204657382
    a90  =       90.000000000000000
    r2   =        1.333030062830918
    d180 =      180.000000000000000
    r3   =        1.085894615830395
    a1   =      121.302559724582778
    d0   =        0.000000000000000

Table 3:  Calculated and experimental spectroscopic parameters of H2CCC and
H2CCCH+ (in MHz). Equilibrium rotational constants calculated at the ae-
CCSD(T)/cc-pwCVQZ level of theory, zero-point vibrational corrections $\Delta
A_{0}$, $\Delta B_{0}$, $\Delta C_{0}$ and centrifugal distortion constants
calculated at the fc-CCSD(T)/ANO1 level. Best theoretical estimates (BE) for
the rotational and centrifugal distortion constants $X$ of H2CCCH+ are
estimated as
$X^{scaled}_{H_{2}CCCH^{+}}=\frac{X^{exp}_{H_{2}CCC}}{X^{calc}_{H_{2}CCC}}\times
X^{calc}_{H_{2}CCCH^{+}}$, i.e., using isoelectronic H2CCC as a calibrator.

|  H2CCC |  H2CCCH+ |  |
---|---|---|---|---
| Experiment | Calculated | Calculated | Scaled | Experiment
Parameter | Vrtilek et al. (1990) | This study | This study | This study, BE | This study
$A_{e}$ | $\cdots$ | 292302 | .539 | 285650 | .316 | $\cdots$ | $\cdots$
$B_{e}$ | $\cdots$ | 10578 | .491 | 9675 | .219 | $\cdots$ | $\cdots$
$C_{e}$ | $\cdots$ | 10209 | .024 | 9358 | .247 | $\cdots$ | $\cdots$
$\Delta A_{0}$ | $\cdots$ | 2357 | .134 | 2331 | .884 | $\cdots$ | $\cdots$
$\Delta B_{0}$ | $\cdots$ |  $-$3 | .807 | 5 | .043 | $\cdots$ | $\cdots$
$\Delta C_{0}$ | $\cdots$ | 10 | .869 | 20 | .431 | $\cdots$ | $\cdots$
$A_{0}$ | 288783 | .(34) | 289945 | .405 | 283318 | .432 | 282182 | .595 | 281856 | .(247)
$B_{0}$ | 10588 | .639(2) | 10582 | .298 | 9670 | .175 | 9675 | .970 | 9675 | .841(1)
$C_{0}$ | 10203 | .966(2) | 10198 | .155 | 9337 | .816 | 9343 | .136 | 9342 | .877(1)
$D_{J}\times 10^{3}$ | 4 | .248(2) | 3 | .722 | 2 | .684 | 3 | .063 | 2 | .942(4)
$D_{JK}$ | 0 | .5164(5) | 0 | .571 | 0 | .479 | 0 | .4328 | 0 | .4388(9)
$D_{K}$ | 23 |  .535a | 21 | .934 | 20 | .678 | [20 | .678] | 20 |  .678a
$d_{1}\times 10^{3}$ |  $-$0 | .153(2) |  $-$0 | .143 |  $-$0 | .096 |  $-$0 | .103 |  $-$0 | .121(3)
$d_{2}\times 10^{3}$ |  $-$0 | .070(1) |  $-$0. | 050 |  $-$0 | .037 |  $-$0 | .051 |  $-$0 |  .051a
555
$a$$a$footnotetext: Kept fixed in the analysis.

##  Appendix B Line parameters of H2CCCH+ and H2CCC

The line parameters of H2CCCH+ and H2CCC have been obtained by fitting a
Gaussian line profile to the observed data. The results are given in Tables 4
and 5, respectively. The observed lines of H2CCCH+ are shown in Fig. 2 and
those of H2CCC in Fig. 4.

Table 4: Line parameters of the observed transitions of H2CCCH+ in TMC-1.

Transition |  $\nu_{obs}$a |  $\int$ $T_{A}^{*}$ dv b |  $\Delta$vc | T${}_{A}^{*}$d |  $\sigma$e | Notes
---|---|---|---|---|---|---
| (MHz) | (mK $\times$ kms-1) | (kms-1) | (mK) | (mK) |
$2_{1,2}-1_{1,1}$ | 37702.627$\pm$0.010 | 0.38$\pm$0.13 | 0.63$\pm$0.20 | 0.56 | 0.13 |
$2_{0,2}-1_{0,1}$ | 38037.044$\pm$0.010 | 0.36$\pm$0.09 | 0.78$\pm$0.17 | 0.43 | 0.13 |
$2_{1,1}-1_{1,0}$ | 38368.581$\pm$0.010 | 0.57$\pm$0.07 | 0.77$\pm$0.11 | 0.70 | 0.13 | A
$4_{0,4}-3_{0,3}$ | 76071.062$\pm$0.020 |  |  |  | 3.40 | B
$4_{1,3}-3_{1,2}$ | 76735.888$\pm$0.020 | 2.28$\pm$0.07 | 0.40$\pm$0.10 | 5.50 | 1.80 | CD
$5_{1,5}-4_{1,4}$ | 94254.055$\pm$0.004 |  |  |  | 2.04 | B
$5_{0,5}-4_{0,4}$ | 95085.988$\pm$0.030 | 3.27$\pm$0.08 | 0.85$\pm$0.23 | 3.60 | 1.40 | BD
$5_{1,4}-4_{1,3}$ | 95918.765$\pm$0.020 | 1.04$\pm$0.20 | 0.40$\pm$0.10 | 3.30 | 0.67 |
$6_{1,6}-5_{1,5}$ | 113103.272$\pm$0.004 |  |  |  | 4.40 | B
$6_{0,6}-5_{0,5}$ | 114099.086$\pm$0.006 |  |  |  | 4.90 | B
666
$a$$a$footnotetext: Observed frequency assuming a vLSR of 5.83 km s-1.
$b$$b$footnotetext: Integrated line intensity in mK $\times$ km s-1.
$c$$c$footnotetext: Line width in km s-1. $d$$d$footnotetext: Antenna
temperature in mK. $e$$e$footnotetext: Root mean square noise of the data.
$A$$A$footnotetext: Blended with a weak feature at 38638.7 MHz.
$B$$B$footnotetext: For undetected lines their rest frequencies correspond to
the observed frequency in the laboratory (Table 1), or to the predicted one
from the molecular constants of Table 2. $C$$C$footnotetext: The feature
appears too strong. $D$$D$footnotetext: Marginal detection.

Figure 4: Observed lines of H2CCC in the 31–115 GHz domain towards TMC-1. The
abscissa corresponds to the velocity of the cloud with respect to the Local
Standard of Rest (vLSR). The ordinate is the antenna temperature corrected for
atmospheric and telescope losses in mK. The spectral resolution is 38 kHz
below 50 GHz and 48 kHz above. Derived line parameters are given in Table 5.
The red lines show the computed synthetic spectra for the lines of H2CCC (see
text). The lines in the left column correspond to para lines of H2CCC, while
those of the central and right columns correspond to transitions of the ortho
species. The adopted frequencies are those of the CDMS catalogue (Müller et
al. 2005) and are given in Table 5. Table 5: Line parameters of the observed
transitions of H2CCC in TMC-1.

Transition |  $\nu_{obs}$a |  $\int$ $T_{A}^{*}$ dv b | vLSRc |  $\Delta$vd | T${}_{A}^{*}$e |  $\sigma$f
---|---|---|---|---|---|---
| (MHz) | (mK $\times$ kms-1) | (kms-1) | (kms-1) | (mK) | (mK)
$2_{1,2}-1_{1,1}$ | 41198.335$\pm$0.002 | 82.7$\pm$0.1 | 5.70$\pm$0.01 | 0.65$\pm$0.01 | 120.4 | 0.1
$2_{0,2}-1_{0,1}$ | 41584.675$\pm$0.001 | 57.8$\pm$0.1 | 5.76$\pm$0.01 | 0.63$\pm$0.01 | 86.9 | 0.2
$2_{1,1}-1_{1,0}$ | 41967.671$\pm$0.002 | 84.6$\pm$0.1 | 5.75$\pm$0.01 | 0.61$\pm$0.01 | 130.8 | 0.2
$4_{1,4}-3_{1,3}$ | 82395.089$\pm$0.003 | 61.9$\pm$1.3 | 5.71$\pm$0.01 | 0.47$\pm$0.01 | 124.8 | 3.0
$4_{0,4}-3_{0,3}$ | 83165.345$\pm$0.003 | 31.9$\pm$0.4 | 5.77$\pm$0.02 | 0.54$\pm$0.01 | 55.5 | 0.9
$4_{1,3}-3_{1,2}$ | 83933.699$\pm$0.003 | 56.6$\pm$0.4 | 5.75$\pm$0.02 | 0.54$\pm$0.01 | 97.9 | 0.9
$5_{1,5}-4_{1,4}$ | 102992.379$\pm$0.004 | 25.4$\pm$1.2 | 5.73$\pm$0.02 | 0.47$\pm$0.03 | 51.1 | 3.1
$5_{0,5}-4_{0,4}$ | 103952.926$\pm$0.003 | 86.7$\pm$1.1 | 5.71$\pm$0.03 | 0.29$\pm$0.05 | 28.0 | 3.6
$5_{1,4}-4_{1,3}$ | 104915.583$\pm$0.003 | 19.0$\pm$1.1 | 5.77$\pm$0.03 | 0.48$\pm$0.03 | 37.4 | 2.9
777
$a$$a$footnotetext: Adopted rest frequencies from the CDMS entry for H2CCC.
$b$$b$footnotetext: Integrated line intensity in mK $\times$ km s-1.
$c$$c$footnotetext: Velocity of the line in km s-1. $d$$d$footnotetext: Line
width in km s-1. $e$$e$footnotetext: Antenna temperature in mK.
$f$$f$footnotetext: Root mean square noise of the data.


# Emergent Classical Spin Liquid Phases in an Ising Lattice via Size Effects

Pratik Brahma1, Sayeef Salahuddin1,2 1 Department of Electrical Engineering
and Computer Sciences, University of California, Berkeley
2 Materials Science Division, Lawrence Berkeley National Laboratory
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We show that a classical spin liquid phase can emerge from an ordered magnetic
state in the two-dimensional frustrated Shastry-Sutherland Ising lattice due
to lateral confinement. Two distinct classical spin liquid states are
stabilized (i) long-range spin-correlated dimers, and (ii) exponentially
decaying spin-correlated disordered states, depending on widths of W=3n, 3n+1
or W=3n+2, n being a positive integer. Stabilization of spin liquids in a
square-triangular lattice moves beyond the conventional geometric paradigm of
kagome, triangular or tetrahedral arrangements of antiferromagnetic ions,
where spin liquids have been discussed conventionally.

††preprint: APS/123-QED

Geometric frustration is an essential tool for stabilizing emergent quantum
phenomena in magnetic materials due to the absence of a unique state that
simultaneously satisfies all interactions among the magnetic ions. The quantum
spin liquid state is an example of a frustration-driven emergent phase
comprising highly correlated spin fluctuations down to zero temperature [1, 2,
3, 4]. This disordered state harbors many exotic quantum properties like
fractional charges, and long-range entanglement, which is relevant for high-
temperature superconductivity [5] and topological quantum computing [6].
Similarly, the classical counterpart defined by classical magnetic moments on
a geometrically frustrated lattice demonstrates finite thermodynamic entropy
per spin down to zero temperature without any quantum entanglement [7, 8, 9].
This is due to the presence of macroscopic ground state degeneracy, which is
considered a possible route to quantum spin liquids [10].

As yet, these exotic phases of matter have been stabilized only in a few
select geometric arrangements: (i) two-dimensional Ising triangular lattice
[11, 12], (ii) two-dimensional kagome lattice [13], (iii) three-dimensional
hyper-hyper kagome lattice [14], and (iv) three-dimensional pyrochlore spin
ice [15, 9]. In this context, geometric confinement may provide a new and
powerful way to induce novel phases. Indeed, in the recent years, confinement
has been utilized for enhanced ferroelectricity in unit cell thick HfO2 and
ZrO2[16, 17], ferromagnetic order in bilayer van der Waals crystal Cr2Ge2Te6
[18],high-temperature superconductivity in monolayer Bi2Sr2CaCu2O8+δ [19].
Here, we show that geometric confinement in a Sastry-Sutherland Ising lattice
can lead to two distinct Spin Liquid phases, that has not been previously
reported for a square-triangular lattice.

We study the classical 2D frustrated Shastry-Sutherland (SS) Ising lattice
(Fig.1a), which is topologically equivalent to the square-triangular
Archimedean lattice [20, 21]. This lattice model and its variations may
explain the presence of several fractional magnetization plateaus in quasi-
two-dimensional magnetic materials, such as the rare-earth tetraborides
($\text{RB}_{4}$, $\text{R}={\text{Tm},\text{Er},\text{Ho}}$ [22, 23, 24,
25]), where each magnetic atom in a layer is located on a 2D lattice that is
topologically equivalent to the SS lattice. Due to the large magnetic moments
of the ions and crystal field effects, these compounds are well described by
the classical Ising version of the SS Hamiltonian given by:

$\displaystyle\mathcal{H}\left(\mathbf{\sigma}\right)=J_{1}\sum_{\langle
i,j\rangle}\sigma_{i}\sigma_{j}+J_{2}\sum_{\langle\langle
i,j\rangle\rangle}\sigma_{i}\sigma_{j}+h\sum_{i}\sigma_{i}$ (1)

$\sigma_{i}\in\\{+1,-1\\}$, $\langle i,j\rangle$ are the nearest neighbor spin
interactions and $\langle\langle i,j\rangle\rangle$ are the diagonal and off-
diagonal row-wise alternating second nearest neighbor spin interactions shown
in Fig.1a. Frustration in the lattice arises from the competition between the
$J_{1}$ and $J_{2}$ antiferromagnetic interacting bonds giving rise to the 1/3
UUD (Up-Up-Down) fractional magnetization state inside a suitable range of
longitudinal magnetic field ($h_{z}$), where each triangular motif of the
lattice contains two up spins and one down spin [20]. Recent studies show the
presence of a narrow spin liquid phase in the bulk quantum SS model [26, 27]
and new quantum phases in the SS magnet SrCu2(BO3)2 under high field and
pressure [28] however, the role of size effects in stabilizing spin liquids in
the SS Ising lattice has yet to be explored.

Figure 1: Emergence of classical spin liquid phase due to size effects: (a)
All neighboring interactions in the Shastry-Sutherland (SS) Ising lattice
which is topologically equivalent to the square-triangular lattice
(Archimedean $3^{2}.4.3.4$ lattice) [20, 21]. (b) Introducing size effects
along the lateral direction in the SS lattice via open boundary conditions.
The ordered UUD spin configuration converts to a disordered classical spin
liquid characterized by the presence of macroscopically degenerate ground
states. The bottom graphs show the general behavior of $N_{G}$ (number of
degenerate ground states), and $C(r)$ (spin-spin correlation) for the ordered
spin state and the disordered spin liquid.  Figure 2: Transformation of
magnetization phase diagram due to size effects along lateral direction: (a)
Magnetization phase diagram in the thermodynamic limit of the Shastry-
Sutherland (SS) Ising lattice. The ground state spin configurations for all
the magnetization phases are shown beside the phase diagram. Only the Dimer
phase is macroscopically degenerate with zero spin correlation among the
dimers (shown in green boxes), while the other spin ground states are finitely
degenerate. (b) The phase diagrams of the laterally confined SS lattice via
open boundary conditions. Depending on the width of the lattice (grouped by
the modulus of 3), different regions of the 1/3 magnetization phase convert to
a classical spin liquid phase (shown as black hashed regions), identified by
its macroscopic degeneracy.

We calculate the ground state energy, ground state spin configuration, spin-
spin correlation, and magnetization phase diagram of the SS Ising lattice
using the total partition function [29]. Since brute force calculation of the
partition function is time-consuming and exponential in the lattice size, we
use the transfer matrix method, which has been used for various spin lattice
models [30, 31]. We assume the size of the lattice to be $W\times L$ and
define variables and calculate the Transfer Matrix from:

$\displaystyle T^{k}_{\left[\sigma_{i}\right],\left[\sigma_{i+1}\right]}$
$\displaystyle=$ $\displaystyle\exp\left(-\beta
E^{k}_{\left[\sigma_{i:i+1}\right]}\right),\quad k\in\\{0,1\\}$ (2)

$\beta=1/T$ is the inverse temperature of the system,
$E^{0}_{\left[\sigma_{i:i+1}\right]}\left(E^{1}_{\left[\sigma_{i:i+1}\right]}\right)$
is the lattice energy of the spins between columns $i$ and $i+1$ when the
diagonal bonds connect the top (bottom) left and bottom (top) right spin-
lattice sites and $T^{k}_{\left[\sigma_{i}\right],\left[\sigma_{i+1}\right]}$
is the transfer matrix whose rows and columns represent all the possible
combinations of spins at column $i$ $\left(\left[\sigma_{i}\right]\right)$ and
column $i+1$ $\left(\left[\sigma_{i+1}\right]\right)$ respectively. We
redefine the total partition function ($Z_{\beta}$) of the spin system as:

$\displaystyle Z_{\beta}$ $\displaystyle=$
$\displaystyle\sum_{\left[\sigma_{0:L-1}\right]}\exp\left(-\beta\mathcal{H}\left(\sigma\right)\right)=\text{Tr}\left[T^{\lfloor
L/2\rfloor}\right]$ (3)

With the partition function $\left(Z_{\beta}\right)$ of the lattice, we then
calculate the expectation of one and two spin operators:

$\displaystyle\langle\sigma_{2i}\sigma_{2i+2j}\rangle$ $\displaystyle=$
$\displaystyle\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\frac{1}{Z_{\beta}}\text{Tr}\biggl{[}T^{i}\tau^{z}T^{j}$
(4) $\displaystyle\tau^{z}T^{\lfloor L/2\rfloor-i-j}\biggr{]}$
$\displaystyle\langle\sigma_{2i}\rangle$ $\displaystyle=$
$\displaystyle\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\frac{1}{Z_{\beta}}\text{Tr}\biggl{[}T^{i}\tau^{z}T^{\lfloor
L/2\rfloor-i}\biggr{]}$ (5)

$\tau^{z}$ is the Pauli-Z matrix and $T=T^{0}T^{1}$. These spin expectation
values can then be used to calculate the ground state energy and the spin
correlation. A more detailed description is provided in the Supplementary
Information

As the transfer matrix for our studied lattice is non-diagonalizable, we
resort to the numerical estimation of the partition function and the required
observables. We choose a large enough length $(L\geq 20$) and low enough
temperature $(\beta\geq 15$) so that the calculated observables converge to
their thermodynamic limits (Fig. S5). To study size effects, we introduce
periodic boundary conditions along the transverse direction (length) and open
boundary conditions along the lateral direction (width) of the lattice.
Subsequently, we use the numerical estimation of the partition function to
obtain phase diagrams by varying the interaction bond ratio, $J_{2}/J_{1}$
from 0 to 3, and the applied magnetic field, $h_{z}$ from -8 to 8 for widths
$W$ varying from 4 to 12. The above parameter ranges are chosen to ensure the
presence of $1/3$ magnetization plateau, according to studies on the SS
lattice [32, 20]. The ground state spin configuration for every magnetization
phase is deduced from the spin-spin correlation to all its nearest neighbors
and the average magnetization spin site calculations. These ground states are
also verified by comparing their energies to the ground state energies
obtained by Monte-Carlo sampling the Boltzmann distribution of the Shastry-
Sutherland Ising lattice [33, 34]. To further characterize the ground states
of the laterally confined SS Ising lattice, we calculate the spatial variation
of spin-spin correlation along the length of the lattice and the number of
degenerate ground states (as shown in Fig.3e and Fig.3f) using the following
equations:

$\displaystyle F$ $\displaystyle=$ $\displaystyle-
k_{b}T\ln\left(Z_{\beta}\right)$ (6) $\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{E_{gs}-F}{T}=k_{b}\ln(N_{G})$ (7)

$k_{b}$ is the Boltzmann constant, $F$ is the free energy, $S$ is the entropy,
$E_{gs}$ is the ground state energy and $N_{G}$ is the number of degenerate
ground states. For calculating $N_{G}$, we choose low enough temperatures such
that the average energy of the system converges to the ground state energy
($\langle E\rangle\approx E_{gs}$) and every spin configuration other than the
ground state has negligible Boltzmann probabilities. Thus, the entropy $S$
takes the form of eq.(7). It should be noted that macroscopically degenerate
systems (the number of degenerate ground states is a macroscopic property)
have non-trivial entropy even in the zero temperature limit [4].

Our main results predict the emergence of a classical spin liquid phase from
an ordered magnetic state when the Shastry-Sutherland (SS) Ising lattice is
confined along the lateral direction. Fig.2a shows the magnetization phase
diagram of the SS lattice in the thermodynamic limit. The ferromagnetic (FM),
antiferromagnetic (AFM), and 1/3 fractional magnetic phases have a finite
number of degenerate ground states while the dimer magnetic phase is
macroscopically degenerate at the zero temperature limit, due to the
uncorrelated antiferromagnetic dimers present along the diagonal bonds of the
SS lattice shown in Fig.2a. But our study focuses on the $1/3$ fractional
phase containing 6 degenerate ground states where each triangle in the lattice
contains two up spins and one down spin (UUD spin configuration) [20]. With
size effects introduced by confining the lattice along the lateral direction,
the ordered 1/3 fractional magnetization phase gets fragmented into slightly
different valued magnetic phases depending upon the width of the lattice. In
contrast, the other magnetization phases (FM, AFM, and Dimer) retain their
shapes and magnetization values as shown in Fig.2. We categorize the phase
diagrams of the laterally confined lattice according to their widths: $W=3n$,
$W=3n+1$, and $W=3n+2$, which also determines the region of the 1/3 fractional
phase that transforms into the classical spin liquid phase, identified by its
finite temperature-entropy product in the zero temperature limit.

Figure 3: Characterisation of classical spin liquid phase: (a) The ordered UUD
(1/3) magnetization spin state in the thermodynamic limit. (b) The number of
degenerate ground states for UUD is independent of the size of the system.
Size effect introduces two types of spin liquid depending on the width of the
lattice: (c) classical spin liquid for widths $3n$ and $3n+1$ where the
correlated dimers in pink boxes are sandwiched between antiferromagnetic rows
and (d) disordered classical spin liquid for width $3n+2$ where each spin is
fractionally spin-spin correlated to all its neighboring spins. The color bar
represents the average magnetization for each spin site. (e), (f) Both spin
states are identified as classical spin liquids via the exponential increase
of ground states with the length of the system. (g) Long-range spin-spin
correlation is present along a row of the lattice in the correlated dimer spin
liquid state. (h) Exponential spin-spin correlation decay in the disordered
spin liquid state.

Below a critical width, in (i) $W=3n$, the region of the 1/3 fractional phase
along the borders of the FM and within a longitudinal field ($h_{z}$) range
transforms into the $1/3+1/3n$ classical spin liquid phase, in (ii) $W=3n+1$,
the region of the 1/3 fractional phase bordering the dimer phase where
$J_{2}\geq 2J_{1}$ transforms into the $1/(3+1/n)$ classical spin liquid phase
and finally, in (iii) $W=3n+2$, the thin region of the 1/3 fractional phase
along the FM and the dimer phase transforms into the $(1+2/n)/(3+2/n)$ and the
$1/(3+2/n)$ classical spin liquid phase respectively. Fig. 2b shows all the
above emergent classical spin liquid phases as the black hashed region.

The ground state spin configurations of the aforementioned classical spin
liquid phases can be grouped into two categories, as shown in Fig.3. Both
$W=3n$ and $W=3n+1$ have the same ground state spin liquid phase described by
rows of correlated dimers sandwiched between antiferromagnetic rows
(Fig.3(c)). Due to the exponential combinations possible among the correlated
dimers, the number of degenerate ground states is a macroscopic property
similar to the dimer phase (Fig.3(e)) [20]. Lattices with small widths
($W=4,6$) contain uncorrelated dimers whereas for larger widths, there exist
persistent long-range correlated dimers (Fig.3(g)). For $W=3n+2$, the
classical spin liquid phase is entirely disordered, characterized by the
average fractional magnetization per site and the fractional spin-spin
correlation among all its nearest and second-nearest neighbors (Fig.3(d)).
Similar to the previous spin liquid phase, the number of degenerate ground
states is a macroscopic property (Fig.3(f)), but in contrast to the long-range
spin-spin correlation present in the previous case, the spin-spin correlation
decays exponentially along the length of the lattice for this particular spin
liquid state (Fig.3(h)). The detailed magnetization phase diagrams and the
ground states of all the magnetization plateaus for widths varying from $4-12$
are given in the supplemental Fig.S1 - S4. The emergence of these exotic
phases is explained by the lateral confinement or size effects introducing
additional frustration in the bulk of the lattice where the spin arrangement
and interaction energy at the edges destroy the long-range spin order in the
bulk and reduces the total energy of the system [35]. Thus, these results
demonstrate the key role of size effects in stabilizing spin liquids.

The dimers in the spin liquid states of lattice widths $3n$ and $3n+1$ are
always sandwiched between antiferromagnetic rows of opposite spin
arrangements. Thus, a lattice structure must have width $3d+1$ to contain d
rows of dimers interleaved with d + 1 antiferromagnetic rows. This spin
arrangement also leads to AFM spin configuration at the edges, which is
favorable for low hz values, stabilizing the dimer spin liquid for $W=3n+1$.
In contrast, for the higher hz values, the spins along both edges favor the FM
configuration to lower the ground state energy compared to the FM and AFM spin
edge configuration in the bulk UUD configuration. Adding the two FM spin rows
at the edges to the previous $d$ rowed dimer structure gets us a total of
$3d+3=3n$, thus leading to the spin liquid stabilization for $W=3n$ at higher
hz values. But for lattices with width, $3n+2$ can only support the $d$-dimer
structure with an FM and an AFM spin configuration at the edges, leading to an
increase in the energy for lower and higher hz values. This leads to the
formation of a new spin liquid with no dimers and exponential decaying spin
correlation, with both the edges having FM spin configurations for higher hz
values and AFM spin configurations for the lower hz values (Fig. S4).

In summary, a classical spin liquid state emerges from the ordered UUD
magnetization state in a laterally confined classical SS Ising lattice. In
particular, two classical spin liquid states are observed depending on the
width of the lattice: (i) correlated dimers sandwiched between
antiferromagnetic rows (widths $3n$ and $3n+1$), and (ii) spin liquid state
with fractional correlation to all its neighboring interacting spins (width
$3n+2$). This novel physical phenomenon originates from spin arrangement at
the edges which destroy the long-range bulk order to reduce the total energy
of the lattice; therefore, size effects and boundary conditions underlie spin
liquid stabilization in laterally confined frustrated 2D lattices. This
approach can pave the way for stabilizing new spin liquid systems through
lateral confinement of quasi-two-dimensional magnetic materials, e.g. rare-
earth tetraborides in which the ionic arrangement is topologically equivalent
to the SS lattice [22, 23, 24, 25].

This work was supported by ASCENT, one of six centers in JUMP, a Semiconductor
Research Corporation (SRC) program sponsored by DARPA.

## References

  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Reports on Progress in Physics 80, 016502 (2016).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Anderson [1987] P. W. Anderson, The resonating valence bond state in La2CuO4 and superconductivity, Science 235, 1196 (1987).
  * Ioffe  _et al._ [2002] L. B. Ioffe, M. V. Feigel’man, A. Ioselevich, D. Ivanov, M. Troyer, and G. Blatter, Topologically protected quantum bits using josephson junction arrays, Nature 415, 503 (2002).
  * Bramwell and Gingras [2001] S. T. Bramwell and M. J. P. Gingras, Spin ice state in frustrated magnetic pyrochlore materials, Science 294, 1495 (2001).
  * Gingras [2011] M. J. P. Gingras, Spin ice, in  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 293–329.
  * Ramirez  _et al._ [1999] A. P. Ramirez, A. Hayashi, R. J. Cava, R. Siddharthan, and B. S. Shastry, Zero-point entropy in ‘spin ice’, Nature 399, 333 (1999).
  * Niggemann  _et al._ [2019] N. Niggemann, M. Hering, and J. Reuther, Classical spiral spin liquids as a possible route to quantum spin liquids, Journal of Physics: Condensed Matter 32, 024001 (2019).
  * Arh  _et al._ [2022] T. Arh, B. Sana, M. Pregelj, P. Khuntia, Z. Jagličić, M. D. Le, P. K. Biswas, P. Manuel, L. Mangin-Thro, A. Ozarowski, and A. Zorko, The ising triangular-lattice antiferromagnet neodymium heptatantalate as a quantum spin liquid candidate, Nature Materials 21, 416 (2022).
  * Li  _et al._ [2015] Y. Li, G. Chen, W. Tong, L. Pi, J. Liu, Z. Yang, X. Wang, and Q. Zhang, Rare-earth triangular lattice spin liquid: A single-crystal study of ${\mathrm{ybmggao}}_{4}$, Phys. Rev. Lett. 115, 167203 (2015).
  * Kermarrec  _et al._ [2021] E. Kermarrec, R. Kumar, G. Bernard, R. Hénaff, P. Mendels, F. Bert, P. L. Paulose, B. K. Hazra, and B. Koteswararao, Classical spin liquid state in the s=5/2 heisenberg kagome antiferromagnet $\text{Li}_{9}\text{Fe}_{3}(\text{P}_{2}\text{O}_{7})_{3}(\text{PO}_{4})_{2}$, Physical review letters 127, 157202 (2021).
  * Chillal  _et al._ [2020] S. Chillal, Y. Iqbal, H. O. Jeschke, J. A. Rodriguez-Rivera, R. Bewley, P. Manuel, D. Khalyavin, P. Steffens, R. Thomale, A. T. M. N. Islam, J. Reuther, and B. Lake, Evidence for a three-dimensional quantum spin liquid in PbCuTe2O6, Nature Communications 11, 2348 (2020).
  * Carrasquilla  _et al._ [2015] J. Carrasquilla, Z. Hao, and R. G. Melko, A two-dimensional spin liquid in quantum kagome ice, Nature Communications 6, 7421 (2015).
  * Cheema  _et al._ [2020] S. S. Cheema, D. Kwon, N. Shanker, R. dos Reis, S. L. Hsu, J. Xiao, H. Zhang, R. Wagner, A. Datar, M. R. McCarter, C. R. Serrao, A. K. Yadav, G. Karbasian, C. H. Hsu, A. J. Tan, L. C. Wang, V. Thakare, X. Zhang, A. Mehta, E. Karapetrova, R. V. Chopdekar, P. Shafer, E. Arenholz, C. Hu, R. Proksch, R. Ramesh, J. Ciston, and S. Salahuddin, Enhanced ferroelectricity in ultrathin films grown directly on silicon, Nature 580, 478 (2020).
  * Cheema  _et al._ [2022] S. S. Cheema, N. Shanker, S.-L. Hsu, Y. Rho, C.-H. Hsu, V. A. Stoica, Z. Zhang, J. W. Freeland, P. Shafer, C. P. Grigoropoulos, J. Ciston, and S. Salahuddin, Emergent ferroelectricity in subnanometer binary oxide films on silicon, Science 376, 648 (2022).
  * Gong  _et al._ [2017] C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, and X. Zhang, Discovery of intrinsic ferromagnetism in two-dimensional van der waals crystals, Nature 546, 265 (2017).
  * Yu  _et al._ [2019] Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang, High-temperature superconductivity in monolayer Bi2Sr2CaCu2O8+δ, Nature 575, 156 (2019).
  * Dublenych [2012] Y. I. Dublenych, Ground states of the ising model on the shastry-sutherland lattice and the origin of the fractional magnetization plateaus in rare-earth-metal tetraborides, Phys. Rev. Lett. 109, 167202 (2012).
  * Their and St-Pierre [2017] T. Their and L. St-Pierre, Stiffness and strength of a semi-regular lattice, Rakenteiden Mekaniikka 50, 137 (2017).
  * Yoshii  _et al._ [2008] S. Yoshii, T. Yamamoto, M. Hagiwara, S. Michimura, A. Shigekawa, F. Iga, T. Takabatake, and K. Kindo, Multistep magnetization plateaus in the shastry-sutherland system tbb4, Phys. Rev. Lett. 101, 087202 (2008).
  * Kim  _et al._ [2009] J. Y. Kim, B. K. Cho, and S. H. Han, Anisotropic magnetic phase diagrams of HoB4 single crystal, Journal of Applied Physics 105, 07E116 (2009).
  * Trinh  _et al._ [2018] J. Trinh, S. Mitra, C. Panagopoulos, T. Kong, P. C. Canfield, and A. P. Ramirez, Degeneracy of the $1/8$ plateau and antiferromagnetic phases in the shastry-sutherland magnet ${\mathrm{tmb}}_{4}$, Phys. Rev. Lett. 121, 167203 (2018).
  * Siemensmeyer  _et al._ [2008] K. Siemensmeyer, E. Wulf, H.-J. Mikeska, K. Flachbart, S. Gabáni, S. Mat’aš, P. Priputen, A. Efdokimova, and N. Shitsevalova, Fractional magnetization plateaus and magnetic order in the shastry-sutherland magnet ${\mathrm{tmb}}_{4}$, Phys. Rev. Lett. 101, 177201 (2008).
  * Yang  _et al._ [2022] J. Yang, A. W. Sandvik, and L. Wang, Quantum criticality and spin liquid phase in the shastry-sutherland model, Phys. Rev. B 105, L060409 (2022).
  * Keleş and Zhao [2022] A. Keleş and E. Zhao, Rise and fall of plaquette order in the shastry-sutherland magnet revealed by pseudofermion functional renormalization group, Phys. Rev. B 105, L041115 (2022).
  * Shi  _et al._ [2022] Z. Shi, S. Dissanayake, P. Corboz, W. Steinhardt, D. Graf, D. M. Silevitch, H. A. Dabkowska, T. F. Rosenbaum, F. Mila, and S. Haravifard, Discovery of quantum phases in the shastry-sutherland compound srcu2(bo3)2 under extreme conditions of field and pressure, Nature Communications 13, 10.1038/s41467-022-30036-w (2022).
  * Schmitz [2017] K. S. Schmitz, Chapter 12 - classical statistical mechanics, in  _Physical Chemistry_, edited by K. S. Schmitz (Elsevier, Boston, 2017) pp. 559–632.
  * Kramers and Wannier [1941] H. A. Kramers and G. H. Wannier, Statistics of the two-dimensional ferromagnet. part i, Phys. Rev. 60, 252 (1941).
  * Onsager [1944] L. Onsager, Crystal statistics. i. a two-dimensional model with an order-disorder transition, Phys. Rev. 65, 117 (1944).
  * Kairys  _et al._ [2020] P. Kairys, A. D. King, I. Ozfidan, K. Boothby, J. Raymond, A. Banerjee, and T. S. Humble, Simulating the shastry-sutherland ising model using quantum annealing, PRX Quantum 1, 020320 (2020).
  * Kotze [2008] J. Kotze, Introduction to monte carlo methods for an ising model of a ferromagnet (2008).
  * Huang  _et al._ [2012] W. C. Huang, L. Huo, G. Tian, H. R. Qian, X. S. Gao, M. H. Qin, and J.-M. Liu, Multi-step magnetization of the ising model on a shastry–sutherland lattice: a monte carlo simulation, Journal of Physics: Condensed Matter 24, 386003 (2012).
  * Marić  _et al._ [2020] V. Marić, S. M. Giampaolo, D. Kuić, and F. Franchini, The frustration of being odd: how boundary conditions can destroy local order, New Journal of Physics 22, 083024 (2020).

Supplemental Materials: Emergent Classical Spin Liquid Phases in an Ising
Lattice via Size Effects

##  I Calculation of Partition Function

Using the total partition function, we calculate the ground state energy,
ground state spin configuration, spin-spin correlation, and magnetization
phase diagram of the SS Ising lattice. We assume the size of the lattice to be
$W\times L$ and define variables (i) $\sigma_{i}^{j}$, the spin site at column
$i$ and row $j$, (ii) $\left[\sigma_{a:b}\right]$, the set of all spin sites
between column $a$ and column $b$. Using the above definitions, we introduce
the following variables:

$\displaystyle E^{H}_{\left[\sigma_{i:i+1}\right]}$ $\displaystyle=$
$\displaystyle J_{1}\sum_{j=0}^{W-1}\sigma_{i}^{j}\sigma_{i+1}^{j}$ (S1)
$\displaystyle E^{V}_{\left[\sigma_{i:i+1}\right]}$ $\displaystyle=$
$\displaystyle
J_{1}\sum_{j=0}^{W-1}\sigma_{i}^{j}\sigma_{i}^{j+1}+\sigma_{i+1}^{j}\sigma_{i+1}^{j+1}$
(S2) $\displaystyle E^{k}_{\left[\sigma_{i:i+1}\right]}$ $\displaystyle=$
$\displaystyle
E^{H}_{\left[\sigma_{i:i+1}\right]}+\frac{1}{2}E^{V}_{\left[\sigma_{i:i+1}\right]}+J_{2}\sum_{j=0}^{\lfloor
W/2\rfloor}\sigma_{i}^{2j}\sigma_{i+1}^{2j+(-1)^{k}}\quad,k\in\\{0,1\\}$ (S3)
$\displaystyle T^{k}_{\left[\sigma_{i}\right],\left[\sigma_{i+1}\right]}$
$\displaystyle=$ $\displaystyle\exp\left(-\beta
E^{k}_{\left[\sigma_{i:i+1}\right]}\right)$ (S4)

$\beta=1/T$ is the inverse temperature of the system,
$E^{0}_{\left[\sigma_{i:i+1}\right]}\left(E^{1}_{\left[\sigma_{i:i+1}\right]}\right)$
is the lattice energy of the spins between columns $i$ and $i+1$ when the
diagonal bonds connect the top (bottom) left and bottom (top) right spin-
lattice sites and $T^{k}_{\left[\sigma_{i}\right],\left[\sigma_{i+1}\right]}$
is the transfer matrix whose rows and columns represent all the possible
combinations of spins at column $i$ $\left(\left[\sigma_{i}\right]\right)$ and
column $i+1$ $\left(\left[\sigma_{i+1}\right]\right)$ respectively. We
redefine the total partition function ($Z_{\beta}$) of the spin system as:

$\displaystyle Z_{\beta}$ $\displaystyle=$
$\displaystyle\sum_{\left[\sigma_{0:L-1}\right]}\exp\left(-\beta\mathcal{H}\left(\sigma\right)\right)$
(S5) $\displaystyle=$
$\displaystyle\sum_{\left[\sigma_{0:L-1}\right]}\exp\biggl{(}-\beta\sum_{i=0}^{\lfloor
L/2\rfloor}E^{0}_{\left[\sigma_{2i:2i+1}\right]}+E^{1}_{\left[\sigma_{2i+1:2i+2}\right]}\biggr{)}$
$\displaystyle Z_{\beta}$ $\displaystyle=$
$\displaystyle\sum_{\left[\sigma_{0:L-1}\right]}\prod_{i=0}^{\lfloor
L/2\rfloor}T^{0}_{\left[\sigma_{2i}\right],\left[\sigma_{2i+1}\right]}T^{1}_{\left[\sigma_{2i+1}\right],\left[\sigma_{2i+2}\right]}$
(S6) $\displaystyle Z_{\beta}$ $\displaystyle=$
$\displaystyle\text{Tr}\left[\left(T^{0}T^{1}\right)^{\lfloor
L/2\rfloor}\right]=\text{Tr}\left[T^{\lfloor L/2\rfloor}\right]$ (S7)

With the partition function $\left(Z_{\beta}\right)$ of the lattice, we then
calculate the expectation of one and two spin operators:

$\displaystyle\langle\sigma_{2i}\sigma_{2i+2j}\rangle$ $\displaystyle=$
$\displaystyle\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\frac{1}{Z_{\beta}}\text{Tr}\biggl{[}T^{i}\tau^{z}T^{j}\tau^{z}T^{\lfloor
L/2\rfloor-i-j}\biggr{]}$ (S8) $\displaystyle\langle\sigma_{2i}\rangle$
$\displaystyle=$
$\displaystyle\lim_{\beta\rightarrow\infty}\lim_{L\rightarrow\infty}\frac{1}{Z_{\beta}}\text{Tr}\biggl{[}T^{i}\tau^{z}T^{\lfloor
L/2\rfloor-i}\biggr{]}$ (S9)

$\tau^{z}$ is the Pauli-Z matrix. These spin expectation values can then be
used to calculate the ground state energy and the spin correlation:

$\displaystyle C(2j)$ $\displaystyle=$
$\displaystyle\langle\sigma_{2i}\sigma_{2i+2j}\rangle-\langle\sigma_{2i}\rangle\langle\sigma_{2i+2j}\rangle$
(S10) $\displaystyle E_{gs}$ $\displaystyle=$ $\displaystyle
J_{1}\sum_{\langle
i,j\rangle}\langle\sigma_{i}\sigma_{j}\rangle+J_{2}\sum_{\langle\langle
i,j\rangle\rangle}\langle\sigma_{i}\sigma_{j}\rangle+h\sum_{i}\langle\sigma_{i}\rangle$
(S11)

$E_{gs}$ is the ground state energy and $C(2j)$ is the spin-spin correlation
between two spins at lattice sites $2j$ distance apart.

##  II Phase Diagrams and Ground State Lattices

Figure S1: Magnetization phase diagrams for widths 4 to 12: The phase diagrams
have been grouped according to modulus 3 of the width of the lattice. The
following phases were not shown in the main Fig.2 due to brevity : (i) M = 1/9
for W = 9, (ii) M = 1/7 for W = 7, (iii) M = 1/11, M = 3/11 for W = 11, (iv) M
= 1/5 for W = 5.  Figure S2: Ground state solutions of all magnetization
phases for width $3n$:  The description of the spin configuration and the
number of degenerate ground states are given below each lattice diagram. The
green boxes represent the antiferromagnetic dimers, the white circle
represents +1 and the black circle represents -1. (a) is the original dimer
phase in the thermodynamic limit with antiferromagnetic edges, (b)-(e) are the
ordered spin lattice states and (f)-(h) are the classical spin liquids. The
correlation among the dimers is written beside each dimer row. Figure S3:
Ground state solutions of all magnetization phases for width $3n+1$: The
description of the spin configuration and the number of degenerate ground
states are given below each lattice diagram. The green boxes represent the
antiferromagnetic dimers, the white circle represents +1 and the black circle
represents -1. (a) is the original dimer phase in the thermodynamic limit with
antiferromagnetic edges, (b)-(e) are the ordered spin lattice states and
(f)-(h) are the classical spin liquids. The correlation among the dimers is
written beside each dimer row. Figure S4: Ground state solutions of all
magnetization phases for width $3n+2$: The description of the spin
configuration and the number of degenerate ground states are given below each
lattice diagram. The green boxes represent the antiferromagnetic dimers, the
white circle represents +1 and the black circle represents -1. (a) is the
original dimer phase in the thermodynamic limit with antiferromagnetic edges,
(b)-(f) are the ordered spin lattice states and (g)-(h) are the disordered
classical spin liquids due to the fractional spin-spin correlations among all
the interacting spins. The color bar represents the average magnetization for
each spin site. Figure S5: Convergence to the thermodynamic value: (a) Energy
per spin converges to the thermodynamic value with an increase in inverse
temperature ($\beta$). Therefore, we choose $\beta>=15$ for all computations.
(b) Energy per spin remains almost constant with an increase in the length of
the lattice. We choose $L>=20$ to ensure thermodynamic convergence. (c) The
magnetization phase converges to the thermodynamic limit with an increase in
inverse temperature ($\beta$). Therefore, we choose $\beta>=15$ for all
computations. (d) The magnetization phase remains almost constant with an
increase in the length of the lattice. We choose $L>=20$ to ensure
thermodynamic convergence.