11institutetext:  Key Laboratory of Optical Astronomy, National Astronomical
Observatories, Chinese Academy of Sciences,
20A Datun Road, Chaoyang District, Beijing 100012,People’s Republic of China
22institutetext: College of Astronomy and Space Sciences, University of
Chinese Academy of Sciences, Beijing 100049, China 33institutetext: Department
of Astronomy, Beijing Normal University, Beijing 100875, People’s Republic of
China 44institutetext: Centro de Estudios de Física del Cosmos de Aragón
(CEFCA), Unidad Asociada al CSIC, Plaza San Juan 1, 44001 Teruel, Spain
55institutetext: Instituto de Astronomia, Geofísica e Ciências Atmosféricas,
Universidade de São Paulo, 05508-090 São Paulo, Brazil 66institutetext:
Departamento de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC Campus,
Camino Bajo del Castillo s/n, E-28692 Villanueva de la Cañada, Madrid, Spain
77institutetext: Observatório Nacional - MCTI (ON), Rua Gal. José Cristino 77,
São Cristóvão, 20921-400 Rio de Janeiro, Brazil 88institutetext: Donostia
International Physics Centre (DIPC), Paseo Manuel de Lardizabal 4, 20018
Donostia-San Sebastián, Spain 99institutetext: IKERBASQUE, Basque Foundation
for Science, 48013, Bilbao, Spain 1010institutetext: University of Michigan,
Department of Astronomy, 1085 South University Ave., Ann Arbor, MI 48109, USA
1111institutetext: University of Alabama, Department of Physics and Astronomy,
Gallalee Hall, Tuscaloosa, AL 35401, USA 1212institutetext: Instituto de
Astrofísica de Canarias, La Laguna, 38205, Tenerife, Spain 1313institutetext:
Departamento de Astrofísica, Universidad de La Laguna, 38206, Tenerife, Spain

# J-PLUS: Support vector regression to measure stellar parameters

Cunshi Wang J-PLUS: Support vector regression to measure stellar parametersJ-
PLUS: Support vector regression to measure stellar parametersJ-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters Yu Bai J-PLUS: Support vector
regression to measure stellar parametersJ-PLUS: Support vector regression to
measure stellar parameters Haibo Yuan J-PLUS: Support vector regression to
measure stellar parametersJ-PLUS: Support vector regression to measure stellar
parameters Jifeng Liu J-PLUS: Support vector regression to measure stellar
parametersJ-PLUS: Support vector regression to measure stellar parametersJ-
PLUS: Support vector regression to measure stellar parametersJ-PLUS: Support
vector regression to measure stellar parameters J.A. Fernández-Ontiveros
J-PLUS: Support vector regression to measure stellar parametersJ-PLUS: Support
vector regression to measure stellar parameters Paula R. T. Coelho J-PLUS:
Support vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters F. Jiménez-Esteban J-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters Carlos Andrés Galarza J-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters R. E. Angulo J-PLUS: Support vector
regression to measure stellar parametersJ-PLUS: Support vector regression to
measure stellar parametersJ-PLUS: Support vector regression to measure stellar
parametersJ-PLUS: Support vector regression to measure stellar parameters A.
J. Cenarro J-PLUS: Support vector regression to measure stellar parametersJ-
PLUS: Support vector regression to measure stellar parameters D. Cristóbal-
Hornillos J-PLUS: Support vector regression to measure stellar parametersJ-
PLUS: Support vector regression to measure stellar parameters R. A. Dupke
J-PLUS: Support vector regression to measure stellar parametersJ-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parametersJ-PLUS: Support vector regression to
measure stellar parametersJ-PLUS: Support vector regression to measure stellar
parametersJ-PLUS: Support vector regression to measure stellar parameters A.
Ederoclite J-PLUS: Support vector regression to measure stellar parametersJ-
PLUS: Support vector regression to measure stellar parameters C. Hernández-
Monteagudo J-PLUS: Support vector regression to measure stellar parametersJ-
PLUS: Support vector regression to measure stellar parametersJ-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters C. López-Sanjuan J-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters A. Marín-Franch J-PLUS: Support
vector regression to measure stellar parametersJ-PLUS: Support vector
regression to measure stellar parameters M. Moles J-PLUS: Support vector
regression to measure stellar parametersJ-PLUS: Support vector regression to
measure stellar parameters L. Sodré Jr J-PLUS: Support vector regression to
measure stellar parametersJ-PLUS: Support vector regression to measure stellar
parameters H. Vázquez Ramió J-PLUS: Support vector regression to measure
stellar parametersJ-PLUS: Support vector regression to measure stellar
parameters J. Varela J-PLUS: Support vector regression to measure stellar
parametersJ-PLUS: Support vector regression to measure stellar parameters

(Received April 29, 2022; accepted April 29, 2022)

###### Abstract

Context. Stellar parameters are among the most important characteristics in
studies of stars which, in traditional methods, are based on atmosphere
models. However, time, cost, and brightness limits restrain the efficiency of
spectral observations. The Javalambre Photometric Local Universe Survey
(J-PLUS) is an observational campaign that aims to obtain photometry in 12
bands. Owing to its characteristics, J-PLUS data have become a valuable
resource for studies of stars. Machine learning provides powerful tools for
efficiently analyzing large data sets, such as the one from J-PLUS, and
enables us to expand the research domain to stellar parameters.

Aims. The main goal of this study is to construct a support vector regression
(SVR) algorithm to estimate stellar parameters of the stars in the first data
release of the J-PLUS observational campaign.

Methods. The training data for the parameters regressions are featured with
12-waveband photometry from J-PLUS and are cross-identified with spectrum-
based catalogs. These catalogs are from the Large Sky Area Multi-Object Fiber
Spectroscopic Telescope, the Apache Point Observatory Galactic Evolution
Experiment, and the Sloan Extension for Galactic Understanding and
Exploration. We then label them with the stellar effective temperature, the
surface gravity, and the metallicity. Ten percent of the sample is held out to
apply a blind test. We develop a new method, a multi-model approach, in order
to fully take into account the uncertainties of both the magnitudes and the
stellar parameters. The method utilizes more than two hundred models to apply
the uncertainty analysis.

Results. We present a catalog of 2,493,424 stars with the root mean square
error of 160K in the effective temperature regression, 0.35 in the surface
gravity regression, and 0.25 in the metallicity regression. We also discuss
the advantages of this multi-model approach and compare it to other machine-
learning methods.

###### Key Words.:

methods: data analysis – techniques: spectroscopic - astronomical databases:
miscellaneous

##  1 Introduction

In modern astronomy, the large amount of raw data produced by the newest
surveys is far beyond our traditional processing capacity, and this
insufficient computing power has become a bottleneck to the rapid development
of astrophysics. Fortunately, the development of computer science has provided
us with new ways to understand and gain knowledge from these raw data. Machine
learning in particular has achieved remarkable success in offering novel
solutions to complex problems based on applications of loss functions,
optimization methods (Ruder 2017), and statistical models (MacKay 2003). Owing
to its nonfunctional algorithms (Cortes & Vapnik 1995; Boser et al. 1992;
Cristianini & Shawe-Taylor 2000; Cover & Hart 1967; Stone 1977; Quinlan 1986;
Breiman 2001), machine learning can reveal potential patterns and important
parameters that are indistinguishable using traditional scientific or
statistical methods.

Stars are the cornerstone of astronomy, and stellar parameters are among the
most crucial characteristics of understanding and characterizing stars. Among
the most powerful and accurate methods of determining stellar parameters is
spectral analysis (Wu et al. 2011; Boeche et al. 2018; Anguiano, B. et al.
2018).

However, spectral observations can be costly. The most powerful spectroscopy
survey — the Large Sky Area Multi-Object Fiber Spectroscopy Telescope (LAMOST)
— has taken spectra of about ten million stars in low resolution and about
three million in medium resolution with a four-meter telescope in the past 8
years. In contrast to spectral observations, photometric observations have
much higher observational efficiency. Compared to LAMOST, the Javalambre
Photometric Local Universe Survey (J-PLUS, Cenarro et al. 2019) observed 223.6
million objects in its five-year observational campaign
111http://J-PLUS.es/datareleases/data_release_dr2. Therefore, many studies
have been working on estimating stellar parameters from photometry-based data
by constructing models of photometric observation results (Bailer-Jones 2011;
Sichevskij 2012; Sichevskij et al. 2014). Their modeling requires highly
precise photometric data from a few surveys, while another mighty tool —
machine learning — can reveal the distribution or pattern hidden inside the
photometry with rougher data.

Machine learning, a cross-disciplinary subject of statistics, optimization,
and computer science, creates algorithms that can process data based on a
well-chosen sample set with features. Different algorithms assign different
models and reveal potential patterns in the sample (MacKay 2003; Shalev-
Shwartz & Ben-David 2014). The algorithms select features from samples and
optimize the parameters based on a given loss function. The loss function
evaluates the difference between the predicted result and the true value.
Machine-learning algorithms have been applied in many different disciplines,
including finance, medical science, and computer vision.

Machine-learning technology has shown the ability to obtain valuable
information from multi-band photometric data. Several studies have made some
efforts to determine stellar parameters from photometric surveys. Bai et al.
(2019) have derived stellar effective temperatures from $Gaia$ second data
release using a random forest (RF) algorithm. Bai et al. (2018) have presented
RF models to categorize objects such as stars, galaxies, and quasi-stellar
objects (QSOs) and classified and obtained the effective temperature of stars.
Lu & Li (2015) developed a scheme to achieve stellar parameters from the Least
Absolute Shrinkage and Selection Operator (LASSO) algorithm and support vector
regression (SVR) models. Yang et al. (2021) designed a cost-sensitive
artificial neural network and achieved stellar parameters from two million
stars from the J-PLUS first data release (DR1). Galarza et al. (2022)
developed the stellar parameters estimation based on the ensemble methods
(SPEEM) pipeline, which is a stack of feature searching, normalization, and
multi-output regressor. The SPEEM is based on RF and extreme gradient boosting
(XGB). These works are all based on machine-learning algorithms which
indicates their powerful ability to reveal potential patterns within data.

J-PLUS has also been used to gain knowledge of objects ranging from our Solar
System to the deep universe, 222http://J-PLUS.es/survey/science such as the
coma cluster (Jiménez-Teja et al. 2019), low metallicity stars (Whitten et al.
2019; Galarza et al. 2022), and galaxy formation (Nogueira-Cavalcante et al.
2019). J-PLUS has 12-waveband photometry, which makes the survey ideal for the
application of machine learning.

In this study, we adopt the SVR algorithm to obtain the stellar parameters,
which include the effective temperature ($T_{\rm eff}$), the surface gravity
(log $g$), and the metallicity ([Fe/H]) of the stars in the J-PLUS DR1
catalog. To construct the training sample, we use data from LAMOST, the Apache
Point Observatory Galactic Evolution Experiment (APOGEE, Zasowski et al.
2013), and the Sloan Extension for Galactic Understanding and Exploration
(SEGUE, Yanny et al. 2009; Sect. 2).

We adopt the SVR (Awad & Khanna 2015; Drucker et al. 1997; Cortes & Vapnik
1995) algorithm and determine the kernel scales (Sect. 3.1). We construct 80
training sets based on the uncertainties of the stellar parameters. Each
parameter of any object in the 80 sets obeys a Gaussian distribution (Sect.
3.2). In all, we construct more than two hundred different models. The blind
test is presented in Sect. 3.3.

The result of applying our method to the entire J-PLUS DR1 catalog is
presented in Sect. 4. In Sect. 5.1 and 5.2, we compare different methods for
constructing the training sample and the inconsistency of stellar parameters
among different pipelines (Sect. 5.3). We discuss the distribution of the
regressed stellar parameters in Sect. 5.4, and compare our result with Yang et
al. (2021) in Sect. 5.5. We present a conclusion in Sect. 6.

##  2 Data

###  2.1 J-PLUS

J-PLUS is conducted by the Observatorio Astrofísico de Javalambre (OAJ,
Teruel, Spain; Cenarro et al. 2014). It uses the 83 cm Javalambre Auxiliary
Survey Telescope (JAST80) and T80Cam, which is a panoramic camera of 9.2k
$\times$ 9.2k pixels that provides a $2\deg^{2}$ field of view (FoV) with a
pixel scale of 0.55 arcsec pix-1 (Marín-Franch et al. 2015). The J-PLUS
contains a 12-passband filter system that comprises five broad ($u,g,r,i,z$)
and seven medium bands from 3000 to 9000 Å. Cenarro et al. (2019) illustrate
the scientific goals and the observational and image reduction strategy of
J-PLUS.

