# The ramp-up of interstellar medium enrichment at z $>$ 4

M. Franco∗1 K. E. K. Coppin1 J. E. Geach1 C. Kobayashi1 S. C. Chapman2,3 C.
Yang4 E. González-Alfonso5 J. S. Spilker6 A. Cooray7 M. J. Michałowski8

###### Abstract

Fluorine is one of the most interesting elements for nuclear and stellar
astrophysics[1, 2]. Fluorine abundance was first measured for stars other than
the Sun in 1992[[1]], then for a handful metal-poor stars[3], which are likely
to have formed in the early Universe. The main production sites of fluorine
are under debate and include asymptotic giant branch (AGB) stars,
$\nu$-process in core-collapse supernovae, and Wolf-Rayet (WR) stars[4, 5, 6,
7, 8, 9, 10]. Due to the difference in the mass and lifetime of progenitor
stars, high redshift observations of fluorine can help constrain the mechanism
of fluorine production in massive galaxies. Here, we report the detection of
HF (S/N = 8) in absorption in a gravitationally lensed dusty star-forming
galaxy at redshift z=4.4 with $N_{\rm HF}$/$N_{\rm{H_{2}}}$ as high as $\sim
2\times 10^{-9}$, indicating a very quick ramp-up of the chemical enrichment
in this high-z galaxy. At z=4.4, AGB stars of a few solar masses are very
unlikely to dominate the enrichment. Instead, we show that WR stars are
required to produce the observed fluorine abundance at this time, with other
production mechanisms becoming important at later times. These observations
therefore provide an insight into the underlying processes driving the ‘ramp-
up’ phase of chemical enrichment alongside rapid stellar mass assembly in a
young massive galaxy.

Centre for Astrophysics Research, School of Physics, Engineering and Computer
Science, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK

Department of Physics and Atmospheric Science, Dalhousie University, Halifax,
NS B3H 4R2, Canada

National Research Council, Herzberg Astronomy and Astrophysics, 5071 West
Saanich Rd, Victoria, BC, V9E 2E7, Canada

European Southern Observatory, Alonso de Córdova 3107, Casilla, 19001
Vitacura, Santiago, Chile

Universidad de Alcalá, Departamento de Física y Matemáticas, Campus
Universitario, 28871 Alcalaá de Henares, Madrid, Spain

Department of Astronomy, University of Texas at Austin, 2515 Speedway, Stop
C1400, Austin, TX 78712, USA

Department of Physics and Astronomy, University of California, Irvine, CA
92697, USA

Astronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz
University, ul. Słoneczna 36, 60-286 Poznań, Poland

We observed the gravitationally lensed galaxy NGP–190387 at z=4.420 at
211.6–230.0 GHz (see Methods) with the Atacama Large Millimetre/submillimetre
Array (ALMA) seeking to detect the millimetre continuum and various molecular
and atomic emission and absorption features in this spectral range. Figure 1
shows 1.36-mm continuum contours (r.m.s = 15 $\mu$Jy beam-1) overlaid on a
K-band image. The arc-like morphology of the millimetre emission and presence
of a counter image with identical emission lines confirm without ambiguity
that NGP–190387 is lensed by a foreground structure. Lensing was first
suggested based on the spatial offset between the millimetre-wave continuum
emission and the sources detected in near-infrared imaging[11]. The Gemini
Near InfraRed Imager and spectrograph (NIRI) imaging, and data taken in the R
and I-bands with the Gemini Multi-Object Spectrograph-North (GMOS-N), reveal a
group of at least three massive galaxies, most likely members of the
foreground lensing structure. The redshift of the lensing group can be
estimated from the K-band Hubble relation[12], suggesting they lie at $z\sim
1$–$1.5$.

We detect a strong absorption profile corresponding unambiguously to the
$J=1\rightarrow 0$ transition of HF at $\nu_{\rm rest}=1232.48$ GHz,
redshifted to $\nu_{\rm obs}=227.39$ GHz at $z=4.420$ (Figure 2). Other
emission features are also present, mostly from H2O. We show the line-averaged
absorption contours in Figure 1. In Figure 2, we present the continuum-
subtracted millimetre-wave spectrum centred at HF(1–0). The continuum flux
density at the frequency of the line is $\mu S_{\nu}=7.2\pm 0.3$ mJy, where
$\mu$ is the magnification factor, which we estimate to be $\mu\approx 5$ (see
Methods). The HF(1–0) line can be described adequately by a single Gaussian
profile with $\text{FWHM}=416\pm 39$ km s-1. We measure an equivalent width of
${\rm EW}=63\pm 10$ km s-1 for the absorption line and a velocity-integrated
line flux of $\mu S\Delta V=-490\pm 80$ mJy km s-1 (uncertainties being
computed by Monte-Carlo simulations). The peak of the Gaussian is slightly
offset from the expected position of the transition for the systemic redshift
(as defined by the emission lines) by $-64\pm 21$ km s-1. Therefore, the
absorption line shows signs of asymmetry (see Figure 2). If we fit the line
with a double Gaussian, the centre of the strongest component is close to zero
velocity, $-23\pm 30$ km s-1, and a weaker component is blue-shifted by
$-277\pm 35$ km s-1. HF(1–0) is a potential tracer of molecular outflows[13,
14], and therefore we could be seeing tentative evidence of a molecular
outflow, as has been detected in a high fraction of similarly luminous objects
at high redshift[15]. However, the reduced $\chi^{2}$ of the double Gaussian
fit is not significantly better than the single Gaussian fit and could be
affected by differential lensing. For these reasons, we use the simple single
Gaussian fit in our analysis.

In Figure 3, we velocity-register and normalise other emission lines detected
during the same observation in ALMA Bands 3, 4 and 6. These include several
water transitions, CO(7–6) and the two atomic carbon transitions [C i].(1–0)
and [C i].(2–1). The emission lines all have a very similar shape reminiscent
of the double horned profile of a rotating disk, although merger dynamics
could also produce this profile. In contrast, the HF(1–0) absorption line has
a profile best described by a single Gaussian with a similar full width at
zero intensity as the emission lines. Given the similarity of the profiles of
the emission line molecular gas tracers, which span a wide range of excitation
conditions (e.g. [C i].(1–0) traces diffuse molecular gas while the high-J CO
and water lines trace denser and warmer material), at first it might seem
surprising that HF(1–0) has a different line profile. A simple explanation is
that in order to be seen in absorption, the HF-laden gas must be in the
foreground of an illuminating continuum source. If, as the emission-line
profiles hint, the galaxy is a rotating disk inclined to the line-of-sight,
then we will preferentially see absorption due to gas at the near-side of the
structure, and thus the full kinematic imprint of the disk rotation will not
be observed in the HF(1–0) absorption line.

Aside from the interest in this element from a chemical evolution perspective,
HF is potentially an excellent tracer of molecular hydrogen in the diffuse
interstellar medium (ISM)[16, 17, 18, 19] and could be effective to calibrate
other tracers over a large range of conditions. Indeed, the main channel for
HF formation is through a direct exothermic reaction of fluorine with
molecular hydrogen: H2 \+ F -¿ HF + H, with fluorine being the only element
that reacts with H2 in this way[20]. Once formed, HF is difficult to destroy
under most interstellar conditions even in low-density environments[16] and is
therefore the dominant gas-phase reservoir of fluorine[16], representing
nearly 100% of gas-phase fluorine over a wide range of temperatures and
densities[21]. Despite the relatively low abundance of fluorine, the
extraordinary stability of HF allows the accumulation of large amounts of this
molecule in the ISM. Since the first detection of HF(2–1) in absorption
towards Sagittarius B2[21] with the Infrared Space Observatory, HF(1–0) has
been observed several times in absorption and in emission in the Milky Way. HF
has also been detected mostly in absorption in a handful of other local
galaxies (see Figure 4-left panel). At higher redshifts, observations are
scarce: the molecule has only been detected quasars; in absorption in PKS
1830$-$211[[22]] at $z=0.9$, in the Cloverleaf[23] at $z=2.6$ and in emission
in BR 1202$-$0725 at $z=4.7^{[}$[14]]. This makes NGP–190387, in addition to
being the most distant detection observed in absorption, the only detection
outside of the local Universe in a non-quasar host.

We observe both HF absorption and [C i].(1–0) emission in this galaxy. While
we must take care to consider the radiative transfer effects alluded to above,
we can use this information to estimate the HF/H2 abundance, since [C i].(1–0)
offers an independent tracer of H2[[24]]. We derive $N_{\rm HF}>5\times
10^{15}$ cm-2 and N${}_{\mathrm{H_{2}}}$ = (2.1 $\pm$ 0.4) $\times 10^{24}$
cm-2 (see Methods). This implies an abundance ratio $N_{\rm HF}/$ $N_{\rm
H_{2}}>$ 2 $\times 10^{-9}$.

We compare our measurement with others from the literature, spanning Galactic
studies to the distant Universe, in Figure 4-left panel. With the caveat that
our HF abundance should be taken as a lower limit, the ratio $N_{\rm
HF}/N_{\rm H_{2}}$ in NGP–190387 may be up to one order of magnitude below
that measured for the Solar neighbourhood[25] and in diffuse clouds within the
Milky Way[17, 26, 27]. It is important to consider what internal processes
could affect the observed abundance ratio. While at low redshift and low H2
column density ($N_{\rm{H}_{2}}$ $<$ 1023 cm-2), the ratio
$N_{\rm{H}_{2}}/N_{\rm{HF}}$ seems to follow a clear trend with a low scatter
(consistent with the trend found in the “diffuse cloud abundance”, 1.8
$\times$ 10-8 [[16]]), at high redshift and high H2 column density, the trend
is less clear. Although it is difficult to draw conclusions with so few robust
detections of HF at high redshift, it is interesting to note a possible
plateau in the HF column density at high ${\rm H_{2}}$ density, with $N_{\rm
HF}\sim 5\times 10^{14}$ cm-2. NGP–190387 could be the first high redshift
example for which an abundance consistent with that found in the diffuse
clouds of the local Universe cannot be excluded. Investigations of a larger
number of high redshift galaxies at several periods of their evolution will be
needed to understand the abundance of fluorine at high H2 column density. It
has been suggested that the $N_{\rm HF}/N_{\rm H_{2}}$ ratio can decrease
sharply in dense environments[18], possibly due to accretion (including
freeze-out) of HF onto dust grains[16]. This phenomenon might be prevalent in
the high density, dusty environment of the ISM in galaxies such as
NGP–190387[[28]]. If the low ratio of $N_{\rm HF}/N_{\rm H_{2}}$ is only due
to depletion, that is, the observed relative abundance of the two components
is only due to evolution of the gas-phase, we can estimate that the fluorine
nuclei in the gas phase to HF molecules in the solid phase ratio, $R=n_{\rm
F,gas}/n_{\rm HF,solid}$[16], is less than 20% in NGP–190387. In this
scenario, the main fraction of fluorine in NGP–190387 would be locked up in
dust. However we think this level of freeze-out unlikely as the Cosmic
microwave background (CMB) temperature at the redshift of the source is
$T_{\rm CMB}\approx$ 15 K, and dust is likely to be in the vicinity of an
intense starburst (see Methods).

A simpler explanation for the relatively low HF abundance relative to H2
compared to Solar is the picture where early galaxies are gas-rich, with high
molecular gas fractions[29] and star-formation rates, but enrichment of the
ISM is still ramping up. The origin of fluorine is still poorly understood and
not well constrained, and so-far, fluorine nucleosynthesis processes have only
been observationally confirmed in AGB stars[2]. The formation of fluorine in
AGB stars during the thermally pulsing phase depends on the stars’ initial
masses and metallicity, with a production that is higher for initial masses
2–4 M⊙[[30]] with a peak at $\sim$3 M⊙[[2]]. AGB stars can produce more than
50% of the total fluorine in the Solar neighbourhood[31]. The lifetimes of 2–4
M⊙ stars are about 0.2–0.9 Gyr at low metallicities. It is expected that the
F/H ratio is much lower than the Solar, since F is produced from low-mass
stars during their AGB phase. If the star formation history of NGP–190387 were
similar to our Milky Way Galaxy, the fluorine abundance would be $\sim 1\%$ of
the Solar neighbourhood abundance at redshift of 4.4, since there is not
enough time to have the enrichment not only from AGB stars but also from
massive stars. However, this galaxy experience an episode of intense star
formation (see Methods). With a timescale of less than 0.1 Gyr, it could be
possible to reproduce the observed F/H with a rapid enrichment from massive
stars via supernovae. If the star formation timescale is longer, then it would
be necessary to have extra sources, which are WR stars, i.e., rotating massive
stars.

