Original Article

Shogo NishiyamalabelATomohiro KaralabelA Brian ThorsbrolabelB Hiromi
SaidalabelC Yohsuke TakamorilabelD Masaaki TakahashilabelE Takayuki
OhgamilabelF Kohei IchikawalabelG, labelH Rainer SchödellabelI

[labelA]Miyagi University of Education, Sendai, Miyagi, Japan
[labelB]Observatoire de la Côte d’Azur, CNRS UMR 7293, BP4229, Laboratoire
Lagrange, F-06304 Nice Cedex 4, France [labelC]Daido University, Naogya,
Aichi, Japan [labelD]National Institute of Technology, Gobo, Wakayama, Japan
[labelE]Aichi University of Education, Kariya, Aichi, Japan [labelF]National
Astronomical Observatory of Japan, Tokyo, Japan [labelG]Tohoku University,
Sendai, Miyagi, Japan [labelH]Waseda University, Tokyo, Japan
[labelI]Instituto de Astrofísica de Andalucía (IAA)-CSIC, Granada, Spain

Shogo Nishiyama, Miyagi University of Education, 149 Aramaki-Aza-Aoba, Aoba-
ku, Sendai, Miyagi 980-0845, Japan<EMAIL_ADDRESS>
# Origin of an Orbiting Star Around the Galactic Supermassive Black Hole

###### Abstract

The tremendous tidal force that is linked to the supermassive black hole
(SMBH) at the center of our galaxy is expected to strongly subdue star
formation in its vicinity. Stars within $1\mbox{${}^{\prime\prime}$}$ from the
SMBH thus likely formed further from the SMBH and migrated to their current
positions. In this study, spectroscopic observations of the star S0-6/S10, one
of the closest (projected distance from the SMBH of $\approx
0\mbox{$.\\!\\!^{\prime\prime}$}3$) late-type stars were conducted. Using
metal absorption lines in the spectra of S0-6, the radial velocity of S0-6
from 2014 to 2021 was measured, and a marginal acceleration was detected,
which indicated that S0-6 is close to the SMBH. The S0-6 spectra were employed
to determine its stellar parameters including temperature, chemical abundances
([M/H], [Fe/H], [$\alpha$/Fe], [Ca/Fe], [Mg/Fe], [Ti/Fe]), and age. As
suggested by the results of this study, S0-6 is very old ($\gtrsim 10$ Gyr)
and has an origin different from that of stars born in the central pc region.

###### keywords:

astrophysics, infrared astronomy, black hole, spectroscopy

##  1 Introduction

Stars within approximately 1′′( $\approx 0.04$ pc at a distance of 8 kpc from
the Earth) from the supermassive black hole (SMBH) at the center of our galaxy
are referred to as “S-stars”111There is no common definition for S-stars.
Here, S-stars are defined as those within $1\mbox{${}^{\prime\prime}$}$ from
Sgr A*, which include both early- and late-type stars.. S-stars are a
collection of both early-type and late-type stars, and the early-type are
B0–B3 main sequence stars with masses of $8$–$14\,M_{\odot}$ [2, 1]. The late-
type stars are G–M red giants with estimated initial masses of
$0.5$–$2\,M_{\odot}$ [3].

S-stars provide a unique test bed for probing the strong gravity around SMBHs.
Owing to the large mass ratio between S-stars and the SMBH, named Sagittarius
A* (Sgr A*; $\approx 4\times 10^{6}M_{\odot}$[4]), the stars can be considered
as test particles that move in a static potential. A particularly interesting
target is the star named S0-2 or S2 (Fig. 1), whose orbital period is $\approx
16$ years [4, 2]. We have been carrying out spectroscopic monitoring
observations since 2014 with the aim to precisely monitor the radial velocity
(RV) of S0-2/S2 around the pericenter passage in 2018 [5]. We were able to
spatially resolve S0-2/S2 and obtain high-resolution near-infrared (NIR)
spectra with the use of the infrared camera and spectrograph (IRCS) [6] on the
Subaru telescope [7], in combination with the Subaru adaptive optics (AO) [8]
and laser guide star (LGS) systems [9]. The observational results and the
obtained orbit of S0-2/S2 are shown in Fig. 2.

Figure 1: $K^{\prime}$-band image of the galactic center region taken with
Subaru/IRCS/AO188 in 2017. North is up, and east is to the left in the image.
The positions of S0-2/S2, Sgr A*, and our target S0-6/S10, are shown by yellow
circles.

The observations and analysis provide not only the orbital parameters for
S0-2/S2 but also the physical parameters for the SMBH. In our study, the mass
and distance of the SMBH were determined to be $(4.23\pm 0.07)\times
10^{6}\,M_{\odot}$ and $8.10\pm 0.07$ kpc, respectively (Saida et al.[10]; see
also [11, 12]). It should be noted that by observing stars that orbit the
SMBH, we can determine the mass of the SMBH with an uncertainty of only a few
percent, which is much smaller than those from observations of the shadow of
the SMBH [13].

Figure 2: Observational results and derived orbit of S0-2/S2. The blue
squares, red diamonds, and green circles respectively denote the results of
the VLT, Keck/Gemini, and Subaru observations. The darkness of the colors
indicates the epoch of observations, with the most recent ones represented by
the darkest points. The results of an orbital fitting are shown by blue and
red lines, which indicate the blue- and red-shifted parts of the orbit,
respectively. (Left) Orbital motion of S0-2 on the sky. The positions are
relative to Sgr A*. S0-2/S2 revolves around Sgr A* in a clockwise direction in
this plot. (top middle) From top to bottom, offset in right ascension (RA),
offset in declination (Dec), and distance from Sgr A*, respectively, as
function of time. (top right) Three-dimensional visualization of the orbit of
S0-2/S2. The direction to the Earth is indicated by the white arrow, and the
position of Sgr A* is depicted by the black cross. (bottom right) Radial
velocity of S0-2 as function of time.

Because of the tremendous tidal force associated with the SMBH, it is highly
unlikely that the S-stars formed at their present position [14]. Therefore,
the stars are generally assumed to form at larger distances from Sgr A* and
were brought to their current location via dynamical processes. Considering
that the late-type S-stars are likely to be old, $\sim 3$–$10$ Gyr [3], they
are expected to be dynamically relaxed. In other words, identifying where they
formed using only kinematic information is impossible. In this case, stellar
properties including age and chemical abundances are crucial in order to
understand their origin.

The SMBH and S-stars are surrounded by a dense stellar system, the nuclear
star cluster (NSC). This is a massive ($\sim 2.5\times 10^{7}M_{\odot}$)
cluster with a half-light radius of $\sim 5$ pc [15, 16]. Recent observations
revealed that, besides a major stellar population with a mean metallicity at
least as high as twice the solar, there exists a subpopulation of stars that
are characterized by sub-solar metallicity [17, 18]. The subsolar population
may be a relic of an old star cluster migrated to the center of our galaxy
[19, 20]. In this paper, “NSC stars” is defined as stars within $\sim 5$ pc
from Sgr A* except S-stars. Habibi et al.[3] obtained spectra of 21 late-type
S-stars within $0.5\,$pc from Sgr A*, and estimated their $T_{\mathrm{eff}}$,
spectral type, age, and initial mass. However, owing to the medium resolution
($\sim 1500$) of the spectra, they could not determine the chemical abundance
of the late-type S-stars. Without metallicity information, stellar ages are
impossible to determine precisely. In this study, we obtain high-resolution
spectra of one of the late-type S-star S0-6 and determine its chemical
abundances and age to investigate its origin. S0-6 is located at $\approx
0\mbox{$.\\!\\!^{\prime\prime}$}3$ from Sgr A* (Fig. 1), which corresponds to
a transverse distance of $\approx 0.012\,$pc or $2400\,$AU. S0-6 is of great
interest because it appears to be one of the rare red giants that are very
close to Sgr A*.

##  2 Observations and data reduction

In this study, NIR spectroscopic observations were conducted using the Subaru
telescope and IRCS in the echelle mode (spectral resolution of $\sim 20,000$).
Spectra were obtained in the $K+$ setting each year from 2014 to 2019 and in a
special setting in the $K$-band ($K_{\mathrm{sp}}$) in 2021. The AO system on
Subaru AO188 was employed in the observations. During the observations in
2014, 2016, and 2017, the LGS system was applied.

Table 1 shows the echelle order numbers and wavelength coverage for the
settings. Although seven echelle orders can be simultaneously observed, we
utilized only three orders, namely, 25, 26 and 27, in the following analysis.
This is mainly due to strong atmospheric absorption lines in other orders.

Table 1: Wavelength coverage for echelle $K+$ and $K_{\mathrm{sp}}$ mode.  Setting | Order | Wavelength coverage
---|---|---
|  | [$\mathrm{\AA}$]
$K+$ | 25 |  $22500$–$23039$
| 26 |  $21634$–$22154$
| 27 |  $20833$–$21335$
$K_{\mathrm{sp}}$ | 25 |  $22324$–$22878$
| 26 |  $21466$–$21998$
| 27 |  $21670$–$21185$

The details of the data reduction are described in Nishiyama et al.[5]. The
first process of data reduction involves flat fielding, bad pixel correction,
and cosmic ray removal. Flat field images were acquired using observations of
a continuum source.

Prior to the extraction of spectra of S0-6, wavelength-calibrated stripe
images, where the $x$\- and $y$-axes of the stripe images are wavelength and
spatial position along the slit, respectively, were constructed. Here, all OH
emission lines along the slit direction are used for wavelength calibration of
the stripe images. The numbers of the OH lines used for the calibration are 6,
7, and 14 for the echelle orders 25, 26, and 27, respectively.

We extracted S0-6 spectra and then performed telluric correction of the S0-6
spectra using spectra of standard stars, which were bright early-A-type
dwarfs. Then, the S0-6 spectra were divided by the standard star spectra. We
conducted continuum fitting for the telluric-corrected spectrum for each
observation run. Finally, we combined all spectra to come up with a combined
spectrum for each observed epoch. To determine the metallicity, we combined
all the RV-corrected spectra obtained in our observations.

We estimated the signal-to-noise (S/N) ratio for S0-6 employing the procedure
described in Fukue et al.[21]. Here, the S/N ratio is defined as the inverse
of the noise for the normalized spectrum, $\sigma$, that is, S/N ratio
$=\sigma^{-1}$. Table 2 shows The S/N ratio for each observing run and each
order. The S/N ratios vary from 9 to 36 and depend not only on the number of
the used data sets but also on the seeing condition of the observations and
used guide star system, LGS or NGS. The S/N ratios calculated for the combined
spectra, composed of all the datasets from 2014 to 2021, were 57.3, 52.5, and
41.8, for echelle orders 27, 26, and 25, respectively.

Table 2: Signal-to-noise ratio of spectra with $R=20,000$.

Date |  | S/N ratio |
---|---|---|---
| order 27 | order 26 | order 25
2014.379 | 23.1 | 22.9 | 21.1
2015.635 | 22.5 | 23.5 | 20.4
2016.381 | 18.9 | 18.2 | 16.4
2017.341 | 22.1 | 20.7 | 19.6
2017.344 | 33.1 | 36.2 | 30.6
2017.347 | 30.1 | 30.2 | 28.9
2017.603 | 11.2 | 13.3 | 9.3
2017.609 | 23.1 | 23.1 | 14.1
2018.087 | 22.1 | 24.5 | 19.6
2021.420 | 28.1 | 12.3 | 22.7
Combined(a) | 57.3 | 52.5 | 41.8

(a) Combined spectra from 2014 to 2021.

##  3 Radial velocities and location of S0-6

Absorption lines in its spectra were used in order to determine the RV for
S0-6. First, we roughly estimated the wavelength shift of the observed spectra
using strong distinctive absorption lines in order 26. Second, the observed
spectra were compared with a wavelength-shifted model spectrum of a late-type
giant with a temperature of $\sim 4,000$ K [3]. If we found an absorption line
that was also identified in the model spectrum, we fitted the observed line
with a Gaussian function and calculated the RV of S0-6 from the wavelength
difference between the observed and model spectra. Finally, the amount of RV
correction was calculated using the IRAF rvcorrect task in order to obtain RVs
in the Local Standard of Rest reference frame.

The time variation of the RV for S0-6 is shown in Fig. 3. The error bars
include both statistical and systematic uncertainties (Table 3). The
statistical uncertainties are taken from the standard errors of the RV
measurements using the absorption lines, where 21 to 45 lines are employed.
The average statistical uncertainty is 1.25 km s-1.

Table 3: Radial velocities and uncertainties for S0-6.

