11institutetext: Nicolaus Copernicus Astronomical Centre, Polish Academy of
Sciences, ul. Bartycka 18, 00-716 Warsaw, Poland 22institutetext: Fesenkov
Astrophysical Institute, 050020, Almaty, Kazakhstan 33institutetext: Zentrum
für Astronomie der Universität Heidelberg, Astronomisches Rechen-Institut,
Mönchhofstr 12-14, 69120 Heidelberg, Germany 44institutetext: Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA
55institutetext: Konkoly Observatory, HUN-REN Research Centre for Astronomy
and Earth Sciences, Konkoly Thege Miklós út 15-17, 1121 Budapest, Hungary

# Star-by-star dynamical evolution of the physical pair of the Collinder 135
and UBC 7 open clusters

Maryna Ishchenko<EMAIL_ADDRESS>Dana A. Kovaleva  Peter Berczik
11223355 Nina V. Kharchenko  Anatoly E. Piskunov  Evgeny Polyachenko  3311
Ekaterina Postnikova  Andreas Just  33 Olga Borodina  44 Chingis Omarov  22
Olexandr Sobodar  11

(Received … / Accepted …)

###### Abstract

Context. In a previous paper using Gaia DR2 data, we demonstrated that the two
closely situated open clusters Collinder 135 and UBC 7 might have formed
together about 50 Myr ago.

Aims. In this work, we performed star-by-star dynamical modelling of the
evolution of the open clusters Collinder 135 and UBC 7 from their supposed
initial state to their present-day state, reproducing observational
distributions of members.

Methods. Modelling of the Collinder 135 and UBC 7 dynamical evolution was done
using the high-order parallel $N$-body code $\varphi-$GPU with up-to-date
stellar evolution. Membership and characteristics of the clusters were
acquired based on Gaia DR3 data.

Results. The comparison of the present-day radial cumulative star count
obtained from the $N$-body simulations with the current observational data
gave us full consistency of the model with observational data, especially in
the central 8 pc, where 80% of the stars reside. The proper motion velocity
components obtained from the $N$-body simulations of the stars are also quite
consistent with the observed distributions and error bars.

Conclusions. These results show that our numerical modelling is able to
reproduce the open clusters’ current complex 6D observed phase-space
distributions with a high level of confidence. Thus, the model demonstrates
that the hypothesis of a common origin of Collinder 135 and UBC 7 complies
with present-day observational data.

###### Key Words.:

open clusters and associations: individual – numerical integration in Galactic
potential: open cluster orbits – binary open star cluster evolution: initial
conditions for open clusters

##  1 Introduction

The data of the space mission Gaia (Gaia Collaboration et al., 2016) have
stimulated investigations into open clusters (hereafter OCs) in the Galaxy.
One of the areas of research that has significantly benefited is related to
the binarity of OCs. While a lack of observational data of binary star
clusters in our Galaxy in comparison with the Magellanic Clouds has been
discussed previously (e.g. Vázquez et al., 2010; de La Fuente Marcos & de La
Fuente Marcos, 2009), recent investigations have provided a number of examples
of pairs of OCs in our Galaxy that suggest a common origin (see, e.g. Qin et
al., 2023; Angelo et al., 2022; Ye et al., 2022; Casado, 2022a; Bisht et al.,
2021; Zhong et al., 2019). Some clusters that were previously considered
single have shown a binary signature upon detailed examination (Alejo et al.,
2020; Dias et al., 2018), including clusters with signs of merging by tidal
capture (Camargo, 2021), according to Gaia data.

There is observational evidence that indicates binary star clusters are
usually coeval and form as a result of hierarchical clustering in parental
molecular clouds (Casado, 2022b; Angelo et al., 2022). These pairs may or may
not be bound and physically interacting. Notably, Angelo et al. (2022) found
four bound and two interacting but unbound pairs of star clusters.
Alternatively, coeval clusters may be located relatively close to each other
due to being born in the same star forming complex but may never interact
(Casado, 2022b, a; Casado & Hendy, 2023), and eventually they may become
disaggregated while maintaining similar Galactic orbits (Coronado et al.,
2022). On the other hand, pairs of closely located OCs with significantly
different ages, which represent bona fide events of an occasional encounter,
have been found (Camargo, 2021; Piatti & Malhan, 2022; Olivares et al., 2023).
Dynamical simulation of such encounters has so far been limited to the
publicly available orbit integration code galpy (Bovy, 2015). The dynamical
simulations of the evolution of initially bound primordial OCs in the Galactic
tidal field by de la Fuente Marcos & de la Fuente Marcos (2010) have
demonstrated that very close pairs, if formed, are merged in a short
timescale, and in other pairs, either the smaller component is disrupted or a
pair is decoupled in the background tidal field due to mass loss, making open
star cluster binarity a transient phenomenon. Priyatikanto et al. (2016) used
$N$-body simulations to explore the evolution of binary open star clusters in
the Galactic tidal field with different initial orbital orientations with
respect to the Galactic plane, which proved to be another important factor
influencing the dynamical fate of such clusters. Darma et al. (2021) carried
out $N$-body simulations to investigate the formation and evolution of binary
star clusters in the Milky Way and in the Large Magellanic Cloud and found
that binary open star clusters can form from stellar aggregates with a variety
of initial conditions.

In Kovaleva et al. (2020), hereafter Paper I, we discussed possible scenarios
of origin and dynamical evolution in a simplified approach and constrained
parameters, within the limits of observational values, that suggest a common
origin of the open star clusters Collinder 135 (hereafter Cr 135) and UBC 7
(see also discussion in Pang et al., 2021). As one of the results of that
work, we identified, using backward orbital integration, a presumable grouping
of initial states in 6D space that corresponds to a common origin of these
clusters. Our previous work was limited only by the investigation of motion of
OC dynamical centres (centre of masses). So, we completely neglected the OC
internal dynamical and stellar evolution, including the mass loss of the
clusters. By overcoming such restrictions in our current investigation, we
make our previous results more robust and reliable.

In this article, our purpose is to reproduce the dynamical evolution of the
system of Cr 135 and UBC 7 using star-by-star simulations. During our
investigation, we especially paid attention to the complex dynamical co-
evolution and cross-penetration of the clusters’ stellar populations. This
behaviour of the modelled stellar systems made our attempts of system
modelling quite challenging.

In comparison with Paper I, we have updated the numerical implementation of
the $\varphi-$GPU with a new stellar evolution block (Banerjee et al., 2020;
Kamlah et al., 2022) combined with our high-order dynamical $N$-body code.
This new implementation allowed us to predict with a good degree of confidence
the evolved clusters’ internal stellar structure and cumulative mass
distribution.

After the publication of Paper I, the Gaia DR3 catalogue (Gaia Collaboration
et al., 2021, 2023) was published. With respect to the previously used Gaia
DR2 catalogue (Gaia Collaboration et al., 2018b), Gaia DR3 provides enhanced
astrometric and photometric data for an enlarged sample of the stars as well
as a significantly larger number of sources with mean radial velocities. Some
systematic effects have been resolved or mitigated regarding Gaia DR2, which
is also beneficial for the study of the population of OCs. These improvements
allowed us to obtain refined characteristics of the system composed of the two
clusters based on the deep census of the probable members of this system.

Our goal is to reproduce the present-day state of OCs and compare it with the
most up-to-date information about their stellar population from Gaia DR3 data.
Such an ambitious investigation has only been possible since last few years
thanks to the Gaia’s unprecedented position and proper motion measurements. As
a possible physical pair of OCs, Cr 135 and UBC 7 give us a unique laboratory
to also study the common processes of the ongoing star formation on larger
scales in our Milky Way disc.

The paper is organised as follows. In Sec. 2, we describe dynamical $N$-body
modelling of evolution of the pair of clusters. In Section 3, we discuss the
data filtering process that resulted in the following view onto the system of
the two OCs, Cr 135 and UBC 7, with Gaia DR3 data. Section 4 is focused on
characterisation of the system based on the obtained characteristics of its
members. In Section 5, we describe results obtained in our $N$-body
simulations in comparison with the observational data presented in Sections 3
and 4. In Section 5, we also discuss kinematics of the surrounding extended
halo. In Section 6, we sum up the main results of the work.

##  2 Cr 135 and UBC 7: Star-by-star N-body modelling

###  2.1 Initial conditions and integration procedure

