# Updated radial velocities and new constraints on the nature of the unseen
source in NGC1850 BH1

S. Saracino1, T. Shenar2,3, S. Kamann1, N. Bastian4,5,1, M. Gieles6,7, C.
Usher8,J. Bodensteiner9, A. Kochoska10, J. A. Orosz11, H. Sana3

1 Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK
2Anton Pannekoek Institute for Astronomy, Science Park 904, 1098XH Amsterdam,
The Netherlands
3 Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Leuven,
Belgium
4 Donostia International Physics Center (DIPC), Paseo Manuel de Lardizabal, 4,
20018, Donostia-San Sebastián, Guipuzkoa, Spain
5 IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain
6 ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain
7 Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),
Martí i Franquès 1, E08028 Barcelona, Spain
8 The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-106 91 Stockholm, Sweden
9 European Organisation for Astronomical Research in the Southern Hemisphere
(ESO), Karl-Schwarzschild-Str. 2, 85748 Garching, Germany
10 Villanova University, Deptartment of Astrophysics and Planetary Sciences,
800 East Lancaster Avenue, Villanova, PA 19085, USA
11 San Diego State University, Department of Astronomy, 5500 Campanile Drive,
San Diego, CA 92182-1221, USA E-mail<EMAIL_ADDRESS>
(Accepted 08 March 2023. Received 14 February 2023; in original form 16
December 2022)

###### Abstract

A black hole candidate orbiting a luminous star in the Large Magellanic Cloud
young cluster NGC 1850 ($\sim 100$ Myr) has recently been reported based on
radial velocity and light curve modelling. Subsequently, an alternative
explanation has been suggested for the system: a bloated post-mass transfer
secondary star (M${}_{\rm initial}\sim 4-5$ M⊙, M${}_{\rm current}\sim 1-2$
M⊙) with a more massive, yet luminous companion (the primary). Upon reanalysis
of the MUSE spectra, we found that the radial velocity variations originally
reported were underestimated ($K_{\rm 2,revised}=176\pm 3$ km/s vs $K_{\rm
2,original}=140\pm 3$ km/s) because of the weighting scheme adopted in the
full-spectrum fitting analysis. The increased radial velocity semi-amplitude
translates into a system mass function larger than previously deduced ($f_{\rm
revised}$=2.83 M⊙vs $f_{\rm original}$=1.42 M⊙). By exploiting the spectral
disentangling technique, we place an upper limit of 10% of a luminous primary
source to the observed optical light in NGC1850 BH1, assuming that the primary
and secondary are the only components contributing to the system. Furthermore,
by analysing archival near-infrared data, we find clues to the presence of an
accretion disk in the system. These constraints support a low-mass post-mass
transfer star but do not provide a definitive answer whether the unseen
component in NGC1850 BH1 is indeed a black hole. These results predict a
scenario where, if a primary luminous source of mass M $\geq 4.7$ M⊙ is
present in the system (given the inclination and secondary mass constraints),
it must be hidden in a optically thick disk to be undetected in the MUSE
spectra.

###### keywords:

globular clusters: individual: NGC 1850 – techniques: imaging spectroscopy,
photometry – techniques: radial velocities – binaries: spectroscopic

††pubyear: 2022††pagerange: Updated radial velocities and new constraints on
the nature of the unseen source in NGC1850 BH1–Updated radial velocities and
new constraints on the nature of the unseen source in NGC1850 BH1

##  1 Introduction

Recently, Saracino et al. (2022) reported the discovery of a black hole (BH)
candidate orbiting a luminous star in the massive young ($\sim 100$ Myr) star
cluster NGC 1850, in the Large Magellanic Cloud. Based on the measured radial
velocity and luminosity variations of the observed source, and its position in
the colour-magnitude diagram (CMD), the authors concluded that the source is a
main-sequence turn-off (MSTO) B-type star ($M\sim 4.9$ M⊙) and that the unseen
companion is an $\sim 11$ M⊙ BH. Furthermore, the authors suggested that the
system is in a semi-detached configuration meaning that the luminous star is
beginning to fill its Roche Lobe (they also studied the case of a detached
configuration). The system does not display obvious emission lines in the
optical region of the spectrum (although the presence of nebular contamination
combined with the low spectral resolution of the MUSE observations makes this
analysis complicated). However, a faint but significant X-ray detection
appears at the position of the source. The lack of a persistent X-ray emission
from NGC1850 BH1, although surprising, does not in itself exclude the presence
of a BH in the system. Low-mass X-ray binaries with both persistent and
transient X-ray emissions are indeed known in the literature (e.g., Cyg X-2,
Orosz & Kuulkers 1999 and V404 Cyg, Casares et al. 1992, respectively),
although neither of them can be directly compared to NGC1850 BH1.

A potential caveat to this discovery is that stars of different masses, which
have undergone different evolutionary paths, can display B-type spectra. As an
alternative explanation for this system, El-Badry & Burdge (2022) and Stevance
et al. (2022) have suggested that NGC1850 BH1 is a post-mass transfer binary
system, with the brighter source a bloated stripped star with a current mass
of $\sim 1-2$ M⊙ (M${}_{\rm initial}\sim 5$ M⊙) and the fainter source a more
massive star that has gained a lot of mass from the companion (M${}_{\rm
current}\sim 2-5$ M⊙). The latter is predicted to be significantly fainter (by
approx. 1-2.3 mag in the optical bands) than a main sequence (MS) star of the
same mass at the age of NGC 1850 due to rejuvenation episodes occurring during
mass transfer (Stevance et al., 2022), but see Wang et al. (2020) for an
alternative discussion on the impact of mass transfer on the luminosity of the
mass gainer.

We note here that there is a precedence for preferring such a configuration,
as previously suggested stellar-mass BH candidates LB-1 (Liu et al., 2019) and
HR 6819 (Rivinius et al., 2020) appear to be best explained instead as post-
mass transfer binary stars with two luminous companions (e.g., Shenar et al.,
2020; Bodensteiner et al., 2020; El-Badry & Quataert, 2021). One important
difference, however, between the LB-1 and HR 6819 systems compared to NGC1850
BH1, is that the former systems contain Be stars, i.e., fast rotating B-type
stars that display prominent emission lines, while no similar emission is
observed in the latter case (see Kamann et al. 2023 for a detailed study of
the sample of Be stars in NGC 1850).

Additionally, El-Badry & Burdge (2022) noted an inconsistency in the Saracino
et al. (2022) interpretation, namely that if the system is in a semi-detached
configuration, then a 5 M⊙ MSTO star would be more luminous than permitted by
the observed photometry. In the detached configuration, its implied radius
would instead be smaller than the Roche radius, and seems inconsistent with
the photometric variability suggesting a (near) Roche filling donor. On the
other hand, in the post-mass transfer model for NGC1850 BH1, we must be
catching the system at a unique time, specifically as it is transferring
across the Hertzsprung-Russell (HR) diagram from a cool bloated star to a hot
sub-dwarf state. The rarity of catching such a system at this time is
highlighted in Stevance et al. (2022), where the authors systematically
explored a large grid of pre-computed binary models (including mass transfer)
and could only find a matching system by significantly expanding the allowed
temperature range of the secondary ($\sim 10,000$ K compared to the observed
$\sim 14,500$ K). The chance of catching the luminous component as it crosses
the HR diagram from cool to hot, directly on the MS is then rather small
(approx. 1% of the lifetime of the system), but in principle easier to be
detected in this stage than in the later subdwarf stage (Bodensteiner et al.,
2020).

Upon further modelling of the NGC1850 BH1 system, we uncovered a systematic
bias in the published radial velocity measurements. This bias, which will be
accurately described in the following Sections, resulted in the
underestimation of the radial velocity semi-amplitude $K_{2}$111To avoid any
confusion in the reader, we specify here that throughout the paper, the
observed star is labelled with index 2 and called secondary, while the unseen
(more massive) object is labelled with index 1 and called primary. of the
visible source which in turn resulted in an underestimated mass function for
the system. In the present work we discuss the updated radial velocity
measurements in Section 2 along with the implications on the estimated orbital
properties of the system, especially the mass function. In Section 3 we
present upper limits to the presence of a luminous primary stellar component
in the system through the technique of spectral disentangling. In Section 4 we
focus on the visible secondary star and investigate a plausible lower limit in
mass for it. In Section 5 we combine these results and discuss the possible
nature of the unseen component based on the new constraints available.
Finally, in Section 6 we present our conclusions.

##  2 Revised radial velocity and mass function

Unlike what was done in Saracino et al. (2022), we present here an alternative
method to derive the relative radial velocities of the system, which relies on
cross-correlation of the observations with a template spectrum (Zucker &
Mazeh, 1994; Shenar et al., 2017; Dsilva et al., 2020). We perform the cross-
correlation in the range 7800-8900 Å, where several hydrogen lines of the
Paschen series are present. As a template, we first use one of the
observations themselves. Once a first set of radial velocities has been
determined, we compute the co-added spectrum, and use it as a new template for
measuring the radial velocities, repeating this process a few times until no
notable change in the radial velocities is observed. We convert the relative
radial velocities to absolute ones by using the systemic velocity of
$V_{0}=253.30\,$km ${\rm s}^{-1}$ measured by Saracino et al. (2022). Using
this new set of radial velocities we find consistent orbital parameters (e.g.
orbital period, eccentricity) to those derived by Saracino et al. (2022),
except for the radial velocity amplitude, which is found to be $K_{2}=175.6\pm
2.6$ km ${\rm s}^{-1}$. The orbital solution thus derived is shown in Figure 1
while the new single-epoch radial velocities are presented in Table 1.