Time | Number of lines | RV | Statistical uncertainty | Total uncertainty(a)
---|---|---|---|---
(year) |  | (km/s) | (km/s) | (km/s)
2014.379 | 31 | 94.96 | 1.27 | 1.46
2015.635 | 31 | 98.18 | 0.90 | 1.15
2016.381 | 26 | 96.07 | 1.52 | 1.68
2017.341 | 31 | 98.72 | 1.16 | 1.37
2017.344 | 39 | 97.72 | 0.81 | 1.09
2017.347 | 44 | 97.21 | 0.85 | 1.12
2017.603 | 21 | 98.90 | 2.00 | 2.13
2017.609 | 28 | 98.01 | 1.53 | 1.69
2018.087 | 45 | 96.93 | 1.04 | 1.27
2021.420 | 34 | 99.60 | 1.37 | 1.55

(a)The systematic uncertainty of $0.72\,$km s-1 is quadratically summed to the
statistical uncertainty.

To understand the long-term stability of the RV measurements, we applied
atmospheric OH emission lines. We extracted spectra at a position close to but
different from S0-6 without subtracting background emission and fit the OH
emission lines to measure wavelength-calibrated OH line wavelengths. By
comparing them with the OH emission wavelengths measured in laboratories, we
obtain a long-term stability of 0.72 km s-1 for the 7-year monitoring
observations. We include this as a systematic uncertainty (Table 3). When the
statistical and systematic uncertainties are summed quadratically, the average
of the total uncertainties in our RV measurements is 1.45 km s-1. The small
uncertainties shown above depict that this study is one of the most accurate
and precise RV measurements that were conducted for stars in the central pc
region of our galaxy.

Figure 3: (Top) Radial velocity of S0-6 as function of time. The blue line
represents the best-fit linear function. (Bottom) Residual from best-fit
linear function shown in the top panel.

Using a linear function, we fitted the RV plot, and the best-fit function is
shown in Fig. 3. The slope of the best-fit function is $0.48\pm 0.20$ km s-1
yr-1. This might suggest a marginal detection of the acceleration of S0-6
along the line of sight. The bottom panel in Fig. 3 represents the residuals
of the data points from the best-fit linear function. The standard deviation
of the data points around the best-fit function is 0.92 km/s, suggesting that
the observed data points are well-fitted with a linear function.

The RV of S0-6 has been also measured for more than 13 years using the Keck
telescope[22]. Acceleration along the line of sight was also detected as
$0.83\pm 0.12$ km s-1 yr-1, and the difference between the results is only
$1.5\sigma$, considering the uncertainties. The positive value of the
acceleration of S0-6 suggests that it is located in front of Sgr A* with its
acceleration vector pointing away from us.

The marginal detection of acceleration allows us to estimate the three-
dimensional distance $r$ of S0-6 from Sgr A*. The equation of motion for S0-6
whose mass is $m_{*}$ is

$m_{*}\boldsymbol{a}=-G\frac{Mm_{*}}{r^{2}}\frac{\boldsymbol{r}}{r},$ (1)

where $\boldsymbol{a}$ is the acceleration of S0-6, $M$ is the mass of Sgr A*,
$\boldsymbol{r}$ is the position vector, and $G$ is the gravitational
constant. We measured the acceleration along the line of sight, $a_{z}$, where
$\boldsymbol{a}=(a_{x},a_{y},a_{z})$ and $|a_{z}|$ is smaller than the three
dimensional acceleration $|\boldsymbol{a}|$. Hence, we obtain the following:

$|a_{z}|<|\boldsymbol{a}|=G\frac{M}{r^{2}}.$ (2)

This enables us to calculate the upper limit of the distance $r$ using the
following equation:

$r<\sqrt{\frac{GM}{|a_{z}|}}.$ (3)

Substituting $G=6.67\times 10^{-11}$ m3kg-1s-2, $M=4.23\times 10^{6}M_{\odot}$
[10], and $|a_{z}|=0.48$ km s-1 yr-1 into equation [3], we acquire $r<0.2$ pc
$\approx 4\times 10^{4}$ AU. This suggests that S0-6 is truly close to Sgr A*.

Note that on the basis of our current data set, the zero-acceleration
hypothesis cannot be ruled out. The significance of the detection of
$\boldsymbol{a}$ is still $2.4\,\sigma$. When fitting the observed RVs using a
function without acceleration (i.e., slope $=0$), we obtain
$\chi^{2}/\mathrm{d.o.f}=7.53/9$. This results in a p-value of 0.58, and we
cannot rule out the zero-acceleration hypothesis for S0-6. Although the
acceleration was detected using the Keck observations ($6.9\,\sigma$)[22],
more observations are needed to conclusively confirm its detection using our
own data.

##  4 Stellar parameters of S0-6

###  4.1 Data sets, StarKit, and calibration

We constrained the stellar parameters of S0-6 by computing synthetic spectra
and comparing them with the observed spectra. To compare the observed and
synthetic spectra, we utilized the StarKit code developed by Kerzendorf & Do
(2015)[23]. StarKit is a spectral fitting framework using Bayesian inference
to determine the best-fit parameters. StarKit simultaneously fits stellar
parameters such as the following: the effective temperature
$T_{\mathrm{eff}}$, the surface gravity $\log g$, the overall metallicity of
all elements [M/H], the $\alpha$ element abundance [$\alpha$/Fe], and the
spectral continuum.

A synthetic spectral library is necessary for StarKit to determine the set of
stellar parameters. Bentley et al.[24] found that the BOSZ grid[25] provides
smaller offsets and scatter than other ones, and we thus decided to utilize
the BOSZ grid. The parameter spaces covered by the BOSZ grid in our analysis
are as follows: $3500\leq T_{\mathrm{eff}}[\mathrm{K}]\leq 35,000$,
$0.0\leq\log g\leq+4.5$, $-2.0\leq\mathrm{[M/H]}\leq+0.75$, and
$-0.25\leq[\alpha/\mathrm{Fe}]\leq+0.5$.

Upon fitting $K$-band stellar spectra with StarKit, we excluded the following
absorption lines and wavelength ranges from the analysis. We did not use
echelle order 25 for the analysis as order 25 includes CO absorption features
at $\gtrsim 20900$ Å, and fitting the CO lines could introduce significant
biases in the determinations of the parameters. The absorption lines due to
Na, Ca, Sc, and V were known to be strong when compared with those for stars
in the galactic disk or the solar neighborhood [28, 26, 27]. Thus, the lines
associated with these four elements were excluded from the fit of the spectra.

Table 4: Data for calibration sample

Name | Spectral type | $T_{\rm{eff}}$ | $\log g$ | [Fe/H] |  [$\alpha$/Fe]
---|---|---|---|---|---
|  | (K) | (dex) | (dex) | (dex)
BD$-01$ 2971 | M5 | $3587\pm 19$ | $0.5^{(\mathrm{a})}$ | $-0.63\pm 0.26$ | $0.023\pm 0.222^{(\mathrm{b})}$
2MASS J19213390+3750202 | K6-7 | $3734\pm 97$ | $1.12\pm 0.10$ | $+0.30\pm 0.08$ | $0.068\pm 0.012^{(\mathrm{c})}$
HD 787 | K4 | $3962\pm 72$ | $1.12\pm 0.31$ | $-0.04\pm 0.13$ | $0.133\pm 0.114^{(\mathrm{b})}$
$\alpha$ Boo (Arcturus, HD124897) | K1.5 | $4296\pm 110$ | $1.66\pm 0.29$ | $-0.53\pm 0.11$ | $0.228\pm 0.015^{(\mathrm{c})}$
$\mu$ Leo (HD 85503) | K2 | $4494\pm 93$ | $2.44\pm 0.24$ | $+0.27\pm 0.15$ | $0.019\pm 0.009^{(\mathrm{c})}$
$\alpha$ Ari (HD 12929) | K2 | $4587\pm 65$ | $2.59\pm 0.18$ | $-0.11\pm 0.06$ | $0.06^{(\mathrm{b,d})}$
$\delta$ Dra (HD 180711) | G9 | $4856\pm 51$ | $2.69\pm 0.27$ | $-0.12\pm 0.10$ | $0.02\pm 0.02^{(\mathrm{e})}$
$\varepsilon$ Leo (HD 84441) | G1 | $5365\pm 76$ | $2.08\pm 0.29$ | $-0.05\pm 0.20$ | –

(a)Only two measurements with the same results are found. (b)[$\alpha$/Fe] $=$
([Mg/Fe]$+$[Ca/Fe]$+$[Ti/Fe])/3. (c)[$\alpha$/Fe] $=$
[$\alpha$/M]$+$[M/H]$-$[Fe/H], where [$\alpha$/M] and [M/H] are from APOGEE.
(d)[X/Fe] $=$ [X/H]$-$[Fe/H] $=$
$A$(X)$-\log(N_{\mathrm{X}}/N_{\mathrm{H}})_{\odot}-$[Fe/H], where $A$(X) is
the abundance of element X and $N_{\mathrm{X}}$ and $N_{\mathrm{H}}$ are the
number of atoms per unit volume for the element X and hydrogen, respectively.
Here, $\log(N_{\mathrm{X}}/N_{\mathrm{H}})_{\odot}$ values for Mg, Ca, and Ti
are taken from Asplund et al.(2009)[30]. (e)Taken from Tautvaišienė et
al.(2020)[31].

The stellar parameters determined by StarKit appear to depict a systematic
offset in comparison to other studies [18, 24]. We analyzed the stellar
parameters for a calibration sample (Table 4) to evaluate the systematic
offsets and the scatter when we use StarKit and IRCS $K$-band spectra. The
equivalent widths of Na, Ca, and CO absorption features in the $K$-band are
known to be sensitive to $\log g$[29]. However, the CO absorption features are
only partly included in our observed spectrum. Thus, we fixed $\log g$ for the
calibration sample to values determined in past studies (Table 4).

Table 5 shows the acquired parameter offsets between the StarKit and
literature values. Here, we assume that the stellar [M/H] values are the same
as the [Fe/H] values. These results suggest that the StarKit and IRCS spectra
tend to underestimate $T_{\mathrm{eff}}$, [M/H], and [$\alpha$/Fe] compared to
the literature values. These results are consistent with the studies by
Bentley et al.[24], where they concluded that StarKit underestimates the three
parameters.

Table 5:  Offsets between StarKit (ST) and literature (lit) values determined
for different parameters of the calibration sample

Order | $\Delta T_{\mathrm{eff}}^{(a)}$ | $\Delta\mathrm{[M/H]}^{(b)}$ |  $\Delta$[$\alpha$/Fe](c)
---|---|---|---
| (K) | (dex) | (dex)
27 | $-199\pm 196$ | $-0.30\pm 0.15$ | $-0.13\pm 0.27$
26 | $-204\pm 142$ | $-0.17\pm 0.06$ | $-0.18\pm 0.06$

(a) $\Delta T_{\mathrm{eff}}=T_{\mathrm{ST}}-T_{\mathrm{lit}}$ (b)
$\Delta\mathrm{[M/H]}=\mathrm{[M/H]_{ST}}-{\mathrm{[Fe/H]_{lit}}}$ (c)
$\Delta[\alpha/\mathrm{Fe}]=[\alpha/\mathrm{Fe}]_{\mathrm{ST}}-[\alpha/\mathrm{Fe}]_{\mathrm{lit}}$

###  4.2 Analysis of S0-6 with StarKit

We conducted spectral fitting of S0-6 using StarKit. The fit was performed for
the spectra of order 27 and 26 separately. Since the wavelength range of our
observed spectra is insensitive to $\log g$, we estimated $\log g$ by using
Yonsei-Yale isochrones [32], in which we find correlations between
$T_{\mathrm{eff}}$, the metallicity, and $\log g$. The late-type stars in the
central 0$.\\!\\!^{\prime\prime}$5 region are likely to be old, 3–10 Gyr [3];
thus we used 5- and 10-Gyr isochrones, the difference of which is negligible.

S0-6 is a late-type giant, and its $T_{\mathrm{eff}}$ was observed to be $\sim
4,000$ K [3]. Giant stars with $T_{\mathrm{eff}}\sim 4000$ are expected to
have a $\log g$ value in the range $\approx 0.8$ to $\approx 2.0$, if they are
not very metal poor ($[\mathrm{Fe/H}]>-1$). [19] For the first attempt to
identify the stellar parameters, we fit the spectra of S0-6 using StarKit, by
constraining $\log g$ to $0.8\leq\log g\leq 2.0$. We then obtained the first
results of the fit for $T_{\mathrm{eff}}$, [M/H], and [$\alpha$/Fe]. On the
basis of the first results and the systematic offset and uncertainties
obtained in §4.1, we applied a stronger constraint on $\log g$ using the
Yonsei-Yale isochrones. With the stronger second constraint on $\log g$, we
fit the spectra with StarKit and achieved a second set of results for
$T_{\mathrm{eff}}$, [M/H], and [$\alpha$/Fe]. By repeating this procedure
several times, we obtained stronger constraints on $\log g$ and more reliable
stellar parameters. Finally, we fit the spectra with a fixed $\log g$. The
fixed values are $\log g=0.9$ for order 27 and $1.2$ for order 26.

