# High-Harmonic Generation with a twist: all-optical characterization of
magic-angle twisted bilayer graphene.

Eduardo B. Molinero<EMAIL_ADDRESS>_Instituto de Ciencia de Materiales
de Madrid (ICMM), Consejo Superior de Investigaciones Científicas (CSIC), Sor
Juana Inés de la Cruz 3, 28049 Madrid, Spain_ Anushree Datta  _Instituto de
Ciencia de Materiales de Madrid (ICMM), Consejo Superior de Investigaciones
Científicas (CSIC), Sor Juana Inés de la Cruz 3, 28049 Madrid, Spain_ M. J.
Calderón  _Instituto de Ciencia de Materiales de Madrid (ICMM), Consejo
Superior de Investigaciones Científicas (CSIC), Sor Juana Inés de la Cruz 3,
28049 Madrid, Spain_ E. Bascones  _Instituto de Ciencia de Materiales de
Madrid (ICMM), Consejo Superior de Investigaciones Científicas (CSIC), Sor
Juana Inés de la Cruz 3, 28049 Madrid, Spain_ Rui E. F. Silva
<EMAIL_ADDRESS><EMAIL_ADDRESS>_Instituto de Ciencia de Materiales
de Madrid (ICMM), Consejo Superior de Investigaciones Científicas (CSIC), Sor
Juana Inés de la Cruz 3, 28049 Madrid, Spain_ _Max Born Institute, Max-Born-
Strasse 2A, D-12489 Berlin, Germany_

###### Abstract

If we stack up two layers of graphene while changing their respective
orientation by some twisting angle, we end up with a system that has striking
differences when compared to single-layer graphene. For a very specific value
of this twist angle, known as _magic angle_ , twisted bilayer graphene
displays a unique phase diagram that cannot be found in other systems.
Recently, high harmonic generation spectroscopy has been successfully applied
to elucidate the electronic properties of quantum materials. The purpose of
the present work is to exploit the nonlinear optical response of magic-angle
twisted bilayer graphene to unveil its electronic properties. We show that the
band structure of magic-angle twisted bilayer graphene is imprinted onto its
high-harmonic spectrum. Specifically, we observe a drastic decrease of
harmonic signal as we approach the magic angle. Our results show that high
harmonic generation can be used as a spectroscopy tool for measuring the twist
angle and also the electronic properties of twisted bilayer graphene, paving
the way for an all-optical characterization of moiré materials.

## Introduction

The discovery of magic-angle twisted bilayer graphene (MATBG) (Cao  _et al._ ,
2018a, b) has sparked a huge increase of attraction into moiré quantum
materials (Andrei  _et al._ , 2021) in the recent years. One of the key
reasons for that increase is that MATBG shows a unique plethora of exotic
phenomena (Andrei and MacDonald, 2020), which includes correlated insulators
(Cao  _et al._ , 2018a; Lu  _et al._ , 2019), unconventional superconductivity
(Cao  _et al._ , 2018b; Yankowitz  _et al._ , 2019; Lu  _et al._ , 2019),
interacting topological phases (Nuckolls  _et al._ , 2020; Choi  _et al._ ,
2021; Xie  _et al._ , 2021a), ferromagnetism (Sharpe  _et al._ , 2019) and
also strange metal behavior (Cao  _et al._ , 2020; Jaoui  _et al._ , 2022).
Such phase diagram has some resemblance with the phenomenology present in
cuprates (Keimer  _et al._ , 2015) and iron pnictides (Fernandes  _et al._ ,
2022) superconductors. For instance, both systems share the presence of
superconducting domes surrounded by correlated insulating states as doping is
varied. Moiré materials have opened completely new avenues in research thanks
to the unprecedented in-situ tunability of their electronic properties through
electronic gates, strain, magnetic field, pressure or substrate alignement
among others. This has for example allowed observing many different correlated
insulating states in a single MATBG sample, including integer and fractional
Chern insulators as well as trivial insulators with or without charge density
modulations (Xie  _et al._ , 2021a).

At the heart of the correlated phases of MATBG it is the existence of two flat
bands (Bistritzer and MacDonald, 2011; Lisi  _et al._ , 2021; Lopes dos Santos
_et al._ , 2007) at the so-called magic angle $\theta_{\text{mag}}\approx
1.1^{\circ}$. As the bandwidth is reduced, the role played by the interactions
increases (Choi  _et al._ , 2019; Datta  _et al._ , ) and rich and complex
physics arises. The nature of these flat bands is far from trivial, as they
are formed by the hybridization of the Dirac cones of each layer of graphene
(Cao  _et al._ , 2018a; San-Jose  _et al._ , 2012). A proper characterization
of the low energy bandstructure of MATBG is key to unveil the origin of the
correlated states. On the other hand, studying the electronic properties of
these flat bands requires of new spectroscopic tools as the large moiré unit
cell and the two dimensionality of MATBG complicate the applicability of many
experimental techniques of common use in other strongly correlated materials.

Figure 1: Schematic figure of the moiré superlattice alongside with the laser
directions. Note that the $\Gamma-M$ direction connects $AA$ to $AA$ zones
while $\Gamma-K$ direction connects $AA$ to $AB/BA$ zones. We note that the
moiré Brillouin zone is not shown at scale.

In parallel, the study of the interaction between ultrashort laser pulses and
condensed matter systems has also experienced a huge growth in relevance
(Kruchinin  _et al._ , 2018). Such short pulses can drive the system’s
electrons into a non-equilibrium excited state within times below their own
cycle. As a consequence of this excitation, the electrons emit coherent
radiation in multiples (up to the hundreds) of the laser frequency, providing
us a tool to unveil their dynamics. Even though this high-harmonic generation
(HHG) had its origins in atomic and molecular physics (Krausz and Ivanov,
2009), it has been recently, and successfully, applied to condensed matter
systems as a spectroscopy tool (Ghimire and Reis, 2018; Goulielmakis and
Brabec, 2022). For instance, it has allowed for the observation of ultrafast
electron-hole dynamics (Vampa  _et al._ , 2015a), insulator-to-metal
transitions (Silva  _et al._ , 2018; Bionta  _et al._ , 2021), topological
transitions (Silva  _et al._ , 2019a; Chacón  _et al._ , 2020), and light-
driven band structure (Uzan-Narovlansky  _et al._ , 2022).

Figure 2:  Band structure of TBG for 3 distinct twist angles alongside its
Fermi-Dirac distribution. a corresponds to $\theta=1.05^{\circ}$, b to
$\theta=1.12^{\circ}$ and c to $\theta=1.22^{\circ}$, showing the appearance
of flat bands at the so-called magic angle $\theta=1.12^{\circ}$. Red color
indicates the band structure associated with the $K$ valley, while blue lines
depicts its $K^{\prime}$ counterpart (see Methods for a discussion on the
model).

Why HHG can be effectively used as a spectroscopy tool, can be understood
thanks to its underlying physical mechanisms. The HHG spectrum in solids comes
from two distinct contributions, namely intraband and interband currents.
Intraband harmonics are generated due to the accelerated motion of the
electron (hole) in the conduction (valence) band after a multiphoton/tunneling
excitation of the electron from the valence band. Afterwards, the
recombination of such pair gives rise to the interband harmonics. While
intraband contributions are generally associated with low order harmonics
(Jürgens  _et al._ , 2020), higher-order ones are caused by the latter
contributions (Vampa  _et al._ , 2014). In the presence of a strong low
frequency laser field, the charge carriers are drifted through the whole
Brillouin zone, exploring sections, and hence physics, that are not commonly
seen by other optical tests, such as second harmonic generation (Boyd, 2008).

In the present work, we aim to bridge the gap between these two distinct
fields by studying HHG in magic-angle twisted bilayer graphene. Remarkably, in
spite of extensive work in twisted bilayer graphene (Bernevig  _et al._ ,
2021a; Song  _et al._ , 2021; Bernevig  _et al._ , 2021b; Lian  _et al._ ,
2021; Bernevig  _et al._ , 2021c; Xie  _et al._ , 2021b), its ultrafast
dynamics have been hardly explored, with very few exceptions (Ikeda, 2020;
Mrudul and Dixit, 2021; Du  _et al._ , 2021; Yang  _et al._ , 2020). However
and to the knowledge of the authors, no previous work has dealt with highly
nonlinear optical response in MATBG. Here we show that the emitted spectrum
reveals the formation of the flat bands at the magic angle.

## results

As it has been discussed earlier, the HHG of solids is extremely dependent on
the electronic band structure of the system (Tancogne-Dejean  _et al._ , 2017;
Silva  _et al._ , 2019b), up to the point that one can reconstruct the band
structure from its high-harmonic spectrum (Vampa  _et al._ , 2015b).
Therefore, the harmonic spectra of twisted bilayer graphene must inevitably
reflect the existence (or absence) of the flat bands.

Figure 3:  a and c, High harmonic spectrum for a set of twist angles. b and d,
Yield for the third (H3), fifth (H5) and seventh (H7) harmonic. a and c
results are obtained for a laser in the $\Gamma-M$ direction, while c and d
are obtained for the same laser but in the $\Gamma-K$ direction. All
calculations have been done for a temperature of $80$ K. One can see the clear
depression of the harmonic signal, for all subfigures, near the magic angle.

The purpose of this work is to show that one can characterize the magic angle
in an all-optical way, offering an alternative to other characterization
techniques, such as scanning tunneling microscopy (Lisi  _et al._ , 2021). We
show that when the twist angle of the TBG gets near $\theta_{\text{mag}}$, the
intensity of the high harmonics (from the third upwards) is reduced by several
orders of magnitude. Moreover, it is shown that the effect is robust against
important variations of the laser field and to temperature changes, providing
a more robust probe of magic angle over the existing ones.To check this
hypothesis, we performed numerical simulations (see Methods for more details
on the model and the time propagation scheme). To test the effect of the flat
bands on the emission spectra, we will scan a set of twist angles, near
$\theta_{\text{mag}}$, ranging from $1.00^{\circ}$ to $1.35^{\circ}$. In Fig.
2, we show the band structure for 3 different angles: the magic angle itself
($\theta=1.12^{\circ}$ in our model), one angle above ($\theta=1.22^{\circ}$)
and one below it ($\theta=1.05^{\circ}$). One can notice how the bands lose
their flatness when we move away from $\theta_{\text{mag}}$.

Figure 3(a) shows the harmonic spectrum of the TBG for a laser in the
$\Gamma-M$ direction (see Fig. 1) as a function of the twist angle. A main
feature emerge: a strong suppression of the intensity of the odd non-linear
harmonics for twist angles near $\theta_{{\rm mag}}$. In order to obtain a
clearer grasp of the magnitude of this decrease, Fig. 3(b) presents the
harmonic yield from the third up to the seventh harmonic. In spite of the
small angle differences, the intensities of each of the three harmonics differ
by several orders of magnitude. This disparity would allow for all-optical
characterization of the magic-angle twisted bilayer graphene, even for twist
angles close to each other; see for instance, the difference between
$1.06^{\circ}$ versus $1.12^{\circ}$.

