# MUSE narrow field mode observations of the central kinematics of M15

Christopher Usher,1 Sebastian Kamann,2 Mark Gieles,3,4 Vincent Hénault-
Brunet,5 Emanuele Dalessandro,6 Eduardo Balbinot,7 Antonio Sollima6
1Department of Astronomy, Oskar Klein Centre, Stockholm University, AlbaNova
University Centre, SE-106 91 Stockholm, Sweden
2Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK
3ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain
4Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),
Martí Franquès 1, E08028 Barcelona, Spain
5Department of Astronomy and Physics, Saint Mary’s University, 923 Robie
Street, Halifax, NS B3H 3C3, Canada
6INAF - Astrophysics and Space Science Observatory Bologna, Via Gobetti 93/3,
I-40129 Bologna, Italy
7Kapteyn Astronomical Institute, University of Groningen, Postbus 800,
NL-9700AV Groningen, the Netherlands email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present observations of the stellar kinematics of the centre of the core
collapsed globular cluster M15 obtained with the MUSE integral field
spectrograph on the VLT operating in narrow field mode. Thanks to the use of
adaptive optics, we obtain a spatial resolution of 0.1 arcsec and are able to
reliably measure the radial velocities of 864 stars within 8 arcsec of the
centre of M15 thus providing the largest sample of radial velocities ever
obtained for the innermost regions of this system. Combined with previous
observations of M15 using MUSE in wide field mode and literature data, we find
that the central kinematics of M15 are complex with the rotation axis of the
core of M15 offset from the rotation axis of the bulk of the cluster. While
this complexity has been suggested by previous work, we confirm it at higher
significance and in more detail.

###### keywords:

stars: kinematics and dynamics - globular clusters: individual: NGC 7078 (M15)

††pubyear: 2019††pagerange: MUSE narrow field mode observations of the central
kinematics of M15–MUSE narrow field mode observations of the central
kinematics of M15

##  1 Introduction

The high stellar densities of globular clusters (GCs) have long proved a
challenge in studying the kinematics of their resolved stars, particularly in
their centres. However, recent improvements in instrumentation and analysis
techniques have allowed the kinematics of large numbers of GC stars to be
studied. Imaging from space using the Hubble Space Telescope has allowed
proper motion measurements of GCs in their centres (e.g. Bellini et al., 2014;
Watkins et al., 2015) while Gaia provides high quality astrometry in cluster
outskirts (e.g. Gaia Collaboration et al., 2018; Bianchini et al., 2018;
Sollima et al., 2019; Vasiliev, 2019). Recent advances in adaptive optics also
provide the potential for useful astrometry in the centres of GCs (e.g.
Dalessandro et al., 2016; Massari et al., 2016).

Multiobject spectroscopy has allowed large samples of radial velocities to be
measured (e.g. Lane et al., 2011; Kimmig et al., 2015; Ferraro et al., 2018)
but conventional fibre or slit spectroscopy is limited to the lower density
outskirts. Integral field spectroscopy (e.g. Lanzoni et al., 2013; Kamann et
al., 2014; Dalgleish et al., 2020) allows spectra to be observed even in dense
fields, especially when used with techniques to reliably extract spectra from
blended stars (Kamann et al., 2013). With wide-field integral field
spectrographs, large numbers of spectra can be obtained at once, with Kamann
et al. (2018a) measuring the radial velocities of over 200 000 stars in 22 GCs
using MUSE on the VLT (Bacon et al., 2010).

M15 (NGC 7078) is an old ($\sim 13$ Gyr, e.g. Dotter et al. 2010; VandenBerg
et al. 2013), metal poor ([Fe/H] $=-2.33$ Carretta et al. 2009) GC with a mass
of $\sim 5\times 10^{5}\,{\rm M}_{\odot}$ (e.g. Dalgleish et al., 2020).
According to the 2010 edition of the Harris (1996) catalgoue, M15 lies 10.4
kpc from the Sun as well as 10.4 kpc from the Galactic centre. Its steep,
cuspy surface brightness profile strongly suggests that it is a core collapse
GC (e.g. Djorgovski & King, 1986).

To date, many groups investigated the kinematics of this GC through radial
velocity (Peterson et al., 1989; Gebhardt et al., 1994, 2000; Gerssen et al.,
2002) and proper motion (Bellini et al., 2014) surveys.

Interest on the central kinematics of M15 (e.g. Newell et al., 1976; Gerssen
et al., 2002; McNamara et al., 2003; Baumgardt et al., 2003; Murphy et al.,
2011; den Brok et al., 2014) initially focused on the possible presence of an
intermediate mass black hole (IMBH) in the cluster. However, it is now
generally considered unlikely that a core collapsed GC could host an IMBH as
the excess kinetic energy that stars gain in interactions with the IMBH
prevent their accumulation in the cluster centre (e.g. Illingworth & King,
1977; Baumgardt et al., 2003; Gieles et al., 2018).

More recent work has focused on the surprising detection of rotation in the
central regions of M15. van den Bosch et al. (2006) found evidence for a fast
spinning, kinematically decoupled core in the inner 4 arcsec of M15. This
result, already suggested by the analysis of Gebhardt et al. (1997), was later
supported by the work of Bianchini et al. (2013) and Kamann et al. (2018a) who
found the core of M15 shows a higher rotation and a different rotation axis
than what is found at large radii. This discovery was surprising given that
ordered rotation should be quickly erased by the short relaxation time scale
($<100$ Myr Harris 1996) in the centre of M15. Kamann et al. (2018a) found a
clear connection between the angular momentum of a GC and its relaxation time
with GCs with shorter relaxation times having lower specific angular momenta
(see also Bianchini et al., 2018; de Boer et al., 2019).

In this paper we combine the techniques of adaptive optics and integral field
spectroscopy, by analysing MUSE observations in the AO-supported narrow field
mode to study the central kinematics of M15. In Section 2 we describe our
observations and data reduction as well as our radial velocity measurements
and the literature datasets we include in our analysis. In Section 3 we
present the spatially resolved inner kinematics of M15 before discussing these
kinematics in 4.

##  2 Observations and data reduction

###  2.1 MUSE narrow field mode

We observed the centre of M15 using MUSE on the VLT in narrow field mode (NFM)
as part of the NFM science verification (program 60.A-9492(A)) on the nights
of 2018 September 13 and 14. We observed four pointings in a 2 by 2 mosaic
centred on the cluster centre. Two of the pointings were observed for $5\times
250$ s; the remaining two $4\times 250$ s. The DIMM seeing ranged from 0.5 to
1.2 arcsec. The observations were all carried out using the GALACSI adaptive
optics module (Arsenault et al., 2012; Ströbele et al., 2012) which provides
MUSE a 7.5 by 7.5 arcsec field-of-view with a spatial sampling of 0.025 arcsec
in NFM. In nominal mode, MUSE provides wavelength coverage between 4800 Å and
9300 Å with a gap between 5780 and 6050 Å due to the notch filter used to
suppress the light from the laser guide stars. MUSE has a near constant
spectral FWHM of $\sim 2.5$ Å across the wavelength range (corresponding to a
resolution of $R\sim 1700$ at the blue end of the wavelength range and $R\sim
3500$ at the red end). The datacubes were reduced using the standard MUSE
pipeline (Weilbacher et al., 2020) in the same manner as in Kamann et al.
(2018b).

Figure 1: Comparison of an $i$-band image created from one of the MUSE NFM
pointings (left) and a HST ACS F814W image (right). The MUSE NFM data shows a
similar spatial resolution to the HST data. Both images measure 7.98 by 7.88
arcsec.

We used PampelMuse (Kamann et al., 2013) to extract the stellar spectra from
the MUSE datacubes. PampelMuse determines the point spread function (PSF) from
bright, isolated stars in the datacube and uses the PSF to extract spectra of
all the sources in its input catalogue. We utilised the Beccari et al. (2019)
HST ACS/HRC photometric catalogue as the input catalogue. During the
extraction, we modelled the PSF using the combination of an outer Moffat and a
central Gaussian profile. The idea behind this hybrid PSF model is that AO
observations typically result in a PSF consisting of a diffraction-limited
core surrounded by a seeing-limited halo. For the Gaussian component, we
obtain a FWHM of $0\aas@@fstack{\prime\prime}08$ at $H_{\rm\alpha}$ and
$0\aas@@fstack{\prime\prime}06$ at the calcium triplet. We compare a $i$-band
image created from one of the MUSE NFM pointings with the HST ACS/WFC F814W
image from the ACS Globular Cluster Treasury Survey (GO-10775 Sarajedini et
al., 2007) in Figure 1. In this wavelength region, the spatial resolution of
the MUSE and HST data is is similar, with the MUSE PSF showing a sharper core
but broader wings than the HST PSF.

Initial radial velocities were obtained by cross correlating the extracted
spectra with synthetic templates from the PHOENIX library of Husser et al.
(2013). In order to find a suitable template for each spectrum, we compared
the photometry of Beccari et al. (2019) to an isochrone from the MIST (Choi et
al., 2016) database, adopting a cluster metallicity of $[{\rm Fe/H}]=-2.37$
(Harris, 1996). The comparison provided us with estimates of the effective
temperature and surface gravity of each star with a MUSE spectrum available,
which were in turn used to select a synthetic template for the cross
correlation.

In the final step of the spectrum analysis, we processed each spectrum with
Spexxy (Husser et al., 2016). The code performs a full spectrum fit against
the same library of synthetic spectra used for the cross correlation to
measure the radial velocity, effective temperature, and metallicity of each
star. Initial guesses for the fit were obtained from the isochrone comparison
described above as well as from the cross correlation.

As detailed in Kamann et al. (2016), the radial velocities obtained with
Spexxy tend to be more precise than those measured via cross correlation,
hence we gave preference to the former. We obtained reliable radial velocities
for 864 stars. This is a significant increase compared to previous studies in
the number of stars with radial velocities in the central few arcsec of M15
(see below). Owing to the large range in brightness of our target stars,
spectra were extracted over a wide range of signal-to-noise ratios, leading to
very different uncertainties of the individual radial velocities. Our median
accuracy is $7.3\,{\rm km\,s^{-1}}$, while the 16th and 84th percentiles of
the uncertainties distribution are $2.6\,{\rm km\,s^{-1}}$ and $11.9\,{\rm
km\,s^{-1}}$, respectively.