Figure 4: Observed spectra (black) and best-fitting model (red) of S0-6, and
residuals between observed and model spectra for order 27 (top) and 26
(bottom). The model parameters for order 27 are $T_{\mathrm{eff}}=3630$ K,
$\log g=0.9$, [M/H]$=-0.69$, and [$\alpha$/Fe] $=-0.23$, and those for order
26 are $T_{\mathrm{eff}}=3790$ K, $\log g=1.2$, [M/H]$=-0.72$, and
[$\alpha$/Fe] $=-0.22$. The BOSZ synthetic spectral library [25] was applied
in order to generate the model spectra. The vertical dashed lines represent
known spectral lines. They are labeled on top of the panels. The grey regions
of the spectra are excluded from the fitting process (see § 4.1).

Table 6 lists the resulting stellar parameters. Fig. 4 presents the spectra of
S0-6 and the best-fitting models for orders 27 and 26. The residuals between
the observed and model spectra are also represented. The standard deviations
of the residuals for orders 27 and 26 are 0.027 and 0.029, respectively.
Considering the systematic offset shown in Table 5, the calibrated parameters
are $T_{\mathrm{eff}}=3749$ K, $\mathrm{[M/H]}=-0.385$,
[$\alpha/\mathrm{Fe}]=-0.097$ (order 27), and $T_{\mathrm{eff}}=3998$ K,
$\mathrm{[M/H]}=-0.546$, [$\alpha/\mathrm{Fe}]=-0.041$ (order 26). Here, we
assume that the systematic offsets are the same for good S/N spectra (the
calibration samples) and for lower-quality spectra (S0-6).

Table 6: Results of spectral fitting  Order |  StarKit results |  | After calibration
---|---|---|---
| $T_{\mathrm{eff}}$ | $\log g$ | [M/H] | [$\alpha$/Fe] |  | $T_{\mathrm{eff}}$ | $\log g$ | [M/H] | [$\alpha$/Fe]
| (K) | fixed | (dex) | (dex) |  | (K) | fixed | (dex) | (dex)
27 | $3630$ | $0.90$ | $-0.685$ | $-0.227$ |  | $3749$ | $0.90$ | $-0.385$ | $-0.097$
26 | $3794$ | $1.2$ | $-0.716$ | $-0.221$ |  | $3998$ | $1.2$ | $-0.546$ | $-0.041$
Average |  |  |  |  |  | $3874$ |  | $-0.466$ | $-0.069$
Table 7: Error budget in stellar parameter determination

| $T_{\mathrm{eff}}$ |  | [M/H] |  | [$\alpha$/Fe]
---|---|---|---|---|---
Order | 27 | 26 |  | 27 | 26 |  | 27 | 26
Internal(a) | $56$ | $163$ |  | $0.20$ | $0.13$ |  | $0$ | $0.019$
Calibration(b) | $196$ | $142$ |  | $0.15$ | $0.06$ |  | $0.27$ | $0.06$
Total in each order(c) | $204$ | $216$ |  | $0.25$ | $0.14$ |  | $0.27$ | $0.063$
Average(d) | $210$ |  | $0.20$ |  | $0.196$
Standard error(e) | $124$ |  | $0.081$ |  | $0.028$
Total uncertainty(f) | $244$ |  | $0.22$ |  | $0.20$
Final result | $3870\pm 240$ |  | $-0.47\pm 0.22$ |  | $-0.07\pm 0.20$

(a)Standard deviation of results for two spectra combined from two subsets.
(b)Uncertainty in calibration, taken from Table 5. (c)Quadratic sum of
internal and calibration uncertainties. (d)Error propagation in averaging
results for orders 27 and 26. (e)Standard error of the results for orders 27
and 26, where the standard deviation of [M/H] for the two orders is divided by
$\sqrt{2}$. (f)Quadratic sum of uncertainties in averaging and standard error.

We summarize the uncertainties of our measurements in Table 7. We divided the
observed data into two subsets to measure the internal uncertainty for each
order. The spectra in the subsets were combined, which resulted in two
different spectra for each order. We performed fitting of the spectra with
StarKit and determined the stellar parameters. The “internal” uncertainties in
Table 7 are the standard errors of the results of the two subsets. To estimate
the total uncertainty in each order, the internal uncertainty and that derived
in calibration were quadratically summed (“Total in each order” in Table 7).
The “Average” uncertainties were determined by error propagation in averaging
the results for the two orders. Besides the uncertainties described above, we
include standard errors derived from the results of the two orders as an
uncertainty. Finally, the uncertainties in “average” and “standard error” are
quadratically summed in order to derive the total uncertainties. The derived
stellar parameters of S0-6 are as follows: $T_{\mathrm{eff}}=3870\pm 240$ K,
[M/H] $=-0.47\pm 0.22$, and [$\alpha$/Fe] $=-0.07\pm 0.20$. For
$T_{\mathrm{eff}}$, we find very good agreement between our result and the one
of Habibi et al.[3], $\approx 4000$–$4100$ K.

###  4.3 Analysis of S0-6 with Scandium lines and SME

We know from earlier studies that the scandium energy level transition with
the electron configuration 3d24s$-$3d4s4p, found near $2.2\,\mu$m, is
sensitive to $T_{\textrm{eff}}$ of stars[27, 33]. This is because, toward 4000
K, neutral scandium is becoming a minority species, which caused the scandium
spectral features to be very sensitive to temperature. This is further
strengthened by the fact that neutral scandium has a strong nuclear spin and
hence features hyperfine structures, and for this energy level transition with
a single s-electron at the lower level, the hyperfine structure is
particularly strongly enhanced. Moreover, this energy-level transition is
metastable, that is, has a long lifetime for spontaneous decay and therefore
has a high electron population count. The combination of a hyperfine structure
and high electron count makes the lines in the spectra particularly strong.
This benefits establishing an empirical relation between the equivalent width
of the scandium lines and $T_{\textrm{eff}}$ [33].

The empirical relation developed by Thorsbro et al.[33] finds the equivalent
width using the IRAF sbands task with a width of 3 Å around the scandium
lines. The linear regression between the equivalent width and the temperatures
results in a relationship with a standard deviation on the order of 50 K.

Based on the scandium lines 21730 Å, 21812 Å, and 21842 Å, we employ the
empirical relation between the equivalent width of the scandium lines and
$T_{\textrm{eff}}$. Assuming that [Sc/Fe] is not likely to vary beyond 0.15
dex [34] and that the variation of 0.1 dex in [Sc/Fe] leads to a variation of
$\sim$100 K [27], the uncertainty is $\sim$150 K, which includes the
statistical standard deviation from the empirical relation. Hence, using the
empirical relation, we found $T_{\textrm{eff}}=3830\pm 150$ K. This agrees
well with the results obtained using StarKit.

We also analyzed the star using SME [35, 36] to determine chemical abundances.
SME interpolates in a grid of one-dimensional MARCS atmosphere models [37].
These are hydrostatic model atmospheres in spherical geometry, computed
assuming LTE, chemical equilibrium, homogeneity, and conservation of the total
flux. The SME code has the advantage that it includes a flexible $\chi^{2}$
minimization tool to find the solution that best fits an observed spectrum in
a pre-specified spectral window.

SME also has a powerful continuum normalization code, normalizing the spectra
for every attempt to fit. In this way, SME can consider possible suppressed
continuum levels that could exist due to the crowded spectra of cool stars. We
found that for orders 26 and 27, we were able to do a good continuum
normalization, but less so for order 25, thus confirming our previous choice
not to use this order for the analysis. Order 25 has a CO-bandhead at the
later part of the order, which extends beyond the edge of the order, causing
it to have no anchor points for the continuum normalization on the higher
wavelength range of the order.

To determine the chemical abundances, we need to have an accurate list of the
atomic and molecular energy-level transitions. We use a list of atomic energy
level transitions from a previous work by Thorsbro et al. [38], where
wavelengths and line strengths (astrophysical $\log gf$-values) have been
updated using the solar center intensity atlas [39]. We also include CN
molecular lines, since the $K$-band is known for crowded CN-lines [40]. Table
8 summarizes the lines used in this work.

Table 8: Lines used for abundance determinations. Wavelength in air (Å) |  $E_{\mathrm{exc}}$ (eV) | $\log gf$
---|---|---
Fe I
20948.086 | $-0.883$ | 6.1190
20991.083 | $-3.019$ | 4.1427
21123.885 | $-0.857$ | 6.0664
21124.505 | $-1.647$ | 5.3342
21756.929 | $-0.715$ | 6.2182
21779.651 | $-4.298$ | 3.6398
Mg I
21059.757 | $-0.709$ | 6.7791
21061.095 | $-0.278$ | 6.7791
Ca I
20937.903 | $-1.357$ | 4.6807
20962.570 | $-0.740$ | 4.6814
20972.529 | $-1.002$ | 4.6807
Ti I
21767.610 | $-2.155$ | 2.5784
21782.997 | $-1.170$ | 1.7489
21897.437 | $-1.476$ | 1.7393
22004.510 | $-1.950$ | 1.7335

Table 9 summarizes the results of the analysis with SME. We found a
metallicity value of [Fe/H] = $-0.4\pm 0.15$. For the $\alpha$ element
abundances, we found [Mg/Fe] = $-0.30\pm 0.15$, [Ca/Fe] = $0.1\pm 0.15$, and
[Ti/Fe] = $-0.4\pm 0.15$, giving an average [$\alpha$/Fe] = $-0.2\pm 0.15$,
based on an average of Mg, Ca, and Ti. These also agree well with the results
from StarKit.

Table 9:  S0-6 parameters using the scandium method and SME

$T_{\mathrm{eff}}$ | [Fe/H] | [Mg/Fe] | [Ca/Fe] | [Ti/Fe] | [$\alpha$/Fe](a)
---|---|---|---|---|---
(K) | (dex) | (dex) | (dex) | (dex) | (dex)
$3830\pm 150$ | $-0.40\pm 0.15$ | $-0.30\pm 0.15$ | $+0.10\pm 0.15$ | $-0.40\pm 0.15$ | $-0.20\pm 0.15$

(a) [$\alpha$/Fe] = ([Ca/Fe] + [Mg/Fe] + [Ti/Fe])/3.

##  5 Discussion

###  5.1 Location and age of S0-6

We have marginally detected an acceleration in the RV for S0-6 (§ 3). This
suggests that S0-6 is accelerated by the SMBH [22]. When we use the
acceleration value and the mass of the SMBH, we obtain an upper limit for the
distance of S0-6 from the SMBH of $\sim 0.2$ pc $\approx 4\times 10^{4}$ AU.
Hence, our result suggests that S0-6 is located close to the SMBH. For the
late-type stars within $\sim 1\mbox{${}^{\prime\prime}$}$ from Sgr A*,
$T_{\mathrm{eff}}$, bolometric magnitude, and spectral types are determined
[3], but no measurement of metallicity has been conducted. If S0-6 is indeed
located at $0.2$ pc from the SMBH, to our knowledge, this study is the first
measurement of the chemical abundances of a late-type S-star.

No clear signal of binarity in observations of S0-6 was found. Fig. 3 shows
that the observed RVs are very well-fitted with a linear function, and no
clear periodic signal is found in our observations. S0-6 was also observed
with Keck/OSIRIS for 13 years [22]. They find no significant periodic signal
in the RVs of S0-6 and provide an upper limit of $0.1M_{\odot}$ on the mass of
a hypothetical companion object.

In order to estimate the age of S0-6, we plotted it on an HR diagram. The
bolometric magnitude $M_{\mathrm{bol}}$ can be calculated to be
$M_{\mathrm{bol}}=K-A_{K}-\mathrm{DM}+\mathrm{BC}_{K}$, where $K$, $A_{K}$,
and DM are the observed $K$-band magnitude, the amount of the interstellar
extinction in the $K$-band, and the distance modulus, respectively. Here, we
use $K=13.95\pm 0.04$ [41], $A_{K}=2.46\pm 0.03$ [42], and
$\mathrm{DM}=14.5\pm 0.18$ [10]. We also used the equation
$\mathrm{BC}_{K}=2.6-(T_{\mathrm{eff}}-3800)/1500$,[43] to determine the
bolometric correction $\mathrm{BC}_{K}$. We finally obtained
$M_{\mathrm{bol}}=-0.50\pm 0.25$ and $-0.47\pm 0.21$ for the StarKit and
scandium method results, respectively.

