††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# Exploring light nuclei production at RHIC and LHC energies with A Multi-
Phase Transport model and a coalescence afterburner

Yoshini Bailung1 Neha Shah2 Ankhi Roy1 1Department of Physics, Indian
Institute of Technology Indore, Simrol, Indore, Madhya Pradesh, India
2Department of Physics, Indian Institute of Technology Patna, Bihta, Patna,
Bihar, India

###### Abstract

In heavy-ion collisions, understanding how light nuclei species are produced
can provide insight into the nature of hadronic interactions in extreme
conditions. It can also shed light on understanding the matter-antimatter
asymmetry and dark matter searches in astrophysical processes. To investigate
the production mechanism of light nuclei such as deuteron, triton, and
helium-3, we use a naive coalescence afterburner coupled to the well-known
$``$A Multi-Phase Transport model” (AMPT). We focus on studying the production
of light nuclei in central Au+Au collisions at different center of mass
energies ($\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, and 200 GeV) and in Pb+Pb
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV, at mid-rapidity. We
generate events with the string melting version of AMPT, and feed the
information of the nucleons with spatial and momentum conditions into the
coalescence afterburner. Our study reports differential and integrated yields
in transverse momentum ($p_{\rm{T}}$) of the light nuclei in different center
of mass energies. We also estimate the coalescence parameters ($B_{A}$) as a
function of $p_{\rm{T}}$ and collision energy for (anti-)deuterons, tritons
and helium-3s for Au+Au and Pb+Pb collisions, which are compared to other
light nuclei production studies. All results are compared with measurements
from the STAR and ALICE experiments.

##  I Introduction

Experiments at the Relativistic Heavy-Ion Collider (RHIC), BNL, USA
STAR:2005gfr  and the Large Hadron Collider (LHC), CERN, Geneva, ALICE:2008ngc
are designed to conduct high energy hadronic and heavy-ion collisions and
examine the ensuing produced particles. Amongst all, measurements of light
(anti-)nuclei and hyper (anti-)nuclei production are of great importance to
understand the nature of hadron-hadron interactions and extract the freeze-out
parameters of the QCD matter produced in these collisions. These studies also
lead towards unfolding the mysteries of matter-antimatter asymmetry and
provide relevant insights in the search for signatures of dark matter
Winkler:2020ltd . In hadronic or heavy-ion collision experiments,
(anti-)nuclei production mechanisms are not well understood. They are believed
to be produced in the final stages of the hadronic evolution, specifically
during chemical and kinetic freeze-out. The typical binding energies of light
nuclei are of $\mathcal{O}(2\ \rm{MeV})$, which is much smaller than the
temperatures during chemical and kinetic freeze-out $\mathcal{O}(100\
\rm{MeV})$.

A substantial amount of studies with theoretical models have described how the
light (anti-)nuclei and hyper (anti-)nuclei are produced in these extreme
conditions. Models such as the thermal model, which consider the system to be
in chemical and kinetic equilibrium, are successful in reproducing the anti-
baryon over baryon and light anti-nuclei over nuclei ratios Andronic:2010qu ;
Cleymans:2011pe ; Steinheimer:2012tb . On the other hand, coalescence models;
based on the ansatz that the formation of light (anti-)nuclei and hyper
(anti-)nuclei occur after the system has reached a kinetic freeze-out
successfully describe nuclei production and the respective particle ratio
measurements from experiments Shah:2015oha ; Zhao:2021dka . In a coalescence
approach, nucleons come close in phase space to form nuclei. The coalescence
parameter ($B_{A}$) represents the probability of production via coalescing
nucleons at kinetic freezeout, and is inversely related to the nucleon
correlation volume (volume of nuclear matter during the time coalescence of
nucleons to nuclei). Often the coalescence mechanism is applied as an
afterburner, after being driven by transport or hydrodynamic model for phase
space distribution of nucleons Sun:2020uoj ; Hillmann:2021zgj ; Zhao:2020irc ;
Zhu:2015voa . Moreover, kinetic processes or potential interactions during the
hadronic evolution also serve as a tangible mechanism for nuclei production.
Recent advancements in n-body transport models Aichelin:1991xy are allowing
detailed and concrete studies in light as well as hyper nuclei cluster
formation via potential interactions Glassel:2021rod .

The first measurements in light (anti-)nuclei production were carried out at
the CERN Inner Storage Ring experiments Alper:1973my ; British-Scandinavian-
MIT:1977tan . Several measurements that followed in RHIC and the LHC reported
light nuclei and hyper nuclei ($d,\bar{d}$, $\rm{{}^{3}H}$,
$\rm{{}^{4}\overline{He}}$, $\rm{{}^{4}_{\Lambda}\overline{He}}$) production
in $p+p$, Au+Au, $p$+Pb and Pb+Pb systems in different center of mass energies
STAR:2019sjh ; ALICE:2015wav ; ALICE:2015oer ; STAR:2011eej ; STAR:2022hbp .
All these studies open a lot of doors to test the model calculations. The
dependence of nuclei production with energy allow researchers to narrow down
the location of the critical point in the QCD phase diagram. The coalescence
parameter ($B_{A}$) is reported in these measurements to have a finite slope
as a function of $p_{\rm{T}}$. This is regarded to as an increase in
coalescence probability, or a decrease in the nucleon correlation volume with
increasing $p_{\rm{T}}$.

This article presents the implementation of a simple transport + coalescence
methodology for producing light nuclei. The coalescence model, which requires
only a small number of tuning parameters, effectively captures light nuclei
($d,\bar{d}$, $\rm{{}^{3}H}$, $\rm{{}^{3}He}$) production for a broad range of
beam energies in heavy-ion collisions. In section II, the AMPT model
Lin:2004en  is described, with its use to obtain the phase space information
of the constituent nucleons. We employ the string melting version of AMPT
(version 1.26t9b-v2.26t9b), which is then coupled to a coalescence
afterburner, from which the light nuclei yields are estimated. The coalescence
model and its implementation is discussed in section III. AMPT also employs
kinetic processes for deuteron and anti-deuteron formation during the hadronic
transport, which gives an opportunity for a comparative study between the two
modes (transport vs coalescence). The study is carried out in various center
of mass energies $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, and 200 GeV for central
(0-10%) Au+Au collisions and $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV for central
(0-10% and 0-20%) Pb+Pb collisions at mid-rapidity. We also compare our
results with some recent predictions from the parton-hadron-quantum-molecular-
dynamics (PHQMD) model Aichelin:2019tnk  with light nuclei clustering
algorithms Glassel:2021rod . The nuclei clustering is done dynamically in
various cluster freezeout times via Simulated Annealing Clustering Algorithm
(SACA) and Minimum Spanning Tree (MST) algorithm Puri:1998te ;
Aichelin:2019tnk . In section IV, the results are presented and compared to
the existing experimental measurements from the STAR and ALICE experiments.

##  II The AMPT model

AMPT Lin:2004en  consists of a hybrid framework of four primary components;
the initial conditions, parton level interactions, a hadronization mechanism,
and the hadronic scatterings. The initial conditions comprise the hard
scatterings between mini-jet partons and soft scattering strings from the
Heavy-Ion Jet INteraction Generator (HIJING) model Wang:1991hta . Zhang Parton
Cascade (ZPC) model Zhang:1997ej  handles the parton evolution and the
scatterings between them, which recombine to form hadrons via the Lund string
fragmentation model Andersson:1983ia . However, in the string melting
Lin:2001zk  version of AMPT, the ZPC handles both minijet partons and the
excited strings during parton evolution, which then goes through hadronization
via a parton coalescence mechanism Molnar:2003ff . The final stage is the
hadronic evolution, which is brought in by A Relativistic Transport (ART)
model Li:1995pra  and additional resonance channels that contribute to the
production of final state hadrons.

The AMPT model is well established in describing various experimental
observables obtained at RHIC and the LHC STAR:2010vob ; ALICE:2013rdo ;
STAR:2013ayu . Several interpretations of nuclei formation are covered using
AMPT. AMPT natively produces (anti-)deuterons via elementary interactions of
$NN\rightarrow d\pi$, $NNN\rightarrow dN$, and $NN\pi\rightarrow d\pi$ during
the hadronic evolution. In this study, we investigate the deuteron production
processes in two different modes. For $``$Coalescence” mode, we turn off the
native deuteron processes in AMPT to ensure that all (anti-)nucleon pairs are
available for the coalescence afterburner. For the $``$Transport” mode, we
explicitly turn on the hadronic rescattering processes and do not implement
the afterburner on this freeze-out hypersurface.
We generate AMPT events for the phase space information of (anti-)nucleons
which are fed to the coalescence afterburner to produce (anti-)deuterons,
tritons and helium-3s. We generate events for central Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, and 200 GeV and central Pb+Pb collisions
at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV. For the Pb+Pb system, an additional
parton mini-jet momentum cutoff at 5 GeV/$c$ is chosen which minimizes the
softening of the (anti-)proton $p_{\rm{T}}$ spectrum. The ALICE results for
Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV is well described
with the application of an additional cut-off on the parton mini-jet momentum
Zhang:2019utb . The Lund String Fragmentation function, which is defined as

$f(z)\propto z^{-1}(1-z)^{a}\exp(-bm_{\perp}^{2}/z)$ (1)

where, $z$ denotes the light cone momentum fraction of a produced hadron with
respect to that of the fragmentation string, $m_{\perp}$ is the transverse
mass of the hadron, and $a$ and $b$ are the Lund Fragmentation parameters. It
can be shown that the mean squared transverse momentum of the particles is
proportional to the string tension ($k$) of fragmentation string ($\langle
p_{\rm T}^{2}\rangle\propto 1/b(2+a)$), which is related to the parameters $a$
and $b$ by

$k\propto\frac{1}{b(2+a)}$ (2)