At $z=4.4$ the contribution to fluorine production by AGB stars is therefore
small, unless (i) the formation redshift is very high ($z>10$; when the
Universe was 450 Myr old, 900 Myr before z=4.4), (ii) a significant fraction
of stars have formed almost instantaneously, and (iii) this galaxy is already
chemically evolved with the stellar metallicity greater than $0.1\,Z_{\odot}$
(see Methods). However, this problem can be solved with the inclusion of WR
stars. In Figure 4-right panel, we show the predictions of state-of-the-art
chemical evolution models with and without WR contribution. The evolution of
the fluorine abundance, $\rm{[F/H]}\sim\log_{10}[N_{\rm HF}/N_{\rm
H_{2}}/(N_{{\rm F},\odot}/N_{{\rm H},\odot})]$, in the ISM of NGP–190387 is
calculated with so-called one-zone model[32], which includes the metal-
dependent nucleosynthesis yields of the full mass range of stars.

In the standard model with nucleosynthesis yields that are consistent with the
chemical evolution of the Milky Way, the first and second peaks of the
fluorine abundance evolution (that we can see at $z\sim 4.4$ and $z\sim 3$ in
the Figure 4-right panel) are caused by core-collapse supernovae (13–50
M⊙[[31]]) and AGB stars (1–8 M⊙[[2]]), respectively (we note that the
contribution from super-AGB stars (8–10M ⊙) is also included in the models but
they make no difference[31]). With only AGB stars, the model fluorine
abundance is too low to explain the observations. In addition to the standard
model, we also show a model including stellar rotation, which predicts a
significant contribution from WR stars even at very low metallicities. In this
case, pre-supernova and core-collapse supernova yields are taken from the
literature[33] (set “R” , for set Recommended[33]) with stellar rotation of
300 km s-1 at 13–120 M⊙, in addition to the same yields for super-AGB and AGB
yields as in the standard model, and considering that all stars with mass
greater than 25 M⊙ fully collapse to a black hole (except for pre-supernova
mass-loss). This contribution from WR stars should be considered as a maximum
as the high rotation velocity is assumed for all massive stars regardless of
metallicity; note that in observations of nearby stars, the peak of rotational
velocity is much smaller[34]. It is also important to note that the fluorine
abundance also depends on the stellar initial mass function (IMF), and a
Kroupa IMF with a slope of $x=1.3$ is adopted for 0.01–120 M⊙ in these models.
A more bottom-heavy IMF would give enhanced AGB contribution but a lower
contribution from supernovae and WR stars, resulting in lower F abundances
than in these models. The flattest IMF with the slope of $x=1.1$ gives a
$\sim$0.4 dex higher abundance, which is also not large enough to explain the
observation unless the star formation timescale is extremely short.

The predicted fluorine abundance with the inclusion of WR stars increases very
quickly to reach the observed value of NGP–190387, and stays high afterwards.
This may add further observational evidence (and the first outside the local
Universe) of the need of WR stars in the nucleosynthesis of fluorine and
further strengthens the hypothesis that WR stars (and/or $\nu$-process[8]) are
a dominant channel of ISM chemical enrichment during the key formation epoch
of massive galaxies in the early Universe. Novae can also be fluorine
producers[35] but do not contribute in our high-redshift galaxy due to the
long timescale.

##  0.1 Galaxy selection

HATLAS J133337.6+241541, hereafter NGP–190387, was discovered[36] in the
H-ATLAS survey[37] (Herschel Astrophysical Terahertz Large Area Survey), which
was undertaken with the Spectral and Photometric Imaging Receiver[38] (SPIRE)
on-board the Herschel Space Observatory[39]. NGP–190387 was identified as a
‘500-$\mu$m riser’[40], with $S_{250}$ $<$ $S_{350}$ $<$ $S_{500}$, which
identifies potential high-redshift ($z>4$) dusty star-forming galaxies
(DSFGs). The redshift ($z=4.420\pm 0.001$) and astrometric position
(13h33m37.47s, +24∘15′39.3′′ J2000) of NGP–190387 were measured unambiguously
using the NOrthern Extended Millimetre Array (NOEMA) through 12CO(5–4) and
(4–3) line emission[11]. This study detects multiple additional emission lines
at the same position and redshift.

##  0.2 Data reduction

The ALMA observations were carried out on 20 March 2018 during the ALMA Cycle
5 campaign (project 2017.1.00510.S). The observations used 46 antennas of the
12 m array, in configuration C43–4 in Band 6. The target was observed for 61.5
min (on-source integration time 42.3 min) to reach a median r.m.s. per 40-km
s-1 channel of 0.14 mJy beam-1 in Band 6. The full width half maximum (FWHM)
synthesised beam size is 0$\overset{\prime\prime}{.}$88 $\times$
0$\overset{\prime\prime}{.}$57 at position angle 29.5°. J1229+0203 and
J1327+2210 were observed as calibrators; the flux density scale should be
accurate to 5 per cent. The data were calibrated using the Common Astronomy
Software Applications (CASA) ALMA pipeline v.5.6.1-8[[41]]. We used CASA task
tclean, employing multi-scale cleaning using a geometric progression of
scales[42], with natural weighting to maximise sensitivity, masking the
central region. Continuum subtraction was performed in UV-space using the
line-free channels with the CASA task uvcontsub. New continuum subtracted
cubes have been created using the same procedure as before. Galaxy sizes
(including error estimations[43]) were obtained using the CASA task imfit in
each channel. The sizes presented in this paper are the deconvolved sizes,
uncorrected for magnification factor unless otherwise stated. For a S/N of
$\sim$3 per channel, the reliable size measurement limit for an interferometre
is approximately half the beam size[44].

##  0.3 Infrared properties

We estimate a rest-frame 8–1000-$\mu$m infrared luminosity, uncorrected for
gravitational amplification, of ${\log_{10}}(\mu L_{\rm IR}/L_{\odot})=13.5\pm
0.2$ from the SED-fitting code CIGALE[45]. We measure the dust temperature by
fitting a modified blackbody model with emissivity $\beta=1.5$[[46, 47]],
$T_{\rm dust}=42\pm 2$ K. Even for the highest plausible magnifications, this
suggests a prodigious star-formation rate, $\text{SFR}\approx 3500/\mu$ M⊙
yr-1, for a Kroupa initial mass function[48, 49].

##  0.4 Details on the computation of N${}_{\rm{H}_{2}}$

To determine $N_{\rm H_{2}}$ we use the [C i].(1–0) emission line and assumed
a C i/H2 conversion factor $X_{\mathrm{C\,{\sc I}}}$:

$N_{\mathrm{H}_{2}}~{}[\mathrm{cm^{-2}}]=X_{\mathrm{CI}}~{}\mathrm{[cm^{-2}~{}(K~{}km\,s^{-1})^{-1}]}\times
W_{\mathrm{CI}}~{}[\mathrm{K~{}km\,s^{-1}}]$ (1)

where $N_{\mathrm{H}_{2}}$ is the $\mathrm{H}_{2}$ column density and
$W_{\mathrm{CI}}$ is the integrated brightness temperature of the [C i].(1–0)
transition. Integrating over the source we measure $W_{\mathrm{CI}}=(6.4\pm
1.2)\times 10^{3}$ K km s-1 (the uncertainties were calculated performing 10
000 Monte-Carlo simulations perturbing the measurement randomly within the
uncertainties). We assume $X_{\mathrm{CI}}=$ $(3.3\pm 0.2)\times 10^{20}$ cm-2
(K km s-1)-1 from the ratio between the H2 mass and the line luminosity
L${}^{\prime}_{\mathrm{C\,{\sc I}(1-0)}}$:

$X_{\mathrm{CI}}~{}\mathrm{[cm^{-2}~{}(K~{}km\,s^{-1})^{-1}]}=\frac{\rm
M_{\rm{H2}}~{}\mathrm{[M_{\odot}]}}{\rm L^{\prime}_{\mathrm{C\,{\sc
I}(1-0)}}~{}\mathrm{[K\ km\ s^{-1}\ pc^{2}]}}\frac{1}{1.6\times
10^{-20}~{}\mathrm{[M_{\odot}\ pc^{-2}\ cm^{2}]}}$ (2)

with $\mu$L${}^{\prime}_{\mathrm{C\,{\sc I}(1-0)}}=(6.0\pm 0.3)\times 10^{10}$
K km s-1 pc2, $\mu$M${}_{\rm{H2}}=(3.2\pm 0.2)\times 10^{11}$ M⊙, using[50,
51]:

$M_{\rm H_{2}}\
[M_{\odot}]=1375.8\frac{d_{L}^{2}}{1+z}\left(\frac{X^{\prime}_{\mathrm{C\,{\sc
I}(1-0)}}}{10^{-5}}\right)^{-1}\left(\frac{A_{\mathrm{10}}}{10^{-7}}\right)^{-1}Q_{10}^{-1}\
S_{\mathrm{C\,{\sc I}}}\Delta v\ [\mathrm{Jy\ km\ s^{-1}}],$ (3)

with the Einstein coefficient A10 = 7.93$\times$10-8s-1 , the [C i]-to-H2
abundance ratio $X^{\prime}_{\mathrm{C\,{\sc I}}}$ = 3$\times 10^{-5}$ [[52]]
and Q10 = 0.6[[53]] the excitation factor. This result has then been divided
by 1.36 to remove the contribution from helium. This ratio is consistent with
that found for local ULIRGs[54] giving N${}_{\mathrm{H_{2}}}$ = $(2.1\pm
0.4)\times 10^{24}$ cm-2. Although for nearby galaxies, atomic hydrogen
dominates the total mass of gas in galaxies, this proportion is reversed for
more distant galaxies due to a strong increase in molecular hydrogen[55, 56,
29] and we adopt $M_{\mathrm{H_{2}}}+M_{\mathrm{H}}\sim M_{\mathrm{H_{2}}}$,
valid at $z>0.4$[[29]].

The [C i].(1–0) surface brightness of the galaxy ($W_{\mathrm{CI}}$) is given
by[57]:

$W_{\mathrm{CI}}\ [\mathrm{K~{}km\,s^{-1}}]=\frac{(1+z)\ c^{2}}{2k\ \nu_{\rm
obs}^{2}\ \Omega_{s}}S_{\mathrm{C\,{\sc I}}}\Delta v=1.38\times
10^{6}\frac{(1+z)\ S_{\mathrm{C\,{\sc I}}}\Delta v\ [\mathrm{Jy\ km\
s^{-1}}]}{\nu_{\rm obs}^{2}[\mathrm{GHz}]\ \Omega_{s}\ [\mathrm{arcsec^{2}}]}$
(4)

where $\nu$ is the observed frequency, $k$ the Boltzmann constant,
$\Omega_{s}$ the source solid angle ($\Omega_{s}$ = $\pi\theta_{\rm
min}\theta_{\rm max}$/ 4$\ln$(2)), and $S_{\mathrm{C\,{\sc I}}}$ $\Delta v$
the velocity integrated flux density. We corrected the observed [C i].(1–0)
flux for the CMB[58] by assuming local thermodynamic equilibrium (LTE) between
the gas and the dust. This represents 24$\pm$2% of the flux that had to be
corrected (added) in the column density calculation above. The size of the [C
i].(1–0) emission line is too uncertain in Band 3 to be measured accurately,
the beam size in Band 3 is 1.39” $\times$ 0.95”. We took advantage of the
higher resolution in Band 6, to estimate the size of the emission line. We
assume that the CO(10-9) emission line traces the same material as [C
i].(1–0)[[59, 60, 24]] (both lines show similar profiles) and its size can
therefore be used to calculate the H2 column density. In addition, because HF
has the same spatial distribution as CO(10-9) (see Details on the computation
of NHF), they are both reasonable tracers of the spatial distribution of the
gas responsible for the emission/absorption. The CO(10-9) emission line being
partially blended with the emission line of water H2O(312-221) (see Figure
2-left panel), and since there is no significant variation in size between
these two components, we have considered the average size of these blended
elements. We measure a source size of $(0.70^{\prime\prime}\pm
0.05^{\prime\prime})\times(0.28^{\prime\prime}\pm 0.04^{\prime\prime})$, the
uncertainties corresponding to the standard deviation of the sample means. For
ULIRGs, there may be a fine structure deficit. For this reason, we also
calculated N${}_{\rm H_{2}}$ from the continuum and water emission lines of
NGP–190387. The results from these different techniques give comparable
results.