Even though the process here described is rather complex, the physical
intuition behind it, is not. At a semiclassical level, the intraband current
can be expressed as (Vampa  _et al._ , 2014)

$\mathbf{j}_{\text{intra}}(t)\propto\sum_{\mathbf{k},n}\nabla_{\mathbf{k}}\varepsilon_{n}(\mathbf{k}(t)),$
(1)

where $\varepsilon_{n}(\mathbf{k}(t))$ is the energy of the $n$-th band and
the crystal momentum is given by $\mathbf{k}(t)=\mathbf{k}+e\mathbf{A}(t)$.
Therefore, the most relevant contribution to the emission spectra is the
_curvature_ of the bands. Hence, when the band starts to flatten, the harmonic
intensity will be reduced. While this semiclassical explanation cannot account
for all the subtle, yet interesting, physics of the system, it provides an
understandable picture of this process.

Consequently, this phenomenon must be robust under important variations of the
laser because it is only based on the existence of the flat bands. This
intuition is confirmed by our numerical simulations. Figure 3(c,d) shows the
same quantities as in 2(a,b) but for a laser in the $\Gamma-K$ direction, see
Fig. 1. Independently of the laser direction, the same physics appears; when
the twist angle gets close to $\theta_{\text{mag}}$, a depression of several
orders of magnitude appears in the emission spectrum.

As seen above, the suppression of the harmonic intensity close to the magic
angle crucially depends on the narrowing of the flat bands but, what happens
when other bands are involved on the dynamics? For instance, when the upper
conduction bands have some population due to an increase in the temperature,
see Fig. 2. Intuitively, one would expect that when upper bands get populated,
the current will have relevant intraband contributions and hence, the decrease
on the emission will be smoothed out, as the harmonic signal will not be
exclusively originated from the flat bands.

To check this hypothesis, we performed a wide temperature scan for a set of
angles near $\theta_{\text{mag}}$. Fig. 4 shows the yield for the 3rd harmonic
in terms of the twist angle and the temperature. For temperatures below $100$
K, the depression on the intensity is heavily located around the magic angle.
As seen from the location of black circles which signify the minimum of the
3rd harmonic, the minimum is located at the magic angle ($1.12^{\circ}$).
However, this depression near $\theta_{\text{mag}}$ disappears when the
temperature is enough to populate the upper bands, and the minimum moves away
from the magic-angle. This effect is still visible for temperatures that are
well above the ones at which the correlated insulators and superconductivity
have been observed in MATBG.

For low frequency laser pulses, the most important contributions to the
current, and therefore to the harmonic signal, come from bands near the Fermi
energy, in the case of TBG, the flat bands. In particular, by choosing laser
frequencies that are smaller than the typical bandwidth of the flat bands in
TBG (here $\omega=1\text{ meV}$), spurious contributions from higher
conductions bands can be avoided. Such small laser frequencies are crucial for
our study and to detect experimentally the drop in the harmonic efficency due
to the flattening of the bands. Laser frequencies and intensities similar to
the ones considered here have been used in high harmonic generation
measurements on a graphene monolayer (Hafez  _et al._ , 2018), making the
experiments proposed experimentally feasible.

Figure 4: Yield of the third harmonic in terms of temperature. Black circles
depicts the minimum of the yield for each twist angle.

## discussion

Our results have shown that the flatness of the bands of MATBG has a strong
impact on its high harmonic spectrum. Their narrowing close to the magic angle
produces a decrease of several orders of magnitude in the emission intensity.
Such discrepancy in the optical response can be used in a fully optical study
of MATBG. This effect is crucially based _only_ on the existence of the flat
bands and it holds for a wide range of laser parameters, even showing no
anisotropy in the laser direction. In addition, we have confirmed that the
suppression is maintained within a relatively wide range of temperatures up
until 100 K, far away from the ultralow temperature regime used in most MATBG
experiments. Nevertheless, we hope that our work can not only be useful as a
spectroscopy tool but also, shed light onto the mechanisms behind electronic
ultrafast dynamics in twisted materials.

## Methods

To model TBG we start from the continuum model (Lopes dos Santos  _et al._ ,
2007; Bistritzer and MacDonald, 2011; San-Jose  _et al._ , 2012) taking a
graphene Fermi velocity of $v_{F}=2.34a$ eV, where $a$ is the lattice constant
of graphene, and a ratio between the interlayer tunneling at AA and AB regions
of $w_{0}/w_{1}=0.78$ (Nam and Koshino, 2017). With these parameters the magic
angle is found at $\theta\approx 1.12^{\circ}$. In the continuum model the two
valleys of graphene are assumed to be uncoupled. The spectrum of TBG equals
the sum of the contributions of each valley. In each of the valleys the model
satisfies the symmetries: $\mathcal{C}_{3}$, $\mathcal{C}_{2}\mathcal{T}$ (the
combined operation involving time reversal symmetry $\mathcal{T}$ and
$\mathcal{C}_{2}$), and a valley preserving mirror symmetry $M$. Here
$\mathcal{C}_{3}$ and $\mathcal{C}_{2}$, involve respectively $120^{\circ}$
and $180^{\circ}$ rotations with respect to an axis perpendicular to the TBG
and $M$ is a layer exchanging two fold rotation around an axis parallel to the
TBG. We then adapt the Wannier function model of TBG with 8 orbitals per
valley (Carr  _et al._ , 2019; Calderón and Bascones, 2020) to approximate the
bandstructure model obtained with the continuum model. The resulting tight
binding model includes hopping up to a radius of 10 moiré lattice constant.
This Wannier model satisfies the symmetries of the continuum model in each
valley. The TBG is assumed to be undoped for all the calculations by setting
the chemical potential for each temperature. Since we are interested in the
optical response of the bilayer, we need to assume some form for the position
operator $\mathbf{r}$, and we restrain ourselves to the diagonal approximation
(nez Azpiroz  _et al._ , 2022; Silva  _et al._ , 2019b).

Inversion and time reversal symmetry are satisfied only when the two valleys
are included. The spectrum of the even harmonics is finite for each valley but
it vanishes when adding the contribution of both valleys.

The time-dependent calculations were performed by solving the semiconductor
Bloch equations (SBE), using Maximally Localized Wannier functions (Marzari
_et al._ , 2012) that are well suited for the model describe above, following
the formalism presented in (Silva  _et al._ , 2019b). During this work, we
have considered the interaction of the twisted bilayer graphene with a
relatively weak $(E_{0}=1\cdot 10^{5}\text{ V/m})$ pulse in the microwave
regime $(\omega=1\text{ meV})$ and a gaussian envelope with a full width at
half-maximum of $23.5$ ps. This optical regime was chosen so to properly
isolate the flat bands, i.e., to avoid transitions between them and the upper
conduction bands. Otherwise, such processes would give relevant interband
contributions to the spectrum. Nevertheless, the physical behavior is robust
under variations of both the amplitude and the frequency of the field. We note
that, even though the field amplitude may seem weak, for the case of MATBG it
will be relatively strong due to the reduced bandwidth of the system. The
reduced density matrix was propagated in the Wannier gauge using a $100\times
100$ Monkhorst–Pack grid and a time step of $2.41$ fs. We checked that the
dynamics were properly converged for all numerical parameters. The dephasing
term (Vampa  _et al._ , 2014; Silva  _et al._ , 2019b) was set to $T_{2}=1000$
fs $\approx\Omega/4$, where $\Omega$ is the laser period. However, we checked
that the same results were obtained if the value of $T_{2}$ was varied.

## Acknowledgments

E. B. M. and R. E. F. S. acknowledge support from the fellowship
LCF/BQ/PR21/11840008 from “La Caixa” Foundation (ID 100010434). This research
was supported by Grant PID2021-122769NB-I00 funded by
MCIN/AEI/10.13039/501100011033.