In Paper I we found that it is highly probable that Cr 135 and UBC 7 form a
physical pair. This conclusion is based on backward-in-time numerical
simulations and 6D observations from the Gaia DR2 catalogue. The resulting
initial positions and velocities for cluster centres of mass (initial
separation and velocity: 10 pc and 0.95 km/s) showed that the probability of a
mutual formation scenario is about 98% (see the numerical model #[53,61], and
Fig. 2 of Paper I). In order to estimate the impact of Gaia DR3 data on new
simulations, we used newly determined present-day cluster observed parameters
determined as described in detail in Sec. 3. We found that there is good
agreement in the newly determined Galactic coordinates, proper motions, and
radial velocities (see Table 1) with the Gaia DR2 parameters derived in Paper
I).

Table 1: Comparison of the main parameters of the OCs based on two Gaia data
releases.

DR | OC |  $l$, |  $b$, |  $\varpi$, |  $\mu_{l}^{*}$, |  $\mu_{b}$, |  $V_{r}$,
---|---|---|---|---|---|---|---
|  | deg | deg | mas | mas/yr | mas/yr | km/s
3 | Cr 135 | 248.75 |  $-$11.10 | 3.38 |  $-$9.99 |  $-$6.40 | 15.2
|  |  $\pm$0.06 |  $\pm$0.04 |  $\pm$0.03 |  $\pm$0.03 |  $\pm$0.03 |  $\pm$1.5
2 | Cr 135 | 248.98 |  $-$11.10 | 3.31 |  $-$9.92 |  $-$6.47 | 17.4
|  |  $\pm$0.06 |  $\pm$0.05 |  $\pm$ 0.02 |  $\pm$ 0.05 |  $\pm$ 0.06 |  $\pm$ 1.3
3 | UBC 7 | 248.62 |  $-$13.45 | 3.61 |  $-$10.38 |  $-$5.90 | 15.8
|  |  $\pm$0.04 |  $\pm$0.05 |  $\pm$0.03 |  $\pm$ 0.03 |  $\pm$ 0.03 |  $\pm$ 2.8
2 | UBC 7 | 248.62 |  $-$13.37 | 3.56 |  $-$10.25 |  $-$5.98 | 16.7
|  |  $\pm$0.04 |  $\pm$ 0.05 |  $\pm$ 0.02 |  $\pm$0.05 |  $\pm$ 0.05 |  $\pm$ 1.5

∗ – here and further $\mu_{l}\equiv\mu_{l}^{0}\cdot\cos b$

To evaluate the possible distinction (DR2 vs. DR3) in the distances between
the two OC centres, we carried out a new lookback integration using the new
DR3 initial data from Table 1 and compared it to our earlier selected best
numerical model #[53,61] (see Fig. 1). As can be seen in the figure, the
agreement is good, and the effect of the different data releases is small.

Figure 1: Evolution of the distance between the centres of Cr 135 and UBC 7
for the DR2 (red) and DR3 (blue) initial conditions.

For the OCs’ dynamical orbital integration including the stellar evolution, we
used the high-order parallel $N$-body code $\varphi-$GPU ,111$N$-body code
$\varphi-$GPU :
https://github.com/berczik/phi-GPU-mole. which is based on the fourth-order
Hermite integration scheme with individual hierarchical block time steps
(Berczik et al., 2011, 2013). The current version of the $\varphi-$GPU code
uses a native GPU support and a direct code access to the GPU using the NVIDIA
native CUDA library. The present code is well tested and was used to obtain
important results in our earlier OC simulations Just et al. (2009);
Shukirgaliyev et al. (2017, 2018, 2021a).

Because of the relatively short evolutionary timescale of our system, we used
the so-called fixed external potential with parameters independent of time.
Similar to Kharchenko et al. (2009); Ishchenko et al. (2023), we used the
three-component (bulge-disc-halo) axisymmetric Plummer-Kuzmin model (Miyamoto
& Nagai, 1975a).

The current code also takes into account up-to-date stellar evolution models.
The most important updates were made to components of stellar evolution and
include updated metallicity dependent stellar winds; updated metallicity
dependent core-collapse supernovae, along with a new fallback prescription and
their remnant masses; updated electron-capture supernovae accretion induced
collapse and merger-induced collapse; remnant masses and natal kicks; and
black hole natal spins, among other updates. (For more details about the new
stellar evolution library, see Kamlah et al. (2022)).

To carry out the star-by-star $N$-body modelling together with the stellar
evolution for both OCs, we assumed the clusters to be in initial dynamical
equilibrium. For the star particles, we used the King model distribution
function King (1966) with three free parameters: the initial total mass, the
initial half-mass radius $r_{\rm hm}$, and the dimensionless central potential
$W_{0}$. Using the Kroupa (Kroupa, 2001) initial mass function with the lower
and upper mass limits equal to 0.1 – 10 $\mathrm{M}_{\odot}$, and the
classical solar value for the initial stellar metallicity, $Z=Z_{\odot}=0.02$
(Grevesse & Sauval, 1998), for both clusters. For the centres of the clusters,
we used the initial coordinates and velocities at the 50 Myr lookback time
integration (see Table 2).

We tested more than 50 runs with various combinations of $M,r_{\rm hm}$, and
$W_{0}$ separately for each cluster using the full $N$-body model calculations
up to the present time. For the initial positions and orbital velocities in
Cartesian Galactic coordinates of the star clusters (SC) centre at the -50 Myr
lookback in time, we used the values obtained during our own set of lookback
integration of the SC centre of mass. The initial masses, radii, and King
profiles were found during the fitting test runs, with the aim of matching the
theoretical and observed parameters. The basic initial parameters of the model
are given in Table 2 separately for each SC.

Table 2: Initial positions, velocities, and parameters for OCs at 50 Myr lookback time. SC |  $X$, pc |  $Y$, pc |  $Z$, pc |  $V_{x}$, km/s |  $V_{y}$, km/s |  $V_{z}$, km/s |  $M$, $\mathrm{M}_{\odot}$ | $N$ |  $r_{\rm hm}$, pc | $W_{0}$
---|---|---|---|---|---|---|---|---|---|---
Cr 135 | –1073.5 | –8415.3 | –20.4 | –230.17 | 34.14 | –5.38 | 230 | 442 | 10 | 3
UBC 7 | –1076.0 | –8417.7 | –12.0 | –230.31 | 34.96 | –5.76 | 200 | 384 | 7 | 11

Figure 2 shows the observed stars of the OCs (light green for Cr 135 and light
blue for UBC 7) overplotted with the present-day positions of all particles
from our basic model according to the initial data in Table 2 (dark green for
Cr 135 and dark blue for UBC 7). The mutual orientation of the OCs is correct,
although both sets of particles (representing Cr 135 and UBC 7) are slightly
misaligned with their prototype clusters (see Table 3). In Table 3, we present
the positions and velocities for the current centres of mass for both OCs. The
identifier Sim indicates the cluster’s centre of mass as obtained from our
basic model, and Obs indicates the centre of mass as obtained from
observations (see Section 3). The small disagreement between the 3D
coordinates and velocities inside the table most probably arises from the
different mutual interactions of the clusters in the backward integration of
the orbit compared to the forward integration of the full OCs. The backward
integration was performed with the present-day masses of the OCs. The forward
integration started with the larger initial cluster masses and includes mass
loss by stellar evolution and the tidal field. The present-day shift of the
cluster centres demonstrates that Cr 135 and UBC 7 interact gravitationally
with each other.

Figure 2: Comparison of observed positions of probable cluster members and of model stars of basic $N$-body simulations. The light green and light blue dots are the members of Cr 135 and UBC 7, respectively. The coordinate frame was obtained by parallel translation of the standard Galactic coordinate system. Table 3: Comparison of the observed and numerical centre-of-mass positions and velocities for Cr 135 and UBC 7 in Cartesian Galactic coordinates for the present day. SC |  |  $X$, pc |  $Y$, pc |  $Z$, pc |  $V_{x}$, km/s |  $V_{y}$, km/s |  $V_{z}$, km/s
---|---|---|---|---|---|---|---
Cr 135 | Sim | –8281.3 | –277.0 | –36.2 | –7.63 | 237.35 | –5.11
| Obs | –8282.1 | –271.0 | –36.2 | –7.38 | 237.94 | –4.85
UBC 7 | Sim | –8283.1 | –239.5 | –45.2 | –7.10 | 237.71 | –4.72
| Obs | –8276.2 | –250.8 | –43.6 | –7.00 | 238.13 | –4.27

###  2.2 Stellar evolution and lifetime