##  0.5 Details on the computation of NHF

For extragalactic sources the HF absorption is generally expected to be
physically associated with the source(s) of far-infrared continuum, unless the
galaxy is seen edge-on or there is a foreground component in the direction of
the continuum source. We note that in the case of NGP–190387, the [C i].(1–0)
profile that is used to estimate N${}_{\rm H_{2}}$ mimics the profiles of the
excited CO and H2O lines (Figure 3), so that there is no evidence for detached
foreground low-excitation gas. The HF(1-0) absorption profile is slightly
blueshifted ($\sim$ 60 km s-1) relative to the emission [C i], CO and H2O
profiles, but this is a rather common characteristic of local ULIRGs where the
excited OH 65$\mu$m absorption, radiatively pumped by the nuclear far-IR
field, is also slightly blueshifted relative to the [C ii] 158$\mu$m emission
line[61]. Furthermore, the velocity coverage of HF(1-0) is the same as that of
the CO and H2O lines.

If the HF is physically associated with the continuum, we should then consider
the excitation effects of the infrared radiation on HF(1–0), and also the fact
that only the HF in front of the continuum will produce absorption. In order
to consider these effects and quantify their impact on the estimation of NHF,
we have modeled the source using the non-LTE, non-local approach[62], where
dust and molecules are assumed to be mixed. Any geometrical effect is ignored
by considering a spherically symmetric source; even if the actual geometry is
a disk, it can be considered in a first approach as an ensemble of these
spherical clouds, and EWHF(1-0) would be the same for the individual clouds
and for the ensemble.

By taking the following parameters in the model: a dust temperature of 50 K,
$\tau_{100}$ = 2, corresponding to N${}_{\rm H_{2}}\sim 1.3\times$1024 cm-2
(consistent with the value found by the [C i].(1–0) line) for gas-to-dust
ratio of 100, Vturb = 240 km s-1 and by excluding the effects of the
collisions, we tested two different values of the HF column density: NHF =
5$\times$1014 cm-2 and NHF = 5$\times$1015 cm-2. The absorption produced in
the NHF = 5$\times$1014 cm-2 model is very weak, with only EWHF(1-0) = 8 km
s-1. In contrast, the NHF = 5$\times 10^{15}$ cm-2 model gives EWHF(1-0) = 66
km s-1, in agreement with the observations. Since a lower value of Vturb (i.e.
a profile broadened by systemic motions rather than by turbulence alone) would
increase the required NHF, the HF column density of $5\times 10^{15}$ cm-2
would be still considered a lower limit; we conservatively estimate an
abundance of $2\times 10^{-9}$ with 0.4 dex of uncertainty (Figure 4).

We also measured the deconvolved size of the HF(1–0) absorption using the CASA
task imfit after subtraction of the continuum. In order to be consistent with
the emission line size measurements, we multiplied the flux density by -1.
Because of the weaker signal of the HF(1–0) absorption line compared to the
emission lines (see Figure 2-left panel), the size measurement could only be
made on fewer channels and may suffer from more uncertainty. We measure a
source size of $(0.59^{\prime\prime}\pm
0.06^{\prime\prime})\times(0.19^{\prime\prime}\pm 0.07^{\prime\prime})$, the
uncertainties corresponding to the standard deviation of the individual sizes
across different channels. The measured size of the HF(1–0) absorption line is
slightly smaller than the emission line region we consider for the gas.
However, as the sizes are consistent within uncertainties, we consider that
the assumption of considering that the H2 gas probes the same region as HF is
justified.

##  0.6 Chemical evolution models

The fluorine evolution depends on the star-formation history (SFH), in
addition to the enrichment source and IMF as already discussed. This galaxy is
extremely faint/not detected across the UV-optical-NIR, and therefore we have
only limited constraints on the SFH from stellar population modelling[45]. We
therefore investigated and modelled all possible SFHs varying i) star
formation efficiency (SFR/Mgas), ii) duration of star burst, iii) the fraction
of additional, older stellar population, iv) the formation epoch of the old
stellar population, v) gas accretion timescale, and vi) the IMF. Based on its
observed properties (infrared luminosity, gas mass, SF timescale), this galaxy
is similar to submillimetre galaxies (SMGs). We therefore found that iii),
iv), and v) are not important for chemical evolution of galaxies. There is no
observational indication that the IMF of this galaxy should be different from
that in the Milky Way, and thus we present the models with the Kroupa IMF in
Figure 4-right panel. We then find that the observed F/H ratio of NGP–190387
can be reproduced (A) if the star formation timescale is shorter than 0.3 Gyr,
or (B) with significant F production from WR stars.

In Figure 4-right panel, we present the models with star formation timescales
from 0.1 to 0.5 Gyr, which is a relative constant between SFR and gas fraction
at a give time (also called depletion timescale)[32]. The gas fraction
increases due to the accretion of pristine gas with an exponential timescale
of 1 Gyr in these models, and almost identical results are obtained with a
longer timescale (e.g., 5 Gyr); the accretion timescale would have an impact
for a much longer star formation timescale such as in our Milky Galaxy. The
high SFR of this galaxy requires some mechanism of a starburst, such as galaxy
merger, and in that case, the duration of starburst is expected to be as short
as 0.1 Gyr. We also assume the formation epoch to be the epoch of
reionization, $z=$ 7.7[[63]], which is most plausible; an earlier formation
epoch (e.g., $z=10$) results in a slightly more old stellar population, and
hence in a slightly higher F abundance at the onset of the starburst. However,
this contribution from old stellar populations is not important because the F
abundance rapidly increases during the starburst as seen in the Figure.
Similarly, increasing the SFR before the starburst (by a factor of 10) does
not affect either. Finally, if the duration of starburst is longer than in the
models, the predicted gas fraction becomes too small for this galaxy at this
redshift.

Our chemical evolution model includes the production of F from all plausible
sources, i.e., core-collapse supernovae, AGB stars (including super-AGB
stars), and WR stars using the most recent nucleosynthesis yields[31], self-
consistently. There is an uncertainty of the F yields due to the nuclear
reaction rates and the description of stellar mass loss, rotation, and
rotation-induced mixing, which are calibrated by the observations of
stars[33]. 19F is synthesized from 14N in He convective shell but can be
destroyed by the reaction 19F($\alpha$, $p$)22Ne, which is uncertain. The
Wolf-Rayet stars might not produce much 19F [[64]], but noted that their mass-
loss rates, and hence F yields, might be underestimated. It is also important
to note that a significant production of F from WR stars is required from the
Galactic chemical evolution models[31, 9, 65] and also the observational of
stars[66] in the Milky Way.

##  0.7 Determination of the magnification factor

The determination of the magnification factor has been done in the UV plane
using the VISILENS[67, 68] package. A detailed description of the lens
modeling will be presented in a companion paper. Using continuum emission,
having first subtracted the emission lines in Band 6 and taking the hypothesis
of a single lens we obtain a magnification factor $\mu=4.9\pm 0.2$.

##  0.8 Freeze-out of the fluorine on dust grains

One hypothesis to explain the low abundance of hydrogen fluoride in NGP–190387
compared to the abundance in the vicinity of the sun is that HF can condense
onto dust grains. We have investigated this possibility by considering that
the HF abundance follows from balancing its adsorption onto dust grains and
desorption[69, 70, 71]. The characteristic time for the HF(1–0) to freeze-out
($\tau_{\mathrm{freeze-out}}$) is defined by[69]:

$\tau_{\mathrm{freeze-out}}(n_{\mathrm{H}},T_{\mathrm{gas}})\ [yr]=6.97\times
10^{9}\left(\frac{m_{\mathrm{HF}}}{T_{\mathrm{gas}}}\right)^{0.5}n_{\mathrm{H}}^{-1},$
(5)

where $T_{\mathrm{gas}}$ is the temperature of the gas, $m_{\mathrm{HF}}$ is
the molecular weight of HF, and $n_{\mathrm{H}}$ is the density of hydrogen
nuclei in m-3. The characteristic desorption time
($\tau_{\mathrm{desorption}}$) is more difficult to estimate because it
depends on the type of grain and is expressed as[69]:

$\tau_{\mathrm{desorption}}(T_{\text{dust }})\ [yr]=\frac{3.17\times
10^{-21}}{\exp\left(-\frac{E_{\mathrm{b},\mathrm{HF}}}{kT_{\text{dust}}}\right)},$
(6)

where $T_{\mathrm{dust}}$ is the temperature of the dust, $k$ is the Boltzmann
constant, and $E_{\mathrm{b},\mathrm{HF}}$ is the binding energy of HF to the
the dust grain surface. HF(1–0) may condense on to the dust grains if
$\tau_{\mathrm{freeze-out}}<\tau_{\mathrm{desorption}}$. By taking the most
conservative values of the parameters possible:
$T_{\mathrm{gas}}\in[35-100K]$, $T_{\mathrm{dust}}\in[35-50K]$ and
$n_{\mathrm{H}}\in[10^{3}-10^{5}$ cm-3], and different values of
$E_{\mathrm{b},\mathrm{HF}}$ corresponding to Hydrogenated crystaline silica,
CO ice on amorphous silica, and CO2 ice on amorphous silica[69], the freeze-
out phenomenon is negligible (the transition is highly dependent on the
temperature and occurs at $\sim$ 20K). On the other hand, one possibility
could lead to an adsorption of HF on the grains, if they are covered by a pure
water layer. This possibility has been studied in the case of NGC 6334 I and
I(N) and has been rejected[69]. While being well aware that the conditions of
the interstellar medium of NGP–190387 are not homogeneous and may suffer from
large local disparities, our conclusion is that freeze-out does not explain
the low HF abundance relative to the Solar abundance.

##  0.9 Data Availability

This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2017.1.00510.S, archived at https://almascience.nrao.edu/alma-
data/archive. Tabulated spectral data used in this study is provided at
https://github.com/maximilienfranco/f21_fluorine_spectrum

##  0.10 Code Availability

The ALMA data are processed using the CASA ALMA pipeline (v.5.6.1-8) available
at https://almascience.nrao.edu/processing/science-pipeline. The lens model
was done using the VISILENS package publicly available at
https://github.com/jspilker/visilens/

##  0.11 Acknowledgements

M.F. is grateful to Yashar Hezaveh for his advice on the lens model. M.F. and
K.E.K.C. acknowledge support from the UK Science and Technology Facilities
Council (STFC) (grant number ST/R000905/1). K.E.K.C. acknowledges support from
a Royal Society Leverhulme Trust Senior Research Fellowship (grant number RSLT
SRF/R1/191013). J.E.G. acknowledges support from a Royal Society University
Research Fellowship. C.K. acknowledges funding from the UK STFC through grant
ST/M000958/1 & ST/ R000905/1. S.C.C. acknowledges the Natural Sciences and
Engineering Research Council of Canada (NSERC). C.Y. acknowledges support from
an ESO Fellowship. J.S. is a NHFP Hubble Fellow supported by NASA Hubble
Fellowship grant no. HF2-51446 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, for NASA, under contract NAS5-26555. M.J.M. acknowledges the
support of the National Science Centre, Poland through the SONATA BIS grant
2018/30/E/ST9/00208. This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2017.1.00510.S. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ.

M.F. reduced and analysed the data; M.F., K.E.K.C., J.E.G., C.K. interpreted
the results and wrote the paper. C.K. created the chemical evolution models.
C.Y. interpreted the results. J.S. helped create the lens model. E.G.-A.
computed the HF column density and contributed to various aspects of the
analysis. S.C.C. provided the data. All other authors contributed to the ALMA
proposals, to the scientific discussion, and provided comments to the
manuscript.

The authors declare that they have no competing interests.