###  2.2 Literature radial velocities

We used the radial velocity catalogue from Kamann et al. (2018a) which used
MUSE in wide field mode (WFM) to observe four pointings covering M15 out to
$\sim 1$ arcmin at the same spectral resolution of the NFM dataset. The MUSE
WFM dataset was reduced and analysed in the same manner as the MUSE NFM data,
providing radial velocities for 10420 stars.

We also utilised the literature radial velocities from Gebhardt et al. (2000),
Gerssen et al. (2002) and Baumgardt & Hilker (2018). Gebhardt et al. (2000)
used an imaging Fabry-Perot spectrograph together with the CFHT Adaptive
Optics Bonnette (Rigaut et al., 1998) to obtain radial velocities of 104 stars
within 20 arcsec of the centre of M15. We note that these observations
required 6.5 hours of exposure time. Gebhardt et al. (2000) combine this
sample with their previous natural seeing Fabry-Perot observations and
literature data to provide a dataset of 1779 radial velocities. We further
note that the astrometry of this catalogue is poor, especially in the
outskirts of M15.

Gerssen et al. (2002) used multiple long slit HST STIS observations to mimic
an integral field spectrograph. They used 25 orbits of HST time to measure
radial velocities for 64 stars.

The Baumgardt & Hilker (2018) sample combines 273 radial velocities from
APOGEE (Abolfathi et al., 2018) with their own radial velocity measurements of
archival VLT FLAMES (Pasquini et al., 2002) and Keck HIRES (Vogt et al., 1994)
spectra. As these spectra were observed using traditional long slit or fibre
spectroscopy, they all lie more than 30 arcsec away from the cluster centre.

We compare the spatial distribution of stars from the first four catalogs in
the inner 15 by 15 arcsec in Figure 2. We plot the radial distributions of
stars from the five radial velocity catalogs in the upper panel of Figure 3.
Our MUSE NFM observations provide significantly more radial velocities within
8 arcsec compared to previous works. The MUSE NFM observations provide 33
stars with radial velocities within 1 arcsec of the centre, 191 within 3
arcsec and 801 within 8 arcsec compared to 7, 57 and 363 stars respectively in
the MUSE WFM observation within the same radii and 14, 57 and 173 stars within
the same radii from the combined Gebhardt et al. (2000) and Gerssen et al.
(2002) samples.

In the lower panel of Figure 3 we show the HST ACS/HRC F435W luminosity
functions (using the photometry from Beccari et al. 2019) for each of the
radial velocity catalogs within 3 arcsec of the centre of M15. Our MUSE NFM
observations provide significantly more radial velocities for fainter stars
than previous work. We are able to measure radial velocities for virtually all
stars in the ACS/HRC input catalogue brighter than F435W$\sim 19$.

Figure 2: Spatial comparison of the stars from different radial velocity
sources. Each panel shows the HST ACS F814W image of the central 15 by 15
arcsec of M15 over plotted by stars with our MUSE NFM measurements (top left
in blue), the HST STIS measurements from Gerssen et al. (2002, top right in
orange), the MUSE WFM measurements from Kamann et al. (2018a, bottom left in
green) and the CFHT Fabry-Perot measurements from Gebhardt et al. (2000,
bottom right in red). None of the velocity measurements from Baumgardt &
Hilker (2018) are within the inner 15 by 15 arcsec. Our MUSE NFM observations
allow us to measure radial velocities for significantly more stars compared to
the MUSE WFM or CFHT FP observations in the same footprint and to measure
radial velocities over a much larger area than the HST STIS observations.
Figure 3: Top: cumulative radial distribution of stars from different radial
velocity sources. Our MUSE NFM observations provide a dramatic increase in the
number of stars with radial velocity measurements within the inner few arcsec
of M15. Bottom: cumulative HST ACS F435W apparent magnitude distribution of
stars from different radial velocity sources within 3 arcsec of the cluster
centre. The black dotted line shows the magnitude distribution of all stars
detected within 3 arcsec in the F435W imaging. The MUSE NFM observations allow
us to measure radial velocities for fainter stars in the centre of M15 than
previous studies. The MUSE NFM observations appear to be as complete as the
HST photometry to F435W magnitude of $\sim 19$.

In Figure 4 we show a comparison of the MUSE NFM velocities and the other
velocities. When combining the velocity catalogues we applied shifts to each
of the other catalogues such that the median velocities of stars in common
were the same as the median for the same stars in the MUSE NFM sample. For the
312 stars in common between the MUSE WFM and MUSE NFM datasets, 68 % of the
velocity differences lie between $-11.5$ and $8.9$ km s-1, for the 54 stars in
common the HST STIS and MUSE NFM, 68 % of the velocity differences lie between
$-12.6$ and $3.8$ km s-1, and for the 126 stars in common the CFHT and MUSE
NFM, 68 % of the velocity differences lie between $-11.0$ and $9.8$ km s-1.
The widths of the distributions of the velocity differences are roughly two
times larger than would be expected from the uncertainties. The median
uncertainties on the MUSE NFM velocities are similar to those of the MUSE WFM
velocities (3.5 versus 3.2 km s-1) but smaller than the HST STIS (2.8 versus
5.5 km s-1) or the CFHT FP (2.0 versus 3.8 km s-1) velocity uncertainties for
stars in common between the studies.

To check whether the uncertainties are systematically larger or smaller for
one sample or another, we used the same Markov chain Monte Carlo kinematics
analysis presented in Section 3 to fit the median velocity and the velocity
dispersion for the stars in common between the other velocity source and the
MUSE NFM using both the NFM velocities and the other sources for each of the
MUSE WFM, HST STIS and CFHT FP datasets. In all cases the fitted kinematics
agree within the uncertainties suggesting that the uncertainties of none of
the samples are systematically different.

We expect that stars in short period binaries would have radial velocities
that would vary between the different studies, inflating the differences in
velocities beyond what is expected from their uncertainties. The binary
fraction is also typically higher in the centres of GCs than in their
outskirts (e.g. Milone et al., 2012; Ji & Bregman, 2015; Giesers et al.,
2019).

Figure 4: Comparison of the MUSE NFM velocities with velocities from other
sources. We show the differences between the MUSE NFM velocities and the MUSE
WFM velocities (top), the HST STIS velocities (middle) and the CFHT FP
velocities (bottom) as a function of the MUSE NFM velocities for the stars in
common between the studies. The scatter between studies is larger than the
uncertainties.

##  3 Inner kinematics of M15

We analysed the spatial resolved kinematics in a similar manner as in Kamann
et al. (2018a, 2020), by deriving 1D radial profiles and 2D Voronoi maps of
the first and second moments of the line-of-sight velocity distribution. The
radial profiles are depicted in Fig. 5. The continuous curves were obtained by
describing the unbinned data using analytical models. More specific, we
adopted the rotation profile expected for systems that underwent violent
relaxation (Lynden-Bell, 1967; Gott, 1973, eq. 1 in Kamann et al. 2020), while
a Plummer (1911) profile (eq. 2 in Kamann et al., 2020) was adopted for the
dispersion. The model parameters were obtained using a maximum-likelihood
approach, where we used the Markov chain Monte Carlo code emcee (Foreman-
Mackey et al., 2013) to sample their posterior distributions. Our model also
included a background component, characterised by a Gaussian velocity
distribution with mean velocity $v_{\rm back}$ and dispersion $\sigma_{\rm
back}$ and a surface stellar density $f_{\rm back}$. Inclusion of the
background component allowed us to determine a posteriori cluster membership
probabilities for all stars in our sample, via

$p_{i}=\left[1+\frac{f_{back}}{1-f_{back}}\sqrt{\frac{\sigma_{back}^{2}+\epsilon_{i}^{2}}{\sigma^{2}+\epsilon_{i}^{2}}}\exp{\left(\frac{(v_{i}-v_{sys})^{2}}{2(\sigma^{2}+\epsilon_{i}^{2})}-\frac{(v_{i}-v_{back})^{2}}{2(\sigma_{back}^{2}+\epsilon_{i}^{2})}\right)}\right]^{-1}\,.$

where $v_{i}$ and $\epsilon_{i}$ are the radial velocity of the $i$-th star
and its associated uncertainty and $\sigma$ is the velocity dispersion of the
bin containing it. We do not try to identify or remove the contribution of
binary stars.

We also analysed the kinematics in radial bins that contained at least 100
stars and covered at least a radial range of 0.2 dex in log radius. From this,
we excluded all stars with membership probabilities $p<0.5$. Further, we fixed
the mean velocity in each bin to the value of $v_{\rm sys}=-106.5\pm 0.1\,{\rm
km\,s^{-1}}$ that was determined in the previous step. Hence, in each bin, we
optimized three parameters: the rotation amplitude, the position angle of the
rotation axis, and the velocity dispersion. To do so, we used the same
maximum-likelihood approach introduced above. The values of the three
parameters in all radial bins are included in Fig. 5.

Second, we used the Cappellari & Copin (2003) code to create Voronoi maps of
the kinematics of M15 with an average number of 100 stars per bin. In each
bin, we used a simple Nelder-Mead optimization in order to obtain the mean
velocity and the velocity dispersion. Stars with membership probabilities
$p<0.5$ were again excluded from the calculation. The resulting Voronoi maps
of the mean velocity and the velocity dispersion are shown in Fig. 6.