The parameters ($a$ and $b$) which best describe the experimental data for
protons AMPT are chosen. The same set of parameters are used for the phase
space information of (anti-)nucleons, in order to coalesce them into the light
nuclei species. The parton scattering cross-section is kept as 3$~{}mb$
Tiwari:2020ult . The centrality of the collisions, 0-10% and 0-20% are taken
into account by considering the collision impact parameter range within 0-4
$fm$ and 0-7 $fm$ respectively. A detailed description of $a$ and $b$ values
for $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, 200 GeV and 2.76 TeV is reported in
Table 1. A total of $\mathcal{O}(10^{5})$ AMPT events are generated for each
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, 200 GeV (Au+Au) and 2.76 TeV (Pb+Pb)
corresponding to both centrality classes.

$\sqrt{s_{\mathrm{NN}}}~{}(\rm{GeV})$ | $a$ | $b~{}(\rm{GeV^{-2}})$
---|---|---
19.6 | 0.55 | 0.15
39 | 0.55 | 0.15
200 | 0.55 | 0.15
2760 | 0.8 | 0.05
Table 1: Parameters ($a$ and $b$) for Lund String Fragmentation function
chosen for Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, 200 GeV
and Pb+Pb collisions $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV corresponding to
0-10% and 0-20% centrality class.

It is worth noting that the present iteration of the string melting AMPT has
certain limitations in precisely characterizing baryons He:2017tla , a key
aspect with regard to the findings presented in this paper. Notably, AMPT
without transport deuterons (under)overestimates the yields of (anti-)protons
at mid-rapidity. A possible explanation in this regard will be mentioned in
Section IV. Additionally, the model implements quark coalescence via separate
conservation of the numbers of meson and (anti-)baryons per event. This
constraint is removed with an upcoming string melting AMPT by assigning a
control parameter for the quarks to form a baryon or a meson, during quark
coalescence. Studies involving the development of an updated string melting
AMPT with the improved quark coalescence prescription cites a better
description of experimental measurements for baryons He:2017tla .

Parameters | (Anti-)Deuteron | Triton/Helium-3
---|---|---
$\Delta r_{max}~{}(fm)$ | 3.15 | 4.18
$\Delta p_{max}$ (MeV/$c$) | 190 | 350
Table 2: Parameters ($\Delta p_{max}$ and $\Delta r_{max}$) used in
coalescence of deuterons, tritons and helium-3.

##  III The Coalescence Model