J-PLUS DR1 covers an area of 1,022 ${\rm deg}^{2}$ on the sky with a magnitude
limit of $\sim 21.5$ for a S/N of $\sim$ 3\. These 12 bands provide a large
sample for us to characterize the spectral energy distribution of the detected
sources (Cenarro et al. 2019). The 12-band magnitudes are adopted as our
training features, which are $u$, $J0378$, $J0395$, $J0410$, $J0430$, $g$,
$J0515$, $r$, $J0660$, $i$, $J0861$, and $z$. We name them mag1 to mag12 for
simplicity.

Yuan (2021) recalibrated the J-PLUS DR1 catalog using stellar color
regression. The method is described in detail in Yuan et al. (2015). The
catalog in Yuan (2021) contains 13,265,168 objects, including 4,126,928
objects with all 12 magnitudes. In Wang et al. 2021, the objects in the
recalibrated J-PLUS catalog are classified using the support vector machine
(SVM) algorithm, which distinguishes them among the classes STAR, GALAXY, and
QSO. Here they chose the 12 J-PLUS magnitudes as features and used their
corresponding uncertainties as weight for the SVM. They provide two catalogs
based on the 12 density contours on the 12 magnitudes of the training sample.
The objects that fall into all the 12 contours are assigned as interpolations,
and the others are assigned as extrapolation. Interpolations have better
classification accuracy than extrapolations. In this study, we use all objects
classified as STAR in these two catalogs for the stellar parameter
regressions.

###  2.2 LAMOST spectra

LAMOST is a northern spectroscopic survey situated at Xinglong Observatory,
China. LAMOST is able to observe 4,000 objects simultaneously with a 20 ${\rm
deg}^{2}$ FoV. The main scientific project of LAMOST aims to understand the
structure of the Milky Way (Deng et al. 2012) and external galaxies. We used
the A, F, G, and K catalogs in LAMOST Data Release 7 (DR7), low resolution
spectra (LRS), and medium resolution spectra (MRS)
333http://dr7.lamost.org/catalogue.

LAMOST LRS have a limiting magnitude of about $20$ in the $g$ band and its S/N
is higher than 6 on dark nights or 15 on bright nights. In MRS, the S/N is
always larger than 10. The parameters are given by the LAMOST Stellar
Parameter pipeline (LASP), and their uncertainties are mainly from the stellar
S/N and chi-square of the best-matched theoretical spectrum (Wu et al. 2014;
Liu et al. 2020). The parameter differences of the LRS data and MRS data are
very small, on average 17.6 K for $T_{\rm eff}$, 0.028 dex for log $g$ ($g$ in
${\rm(cms^{-2})}$), and 0.084 dex for [Fe/H], according to our test. The
internal uncertainties of LASP are estimated to be $\Delta T_{\rm eff}\sim
80\rm K$, $\Delta{\rm log}g\sim 0.15\rm dex$, and $\Delta{\rm[Fe/H]}\sim
0.09\rm dex$ (Wang et al. 2020). We cross-matched these catalogs with J-PLUS
DR1 to within one arcsec using the Tool for OPerations on Catalogues And
Tables (TOPCAT, Taylor 2005). There are 216,114 and 25,170 cross-matched stars
in the LRS and MRS catalogs, respectively.

###  2.3 APOGEE

APOGEE (Zasowski et al. 2013) has observed about 150,000 stars in the Milky
Way and obtained precise stellar information, including stellar atmospheric
parameters and radial velocities. The APOGEE Stellar Parameter and Chemical
Abundance Pipeline (ASPCAP, García Pérez et al. 2016) uses a chi-square
minimization to determine the stellar parameters. ASPCAP has high precision in
stellar parameters ($\Delta T_{\rm eff}\sim 2\%$, $\Delta{\rm log}g\sim 0.1\rm
dex$, and $\Delta{\rm[Fe/H]}\sim 0.05\rm dex$). Jönsson et al. (2018) compared
the spectra and determined the uncertainties of the stellar parameters
($\Delta T_{\rm eff}\sim 100\rm K$, $\Delta{\rm log}g\sim 0.2dex$, and
$\Delta{\rm[Fe/H]}\sim 0.1dex$).

We adopted the APOGEE catalog to enlarge our training sample. Using the Sloan
Digital Sky Survey (SDSS) Catalog Archive Server
Jobs444http://skyserver.sdss.org/casjobs/, we extracted table aspcapStar and
find 12,931 stars that satisfy our one arcsec cross-match tolerance.

###  2.4 SEGUE

SEGUE (Yanny et al. 2009) is designed to obtain images in the $u,g,r,i$, and
$z$ wavebands for 3,500 square degrees of the sky located primarily at low
galactic latitudes ($|b|<35^{\circ}$). It delivers observations of about
240,000 stars with a $g$ band magnitude between 14.0 and 20.3 mag and
moderate-resolution spectra from 3,900 to 9,000 Å. The stellar parameters are
presented in SDSS Data Release 7 for these Milky Way stars with S/N greater
than ten. The SEGUE Stellar Parameter Pipeline (SSPP; Lee et al. 2008a, b;
Allende Prieto et al. 2008) developed a multi-method technology to calculate
the stellar parameters which includes non-linear regression models, a minimal
distance of observed spectra and grids of synthetic spectra, and the
correlations between spectral lines and color relations. The average
uncertainties of SSPP are $\Delta T_{\rm eff}\sim 130\rm K$, $\Delta{\rm
log}g\sim 0.21\rm dex$, and $\Delta{\rm[Fe/H]}\sim 0.11\rm dex$. One arcsec
tolerance cross-match yields 25,487 stars.

###  2.5 RAVE

Bai et al. (2018) adopted the RAdial Velocity Experiment (RAVE, Steinmetz et
al. 2020) catalog for their training sample. Our one arcsec tolerance cross-
matching returns only 70 stars. Moreover, the RAVE catalog includes [m/Fe] not
[Fe/H]. The conversion from [m/Fe] to [Fe/H] could import bias or deviation
when using an empirical formula. Therefore, the RAVE catalog is not included
in our training sample.

###  2.6 Normalization

The training sample was constructed with J-PLUS DR1 and three spectral
surveys, 12,931 stars from APOGEE, 216,114 from LAMOST LRS, 25,170 from LAMOST
MRS, and 25,487 from SSPP. After removing the stars with missing photometric
observations in one or more bands, there are 279,702 stars left.

Prior to training our model, we centered all parameters $x_{i}$ in our
training data (both the targets and input magnitudes) by subtracting their
mean: $x_{i,t}=x_{i}-{\rm E}(x_{i})$ 555E stands for mathematical
expectation.. This normalization reduces the upper and lower bound of the
parameter distribution and makes them all zero-mean distributed, which can
accelerate the training speed (Shalev-Shwartz & Ben-David 2014).

During the prediction procedure, the query input magnitudes are centered in
the same fashion, using the respective means from the training data, before
being fed to the model. The model output $x_{i,p}^{*}$ is then returned to the
true target parameter space following: $x_{i,p}=x_{i,p}^{*}+{\rm E}(x_{i})$ to
retrieve the prediction $x_{i,p}$.

We kept the parameters of the same stars from different spectroscopy surveys.
For example, star A appears in both LAMOST and APOGEE, and the corresponding
effective temperatures are $\rm Teff_{1}$ and $\rm Teff_{2}$. In our sample
set, there would be a star $\rm A_{1}$ with effective temperature $\rm
Teff_{1}$ and another star $\rm A_{2}$ with $\rm Teff_{2}$ . The stellar
parameters are only merged when all parameters are the same.

###  2.7 Contours

The prior distributions of the parameters in the training sample set are shown
in Appendix A. The effective temperature is mostly distributed from $5000\sim
6000K$, with log $g$ about $4\sim 4.8$ and a metallicity similar to that of
our Sun, which means that most stars are in their main sequence stage. There
are also some giants with lower surface gravity ($2\sim 3$), lower metallicity
(¡-1), and lower effective temperature ($4000\sim 5000K$).

However, the LAMOST stellar parameters are extracted from the spectra whose
S/N is higher than 10. This criterion makes our LAMOST sample brighter than
the limit magnitude of J-PLUS stars. Such disagreement would decrease the
accuracy of prediction. To solve this problem, we applied the method to
control the interpolation and extrapolation in Wang et al. (2021). This method
could increase the precision in prediction. To quantify the density of the
sample space, we considered 12 sub-space combinations, which are (mag1, mag2,
mag3), … , (mag10, mag11, mag12), (mag11, mag12, mag1), and (mag12, mag1,
mag2), and calculated the 95% density contours to make clear how they cycle
through the parameter space. The sample-dense space can be approached by the
intersection space inside all these twelve density contours. Stars situated in
all contours are assigned as interpolations.

##  3 Methodology

###  3.1 Support vector regression

The SVR (Awad & Khanna 2015; Drucker et al. 1997; Cortes & Vapnik 1995)
algorithm is a regression method based on the SVM (SVM, Cortes & Vapnik 1995;
Boser et al. 1992; Cristianini & Shawe-Taylor 2000; Shalev-Shwartz & Ben-David
2014). The data located inside the margin given by SVR algorithm are not
involved in the calculation. These data gave us the flexibility to define the
tolerance in regression.

For a nonlinear regression problem, we first embedded the sample space into a
new feature space with a higher dimension and transformed it into a linear
regression problem. To transform the problem to the higher dimensional feature
space, the SVR algorithm makes use of the kernel trick, which represents the
inner product of the image of objects mapped from the sample space to the
feature space with a kernel function. This process can accelerate the
calculation by doing calculations in a lower dimension. In the feature space,
the algorithm then fits the data linearly.

Similar to the SVM algorithm, the SVR also has a strip area (called a tube).
The difference between the two algorithms is that SVM maximizes the minimal
distance between each sample to the strip, while the SVR considers how to
minimize the width of the strip to put all samples in. The width of the strip
area is called the margin. The vectors from the final line to the samples that
finally determined the strip area are called the support vector. Details on
the algorithm are given in Smola & Schölkopf (2004).

The root mean square error (RMSE) is one of the most useful standards for
evaluating the efficiency of regression. The RMSE of $\theta$ is given by
$\sqrt{{\bf E}((\hat{\theta}-\theta)^{2})}$, where $\hat{\theta}$ is the
estimate of $\theta$. In regression, the RMSE is given by

${\rm RMSE}=\sqrt{\frac{1}{n}\sum_{t=1}^{n}(\hat{y_{t}}-y_{t})^{2},}$

where y is the stellar parameter.

In the pre-training, we tested both magnitudes and colors as the input
features with the same preprocessing. As a result, the RMSE of the magnitudes
was lower than that of the mag(n-1)$-$mag(n) color under the default setting
of SVR training (with Gaussian kernel 0.83). Then, we constructed independent
models with different Gaussian kernel scales for each stellar parameter. We
present the RMSE as a function of different kernel scales in Fig. 1. We used
the kernel scales with the lowest RMSEs for the training. We also tested the
mag(n-1)$-$mag(n) and mag1$-$mag(n) color(Yang et al. 2021), and optimized the
kernel scale with the Bayesian optimizer. Table 1 presents the adopted kernel
scales and the corresponding RMSEs, which are based on ten-fold validation. We
adopted 0.8 as the kernel scale of $T_{\rm eff}$, 0.325 as the kernel scale of
log $g$, and 0.45 as the kernel scale of [Fe/H].

Table 1: Adopted kernel scales and RMSEs.

| Kernel scale | Magnitude | Kernel Scale C1 | Color1 | Kernel Scale C2 | Color2
---|---|---|---|---|---|---
$T_{\rm eff}$ | 0.8 | 144.059 | 1.01 | 165.71 | 5.77 | 195.77
log $g$ | 0.325 | 0.310 | 1.646 | 0.312 | 0.415 | 0.316
$\rm[Fe/H]$ | 0.45 | 0.221 | 1.944 | 0.233 | 0.440 | 0.233
666The ”Magnitude” column lists the RMSE of using magnitudes as features,
while the column of ”Color1” is the mag(n-1)$-$mag(n) color and ”Color2”,
presents the optimized mag1$-$mag(n) color.

The ’Kernel Scale’ columns match the magnitudes, color1, and color2,
respectively.