Figure 5: Stellar kinematics of M15. The panels show the position angle (top),
rotation velocity (middle) and velocity dispersion (bottom), all as a function
of clustercentric radius. The solid lines show the results of fitting our
unbinned data with simple analytical models, while the green circles show the
results from radial binning of our dataset. The blue dashed and dotted lines
show the core and half light radii respectively from the 2010 edition of the
Harris (1996) catalogue. The MUSE NFM data allows the inner kinematics to be
studied in great detail, revealing the strong twist in rotation in the centre
of M15 and a double peaked rotation curve that rises towards the centre.
Figure 6: Voronoi maps of the line-of-sight heliocentric velocity (left) and
velocity dispersion (right) for the centre of M15. Each bin contains $\sim
100$ stars. The colour scale of the left panel is centred on the cluster
systemic velocity of $-106.5~{}{\rm km\,s^{-1}}$. The solid green circle shows
the core radius while the dashed green line shows the half light radius. A
clear 90 degree twist is seen in the rotation field between the core radius
and the half light radius. However the velocity field is more complicated than
a twisted disc with the negative velocities below and to the left of the
centre not having corresponding features above and to the right of the centre
that would be expected if the field is symmetric. The velocity dispersion map
shows that the peak is extended in a north-south direction.

We see clear evidence for a kinematically distinct core in M15 in both the
radial profiles and the Voronoi bins. While at large radii, the measured
rotation profile is largely consistent with the analytic profile (cf. Fig. 5,
upper left panel), this is not true for radii $r\lesssim 5\arcsec$, where the
rotation measurements obtained via radially binning the data diverge from the
fitted profile. Rather than showing a single peak, the rotation amplitude
shows a minimum at $\sim 20$ arcsec with two peaks at $\sim 2$ and $\sim 100$
arcsec. Since the analytic profile is fit to the individual stars and not the
binned data, the profile is a better fit to the bins at large radii due to the
larger number of stars in this range. Furthermore, the axis of rotation
clearly twists by $\sim 90$ deg at roughly 10 arcsec of the centre (cf. Fig.
5, right panel). This twist has been previously identified by van den Bosch et
al. (2006), Bianchini et al. (2013) and Kamann et al. (2018a). Another twist
may be present even closer to the centre with the inner-most bin showing a
rotation axis $\sim 180$ deg different from the next further out bin.

The Voronoi map shows a velocity field more complex and asymmetric than
expected from a simple twisted disc model. For example, a region with negative
velocity is seen south-west of the centre that does not have a corresponding
region with positive velocity to the north-east. The velocity dispersion map
shows an elongation of the central velocity peak in the north-south direction
such that the velocity dispersion north and south of the centre is higher than
east or west of the centre at the same distance.

##  4 Discussion

Kinematically distinct cores (KDCs) are found in the centres of some massive
early-type galaxies (e.g. 19 out 260 galaxies in the ATLAS${}^{\text{3D}}$
sample, Krajnović et al. 2011). In galaxies, KDCs are usually thought to form
due to mergers (e.g. Kormendy, 1984; Franx & Illingworth, 1988; Balcells &
Quinn, 1990; Hernquist & Barnes, 1991; Hoffman et al., 2010; Bois et al.,
2011). However, mergers of GC are thought to be relatively rare (Khoperskov et
al., 2018). We note that M15 does not show any spread in Fe (Carretta et al.,
2009) which might be expected if M15 is a merger remnant.

Alternatively, a KDC may be due to the projection effects of a triaxial system
with a different distribution of orbits in the centre than in the outskirts
(Statler, 1991; van den Bosch et al., 2008). Core collapse could conceivably
cause the orbits in the centre of the cluster to be different than in the
outskirts.

Kamann et al. (2018a) found a correlation between rotation and relaxation
times at the half light radii, suggesting that GCs loose angular momentum via
two-body interactions as they dynamically evolve (see Tiongco et al., 2017,
for a theoretical modeling of this process). Due to their significantly larger
sizes, in galaxies two-body relaxation timescales are longer than the age of
the Universe so KDCs are expected to persist in them till the present day. Due
to the short two-body relaxation time scale of the core of M15 ($<100$ Myr,
Harris 1996), any ordered rotation should be quickly erased. Thus it is
unlikely that the complex kinematics of centre M15 are a relic of its
formation.

Beccari et al. (2019) used the same HST ACS/HRC imaging used in this work to
study the blue straggler population of M15. They found evidence of multiple
generations of blue stragglers formed by stellar collisions. Beccari et al.
suggest that these generations of blue stragglers formed during brief periods
of increased stellar density caused by core collapse and by the oscillations
that occur after core collapse (e.g. Bettwieser & Sugimoto, 1984; Goodman,
1987; Cohn et al., 1989; Takahashi & Inagaki, 1991; Breen & Heggie, 2012). The
complicated velocity field of the centre M15 may also be due to these post
core collapse oscillations. The post core collapse model of Grabhorn et al.
(1992) reproduces the observed large core radius, cuspy light profile and high
central velocity peak of M15. Grabhorn et al. (1992) also find that the post
core collapse oscillations only produce significant temporal variations in
surface brightness and in kinematics within $\sim 1$ arcsec. This would neatly
explain how the central kinematics of M15 are disconnected from the outer
kinematics. Detailed dynamical modelling will be required to test whether the
observed kinematics can be explained by post core collapse oscillations.

Mapelli et al. (2005) propose that a central binary with massive components
(i.e. two massive stellar-mass BHs $\gtrsim 20\,{\rm M}_{\odot}$ or two IMBHs)
could transfer angular momentum and energy to nearby stars, causing them to
co-rotate with the binary and have high velocity dispersion. These predictions
are very similar to what we see in the inner few arc seconds and would also
naturally explain the misalignment between the inner and outer rotation, but
the collapsed core argues strongly against the presence of BH binaries of any
mass (Hurley, 2007; Heggie et al., 2007; Breen & Heggie, 2013). It could be
that M15 only recently ejected its last BH binary, and still shows features of
the co-rotation imprinted by this binary. For this to work, angular momentum
needs to be retained during core collapse and some time thereafter.
Alternatively, the mechanism proposed by Mapelli et al. (2005) also works with
lower mass binaries (e.g. a recently formed binary with non-BH components,
such as neutrons stars or white dwarfs). Dynamical modelling is required to
test these ideas and asses whether they can explain the observations.

We note that Boberg et al. (2017) discussed the case of another cluster, M53,
which also exhibits strong internal rotation in its central regions with a
complex radial profile and a radial variation of the position angle of the
rotation axis. Their data however do not probe the innermost few arcseconds
that we are able to study in M15 with MUSE, and M53 is not a core collapse
cluster. Kamann et al. (2018a) also found evidence for enhanced rotation in
the centre and twists in the position angle with radius in five of their 25
GCs (NGC 362, M5, M10, M15 and M62). All five of these GCs are relatively
compact (King profile concentrations $c>1.7$) but only three (NGC 362, M15 and
M62) are listed as possible core collapse GCs in the 2010 edition of the
Harris (1996) catalogue.

##  5 Conclusion

Using the new narrow field mode of MUSE we have studied the central kinematics
of the core collapse GC M15 in detail. The combination of AO delivering HST-
like spatial resolution with MUSE’s large field-of-view and analysis
techniques to extract spectra from crowded fields (i.e. PampelMuse) allows for
a significantly larger sample of stars (864) with radial velocities compared
to earlier studies with AO (Gebhardt et al., 2000, 104 stars) and HST (Gerssen
et al., 2002, 64 stars). The improvements in instrumentation and analysis
techniques over the last two decades are highlighted by the fact our
significantly larger sample required a significantly smaller amount of
telescope time (1.3 hours of exposure time) compared to previous studies (6.5
hours for Gebhardt et al. 2000; 25 HST orbits or 16.2 hours of exposure time
for Gerssen et al. 2002).

We see clear evidence for complex central kinematics that are distinct from
those of the outer parts of the cluster. Due to the short two-body relaxation
time scale, this KDC is unlikely to have survived from the formation of M15.
Nor is it likely that this KDC shares a similar merger origin to the much
larger KDCs observed in some massive early-type galaxies. In a future work we
will use dynamical models to fit this kinematic data set. This will allow us
to explore possible hypotheses for the peculiar kinematics of this cluster,
including the transfer of angular momentum from black hole binaries and post
core-collapse oscillations. The dynamical models will further allow us to
place strong constraints on the presence of any compact system of stellar
remnants (e.g. Gieles et al., 2018; Hénault-Brunet et al., 2019; Hénault-
Brunet et al., 2020).

## Acknowledgements

We thank the anonymous referee for their helpful suggestions which improved
the manuscript. CU acknowledges the support of the Swedish Research Council,
Vetenskapsrådet. SK gratefully acknowledges funding from UKRI in the form of a
Future Leaders Fellowship (grant no. MR/T022868/1) and financial support from
the European Research Council (ERC-CoG-646928, Multi-Pop). MG acknowledges
support from the Ministry of Science and Innovation through a Europa
Excelencia grant (EUR2020-112157). VHB acknowledges the support of the Natural
Sciences and Engineering Research Council of Canada (NSERC) through grant
RGPIN-2020-05990. EB acknowledges financial support from a Vici grant by the
Netherlands Organisation for Scientific Research (NWO).

Based on observations collected at the European Organisation for Astronomical
Research in the Southern Hemisphere under ESO programme 60.A-9492(A).

This work made use of numpy (van der Walt et al., 2011), scipy (Jones et al.,
01), matplotlib (Hunter, 2007) and corner (Foreman-Mackey, 2016) as well as
astropy, a community-developed core Python package for astronomy (Astropy
Collaboration et al., 2013).

## Data availability