S0-6 is plotted on an HR diagram in Fig. 5. The result of Habibi et al.[3] is
also plotted, which agrees well with our results. When S0-6 is compared with
the theoretical $Z=0.3Z_{\odot}$ isochrones[44] (solid lines in Fig. 5), it is
located to the right of the 10 Gyr isochrone, which suggests that S0-6 is an
old star, as old as $\gtrsim 10$ Gyr.

Figure 5: HR diagram showing $T_{\mathrm{eff}}$ and $M_{\mathrm{bol}}$ of S0-6
obtained by our study and Habibi et al. [3] (black circle). The result of the
StarKit, and that of the scandium method and SME are represented by the red
circle and orange square, respectively. Overlaid are $Z=0.3Z_{\odot}$
theoretical isochrones [44] for 0.5, 1, 5, and 10 Gyr (solid lines, from left
to right). The dashed line is a theoretical 10 Gyr isochrone for
$Z=Z_{\odot}$.

To further confirm the reliability of our $T_{\mathrm{eff}}$ determination, we
estimated it using the line-depth-ratio (LDR) method. In the LDR method,
ratios of the depths of two stellar absorption lines are employed to determine
$T_{\mathrm{eff}}$. Using well-studied bright stars such as $\alpha$ Boo and
$\mu$ Leo, we found good LDR–$T_{\mathrm{eff}}$ relations for six line pairs
in echelle order 26 (Nishiyama et al, in preparation for publication). The
line pairs are shown in Table 10. The mean and standard deviation of
$T_{\mathrm{eff}}$ using the six LDR-$T_{\mathrm{eff}}$ relations is $4100\pm
90$ K, where we show only the statistical uncertainty. This is slightly higher
than $T_{\mathrm{eff}}$ derived by StarKit and the scandium method but
consistent with them within the uncertainties. Additionally,
$T_{\mathrm{eff}}$ derived by StarKit for echelle order 26 is 3998 K, which
agrees very well with that derived by the LDR method.

Table 10:  LDR pairs and derived $T_{\mathrm{eff}}$

Element | Wavelength | Element | Wavelength |  $T_{\mathrm{eff}}$(a)
---|---|---|---|---
| (Å) |  | (Å) | (K)
Ti I | 22010.51 | Na I | 22089.69 | 4105
Ti I | 22010.51 | Na I | 22062.42 | 4134
Sc I | 22058.00 | Na I | 22089.69 | 3977
Sc I | 22058.00 | Na I | 22062.42 | 4012
Ti I | 21903.35 | Na I | 22089.69 | 4166
Ti I | 21903.35 | Na I | 22062.42 | 4190

(a) Mean and standard deviation of six $T_{\mathrm{eff}}$ is $4100\pm 90$ K.

Similar to S0-6, NSC stars distributed to the right and below the oldest
isochrone were found [45, 46, 3]. Chen et al. [47] claimed that this
discrepancy can be explained by assuming younger stars with higher
metallicity. In this study, we found a similar discrepancy for an old and low-
metal star, which suggests a previously unknown systematic uncertainty in
observations and/or theoretical models.

###  5.2 Chemical abundance of S0-6 and stars in our galaxy

Recent observations [18, 19, 20] found that there is a metal-poor stellar
population ($\sim 7\,$%) in the NSC, whose mean metallicity is
$[\mathrm{M/H}]\sim-0.5$. Although Do et al. did not find stars with
$[\mathrm{M/H}]\sim-0.5$ [17], Rich et al. found stars with
$[\mathrm{Fe/H}]\sim-0.5$ [48]. These suggest that stars with similar [M/H]
values to S0-6 are not rare in the NSC, and S0-6 may have the same origin as
that of the NSC metal-poor population. An interesting aspect of S0-6 is its
closeness to Sgr A*, because strong tidal force associated with the SMBH may
provide an observable impact on the nature of S0-6.

To examine the origin of S0-6 in more detail, we analyze its position in an
[$\alpha$/Fe] vs. [M/H] diagram, which is useful for clarifying the chemical
evolution of stars and stellar systems. The top panel in Fig. 6 plots S0-6 and
stars in our galaxy. It is clear that most of the stars in the galaxy with
[M/H]$\lesssim 0.0$ show positive [$\alpha$/Fe]. We can observe two sequences
in the diagram for disk stars (green circles) and bulge stars (yellow
circles). They overlap, but the bulge stars tend to have larger [$\alpha$/Fe]
than the disk stars. At [Fe/H] $\sim-0.5$, the [$\alpha$/Fe] value for S0-6 is
slightly smaller than that for the bulge stars, whereas it is consistent with
the disk stars considering the observational uncertainty. There is a metal-
poor population in the bulge, and for some elements, these metal-poor stars
show different distributions in the diagrams from the other bulge stars and
disk stars [52]. Although the sequence is not clear in Fig. 6, S0-6 might be
on the bulge metal-poor population sequence. We currently do not have enough
samples to determine the distribution of late-type NSC stars in the diagram,
because of the difficulty in observing such NSC stars. A pioneering study to
determine the $\alpha$-element abundance of the late-type stars in the NSC
region, a few pc from Sgr A*, was recently conducted by Thorsbro et al.[33],
in which they determined [Si/Fe] for 15 stars (Fig. 6). All of the stars
measured by Thorsbro et al.[33] show larger [M/H] and [$\alpha$/Fe] than S0-6.

The [M/H] and [$\alpha$/Fe] values were measured for two late-type stars
within 1 pc from Sgr A* by Bentley et al.[24]. They employed a method very
similar to ours and found low-metal stars, whereas their [$\alpha$/Fe] values
are positive. One of the stars, NE1-1-003, shows a very similar chemical
abundances to S0-6, [M/H]$=-0.59\pm 0.11$ and [$\alpha$/Fe]$=0.05\pm 0.15$.
Although the abundances of S0-6 are likely to be different from a large part
of the stars in the NSC region, there might exist low [M/H] and low
[$\alpha$/Fe] abundance stars in the region, for example, S0-6 and NE1-1-003.

Figure 6:  (Top) [$\alpha$/Fe] vs. [M/H] plot for stars and stellar systems in
our galaxy. S0-6 is plotted as the red filled circle (StarKit result) and red
square (SME result). The abundances of globular clusters are from Pritzl et
al.[49] (cyan circles) assuming [Fe/H] $=$ [M/H]. The abundances of disk
stars, bulge stars, and bulge metal-poor stars are from Bensby et al.[50]
(green circles), Bensby et al.[51](yellow circles), and Lucey et al.[52] (gold
circles), respectively, assuming [Fe/H] $=$ [M/H] and [$\alpha$/Fe]
$=($[Mg/Fe]$+$[Si/Fe]$+$[Ca/Fe]$)/3$. The abundance of NSC stars are from
Thorsbro et al.[33] (magenta circles) assuming [Fe/H] $=$ [M/H] and
[$\alpha$/Fe] $=$ [Si/Fe], and from Bentley et al.[24] (magenta squares).
(Bottom) [$\alpha$/Fe] vs. [M/H] plot for stars in nearby dwarf galaxies. S0-6
is plotted by the red-filled circle. The abundance of the Sculptor dwarf
galaxy is from Hill et al.[53] (olive squares), and those for the Sgr dwarf
galaxy are from Monaco et al.[54] (light green circles), assuming [Fe/H] $=$
[M/H] and [$\alpha$/Fe] = ([Mg/Fe] + [Ca/Fe])/2. The abundance of LMC is from
Van der Swaelmen et al.[55] (dark green triangle) assuming [Fe/H] $=$ [M/H]
and [$\alpha$/Fe] = ([Mg/Fe] + [Si/Fe] + [Ca/Fe])/3.

###  5.3 Chemical abundance of S0-6 and stars in dwarf galaxies

Generally, stars in dwarf galaxies experience a different chemical evolution
from those in our galaxy, and thus they show a different distribution in the
[$\alpha$/Fe] vs. [M/H] diagram. In the bottom panel in Fig. 6, stellar
abundances for stars in three dwarf galaxies, the Sagittarius dwarf spheroidal
(Sgr) galaxy, the Sculptor dwarf galaxy, and the Large Magellanic Cloud (LMC)
are shown and compared with S0-6.

As can be seen, [M/H] of S0-6 is higher than that for the highest stellar
population in the Sculptor dwarf galaxy. On the other hand, the position of
S0-6 lies within the distributions of the stars of the LMC and the Sgr galaxy.
Even so, we do not conclude that S0-6 was born in the LMC or the Sgr galaxy.
The LMC is located at a distance of $\approx 50$ kpc and is orbiting around
our galaxy. A close encounter between the LMC and our galaxy likely occurred
$\sim 1.5$ Gyr ago, but the LMC has not been clearly tidally disrupted. The
Sgr galaxy is experiencing strong and disruptive tidal interactions with our
galaxy and has been almost totally disrupted, but has not yet been fully
merged. Hence it is difficult to imagine that stars born in the LMC or Sgr
galaxy would have migrated to the very center of our galaxy.

A more natural explanation for the origin might be that S0-6 was born in a
dwarf galaxy and then migrated to the center of our galaxy. Recently, a minor,
metal-poor population was found in the Galactic NSC with kinematics distinct
from the major, metal-rich population [18, 19], The metal-poor NSC population
shows kinematics distinct from the metal-rich population [18, 19], and its
mass fraction was estimated to be $\sim 7$–$18\,\%$ [19, 56]. It might be
natural to assume that S0-6 was born in the same cluster as the metal-poor
population, migrated to the center of our galaxy as a member of the cluster,
and reached the central pc. One of the scenarios to explain the low chemical
abundances of S0-6 is that S0-6 was born as a star in an NSC of a dwarf
galaxy, rather than a globular cluster of our galaxy.

##  6 Conclusions

In this paper, we report the results of our NIR spectroscopic monitoring
observations of the late-type S-star S0-6/S10 using the Subaru telescope and
IRCS. Our main results are as follows. We measured the RV of S0-6 with the
average uncertainty of $1.45\,$km s-1. This is one of the most precise RV
measurements of a star in the central $1\mbox{${}^{\prime\prime}$}$ region of
our galaxy. We marginally detected an acceleration in the RV of S0-6 of
$0.48\pm 0.20$ km s-1 yr-1 from $2014$ to $2021$, and obtained an upper limit
for the distance of S0-6 from the SMBH of $0.2\,\mathrm{pc}\approx 4\times
10^{4}\,$AU. We determined the stellar parameters for S0-6 using the StarKit
code and the observed spectra. The parameters are $T_{\mathrm{eff}}=3870\pm
240$ K, [M/H] $=-0.47\pm 0.22$, [$\alpha$/Fe] $=-0.07\pm 0.20$. We also
determined stellar parameters using the scandium method and SME, and obtained
$T_{\mathrm{eff}}=3830\pm 150$ K, [Fe/H] $=-0.40\pm 0.15$, which are
consistent with the StarKit results, and both results suggest the very old age
of S0-6, $\gtrsim 10$ Gyr. The abundance of other elements are: [Ca/Fe]
$=+0.10\pm 0.15$, [Mg/Fe]$=-0.30\pm 0.15$, and [Ti/Fe] $=-0.40\pm 0.15$. The
chemical abundances and age suggest that S0-6 experienced a different chemical
evolution from other stars in the center of our galaxy.

##  7 Acknowledgements

We wish to thank the Subaru Telescope staff for the support provided for our
observations. We thank Tuan Do, Rory Bentley, and Anja Feldmeier-Krause for
their kind assistance with the analysis using Starkit. This work was supported
by JSPS KAKENHI, Grant in Aid for Challenging Exploratory Research (Grant
Number 18K18760, 20H00178), Grant-in-Aid for Scientific Research(A) (Grant
Number 20H00178, 19H00695, 18H03720) and Grant-in-Aid for Scientific
Research(B) (50640843). This work was supported by the Tohoku Initiative for
Fostering GlobalResearchers for Interdisciplinary Sciences (TI-FRIS) of
MEXT’sStrategic Professional Development Program for Young Researchers. RS
acknowledges financial support from the State Agency for Research of the
Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the
Instituto de Astrofísica de Andalucía (SEV-2017-0709) and from grant
EUR2022-134031 funded by MCIN/AEI/10.13039/501100011033 and by the European
Union NextGenerationEU/PRTR. BT acknowledges the financial support from the
Wenner-Gren Foundation (WGF2022-0041).This research is based on data collected
at the Subaru Telescope,which is operated by the National Astronomical
Observatory of Japan.We are honored and grateful for the opportunity to
observe theUniverse from Maunakea, which has cultural, historical, and
naturalsignificance in Hawaii.