Figure 1: RMSEs of the stellar parameters vs. the kernel scale. The blue line
is the effective temperature. The surface gravity is the red solid line with
the star markers, and the metallicity is the dotted line with the cross
markers.

###  3.2 Data enhancement

We duplicated the sample into 80 sets to enhance the cardinality of data for
each stellar parameter. In image recognition, the method is usually carried
out by rotating or mirroring the image and producing more samples for the
training. We applied this method here by generating random training sets to
fully use the uncertainties of both the magnitudes and the stellar parameters.

For each star in our training sample, after the centering process, we
generated 80 stars with its 12 magnitudes, 3 stellar parameters, and 15
corresponding uncertainties. These 80 stars have Gaussian distributions of all
the magnitudes and the stellar parameters following
$\widetilde{x_{i,t}}=x_{i,t}+\sigma(x_{i,t})e_{i}$, where $\sigma(x_{i,t})$ is
the uncertainty of $x_{i,t}$, and $\sigma_{i}\sim{\rm N}(0,1)$. For example,
the mean temperature is 6,000K and the $\sigma$ is 100K for the star with
$T_{\rm eff}=6000\pm 100$K. The simulation process did not change the original
distribution or introduce new errors (Appendix A).

Each of the 80 constructed samples has 251,732 stars. We trained 80 models
using the SVR algorithm with the same kernel scale in Sect. 3.1 for each
stellar parameter, which resulted in a total of 240 different models (same
scheme but different sample set) for $T_{\rm eff}$, log $g,$ and [Fe/H]. We
then used the corresponding 80 regression models to retrieve a distribution of
80 predictions for a given query observation. Then, we used Gaussian functions
to fit these predicted distributions to obtain the centers and standard
deviations.

The data enhancement method takes the uncertainties of stellar parameters into
account, which are not considered in other studies. The absence of these
uncertainties could cause unexpected errors in the prediction since the
spectral precision is not involved in the model construction.

###  3.3 Model validation

We applied model validation to illustrate the effectiveness and to avoid
potential overfitting. A widespread method for model validation is blind
tests, which can reveal potential overfitting and quantify a model’s ability
to generalize to new, previously unseen data. We reserved 10%, 27,970 stars as
our blind test sample, producing 80 sets of tests for each stellar parameter
according to our data enhancement procedure. Similar to the construction of
the sample set, we did not use any criterion to preselected data in order to
avoid the selection effect.

We provide the RMSE and normalized RMSE (NRMSE) of our blind test for both
interpolation and extrapolation in Table 2, where the ${\rm
NRMSE}=\frac{\sqrt{\sum_{i=1}^{n}(\hat{x_{i}}-x_{i})^{2}}}{max(x_{i})-min(x_{i})}$.
The prediction has the best result in effective temperature and the least in
surface gravity. The extrapolation contains 445 stars, which suffered from a
large deviation in the blind test. The result shows that the training model is
better restricted in the space in which the samples are densely distributed.

Table 2: Blind test RMSEs

| BT RMSE | BT NRMSE | Inter RMSE | Inter NRMSE | Extra RMSE | Extra NRMSE
---|---|---|---|---|---|---
$T_{\rm eff}$ | 159.6 | 0.0283 | 150.2 | 0.0266 | 454.0 | 0.0805
log $g$ | 0.3453 | 0.0677 | 0.3380 | 0.0663 | 0.6543 | 0.1283
$\rm[Fe/H]$ | 0.2503 | 0.0484 | 0.2428 | 0.0470 | 0.5402 | 0.1045

777Inter stands for the interpolation and Extra stands for the extrapolation.

Figure 2: Two examples of the $T_{\rm eff}$ distributions. The left panel is
the distribution of $T_{\rm eff}$ of 26109-16130 in the blind test. The right
panel is the distribution of 26016-295 in the prediction. The red lines are
Gaussian fits.

Figure 2 shows the $T_{\rm eff}$ distribution of one star (26109-16130) in our
blind tests. After applying the Lilliefors test (Lilliefors 1967), the
$p$-value of the distribution is 0.375, greater than 0.05, implying that it
may follow a Gaussian distribution. However, this may not generally be the
case (Section 5.4). The results of the blind tests for log $g$ and [Fe/H] are
presented in Appendix B.

##  4 Result

We then applied 240 stellar parameter models to the interpolation and
extrapolation catalogs from Wang et al. (2021), separating our predictions
into two catalogs. There are 2,493,424 objects in the interpolation catalog
and 233,924 objects in the extrapolation catalog. In the classification, the
authors restrained the 12-band magnitudes and constructed 12 contours based on
the sample set distribution. All classified objects were assigned to the area
inside or outside the contours and categorized as either interpolation or
extrapolation. The distribution of inside objects is similar to the samples
and therefore they have higher accuracy. In short, the stars in the
interpolation catalog have a higher classification confidence than those in
the extrapolation catalog.

We also applied this method in the prediction. There are 1,898,154 stars
situated in our training set contours, and 595,270 stars are categorized as
extrapolation in the classification interpolation catalog. In the
classification extrapolation catalog, there are 13,274 interpolations and
220,650 extrapolations. We then input the 12-band magnitudes of each star into
the $3\times 80$ models and obtained the stellar parameters of the
interpolation and extrapolation catalogs. An example of the Gaussian fit is
given in the right panel of Fig. 2. Table 3 is a stellar parameter catalog of
the interpolation catalog.

Table 3: Parameter catalog of the interpolation

ID | R.A. (deg) | Dec. (deg) | Teff (K) |  $\Delta$Teff (K) | log $g$ (dex) |  $\Delta$log $g$ (dex) | [Fe/H] (dex) |  $\Delta$[Fe/H] (dex)
---|---|---|---|---|---|---|---|---
26016-5* | 255.4683 | 22.7600 | 6831 | 101 | 3.71 | 0.14 |  $-$1.20 | 0.10
26016-15 | 255.3051 | 22.7608 | 5950 | 15 | 4.15 | 0.07 |  $-$0.36 | 0.04
26016-16* | 255.3675 | 22.7611 | 5659 | 127 | 3.89 | 0.08 |  $-$0.79 | 0.10
26016-50* | 255.1639 | 22.7644 | 6034 | 80 | 3.71 | 0.13 |  $-$1.19 | 0.17
26016-55* | 255.6631 | 22.7617 | 4061 | 120 | 3.93 | 0.06 |  $-$0.67 | 0.09
26016-64* | 254.7562 | 22.7650 | 6175 | 58 | 3.73 | 0.13 |  $-$1.18 | 0.15
26016-68 | 255.7539 | 22.7623 | 6467 | 13 | 3.59 | 0.07 |  $-$1.63 | 0.06
26016-70 | 255.3743 | 22.7649 | 6137 | 41 | 3.54 | 0.14 |  $-$1.62 | 0.10
26016-74 | 255.5894 | 22.7630 | 6162 | 14 | 3.87 | 0.06 |  $-$1.04 | 0.05
26016-84* | 254.8371 | 22.7663 | 5532 | 151 | 3.84 | 0.04 |  $-$0.76 | 0.09

888ID stands for the J-PLUS ID. RA and Dec. are the right ascension and the
declination in degree (J2000). The extrapolation stars have an asterisk on
their ID, which indicates that their stellar parameters are not very reliable.
The $\Delta$ shows the uncertainty.

Both median and mean values of the uncertainties are presented in Table 4. The
difference in mean versus median shows the existence of a few stars with a
large residual that raise the mean in both the interpolation and extrapolation
catalogs. These uncertainties were decreased by about $3\%\sim 5\%$ by
applying the contours to constrain the prediction sample. The uncertainties of
the prediction are smaller than the residuals among different pipelines,
implying that our regressions have fairly good precisions (see Sect. 5.3).
This also implies that the uncertainties are mainly generated from the
pipeline difference and observation uncertainties.

Table 4: Uncertainties of the classification interpolation and the
extrapolation catalogs

| Inter median | Inter mean | Extra median | Extra mean
---|---|---|---|---
$T_{\rm eff}$ (K) | 14.0 | 34.7 | 59.5 | 71.6
log $g$ (dex) | 0.0689 | 0.0780 | 0.0268 | 0.0515
$\rm[Fe/H]$ (dex) | 0.0566 | 0.0677 | 0.0187 | 0.0477
999Inter stands for the classification interpolation catalog and Extra stands
for the classification extrapolation catalog. We present both median and mean
values in this table.

##  5 Discussion

###  5.1 Training and blind test

We tested three methods to regress the parameters. In the first method, we
trained a single SVR model using the magnitude uncertainties as weights, but
excluding the stellar parameter uncertainties. This method aims at making a
comparison with our multi-model simulation. The weight is given by
$w_{i}=lg\sum_{j=1}^{1}2(\frac{1}{e_{ij}^{2}})$ for an object $x_{i},$ where
the $e_{ij}$ shows the magnitude uncertainty of mag$j$. We adopted a
logarithmic scale to reduce potential steep cliffs in the feature space (Wang
et al. 2021). In the second and third tests, we randomly constructed 80
training sets for each stellar parameter. The difference between the two tests
is whether LAMOST MRS were included in the training sample. We constructed the
second test to show the generalization ability of extrapolating data. The
first and second tests use LAMOST MRS as a blind test, while the third one
includes all the stars in the training sample and holds $10\%$ for the blind
test. The third method is the final model we applied. It has lower RMSEs than
the single model. We present the results of these three tests in Table 5.

Table 5: RMSEs of the blind tests in the three tests  Method |  $T_{\rm eff}$ (K) | log $g$ (dex) | [Fe/H] (dex)
---|---|---|---
Single Model | 239.9 | 0.373 | 0.296
Generalization Ability | 239.6 | 0.400 | 0.313
Final Model | 159.6 | 0.345 | 0.250

###  5.2 Multi-model approach

Table 6: Gaussian fits of the residuals for our blind tests and other works.

Parameter |  $\mu$(BT) |  $\Delta$(BT) |  $\Delta$(Train) |  $\Delta$(Bai) |  $\Delta$(Yang) |  $\mu$(Galarza) |  $\Delta$(Galarza)
---|---|---|---|---|---|---|---
$T_{\rm eff}$ (K) | 2.3125 | 159.5216 | 84.6750 | 100 $\sim$ 200 |  $\sim$ 55 | 41 | 61
log $g$ (dex) | 0.0520 | 0.3413 | 0.1200 |  | 0.15 | 0.11 | 0.17
$\rm[Fe/H]$ (dex) | 0.0366 | 0.2476 | 0.0660 |  | 0.07 | 0.09 | 0.14

101010BT stands for the blind test. Bai, Yang, and Galarza present the
uncertainties from Bai et al. 2018, Yang et al. 2021, and Galarza et al. 2022,
respectively. Bai et al. 2018 only contains effective temperature regression.
$\mu$ is the center of the Gaussian fit and $\Delta$ is the standard
deviation. We present the median values of the training uncertainties.

Bai et al. (2018) constructed a single model to train the stellar parameters
and enlarge the sample to constrain the discrepancy caused by different
pipelines. The uncertainty of the stellar parameter is not included in their
calculation. Yang et al. (2021) also constructed one deep learning artificial
neural network (ANN) — $CSNet$ — and gained stellar parameters with a
precision of $\sim 55\rm K$ in effective temperature, 0.15 dex in log $g,$ and
0.07 dex in [Fe/H]. Galarza et al. (2022) presented a multi-regressior
pipeline, SPEEM, based on three single XGB algorithms and validated the result
with SSPP. The errors for each stellar parameter ($T_{\rm eff}$, log $g$,
[Fe/H]) are $41\pm 61$ K, $0.11\pm 0.17$ dex, and $0.09\pm 0.14$ dex,
respectively. More information is presented in Table 6. The precision of all
these single models is similar to our results.

In our study, the RMSE of a single model is similar to the RMSEs of 80 models,
indicating that the multi-model approach performs similarly to the weight-
controlled model. We can incorporate the uncertainties of both features and
target parameters, which is the main advantage of our approach. Furthermore,
if a certain photometric measurement suffers a higher uncertainty, its
contribution to the regression would decrease significantly, independently of
the uncertainties in the remaining 11 measurements. This could introduce
biases into the models and result in inaccurate predictions. Such bias is
avoided in our results, since all the uncertainties of the 12-band magnitudes
are included in the model construction.

We also find that the uncertainties of the stellar parameters are larger than
the blind test’s RMSEs. If these uncertainties are not considered in the
calculation, the error from the spectral fitting pipeline would not propagate
to the final regression. The uncertainties of the prediction would be highly
underestimated in this case.