In Fig. 3, we present the evolution of the tidal radius and bound particle
numbers according to our basic $N$-body model. The radii of both clusters show
a slight increase during the first several tens of million years, after which
the radius of Cr 135 starts to sharply decrease, and after approximately 110
million years of evolution, the core falls apart (technically, the number of
bound particles decreases to zero). At this time, the separation between the
two OCs is about 50 pc. UBC 7 has a more gradual decrease in tidal radius in
contrast to Cr 135, and it falls apart approximately 50 million years after
the dissolution of Cr 135. The differences between the dynamical behaviour of
our clusters can be explained as being a result of the concentration
parameters $W_{0}$ = 3 and $W_{0}$ = 11 of the different initial King
profiles. Thus, for Cr 135 the profile indicates a looser and more unstable
core. In UBC 7, the core is more concentrated, and therefore, despite the
smaller number of stars, it demonstrates a greater ability to survive.

Figure 3: Evolution of the tidal radius (upper panel) and bound particle
number (bottom panel). The green line represents Cr 135, and the blue line
indicates UBC 7.

In Fig. 4, we show the orbital evolution up to $\sim$150 Myr from the
formation of the OCs. As can be seen from Fig. 4, during the time of
evolution, the clusters follow mildly eccentric clockwise orbits around the
Galactic centre, as seen from the Galactic north pole, and simultaneously
rotate around each other. The orbits have considerable vertical oscillations
in the $z$-direction between $\pm 80$ pc.

Figure 4: Orbital evolution of Cr 135 and UBC 7. The colour-coding represents
the time of integration. The filled and open black circles indicate the
present positions of the clusters. Th black crosses indicate the positions
where the clusters dissolve.

In Fig. 5, we present the Hertzsprung-Russell diagram (HRD) for three moments
of evolution (50, 100, and 150 Myr in forward integration) for each OC
starting from 0 Myr. Due to the low number of available stars and the short
period of the OCs’ lifetime, we mainly have main-sequence stars. When
analysing the type of stars for the present day (which have an age of 50 Myr
according to the $N$-body simulation), $\sim 84\%$ of them are low-mass stars
(m ¡ 0.8 M⊙), and only $\sim 16\%$ of them are higher-mass main-sequence stars
for both OCs. At the end of the forward integration, low-mass stars were still
dominant, and there were only a few carbon/oxygen and oxygen/neon white dwarf
stars.

Figure 5: Hertzsprung-Russell diagrams of Cr 135 and UBC 7 (including all
unbound stars). Red dots correspond to an age of 0 Myr, light blue to 50 Myr,
green to 100 Myr, and black to 150 Myr. The crosses and the stars show two
types of white dwarfs from the $N$-body simulation.

##  3 Cluster membership with Gaia DR3 data

Our work in Paper I is based on Gaia DR2 data and aimed at selecting the most
probable members of Cr 135 and UBC 7 in order to obtain their global space and
kinematic parameters. We incidentally discovered an extended corona of stars
with an almost constant projected density inside a 20-pc radius around the
centre of mass of both clusters. In the present investigation, our aim is to
use the wealth of the Gaia DR3 data to explore stellar membership and the
structure of the whole system in detail.

Recent investigations have proved the existence of extended features connected
with groups of OC (see, for example, Beccari et al. (2020); Zhong et al.
(2019); Pang et al. (2020)). On the other hand, the discussed Cr 135 and UBC 7
are situated in the Vela-Puppis star-formation region (Galactic longitude
$245\degr\lesssim l\lesssim 265\degr$ and latitude $-15^{\circ}\lesssim
b\lesssim-5^{\circ}$), which hosts several young OCs with ages between
approximately 15 and 50 Myr (Cantat-Gaudin et al., 2019b, a; Beccari et al.,
2020). Of these clusters, Cr 140, NGC 2451B, and NGC 2547 are situated at
about 100 pc from the investigated pair. The nearest neighbour to Cr 135 and
UBC 7 is the evolved open cluster Alessi 3 (625 Myr, Hunt & Reffert, 2023).
Thus, in order to reveal the size and characteristics of the coronae and check
whether it is due to contamination from neighbouring stellar groups, we
considered a very wide region surrounding Cr 135 and UBC 7. This required
refining of the data selection principles used in Paper I.

We performed an intensive search for members related to OCs Cr 135 and UBC 7
among all Gaia DR3 sources satisfying the following conditions:

$2.5<\varpi\ {\rm[mas]}<5$,

$220\degr<l<280\degr$,

$-25\degr<b<1\degr$,

$1.2<\mu_{\rm\alpha}^{*}\ {\rm[mas/yr]}<11.8$, and

$-15<\mu_{\rm\delta}\ {\rm[mas/yr]}<-5$.

The linear size of a region was determined based on previous research in order
to include regions that are most probably free of members of central clusters.
These external regions are useful for the estimation of random background
contamination from the procedure. We applied parallax zero-point correction
(Lindegren et al., 2021) with the algorithm provided by Gaia.222 Gaia EDR3
source code from
https://www.cosmos.esa.int/web/gaia/edr3-code.

Figure 6 presents the primary selection containing 116 975 sources. Their
colour-coded parallax ranges from yellow, for the most distant stars of the
sample at about 400 pc from the observer, to dark violet, for foreground stars
at 250 pc from the Sun. Even with this rough selection, some of the OCs of a
region became noticeable, Cr 135 and UBC 7 being the nearest two. At a
distance of about 280 to 300 pc from the Sun where they are situated, 1 deg
$\approx$ 5 pc, so the linear size of the rectangle in the frontal plane is
about 300 pc in $l$ and about 75 pc in $b$ (Fig.6, left panel). One can see an
enhanced density of sources towards the Galactic plane in $b$ and towards the
Vela Molecular Ridge and the Galactic centre in $l$. As we applied no filter
to the quality of astrometric and photometric solutions, part of the faint
stars may belong to the background and may reflect the stellar density of
distant regions (Lindegren et al., 2018).

Figure 6: A selection of Gaia DR3 data in the broad vicinity of the Cr 135 and
UBC 7. Cr 135, UBC 7, and their closest neighbouring clusters according to
Cantat-Gaudin et al. (2020), left panel. Refined sample of $60\%$ of the
members of Cr 135 and UBC 7 over the same sky area, right panel. The colour-
coding indicates the MP according to Eq. 1. The circles in both panels
represent the half-number radii of the clusters according to Cantat-Gaudin et
al. (2020). See discussion in the text.

We kept as a basis the procedure of probabilistic selection of cluster members
described by Kharchenko et al. (2012). Each star of the considered dataset was
assigned a probability $P$ of sharing age, space, and kinematics of a pair of
clusters, namely Cr 135 and UBC 7:

$P=\min(P_{\rm kin},P_{\rm\varpi},P_{\rm ph}).$ (1)

where $P_{\rm kin},P_{\rm\varpi},P_{\rm ph}$ are the probabilities based on
kinematics, parallax, and photometric data. The procedure of calculating the
probabilities is described in detail in Appendix A.

After each source in our dataset was assigned a membership probability (MP)
according to Eq. (1), we selected stars with $P>0.6$ (so-called
$60\%$-members) for further consideration. In Fig. 6 (right panel), one can
observe two dense central groups (i.e. the cores of Cr 135 and UBC 7) and an
extended group of stars with a high MP resembling a halo surrounding the
cores. Also present are sparse, chaotically distributed sources of lower
probability that dominate the periphery, which presumably are contaminating
foreground and background stars, and an elongated dense structure between
$250\degr<l<263\degr$ and $-11\degr<b<-7\degr$ formed by members with an
intermediate MP. After investigation, we found that this elongated stellar
group has a systematically different space motion and does not belong to the
system we are focused on. This group may in part be associated with a group of
stars identified as a loosely spread cluster [BBJ2018] 6 in Beccari et al.
(2018), but due to its spatial scale, it might be referred to as a stream
instead. For the purposes of the present work, we identified and removed from
consideration sources with MPs to the [BBJ2018] 6 group larger than MPs to the
group of Cr 135 and UBC 7. The remaining dataset represented in Fig. 7
consists of 1714 stars. One may observe that the probable members of the group
composed of the two clusters still form a distinct asymmetrical halo around
the location of the dense clusters.

Random contamination can be estimated based on the number of probable members
in the outer parts of the field, which are approximately at $l>270\degr$,
$l<235\degr$, $b>-4\degr$, $b<-20\degr$. For this purpose, as well as to set a
distinction between the halo and the background, we relied on the density
distribution of the stars in the tangential plane. For each source, we
estimated its weight over the background $W_{i}$ as the residual:

$W_{i}=\rho(l_{i},b_{i})-\rho_{\rm r}(l_{i},b_{i}),$ (2)