Correspondence and requests for materials should be addressed to M. Franco
(email: m.franco@herts.ac.uk).

Figure 1: ALMA contours of NGP–190387. The contours (4 to 16$\sigma$ in steps
of 2$\sigma$) show the 1.36-mm continuum (cyan) and HF absorption line
(yellow). The background image is a K-band image from Gemini/NIRI. The ALMA
synthesised beam is shown in the lower left-hand corner. Figure 2: ALMA
spectrum of NGP–190387. a: Spectrum in Band 6 without continuum subtraction.
The position of emission and absorption lines are indicated by vertical lines
and labelled. For comparison, the Arp 220 spectrum[72], shifted and scaled by
arbitrary units is also displayed in gray. The frequency gap is due to the
separation of the ALMA spectral windows. By inspection, the plot of the flux
density without continuum subtraction shows that the absorption line is not
saturated. By integrating the S/N over the absorption line, we obtain an S/N =
8. b: Observed HF(1-0) line profile ($\nu_{\rm rest}=1232.48$ GHz[[73]]),
continuum subtracted. The spectral resolution is 40 km s-1, with the velocity
scale centered on the HF(1–0) transition at $z=4.420$. The single and double
component Gaussian fits are displayed in black and red respectively, with
dashed lines indicating the individual components of the double Gaussian fit.
The error bars show the uncertainties on the flux measurement given by imfit.
The bottom panel shows the residuals (data - model) in units of mJy in black
for the single Gaussian fit and in red for the double Gaussian fit. In theses
two panels, the flux density is not corrected by the magnification factor
($\mu\approx 5$). Figure 3: Overlay of the different ALMA lines. We displayed
the different lines present in the three ALMA bands, omitting CO(10–9) and
H2O(312-221) due to blending and H2O+ lines because of the low signal-to-
noise. The pattern observed for the HF line at a velocity $>$ 500 km s-1
corresponds to the line H2O(220-211) (see Figure 2, left panel). For
comparison, we also display the HF(1–0) absorption line, normalised and
multiplied by (-1) with a dashed black line. Figure 4: Evolution of fluorine
abundance. a: HF column density plotted versus H2 column density. We compare
our observations to literature data [28, 17, 74, 75, 72, 76, 23, 18, 77, 78,
27, 26, 13, 22, 69, 79, 80, 19, 14] of HF in absorption and emission (red
squares) color-coded as a function of redshift. Upper/Lower limits (3$\sigma$)
not associated with a data point imply that the HF line has not been detected,
while a limit associated with a data point indicates that the method used by
the authors may over/underestimate the abundance of the chemical element. b:
Top-panel: Contribution of SNe (blue), AGB (orange) and WR stars (green) to
the total amount of fluorine (SNe + AGB + WR) predicted by the different
chemical evolution models: timescale ($\tau$) = 0.1 Gyr (dotted line), 0.3 Gyr
(dashed line), 0.5 Gyr (solid line). The redshift of NGP–190387 is indicated
by a vertical dashed gray line. Bottom-panel: Fluorine abundance of NGP–190387
put into perspective with different chemical evolution models. For the
different models the blue, orange, and green lines take into account the
production of fluorine by the SNe only, SNe and AGB stars, and SNe, AGB and WR
stars, respectively. With a timescale longer than 0.1 Gyr, it would be
necessary to have WR stars to reproduce the measured fluorine abundance. This
shows that Fluorine is already present in the interstellar medium at very high
redshift and WR stars could be the main Fluorine producer during the key
formation epoch of massive galaxies in the early Universe. In both panels, the
Solar neighbourhood number ratio of 3.63 $\times$10-8 [[25]] is indicated with
a dashed line and NGP–190387 is displayed as a large circle outlined in black.

## References

  * [1] Jorissen, A., Smith, V. V. & Lambert, D. L.  Fluorine in red giant stars : evidence for nucleosynthesis.  _Astron. Astrophys._ 261, 164–187 (1992).
  * [2] Kobayashi, C., Karakas, A. I. & Umeda, H.  The evolution of isotope ratios in the Milky Way Galaxy.  _Mon. Not. R. Astron. Soc._ 414, 3231–3250 (2011).
  * [3] Cunha, K., Smith, V. V., Lambert, D. L. & Hinkle, K. H.  Fluorine Abundances in the Large Magellanic Cloud and $\omega$ Centauri: Evidence for Neutrino Nucleosynthesis?  _Astron. J._ 126, 1305–1311 (2003).
  * [4] Forestini, M., Goriely, S., Jorissen, A. & Arnould, M.  Fluorine production in thermal pulses on the asymptotic giant branch.  _Astron. Astrophys._ 261, 157–163 (1992).
  * [5] Woosley, S. E. & Haxton, W. C.  Supernova neutrinos, neutral currents and the origin of fluorine.  _Nature_ 334, 45–47 (1988).
  * [6] Meynet, G. & Arnould, M.  Synthesis of 19F in Wolf-Rayet stars.  _Astron. Astrophys._ 355, 176–180 (2000).
  * [7] Renda, A. _et al._ On the origin of fluorine in the Milky Way.  _Mon. Not. R. Astron. Soc._ 354, 575–580 (2004).
  * [8] Kobayashi, K. _et al._ Evolution of Fluorine in the Galaxy with the $\nu$-process.  _Astrophys. Journal Letter_ 739, L57 (2011).
  * [9] Spitoni, E., Matteucci, F., Jönsson, H., Ryde, N. & Romano, D.  Fluorine in the solar neighborhood: Chemical evolution models.  _Astron. Astrophys._ 612, A16 (2018).
  * [10] Grisoni, V. _et al._ Fluorine in the solar neighbourhood: modelling the Galactic thick and thin discs.  _Mon. Not. R. Astron. Soc._ 498, 1252–1258 (2020).
  * [11] Fudamoto, Y. _et al._ The most distant, luminous, dusty star-forming galaxies: redshifts from NOEMA and ALMA spectral scans.  _Mon. Not. R. Astron. Soc._ 472, 2028–2041 (2017).
  * [12] Rocca-Volmerange, B., Le Borgne, D., De Breuck, C., Fioc, M. & Moy, E.  The radio galaxy K-z relation: The 1012 M⊙ mass limit. Masses of galaxies from the LK luminosity, up to z $>$ 4\.  _Astron. Astrophys._ 415, 931–940 (2004).
  * [13] Monje, R. R. _et al._ Hydrogen Fluoride toward Luminous Nearby Galaxies: NGC 253 and NGC 4945.  _Astrophys. J._ 785, 22 (2014).
  * [14] Lehnert, M. D. _et al._ Etching glass in the early Universe: Luminous HF and H2O emission in a QSO-SMG pair at z = 4.7.  _Astron. Astrophys._ 641, A124 (2020).
  * [15] Spilker, J. S. _et al._ Fast molecular outflow from a dusty star-forming galaxy in the early Universe.  _Science_ 361, 1016–1019 (2018).
  * [16] Neufeld, D. A., Wolfire, M. G. & Schilke, P.  The Chemistry of Fluorine-bearing Molecules in Diffuse and Dense Interstellar Gas Clouds.  _Astrophys. J._ 628, 260–274 (2005).
  * [17] Sonnentrucker, P. _et al._ Detection of hydrogen fluoride absorption in diffuse molecular clouds with Herschel/HIFI: an ubiquitous tracer of molecular gas.  _Astron. Astrophys._ 521, L12 (2010).
  * [18] Emprechtinger, M. _et al._ Hydrogen Fluoride in High-mass Star-forming Regions.  _Astrophys. J._ 756, 136 (2012).
  * [19] Kavak, Ü., van der Tak, F. F. S., Tielens, A. G. G. M. & Shipman, R. F.  Origin of hydrogen fluoride emission in the Orion Bar. An excellent tracer for CO-dark H2 gas clouds.  _Astron. Astrophys._ 631, A117 (2019).
  * [20] Gerin, M., Neufeld, D. A. & Goicoechea, J. R.  Interstellar Hydrides.  _Annu. Rev. Astron. Astrophys._ 54, 181–225 (2016).
  * [21] Neufeld, D. A., Zmuidzinas, J., Schilke, P. & Phillips, T. G.  Discovery of Interstellar Hydrogen Fluoride 1.  _Astrophys. J._ 488, L141–L144 (1997).
  * [22] Kawaguchi, K. _et al._ Detection of HF Toward PKS 1830-211, Search for Interstellar H2F+, and Laboratory Study of H2F+ and H2Cl+ Dissociative Recombination.  _Astrophys. J._ 822, 115 (2016).
  * [23] Monje, R. R. _et al._ Discovery of Hydrogen Fluoride in the Cloverleaf Quasar at z = 2.56.  _Astrophys. J. Lett._ 742, L21 (2011).
  * [24] Papadopoulos, P. P., Thi, W. F. & Viti, S.  CI lines as tracers of molecular gas, and their prospects at high redshifts.  _Mon. Not. R. Astron. Soc._ 351, 147–160 (2004).
  * [25] Asplund, M., Grevesse, N., Sauval, A. J. & Scott, P.  The Chemical Composition of the Sun.  _Annu. Rev. Astron. Astrophys._ 47, 481–522 (2009).
  * [26] Indriolo, N., Neufeld, D. A., Seifahrt, A. & Richter, M. J.  Direct Determination of the HF/H2 Abundance Ratio in Interstellar Gas.  _Astrophys. J._ 764, 188 (2013).
  * [27] Pereira-Santaella, M. _et al._ Herschel/SPIRE Submillimeter Spectra of Local Active Galaxies.  _Astrophys. J._ 768, 55 (2013).
  * [28] Neufeld, D. A. _et al._ Strong absorption by interstellar hydrogen fluoride: Herschel/HIFI observations of the sight-line to G10.6-0.4 (W31C).  _Astron. Astrophys._ 518, L108 (2010).
  * [29] Tacconi, L. J. _et al._ PHIBSS: Unified Scaling Relations of Gas Depletion Time and Molecular Gas Fractions.  _Astrophys. J._ 853, 179 (2018).
  * [30] Lugaro, M. _et al._ Reaction Rate Uncertainties and the Production of 19F in Asymptotic Giant Branch Stars.  _Astrophys. J._ 615, 934–946 (2004).
  * [31] Kobayashi, C., Karakas, A. I. & Lugaro, M.  The Origin of Elements from Carbon to Uranium.  _Astrophys. J._ 900, 179 (2020).
  * [32] Kobayashi, C., Tsujimoto, T. & Nomoto, K.  The History of the Cosmic Supernova Rate Derived from the Evolution of the Host Galaxies.  _Astrophys. J._ 539, 26–38 (2000).
  * [33] Limongi, M. & Chieffi, A.  Presupernova Evolution and Explosive Nucleosynthesis of Rotating Massive Stars in the Metallicity Range -3 $\leq$ [Fe/H] $\leq$ 0\.  _Astrophys. J. Supp. Series_ 237, 13 (2018).
  * [34] Ramírez-Agudelo, O. H. _et al._ The VLT-FLAMES Tarantula Survey. XII. Rotational velocities of the single O-type stars.  _Astron. Astrophys._ 560, A29 (2013).
  * [35] Wiescher, M., Gorres, J., Thielemann, F. K. & Ritter, H.  Explosive hydrogen burning in novae.  _Astron. Astrophys._ 160, 56–72 (1986).
  * [36] Ivison, R. J. _et al._ The Space Density of Luminous Dusty Star-forming Galaxies at z > 4: SCUBA-2 and LABOCA Imaging of Ultrared Galaxies from Herschel-ATLAS.  _Astrophys. J._ 832, 78 (2016).
  * [37] Eales, S. _et al._ The Herschel ATLAS.  _Publ. Astron. Soc. Pacific_ 122, 499 (2010).
  * [38] Griffin, M. J. _et al._ The Herschel-SPIRE instrument and its in-flight performance.  _Astron. Astrophys._ 518, L3 (2010).
  * [39] Pilbratt, G. L. _et al._ Herschel Space Observatory. An ESA facility for far-infrared and submillimetre astronomy.  _Astron. Astrophys._ 518, L1 (2010).
  * [40] Cox, P. _et al._ Gas and Dust in a Submillimeter Galaxy at z = 4.24 from the Herschel Atlas.  _Astrophys. J._ 740, 63 (2011).
  * [41] McMullin, J. P., Waters, B., Schiebel, D., Young, W. & Golap, K.  CASA Architecture and Applications.  In Shaw, R. A., Hill, F. & Bell, D. J. (eds.) _Astronomical Data Analysis Software and Systems XVI_ , vol. 376 of _Astron. Soc. Pac. Conf. Series_ , 127 (2007).
  * [42] Cornwell, T. J.  Multiscale CLEAN Deconvolution of Radio Synthesis Images.  _IEEE J. of Selected Topics in Signal Proc._ 2, 793–801 (2008).
  * [43] Condon, J. J.  Errors in Elliptical Gaussian FITS.  _Publ. Astron. Soc. Pacific_ 109, 166–172 (1997).
  * [44] Martí-Vidal, I., Pérez-Torres, M. A. & Lobanov, A. P.  Over-resolution of compact sources in interferometric observations.  _Astron. Astrophys._ 541, A135 (2012).
  * [45] Boquien, M. _et al._ CIGALE: a python Code Investigating GALaxy Emission.  _Astron. Astrophys._ 622, A103 (2019).
  * [46] Kovács, A. _et al._ SHARC-2 350 $\mu$m Observations of Distant Submillimeter-selected Galaxies.  _Astrophys. J._ 650, 592–603 (2006).
  * [47] Gordon, K. D. _et al._ Determining dust temperatures and masses in the Herschel era: The importance of observations longward of 200 micron.  _Astron. Astrophys._ 518, L89 (2010).
  * [48] Kennicutt, J., Robert C.  The Global Schmidt Law in Star-forming Galaxies.  _Astrophys. J._ 498, 541–552 (1998).
  * [49] Kroupa, P.  On the variation of the initial mass function.  _Mon. Not. R. Astron. Soc._ 322, 231–246 (2001).
  * [50] Papadopoulos, P. P. & Greve, T. R.  C I Emission in Ultraluminous Infrared Galaxies as a Molecular Gas Mass Tracer.  _Astrophys. J._ 615, L29–L32 (2004).
  * [51] Wagg, J., Wilner, D. J., Neri, R., Downes, D. & Wiklind, T.  Atomic Carbon in APM 08279+5255 at z = 3.91.  _Astrophys. J._ 651, 46–50 (2006).
  * [52] Weiß, A., Henkel, C., Downes, D. & Walter, F.  Gas and dust in the Cloverleaf quasar at redshift 2.5.  _Astron. Astrophys._ 409, L41–L45 (2003).
  * [53] Bothwell, M. S. _et al._ ALMA observations of atomic carbon in z $\sim$ 4 dusty star-forming galaxies.  _Mon. Not. R. Astron. Soc._ 466, 2825–2841 (2017).
  * [54] Jiao, Q. _et al._ Neutral Carbon Emission in Luminous Infrared Galaxies: The [C I] Lines as Total Molecular Gas Tracers.  _Astrophys. J. Lett._ 840, L18 (2017).
  * [55] Daddi, E. _et al._ Different Star Formation Laws for Disks Versus Starbursts at Low and High Redshifts.  _Astrophys. J. Lett._ 714, L118–L122 (2010).
  * [56] Lagos, C. d. P. _et al._ Molecular hydrogen abundances of galaxies in the EAGLE simulations.  _Mon. Not. R. Astron. Soc._ 452, 3815–3837 (2015).
  * [57] Solomon, P. M. & Vanden Bout, P. A.  Molecular Gas at High Redshift.  _Annu. Rev. Astron. Astrophys._ 43, 677–725 (2005).
  * [58] da Cunha, E. _et al._ On the Effect of the Cosmic Microwave Background in High-redshift (Sub-)millimeter Observations.  _Astrophys. J._ 766, 13 (2013).
  * [59] Ojha, R. _et al._ AST/RO Observations of Atomic Carbon near the Galactic Center.  _Astrophys. J._ 548, 253–257 (2001).
  * [60] Ikeda, M., Oka, T., Tatematsu, K., Sekimoto, Y. & Yamamoto, S.  The Distribution of Atomic Carbon in the Orion Giant Molecular Cloud 1.  _Astrophys. J. Supp. Series_ 139, 467–485 (2002).
  * [61] González-Alfonso, E. _et al._ Molecular Outflows in Local ULIRGs: Energetics from Multitransition OH Analysis.  _Astrophys. J._ 836, 11 (2017).
  * [62] González-Alfonso, E., Fischer, J., Aalto, S. & Falstad, N.  Modeling the H2O submillimeter emission in extragalactic sources.  _Astron. Astrophys._ 567, A91 (2014).
  * [63] Planck Collaboration _et al._ Planck 2018 results. VI. Cosmological parameters.  _Astron. Astrophys._ 641, A6 (2020).
  * [64] Palacios, A., Arnould, M., & Meynet, G.  The thermonuclear production of 19F by Wolf-Rayet stars revisited.  _Astron. Astrophys._ 443, 243 (2005).
  * [65] Prantzos, N., Abia, C., Limongi, M., Chieffi, A. & Cristallo, S.  Chemical evolution with rotating massive star yields \- I. The solar neighbourhood and the s-process elements.  _Mon. Not. R. Astron. Soc._ 476, 3432–3459 (2018).
  * [66] Jönsson, H. _et al._ Chemical evolution of fluorine in the bulge. High-resolution K-band spectra of giants in three fields.  _Astron. Astrophys._ 564, A122 (2014).
  * [67] Hezaveh, Y. D. _et al._ ALMA Observations of SPT-discovered, Strongly Lensed, Dusty, Star-forming Galaxies.  _Astrophys. J._ 767, 132 (2013).
  * [68] Spilker, J. S. _et al._ ALMA Imaging and Gravitational Lens Models of South Pole Telescope—Selected Dusty, Star-Forming Galaxies at High Redshifts.  _Astrophys. J._ 826, 112 (2016).
  * [69] van der Wiel, M. H. D., Naylor, D. A., Makiwa, G., Satta, M. & Abergel, A.  Three-dimensional distribution of hydrogen fluoride gas toward NGC 6334 I and I(N).  _Astron. Astrophys._ 593, A37 (2016).
  * [70] Rodgers, S. D. & Charnley, S. B.  Chemical Evolution in Protostellar Envelopes: Cocoon Chemistry.  _Astrophys. J._ 585, 355–371 (2003).
  * [71] Jørgensen, J. K., Schöier, F. L. & van Dishoeck, E. F.  Molecular freeze-out as a tracer of the thermal and dynamical evolution of pre- and protostellar cores.  _Astron. Astrophys._ 435, 177–182 (2005).
  * [72] Rangwala, N. _et al._ Observations of Arp 220 Using Herschel-SPIRE: An Unprecedented View of the Molecular Gas in an Extreme Star Formation Environment.  _Astrophys. J._ 743, 94 (2011).
  * [73] Pickett, H. M. _et al._ Submillimeter, millimeter and microwave spectral line catalog.  _J. Quant. Spec. Radiat. Transf._ 60, 883–890 (1998).
  * [74] Phillips, T. G. _et al._ Herschel observations of EXtra-Ordinary Sources (HEXOS): Detection of hydrogen fluoride in absorption towards Orion KL.  _Astron. Astrophys._ 518, L109 (2010).
  * [75] van der Werf, P. P. _et al._ Black hole accretion and star formation as drivers of gas excitation and chemistry in Markarian 231.  _Astron. Astrophys._ 518, L42 (2010).
  * [76] Agúndez, M. _et al._ HIFI detection of hydrogen fluoride in the carbon star envelope IRC +10216.  _Astron. Astrophys._ 533, L6 (2011).
  * [77] van der Tak, F. F. S. _et al._ Detection of HF emission from the Orion Bar.  _Astron. Astrophys._ 537, L10 (2012).
  * [78] Kamenetzky, J. _et al._ Herschel-SPIRE Imaging Spectroscopy of Molecular Gas in M82.  _Astrophys. J._ 753, 70 (2012).
  * [79] Lu, N. _et al._ A Herschel Space Observatory Spectral Line Survey of Local Luminous Infrared Galaxies from 194 to 671 Microns.  _Astrophys. J. Supp. Series_ 230, 1 (2017).
  * [80] Pérez-Beaupuits, J. P. _et al._ A Thorough View of the Nuclear Region of NGC 253: Combined Herschel, SOFIA, and APEX Data Set.  _Astrophys. J._ 860, 23 (2018).