###  5.3 Pipelines

Different pipelines result in different stellar parameters (Fig. LABEL:Fig56,
see more in Appendix C), and we adopted all parameters from all pipelines. We
selected the stars that have parameters from multiple spectroscopy surveys and
show the differences and variances in Table 7. In Table 6 and Table 7, the
differences in the blind test are smaller than those caused by pipelines. When
the stellar parameters from different pipelines are included in our training,
the regression models are unbiased among these different pipelines and can
thus be used to provide general predictions from different surveys. Including
a larger variety of surveys in the training data will improve the
generalization ability.

Table 7: Differences among pipelines

Stellar parameter | pipelines | overlap | $\mu$ | $\sigma$
---|---|---|---|---
| L$-$A | 6,956 | 27.8 | 143.2
$T_{\rm eff}$ | L$-$S | 3,318 |  $-$88.4 | 143.1
| S$-$A | 64 |  $-$112.3 | 149.7
| L$-$A |  | 0.048 | 0.191
log $g$ | L$-$S |  | 0.206 | 0.353
| S$-$A |  |  $-$0.150 | 0.228
| L$-$A |  | 0.008 | 0.099
${\rm[Fe/H]}$ | L$-$S |  | 0.068 | 0.184
| S$-$A |  |  $-$0.111 | 0.144
111111The mean and standard variance of the difference of pipelines. In this
table, the ”L” stands for LAMOST, while ”S” for SSPP and ”A” for APOGEE.
Overlap stands for the number of cross-match overlapping between these
catalogs. The unit of effective temperature is in K, and dex for log $g$ and
[Fe/H].

Other studies usually adopt data from one survey to avoid biases among
different pipelines (Bu et al. 2020; Lu & Li 2015). Another method is to set a
priority for each catalog. Some studies have concluded that a larger sample
size can decrease the bias caused by pipelines, for example Bai et al. (2018).

Although a single catalog and priority-based methods do not suffer the
systematic error caused by pipelines, the models built from them probably do
not perform well for other surveys. Our training sample contains different
surveys, and their biases propagate to the final regression models. Therefore,
our models would have a wider application.

One catalog-based method may be good for regression of the same catalog, but
J-PLUS is a more general catalog with a large sky coverage and many
photometric wavebands. J-PLUS requires models from hybrid samples for general
applications.

There may not exist a regressor that can properly fit the data for all the
surveys; this is analogous to the no-free-lunch theorem, which says that a
machine-learning algorithm that can solve every problem does not exist
(Shalev-Shwartz & Ben-David 2014). This implies that the prediction would be
biased if we use single catalog or priority-based methods to regress the data
in another catalog. In machine learning, a larger amount of data does not
always work, but higher diversity does (Wang et al. 2021).

###  5.4 Distribution of stellar parameters

While our data augmentation process employs Gaussian-based uncertainty
sampling, the predicted stellar parameters may not necessarily follow a
Gaussian distribution. The SVR embeds the samples to a higher-dimensional
feature space by using a nonlinear embeding mapping. The nonlinear mapping may
not hold the distribution of sample to the feature space.

To quantify whether our approach preserves the shape of the distribution, we
applied the Lilliefors test (Lilliefors 1967) to the blind test samples. About
$19.01\%$ of the $T_{\rm eff}$, $19.49\%$ for log $g,$ and $19.39\%$ for
[Fe/H] pass the test with a significance level of 0.05. When the level
decreases to 0.01, the number increases to $39.94\%$ of the $T_{\rm eff}$,
$41.40\%$ for log $g$ and $40.54\%$ for [Fe/H]. Given the low number
statistics in our predicted distributions of only 80 samples, it is possible
that this test underestimates the true fraction of Gaussian distributed
predictions. To estimate the influence of the number of samples, we
constructed 27,970 sets (which was the same size as the blind test size) of
random numbers based on $N(100,25)$ and find that about $94.87\%$ passed the
test with a significance level of 0.05. When we change the number of sets to
800, the result increases to $95.11\%$. The difference is less than $1\%$.
Therefore, there is no evidence that the resulted stellar parameters follow a
Gaussian distribution.

###  5.5 Comparison with Yang et al. (2021)

To determine the robustness of our predictions, we performed a comparison with
the stellar parameter catalog from Yang et al. (2021) 121212available at
http://www.j-plus.es/ancillarydata/dr1_stellar_parameters_elemental_abundances.
We selected the reliable stars in both catalogs, which were stars in the
interpolation catalog and the stars with stellar parameter FLAG=0 in Yang et
al. (2021). The cross-match yielded 2,008,654 stars. The average stellar
parameter differences are -247.9764 K for $T_{\rm eff}$, -0.1984 dex for log
$g$, and -0.0998 dex for [Fe/H]. There are some extreme values in Yang et al.
(2021), and some of them are unreliable. For example, the maximum and minimum
value of surface gravity are 882.44 dex and -237.24 dex and 460.39 dex and
-567.16 dex for [Fe/H]. These stars may be situated at the edge of the feature
space and, thus, are subject to overfitting by the ANN model.

The training sample ranges of Yang et al. (2021) are about 4000 $\sim$ 7500 K
for $T_{\rm eff}$, 0 $\sim$ 5 dex for log $g$, and -3 $\sim$ 1 dex for [Fe/H].
We restricted their results to these ranges and obtain 1,663,053 stars. The
average differences decrease to -63.6817 K for $T_{\rm eff}$, 0.0886 dex for
log $g$, and -0.1069 dex for [Fe/H] (Fig. 4). These differences are similar to
our blind test uncertainties. The main difference between our work and Yang et
al. (2021) is the control of the extrapolations. These extreme values might
fall into the extrapolation category and suffer from overfitting. The training
sample of Yang et al. (2021) is not an official LAMOST release, and the
training sample difference may also be a minor reason for such a stellar
parameter difference.

Figure 4: Comparison with Yang et al. (2021). The three panels are $T_{\rm
eff}$, log $g$, and [Fe/H], from top to bottom. Only the density contours are
shown in these figures.

##  6 Conclusion

In this work, we predicted stellar parameters, the effective temperature
$T_{\rm eff}$, the surface gravity log $g,$ and the metallicity [Fe/H] for
J-PLUS DR1 by using $3\times 80$ SVR models. We chose stars from spectra-based
surveys (LAMOST, APOGEE, and SEGUE) to construct our training sample in order
to improve the reliability of the sample and gained 279,702 stars from the
cross-match. We normalized the features of these stars by subtracting their
average in order to accelerate the calculation. We held out ten percent of the
set for the blind test and the remaining 251,732 stars were used for training.
We used 12 three-dimensional density contours to distinguish the prediction
reliability for a given star. We presented two catalogs that correspond to the
classification catalogs of Wang et al. (2021). For the classification
interpolation catalog, we present 1,898,154 stars with interpolated stellar
parameters and 595,270 stars with extrapolations in our regression. For their
classification extrapolation catalog, on the other hand, we present 13,274
interpolated and 220,650 extrapolated predictions with our SVR. Regarding the
construction of the models, multi-model simulations give better results than
the single model. The RMSEs of our validations are 159.6 for $T_{\rm eff}$,
0.3453 for log $g,$ and 0.2502 for $[Fe/H]$. Using different catalogs as a
sample set may increase the generalization ability of our model. Lastly, we
compared our results to Yang et al. (2021) and find a decent agreement with
their work, with discrepancies comparable to the RMSEs that we found in our
blind test.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(NSFC) through grants NSFC-11988101/11973054/11933004 and the National
Programs on Key Research and Development Project (grant No.2019YFA0405504 and
2019YFA0405000). Strategic Priority Program of the Chinese Academy of Sciences
under grant number XDB41000000. H, B, Yuan, acknowledge funding from National
Natural Science Foundation of China (NSFC) through grants NSFC-12173007. JAFO
acknowledges the financial support from the Spanish Ministry of Science and
Innovation and the European Union – NextGenerationEU through the Recovery and
Resilience Facility project ICTS-MRR-2021-03-CEFCA. F.J.E. acknowledges
financial support by the Spanish grant PGC2018-101950-B-I00 and MDM-2017-0737
at Centro de Astrobiología (CSIC-INTA), Unidad de Excelencia María de Maeztu.
We acknowledge the science research grants from the China Manned Space Project
with NO.CMS-CSST-2021-B07.
Based on observations made with the JAST80 telescope at the Observatorio
Astrofísico de Javalambre (OAJ), in Teruel, owned, managed, and operated by
the Centro de Estudios de Física del Cosmos de Aragón (CEFCA). We acknowledge
the OAJ Data Processing and Archiving Unit (UPAD) for reducing the OAJ data
used in this work. Funding for the J-PLUS Project has been provided by the
Governments of Spain and Aragón through the Fondo de Inversiones de Teruel;
the Aragón Government through the Research Groups E96, E103, and E16_17R; the
Spanish Ministry of Science, Innovation and Universities (MCIU/AEI/FEDER, UE)
with grants PGC2018-097585-B-C21 and PGC2018-097585-B-C22; the Spanish
Ministry of Economy and Competitiveness (MINECO) under AYA2015-66211-C2-1-P,
AYA2015-66211-C2-2, AYA2012-30789, and ICTS-2009-14; and European FEDER
funding (FCDD10-4E-867, FCDD13-4E-2685). The Brazilian agencies FINEP, FAPESP,
and the National Observatory of Brazil have also contributed to this project.
Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic
Telescope LAMOST) is a National Major Scientific Project built by the Chinese
Academy of Sciences. Funding for the project has been provided by the National
Development and Reform Commission. LAMOST is operated and managed by the
National Astronomical Observatories, Chinese Academy of Sciences. Funding for
the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan
Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High-Performance Computing at the University of Utah. The SDSS
website is http://www.sdss.org/. SDSS-IV is managed by the Astrophysical
Research Consortium for the Participating Institutions of the SDSS
Collaboration including the Brazilian Participation Group, the Carnegie
Institution for Science, Carnegie Mellon University, the Chilean Participation
Group, the French Participation Group, Harvard-Smithsonian Center for
Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins
University, Kavli Institute for the Physics and Mathematics of the Universe
(IPMU)/University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz
Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie
(MPIA Heidelberg), Max- Planck-Institut für Astrophysik (MPA Garching), Max-
Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatário Nacional/MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University.