The process of light (anti-)nuclei and hyper (anti-)nuclei production in
heavy-ion collisions can be explained by the coalescence mechanism of
respective nucleons or baryons Steinheimer:2012tb . The model assumes that
nuclei production begins at kinetic freeze-out during the later stages of a
collision. In the limit of no-baryon number transport at mid-rapidity, the
coalescence probability is described by the coalescence parameter $`B_{A}$’,
as

$\displaystyle E_{A}\frac{d^{3}N_{A}}{d^{3}p_{A}}$ $\displaystyle=$
$\displaystyle
B_{A}\left(E_{p}\frac{d^{3}N_{p}}{d^{3}p_{p}}\right)^{Z}\left(E_{n}\frac{d^{3}N_{n}}{d^{3}p_{n}}\right)^{A-Z}$
$\displaystyle\approx$ $\displaystyle
B_{A}\left(E_{p}\frac{d^{3}N_{p}}{d^{3}p_{p}}\right)^{A}$

where $E\frac{d^{3}N}{d^{3}p}$ is the invariant yield, $A$ and $Z$ are the
atomic mass number and the atomic number of (anti-)nucleons or (anti-)nuclei.
To calculate $B_{A}$, the proton and light nuclei transverse momentum are
related as $p_{\mathrm{T}A}=A\cdot p_{\mathrm{T}p}$. The coalescence
conditions for two nucleons to form a $d,\bar{d}$, $\rm{{}^{3}H}$,
$\rm{{}^{3}He}$ are characterized by their phase space distribution. In this
case, we implement a box description of coalescence where we define conditions
based on the relative distance $(\Delta r)$ and momentum $(\Delta p)$ of the
nucleon pairs. The coalescence conditions are met when the relative momentum
between the nucleons are less than a critical value ($\Delta p^{max}$), and
the distance between the nucleons is less than twice the size of the nuclear
force radius ($\Delta r^{max}<2R_{0}$). $R_{0}$ is taken to be $\sim$1.6 fm
for deuterons and $\sim$2.1 fm for triton and helium-3 Hillmann:2021zgj . A
summary of the coalescence parameters used for (anti-)deuteron, triton and
helium-3 is tabulated in Table 2. A step by step procedure to achieve this
phase space coalescence is described below:

  1. 1.

The (anti-)nucleon states, which is $pn(\bar{p}\bar{n})$ for $d(\bar{d})$,
$pnn$ for $\rm{{}^{3}{H}}$ and $ppn$ for $\rm{{}^{3}{He}}$ are boosted to the
center of mass frame.

  2. 2.

If their relative momentum $\Delta p=\mid\vec{p}_{1}-\vec{p}_{2}\mid<\Delta
p^{max}$ , their momenta are combined ($\vec{p}_{d}=\vec{p}_{p}+\vec{p}_{n}$
or
$\vec{p}_{\rm{{}^{3}H},\rm{{}^{3}He}}=\vec{p}_{p}+\vec{p}_{n,p}+\vec{p}_{n}$)

  3. 3.

Among the (anti-)nucleon pairs, the position(s) of the (anti-)nucleon(s) at an
earlier time are reset to the (anti-)nucleon at the later freeze-out time.
This is done to ensure that the pair (triplet) of nucleons do not interact
between these two times.

  4. 4.

The relative distances between the pairs are checked so that the condition
$\Delta r=\mid\vec{r}_{1}-\vec{r}_{2}\mid<\Delta r^{max}$ is satisfied.

  5. 5.

The spin-isospin coupling probabilities are taken into account (3/8 for a
deuteron and 1/12 for triton and helium-3) Hillmann:2021zgj ; Xia:2014rua .
The finally chosen deuteron, triton and helium-3 candidates that fulfill the
coalescence conditions are then removed from the phase space distribution.

##  IV Result and Discussion

Figure 1: (Color online) Transverse momentum spectra of protons (left) and
anti-protons (right) in central Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$
= 19.6 GeV, 39 GeV ($|y|<0.1$), 200 GeV ($|y|<0.5$) and central Pb+Pb
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV ($|y|<0.5$) from
$``$Coalescence” or $``$Transport OFF” (solid lines) and $``$Transport ON”
(dashed lines) mode. Experimental results (symbols) are taken from STAR and
ALICE experiments STAR:2007zea ; STAR:2017sal ; ALICE:2013mez . The lower
panels show the ratio of the experimental results to the model predictions.

Figure 2: (Color online) Transverse momentum spectra of deuterons (left) and
anti-deuterons (right) in central Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 GeV, 39 GeV, 200 GeV ($|y|<0.3$) and central
Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV ($|y|<0.5$) from
$``$Coalescence” (solid lines) and $``$Transport” (dashed lines) mode.
Experimental results (symbols) are taken from the STAR and ALICE measurements
STAR:2019sjh ; ALICE:2015wav . The lower panels show the ratio of the
experimental results to the model predictions. Figure 3: (Color online)
Transverse momentum spectra of tritons in central Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 GeV, 39 GeV, 200 GeV ($|y|<0.5$) and
helium-3s in central Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV
($|y|<0.5$) from the $``$Coalescence” mode (solid lines). Experimental results
(symbols) are taken from the STAR and ALICE measurements STAR:2022hbp ;
ALICE:2015wav . The lower panels show the ratio of the experimental results to
the model predictions. Figure 4: (Color Online) Particle ratios for anti-
deuterons over deuterons as a function of center of mass energies from
$``$Coalescence” and $``$Transport” modes. Experimental measurment of
$\bar{d}/d$ ratios from the STAR and ALICE experiments for the central Au+Au
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 GeV, 39 GeV, 200 GeV and
central Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV are also
shown STAR:2019sjh ALICE:2015wav

The transverse momentum spectra ($p_{\rm{T}}$) for the (anti-)proton and their
comparisons to the STAR and ALICE measurements STAR:2007zea ; STAR:2017sal ;
ALICE:2013mez  at mid-rapidity for $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39
($|y|<0.1$), 200 GeV ($|y|<0.5$) (Au+Au, 0-10%) and 2.76 TeV (Pb+Pb, 0-10%,
$|y|<0.5$) respectively, are shown in Figure 1. The $``$Coalescence”
(transport deuterons are OFF) results overestimate the data for protons at
$p_{\rm{T}}$ $<1$ GeV/$c$ for $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV and
underestimate at $p_{\rm{T}}$ $\sim$ 2 GeV/$c$. In this aspect, the
$``$Transport” (transport deuterons are ON) mode provides a better description
of the $p_{\rm{T}}$ spectra for $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV
STAR:2017sal . Both the modes show good agreement with the data for protons at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au) and $\sqrt{s_{{}_{\rm{NN}}}}$ =
2.76 TeV (Pb+Pb) STAR:2007zea ; ALICE:2013mez . However, the $``$Coalescence”
mode largely underestimates the data for Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV, but describe the experimental
measurements at $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV and 2.76 TeV. The anti-
protons from $``$Transport” mode provides a slightly harder $p_{\rm{T}}$
spectra for $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV when compared to data,
and performs similar to the $``$Coalescence” mode for
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV and 2.76 TeV.

Figure 2 displays the (anti-)deuteron $p_{\rm{T}}$ spectra for
$``$Coalescence” and $``$Transport” modes, which are compared to experimental
measurements at midrapidity for Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$
= 19.6, 39, 200 GeV ($|y|<0.3$), as well as Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV ($|y|<0.5$). While the deuterons produced
via $``$Coalescence” accurately reproduce the experimental data at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au) STAR:2019sjh  and 2.76 TeV (Pb+Pb)
ALICE:2015wav  at all $p_{\rm{T}}$ values, they tend to overestimate the data
at low $p_{\rm{T}}$ and underestimate it at high $p_{\rm{T}}$ in Au+Au
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV STAR:2019sjh . This
disagreement is attributed to the softer proton $p_{\rm{T}}$ spectra from
AMPT. The $``$Coalescence” anti-deuteron spectra show minor deviations across
$p_{\rm{T}}$, and fit well with the experimental measurements for Au+Au and
Pb+Pb collisions. STAR:2019sjh ; ALICE:2015wav . In contrast, the deuterons
and anti-deuteron yields from $``$Transport” mode significantly deviates from
the experimental results. In particular, the anti-deuteron yields are largely
overestimated for Au+Au collisions. While, the $``$Transport” deuteron spectra
for Au+Au collision at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV sits close
to the $``$Coalescence” spectra, the yields are overestimated for
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV. In Pb+Pb collisions, the $``$Transport”
(anti-)deuteron yields are underestimated, with a visible drop at low
$p_{\rm{T}}$.

The results shown in Figure 1 and 2 clearly demonstrate a trade-off between
(anti-)proton and (anti-)deuteron production in the AMPT model. This trade-off
likely contributes to the observed differences in (anti-)proton yields between
the $``$Coalescence” and $``$Transport” modes, especially for anti-proton
yields in Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 GeV and 39 GeV.
The overestimation of anti-deuteron yields when using the $``$Transport” mode
suggests that there might be an inherent imbalance in the cross-sections
governing anti-deuteron production and dissociation processes. This imbalance
directly impacts the anti-proton yields when these processes are turned OFF.
In the $``$Coalescence” mode, where these kinetic processes are absent, the
production of anti-protons is underestimated, particularly at collision
energies where the net proton number is non-zero.

The triton and helium-3 $p_{\rm{T}}$-spectra results from $``$Coalescence”
mode are presented and compared to experimental results in Figure 3 for
central (0-10%) Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, 200
GeV STAR:2022hbp  and central (0-20%) Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV ALICE:2015wav  at mid-rapidity
($|y|<0.5$). For tritons in Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ =
19.6 and 39 GeV, the model exhibits a deviation across $p_{\rm{T}}$ that
follows a similar pattern to that observed for deuterons. The model slightly
overestimates the experimental data at low $p_{\rm{T}}$ and underestimates at
larger values. However, the behaviour changes for tritons in Au+Au collisions
at $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV, giving a harder spectra when compared
to the data. The helium-3 spectra for Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV is well described, with a slight
overestimation for $p_{\rm{T}}$ $<3.5$ GeV/$c$.

The $p_{\rm{T}}$-integrated yields for $d,\bar{d}$ and $\rm{{}^{3}H}$,
$\rm{{}^{3}He}$ are calculated for central Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39, 200 GeV and Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV at mid-rapidity via $``$Coalescence” and
$``$Transport” modes. Figure 4 shows the $\bar{d}/d$ ratios as a function of
center of mass energy compared to experimental measurements. Both modes are
able to qualitatively describe the trend of energy dependence of $d$ and
$\bar{d}$ production successfully at the respective RHIC and LHC
energiesSTAR:2019sjh . Quantitatively, the $``$Transport” mode systematically
overestimates the experimental data for Au+Au collisions from 19.6 to 200 GeV,
but largely underestimates the value for Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV. The $``$Coalescence” mode underestimates
the ratios for Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39
GeV, and slightly overestimates for $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV. The
value at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV is slightly underestimated.

Figure 5: (Color Online) Particle ratios for deuterons over protons and anti-
deuterons over anti-protons as a function of center of mass energies from
$``$Coalescence” and $``$Transport” modes. Solid lines are predictions from
the PHQMD model Glassel:2021rod . Experimental measurement of $d/p$ and
$\bar{d}/\bar{p}$ ratios from the STAR and ALICE experiments for the central
Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 GeV, 39 GeV, 200 GeV and
central Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV are also
shown STAR:2019sjh ; ALICE:2015wav .

The $d(\bar{d})$ over $p(\bar{p})$ ratios are presented in Figure 5 as a
function of center of mass energy. Both modes capture the trends for $d/p$ and
$\bar{d}/\bar{p}$ in Au+Au collisions. However, there are substantial
quantitative discrepancies between the model predictions and the experimental
data. In particular, the $``$Transport” mode significantly overestimates the
yield ratios compared to the experimental measurements, indicating a large
deviation from the data STAR:2019sjh . The $``$Coalescence” mode overestimates
the values, especially at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV, but
gives a good description for $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV and 2.76 TeV
ALICE:2015wav . For completeness, we append additional model comparisons from
Glassel:2021rod  in Figure 5. The predictions from the light nuclei cluster
production with a PHQMD approach shows a comprehensive study with different
cluster freezeout times or $``$physical” times. The $d/p$ ratio is well
described for t $=60$ fm/$c$ for Au+Au collisions. The $\bar{d}/\bar{p}$
ratios for t $=70$ fm/$c$ are closer to data and consistent to our
$``$Coalescence” predictions. However, this is accredited to an over-predicted
$\bar{d}$ yield by the model.

In Figure 6, the $\rm{{}^{3}H}$, $\rm{{}^{3}He}$ $p_{\rm{T}}$-integrated
yields over $p$ and $d$ yields for central (0-10%) Au+Au STAR:2022hbp  and
Pb+Pb collisions with $``$Coalescence” mode are presented ALICE:2015wav . All
the ratios show the dependence with center of mass energy, and capture the
trend of the experimental results. The $\rm{{}^{3}H}/p$ ratios are slightly
overestimated for Au+Au collisions, but the $\rm{{}^{3}He}/p$ ratio for Pb+Pb
collisions is well described. The $\rm{{}^{3}H}/d$ and $\rm{{}^{3}He}/d$
values at 0-10% are also closely compatible with the experimental measurements
STAR:2022hbp ; ALICE:2015wav .

Figure 6: (Color Online) Particle ratios from the $``$Coalescence” mode for
triton over proton and deuteron for 0-10% central Au+Au and Pb+Pb collisions,
helium-3 over deuteron for 0-10% central Pb+Pb collisions and helium-3 over
proton for 0-20% central Pb+Pb collisions. Experimental measurement of
${}^{3}{H}/p$, ${}^{3}{H}/d$, ${}^{3}{He}/p$, ${}^{3}{He}/d$ ratios from the
STAR and ALICE experiments for the central Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 GeV, 39 GeV, 200 GeV and central Pb+Pb
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV are also shown STAR:2022hbp
; ALICE:2015wav .

It can be concluded that while there are some discrepancies between the
model’s $p_{\rm{T}}$-differential and $p_{\rm{T}}$-integrated yields of light
nuclei and experimental measurements from a quantitative standpoint, the model
is still able to provide a reliable qualitative description of light nuclei
yields as a function of $p_{\rm{T}}$ and center of mass energy. These findings
demonstrate the effectiveness of a simple phase space coalescence approach in
qualitatively describing light nuclei production across a wide range of
collision energies and systems. The model’s ability to offer a comprehensive
explanation of experimental data with only a few tuning parameters is
particularly noteworthy. Moving forward, these results can be further refined
by leveraging the enhanced capabilities of the updated string melting AMPT
model.

The $B_{2}$ values as a function of $p_{\rm{T}}$ are calculated for central
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 39 GeV (Au+Au) and 2.76 TeV (Pb+Pb)
using equation III. Figure 7 (left) shows the results from $``$Coalescence”
and $``$Transport” modes compared to experimental measurements STAR:2019sjh
ALICE:2015wav . The $``$Coalescence” mode is promising in reproducing the
increasing $p_{\rm{T}}$ dependence of deuteron $B_{2}$ for central Au+Au
collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 39 GeV and 2.76 TeV. For deuterons,
the model predictions are compatible with the experimental results for Au+Au
as well as Pb+Pb collisions. The anti-deuterons at $\sqrt{s_{{}_{\rm{NN}}}}$ =
39 GeV (Au+Au) shows an overestimated $B_{2}$ vs $p_{\rm{T}}$ with a finite
slope. The $``$Transport” mode is unable to reproduce the increasing
dependence of $B_{2}$ with $p_{\rm{T}}$ for Au+Au collisions. Both deuteron
and anti-deuteron results are largely overestimated by $``$Transport” mode for
$\sqrt{s_{{}_{\rm{NN}}}}$ = 39 GeV (Au+Au). For Pb+Pb collisions, the
$``$Transport” mode is able to reproduce the increasing trend with
$p_{\rm{T}}$. However, the results are underestimated in comparison to the
experimental measurementsALICE:2015wav . The highlight of these results is the
ability of the $``$Coalescence” mode to reproduce the increasing trend of
$B_{2}$ with $p_{\rm{T}}$. The PHQMD model results for deuteron $B_{2}$ for
Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 39 GeV also successfully
reproduces the increasing trend with $p_{\rm{T}}$, for cluster freezeout time
t $=70$ fm/$c$. A lower freezeout time overestimates the values. Moreover, the
simple coalescence model outperforms in quantitative estimation of the $B_{2}$
for deuterons.

Furthermore, we calculate the $B_{3}$ values as a function of $p_{\rm{T}}$ for
central collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au) and 2.76 TeV
(Pb+Pb), and we present the results from the $``$Coalescence” mode in Figure 7
(right). The predictions are fairly close to the experimental data for 0-20%
central Pb+Pb collisions ALICE:2015wav  and 0-10% central Au+Au collisions
Zhang:2020ewj , respectively. The trend of $B_{3}$ as a function of
$p_{\rm{T}}$ is quantitatively reproduced for helium-3s with good agreement
with the experimental results. On the other hand, the data for triton $B_{3}$
shows a flat or decreasing dependence with $p_{\rm{T}}$, whereas the model
shows somewhat of an increasing trend.

Figure 8 reports the $B_{A}$ as a function of center of mass energies for
(anti-)deuterons, tritons and helium-3 for central Au+Au and Pb+Pb collisions.
The values are obtained for $p_{\rm{T}}/A$ = 0.65 GeV/$c$, which is 1.3
GeV/$c$ for deuterons and 1.95 GeV/$c$ for triton and helium-3. The $B_{2}$ of
deuterons from $``$Coalescence” shows a good match within uncertainties with
the experimental results STAR:2019sjh  across the $\sqrt{s_{{}_{\rm{NN}}}}$.
The points for anti-deuterons at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV
are systematically overestimated and appear higher than the $B_{2}$ for
deuteron. For $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au), the anti-deuteron
$B_{2}$ falls within the uncertainties of the experimental results, nearly
superimposed on the $B_{2}$ of the deuteron. The $``$Transport” mode,
overestimates the coalescence parameter for both deuteron and anti-deuteron
for the Au+Au collisions. The same is underestimated for Pb+Pb collisions. The
PHQMD predictions also show the deuteron $B_{2}$ as a function of collision
energy for Au+Au collisions. Out of the two freezeout times, t $=70$ fm/$c$
shows a better agreement to the experimental measurements. The $B_{3}$ for
tritons in Au+Au collisions are also shown and compared to experimental
measurements from STAR Zhang:2020ewj . The results at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV is accurately described by the model,
whereas at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV the values are slightly
underestimated.

Figure 7: (Color Online) (Left) The $B_{2}$ parameter versus $p_{\rm{T}}$ of
deuterons and anti-deuterons in central (0-10%) Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 39 GeV and central Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV from $``$Coalescence” (solid lines) and
$``$Transport” (dashed lines) modes. Predictions from the PHQMD model are
superimposed Glassel:2021rod . (Right) The $B_{3}$ parameter versus
$p_{\rm{T}}$ for triton and helium-3 in central (0-10%) Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV and central (0-20%) Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV from $``$Coalescence” (solid lines) mode.
The results are compared to experimental measurements (symbols) from STAR and
ALICE STAR:2019sjh ; Zhang:2020ewj ; ALICE:2015wav . Figure 8: Coalescence
parameter $B_{2}$ and $\sqrt{B_{3}}$ as a function of
$\sqrt{s_{{}_{\rm{NN}}}}$ in central (0-10%) Au+Au and Pb+Pb collisions for
$``$Coalescence” and $``$Transport” modes. Solid lines are predictions from
the PHQMD model Glassel:2021rod . The results are compared with measurements
from the STAR and ALICE experiments for $p_{\rm{T}}/A$ = 0.65 GeV/$c$.
STAR:2019sjh ; Zhang:2020ewj ; ALICE:2015wav .

##  V Conclusions

This article presents the results of a systematic study on the production of
light nuclei in central Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6,
39, 200 GeV, as well as Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76
TeV using the AMPT driven coalescence afterburner model. The coalescence model
successfully provides a comprehensive description of all light nuclei species
across various collision energies in heavy-ion collisions. We also test the
feasibility of (anti-)deuteron production processes implemented during the
hadronic evolution in AMPT. Furthermore, we make a comparison of nuclei
cluster production with a PHQMD approach with our calculations. The PHQMD
approach considers deuteron production via potential interactions during the
hadronic evolution. Overall, our phase space coalescence model out-performs
the kinetic implementation of deuteron production in AMPT and performs closely
with the PHQMD results.

The string melting AMPT model under the $``$Coalescence” mode struggles to
quantitatively reproduce the experimental measurements of (anti-)proton
$p_{\rm{T}}$ spectra for Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6,
39 GeV but provides a good description for Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV and Pb+Pb collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV. The $``$Transport” mode shows a better
agreement of the (anti-)proton $p_{\rm{T}}$ spectra in all collision energies.
While the predictions from $``$Coalescence” mode deviate from the data for the
$p_{\rm{T}}$ differential yields of deuterons for Au+Au collisions at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6, 39 GeV, it gives a reasonable description
for the deuterons in central Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ =
200 GeV and Pb+Pb collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV. A
similar outcome is seen for tritons and helium-3s, with a slight deviation for
Au+Au collisions at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV, with a good
description at $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au) and 2.76 TeV
(Pb+Pb). The $``$Transport” mode largely overestimate the deuteron and anti-
deuteron yields as a function of $p_{\rm{T}}$ across all energies.

The $``$Coalescence” mode qualitatively describes the various particle ratios
as a function of collision energy. The trend of the yields with center of mass
energies is successfully described by the model predictions. The
$``$Coalescence” results are fairly compatible with the $d/p$,
$\bar{d}/\bar{p}$ ratios predicted by the PHQMD model. The $``$Transport”
model overestimates these ratios. The $\bar{d}/d$ ratios are qualitatively
reproduced by both $``$Coalescence” and $``$Transport” mode for Au+Au
collisions.

The coalescence parameter $B_{2}$ and $B_{3}$ as a function of $p_{\rm{T}}$
are estimated for RHIC and LHC energies. The results from $``$Coalescence”
mode describe the trends of the coalescence parameter fairly well, reproducing
the increasing dependence of $B_{2}$ with $p_{\rm{T}}$ from the experimental
measurements. The results are better when compared to the PHQMD model for
deuteron $B_{2}$. The $B_{3}$ vs $p_{\rm{T}}$ for helium-3 at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 2.76 TeV (Pb+Pb) is well described by the model,
however struggles to reproduce the trend for tritons at
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au). The description of the
$p_{\rm{T}}$-spectra of (anti-)protons and light-nuclei by the model may vary
at low and high $p_{\rm{T}}$, but these effects are diminished for the
predictions of $B_{A}$. The coalescence parameter $B_{A}$ as a function of
$\sqrt{s_{{}_{\rm{NN}}}}$ is also well described by the $``$Coalescence” mode
for the deuterons for the RHIC and LHC energies. The $``$Coalescence” mode
overestimates the anti-deuteron $B_{2}$ at $\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6
and 39 GeV, but sits well for $\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV. The $B_{3}$
values for triton and helium-3 are well reproduced for
$\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV, and underestimates for
$\sqrt{s_{{}_{\rm{NN}}}}$ = 19.6 and 39 GeV.

In conclusion, the results show the effectiveness of a simple coalescence
model and its ability to describe light nuclei production across different
center of mass energies in heavy-ion collisions. We broadly observe that the
model provides better predictions for higher collision energies
($\sqrt{s_{{}_{\rm{NN}}}}$ = 200 GeV (Au+Au) and 2.76 TeV (Pb+Pb)). The model
successfully covers various aspects of light nuclei production while
accurately predicting a wide range of experimental observables. This study can
be extended to hyper-nuclei production, as well as exotic states that are
measured in heavy-ion experiments.

##  VI Acknowledgements

Author Y.B is thankful to Z.W. Lin for his helpful suggestions to this work,
and to all the authors of AMPT. Y. B is grateful to Sudhir P. Rode for
carefully reading the manuscript, and to Sumit Kundu for his contribution to
the data generating process. Y.B is also thankful to Ravindra Singh and
Swapnesh Khade for the fruitful discussions.

## References

  * (1) J. Adams et al. [STAR], Nucl. Phys. A 757 (2005), 102-183 doi:10.1016/j.nuclphysa.2005.03.085 [arXiv:nucl-ex/0501009 [nucl-ex]].
  * (2) K. Aamodt et al. [ALICE], JINST 3 (2008), S08002 doi:10.1088/1748-0221/3/08/S08002
  * (3) M. W. Winkler and T. Linden, Phys. Rev. Lett. 126 (2021) no.10, 101101 doi:10.1103/PhysRevLett.126.101101 [arXiv:2006.16251 [hep-ph]].
  * (4) A. Andronic, P. Braun-Munzinger, J. Stachel and H. Stocker, Phys. Lett. B 697 (2011), 203-207 doi:10.1016/j.physletb.2011.01.053 [arXiv:1010.2995 [nucl-th]].
  * (5) J. Cleymans, S. Kabana, I. Kraus, H. Oeschler, K. Redlich and N. Sharma, Phys. Rev. C 84 (2011), 054916 doi:10.1103/PhysRevC.84.054916 [arXiv:1105.3719 [hep-ph]].
  * (6) J. Steinheimer, K. Gudima, A. Botvina, I. Mishustin, M. Bleicher and H. Stocker, Phys. Lett. B 714 (2012), 85-91 doi:10.1016/j.physletb.2012.06.069 [arXiv:1203.2547 [nucl-th]].
  * (7) N. Shah, Y. G. Ma, J. H. Chen and S. Zhang, Phys. Lett. B 754 (2016), 6-10 doi:10.1016/j.physletb.2016.01.005 [arXiv:1511.08266 [nucl-ex]].
  * (8) W. Zhao, K. j. Sun, C. M. Ko and X. Luo, Phys. Lett. B 820 (2021), 136571 doi:10.1016/j.physletb.2021.136571 [arXiv:2105.14204 [nucl-th]].
  * (9) K. J. Sun and C. M. Ko, Phys. Rev. C 103 (2021) no.6, 064909 doi:10.1103/PhysRevC.103.064909 [arXiv:2005.00182 [nucl-th]].
  * (10) P. Hillmann, K. Käfer, J. Steinheimer, V. Vovchenko and M. Bleicher, J. Phys. G 49 (2022) no.5, 055107 doi:10.1088/1361-6471/ac5dfc [arXiv:2109.05972 [hep-ph]].
  * (11) W. Zhao, C. Shen, C. M. Ko, Q. Liu and H. Song, Phys. Rev. C 102 (2020) no.4, 044912 doi:10.1103/PhysRevC.102.044912 [arXiv:2009.06959 [nucl-th]].
  * (12) L. Zhu, C. M. Ko and X. Yin, Phys. Rev. C 92 (2015) no.6, 064911 doi:10.1103/PhysRevC.92.064911 [arXiv:1510.03568 [nucl-th]].
  * (13) J. Aichelin, Phys. Rept. 202 (1991), 233-360 doi:10.1016/0370-1573(91)90094-3
  * (14) S. Gläßel, V. Kireyeu, V. Voronyuk, J. Aichelin, C. Blume, E. Bratkovskaya, G. Coci, V. Kolesnikov and M. Winn, Phys. Rev. C 105 (2022) no.1, 014908 doi:10.1103/PhysRevC.105.014908 [arXiv:2106.14839 [nucl-th]].
  * (15) B. Alper, H. B$\gimel$gild, P. Booth, F. Bulos, L. J. Carroll, G. von Dardel, G. Damgaard, B. Duff, F. Heymann and J. N. Jackson, et al. Phys. Lett. B 46 (1973), 265-268 doi:10.1016/0370-2693(73)90700-4
  * (16) S. Henning et al. [British-Scandinavian-MIT], Lett. Nuovo Cim. 21 (1978), 189 doi:10.1007/BF02822248
  * (17) J. Adam et al. [STAR], Phys. Rev. C 99 (2019) no.6, 064905 doi:10.1103/PhysRevC.99.064905 [arXiv:1903.11778 [nucl-ex]].
  * (18) J. Adam et al. [ALICE], Phys. Rev. C 93 (2016) no.2, 024917 doi:10.1103/PhysRevC.93.024917 [arXiv:1506.08951 [nucl-ex]].
  * (19) J. Adam et al. [ALICE], Phys. Lett. B 754 (2016), 360-372 doi:10.1016/j.physletb.2016.01.040 [arXiv:1506.08453 [nucl-ex]].
  * (20) H. Agakishiev et al. [STAR], Nature 473 (2011), 353 [erratum: Nature 475 (2011), 412] doi:10.1038/nature10079 [arXiv:1103.3312 [nucl-ex]].
  * (21) [STAR], [arXiv:2209.08058 [nucl-ex]].
  * (22) J. Aichelin, E. Bratkovskaya, A. Le Fèvre, V. Kireyeu, V. Kolesnikov, Y. Leifels, V. Voronyuk and G. Coci, Phys. Rev. C 101 (2020) no.4, 044905 doi:10.1103/PhysRevC.101.044905 [arXiv:1907.03860 [nucl-th]].
  * (23) R. K. Puri and J. Aichelin, J. Comput. Phys. 162 (2000), 245-266 doi:10.1006/jcph.2000.6534 [arXiv:nucl-th/9811018 [nucl-th]].
  * (24) Z. W. Lin, C. M. Ko, B. A. Li, B. Zhang and S. Pal, Phys. Rev. C 72 (2005), 064901 doi:10.1103/PhysRevC.72.064901 [arXiv:nucl-th/0411110 [nucl-th]].
  * (25) X. N. Wang and M. Gyulassy, Phys. Rev. D 44 (1991), 3501-3516 doi:10.1103/PhysRevD.44.3501
  * (26) B. Zhang, Comput. Phys. Commun. 109 (1998), 193-206 doi:10.1016/S0010-4655(98)00010-1 [arXiv:nucl-th/9709009 [nucl-th]].
  * (27) B. Andersson, G. Gustafson, G. Ingelman and T. Sjostrand, Phys. Rept. 97 (1983), 31-145 doi:10.1016/0370-1573(83)90080-7
  * (28) Z. w. Lin and C. M. Ko, Phys. Rev. C 65 (2002), 034904 doi:10.1103/PhysRevC.65.034904 [arXiv:nucl-th/0108039 [nucl-th]].
  * (29) D. Molnar and S. A. Voloshin, Phys. Rev. Lett. 91 (2003), 092301 doi:10.1103/PhysRevLett.91.092301 [arXiv:nucl-th/0302014 [nucl-th]].
  * (30) B. A. Li and C. M. Ko, Phys. Rev. C 52 (1995), 2037-2063 doi:10.1103/PhysRevC.52.2037 [arXiv:nucl-th/9505016 [nucl-th]].
  * (31) M. M. Aggarwal et al. [STAR], [arXiv:1007.2613 [nucl-ex]].
  * (32) B. B. Abelev et al. [ALICE], Phys. Lett. B 727 (2013), 371-380 doi:10.1016/j.physletb.2013.10.054 [arXiv:1307.1094 [nucl-ex]].
  * (33) L. Adamczyk et al. [STAR], Phys. Rev. C 88 (2013), 014902 doi:10.1103/PhysRevC.88.014902 [arXiv:1301.2348 [nucl-ex]].
  * (34) C. Zhang, L. Zheng, F. Liu, S. Shi and Z. W. Lin, Phys. Rev. C 99 (2019) no.6, 064906 doi:10.1103/PhysRevC.99.064906 [arXiv:1903.03292 [nucl-th]].
  * (35) K. Tiwari and M. Nasim, Nucl. Phys. A 999 (2020), 121751 doi:10.1016/j.nuclphysa.2020.121751 [arXiv:2002.01201 [nucl-th]].
  * (36) Y. He and Z. W. Lin, Phys. Rev. C 96 (2017) no.1, 014910 doi:10.1103/PhysRevC.96.014910 [arXiv:1703.02673 [nucl-th]].
  * (37) Y. Xia, J. Xu, B. A. Li and W. Q. Shen, Nucl. Phys. A 955 (2016), 41-57 doi:10.1016/j.nuclphysa.2016.06.001 [arXiv:1411.3057 [nucl-th]].
  * (38) B. I. Abelev et al. [STAR], Phys. Lett. B 655 (2007), 104-113 doi:10.1016/j.physletb.2007.06.035 [arXiv:nucl-ex/0703040 [nucl-ex]].
  * (39) L. Adamczyk et al. [STAR], Phys. Rev. C 96 (2017) no.4, 044904 doi:10.1103/PhysRevC.96.044904 [arXiv:1701.07065 [nucl-ex]].
  * (40) B. Abelev et al. [ALICE], Phys. Rev. C 88 (2013), 044910 doi:10.1103/PhysRevC.88.044910 [arXiv:1303.0737 [hep-ex]].
  * (41) D. Zhang [STAR], Nucl. Phys. A 1005 (2021), 121825 doi:10.1016/j.nuclphysa.2020.121825 [arXiv:2002.10677 [nucl-ex]].


††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work

# Tunable spin polarization and electronic structure of bottom-up synthesized
MoSi2N4 materials

Rajibul Islam  International Research Centre MagTop, Institute of Physics,
Polish Academy of Sciences, Aleja Lotnikow 32/46, PL-02668 Warsaw, Poland
Barun Ghosh  Department of Physics, Indian Institute of Technology, Kanpur
208016, India  Department of Physics, Northeastern University, Boston,
Massachusetts 02115, USA  Carmine Autieri  International Research Centre
MagTop, Institute of Physics, Polish Academy of Sciences, Aleja Lotnikow
32/46, PL-02668 Warsaw, Poland  Sugata Chowdhury  Department of Physics and
Astronomy, Howard University, Washington, D.C. 20059, USA  Arun Bansil
Department of Physics, Northeastern University, Boston, Massachusetts 02115,
USA  Amit Agarwal  Department of Physics, Indian Institute of Technology,
Kanpur 208016, India  Bahadur Singh<EMAIL_ADDRESS>Department of
Condensed Matter Physics and Materials Science, Tata Institute of Fundamental
Research, Mumbai 400005, India

###### Abstract

Manipulation of spin-polarized electronic states of two-dimensional (2D)
materials under ambient conditions is necessary for developing new quantum
devices with small physical dimensions. Here, we explore spin-dependent
electronic structures of ultra-thin films of recently introduced 2D synthetic
materials MSi2Z4 (M = Mo or W and Z = N or As) using first-principles
modeling. Stacking of MSi2Z4 monolayers is found to generate dynamically
stable bilayer and bulk materials with thickness-dependent properties. When
spin-orbit coupling (SOC) is included in the computations, MSi2P4 monolayers
display indirect bandgaps and large spin-split states at the $K$ and
$K^{\prime}$ symmetry points at the corners of the Brillouin zone with nearly
100% spin-polarization. The spins are locked in opposite directions along an
out-of-the-plane direction at $K$ and $K^{\prime}$, leading to spin-valley
coupling effects. As expected, spin-polarization due to the presence of
inversion symmetry in the pristine bilayers is absent, but it can be induced
via an external out-of-plane electric field much like the case of Mo(W)S2
bilayers. A transition from indirect to direct bandgap can be driven by
replacing N by As in MSi2(N, As)4 monolayers. Our study indicates that the
MSi2Z4 materials can provide a viable alternative to the MoS2 class of 2D
materials for valleytronic and optoelectronic applications.

Introduction. Since the isolation of two-dimensional (2D) graphene from its
parent graphite in 2004 Novoselov  _et al._ (2004, 2005); Zhang  _et al._
(2005), a variety of atomically-thin materials have been exfoliated from bulk
layered compounds with electronic states covering insulators to semiconductors
to semimetals/metals. Prominent examples include hexagonal boron nitride Lin
_et al._ (2010), 2D transition metal dichalcogenides (TMDs) Mak  _et al._ ;
Wang  _et al._ (2012); Chhowalla  _et al._ (2013); Huang  _et al._ (2013); Tan
and Zhang (2015); Yuan  _et al._ (2013); Chang  _et al._ (2014), phosphorene
Li  _et al._ (2014); Liu  _et al._ (2015), and MXenes Naguib  _et al._ (2014),
among other materials Novoselov  _et al._ (2016). These 2D materials offer
exciting opportunities for exploring novel electronic, excitonic, correlated,
and topological states under 2D charge confinement for spintronic,
valleytronic, and optoelectronic applications and developing materials
platforms for high-density devices with minimal physical dimensions. Stacking,
twisting, and straining of such 2D layers to form moire superlattices and
heterostructures brings unprecedented possibilities for tailoring properties
Novoselov and Neto (2012); Heine (2015); Gong  _et al._ (2014); Wang  _et al._
(2017); Cao  _et al._ (2018); Novoselov  _et al._ (2016); Chen  _et al._
(2018); Vargas  _et al._ (2017). A common approach for obtaining 2D materials
is exfoliation from appropriate 3D layered materials using a top-to-bottom
approach. Finding new 2D materials without parental analogs can provide a new
paradigm for engineering states with diverse functionalities and offer new
pathways for designing synthetic materials with desirable properties Novoselov
and Neto (2012); Heine (2015); Gong  _et al._ (2014); Wang  _et al._ (2017);
Cao  _et al._ (2018); Novoselov  _et al._ (2016); Chen  _et al._ (2018);
Vargas  _et al._ (2017); Novoselov  _et al._ (2016).

Among the methods of growing materials in a bottom-up approach is the use of a
substrate with strong adatom adhesion. This method has shown success in
synthesizing atomically-thin films such as silicene Fleurence  _et al._
(2012), germanene Dávila  _et al._ (2014), bismuthene Cheng  _et al._ (2014),
and borophene Feng  _et al._ (2017). The stability and morphology of such
materials are, however, strongly dependent on growth conditions due to the
presence of dangling bonds of adatoms that either reorganize to generate
complicated surface morphologies or get oxidized when exposed to air
Grazianetti  _et al._ (2015). An alternate route proposed recently involves
passivation of the high-energy surfaces of materials with elements that can
generate synthetic layered 2D materials Hong  _et al._ (2020); Novoselov
(2020). By passivating non-layered molybdenum nitride with elemental silicon
during chemical vapor deposition growth, large area (15 mm$\times$ 15 mm)
layered 2D MoSi2N4 materials were synthesized. Importantly, MoSi2N4 shows
remarkable properties such as stability under ambient conditions, a
semiconducting behavior, and high mobility of 270/1200 cm2V-1s-1, which is
better than that of the widely used MoS2 class of 2D materials Hong  _et al._
(2020); Novoselov (2020); Bertolazzi  _et al._ (2011); Cai  _et al._ (2014).
MoSi2N4 and its derivative monolayers host gapped states in a pair of valleys
located at the corners of the hexagonal Brillouin zone (BZ) Li  _et al._
(2020); Ai  _et al._ (2021); Yang  _et al._ (2021). Due to the breaking of the
spatial inversion symmetry, the spin states in these monolayers become
separated in energy and give rise to unique spin-valley couplings in the
vicinity of the Fermi level and valley-contrasting Berry curvatures and
orbital magnetic moments, which could potentially enable wide-ranging
valleytronics and optoelectronics applications Gunawan  _et al._ (2006); Xiao
_et al._ (2007); Gunlycke and White (2011); Xiao  _et al._ (2012); Cai  _et
al._ (2013); Xu  _et al._ (2014); Autieri  _et al._ (2017)

Motivated by the new possibilities offered by a bottom-up approach and the
excellent stability of MoSi2N4 2D materials under ambient conditions, here we
report layer-dependent stability and valleytronic properties of MSi2Z4 (M = Mo
or W, and Z = N or As) materials. Using density-functional-perturbation theory
based first-principles modeling, we show that the MoSi2N4 materials are
dynamically stable up to the bulk limit. The monolayers are found to exhibit
large spin-split states at the BZ corners $K$ and $K^{\prime}$ with nearly
100% spin-polarization, similar to MoS2 materials class. As expected, the
spin-splitting is zero in the bilayer films due to the presence of inversion
symmetry. However, spin-splitting can be switched on and manipulated in the
bilayers via an out-of-plane electric field. An indirect to direct bandgap
transition in MSi2Z4 is driven by the replacement of N by As. In addition to
highlighting unique thickness-dependent properties of MSi2Z4, our study
highlights the value of a bottom-up approach for synthesizing viable 3D bulk
materials based on synthetic 2D vdW materials.

Methods. Electronic structure calculations were performed within the density
functional theory (DFT) framework using the Vienna ab-initio simulation
package (VASP) Kresse and Furthmüller (1996); Kresse and Joubert (1999). The
projector augmented wave (PAW) pseudopotentials were used with generalized
gradient approximation Perdew  _et al._ (1996) to consider exchange-
correlation effects. A plane-wave cutoff of 500 eV was used throughout the
calculations. The surface BZ integration was performed using a $10\times
10\times 1$ Monkhorst-pack $k-$grid. The spin-orbit coupling (SOC) was
considered self-consistently. The structural parameters were optimized until
the residual forces on each atom became less than 10-4 eV/Å, and these
optimized parameters were used in all our calculations. An energy tolerance of
10-8 eV was used. The thin-film calculations were performed using slab
geometry with a vacuum of 20 Å to eliminate the spurious interaction between
the periodically repeated 2D layers. The phonon dispersion curves were
calculated within the density functional perturbation theory (DFPT) framework
using PHONOPY code Togo and Tanaka (2015) with a $4\times 4\times 1$
supercell. The PyProcar Herath  _et al._ (2020) and Pymatgen Ong  _et al._
(2013) packages were used for band structure illustrations.

Crystal structure and dynamical stability of MoSi2N4. The monolayer MoSi2N4
crystallizes in the hexagonal lattice with space group $D^{1}_{3h}$
($P\bar{6}m2$, No. 187). It consists of strongly bonded seven sublayers
stacking in the order of N-Si-N-Mo-N-Si-N that can be viewed as a sandwich of
MoN2 layer by two Si-N bilayers (see Fig. 1(a)-(e)). This monolayer structure
preserves a trigonal $C_{3v}$ and $M_{z}(z\rightarrow-z)$ mirror plane
symmetries and breaks inversion symmetry. The monolayers are stacked by weak
van der Waals forces in -A-B-A- order to realize a 2H bilayer structure in a
way similar to MoS2. Unlike monolayers, the bilayer MoSi2N4 realizes a high-
symmetry group $D^{4}_{6h}$ ($P6_{3}/mmc$, No. 194) Hong  _et al._ (2020);
Zhong  _et al._ (2021), restoring the spatial center of inversion which is
marked by the red dot in Fig. 1(b). The equilibrium interlayer distance
($d_{0}$) between Mo1 and Mo2 sublayers in the bilayer is found to be 10.65 Å.
Notably, the 2H-bilayer structure can be repeated to realize the bulk of
MoSi2N4 materials in a way similar to transition metal dichalcogenides. The
optimized structural parameters and Wyckoff positions for bulk MSi2Z4 are
listed in Table 1.

Figure 1: Atomic arrangements of (a) four and (b) two layers of MoSi2N4 with
AB stacking. The dashed box identifies the bulk unit cell of 2H-MoSi2N4. The
red dot at the middle of the van der Waals gap in (b) marks the spatial center
of inversion, which is absent in the monolayer structure. (c) Top view of
monolayer MoSi2N4. (d) Mo-N trigonal and (e) Si-N tetrahedral local
coordination structures of MoSi2N4 monolayers. The calculated phonon
dispersion of (f) monolayer (1ML), (g) bilayer (2ML), and (h) bulk MoSi2N4.

To showcase the stability of the monolayer and multilayer MoSi2N4, we present
the calculated phonon dispersion curves for monolayer, bilayer, and bulk
structures in Figs. 1(f), (g), and (h), respectively. The absence of the
imaginary phonon frequencies in the whole hexagonal BZ confirms the dynamical
stability of these structures. Importantly, the bulk phonon spectrum lacks
imaginary phonon frequencies similar to monolayer and bilayer MoSi2N4. This
infers that the stable 3D bulk of MoSi2N4 can be realized experimentally.

Table 1: Calculated lattice parameters for 2H-bulk MoSi2N4, MoSi2As4, WSi2P4, and WSi2As4. $a$ and $c$ are the hexagonal lattice constants and $u_{Si}$ , $u_{N/As}$, and $v_{N/As}$ are the internal parameter associated with Wyckoff position $4e~{}(0,0,u_{Si})$, $4f~{}(\frac{1}{3},\frac{2}{3},u_{N/As})$, and $4e~{}(0,0,v_{N/As})$, respectively. The subscripts identify the atoms. | a (Å) | c (Å) | $u_{Si}$ | $u_{N/As}$ | $v_{N/As}$ | bandgap (eV)
---|---|---|---|---|---|---
MoSi2N4 | 2.910 | 21.311 | 0.1095 | 0.1915 | 0.0859 | 1.655
MoSi2As4 | 3.622 | 27.617 | 0.1106 | 0.1960 | 0.0703 | 0.508
WSi2N4 | 2.914 | 21.439 | 0.1099 | 0.1914 | 0.0865 | 1.970
WSi2As4 | 3.628 | 27.940 | 0.1121 | 0.1967 | 0.0723 | 0.207

Electronic and spin structure of monolayer MoSi2N4. We present the orbital-
resolved band structure of monolayer MoSi2N4 without SOC in Fig 2(a). An
indirect bandgap of $\sim 1.778$ eV is obtained between the valence band
maxima (VBM) and conduction band minima (CBM) which are located at zone-center
$\Gamma$ and zone-boundary $K/K^{\prime}$ points, respectively. Notably, the
energy difference of top of valence bands, $\Delta_{\Gamma K}$, at $\Gamma$
and $K/K^{\prime}$ points is $~{}322$ meV. It can be tuned by strain to
realize a direct band gap at $K/K^{\prime}$ point Li  _et al._ (2020). The VBM
and CBM consist of $d_{z^{2}}$ states of the Mo atoms. All the states remain
two-fold spin degenerate without SOC as seen in Fig. 2(a). When SOC is
included in the computations, the top valence bands display a large spin-
splitting of $\sim$129 meV at $K$ due to the broken spatial inversion
symmetry. Since $K$ is not a time-reversal invariant momentum (TRIM) point,
the spin-split states are not two-fold degenerate at $K$. In contrast, the
bands at $\Gamma$ and $M$ point remain two-fold spin degenerate as they are
TRIM points [see Figs. 2(b) and (c)]. Despite the spin-resolved states, the
indirect nature of the monolayer bandgap remains preserved with a value of
$\sim 1.774$ eV with SOC.

An analysis of spin-polarization reveals that the two spin-split states at $K$
have nearly 100% out-of-plane ($S_{z}$) spin-polarization. This can be
attributed to the presence of the horizontal mirror plane $M_{z}$ in monolayer
MoSi2N4 that ensures the dissolution of the $S_{x}$ and $S_{y}$ spin
components. The spin-split states at the $K$ and $K^{\prime}$ have opposite
spin-polarization consistent with the fact that they are Kramers pair obeying
time-reversal symmetry constraint $E(\vec{k},\uparrow)$ =
$E(-\vec{k},\downarrow)$. Figure 2(d) shows the evolution of spin-polarization
of the top valence band with momentum away from the $K$ point. The spin-
polarization decays slightly to 99.88% for $\Delta k=$0.553 Å-1 ( 38% of the
$\Gamma-K$ distance) change in the momentum, demonstrating its robustness. The
spin texture of the top valence band in the hexagonal BZ is shown
schematically in Fig. 2(f). Such spin behavior reveals a spin-valley locking
in MoSi2N4 monolayer similar to the previously observed spin-valley locking
for MoS2 series of TMDs Chang  _et al._ (2014).

Figure 2: Orbital-resolved band structure of monolayer MoSi2N4 (a) without and
(b) with spin-orbit coupling (SOC). The spin-resolved bands around $K$ in
$\Gamma-K$ directions for (c) $E_{z}$ = 0 eV/Å and (d) $E_{z}$ = 0.03 eV/Å
with SOC. Colorbar in (d) denotes the strength of the spin-polarization. (e)
Spin-polarization decay profile of topmost valence band near the $K$ point.
Large spin polarization ($>99.9\%$) persists over a wide momentum range along
the $\Gamma-K$ direction. (f) Schematic representation of spin-valley locking
in monolayer MoSi2N4. Red and blue represent out-of-the-plane and inward spin
polarizations, respectively.

We emphasize that MoSi2N4 monolayer shows a Zeeman-type out-of-the-plane spin
polarization in the vicinity of $K$ points due to its structural constraints.
The spin-polarization is an intrinsic property tied to the crystal structure
and thus cannot be destroyed or manipulated with an out-of-plane electric
field $E_{z}$. We have verified this by calculating the band structure and
spin texture by applying a vertical electric field to the monolayer plane.
Figure 2(d) shows the results for an electric field $E_{z}=0.03$ eV/Å. Both
the spin-splitting and spin-polarization features are retained with the
electric field, indicating the robustness of these states.

Electric-field tunable spin structure of bilayer MoSi2N4. Figure 3(a) shows
the calculated band structure of bilayer MoSi2N4. Similar to the monolayer, it
is an indirect bandgap semiconductor with VBM and CBM located at the $\Gamma$
and $K/K^{\prime}$ points, respectively. However, in contrast to the
monolayer, the inversion symmetry is restored in the bilayer and as a
consequence all the bands become two-fold spin-degenerate in the full BZ. A
small splitting at the $\Gamma$ point appears due to interlayer interaction
between two MoSi2N4 layers whereas the splitting at the $K/K^{\prime}$ points
arises due to the SOC. Importantly, the inversion symmetry of bilayer MoSi2N4
can be broken by an out-of-plane external electric field $E_{z}$. It lifts the
spin-degeneracy at non-TRIM $K/K^{\prime}$ points, allowing the manipulation
of spin-splitting at the top of the valence bands with $E_{z}$.

Figure 3: (a) Band structure of bilayer MoSi2N4 in the absence of external
electric field (Ez = 0 eV/Å). (b) Same as (a) but for electric field
$E_{z}=0.03$ eV/Å. The spin-splitting in the band structure is evident due to
the breaking of inversion symmetry by an external electric field. (c)
Evolution of the top four valence bands around the $K$ point for varying
external electric field strengths. Color scale gives the strength of spin-
polarization in each band whereas the numbers 1 and 2 identify spin-split
doublet of first and second layers of bilayer MoSi2N4. (d) The degree of spin-
splitting at $K$ point as a function of electric field $E_{z}$. Blue (red)
markers show the intra (inter)-layer spin-splitting $\Delta_{intra}$
($\Delta_{inter}$). (e) Schematic description of the electric field effect on
the bilayer band structure. The spin-split bands of top and bottom layers in
the absence of external electric field are spin-degenerate due to the presence
of inversion symmetry (middle row). These states lead to spin-split doublet in
the presence of an external electric field (top and bottom rows). The
direction of the applied electric field determines the spin-polarization of
the states near the Fermi level.

Figure 3(b) shows the spin-resolved bilayer band structure for $E_{z}=0.03$
eV/Å. The spin-split states are resolved at $K$ and $K^{\prime}$ points with
opposite spin-polarization for top bands. There are four spin-polarized
valence bands near the Fermi level, two of which come from the first layer
whereas the other two come from the second layer in the bilayer system. The
evolution of these four states with electric field $E_{z}$ is shown in Fig.
3(c). To quantify the change in the spin-split states, we define two energy
parameters $\Delta_{intra}$ and $\Delta_{inter}$. While $\Delta_{intra}$ is
the energy difference between first (second) layer spin-up and first (second)
layer spin-down states, $\Delta_{inter}$ is defined as the energy difference
between first-layer spin-up and second-layer spin-down state. $\Delta_{intra}$
thus gives the spin-splitting due to SOC whereas $\Delta_{inter}$ shows the
splitting due to potential difference between two layers caused by external
field $E_{z}$. When $E_{z}=0.01$ eV/Å, the spin-split doublet states from
second layer are slightly lower than the first layer, making $\Delta_{inter}$
lower than $\Delta_{intra}$. The two top-most valence states are thus
constituted by states from two separate layers. When $E_{z}$ is larger than a
critical value, $\Delta_{inter}$ becomes larger than $\Delta_{intra}$ and
thus, two topmost valence states come from the same layer. The variation of
$\Delta_{intra}$ and $\Delta_{inter}$ with electric field is shown in Fig.
3(d). We find that $\Delta_{inter}$ varies linearly with the electric field
while $\Delta_{intra}$ shows a negligible electric field dependence. A
crossover between $\Delta_{intra}$ and $\Delta_{inter}$ is observed for an
electric field of $E_{z}=0.012$ eV/Å. Notably, the spin polarization of the
top most valence states at the $K/K^{\prime}$ points remains nearly 100% in
the presence $E_{z}$.

Figure 4: (a) Orbital-resolved band structure of bulk MoSi2N4 in the bulk
hexagonal Brillouin zone. (b) Band gap and (c) average spin-polarization of
the top layer as a function of the number of layers of MoSi2N4.

We find that the applied electric field $E_{z}$ in bilayer changes the
splitting ($\Delta_{inter}$) between states coming from different layers. It
shows the negligible effect on the spin-splitting in states coming from the
same layer as well as on the degree of spin-polarization of the states. The
sign of spin-polarization of states at $K/K^{\prime}$ points is however
electric-field direction-dependent. The evolution of bilayer states at the $K$
points under positive and negative field direction is shown schematically in
Fig 3(e). These results provide a clear pathway for manipulating the spin-
states in bilayer MoSi2N4. Such electric-field evolution of bilayer states
remains robust for all investigated MSi2Z4 materials. Interestingly, the
values of the electric field required to manipulate these states are much
lower than MoS2 Chang  _et al._ (2014).

Layer-dependent states and spin polarization. We now turn to discuss the
evolution of bandgap and spin-polarization of multilayer MoSi2N4. Figure 4(a)
shows the calculated 2H bulk band structure using our optimized lattice
parameters (Table 1). It has an indirect bandgap of $~{}1.655$ eV. The CBM at
$K$ and VBM at $\Gamma$ consist of Mo $d_{z^{2}}$ states of Mo atoms similar
to the monolayer and bilayer systems. The bands along the $\Gamma-A$ direction
remain weakly dispersive showing a weak interlayer coupling. However, the SOC
split states can be seen at the $K$ and $H$ points. The evolution of bandgap
as a function of the layer thickness is shown in Fig. 4(b). The bandgap
reduces slightly with the number of MoSi2N4 layers and converges to the bulk
value of $~{}1.655$ eV at eight layers thin-film. Such a weak reduction of
bandgap with the number of layers indicates that weak Van der Waal’s coupling
dominates the interlayer interactions in the MoSi2N4.

Figure 4(c) shows the evolution of the degree of spin-polarization of valence
state as a function of the number of layers. Since the system with an even
number of MoSi2N4 layers are inversion symmetric, the degree of spin-
polarization of their states is zero. Adding a layer to such an even layer
stack breaks inversion symmetry and induces a net spin-polarization in the
system. We find the spin-polarization of odd-layered systems varies as
$~{}1/N$, where $N$ is the number of odd layers.

Figure 5: Calculated phonon spectrum of (a) monolayer (1ML), (b) bilayer
(2ML), (c) and bulk MoSi2As4. The orbital-decomposed band structure of (d)
monolayer, (e) bilayer, and (f) bulk MoSi2As4.

Band structure of MSi2Z4 materials. We now explore the dynamical stability and
band-structure of other MSi2Z4 thin films similar to MoSi2N4. Figure 5(a)-(c)
shows the phonon spectrum of monolayer, bilayer, and bulk MoSi2As4. The phonon
dispersion curves are stable without the presence of any imaginary branch in
the whole BZ, showing the dynamical stability of MoSi2As4 layers. The band
structures of MoSi2As4 films and bulk are displayed in Fig. 5(d)-(f). In
contrast to MoSi2N4 layers, the monolayer MoSi2As4 is a direct bandgap
semiconductor with a value of $~{}0.508$ eV at the $K/K^{\prime}$ point. The
direct bandgap nature stays preserved as one increases the thickness from
monolayer to bulk. The interlayer coupling between MoSi2As4 is large as
compared to MoSi2N4 and thus, a direct bandgap location changes from $K$ to
$H$ point in the bulk limit [Fig. 5(f)]. Furthermore, the spin-split states
with nearly 100% spin-polarization remain preserved in MoSi2As4 monolayers. We
present phonon spectrum and orbital resolved bandstructure of WSi2N4 and
WSi2As4 in Supplementary Materials (SMs). Both of these systems are stable up
to the bulk limit and hold the spin-polarized states similar to MoSi2N4 and
MoSi2As4. However, the presence of W atoms enhances the spin-splitting at the
$K/K^{\prime}$ points in these materials owing to their inherent large SOC.

Conclusion. Using first-principles modeling, we have carried out a systematic
thickness-dependent investigation of dynamical stabilities and electronic and
spin-polarization properties of the MSi2Z4 (M = Mo or W and Z = N or As)
compounds. These materials are found to be dynamically stable from the
monolayer to the bulk limit, indicating that multilayer films and bulk of such
bottom-up synthetic 2D vdW materials should be possible to realize
experimentally. Our analysis reveals that the monolayers host two nearly 100%
out-of-the-plane spin-polarized states at the $K$ points in the BZ with
Zeeman-type spin splittings. The spin-polarization is reversed at the
$K^{\prime}$ points while the degree of spin-polarization remains preserved.
The spin-polarization of the states in the bilayers, which is zero due to the
restoration of the inversion symmetry in the pristine bilayers, can be
switched on and manipulated using an external electric field. MoSi2N4 and
WSi2N4 exhibit a robust indirect bandgap in the visible and infrared from the
monolayer to the bulk limit. In contrast, MoSi2As4 and WSi2As4 monolayers
display a direct bandgap at the $K/K^{\prime}$-point, which is preserved from
the monolayer to the bulk. Our study provides insight into the bandgap, spin-
polarization, and spin-valley locking of electronic states in MSi2Z4
materials, and indicates that these materials could provide a viable materials
platform as an alternative to the MoS2 materials in common current use for
spintronic, valleytronic and optoelectronic applications.

## Acknowledgments

We thank Tomasz Dietl for valuable discussions. The work is supported by the
Foundation for Polish Science through the international research agendas
program co-financed by the European union within the smart growth operational
programme. We acknowledge the access to the computing facilities of the
Interdisciplinary Center of Modeling at the University of Warsaw, Grant No.
G75-10, GB84-0 and GB84-7. R. I and C. A. acknowledge support by Narodowe
Centrum Nauki (NCN, National Science Centre, Poland) Project
No.2020/37/N/ST3/02338. We acknowledge IIT Kanpur, Science Engineering and
Research Board (SERB) and the Department of Science and Technology (DST) for
financial support. The work at Northeastern University was supported by the
Air Force Office of Scientific Research under award number FA955-20-1-0322,
and benefited from the computational resources of Northeastern University’s
Advanced Scientific Computation Center (ASCC) and the Discovery Cluster.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. Grigorieva, S. Dubonos, Firsov, and AA, Nature 438, 197 (2005).
  * Zhang  _et al._ (2005) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005).
  * Lin  _et al._ (2010) Y. Lin, T. V. Williams, and J. W. Connell, J. Phys. Chem. Lett. 1, 277 (2010).
  * (5) K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 136805.
  * Wang  _et al._ (2012) Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, Nat. Nanotechnol. 7, 699 (2012).
  * Chhowalla  _et al._ (2013) M. Chhowalla, H. S. Shin, G. Eda, L.-J. Li, K. P. Loh, and H. Zhang, Nat. Chem. 5, 263 (2013).
  * Huang  _et al._ (2013) X. Huang, Z. Zeng, and H. Zhang, Chem. Soc. Rev. 42, 1934 (2013).
  * Tan and Zhang (2015) C. Tan and H. Zhang, Chem. Soc. Rev. 44, 2713 (2015).
  * Yuan  _et al._ (2013) H. Yuan, M. S. Bahramy, K. Morimoto, S. Wu, K. Nomura, B.-J. Yang, H. Shimotani, R. Suzuki, M. Toh, C. Kloc, X. Xu, R. Arita, N. Nagaosa, and Y. Iwasa, Nat. Phys. 9, 563 (2013).
  * Chang  _et al._ (2014) T.-R. Chang, H. Lin, H.-T. Jeng, and A. Bansil, Sci. Rep. 4, 6270 (2014).
  * Li  _et al._ (2014) L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nat. Nanotechnol. 9, 372 (2014).
  * Liu  _et al._ (2015) H. Liu, Y. Du, Y. Deng, and D. Y. Peide, Chem. Soc. Rev. 44, 2732 (2015).
  * Naguib  _et al._ (2014) M. Naguib, V. N. Mochalin, M. W. Barsoum, and Y. Gogotsi, Adv. Mater. 26, 992 (2014).
  * Novoselov  _et al._ (2016) K. Novoselov, O. A. Mishchenko, O. A. Carvalho, and A. C. Neto, Science 353 (2016).
  * Novoselov and Neto (2012) K. S. Novoselov and A. H. C. Neto, Phys. Scr. T146, 014006 (2012).
  * Heine (2015) T. Heine, Acc. Chem. Res. 48, 65 (2015).
  * Gong  _et al._ (2014) Y. Gong, J. Lin, X. Wang, G. Shi, S. Lei, Z. Lin, X. Zou, G. Ye, R. Vajtai, B. I. Yakobson, _et al._ , Nat. Mater. 13, 1135 (2014).
  * Wang  _et al._ (2017) J. Wang, F. Ma, and M. Sun, RSC Adv. 7, 16801 (2017).
  * Cao  _et al._ (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * Chen  _et al._ (2018) C. Chen, B. Singh, H. Lin, and V. M. Pereira, Phys. Rev. Lett. 121, 226602 (2018).
  * Vargas  _et al._ (2017) A. Vargas, F. Liu, C. Lane, D. Rubin, I. Bilgin, Z. Hennighausen, M. DeCapua, A. Bansil, and S. Kar, Sci. Adv. 3, e1601741 (2017).
  * Fleurence  _et al._ (2012) A. Fleurence, R. Friedlein, T. Ozaki, H. Kawai, Y. Wang, and Y. Yamada-Takamura, Phys. Rev. Lett. 108, 245501 (2012).
  * Dávila  _et al._ (2014) M. E. Dávila, L. Xian, S. Cahangirov, A. Rubio, and G. L. Lay, New J. Phys. 16, 095002 (2014).
  * Cheng  _et al._ (2014) L. Cheng, H. Liu, X. Tan, J. Zhang, J. Wei, H. Lv, J. Shi, and X. Tang, J. Phys. Chem. C 118, 904 (2014).
  * Feng  _et al._ (2017) B. Feng, O. Sugino, R.-Y. Liu, J. Zhang, R. Yukawa, M. Kawamura, T. Iimori, H. Kim, Y. Hasegawa, H. Li, L. Chen, K. Wu, H. Kumigashira, F. Komori, T.-C. Chiang, S. Meng, and I. Matsuda, Phys. Rev. Lett. 118, 096401 (2017).
  * Grazianetti  _et al._ (2015) C. Grazianetti, D. Chiappe, E. Cinquanta, M. Fanciulli, and A. Molle, J. Condens. Matter Phys. 27, 255005 (2015).
  * Hong  _et al._ (2020) Y.-L. Hong, Z. Liu, L. Wang, T. Zhou, W. Ma, C. Xu, S. Feng, L. Chen, M.-L. Chen, D.-M. Sun, _et al._ , Science 369, 670 (2020).
  * Novoselov (2020) K. S. Novoselov, Natl. Sci. Rev. 7, 1842 (2020).
  * Bertolazzi  _et al._ (2011) S. Bertolazzi, J. Brivio, and A. Kis, ACS nano 5, 9703 (2011).
  * Cai  _et al._ (2014) Y. Cai, G. Zhang, and Y.-W. Zhang, J. Am. Chem. Soc. 136, 6269 (2014).
  * Li  _et al._ (2020) S. Li, W. Wu, X. Feng, S. Guan, W. Feng, Y. Yao, and S. A. Yang, Phys. Rev. B 102, 235435 (2020).
  * Ai  _et al._ (2021) H. Ai, D. Liu, J. Geng, S. Wang, K. H. Lo, and H. Pan, Phys. Chem. Chem. Phys. 23, 3144 (2021).
  * Yang  _et al._ (2021) C. Yang, Z. Song, X. Sun, and J. Lu, Phys. Rev. B 103, 035308 (2021).
  * Gunawan  _et al._ (2006) O. Gunawan, Y. P. Shkolnikov, K. Vakili, T. Gokmen, E. P. De Poortere, and M. Shayegan, Phys. Rev. Lett. 97, 186404 (2006).
  * Xiao  _et al._ (2007) D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007).
  * Gunlycke and White (2011) D. Gunlycke and C. T. White, Phys. Rev. Lett. 106, 136806 (2011).
  * Xiao  _et al._ (2012) D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. 108, 196802 (2012).
  * Cai  _et al._ (2013) T. Cai, S. A. Yang, X. Li, F. Zhang, J. Shi, W. Yao, and Q. Niu, Phys. Rev. B 88, 115140 (2013).
  * Xu  _et al._ (2014) X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Nat. Phys. 10, 343 (2014).
  * Autieri  _et al._ (2017) C. Autieri, A. Bouhon, and B. Sanyal, Philos. Mag. 97, 3381 (2017).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Togo and Tanaka (2015) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Herath  _et al._ (2020) U. Herath, P. Tavadze, X. He, E. Bousquet, S. Singh, F. Munoz, and A. H. Romero, Comput. Phys. Commun. 251, 107080 (2020).
  * Ong  _et al._ (2013) S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson, and G. Ceder, Comput. Mater. Sci. 68, 314 (2013).
  * Zhong  _et al._ (2021) H. Zhong, W. Xiong, P. Lv, J. Yu, and S. Yuan, Phys. Rev. B 103, 085124 (2021).