The raw MUSE data are available from the ESO archive
(https://archive.eso.org/eso/eso_archive_main.html). The reduced data
underlying this work is available from the authors upon reasonable request.

## References

  * Abolfathi et al. (2018) Abolfathi B., et al., 2018, ApJS, 235, 42
  * Arsenault et al. (2012) Arsenault R., et al., 2012, ESO adaptive optics facility progress report. p. 84470J, doi:10.1117/12.926074
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Bacon et al. (2010) Bacon R., et al., 2010, The MUSE second-generation VLT instrument. p. 773508, doi:10.1117/12.856027
  * Balcells & Quinn (1990) Balcells M., Quinn P. J., 1990, ApJ, 361, 381
  * Baumgardt & Hilker (2018) Baumgardt H., Hilker M., 2018, MNRAS, 478, 1520
  * Baumgardt et al. (2003) Baumgardt H., Hut P., Makino J., McMillan S., Portegies Zwart S., 2003, ApJ, 582, L21
  * Beccari et al. (2019) Beccari G., et al., 2019, ApJ, 876, 87
  * Bellini et al. (2014) Bellini A., et al., 2014, ApJ, 797, 115
  * Bettwieser & Sugimoto (1984) Bettwieser E., Sugimoto D., 1984, MNRAS, 208, 493
  * Bianchini et al. (2013) Bianchini P., Varri A. L., Bertin G., Zocchi A., 2013, ApJ, 772, 67
  * Bianchini et al. (2018) Bianchini P., van der Marel R. P., del Pino A., Watkins L. L., Bellini A., Fardal M. A., Libralato M., Sills A., 2018, MNRAS, 481, 2125
  * Boberg et al. (2017) Boberg O. M., Vesperini E., Friel E. D., Tiongco M. A., Varri A. L., 2017, ApJ, 841, 114
  * Bois et al. (2011) Bois M., et al., 2011, MNRAS, 416, 1654
  * Breen & Heggie (2012) Breen P. G., Heggie D. C., 2012, MNRAS, 425, 2493
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
  * Cappellari & Copin (2003) Cappellari M., Copin Y., 2003, MNRAS, 342, 345
  * Carretta et al. (2009) Carretta E., Bragaglia A., Gratton R., D’Orazi V., Lucatello S., 2009, A&A, 508, 695
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Cohn et al. (1989) Cohn H., Hut P., Wise M., 1989, ApJ, 342, 814
  * Dalessandro et al. (2016) Dalessandro E., et al., 2016, ApJ, 833, 111
  * Dalgleish et al. (2020) Dalgleish H., et al., 2020, MNRAS, 492, 3859
  * Djorgovski & King (1986) Djorgovski S., King I. R., 1986, ApJ, 305, L61
  * Dotter et al. (2010) Dotter A., et al., 2010, ApJ, 708, 698
  * Ferraro et al. (2018) Ferraro F. R., et al., 2018, ApJ, 860, 50
  * Foreman-Mackey (2016) Foreman-Mackey D., 2016, The Journal of Open Source Software, 24
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Franx & Illingworth (1988) Franx M., Illingworth G. D., 1988, ApJ, 327, L55
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A12
  * Gebhardt et al. (1994) Gebhardt K., Pryor C., Williams T. B., Hesser J. E., 1994, AJ, 107, 2067
  * Gebhardt et al. (1997) Gebhardt K., Pryor C., Williams T. B., Hesser J. E., Stetson P. B., 1997, AJ, 113, 1026
  * Gebhardt et al. (2000) Gebhardt K., Pryor C., O’Connell R. D., Williams T. B., Hesser J. E., 2000, AJ, 119, 1268
  * Gerssen et al. (2002) Gerssen J., van der Marel R. P., Gebhardt K., Guhathakurta P., Peterson R. C., Pryor C., 2002, AJ, 124, 3270
  * Gieles et al. (2018) Gieles M., Balbinot E., Yaaqib R. I. S. M., Hénault-Brunet V., Zocchi A., Peuten M., Jonker P. G., 2018, MNRAS, 473, 4832
  * Giesers et al. (2019) Giesers B., et al., 2019, A&A, 632, A3
  * Goodman (1987) Goodman J., 1987, ApJ, 313, 576
  * Gott (1973) Gott III R. J., 1973, ApJ, 186, 481
  * Grabhorn et al. (1992) Grabhorn R. P., Cohn H. N., Lugger P. M., Murphy B. W., 1992, ApJ, 392, 86
  * Harris (1996) Harris W. E., 1996, AJ, 112, 1487
  * Heggie et al. (2007) Heggie D. C., Hut P., Mineshige S., Makino J., Baumgardt H., 2007, PASJ, 59, L11
  * Hénault-Brunet et al. (2019) Hénault-Brunet V., Gieles M., Sollima A., Watkins L. L., Zocchi A., Claydon I., Pancino E., Baumgardt H., 2019, MNRAS, 483, 1400
  * Hénault-Brunet et al. (2020) Hénault-Brunet V., Gieles M., Strader J., Peuten M., Balbinot E., Douglas K. E. K., 2020, MNRAS, 491, 113
  * Hernquist & Barnes (1991) Hernquist L., Barnes J. E., 1991, Nature, 354, 210
  * Hoffman et al. (2010) Hoffman L., Cox T. J., Dutta S., Hernquist L., 2010, ApJ, 723, 818
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Hurley (2007) Hurley J. R., 2007, MNRAS, 379, 93
  * Husser et al. (2013) Husser T. O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6
  * Husser et al. (2016) Husser T.-O., et al., 2016, A&A, 588, A148
  * Illingworth & King (1977) Illingworth G., King I. R., 1977, ApJ, 218, L109
  * Ji & Bregman (2015) Ji J., Bregman J. N., 2015, ApJ, 807, 32
  * Jones et al. (01 ) Jones E., Oliphant T., Peterson P., et al., 2001–, SciPy: Open source scientific tools for Python, http://www.scipy.org/
  * Kamann et al. (2013) Kamann S., Wisotzki L., Roth M. M., 2013, A&A, 549, A71
  * Kamann et al. (2014) Kamann S., Wisotzki L., Roth M. M., Gerssen J., Husser T. O., Sandin C., Weilbacher P., 2014, A&A, 566, A58
  * Kamann et al. (2016) Kamann S., et al., 2016, A&A, 588, A149
  * Kamann et al. (2018a) Kamann S., et al., 2018a, MNRAS, 473, 5591
  * Kamann et al. (2018b) Kamann S., et al., 2018b, MNRAS, 480, 1689
  * Kamann et al. (2020) Kamann S., et al., 2020, MNRAS, 492, 966
  * Khoperskov et al. (2018) Khoperskov S., Mastrobuono-Battisti A., Di Matteo P., Haywood M., 2018, A&A, 620, A154
  * Kimmig et al. (2015) Kimmig B., Seth A., Ivans I. I., Strader J., Caldwell N., Anderton T., Gregersen D., 2015, AJ, 149, 53
  * Kormendy (1984) Kormendy J., 1984, ApJ, 287, 577
  * Krajnović et al. (2011) Krajnović D., et al., 2011, MNRAS, 414, 2923
  * Lane et al. (2011) Lane R. R., Kiss L. L., Lewis G. F., Ibata R. A., Siebert A., Bedding T. R., Székely P., Szabó G. M., 2011, A&A, 530, A31
  * Lanzoni et al. (2013) Lanzoni B., et al., 2013, ApJ, 769, 107
  * Lynden-Bell (1967) Lynden-Bell D., 1967, MNRAS, 136, 101
  * Mapelli et al. (2005) Mapelli M., Colpi M., Possenti A., Sigurdsson S., 2005, MNRAS, 364, 1315
  * Massari et al. (2016) Massari D., et al., 2016, A&A, 595, L2
  * McNamara et al. (2003) McNamara B. J., Harrison T. E., Anderson J., 2003, ApJ, 595, 187
  * Milone et al. (2012) Milone A. P., et al., 2012, A&A, 540, A16
  * Murphy et al. (2011) Murphy B. W., Cohn H. N., Lugger P. M., 2011, ApJ, 732, 67
  * Newell et al. (1976) Newell B., Da Costa G. S., Norris J., 1976, ApJ, 208, L55
  * Pasquini et al. (2002) Pasquini L., et al., 2002, The Messenger, 110, 1
  * Peterson et al. (1989) Peterson R. C., Seitzer P., Cudworth K. M., 1989, ApJ, 347, 251
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Rigaut et al. (1998) Rigaut F., et al., 1998, PASP, 110, 152
  * Sarajedini et al. (2007) Sarajedini A., et al., 2007, AJ, 133, 1658
  * Sollima et al. (2019) Sollima A., Baumgardt H., Hilker M., 2019, MNRAS, 485, 1460
  * Statler (1991) Statler T. S., 1991, AJ, 102, 882
  * Ströbele et al. (2012) Ströbele S., et al., 2012, GALACSI system design and analysis. p. 844737, doi:10.1117/12.926110
  * Takahashi & Inagaki (1991) Takahashi K., Inagaki S., 1991, PASJ, 43, 589
  * Tiongco et al. (2017) Tiongco M. A., Vesperini E., Varri A. L., 2017, MNRAS, 469, 683
  * VandenBerg et al. (2013) VandenBerg D. A., Brogaard K., Leaman R., Casagrande L., 2013, ApJ, 775, 134
  * Vasiliev (2019) Vasiliev E., 2019, MNRAS, 489, 623
  * Vogt et al. (1994) Vogt S. S., et al., 1994, HIRES: the high-resolution echelle spectrometer on the Keck 10-m Telescope. p. 362, doi:10.1117/12.176725
  * Watkins et al. (2015) Watkins L. L., van der Marel R. P., Bellini A., Anderson J., 2015, ApJ, 803, 29
  * Weilbacher et al. (2020) Weilbacher P. M., et al., 2020, A&A, 641, A28
  * de Boer et al. (2019) de Boer T. J. L., Gieles M., Balbinot E., Hénault-Brunet V., Sollima A., Watkins L. L., Claydon I., 2019, MNRAS, 485, 4906
  * den Brok et al. (2014) den Brok M., van de Ven G., van den Bosch R., Watkins L., 2014, MNRAS, 438, 487
  * van den Bosch et al. (2006) van den Bosch R., de Zeeuw T., Gebhardt K., Noyola E., van de Ven G., 2006, ApJ, 641, 852
  * van den Bosch et al. (2008) van den Bosch R. C. E., van de Ven G., Verolme E. K., Cappellari M., de Zeeuw P. T., 2008, MNRAS, 385, 647
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science & Engineering, 13, 22


# More on Projected Type Iteration Method and Linear Complementarity Problem

Bharat Kumara,1, Deepmalaa,2 and A.K. Das b,3
_a Mathematics Discipline,_
_PDPM-Indian Institute of Information Technology, Design and Manufacturing,_
_Jabalpur - 482005 (MP), India_
_b Indian Statistical Institute, 203 B.T. Road, _
_Kolkata - 700108, India_
_1<EMAIL_ADDRESS>, 2Email<EMAIL_ADDRESS>
_3 Email<EMAIL_ADDRESS>
###### Abstract

In this article, we establish a class of new projected type iteration methods
based on matrix spitting for solving the linear complementarity problem. Also,
we provide a sufficient condition for the convergence analysis when the system
matrix is an $H_{+}$-matrix. We show the efficiency of the proposed method by
using two numerical examples for different parameters.

Keywords. Iterative method, Linear complementarity problem, $H_{+}$-matrix,
$P$-matrix, Matrix splitting, Convergence.

Mathematics Subject Classification. 90C33, 65F10, 65F50.

##  1 Introduction

The LCP frequently appears in an extensive range of applications that include
scientific computing and engineering, such as the free boundary problem and
the Nash equilibrium point of the bimatrix game; the American option pricing
problem; mathematical economics; operations research; control theory;
optimization theory; stochastic optimal control; economics; and elasticity
theory. For details see [5], [27], [24] [16], [10], [26], [7], [15], [18] and
[21].
Assuming $\mathcal{A}\in\mathcal{R}^{n\times n}$ and a vector
$\,\sigma\,\in\,\mathcal{R}^{n}.$ The linear complementarity problem denoted
as LCP$(\sigma,\mathcal{A})$ is to find the solution
$\lambda\in\mathcal{R}^{n}$ to the following system

$\displaystyle\lambda\geq 0,~{}~{}~{}~{}\mathcal{A}\lambda+\sigma\geq
0,~{}~{}~{}~{}\lambda^{T}(\mathcal{A}\lambda+\sigma)=0$ (1)

The methods for solving linear complementarity problems are divided into two
categories: the pivoting method [6] [8], [14] and iterative method [25], [13],
[17], [20] and [19]. Lemke and Howson [22] introduced the complementary pivot
method, but some matrices are not processable by this method as well as by
Lemke’s Method. The linear complementarity problem can be solved in a number
of ways by an iterative process; namely, the projected type methods [3], [13],
[25], the modulus method [2], [9] and the modulus based matrix splitting
iterative methods [23] and [28].
Fang proposed a general fixed point method (GFP) [11] assuming the case where
$\Omega=\omega A_{\mathcal{D}}^{-1}$ with $\omega\textgreater 0$ and
$A_{\mathcal{D}}$ is the diagonal matrix of $\mathcal{A}$. The GFP approach
takes less iterations than the modulus-based successive over-relaxation (MSOR)
[2] iteration method. However, the GFP approach calculates the numerical
solution component by component of vectors, which takes a long time.
In this article, we present a class of new projected type iteration methods by
using the ideas of Xi [11] and Ali [1]. Also, we show that the fixed point
equation and the linear complementarity problem are equivalent, discuss
convergence conditions and provide a convergence domain for our proposed
method.
The article is organized as follows: some required definitions, notations and
well-known lemmas are given in Section 2, which will be used for the
discussions in the remaining sections of this work. New projected type
iteration methods are constructed in Section 3 with the help of the new
equivalent fixed point form of the LCP$(\sigma,\mathcal{A})$. In Section 4, we
establish the convergence domain of our proposed method. A numerical
comparison between the proposed methods and modulus-based matrix splitting
iteration methods, introduced by Bai [2], is illustrated in Section 5. Section
6 contains the conclusion of the article.

##  2 Preliminaries

In this section, we provide an overview of various essential notations,
definitions, and foundational results.
Suppose $\mathcal{A}=({a}_{ij})\in\mathcal{R}^{n\times n}$ and
$\mathcal{B}=({b}_{ij})\in\mathcal{R}^{n\times n}$ are real square matrices.
For $\mathcal{A}=({a}_{ij})\in\mathcal{R}^{n\times n}$ and
$\mathcal{B}=({b}_{ij})\in\mathcal{R}^{n\times n}$, $\mathcal{A}\geq$
$(\textgreater)$ $\mathcal{B}$ means ${a}_{ij}\geq(\textgreater)$ ${b}_{ij}$
for all $i,j$.

######  Definition 2.1.

[11] Let $\mathcal{A}=({a}_{ij})\in\mathcal{R}^{n\times n}$. Then
$|\mathcal{A}|=({c}_{ij})$ is defined by ${c}_{ij}=|{a}_{ij}|$ $\forall~{}i,j$
and $|\mathcal{A}|$ represent that ${a}_{ij}\geq 0$ $\forall~{}i,j$.

######  Definition 2.2.

[11] Let $\mathcal{A},\mathcal{B}\in\mathcal{R}^{n\times n}$. Then
$|\mathcal{A}+\mathcal{B}|\leq|\mathcal{A}|+|\mathcal{B}|$ and
$|\mathcal{A}\mathcal{B}|\leq|\mathcal{A}||\mathcal{B}|$. Moreover
${x},{y}\in\mathcal{R}^{n}$ then $|{x}+{y}|\leq|{x}|+|{y}|$ and
$||{x}|-|{y}||\leq|{x}-{y}|$.

######  Definition 2.3.

[8] Let $\mathcal{A}\in\mathcal{R}^{n\times n}$. $\mathcal{A}$ is said to be a
$P$-matrix if all its principle minors are positive i.e.
det$({\mathcal{A}}_{\gamma\gamma})~{}\textgreater~{}0$ for all
$\gamma\subseteq\\{1,2,\ldots,n\\}$.

######  Definition 2.4.

[11] Suppose $\mathcal{A}\in\mathcal{R}^{n\times n}$. Then its comparison
matrix is defined as $\langle{a}_{ij}\rangle=|{a}_{ij}|$ if $i=j$ and
$\langle{a}_{ij}\rangle=-|{a}_{ij}|$ if $i\neq j$.

######  Definition 2.5.

[12] Suppose $\mathcal{A}\in\mathcal{R}^{n\times n}$. $\mathcal{A}$ is said to
be a $Z$-matrix if all of its non-diagonal elements are less than or equal to
zero; $\mathcal{A}$ is said to be an $M$-matrix if $\mathcal{A}^{-1}\geq 0$ as
well as $Z$-matrix; $\mathcal{A}$ is said to be an $H$-matrix if
$\langle\mathcal{A}\rangle$ is an $M$-matrix; $\mathcal{A}$ is an
$H_{+}$-matrix if it is an $H$-matrix with
${a}_{ii}~{}\textgreater~{}0~{}\forall~{}i\in\\{1,2,\ldots,n\\}$.

######  Definition 2.6.

[12] Suppose $\mathcal{A}\in\mathcal{R}^{n\times n}$. The splitting
$\mathcal{A}=\mathcal{M}-\mathcal{N}$ is called an $M$-splitting if
$\mathcal{M}$ is a nonsingular $M$-matrix and $\mathcal{N}\geq 0$; an
$H$-splitting if $\langle\mathcal{M}\rangle-|\mathcal{N}|$ is an $M$-matrix;
an $H$-compatible splitting if
$\langle\mathcal{A}\rangle=\langle\mathcal{M}\rangle-|\mathcal{N}|$.

######  Lemma 2.1.

[1] Let ${x},y\in\mathcal{R}^{n}$. ${x}\geq 0$, ${y}\geq 0$, $x^{T}{y}=0$ if
and only if ${x}+{y}=|{x}-{y}|$.

######  Lemma 2.2.

[12] Suppose $\mathcal{A},B_{1}\in\mathcal{R}^{n\times n}$. If $\mathcal{A}$
and $\mathcal{B}$ are $M$ and $Z$-matrices respectively with
$\mathcal{A}\leq\mathcal{B}$ then $\mathcal{B}$ is an $M$-matrix. If
$\mathcal{A}$ is an $H$-matrix then
$|\mathcal{A}^{-1}|\leq\langle\mathcal{A}\rangle^{-1}$. If
$\mathcal{A}\leq\mathcal{B}$, then $\rho(\mathcal{A})\leq\rho(\mathcal{B})$.

######  Lemma 2.3.

[11] Let $\mathcal{A}\in\mathcal{R}^{n\times n}$ be an $M$-matrix and
$\mathcal{A}=\mathcal{M}-\mathcal{N}$ be an $M$-splitting. Let $\rho$ be the
spectral radius, then $\rho(\mathcal{M}^{-1}\mathcal{N})$ $\textless$ $1$.

######  Lemma 2.4.

[4] Suppose $\mathcal{A}\in\mathcal{R}^{n\times n}$ with splitting
$\mathcal{A}=\mathcal{M}-\mathcal{N}$. Let splitting be an $H$-compatible of
an $H$-matrix, then it is an $H$-splitting but converse is not true.

######  Lemma 2.5.

[12] Suppose $\mathcal{A}\geq 0\in\mathcal{R}^{n\times n}$, if there exists
$v~{}\textgreater~{}0\in\mathcal{R}^{n}$ and a scalar
$\alpha_{1}~{}\textgreater~{}0$ such that $\mathcal{A}v\leq\alpha_{1}v$ then
$\rho(\mathcal{A})\leq\alpha_{1}$. Moreover, if $\mathcal{A}v~{}\textless~{}v$
then $\rho(\mathcal{A})\textless 1$.

##  3 Main results

For a given vector $\zeta\in\mathcal{R}^{n}$, we indicate the vector
$\zeta_{+}=$max$\\{0,\zeta\\}$ and matrix
$\mathcal{A}=(\mathcal{M}+I+D_{\mathcal{A}})-(\mathcal{N}+I+D_{\mathcal{A}})$,
where $D_{\mathcal{A}}$ is diagonal matrix of $\mathcal{A}$. In the following
result, we convert the LCP $(\sigma,\mathcal{A})$ into a fixed point
formulation.

######  Theorem 3.1.

Let $\mathcal{A}\in\mathcal{R}^{n\times n}$ with the splitting
$\mathcal{A}=(\mathcal{M}+I+D_{\mathcal{A}})-(\mathcal{N}+I+D_{\mathcal{A}})$.
Let $\lambda=\zeta_{+}$, then equivalent formulation of the
LCP$(\sigma,\mathcal{A})$ in form of fixed point equation is

$\displaystyle\zeta_{+}=(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}[(\mathcal{N}+I+D_{\mathcal{A}})\zeta_{+}+|(\mathcal{A}-I)\zeta_{+}+\sigma|-\sigma]$
(2)

###### Proof.

We have $\lambda=\zeta_{+}\geq 0$ and $\mathcal{A}\lambda+\sigma\geq 0$, from
Lemma 2.1,

$\begin{split}(\mathcal{A}\zeta_{+}+\sigma+\zeta_{+})&=|\mathcal{A}\zeta_{+}+\sigma-\zeta_{+}|\\\
(I+\mathcal{A})\zeta_{+}&=|(\mathcal{A}-I)\zeta_{+}+\sigma|-\sigma\\\
(\mathcal{M}+2I+D_{\mathcal{A}})\zeta_{+}&=(\mathcal{N}+I+D_{\mathcal{A}})\zeta_{+}+|(\mathcal{A}-I)\zeta_{+}+\sigma|-\sigma,\\\
\end{split}$

the above equation can be rewritten as,

$\zeta_{+}=(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}[(\mathcal{N}+I+D_{\mathcal{A}})\zeta_{+}+|(\mathcal{A}-I)\zeta_{+}+\sigma|-\sigma]$
(3)

∎

In the following, Based on Equation (2), we propose an iteration method which
is known as Method 3.1 to solve the LCP$(\sigma,\mathcal{A})$.

######  Method 3.1.

Let
$\mathcal{A}=(\mathcal{M}+I+D_{\mathcal{A}})-(\mathcal{N}+I+D_{\mathcal{A}})$
be a splitting of the matrix $\mathcal{A}\in\mathcal{R}^{n\times n}$ and the
matrix $(\mathcal{M}+2I+D_{\mathcal{A}})$ be the nonsingular. Then we use the
following equation for Method 3.1 is

$\zeta^{(\eta+1)}_{+}=(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}[(\mathcal{N}+I+D_{\mathcal{A}})\zeta^{(\eta)}_{+}+|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|-\sigma]$
(4)

Let Residual be the Euclidean norm of the error vector, which is defined as
follows:

$Res(\lambda^{(\eta)})=|min(\lambda^{(\eta)},\mathcal{A}\lambda^{(\eta)}+\sigma)|_{2}.$

Consider a nonnegative initial vector $\lambda^{(0)}\in\mathcal{R}^{n}$. The
iteration process continues until the iteration sequence
$\\{\lambda^{(\eta)}\\}_{\eta=0}^{+\infty}\subset\mathcal{R}^{n}$ converges.
For $\eta=0,1,2,\ldots$, the iterative process continues until the iterative
sequence $\lambda^{(\eta+1)}\in\mathcal{R}^{n}$ converges. The iteration
process stops if $Res(\lambda^{(\eta)})$ $\textless$ $\epsilon$. For computing
$\lambda^{(\eta+1)}$ we use the following steps.
Step 1: Given an initial vector $\zeta^{(0)}\in\mathcal{R}^{n}$,
$\epsilon~{}\textgreater~{}0$ and set $\eta=0$.
Step 2: Using the following scheme, generate the sequence $\lambda^{(\eta)}$:

$\zeta^{(\eta+1)}_{+}=(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}[(\mathcal{N}+I+D_{\mathcal{A}})\zeta^{(\eta)}_{+}+|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|-\sigma],$

and set $\lambda^{(\eta+1)}=\zeta_{+}^{(\eta+1)}$, where
$\zeta_{+}^{(\eta+1)}$ is the $(\eta+1)^{th}$ approximate solution of Equation
(3) .
Step 3: If $Res(\lambda^{(\eta)})$ $\textless$ $\epsilon$ then stop;
otherwise, set $\eta=\eta+1$ and return to step 2.

Moreover, Method 3.1 provides a general structure for solving
LCP$(\sigma,\mathcal{A})$. We obtain a class of new projected type iteration
relaxation methods using matrix splitting. We express the system matrix
$\mathcal{A}=(\mathcal{M}+I+D_{\mathcal{A}})-(\mathcal{N}+I+D_{\mathcal{A}})$.
Then

  1. 1.

when $\mathcal{M}=D_{\mathcal{A}}-L_{\mathcal{A}}$ and
$\mathcal{N}=U_{\mathcal{A}}$, Equation (4) gives the new projected type Gauss
Seidel iteration (NPGS) method

$\displaystyle\begin{split}\zeta^{(\eta+1)}_{+}&=(D_{\mathcal{A}}-L_{\mathcal{A}}+2I+D_{\mathcal{A}})^{-1}[(U_{\mathcal{A}}+I+D_{\mathcal{A}})\zeta^{(\eta)}_{+}\\\
&+|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|-\sigma].\end{split}$

  2. 2.

when $\mathcal{M}=(\frac{1}{\alpha_{1}}D_{\mathcal{A}}-L_{\mathcal{A}})$ and
$\mathcal{N}=(\frac{1}{\alpha_{1}}-1)D_{\mathcal{A}}+U_{\mathcal{A}}$,
Equation (4) gives the new projected type successive overrelaxation iteration
(NPSOR) method

$\displaystyle\begin{split}\zeta^{(\eta+1)}_{+}&=(D_{\mathcal{A}}-\alpha_{1}L_{\mathcal{A}}+\alpha_{1}(2I+D_{\mathcal{A}}))^{-1}[((1-\alpha_{1})D_{\mathcal{A}}+U_{\mathcal{A}}\\\
&+\alpha_{1}I+D_{\mathcal{A}})\zeta^{(\eta)}_{+}+\alpha_{1}|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|-\alpha_{1}\sigma].\end{split}$

  3. 3.

when
$\mathcal{M}=(\frac{1}{\alpha_{1}})(D_{\mathcal{A}}-\beta_{1}L_{\mathcal{A}})$
and
$\mathcal{N}=(\frac{1}{\alpha_{1}})[(1-\alpha_{1})D_{\mathcal{A}}+(\alpha_{1}-\beta_{1})L_{\mathcal{A}}+\alpha_{1}U_{\mathcal{A}}]$,
Equation (4) gives the new projected type accelerated overrelaxation iteration
(NPAOR) method

$\displaystyle\begin{split}\zeta^{(\eta+1)}_{+}&=(D_{\mathcal{A}}-\beta_{1}L_{\mathcal{A}}+\alpha_{1}(2I+D_{\mathcal{A}}))^{-1}[((1-\alpha_{1})D_{\mathcal{A}}+(\alpha_{1}-\beta_{1})U_{\mathcal{A}}\\\
&+\alpha_{1}I+D_{\mathcal{A}})\zeta^{(\eta)}_{+}+\alpha_{1}|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|-\alpha_{1}\sigma].\end{split}$

When $(\alpha_{1},\beta_{1})$ takes the values $(\alpha_{1},\alpha_{1})$,
$(1,1)$ and $(1,0)$, the NPAOR method transforms into the new projected type
successive overrelaxation (NPSOR), new projected type Gauss-Seidel (NPGS) and
new projected type Jacobi (NPJ) methods respectively.

##  4 Convergence analysis

In the following, we present the convergence condition when the system matrix
$\mathcal{A}$ of LCP$(\sigma,\mathcal{A})$ is a $P$-matrix.

######  Theorem 4.1.

Let $\mathcal{A}\in\mathcal{R}^{n\times n}$ be a $P$-matrix and
$\zeta_{+}^{*}$ be the solution of Equation $(\ref{eq2})$. Let
$\rho(|(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}|(|\mathcal{N}+I+D_{\mathcal{A}}|+|\mathcal{A}-I|))~{}\textless~{}1$.
Then the sequence $\\{\zeta_{+}^{(\eta)}\\}^{+\infty}_{\eta=1}$ generated by
Method 3.1 converges to the solution $\zeta_{+}^{*}$ for any initial vector
$\zeta^{(0)}\in\mathcal{R}^{n}$.

###### Proof.

Let $\zeta^{*}_{+}$ be the solution of Equation $(\ref{eq2})$, then error is

$\begin{split}(\mathcal{M}+2I+D_{\mathcal{A}})(\zeta^{(\eta+1)}_{+}-\zeta^{*}_{+})&=(\mathcal{N}+I+D_{\mathcal{A}})(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})+|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|\\\
&-|(\mathcal{A}-I)\zeta^{*}_{+}+\sigma|\\\
|(\mathcal{M}+2I+D_{\mathcal{A}})(\zeta^{(\eta+1)}_{+}-\zeta^{*}_{+})|&=|(\mathcal{N}+I+D_{\mathcal{A}})(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})+|(\mathcal{A}-I)\zeta^{(\eta)}_{+}+\sigma|\\\
&-|(\mathcal{A}-I)\zeta^{*}_{+}+\sigma||\\\ \end{split}$
$\begin{split}&\leq|(\mathcal{N}+I+D_{\mathcal{A}})(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})+|(\mathcal{A}-I)(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})||\\\
&\leq|(\mathcal{N}+I+D_{\mathcal{A}})||(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})|+|(\mathcal{A}-I)||(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})|\\\
|(\zeta^{(\eta+1)}_{+}-\zeta^{*}_{+})|&\leq|(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}|[|(\mathcal{N}+I+D_{\mathcal{A}})|+|(\mathcal{A}-I)|]|(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})|\\\
|(\zeta^{(\eta+1)}_{+}-\zeta^{*}_{+})|&\textless|(\zeta^{(\eta)}_{+}-\zeta^{*}_{+})|.\end{split}$