# Study of the $\Omega_{ccc}\Omega_{ccc}$ and $\Omega_{bbb}\Omega_{bbb}$
dibaryons in constituent quark model

Pablo Martín-Higueras<EMAIL_ADDRESS>Departamento de Ciencias
Integradas y Centro de Estudios Avanzados en Física, Matemática y Computación,
Universidad de Huelva, 21071 Huelva, Spain  David R. Entem<EMAIL_ADDRESS>Grupo de Física Nuclear, Universidad de Salamanca, E-37008 Salamanca, Spain
Instituto Universitario de Física Fundamental y Matemáticas (IUFFyM),
Universidad de Salamanca, E-37008 Salamanca, Spain  Pablo G. Ortega
<EMAIL_ADDRESS>Instituto Universitario de Física Fundamental y Matemáticas
(IUFFyM), Universidad de Salamanca, E-37008 Salamanca, Spain  Jorge Segovia
<EMAIL_ADDRESS>Departamento de Sistemas Físicos, Químicos y Naturales,
Universidad Pablo de Olavide, E-41013 Sevilla, Spain  Francisco Fernández
<EMAIL_ADDRESS>Grupo de Física Nuclear, Universidad de Salamanca, E-37008
Salamanca, Spain  Instituto Universitario de Física Fundamental y Matemáticas
(IUFFyM), Universidad de Salamanca, E-37008 Salamanca, Spain

###### Abstract

Dibaryons are the simplest system in which the baryon-baryon interaction, and
hence the underlying quark-quark interaction, can be studied in a clear way.
Although the only dibaryon known today is the deuteron (and possibly the
$d^{*}$), fully heavy dibaryons are good candidates for bound states because
in such systems the kinetic energy is small and the high symmetry of the wave
function favours binding. In this study, the possible existence of
$\Omega_{ccc}\Omega_{ccc}$ and $\Omega_{bbb}\Omega_{bbb}$ dibaryons is
investigated in the framework of a constituent quark model that satisfactorily
describes the deuteron, the $d^{*}(2380)$ and the $NN$ interaction.
$J^{P}=0^{+}$ candidates are found in both systems with binding energies of
the order of MeV.

Potential models, Quark models, Bottom charmed mesons, Exotic mesons

###### pacs:

12.39.Pn, 14.40.Lb, 14.40.Rt

##  I Introduction

Understanding the nucleon-nucleon interaction has been one of the priority
problems in Nuclear Physics since Yukawa’s one pion exchange theory. The
subsequent development of QCD paved the way to describe the strong
interactions in terms of quark degrees of freedom and facilitate to enlarge
the field to other flavors like charm en bottom.

Dibaryons are the simplest systems in which these studies can be addressed in
a transparent way. Until recently, the only well-established bound state of
two baryons was the deuteron. Then, in 2011, another unstable light dibaryon,
the $d^{*}(2380)$, was reported by the WASA-at-COSY Collaboration Adlarson
_et al._ (2011) from the double pionic fusion reaction $pn\to
d\pi^{0}\pi^{0}$. This resonance can be described as a nonstrange
$\Delta\Delta$ dibaryon with $I(J^{P})=0(3^{+})$. In 1989, Goldman noted that
due to the special symmetry of such a state, any model based on confinement
and gluon exchange should predict it Goldman  _et al._ (1989). The long
history of the search for dibaryons in the light quark sector can be found in
Ref Clement (2017).

It is well known that the binding of the deuteron is due to the coupling of
the ${}^{3}S_{1}$ and ${}^{3}D_{1}$ partial waves by one-pion exchange tensor
interactions. Similarly, the binding of the $d^{*}(2380)$ can be explained in
terms of Goldstone-boson exchanges Huang  _et al._ (2014). These two systems
then prove that the interaction binding these dibaryons arises from QCD chiral
symmetry breaking in the light quark sector.