In order to understand the discrepancy in $K_{2}$ values derived above and
reported in Saracino et al. (2022), who originally found $K_{2}$ = $140.4\pm
3.3$ km ${\rm s}^{-1}$, we performed additional full-spectrum fitting analyses
using the Spexxy code (Husser et al., 2016), which was used to measure the
velocities of the visible star in Saracino et al. (2022). We found that for
this particular star, the weighting scheme used by Spexxy has a significant
impact on the measured velocities. By default, Spexxy weighs the spectral
pixels by the inverse of their uncertainties during the fitting. If we switch
to a more physically motivated inverse-variance weighting scheme, we get
radial velocities consistent with those shown in Figure 1. Using both
weighting schemes, we then performed an analysis with Spexxy where we only
used the spectral range with $\lambda>7\,800\,\text{\AA}$, in effect using the
Paschen series as the only spectral lines in the fit. We found that either
weighting scheme (as well as using no weighting at all) resulted again in a
velocity curve consistent with the one shown in Figure 1. We repeated the
fitting with pPXF code (Cappellari & Emsellem, 2004; Cappellari, 2012) for
both the entire wavelength range and $\lambda>7\,800\,\text{\AA}$, finding
identical results as with Spexxy.

Given the high effective temperature of the observed star, its flux is much
higher in the blue part of the MUSE spectral range than in the red (see Figure
2 in Saracino et al., 2022). As a consequence, when using the inverse
uncertainties as weights, the blue part (with the strong H$\beta$ and
H$\alpha$ lines) has a larger impact on the fit than the red part (containing
the Paschen series). Ideally, this over-weighting of the blue part would not
affect the kinematics, as all lines are shifted by the same radial velocity.
In the case of NGC 1850, however, the strong nebular emission, associated with
the 5 Myr old cluster NGC 1850B, represents an additional complication for the
data analysis. As both H$\beta$ and H$\alpha$ are particularly strong in the
nebular line spectrum, it is conceivable that residual nebular emission
contaminates the two line profiles sufficiently so that the velocity estimates
based on them are biased towards the cluster mean. An alternative, more
physical explanation, could be that accretion onto the unseen companion
creates some H$\alpha$ (and potentially H$\beta$) emission that partially
fills up the absorption line. However, as the emission component would follow
the motion of the unseen primary and hence appear with a phase shift of 180∘
in velocity space, one would expect an over- rather than an under-estimation
of the $K_{2}$ value derived using the blue part of the MUSE spectral range
(see the discussion in Abdul-Masih et al. (2020) for the LB-1 system).

When comparing the velocity curve shown in Figure 1 to the one depicted in
Figure 5 of Saracino et al. (2022), one can see that the scatter of the
individual velocity measurements around the Keplerian model is significantly
reduced when the revised measurements are used (reduced $\chi^{2}$ is 0.52).
Given this improvement, and the potential issues regarding the usage of the
H$\beta$ and H$\alpha$ lines, we give preference to the new result. The
difficulty in determining $K_{2}$ discussed here highlights the need to study
the system at higher spectral resolution over a broad wavelength range. This
would allow: 1) for a better cleaning of the nebular emission lines (which
would be significantly narrower in high-resolution data), also thanks to
strong metal lines in the blue that do not appear in the nebular emission; 2)
to add stricter constraints from the spectral disentangling technique; and 3)
to increase the chances of finding a potential emission-line contribution from
accretion onto the unseen companion.

###  2.1 Mass Function

A change in the derived semi-amplitude velocity $K_{2}$ of the visible source
in NGC1850 BH1 has a direct effect on the mass function of the system, even if
all other orbital parameters (e.g. period, eccentricity) stay the same. In
fact, based on the formula of the binary mass function (Remillard &
McClintock, 2006), that can be expressed in terms of observational quantities
as:

$f=\frac{P_{\rm orb}K_{\rm 2}^{\rm 3}(1-e^{2})^{3/2}}{2\pi G},$ (1)

which does not make any assumptions on the mass of the visible source, we
obtain $f=2.83^{+0.14}_{-0.12}$ M⊙, significantly higher than $f=1.42$ M⊙ as
derived in Saracino et al. (2022). Orbital period $P_{\rm orb}$ and
eccentricity $e$ are the same as in Table 2 of Saracino et al. (2022). This
implies that, regardless of the mass of the visible star (a normal MS star vs
a bloated star), the unseen primary companion is substantially more massive
than previously predicted. All the revised and relevant properties of NGC1850
BH1 are listed in Table 2, to provide the reader with a clearer reference.

By using Kepler’s third law, the binary mass function can also be written in
the form:

$f=\frac{M_{\rm 1}^{3}\sin(i)^{3}}{(M_{\rm 1}+M_{\rm 2})^{2}},$ (2)

where $M_{\rm 1}$ and $M_{\rm 2}$ are the masses of the primary unseen
component and the secondary visible star, respectively, and $i$ the
inclination of the system with respect to the line of sight. This formula
suggests that once the mass of the visible star and the inclination of the
system are known, the mass of the unseen companion can be determined.
Unfortunately, these two additional quantities are uncertain in the case of
NGC1850 BH1. In Section 4 we will define an alternative way to put constraints
on the mass of the unseen source.

Figure 1: The revised MUSE radial velocity curve of the luminous secondary star in NGC1850 BH1, phase-folded using its orbital period of P = 5.04 d, is shown as black dots. The red solid line represents the new best-fitting model (reduced $\chi^{2}=0.52$, RMS = 13.4 km ${\rm s}^{-1}$). The bottom panel shows the residuals of the comparison between the observed radial velocities and the best fit model. Table 1: Updated radial velocities using the method outlined in Section 2. Time (MJD) | VR (km/s) |  $\sigma$ VR (km/s)
---|---|---
58497.08534751 | 379.44 | 10.12
58498.15614836 | 407.7 | 14.1
58498.17027805 | 404.9 | 25.9
58550.02867354 | 115.18 | 14.15
58550.04180262 | 155.2 | 27.7
58553.01788807 | 418.4 | 10.3
58553.03123202 | 442.5 | 38.0
58556.01231551 | 100.4 | 12.4
58556.02541888 | 100.5 | 27.7
59174.27808058 | 252.7 | 14.0
59174.32203504 | 255.1 | 27.1
59175.16932389 | 103.94 | 12.04
59175.29610617 | 81.6 | 17.7
59176.13916114 | 116.8 | 13.7
59176.30463698 | 149.7 | 14.8
59177.30703879 | 365.2 | 25.1
59190.19805799 | 110.6 | 14.3
59201.25318966 | 120.9 | 13.8
59203.14316930 | 428.29 | 11.16
59251.14817479 | 74.4 | 14.8
Table 2: Revised properties of NGC1850 BH1 Period $P_{\rm orb}$ |  $5.0402\pm 0.0004$ d
---|---
Velocity semi-amplitude $K_{2}$ |  $175.6\pm 2.6$ km/s
Barycentric radial velocity $v_{0}$ |  $253.30^{+2.59}_{-2.44}$ km/s
Mass function $f$ |  $2.83^{+0.14}_{-0.12}$ M⊙
Eccentricity $e$ | $0.029^{+0.010}_{-0.014}$

##  3 Spectral disentangling

Based on the newly derived mass function, which points towards a rather
massive unseen companion, we used the MUSE spectra available to set an upper
limit on how much this object actually contributes light to the total flux of
the system. In fact, if the unseen companion is a massive star as suggested by
El-Badry & Burdge (2022) and Stevance et al. (2022), it is rather luminous, so
it is expected to contribute significantly to the total flux of the system
(but see the discussion about the rejuvenation factor in Section 4). If the
unseen companion is instead a compact object (such as a BH) as suggested by
Saracino et al. (2022), it does not contribute to the light of the system at
all if there is no accretion disk around it, regardless of its mass. To do
this test, we employed the shift-and-add spectral disentangling technique
(Marchenko et al., 1998; González & Levato, 2006; Shenar et al., 2020, 2022),
which was successfully used to uncover hidden companions in other SB1 binaries
(e.g., LB-1, Shenar et al. 2020; HR 6819, Bodensteiner et al. 2020; 28 O-type
binaries, Shenar et al. 2022), which have companions contributing as little as
$\approx 1-2\%$ to the visual flux.

Briefly, spectral disentangling is the separation of composite spectra into
the component spectra of multiple systems, usually performed simultaneously to
the derivation of the orbital parameters (Hadrava, 1995; Bagnuolo & Gies,
1991; Mahy et al., 2012). For given orbital elements, the shift-and-add method
relies on an iterative procedure that uses the disentangled spectra obtained
in the $j^{\rm th}$ iteration, to calculate the disentangled spectra for the
$j+1$ iteration through consecutive shifting-and-adding of the spectra. By
minimizing $\chi^{2}$ between the added disentangled spectra and the
observations, one can derive the orbital elements; we refer to Shenar et al.
(2020, 2022) for details. Here, we fix the orbital parameters to those given
in Table 2 of Saracino et al. (2022), except for the radial velocity
amplitudes $K_{1},K_{2}$, which are used to minimise $\chi^{2}$. We note that
the light ratio of the components cannot be derived from the disentangling
procedure. The adopted light ratio only impacts the final scaling of the
spectra.