## References

  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Andrei  _et al._ (2021) E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, Nature Reviews Materials 6, 201 (2021).
  * Andrei and MacDonald (2020) E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020).
  * Lu  _et al._ (2019) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019).
  * Yankowitz  _et al._ (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Nuckolls  _et al._ (2020) K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 588, 610 (2020).
  * Choi  _et al._ (2021) Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Nature 589, 536 (2021).
  * Xie  _et al._ (2021a) Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Nature 600, 439 (2021a).
  * Sharpe  _et al._ (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019).
  * Cao  _et al._ (2020) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * Jaoui  _et al._ (2022) A. Jaoui, I. Das, G. Di Battista, J. Díez-Mérida, X. Lu, K. Watanabe, T. Taniguchi, H. Ishizuka, L. Levitov, and D. K. Efetov, Nature Physics 18, 633 (2022).
  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Fernandes  _et al._ (2022) R. M. Fernandes, A. I. Coldea, H. Ding, I. R. Fisher, P. J. Hirschfeld, and G. Kotliar, Nature 601, 35 (2022).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Lisi  _et al._ (2021) S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov, J. R. Duran, F. Margot, I. Cucchi, E. Cappelli, A. Hunter, A. Tamai, V. Kandyba, A. Giampietri, A. Barinov, J. Jobst, V. Stalman, M. Leeuwenhoek, K. Watanabe, T. Taniguchi, L. Rademaker, S. J. van der Molen, M. P. Allan, D. K. Efetov, and F. Baumberger, Nature Physics 17, 189 (2021).
  * Lopes dos Santos  _et al._ (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Choi  _et al._ (2019) Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. Nadj-Perge, Nature Physics 15, 1174 (2019).
  * (19) A. Datta, M. Calderón, A. Camjayi, and E. Bascones, “Heavy quasiparticles and cascades without symmetry breaking in twisted bilayer graphene,” ArXiv:2301.13024.
  * San-Jose  _et al._ (2012) P. San-Jose, J. González, and F. Guinea, Phys. Rev. Lett. 108, 216802 (2012).
  * Kruchinin  _et al._ (2018) S. Y. Kruchinin, F. Krausz, and V. S. Yakovlev, Rev. Mod. Phys. 90, 021002 (2018).
  * Krausz and Ivanov (2009) F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163 (2009).
  * Ghimire and Reis (2018) S. Ghimire and D. A. Reis, Nature Physics 15, 10 (2018).
  * Goulielmakis and Brabec (2022) E. Goulielmakis and T. Brabec, Nature Photonics 16, 411 (2022).
  * Vampa  _et al._ (2015a) G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, and P. B. Corkum, Nature 522, 462 (2015a).
  * Silva  _et al._ (2018) R. E. F. Silva, I. V. Blinov, A. N. Rubtsov, O. Smirnova, and M. Ivanov, Nature Photonics 12, 266 (2018).
  * Bionta  _et al._ (2021) M. R. Bionta, E. Haddad, A. Leblanc, V. Gruson, P. Lassonde, H. Ibrahim, J. Chaillou, N. Émond, M. R. Otto, A. Jiménez-Galán, R. E. F. Silva, M. Ivanov, B. J. Siwick, M. Chaker, and F. m. c. Légaré, Phys. Rev. Research 3, 023250 (2021).
  * Silva  _et al._ (2019a) R. E. F. Silva, Á. Jiménez-Galán, B. Amorim, O. Smirnova, and M. Ivanov, Nature Photonics 13, 849 (2019a).
  * Chacón  _et al._ (2020) A. Chacón, D. Kim, W. Zhu, S. P. Kelly, A. Dauphin, E. Pisanty, A. S. Maxwell, A. Picón, M. F. Ciappina, D. E. Kim, C. Ticknor, A. Saxena, and M. Lewenstein, Phys. Rev. B 102, 134115 (2020).
  * Uzan-Narovlansky  _et al._ (2022) A. J. Uzan-Narovlansky, Á. Jiménez-Galán, G. Orenstein, R. E. F. Silva, T. Arusi-Parpar, S. Shames, B. D. Bruner, B. Yan, O. Smirnova, M. Ivanov, and N. Dudovich, Nature Photonics 16, 428 (2022).
  * Jürgens  _et al._ (2020) P. Jürgens, B. Liewehr, B. Kruse, C. Peltz, D. Engel, A. Husakou, T. Witting, M. Ivanov, M. J. J. Vrakking, T. Fennel, and A. Mermillod-Blondin, Nature Physics 16, 1035 (2020).
  * Vampa  _et al._ (2014) G. Vampa, C. R. McDonald, G. Orlando, D. D. Klug, P. B. Corkum, and T. Brabec, Phys. Rev. Lett. 113, 073901 (2014).
  * Boyd (2008) R. W. Boyd,  _Nonlinear optics_ (Academic press, 2008).
  * Bernevig  _et al._ (2021a) B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205411 (2021a).
  * Song  _et al._ (2021) Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 205412 (2021).
  * Bernevig  _et al._ (2021b) B. A. Bernevig, Z.-D. Song, N. Regnault, and B. Lian, Phys. Rev. B 103, 205413 (2021b).
  * Lian  _et al._ (2021) B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 103, 205414 (2021).
  * Bernevig  _et al._ (2021c) B. A. Bernevig, B. Lian, A. Cowsik, F. Xie, N. Regnault, and Z.-D. Song, Phys. Rev. B 103, 205415 (2021c).
  * Xie  _et al._ (2021b) F. Xie, A. Cowsik, Z.-D. Song, B. Lian, B. A. Bernevig, and N. Regnault, Phys. Rev. B 103, 205416 (2021b).
  * Ikeda (2020) T. N. Ikeda, Phys. Rev. Research 2, 032015 (2020).
  * Mrudul and Dixit (2021) M. S. Mrudul and G. Dixit, Phys. Rev. B 103, 094308 (2021).
  * Du  _et al._ (2021) M. Du, C. Liu, Z. Zeng, and R. Li, Phys. Rev. A 104, 033113 (2021).
  * Yang  _et al._ (2020) F. Yang, W. Song, F. Meng, F. Luo, S. Lou, S. Lin, Z. Gong, J. Cao, E. S. Barnard, E. Chan, L. Yang, and J. Yao, Matter 3, 1361 (2020).
  * Tancogne-Dejean  _et al._ (2017) N. Tancogne-Dejean, O. D. Mücke, F. X. Kärtner, and A. Rubio, Phys. Rev. Lett. 118, 087403 (2017).
  * Silva  _et al._ (2019b) R. E. F. Silva, F. Martín, and M. Ivanov, Phys. Rev. B 100, 195201 (2019b).
  * Vampa  _et al._ (2015b) G. Vampa, T. J. Hammond, N. Thiré, B. E. Schmidt, F. Légaré, C. R. McDonald, T. Brabec, D. D. Klug, and P. B. Corkum, Phys. Rev. Lett. 115, 193603 (2015b).
  * Hafez  _et al._ (2018) H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Müllen, M. Bonn, M. Gensch, and D. Turchinovich, Nature 561, 507 (2018).
  * Nam and Koshino (2017) N. N. T. Nam and M. Koshino, Phys. Rev. B 96, 075311 (2017).
  * Carr  _et al._ (2019) S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, Phys. Rev. Research 1, 033072 (2019).
  * Calderón and Bascones (2020) M. J. Calderón and E. Bascones, Phys. Rev. B 102, 155149 (2020).
  * nez Azpiroz  _et al._ (2022) J. I. nez Azpiroz, F. de Juan, and I. Souza, SciPost Phys. 12, 070 (2022).
  * Marzari  _et al._ (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Reviews of Modern Physics 84, 1419 (2012).


# NEATH III: a molecular line survey of a simulated star-forming cloud

F. D. Priestley1, P. C. Clark1, S. C. O. Glover2, S. E. Ragan1, O. Fehér1, L.
R. Prole3, R. S. Klessen2,4
1School of Physics and Astronomy, Cardiff University, Queen’s Buildings, The
Parade, Cardiff CF24 3AA, UK
2Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Straße 2, D-69120 Heidelberg, Germany
3Centre for Astrophysics and Space Science Maynooth, Department of Theoretical
Physics, Maynooth University, W23 F2H6 Maynooth, Ireland
4Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches
Rechnen, Im Neuenheimer Feld 205, D-69120 Heidelberg, Germany
Email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present synthetic line observations of a simulated molecular cloud,
utilising a self-consistent treatment of the dynamics and time-dependent
chemical evolution. We investigate line emission from the three most common CO
isotopologues (12CO, 13CO, C18O) and six supposed tracers of dense gas (NH3,
HCN, N2H+, HCO+, CS, HNC). Our simulation produces a range of line intensities
consistent with that observed in real molecular clouds. The HCN-to-CO
intensity ratio is relatively invariant with column density, making HCN (and
chemically-similar species such as CS) a poor tracer of high-density material
in the cloud. The ratio of N2H+ to HCN or CO, on the other hand, is highly
selective of regions with densities above $10^{22}\,{\rm cm}^{-2}$, and the
N2H+ line is a very good tracer of the dynamics of high volume density
($>10^{4}\,{\rm cm}^{-3}$) material. Focusing on cores formed within the
simulated cloud, we find good agreement with the line intensities of an
observational sample of prestellar cores, including reproducing observed CS
line intensities with an undepleted elemental abundance of sulphur. However,
agreement between cores formed in the simulation, and models of isolated cores
which have otherwise-comparable properties, is poor. The formation from and
interaction with the large-scale environment has a significant impact on the
line emission properties of the cores, making isolated models unsuitable for
interpreting observational data.

###### keywords:

astrochemistry – stars: formation – ISM: molecules – ISM: clouds

††pubyear: 2024††pagerange: NEATH III: a molecular line survey of a simulated
star-forming cloud–B

##  1 Introduction

Star formation occurs in molecular clouds, where high gas densities and strong
shielding from the external ultraviolet (UV) radiation field result in most
hydrogen existing in the form of H2 (Bergin & Tafalla, 2007). Temperatures in
molecular clouds are too low to excite any significant line emission from H2,
making the bulk of the star-forming gas essentially invisible. Studies of star
formation therefore rely on indirect tracers of the gas mass, such as far-
infrared (far-IR) thermal emission from dust grains (e.g. Könyves et al.,
2015) or extinction measurements of background stars (Zucker et al., 2021). Of
particular utility are rotational emission lines from molecular species at
millimetre wavelengths. Line emission provides information on the line-of-
sight velocity of the gas, crucial for assessing its dynamics and stability,
and the large number of lines detectable with modern observational facilities
(e.g. Pety et al., 2017; Kauffmann et al., 2017; Tafalla et al., 2021) can be
used to estimate properties such as volume density and temperature, exploiting
the fact that the excitation of each line has a unique response to the local
physical conditions (Shirley, 2015).

Analysis of molecular line data is complicated by the fact that the observed
emission depends on a complex combination of the physical structure of the
cloud, radiative transfer effects, and the chemical composition of the gas.
Disentangling all of these factors to accurately reconstruct detailed cloud
properties from observational data is not practically feasible. Many studies
have therefore investigated forward modelling the problem and producing
synthetic line observations of simulated clouds, which can be used to help
inform a physical interpretation of real observational data. To date, most of
these studies have focused on lines from important coolant species, in
particular CO (e.g. Peñaloza et al., 2017; Seifried et al., 2017; Peñaloza et
al., 2018; Clarke et al., 2018; Clark et al., 2019), as its chemical evolution
is often followed self-consistently in modern hydrodynamical simulations.
However, due to its high abundance and correspondingly high line optical
depths, CO is a poor tracer of the high-density material where star formation
actually occurs in molecular clouds (Clark et al., 2019; Priestley et al.,
2023a). Isotopologues of CO, such as C18O, have much lower optical depths and
hence are better tracers of dense gas, but they are also faint and hard to
observe in distant clouds or extragalactic systems.

Studies investigating lines from rarer molecules such as HCN and N2H+, which
can probe denser gas more effectively, generally assume that the molecular
abundances are constant, or can be specified as a function of the local
physical properties (e.g. Offner et al., 2008; Smith et al., 2012, 2013; Chira
et al., 2014; Jones et al., 2023; Jensen et al., 2023). This approach ignores
the drastic variations in abundance driven by differing evolutionary
histories, which become particularly important at these high densities
(Priestley et al., 2023c). Obtaining reliable line emission properties from
simulated molecular clouds requires the physical and chemical evolution of the
clouds to be treated self-consistently.

In this paper, we achieve this by performing radiative transfer simulations of
a cloud for which the molecular composition has been obtained from a full
time-dependent chemical network. The abundances at each point in the cloud are
thus consistent with the physical evolution of the material that ended up at
that point, making the resulting line emission data genuinely comparable to
that from real molecular clouds with their own complex formation histories. We
use this to assess the degree of correspondance between our simulated cloud
and reality, the limitations of smaller-scale models in interpreting line
emission, and the accuracy of commonly-used observational diagnostics of the
physical structure of clouds.

##  2 Method

Table 1: Elemental abundances used in the chemical modelling. Element | Abundance |  | Element | Abundance
---|---|---|---|---
C | $1.4\times 10^{-4}$ |  | S | $1.2\times 10^{-5}$
N | $7.6\times 10^{-5}$ |  | Si | $1.5\times 10^{-7}$
O | $3.2\times 10^{-4}$ |  | Mg | $1.4\times 10^{-7}$
Table 2: Molecules investigated, collisional partners, and sources of collisional rate data. Molecule | Partners | Reference
---|---|---
12CO | p-H2, o-H2 | Yang et al. (2010)
13CO | p-H2, o-H2 | Yang et al. (2010)
C18O | p-H2, o-H2 | Yang et al. (2010)
p-NH3 | p-H2 | Loreau et al. (2023)
HCN | p-H2, e- | Faure et al. (2007)
Magalhães et al. (2018)
N2H+ | p-H2 | Lique et al. (2015)
HCO+ | p-H2, o-H2 | Denis-Alpizar et al. (2020)
CS | H2 | Lique et al. (2006)
HNC | H2 | Dumouchel et al. (2010)

###  2.1 Hydrodynamical simulation

We simulate the formation and subsequent evolution of a molecular cloud via
the collision of two spherical, initially-atomic gas clouds. The simulations
are performed using the magnetohydrodynamic (MHD) moving-mesh code arepo
(Springel, 2010; Pakmor et al., 2011), modified to properly capture the
thermodynamics of the gas and dust (Glover & Mac Low, 2007; Glover & Clark,
2012). This includes a simplified chemical network for H2 and CO (based on
Gong et al. 2017, with some additions and modifications described in Hunter et
al. 2023), and a self-consistent treatment of shielding from the background UV
radiation field (Clark et al., 2012a).

The initial conditions of the simulation are two spherical $10^{4}\,{\rm
M}_{\odot}$ clouds with radii $R=19\,{\rm pc}$, giving an initial volume
density of hydrogen nuclei $n_{\rm H}=10\,{\rm cm}^{-3}$. The gas and dust
temperatures are initialised to $300\,{\rm K}$ and $15\,{\rm K}$ respectively.
The clouds are displaced in the $\pm x$ direction by $R$, so they are just in
contact at the outset, and given a bulk motion along the $x$-axis of $\mp
7\,{\rm km\,s^{-1}}$. Each cloud also has a virialised turbulent velocity
field with 3D velocity dispersion $\sigma=0.95\,{\rm km\,s^{-1}}$, and a
$3\,{\rm\mu G}$ magnetic field is present parallel to the collision axis. In
collapsing high-density regions of the simulation, sink particles are
introduced (following Tress et al. 2020; see also Prole et al. 2022), with a
threshold density of $2\times 10^{-16}\,{\rm g\,cm}^{-3}$ and a formation
radius of $9\times 10^{-4}\,{\rm pc}$. The simulation is run for $5.53\,{\rm
Myr}$, at which point sink particles have accreted $102\,{\rm M}_{\odot}$ of
material ($0.5\%$ of the original cloud masses); beyond this point, feedback
effects from newly-formed stars (which are not modelled) are likely to become
important (Whitworth, 1979; Walch et al., 2012).

As the mesh cells can be created and destroyed, and so do not correspond to
coherent parcels of gas, we follow the physical evolution of representative
gas parcels using Monte Carlo tracer particles (Genel et al., 2013), which
provide the input for the chemical modelling described below. Global
properties, shared between the MHD and chemical models, are the UV radiation
field strength of $1.7$ times the Habing (1968) field, the cosmic ray
ionisation rate111Our simulations do not include attenuation of cosmic rays,
so the ionisation rate is constant with density. The observed rate in high-
density regions of molecular clouds (Pineda et al., 2024; Redaelli et al.,
2024) is typically somewhat lower than the value we have chosen
(representative of the diffuse interstellar medium; Indriolo et al. 2015),
which may be an overestimate for the corresponding regions in the simulation.
of $10^{-16}\,{\rm s^{-1}}$ per H atom, the dust-to-gas mass ratio of $0.01$,
and the elemental abundances of carbon ($1.4\times 10^{-4}$) and oxygen
($3.2\times 10^{-4}$) from Sembach et al. (2000). The ‘metal’ abundance in the
MHD simulation, representing heavier elements such as silicon, is set to
$10^{-7}$, corresponding to the high levels of depletion found in the dense
interstellar medium (Jenkins, 2009).

###  2.2 Chemical modelling

The tracer particles in the MHD simulation record the density, gas
temperature, and effective shielding column densities at intervals of
$44\,{\rm kyr}$. These particle histories are used to model the chemical
evolution under the NEATH framework222https://fpriestley.github.io/neath/
(Priestley et al., 2023c), which is designed so that the returned H2 and CO
abundances are consistent with those of the underlying MHD chemical network.
We use a modified version of the time-dependent gas-grain code uclchem
(Holdship et al., 2017) to evolve the UMIST12 (McElroy et al., 2013) reaction
network, with the same UV and cosmic ray values as the MHD simulation. The
assumed elemental abundances are listed in Table 1. We take carbon, nitrogen
and oxygen abundances from Sembach et al. (2000), and deplete silicon and
magnesium by factors of $100$ from their warm neutral medium (WNM) values,
again representing depletion into refractory dust grains. Unlike most
astrochemical models of star-forming clouds and cores, we do not deplete
sulphur from its WNM value, which is close to the undepleted Solar value; we
return to this point in Section 4.2.

We follow the chemical evolution of $10^{5}$ tracer particles chosen from the
central $16.2\,{\rm pc}$ of the computational domain at the simulation’s end,
which contains virtually all of the molecular material. Particles are selected
randomly to evenly sample densities between $10-10^{6}\,{\rm cm}^{-3}$, the
highest density for which the chemical evolution is converged (Priestley et
al., 2023c). Lower-density material has a negligible molecular content, and
higher-density material makes up a negligible fraction of the cloud mass. We
confirm that this number and selection of particles is sufficient for our
purposes in Appendix A.

###  2.3 Radiative transfer

We perform radiative transfer modelling of our simulated cloud using radmc3d
(Dullemond et al., 2012). The unstructured Voronoi mesh of the MHD simulation
is interpolated onto a cubic adaptive mesh, with side length $32.4\,{\rm pc}$
and a base resoluton of $20^{3}$. Any cell which contains more than one
Voronoi sampling point is refined into eight subcells, and this process is
repeated until all cells hold at most one sampling point each. Cells are then
assigned physical properties (density, velocity, gas and dust temperatures)
from the nearest Voronoi sampling point, and molecular abundances from the
nearest post-processed tracer particle. We assess the accuracy of the
interpolation in Appendix A. The output position-position-velocity (PPV) cubes
have a spatial resolution of $0.06\,{\rm pc}$, and a velocity resolution of
$0.03\,{\rm km\,s^{-1}}$ for lines without hyperfine structure (lines with
hyperfine structure are discussed below).

We use collisional rate data taken from the LAMDA (Schöier et al., 2005) and
EMAA databases, for the molecules and collision partners listed in Table 2.
Where collisional data for both ortho- and para-H2 are available, we assume an
ortho:para ratio of $3:1$, and we assume para-NH3 makes up half the total
abundance of NH3. As our reaction network does not include separate isotopic
chemistry, we assume, where necessary, that the ratio of molecular
isotopologues is equal to the underlying isotope ratio, namely 12C/13C$=77$
and 12O/18O$=550$ (Wilson & Rood, 1994). Deviations from these ratios under
the conditions we investigate are unlikely to be larger than a factor of a few
(Szűcs et al., 2014). We assume thick ice mantle dust opacities from Ossenkopf
& Henning (1994) for wavelengths $\lambda>1\,{\rm\mu m}$, and Mathis et al.
(1983) opacities at shorter wavelengths (see Clark et al. 2012b).

We produce PPV line intensity cubes for the $J=1-0$ rotational transitions of
12CO, 13CO, C18O, HCN, N2H+, HCO+ and HNC, and the $J=2-1$ transition of CS.
These lines all fall within the $3\,{\rm mm}$ observational window, and are
typically among the brightest in star-forming regions, leading to plentiful
observational data as a point of comparison (Pety et al., 2017; Kauffmann et
al., 2017; Barnes et al., 2020; Tafalla et al., 2021, 2023; Hacar et al.,
2024). We also investigate the para-NH3 $(1,1)$ inversion transition, for
which similarly rich observational data is available (Ragan et al., 2011;
Ragan et al., 2012b; Friesen et al., 2017; Fehér et al., 2022).

The NH3, HCN and N2H+ transitions involve multiple hyperfine components, and
we use larger velocity resolutions of $0.25$, $0.1$ and $0.1\,{\rm
km\,s^{-1}}$ respectively to capture all the structure without increasing the
computational cost of the radiative transfer. While radmc3d does account for
the optical depth effects of multiple overlapping hyperfine transitions, this
is not included when determining level populations, meaning that the hyperfine
levels are not radiatively coupled with each other. The distribution of flux
between the individual components may therefore be calculated incorrectly if
this coupling is significant, although we expect the integrated intensities
over all hyperfine components to be valid.

Figure 1: Column density maps of the cloud seen face-on (left) and edge-on
(right).

Figure 2: Integrated line intensity maps for the cloud seen face-on.

Figure 3: Integrated line intensity versus column density for the cloud seen
face-on (blue) and edge-on (red). Crosses show the median pixel values, error
bars the 16th/84th percentiles. Observational data from the Perseus molecular
cloud (Friesen et al., 2017; Tafalla et al., 2021) are shown in black.

Figure 4: Maps of the ratio of integrated line intensities for HCN and 12CO
(left), N2H+ and 12CO (centre), and N2H+ and HCN (right), for the cloud seen
face-on. Pixels with integrated intensities of either line below $0.2\,{\rm
K}\,{\rm km\,s^{-1}}$ have been masked.

Figure 5: Ratio of integrated line intensities versus column density for HCN
and 12CO (left), N2H+ and 12CO (centre), and N2H+ and HCN (right), for the
cloud seen face-on (blue) and edge-on (red). Crosses show the median pixel
values, error bars the 16th/84th percentiles.

Figure 6: Cloud-averaged line profiles for 12CO (blue), HCN (red) and N2H+
(green) seen face-on, compared to the velocity distribution of all cloud mass
(black, left) and mass above a density of $10^{4}\,{\rm cm}^{-3}$ (black,
right). All quantities have been normalised by their peak values. Figure 7:
Positions of cores selected from the face-on column density map, highlighted
with white circles. Numbers correspond to the labelling in Figure 8. For the
purpose of clarity, the circles are three times larger than the $0.1\,{\rm
pc}$ extraction regions.

Figure 8: Line profiles of C18O (blue), HCN (red) and N2H+ (green) for the
cores identified from the simulated cloud in Figure 7. Text labels indicate
the column density of the $0.1\,{\rm pc}$ extraction region.

Figure 9: Distribution of the peak line intensities (left) and FWHMs (right)
of CS versus N2H+. Blue crosses show the cores identified from the simulated
cloud in Figure 7. Green circles show observational data taken from Lee et al.
(1999) and Tafalla et al. (2002). Orange triangles show the results of
initially-spherical core models from Priestley et al. (2022).

##  3 Results

###  3.1 Cloud-scale line properties

Figure 1 shows column density maps in total hydrogen nuclei of the simulated
cloud, viewed both parallel (‘face-on’) and perpendicular (‘edge-on’) to the
collision axis. The collision between the initial atomic clouds forms a layer
of enhanced density at their interface, which subsequently fragments into a
network of filaments and cores. The initial turbulent velocity field causes
deviations from spherical symmetry, such as the ‘bend’ in the cloud when
viewed edge-on. Figure 2 shows integrated intensity maps of lines from the
molecules listed in Table 2, for the face-on cloud orientation (line emission
maps of the edge-on cloud show the same qualitative behaviour).

The CO isotopologues, with their relatively high abundances and low critical
excitation densities, show widespread emission throughout the cloud, although
this is less pronounced for the rarest (C18O). This behaviour is also seen for
several species commonly thought of as tracers of high-density material, most
notably HCN, CS and HNC, as is found observationally (Pety et al., 2017;
Kauffmann et al., 2017; Evans et al., 2020). The line emission maps of NH3 and
HCO+ more closely track the actual column density of material shown in Figure
1, while N2H+ is the only species investigated which selectively traces the
densest regions of the cloud. These species’ tendency to trace higher density
gas than CO is due to chemical differences causing their abundances to peak at
higher volume densities, whereas species such as CS and HCN are chemically
quite similar to CO (Priestley et al., 2023c).

Figure 3 shows the relationships between integrated line intensity and column
density, compared to observations of the Perseus molecular cloud. For all
lines except NH3, the data are the pointed observations taken by Tafalla et
al. (2021). For NH3, we use maps of the NGC 1333 subregion from Friesen et al.
(2017), which have been re-reduced (Pineda et al. in prep.) and matched to the
Herschel-derived column density333Converted from H2 column to total hydrogen
column assuming all hydrogen is in the form of H2. maps of Singh & Martin
(2022). The noise level of the NH3 observations is $\sim 0.5\,{\rm K}$, so the
apparently-constant intensity at around this value below a column density of
$10^{22}\,{\rm cm}^{-2}$ should be interpreted as an upper limit. We note that
the commonly-used approach of fitting single-temperature modified blackbodies
to far-IR data may introduce an observational bias in the measured column
densities, but in Appendix B we demonstrate that this effect is modest (a
factor of a few at the highest column density values), so we use the true
column densities from the simulation throughout the rest of this paper.

The predicted relationships between column density and line intensity for the
three CO isotopologues are in general agreement with those observed, although
the 12CO line is somewhat weaker than observed. A possible explanation for
this difference is that our assumed radiation field strength of $1.7$ Habing
may be too low to properly represent the Perseus molecular cloud; Tafalla et
al. (2023) find that regions with more active star formation produce more 12CO
emission at a given column density, which they attribute to higher gas
temperatures caused by a stronger local UV field. The predicted strengths of
the 13CO and C18O lines, which originate from further within the cloud, are in
better agreement with the data, suggesting that this effect, if relevant, is
only important in the outer, less-shielded regions.

Alternatively, the difference in 12CO intensities may be due to systematic
differences in the velocity dispersions of our simulated cloud and Perseus.
Tafalla et al. (2021) find typical 12CO full-widths at half-maximum (FWHMs) of
$3-4\,{\rm km\,s^{-1}}$ in Perseus, with little column density variation. Our
simulated cloud also has a flat FWHM-column density relationship, but at a
lower value of $\sim 1\,{\rm km\,s^{-1}}$. If the velocity dispersion were to
be increased to match the Perseus FWHM values, the increase in integrated
intensity could plausibly reach the factor of $2-3$ required to bring our
simulated data into agreement with the Perseus observations in Figure 3. For
the rarer CO isotopologues, the differences between observed and simulated
FWHMs become less severe ($2$ versus $1\,{\rm km\,s^{-1}}$ for 13CO, and $1$
versus $0.7\,{\rm km\,s^{-1}}$ for C18O) and the line intensities are in
correspondingly better agreement. This again suggests that discrepancies
between our simulated cloud and the Perseus data primarily originate in the
outer, lower-density regions of the cloud.

Line emission from molecules with similar chemical behaviour to CO (HCN, CS,
HNC; Priestley et al. 2023c) differs from the Perseus data in our simulated
observations. The intensities rise sharply with column density up to $\sim
10^{22}\,{\rm cm}^{-2}$ and then saturate at a near-constant value, whereas
the observed behaviour is a continuous linear rise between $10^{21-23}\,{\rm
cm}^{-2}$. For species which preferentially trace denser gas (NH3, N2H+,
HCO+), the simulation is in somewhat better agreement with the data. As with
the CO lines, this suggests that while our simulated cloud differs from
Perseus in some respects along moderate-density sightlines ($\sim
10^{21}\,{\rm cm}^{-2}$), we are accurately capturing the physical and
chemical properties of the densest regions. We note that while the linear
relationship between line intensity and column density has also been observed
in the Orion A and California clouds (Tafalla et al., 2023), similar studies
of Orion B (Pety et al., 2017; Santa-Maria et al., 2023) and W49 (Barnes et
al., 2020) find plateau-like behaviour for HCN and other molecules, which more
closely resembles our simulated cloud. We discuss interpretations of the
unexpected linear relationship for optically-thick lines in Section 4.3.

###  3.2 Tracers of dense gas

Figure 4 shows maps of the intensity ratios of the 12CO, HCN and N2H+ lines
seen face-on, and Figure 5 shows how their values vary with column density.
Edge-on maps, as with the total intensity maps, are qualitatively similar to
their face-on equivalents. The HCN/12CO ratio is frequently used as an
indicator of the fraction of dense gas (with the exact interpretation of
‘dense’ varying) on extragalactic scales, but in our simulated cloud, its
value varies by only a factor of two over two orders of magnitude in column
density, suggesting it is a poor indicator of genuinely-dense material. Ratios
of other lines, such as CS and HCO+, also show relatively little variation
over the entire range of column densities where they are detectable (above
$3\times 10^{21}\,{\rm cm}^{-2}$). As argued by several observational studies
(Pety et al., 2017; Kauffmann et al., 2017; Tafalla et al., 2021), ratios
involving N2H+ are much more sensitive to the presense of material above a
column density of $10^{22}\,{\rm cm}^{-2}$, thought to be related to the onset
of star formation (Lada et al., 2010), although effects not considered in our
simulation such as protostellar heating may complicate this interpretation
(Fehér et al., 2024).

Of particular note, the N2H+/HCN ratio in our cloud rises by a factor of ten
over a relatively small range in column density, making it a sensitive probe
of high-density material. This line ratio has become observationally
accessible on extragalactic scales in recent years, with a typical value of
$\sim 0.2$ (Jiménez-Donaire et al., 2023; Stuber et al., 2023). On the scale
of our simulated cloud, this value corresponds closely to the $10^{22}\,{\rm
cm}^{-2}$ star formation threshold found by Lada et al. (2010), suggesting
that a significant fraction of the $\,{\rm kpc}$-scale beam area in these
extragalactic studies is made up of actively star-forming gas.

The inadequacy of HCN as a dense gas tracer is also apparent in velocity-
resolved data. Figure 6 shows the face-on line profiles for 12CO, the central
HCN hyperfine component, and the isolated N2H+ component, averaged over the
central $16.2\,{\rm pc}$ region, as compared to the actual line-of-sight
velocity distribution of the mass in this region. Both the 12CO and HCN lines
are quite accurate tracers of the whole-cloud dynamics, and in fact have
almost indistinguishable profile shapes. The N2H+ profile is significantly
narrower, being a much better (although still imperfect) tracer of the
dynamics of the high-density ($>10^{4}\,{\rm cm}^{-3}$) material.

The difference in the velocity dispersions measured via different tracers has
been noted previously in observations (Goodman et al., 1998; Pineda et al.,
2010; Ragan et al., 2015; Peretto et al., 2023), and appears to represent a
real reduction in the magnitude of turbulent motions in the densest material,
possibly due to dissipation in accretion shocks (Klessen et al., 2005;
Priestley et al., 2023b). This means that the observational definition of
velocity dispersion, used to determine quantities such as virial stability
(e.g. Rigby et al., 2024) and magnetic field strength (Wang et al., 2024), is
tracer-dependent in a highly non-trivial manner. A structure’s virial ratio
determined via the HCN line will be different to the one calculated using the
N2H+ line, and the ‘correct’ value to use will depend entirely on the
underlying science question. The association between a line’s velocity
structure and the underlying gas dynamics should be made with considerable
caution.

###  3.3 Core-scale line properties

To investigate the line emission properties on the scales of prestellar cores,
we identify nine peaks in the face-on column density map, shown in Figure 7,
and extract line profiles from $0.1\,{\rm pc}$ apertures around these peaks,
approximating the effects of beam size for single-dish observations of nearby
clouds (e.g. Lee et al., 1999; Tafalla et al., 2002). Figure 8 shows the
resulting line profiles of three species commonly used to trace core-scale
emission: C18O, HCN, and N2H+ (we again focus on the central HCN and isolated
N2H+ hyperfine components). The HCN line is in most cases substantially
broader than the other two, and often shows evidence of multiple velocity
components which are not visible in the other two lines. Rather than
selectively tracing the densest ‘core’ material, HCN emission also originates
from gas with significantly lower densities ($\ll 10^{4}\,{\rm cm}^{-3}$;
Evans et al. 2020; Jones et al. 2023), making it a poor tracer of the cores’
dynamics. Several cores also show evidence of multiple velocity components in
the N2H+ profiles, suggestive of line-of-sight confusion between physically-
distinct regions of dense gas. These features are rarely seen in observational
samples of nearby cores with comparable ($\sim 0.1\,{\rm pc}$) spatial
resolution (Lee et al., 1999; Tafalla et al., 2002), although complex multi-
component N2H+ profiles on similar spatial scales do occur in more distant
infrared-dark clouds (Rigby et al., 2024).

Figure 9 shows the peak line intensity of CS versus that of N2H+ for the
simulated cloud core sample, compared to observational data from Lee et al.
(1999) and Tafalla et al. (2002). The ratio of these lines was identified by
Yin et al. (2021) as a diagnostic of the degree of magnetic support, based on
simulations of isolated, initially-spherical cores. Priestley et al. (2022)
found that model cores with supercritical mass-to-flux ratios (Mouschovias &
Spitzer, 1976), also shown in Figure 9, were incapable of matching the
observed distribution of this ratio, as in the absence of magnetic support
collapse proceeds too rapidly to significantly deplete CS from the gas phase,
and its $J=2-1$ line is subsequently stronger than observed. By contrast, the
cores formed in our simulation agree well with the observational data, despite
the initial mass-to-flux ratio of the simulation being $2.4$ times the
critical value (i.e. moderately supercritical). The complex formation
histories of the cores in our simulation, involving material being repeatedly
cycled in and out of high-density phases before undergoing runaway
gravitational collapse (Priestley et al., 2023c), results in them having very
different chemical properties to equivalent models of isolated cores, which
instead experience a monotonic increase in density with time.

Another problematic aspect of supercritical core models is their predicted
line widths, which are inevitably much broader than the observed subsonic
values due to the supersonic infall velocities of unimpeded gravitational
collape (Larson, 1969). Figure 9 shows the FWHMs444In all cases measured
directly from the line profiles, rather than being the FWHM of a Gaussian fit
to the lines. of the CS and N2H+ lines of the cores from our simulation and
from the Priestley et al. (2022) supercritical sample, compared to
observational data from Tafalla et al. (2002). In this case, forming cores
self-consistently does not improve the agreement with observations: simulated
cores still have much broader lines ($>0.5\,{\rm km\,s^{-1}}$) than real ones.
Previous studies (Offner et al., 2008; Smith et al., 2012; Priestley et al.,
2023b) have reproduced subsonic line widths in cores and filaments formed via
supersonic converging flows, which efficiently dissipate the initial kinetic
energy (Klessen et al., 2005), so the discrepancy here may be due to our
initial turbulent velocity field being somewhat subsonic, even though the
cloud collision itself is highly supersonic. Alternatively, simulations may
still require a substantial degree of magnetic support to reproduce the
observed properties of real objects, but it is clearly impossible to make this
argument based on models of isolated spheres, which are far too idealised
compared to the complex formation histories of cores in turbulent molecular
clouds.

##  4 Discussion

Figure 10: Distribution of arepo cells by volume density and column density as
seen face-on (left) and edge-on (right). The dashed black lines show the
relationship proposed by Tafalla et al. (2021) to explain linear correlations
between line intensities and column density.

Figure 11: Line-of-sight velocity dispersion (left) and HCN line intensity
(right) versus column density for the colliding (blue) and non-colliding
(cyan) simulations. Crosses show the median pixel values, error bars the
16th/84th percentiles. Observational data from the Perseus molecular cloud for
HCN (Tafalla et al., 2021) are shown in black.

###  4.1 The (in)accuracy of isolated core models

Simulations of isolated prestellar cores are widely used to interpret
molecular line observations (Keto et al., 2015; Sipilä & Caselli, 2018; Young
et al., 2019; Sipilä et al., 2022; Tritsis et al., 2023; Redaelli et al.,
2024), under the assumption that cores are sufficiently decoupled from their
environment to evolve more-or-less independently. This does not appear to be
the case; Wurster & Rowan (2024) find that resimulating the regions which form
cores within a larger-scale cloud produces significantly different results to
the original simulation, due to the chaotic nature of star formation, while
several observational studies (Peretto et al., 2020; Rigby et al., 2021;
Anderson et al., 2021; Redaelli et al., 2022) have suggested that cores
continue to accrete mass from their surroundings throughout their lifetimes.

The prior evolutionary history of the material making up a core, and
continuing accretion onto it, can significantly alter its chemical composition
(Priestley et al., 2023c), with a correspondingly large effect on the
predicted line emission properties. The peak line intensities of our core
samples in Figure 9 resemble those of magnetically-supported isolated models
from Priestley et al. (2022), despite the opposite being the case in the
simulation, while the subsonic linewidths of observed cores may be a result of
their formation from cloud-scale turbulent motions (Offner et al., 2008; Smith
et al., 2012; Priestley et al., 2023b) rather than the isolated quasi-
equilibrium scenario adopted by Keto et al. (2015) and subsequent other
authors. Taking these highly-idealised models as representative of real cores
risks making fundamental errors in the physical interpretation of molecular
line data.

We note that a similar argument could be made regarding our idealised setup of
two colliding, uniform-density spheres. Real molecular clouds are rarely
spherical (Faustino Vieira et al., 2024), and may be better characterised by a
continuous injection of turbulence (Zhou et al., 2022), rather than the
decaying velocity field we employ. However, the generally-good agreement
between the predicted line emission from our simulations and observational
data, on both cloud and core scales, suggests that they are accurately
capturing the physical-chemical structure of real molecular clouds (see also
Priestley et al. 2023a), despite their idealised initial conditions.
Simulations on the galactic disc scales necessary to self-consistently follow
cloud formation (e.g. Panessa et al., 2023) typically lack the resolution to
investigate chemistry in detail in pre- and protostellar cores. These objects
are likely to be sufficiently decoupled from their galactic environment for
cloud-scale simulations to be an entirely adequate representation of their
formation and evolution.

###  4.2 Sulphur depletion in molecular gas

It is common practice in astrochemical modelling to reduce the elemental gas-
phase sulphur abundance by a factor of ten or more from its Solar value, as
otherwise many sulphur-bearing molecules are predicted to be far more
abundant, and their line emission far stronger, than is actually observed
(Fuente et al., 2019; Navarro-Almaida et al., 2020; Fuente et al., 2023). The
assumption is that most of the sulphur is locked up in some chemically-
unavailable solid phase, but direct measurements of sulphur depletion in the
interstellar medium find a near-Solar gas-phase abundance even along the
densest sightlines (Jenkins, 2009), suggesting that only a modest fraction of
the total sulphur is in solid form.

Our simulation, using an undepleted sulphur abundance (Table 1), successfully
reproduces the observed peak intensities of CS in cores (Figure 9), often one
of the more problematic species in astrochemical modelling (Navarro-Almaida et
al., 2020). As argued by Hily-Blant et al. (2022), the main sulphur reservoir
in molecular clouds may be atomic S, rather than sulphur-bearing molecules
(Holdship et al., 2019; Kushwahaa et al., 2023) and ices (Boogert et al.,
2015; McClure et al., 2023), which make up only a few percent of the Solar
abundance when combined. It was speculated in Priestley et al. (2022) that the
sulphur abundance problem may be resolved by more realistic physical models of
molecular gas, rather than an improved understanding of its chemistry. We will
investigate this possibility, and the detailed behaviour of sulphur-bearing
species in our chemical model, in future work.

###  4.3 Linear intensity scaling in optically-thick lines

The linear relationship between line intensity and column density for
molecules such as HCN in Perseus (Tafalla et al., 2021) and some other
molecular clouds (Tafalla et al., 2023) is unexpected, as these lines are
optically thick (typical line-centre opacities $>10$ for column densities
$>10^{21}\,{\rm cm}^{-2}$). The amount of emitting material along the line-of-
sight should then have little impact on the strength of the observed emission.
Explanations for this behaviour (Tafalla et al., 2021; Dame & Lada, 2023)
typically rely on the presence of a positive correlation between volume and
column densities; if the characteristic volume density is larger along higher-
column sightlines, and the level populations are subthermally-excited, then an
increasing column density can produce an increase in line intensity even in
the optically-thick limit.

Simulated molecular clouds typically do show such a relationship (Clark &
Glover, 2014; Bisbas et al., 2019), but not one capable of reproducing the
linear behaviour observed by Tafalla et al. (2021). Figure 10 shows the
distribution of the volume and line-of-sight column densities of arepo cells
in our simulation. The correlation between column and volume density is much
flatter than the $N_{\rm H}\propto n_{\rm H}^{4/3}$ proposed by Tafalla et al.
(2021), and more significantly, departs severely from a one-to-one
relationship; a given line-of-sight column density can comprise gas covering
several orders of magnitude in volume density. It seems probable that at least
the latter behaviour is shared with real molecular clouds, making the
excitation temperature argument for the linear intensity scaling unviable;
there should be little, if any, correlation between the excitation temperature
of a transition, a complicated mass- and optical depth-weighted average over
material of all densities along the line of sight, and the observed column
density. We note that despite the face-on and edge-on cloud orientations
having very different volume density distributions for line-of-sight columns
above $10^{22}\,{\rm cm}^{-2}$, the resulting line intensities in this regime
in Figure 3 are effectively indistinguishable, suggesting that the
characteristic volume density of a line of sight has a negligible effect on
the resulting line emission.

An alternative explanation for the increasing intensity of optically-thick
lines is if the line-of-sight velocity dispersion is higher along high-density
sightlines, providing more ‘bandwidth’ available for emission (Whitworth &
Jaffa, 2018). Simulations in Priestley et al. (2023a) where the velocity
dispersion increases with column density produced linearly-increasing HCN line
intensities, whereas simulations without a rise in velocity dispersion
resulted in an intensity plateau similar to those in Figure 3. Our simulation
has an almost-constant velocity dispersion over the range of column densities
where HCN is detectable (Figure 11), so no additional velocity channels become
available for emission, and once the HCN line becomes optically thick, its
intensity saturates.

The distinction between rising and flat velocity dispersions in Priestley et
al. (2023a) was attributed to simulations of isolated versus colliding clouds.
However, we do not reproduce this behaviour here. Figure 11 shows the velocity
dispersion and HCN intensity relationships for a simulation where we have
reduced the collision velocity from $7$ to $0.5\,{\rm km\,s^{-1}}$; the
subsonic collision velocity (also below the average velocity of the initial
turbulent field) results in behaviour comparable to an isolated, turbulent
cloud as studied in Priestley et al. (2023a). Rather than increasing, the
velocity dispersion actually declines over the range of column densities
relevant for HCN emission, resulting in line intensities which are
indistinguishable from the colliding case. We attribute this to our use of a
lower initial density compared to Priestley et al. (2023a) ($10$ versus
$250\,{\rm cm}^{-3}$), so that the initial turbulent velocity field is
dissipated before the cloud material has had time to collapse to the higher
densities where HCN emission can be produced. Establishing the relevance of
this distinction to the interpretation of observational data is left to future
work.

##  5 Conclusions

We have combined an MHD simulation of the formation and dynamical evolution of
a molecular cloud with time-dependent chemistry and radiative transfer
modelling, producing self-consistent synthetic observations in the $J=1-0$
rotational transitions of 12CO, 13CO, C18O, HCN, N2H+, HCO+ and HNC, the
$J=2-1$ transition of CS, and the $(1,1)$ inversion transition of para-NH3.
Our main results are as follows:

  * •

Our simulated cloud reproduces the observed range of intensities for most
lines well, with better agreement for species which preferentially trace high-
density regions. This suggests that despite the idealised initial conditions
and modelling assumptions, the simulation has produced a molecular cloud with
similar structural properties to real ones.

  * •

HCN, along with many other molecules often described as dense-gas tracers, is
actually a poor tracer of dense gas, being detectable down to low column
density. The velocity information derived from the HCN line profile mostly
traces the bulk dynamics of the cloud. N2H+ is a much more selective tracer of
high-density ($>10^{22}\,{\rm cm}^{-2}$) regions, and its line profile shape
is strongly correlated with the dynamics of material above $10^{4}\,{\rm
cm}^{-3}$.

  * •

Cores formed from the large-scale cloud dynamics in the simulation agree well
with the distribution of line intensities from observational samples, although
their line widths are somewhat broader. They do not agree with the predictions
of an equivalent sample of isolated core models, which neglect the formation
of the cores themselves from their more diffuse environment. Using these
isolated core models to interpret molecular line observations is likely to
result in erroneous conclusions.

  * •

We obtain CS line intensities in good agreement with observed values using an
undepleted elemental abundance of sulphur. This suggests that the common
modelling practice of reducing this quantity by several orders of magnitude to
match observations of sulphur-bearing species may be a sign of the limitations
of the physical models employed, rather than genuine evidence for sulphur
depletion in molecular clouds.

Our work demonstrates the importance of simulating the large-scale dynamics of
molecular clouds when investigating their line emission properties, even when
the subject of interest is on much smaller scales. This serves as a starting
point for future work to fully exploit the ability of molecular line emission
to advance our understanding of star formation.

## Acknowledgements

We are grateful to Rachel Friesen for providing the GAS data, and to Nicolas
Peretto for a useful discussion about the comparison with observations. FDP,
PCC, SER and OF acknowledge the support of a consolidated grant (ST/K00926/1)
from the UK Science and Technology Facilities Council (STFC). SCOG and RSK
acknowledge funding from the European Research Council (ERC) via the ERC
Synergy Grant “ECOGAL-Understanding our Galactic ecosystem: From the disk of
the Milky Way to the formation sites of stars and planets” (project ID
855130), from the Heidelberg Cluster of Excellence (EXC 2181 - 390900948)
“STRUCTURES: A unifying approach to emergent phenomena in the physical world,
mathematics, and complex data”, funded by the German Excellence Strategy, and
from the German Ministry for Economic Affairs and Climate Action in project
“MAINN” (funding ID 50OO2206). The team in Heidelberg also thanks for
computing resources provided by The Länd through bwHPC and DFG through grant
INST 35/1134-1 FUGG and for data storage at SDS@hd through grant INST
35/1314-1 FUGG. LRP acknowledges support from the Irish Research Council
Laureate programme under grant number IRCLA/2022/1165. This research was
undertaken using the supercomputing facilities at Cardiff University operated
by Advanced Research Computing at Cardiff (ARCCA) on behalf of the Cardiff
Supercomputing Facility and the Supercomputing Wales (SCW) project. We
acknowledge the support of the latter, which is part-funded by the European
Regional Development Fund (ERDF) via the Welsh Government. This research has
made use of spectroscopic and collisional data from the EMAA database
(https://emaa.osug.fr and https://dx.doi.org/10.17178/EMAA). EMAA is supported
by the Observatoire des Sciences de l’Univers de Grenoble (OSUG)

## Data Availability

The data underlying this article will be shared on request.

## References

  * Anderson et al. (2021) Anderson M., et al., 2021, MNRAS, 508, 2964
  * Barnes et al. (2020) Barnes A. T., et al., 2020, MNRAS, 497, 1972
  * Bergin & Tafalla (2007) Bergin E. A., Tafalla M., 2007, ARA&A, 45, 339
  * Bisbas et al. (2019) Bisbas T. G., Schruba A., van Dishoeck E. F., 2019, MNRAS, 485, 3097
  * Boogert et al. (2015) Boogert A. C. A., Gerakines P. A., Whittet D. C. B., 2015, ARA&A, 53, 541
  * Chira et al. (2014) Chira R.-A., Smith R. J., Klessen R. S., Stutz A. M., Shetty R., 2014, MNRAS, 444, 874
  * Clark & Glover (2014) Clark P. C., Glover S. C. O., 2014, MNRAS, 444, 2396
  * Clark et al. (2012a) Clark P. C., Glover S. C. O., Klessen R. S., 2012a, MNRAS, 420, 745
  * Clark et al. (2012b) Clark P. C., Glover S. C. O., Klessen R. S., Bonnell I. A., 2012b, MNRAS, 424, 2599
  * Clark et al. (2019) Clark P. C., Glover S. C. O., Ragan S. E., Duarte-Cabral A., 2019, MNRAS, 486, 4622
  * Clarke et al. (2018) Clarke S. D., Whitworth A. P., Spowage R. L., Duarte-Cabral A., Suri S. T., Jaffa S. E., Walch S., Clark P. C., 2018, MNRAS, 479, 1722
  * Dame & Lada (2023) Dame T. M., Lada C. J., 2023, ApJ, 944, 197
  * Denis-Alpizar et al. (2020) Denis-Alpizar O., Stoecklin T., Dutrey A., Guilloteau S., 2020, MNRAS, 497, 4276
  * Dullemond et al. (2012) Dullemond C. P., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library, record ascl:1202.015 (ascl:1202.015)
  * Dumouchel et al. (2010) Dumouchel F., Faure A., Lique F., 2010, MNRAS, 406, 2488
  * Evans et al. (2020) Evans Neal J. I., Kim K.-T., Wu J., Chao Z., Heyer M., Liu T., Nguyen-Lu’o’ng Q., Kauffmann J., 2020, ApJ, 894, 103
  * Faure et al. (2007) Faure A., Varambhia H. N., Stoecklin T., Tennyson J., 2007, MNRAS, 382, 840
  * Faustino Vieira et al. (2024) Faustino Vieira H., Duarte-Cabral A., Davis T. A., Peretto N., Smith M. W. L., Querejeta M., Colombo D., Anderson M., 2024, MNRAS, 527, 3639
  * Fehér et al. (2022) Fehér O., et al., 2022, ApJS, 258, 17
  * Fehér et al. (2024) Fehér O., Ragan S. E., Priestley F. D., Clark P. C., Moore T. J. T., 2024, MNRAS, 530, 1311
  * Friesen et al. (2017) Friesen R. K., et al., 2017, ApJ, 843, 63
  * Fuente et al. (2019) Fuente A., et al., 2019, A&A, 624, A105
  * Fuente et al. (2023) Fuente A., et al., 2023, A&A, 670, A114
  * Genel et al. (2013) Genel S., Vogelsberger M., Nelson D., Sijacki D., Springel V., Hernquist L., 2013, MNRAS, 435, 1426
  * Glover & Clark (2012) Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 116
  * Glover & Mac Low (2007) Glover S. C. O., Mac Low M.-M., 2007, ApJS, 169, 239
  * Gong et al. (2017) Gong M., Ostriker E. C., Wolfire M. G., 2017, ApJ, 843, 38
  * Goodman et al. (1998) Goodman A. A., Barranco J. A., Wilner D. J., Heyer M. H., 1998, ApJ, 504, 223
  * Habing (1968) Habing H. J., 1968, Bull. Astron. Inst. Netherlands, 19, 421
  * Hacar et al. (2024) Hacar A., et al., 2024, arXiv e-prints, p. arXiv:2403.08091
  * Hily-Blant et al. (2022) Hily-Blant P., Pineau des Forêts G., Faure A., Lique F., 2022, A&A, 658, A168
  * Holdship et al. (2017) Holdship J., Viti S., Jiménez-Serra I., Makrymallis A., Priestley F., 2017, AJ, 154, 38
  * Holdship et al. (2019) Holdship J., et al., 2019, ApJ, 878, 64
  * Hunter et al. (2023) Hunter G. H., Clark P. C., Glover S. C. O., Klessen R. S., 2023, MNRAS, 519, 4152
  * Indriolo et al. (2015) Indriolo N., et al., 2015, ApJ, 800, 40
  * Jenkins (2009) Jenkins E. B., 2009, ApJ, 700, 1299
  * Jensen et al. (2023) Jensen S. S., Spezzano S., Caselli P., Grassi T., Haugbølle T., 2023, A&A, 675, A34
  * Jiménez-Donaire et al. (2023) Jiménez-Donaire M. J., et al., 2023, A&A, 676, L11
  * Jones et al. (2023) Jones G. H., Clark P. C., Glover S. C. O., Hacar A., 2023, MNRAS, 520, 1005
  * Kauffmann et al. (2017) Kauffmann J., Goldsmith P. F., Melnick G., Tolls V., Guzman A., Menten K. M., 2017, A&A, 605, L5
  * Keto et al. (2015) Keto E., Caselli P., Rawlings J., 2015, MNRAS, 446, 3731
  * Klessen et al. (2005) Klessen R. S., Ballesteros-Paredes J., Vázquez-Semadeni E., Durán-Rojas C., 2005, ApJ, 620, 786
  * Könyves et al. (2015) Könyves V., et al., 2015, A&A, 584, A91
  * Kushwahaa et al. (2023) Kushwahaa T., Drozdovskaya M. N., Tychoniec Ł., Tabone B., 2023, A&A, 672, A122
  * Lada et al. (2010) Lada C. J., Lombardi M., Alves J. F., 2010, ApJ, 724, 687
  * Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  * Lee et al. (1999) Lee C. W., Myers P. C., Tafalla M., 1999, ApJ, 526, 788
  * Lique et al. (2006) Lique F., Spielfiedel A., Cernicharo J., 2006, A&A, 451, 1125
  * Lique et al. (2015) Lique F., Daniel F., Pagani L., Feautrier N., 2015, MNRAS, 446, 1245
  * Loreau et al. (2023) Loreau J., Faure A., Lique F., Demes S., Dagdigian P. J., 2023, MNRAS, 526, 3213
  * Magalhães et al. (2018) Magalhães V. S., Hily-Blant P., Faure A., Hernandez-Vera M., Lique F., 2018, A&A, 615, A52
  * Mathis et al. (1983) Mathis J. S., Mezger P. G., Panagia N., 1983, A&A, 128, 212
  * McClure et al. (2023) McClure M. K., et al., 2023, Nature Astronomy, 7, 431
  * McElroy et al. (2013) McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. J., 2013, A&A, 550, A36
  * Mouschovias & Spitzer (1976) Mouschovias T. C., Spitzer L. J., 1976, ApJ, 210, 326
  * Navarro-Almaida et al. (2020) Navarro-Almaida D., et al., 2020, A&A, 637, A39
  * Offner et al. (2008) Offner S. S. R., Krumholz M. R., Klein R. I., McKee C. F., 2008, AJ, 136, 404
  * Ossenkopf & Henning (1994) Ossenkopf V., Henning T., 1994, A&A, 291, 943
  * Pakmor et al. (2011) Pakmor R., Bauer A., Springel V., 2011, MNRAS, 418, 1392
  * Panessa et al. (2023) Panessa M., Seifried D., Walch S., Gaches B., Barnes A. T., Bigiel F., Neumann L., 2023, MNRAS, 523, 6138
  * Peñaloza et al. (2017) Peñaloza C. H., Clark P. C., Glover S. C. O., Shetty R., Klessen R. S., 2017, MNRAS, 465, 2277
  * Peñaloza et al. (2018) Peñaloza C. H., Clark P. C., Glover S. C. O., Klessen R. S., 2018, MNRAS, 475, 1508
  * Peretto et al. (2020) Peretto N., et al., 2020, MNRAS, 496, 3482
  * Peretto et al. (2023) Peretto N., Rigby A. J., Louvet F., Fuller G. A., Traficante A., Gaudel M., 2023, MNRAS, 525, 2935
  * Pety et al. (2017) Pety J., et al., 2017, A&A, 599, A98
  * Pineda et al. (2010) Pineda J. E., Goodman A. A., Arce H. G., Caselli P., Foster J. B., Myers P. C., Rosolowsky E. W., 2010, ApJ, 712, L116
  * Pineda et al. (2024) Pineda J. E., et al., 2024, arXiv e-prints, p. arXiv:2402.16202
  * Priestley et al. (2022) Priestley F. D., Yin C., Wurster J., 2022, MNRAS, 515, 5689
  * Priestley et al. (2023a) Priestley F. D., Clark P. C., Whitworth A. P., 2023a, MNRAS, 519, 6392
  * Priestley et al. (2023b) Priestley F. D., Arzoumanian D., Whitworth A. P., 2023b, MNRAS, 522, 3890
  * Priestley et al. (2023c) Priestley F. D., Clark P. C., Glover S. C. O., Ragan S. E., Fehér O., Prole L. R., Klessen R. S., 2023c, MNRAS, 524, 5971
  * Prole et al. (2022) Prole L. R., Clark P. C., Klessen R. S., Glover S. C. O., 2022, MNRAS, 510, 4019
  * Ragan et al. (2011) Ragan S. E., Bergin E. A., Wilner D., 2011, ApJ, 736, 163
  * Ragan et al. (2012a) Ragan S., et al., 2012a, A&A, 547, A49
  * Ragan et al. (2012b) Ragan S. E., Heitsch F., Bergin E. A., Wilner D., 2012b, ApJ, 746, 174
  * Ragan et al. (2015) Ragan S. E., Henning T., Beuther H., Linz H., Zahorecz S., 2015, A&A, 573, A119
  * Redaelli et al. (2022) Redaelli E., Chacón-Tanarro A., Caselli P., Tafalla M., Pineda J. E., Spezzano S., Sipilä O., 2022, ApJ, 941, 168
  * Redaelli et al. (2024) Redaelli E., Bovino S., Lupi A., Grassi T., Gaete-Espinoza D., Sabatini G., Caselli P., 2024, A&A, 685, A67
  * Rigby et al. (2021) Rigby A. J., et al., 2021, MNRAS, 502, 4576
  * Rigby et al. (2024) Rigby A. J., et al., 2024, MNRAS, 528, 1172
  * Santa-Maria et al. (2023) Santa-Maria M. G., et al., 2023, A&A, 679, A4
  * Schöier et al. (2005) Schöier F. L., van der Tak F. F. S., van Dishoeck E. F., Black J. H., 2005, A&A, 432, 369
  * Seifried et al. (2017) Seifried D., Sánchez-Monge Á., Suri S., Walch S., 2017, MNRAS, 467, 4467
  * Sembach et al. (2000) Sembach K. R., Howk J. C., Ryans R. S. I., Keenan F. P., 2000, ApJ, 528, 310
  * Shetty et al. (2009) Shetty R., Kauffmann J., Schnee S., Goodman A. A., Ercolano B., 2009, ApJ, 696, 2234
  * Shirley (2015) Shirley Y. L., 2015, PASP, 127, 299
  * Singh & Martin (2022) Singh A., Martin P. G., 2022, ApJ, 941, 135
  * Sipilä & Caselli (2018) Sipilä O., Caselli P., 2018, A&A, 615, A15
  * Sipilä et al. (2022) Sipilä O., Caselli P., Redaelli E., Spezzano S., 2022, A&A, 668, A131
  * Smith et al. (2012) Smith R. J., Shetty R., Stutz A. M., Klessen R. S., 2012, ApJ, 750, 64
  * Smith et al. (2013) Smith R. J., Shetty R., Beuther H., Klessen R. S., Bonnell I. A., 2013, ApJ, 771, 24
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Stuber et al. (2023) Stuber S. K., et al., 2023, A&A, 680, L20
  * Szűcs et al. (2014) Szűcs L., Glover S. C. O., Klessen R. S., 2014, MNRAS, 445, 4055
  * Tafalla et al. (2002) Tafalla M., Myers P. C., Caselli P., Walmsley C. M., Comito C., 2002, ApJ, 569, 815
  * Tafalla et al. (2021) Tafalla M., Usero A., Hacar A., 2021, A&A, 646, A97
  * Tafalla et al. (2023) Tafalla M., Usero A., Hacar A., 2023, A&A, 679, A112
  * Tress et al. (2020) Tress R. G., Smith R. J., Sormani M. C., Glover S. C. O., Klessen R. S., Mac Low M.-M., Clark P. C., 2020, MNRAS, 492, 2973
  * Tritsis et al. (2023) Tritsis A., Basu S., Federrath C., 2023, MNRAS, 521, 5087
  * Walch et al. (2012) Walch S. K., Whitworth A. P., Bisbas T., Wünsch R., Hubber D., 2012, MNRAS, 427, 625
  * Wang et al. (2024) Wang J.-W., et al., 2024, ApJ, 962, 136
  * Whitworth (1979) Whitworth A., 1979, MNRAS, 186, 59
  * Whitworth & Jaffa (2018) Whitworth A. P., Jaffa S. E., 2018, A&A, 611, A20
  * Wilson & Rood (1994) Wilson T. L., Rood R., 1994, ARA&A, 32, 191
  * Wurster & Rowan (2024) Wurster J., Rowan C., 2024, MNRAS, 528, 2257
  * Yang et al. (2010) Yang B., Stancil P. C., Balakrishnan N., Forrey R. C., 2010, ApJ, 718, 1062
  * Yin et al. (2021) Yin C., Priestley F. D., Wurster J., 2021, MNRAS, 504, 2381
  * Young et al. (2019) Young A. K., Bate M. R., Harries T. J., Acreman D. M., 2019, MNRAS, 487, 2853
  * Zhou et al. (2022) Zhou J.-X., Li G.-X., Chen B.-Q., 2022, MNRAS, 513, 638
  * Zucker et al. (2021) Zucker C., et al., 2021, ApJ, 919, 35

##  Appendix A Resolution tests

Figure 12: Column density of the radmc3d grid versus the original arepo
Voronoi mesh. Crosses show the median pixel values, error bars the 16th/84th
percentiles. The dashed red line shows the expected one-to-one relationship.

Figure 13: Integrated 12CO (left), 13CO (centre) and C18O (right) intensities
using the NEATH tracer particle abundances versus the AREPO cell values.
Crosses show the median pixel values, error bars the 16th/84th percentiles.
The dashed red lines shows the expected one-to-one relationship.

The MHD, chemical and radiative transfer simulations in this paper each have
mutually-incompatible concepts of resolution, which results in some losses
when properties are mapped between them. The cubic adaptive mesh used by
radmc3d loses detail when compared to the unstructured Voronoi mesh used by
arepo, and with only $10^{5}$ post-processed tracer particles compared to
several million grid cells, the information on molecular abundances is
inevitably ‘smeared out’ to some extent. We demonstrate below that neither of
these effects is large enough to substantially alter our results.

Figure 12 shows the relationship between the true column density obtained from
the arepo Voronoi mesh, and the column density from the radmc3d grid used for
the radiative transfer. We determine the latter by fixing the dust temperature
to $10\,{\rm K}$ in all cells, and producing an image of the cloud at
$850\,{\rm\mu m}$, where the dust optical depths are $\lesssim 10^{-4}$. As
the dust opacity and gas-to-dust ratio are both constant, the intensity is
directly proportional to the column density of the radmc3d grid. After
converting from intensity to $N_{\rm H}$, the radmc3d column densities are
almost identical to the true values, with the exception of the highest end of
the range near $10^{23}\,{\rm cm}^{-2}$. Even here, the difference is less
than a factor of two, and these column densities represent a small fraction of
the cloud area, so we consider the mapping between the arepo and radmc3d grids
to be robust.

While the radmc3d grid appears to accurately represent the underlying
structure of the MHD simulation, the molecular abundances used in the
radiative transfer are taken from $10^{5}$ post-proessed tracer particles,
compared to $\sim 20$ million cells in the radmc3d grid. We assess the
accuracy of this approach using CO emission lines, as the CO abundance for
every cell can be taken directly from the arepo simulation as a point of
comparison to the NEATH tracer particle abundance.

Figure 13 shows the relationship between integrated intensity of the three CO
isopologues for the two abundance determinations. All three lines show a
linear correlation in intensity between the NEATH and arepo values, albeit
with some scatter. The two rarer isotopologues tend to have somewhat lower
line intensities for the NEATH abundances than when using the arepo values,
but this can be attributed to the neglect of freeze-out in the arepo chemical
model, which severely reduces the CO abundance at high volume densities
(Priestley et al., 2023c). This effect is not visible for the more-abundant
12CO line, which is too optically-thick to probe the densities where freeze-
out occurs. The otherwise-good correlation for all isotopologues suggests that
the number of tracer particles is adequate to correctly capture the molecular
structure of the cloud.

##  Appendix B Dust-derived column densities

Figure 14: Left: Column density map of the cloud derived from the far-IR dust
emission (compare Figure 1). Right: True column density versus dust-derived
column density. Crosses show the median pixel values, error bars the 16th/84th
percentiles. The dashed red line shows the expected one-to-one relationship.

Figure 14 shows the column density of the cloud inferred from thermal dust
emission, the most common means of assessing molecular cloud structures
observationally. We produce images of the cloud at wavelengths of $100$,
$160$, $250$, $350$ and $500\,{\rm\mu m}$, corresponding to the central
wavelengths of the Herschel filters generally used to study cold molecular gas
(e.g. Ragan et al., 2012a; Könyves et al., 2015). We fit the resulting pixel-
by-pixel spectral energy distributions (SEDs) with a single-temperature
modified blackbody, using the same dust opacity as in the radiative transfer
calculation, and the same dust-to-gas ratio to convert dust mass into total
column density. Any differences between the two column densities are therefore
due to temperature variations along the line of sight (Shetty et al., 2009).

The dust-derived column densities are generally in good agreement with the
true values. However, this breaks down at the highest values, with peaks in
the true column density images less prominent in the dust-column maps. The
correlation between the column densities derived from dust emission and the
true values, shown in Figure 14, are nearly identical up to true columns of
$10^{22}\,{\rm cm}^{-2}$. Beyond this point, the column densities derived from
dust emission are lower than the true values by factors of up to $2-3$, as an
increasing proportion of the mass is made up of colder, less emissive grains,
whereas the emission is still dominated by the highest grain temperatures. The
luminosity-weighted temperature of the blackbody fit is therefore higher than
the mass-weighted temperature needed to retrieve the physical column density,
which is correspondingly underestimated. As we are mainly concerned with the
line emission properties of the cloud, we use the true simulation column
densities for simplicity, but note that this effect could become important for
the highest-density sightlines.