Hence $\zeta^{(\eta)}_{+}$ converges to the solution $\zeta^{*}_{+}$. ∎

Now we discuss the convergence conditions for Method 3.1 when the system
matrix $\mathcal{A}$ of LCP$(\sigma,\mathcal{A})$ is an $H_{+}$-matrix.

######  Theorem 4.2.

Let $\mathcal{A}\in\mathcal{R}^{n\times n}$ be an $H_{+}$-matrix and
$\mathcal{A}=\mathcal{M}-\mathcal{N}=(\mathcal{M}+I+D_{\mathcal{A}})-(\mathcal{N}+I+D_{\mathcal{A}})$
be an $H$-compatible splitting of the matrix $\mathcal{A}$, such that
$\langle\mathcal{A}\rangle=\langle\mathcal{M}\rangle-|\mathcal{N}|$=
$\langle\mathcal{M}+I+D_{\mathcal{A}}\rangle-|\mathcal{N}+I+D_{\mathcal{A}}|$
and either one of the following conditions hold:
(1) $D_{\mathcal{A}}~{}\geq~{}I$ and
$\langle\mathcal{A}\rangle+2I-D_{\mathcal{A}}-|B|$ is an $M$\- matrix, where
$B=L_{\mathcal{A}}+U_{\mathcal{A}}$;
(2) $D_{\mathcal{A}}~{}\textless~{}I$.
Then the sequence $\\{\zeta_{+}^{(\eta)}\\}^{+\infty}_{\eta=1}$ generated by
Method 3.1 converges to the solution $\zeta_{+}^{*}$ for any initial vector
$\zeta^{(0)}\in\mathcal{R}^{n}$.