In Figure 2, we show the reduced $\chi^{2}$($K_{1},K_{2}$) map obtained when
disentangling the four Paschen members (members 8 - 11) in the region
$8570-8910\text{\AA}$. Evidently, $K_{2}$ can be reasonably well constrained
and is consistent with the radial velocity measurements in Table 2 to within
1$\sigma$. In contrast, the value $K_{1}$ is poorly constrained, virtually
ranging across the entire plausible range of values. We note that
disentangling generally yields much larger formal errors compared with
standard radial velocity fitting due to the freedom in varying each pixel in
the disentangled spectrum. In the figure a slight correlation between $K_{\rm
1}$, $K_{\rm 2}$ is observed. This may possibly indicate that there is some
contributing signal from a primary star or disk in NGC 1850 BH1, although this
contribution is too small to actually be extracted from the noise using the
MUSE data. The presence of a putative accretion disk or a luminous primary
(see the discussion in Section 5) could provide the light signal to explain
such a trend in the residual map. Alternatively, this correlation could be a
spurious result caused by contaminants (e.g., nebular contamination,
tellurics, uncertain normalisation).

Figure 3 shows the disentangled spectra for $K_{2}=175.6\,$km ${\rm s}^{-1}$
and $K_{1}=41\,$km ${\rm s}^{-1}$ (i.e., assuming the primary is roughly four
times more massive than the luminous secondary). The shifted spectra and their
sum are compared to the observations at radial velocity extremes (conjunction
phases). Generally, the disentangled spectra of the primary appear close to
flat, with the possible exception of H i $\lambda 8750$. We note that the
results depend weakly on the adopted value of $K_{1}$ (values in the range
$0.25\,K_{2}<K_{1}<4\,K_{2}$ were considered). In all cases, the features seen
in the spectrum of the primary are comparable to the noise level of the
disentangled spectrum.

In Figure 4, we show the disentangled spectra of a few neighbouring Paschen
lines calculated for $K_{2}=175.6\,$km ${\rm s}^{-1}$ and $K_{1}=41\,$km ${\rm
s}^{-1}$, assuming a low light contribution for the primary of $l_{1}=0.1$
(i.e., the intrinsic spectrum is multiplied by a factor of 10). The features
that are observed in the disentangle spectra are again of level of the noise,
and generally do not overlap with spectral lines. Such features can easily
result from non-Gaussian noise, imperfect normalisation, tellurics, or other
contaminants. The results imply that, if a non-degenerate companion is
present, it must be rather faint. This is corroborated by the simulations
below (Section 3.1), where this statement is further quantified.

Figure 2: $\chi^{2}$($K_{1},K_{2}$) from disentangling the spectra in the
wavelength region from 8570 to 8910 Å. A slight correlation between $K_{\rm
1}$ and $K_{\rm 2}$ is observed (see the dashed green line).
---
Figure 3:  Disentangled spectra of the primary and secondary for the H i
$\lambda 8863,8750,8665$ lines (top, middle, and bottom panels) and their sum,
compared to observations at radial velocity extremes (left and right panels).
The spectra are calculated for $K_{2}=175.6\,$km ${\rm s}^{-1}$ and
$K_{1}\approx K_{2}/4\,$ (41 km ${\rm s}^{-1}$). The spectra are not scaled by
the light ratio in this figure. The results depend weakly on $K_{1}$. The
spectrum of the primary appears featureless, with the possible exception of
the H i $\lambda 8750$.

###  3.1 Simulations

Figure 4: Disentangled spectra of the primary and secondary in NGC1850 BH1,
obtained for $K_{2}=175.6\,$km ${\rm s}^{-1}$and $K_{1}=41\,$km ${\rm
s}^{-1}$, and assuming a light contribution of $l_{1}$ = 10% for the unseen
primary.  Figure 5: The disentangled spectra of the bright secondary (top) and
faint primary (bottom) of our simulated binary, compared with the input
templates. The disentangling was performed using the same input orbital
parameters as used for the simulation, but for $K_{1}=20,50$, and $150\,$km
${\rm s}^{-1}$(see legend), illustrating the minor impact of $K_{1}$ on the
spectral appearance of the secondary.

To test the validity of our method and explore the sensitivity down to which
we could detect a hidden MS companion, we simulate a mock binary that mimics
the orbit and observational data set of NGC1850 BH1, but contains a non-
degenerate companion. For the simulation, we use the co-added spectrum as a
template for the luminous secondary. For the mock spectrum of the unseen
primary, we use the grids computed with the TLUSTY model atmosphere code
(Hubeny & Lanz, 1995; Lanz & Hubeny, 2003, 2007). We use a model with $T_{\rm
eff}=20,000\,$K, $\log g=4.0\,$[cgs], and assume that it moves with $K_{\rm
1,true}=50\,$km ${\rm s}^{-1}$. To be conservative, we convolve the emergent
spectrum of the model with $\varv\sin i=300\,$km ${\rm s}^{-1}$ and a
macroturbulent velocity of $\varv_{\rm mac}=30\,$km ${\rm s}^{-1}$. Finally,
motivated by the results in Section 3, we adopt a low light contribution for
the primary of $l_{1}=0.1$. The mock observations use the exact S/N values and
phases of the original spectra, and are degraded to the MUSE resolution and
sampling.

We then attempted to derive $K_{1}$ through $\chi^{2}$ minimisation. However,
the $\chi^{2}$ map is virtually flat, implying that $K_{1}$ cannot be
retrieved. Given the low resolution (compared to the secondary’s radial
velocities), the intrinsic and rotational broadening of the lines, and the
modest S/N, it is not surprising that $K_{1}$ cannot be retrieved. However,
disentangling can still be performed, assuming various values of $K_{1}$. In
truth, the $K_{1}$ value has a very small impact on the spectral appearance of
the disentangled spectra, as long as it is varied in a plausible range. To
illustrate this, in Figure 5, we show the results from three disentangling
experiments of the mock data, varying $K_{1}$ between 20 and 150 km ${\rm
s}^{-1}$. The spectra are virtually indistinguishable. This seeming
independence on $K_{1}$ is the result of the broad profiles of the simulated
primary and the low spectral resolution of the data.

Evidently, while we cannot retrieve $K_{1}$, the method yields a spectrum for
the hidden primary that matches the original template reasonably well. Some
differences are apparent for the primary, which are intrinsic to the method.
Since the lines are constantly blended, the disentangling procedure is bound
to have some cross-contamination between the stars. However, we note that the
differences are boosted by a factor of ten due to the faintness of the
primary, such that the deviations seen in Figure 5 amount to deviations of the
order of a few percent with respect to the mock observations. The exact
sensitivity down to which we could detect companions is difficult to
establish, since it depends on the stellar parameters, rotation, and light
contribution of the primary. However, the experiment described here
illustrates that we would very likely be able to detect companions
contributing more than $\approx 5-10\%$ to the light.

We visually illustrate how faint the primary (unseen) star must be in terms of
magnitude to be undetectable in the MUSE spectra, based on the results of the
spectral disentangling. In the left panel of Figure 6 we present the CMD of
NGC 1850 where a MIST isochrone (Choi et al., 2016) of the appropriate age is
overplotted to guide the eye. A red star indicates the position of NGC1850 BH1
(F438W = 16.7, F814W = 16.6) while the red solid line indicates the magnitude
level of a MS star (F438W = 19.2, F814W = 19.1) corresponding exactly to 10%
of the brightness of NGC1850 BH1, the limit set by the spectral disentangling.

The conclusions of these tests are twofold: First, if a non-degenerate stellar
companion is present in the binary, as suggested by El-Badry & Burdge (2022),
it is fainter than $\approx 10\%$ in the visual. Second, even if a non-
degenerate companion is present, it cannot significantly contaminate the
spectrum due to its faintness. The stellar parameters determined for the
luminous secondary using the co-added observations (which are virtually
identical to the disentangled spectrum) should therefore represent the
secondary well, unless its light is diluted by an additional sources (e.g.,
excess emission stemming from a disk). In fact, if there should be an
accretion disk in NGC1850 BH1 orbiting around the invisible source, then the
results presented above could change to account for this additional component.
An extensive discussion of this aspect will be provided in Section 5.

Figure 6: Left panel: (F336W-F438W, F438W) CMD of NGC 1850, with a MIST
isochrone of 100 Myr overplotted. NGC1850 BH1 is represented by a red star in
the figure. The red solid line shows the magnitude level corresponding to a
brightness of only 10% of our target (the limit derived by the spectral
disentangling). Right panel: The CMD of NGC 1850 with the same MIST isochrone
overplotted. The closed vs open blue dots indicate the magnitude level of a
standard (non-interacting) MS star with a mass of 5 M⊙(i.e. the mass of a star
as bright as NGC1850 BH1) and 4.7 M⊙(i.e. the minimum mass for the primary
star in NGC1850 BH1, see the text), respectively. The closed vs open green
dots instead show how the same two stars look like when the rejuvenation
factor due to mass transfer in the binary is applied. They appear $\approx$
1.5 mag fainter, but still above (by 1 mag and 0.4 mag, respectively) the red
solid line which defines the 10% brightness limit in F438W.

##  4 Minimum mass of the secondary and its implications for the primary