Another interesting system is the fully heavy dibaryon. In such a system the
relativistic effects are negligible and the kinetic energy is small. As
originally pointed out by Bjorken Bjorken (1985), the triply-charmed baryon
$\Omega_{ccc}$ is stable against strong interactions. This fact opens the
possibility to study systems like $\Omega_{ccc}\Omega_{ccc}$ or
$\Omega_{bbb}\Omega_{bbb}$. Moreover, in contrast to the deuteron and the
$d^{*}$ case, the latter systems provide an ideal scenario to explore the
baryon-baryon interaction in an environment free of chiral dynamics.

In this work we will focus on the study of the fully heavy dibaryons. Two
recent Lattice QCD calculations have explored these systems: Ref. Lyu  _et
al._ (2021) showed that $\Omega_{ccc}\Omega_{ccc}$ is loosely bound by
$5.68(0.77)$ MeV, while Ref. Mathur  _et al._ (2023) found a very deep
$\Omega_{bbb}\Omega_{bbb}$ state with a binding energy of $81_{-16}^{+14}$
MeV. These conclusions are confirmed by several quark model calculations, but
are contradicted by others. For example, Huang et al. Huang  _et al._ (2022),
using a constituent quark model based on the one-gluon exchange interaction
and the resonating group method, studied the possible bound states of the
$\Omega_{ccc}\Omega_{ccc}$ and $\Omega_{bbb}\Omega_{bbb}$, among others. They
found a $J^{P}=0^{+}$ bound state for the $\Omega_{ccc}\Omega_{ccc}$ system
with a binding energy of $2.5$ MeV and another $\Omega_{bbb}\Omega_{bbb}$
state bound by $0.9$ MeV, contrary to naive espectations. Deng Deng (2023)
performed a study of the di-$\Delta^{++}$, di-$\Omega_{ccc}$ and
di-$\Omega_{bbb}$ systems using a naive one-gluon exchange quark model and a
chiral quark model including $\pi$ and $\sigma$ exchanges between quarks.
Obviously, in this case these parts of the interaction apply only to the light
quarks, but the set of parameters is different in the two models. Both studies
predict very shallow di-$\Omega_{ccc}$ and di-$\Omega_{bbb}$ states with
binding energies around $1$ MeV. Using a different model, namely QCD sum
rules, Wang Wang (2022) found for each di-$\Omega_{ccc}$ and di-$\Omega_{bbb}$
systems two $J^{P}=0^{+}$ and $J^{P}=1^{-}$ states that are slightly below
their respective thresholds.

On the other hand, several studies within the quark model have ruled out the
existence of fully heavy dibaryons. In Ref Richard  _et al._ (2020) the
authors investigated the existence of $bbbccc$ dibaryons and extrapolated
their results to the properties of the $bbbbbb$ and $cccccc$ systems. They
found no bound states for $\Omega_{ccc}\Omega_{ccc}$ or
$\Omega_{bbb}\Omega_{bbb}$ combinations. On the other hand, Alcaraz-Peregrina
et al. Alcaraz-Pelegrina and Gordillo (2022) used the Difussion Montecarlo
technique to describe fully heavy compact six-quark arrangements. They found
that all the hexaquarks have smaller masses than those of their constituents,
i.e., all the hexaquarks are bound systems. However, their masses are also
larger than those of any pair of baryons into which they can be divided. This
means that each hexaquark is unstable with respect to its splitting into two
baryons. Finally, two more calculations, in the framework of the constituent
quark model Lü  _et al._ (2022) or the extended chromomagnetic model Weng and
Zhu (2024), showed that all the fully heavy dibaryons lie above their
corresponding baryon-baryon thresholds.

In view of this controversial situation, since different approaches lead to
quite different conclusions, we will study the possible existence of
$\Omega_{ccc}\Omega_{ccc}$ and $\Omega_{bbb}\Omega_{bbb}$ dibaryons using the
constituent quark model of Ref. Fernandez  _et al._ (1993) and its extension
to the heavy quark sector Vijande  _et al._ (2005); Segovia  _et al._ (2008),
which has been able to describe a large variety of hadronic phenomenology. In
particular, the model reproduces the properties of the deuteron Valcarce  _et
al._ (1994); Juliá-Díaz  _et al._ (2002) and predicts the existence of the
$d^{*}(2380)$ as a $\Delta\Delta$ dibaryon Garcilazo  _et al._ (1997);
Valcarce  _et al._ (2001). Although the binding energy of the $d^{*}$
predicted in the latter references is smaller than the experimental value, it
is laso worth mentioning that the these studies were performed without
coupling to the $NN$ channel.

The paper is structured as follows. In Sec. II we describe the main aspects of
our theoretical model, giving details about the wave functions used to
describe $\Omega_{ccc}$ ($\Omega_{bbb}$) baryons and the way we derive the
$\Omega_{ccc}\Omega_{ccc}$ interaction using the Resonating Group Method
(RGM). Section IV is devoted to presenting our results for the possible
dibaryons. Finally, we summarize and give some conclusions in Sec. V.

##  II Theoretical formalism

###  II.1 The constituent quark model

Our theoretical framework is a QCD-inspired constituent quark model (CQM)
proposed in Ref. Vijande  _et al._ (2005) and extended to the heavy quark
sector in Ref. Segovia  _et al._ (2008). The main pieces of the model are the
constituent light quark masses and Goldstone-boson exchanges, which appears as
consequences of spontaneous chiral symmetry breaking of the QCD Lagrangian
together with perturbative one-gluon exchange (OGE) and nonperturbative color
confining interactions.

Following Diakonov Diakonov (2003), a simple Lagrangian invariant under chiral
transformations can be written as

${\mathcal{L}}=\bar{\psi}(i\,{/\\!\\!\\!\partial}-M(q^{2})U^{\gamma_{5}})\,\psi\,,$
(1)

where $M(q^{2})$ is the dynamical (constituent) quark mass and
$U^{\gamma_{5}}=e^{i\lambda_{a}\phi^{a}\gamma_{5}/f_{\pi}}$ is the matrix of
Goldstone-boson fields that can be expanded as

$U^{\gamma_{5}}=1+\frac{i}{f_{\pi}}\gamma^{5}\lambda^{a}\pi^{a}-\frac{1}{2f_{\pi}^{2}}\pi^{a}\pi^{a}+\ldots$
(2)

The first term of the expansion generates the constituent quark mass, while
the second term gives rise to a one-boson exchange interaction between quarks.
The main contribution of the third term comes from the two-pion exchange which
has been simulated by means of a scalar-meson exchange potential.

In the heavy quark sector, chiral symmetry is explicitly broken and Goldstone-
boson exchange does not occur. However, the full interaction constrains the
model parameters through the light-meson phenomenology Vijande  _et al._
(2005); Segovia  _et al._ (2008). Thus, OGE and confinement are the only
remaining interactions between the heavy quarks.

The OGE potential is generated from the vertex Lagrangian

${\mathcal{L}}_{qqg}=i\sqrt{4\pi\alpha_{s}}\,\bar{\psi}\gamma_{\mu}G^{\mu}_{c}\lambda^{c}\psi,$
(3)

where $\lambda^{c}$ are the $SU(3)$ colour matrices, $G^{\mu}_{c}$ is the
gluon field and $\alpha_{s}$ is the strong coupling constant. The scale
dependence of $\alpha_{s}$ allows a consistent description of light, strange
and heavy mesons. Its explicit expression can be found in, e.g., Ref. Vijande
_et al._ (2005),

$\alpha_{s}(\mu)=\frac{\alpha_{0}}{\ln\left(\frac{\mu^{2}+\mu_{0}^{2}}{\Lambda_{0}^{2}}\right)}$
(4)

Regarding the confinement potential, it is well known that multi-gluon
exchanges produce an attractive linearly rising potential proportional to the
distance between infinite-heavy quarks Mateu  _et al._ (2019). However, sea
quarks are also important components of the strong interaction dynamics that
contribute to the screening of the rising potential at low momenta and
eventually to the breaking of the quark-antiquark binding string Bali  _et
al._ (2005). Our model tries to mimic this behaviour with a screening
potential at high distances.

Then, the full interaction between heavy quarks is given by

$\displaystyle V_{ij}(r)$
$\displaystyle=\bigg{[}-a_{c}(1-e^{-\mu_{c}r})+\Delta+\frac{\alpha_{s}(\mu)}{4}\frac{1}{r}\bigg{]}(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j})$
$\displaystyle V_{ij}^{S}(r)$
$\displaystyle=-\frac{\alpha_{s}(\mu)}{4}\frac{1}{6m_{i}m_{j}}\frac{e^{-r/r_{0}(\mu)}}{rr_{0}^{2}(\mu)}(\vec{\sigma}_{i}\cdot\vec{\sigma}_{j})(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j})$
$\displaystyle V_{ij}^{T}(r)$
$\displaystyle=-\frac{1}{16}\frac{\alpha_{s}(\mu)}{m_{i}m_{j}}S_{ij}\,(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j})\times$
$\displaystyle\times\bigg{[}\frac{1}{r^{3}}-\frac{e^{-r/r_{g}(\mu)}}{r}\bigg{(}\frac{1}{r^{2}}+\frac{1}{3r_{g}^{2}(\mu)}+\frac{1}{rr_{g}(\mu)}\bigg{)}\bigg{]}$
(5)

where $r_{0}(\mu)=\hat{r}_{0}\frac{m_{n}}{2\mu}$ with $\mu$ the reduced mass
of the $(ij)$ heavy quark pair, $\vec{\lambda}$ are the colour matrices,
$\vec{\sigma}$ the spin matrices and
$S_{ij}=3(\vec{\sigma}_{i}\cdot\hat{r})(\vec{\sigma}_{j}\cdot\hat{r})-(\vec{\sigma}_{i}\cdot\vec{\sigma}_{j})$
the tensor operator of the $(ij)$ pair with $\vec{r}$ their relative position.

All the parameters of the model are given in Table 1. We have not included the
spin-orbit interaction parts coming from the one-gluon exchange and
confinement because they should give small contributions in this calculation.
For the same reason, the spin-tensor terms are neglected in the calculation of
the $\Omega_{ccc}$ ($\Omega_{bbb}$) masses, but are included in the
$\Omega_{ccc}$$\Omega_{ccc}$ ($\Omega_{bbb}$$\Omega_{bbb}$) interaction.

Quark masses (MeV) | $m_{c}$ | 1763
---|---|---
| $m_{b}$ | 5110
OGE |  $\hat{r}_{0}$ (fm) | 0.181
| $\alpha_{0}$ | 2.118
|  $\Lambda_{0}$ (fm-1) | 0.113
|  $\mu_{0}$ (MeV) | 36.976
Confinement |  $a_{c}$ (MeV) | 507.4
|  $\mu_{c}$ (fm-1) | 0.576
|  $\Delta$ (MeV) | 184.432
Table 1:  Parameters for the quark-quark interaction.

###  II.2 The wave function of the $\Omega_{ccc}$($\Omega_{bbb}$)

A precise definition of the wave functions of the $\Omega_{ccc}$ and
$\Omega_{bbb}$ baryons (henceforth $\Omega_{QQQ}$) is an essential part of the
calculation, because it defines the size of the baryon, which is important for
the baryon-baryon interaction.

Once we know the quark-quark interaction, the $\Omega_{QQQ}$ wave function can
be calculated by solving the Schrödinger equation with the Gaussian Expansion
Method (GEM) Hiyama  _et al._ (2003). In the GEM framework one makes an
expansion in gaussian wave functions but instead of using only one set of
Jacobi coordinates, one includes the lowest orbital angular momentum wave
functions using the three sets of possible Jacobi coordinates. The reason to
use different sets is that lowest angular momentum wave functions in one set
generates higher angular momentum wave functions in the other sets, making a
very numerically efficient way to include such high angular momentum
components.

However the wave function given by GEM would be quite complicate and would
make the calculation of the dibaryon interaction slow. Alternatively, for the
calculation of the dibaryon interaction (that would be justified later), the
following orbital wave function can be used

$\displaystyle\phi(\vec{p}_{\xi_{1}},\vec{p}_{\xi_{2}})$
$\displaystyle=\bigg{[}\frac{2b^{2}}{\pi}\bigg{]}^{3/4}e^{-b^{2}p_{\xi_{1}}^{2}}\bigg{[}\frac{3b^{2}}{2\pi}\bigg{]}^{3/4}e^{-\frac{3b^{2}}{4}p_{\xi_{1}}^{2}}$
(6)