## References

  * [1] Habibi, M., Gillessen, S., Martins, F., Eisenhauer, F., Plewa, P. M., Pfuhl, O., et al. (2017) Twelve Years of Spectroscopic Monitoring in the Galactic Center: The Closest Look at S-stars near the Black Hole. ApJ. 847, id.120.
  * [2] Ghez, A. M., Duchêne, G., Matthews, K., Hornstein, S. D., Tanner, A., Larkin, J., et al.. (2003) The First Measurement of Spectral Lines in a Short-Period Star Bound to the Galaxy’s Central Black Hole: A Paradox of Youth. ApJ. 586, L127–L131.
  * [3] Habibi, M., Gillessen, S., Pfuhl, O., Eisenhauer, F., Plewa, P. M., von Fellenberg, S., et al. (2019) Spectroscopic Detection of a Cusp of Late-type Stars around the Central Black Hole in the Milky Way. ApJ. 872, id.L15
  * [4] Schödel, R., Ott, T., Genzel, R., Hofmann, R., Lehnert, M., Eckart, A., et al. (2002) A star in a 15.2-year orbit around the supermassive black hole at the centre of the Milky Way. Nature. 419, 694–696.
  * [5] Nishiyama, S., Saida, H., Takamori, Y., Takahashi, M., Schödel, R., Najarro, F., et al. (2018) Radial velocity measurements of an orbiting star around Sgr A*. PASJ. 70, id.74.
  * [6] Kobayashi, N., Tokunaga, A. T., Terada, H., Goto, M., Weber, M., Potter, R., et al. (2000) IRCS: infrared camera and spectrograph for the Subaru Telescope. SPIE Conference Series. 4008, 1056–1066.
  * [7] Iye, M., Karoji, H., Ando, H., Kaifu, N., Kodaira, K., Aoki, K., et al. (2004) Current Performance and On-Going Improvements of the 8.2 m Subaru Telescope. PASJ. 56, 381–397.
  * [8] Hayano, Y., Takami, H., Oya, S., Hattori, M., Saito, Y., Watanabe, M., et al. (2010) Commissioning status of Subaru laser guide star adaptive optics system. SPIE Conference Series. 7736, id. 77360N.
  * [9] Minowa, Y., Hayano, Y., Terada, H., Pyo, T-S., Oya, S., Hattori, M., et al. (2012) Subaru laser guide adaptive optics system: performance and science operation. 8447, id.84471F.
  * [10] Saida, H., Nishiyama, S., Ohgami, T., Takamori, Y., Takahashi, M., Minowa, Y., et al. (2019) A significant feature in the general relativistic time evolution of the redshift of photons coming from a star orbiting Sgr A*. PASJ. 71, id.126.
  * [11] GRAVITY Collaboration. (2019) A geometric distance measurement to the Galactic center black hole with 0.3% uncertainty. A&A. 625, id.L10.
  * [12] Do, T., Hees, A., Ghez, A., Martinez, G. D., Chu, D. S., Jia, S., et al. (2019b) Relativistic redshift of the star S0-2 orbiting the Galactic Center supermassive black hole. Science. 365, 664–668.
  * [13] Event Horizon Telescope Collaboration. (2022) First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way. ApJ. 930, id.L12.
  * [14] Mapelli, M. and Gualandris, A. (2016) Star Formation and Dynamics in the Galactic Centre. Lecture Notes in Physics, Berlin Springer Verlag. 905, 205.
  * [15] Schödel, R., Feldmeier, A., Neumayer, N., Meyer, L., and Yelda, S. (2014) The nuclear cluster of the Milky Way: our primary testbed for the interaction of a dense star cluster with a massive black hole. Classical and Quantum Gravity, 31, id.244007
  * [16] Neumayer, N., Seth, A. and Böker, T. (2020) Nuclear star clusters. A&A Rev.. 28, eid.4.
  * [17] Do, T., Kerzendorf, W., Winsor, N., Støstad, M., Morris, M. R., Lu, J. R., et al. (2015) Discovery of Low-metallicity Stars in the Central Parsec of the Milky Way. ApJ. 809, id.143.
  * [18] Feldmeier-Krause, A., Kerzendorf, W., Neumayer, N., Schödel, R., Nogueras-Lara, F., Do, T., et al. (2017) KMOS view of the Galactic Centre – II. Metallicity distribution of late-type stars. MNRAS. 464, 194–209.
  * [19] Feldmeier-Krause, A., Kerzendorf, W., Do, T., Nogueras-Lara, F., Neumayer, N., Walcher, C. J., et al. (2020) Asymmetric spatial distribution of subsolar metallicity stars in the Milky Way nuclear star cluster. MNRAS. 494, 396–410.
  * [20] Do, T., David Martinez, G., Kerzendorf, W., Feldmeier-Krause, A., Arca Sedda, M., Neumayer, N., et al. (2020) Revealing the Formation of the Milky Way Nuclear Star Cluster via Chemo-dynamical Modeling. ApJ. 901, id.L28.
  * [21] Fukue, K., Matsunaga, N., Yamamoto, R., Kondo, S., Kobayashi, N., Ikeda, Y., et al. (2015) Line-depth Ratios in H-band Spectra to Determine Effective Temperatures of G- and K-type Giants and Supergiants. ApJ. 812, id.64.
  * [22] Chu, D. S., Do, T., Ghez, A., Gautam, A. K.., Ciurlo, A., O’neil, K. K., et al. (2023) Evidence of a Decreased Binary Fraction for Massive Stars within 20 milliparsecs of the Supermassive Black Hole at the Galactic Center. ApJ. 948, id.94.
  * [23] Kerzendorf, W. and Do, T. (2015) starkit: First real release.
  * [24] Bentley, R. O., Do, T., Kerzendorf, W., Chu, D. S., Chen, Z., Konopacky, Q. et al. (2022) Measuring the $\alpha$-abundance of Subsolar-metallicity Stars in the Milky Way’s Central Half-parsec: Testing Globular Cluster and Dwarf Galaxy Infall Scenarios. ApJ. 925, id.77.
  * [25] Bohlin, R. C., Mészáros, S., Fleming, S. W., Gordon, Karl D., Koekemoer, Anton M. and Kovács, J. (2017) A New Stellar Atmosphere Grid and Comparisons with HST/STIS CALSPEC Flux Distributions. AJ. 153, id.234.
  * [26] Do, T., Kerzendorf, W., Konopacky, Q., Marcinik, J. M., Ghez, A., Lu, J. R. et al. (2018) Super-solar Metallicity Stars in the Galactic Center Nuclear Star Cluster: Unusual Sc, V, and Y Abundances. ApJ. 855, id.L5.
  * [27] Thorsbro, B., Ryde, N., Schultheis, M., Hartman, H., Rich, R. M., Lomaeva, M., et al. (2018) Evidence against Anomalous Compositions for Giants in the Galactic Nuclear Star Cluster. ApJ. 866, id.52.
  * [28] Cunha, K., Sellgren, K., Smith, V. V., Ramirez, S. V., Blum, R. D. and Terndrup, D. M. (2007) Chemical Abundances of Luminous Cool Stars in the Galactic Center from High-Resolution Infrared Spectroscopy. ApJ. 669, 1011
  * [29] Ramirez, S. V., Depoy, D. L., Frogel, Jay A., Sellgren, K., Blum, R. D. (1997) Luminosity and Temperature from Near-Infrared Spectra of Late-Type Giant Stars. AJ. 113, 1411–1420
  * [30] Asplund, M., Grevesse, N., Sauval, A. J. and Scott, P. (2009) The Chemical Composition of the Sun. ARA&A. 47, 481–522.
  * [31] Tautvaišienė, G., Mikolaitis, Š., Drazdauskas, A., Minkevičiūtė, R., Kjeldsen, H., Brogaard, K., et al. (2020) Chemical Composition of Bright Stars in the Continuous Viewing Zone of the TESS Space Mission. ApJS. 248, id.19.
  * [32] Demarque, P., Woo, J.-H., Kim, Y.-C. and Yi, S. K. (2004) Y2 Isochrones with an Improved Core Overshoot Treatment. ApJS. 155, 667–674.
  * [33] Thorsbro, B., Ryde, N., Rich, R. M., Schultheis, M., Renaud, F., Spitoni, E., et al. (2020) Detailed Abundances in the Galactic Center: Evidence of a Metal-rich Alpha-enhanced Stellar Population. ApJ. 894, id.26.
  * [34] Battistini, C., and Bensby, T. (2015) The origin and evolution of the odd-Z iron-peak elements Sc, V, Mn, and Co in the Milky Way stellar disk. A&A. 577, A9.
  * [35] Valenti, J. A., and Piskunov, N. (1996) Spectroscopy made easy: A new tool for fitting observations with synthetic spectra. A&AS. 118, 595-603.
  * [36] Valenti, J. A. and Piskunov, N. (2012) SME: Spectroscopy Made Easy. Astrophysics Source Code Library, record ascl:1202.013
  * [37] Gustafsson, B., Edvardsson, B., Eriksson, K., Jørgensen, U. G., Nordlund, Å. and Plez, B. (2008) A grid of MARCS model atmospheres for late-type stars. I. Methods and general properties. A&A. 486, 951–970.
  * [38] Thorsbro, B., Ryde, N., Rich, R. M., Schultheis, M. Fritz, T. K., and Origlia, L. (2017) Developing an astrophysical line list for Keck/Nirspec observations of red giants in the Galactic centre. Proceedings of the International Astronomical Union. 334. 372–373.
  * [39] Livingston, W., and Wallace, L. (1991) An atlas of the solar spectrum in the infrared from 1850 to 9000 cm${}^{-}{1}$ (1.1 to 5.4 micrometer). NSO Technical Report, Tucson: National Solar Observatory, National Optical Astronomy Observatory.
  * [40] Sneden, C., Lucatello, S., Ram, R. S., Brooke, J. S. A., and Bernath, P. (2014) Line Lists for the A 2$\Pi$-X 2$\Sigma$+ (Red) and B 2$\Sigma$+-X 2$\Sigma$+ (Violet) Systems of CN, 13C14N, and 12C15N, and Application to Astronomical Spectra. ApJS. 214. id26.
  * [41] Gautam, A. K., Do, T., Ghez, A. M., Morris, M. R., Martinez, G. D., Hosek, M. W. Jr., et al. (2019) An Adaptive Optics Survey of Stellar Variability at the Galactic Center. ApJ. 871, id.103.
  * [42] Schödel, R., Najarro, F., Muzic, K. and Eckart, A. (2010) Peering through the veil: near-infrared photometry and extinction for the Galactic nuclear star cluster. Accurate near infrared H, Ks, and L’ photometry and the near-infrared extinction-law toward the central parsec of the Galaxy. A&A. 511, id.A18.
  * [43] Blum, R. D., Ramírez, S. V., Sellgren, K. and Olsen, K. (2003) Really Cool Stars and the Star Formation History at the Galactic Center. ApJ. 597, 323–346.
  * [44] Marigo, P., Girardi, L., Bressan, A., Groenewegen, M. A. T., Silva, L. and Granato, G. L. (2008) Evolution of asymptotic giant branch stars. II. Optical to far-infrared isochrones with improved TP-AGB models. A&A. 482, 883–905.
  * [45] Pfuhl, O., Fritz, T. K., Zilka, M., Maness, H., Eisenhauer, F., Genzel, R., et al. (2011) The Star Formation History of the Milky Way’s Nuclear Star Cluster. ApJ. 741, id.108.
  * [46] Nishiyama, S., Schödel, R., Yoshikawa, T., Nagata, T., Minowa, Y. and Tamura, M. (2016) Spectroscopically identified intermediate age stars at 0.5-3 pc distance from Sagittarius A*. A&A. 588, id.A49.
  * [47] Chen, Z., Do, T., Ghez, A. M., Hosek, M. W., Feldmeier-Krause, A., Chu, D. S., et al. (2023) The Star Formation History of the Milky Way’s Nuclear Star Cluster. ApJ. 944, id.79.
  * [48] Rich, R. M., Ryde, N., Thorsbro, B., Fritz, T. K., Schultheis, M., Origlia, L., et al. (2017) Detailed Abundances for the Old Population near the Galactic Center. I. Metallicity Distribution of the Nuclear Star Cluster. AJ. 154, id.239.
  * [49] Pritzl, B. J., Venn, K. A. and Irwin, M. (2005) A Comparison of Elemental Abundance Ratios in Globular Clusters, Field Stars, and Dwarf Spheroidal Galaxies. AJ. 130, 2140–2165.
  * [50] Bensby, T., Feltzing, S. and Oey, M. S. (2014) Exploring the Milky Way stellar disk. A detailed elemental abundance study of 714 F and G dwarf stars in the solar neighbourhood. A&A. 562, id.A71.
  * [51] Bensby, T., Feltzing, S., Gould, A., Yee, J. C., Johnson, J. A., Asplund, M., et al. (2017) Chemical evolution of the Galactic bulge as traced by microlensed dwarf and subgiant stars. VI. Age and abundance structure of the stellar populations in the central sub-kpc of the Milky Way. A&A. 605, id.A89.
  * [52] Lucey, M., Hawkins, K., Ness, M., Nelson, T., Debattista, V. P.., Luna, A., et al. (2022) The COMBS Survey - III. The chemodynamical origins of metal-poor bulge stars. MNRAS. 509, 122–144.
  * [53] Hill, V., Skúladóttir, Á., Tolstoy, E., Venn, K. A., Shetrone, M. D., Jablonka, P., et al. (2019) VLT/FLAMES high-resolution chemical abundances in Sculptor: a textbook dwarf spheroidal galaxy. A&A. 626, id.A15.
  * [54] Monaco, L., Bellazzini, M., Bonifacio, P., Ferraro, F. R., Marconi, G., Pancino, E., et al. (2005) The Ital-FLAMES survey of the Sagittarius dwarf spheroidal galaxy. I. Chemical abundances of bright RGB stars. A&A. 441, 141–151.
  * [55] Van der Swaelmen, M., Hill, V., Primas, F. and Cole, A. A. (2013) Chemical abundances in LMC stellar populations. II. The bar sample. A&A. 560, id.A44.
  * [56] Dong, H., Schödel, R., Williams, B. F., Nogueras-Lara, F., Gallego-Cano, E., Gallego-Calvente, T., et al. (2017) Near-infrared variability study of the central 2.3 × 2.3 arcmin2 of the Galactic Centre - II. Identification of RR Lyrae stars in the Milky Way nuclear star cluster”. MNRAS. 471, 3617–3631.