###### Proof.

Let
$\mathcal{A}=\mathcal{M}-\mathcal{N}=(\mathcal{M}+I+D_{\mathcal{A}})-(\mathcal{N}+I+D_{\mathcal{A}})$
and it holds that
$\langle\mathcal{A}\rangle\leq\langle\mathcal{M}+I+D_{\mathcal{A}}\rangle\leq
diag(\mathcal{M}+I+D_{\mathcal{A}})$, hence $(\mathcal{M}+I+D_{\mathcal{A}})$
is an $H_{+}$-matrix and it holds that

$|(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}|\leq(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}.$

Let
$T=|(\mathcal{M}+2I+D_{\mathcal{A}})^{-1}|(|(\mathcal{N}+I+D_{\mathcal{A}})|+|(\mathcal{A}-I)|)$.
Then

$\begin{split}T&=|(2I+\mathcal{M}+D_{\mathcal{A}})^{-1}|[|\mathcal{N}+I+D_{\mathcal{A}}|+|(\mathcal{A}-I)|]\\\
&\leq(2I+\langle\mathcal{M}\rangle+D_{\mathcal{A}})^{-1}[|\mathcal{N}+I+D_{\mathcal{A}}|+|(\mathcal{A}-I)|]\\\
&\leq(2I+\langle\mathcal{M}\rangle+D_{\mathcal{A}})^{-1}[|\mathcal{N}+I+D_{\mathcal{A}}|+|(D_{\mathcal{A}}-I)-(L_{\mathcal{A}}+U_{\mathcal{A}})|]\\\
&\leq(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})\\\
&+|\mathcal{N}+I+D_{\mathcal{A}}|+|D_{\mathcal{A}}-I|+|L_{\mathcal{A}}+U_{\mathcal{A}}|].\\\
\end{split}$