An additional hint on what is the nature of the invisible source can actually
be derived from the analysis of the visible secondary star. The two
alternative scenarios that have been proposed to explain NGC1850 BH1 assumed
very different masses for the visible source (a normal 5 ${\rm M}_{\odot}$ MS
star vs a 1-2 ${\rm M}_{\odot}$ bloated-stripped star). El-Badry & Burdge
(2022) argued that the secondary star, if it is indeed filling its Roche lobe,
must be less massive than the 5 ${\rm M}_{\odot}$ adopted by Saracino et al.
(2022). According to Eggleton (1983), who derived a formula for the mean
density of Roche-lobe filling stars, a 5 $M_{\odot}$ star would be too large
(and hence too luminous) to satisfy the photometric constraints. Therefore, we
are inclined to adopt the scenario in which the secondary is a low-mass post
mass transfer star. Unfortunately, deriving its actual current mass is not
possible with the available data, because it is unknown if and how much other
sources contribute to the observed magnitudes. Because of this limitation, in
this work we define a physically motivated minimum mass for the secondary star
and adopt this value in the subsequent analysis. This lower limit directly
translates into a minimum mass for the unseen primary star as well, once the
inclination of the system is known or can be set to a reasonable value.

From the OGLE light curves available for NGC1850 BH1 and presented in Saracino
et al. (2022), the system does not show any evidence for total or partial
eclipses. This might be the case for two reasons: 1) the binary system is made
of two luminous sources but it has an inclination such that one source never
obscures the other; 2) the unseen source is a dark object (e.g. a BH) so it
does not produce eclipses regardless of how the system is inclined. We note
here that if the BH is surrounded by an accretion disk, eclipses are still
expected unless the system has a geometric configuration similar to that of
1).

According to Beech (1989), in a binary system, the geometric condition for
eclipses not to occur is the following:

${\rm\cos}(i)>(R_{\rm 1}+R_{\rm 2})/a,$ (3)

where $a$ and $i$ are the semi-major axis and the inclination of the system,
and $R_{\rm 1}$ and $R_{\rm 2}$ the radii of the primary and secondary stars,
respectively. Based on the observational constraints derived from the
modelling of the radial velocity curve and the constraints on the radius of
the luminous (secondary) star (4.9 ${\rm R}_{\odot}\leq R_{\rm 2}\leq 6.5{\rm
R}_{\odot}$) imposed by observational uncertainties and the possibility that a
second luminous star could contribute to the observed photometry, the lack of
eclipses in the OGLE light curves places a limit of $i\leqslant 67^{\circ}$ on
the inclination of NGC1850 BH1 (see also El-Badry & Burdge 2022).

In Figure 7 we show, as red solid lines, the mass of the primary (unseen)
component as a function of the luminous secondary component, based on the
newly measured binary mass function and by adopting equation (2) above, for
two different inclinations: when the binary is seen edge-on ($i=90^{\circ}$)
and when the inclination is $i=67^{\circ}$, as labelled in the plot. The red
shaded area is the region in this parameter space where eclipses are expected
to occur, while the white area above the red area is where no eclipses are
expected to be observed. In other words, the red solid line at $i=67^{\circ}$
sets the lower limit to the mass of the primary (as a function of the
secondary) if the primary is a star or a BH with an accretion disk. If the
primary is instead a dark compact object not surrounded by an accretion disk,
the reference line to be considered is the red solid line at $i=90^{\circ}$.

El-Badry & Burdge (2022) estimate the current mass of the luminous component
(a bloated stripped star in their model) to be $\sim 1-2$ M⊙. In their model,
the secondary (initially most massive) component would have recently left the
MS and began expanding. During this expansion, the Roche Lobe of the star
would have been filled, and mass transfer onto the primary companion would
have followed. At the age of NGC 1850 ($\sim 100$ Myr), this implies that the
initial mass of the secondary would have been $\sim 5$ M⊙. This is consistent
with what Götberg et al. (2018) found in their binary interaction models: when
assuming an initial mass of $\sim 5$ M⊙ for the secondary star, after the mass
transfer, the final mass of the stripped star turns out to be of $\sim 1-2$
M⊙. Although a mass range between 1 and 2 M⊙ is in agreement with current
binary evolution models (e.g. Götberg et al. 2018), it is worth mentioning
here that there is not yet enough information about the binary system to allow
us to assign a value to the present mass of the visible star. What we do here
instead is test how massive the primary would have to be under the assumption
of a given secondary mass (and in particular the value of $1$ M⊙ proposed by
El-Badry & Burdge 2022).

Indeed, by assuming a current mass of the stripped star of $\sim 1-2$ M⊙, and
by applying the no-eclipse condition, the mass of the unseen component is $>5$
M⊙. In particular, for a $1$ M⊙ secondary mass, the primary has a minimum mass
of $M_{1}=5.17$ M⊙, as also highlighted with a black dashed line in Figure 7.
An accretor star of this mass would be as luminous as the visible star itself,
so its contribution should be clearly detectable in the MUSE spectra according
to the brightness limit set by the spectral disentangling (see Figure 6, left
panel).

Based on the BPASS models (Eldridge & Stanway, 2016), Stevance et al. (2022)
pointed out that in a post-mass transfer system, the star that has gained a
considerable amount of mass from the companion does not look like a standard
(non-interacting) MS star of the same mass and age (in terms of brightness),
but instead experiences an episode of rejuvenation, so that in the end it
looks much fainter (by up to approx. 2.3 mag in the optical filters). The
uncertainties associated with this process are quite large and different
binary models tend to predict different scenarios. For example, by using the
MESA binary evolution models (Paxton et al., 2015), Wang et al. (2020)
recently found that stars who gained a significant amount of mass from their
companions in binaries are systematically brighter and more rapidly rotating
than they were pre-interaction. This shows that the rejuvenation factor in
binary models is still an open question.

While the real factor is somewhat uncertain, we decide to be conservative here
and adopt a value of 1.5 mag (fainter) for the rejuvenation factor in the
analysis hereafter. Based on this assumption, a primary of $M_{1}=5.17$ M⊙ is
significantly fainter than expected from stellar evolution (F438W$\sim$18.2 vs
F438W$\sim$16.7) but well detectable in the spectra, as it is still one
magnitude brighter than the 10% brightness limit of F438W = 19.2 imposed by
the spectral disentangling222Since the 10% flux limit was determined using
Paschen lines ($\lambda>7\,800\,\text{\AA}$), we verified that the same
behavior is also observed using the F814W filter. In fact, the 10% limit
corresponds to F814W = 19.1, while a primary as bright as the visible star in
NGC1850 BH1 would appear to be F814W = 16.6 + 1.5 = 18.1, still one magnitude
above the limit.. This is shown in the right panel of Figure 6, where the
position of a standard (non-interacting) MS star as bright as NGC1850 BH1 is
shown as a closed blue dot, while the same rejuvenated star as a closed green
square, overplotted on the CMD of NGC 1850. The impact of the rejuvenation
factor on the brightness of the primary (unseen) source will be illustrated
more clearly later in the Section, when we define the minimum mass that a
primary mass can assume.

Figure 7: Secondary mass $M_{2}$ vs Primary Mass $M_{1}$ for NGC1850 BH1. The
red shaded area defines the region where eclipses are observed, while for
$i=67^{\circ}$ or below (top white area) no eclipses are observed. Our target,
NGC1850 BH1, does not show eclipses and assuming $M_{2}=1~{}$M⊙ as suggested
by El-Badry & Burdge (2022), we obtain a minimum mass $M_{1}=5.17~{}$M⊙ for
the primary (black dashed line). This is also the case for the lowest
plausible mass for the secondary star that we set to $M_{2}=0.65~{}$M⊙ based
on what also suggested by El-Badry & Burdge (2022) (dotted line in the plot).
The grey shaded area sets the physically motivated lower limit to the mass of
the secondary star (see text for details). The white region on the bottom is
excluded by the measured binary mass function.

Both Stevance et al. (2022) and El-Badry & Burdge (2022) explored evolutionary
scenarios that could produce the binary NGC1850 BH1 using BPASS (Eldridge &
Stanway, 2016) and MESA (Paxton et al., 2015), respectively. They could
reproduce the observational properties of the binary system by assuming that
the luminous star was a bloated stripped star of $\sim 1$ M⊙. While it seems
very unlikely that the secondary star is significantly less massive than
predicted in their models, we want to be conservative here and set the minimum
mass for the visible star in NGC1850 BH1 at $M_{2}=0.65\,{\rm M}_{\odot}$,
which is the lower limit that El-Badry & Burdge (2022) derived in their work
on the basis of the minimum possible radius ($R_{\rm 2}=4.9{\rm R}_{\odot}$)
this star can assume given its temperature, its CMD position, and the possible
presence of a primary luminous source contributing to the photometry.
Moreover, since the observed spectrum of the luminous source contains
prominent hydrogen lines and easily resembles that of a normal B-type star, it
is reasonable to infer that this source has not yet been completely stripped
of its entire hydrogen envelope, hence its mass cannot be very low.