where $\rho(l,b)$ is the stellar density distribution of the refined dataset
convolved by the 2D Gaussian kernel of the width of 2°and $\rho_{\rm r}(l,b)$
is the similarly convolved random flat distribution of the same number of
stars over the same field. We considered the 293 sources with $W_{i}\leq 0$ to
be background contamination stars, and we treated the 1421 stars with a
positive $W_{i}$ as those composing the clusters and the halo. From hereon, we
refer to this dataset as Sample I (available via CDS).

Our purpose was to recover the most probable distribution of sources in Sample
I between three substructures: Cr 135, UBC 7, and their extended halo. In
space, these groups appear to be separated, but there is no distinct border
between them. Given that the separation between the centres of the clusters is
just 22 pc, peripheral stars of the clusters are mixed and may only be
conditionally attributed to one of them. Due to the similarity of their mean
kinematic characteristics, we attributed each source to one of these
components based solely on its maximum likelihood location in space with
respect to the clusters as follows.

Using the mean position of the reliable members from Table 1, we suggested for
the clusters a Gaussian 3D PDF $P^{sc}_{k}$ over a Cartesian distance to the
centre with coordinates based on Gaia DR3
($l^{0}_{k},b^{0}_{k},\varpi^{0}_{k}$) from Table 1, where $k=1$ refers to Cr
135 and $k=2$ is for UBC 7:

$\Pi(d_{k})=\exp\left[\frac{\left(d^{i}_{k}\right)^{2}}{2\sigma_{d^{k}}^{2}}\right]\,,$
(3)

where $\sigma_{d^{k}}$ is based on the distribution of the reliable cluster
members by Gaia DR3 (Table 1). We considered the distance $d$ between the
centre of an OC ($l^{0}_{k},b^{0}_{k},\varpi^{0}_{k}$) and the $i^{th}$ source
of Sample I ($l^{i},b^{i},\varpi^{i}$) to be a function of $\varpi^{i}_{m}$
with a Gaussian PDF:

$\Pi(\varpi^{i})=\exp\left[-\frac{\left(\varpi^{i}-\varpi^{i}_{m}\right)^{2}}{2\sigma^{2}_{\varpi^{i}}}\right]\,,$
(4)

where the standard deviation $\sigma^{i}_{\varpi}\equiv\epsilon_{\varpi}^{i}$.
We simulated numerically the resulting PDF for the distance of each source of
Sample I to the cluster centre
$\Pi(d_{k},\varpi^{i})=\Pi(d_{k})\cdot\Pi(\varpi^{i}_{k})$, and we compared
the maximum values of PDFs for $k=1,2$:

if $\max(\Pi(d_{1},\varpi^{i}))\geq\max(\Pi(d_{2},\varpi^{i}))\land
d_{1}^{i}[\max(\Pi(d_{1},\varpi^{i}))]\leq 15$ pc, then

the source $i\in$ Cr 135,

else

if $\max(\Pi(d_{2},\varpi^{i}))\geq\max(\Pi(d_{1},\varpi^{i}))\land
d_{2}^{i}(\max(\Pi(d_{2},\varpi^{i})))\leq 15$ pc, then

the source $i\in$ UBC 7,

else

the source $i\in$ halo.

end if

end if

The border between the cluster and halo members was selected to be equal to 15
pc, tentatively, based on the results of the analysis of the spatial
distribution of Sample I. As a result, we related 294 stars to Cr 135, 243
stars to UBC 7, and 884 stars to the extended halo (see Fig. 7).

##  4 Estimation of the parameters of the clusters

The system includes two quasi-ellipsoidal dense stellar groups that we refer
to as clusters, which have similar but distinct mean space motions. They are
surrounded by an extended stellar halo of irregular form and have an
approximately constant stellar density in projection onto the $lb$-plane (Fig.
7, upper panel). This structure extends for more than 80 pc towards the
Galactic anti-centre and half as much in the opposite direction. The
dimensions towards and from the Galactic plane are comparable with the size of
the clusters. In space, the structure is aligned along the direction joining
the centres of Cr 135 and UBC 7, becoming closer to the Galactic plane the
farther away it is from the Sun (see Fig. 7, bottom panel).

The similarity of the populations of two clusters and of the population of the
extended halo is illustrated by their distribution over the colour-magnitude
diagram in Fig. 8. In the figure, apparent magnitudes and reddened colours
have been reduced to the absolute scale using individual parallaxes of the
stars and mean reddening of the central probable members as in Paper I. The
Green line represents the Padova isochrone for 50 Myr for single stars, and
the magenta line is for the same isochrone for the unresolved binary stars
with identical components. The red line is for the zero-age main sequence
built in Paper I as a hot envelope of the related set of Padova isochrones of
different ages. The probable members of the components of the system are
represented by the same colours as in Fig. 7. For reference, grey dots
represent each tenth background source of the same spatial region (as in Fig.
6). One may notice not only the similar age of the OCs and the halo but also
the similar density distribution of probable members along the isochrone,
which suggests a common origin.

The characteristics of the components of a system composed of the two OCs and
their halo based on Sample I data selection are listed in Table 4. This
includes the number of probable members as defined in Section 3, $N$, and the
half-number radii of the clusters and halo, $R_{HN}$.

Figure 7: Components of the OC system. The probable members of Cr 135 are
shown with red dots, probable members of UBC 7 are indicated with blue dots,
and extended halo members are shown with yellow dots. The upper panel shows
the $lb$ Galactic coordinates sky plane; the bottom panel shows the
$\varpi,\,b$ plane. Figure 8:  Components of the OC system at the colour-
magnitude diagram. Probable members of Cr 135 are shown with red dots,
probable members of UBC 7 are indicated with blue dots (left panel), and
probable members of the extended halo are shown with yellow dots (right
panel). Grey dots are for background stars. The green line represents the
isochrone for 50 Myr for single stars, and the magenta line is for the same
isochrone for unresolved binary stars with identical components. The red line
is for the zero-age main sequence.

One may expect that the halo contains significant contamination and a random
background. Its numerical estimate based on the surface density of a random
background in the neighbouring regions led to $0.25\pm 0.05$ contaminating
stars per square degree. This decreased the expected actual population of the
halo by 60 to 120 stars, taking into account the vagueness of its borders.

###  4.1 Counted masses

The estimate of the present-day mass of the system may be assessed based on
masses of visible probable members as $\approx 900M_{\odot}$ (see Table 4,
$M^{\Sigma}$). More than half of the mass lies in the halo. There are a few
probable causes of it being underestimated, such as hidden mass in low-mass
stars beyond the photometric limit, members with an MP less than $60\%$, and
unresolved binarity. As well, there is likely random contamination that causes
present-day mass overestimation. The effect of random contamination is mostly
important for the halo due to its large extent. We provide a numerical
estimate for the expected contamination based on the surface density of the
random background in the neighbouring regions as obtained in Sec. 4. Filtering
out the expected background reduced the mass of the halo by 40 to
$70M_{\odot}$. For the clusters, the halo itself should be considered a
contaminating background; the estimate of the mean density for it is $3.6\pm
0.2$ stars per square degree. This provides an upper estimate for a
contaminating mass of $\approx 27M_{\odot}$ for Cr 135 and $\approx
17M_{\odot}$ for UBC 7.

A possible source of mass underestimation lies within the unresolved stellar
multiplicity and, primarily, binarity (see discussion in Borodina et al.,
2019, 2021). The effect depends on the fraction of binary to multiple stars
and their distribution over the mass fraction $q=m_{B}/m_{A}$, where $m_{A}$
and $m_{B}$ are respectively masses of the most massive and least massive
components of a binary star.We estimated the fraction of binary stars
$\alpha=N_{bin}/N_{tot}$ with mass ratio of components $q>0.4$ for the system
directly from the colour-magnitude diagram (Borodina et al., 2019; Borodina &
Kovaleva, 2020) as $\alpha=23\%\pm 12\%$. The result is consistent with the
binary fraction estimates carried out by Pang et al. (2023). The estimated
fraction of binary stars leads to the increase of the expected mass of the
system by a factor of 1.1.

The complete mass of all stellar members of the system includes low-mass stars
beyond the photometric completeness limit. We estimated the total hidden mass,
including the unresolved components of binaries, based on the suggestion of a
Kroupa-like IMF (Kroupa, 2001) of the system, with normalisation to the number
of stars with a visible G magnitude $phot\\_g\\_mean\\_mag\\_corr<17.5$. This
provided an estimate of the total mass, including unseen stars with masses
down to $0.09\mathrm{M}_{\odot}$, that is, 1.3 to 1.4 times larger than the
straightforward sum of masses of observed sources. These suggestions provide a
scale of reliability of counted mass and resulted in the asymmetrical
uncertainties indicated in Table 4.