## References

  * Allende Prieto et al. (2008) Allende Prieto, C., Sivarani, T., Beers, T. C., et al. 2008, AJ, 136, 2070
  * Anguiano, B. et al. (2018) Anguiano, B., Majewski, S. R., Allende-Prieto, C., et al. 2018, A&A, 620, A76
  * Awad & Khanna (2015) Awad, M. & Khanna, R. 2015, Support Vector Regression (Berkeley, CA: Apress), 67–80
  * Bai et al. (2019) Bai, Y., Liu, J., Bai, Z., Wang, S., & Fan, D. 2019, The Astronomical Journal, 158, 93
  * Bai et al. (2018) Bai, Y., Liu, J., Wang, S., & Yang, F. 2018, The Astronomical Journal, 157, 9
  * Bailer-Jones (2011) Bailer-Jones, C. A. L. 2011, Monthly Notices of the Royal Astronomical Society, 411, 435
  * Boeche et al. (2018) Boeche, C., Smith, M. C., Grebel, E. K., et al. 2018, The Astronomical Journal, 155, 181
  * Boser et al. (1992) Boser, B. E., Guyon, I. M., & Vapnik, V. N. 1992, in Proceedings of the Fifth Annual Workshop on Computational Learning Theory, COLT ’92 (New York, NY, USA: Association for Computing Machinery), 144–152
  * Breiman (2001) Breiman, L. 2001, Statist. Sci., 16, 199
  * Bu et al. (2020) Bu, Y., Kumar, Y. B., Xie, J., et al. 2020, The Astrophysical Journal Supplement Series, 249, 7
  * Cenarro et al. (2019) Cenarro, A. J., Moles, M., Cristóbal-Hornillos, D., et al. 2019, Astronomy & Astrophysics, 622, A176
  * Cenarro et al. (2014) Cenarro, A. J., Moles, M., Marín-Franch, A., et al. 2014, in Proc. SPIE, Vol. 9149, Observatory Operations: Strategies, Processes, and Systems V, 91491I
  * Cortes & Vapnik (1995) Cortes, C. & Vapnik, V. 1995, Machine Learning, 20, 273
  * Cover & Hart (1967) Cover, T. & Hart, P. 1967, IEEE Trans. Inf. Theory, 13, 21
  * Cristianini & Shawe-Taylor (2000) Cristianini, N. & Shawe-Taylor, J. 2000, An Introduction to Support Vector Machines and Other Kernel-based Learning Methods (Cambridge University Press)
  * Deng et al. (2012) Deng, L.-C., Newberg, H. J., Liu, C., et al. 2012, Research in Astronomy and Astrophysics, 12, 735
  * Drucker et al. (1997) Drucker, H., Burges, C. J., Kaufman, L., et al. 1997, Advances in neural information processing systems, 9, 155
  * Galarza et al. (2022) Galarza, C. A., Daflon, S., Placco, V. M., et al. 2022, A&A, 657, A35
  * García Pérez et al. (2016) García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144
  * Jiménez-Teja et al. (2019) Jiménez-Teja, Y., Dupke, R. A., Lopes de Oliveira, R., et al. 2019, Astronomy & Astrophysics, 622, A183
  * Jönsson et al. (2018) Jönsson, H., Prieto, C. A., Holtzman, J. A., et al. 2018, The Astronomical Journal, 156, 126
  * Lee et al. (2008a) Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008a, AJ, 136, 2022
  * Lee et al. (2008b) Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008b, AJ, 136, 2050
  * Lilliefors (1967) Lilliefors, H. W. 1967, Journal of the American Statistical Association, 62, 399
  * Liu et al. (2020) Liu, C., Fu, J., Shi, J., et al. 2020, LAMOST Medium-Resolution Spectroscopic Survey (LAMOST-MRS): Scientific goals and survey plan
  * Lu & Li (2015) Lu, Y. & Li, X. 2015, Monthly Notices of the Royal Astronomical Society, 452, 1394
  * MacKay (2003) MacKay, D. J. 2003, Information theory, inference and learning algorithms (Cambridge university press)
  * Marín-Franch et al. (2015) Marín-Franch, A., Taylor, K., Cenarro, J., Cristobal-Hornillos, D., & Moles, M. 2015, in IAU General Assembly, Vol. 29, 2257381
  * Nogueira-Cavalcante et al. (2019) Nogueira-Cavalcante, J. P., Dupke, R., Coelho, P., et al. 2019, Astronomy & Astrophysics, 630, A88
  * Quinlan (1986) Quinlan, J. R. 1986, Machine Learning, 1, 81
  * Ruder (2017) Ruder, S. 2017, An overview of gradient descent optimization algorithms
  * Shalev-Shwartz & Ben-David (2014) Shalev-Shwartz, S. & Ben-David, S. 2014, Understanding Machine Learning: From Theory to Algorithms (Cambridge University Press)
  * Sichevskij (2012) Sichevskij, S. G. 2012, Astronomy Reports, 56, 710
  * Sichevskij et al. (2014) Sichevskij, S. G., Mironov, A. V., & Malkov, O. Y. 2014, Astrophysical Bulletin, 69, 160
  * Smola & Schölkopf (2004) Smola, A. J. & Schölkopf, B. 2004, Statistics and Computing, 14, 199
  * Steinmetz et al. (2020) Steinmetz, M., Matijevič, G., Enke, H., et al. 2020, The Astronomical Journal, 160, 82
  * Stone (1977) Stone, C. J. 1977, The Annals of Statistics, 5, 595
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Wang et al. (2021) Wang, C., Bai, Y., Yuan, H., Wang, S., & Liu, J. 2021, Machine Learning Applied to STAR-GALAXY-QSO Classification of The Javalambre-Photometric Local Universe Survey
  * Wang et al. (2020) Wang, J., Fu, J.-N., Zong, W., et al. 2020, The Astrophysical Journal Supplement Series, 251, 27
  * Whitten et al. (2019) Whitten, D. D., Placco, V. M., Beers, T. C., et al. 2019, Astronomy & Astrophysics, 622, A182
  * Wu et al. (2014) Wu, Y., Luo, A., Du, B., Zhao, Y., & Yuan, H. 2014, Automatic stellar spectral parameterization pipeline for LAMOST survey
  * Wu et al. (2011) Wu, Y., Luo, A. L., Li, H.-N., et al. 2011, Research in Astronomy and Astrophysics, 11, 924
  * Yang et al. (2021) Yang, L., Yuan, H., Xiang, M., et al. 2021, arXiv e-prints, arXiv:2112.07304
  * Yanny et al. (2009) Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, The Astronomical Journal, 137, 4377
  * Yuan (2021) Yuan, H. 2021, In preparation and private communication
  * Yuan et al. (2015) Yuan, H., Liu, X., Xiang, M., et al. 2015, The Astrophysical Journal, 799, 133
  * Zasowski et al. (2013) Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, The Astronomical Journal, 146, 81

##  Appendix A Sample

We present the distribution of a sample set in this section. The standard
sample set distributions are shown in red blocks. We then split the sample set
into 80 new randomly constructed sample sets to perform our simulation. The
distribution contains all constructed objects of the 80 sets and is presented
in the white block. The overlapping part of these blocks will fade to pink,
and we can see that these figures are nearly all in pink.

Figure 5: Stellar parameter distributions of the sample set. The red block
shows the distribution of standard sample set. The white block shows the
simulation distribution. The overlap is shown in pink (red+white). The y axis
shows the probability.

Figure 6: Magnitude distributions of the sample set. The yellow block shows
the J-PLUS distribution.

##  Appendix B Blind test result

We plot the results of the blind test, showing the differences between the
predicted and the true parameters. We fit the residuals with Gaussian
functions.

Figure 7: Results of the blind test validation. The upper two panels are
$T_{\rm eff}$. The middle two panels are log $g$, and the lower two panels are
[Fe/H].

##  Appendix C Pipeline differences

We present the differences among the pipelines. The residuals are fitted by
Gaussian functions, and the means and variances are shown in Table 7.

Figure 8: Pipeline differences of $T_{\rm eff}$. The upper two panels show the
difference between LAMOST and APOGEE. The middle two are the difference
between LAMOST and SSPP, and the lower two are the difference between SSPP and
APOGEE.

Figure 9: Pipeline difference of log $g$.

Figure 10: Pipeline difference of [Fe/H].


# Femtoscopy analysis of ultra-soft pion trap at LHC energies

W. Rzesa1, G. Kornakov1, A.R. Kisiel1, Yu.M. Sinyukov1,2, V.M. Shapoval2 1
Warsaw University of Technology, Faculty of Physics, ul. Koszykowa 75, 00-662,
Warsaw, Poland
2 Bogolyubov Institute for Theoretical Physics, 14b Metrolohichna street, Kyiv
03143, Ukraine

###### Abstract

Femtoscopy studies of pion radiation in heavy-ion collisions have been
conducted extensively at all available collider energies, both theoretically
and experimentally. In all these studies a special interest is given to
$m_{T}$ dependency of pion femtoscopy radii, usually approximated by a power-
law function at transverse momenta above 200 MeV/$c$. However, the radii
behaviour has been much less explored for the ultra-soft pions, possessing the
transverse momentum comparable to or lower than the pion mass. For many
experimental setups this region is difficult to measure. In this work we
present theoretical calculations of pion emission in the ultra-soft region in
the two hybrid models — iHKM and LHYQUID+THERMINATOR2. Along with the particle
transverse momentum spectra, we present the calculated femtoscopy radii, both
in one-dimensional and three-dimensional representations. We investigate the
radii dependence on pair $m_{T}$ and observe, in particular, a departure from
the power-law behaviour at ultra-soft momenta, potentially reflecting a
decoupling of such slow pions from the rest of collectively expanding system.
We provide the theoretical interpretation of this result and discuss its
significance, in particular, for the ongoing non-identical particle femtoscopy
analysis for pairs consisting of a pion and a baryon (or of a pion and a
charmed meson).

##  I INTRODUCTION

The interferometry of particles emitted from the small radiating source
enables one to study the evolving geometry of that source. The particle
interferometry has originated from the experimental paper [1] and theoretical
works [2, 3]. In the latter it was also found that the results for the quantum
statistics correlation effect for two close sources of identical particle
emission in collision physics are formally similar to those for classical
intensity correlations of the light radiation coming from different parts of
the star, known in astronomy as the Hanbury Brown – Twiss (HBT) effect [4, 5,
6]. The similarity, as well as the difference between the quantum correlations
in identical-particle pairs and the classical HBT effect are analyzed, e.g.,
in [7]. Note, that the correlation femtoscopy method, initially developed for
finding the overall size of identical quantum particles’ source, now is
considered to measure the homogeneity lengths [8] of such femto-small sources.
This general interpretation is especially important for the studies of the
space-time structure of expanding sources, such as those created in
proton/antiproton and heavy-ion collisions [1, 2, 3, 9, 10, 11, 12, 13, 14].

The homogeneity lengths, defined as a result of correlation femtoscopy
analysis, reflect the spatial dimensions of the region within the entire
strongly interacting system formed in the collision, from which particles are
emitted with similar velocities (having close, nearly coinciding values and
directions) [8, 15]. Gaussian parameterization for the particle emission
source function (SF) and the two-particle correlation function (CF) are
usually used in femtoscopy studies, although the realistic source shape
differs from a perfect Gaussian — e.g., resonance decay contributions cause
the exponential behaviour near the peak of the correlation function. The
conventional use of the Gaussian shape allows one to standardize the
description of experimental data and easily compare the results of different
femtoscopy measurements, as well as to interpret the obtained radii as the
homogeneity lengths. It is also well-motivated experimentally in heavy-ion
collisions, where correlation shapes for pions in 3D analyses are universally
observed to be Gaussian in a wide range of centralities and pair transverse
momenta. The study of the dependence of the radii on the particle species,
collision type and collision energy is the main objective of the femtoscopic
analysis of heavy-ion collisions.

In experimental studies, one clearly observes certain types of universal
scaling behaviour for the measured femtoscopy scales. For a given colliding
system and collision energy a linear scaling of the femtoscopic radii is
universally observed versus cube root of the final state mean particle
multiplicity $\langle dN_{ch}/d\eta\rangle^{1/3}$ [16, 17, 18, 19, 20, 21,
22]. A similar scaling across collision energies and colliding systems is only
approximate [23, 24]. The radii versus the pair transverse momentum exhibit a
power-law-like scaling in pair transverse mass $m_{T}$ ($m_{T}$ =
$\sqrt{k_{T}^{2}+m^{2}}$, where $k_{T}$ is pair mean transverse momentum,
$k_{T}=|\textbf{p}_{T1}+\textbf{p}_{T2}|$/2, and $m$ is particle mass) [25,
26, 27, 28, 29, 16, 17, 30]. In heavy-ion collisions both these scalings are
predicted by hydrodynamic models [31]. Specifically, the $m_{T}$ scaling is
explained as a direct consequence and one of the main signatures of the
collective radial flow of the system.

The observed $m_{T}$ power-law scaling means a decrease of the radii with
growing $m_{T}$ that is a signature of hydrodynamic collectivity, typical for
all particle species affected by the same flow field. Such dependence can be
helpful for the prediction of the source size at certain $m_{T}$ through the
interpolation between radii measured at different $m_{T}$ values. Nonetheless,
the $m_{T}$ range of experimental radii measurement depends on the acceptance
of the detectors used in the study, and the very-low-$m_{T}$ region has not
been reached yet via such measurements. Therefore, the femto-scale predictions
in that region can come only from the extrapolation, which following the
commonly used power-law function, shows a very rapid growth of the radii with
decreasing $m_{T}$. The theoretical studies, by contrast, do not face such
problems and so, simulations within a realistic collision model can help
estimate the possible character of the femtoscopy radii behaviour in ultra-
soft momentum region.

The knowledge of femtoscopic radii for ultra-soft pions is becoming
increasingly important in the context of rapidly developing experimental
analyses of femtoscopic correlations of non-identical particles [32, 33]. In
the case of such pairs the correlation occurs between two particles with
similar velocity. If the pair in question contains a pion and a relatively
massive particle (e.g. a deuteron or a charmed meson), then the pion in the
pair needs to be ultra-soft in order to be correlated. The femtoscopic radius
of the pair is directly related to the size of emission regions of both
particles in the pair. Therefore, in order to study the non-identical
correlations between pions and massive particles, the knowledge of femtoscopic
radii of ultra-soft pions becomes essential.