By imposing a lower limit for the secondary mass of 0.65 M⊙, we consider any
mass below this threshold as non-physical and present it as a grey shaded area
in Figure 7. The dotted line in the Figure instead shows that a secondary mass
of $M_{2}=0.65~{}{\rm M}_{\odot}$ directly translates into a mass of
$M_{1}=4.71~{}{\rm M}_{\odot}$ for the primary (unseen) component based on the
mass function of the system. Assuming it is a normal MS star (F438W $\sim$
17.3, F814W $\sim$ 17.25), this mass corresponds to a brightness of $\sim 56$%
of that observed for NGC1850 BH1 in the visual (F438W). If we take into
account the rejuvenation factor, which makes this star appear fainter by 1.5
mag, we deduce for it a magnitude F438W = 17.3 + 1.5 = 18.8 (F814W = 18.75),
which has $\sim$ 14.5% of the brightness of NGC1850 BH1. It would be barely
but still observable in the spectrum, given the 10% limit found by the
disentangling. To illustrate this, the right panel of Figure 6 shows the CMD
of NGC 1850, with the MIST isochrone appropriate for the cluster. The open
blue dot in the Figure represents the F438W magnitude of a primary (unseen)
star with the minimum allowable mass, 4.71 ${\rm M}_{\odot}$ (assuming it to
be a MS star). The open green square instead shows its position once the
rejuvenation factor of 1.5 mag is applied, due to mass transfer in the binary.
Even for the lowest possible primary mass (hence faintest), this star is
expected to be visible, as it is still brighter (by 0.4 and 0.35 mag in F438W
and F814W, respectively) than the brightness limit set by the disentangling
(here shown as a red solid line). More massive, rejuvenated, primary stars
would appear even brighter, thus even easier to detect in the spectra.

From this analysis we can conclude that, if we want to have a system with two
luminous stars, for any reasonable secondary (present-day) mass for the
visible star, we must invoke the presence of an additional component in the
system which by shielding the light of this massive luminous source causes it
to contribute very little to the total light.

##  5 The unseen source in NGC1850 BH1

The combination of i) the measured high binary mass function
($f=2.83^{+0.14}_{-0.12}$ M⊙), ii) the lack of observed eclipses in the OGLE
optical light curves, and iii) the lower limit of $M_{\rm 2}=0.65\,{\rm
M}_{\odot}$ set to the mass of the visible star results in the primary
(unseen) companion of NGC1850 BH1 to have a mass $M_{\rm 1}>4.71$ M⊙. A
standard (non-accreting) MS star of such a mass would be more than half as
bright as the observed system, hence its contribution would be easily
detectable in the MUSE spectra of the visible companion. However, when the
rejuvenation factor of 1.5 mag is taken into account, the brightness of this
source drops from 56% to 14.5% brighter than NGC1850 BH1 (see Figure 7, right
panel). It is important to note that, although much fainter, it does not go
below the 10% brightness constraint we derived from the spectral
disentangling. In other words, the contribution of the primary in this
configuration is expected to be faint but still detectable in the spectra.

To reconcile this result with the fact that we do not observe any contribution
from the primary in the analysis of the MUSE spectra, there are two viable
possibilities to explore: First, the unseen primary star is a non-luminous
compact object, and since its minimum mass is higher than that of any possible
neutron star (${\rm M}$$\sim$3 ${\rm M}_{\odot}$; Lattimer & Prakash 2001), it
is a BH. Second, the unseen primary star is a rather massive luminous source
enshrouded in a thick accretion disk which absorbs part of its optical light,
making it undetectable. Which of these two possibilities has to be preferred
is unclear to date but the scope of this paper is to present the current
knowledge about NGC1850 BH1 and suggest possible ways to distinguish one or
the other scenario with further observations.

As already noted by Saracino et al. (2022), the NGC1850 BH1 system belongs to
a class of objects called Double Period Variables (DPVs, Mennickent et al.
2003), i.e. it shows two periodicities where one is about 33 times longer than
the other. There is not much literature on DPVs, and the origin of the longest
periodicity is still unknown, however, the general consensus is that these
systems are semi-detached (one of the two components fills its Roche Lobe) and
made up of two stars, one of which (the gainer) is typically a B-type star
surrounded by an accretion disk. In one specific case, HD 170582, it has been
suggested that the gainer is bright and massive and should contribute nearly
50% of the total system light, but since it is encased in an optically thick
disk that almost completely obscures it, it only contributes for about 10%,
thus becoming barely detectable (Mennickent et al., 2015). For sake of
completeness, it is worth mentioning that the accretion disk ($\sim 21\,{\rm
R}_{\odot}$) and the semi-major axis ($\sim 61\,{\rm R}_{\odot}$) deduced for
the system HD 170582 are much larger than allowed for the configuration of
NGC1850 BH1. This is an example but an in-depth comparison of NGC1850 BH1 with
the properties of other DPV systems is beyond the scope of this document.

By analyzing archival near-infrared HST/WFC3 data of NGC1850 BH1 (from F105W
to F160W), we measured a 2$\sigma$ excess both in F140W and F160W compared to
other cluster members in a similar position in the CMD, which seems to support
the presence of a third component in the system, namely an accretion disk.
Figure 8 shows an optical/near-infrared (F438W-F160W, F438W) CMD of NGC 1850,
where the position of NGC1850 BH1 is presented as a red star. As a comparison,
we highlight in green and yellow respectively a sample of Be and shell stars
of NGC 1850 studied in Kamann et al. (2023). Shell stars are Be stars (i.e.
rapidly rotating B stars) observed (partially) through their disks (Rivinius
et al., 2006, 2013). As shown in the Figure, both Be and shell stars exhibit
near-infrared excesses, i.e. they are systematically located on redder colors
than normal stars of similar magnitudes. This excess is believed to be mainly
caused by emission from their disks. NGC1850 BH1 shares a similar colour,
hence a similar near-infrared excess, with many of the shell stars in the
cluster so, although the nature of NGC1850 BH1 is very different from that of
shell stars, an analogy between these sources can still be drawn. The observed
excess of the binary provides further support for the existence of a disk in
the system. Additional constraints (e.g. the expected slow rotation of the
secondary (luminous) star for a synchronized binary) suggest that, if present,
the disk is around the primary (unseen) star.

If the evidence for an accretion disk around the unseen source in NGC1850 BH1
is confirmed by further studies, this will unfortunately still not provide a
final answer as to what is the unseen component. Given the high probability
that NGC1850 BH1 is a post-mass transfer system, based on recent findings and
reasoning, it would be equally plausible to have a disk around a BH or a
massive luminous star. The only way to effectively discriminate between the
two scenarios is to measure the temperature of the putative disk itself as it
is expected to be very different in the two configurations. In particular, if
a 5 M⊙ luminous star (and $T\textsubscript{eff}\sim$ 15,000K) is enclosed in
an optically thick disk, the light it emits is almost completely absorbed by
the disk at optical wavelengths and re-emitted by it at infrared wavelengths.
The system thus becomes particularly bright in the near and mid-infrared,
given the lower extinction of starlight and the added contribution of the
disk, which is expected to be much cooler than the star itself
($T\textsubscript{eff}\ll$ 15,000K). Alternatively, if a BH is part of the
system, the properties of the accretion disk are significantly different, with
a temperature likely higher than 15,000K near the inner edge, but decreasing
with radius as predicted by Shakura & Sunyaev (1973). Near and mid infrared
observations of the NGC1850 BH1 system, as those recently secured with the new
ERIS/NIX imager at the VLT (Davies et al., 2018), will help to investigate
this aspect in more detail.

Figure 8: HST/WFC3 near-infrared CMD of NGC 1850. All stars in the cluster are
shown as grey dots. Our target, NGC1850 BH1, is instead highlighted as a red
star, along with a sample of Be and shell stars spectroscopically identified
in the cluster by Kamann et al. (2023) and presented as green and yellow dots,
respectively. Both Be and shell stars show a near-infrared excess due to disk
emission. Although the nature of NGC1850 BH1 is very different, this system
exhibits a similar color to many of the shell stars, supporting the idea that
an accretion disk is also present in this system.

##  6 Conclusions

The main results of this study can be summarised as follows:

  * •

We have discovered a systematic bias in the measured radial velocities which
has been traced to the weighting scheme adopted in the SPEXXY results. We have
provided updated radial velocity measurements for each epoch. Based on the new
modelling of the radial velocity curve, we have updated the radial velocity
semi-amplitude to $K_{2}=175.6\pm 2.6$ km ${\rm s}^{-1}$.

  * •

The increased (by 20%) semi-amplitude velocity thus derived has significantly
increased the mass function of the system to $f=2.83^{+0.14}_{-0.12}$ M⊙.

  * •

From spectral disentangling we find that only one source is significantly
contributing to the spectrum, i.e., any possible stellar secondary contributes
at most $10$% to the optical flux of the system.

  * •

The secondary (visible) star is most likely a low-mass post-mass transfer
star, but the information available so far does not allow us to assign a value
to the present mass of this binary component. Indeed, it is unknown if and how
much (other) sources contribute to the observed magnitudes.

  * •

Based on the new binary mass function, lack of observed eclipses in the light
curves of NGC1850 BH1, and constraints on the luminosity of the system’s
components, there are two viable possibilities: the unseen component in
NGC1850 BH1 is: 1) a BH, its mass being $>3$ M⊙ with the possible addition of
an accretion disk; or 2) a bright, rejuvenated star with a minimum mass of
$M_{1}\sim 4.7$ M⊙, enshrouded in an optically thick disk that partially
absorbs its light so that it is undetectable in the currently available
spectra.

  * •

NGC1850 BH1 is a DPV and appears to show an excess in the near-infrared, which
can be interpreted as evidence for the presence of a disk in the system.
However, both scenarios are still equally likely. Constraining the properties
of the disk (e.g. temperature, size) will be one good way to shed more light
on the nature of the invisible source.

  * •