In Table 4, $M^{\Sigma}$ is the apparent counted masses of the
$60\%$-probability members. The indicated confidence limits in the table
represent an extreme possible increase and decrease of this value due to
taking into account hidden members or filtering out a maximum possible random
contamination. The mean mass of a star is shown by $\overline{M_{i}}$. The
mean mass of the probable members of Cr 135 $\overline{M_{i}}$ exceeds the
mean mass of the probable members of UBC 7 and that of the halo (which are
similar; see Table 4). This may be a result of mass segregation.

While the halo as a whole is more massive than both clusters together, its
density is much lower, below $0.001\mathrm{M}_{\odot}/pc^{3}$. The mean
stellar density within 15 pc from the centre of Cr 135 and UBC 7 is 20 to 30
times larger, and in the centre of the cores of the clusters, it is up to
around $3.5\mathrm{M}_{\odot}/pc^{3}$.

Figure 9: Fits of the King profile according to Eq. 5 for the observed radial
density profiles, based on the confident members of the clusters. The curves
present distributions including the residual background stars (black dots on
the yellow lines), the distributions corrected for the background (circles
with statistical error bars), and the fitted King profiles (solid red curves).
The vertical dashed lines show the core and tidal radii ($r_{\rm c}$ and
$r_{\rm t}$). The colour bands show errors in the radii.

###  4.2 Tidal parameters of substructures of the Cr 135 and UBC 7

The spatial structure of an OC can be standardised with the help of a well-
known empirical model proposed by King (1962). The model parameterises the
radial surface density profile of a cluster in terms of its core and tidal
(total) radii $r_{c}$ and $r_{t}$. In the integrated form, the profile can be
represented with the following equation:

$\begin{split}n(r)=\pi\,r_{c}^{2}k\,&\left\\{\ln[1+(r/r_{c})^{2}]-4\frac{\left[1+(r/r_{c})^{2}\right]^{1/2}-1}{\left[1+(r_{t}/r_{c})^{2}\right]^{1/2}}+\right.\\\
&\left.+\frac{(r/r_{c})^{2}}{1+(r_{t}/r_{c})^{2}}\right\\}\,,\end{split}$ (5)

where $n(r)$ is the number of stars projected onto the ideal (tangential)
plane within a circle with radius $r$ and $k$ is a normalising coefficient.

The observed profiles were counted based on the sources of Sample I with
brightness limitation $phot\\_g\\_mean\\_mag\\_corr<18.0$. In order to avoid
mutual contamination, we evaluated the residual background level produced by
the neighbouring population via estimation of the average surface density
produced by the sample stars in the neighbouring circular area outside the
visual limit of the respective cluster. The fit and determination of the
parameters was carried out using the MPFIT procedure from the Markwardt (2009)
IDL library.

Figure 9 shows the empirical profile fits for our working samples. Figure 10
illustrates how reasonably computed tidal parameters outline the general
structure of our pair of star clusters. As one can see, the core radii clearly
mark off the densest parts of the star clusters. In 3D, both clusters are
touching at their tidal limits. We regard this as evidence of the physical
relation of both clusters residing in the common potential well. In fact it is
clear that the bodies that are currently recognised as separate OCs themselves
are rather double cluster cores immersed in the common halo.

In Table 4 we show the derived profile parameters. The fit procedure works
well in our case, and the profile correction for residual outer background
plays a critical role for the correct determination of the tidal radii.

Table 4: Global population characteristics and tidal parameters of the system. Parameter | Cr 135 | UBC 7 | Halo
---|---|---|---
$N$ | 294 | 243 | 884
$R_{HN}$, pc |  6.4 $\pm$ 0.2 |  5.8 $\pm$ 0.2 |  42.8 $\pm$ 2.1
$\overline{M_{i}},\,\mathrm{M}_{\odot}$ |  0.74 $\pm$ 0.04 |  0.60 $\pm$ 0.04 |  0.63 $\pm$ 0.02
$M^{\Sigma},\,\mathrm{M}_{\odot}$ |  218${}^{+87}\rm_{-27}$ |  146${}^{+58}\rm_{-17}$ |  558${}^{+220}\rm_{-70}$
$r_{\rm c}$, pc |  4.21 $\pm$ 0.75 |  3.28 $\pm$ 0.72 | –
$r_{\rm t}$, pc |  15.6 $\pm$ 2.3 |  11.0 $\pm$ 1.9 | –
$m_{\rm t},\,\mathrm{M}_{\odot}$ |  770 $\pm$ 340 |  270 $\pm$ 140 | –

In order to understand whether one can observe a footprint of cluster
gravitational interaction in the observed structure of the pair of star
clusters, we tried to compare tidal parameters for the component profiles
counted in sectors opened at different directions. As we found no systematic
difference is observed, the derived tidal parameters vary within the fit
errors only.

Knowledge regarding cluster structure tidal parameters can be used for
independent mass estimates of an OC. To determine the tidal mass $m_{t}$, we
followed King (1962), Ernst et al. (2010) and Just et al. (2023) and used a
condition for the balance of gravitational forces between the Galaxy and the
OC:

$m_{t}=\frac{r_{J}^{3}}{G}(4-\beta^{2})\Omega^{2},$ (6)

where $r_{J}$ is the Jacobi radius (distance from the cluster centre to
Lagrange points $L_{1},$ or $L_{2}$, where its own gravity is equal to the
Galaxy field), $\varOmega$ is the angular velocity of rotation of the Galaxy
at the cluster’s galactocentric distance $R$, $\beta$ is the ratio of
epicyclic $\kappa$ and rotational $\varOmega$ frequencies, and $G$ is the
gravitational constant.

The relation between Jacoby and tidal radii for realistic clusters in the case
of the three-component Plummer-Kuzmin model of Galactic potential (Miyamoto &
Nagai, 1975b) was studied using $N$-body calculations by Ernst et al. (2010).
They showed that the ratio $r_{t}/r_{J}$ depends on the cluster position on
the sky (mainly on Galactic latitude) and on their age. Its typical value at
low latitudes of $|b|<20\degr$ computed by Ernst et al. (2010) for age of
cluster model of 0.62 Gyr is lower than $r_{t}/r_{J}\approx 1.19$. For the
location of our clusters at $R\approx 8.29$ kpc, the rotation parameters
computed from the Plummer-Kuzmin model are equal to $\beta\approx 1.45$,
$\varOmega\approx 28.20$ km/s/kpc. The respective values of the tidal masses
are shown in Table 4, and the formal uncertainty was estimated via propagation
of the error in $r_{t}$. Since our clusters are much younger than the lower
age limit of the Ernst et al. (2010) model and since we observed that the
ratio tends to decrease at younger ages, we believe that the found values
manifest the upper limits of the respective tidal masses.

Figure 10: Open clusters Cr 135 and UBC 7 in the tangential plane contacting
the celestial sphere at the common centre of the mass of the clusters. The red
dots show probable members of Cr 135, and blue dots show probable members of
UBC 7. The circles show the core and tidal radii. Black dots show the centres
of the clusters. The “+” sign shows the centre of mass, which is assumed to be
at 292 pc from the Sun. Linear coordinates $\xi$ and $\eta$ are parallel to
equatorial ones.

The mass estimate obtained for UBC 7 via correction of counted mass for hidden
mass ($\approx 204M_{\odot}$; see Table 4) agrees rather well with the tidal
mass ($270\pm 140M_{\odot}$). The tidal mass for Cr 135 ($770\pm
340M_{\odot}$) exceeds the estimate obtained from the counted mass even after
correcting for the expected hidden mass ($\approx 305M_{\odot}$) and taking
into account the large uncertainty. This may be related to difficulties with
the treatment of a populated halo. In addition, as noted previously, the found
values probably manifest the upper limits of respective tidal masses.

One may also compare present-day tidal radii with the results of the dynamical
evolution modelling (Fig. 3 for the present time $T=50$ Myr). The model
predicts very similar tidal radii for the two OCs, about 12 pc for Cr 135,
which is slightly less than for UBC 7. Taking into account the errors, the
results match the tidal radii obtained from observations quite satisfactorily,
especially for UBC 7. For Cr 135, the estimated tidal radius ($r_{t}=15.6\pm
2.6$ pc; Table 4) is larger than what the dynamical evolution model predicts.

##  5 Discussion