where $p_{\xi_{i}}$ are the Jacobi coordinates defined as

$\displaystyle\vec{p}_{\xi_{1}}$
$\displaystyle=\frac{1}{2}(\vec{p}_{1}-\vec{p}_{2})$
$\displaystyle\vec{p}_{\xi_{2}}$
$\displaystyle=\frac{2}{3}\vec{p}_{3}-\frac{1}{3}(\vec{p}_{1}+\vec{p}_{2})$
(7)

In the notation we use for the baryon calculation this corresponds to mode 3
of the GEM basis using only one gaussian with angular momentum zero and the
parameters $\nu=\frac{1}{4b^{2}}$ and $\lambda=\frac{1}{3b^{2}}$. Notice that
fixing the relation between the parameters of the gaussians $\nu$ and
$\lambda$ to these values ($\nu=\frac{3}{4}\lambda$), the orbital wave
functions is totally symetric, necessary to get a totally antisymetric wave
function for the baryon of lowest energy. The spin wave function has to be
also symmetric and implies $S=\frac{3}{2}$ and the color wave function will be
a color singlet.

So our wave function for the baryon is

$\displaystyle\psi_{B}$
$\displaystyle=\phi(\vec{p}_{\xi_{1}},\vec{p}_{\xi_{2}})\chi_{B}\xi_{c}[1^{3}]$
(8)

with $\chi_{B}=((\frac{1}{2}\frac{1}{2})1\frac{1}{2})\frac{3}{2})$ the spin
wave function and $\xi_{c}[1^{3}]$ a singlet color wave function.

Using the wave function of Eq. (6), the kinetic energy is given by

$\displaystyle T$
$\displaystyle=\langle\psi_{B}|\frac{p_{\xi_{1}}^{2}}{m}+\frac{3p_{\xi_{2}}^{2}}{4m}|\psi_{B}\rangle=\frac{3}{2mb^{2}}$
(9)

For the interaction energy we can evaluate
$\langle\psi_{B}|V_{12}|\psi_{B}\rangle$ and multiply by 3, since we have 3
interactions between equivalent quarks. It is easier to evaluate it in
coordinate space. The wave function in coordinate space is

$\displaystyle\phi_{B}(r_{3},R_{3})=\bigg{[}\frac{1}{2\pi
b^{2}}\bigg{]}^{3/4}e^{-\frac{r_{3}^{2}}{4b^{2}}}\bigg{[}\frac{2}{3\pi
b^{2}}\bigg{]}^{3/4}e^{-\frac{R_{3}^{2}}{3b^{2}}}$ (10)

and so

$\displaystyle\langle\psi_{B}|V_{12}|\psi_{B}\rangle$
$\displaystyle=4\pi\bigg{[}\frac{1}{2\pi
b^{2}}\bigg{]}^{3/2}\int_{0}^{\infty}r_{3}^{2}dr_{3}e^{-\frac{r_{3}^{2}}{2b^{2}}}V(r_{3})$
(11)

The mean value distance between quarks is given by

$\displaystyle\sqrt{\langle r_{ij}^{2}\rangle}$ $\displaystyle=\sqrt{3}b$ (12)

and the mass is given by

$\displaystyle M=3m_{b}+T+3\langle\psi_{B}|V_{12}|\psi_{B}\rangle$ (13)

Finally, the value of the $b$ parameter is obtained by minimizing the mass

$\displaystyle\frac{\partial M}{\partial b}=0$ (14)

| $\Omega_{ccc}$ | $\Omega_{bbb}$
---|---|---
| $\partial M/\partial b$ | GEM | $\partial M/\partial b$ | GEM
$M$ [MeV] | 4810.9 | 4798.6 | 14413.8 | 14396.9
$\sqrt{\langle r_{ij}^{2}\rangle}$ [fm] | 0.4360 | 0.4432 | 0.2716 | 0.2762
$\langle T\rangle$ [MeV] | 522.9 | 522.8 | 465.0 | 471.5
$\langle V\rangle$ [MeV] | -333.7 | -337.7 | -460.4 | -468.2
Table 2:  Parameters for the $\Omega_{ccc}$ and $\Omega_{bbb}$ baryons
obtained from the mass minimization procedure ($\frac{\partial M}{\partial
b}=0$) and the Gaussian Expansion Method.

Although the GEM method provides a more complete description of the wave
function as mentioned before, the calculation is simplified if we use the
analytical wave function of Eq. (6). In Table 2 we show the results of the
$\Omega_{ccc}$ and $\Omega_{bbb}$ wave functions using the mass minimization
procedure and compare with the GEM solution to justify the use of the simple
wave function given by Eq. (6). We see that we get a resonable agreement for
the sizes and energies in both cases, although the agreement is better in the
beauty sector. The minimal values for $b$ are given by $b_{\rm min}=0.15679$
fm for the $\Omega_{bbb}$ and $b_{\rm min}=0.25172$ fm for the $\Omega_{ccc}$.

##  III The $\Omega_{QQQ}\Omega_{QQQ}$ interaction

The system under study has six identical quarks. Then, the baryon-baryon total
wave function must be fully antisymmetric. As the wave function of the baryons
is already antisymmetric, the antisymmetrizer operator is just given by,

$\displaystyle{\mathcal{A}}=1-9P_{36}$ (15)

In order to obtain the effective baryon-baryon interaction from the underlying
quark dynamics we use the Resonating Group Method (RGM) Wheeler (1937); Tang
_et al._ (1978). Then, we need to solve the projected Schrödinger equation,

$\displaystyle 0=\bigg{(}\frac{p^{\prime
2}}{2\mu_{\Omega\Omega}}-E\bigg{)}\chi(\vec{P}^{\prime})+\int$
$\displaystyle\bigg{(}\,^{\rm RGM}V_{D}(\vec{P}^{\prime},\vec{P}_{i})+$
$\displaystyle+\,^{\rm
RGM}K(\vec{P}^{\prime},\vec{P}_{i})\bigg{)}\chi(\vec{P}_{i})d^{3}P_{i}$ (16)

where $\vec{P}^{\prime}$ ($\vec{P}_{i}$) is the relative
$\Omega_{QQQ}-\Omega_{QQQ}$ final (initial) momentum, $E=E_{T}-2M_{\Omega}$
the relative energy of the system with respect to the threshold, ${}^{\rm
RGM}V_{D}(\vec{P}^{\prime},\vec{P}_{i})$ is the direct kernel and ${}^{\rm
RGM}K(\vec{P}^{\prime},\vec{P}_{i})$ is the exchange kernel and
$\mu_{\Omega\Omega}$ is the reduced mass of two $\Omega_{QQQ}$ baryons.

Here, $M_{\Omega}$ is

$\displaystyle M_{\Omega}$
$\displaystyle=3m_{b}+\frac{3}{2m_{Q}b^{2}}+3E_{int}$ (17) $\displaystyle
E_{int}$ $\displaystyle=\langle V_{ij}\rangle=\int
d^{3}qe^{-\frac{q^{2}b^{2}}{2}}\langle V_{ij}(q)\rangle$ (18)

The direct term will be zero in the present model since the color
coefficients, $(\vec{\lambda}_{i}\cdot\vec{\lambda}_{j})$, are zero between
color singlets.

Then, the full interaction is driven by exchange diagrams, which take into
account the quark rearrangement between baryons. The exchange kernel can be
written as,

$\,{}^{\rm RGM}K(\vec{P}^{\prime},\vec{P}_{i})=$ $\,{}^{\rm
RGM}T(\vec{P}^{\prime},\vec{P}_{i})+^{\rm
RGM}V_{ij\,E}(\vec{P}^{\prime},\vec{P}_{i})-$ $\displaystyle-E_{T}\,^{\rm
RGM}N(\vec{P}^{\prime},\vec{P}_{i})$ (19)

where ${}^{\rm RGM}T(\vec{P}^{\prime},\vec{P}_{i})$ is the exchange kinetic
term, ${}^{\rm RGM}N(\vec{P}^{\prime},\vec{P}_{i})$ is a normalization term
and ${}^{\rm RGM}V_{ij\,E}(\vec{P}^{\prime},\vec{P}_{i})$ is the exchange
potential (for explicit expressions see, e.g., Refs. Ortega  _et al._ (2022);
Fernández  _et al._ (2020)).

##  IV Results

Let us first study the $\Omega_{bbb}\Omega_{bbb}$ system. One of the states in
$S$ wave is the $J^{P}=0^{+}$, which corresponds to the ${}^{1}S_{0}$ and
${}^{5}D_{0}$ partial waves. As in the case of the deuteron, $S$ and $D$ waves
are mixed. We first calculate the binding energy considering the parameter $b$
and the reduced mass given by the minimization procedure. Without tensor
interactions they are decoupled and only a bound state appears in the
${}^{1}S_{0}$ partial wave. The binding energy of this state is $E=-1.9859$
MeV. The ${}^{5}D_{0}$ partial wave is not bound. If we include the tensor
interaction of OGE the partial waves are coupled and the binding energy
increases very slightly to $E=-1.9876$ MeV. The probability of the $D$ wave is
only $6.6\cdot 10^{-4}\%$. As in the deuteron, the binding energy has a
sizeable cancellation between the kinetic and interaction parts. The mean
value of the kinetic energy is $\langle T\rangle=11.2$ MeV, while for the
interaction we have $\langle V\rangle=-13.2$ MeV. The confinement interaction
dominates and gives the needed attraction to bind the system. If we exclude
the OGE we get $E=-7.3698$ MeV with $\langle T\rangle=20.9$ MeV and $\langle
V\rangle=-28.3$ MeV.

The potential for the ${}^{1}S_{0}$ partial wave is given in Fig. 1. The
relative wave functions are shown in Fig. 2.

Figure 1:  The potential $V(p^{\prime},p)$ in the ${}^{1}S_{0}$ partial wave
for the $\Omega_{bbb}\Omega_{bbb}$ interaction.

Figure 2:  The relative wave functions in the ${}^{1}S_{0}$ and ${}^{5}D_{0}$
partial waves for the $\Omega_{bbb}\Omega_{bbb}$ dibaryon.

If we consider $b=\frac{\sqrt{r_{ij}^{2}}}{\sqrt{3}}$ and the reduced mass
given by the GEM calculation we get a binding energy in the coupled case of
$E=-1.81$ MeV. With $b$ given by the minimization procedure and the reduced
mass is given by GEM we get $=-1.9754$ MeV. The effect of the different
reduced mass is very small and dominates the effect of the different $b$
parameters. In principle with only one gaussian one should use the value given
by the minimization procedure, but this gives us a feeling of the uncertainty
due to the simplification of the wave function. Although the binding energy
varies a little bit, in both cases the system is bounded.

Another possible state is the $J^{P}=2^{+}$ which includes the ${}^{5}S_{2}$,
${}^{1}D_{2}$, ${}^{5}D_{2}$ and ${}^{5}G_{2}$ partial waves. None of them are
bound. One could expect a bound state for the ${}^{5}S_{2}$ partial wave, but
in this case one can see that the potential coming from the
$\vec{\lambda}_{i}\cdot\vec{\lambda}_{j}$ have opposite sign for $S=2$ with
respect to $S=0$. So if we have attraction for the $S=0$, this implies
repulsion for $S=2$. Higher partial waves are more difficult to bind.

Antisymmetry implies $L+S=$ even and parity is given by $P=(-1)^{L}$. So for
$P=+$, the spin $S$ has to be even. This means that $1^{+}$ and $3^{+}$ can be
only in $D$ or $G$ waves, which will be difficult to bind as it was seen for
the $0^{+}$ and $2^{+}$ states. In more detail:

  * •

We start with the $1^{+}$ state and include ${}^{5}D_{1}$ which is the only
partial wave. It should be the same as the ${}^{5}D_{2}$ partial wave with the
exception of the contribution of the OGE tensor interaction. It does not bind.

  * •

For the $3^{+}$ state we have the ${}^{5}D_{3}$ and ${}^{5}G_{3}$ partial
waves and they do not bind.