A scenario in which the primary (unseen) component in NGC1850 BH1 is a BH
faces substantial issues regarding its evolutionary history, if we assume a
binary origin for it as in the BPASS and MESA models, i.e. the initial period
of the binary would be lower than physically allowed given the size of the
individual components. A possible caveat of these models, however, is that
they only consider isolated binaries and do not include hierarchical
triples/quadruples nor the effect of dynamical interactions in clusters. This
might instead be appropriate for NGC1850 BH1 which belongs to NGC 1850. In
conclusion, since neither the exact configuration of the binary (in terms of
$M_{1}$, $M_{2}$, mass ratio etc.) nor its evolutionary history are known, we
unfortunately cannot draw any definitive conclusion on this aspect.

In a future study we will present detailed modelling of the OGLE light curves
of NGC1850 BH1, which also includes the presence of an accretion disk. This
will provide stricter constraints on the nature of both the luminous secondary
star and the unseen primary companion in the system, sensibly limiting the
parameter space we can move in. Moreover, this work clearly shows the urgent
need for further and more detailed studies of this peculiar binary system,
NGC1850 BH1. They would be helpful to investigate a few important but still
unknown aspects: First, high resolution spectroscopy with a wide wavelength
coverage will be essential 1) to apply the disentangling technique in order to
be able to detect companions contributing as little as $\approx 1-2\%$ to the
visual flux of the system; 2) to study the properties of the luminous
(secondary) component (e.g. surface gravity, rotational velocity, chemical
abundances); 3) to assess whether the putative disk, if present, dilutes the
companion at all bands in a similar way; 4) to place unprecedented constraints
on the rejuvenation episodes that occur in binary systems when one of the two
sources gains a significant fraction of mass from the companion. The
rejuvenation factor is a very uncertain parameter and limiting its allowed
range would be a great achievement for future binary evolution studies.
Secondly, near-infrared high resolution photometry will be important to
investigate the detailed properties of the putative accretion disk in the
system, for example in terms of radius and temperature.

Those mentioned above are essential steps in deciphering the properties of the
unseen source in NGC1850 BH1.

## Acknowledgements

We thank the anonymous referee for the careful reading and analysis of the
paper. We thank Kareem El-Badry and Héloïse F. Stevance for the insightful
discussions on the system. We are also thankful to Matti Dorsch and Ulrich
Heber for their contributions to the analysis of the MUSE spectra. SS
acknowledges funding from STFC under the grant no. R276234. TS acknowledges
support from the European Union’s Horizon 2020 under the Marie Skłodowska-
Curie grant agreement No 101024605. SK acknowledges funding from UKRI in the
form of a Future Leaders Fellowship (grant no. MR/T022868/1). MG acknowledges
support from the Ministry of Science and Innovation (EUR2020-112157,
PID2021-125485NB-C22, CEX2019-000918-M funded by
MCIN/AEI/10.13039/501100011033) and from AGAUR (SGR-2021-01069). CU
acknowledges the support of the Swedish Research Council, Vetenskapsrådet.
This research has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(H.S., grant agreement 772225: MULTIPLES).

## Data Availability

The data underlying this work are already publicly available to the community.
The updated radial velocity measurements are instead listed in Table 1.

## References

  * Abdul-Masih et al. (2020) Abdul-Masih M., et al., 2020, Nature, 580, E11
  * Bagnuolo & Gies (1991) Bagnuolo William G. J., Gies D. R., 1991, ApJ, 376, 266
  * Beech (1989) Beech M., 1989, Ap&SS, 152, 329
  * Bodensteiner et al. (2020) Bodensteiner J., et al., 2020, A&A, 641, A43
  * Cappellari (2012) Cappellari M., 2012, pPXF: Penalized Pixel-Fitting stellar kinematics extraction, Astrophysics Source Code Library, record ascl:1210.002 (ascl:1210.002)
  * Cappellari & Emsellem (2004) Cappellari M., Emsellem E., 2004, PASP, 116, 138
  * Casares et al. (1992) Casares J., Charles P. A., Naylor T., 1992, Nature, 355, 614
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Davies et al. (2018) Davies R., et al., 2018, in Evans C. J., Simard L., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII. p. 1070209 (arXiv:1807.05089), doi:10.1117/12.2311480
  * Dsilva et al. (2020) Dsilva K., Shenar T., Sana H., Marchant P., 2020, A&A, 641, A26
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * El-Badry & Burdge (2022) El-Badry K., Burdge K. B., 2022, MNRAS, 511, 24
  * El-Badry & Quataert (2021) El-Badry K., Quataert E., 2021, MNRAS, 502, 3436
  * Eldridge & Stanway (2016) Eldridge J. J., Stanway E. R., 2016, MNRAS, 462, 3302
  * González & Levato (2006) González J. F., Levato H., 2006, A&A, 448, 283
  * Götberg et al. (2018) Götberg Y., de Mink S. E., Groh J. H., Kupfer T., Crowther P. A., Zapartas E., Renzo M., 2018, A&A, 615, A78
  * Hadrava (1995) Hadrava P., 1995, A&AS, 114, 393
  * Hubeny & Lanz (1995) Hubeny I., Lanz T., 1995, ApJ, 439, 875
  * Husser et al. (2016) Husser T.-O., et al., 2016, A&A, 588, A148
  * Kamann et al. (2023) Kamann S., et al., 2023, MNRAS, 518, 1505
  * Lanz & Hubeny (2003) Lanz T., Hubeny I., 2003, ApJS, 146, 417
  * Lanz & Hubeny (2007) Lanz T., Hubeny I., 2007, ApJS, 169, 83
  * Lattimer & Prakash (2001) Lattimer J. M., Prakash M., 2001, ApJ, 550, 426
  * Liu et al. (2019) Liu J., et al., 2019, Nature, 575, 618
  * Mahy et al. (2012) Mahy L., Gosset E., Sana H., Damerdji Y., De Becker M., Rauw G., Nitschelm C., 2012, A&A, 540, A97
  * Marchenko et al. (1998) Marchenko S. V., Moffat A. F. J., Eenens P. R. J., 1998, PASP, 110, 1416
  * Mennickent et al. (2003) Mennickent R., Pietrzyński G., Diaz M., Gieren W., 2003, A&A, 399, L47
  * Mennickent et al. (2015) Mennickent R. E., Djurašević G., Cabezas M., Cséki A., Rosales J. G., Niemczura E., Araya I., Curé M., 2015, MNRAS, 448, 1137
  * Orosz & Kuulkers (1999) Orosz J. A., Kuulkers E., 1999, MNRAS, 305, 132
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Rivinius et al. (2006) Rivinius T., Štefl S., Baade D., 2006, A&A, 459, 137
  * Rivinius et al. (2013) Rivinius T., Carciofi A. C., Martayan C., 2013, A&ARv, 21, 69
  * Rivinius et al. (2020) Rivinius T., Baade D., Hadrava P., Heida M., Klement R., 2020, A&A, 637, L3
  * Saracino et al. (2022) Saracino S., et al., 2022, MNRAS, 511, 2914
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shenar et al. (2017) Shenar T., et al., 2017, A&A, 606, A91
  * Shenar et al. (2020) Shenar T., et al., 2020, A&A, 639, L6
  * Shenar et al. (2022) Shenar T., et al., 2022, A&A, 665, A148
  * Stevance et al. (2022) Stevance H. F., Parsons S. G., Eldridge J. J., 2022, MNRAS, 511, L77
  * Wang et al. (2020) Wang C., Langer N., Schootemeijer A., Castro N., Adscheid S., Marchant P., Hastings B., 2020, ApJ, 888, L12
  * Zucker & Mazeh (1994) Zucker S., Mazeh T., 1994, ApJ, 420, 806


# Bijections between Variants of
Dyck Paths and Integer Compositions

Manosij Ghosh Dastidar TU Wien
Wien, AustriaInstitute of Discrete Mathematics and Geometry
<EMAIL_ADDRESS>TU Wien
Wien, AustriaInstitute of Discrete Mathematics and Geometry Michael Wallner TU
Wien
Wien, AustriaInstitute of Discrete Mathematics and Geometry
<EMAIL_ADDRESS>dmg.tuwien.ac.at/mwallner

###### Abstract

We give bijective results between several variants of lattice paths of length
$2n$ (or $2n-2$) and integer compositions of $n$, all enumerated by the
seemingly innocuous formula $4^{n-1}$. These associations lead us to make new
connections between these objects, such as congruence results.

Keywords:  Integer compositions, lattice paths, Dyck paths, bijections

##  1 Introduction

We explore several links between different variants of integer compositions
and generalizations of Dyck paths. Let us first introduce these objects.
First, an _integer composition_ of a nonnegative integer $n$ is a tuple
$(n_{1},\dots,n_{k})$ of nonnegative integers such that $n=n_{1}+\dots+n_{k}$.
Note that _integer partitions_ are integer compositions such that $n_{1}\geq
n_{2}\geq\dots\geq n_{k}$, (equivalently, the order of summands is not
significant). Second, a _Dyck path_ is a sequence of steps up
$\mathtt{u}=(1,1)$ and down $\mathtt{d}=(1,-1)$ that starts at the origin,
ends on the $x$-axis, and never crosses the $x$-axis. All classes of paths we
consider will start at the origin and consist of steps $\mathtt{u}$ and
$\mathtt{d}$, but the constraints will differ. Natural classes are _Dyck
walks_ that have no constraints, i.e., they may end anywhere and go below the
$x$-axis, and _Dyck bridges_ (also known as _grand Dyck paths_) that have to
end on the $x$-axis but may go below it.