###  5.1 Cumulative star number and proper motion distributions inside the
modelled clusters

In this section, we present a comparison of the present-day radial cumulative
star count obtained from our basic model with initial conditions presented in
Table 2 with the current observational data (see Section 4.2, Fig. 9 with
Table 4). In Fig. 11, the total cumulative number distribution (CND) of
cluster stars excluding the stellar background but including the observational
errors in number count are presented as a grey zone. For Cr 135 and UBC 7, we
used the fit range within $0<r<20$ pc. We set such a large distance limit in
order to exceed the OCs’ current Jacobi radii (see Table 4). We present our
basic model as a thick black line in these plots, which was obtained directly
from $N$-body simulations. Because of the magnitude limit and incompleteness
of the faint sources in the observational data due to the Gaia satellite
specifications, we selected from the numerical models only the stars that are
within the specific observed stellar mass range: from 0.28 to 4.0
$\mathrm{M}_{\odot}$. In this case, our model (black line) has been adopted as
an observed CND of stars for both clusters.

After the physical parameters of the model clusters were fixed, we carried out
an additional set of 2 $\times$ 9 = 18 random numerical realisations of the
same physical model. In the first nine realisations, we used different samples
of the initial mass function (IMF) of the clusters, while the initial
positions and velocities of the stars were fixed. In the second set of nine
realisations, we randomised the initial positions and velocities (PVs) of the
stars, while the IMF distribution was fixed. This enabled us to estimate the
possible uncertainties directly connected with the randomisation of the
initial individual stellar masses and PVs phase-space distributions.

Figure 11: Cumulative number distribution of stars compared against the
projected radius for Cr 135 (left panels) and for UBC 7 (right panels). The
coloured lines represent different $N$-body randomisations. The upper panel
presents the IMF sampling, and the bottom panels show the PV sampling. The
shaded grey regions show the observational data with uncertainties.

In Fig. 11, all the coloured lines are inside the grey area, especially in the
central 8 pc, which is where 80% of the stars reside. However, our plots show
that half of the coloured lines are on the edge of the grey area, around the
10-11 pc. We interpret this as an indication of a somewhat smaller
concentration of the real initial OCs’ PV model distribution compared to our
initially selected numerical values. Taking into account the OCs’ centres of
mass uncertainties obtained from the 18 samples, we got the next standard
deviation in the three Cartesian Galactic coordinates for Cr 135 ($\pm 2.70$,
$\pm 2.50$, and $\pm 1.86$ pc) and for UBC 7 ($\pm 1.85$, $\pm 1.41$, and $\pm
0.96$ pc). As can be seen, we got slightly larger deviations for Cr 135 as
compared to UBC 7.

In Fig. 12, we present the same CNDs but compare them against the residual
proper motion magnitude inside the clusters. As can be seen from the plots,
the proper motion velocity components of the stars are also quite consistent
with the observed distributions and error bars. Together with the distance
CNDs, these results show that our numerical modelling is able to reproduce the
currently observed 6D phase-space distributions with a high level of
confidence.

Figure 12: Same as in Fig. 11 but compared against the residual magnitude of
the proper motion.

###  5.2 Origin of the extended halo

In Section 3, we pointed out the presence of a halo surrounding both clusters
with a mean density of $3.6\pm 0.2$ stars per square degree. The halo stars
are not in our current $N$-body simulations because we started with two King
spheres in virial equilibrium that are accurate according to mutual attraction
and the Galactic tidal field.

An explanation for the halo can be found in the cluster formation scenario
with low star-formation efficiency (SFE; see Shukirgaliyev et al., 2021b, and
references therein). It suggests that most of the gas is expelled abruptly
from the star-forming region due to stellar feedback from massive OB stars. A
sudden loss of approximate equilibrium in the gas-stars system leads to so-
called violent relaxation followed by star ejection.

To find evidence in favour of this scenario, we evaluated the net residual
velocities of the stars with respect to the centre of the pair of star
clusters in the line-of-sight and tangential directions. The obtained net
line-of-sight velocities of the receding and approaching stars are shown in
Fig. 13. In the figure, the mean residual radial and tangential velocities of
the cluster members (red, blue, and green lines) are compared with low-SFE
cluster values (grey line), which were obtained using data by Shukirgaliyev et
al. (2021b). In this paper, we used the Dehnen $\gamma=0$ initial stellar
profile with SFE $=$ 0.05 after 50 Myr of evolution, which is presented for
the first time. This choice of the model was determined by the bound mass of
the numerical cluster, which should be comparable to the mass of our pair of
star clusters. We calculated the runaway velocities of stars from the centre
of mass of the pair of clusters with the projection on the sky plane. The mean
value of the velocities is shown with the grey line in Fig. 13, and the shadow
represents its standard deviation. The agreement between the net tangential
observed velocity and the mean radial velocity in the simulations is
excellent, supporting the scenario of low-SFE formation of the pair of star
clusters.

Figure 13: Outwards motion of the corona stars with respect to the pair of
star clusters’ centre of mass. Smoothed magnitudes $V_{\rm r}$ of residual
line-of-sight velocities of receding (approaching) stars from Gaia DR3 are
shown by red (blue) lines. The mean residual tangential velocities in the
picture plane $V_{\rm t}$ are shown by the green line. The grey line shows the
mean tangential components of the runaway velocities obtained in $N$-body
simulations for the Dehnen $\gamma=0$ initial stellar profile with SFE $=$
0.05 after 50 Myr of evolution (Shukirgaliyev et al., 2021b), where the shadow
represents the uncertainty.

##  6 Conclusions

We performed star-by-star dynamical modelling of the evolution of the two OCs
Collinder 135 and UBC 7 from their supposed initial state (mean coordinates
and velocities) to their present-day state. The models were compared with the
data obtained based on the performed census of membership of the two clusters
with Gaia DR3 data. In particular, we used CNDs of the cluster members over
the distance from their centres from observations to adjust the initial
parameters of the clusters, such as the number of stars and their
concentration parameter in the model. We selected these parameters in order to
reproduce well the present-day state and distribution of the members of Cr 135
and UBC 7 in 6D space (coordinates, velocities).

Our research supports the potentially shared origin of the Cr 135 and UBC 7
OCs. By integrating the dynamical co-evolution of these clusters within an
external Galactic potential, we obtained results that align closely with the
internal stellar density and velocity distributions of the clusters. In our
parameter-fitting procedure for the clusters, we focused solely on the
resulting stellar CNDs (Fig. 11) at 50 Myr post cluster formation. The model’s
resulting profiles for the residual proper motions and the Gaia DR3 population
of the clusters show a strong agreement without the need for any additional
fitting, Fig. 12.

The star-by-star dynamical evolution model shows the future of the stellar
populations of Cr 135 and UBC 7 (including the existence of the stellar
remnants, such as C/O and O/N white dwarfs). Due to the limited number of
stars (only a few hundred) in both clusters, their 50 Myr-evolved stellar
populations have a lack of neutron stars and black holes.

Based on our comprehensive dynamical modelling with the up-to-date stellar
evolution prescription, we are able to predict the future $\sim$100 Myr
evolution of both clusters with a good degree confidence. In our models, the
UBC 7 cluster lives significantly longer compared to the Cr 135 cluster, Fig.
3. This is not surprising given the fact that, according to our numerical
modelling, UBC 7 is denser (see Table 2) than Cr 135.

During our investigation, we found strong evidence of the presence of the
common stellar halo in the Cr 135 and UBC 7 OC system, and we estimated its
characteristics. The members of the halo demonstrate net residual velocities
of stars with respect to the centre of the pair of star clusters that are
compliant with the mean radial velocity of particles of $N$-body simulations
of a low-SFE cluster stimulated by violent relaxation after gas loss.

Detailed common dynamical modelling of the OCs and the extended halo system
was beyond of the scope of our current paper. Such modelling, however, will
definitely be of interest and valuable for future investigations. The presence
of such an extended halo can also have an impact on the clusters’ cumulative
mass distributions and internal proper motions.

###### Acknowledgements.

The authors thank the anonymous referee for a very constructive report and
suggestions that helped significantly improve the quality of the manuscript.
The work of MI, PB and OS was partially supported under the special program
“Long-term program of support of the Ukrainian research teams at the PAS
Polish Academy of Sciences carried out in collaboration with the US National
Academy of Sciences with the financial support of external partner”,
PAN.BFB.S.BWZ.329.022.2023. The work of PB was also supported by the
Volkswagen Foundation under the special stipend No. 9D154. The work of CO was
funded by the Science Committee of the Ministry of Science and Higher
Education of the Republic of Kazakhstan, Grant No. BR20280974. This research
is funded by the Science Committee of the Ministry of Science and Higher
Education of the Republic of Kazakhstan, Grant No. BR21881880. This work has
made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. The use of TOPCAT, an
interactive graphical viewer and editor for tabular data (Taylor, 2005), is
acknowledged.