Several previous theoretical studies of the ultra-low momentum pions were
motivated by the experimental evidence of enhanced production [34, 35]. These
observations could be addressed by non-equilibrium models [36]. In addition to
the semi-classical models, pure quantum effects such as pion condensation and
Bose-Einstein enhancement were also explored [37]. However, none of these
studies addressed the pion femtoscopic radii nor their scaling at very low
transverse momenta.

Thus, in this paper, we focus on studying the pion femtoscopic radii
dependence on pair $m_{T}$, including the region of very low transverse
masses, in Pb-Pb collisions at the LHC energy $\sqrt{s_{NN}}$ = 5.02 TeV
simulated within the integrated HydroKinetic Model (iHKM) [31, 38, 39] and the
hydrodynamics model [40] coupled to the statistical hadronisation code
THERMINATOR2 (LHYQUID + THERMINATOR2) [41]– (LQTH).

##  II  MODELS’ DESCRIPTION

The collisions of heavy ions at ultra-relativistic energies allow to create
quark-gluon plasma (QGP) — a new state of strongly interacting matter, where
quarks and gluons are no longer confined within individual nucleons. The
evolution of such matter can be successfully described in complicated models
realized in the form of codes for labor-consuming computer simulations. Such
codes typically include modules describing the initial state formation,
relativistic hydrodynamics expansion of liquid-like QGP at the intermediate
stage of the system’s evolution, its subsequent particlization (turning into a
set of hadrons) and, finally, the hadron-resonance gas expansion. As a result
of typical calculation, one obtains from the model a set of created hadrons,
characterized by the space-time points of their last collision and 4-momenta.
Based on these data one can construct different observables, like spectra,
correlation functions etc.

The data analyzed in this study were generated using two such models, iHKM and
LQTH, each having its own characteristics.

The iHKM is one of the most complete models, describing all the essential
phases of the matter evolution in the course of a relativistic A+A collision
(see [31, 38] for details). The initial conditions (IC) for each simulation
include the energy density spatial distribution at the starting time
$\tau_{0}$, usually close to 0.1 fm/$c$ (for the high-energy collisions we use
GLISSANDO code [42] to generate it), and anisotropic momentum distribution,
inspired by the Color Glass Condensate model. These IC correspond to the very
initial non-equilibrium partonic state right after the two nuclei collision.
At the next stage of the system’s evolution, it gradually thermalizes and
approaches to a nearly hydrodynamical state. Then (starting from the
thermalization time $\tau_{th}\approx 1$ fm/$c$) follows a continuous medium
expansion described in the Israel-Stewart viscous hydrodynamics formalism. As
the matter expands and cools down, it eventually reaches the “particlization
temperature” $T_{p}\approx 160$ MeV (depending on the QGP equation of state
used at the hydrodynamics stage), when one switches to the description of the
system in terms of hadrons (here the Cornelius routine [43] is applied).
Produced hadrons are then fed to the UrQMD hadron cascade code [44, 45],
performing multiple hadronic re-scatterings and resonance decays taking place
at this final “afterburner” phase of the collision. The model is calibrated
based on the experimental mean charged particle multiplicity and pion $p_{T}$
spectrum slope in the most central collisions of a given type. In this paper,
we use the iHKM tuning, used for the LHC $5.02A$ TeV Pb-Pb collisions
simulation described in detail in [39].

The LQTH approach assumes the ideal (3+1)D hydrodynamics evolution of the
created fireball (however, bulk and shear viscosity can be implemented as
well) until the single (both chemical and kinetic) freeze-out at certain
temperature ($T_{p}\approx 140$ MeV), where statistical hadronization takes
place, and the system’s chemical composition and momentum spectra do not
change anymore, except for the resonance decays. As well as the iHKM, for the
current study the model is properly calibrated to describe the
$\sqrt{s_{NN}}=5.02$ TeV Pb-Pb collisions at the LHC.

Both models describe well the experimental results for pion, kaon and proton
transverse momentum spectra, including the region of very soft momenta (see
Fig. 1). Note, that such a description is achieved without the implementation
of additional specific mechanisms of the soft pion emission, like the Bose-
Einstein condensation, discussed, e.g., in [37]. Also, the two utilized models
were previously used to successfully describe the LHC data on $2.76A$ TeV Pb-
Pb collisions [38, 39, 46, 47].

Figure 1:  Pion, kaon and proton transverse momentum spectra calculated in the
iHKM (solid lines) and in the LQTH models (dashed lines) for the most central
$\sqrt{s_{NN}}=5.02$ TeV Pb-Pb collisions at the LHC [48].

##  III  CONSTRUCTION OF THE CORRELATION FUNCTION

The femtoscopic radii constituting the main subject of this study are
typically defined from the Gaussian fits to the correlation functions
depending on pair relative momentum $q=p_{1}-p_{2}$. Here we consider the
correlations in pairs of identical charged pions ($\pi^{+}$ and $\pi^{-}$),
performing 1D and 3D analysis. The 1D correlation of identical particles is
calculated as a function of $q_{\mathrm{inv}}=\sqrt{q_{0}^{2}-\textbf{q}^{2}}$
in the Pair Rest Frame (PRF) and requires minimal statistics in contrast to
the 3D study. However, the radius $R_{\mathrm{inv}}$, in this case, is the
only length extracted, so that the emission source is assumed to be a
spherically symmetric Gaussian one. The 3D study is based on relative momentum
components $q_{\mathrm{out}}$, $q_{\mathrm{side}}$, and $q_{\mathrm{long}}$
calculated in the Longitudinally Co-Moving System (LCMS), where the
longitudinal direction is along the beam axis, the outward direction is along
the pair transverse momentum and the sideward direction is perpendicular to
the other two. Accordingly, the three radii $R_{\mathrm{out}}$,
$R_{\mathrm{side}}$, and $R_{\mathrm{long}}$ are extracted from the CF fit in
this case.

The correlation function is generally defined as

$C(p_{1},p_{2})=\frac{P_{12}(p_{1},p_{2})}{P_{1}(p_{1})P_{2}(p_{2})},$ (1)

and can be understood as the ratio of conditional probability to detect a pair
of particles with the specific momentum values, $P_{12}(p_{1},p_{2})$, to the
probability of finding them with such momenta independently,
$P_{1}(p_{1})P_{2}(p_{2})$.

In experimental analysis, the procedure of obtaining the correlation function
consists of making pairs of particles coming from the same event to build the
correlated distribution (numerator), and for the background distribution
(denominator) pairs are created using particles coming from different events
(meaning that they retain all aspects of experimental acceptance, while not
being correlated due to the mutual interaction).

In simulation studies, one usually fills two histograms (1D or 3D, depending
on the analysis type) representing the corresponding particle pairs’ relative
momentum $q$ distribution — one histogram for the numerator and another one
for the denominator in Eq. (1) — based on the generated model output (set of
hadrons). However, since in all currently available heavy-ion collision models
the quantum statistical effects (including the Bose-Einstein symmetrization of
the two identical bosons’ wave-function) are not implemented directly, on a
microscopic level, one usually introduces a weight factor

$w=1+\cos{(qr)},$ (2)

coming from the Bose-Einstein interference at the so-called “afterburner”
stage (to simplify the presentation, Coulomb and strong final-state
interactions are not considered here). This factor $w$ is taken as a weight
for the pairs entering the correlated distribution histogram (a numerator
one). For the background distribution (denominator) the weight is equal to 1.
In this study, we also use the rapidity cut, $|y|<1$, during the identical
charged pion pairs selection to reproduce the acceptance of the ALICE
detector.

To obtain a fitting formula for the CF, which will allow us to define the
desired femtoscopy radii, one can consider the 7D (if particles are on the
mass shell) particle emission functions $S_{i}(x,p)$ for each particle species
$i$ and derive the approximate expressions for single-particle and two-
particle momentum spectra, and then using Eq. (1) obtain the well known
Bertsch-Pratt [49, 50] representation for the correlation function of two
identical bosons in LCMS:

$\displaystyle
C(\textbf{k,q})=1+\lambda_{3D}(\textbf{k})\,\mathrm{exp}(-R^{2}_{\mathrm{out}}(\textbf{k})q^{2}_{\mathrm{out}}-$
$\displaystyle
R^{2}_{\mathrm{side}}(\textbf{k})q^{2}_{\mathrm{side}}-R^{2}_{\mathrm{long}}(\textbf{k})q^{2}_{\mathrm{long}}),$
(3)

where femto-radii $R_{i}$ depend on mean pair momenta k, and $\lambda_{3D}$
are referred to as correlation strength factors.

The invariant one-dimensional form of the correlation function Gaussian
parameterization is

$C(k,q_{\mathrm{inv}})=1+\lambda_{1D}(k)\exp{(-R^{2}_{\mathrm{inv}}(k)q^{2}_{\mathrm{inv}})}.$
(4)

##  IV  PION FEMTOSCOPY SCALES

This section presents the results on pion femto-radii obtained as a result of
correlation function fitting for identical pions in 1D and 3D studies using
the two models (iHKM and LQTH). The correlation functions were calculated
separately for the six centrality classes
($0-5\%,5-10\%,10-20\%,20-30\%,30-40\%$ and $40-50\%$) and the ten $k_{T}$
ranges
$(0.00-0.05,0.05-0.10,0.10-0.15,0.15-0.20,0.20-0.25,0.25-0.30,0.30-0.35,0.35-0.45,0.45-0.60,0.60-0.80)$
GeV/$c$.

Figure 2:  Femtoscopic radii of charged identical pions as functions of
$\langle dN_{ch}/d\eta\rangle^{1/3}$ calculated using iHKM (full markers) and
LQTH (empty markers) models, for 5 $k_{T}$ ranges (from top to bottom:
$0.00-0.05,0.10-0.15,0.20-0.25,0.30-0.35,0.45-0.60$ GeV/$c$). Different panels
from top to bottom correspond to 1D $R_{\mathrm{inv}}$, and 3D
$R_{\mathrm{out}}$, $R_{\mathrm{side}}$, and $R_{\mathrm{long}}$ radii
respectively. Lines correspond to linear fits to the radii dependencies.
Symbols are slightly shifted in the x direction for visibility. Figure 3:
Femtoscopic radii of charged identical pions as functions of $m_{T}$
calculated using iHKM (full markers) and LQTH (empty markers) models, for the
three centrality classes ($0-5\%$ — blue, $10-20\%$ — red, $30-40\%$ — green).
Different panels from top to bottom correspond to 1D $R_{\mathrm{inv}}$, and
3D $R_{\mathrm{out}}$, $R_{\mathrm{side}}$, and $R_{\mathrm{long}}$ radii
respectively.  Figure 4:  First-order derivative of the charged identical
pions femtoscopic radii $m_{T}$ dependency, $dR/dm_{T}$, in the HKM (left) and
the LQTH (right) models, for the three collision centralities (from top to
bottom: $0-5\%,10-20\%,30-40\%$). Each panel contains derivatives of the three
radii from the 3D study ($R_{\mathrm{out}}$ — blue, $R_{\mathrm{side}}$ — red,
$R_{\mathrm{long}}$ — orange) and $R_{\mathrm{inv}}$ from the 1D study
(green).

In Fig.2 all the radii for identical pion pairs, corresponding to 1D and 3D
calculations in the two models are demonstrated as functions of $\langle
dN_{ch}/d\eta\rangle^{1/3}$. The results for the five of all ten $m_{T}$ bins
are shown for better visibility since the radii in all bins follow a monotonic
trend. The radii values vary in the range from about 15 fm
($R_{\mathrm{long}}$ for the lowest $m_{T}$ and the highest multiplicity in
LQTH) to about 2 fm ($R_{\mathrm{side}}$ for high $m_{T}$ and low multiplicity
in iHKM). All the radii universally grow linearly with the cube root of final-
state multiplicity. The multiplicity dependencies for various $m_{T}$ ranges
generally show $m_{T}$ ordering, except for 1D and transverse 3D results from
the LQTH model, where the first bins of $m_{T}$ overlap. This is because in
the LQTH model, radii of the first $m_{T}$ bins take similar values. Fig. 3,
where femto-radii are shown as functions of $m_{T}$ presents it more clearly.
All the measured radii decrease with increasing $m_{T}$. However, in the case
of the LQTH model, we can observe a small plateau for the 1D and the
transverse 3D radii at low $m_{T}$. The iHKM results, however, mostly fall
monotonically in the lowest $m_{T}$ region. The Figure 3 also exhibits an
expected ordering with respect to analysed centrality — the more peripheral
collision, the smaller the radii. In both figures 2 and 3, the difference in
femto-scales between the two models is also evident.