Case 1. Suppose $D_{\mathcal{A}}~{}\geq~{}I$ and
$\langle\mathcal{A}\rangle+2I-D_{\mathcal{A}}-|B|$ is an $M$-matrix then

$\begin{split}T&\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})-|\mathcal{N}+I+D_{\mathcal{A}}|\\\
&-D_{\mathcal{A}}+I-|L_{\mathcal{A}}+U_{\mathcal{A}}|]\\\ &\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[I+(\langle\mathcal{M}\rangle+I+D_{\mathcal{A}})-|\mathcal{N}+I+D_{\mathcal{A}}|\\\
&-D_{\mathcal{A}}+I-|L_{\mathcal{A}}+U_{\mathcal{A}}|]\\\ &\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[(\langle\mathcal{M}\rangle+I+D_{\mathcal{A}})-|\mathcal{N}+I+D_{\mathcal{A}}|+2I-D_{\mathcal{A}}\\\
&-|L_{\mathcal{A}}+U_{\mathcal{A}}|]\\\ &\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}(\langle\mathcal{A}\rangle+2I-D_{\mathcal{A}}-|B|).\\\
\end{split}$

Since $\langle\mathcal{A}\rangle+2I-D_{\mathcal{A}}-|B|$ is an $M$-matrix,
there exists a positive vector $v~{}\textgreater~{}0$ such that

$(\langle\mathcal{A}\rangle+2I-D_{\mathcal{A}}-|B|)v~{}\textgreater~{}0.$

Therefore,

$Tv\leq(I-2(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}(\langle\mathcal{A}\rangle+2I-D_{\mathcal{A}}-|B|)v~{}\textless~{}v$
$Tv~{}\textless~{}v.$

Case 2. Suppose $D_{\mathcal{A}}\textless I$ then

$\begin{split}T&\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})-|\mathcal{N}+I+D_{\mathcal{A}}|\\\
&+D_{\mathcal{A}}-I-|L_{\mathcal{A}}+U_{\mathcal{A}}|]\\\ &\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[(\langle\mathcal{M}\rangle+I+D_{\mathcal{A}})-|\mathcal{N}+I+D_{\mathcal{A}}|\\\
&+D_{\mathcal{A}}-|L_{\mathcal{A}}+U_{\mathcal{A}}|]\\\ &\leq
I-(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}[(\langle\mathcal{M}\rangle+I+D_{\mathcal{A}})-|\mathcal{N}+I+D_{\mathcal{A}}|+D_{\mathcal{A}}\\\
&-|L_{\mathcal{A}}+U_{\mathcal{A}}|]\\\ &\leq
I-2(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}\langle\mathcal{A}\rangle.\\\
\end{split}$

Since $\langle\mathcal{A}\rangle$ is an $M$-matrix, then there exists a
positive vector $v~{}\textgreater~{}0$ such that

$(\langle\mathcal{A}\rangle)v~{}\textgreater~{}0.$

Therefore,