## References

  * Alejo et al. (2020) Alejo, A. D., González, J. F., & Veramendi, M. E. 2020, A&A, 633, A146
  * Angelo et al. (2022) Angelo, M. S., Santos, J. F. C., Maia, F. F. S., & Corradi, W. J. B. 2022, MNRAS, 510, 5695
  * Banerjee et al. (2020) Banerjee, S., Belczynski, K., Fryer, C. L., et al. 2020, A&A, 639, A41
  * Beccari et al. (2020) Beccari, G., Boffin, H. M. J., & Jerabkova, T. 2020, MNRAS, 491, 2205
  * Beccari et al. (2018) Beccari, G., Boffin, H. M. J., Jerabkova, T., et al. 2018, MNRAS, 481, L11
  * Berczik et al. (2011) Berczik, P., Nitadori, K., Zhong, S., et al. 2011, in International conference on High Performance Computing, 8–18
  * Berczik et al. (2013) Berczik, P., Spurzem, R., Wang, L., Zhong, S., & Huang, S. 2013, in Third International Conference ”High Performance Computing”, HPC-UA 2013, p. 52-59, 52–59
  * Bisht et al. (2021) Bisht, D., Zhu, Q., Yadav, R. K. S., et al. 2021, MNRAS, 503, 5929
  * Borodina et al. (2021) Borodina, O. I., Carraro, G., Seleznev, A. F., & Danilov, V. M. 2021, ApJ, 908, 60
  * Borodina & Kovaleva (2020) Borodina, O. I. & Kovaleva, D. A. 2020, INASAN Science Reports, 5, 351
  * Borodina et al. (2019) Borodina, O. I., Seleznev, A. F., Carraro, G., & Danilov, V. M. 2019, ApJ, 874, 127
  * Bovy (2015) Bovy, J. 2015, ApJS, 216, 29
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Camargo (2021) Camargo, D. 2021, ApJ, 923, 21
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1
  * Cantat-Gaudin et al. (2019a) Cantat-Gaudin, T., Jordi, C., Wright, N. J., et al. 2019a, A&A, 626, A17
  * Cantat-Gaudin et al. (2019b) Cantat-Gaudin, T., Mapelli, M., Balaguer-Núñez, L., et al. 2019b, A&A, 621, A115
  * Casado (2022a) Casado, J. 2022a, Universe, 8, 368
  * Casado (2022b) Casado, J. 2022b, Universe, 8, 113
  * Casado & Hendy (2023) Casado, J. & Hendy, Y. 2023, MNRAS, 521, 1399
  * Coronado et al. (2022) Coronado, J., Fürnkranz, V., & Rix, H.-W. 2022, ApJ, 928, 70
  * Darma et al. (2021) Darma, R., Arifyanto, M. I., & Kouwenhoven, M. B. N. 2021, MNRAS, 506, 4603
  * de La Fuente Marcos & de La Fuente Marcos (2009) de La Fuente Marcos, R. & de La Fuente Marcos, C. 2009, A&A, 500, L13
  * de la Fuente Marcos & de la Fuente Marcos (2010) de la Fuente Marcos, R. & de la Fuente Marcos, C. 2010, ApJ, 719, 104
  * Dias et al. (2018) Dias, W. S., Monteiro, H., Lépine, J. R. D., et al. 2018, MNRAS, 481, 3887
  * Ernst et al. (2010) Ernst, A., Just, A., Berczik, P., & Petrov, M. I. 2010, A&A, 524, A62
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, A&A, 616, A10
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1
  * Grevesse & Sauval (1998) Grevesse, N. & Sauval, A. J. 1998, Space Sci. Rev., 85, 161
  * Hunt & Reffert (2023) Hunt, E. L. & Reffert, S. 2023, A&A, 673, A114
  * Ishchenko et al. (2023) Ishchenko, M., Sobolenko, M., Berczik, P., et al. 2023, A&A, 678, A69
  * Just et al. (2009) Just, A., Berczik, P., Petrov, M. I., & Ernst, A. 2009, MNRAS, 392, 969
  * Just et al. (2023) Just, A., Piskunov, A. E., Klos, J. H., Kovaleva, D. A., & Polyachenko, E. V. 2023, A&A, 672, A187
  * Kamlah et al. (2022) Kamlah, A. W. H., Leveque, A., Spurzem, R., et al. 2022, MNRAS, 511, 4060
  * Kharchenko et al. (2009) Kharchenko, N. V., Berczik, P., Petrov, M. I., et al. 2009, A&A, 495, 807
  * Kharchenko et al. (2012) Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R.-D. 2012, A&A, 543, A156
  * King (1962) King, I. 1962, AJ, 67, 471
  * King (1966) King, I. R. 1966, AJ, 71, 64
  * Kovaleva et al. (2020) Kovaleva, D. A., Ishchenko, M., Postnikova, E., et al. 2020, A&A, 642, L4 (Paper 1)
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4
  * Lindegren et al. (2018) Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2
  * Maíz Apellániz & Weiler (2018) Maíz Apellániz, J. & Weiler, M. 2018, A&A, 619, A180
  * Markwardt (2009) Markwardt, C. B. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler, 251
  * Miyamoto & Nagai (1975a) Miyamoto, M. & Nagai, R. 1975a, PASJ, 27, 533
  * Miyamoto & Nagai (1975b) Miyamoto, M. & Nagai, R. 1975b, PASJ, 27, 533
  * Olivares et al. (2023) Olivares, J., Lodieu, N., Béjar, V. J. S., et al. 2023, A&A, 675, A28
  * Pang et al. (2020) Pang, X., Li, Y., Tang, S.-Y., Pasquato, M., & Kouwenhoven, M. B. N. 2020, ApJ, 900, L4
  * Pang et al. (2023) Pang, X., Wang, Y., Tang, S.-Y., et al. 2023, AJ, 166, 110
  * Pang et al. (2021) Pang, X., Yu, Z., Tang, S.-Y., et al. 2021, ApJ, 923, 20
  * Piatti & Malhan (2022) Piatti, A. E. & Malhan, K. 2022, MNRAS, 511, L1
  * Priyatikanto et al. (2016) Priyatikanto, R., Kouwenhoven, M. B. N., Arifyanto, M. I., Wulandari, H. R. T., & Siregar, S. 2016, MNRAS, 457, 1339
  * Qin et al. (2023) Qin, S., Zhong, J., Tang, T., & Chen, L. 2023, ApJS, 265, 12
  * Röser et al. (2011) Röser, S., Schilbach, E., Piskunov, A. E., Kharchenko, N. V., & Scholz, R. D. 2011, A&A, 531, A92
  * Shukirgaliyev et al. (2021a) Shukirgaliyev, B., Otebay, A., Sobolenko, M., et al. 2021a, A&A, 654, A53
  * Shukirgaliyev et al. (2021b) Shukirgaliyev, B., Otebay, A., Sobolenko, M., et al. 2021b, A&A, 654, A53
  * Shukirgaliyev et al. (2017) Shukirgaliyev, B., Parmentier, G., Berczik, P., & Just, A. 2017, A&A, 605, A119
  * Shukirgaliyev et al. (2018) Shukirgaliyev, B., Parmentier, G., Just, A., & Berczik, P. 2018, ApJ, 863, 171
  * Spina et al. (2021) Spina, L., Ting, Y. S., De Silva, G. M., et al. 2021, MNRAS, 503, 3279
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Vázquez et al. (2010) Vázquez, R. A., Moitinho, A., Carraro, G., & Dias, W. S. 2010, A&A, 511, A38
  * Ye et al. (2022) Ye, X., Zhao, J., Oswalt, T. D., Yang, Y., & Zhao, G. 2022, AJ, 164, 132
  * Zhong et al. (2019) Zhong, J., Chen, L., Kouwenhoven, M. B. N., et al. 2019, A&A, 624, A34

##  Appendix A Calculation of probabilities of membership to the group of Cr
135 and UBC 7

We adjusted our approach to selecting probable members of a system based on
data obtained in our previous investigation. That is, in Paper I, we found the
number of stars in a close vicinity to the considered pair of OCs with space
parameters suitable to one cluster and kinematic parameters of another
cluster. These stars were not included in the previous strict selection.
However, we interpreted these stars as being probable members of a system and
further computed their joint probability based on parallax and their joint
kinematic probability by taking into account the parameters of both clusters.