One can see from Fig. 3 that the femto-radii in the iHKM model decrease faster
than in the LQTH one. This is likely due to a very early, at $\tau_{0}=0.1$
fm/$c$, expansion starts in the iHKM, leading to stronger velocity gradients
at the final stage. This is the main reason for the reduction of the femto-
radii: the homogeneity lengths are reduced with the growing flow intensity.
So, at the intermediate and high $k_{T}$ the iHKM radii are smaller, than
those in the LQTH model with a later start of the expansion ($\tau_{0}\approx
0.6$ fm/$c$). As for the radii values at ultra-soft transverse momenta, the
situation becomes even more complicated. At very small $k_{T}$ the transverse
radii are larger in the iHKM approach (the mentioned gradient could be the
reason), while the longitudinal ones are larger in the LQTH model. The physics
of pion production in the very-low-$k_{T}$ region appears to be non-trivial
and worthy of further study and better understanding.

Figure 5:  Pion femtoscopic radii dependencies on $m_{T}$ calculated in iHKM
(left) and LQTH (right) models, for the three centrality classes (from top to
bottom: $0-5\%,10-20\%,30-40\%$). Each panel contains results for the three
radii from the 3D study ($R_{\mathrm{out}}$ — squares, $R_{\mathrm{side}}$ —
crosses, $R_{\mathrm{long}}$ — tilted crosses) and for $R_{\mathrm{inv}}$ from
the 1D study (cicles). Lines correspond to power-law fits $am_{T}^{-b}$ to
radii dependencies in the range above 0.25 GeV/$c$.

To better investigate the character of the radii $m_{T}$ dependencies, they
were differentiated, and the obtained derivatives are presented in Fig. 4 (for
the three centrality classes). The derivative has been extracted in the range
$\pm 0.07$ GeV/$c$ around each point in the figure. Most of the functions show
similar behaviour, with a clear change in the derivative value at certain
momenta. The iHKM model shows a constant derivative at low momentum
(especially for the 1D study), then a systematic rise and again a higher
constant value at large momenta. Flat parts of the curves correspond to
constant slopes of the radii decrease approximated by a linear function. The
rising/falling parts indicate a change in the radii dependency slope. The
low-$m_{T}$ flat part ends around 0.25 GeV/$c$. The LQTH model, in turn, shows
for the 1D and the 3D transverse radii derivatives a constant or decreasing
behaviour up to 0.3 GeV/$c$, then the dependence starts rising to end up with
a constant value at large momentum. The initial decrease reflects the change
of a nearly constant $R(m_{T})$ dependence to a to decreasing one, while the
next rising $dR/dm_{T}$ part is connected with more gentle radii change at
growing $m_{T}$. The $long$ direction, by contrast, shows growing behaviour
ending up with a plateau without any clear flat beginning.

To interpret above mentioned results, one can use a very recent iHKM femto-
analysis [51] as follows. When the system is just formed (at times near 0.1
fm/$c$), huge gradients of density in the transverse direction take place,
since the system is essentially finite and borders with vacuum. The gradient
is not equal locally along the radial directions: it is more strong at the
periphery, and less strong in the center, where soft hadrons mostly come from.
In the vicinity of the geometrical center of the system, its decay into free
particles happens at significantly later proper times than for the most of
other parts of the system. Such a difference in proper times of spectra
formation can be up to 5 fm/$c$ [51]. Thus, in the context of the present
study, one can expect that maximal formation times (about $15$ fm/$c$) should
be typical for the ultra-soft pions with transverse momenta less than 0.3
GeV/$c$, that are emitted from the central region of the fireball. Whereas for
the pions radiated from other parts of the decaying expanding system and
having higher transverse momenta, $0.45-2$ GeV/$c$, the (proper) time of
maximal emission is close to 10 fm/$c$ [51]. The situation looks like a pion
trap formed in the center of the created quark-gluon — hadron system. The
hadrons stay together (cannot leave the system) for a longer time in the
system’s center because of the following:

  * •

firstly, very low (close to zero) transverse collective velocity;

  * •

secondly, smaller density gradient in the center during almost the entire
duration of the evolution, as compared to non-central and periphery parts.

  * •

thirdly, the initially highest density in the geometrical center in central
and semi-central nucleus-nucleus collisions. Therefore, for this high-mass-
density region, it is difficult to expand because of the relatively small
transverse pressure gradient in the center.

Note, that we see the mentioned effects of relative flattening of the
transverse radii behaviour in the ultra-soft momentum region in “hydro plus
hadronic cascade” models, where the dynamics is built in quasi-classical
approximation. The Bose-Einstein quantum statistics effect is taken into
account in the very final stage only, but even this way, using the femtoscopy
method, one can see specific (delayed freeze-out) behaviour of emission in the
ultra-low-momentum region (see also Fig. 1 from [51]). In particular, as we
already mentioned, the femtoscopy analysis brings us the signal about a long-
lived ultra-soft pion trap formed in the central part of the system created in
ultra-relativistic heavy ion collision.

The dependencies of the radii on $m_{T}$, together with the power-law fits,
are presented for the three centrality classes in Fig. 5. The fitting was
performed outside the plateau region observed in Fig. 4 (above 0.25 GeV/$c$).
As one can see from this figure, the power-law function in the form
$am_{T}^{-b}$ works well, especially for all three $R_{i}$ in the 3D case in
the iHKM model and for the $long$ direction of the LQTH model. As for the
transverse radii from the LQTH, the first 2-3 corresponding points are mostly
below expectations. In both models, the behaviour most different from the
power-law in the low-$m_{T}$ region is observed for the 1D $R_{\mathrm{inv}}$
points, which lie, up to 0.25 GeV/$c$, below the extrapolated line.

Figure 6:  The correlation strength $\lambda$ parameters from 1D (black) and
3D (red) fits to the two-pion CF functions for the ten $m_{T}$ bins. The iHKM
(top) and the LQTH (bottom) results are shown for $c=10-20\%$ Pb-Pb collisions
at $\sqrt{s_{NN}}=5.02$ TeV.

In Fig. 5 the fall of 1D radii is less prominent than the fall for radii from
3D studies. This can be connected with the shape of the respective correlation
functions. The 1D functions are less Gaussian than the 3D functions, and
therefore the fitting has to compensate for this by reducing the correlation
strength parameter $\lambda$. The $m_{T}$ dependencies of lambda in the two
applied models are shown in Fig. 6. In the low $m_{T}$ region the lambda
behaves similarly in both cases, whereas with growing $m_{T}$ it tends to go
down in the 1D case.

Finally, we analyzed the femto-radii dependence on pair velocity, see Fig. 7
for different centrality intervals and for both used model approaches. In all
cases, we notice a slow radii falling at low velocities. All the distributions
were fitted with the analytic formula

$R(\beta)=a+\frac{b}{e^{(\beta-c)/d}+1},$ (5)

where $\beta=k_{T}/m_{T}$ is pair velocity, and $a$, $b$, $c$, $d$ are free
parameters of the fit. As particle velocity decreases, the radii seem to reach
a saturated value, similar in a very broad range of $\beta$. Therefore such
value could be used to determine an extrapolated value of the radii for ultra-
soft pions, instead of the historically used dependence based on power-law
$m_{T}$ scaling.

Figure 7:  The identical pion femtoscopy radii dependencies on pair velocity
$\beta$ in the iHKM (full markers) and the LQTH (empty markers) models. The
results for the six centrality classes are shown. The lines represent fits to
the radii by the formula (5).

##  V SUMMARY

In this paper, the model study of femtoscopic radii obtained from simulated
data of iHKM and LQTH were shown. The results point out that the power-law
character of the femtoscopic $m_{T}$-scaling is no longer valid for 1D studies
in the low $m_{T}$ region (below 250 MeV$/c$). In the case of 3D studies, the
radii are better described by the power-law dependence than 1D radii, however,
the final conclusions are model-dependent. The iHKM model follows a power-law
description in all 3 directions. On the other hand, for the LQTH model power-
law description holds only in the longitudinal direction. Both models show
that in low-$m_{T}$ region radii distribution tends to be more linear than
power-law or even flat in some cases. This results in smaller radii than
expected from the extrapolation of the power-law scaling from intermediate and
high $m_{\rm{T}}$.

This modification of the behaviour is motivated by the capacity of the ultra-
soft pions to decouple from the energetic and dense core of the expanding
system, resulting in longer emission times but from the centre of the created
system. This can be interpreted as both spatial and kinematic trapping of the
ultra-soft pions by the more dense part of the created system. Such pions
cannot escape this phase-space region until this confining mechanism is
released. This finding is fundamental for addressing experimentally the
expected pion radii in non-identical femtoscopic studies when the second
particle has a large mass difference with the pion mass, $e.g.$,
pion–deuteron, pion–omega, or pion–charm hadron pairs.

The observations of the non-monotonic behaviour of ultra-soft pions obtained
from semi-classical models should also be explored from the quantum
perspective. The role of confinement and its connection to the direct photon
puzzle [52, 53, 54, 55, 56] might as well influence the pion homogeneity
lengths and should be studied further.

##  VI ACKNOWLEDGMENTS

This work is funded by the Research University – Excellence Initiative of
Warsaw University of Technology via the strategic funds of the Priority
Research Centre of High Energy Physics and Experimental Techniques, the IDUB
POSTDOC programme, the IDUB YOUNG-PW programme, Scientific Council of the
discipline grant programme and by the Polish National Science Centre under
agreements no. 2022/45/B/ST2/02029, and no. 2023/49/N/ST2/03525.