# Hyperparameter optimization of orthogonal functions in the numerical
solution of differential equations

Alireza Afzal Aghaei<EMAIL_ADDRESS>Department of Computer and
Data Science, Faculty of Mathematical Science , Shahid Beheshti University,
Tehran, Iran.  Kourosh Parand<EMAIL_ADDRESS>Department of Computer and
Data Science, Faculty of Mathematical Science , Shahid Beheshti University,
Tehran, Iran.  Department of Cognitive Modeling, Institute for Cognitive and
Brain Sciences , Shahid Beheshti University, Tehran, Iran.

###### Abstract

This paper considers the hyperparameter optimization problem of mathematical
techniques that arise in the numerical solution of differential and integral
equations. The well-known approaches grid and random search, in a parallel
algorithm manner, are developed to find the optimal set of hyperparameters.
Employing rational Jacobi functions, we ran these algorithms on two nonlinear
benchmark differential equations on the semi-infinite domain. The
configurations contain different rational mappings along with their length
scale parameter and the Jacobi functions parameters. These trials are
configured on the collocation Least-Squares Support Vector Regression (CLS-
SVR), a novel numerical simulation approach based on spectral methods. In
addition, we have addressed the sensitivity of these hyperparameters on the
numerical stability and convergence of the CLS-SVR model. The experiments show
that this technique can effectively improve state-of-the-art results.

Keywords— Nonlinear differential equations, Hyperparameter optimization,
Jacobi polynomials, Machine learning

##  1 Introduction

Estimating the unknown dynamics of a given physical system is an essential
task in science and engineering. Mathematicians simulate these systems after
expressing them in functional equations such as differential and integral
equations. The failure of the analytical approaches in solving these problems,
which mainly contain nonlinear terms, has led researchers to develop various
numerical techniques. Finite Elements (FEM), Finite Volume (FVM), Meshless,
and Spectral methods are some well-known efforts. These techniques consider a
simple mathematical model defined by some parameters and hyperparameters. The
former are internal variables learned during the training process, whereas the
latter are external configurations that should be optimized to find the best
estimator. These terms are initially back in the machine learning literature.
For example, in a supervised machine learning task, the Support Vector Machine
(SVM) algorithm considers a hyperplane for separating the given data into
different classes. The vector that defines the hyperplane is a parameter,
whereas the kernel function is a hyperparameter. Likewise, a linear
combination of unknown weights and basis functions is usually considered for
approximating the solution to a differential equation. Here the unknown
weights are parameters, while the basis functions are a hyperparameter.
Choosing the best combination of hyperparameters is critical and may
significantly affect the result. Scientists take advantage of prior knowledge
to choose sub-optimal ones. As an example, it is common to use logarithmic-
spaced values for the regularization parameter of a machine learning model. In
a mathematics setting, the intrinsic nature of the problem is utilized to
choose appropriate hyperparameters. For instance, the periodic problems are
usually approximated by the Fourier basis functions, and the rational
functions may approximate the problems that reach their steady state, and so
on.

Regarding the development of artificial intelligence, hyperparameter
optimization or simply tuning and finding the optimal hyperparameters, has
evolved. Various approaches and techniques are proposed to handle this issue.
Exhaustive search [1], probabilistic algorithms [2], gradient-based [3] and
meta-heuristic algorithms [4] are some of these efforts. In a similar manner,
mathematicians tried to develop some routines to find optimal hyperparameters
that arises in the numerical simulation techniques. To name a few, Boyd [5]
proved some theorems for rational Chebyshev functions, Sanyasiraju et al. [6]
used a local optimization algorithm for Radial Basis Functions (RBFs),
Cavoretto et al. [7] combined leave-one-out cross-validation with univariate
global optimization techniques for RBFs, Tanguy et al. [8] presented a
quadratic minimization procedure for optimal placement of poles in rational
approximations by Müntz–Laguerre functions and Mi at al. [9] developed two
algorithms for the adaptive selection of continuous rational orthogonal basis
functions. However, it is valuable to provide methods for solving problems in
a broader range rather than just specific cases.

A wide range of physical problems involving initial and boundary value
problems are defined on the infinite or semi-infinite domain. Approximating an
accurate solution on the corresponding domain is essential. These problems are
mostly solved with orthogonal polynomials or rational functions with
completeness and orthogonality properties on the problem domain. The Hermite,
Laguerre, and rational Jacobi functions are the most widely-used functions
applied to these problems. Computational and mathematical properties of the
rational functions encouraged scientists to choose them as the basis functions
[10, 11, 12, 13, 14, 15, 16, 17, 18]. However, these functions suffer of
hyperparameters such as rational mapping and length scale parameter. These can
disturb the numerical solution in some cases [19].

In this research, we develop some algorithms and investigate the applications
of machine learning algorithms to optimize the hyperparameters that appear
during the numerical simulation of differential equations arising on semi-
infinite domains. In the continuation of the article, we will discuss the
preliminaries (section 2), the proposed method (section 3), and the state-of-
the-art numerical results (section 4). Finally, the concluding remarks will be
discussed in the last section.

##  2 Preliminaries

In this section, we explain some prerequisites needed in the following
sections. To do so, we first explain the Jacobi polynomials, then the CLS-SVR
method will be recalled. The hyperparameter optimization techniques used in
the rest of the work will be discussed.

###  2.1 Orthogonal Polynomials

Hermite, Laguerre, and Jacobi polynomials are the most well-known orthogonal
polynomials defined on infinite, semi-infinite, and finite intervals. They can
be used to approximate functions in corresponding domains. However,
approximating rational functions may not be very accurate using polynomials.
Therefore, researchers proposed various rational mappings to handle this
issue. In general, Jacobi polynomials with hyperparameters $\alpha,\beta$ are
defined on the interval $[-1,1]$ by the recursive expression

$\displaystyle
J_{i}^{\alpha,\beta}(x)=-\frac{(\alpha+i-1)(\beta+i-1)(\alpha+\beta+2i)}{i(\alpha+\beta+i)(\alpha+\beta+2i-2)}J_{i-2}^{\alpha,\beta}(x)$
(1)
$\displaystyle+\frac{(\alpha+\beta+2i-1)\left\\{\alpha^{2}-\beta^{2}+x(\alpha+\beta+2i)(\alpha+\beta+2i-2)\right\\}}{2i(\alpha+\beta+i)(\alpha+\beta+2i-2)}$
$\displaystyle\quad\times J_{i-1}^{\alpha,\beta}(x),\quad i=2,3,\ldots,$

where

$J_{0}^{\alpha,\beta}(x)=1,\quad
J_{1}^{\alpha,\beta}(x)=\frac{\alpha+\beta+2}{2}x+\frac{\alpha-\beta}{2}$

Their orthogonality is defined using the $L2$ inner product:

${\displaystyle\int_{-1}^{1}(1-x)^{\alpha}(1+x)^{\beta}J_{m}^{(\alpha,\beta)}(x)J_{n}^{(\alpha,\beta)}(x)\,dx={\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}}{\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)n!}}\delta_{nm}}.$

The Gegenbauer, Chebyshev, and Legendre polynomials are some special cases of
Jacobi polynomials. For Legendre polynomials, the equation (1) with
$\alpha=\beta=0$ reduces to:

$\displaystyle P_{0}(x)=1,\quad P_{1}(x)=x,$ (2)
$\displaystyle(n+1)P_{n+1}(x)=(2n+1)xP_{n}(x)-nP_{n-1}(x),\quad n\geq 1,$

and for Chebyshev with $\alpha=\beta=-\nicefrac{{1}}{{2}}$ we have

$\displaystyle T_{0}(x)=1,\quad T_{1}(x)=x,$ (3) $\displaystyle
T_{n+1}(x)=2x\,T_{n}(x)-T_{n-1}(x),\quad n\geq 1.$

Although these polynomials are defined for the bounded domain, researchers
have used some nonlinear maps $\phi$ with the property
$\phi:[-1,1]\to[0,\infty)$ to transform the orthogonality into the semi-
infinite domain. To our best knowledge, these three mappings are the most used
functions among researchers [20]:

  * •

Algebraic mapping:
$\displaystyle\phi(x)=\nicefrac{{(x-\theta)}}{{(x+\theta)}}$.

  * •

Exponential mapping:
$\displaystyle\phi(x)=1-2{\textrm{e}^{-\displaystyle\nicefrac{{x}}{{\theta}}}}$.

  * •

Logarithmic mapping:
$\displaystyle\phi(x)=2\tanh\left(\nicefrac{{x}}{{\theta}}\right)-1$.

The rational Jacobi functions is defined by a simple transformation
$\tilde{J}(x)=J(\phi(x))$ where $\phi(x)$ is a nonlinear rational mapping.
Therefore, orthogonality takes the form

${\displaystyle\int_{0}^{\infty}\tilde{J}_{m}^{(\alpha,\beta)}(x)\tilde{J}_{n}^{(\alpha,\beta)}(x)w(x)\,dx={\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}}{\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)n!}}\delta_{nm}},$
(4)

where $w(x)$ is the corresponding weight function of the inner product.

###  2.2 Collocation Least-Squares Support Vector Regression

Collocation Least-Squares Support Vector Regression (CLS-SVR) is a novel
formulation of support vector machines for solving functional equations [21,
22, 23]. In this machine learning model, the unknown solution of a
differential equation is approximated by a linear combination of unknown
coefficients and some known basis functions. The basis functions, also known
as feature maps, transform the input data into a nonlinear space in which we
hope the approximation accuracy will increase. In the following, we recall the
formulation of CLS-SVR which is based on the paper [24].

Suppose that we want to solve an arbitrary differential equation in the form
of $\mathcal{N}(u)=f$. Approximating the solution by $m$ basis functions, we
have:

$\tilde{u}(x)=w^{T}\varphi(x)+b=\sum_{i=1}^{m}w_{i}\varphi_{i}(x)+b.$

The primal optimization form of CLS-SVR takes the form:

$\displaystyle\min\limits_{w,e}$
$\displaystyle\frac{1}{2}w^{T}w+\frac{\gamma}{2}e^{T}e$ (5) s.t.
$\displaystyle\mathcal{N}(\tilde{u})(x_{k})-f(x_{k})=e_{k},\quad
k=1,\ldots,n,$ $\displaystyle\tilde{u}(c_{j})=u_{j},\quad j=1,\ldots,d.$

where $n$ is the number of training points, $e_{k}$ is the value of the
residual function at $k$-th training point, and $\tilde{u}(c_{j})=u_{j}$ is
the set of initial and boundary conditions. The regularization hyperparameter
$\gamma$ controls the fluctuations of the learned solution and reduces the
overfitting on the training data. However, this hyperparameter may not be
optimized in the case of problems with unique solutions.

The dual form of this optimization problem leads to a linear or nonlinear
system of equations. This system can be obtained using the Lagrangian function
for (5):

$\displaystyle\mathscr{L}(w,e,\alpha)=\frac{1}{2}w^{T}w+\frac{\gamma}{2}e^{T}e-\sum_{k=1}^{n}\alpha_{k}\left[\mathcal{N}(\tilde{u})(x_{k})-f(x_{k})-e_{k}\right]-\sum_{j=1}^{d}\beta_{j}\left[\tilde{u}(c_{i})-u_{i}\right].$

The saddle point of this function satisfies the solution of the problem:

$\left\\{\frac{\partial\mathscr{L}}{\partial
w_{i}}=0,\frac{\partial\mathscr{L}}{\partial
e_{k}}=0,\frac{\partial\mathscr{L}}{\partial\alpha_{k}}=0,\frac{\partial\mathscr{L}}{\partial\beta_{j}}=0\right\\},$

where $i=1,\cdots,m$, $j=1,\cdots,d$, and $k=1,\cdots,n$. After solving this
system, we can use the obtained weights $w$ to approximate the equation or use
the dual variables in a kernel method sense.

###  2.3 Hyperparameter optimization

The problem of hyperparameter optimization is usually expressed as a non-
convex minimization or maximization problem. Various algorithms proposed to
solve this problem can be categorized into gradient-based and gradient-free
sets. The former uses the gradient of a loss function to find an optimum
value. At the same time, by generating some candidates, the latter tries to
approximate the search space to find a local or possibly global optimum. Grid
search, random search, Bayesian optimization, and meta-heuristic algorithms
are some examples of gradient-free methods. Grid and random search have been
widely used because of their simplicity and acceptable performance. Moreover,
they can be massively run in parallel. Briefly, the grid search seeks the
Cartesian product of given parameter sets, whereas the random search samples a
fixed number of parameters from user-defined distributions. For categorical
cases, the hyperparameter is chosen uniformly. In both methods, the best
parameter set would be measured on a test dataset usually generated by the
K-Fold cross-validation algorithm. Figure 1 compares these two algorithms.

A key difference between grid and random search is their computational
complexity. The grid search seeks all possible values in $O(n^{k})$ time,
having $k$ hyperparameters and $n$ different values for each of them, while
for a random search, the user can define a budget based on the available
resources so that the algorithm will run on the time complexity of $O(n)$. For
a more precise comparison, we refer to the authors [25, 26].

Figure 1: . 1: A comparison between random and grid searches for a simple
search space with only two different parameters. It can be seen that grid
search may fail to explore the search space efficiently. The figure is adapted
from Bergstra et al. [1].

##  3 Method and Results

In this section, we explain the proposed algorithm and then provide examples
to show the methods’ efficiency. Here, we focus on the hyperparameter
optimization of orthogonal rational Jacobi functions, whereas the presented
algorithm can be easily extended to other mathematical models.

The optimal set of hyperparameters should be found on a test dataset to
prevent the overfitting problem. In machine learning, this is easily done
using cross-validation techniques. Nevertheless, in the mathematical
equations, there is no data to split into train and test sets. However, there
are three alternative options to handle this issue. The first option is to use
a set of new nodes (different from training points) in the domain and
calculate the prediction error on these nodes. The second is employing some
physical properties of the model. The last one is to use some criteria found
by mathematical techniques, including analytical or accurate numerical ones.
The last option can be seen as a composition of previous ones. Most of the
criteria measures used by the authors have important physical meanings which
are good enough to be accuracy test criteria. Therefore, we use the absolute
difference value between the exact or the state-of-the-art approximations and
the predicted one given by CLS-SVR. The proposed grid and random search
algorithms are presented in algorithms 1 and 2.

We will obtain the optimal numerical results for some well-known benchmark
nonlinear differential equations using the proposed method in the following
sections. The configuration used to find an accurate numerical solution for
these problems is reported in tables 2 and 2. It is seen that the grid search
seeks $600$ different configurations to find the optimal set. To have a fair
comparison, we set the maximum iterations for random search the same as the
number of grid nodes. In addition, the roots of Legendre and Chebyshev
polynomials are utilized as the training data.

Data: The differential equation As $ODE$

1 $S_{i}\leftarrow\text{Set of all desired values for }i\text{-th
hyperparameter}$

2 $CriteriaList\leftarrow[]$

3 $Parameter\\_Set\leftarrow CartesianProduct(\\{S_{i}:\forall i\\})$

4 for  _$param\\_set\textbf{ in }Parameter\\_Set\textit{ in parallel}$_ do

5 $\tilde{u}\leftarrow CLS-SVR(ODE)$

6 Compute Criteria for $\tilde{u}$

7 Push Criteria to $CriteriaList$

8

9 end for

Result: $parameter\\_set$ associated with best criteria

Algorithm 1 Grid search algorithm

Data: The differential equation As $ODE$

1 $S_{i}\leftarrow\text{A suitable distribution for }i\text{-th
hyperparameter}$

2 $CriteriaList\leftarrow[]$

3

4 for  _$iter\textbf{ from }1\textbf{ to }MAX\\_ITER\textit{ in parallel}$_ do

5 $\tilde{u}\leftarrow CLS-SVR(ODE)$

6 Compute Criteria for $\tilde{u}$

7 Push Criteria to $CriteriaList$

8

9 end for

Result: $parameter\\_set$ associated with best criteria

Algorithm 2 Random search algorithm

Parameter | Values
---|---
Kernel | {Legendre, Chebyshev}
Mapping | {Algebraic, Exponential, Logarithmic}
$\theta$ | $\\{0.1,0.2,\cdots,9.9,10\\}$
Table 1: The search space of the grid search

Parameter | Values
---|---
Kernel | {Legendre, Chebyshev}
Mapping | {Algebraic, Exponential, Logarithmic}
$\theta$ | UniformDistribution(0, 10)
Table 2: The search space of the random search

###  3.1 Volterra’s population model

Volterra’s population model is a nonlinear integro-differential equation that
describes population growth of a species within a closed system [24]. It is
common to express this problem into an equivalent nonlinear differential
equation [16]:

$\displaystyle\kappa{u^{\prime\prime}}(x)={u^{\prime}}(x)-{u^{\prime}}(x)^{2}-{u}(x){u^{\prime}}(x),$
(6) $\displaystyle{u}(0)=0,{u^{\prime}}(0)=0.1.$

Here $\kappa$ is a non-dimensional parameter. The criteria for the prediction
correctness of this problem is the maximum value of the approximated function.
TeBeest [27] showed that the maximum peak is:

$\displaystyle u_{max}=1+\kappa\ln({\frac{\kappa}{1+\kappa+{u(0)}}}).$ (7)

Considering (7) as the exact value for (6), the absolute error for this
equations is defined as

$\displaystyle\text{Absolute Error}\vcentcolon=\mid
u_{max}-\max_{t\in(0,\infty)}\tilde{u}(t)\mid.$ (8)

To find a reasonable range for the length scale parameter and the effect of
the other hyperparameters, we first ran a sensitivity analysis on these
hyperparameters on a large domain. The results are reported in the figure 2.
It can be seen that the large values for the length scale do not yield a good
approximation. The maximum reasonable value for this task is about $10$. In
addition, the basis functions do not impose any significant difference on this
parameter. Moreover, the nonlinear mappings can affect the computational
procedure of the optimization problem (5). This issue has resulted a
discontinuity in these figures. Figure 3 plots some of the successfully
learned approximated solutions. It is seen that the equation may not be
accurately simulated using improper hyperparameters. After the sensitivity
analysis, we simulate this problem with five different most used values for
the non-dimensional parameter $\kappa$ using the grid and random search
algorithms. As reported in tables 2 and 2, the interval $(0,10]$ is chosen for
the length scale parameter. Tables 3 and 4 report the best-obtained
hyperparameters. From there, it can be inferred that algebraic mapping is the
best choice for small values of $\kappa$ while for the larger values, the
exponential mapping could obtain better approximations. Moreover, the kernel
function and its nonlinear mapping can result in a bit different accuracy.

To show the effectiveness of the proposed algorithm, we have compared the
absolute error found by various authors in table 5. Rational Legendre (RL) and
Rational Chebyshev (RC) pseudospectral method [17], Sinc collocation method
(SCM) [16], and Fractional Rational Legendre (FRL) [24] are compared to each
other. The length scale values used by Parand et al. [24] are optimized in a
similar manner proposed by Boyd [19].

(a) Legendre rational functions

(b) Chebyshev rational functions

Figure 2: The effect of length scale parameter $\theta$ on a large domain for
Volterra’s population model with $m=25$ and $\kappa=0.5$. Figure 3: A set of
learned solutions with different values for length scale $\theta$ for
Volterra’s population model with $\kappa=0.5$.

$\kappa$ | Kernel | Mapping | $\theta$ |  Exact [27] | Approximate | Error
---|---|---|---|---|---|---
$0.02$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $0.10000$ | $0.92342717207022$ | $0.92342711545307$ | $5.662\times 10^{-08}$
$0.04$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $0.70000$ | $0.87371998300000$ | $0.87371998000000$ | $3.090\times 10^{-09}$
$0.10$ | Legendre | $\displaystyle 1-2\exp(\nicefrac{{-x}}{{\theta}})$ | $1.00000$ | $0.76974149100000$ | $0.76974149100000$ | $7.890\times 10^{-12}$
$0.20$ | Chebyshev | $\displaystyle 1-2\exp(\nicefrac{{-x}}{{\theta}})$ | $1.90000$ | $0.65905038200000$ | $0.65905038200000$ | $4.660\times 10^{-11}$
$0.50$ | Chebyshev | $\displaystyle 1-2\exp(\nicefrac{{-x}}{{\theta}})$ | $3.80000$ | $0.48519029140942$ | $0.48519029141289$ | $3.468\times 10^{-12}$

Table 3: The obtained results for Volterra population equation using a grid
search algorithm.

$\kappa$ | Kernel | Mapping | $\theta$ |  Exact [27] | Approximate | Error
---|---|---|---|---|---|---
$0.02$ | Chebyshev | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $0.539501186666072$ | $0.92342717207022$ | $0.92342733722160$ | $1.652\times 10^{-07}$
$0.04$ | Chebyshev | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $0.318328463774207$ | $0.87371998315400$ | $0.87371998508417$ | $1.930\times 10^{-09}$
$0.10$ | Chebyshev | $\displaystyle 1-2\exp(\nicefrac{{-x}}{{\theta}})$ | $1.626117351946306$ | $0.76974149070060$ | $0.76974149073275$ | $3.216\times 10^{-11}$
$0.20$ | Chebyshev | $\displaystyle 1-2\exp(\nicefrac{{-x}}{{\theta}})$ | $2.510838579760311$ | $0.65905038155232$ | $0.65905038153414$ | $1.818\times 10^{-11}$
$0.50$ | Legendre | $\displaystyle 2\tanh(\nicefrac{{x}}{{\theta}})-1$ | $6.797026768536748$ | $0.48519029140942$ | $0.48519029141142$ | $2.007\times 10^{-12}$

Table 4: The obtained results for Volterra population equation using a random search algorithm. $\kappa$ |  RL [17] |  RC [17] |  FRL [24] |  SCM [16] | Presented Method
---|---|---|---|---|---
$m$ | $50$ | $50$ | $40$ | $35$ | $40$
$0.02$ | $3.72\times 10^{-07}$ | $7.51\times 10^{-07}$ | $6.33\times 10^{-06}$ | $7.00\times 10^{-08}$ | $5.66\times 10^{-08}$
$0.04$ | $1.43\times 10^{-08}$ | $5.27\times 10^{-08}$ | $6.75\times 10^{-08}$ | $3.00\times 10^{-08}$ | $1.93\times 10^{-09}$
$0.10$ | $1.07\times 10^{-10}$ | $2.13\times 10^{-10}$ | $2.73\times 10^{-08}$ | $8.00\times 10^{-08}$ | $7.89\times 10^{-12}$
$0.20$ | $3.53\times 10^{-11}$ | $2.33\times 10^{-10}$ | $8.57\times 10^{-10}$ | $3.30\times 10^{-07}$ | $1.82\times 10^{-11}$
$0.50$ | $2.44\times 10^{-09}$ | $4.87\times 10^{-09}$ | $2.73\times 10^{-10}$ | $3.40\times 10^{-07}$ | $2.01\times 10^{-12}$
Table 5: A comparison among some mathematical methods solved Volterra’s
population model on the semi-infinite domain.

###  3.2 Kidder equation

The unsteady isothermal flow of a gas through a micro-nano porous medium can
be expressed as a nonlinear differential equation [18]. This equation which is
defined on a semi-infinite domain is modeled as:

$\displaystyle u^{\prime\prime}(x)+\frac{2x}{\sqrt{1-\kappa
u(x)}}u^{\prime}(x)=0,$ (9) $\displaystyle u(0)=1,u(\infty)=0.$

The initial slope of the approximated solution is an essential measure of the
accuracy of the problem. Up to now, no exact solution or initial slope has
been found for this problem. However, some researchers developed advanced
techniques to find an accurate solution. In this research, we use the values
obtained by Parand et al. [20] as an exact approximation. This paper has
obtained the exact initial slope up to $38$ digits of precision in a machine
or software that supports arbitrary precision arithmetic. Here, we are just
focusing on the problem of choosing the best hyperparameters. Thus, we compare
our results on a small number of basis and an almost similar number of basis
functions. Furthermore, they have utilized the Quasi-Linearization Method
(QLM) to convert the problem of approximating solutions for nonlinear Ordinary
Differential Equations (ODEs) into a sequence of dependent linear ODEs. This
increases the computational complexity to the number of QLM iterations. Here
we have just solved the original nonlinear problem, which is more
computationally efficient.