###  A.1 Kinematic probability

As in Paper I, we found a halo of stars around the Cr 135 and UBC 7 pair with
a location and kinematics combination that cannot be attributed to only one of
the two clusters but rather to their system. This is why for the selection of
space of coordinates and velocities, we did not distinguish between the two
clusters but rather chose stars with both velocities and parallaxes
sufficiently close to those of either of the two clusters: Cr 135 or UBC 7. We
describe the details of the selection process below.

Provided that our intention was to search for coeval stars probably related to
the pair of clusters in an extended region of sky, the kinematical selection
principles needed to be respectively changed from 2D to 3D velocities to
account for the projection effect. We used the modification of the convergent
point method described by Röser et al. (2011). Of the initial dataset
including 116975 stars, radial velocities were available only for 4575
sources. Thus, we selected probable members related to the system based on the
components of their tangential velocity.

We expected that the members related to Cr 135 and UBC 7 follow certain
distributions with respect to the mean 6D phase-space parameters cited in
Table 1 (based on Gaia DR3). The respective space velocities were obtained as

$\begin{bmatrix}U^{i}\\\ V^{i}\\\
W^{i}\end{bmatrix}=V_{r}^{i}\cdot\begin{bmatrix}\cos l^{i}\cos b^{i}\\\ \sin
l^{i}\cos b^{i}\\\ \sin
b^{i}\end{bmatrix}+\frac{\kappa\mu_{l}^{i}}{\varpi^{i}}\cdot\begin{bmatrix}-\sin
l^{i}\\\ \cos l^{i}\\\
0\end{bmatrix}+\frac{\kappa\mu_{b}^{i}}{\varpi^{i}}\cdot\begin{bmatrix}-\cos
l^{i}\sin b^{i}\\\ -\sin l^{i}\sin b^{i}\\\ \cos b^{i}\\\ \end{bmatrix},$ (7)

where $\kappa=4.74047$ is the transformation factor from 1 mas/yr at 1 kpc to
1 km/s; the indices are $i=1$ for Cr 135 and $i=2$ for UBC 7; the mean values
for coordinates, parallaxes, proper motions, and radial velocities are from
Table 1.

The kinematic probability $P^{k}_{\textrm{kin}}$ for a $k^{\textrm{th}}$ star
to belong to the investigated system is defined as

$P^{k}_{\textrm{kin}}=\max\cases{\exp}\left\\{-\frac{1}{4}\left[\left(\frac{\mu^{k}_{l}-\mu_{l}^{exp_{1}^{k}}}{\varepsilon_{\mu_{l}}}\right)^{2}+\left(\frac{\mu^{k}_{b}-\mu_{b}^{exp_{1}^{k}}}{\varepsilon_{\mu_{b}}}\right)^{2}\right]\right\\},\\\
\exp\left\\{-\frac{1}{4}\left[\left(\frac{\mu^{k}_{l}-\mu_{l}^{exp_{2}^{k}}}{\varepsilon_{\mu_{l}}}\right)^{2}+\left(\frac{\mu^{k}_{b}-\mu_{b}^{exp_{2}^{k}}}{\varepsilon_{\mu_{b}}}\right)^{2}\right]\right\\}\par{}$
(8)

where

$\displaystyle\mu_{l}^{exp_{i}^{k}}=$ $\displaystyle(-\sin l_{k}\cdot
U_{i}+\cos l_{k}\cdot V_{i})/(\kappa/\varpi_{k});$
$\displaystyle\mu_{b}^{exp_{i}^{k}}=$ $\displaystyle(-\cos l_{k}\sin
b_{k}\cdot U_{i}-\sin l_{k}\sin b_{k}\cdot V_{i}+\cos b_{k}\cdot W_{i})/$ (9)
$\displaystyle(\kappa/\varpi_{k})$

are values of proper motion in $(l,b)$ for given coordinates $l_{k},b_{k}$,
and $\varpi_{k}$ expected for the stars sharing space motion
$(U,V,W)_{1},(U,V,W)_{2}$ of Cr 135 and UBC 7, respectively. Scattering
parameters of the proper motion $\varepsilon_{\mu_{l}},\varepsilon_{\mu_{b}}$
were set to 1.0 mas/yr (about 2 km/s in tangential velocity).

###  A.2 Parallax probability

In the tangential plane, we limited our search of members by the coordinates
selected in Sec. 3, but we did not discriminate between the stars by their
position within this scope. The approach for selection by parallaxes is as
follows: We used the probability distribution based on the mean parallaxes of
Cr 135 and UBC 7 from Table 1 and the notion of mean dispersion of parallaxes
of their core members. We intentionally restricted this value to avoid
contamination with peripheral members of neighbouring clusters (primarily NGC
2451B). The resulting formula for the probability of membership based on
parallax value and parallax error is as follows:

$P^{k}_{\varpi}=\cases{1},3.23<\varpi_{k}<3.66\\\
\exp\left\\{-\frac{1}{2}\left(\frac{\varpi^{k}-3.23}{\varepsilon_{\varpi}^{k}}\right)^{2}\right\\},\varpi_{k}\leq
3.23\\\
\exp\left\\{-\frac{1}{2}\left(\frac{\varpi^{k}-3.66}{\varepsilon_{\varpi}^{k}}\right)^{2}\right\\},\varpi_{k}\geq
3.66\,.{}$ (10)

The selected scheme assigns $P^{k}_{\varpi}=1$ to all stars in the line-of-
sight between the centres of the pair of clusters and within about 7.5 pc (in
naive approach $d^{k}[pc]=1000/\varpi_{k}$) to and from each centre. The wings
of the probability distribution were defined by

$\varepsilon_{\varpi}^{k}=\sqrt{\sigma_{\varpi}^{2}+\sigma_{\varpi_{k}}^{2}}\,,$

where $\sigma_{\varpi}^{2}=0.25$ mas is the approximate mean dispersion of
parallaxes of the joint ensemble of reliable members of Cr 135 and UBC 7, and
$\sigma_{\varpi_{k}}^{2}$ is the nominal Gaia EDR3 parallax error for each
source.

###  A.3 Photometric probability

Isochrones for 50 Myr (according to Paper I) for Gaia DR3 passbands from Maíz
Apellániz & Weiler (2018) were obtained from the Padova webserver CMD3.7
333http://stev.oapd.inaf.it/cmd based on the calculations by Bressan et al.
(2012) for solar metallicity $Z=0.0152$. A recent investigation by Spina et
al. (2021) has confirmed close to solar abundance for Cr 135 ($[M/H]=-0.011$).
We took into account the relation between $A_{G}/E({\rm BP-RP})$ and $(BP-
RP)_{0}$ of a star and $A_{0}$ (interstellar extinction at $\lambda=550$ nm)
from Gaia Collaboration et al. (2018a) to adjust isochrones for photometric
probability calculation at a fixed value of $A_{0}=0.1^{m}$, Paper I.

The photometric probability $P^{k}_{\textrm{ph}}$ was computed in the CMD
($G,BP-RP$) as a maximum of probabilities and was calculated separately for Cr
135 (cluster 1) and UBC 7 (cluster 2) as follows. The photometric probability
was set as equal to one for the stars occupying the CMD in an area between the
isochrones for single stars and for unresolved binaries of equal mass. For all
the other stars, the photometric probability depended on a difference between
magnitudes:

$\Delta G^{i}=\min\left(|G^{i}-G^{b}|,|G^{i}-G^{f}|\right),$ (11)

where $G^{b},G^{f}$ are respectively brighter and fainter values of magnitude
limiting the domain between the isochrones for single stars and for unresolved
binaries of equal mass for a given colour $BP-RP$:

$\begin{split}P^{k}_{\textrm{ph}}=\max\left(\exp\left\\{-\frac{1}{2}\left(\frac{\Delta
G_{1}^{k}}{\varepsilon^{k}_{G_{1}}}\right)^{2}\right\\},\exp\left\\{-\frac{1}{2}\left(\frac{\Delta
G_{2}^{k}}{\varepsilon^{k}_{G_{2}}}\right)^{2}\right\\}\right).\end{split}$
(12)

The value of $\varepsilon^{i}_{\rm ph}$ was defined from the mean dispersion
of parallaxes of the joint ensemble of reliable members of Cr 135 and UBC 7,
which is equal to $0.25$ mas, and the individual photometric error
$\sigma^{i}_{G}$ estimated by the flux and flux error in G.