## References

  * Goldhaber  _et al._ [1960] G. Goldhaber, S. Goldhaber, W. Lee, and A. Pais, Influence of Bose-Einstein Statistics on the Antiproton-Proton Annihilation Process, Phys. Rev. 120, 300 (1960).
  * Kopylov and Podgoretsky [1972] G. I. Kopylov and M. I. Podgoretsky, Correlations of identical particles emitted by highly excited nuclei, Sov. J. Nucl. Phys. 15, 219 (1972).
  * Kopylov [1974] G. I. Kopylov, Like particle correlations as a tool to study the multiple production mechanism, Phys. Lett. B 50, 472 (1974).
  * Hanbury Brown and Twiss [1954] R. Hanbury Brown and R. Q. Twiss, LXXIV. A new type of interferometer for use in radio astronomy, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45, 663 (1954), https://doi.org/10.1080/14786440708520475 .
  * Twiss and Hanbury Brown [1956] R. Q. Twiss and R. Hanbury Brown, Correlation between photons in two coherent beams of light, Nature 177, 27 (1956).
  * Twiss  _et al._ [1957] R. Q. Twiss, A. G. Little, and R. Hanbury Brown, Correlation between photons, in coherent beams of light, detected by a coincidence counting technique, Nature 180, 324 (1957).
  * Sinyukov and Shapoval [2013] Yu. M. Sinyukov and V. M. Shapoval, Correlation femtoscopy of small systems, Phys. Rev. D 87, 094024 (2013).
  * Akkelin and Sinyukov [1995] S. V. Akkelin and Yu. M. Sinyukov, The HBT-interferometry of expanding sources, Phys. Lett. B 356, 525 (1995).
  * Koonin [1977] S. E. Koonin, Proton pictures of high-energy nuclear collisions, Phys. Lett. B 70, 43 (1977).
  * Gyulassy  _et al._ [1979] M. Gyulassy, S. K. Kauffmann, and L. W. Wilson, Pion interferometry of nuclear collisions. I. Theory, Phys. Rev. C 20, 2267 (1979).
  * Zajc  _et al._ [1984] W. A. Zajc, J. A. Bistirlich, R. R. Bossingham, H. R. Bowman, C. W. Clawson, K. M. Crowe, K. A. Frankel, J. G. Ingersoll, J. M. Kurck, C. J. Martoff, _et al._ , Two-pion correlations in heavy ion collisions, Phys. Rev. C 29, 2173 (1984).
  * Makhlin and Sinyukov [1987] A. N. Makhlin and Yu. M. Sinyukov, The pion interferometry theory for the hydrodynamic stage of multiple processes, Sov. J. Nucl. Phys 46, 345 (1987).
  * Makhlin and Sinyukov [1988] A. N. Makhlin and Yu. M. Sinyukov, Hydrodynamics of hadron matter under pion interferometric microscope, Z. Phys. C 39, 69 (1988).
  * Gersch [1987] H. U. Gersch, Coulomb final state interaction in identical boson interferometry, Zeitschrift für Physik A Atomic Nuclei 327, 115 (1987).
  * Sinyukov [1994] Yu. M. Sinyukov, Spectra and correlations in locally equilibrium hadron and quark-gluon systems, Nucl. Phys. A 566, 589 (1994).
  * Aamodt  _et al._ [2011a] K. Aamodt, A. Abrahantes, D. Adamova, M. Aggarwal, G. Aglieri Rinella, A. Agocs, S. Salazar, Z. Ahammed, N. Ahmad, A. Masoodi, _et al._ (ALICE Collaboration), Two-pion Bose-Einstein correlations in central Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Lett. B 696, 328 (2011a).
  * Adam  _et al._ [2015a] J. Adam, D. Adamová, M. M. Aggarwal, G. Aglieri Rinella, M. Agnello, N. Agrawal, Z. Ahammed, S. U. Ahn, I. Aimo, S. Aiola, _et al._ (ALICE Collaboration), One-dimensional pion, kaon, and proton femtoscopy in Pb-Pb collisions at $\sqrt{{s}_{\mathit{NN}}}=2.76$ TeV, Phys. Rev. C 92, 054908 (2015a).
  * Aamodt  _et al._ [2010] K. Aamodt _et al._ (ALICE Collaboration), Two-pion Bose-Einstein correlations in $pp$ collisions at $\sqrt{s}=900$ GeV, Phys. Rev. D 82, 052001 (2010), arXiv:1007.0516 [hep-ex] .
  * Adam  _et al._ [2015b] J. Adam _et al._ (ALICE Collaboration), Two-pion femtoscopy in p-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, Phys. Rev. C 91, 034906 (2015b), arXiv:1502.00559 [nucl-ex] .
  * Aamodt  _et al._ [2011b] K. Aamodt _et al._ (ALICE Collaboration), Femtoscopy of $pp$ collisions at $\sqrt{s}=0.9$ and 7 TeV at the LHC with two-pion Bose-Einstein correlations, Phys. Rev. D 84, 112004 (2011b), arXiv:1101.3665 [hep-ex] .
  * Adam  _et al._ [2016a] J. Adam _et al._ (ALICE Collaboration), Centrality dependence of pion freeze-out radii in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. C 93, 024905 (2016a), arXiv:1507.06842 [nucl-ex] .
  * Acharya  _et al._ [2021] S. Acharya _et al._ (ALICE Collaboration), Pion-kaon femtoscopy and the lifetime of the hadronic phase in Pb$-$Pb collisions at $\sqrt{s_{\rm{NN}}}$ = 2.76 TeV, Phys. Lett. B 813, 136030 (2021), arXiv:2007.08315 [nucl-ex] .
  * Akkelin and Sinyukov [2004] S. V. Akkelin and Yu. M. Sinyukov, Phase-space densities and effects of resonance decays in a hydrodynamic approach to heavy ion collisions, Phys. Rev. C 70, 064901 (2004).
  * Shapoval  _et al._ [2013] V. M. Shapoval, P. Braun-Munzinger, Iu. A. Karpenko, and Yu. M. Sinyukov, Femtoscopic scales in p+p and p+Pb collisions in view of the uncertainty principle, Phys. Lett. B 725, 139 (2013).
  * Adler  _et al._ [2001] C. Adler _et al._ (STAR Collaboration), Pion interferometry of $\sqrt{s_{NN}}=130$ GeV Au+Au collisions at RHIC, Phys. Rev. Lett. 87, 082301 (2001), arXiv:nucl-ex/0107008 .
  * Adler  _et al._ [2004] S. S. Adler _et al._ (PHENIX Collaboration), Bose-Einstein Correlations of Charged Pion Pairs in $\mathrm{A}\mathrm{u}+\mathrm{A}\mathrm{u}$ Collisions at $\sqrt{{s}_{NN}}=200\text{ }\text{ }\mathrm{G}\mathrm{e}\mathrm{V}$, Phys. Rev. Lett. 93, 152302 (2004).
  * Adams  _et al._ [2005] J. Adams _et al._ (STAR Collaboration), Pion interferometry in Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV, Phys. Rev. C 71, 044906 (2005), arXiv:nucl-ex/0411036 .
  * Lisa  _et al._ [2005] M. A. Lisa, S. Pratt, R. Soltz, and U. Wiedemann, Femtoscopy in relativistic heavy ion collisions: Two decades of progress, Annual Review of Nuclear and Particle Science 55, 357 (2005), https://doi.org/10.1146/annurev.nucl.55.090704.151533 .
  * Abelev  _et al._ [2009] B. Abelev, M. Aggarwal, Z. Ahammed, B. Anderson, D. Arkhipkin, G. Averichev, J. Balewski, O. Barannikova, L. Barnby, J. Baudot, _et al._ (STAR Collaboration), Pion interferometry in Au+Au and Cu+Cu collisions at $\sqrt{s_{NN}}=62.4$ and 200 GeV, Phys. Rev. C 80, 024905 (2009).
  * Acharya  _et al._ [2020a] S. Acharya, D. Adamová, A. Adler, J. Adolfsson, M. M. Aggarwal, G. Aglieri Rinella, M. Agnello, N. Agrawal, Z. Ahammed, S. Ahmad, _et al._ (ALICE Collaboration), Search for a common baryon source in high-multiplicity pp collisions at the LHC, Phys. Lett. B 811, 135849 (2020a).
  * Naboka  _et al._ [2015] V. Yu. Naboka, S. V. Akkelin, Iu. A. Karpenko, and Yu. M. Sinyukov, Initialization of hydrodynamics in relativistic heavy ion collisions with an energy-momentum transport model, Phys. Rev. C 91, 014906 (2015).
  * Lednicky [2009] R. Lednicky, Finite-size effects on two-particle production in continuous and discrete spectrum, Phys. Part. Nucl. 40, 307 (2009), arXiv:nucl-th/0501065 .
  * Kisiel [2010] A. Kisiel, Non-identical particle femtoscopy at s(NN)**(1/2) = 200-AGeV in hydrodynamics with statistical hadronization, Phys. Rev. C 81, 064906 (2010), arXiv:0909.5349 [nucl-th] .
  * Abelev  _et al._ [2012] B. Abelev _et al._ (ALICE), Pion, Kaon, and Proton Production in Central Pb–Pb Collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. Lett. 109, 252301 (2012), arXiv:1208.1974 [hep-ex] .
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE), Centrality dependence of $\pi$, K, p production in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV, Phys. Rev. C 88, 044910 (2013), arXiv:1303.0737 [hep-ex] .
  * Begun  _et al._ [2014] V. Begun, W. Florkowski, and M. Rybczynski, Explanation of hadron transverse-momentum spectra in heavy-ion collisions at $\sqrt{s}_{NN}=$ 2.76 TeV within chemical non-equilibrium statistical hadronization model, Phys. Rev. C 90, 014906 (2014), arXiv:1312.1487 [nucl-th] .
  * Begun and Florkowski [2015] V. Begun and W. Florkowski, Bose-Einstein condensation of pions in heavy-ion collisions at the CERN Large Hadron Collider (LHC) energies, Phys. Rev. C 91, 054909 (2015), arXiv:1503.04040 [nucl-th] .
  * Naboka  _et al._ [2016] V. Yu. Naboka, Iu. A. Karpenko, and Yu. M. Sinyukov, Thermalization, evolution, and observables at energies available at the CERN Large Hadron Collider in an integrated hydrokinetic model of A+A collisions, Phys. Rev. C 93, 024902 (2016).
  * Shapoval and Sinyukov [2019] V. M. Shapoval and Yu. M. Sinyukov, Bulk observables in $\mathrm{Pb}+\mathrm{Pb}$ collisions at $\sqrt{{s}_{NN}}=5.02\phantom{\rule{1.60004pt}{0.0pt}}\mathrm{TeV}$ at the CERN Large Hadron Collider within the integrated hydrokinetic model, Phys. Rev. C 100, 044905 (2019).
  * Bożek and Wyskiel-Piekarska [2012] P. Bożek and I. Wyskiel-Piekarska, Particle spectra in Pb-Pb collisions at $\sqrt{{\mathbf{s}}_{\mathbf{NN}}}=\mathbf{2}.\mathbf{76}$ TeV, Phys. Rev. C 85, 064915 (2012).
  * Chojnacki  _et al._ [2012] M. Chojnacki, A. Kisiel, W. Florkowski, and W. Broniowski, THERMINATOR 2: THERMal heavy IoN generATOR 2, Comput. Phys. Commun. 183, 746 (2012).
  * Bozek  _et al._ [2019] P. Bozek, W. Broniowski, M. Rybczynski, and G. Stefanek, GLISSANDO 3: GLauber Initial-State Simulation AND mOre…, ver. 3., Comput. Phys. Commun. 245, 106850 (2019).
  * Huovinen and Petersen [2012] P. Huovinen and H. Petersen, Particlization in hybrid models, Eur. Phys. J. A 48, 171 (2012).
  * Bass  _et al._ [1998] S. A. Bass _et al._ , Microscopic Models for Ultrarelativistic Heavy Ion Collisions, Prog. Part. Nucl. Phys. 41, 225–370 (1998).
  * Bleicher  _et al._ [1999] M. Bleicher _et al._ , Relativistic Hadron-Hadron Collisions in the Ultra-Relativistic Quantum Molecular Dynamics Model, J. Phys. G: Nucl. Part. Phys. 25, 1859–1896 (1999).
  * Kisiel  _et al._ [2014] A. Kisiel, M. Gałażyn, and P. Bożek, Pion, kaon, and proton femtoscopy in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV modeled in (3+1)D hydrodynamics, Phys. Rev. C 90, 064914 (2014).
  * Kisiel [2018] A. Kisiel, Pion-kaon femtoscopy in Pb–Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV modeled in (3+1)D hydrodynamics coupled to Therminator 2 and the effect of delayed kaon emission, Phys. Rev. C 98, 044909 (2018).
  * Acharya  _et al._ [2020b] S. Acharya, D. Adamová, S. P. Adhya, A. Adler, J. Adolfsson, M. M. Aggarwal, G. Aglieri Rinella, M. Agnello, N. Agrawal, _et al._ (ALICE Collaboration), Production of charged pions, kaons, and (anti-)protons in pb-pb and inelastic $pp$ collisions at $\sqrt{{s}_{NN}}=5.02$ tev, Phys. Rev. C 101, 044907 (2020b).
  * Pratt [1986] S. Pratt, Pion interferometry of quark-gluon plasma, Phys. Rev. D 33, 1314 (1986).
  * Bertsch and Brown [1989] G. Bertsch and G. E. Brown, Temporal development of the plasma phase transition, Phys. Rev. C 40, 1830 (1989).
  * Sinyukov  _et al._ [2023] Yu. M. Sinyukov, V. M. Shapoval, and M. D. Adzhymambetov, Space-time structure of particle emission and femtoscopy scales in ultrarelativistic heavy-ion collisions, Universe 9, 433 (2023), arXiv:2310.16233 [nucl-th] .
  * Paquet  _et al._ [2016] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Production of photons in relativistic heavy-ion collisions, Phys. Rev. C 93, 044906 (2016).
  * Adare  _et al._ [2015] A. Adare, S. Afanasiev, C. Aidala, N. Ajitanand, Y. Akiba, R. Akimoto, H. Al-Bataineh, H. Al-Ta’ani, J. Alexander, A. Angerami, _et al._ (PHENIX Collaboration), Centrality dependence of low-momentum direct-photon production in Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV, Phys. Rev. C 91, 064904 (2015).
  * Adam  _et al._ [2016b] J. Adam, D. Adamová, M. M. Aggarwal, G. Aglieri Rinella, M. Agnello, N. Agrawal, Z. Ahammed, S. U. Ahn, S. Aiola, A. Akindinov, _et al._ (ALICE Collaboration), Direct photon production in Pb–Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Lett. B 754, 235 (2016b).
  * Naboka  _et al._ [2020] V. Yu. Naboka, Yu. M. Sinyukov, and G. M. Zinovjev, Photon spectra and anisotropic flow in heavy ion collisions at the top RHIC energy within the integrated hydrokinetic model with photon hadronization emission, Nucl. Phys. A 1000, 121843 (2020).
  * Fujii  _et al._ [2017] H. Fujii, K. Itakura, and C. Nonaka, Photon emission at hadronization, Nucl. Phys. A 967, 704 (2017).