Our results reveal a series of bijections, shown in Figure 1, connecting these
structures through a common enumeration formula

$4^{n-1}.$

This gives the corresponding integer sequence A000302 in the OEIS111The On-
Line Encyclopedia of Integer Sequences: http://oeis.org/ many new
combinatorial interpretations. Additionally, these bijections often map
natural statistics onto each other, such as the height of peaks and the number
of crossings of the $x$-axis.

ptPairs of compositions of $n$$3$-compositions of $n$[2]$2$-colored Dyck
bridges of length $2n-2$Unconstrained Dyck walks of length $2n-2$Left-to-right
max. in Dyck bridges of length $2n$Dyck paths with height-labelled peaks of
length $2n$Prop. 2.2Prop. 2.1Prop. 2.3Prop. 3.3Prop. 3.4

Figure 1: Bijections proved in this paper of classes of paths and integer
compositions that are all enumerated by $4^{n-1}$.

##  2 Bijections involving integer compositions

We start with the simple initial bijection, connecting pairs of compositions
and Dyck walks.

######  Proposition 2.1.

There exists a natural bijection between pairs of compositions of $n$ and Dyck
walks of length $2n-2$.

###### Proof.

Let a pair $(A,B)$ of two compositions of $n$ be given. First, we convert each
composition to a binary sequence: for each element $k$ in a composition,
append $k-1$ zeros followed by a $1$. By construction, both of these sequences
have to end in $1$. So we remove these ones and then concatenate the binary
sequences, with $A$’s sequence coming first. Finally, after replacing each $0$
by an up step $\mathtt{u}$ and each $1$ by a down step $\mathtt{d}$ the claim
follows. For the reverse direction cut the walk in the middle into two parts,
and re-add the ones. ∎

A $k$-composition is an integer compositions whose parts come in $k$ different
colors with the restriction that the last part of the composition is of the
first color; see [2]. We will consider only $3$-compositions.

######  Proposition 2.2.

There exists a natural bijection between $3$-compositions of $n$ and pairs of
compositions of $n$.

###### Proof.

By definition, the parts of $3$-compositions have three labels $1$, $2$, and
$3$. Anticipating the result, we introduce a notion of left and right: remove
the labels of color $1$, use label $L$ for color $2$, and label $R$ for color
$3$.

Now, we describe a map from $3$-compositions of $n$ to pairs of compositions
of $n$. First, we create two identical copies. In the first copy, we remove
the labels $R$ and add the parts labeled by $L$ to the next part. If the next
part has also a label $L$, then the addition continues to the next part, etc.
This gives a composition $A$ without any labels. Similarly, in the second
copy, we remove the labels $L$ and add the parts labeled by $R$ to the next
part. Again, if the next part has also a label $R$, then the addition
continues, and we get a composition $B$ without any labels. Observe that the
size of both compositions has not changed. Therefore, $(A,B)$ is a pair of
compositions of $n$.

To prove that this map is in fact a bijection, let us consider an arbitrary
pair $(A,B)$ of compositions of $n$. The key statistic to consider is the _run
of identical parts_ : Let $A=(a_{1},a_{2},\dots,a_{\ell_{A}})$ and
$B=(b_{1},b_{2},\dots,b_{\ell_{B}})$. A run is a sequence of maximal length
such that $a_{1}=b_{1}$, $a_{2}=b_{2}$, and so on. If $a_{1}\neq b_{2}$ we say
the run has length $0$.

For the inverse map, we will describe a recursive algorithm, which reduces the
sizes of $A$ and $B$ in the pair $(A,B)$ step-by-step and builds a
$3$-composition $C$. We start with an empty $3$-composition $C$. Depending on
the first parts of $A$ and $B$ we distinguish three cases:

  1. 1.

If $a_{1}=b_{1}$ then we attach $a_{1}$ with label $1$ to $C$, remove $a_{1}$
from $A$, and remove $b_{1}$ from $B$.

  2. 2.

If $a_{1}>b_{1}$ then we attach $b_{1}$ with color $2$ to $C$, replace $a_{1}$
in $A$ by $a_{1}-b_{1}$, and remove $b_{1}$ from $B$.

  3. 3.

If $a_{1}<b_{1}$ then we attach $a_{1}$ with color $3$ to $C$, remove $a_{1}$
from $A$, and replace $b_{1}$ in $B$ by $b_{1}-a_{1}$.

We repeat this process with the new values of $A$, $B$, and $C$. As the sizes
of $A$ and $B$ decrease in each step by at least one, this process terminates.
Moreover, note that in each step both parts decrease by the same size. Hence,
in the last step the process ends with case one where both parts are equal,
and therefore the final part gets label $1$, as required in the definition of
$3$-compositions. ∎

######  Example 2.3.

Consider the $3$-composition $6_{1}+1_{2}+4_{3}+2_{1}$ of $n=13$. First, we
remove the labels of color $1$, use label $L$ for color $2$, and label $R$ for
color $3$. Second, we create two identical copies. In the first copy, we
remove the labels $R$ and add the parts labeled by $L$ to the next part. If
the next part has also a label $L$, then the addition to the next part
continues, etc. Similarly, in the second copy, we remove the labels $L$ and
add the parts labeled by $R$ to the next part. This gives a pair of
compositions of $n$ without any labels, and we have shown in Proposition 2.2
that this is in fact a bijection.

$6_{1}+1_{2}+4_{3}+2_{1}$$(6+1_{L}+4_{R}+2)$,
$(6+1_{L}+4_{R}+2)$$(6+1_{L}+4+2)$, $(6+1+4_{R}+2)$$(6,1_{L}+4,2)$,
$(6,1,4_{R}+2)$$(6,5,2)$, $(6,1,6)$

##  3 Bijections involving Dyck paths

Let us now consider more complicated classes of Dyck paths. All of them use
the concept of a _peak_ , which is a consecutive pattern
$\mathtt{u}\mathtt{d}$. The first class we consider are _Dyck paths with a
marked peak_ , which are classical Dyck paths enriched by a marker on a
distinguished peak. Two such paths are different, if the underlying paths
differ, or, if the paths are the same then two different peaks are marked.
Therefore, the number of these paths is equal to the number of peaks in all
Dyck paths, whose enumeration is well-known; see, e.g., [4, Section 6.1].

######  Theorem 3.1.

There is an explicit bijection between Dyck paths with a marked peak of height
$h$ and Dyck bridges starting with a $\mathtt{d}$ step and $h-1$ crossings of
the $x$-axis preserving the length. Therefore, the number of peaks in all Dyck
paths of length $2n$ is equal to $\binom{2n-1}{n}$; see A001700.

###### Proof.

Let $D$ be a Dyck path with marked peak at height $h$. Using this peak, we
decompose the path $D$ into a left part $L$ from the origin to this peak and a
right part $R$ from this peak to the end: $D=LR$ such that $L$ ends with
$\mathtt{u}$ and $R$ starts with $\mathtt{d}$. In $L$ we perform a last-
passage decomposition, cutting at the $\mathtt{u}$ leaving a certain altitude
for the last time; while in $R$ we perform a first-passage decomposition,
cutting at the $\mathtt{d}$ bringing us down to a new altitude for the first
time; see Figure 2. More formally, we have

$\displaystyle L$
$\displaystyle=L_{1}\mathtt{u}L_{2}\mathtt{u}\dots\mathtt{u}L_{h}\mathtt{u},$
$\displaystyle R$
$\displaystyle=\mathtt{d}R_{h}\mathtt{d}R_{h-1}\mathtt{d}\dots\mathtt{d}R_{1},$
(1)

where $L_{i}$ and $R_{i}$ for $i=1,\dots,h$ are Dyck paths. Now, we pair the
paths $L_{i}\mathtt{u}$ and $\mathtt{d}R_{i}$ with the same index and map them
as follows to non-empty Dyck paths:

$\displaystyle D_{i}=\mathtt{u}L_{i}\mathtt{d}R_{i}.$

Then we concatenate these parts, after mapping each part with odd index to its
image obtained by any fixed bijection $\varphi$ between Dyck paths and
negative Dyck paths. This gives the Dyck bridge

$\displaystyle\varphi(D_{1})D_{2}\varphi(D_{3})D_{4}\dots\varphi(D_{h-1})D_{h}$
(2)

when $h$ is even. For odd $h$ it ends with $\varphi(D_{h})$. This bridge
starts with a down step $\mathtt{d}$ and crosses the $x$-axis $h-1$ times, as
claimed. Second, let a Dyck bridge starting with a $\mathtt{d}$ step be given.
We cut at each crossing of the $x$-axis and recover the components $D_{i}$ and
$\varphi(D_{i})$. Hence, it is straightforward to recover the components
$L_{i}$ and $R_{i}$ and to rebuild the Dyck path $D$ with marked peak.

Finally, bridges of length $2n$ are counted by $\binom{2n}{n}$, as there is an
equal number of up and down steps. Since half of them start with a down step,
we get

$\frac{1}{2}\binom{2n}{n}=\binom{2n-1}{n}.\qed$

xyLR

xy

Figure 2: A Dyck path with a marked peak (red dot) at height 6 and its image
under the bijection from Theorem 3.1 given by a Dyck bridge starting with a
$\mathtt{d}$ step and $5=6-1$ crossings (red dots). The black steps are used
in the last-passage (resp., first-passage) decomposition in the proof.