Same as the previous example, we first analyze the effect of the
hyperparameters. Figure 4 demonstrates the absolute error obtained by
different sets of hyperparameters. It is seen that the Legendre kernel can
reach better results in comparison to the Chebyshev functions. Furthermore,
some of the hyperparameters lead to numerical issues with the Legendre basis
functions, and hence the plot is discontinuous. Likewise in the previous
example, the interval $(0,10]$ for the length scale contains the best results,
therefore we focus on this range in the next experiments. Tables 6 and 7
reported the best results obtained by the grid and random search,
respectively. From there, it can be seen that the Legendre with algebraic
mapping is the best-obtained hyperparameters in all of the configurations for
different parameters. Furthermore, the random search algorithm overcame the
grid search in all of the experiments.

A comparison is made between the presented method and other related works in
table 8. Some of the simulated approximations are plotted in the figure 5.

(a) Legendre rational functions

(b) Chebyshev rational functions

Figure 4: The effect of length scale parameter $\theta$ on a large domain for
Kidder equation with $m=25$ Figure 5: A set of learned solutions with
different values for length scale $\theta$ for the Kidder equation with
$\kappa=0.5$.

$\kappa$ | Kernel | Mapping | $\theta$ |  Exact [20] | Approximate | Error
---|---|---|---|---|---|---
$0.10$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $4.80000$ | $-1.13900720617830$ | $-1.13900783881046$ | $6.326\times 10^{-07}$
$0.30$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $4.80000$ | $-1.16294145829591$ | $-1.16294198206589$ | $5.238\times 10^{-07}$
$0.50$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $4.80000$ | $-1.19179064971942$ | $-1.19179149343671$ | $8.437\times 10^{-07}$
$0.90$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $6.20000$ | $-1.28188132220336$ | $-1.28188069301075$ | $6.292\times 10^{-07}$

Table 6: The obtained results for the Kidder equation using a grid search
algorithm.

$\kappa$ | Kernel | Mapping | $\theta$ |  Exact [20] | Approximate | Error
---|---|---|---|---|---|---
$0.10$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $0.975404049994095$ | $-1.13900720617830$ | $-1.13900677991113$ | $4.263\times 10^{-07}$
$0.30$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $1.576130816775483$ | $-1.16294145829591$ | $-1.16294118714719$ | $2.711\times 10^{-07}$
$0.50$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $2.760250769985784$ | $-1.19179064971942$ | $-1.19179065703565$ | $7.316\times 10^{-09}$
$0.90$ | Legendre | $\displaystyle\nicefrac{{(x-\theta)}}{{(x+\theta)}}$ | $4.693906410582058$ | $-1.28188132220336$ | $-1.28188182747538$ | $5.053\times 10^{-07}$

Table 7: The obtained result for the Kidder equation using a random search
algorithm.

$\kappa$ |  CPA [28] |  Sinc [29] |  RL [29] |  HFC [30] |  Bessel [31] |  RJ [32] | Presented method
---|---|---|---|---|---|---|---
$m$ | $-$ | 32 | 32 | 20 | 26 | 20 | 25
$0.1$ | $2.94\times 10^{-04}$ | $-$ | $-$ | $-$ | $2.61\times 10^{-06}$ | $1.29\times 10^{-04}$ | $4.26\times 10^{-07}$
$0.5$ | $9.15\times 10^{-04}$ | $3.10\times 10^{-03}$ | $3.10\times 10^{-03}$ | $7.98\times 10^{-03}$ | $3.42\times 10^{-05}$ | $1.17\times 10^{-04}$ | $7.32\times 10^{-09}$
$0.9$ | $1.72\times 10^{-02}$ | $-$ | $-$ | $-$ | $3.65\times 10^{-06}$ | $1.02\times 10^{-04}$ | $5.05\times 10^{-07}$

Table 8: Comparison among various approximate initial slopes for the Kidder
equation (9). The values reported in [18] are assumed as the exact values.

##  4 Conclusion

This paper developed two machine learning techniques for increasing the
accuracy of the numerical simulation of functional equations. The presented
algorithms 1 and 2 are general tuning algorithms capable of hyperparameter
optimization of various mathematical approaches such as spectral methods,
RBFs, and wavelets. However, in this research, we have just focused on the
spectral method for approximating the solution of nonlinear differential
equations on the semi-infinite domain. To do so, we configured the search
space to the hyperparameters of rational Jacobi functions, including basis
functions, nonlinear rational mappings, and the length scale parameter.
Finally, in the numerical results, various experiments have been conducted to
measure the search capability of the proposed algorithms. We discussed the
role of the length scale parameter on the stability and convergence of the
method. In addition, some comparisons among related works were carried out to
show the superiority of these algorithms over traditional mathematical
procedures. This property, along with the small computation complexity handled
with parallel programming approaches, made this process efficient and easy to
use for other researchers. However, modern gradient-free global optimization
techniques, such as Bayesian optimization and Tree-structured Parzen
estimator, can be developed to get better approximations.

## References

  * [1] J. Bergstra and Y. Bengio, “Random search for hyper-parameter optimization,” Journal of Machine Learning Research, vol. 13, no. 10, pp. 281–305, 2012\.
  * [2] J. Snoek, H. Larochelle, and R. P. Adams, “Practical bayesian optimization of machine learning algorithms,” in Advances in Neural Information Processing Systems (F. Pereira, C. Burges, L. Bottou, and K. Weinberger, eds.), vol. 25, Curran Associates, Inc., 2012.
  * [3] D. Maclaurin, D. Duvenaud, and R. Adams, “Gradient-based hyperparameter optimization through reversible learning,” in International conference on machine learning, pp. 2113–2122, PMLR, 2015.
  * [4] F. Itano, M. A. d. A. de Sousa, and E. Del-Moral-Hernandez, “Extending mlp ann hyper-parameters optimization by using genetic algorithm,” in 2018 International joint conference on neural networks (IJCNN), pp. 1–8, IEEE, 2018\.
  * [5] J. P. Boyd, “The optimization of convergence for chebyshev polynomial methods in an unbounded domain,” Journal of computational physics, vol. 45, no. 1, pp. 43–79, 1982.
  * [6] Y. Sanyasiraju and C. Satyanarayana, “On optimization of the rbf shape parameter in a grid-free local scheme for convection dominated problems over non-uniform centers,” Applied Mathematical Modelling, vol. 37, no. 12-13, pp. 7245–7272, 2013.
  * [7] R. Cavoretto, A. De Rossi, M. S. Mukhametzhanov, and Y. D. Sergeyev, “On the search of the shape parameter in radial basis functions using univariate global optimization methods,” Journal of Global Optimization, vol. 79, no. 2, pp. 305–327, 2021.
  * [8] N. Tanguy, N. Iassamen, M. Telescu, and P. Cloastre, “Parameter optimization of orthonormal basis functions for efficient rational approximations,” Applied Mathematical Modelling, vol. 39, no. 16, pp. 4963–4970, 2015.
  * [9] W. Mi and W. X. Zheng, “Adaptive rational orthogonal basis functions for identification of continuous-time systems,” IEEE Transactions on Automatic Control, vol. 66, no. 4, pp. 1809–1816, 2021.
  * [10] T. Tajvidi, M. Razzaghi, and M. Dehghan, “Modified rational legendre approach to laminar viscous flow over a semi-infinite flat plate,” Chaos, Solitons & Fractals, vol. 35, no. 1, pp. 59–66, 2008.
  * [11] M. Khader and M. Adel, “Numerical approach for solving the riccati and logistic equations via qlm-rational legendre collocation method,” Computational and Applied Mathematics, vol. 39, no. 3, pp. 1–9, 2020.
  * [12] A. Saadatmandi and Z. Sanatkar, “Collocation method based on rational legendre functions for solving the magneto-hydrodynamic flow over a nonlinear stretching sheet,” Applied Mathematics and Computation, vol. 323, pp. 193–203, 2018.
  * [13] M. A. Abd El Salam, M. A. Ramadan, M. A. Nassar, P. Agarwal, and Y.-M. Chu, “Matrix computational collocation approach based on rational chebyshev functions for nonlinear differential equations,” Advances in Difference Equations, vol. 2021, no. 1, pp. 1–17, 2021.
  * [14] S. Deniz and M. Sezer, “Rational chebyshev collocation method for solving nonlinear heat transfer equations,” International Communications in Heat and Mass Transfer, vol. 114, p. 104595, 2020.
  * [15] X. Zhang and J. P. Boyd, “Revisiting the thomas–fermi equation: Accelerating rational chebyshev series through coordinate transformations,” Applied Numerical Mathematics, vol. 135, pp. 186–205, 2019.
  * [16] K. Parand, Z. Delafkar, N. Pakniat, A. Pirkhedri, and M. K. Haji, “Collocation method using sinc and rational legendre functions for solving volterra’s population model,” Communications in Nonlinear Science and Numerical Simulation, vol. 16, no. 4, pp. 1811–1819, 2011.
  * [17] M. Dehghan and M. Shahini, “Rational pseudospectral approximation to the solution of a nonlinear integro-differential equation arising in modeling of the population growth,” Applied Mathematical Modelling, vol. 39, no. 18, pp. 5521–5530, 2015.
  * [18] K. Parand, S. Latifi, M. Delkhosh, and M. M. Moayeri, “Generalized lagrangian jacobi gauss collocation method for solving unsteady isothermal gas through a micro-nano porous medium,” The European Physical Journal Plus, vol. 133, no. 1, pp. 1–12, 2018.
  * [19] J. P. Boyd, Chebyshev and Fourier spectral methods.  Courier Corporation, 2001.
  * [20] K. Parand and M. Delkhosh, “An Accurate Numerical Method for Solving Unsteady Isothermal Flow of a Gas Through a Semi-Infinite Porous Medium,” Journal of Computational and Nonlinear Dynamics, vol. 13, p. 011007, 10 2017\.
  * [21] A. G. Khoee, K. M. Mohammadi, M. Jani, and K. Parand, “A least squares support vector regression for anisotropic diffusion filtering,” arXiv preprint arXiv:2202.00595, 2022.
  * [22] P. Ahadian and K. Parand, “Support vector regression for the temperature-stimulated drug release,” Chaos, Solitons & Fractals, vol. 165, p. 112871, 2022.
  * [23] J. A. Rad, K. Parand, and S. Chakraverty, Learning with Fractional Orthogonal Kernel Classifiers in Support Vector Machines: Theory, Algorithms and Applications.  Springer Singapore, 2023.
  * [24] K. Parand, A. Aghaei, M. Jani, and A. Ghodsi, “Parallel ls-svm for the numerical simulation of fractional volterra’s population model,” Alexandria Engineering Journal, vol. 60, no. 6, pp. 5637–5647, 2021.
  * [25] L. Yang and A. Shami, “On hyperparameter optimization of machine learning algorithms: Theory and practice,” Neurocomputing, vol. 415, pp. 295–316, 2020.
  * [26] T. Yu and H. Zhu, “Hyper-parameter optimization: A review of algorithms and applications,” arXiv preprint arXiv:2003.05689, 2020.
  * [27] K. G. TeBeest, “Classroom note: numerical and analytical solutions of volterra’s population model,” SIAM review, vol. 39, no. 3, pp. 484–493, 1997.
  * [28] R. Iacono and J. P. Boyd, “The kidder equation,” Studies in Applied Mathematics, vol. 135, no. 1, pp. 63–85, 2015.
  * [29] A. Rezaei, K. Parand, and A. Pirkhedri, “Numerical study on gas flow through a micro-nano porous media based on special functions,” Journal of Computational and Theoretical Nanoscience, vol. 8, no. 2, pp. 282–288, 2011\.
  * [30] J. Rad, S. Ghaderi, and K. Parand, “Numerical and analytical solution of gas flow through a micro-nano porous media: A comparison,” Journal of Computational and Theoretical Nanoscience, vol. 8, no. 10, pp. 2033–2041, 2011\.
  * [31] K. Parand and M. Nikarya, “Solving the unsteady isothermal gas through a micro-nano porous medium via bessel function collocation method,” Journal of Computational and Theoretical Nanoscience, vol. 11, no. 1, pp. 131–136, 2014.
  * [32] K. Parand, P. Mazaheri, M. Delkhosh, and A. Ghaderi, “New numerical solutions for solving kidder equation by using the rational jacobi functions,” SeMA Journal, vol. 74, no. 4, pp. 569–583, 2017.