We give in Fig. 3 the Fredholm determinant for the 4 different $J^{+}$ quantum
numbers, where we can see that only the $0^{+}$ channel binds.

Figure 3:  The Fredholm determinant of the $\Omega_{bbb}\Omega_{bbb}$ system
for the $0^{+}$, $1^{+}$, $2^{+}$ and $3^{+}$ channels. Only the $0^{+}$
crosses the zero.

Regarding possible $P=-$ states, this would imply odd partial waves and odd
total spin. We have analyzed the $J^{P}=\\{0^{-},1^{-},2^{-},3^{-}\\}$,
finding no additional bound states. Again, in Fig. 4 the Fredholm determinant
for the 4 different $J^{-}$ quantum numbers is shown, where we can see that no
bound state is predicted.

Figure 4:  The Fredholm determinant of the $\Omega_{bbb}\Omega_{bbb}$ system
for the $0^{-}$, $1^{-}$, $2^{-}$ and $3^{-}$ channels. None of them crosses
the zero.

Concerning the $\Omega_{ccc}\Omega_{ccc}$ system, the situation is similar to
the $\Omega_{bbb}\Omega_{bbb}$ system, and we only find a bound state in the
$0^{+}$ channel. The binding energy is $E=-0.7104$ MeV with a $D$-state
probability of $1.7\cdot 10^{-3}\%$. The mean values of kinetic and
interaction terms are $\langle T\rangle=7.46$ MeV and $\langle V\rangle=-8.17$
MeV. In this case we used the $b$ parameter from the minimization procedure
and the reduced mass from the GEM. Using the reduced mass from the
minimization parameter the binding energy changes to $E=-0.7288$ MeV and both
parameters from the GEM to $E=-0.62$ MeV.

###  IV.1 Dependence on the model parameters

We analyze the dependence on the parameters of the model for the $J^{P}=0^{+}$
state to see in which parameter space region the system will not bind. In all
cases we use the minimization procedure to obtain $b$ and
$\mu_{\Omega\Omega}$.

The dependence on the quark mass $m_{q}$ is shown in Fig. 5. Notice that some
of the parameters of the potential depends on $m_{q}$ since we use scale
dependent parameters. We see that the system binds reducing the quark mass up
to $m_{q}\sim 800-900$ MeV.

Figure 5:  Quark mass dependence of $b$, $M_{QQQ}$ and $E$. We show by a dot
the result of the Chiral Quark Model.

Our model has an effective string tension given by

$\displaystyle\sigma=\frac{8}{3}a_{c}\mu_{c}=0.1537\,{\rm GeV}^{2}$ (20)

We plot the parameters $b$, $M_{QQQ}$ and $E$ as a function of the string
tension in Fig. 6. We vary the value of $\mu_{c}$ from 0.15 to 0.85 fm-1 and
leave $a_{c}$ unchanged so the saturation energy does not change.

We see that for higher string tension values (our value is lower than some
determinations) the binding energy will increase.

Figure 6:  Dependence of $b$, $M_{QQQ}$ and $E$ on the effective string
tension $\sigma=\frac{8}{3}a_{c}\mu_{c}$ for $m_{q}=m_{c}$ (left) and
$m_{q}=m_{b}$ (right). We show by a dot the result of the Chiral Quark Model.

Our confinement effective potential is

$\displaystyle V(r)\sim\sigma\frac{1-e^{-\mu_{c}r}}{\mu_{c}}$ (21)

We vary the value of $\mu_{c}$ from 0.15 to 0.85 fm-1 and change $a_{c}$ so
that $\sigma$ does not change. This is the same interval we used when we
changed the string tension $\sigma$. The saturation energy changes as
$\frac{\sigma}{\mu_{c}}$. Notice that the interaction region is
$\sim\sqrt{3}b$ so if $x\equiv\sqrt{3}b\mu_{c}\ll 1$ the potential in the
interacting region is basically linear. In this calculation we got $x=0.062$
to $x=0.39$ in the charm sector and $x=0.039$ to $x=0.24$ in the bottom
sector. For $\mu_{c}\to 0$ the potential becomes more linear in the
interaction region.

The results varying the saturation are shown in Fig. 7. We see that the
dependence on the saturation point of the properties of the $\Omega_{QQQ}$,
$b$ and $M_{QQQ}$, is smaller than on the string tension $\sigma$ as one would
expect. For the binding energy of the $\Omega_{QQQ}$ dibaryon we see also an
smaller dependence.

Notice that when $\mu_{c}\to 0$ the binding energy increases, so a linear
confinement potential should give more binding.

Figure 7:  Dependence of $b$, $M_{QQQ}$ and $E$ on the saturation parameter
$\mu_{c}$ for $m_{q}=m_{c}$ (left) and $m_{q}=m_{b}$ (right). We show by a dot
the result of the Chiral Quark Model.

Finally we can study the dependence of the binding energy on the size of the
baryon. For that we keep all the parameters unchanged and only vary the
parameter $b$ on the RGM calculation. Results are shown in Fig. 8. With bigger
sizes we get less binding but it has to be increased much more than the
difference between the sizes of the variational and GEM calculation, which
shows that using the exact wave function the system will still bind. This
argument is more robust for the bottom sector but it should also work in the
charm.

The result should be seen as an upper bound of the binding energy, since we
are using a variational calculation. Also other channels may be involved, but
since we are considering the lower energy channel, including more channels
will provide more attraction. We can conclude that the Chiral Quark Model
binds the $\Omega_{QQQ}\Omega_{QQQ}$ system in both cases, when $Q$ is a
bottom or a charm quark. These molecular states are analogs of two-atom
molecules, where the direct interaction is zero for neutral atoms, as in our
model for colorless objects.

Figure 8:  Dependence of $E$ on the size of the baryon given by $b$. Charm
sector on the left and bottom sector on the right. We show by a dot the result
of the Chiral Quark Model.

##  V Summary

In this work we have studied the possible existence of fully-heavy dibaryons
in the charm and bottom sectors. The main conclusion we found is that, using a
wave function which minimizes the mass of the $\Omega_{ccc}$ ($\Omega_{bbb}$)
baryons, the six $c$ quarks or the six $b$ quarks can form bound states with
$J^{P}=0^{+}$ quantum numbers. The binding energy of the charm dibaryon is
$E_{b}=-0.71$ MeV, while in the bottom case the binding energy is slightly
higher, $E_{b}=-1.98$ MeV, which is reasonable due to the highest mass of the
bottom quark. The $J^{P}=0^{+}$ state corresponds to the coupling of
${}^{1}S_{0}$ and ${}^{5}D_{0}$ partial waves, but with a very small
${}^{5}D_{0}$ component. No further bound states are found in other partial
waves.

###### Acknowledgements.

This work has been partially funded by EU Horizon 2020 research and innovation
program, STRONG-2020 project, under grant agreement no. 824093; Ministerio
Español de Ciencia e Innovación under grant nos. PID2022-141910NB-I00 and
PID2022-140440NB-C22; and Junta de Andalucía under contract Nos. PAIDI FQM-370
and PCI+D+i under the title: ”Tecnologías avanzadas para la exploración del
universo y sus componentes” (Code AST22-0001).

## References

  * Adlarson  _et al._ (2011) P. Adlarson _et al._ (WASA-at-COSY), Phys. Rev. Lett. 106, 242302 (2011), arXiv:1104.0123 [nucl-ex] .
  * Goldman  _et al._ (1989) J. T. Goldman, K. Maltman, G. J. Stephenson, Jr., K. E. Schmidt, and F. Wang, Phys. Rev. C 39, 1889 (1989).
  * Clement (2017) H. Clement, Prog. Part. Nucl. Phys. 93, 195 (2017), arXiv:1610.05591 [nucl-ex] .
  * Huang  _et al._ (2014) H. Huang, J. Ping, and F. Wang, Phys. Rev. C 89, 034001 (2014).
  * Bjorken (1985) J. D. Bjorken, AIP Conf. Proc. 132, 390 (1985).
  * Lyu  _et al._ (2021) Y. Lyu, H. Tong, T. Sugiura, S. Aoki, T. Doi, T. Hatsuda, J. Meng, and T. Miyamoto, Physical Review Letters 127 (2021), 10.1103/physrevlett.127.072003.
  * Mathur  _et al._ (2023) N. Mathur, M. Padmanath, and D. Chakraborty, Physical Review Letters 130 (2023), 10.1103/physrevlett.130.111901.
  * Huang  _et al._ (2022) H. Huang, J. Ping, X. Zhu, and F. Wang, Eur. Phys. J. C 82, 805 (2022), arXiv:2011.00513 [hep-ph] .
  * Deng (2023) C.-R. Deng, Physical Review D 108 (2023), 10.1103/physrevd.108.054037.
  * Wang (2022) Z.-G. Wang, International Journal of Modern Physics A 37 (2022), 10.1142/s0217751x22501664.
  * Richard  _et al._ (2020) J.-M. Richard, A. Valcarce, and J. Vijande, Phys. Rev. Lett. 124, 212001 (2020), arXiv:2005.06894 [hep-ph] .
  * Alcaraz-Pelegrina and Gordillo (2022) J. Alcaraz-Pelegrina and M. Gordillo, Physical Review D 106 (2022), 10.1103/physrevd.106.114028.
  * Lü  _et al._ (2022) Q.-F. Lü, D.-Y. Chen, and Y.-B. Dong, (2022), arXiv:2208.03041 [hep-ph] .
  * Weng and Zhu (2024) X.-Z. Weng and S.-L. Zhu, The European Physical Journal C 84 (2024), 10.1140/epjc/s10052-024-12395-4.
  * Fernandez  _et al._ (1993) F. Fernandez, A. Valcarce, U. Straub, and A. Faessler, J. Phys. G 19, 2013 (1993).
  * Vijande  _et al._ (2005) J. Vijande, F. Fernandez, and A. Valcarce, J. Phys. G31, 481 (2005), arXiv:hep-ph/0411299 [hep-ph] .
  * Segovia  _et al._ (2008) J. Segovia, A. M. Yasser, D. R. Entem, and F. Fernández, Phys. Rev. D 78, 114033 (2008).
  * Valcarce  _et al._ (1994) A. Valcarce, A. Buchmann, F. Fernández, and A. Faessler, Phys. Rev. C 50, 2246 (1994).
  * Juliá-Díaz  _et al._ (2002) B. Juliá-Díaz, D. R. Entem, A. Valcarce, and F. Fernández, Phys. Rev. C 66, 047002 (2002).
  * Garcilazo  _et al._ (1997) H. Garcilazo, F. Fernández, A. Valcarce, and R. D. Mota, Phys. Rev. C 56, 84 (1997).
  * Valcarce  _et al._ (2001) A. Valcarce, H. Garcilazo, R. D. Mota, and F. Fernandez, J. Phys. G 27, L1 (2001).
  * Diakonov (2003) D. Diakonov, Prog. Part. Nucl. Phys. 51, 173 (2003), arXiv:hep-ph/0212026 [hep-ph] .
  * Mateu  _et al._ (2019) V. Mateu, P. G. Ortega, D. R. Entem, and F. Fernández, Eur. Phys. J. C 79, 323 (2019), arXiv:1811.01982 [hep-ph] .
  * Bali  _et al._ (2005) G. S. Bali, H. Neff, T. Duessel, T. Lippert, and K. Schilling (SESAM), Phys. Rev. D71, 114513 (2005), arXiv:hep-lat/0505012 [hep-lat] .
  * Hiyama  _et al._ (2003) E. Hiyama, Y. Kino, and M. Kamimura, Prog. Part. Nucl. Phys. 51, 223 (2003).
  * Wheeler (1937) J. A. Wheeler, Phys. Rev. 52, 1083 (1937).
  * Tang  _et al._ (1978) Y. C. Tang, M. Lemere, and D. R. Thompson, Phys. Rept. 47, 167 (1978).
  * Ortega  _et al._ (2022) P. G. Ortega, D. R. Entem, and F. Fernández, “Quark Models for Baryon-Baryon Interactions,” in  _Handbook of Nuclear Physics_, edited by I. Tanihata, H. Toki, and T. Kajino (2022) pp. 1–26.
  * Fernández  _et al._ (2020) F. Fernández, P. G. Ortega, and D. R. Entem, Front. in Phys. 7, 233 (2020).