We return now to the bijections of Figure 1 and we connect our results with
_$2$ -colored bridges_; see [3, Section 6.4]. They are defined as the
concatenation of two bridges such that the first bridge is colored in color
$1$ and the second one in color $2$. Note that contrary to [3], we allow each
part to be empty. Hence, it is easy to see that its generating function is
equal to the square of the generating function
$B(z)=\frac{1}{\sqrt{1-4z^{2}}}$ of bridges:

$B(z)^{2}=\frac{1}{1-4z^{2}}.$

######  Proposition 3.2.

There is an explicit bijection between $2$-colored Dyck bridges and
unconstrained Dyck walks of the same length $2n$.

###### Proof.

Recall the classical notion of a Dyck meander, which is defined as the prefix
of a Dyck path, i.e., a path that starts at $0$, never goes below the
$x$-axis, but does not necessarily end on the $x$-axis. We will use repeatedly
that there is an explicit bijection between Dyck bridges and Dyck meanders of
the same length $2n$; see [7]. Sometimes, it will be necessary to transform a
meander further into a negative meander, by flipping all steps, i.e.,
exchanging $\mathtt{u}$ by $\mathtt{d}$ and vice versa.

We distinguish four cases. First, the first and second bridges are non-empty.
The idea is that the change in color corresponds to the last crossing of the
$x$-axis. For this purpose we transform the second bridge into a meander or
negative meander and attach it to the first bridge such that the attached
meander continues on the other side of the $x$-axis. We can easily reverse
this procedure by cutting at the last crossing of the $x$-axis. All the other
cases will have no crossings. Second, if the first bridge is non-empty and the
second one is empty, we transform the first bridge into a meander. Third, if
the first bridge is empty and the second on is non-empty, we transform the
second bridge into a negative meander. Finally, if both bridges are empty, we
map them to the empty walk. ∎

We continue, with _Dyck bridges with marked strict left-to-right maximum_. A
_strict left-to-right maximum_ is any peak $\mathtt{u}\mathtt{d}$ that has a
greater height than all peaks to its left. We called it marked in the previous
sense, when it is attached with a distinguished marker.

######  Proposition 3.3.

There is an explicit bijection between Dyck bridges of length $2n$ with marked
strict left-to-right maximum at height $h$ and $2$-colored Dyck bridges of
length $2n-2$ with $h-1$ crossings of the $x$-axis in color $1$.

###### Proof.

Let us start with a Dyck bridge with marked strict left-to-right maximum of
length $2n$. Then, we cut the bridge at the first return to the $x$-axis after
this maximum. The second part to the right is a bridge, which we give color
$2$. Onto the first part we apply a similar idea as in the bijection of
Theorem 3.1. As before, we cut the path at the marked left-to-right maximum
into a left and right part given by $LR$, such that $L$ ends with $\mathtt{u}$
and $R$ starts with $\mathtt{d}$. Now, we decompose it similar to (1) into

$\displaystyle L$
$\displaystyle=\varphi(L_{1})\mathtt{u}\varphi(L_{2})\mathtt{u}\dots\mathtt{u}\varphi(L_{h})\mathtt{u},$
$\displaystyle R$
$\displaystyle=\mathtt{d}R_{h}\mathtt{d}R_{h-1}\mathtt{d}\dots
R_{2}\mathtt{d},$

where $h$ is the height of the peak, and $L_{i}$ and $R_{i}$ are Dyck paths.
Note that in this case $R$ ends with a $\mathtt{d}$ step and contains only
$h-1$ Dyck paths $R_{i}$. As in the proof of Theorem 3.1 we form Dyck paths
$D_{i}=\mathtt{u}L_{i}\mathtt{d}R_{i}$ for $i=2,\dots,h$. Finally, we remove
the two steps of the marked peaks, and get the following bridges with two
steps less (compare with Equation (2)):

$\displaystyle\varphi(L_{1})D_{2}\varphi(D_{3})D_{4}\dots\varphi(D_{h-1})D_{h},$

when $h$ is even. For odd $h$ it ends with $\varphi(D_{h})$.

The mapping may be reversed by repeating the aforementioned steps in reverse
order. ∎

We end the bijections shown in Figure 1 by the following link to _Dyck paths
with height-labeled peak_ that are Dyck paths in which one peak is associated
with a label from $\\{1,2,\dots,h\\}$, where $h$ is the height of the specific
peak.

######  Proposition 3.4.

There is an explicit bijection between Dyck paths with height-labeled peak
with label $\mu$ at height $h$ and Dyck bridges of the same length with marked
strict left-to-right maximum at height $\mu$ and $h-\mu$ crossings of the
$x$-axis after this maximum.

###### Proof.

This bijection follows directly from the one described in the proof of Theorem
3.1, whose notation we will use here. The difference is that here we
concatenate the (positive and negative) Dyck paths differently; see Figure 3.

Let a Dyck path with height-labeled peak be given. Let $h$ be the height of
this peak and $\mu\in\\{1,\dots,h\\}$ be its label. First, we apply the
bijection $\varphi$ onto all parts $L_{i}$ in (1). From that we get the
following bridge in which the height-labeled peak is now a left-to-right
maximum (underlined):

$\displaystyle\varphi(L_{1})\mathtt{u}\varphi(L_{2})\mathtt{u}\dots\mathtt{u}\varphi(L_{h})\,\underline{\mathtt{u}\mathtt{d}}\,R_{h}\mathtt{d}R_{h-1}\mathtt{d}\dots\mathtt{d}R_{1}.$

Next, we transform this bridge, such that in the end the height-label $\mu$
constitutes the height of the left-to-right maximum. For this purpose, we
create and concatenate the paths $D_{i}$ and $\varphi(D_{i})$ in an
alternating fashion at the end:

$\displaystyle\varphi(L_{1})\mathtt{u}\varphi(L_{2})\mathtt{u}\dots\mathtt{u}\varphi(L_{\mu})\,\underline{\mathtt{u}\mathtt{d}}\,R_{\mu}\mathtt{d}R_{\mu-1}\mathtt{d}\dots\mathtt{d}R_{1}\,\varphi(D_{\mu+1})D_{\mu+2}\varphi(D_{\mu+3})\dots
D_{h},$ (3)

when $h-\mu$ is even. Otherwise, the last $D_{h}$ is replaced by
$\varphi(D_{h})$.

For the reverse direction, let a Dyck bridge with marked left-to-right maximum
be given. It is then straightforward to decompose it into (3) and to reverse
the steps above to build a Dyck path. The left-to-right maximum becomes the
height-labeled peak, labeled by its current height. Observe that the height-
labeled peak is lifted by the number of crossings of the $x$-axis to the left
of this peak. ∎

xyLR

xy

Figure 3: A Dyck path with a height-labeled peak (red dot) with label $4$ at
height $6$ and its image under the bijection from Proposition 3.4 given by a
Dyck bridge with marked strict left-to-right maximum (red dot) at height $4$
and $2=6-4$ crossings (red squares).

##  4 Conclusion and Outlook

Having established links between the internal structures of Dyck paths and
integer compositions, it is only natural to ask whether important theorems
from one subject can be transported to the other. When we talk about integer
compositions or partitions we are keen to see arithmetic properties in those
structures. In the long version of this work [5], we give further bijective
links and we show that such arithmetic results also exist in lattice paths. In
particular, we are pleased to note that:

######  Theorem 4.1 ([5, Theorem 3.8]).

Let $D_{r}(n)$ be the number of Dyck paths with semi-length $n$ and with
_exactly_ $r$ peaks for every reached height. Then $D_{r}(n)\equiv
0\pmod{r+1}$ for $n>r$.

In the opposite direction we also want to see if important theorems in integer
compositions and partitions can “generate” theorems in the world of lattice
paths. In 2020, Kim, Kim, and Lovejoy [6] observed the phenomenon of _parity
bias_ in partitions, where they showed that: if $p_{o}(n)$ denotes the number
of partitions of $n$ with more odd parts than even parts and if $p_{e}(n)$
denotes the number of partitions of $n$ with more even parts than odd parts,
then $p_{o}(n)>p_{e}(n)$.

In a subsequent article we will show that the analogous theorem is true even
for integer compositions. Furthermore, we will demonstrate that even for Dyck
paths a similar result holds, when segregating paths with respect to whether
they have more peaks at odd or even heights.

## References

  * [1]
  * [2] George E. Andrews (2007): _The Theory of Compositions, IV: Multicompositions._ The Mathematics student 76, p. 25.
  * [3] Cyril Banderier, Markus Kuba & Michael Wallner (2024): _Phase transitions of composition schemes: Mittag-Leffler and mixed Poisson distributions_.  To appear in Ann. Appl. Probab. arXiv:2103.03751.
  * [4] Emeric Deutsch (1999): _Dyck path enumeration_.  Discrete Math. 204(1-3), pp. 167–202, 10.1016/S0012-365X(98)00371-9.
  * [5] Manosij Ghosh Dastidar & Michael Wallner (2024): _Bijections and congruences involving lattice paths and integer compositions_.  arXiv:2402.17849.
  * [6] Byungchan Kim, Eunmi Kim & Jeremy Lovejoy (2020): _Parity bias in partitions_.  Eur. J. Comb. 89, p. 18, 10.1016/j.ejc.2020.103159.  Id/No 103159.
  * [7] Philippe Marchal (2003): _Constructing a sequence of random walks strongly converging to Brownian motion_.  In: Proceedings of Discrete Random Walks (DRW’03), Discr. Math. Theor. Comput. Sci., pp. 181–190, 10.46298/dmtcs.3335.