$Tv\leq(I-2(\langle\mathcal{M}\rangle+2I+D_{\mathcal{A}})^{-1}(\langle\mathcal{A}\rangle)v~{}\textless~{}v$

This implies that

$Tv~{}\textless~{}v.$

From Cases 1 and 2 and based on Lemma 2.5, we obtain that $\rho(T)\textless
1$. Therefore based on Theorem 4.1, the iteration sequence
$\\{\zeta_{+}^{(\eta)}\\}^{+\infty}_{\eta=1}$ generated by Method 3.1
converges to $\zeta^{*}_{+}$ for any initial vector $\zeta^{(0)}.$ ∎

##  5 Numerical examples

In this section, two numerical examples are given in this part to demonstrate
the effectiveness of our proposed method and use some notation as number of
iteration steps (denoted by IT), CPU time in seconds (denoted by CPU). Let
$\zeta^{(0)}=(1,0,\ldots 1,0,\ldots)^{T}\in\mathcal{R}^{n}$ be an initial
vector and set $\epsilon=10^{-5}$. We consider the LCP$(\sigma,\mathcal{A})$
which has always unique solution and define
$\sigma=-\mathcal{A}{\lambda}^{*}$, where
${\lambda}^{*}=(1,2,1,\cdots,1,2)^{T}\in\mathcal{R}^{n}.$ The proposed new
projected Gauss Seidel iteration method (NPGS) and the new projected
successive over relaxation iteration method (NPSOR) are compared with the
modulus based Gauss Seidel (MGS) Method and the modulus based successive over
relaxation (MSOR) method [2] respectively, which are effective in solving
LCP$(\sigma,\mathcal{A})$ and set $\Omega=\frac{1}{2\alpha}D_{\mathcal{A}}$
for MGS and MSOR methods. Matlab version 2021a on an Acer Desktop (Intel(R)
Core(TM) i7-8700 CPU @ 3.2 GHz, 3.19 GHz, 16.00GB RAM) is used for all
calculations. Table 1 and Table 2 list the numerical results for new projected
type iteration matrix splitting Method 3.1 (NPGS, NPSOR) and modulus based
matrix splitting Methods (MGS, MSOR).

######  Example 5.1.

The system matrix $\mathcal{A}$ are generated by
$\mathcal{A}=P_{1}+\delta_{1}I$, where $\delta_{1}$ are nonnegative real
parameter and

$P_{1}=\begin{bmatrix}L_{1}&-I_{1}&0&\ldots&0\\\
-I_{1}&L_{1}&-I_{1}&\ldots&0\\\ 0&-I_{1}&L_{1}&-I_{1}&0\\\
0&\ldots&I_{1}&\ddots&-I_{1}\\\ 0&\ldots&0&-I_{1}&L_{1}\\\
\end{bmatrix}$$\in\mathcal{R}^{n\times n}$,
$L_{1}=\begin{bmatrix}4&-1&\ldots&\ldots&0\\\ -1&4&-1&\ldots&0\\\
0&-1&4&-1&0\\\ 0&\ldots&-1&\ddots&-1\\\ 0&\ldots&\ldots&-1&4\\\ \end{bmatrix}$

$\in\mathcal{R}^{m\times m}$, where $I_{1}$ is the identity matrix of order
$m$.

Table 1: Results for MGS and MSOR methods and NPGS and NPSOR methods, when $\delta_{1}=4$. | n | ${100}$ | ${900}$ | ${2500}$ | ${3600}$ | ${6400}$ | ${10000}$
---|---|---|---|---|---|---|---
MGS | IT | 36 | 40 | 41 | 41 | 42 | 42
$\alpha=1$ | CPU | 0.0030 | 0.0254 | 0.2550 | 0.6083 | 1.8468 | 2.7943
| Res | 9.7e-06 | 8.0e-06 | 7.9e-06 | 8.9e-06 | 7.4e-06 | 8.4e-06
NPGS | IT | 21 | 23 | 254 | 24 | 25 | 25
$\alpha_{1}=1$ | CPU | 0.0021 | 0.0035 | 0.0175 | 0.0636 | 0.1725 | 0.3785
| Res | 5.2e-06 | 7.1e-06 | 6.5e-06 | 8.0e-06 | 5.5e-06 | 7.0e-06
MSOR | IT | 15 | 17 | 18 | 18 | 18 | 19
$\alpha=0.85$ | CPU | 0.0024 | 0.0044 | 0.0134 | 0.0462 | 0.1118 | 0.2246
| Res | 9.5e-06 | 7.6e-06 | 5.2e-06 | 6.5e-06 | 8.9e-06 | 4.3e-06
NPSOR | IT | 15 | 16 | 17 | 17 | 17 | 17
$\alpha_{1}=1.7$ | CPU | 0.0019 | 0.0031 | 0.0138 | 0.0449 | 0.1108 | 0.2217
| Res | 5.9e-06 | 6.8e-06 | 4.2e-06 | 4.9e-06 | 6.3e-06 | 7.7e-06

######  Example 5.2.

The system matrix $\mathcal{A}\in\mathcal{R}^{n\times n}$ is generated by
$\mathcal{A}=P_{1}+\delta_{1}I$, where $\delta_{1}$ are nonnegative real
parameter and

$P_{1}=\begin{bmatrix}L_{1}&-0.5I_{1}&0&\ldots&0\\\
-1.5I_{1}&L_{1}&-0.5I_{1}&\ldots&0\\\ 0&-1.5I_{1}&L_{1}&-0.5I_{1}&0\\\
0&\ldots&-1.5I_{1}&\ddots&-0.5I_{1}\\\
0&\ldots&0&-1.5I_{1}&L_{1}\end{bmatrix}$,
$L_{1}=\begin{bmatrix}4&-1&\ldots&\ldots&0\\\ -1&4&-1&\ldots&0\\\
0&-1&4&-1&0\\\ 0&\ldots&-1&\ddots&-1\\\ 0&\ldots&\ldots&-1&4\\\ \end{bmatrix}$

$\in\mathcal{R}^{m\times m}$, where $I_{1}$ is the identity matrix of order
$m$.

Table 2: Results for MGS and MSOR methods and NPGS and NPSOR methods, when $\delta_{1}=4$. | n | ${100}$ | ${400}$ | ${900}$ | ${1600}$ | ${2500}$ | ${3600}$
---|---|---|---|---|---|---|---
MGS | IT | 24 | 26 | 26 | 26 | 27 | 27
$\alpha=1$ | CPU | 0.0038 | 0.0043 | 0.0203 | 0.0691 | 0.1978 | 0.3915
| Res | 9.0e-06 | 6.5e-06 | 8.7e-06 | 9.6e-06 | 6.6e-06 | 7.4e-06
NPGS | IT | 18 | 21 | 22 | 22 | 22 | 23
$\alpha_{1}=1$ | CPU | 0.0019 | 0.0033 | 0.0138 | 0.0469 | 0.1260 | 0.2637
| Res | 9.8e-06 | 5.5e-06 | 4.8e-06 | 5.9e-06 | 8.0e-06 | 4.9e-06
MSOR | IT | 14 | 14 | 15 | 15 | 15 | 15
$\alpha=0.88$ | CPU | 0.0025 | 0.0037 | 0.0121 | 0.0385 | 0.0932 | 0.1823
| Res | 3.8e-06 | 8.2e-06 | 4.0e-06 | 4.6e-06 | 5.50e-06 | 6.3e-06
NPSOR | IT | 12 | 13 | 14 | 14 | 14 | 14
$\alpha_{1}=1.7$ | CPU | 0.0017 | 0.0030 | 0.0112 | 0.0356 | 0.0921 | 0.1865
| Res | 5.5e-06 | 6.6e-06 | 3.6e-06 | 4.4e-06 | 5.9e-06 | 7.5e-06

From Table 1 and Table 2, we can observe that the number of iteration steps
required for our proposed NPGS and NPSOR methods is less than the MGS and MSOR
methods.

##  6 Conclusion

In this article, we introduce a class of new projected-type iteration methods
based on matrix splitting for solving the linear complementarity problem LCP
$(\sigma,\mathcal{A})$. During the iteration process, the large and sparse
structure of $\mathcal{A}$ is maintained by these iterative forms. Moreover,
the sufficient conditions for convergence for $H_{+}$ matrix or $P$-matrix are
presented. Finally, two numerical examples are provided to demonstrate the
effectiveness of the proposed methods.
Conflict of interest The authors declare that there is no conflicts of
interest.
Acknowledgment. The first author is thankful to the University Grants
Commission (UGC), Government of India, under the JRF fellowship programme no.
1068/(CSIR-UGC NET DEC. 2017).

## References

  * [1] R Ali, I Khan, A Ali, and A Mohamed.  Two new generalized iteration methods for solving absolute value equations using m-matrix.  AIMS mathematics, 7(5):8176–8187, 2022.
  * [2] ZZ Bai.  Modulus-based matrix splitting iteration methods for linear complementarity problems.  Numerical Linear Algebra with Applications, 17(6):917–933, 2010\.
  * [3] ZZ Bai and DJ Evans.  Matrix multisplitting relaxation methods for linear complementarity problems.  International Journal of Computer Mathematics, 63(3-4):309–326, 1997.
  * [4] A Berman and RJ Plemmons.  Nonnegative matrices in the mathematical sciences.  SIAM, 1994.
  * [5] RW Cottle, JS Pang, and RE Stone.  The linear complementarity problem. academic pr ess.  Inc., Boston, MA, 5, 1992.
  * [6] AK Das.  Properties of some matrix classes based on principal pivot transform.  Annals of Operations Research, 243:375–382, 2016.
  * [7] AK Das, Deepmala, and R Jana.  Some aspects on solving transportation problem.  Yugoslav Journal of Operations Research, 30(1):45–57, 2020.
  * [8] AK Das and R Jana.  On generalized positive subdefinite matrices and interior point algorithm.  In Operations Research and Optimization: FOTA 2016, Kolkata, India, November 24-26 1, pages 3–16. Springer, 2018.
  * [9] JL Dong and MQ Jiang.  A modified modulus method for symmetric positive-definite linear complementarity problems.  Numerical Linear Algebra with Applications, 16(2):129–143, 2009\.
  * [10] A Dutta, R Jana, and AK Das.  On column competent matrices and linear complementarity problem.  In Proceedings of the Seventh International Conference on Mathematics and Computing: ICMC 2021, pages 615–625. Springer, 2022.
  * [11] XM Fang.  General fixed-point method for solving the linear complementarity problem.  AIMS Mathematics, 6(11):11904–11920, 2021.
  * [12] A Frommer and DB Szyld.  H-splittings and two-stage iterative methods.  Numerische Mathematik, 63:345–356, 1992.
  * [13] A Hadjidimos and LL Zhang.  Comparison of three classes of algorithms for the solution of the linear complementarity problem with an h+-matrix.  Journal of Computational and Applied Mathematics, 336:175–191, 2018\.
  * [14] R Jana, AK Das, and A Dutta.  On hidden z-matrix and interior point algorithm.  Opsearch, 56:1108–1116, 2019.
  * [15] R Jana, AK Das, and S Sinha.  On processability of lemke’s algorithm.  Applications and Applied Mathematics: An International Journal (AAM), 13(2):31, 2018.
  * [16] R Jana, A Dutta, and AK Das.  More on hidden z-matrices and linear complementarity problem.  Linear and Multilinear Algebra, 69(6):1151–1160, 2021.
  * [17] B Kumar, Deepmala, and AK Das.  On general fixed point method based on matrix splitting for solving linear complementarity problem.  Journal of Numerical Analysis and Approximation Theory, 51(2):189–200, 2022.
  * [18] B Kumar, Deepmala, and AK Das.  More on modulus based iterative method for solving implicit complementarity problem.  arXiv preprint arXiv:2303.12519, 2023.
  * [19] B Kumar, Deepmala, and AK Das.  Projected fixed point iterative method for large and sparse horizontal linear complementarity problem.  Indian Journal of Pure and Applied Mathematics, pages 1–10, 2023\.
  * [20] B Kumar, A Deepmala, Dutta, and AK Das.  More on matrix splitting modulus-based iterative methods for solving linear complementarity problem.  OPSEARCH, pages 1–18, 2023.
  * [21] B Kumar, A Dutta, Deepmala, and AK Das.  Error bound for the linear complementarity problem using plus function.  arXiv preprint arXiv:2209.00377, 2022.
  * [22] CE Lemke and Joseph T Howson, Jr.  Equilibrium points of bimatrix games.  Journal of the Society for industrial and Applied Mathematics, 12(2):413–423, 1964.
  * [23] S Liu, H Zheng, and W Li.  A general accelerated modulus-based matrix splitting iteration method for solving linear complementarity problems.  Calcolo, 53:189–199, 2016.
  * [24] KG Murty and FT Yu.  Linear complementarity, linear and nonlinear programming, volume 3.  Citeseer, 1988.
  * [25] HS Najafi and SA Edalatpanah.  Modification of iterative methods for solving linear complementarity problems.  Engineering computations, 30(7):910–923, 2013.
  * [26] SK Neogy and AK Das.  On singular n0-matrices and the class q.  Linear algebra and its applications, 434(3):813–819, 2011.
  * [27] SK Neogy, AK Das, S Sinha, and A Gupta.  On a mixture class of stochastic game with ordered field property.  In Mathematical programming and game theory for decision making, pages 451–477. World Scientific, 2008.
  * [28] H Zheng, W Li, and S Vong.  A relaxation modulus-based matrix splitting iteration method for solving linear complementarity problems.  Numerical Algorithms, 74:137–152, 2017.