11institutetext: Department of Astronomy, School of Science, The University of
Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan
11email<EMAIL_ADDRESS>22institutetext: RIKEN Center for
Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo
648-0047, Japan 33institutetext: European Space Agency (ESA), European Space
Research and Technology Centre (ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, The
Netherlands,33email<EMAIL_ADDRESS>
# The impact of massive stars and black holes on the fate of open star
clusters and their tidal streams

Long Wang 1122 Tereza Jerabkova 33

(Received June 15, 2021; accepted –, –)

###### Abstract

Context. To investigate how the content of massive OB stars affects the long-
term evolution of young open clusters and their tidal streams, and how such an
effect influences the constraint of initial conditions by looking at the
present-day observations.

Aims. OB stars are typically in binaries, have a strong wind mass loss during
the first few Myr, and many become black holes. These affect the dynamical
evolution of an open star cluster and impact its dissolution in a given
Galactic potential. We investigate the correlation between the mass of OB
stars and the observational properties of open clusters. Hyades-like star
clusters are well represented in the Solar neighborhood and thus allow
comparisons with observational data.

Methods. We perform a large number of star-by-star numerical $N$-body
simulations of Hyades-like star clusters by using the high-performance
$N$-body code petar combined with galpy.

Results.  We find that OB stars and black holes have a major effect on star
cluster evolution. Star clusters with the same initial conditions, but a
different initial content of OB stars, follow very different evolutionary
paths. Thus, the initial total mass and radius of an observed star cluster
cannot be unambiguously determined unless the initial content of OB stars is
known. We show that the stellar counts in the corresponding tidal tails, that
can be identified in the Gaia data, help to resolve this issues. We thus
emphasise the importance of exploring not only star-clusters, but also their
corresponding tidal tails. These findings are relevant for studies of the
formation of massive stars.

###### Key Words.:

star clusters – star formation – numerical methods

##  1 Introduction

Young open star clusters are important stellar systems for understanding star
formation and stellar evolution. Due to the stellar and dynamical evolution,
star clusters are evaporating their stars and are in the process of becoming
part of the Galactic field stellar population. Thus, for many star clusters a
large fraction of their stars are not anymore located in the clusters
themselves, but in the large-scaled tidal tails shaped by the Galactic tides
and internal cluster processes (Baumgardt & Makino, 2003; Heggie & Hut, 2003;
Chumak & Rastorguev, 2006; Küpper et al., 2008, 2010; Dinnbier & Kroupa,
2020a, b; Jerabkova et al., 2021).

The Gaia survey, with its currently most up to date third data release EDR3
(Gaia Collaboration et al., 2021), provides rich astrometric and photometric
data allowing to study kinematic properties of stars in the solar neighborhood
in unprecedented detail. A new era comes as more and more new open clusters
(e.g. Cantat-Gaudin et al., 2018, 2019; Castro-Ginard et al., 2018, 2020; Sim
et al., 2019; Liu & Pang, 2019; Ferreira et al., 2020; He et al., 2021), large
scale co-eval relic filaments (Jerabkova et al., 2019; Beccari et al., 2020)
and extended stellar streams (Kounkel & Covey, 2019) have been discovered
recently using Gaia. With the Gaia data, for the first time, it has become
possible to extract tidal tails of open star clusters, a challenging task due
to the difficulty to distinguish tail stars from the Galactic field population
(e.g. Röser & Schilbach, 2019a; Röser, Schilbach, & Goldman, 2019b; Fürnkranz
et al., 2019; Meingast, Alves & Fürnkranz, 2019; Tang et al., 2019; Zhang et
al., 2020; Bhattacharya et al., 2021; Jerabkova et al., 2021). The kinematic
properties of Galactic young star clusters (e.g. Kuhn et al., 2019; Monteiro &
Dias, 2019; Monteiro et al., 2020; Zhong et al., 2020; Angelo et al., 2021;
Dias et al., 2021; Godoy-Rivera et al., 2021) and their extended structures
(e.g. Pang et al., 2020, 2021; Meingast, Alves, & Rottensteiner, 2021) have
also been significantly improved.

With data on the tidal tails of open clusters becoming available, the initial
conditions of these can be better constrained. These in turn constrain star
formation theories: Do all stars form in initially gravitationally bound
clusters? What are the radii and masses of these, and is the stellar and
binary population in these always the same? Inferring the initial conditions
of star clusters constitutes an important astrophysical problem, because the
observed coeval stellar population can then be matched to the physical
conditions of its birth. However, star cluster formation starting from giant
molecular clouds involves a multi-scale complexity. The formation of a star
cluster is not an isolated event, but is embedded in a large-scale
distribution of molecular could. For example, the Orion molecular clouds has a
length scale of hundreds of pc where several star-forming regions exist. (e.g.
Bally et al., 1987; Genzel & Stutzki, 1989; Carpenter, 2000; Dame et al.,
2001; Kounkel et al., 2018; Jerabkova et al., 2019; Beccari et al., 2020).

Meanwhile, in the dense region where star clusters form, the wind, the
radiation and the supernovae feedback from massive OB stars significantly
affects the star formation. They can remove the gas from the cluster-forming
region in a short timescale ($\approx 1\,$Myr) (e.g. Kroupa, Aarseth, &
Hurley, 2001; Goodwin & Bastian, 2006; Baumgardt & Kroupa, 2007; Dinnbier &
Walch, 2020; González-Samaniego & Vazquez-Semadeni, 2020; Fujii et al.,
2021b). If the star formation efficiency is low, the gas expulsion can
significantly change the gravitational potential and even cause the immediate
destruction of the system. The feedback can also affect the surrounding large-
scale molecular cloud, such as the formation of the Orion–Eridanus Superbubble
(e.g., Blaauw, 1964; Brown et al., 1994; Madsen et al., 2006; O’dell et al.,
1967; O’Dell et al., 2011; Ochsendorf et al., 2015; Kounkel, 2020; Großschedl
et al., 2021).

It is challenging to fully understand this complexity in star formation.
However, during the post-gas-expulsion phase, if a star cluster was still
bound, the stellar dynamics would re-virialise the cluster and would result in
a roughly spherical stellar system as observed today. In this work, we
strengthen that the massive OB stars can also significantly influence the
long-term dynamical evolution of the post-gas-expulsion star clusters. The
wind mass loss of OB stars after gas expulsion continues to reduce the
gravitational potential of the star cluster. Meanwhile, the death of OB stars
can leave compact objects like black holes (BHs). By being much more massive
than stars, they are the most likely candidates to end up in dynamically
formed binaries. The interaction between these binaries and surrounding stars
is a heating mechanism that influences the long-term evolution of the system
(e.g. Binney & Tremaine, 1987; Spitzer, 1987; Mackey et al., 2008; Breen &
Heggie, 2013).

The number of OB stars formed is related to the fundamental nature of the
initial mass function (IMF) in a star-forming region (Kroupa et al., 2013). If
star formation were to be stochastic, the number of OB stars in a gas cloud
would be described by randomly sampling from the IMF. Thus, even if a group of
low-mass open clusters would have the same initial properties (i.e. the same
masses, radii), they would show noticeable differences after hundreds of Myr
of dynamical evolution due to the large variation of the number of OB stars.
In contrast, if star formation were to be highly self-regulated, a strong
correlation of the number of OB stars and the total mass of the cluster-
forming gas cloud would result, and the stochastic effect would be
significantly suppressed (e.g., Weidner & Kroupa, 2006; Weidner, Kroupa, &
Pflamm-Altenburg, 2013; Kroupa et al., 2013; Yan, Jerabkova, & Kroupa, 2017).

By numerically modelling the evolution of Hyades-like open clusters, we can
quantify the cluster bulk properties (mass, radius) as a function of time. For
low-mass open star clusters which only contain thousands of stars or fewer,
the stochastic effect from the randomization of initial masses, positions and
velocities of stars cannot be ignored. The evolution of cluster, such as the
density, morphology and the distribution of escaping stars, diverge
increasingly with time due to the inherent chaotic nature of the system
(Heggie, 1988; Goodman, Heggie, & Hut, 1993; Heggie, 1996). By allowing, for a
model of fixed initial mass, the IMF, positions and velocities of stars to be
randomly sampled, we study the degeneracy of initial configurations of post-
gas-expulsion, expanded re-virialised Hyades-like open clusters (the most
comprehensively observed on) in order to constrain the maximum possible
variation of such initial conditions that yield a comparable present-day
configuration. The maximum variation is obtained by treating the IMF as a
probabilistic distribution function rather than an optimally sampled
distribution function (Kroupa et al., 2013). We point out for the first time
that by combining the present-day properties of the open cluster (its mass,
radius, stellar population) with the population of stars in its tidal tails
and their extend, the degeneracy can be broken: the post-gas-expulsion initial
conditions for an observed open cluster can be inferred uniquely if its
astrophysical age, its bulk properties and the stellar population in its tidal
tails are known. In particular, the initial content of OB stars can be
determined. Such constraints can then, in the future, be linked to the pre-
gas-expulsion birth conditions by applying hydrodynamical modelling of the
forming embedded cluster (e.g. Hirai et al., 2021; Fujii et al., 2021a, b).

First, in Section 2, we describe the numerical $N$-body code, petar equipped
with galpy, and the initial condition for modelling the formation and
evolution of a Hyades-like star cluster and its tidal stream. Then, in Section
3, we analyze the observational features and dynamical evolution of our models
in detail. Especially, we focus on how the initial mass of OB stars affects
the present-day properties of the cluster and its tidal stream. Finally, we
summarize and discuss our findings in Section 4.

##  2 Method

###  2.1 The $N$-body code petar

In this work, we use the $N$-body code petar (Wang et al., 2020) to perform
the numerical simulations of the star clusters. petar is a high-performance
$N$-body code that combines three algorithms:

  * •

the Barnes-Hut tree (Barnes & Hut, 1986) for the long-range interactions;

  * •

the fourth-order Hermite integrator with block time steps for middle-range
interactions (e.g., Aarseth, 2003);

  * •

and the slow-down algorithmic regularization method (SDAR; Wang, Nitadori &
Makino, 2020) for short-range interactions inside multiple systems, such as
hyperbolic encounters, binaries and hierarchical few-body systems.

The three algorithms are combined via the particle-tree particle-particle
method (Oshino, Funato, & Makino, 2011). The framework for developing parallel
particle simulation codes (fdps) is used to achieve a high performance with
the multi-process parallel computing (Iwasawa et al., 2016, 2020; Namekata et
al., 2018).

Although the Barnes-Hut tree introduces an approximation of the force
calculation to reduce the computational cost for the massive stellar systems
with $>10^{5}$ stars, petar can also well simulate low-mass systems when the
accuracy parameters are set properly. The high performance and the low
computational cost allow us to carry out a large ensemble of $N$-body
simulations of low-mass clusters.

###  2.2 The stellar evolution package sse/bse

The recently updated version of single and binary stellar evolution packages,
sse and bse, are used to simulate the wind mass loss, the type changes of
stars, the mass transfer and the mergers of binaries (Hurley, Pols, & Tout,
2000; Hurley, Tout, & Pols, 2002; Banerjee et al., 2020). These population
synthesis codes use the semi-empirical stellar wind prescriptions from
Belczynski et al. (2010). For the formation of compact objects, we adopt the
“rapid” supernova model for the remnant formation and material fallback from
Fryer et al. (2012), along with the pulsation pair-instability supernova
(PPSN; Belczynski et al., 2016). We apply the solar metallicity ($Z=0.02$) in
this study.

###  2.3 The galactic potential package galpy

The galpy code (Bovy, 2015) has been implemented as a plugin to petar via the
c programme language interface. The petar-galpy combination can simulate the
formation and evolution of massive star clusters with their tidal streams in a
variety types of galactic potentials (see Appendix A). In addition, we also
implement the petar-Gaia Python tool for the data analysis and the comparison
with observational data (see Appendix B).

The Python interface of galpy provides more functions than the c interface.
However, the latter has a much better computing performance. In this work,
using petar we can efficiently generate a large amount of simulations for
dynamically evolving open star clusters in the Milky-way Potential. One
simulation listed in the Section 2.4 needs a computing wall-clock time of 7
min to 4 hours with one CPU core. The 4500 models in this study cost about 4
days computing to finish by using a desktop computer with a AMD Ryzen
Threadripper 3990x CPU (64 cores).

In our models, we setup the Galactic potential by using the MWPotential2014
from galpy. This is a simplified static Galactic potential model including a
bulge, a disk and a dark-matter halo. Since exactly reproducing the observed
structure along the tidal stream of Hyades is not the purpose of this work,
this approximated potential is sufficient for our study. The parameters of
MWPotential2014 are listed as following:

$\begin{array}[]{lp{0.8\linewidth}}Bulge&a power-law density profile with an
exponentially cut-off\\\ &\-- power-law exponent: -1.8\\\ &\-- cut-off radius:
1.9 kpc\\\ &\-- mass: $5\times 10^{9}M_{\odot}$\\\
Disk&\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Miyamoto1975}{\@@citephrase{(}}{\@@citephrase{)}}} disk\\\
&\-- scale length: 3.0\leavevmode\nobreak\ kpc\\\ &\-- scale height:
280\leavevmode\nobreak\ pc\\\ &\-- mass: $6.8\times 10^{10}M_{\odot}$\\\
Halo&NFW profile
\cite[citep]{(\@@bibref{AuthorsPhrase1Year}{NFW1995}{\@@citephrase{, }}{})}\\\
&\-- scale radius: 16\leavevmode\nobreak\ kpc\\\ &\-- local dark-matter
density: 0.008 $M_{\odot}pc^{-3}$\\\ \end{array}$

The MWPotential2014 assumes a solar distance to the Galactic center to be 8
kpc and the solar velocity to be 220 km s-1.

###  2.4 Initial conditions

Figure 1: The orbit of Hyades-like cluster in the Galactocentric reference
frame calculated by using galpy. The upper panel shows the orbit in the
$x$-$y$ plane and the lower panel shows the orbit in the $r$-$z$ plane, where
$r$ is the radial distance projected in the x-y plane. The symbol of star
represents the present-day position.

####  2.4.1 The real Hyades cluster

The Hyades is the closest star cluster to the Sun, with a distance of 45 pc.
Röser et al. (2011); Röser, Schilbach, & Goldman (2019b) constrained the
Hyades’ half-mass and the tidal radii to be 4.1 pc and 9 pc (up to 10.4 pc),
respectively, with the total stellar mass of the Hyades being 435 $M_{\odot}$
(up to 509 $M_{\odot}$) and the bound mass being 275 $M_{\odot}$ (up to 322
$M_{\odot}$). The upper estimated are correction based on assumptions about
unresolved binaries, for more details see Röser et al. (2011). See also Reino
et al. (2018) for constraints on the Hyades cluster membership and velocity
dispersion and relevant discussion in Jerabkova et al. (2021).

To find the initial position and velocity of the Hyades cluster in the Galaxy,
we reverse the velocity of the center of Hyades at the present day and
integrate back 648 Myr by using the time-symmetric integrator in galpy to find
the initial position (see Figure 1). Then, we reverse the velocity 648 Myr ago
to obtain the initial velocity. Since the MWPotential2014 is static, if we
start from this initial coordinate and integrate forwards to the present day,
we can recover the correct position and velocity with a small numerical error.
The evaluated initial and the observed present-day coordinates of the Hyades
in the International Celestial Reference System (ICRS) frame are listed in
Table 1 (see Jerabkova et al., 2021; Gaia Collaboration et al., 2018, for more
details).

The stellar wind mass loss and supernovae can break the conservation of the
momentum. Thus, in the $N$-body simulation using petar, the center of the star
cluster may not exactly follow the orbit obtained from the result using galpy.
But the offset at the present day is small so that it does not affect our
analysis. We can also correct the offset by rotating the coordinates of the
stellar system in the Galactrocentric frame111based on
coordinates.basecoordinateframe from the astropy Python module (see Astropy
Collaboration et al., 2013, 2018).

Table 1: The initial and present-day coordinates of Hyades

Parameter | 648 Myr ago | present-day
---|---|---
RA[$\deg$] | 308.488 | 67.985
Dec[$\deg$] | 41.217 | 17.012
Distance[pc] | 4.575 | 47.501
pm(RA,cos(Dec))[mas yr-1] | -0.904 | 101.005
pm(Dec) [mas yr-1] | -0.740 | -28.490
radial velocity [km s-1] | 431.31 | 39.96
222The evaluated coordinate of the Hyades 648 Myr ago by using galpy and the
observed present-day one in the ICRS frame (see Jerabkova et al., 2021; Gaia
Collaboration et al., 2018, for more details)

####  2.4.2 The models

Figure 2: Sketch of the formation and evolution of a star cluster (from left
to right). The left: the gas embedded phase where star formation occurred
along the filaments in a few Myr, which were influenced by the gas flow in the
large scale, the OB feedback and the stellar dynamics. The right: after the
gas expulsion, the bound star cluster re-virialised to a roughly spherical
system, which is the initial condition in our models. With the three black
solid arrows we show that there are a number of evolutionary tracks a star
cluster can take. In our case we focus on how the content of massive stars can
affect the evolution of the cluster with otherwise identical initial
conditions. We show that while this causes degeneracy making it impossible to
uniquely infer the initial conditions of an observed cluster, the properties
of its tidal tails help to break the degeneracy and thus constrains the
initial conditions, as marked by the gray solid arrow. The ultimate question
(indicated by the gray dotted arrow) is then what constraints can be placed on
the gas-embedded phase of the star cluster.

Fig. 2 shows a sketch describing the phase transition for the evolution of a
star cluster. The gray-shaded panel shows the gas-embedded phase involving the
following general steps: i) initial convergent gas flow towards a density
maximum, ii) on-set of star-formation along the filaments and predominantly at
their intersection comprising the bulk of the later stellar body of the
embedded star cluster, iii) star cluster formation being affected by dynamical
processes (mergers, ejections) and stellar feedback as exemplified by the
Orion Nebula Cluster (Kroupa et al., 2018), iv) gas expulsion which gradually
begins with the on-set of star formation and results in a largely gas-free
star cluster. Phases (i)-(iii) are encompassed in the birth conditions as
deduced by Marks & Kroupa (2012). All of these processes are contributing to
the star cluster’s evolutionary path in a complex way as depicted by the
bottom solid arrow pointing to the gas-free star cluster. This point marks the
start of our $N$-body simulations, where we assume the embedded pre-gas-
expulsion cluster (gray-shaded panel) has expanded and re-virialised after the
residual gas has been removed (Banerjee & Kroupa, 2017). Thus, we do not aim
to simulate the formation of the cluster but the evolution starts from the
expanded gas-free state. The long-term evolution of the star cluster and its
tidal stream is controlled by the net effect from the stellar evolution,
stellar dynamics and the Galactic potential. In this work, we focus on how the
content of massive stars can affect the evolution, which causes the degeneracy
of the evolution track and brings the challenge to identify the initial
condition from the present-day observational property.

Observationally, gas-free star clusters are well resembled by Plummer-like
profiles (Plummer, 1911; Röser et al., 2011; Röser & Schilbach, 2019a) and
thus our initial set-up is empirically motivated. Using the initial coordinate
described in Section 2.4.1 in the Milky-way potential, we generate a grid of
star-cluster models by varying the initial total masses $M_{0}$ and the
initial half-mass radii $R_{\mathrm{h},0}$ (see Table 3). We combine each of
$M_{0}$ and $R_{\mathrm{h},0}$, and thus, the total number of model sets is
15. The reference name of each set is described in Table 3. To investigate the
stochastic effect, it is necessary to carry out a large number of simulations
with different random seeds. Thus, for each set, we generate 300 models by
randomly sampling the IMF from Kroupa (2001) with the mass range of
$0.08-150\leavevmode\nobreak\ M_{\odot}$ and the positions and velocities of
stars from the Plummer density profile (Aarseth, Henon, & Wielen, 1974). The
stochastic effect also naturally results in a different initial mass of OB
stars ($M_{\mathrm{OB}}$), where we assume OB stars have the zero-age mass
sequence mass $>8M_{\odot}$.

Table 2: The initial parameters of $N$-body models

$M_{0}[M_{\odot}]$ | 800, 1000, 1200, 1400, 1600
---|---
$R_{\mathrm{h},0}[pc]$ | 0.5, 1, 2
333The initial total masses $M_{0}$ and the initial half-mass radii
$R_{\mathrm{h},0}$ of the star cluster model sets. The reference name of each
set combines the values of $M_{0}$ and $R_{\mathrm{h},0}$ without the decimal
point, such as M800R05, M1000R1, M1600R2.

##  3 Results

###  3.1 The variation of $M_{\mathrm{OB}}$ due to randomly sampling of IMF

Figure 3: The mass of the heaviest star ($m_{\mathrm{max}}$) and the total
mass of OB stars ($M_{\mathrm{OB}}$) by randomly sampling the Kroupa (2001)
IMF in each of the 300 Hyades-like models. The initial total mass of each
model is $1600\leavevmode\nobreak\ M_{\odot}$.

The randomly sampling of IMF results in a large dispersion of the mass of the
heaviest star ($m_{\mathrm{max}}$) and $M_{\mathrm{OB}}$ model by model. Fig.
3 shows the $m_{\mathrm{max}}-M_{\mathrm{OB}}$ relation from all 300 models in
the M1600R2 set. $M_{\mathrm{OB}}$ varies from $150\leavevmode\nobreak\
M_{\odot}$ to $600M_{\odot}$ and $m_{\mathrm{max}}$ varies from 20 $M_{\odot}$
to $140\leavevmode\nobreak\ M_{\odot}$. The maximum $M_{\mathrm{OB}}$
comprises $37.5$ per cent of the initial mass of the cluster. Thus, it is
expected that the dynamical impact from the OB stars can differ significant
model by model.

###  3.2 The stellar winds from OB stars

Figure 4: The comparison of the number density map for the positions and
proper motions of stars at the age of 648 Myr in the two models with the
minimum and the maximum $M_{\mathrm{OB}}$ from the M800R2 set (OB-min and OB-
max models).

To investigate how $M_{\mathrm{OB}}$ affects the present-day morphology, we
pick up the two models with the minimum and the maximum $M_{\mathrm{OB}}$ in
the M1600R2 set (hereafter named as “OB-min” and ”OB-max” models,
respectively). Fig. 4 compares the present-day positions and proper motions of
stars in the ICRS frame of the two models. The OB-min model in the left panels
has $M_{\mathrm{OB}}\approx 154M_{\odot}$ and $m_{\mathrm{max}}\approx
20M_{\odot}$, while the OB-max model in the right panel has a much larger
$M_{\mathrm{OB}}$ and its $m_{\mathrm{max}}$ is comparable to
$M_{\mathrm{OB}}$ of the OB-min model. Their present-day morphology looks
completely different: the host star cluster still exists in the OB-min model
while it has already been disrupted in the OB-max model. This suggests the
strong impact from $M_{\mathrm{OB}}$ on the long-term dynamical evolution.

The proper motions of the two branches of tidal tails also show different
distributions (the upper regions of the lower panels in Fig. 4). This may be
caused by the different offsets of the central position and velocity of the
two models. In the OB-min model, we can correct the central coordinate to be
the exact observed one. But it is difficult for the OB-max model where the
center cannot be determined.

Figure 5: The comparison for the evolution of the global properties for the
OB-min and OB-max models shown in Fig. 4. From the top to the bottom: tidal
radii ($R_{\mathrm{tid}}$), the total masses inside $R_{\mathrm{tid}}$
($M_{\mathrm{tid}}$), the half-mass radii ($R_{\mathrm{h}}$) and the core
radii ($R_{\mathrm{c}}$). The vertical line refers to 30 Myr.

To investigate when the difference started to develop, we compare the
evolution of several global parameters for the two models including the tidal
radii ($R_{\mathrm{tid}}$), the total masses inside $R_{\mathrm{tid}}$
($M_{\mathrm{tid}}$), the half-mass radii ($R_{\mathrm{h}}$) and the core
radii ($R_{\mathrm{c}}$). To calculate $R_{\mathrm{c}}$, we apply the method
defined by Casertano & Hut (1985) as

$R_{\mathrm{c}}=\sqrt{\frac{\sum_{i}\rho_{i}^{2}r_{i}^{\prime
2}}{\sum_{i}\rho_{i}^{2}}},$ (1)

where $\rho_{i}$ is the local density of object $i$ estimated by counting 6
nearest neighbors, and $r_{i}^{\prime}$ is the distance to the central
position of the system. The center is estimated by the density weighted
average:

$\mathbf{r}_{\mathrm{cm}}=\frac{\sum_{i}\rho_{i}\mathbf{r}_{i}}{\sum_{I}\rho_{i}}$
(2)

where $\mathbf{r}_{i}$ is the position vector of each star or compact object
in the coordinate system used in the simulation. During the simulation, the
origin point of the coordinate system follows the motion of the potential
weighted center of the star cluster. Notice that the definition of
$R_{\mathrm{c}}$ from observation is different. Especially, compact objects
cannot be detected by observations. For massive star clusters where many BHs
exist, $R_{\mathrm{c}}$ in $N$-body models can be very different from that of
observations. In this study, only a few BHs can form and their impact on the
calculation of $R_{\mathrm{c}}$ is neglectable. We have confirmed this in our
models.

In the Galactic potential, there is no precise definition of
$R_{\mathrm{tid}}$. But we can obtain the rough estimation via

$R_{\mathrm{tid}}=R_{\mathrm{gal}}\left(\frac{M_{\mathrm{tid}}}{3M_{\mathrm{g}al}}\right)^{\frac{1}{3}}$
(3)

(Binney & Tremaine, 1987), where $R_{\mathrm{gal}}$ is the distance to the
Galactic center and $M_{\mathrm{gal}}$ is the effective mass of the Galaxy
(the mass enclosed by $R_{\mathrm{gal}}$). In our analysis, we first calculate
the (positive) Galactic potential at the cluster center ($P_{\mathrm{gal}}$)
and then approximate $M_{\mathrm{gal}}$ via

$M_{\mathrm{gal}}=\frac{P_{\mathrm{gal}}r_{\mathrm{gal}}}{G}$ (4)

where $G$ is gravitational constant. At the beginning, $M_{\mathrm{tid}}$ is
unknown, thus the total mass of all objects is used as the starting
$M_{\mathrm{tid}}$. Then a few iterations are necessary to obtain the
consistent $R_{\mathrm{tid}}$ and $M_{\mathrm{tid}}$. We stop the iteration
when the difference between the new and the old $R_{\mathrm{tid}}$ is less
then 1 percent. When the cluster reaches the disruption phase, the iteration
may not converge to a positive value of $R_{\mathrm{tid}}$, i.e.,
$R_{\mathrm{tid}}$ becomes zero. Once $R_{\mathrm{tid}}$ is determined, the
half-mass radii ($R_{\mathrm{h}}$) is calculated by counting all stars inside
$R_{\mathrm{tid}}$.

The comparison of these parameters for the two models is shown in Fig. 5. The
difference already appears during the first 30 Myr. $R_{\mathrm{tid}}$ and
$M_{\mathrm{tid}}$ of the OB-max model fast decrease due to the strong wind
mass loss of OB stars. More than half of the initial mass has lost during this
period. As the gravitational potential changes, $R_{\mathrm{h}}$ and
$R_{\mathrm{c}}$ of the OB-max model also increases significantly at the
beginning. After approximately 400 Myr, $R_{\mathrm{h}}$ decreases while
$R_{\mathrm{c}}$ significantly increases. This indicates that the star cluster
loses (virial) equilibrium, and thus, it rapidly expands and approaches to the
disruption (Fukushige & Heggie, 1995). This is also reflected by the
oscillation of $R_{\mathrm{tid}}$, where zero appears sometimes.

In contrast, the OB-min model has a smooth evolution and the star cluster
still survive at 648 Myr with more than half of the initial mass being inside
$R_{\mathrm{tid}}$. The evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ is
flat until $648$ Myr.

###  3.3 The general trend depending on $M_{\mathrm{OB}}$

Figure 6: The comparison for the evolution of the global properties for all
models with different $M_{0}$ (represented by colors) and different
$R_{\mathrm{h,0}}$ (represented by markers and brightness of colors). From the
top to the bottom are: $R_{\mathrm{tid}}$, $M_{\mathrm{tid}}$,
$R_{\mathrm{h}}$ and $R_{\mathrm{c}}$. The 300 models in each set are
collected into the bins of $M_{\mathrm{OB}}$. The error bars represent the
standard deviation. The observed values of Hyades are shown as horizontal
lines or regions (with uncertainties).

The OB-min and OB-max are two specific cases. To investigate the general
trend, we separate the 300 models in each model set into five bins of
$M_{\mathrm{OB}}$. Each bin contains the same number of models (60). Then, for
each model in the bin, we obtain $R_{\mathrm{tid}}$ and $M_{\mathrm{tid}}$ at
648 Myr and the average values of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ from
the snapshots at 632, 640 and 648 Myr. The average can reduces the fluctuation
that appears in the evolution of $R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ as
shown in Fig. 5. Finally, we calculate the statistic average and the standard
deviation of each parameter (including $M_{\mathrm{OB}}$) inside each bin. The
result is shown in Fig. 6. There is a pronounced trend that a larger
$M_{\mathrm{OB}}$ results in smaller $R_{\mathrm{tid}}$ and $M_{\mathrm{tid}}$
(faster mass loss).

The reference data from Hyades are shown as grey horizontal lines and regions
(with uncertainties). Most of models have slightly lower $R_{\mathrm{tid}}$,
$R_{\mathrm{h}}$ and $R_{\mathrm{c}}$ compared to those of the observational
data. There is a degeneracy between $M_{\mathrm{OB}}$ and $M_{0}$ that many
combinations of the two parameters can result in the similar present-day
property of star cluster. For example, all five $M_{0}$ set (different colors)
can result in the same $R_{\mathrm{tid}}$ ($4-8$ pc), $M_{\mathrm{tid}}$
($100-300M_{\odot}$), $R_{\mathrm{h}}$ ($2-3$ pc) and $R_{\mathrm{c}}$ ($1-2$
pc). This suggests that it is hard to determine the initial condition by only
checking these four parameters.

The impact of $R_{\mathrm{h,0}}$ is not significant. Generally, smaller
$R_{\mathrm{h,0}}$ result in larger $R_{\mathrm{tid}}$ ($M_{\mathrm{tid}}$).
But the scatter is large and the effect is less than that from $M_{0}$ and
$M_{\mathrm{OB}}$.

###  3.4 The long-term dynamical impact from black hole binaries

Figure 7: The evolution of the global properties, from the top to the bottom:
$R_{\mathrm{tid}}$, $M_{\mathrm{tid}}$, $R_{\mathrm{h}}$, $R_{\mathrm{c}}$ and
$a_{\mathrm{BBH}}$, for the models with a similar property at 32 Myr to those
in the OB-max model. The columns separates the models with different $M_{0}$.
For each $M_{0}$, a brighter color indicate a higher total mass of BHs
($M_{\mathrm{BH}}$). The half-mass relaxation time ($T_{\mathrm{rh}}$) at 32
Myr are shown in the second row of panels. The values of $M_{\mathrm{BH}}$ and
the corresponding numbers of BHs (in brackets) are shown in the last row of
panels.

The major impact from stellar wind mass loss of OB stars appears at the first
32 Myr as shown in Fig. 5. Thus, if a model has a different initial condition
but shows similar properties (e.g. mass, density and $R_{\mathrm{h}}$) at 32
Myr as those in the OB-max model, they may have identical evolution after
that. To investigate this, we measure the difference between models as

$\Delta=\left(\frac{M_{\mathrm{tid}}-M_{\mathrm{tid,ref}}}{M_{\mathrm{tid,ref}}}\right)^{2}+\left(\frac{R_{\mathrm{h}}-R_{\mathrm{h,ref}}}{R_{\mathrm{h,ref}}}\right)^{2}+\left(\frac{N_{\mathrm{rh}}-N_{\mathrm{rh,ref}}}{N_{\mathrm{rh,ref}}}\right)^{2}$
(5)

where $N_{\mathrm{rh}}$ is the number of stars within $R_{\mathrm{h}}$ and the
suffix “ref” refers to the parameters of the OB-max model. The values of
$M_{\mathrm{tid}}$ at 32 Myr are used for comparison. $R_{\mathrm{h}}$ and
$M_{\mathrm{rh}}$ have large fluctuations as shown in Fig. 5. Thus, their
averaged values at 24, 32 and 40 Myr are compared.

Baumgardt & Makino (2003) showed that the dissolution time of a star cluster
depends the strength of tidal field ($R_{\mathrm{tid}}/R_{\mathrm{h}}$) and
the relaxation time. The formula of the half-mass relaxation time for one
component system can be described as (Spitzer, 1987)

$T_{\mathrm{rh}}\approx
0.138\frac{N^{1/2}R_{\mathrm{h}}^{3/2}}{m^{1/2}G^{1/2}\ln\Lambda},$ (6)

where $\Lambda=0.02N$, $m$ is the mass of star (we use the average mass
instead) and $G$ is the gravitational constant. The factor of $0.02$ follows
the measurement from Giersz & Heggie (1996). The three parameters in Equation
5 can reflect the tidal effect and $T_{\mathrm{rh}}$.

We select models with first five smallest $\Delta$ for each $M_{0}$ set. Fig.
7 shows the comparison for the dynamical evolution of these models. The models
with a larger common $M_{0}$ tend to have a more divergent long-term
evolution. For the M800 set, all models except one have an identical
evolution, while the M1600 models show large differences of $R_{\mathrm{tid}}$
and $M_{\mathrm{tid}}$ after 200 Myr. A part of M1600 models are still bound
at 648 Myr while others have been disrupted. Such a trend depending on $M_{0}$
suggests that the stellar wind mass loss at the beginning is not the only
reason for the divergent evolution like that of the OB-min and OB-max models.

The values of $T_{\mathrm{rh}}$ are shown in the second row of panels in Fig.
7. Generally, the model with a short $T_{\mathrm{rh}}$ loses mass faster, but
it is not the case for the two models in the M1600 set with
$T_{\mathrm{rh}}=231.9$ Myr and $255.4$ Myr (hereafter named as R1 and R2
models, respectively). The R2 model has a large $T_{\mathrm{rh}}$ but
dissolute faster. In the meantime, some models in the M1600 set have larger
$T_{\mathrm{rh}}$ than those of the M800 set, but the former has already been
disrupted at the present day while the latter still survives. Thus, the
different evolution cannot be explained by $T_{\mathrm{rh}}$.

Instead, we find that the existence of BH binaries (BBHs) play an important
role. The total masses ($M_{\mathrm{BH}}$), the number of BHs at 32 Myr and
the evolution of semi-major axies of BBHs are shown in the bottom panels in
Fig. 7. The R2 model in the M1600 set has 4 BHs with
$M_{\mathrm{BH}}=36M_{\odot}$ while the R1 model has no BH. This can explain
why the R2 model loses mass faster. The decrease of $a_{\mathrm{BBH}}$ is
found in all models with BBHs in Fig. 7, which indicates the dynamical effect
of binary heating (e.g., Binney & Tremaine, 1987; Spitzer, 1987; Mackey et
al., 2008; Breen & Heggie, 2013; Wang, 2020). When light stars have close
encounters with the BBH, they gain a large kinetic energy and are easily
kicked out from the center of a star cluster while the separation of the BBH
shrinks. Such a process can drive the expansion of the stellar halo and
accelerate the escaper generation. Thus, a star cluster containing more
massive BBHs tends to lose mass faster. In summary, $M_{\mathrm{OB}}$ not only
affect the strength of stellar wind at the early phase of the star cluster,
but also influences the formation of BBHs and the long-term evolution of the
system.

###  3.5 The structure of tidal streams

Figure 8: The comparison for the morphology, the density and the average mass
of stars for the five models with different $M_{0}$ and the minimum $\Delta$
referring to the OB-min model in the M800R2 set. For each model, the three
sub-panels show the number density map of stars in the spherical coordinate
system of the Galactocentric reference frame, the histogram of star counting
($N_{\mathrm{s}}$) and the average mass of stars ($\langle
m_{\mathrm{s}}\rangle$) along the longitude, respectively. The compacted
objects (WDs, NSs and BHs) are excluded.

It is difficult to reconstruct the initial condition from the present-day
property of a star cluster because of the random value of $M_{\mathrm{OB}}$.
However, the tidal stream, which records the history of escapers, may provide
an additional constraint. To investigate this, we select the OB-min model in
the M800R2 set as the reference, and then, for each $M_{0}$, we find the most
similar model with the minimum $\Delta$ at the present day. Totally, five
models with different $M_{0}$ are collected.

The morphology, the number density distribution and the average mass
distribution of these (five) models at the present day are shown in Fig. 8.
For a better comparison, we correct the central position of each model to be
located at the observed position of the Hyades. The morphology (density map)
is shown in the Galactocentric reference frame with the spherical coordinate
system. In the longitude–latitude plot, the stream distributes in a regular
shape. Thus, it is easy to count the number of stars along the stream (the
longitude) and detect the sub-structures like the overdensity generated by the
epicyclic motion (Küpper et al., 2008).

For all models, the lengths of the streams are similar, and the overdensities
appear at the 179 $\deg$ and 181 $\deg$ along the longitude. Küpper et al.
(2008) showed that for a point-mass potential with a circular orbit, the
distance of the overdensity to the center can be described as

$R_{\mathrm{den}}=R_{\mathrm{tid}}\frac{4\pi\Omega(4\Omega^{2}-\kappa^{2})}{\kappa^{3}}$
(7)

where $\kappa$ and $\Omega$ are epicyclic and circular frequencies,
respectively. Based on the observational data of the solar neighborhood,
$\kappa/\Omega\approx 1.35$ (Binney & Tremaine, 1987). Thus,
$R_{\mathrm{den}}\approx 3.54\pi R_{\mathrm{tid}}$. Since $R_{\mathrm{tid}}$
evolves, we roughly estimate the range of $R_{\mathrm{den}}$ (the grey region
in Fig. 8) by using the minimum and the maximum values of $R_{\mathrm{tid}}$
during the simulation. To calculate the longitude of $R_{\mathrm{den}}$,
$R_{\mathrm{gal}}$ is needed. From Figure 1, $R_{\mathrm{gal}}$ varies between
$7$ to $9$ kpc. We simply use the present-day value to obtain a rough
estimation. The predicted positions of the overdensities are consistent with
those of the models. Since $R_{\mathrm{den}}$ is independent of $M_{0}$, the
morphology of the tidal streams does not show a pronounced difference.

We also analyze the average mass ($\langle m_{\mathrm{s}}\rangle$) along the
longitude. There is a clear sign of mass segregation that $\langle
m_{\mathrm{s}}\rangle$ is much larger inside the star cluster compared to that
along the stream. But there is also no pronounced difference depending on
$M_{0}$.

Although the structures and the mass functions of the five tidal streams are
similar, we find that the star counts inside and outside $R_{\mathrm{tid}}$
strongly depend on $M_{0}$. Table 4 lists the total masses and the numbers of
stars with apparent G-band magnitudes $<20$ inside and outside
$R_{\mathrm{tid}}$ for these 5 models, respectively. All models have a similar
$M_{\mathrm{tid}}$ and $N_{\mathrm{tid}}$ while the low-$M_{0}$ model has a
small $M_{\mathrm{out}}$ and $N_{\mathrm{out}}$. As the high-$M_{0}$ model has
a large $M_{\mathrm{OB}}$, the dissolution of the star cluster is faster
because of stronger stellar wind mass loss and dynamical heating from BBHs.
Thus, the fraction of the number of stars inside $R_{\mathrm{tid}}$
($f_{\mathrm{tid}}=N_{\mathrm{tid}}/N$ where $N$ is the total number of stars)
can help to estimate $M_{0}$ and $M_{\mathrm{OB}}$ for an observed star
cluster.

From the Gaia EDR3 catalogue, 862 objects with a similar magnitude cutoff were
detected, 293 of which belong to the leading tail and 166 to the trailing
tail. All $N_{\mathrm{out}}$ in Table 4 are significantly higher than the
observed counts. But it is difficult to make a solid conclusion based on the
direct comparison of the $N$-body model to the observation. For Gaia data, a
strict selection is needed to remove contamination and obtain signal in the
tails, which affects completeness. With a looser constraint, the contamination
is higher. In addition, an unresolved binary is counted as one object and
large fraction of these are filtered out by quality cuts applied to the Gaia
data (Jerabkova et al., 2021). Due to the overlap of objects along the line-
of-sight and a higher fraction of binaries inside the star cluster, the
numbers of stars might be underestimated. Meanwhile, the background
contamination is expected to be higher along the tails. A better solution is
to use the $N$-body model as the reference to create mock observations
considering the Gaia uncertainty and background contamination. Then compare
the model with the observation to obtain a better constraint. For such a
purpose, the assumption of primordial binaries is also necessary. This will be
our future effort.

This result suggests one strategy to constrain the initial condition of a star
cluster. First, we find a group of models with different $M_{0}$ and
$M_{\mathrm{OB}}$ that can result in the same $R_{\mathrm{tid}}$,
$M_{\mathrm{tid}}$ and $R_{\mathrm{h}}$. Then, we use $f_{\mathrm{tid}}$ to
determine which combination of $M_{0}$ and $M_{\mathrm{OB}}$ best fits the
observational data. Because of the high performance of petar, it is
practically possible to generate a grid of star clusters models to achieve
this goal.

Table 3: The parameters of the model set shown in Fig. 8

Mode set | $M_{\mathrm{OB}}$ | $M_{\mathrm{tid}}$ | $M_{\mathrm{out}}$ | $N_{\mathrm{tid}}$ | $N_{\mathrm{out}}$
---|---|---|---|---|---
M800R2 | 27.0 | 196.3 | 471.9 | 186 | 735
M1000R1 | 207.2 | 196.5 | 518.2 | 184 | 761
M1200R05 | 165.5 | 196.1 | 658.6 | 173 | 958
M1400R05 | 348.8 | 198.5 | 727.1 | 202 | 973
M1600R05 | 456.9 | 196.2 | 857.0 | 180 | 1208
444The corresponding model set, the initial OB star mass ($M_{\mathrm{OB}}$),
the total mass inside $R_{\mathrm{tid}}$ ($M_{\mathrm{tid}}$) and outside
$R_{\mathrm{tid}}$ ($M_{\mathrm{out}}$), the number counts for stars with
apparent G-band magnitude $<20$ inside $R_{\mathrm{tid}}$ ($N_{\mathrm{tid}}$)
and outside $R_{\mathrm{tid}}$ ($N_{\mathrm{out}}$) for the five models shown
in Fig. 8. Their $R_{\mathrm{tid}}\approx 6.5$ pc and $R_{\mathrm{h}}\approx
3.0$ pc at 648 Myr. The mass unit is $M_{\odot}$.

##  4 Discussion and conclusion

In this work, we investigate how the content of OB stars affects the long-term
evolution of Hyades-like open star clusters and their tidal streams. The
models are designed to have the initial condition, including the shape of the
IMF and the density profile, to be identical. The random number of OB stars
can result in a significantly different dynamical evolution of the system. As
the example in Fig. 4 shows, the OB-min model with a small $M_{\mathrm{OB}}$
can survive as a star cluster until the present day, while the OB-max model
with a large $M_{\mathrm{OB}}$ has already dissolved 200 Myr ago.

The stellar wind mass loss of OB stars during the first 32 Myr and the BBH
dynamical heating both accelerate the mass loss of star clusters (Fig. 7) and
cause the difference. Thus, assuming the stochastic variation of
$M_{\mathrm{OB}}$ being the correct physical model, it is not possible to
constrain the initial condition (e.g. $M_{0}$ and $M_{\mathrm{OB}}$) of an
observed star cluster confidently by only checking the present-day properties
inside $R_{\mathrm{tid}}$ (e.g., $M_{\mathrm{tid}}$, $R_{\mathrm{h}}$ and
$R_{\mathrm{c}}$). Meanwhile, we find that the morphology (overdensity) and
average mass distribution along the tidal streams are independent on $M_{0}$
(see Fig. 8). However, the star counts inside and outside $R_{\mathrm{tid}}$
($f_{\mathrm{tid}}$; Table 4) can help to remove the degeneracy and help to
constrain the initial condition.

There are a few aspects that are not included in our analysis. We did not
consider the influence from the uncertainty of the ages in the observation. We
expect that the age may affect the determination of $M_{0}$ by using
$f_{\mathrm{tid}}$, but it may also influence the length of the tidal stream
and the position of the overdensity. Thus, we may disentangle the impacts from
the age and the $M_{\mathrm{OB}}$.

There are no primordial binaries in our models. Observations have shown that
massive stars in open clusters are most likely all in multiple systems (Sana
et al., 2012; Duchêne & Kraus, 2013; Moe & Di Stefano, 2017). The dynamical
interactions of these multiple systems can result in ejections and mergers
(e.g. Wang, Kroupa & Jerabkova, 2019). As a result, the numbers of retained OB
stars and formed BHs decrease. In addition, the primordial binaries can also
attend the binary heating at the early phase that affects the timescale of
core collapse (e.g. Heggie, Trenti, & Hut, 2006). All of these subsequently
affect the long-term dynamical evolution of clusters. In the future work, we
will also consider the impact from primordial binaries.

The gas expulsion process can also significantly affect the dynamics of the
star cluster during the gas embedded phase (e.g., Wang, Kroupa & Jerabkova,
2019; Fujii et al., 2021a, b) and the formation of tidal streams (Dinnbier &
Kroupa, 2020a, b). The UV radiation, stellar winds and supernovae from OB
stars all can drive the gas expulsion. Thus, the gas expulsion is also
sensitive to the stochastic nature of $M_{\mathrm{OB}}$. A stronger feedback
because of more OB stars may quench the star formation earlier, and thus, the
star formation efficiency becomes lower. As a result, the damage to the
stellar system is also stronger and a larger number of escapers is expected to
appear during the first few Myr. By considering gas expulsion, we expect an
even larger scatter of the evolution tracks due to the variation of
$M_{\mathrm{OB}}$.

As a caveat it is noted that we applied a simple Milky Way potential
(MWPotential2014). There is no time-dependent evolution of the potential, no
bar and no spiral arm. Since we are not aiming at obtaining the precise
morphology of the tidal stream for the Hyades, missing these components does
not significantly influence the major conclusion. Since galpy supports to
build up a more realistic Galactic potential, we will investigate such an
effect in the future work.

We assumed in this work that the formation of OB stars is purely statistical
implying a large variation of $M_{\mathrm{OB}}$. But if the star formation has
some degree of self-regulation such that the total number and mass of OB stars
depends on the initial mass of the cluster-forming gas cloud, the scatter of
$M_{\mathrm{OB}}$ would be much smaller than the case studied here. Then,
constraining the initial conditions of an observed open cluster would be much
easier. For example, if the star formation is highly self-regulated such that
the stars formed can be described by optimally sampling the IMF (Kroupa et
al., 2013), $M_{\mathrm{OB}}\approx 289.81M_{\odot}$ and
$m_{\mathrm{max}}\approx 53M_{\odot}$ for $M_{0}=1600\leavevmode\nobreak\
M_{\odot}$ and only a tiny scatter of $M_{\mathrm{OB}}$ is allowed. Therefore,
by invoking stochastic sampling from the IMF, our work provides the upper
limit on the uncertainty in constraining the initial conditions of a given
open star cluster.

###### Acknowledgements.

L.W. thanks the financial support from JSPS International Research Fellow
(School of Science, The university of Tokyo). We thank Jo Bovy for the help on
the implementation of galpy interface in petar and Qi Shu for the help on the
debugging of the code. TJ acknowledges support through the European Space
Agency fellowship programme.

## References

  * Aarseth, Henon, & Wielen (1974) Aarseth S. J., Henon M., Wielen R., 1974, A&A, 37, 183
  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations, Cambridge University Press
  * Angelo et al. (2021) Angelo, M. S., Corradi, W. J. B., Santos, J. F. C., et al. 2021, MNRAS, 500, 4338. doi:10.1093/mnras/staa3192
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33. doi:10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123. doi:10.3847/1538-3881/aabc4f
  * Bally et al. (1987) Bally, J., Langer, W. D., Stark, A. A., et al. 1987, ApJ, 312, L45. doi:10.1086/184817
  * Banerjee & Kroupa (2017) Banerjee S., Kroupa P., 2017, A&A, 597, A28. doi:10.1051/0004-6361/201526928
  * Banerjee et al. (2020) Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2020, A&A, 639, A41. doi:10.1051/0004-6361/201935332
  * Baumgardt & Makino (2003) Baumgardt H., Makino J., 2003, MNRAS, 340, 227
  * Baumgardt & Kroupa (2007) Baumgardt H., Kroupa P., 2007, MNRAS, 380, 1589. doi:10.1111/j.1365-2966.2007.12209.x
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Natur, 324, 446
  * Beccari et al. (2020) Beccari, G., Boffin, H. M. J., & Jerabkova, T. 2020, MNRAS, 491, 2205. doi:10.1093/mnras/stz3195
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217. doi:10.1088/0004-637X/714/2/1217
  * Belczynski et al. (2016) Belczynski K., et al., 2016, A&A, 594, A97
  * Bhattacharya et al. (2021) Bhattacharya, S., Agarwal, M., Rao, K. K., et al. 2021, MNRAS, 505, 1607. doi:10.1093/mnras/stab1404
  * Binney & Tremaine (1987) Binney J., Tremaine S., 1987, Princeton, NJ, Princeton University Press
  * Blaauw (1964) Blaauw, A. 1964, ARA&A, 2, 213. doi:10.1146/annurev.aa.02.090164.001241
  * Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  * Brown et al. (1994) Brown, A. G. A., de Geus, E. J., & de Zeeuw, P. T. 1994, A&A, 289, 101
  * Breen & Heggie (2013) Breen P. G., Heggie D. C., 2013, MNRAS, 432, 2779
  * Carpenter (2000) Carpenter, J. M. 2000, AJ, 120, 3139. doi:10.1086/316845
  * Cantat-Gaudin et al. (2018) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93. doi:10.1051/0004-6361/201833476
  * Cantat-Gaudin et al. (2019) Cantat-Gaudin, T., Krone-Martins, A., Sedaghat, N., et al. 2019, A&A, 624, A126. doi:10.1051/0004-6361/201834453
  * Casertano & Hut (1985) Casertano S., Hut P., 1985, ApJ, 298, 80
  * Chumak & Rastorguev (2006) Chumak Y. O., Rastorguev A. S., 2006, AstL, 32, 157. doi:10.1134/S1063773706030030
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792. doi:10.1086/318388
  * de Bruijne (1999) de Bruijne, J. H. J. 1999, MNRAS, 306, 381. doi:10.1046/j.1365-8711.1999.02643.x
  * Dias et al. (2021) Dias, W. S., Monteiro, H., Moitinho, A., et al. 2021, MNRAS, 504, 356. doi:10.1093/mnras/stab770
  * Dinnbier & Kroupa (2020a) Dinnbier F., Kroupa P., 2020a, A&A, 640, A84. doi:10.1051/0004-6361/201936570
  * Dinnbier & Kroupa (2020b) Dinnbier F., Kroupa P., 2020b, A&A, 640, A85. doi:10.1051/0004-6361/201936572
  * Dinnbier & Walch (2020) Dinnbier F., Walch S., 2020, MNRAS, 499, 748. doi:10.1093/mnras/staa2560
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Fukushige & Heggie (1995) Fukushige T., Heggie D. C., 1995, MNRAS, 276, 206
  * Fürnkranz et al. (2019) Fürnkranz, V., Meingast, S., & Alves, J. 2019, A&A, 624, L11. doi:10.1051/0004-6361/201935293
  * Ferreira et al. (2020) Ferreira, F. A., Corradi, W. J. B., Maia, F. F. S., et al. 2020, MNRAS, 496, 2021. doi:10.1093/mnras/staa1684
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91. doi:10.1088/0004-637X/749/1/91
  * Fujii et al. (2021a) Fujii, M. S., Saitoh, T. R., Wang, L., et al. 2021a, PASJ. doi:10.1093/pasj/psab037
  * Fujii et al. (2021b) Fujii, M. S., Saitoh, T. R., Hirai, Y., et al. 2021b, PASJ. doi:10.1093/pasj/psab061
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018, A&A, 616, A10. doi:10.1051/0004-6361/201832843
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1. doi:10.1051/0004-6361/202039657
  * Genzel & Stutzki (1989) Genzel, R. & Stutzki, J. 1989, ARA&A, 27, 41. doi:10.1146/annurev.aa.27.090189.000353
  * Giersz & Heggie (1996) Giersz M., Heggie D. C., 1996, MNRAS, 279, 1037
  * Castro-Ginard et al. (2018) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2018, A&A, 618, A59. doi:10.1051/0004-6361/201833390
  * Castro-Ginard et al. (2020) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2020, A&A, 635, A45. doi:10.1051/0004-6361/201937386
  * Godoy-Rivera et al. (2021) Godoy-Rivera, D., Pinsonneault, M. H., & Rebull, L. M. 2021, arXiv:2101.01183
  * Goodwin & Bastian (2006) Goodwin S. P., Bastian N., 2006, MNRAS, 373, 752. doi:10.1111/j.1365-2966.2006.11078.x
  * Goodman, Heggie, & Hut (1993) Goodman J., Heggie D. C., Hut P., 1993, ApJ, 415, 715. doi:10.1086/173196
  * González-Samaniego & Vazquez-Semadeni (2020) González-Samaniego A., Vazquez-Semadeni E., 2020, MNRAS, 499, 668. doi:10.1093/mnras/staa2921
  * Großschedl et al. (2021) Großschedl, J. E., Alves, J., Meingast, S., et al. 2021, A&A, 647, A91. doi:10.1051/0004-6361/202038913
  * He et al. (2021) He, Z.-H., Xu, Y., Hao, C.-J., et al. 2021, Research in Astronomy and Astrophysics, 21, 093. doi:10.1088/1674-4527/21/4/93
  * Heggie (1996) Heggie D. C., 1996, IAUS, 174, 131
  * Heggie (1988) Heggie D. C., 1988, ltdb.conf, 329
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, gmbp.book
  * Heggie, Trenti, & Hut (2006) Heggie D. C., Trenti M., Hut P., 2006, MNRAS, 368, 677. doi:10.1111/j.1365-2966.2006.10122.x
  * Hirai et al. (2021) Hirai, Y., Fujii, M. S., & Saitoh, T. R. 2021, PASJ. doi:10.1093/pasj/psab038
  * Hurley, Pols, & Tout (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley, Tout, & Pols (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Iwasawa et al. (2016) Iwasawa M., Tanikawa A., Hosono N., Nitadori K., Muranushi T., Makino J., 2016, PASJ, 68, 54
  * Iwasawa et al. (2020) Iwasawa M., Namekata D., Nitadori K., Nomura K., Wang L., Tsubouchi M., Makino J., 2020, PASJ, 72, 13
  * Jerabkova et al. (2019) Jerabkova, T., Boffin, H. M. J., Beccari, G., et al. 2019, MNRAS, 489, 4418. doi:10.1093/mnras/stz2315
  * Jerabkova et al. (2021) Jerabkova, T., Boffin, H. M. J., Beccari, G., et al. 2021, A&A, 647, A137. doi:10.1051/0004-6361/202039949
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa, Aarseth, & Hurley (2001) Kroupa P., Aarseth S., Hurley J., 2001, MNRAS, 321, 699. doi:10.1046/j.1365-8711.2001.04050.x
  * Kroupa et al. (2013) Kroupa P., Weidner C., Pflamm-Altenburg J., Thies I., Dabringhausen J., Marks M., Maschberger T., 2013, pss5.book, 115. doi:10.1007/978-94-007-5612-0_4
  * Kroupa et al. (2018) Kroupa P., Jeřábková T., Dinnbier F., Beccari G., Yan Z., 2018, A&A, 612, A74. doi:10.1051/0004-6361/201732151
  * Küpper et al. (2008) Küpper, A. H. W., MacLeod, A., & Heggie, D. C. 2008, MNRAS, 387, 1248. doi:10.1111/j.1365-2966.2008.13323.x
  * Küpper et al. (2010) Küpper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010, MNRAS, 401, 105. doi:10.1111/j.1365-2966.2009.15690.x
  * Kuhn et al. (2019) Kuhn, M. A., Hillenbrand, L. A., Sills, A., et al. 2019, ApJ, 870, 32. doi:10.3847/1538-4357/aaef8c
  * Kounkel et al. (2018) Kounkel, M., Covey, K., Suárez, G., et al. 2018, AJ, 156, 84. doi: 10.3847/1538-3881/aad1f1
  * Kounkel & Covey (2019) Kounkel, M. & Covey, K. 2019, AJ, 158, 122. doi:10.3847/1538-3881/ab339a
  * Kounkel (2020) Kounkel, M. 2020, ApJ, 902, 122. doi:10.3847/1538-4357/abb6e8
  * Liu & Pang (2019) Liu L., Pang X., 2019, ApJS, 245, 32. doi:10.3847/1538-4365/ab530a
  * Mackey et al. (2008) Mackey A. D., Wilkinson M. I., Davies M. B., Gilmore G. F., 2008, MNRAS, 386, 65. doi:10.1111/j.1365-2966.2008.13052.x
  * Madsen et al. (2006) Madsen, G. J., Reynolds, R. J., & Haffner, L. M. 2006, ApJ, 652, 401. doi:10.1086/508441
  * Marks & Kroupa (2012) Marks M., Kroupa P., 2012, A&A, 543, A8. doi:10.1051/0004-6361/201118231
  * Meingast, Alves & Fürnkranz (2019) Meingast S., Alves J., Fürnkranz V., 2019, A&A, 622, L13
  * Meingast, Alves, & Rottensteiner (2021) Meingast S., Alves J., Rottensteiner A., 2021, A&A, 645, A84. doi:10.1051/0004-6361/202038610
  * Miyamoto & Nagai (1975) Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15. doi:10.3847/1538-4365/aa6fb6
  * Monteiro & Dias (2019) Monteiro, H. & Dias, W. S. 2019, MNRAS, 487, 2385. doi:10.1093/mnras/stz1455
  * Monteiro et al. (2020) Monteiro, H., Dias, W. S., Moitinho, A., et al. 2020, MNRAS, 499, 1874. doi:10.1093/mnras/staa2983
  * Navarro et al. (1995) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1995, MNRAS, 275, 720. doi:10.1093/mnras/275.3.720
  * Namekata et al. (2018) Namekata D., Iwasawa M., Nitadori K., Tanikawa A., Muranushi T., Wang L., Hosono N., et al., 2018, PASJ, 70, 70. doi:10.1093/pasj/psy062
  * Ochsendorf et al. (2015) Ochsendorf, B. B., Brown, A. G. A., Bally, J., et al. 2015, ApJ, 808, 111. doi:10.1088/0004-637X/808/2/111
  * Oshino, Funato, & Makino (2011) Oshino S., Funato Y., Makino J., 2011, PASJ, 63, 881
  * O’dell et al. (1967) O’dell, C. R., York, D. G., & Henize, K. G. 1967, ApJ, 150, 835. doi:10.1086/149386
  * O’Dell et al. (2011) O’Dell, C. R., Ferland, G. J., Porter, R. L., et al. 2011, ApJ, 733, 9. doi:10.1088/0004-637X/733/1/9
  * Pang et al. (2020) Pang X., Li Y., Tang S.-Y., Pasquato M., Kouwenhoven M. B. N., 2020, ApJL, 900, L4. doi:10.3847/2041-8213/abad28
  * Pang et al. (2021) Pang X., Li Y., Yu Z., Tang S.-Y., Dinnbier F., Kroupa P., Pasquato M., et al., 2021, ApJ, 912, 162. doi:10.3847/1538-4357/abeaac
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Röser & Schilbach (2019a) Röser S., Schilbach E., 2019, A&A, 627, A4. doi:10.1051/0004-6361/201935502
  * Röser, Schilbach, & Goldman (2019b) Röser S., Schilbach E., Goldman B., 2019, A&A, 621, L2. doi:10.1051/0004-6361/201834608
  * Röser et al. (2011) Röser S., Schilbach E., Piskunov A. E., Kharchenko N. V., Scholz R.-D., 2011, A&A, 531, A92. doi:10.1051/0004-6361/201116948
  * Sim et al. (2019) Sim, G., Lee, S. H., Ann, H. B., et al. 2019, Journal of Korean Astronomical Society, 52, 145. doi:10.5303/JKAS.2019.52.5.145
  * Strömberg (1939) Strömberg, G. 1939, Popular Astronomy, 47, 172
  * Sana et al. (2012) Sana H., et al., 2012, Sci, 337, 444
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters, Princeton University Press
  * Tang et al. (2019) Tang S.-Y., Pang X., Yuan Z., Chen W. P., Hong J., Goldman B., Just A., et al., 2019, ApJ, 877, 12. doi:10.3847/1538-4357/ab13b0
  * van Leeuwen (2009) van Leeuwen, F. 2009, A&A, 497, 209. doi:10.1051/0004-6361/200811382
  * Wang, Kroupa & Jerabkova (2019) Wang L., Kroupa P., Jerabkova T., 2019, MNRAS, 484, 1843
  * Wang (2020) Wang L., 2020a, MNRAS, 491, 2413
  * Wang, Nitadori & Makino (2020) Wang L., Nitadori K., Makino J., 2020b, MNRAS, 493, 3398
  * Wang et al. (2020) Wang L., Iwasawa M., Nitadori K., Makino J., 2020c, MNRAS, 497, 536
  * Weidner & Kroupa (2006) Weidner C., Kroupa P., 2006, MNRAS, 365, 1333
  * Weidner, Kroupa, & Pflamm-Altenburg (2013) Weidner C., Kroupa P., Pflamm-Altenburg J., 2013, MNRAS, 434, 84
  * Yan, Jerabkova, & Kroupa (2017) Yan Z., Jerabkova T., Kroupa P., 2017, A&A, 607, A126
  * Zhang et al. (2020) Zhang, Y., Tang, S.-Y., Chen, W. P., et al. 2020, ApJ, 889, 99. doi:10.3847/1538-4357/ab63d4
  * Zhong et al. (2020) Zhong J., Chen L., Wu D., Li L., Bai L., Hou J., 2020, A&A, 640, A127. doi:10.1051/0004-6361/201937131
  * Reino et al. (2018) Reino, S., de Bruijne, J., Zari, E., et al. 2018, MNRAS, 477, 3197. doi:10.1093/mnras/sty793

##  Appendix A petar \- galpy

galpy (Bovy, 2015) is the state-of-the-art code to modelling the orbits of
stars in different kinds of galactic potentials. The online document provides
rich information about the usage of the
code555https://docs.galpy.org/en/v1.7.0/. The code is written mainly in the
Python programme language while it also provides the c programme language for
coupling to the $N$-body codes with a high computing performance. A interface
has been implemented in petar so that it is possible to trace the formation
and evolution of tidal streams. Table 4 lists all names of galpy potentials
(version 1.7.0) that can be accessed by petar. The MWPotential and
MWPotential2014 are wrapped potentials (see Bovy, 2015). The petar.galpy.help
tool provides the basic description of each potential type. It is recommended
to use the Python interface of galpy to configure the potential, and then,
save the parameters and use petar to read them to initialize the potential in
$N$-body simulations.

Table 4: The supported potentials in the petar-galpy interface Burkert | DehnenBar
---|---
DehnenCoreSpherical | DehnenSpherical
DiskSCF | DoubleExponentialDisk
FlattenedPower | Hernquist
HomogeneousSphere | Isochrone
Jaffe | Kepler
King | KuzminDisk
KuzminKutuzovStaeckel | LogarithmicHalo
MN3ExponentialDisk | MWPotential
MWPotential2014 | MiyamotoNagai
NFW | PerfectEllipsoid
Plummer | PowerSpherical
PowerSpherical | PowerTriaxial
PseudoIsothermal | SCF
SoftenedNeedleBar | SpiralArms
TriaxialGaussian | TriaxialHernquist
TriaxialJaffe | TriaxialNFW

When galpy is switched on, petar uses the Galactic center as the reference
frame without rotation. To keep the high digital precision for positions and
velocities of stars, the coordinate origin of the stellar system is not the
Galactic center but follows the motion of the center of the cluster. The
position and velocity of center referring to the Galactic center are saved as
the offset values in the header of the snapshots. Using these offset values,
users can easily obtain the coordinate frame referring to the Galactic center.
The default way to determine the center is to use the long-range potential
weighted average of positions and velocities of all objects. This may not be
the best way. It is challenging to determine the center when the star cluster
is close to the complete disruption. Thus, users may need to redetermine the
center in the late phase of the evolution by using the snapshots.

The petar-galpy interface supports to add two types of potentials: the co-
moving potentials and the fixed central potentials. The former co-moves with
the center of the stellar system. The latter has a fixed position referring to
the Galactic center. Either type can combine arbitrary time-dependent
potentials listed in Table 4. Thus, it is flexible to study a variety topics.
For example, users can setup a star cluster move in the Milky-way potential,
with a gas halo or a dark matter halo surrounding the cluster. The gas
potential can be reduced by time to represent the gas expulsion. The details
to setup the potentials can be found in the online README of petar666GitHub
page: https://github.com/lwang-astro/PeTar and from the commander ’petar -h’.

##  Appendix B Comparison with observation

To construct a bridge to connect the $N$-body simulations and observations, we
have implemented a group of data analysis tool written in Python3 (see online
manual for details), Here we briefly describe the functions of this tool.

###  B.1 The transformation of reference frames and coordinate systems

petar equipped with galpy uses the Galactocentric reference frame and the
Cartesian coordinate system in the simulation. The result cannot be directly
compared with the observational data. The petar data analysis tool includes an
interface to transform the snapshot from a $N$-body simulation to the data
type of astropy.coordinates.skycoord (Astropy Collaboration et al., 2013,
2018). skycoord is a flexible data type that supports an easy way to transform
reference freams and coordinates systems for a group of data with 3D positions
and velocities. Fig. 4 and 8 show the snapshot of $N$-body simulations in the
ICRS frame (RA, Dec, proper motions) and in the Galactocentric frame
(longitude and lattitude), respectively. These can be easily generated by
using this interface.

###  B.2 Convergent point check for proper motions

One of the very useful way of identification of (nearly) co-moving stars in
the Gaia catalog, that is using only proper motions and not the mostly absent
radial velocity values, is so-called convergent point method. Originally, the
convergent point (CP) method has been used to constrain distance to nearby
star-clusters (Strömberg, 1939; de Bruijne, 1999; van Leeuwen, 2009, and
references therein). The innovative usage of this method to identify co-moving
stars instead has been first applied by (Röser, Schilbach, & Goldman, 2019b)
and later discussed in detail by (Jerabkova et al., 2021). The CP method is
correcting measured proper motions for on-the-sky projection effects that can
be substantial for close-by extended objects. For more details we refer reader
to Jerabkova et al. (2021); van Leeuwen (2009); Röser, Schilbach, & Goldman
(2019b).

The mentioned PeTar script computed the CP diagram and allows thus compare
simulations with the Gaia and offers thus more possiblities how to search for
co-moving stars in the Gaia catalog.

###  B.3 Mock observational errors

In the current version of the tool it is possible to generate expected Gaia
EDR3 uncertainties based on stellar magnitudes computed based on respective
stellar mass. The detail description will be provided in the upcoming work
with demonstration cases.


#  High-ordered spectral characterizations of graphs

Lixiang Chen  Lizhu Sun  Changjiang Bu<EMAIL_ADDRESS>College of
Mathematical Sciences, Harbin Engineering University, Harbin 150001, PR China

###### Abstract

The spectrum of the $k$-power hypergraph of a graph $G$ is called the
$k$-ordered spectrum of $G$. If graphs $G_{1}$ and $G_{2}$ have same
$k$-ordered spectrum for all positive integer $k\geq 2$, $G_{1}$ and $G_{2}$
are said to be high-ordered cospectral. If all graphs who are high-ordered
cospectral with the graph $G$ are isomorphic to $G$, we say that $G$ is
determined by the high-ordered spectrum. In this paper, we use the high-
ordered spectrum of graphs to study graph isomorphism and show that all
Smith’s graphs are determined by the high-ordered spectrum. And we give
infinitely many pairs of trees with same spectrum but different high-ordered
spectrum by high-ordered cospectral invariants of trees, it means that we can
determine that these cospectral trees are not isomorphism by the high-ordered
spectrum.

###### keywords:

graph isomorphism, high-ordered spectrum, spectral characterization,
cospectral invariants.
_AMS classification(2020):_ 05C65, 05C50.

##  1 Introduction

If there was a “Holy Grail” in graph theory, it would be a practical test for
graph isomorphism [17]. In 1956, Günthard and Primas used the spectrum of
graphs to determine whether the graphs are isomorphic [11]. In 1957, Collatz
and Sinogowitz presented a pair of non-isomorphic cospectral tree [24].
Subsequently, many non-isomorphic cospectral graphs (i.e., graphs having the
same spectrum) were found [8, 10, 19]. For instance, two graphs shown in
Figure 1 are non-isomorphic cospectral graphs given by Cvetković in 1971 [8].
This pair of graphs is usually called the Saltire pair because the two graphs
superposed give the Scottish flag: Saltire [22].

(a) $G_{1}$

(b) $G_{2}$

Figure 1: The Saltire pair

If all graphs who are cospectral with a graph $G$ are isomorphic to $G$, the
graph $G$ is said to be determined by the spectrum (DS for short). Until now,
the known DS graphs are very special and most of the techniques involved in
proving them to be DS cannot be applied to general graphs [9, 22, 23, 25].

In this paper, we use the high-ordered spectrum of graphs to study the high-
ordered spectral characterizations of graphs. Since two graphs are isomorphic
if and only if their power hypergraphs are isomorphic and the spectrum of the
power hypergraph contains more graph structure information, the high-ordered
spectrum of graphs can determine more graphs than the spectrum of graphs. For
example, the Saltire pair can not be determined by the spectrum but can be
determined by the high-ordered spectrum (shown in Section 3).

Let the $k$-power hypergraph $G^{(k)}$ be the $k$-uniform hypergraph that is
obtained by adding $k-2$ new vertices to each edge of $G$. The spectrum of
$G^{(k)}$ is called the $k$-ordered spectrum of $G$. If graphs $G_{1}$ and
$G_{2}$ have same $k$-ordered spectrum for all positive integer $k\geq 2$,
$G_{1}$ and $G_{2}$ are said to be _high-ordered cospectral_. If all graphs
who are high-ordered cospectral with the graph $G$ are isomorphic to $G$, we
say that $G$ is _determined by the high-ordered spectrum_ (DHS for short). If
graphs which have same $k$-ordered spectrum with the graph $G$ are isomorphic
to $G$, we say that $G$ is determined by the $k$-ordered spectrum. Obviously,
the graphs determined by the $k$-ordered spectrum are DHS. The graphs
determined by the spectrum ($2$-ordered spectrum) are DHS.

In 1970, Smith classified all connected graphs with spectral radius at most
$2$ [20], which usually are called “Smith’s graphs”. In 2009, Van Dam and
Haemers showed that not all Smith’s graphs are determined by the spectrum
[23]. We show that all Smith’s graphs are determined by the high-ordered
spectrum in this paper. In 1973, Schwenk gave a useful method to construct
cospectral trees and proved his famous conclusion:“Almost all trees are not
DS” [19]. We show that the high-ordered spectrum of infinitely many pairs of
cospectral trees constructed by Schwenk’s method are different.

This paper is organized as follows. In Section 2, we introduce the spectrum of
hypergraph and some lemmas. In Section 3, we show that Smith’s graphs are
determined by the high-ordered spectrum. In Section 4, we give some high-
ordered cospectral invariants about the number of some subtrees. As their
applications, we give infinitely many pairs of cospectral trees with different
high-ordered spectrum. In Section 5, conclusions and future work are given.

##  2 Preliminaries

In this section, we introduce the spectrum of hypergraph and some lemmas which
are important to the subsequent results in the paper. For a positive integer
$n$, let $\left[n\right]=\left\\{{1,\ldots,n}\right\\}$. A $k$-order
$n$-dimension complex tensor $T=\left({{t_{{i_{1}}\cdots{i_{k}}}}}\right)$ is
a multidimensional array with $n^{k}$ entries on complex number field
$\mathbb{C}$, where ${i_{j}}\in\left[n\right]$, $j=1,\ldots,k$. Denote the set
of $n$-dimension complex vectors by $\mathbb{C}^{n}$. For
$x={\left({{x_{1}},\ldots,{x_{n}}}\right)^{\rm{T}}}\in{\mathbb{C}^{n}}$,
$T{x^{k-1}}$ denotes a vector in $\mathbb{C}^{n}$ whose $i$-th component is

${\left({T{x^{k-1}}}\right)_{i}}=\sum\limits_{{i_{2}},\ldots,{i_{k}}=1}^{n}{{t_{i{i_{2}}\cdots{i_{k}}}}{x_{{i_{2}}}}\cdots{x_{{i_{k}}}}}.$

If there exist $\lambda\in\mathbb{C}$ and a nonzero vector
$x=(x_{1},x_{2},\ldots,x_{n})^{\mathrm{T}}\in{\mathbb{C}^{n}}$ such that
$T{x^{k-1}}=\lambda{x^{\left[{k-1}\right]}}$, then $\lambda$ is called an
_eigenvalue_ of $T$ and $x$ is an _eigenvector_ of $T$ corresponding to
$\lambda$, where
${x^{\left[{k-1}\right]}}={\left({x_{1}^{k-1},\ldots,x_{n}^{k-1}}\right)^{\rm{T}}}$
[16, 18].

A hypergraph $H=(V(H),E(H))$ is called _$k$ -uniform_ if each edge of $H$
contains exactly $k$ vertices. For a $k$-uniform hypergraph $H$ with $n$
vertices, its _adjacency tensor_ ${A}_{H}=(a_{i_{1}i_{2}\ldots i_{k}})$ is a
$k$-order $n$-dimension tensor, where

${a_{{i_{1}}{i_{2}}\ldots{i_{k}}}}=\left\\{\begin{array}[]{l}\frac{1}{{\left({k-1}\right)!}},{\kern
37.0pt}\mathrm{if}{\kern
2.0pt}{\left\\{{{i_{1}},{i_{2}},\ldots,{i_{k}}}\right\\}\in{E(H)}},\\\
0,{\kern 57.0pt}\mathrm{otherwise}.\\\ \end{array}\right.$

All the eigenvalues of ${A}_{H}$ are called the spectrum of the hypergraph $H$
[6]. When $k=2$, ${A}_{H}$ is the adjacency matrix of the graph $H$.

The $k$-ordered spectrum of graph $G$ is the spectrum of power hypergraph
$G^{(k)}$. In [26], Zhou et al. show that eigenvalues of $G^{(k)}$ can be
obtained from the eigenvalues of $G$. In [2], the authors gave all the
eigenvalues of $G^{(k)}$ without counting multiplicity as following.

######  Lemma 2.1.

[2] Let $G^{(k)}$ be the $k$-power hypergraph of a graph $G$.
(1) When $k=3$, $\lambda$ is an eigenvalue of $G^{(3)}$ if and only if there
is an induced subgraph with an eigenvalue $\beta$ such that
$\beta^{2}=\lambda^{k}$.
(2) When $k>3$, $\lambda$ is an eigenvalue of $G^{(k)}$ if and only if there
is a subgraph with an eigenvalue $\beta$ such that $\beta^{2}=\lambda^{k}$.

In [3], the authors gave the characteristic polynomial and then give all
distinct eigenvalues of power hyperpaths.

######  Lemma 2.2.

[3] All distinct elements in the set

$\left\\{{\left({2\cos\frac{\pi}{{j+1}}t}\right)^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}}:j\in[n],t\in[j],\theta\in[k]\right\\}$

are all distinct eigenvalues of the $k$-uniform hyperpath $P^{(k)}_{n}$, where
$\mathbf{i}=\sqrt{-1}$.

Since the connected subgraphs of the cycle $C_{n}$ are $C_{n}$ and $P_{j}$ for
$j\in[n-1]$, we obtain all distinct eigenvalues of $C^{(k)}_{n}$ from all
eigenvalues of $C_{n}$ and $P_{j}$ by Lemma 2.1.

######  Lemma 2.3.

All distinct elements in the set

$\left\\{{\left({2\cos\frac{\pi}{{j+1}}t}\right)^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}},{\left({2\cos\frac{2\pi}{{n}}r}\right)^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}}:j\in[n-1],t\in[j],r\in[n],\theta\in[k]\right\\}$

are all distinct eigenvalues of the $k$-power hypercycle $C^{(k)}_{n}$ (for
$k>3$), where $\mathbf{i}=\sqrt{-1}$.

###### Proof.

The connected subgraphs of cycle $C_{n}$ are $C_{n}$ and $P_{j}$ for
$j\in[n-1]$. By Lemma 2.1 and 2.2, we know that all distinct numbers in the
set

$\left\\{{\left({2\cos\frac{\pi}{{j+1}}t}\right)^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}}:j\in[n-1],t\in[j],\theta\in[k]\right\\}\bigcup\left\\{{\lambda^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}}:\lambda\in\sigma(C_{n}),\theta\in[k]\right\\}$

are all distinct eigenvalues of $C^{(k)}_{n}$, where $\sigma(C_{n})$ is the
set of all eigenvalues of $C_{n}$. Since $\sigma(C_{n})=\\{2\cos{\frac{2\pi
r}{n}}:r=1,2,\ldots,n\\}$ ( Page 72 in [9]), we know that all distinct numbers
in the set

$\left\\{{\left({2\cos\frac{\pi}{{j+1}}t}\right)^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}},{\left({2\cos\frac{2\pi}{{n}}r}\right)^{\frac{2}{k}}}{e^{\mathbf{i}\frac{{2\pi}}{k}\theta}}:j\in[n-1],t\in[j],r\in[n],\theta\in[k]\right\\}$

are all distinct eigenvalues of $C^{(k)}_{n}$. ∎

The $d$-th order spectral moment $\mathrm{S}_{d}(G)$ of a graph $G$ is the sum
of $d$-th power of all eigenvalues of $G$. Two graphs are cospectral if and
only if their $d$-th order spectral moments are equal [23]. The spectral
moment of graph is an important parameter in the topic of graph spectral
characterizations [1, 7, 15, 22]. Let $\widehat{G}$ be a connected subgraph of
$G$. We use $N_{G}(\widehat{G})$ to denote the number of subgraphs of $G$
isomorphic to $\widehat{G}$, $c_{d}(\widehat{G})$ to denote the number of
closed walks with length $d$ in the graph $\widehat{G}$ running through all
the edges at least once. The $d$-th order spectral moment $\mathrm{S}_{d}(G)$
can be represented as a linear combination of the number of connected
subgraphs [1, 7, 9], i.e., ${\mathrm{S}_{d}}(G)=\sum\nolimits_{\widehat{G}\in
G(d)}{{c_{d}}(\widehat{G}){N_{G}}(\widehat{G})}$, where $G(d)$ is the set of
connected subgraphs of $G$ with at most $d$ edges. The coefficient
$c_{d}(\widehat{G})$ is called the $d$-th order spectral moment coefficient of
$\widehat{G}$. In [4], the authors gave a formula for the $d$-th order
spectral moment coefficients of trees.

######  Lemma 2.4.

[4] The $d$-th order spectral moment of the tree $T$ is

${\mathrm{S}_{d}}\left(T\right)=\left\\{\begin{array}[]{l}\sum\limits_{m=1}^{{\frac{d}{2}}}{{\sum\limits_{\widehat{T}\in{\mathbf{T}({m})}}{c_{d}(\widehat{T})}}N_{T}(\widehat{T})},{\kern
16.0pt}2\mid d,\\\ 0,{\kern 120.0pt}2\nmid d,\\\ \end{array}\right.$

where ${\mathbf{T}({m})}$ is the set of subtrees of $T$ with $m$ edges. The
$d$-th order spectral moment coefficient of the subtree $\widehat{T}$ is

$\displaystyle
c_{d}({\widehat{T}})=\left\\{\begin{array}[]{l}{d}{\sum\limits_{{\sum_{e\in
E({\widehat{T}})}}w(e)=\frac{d}{2}}}\left({\prod\limits_{e\in
E({\widehat{T}})}}w{(e)}{\prod_{v\in
V({\widehat{T}})}}\frac{{({d_{v}}-1)!}}{{{r_{v}}}}\right),{\kern 16.0pt}2\mid
d,\\\ 0,{\kern 228.0pt}2\nmid d,\\\ \end{array}\right.$ (2.3)

where $w(e)$ is a positive integer weight of $e\in E(\widehat{T})$,
${d_{v}}={\sum_{e\in E_{v}(\widehat{{T}})}}{w(e)}$,
${r_{v}}=\prod\nolimits_{{e\in E_{v}(\widehat{{T}})}}{{w(e)}!}$.

The $d$-th order spectral moment $\mathrm{S}_{d}(H)$ of a $k$-uniform
hypergraph $H$ is the sum of $d$-th power of all eigenvalues of $H$, i.e.,
$\mathrm{S}_{d}(H)=\sum_{\lambda\in\sigma(H)}{\lambda^{d}}$, where $\sigma(H)$
is the spectrum of $H$. Two hypergraphs are cospectral if and only if their
$d$-th order spectral moments are equal for all $d\geq 1$ [5]. Clark and
Cooper expressed the characteristic polynomial coefficients of a uniform
hypergraph $H$ by the spectral moments of $H$ and gave the “Harary-Sachs
Theorem” of hypergraphs [5]. In [4], the authors expressed the spectral moment
of power hypertree by the number of subtrees as following.

######  Lemma 2.5.

[4] Let $T^{(k)}$ be the $k$-power hypertree of a tree $T$. Let
$c_{i}({\widehat{T}})$ denote the $i$-th order spectral moment coefficient of
the subtree ${\widehat{T}}$. Then the $d$-th spectral moment of $T^{(k)}$ is

$\displaystyle{\mathrm{S}_{d}}\left(T^{(k)}\right)=\left\\{\begin{array}[]{l}\sum\limits_{m=1}^{{\frac{d}{k}}}\frac{1}{2}{{(k-1)}^{\left({|E(T)|-m}\right)(k-1)}}{k^{m(k-2)+1}}{\sum\limits_{\widehat{T}\in{\mathbf{T}(m)}}{c_{\frac{2d}{k}}(\widehat{T})}}N_{T}(\widehat{T}),{\kern
6.0pt}k\mid d,\\\ 0,{\kern 280.0pt}k\nmid d,\\\ \end{array}\right.$

where $\mathbf{T}(m)$ is the set of subtrees of $T$ with $m$ edges.

##  3 All Smith’s graphs are DHS

Connected graphs with spectral radius at most $2$, as shown as in Figure 2,
usually are called “Smith’s graphs”, since Smith classified all connected
graphs with spectral radius at most $2$ in 1970 [20]. In 2009, Van Dam and
Haemers showed that not all Smith’s graphs are determined by the spectrum
[23]. In this section, we show that all Smith’s graphs are DHS.

Smith’s graphs are determined by the spectrum except for the graphs
$\widetilde{D}_{n}$ and $\widetilde{E}_{6}$ [23], then all Smith’s graphs are
DHS if and only if $\widetilde{D}_{n}$ and $\widetilde{E}_{6}$ are DHS. In
order to prove that $\widetilde{D}_{n}$ and $\widetilde{E}_{6}$ are DHS, we
only need to give all non-isomorphic cospectral graphs of $\widetilde{D}_{n}$
(resp. $\widetilde{E}_{6}$) and then prove that these cospectral graphs are
not high-ordered cospectral with $\widetilde{D}_{n}$ (resp.
$\widetilde{E}_{6}$).

(a) $P_{n}$

(b) $D_{n}$

(c) $E_{6}$

(d) $E_{7}$

(e) $E_{8}$

(f) $C_{n}$

(g) $\widetilde{D}_{n}$

(h) $\widetilde{E}_{6}$

(i) $\widetilde{E}_{7}$

(j) $\widetilde{E}_{8}$

Figure 2: Smith’s graphs

We use cospectral invariants of graphs to give all non-isomorphic cospectral
graphs of graph $\widetilde{D}_{n}$ (resp. $\widetilde{E}_{6}$). In [7],
Cvetković and Rowlinson gave the following cospectral invariant of graphs
about the number of subgraphs.

######  Lemma 3.6.

[7] If graphs $G$ and $G^{*}$ are cospectral, then
$N_{P_{3}}(G)+2N_{C_{4}}(G)=N_{P_{3}}(G^{*})+2N_{C_{4}}(G^{*})$.

We give all non-isomorphic cospectral graphs of $\widetilde{D}_{n}$ (resp.
$\widetilde{E}_{6}$) by the above cospectral invariant.

######  Lemma 3.7.

Graph $G$ is a non-isomorphic cospectral graphs of $\widetilde{D}_{n}$ (resp.
$\widetilde{E}_{6}$) if and only if $G=C_{4}+P_{n}$ (resp. $C_{6}+K_{1}$).

###### Proof.

It is well know that $C_{4}+P_{n}$ is a non-isomorphic cospectral graphs of
$\widetilde{D}_{n}$ ( Page 77 in [9]). Let $G$ be a non-isomorphic cospectral
graph of $\widetilde{D}_{n}$. We will show that $G=C_{4}+P_{n}$.

Since $G$ is cospectral with $\widetilde{D}_{n}$ and $G$ is not a Smith’s
graph, we know that the spectral radius of $G$ is equal to $2$ and $G$ is not
connected. Then the spectral radii of every connected components of $G$ are at
most $2$. Since connected graphs with spectral radius at most $2$ are cycles
or trees, every connected components of $G$ are cycles or trees. From
$|V(G)|=|V(\widetilde{D}_{n})|=n+4$ and $|E(G)|=|E(\widetilde{D}_{n})|=n+3$,
we get $|E(G)|=|V(G)|-1$. Then only one connected component of $G$ is a tree
and the other connected components are cycles. The algebraic multiplicity of
the eigenvalue $\lambda=2$ of $G$ (or $\widetilde{D}_{n}$) is 1 ( Page 77 in
[9]), it yields that only one connected component of $G$ is cycle and the
spectral radii of the other connected components are less than $2$. Then the
graph $G$ contains exactly two connected components, one of them is a cycle
and the other is a tree with spectral radius less than $2$.

We know that $G=C_{s}+T_{n+4-s}$, where $C_{s}$ is a cycle with $s$ vertices
and $T_{n+4-s}$ is a tree whose spectral radius is less than 2. Since $G$ is
cospectral with $\widetilde{D}_{n}$, we have
$N_{P_{3}}(G)+2N_{C_{4}}(G)=N_{P_{3}}(\widetilde{D}_{n})+2N_{C_{4}}(\widetilde{D}_{n})$
by Lemma 3.6. We know that
$N_{P_{3}}(\widetilde{D}_{n})+2N_{C_{4}}(\widetilde{D}_{n})=n+4$ from Figure
2. Then
$N_{P_{3}}(C_{s})+N_{P_{3}}(T_{n+4-s})+2N_{C_{4}}(C_{s})+2N_{C_{4}}(T_{n+4-s})=n+4$.
We have
$2N_{C_{4}}(C_{s})=n+4-N_{P_{3}}(C_{s})-N_{P_{3}}(T_{n+4-s})=n+4-s-N_{P_{3}}(T_{n+4-s})$.
Since the spectral radius of $T_{n+4-s}$ is less than 2, we know that
$T_{n+4-s}$ is isomorphic to one of $P_{n+4-s}$, $D_{n+2-n}$, $E_{6}$, $E_{7}$
and $E_{8}$. If $T_{n+4-s}=P_{n+4-s}$, we have $N_{P_{3}}(T_{n+4-s})=n+2-s$
when $n-s\geq-1$ and $N_{P_{3}}(T_{n+4-s})=0$ when $-3\leq n-s<-1$. When
$n-s\geq-1$, we have $N_{C_{4}}(C_{s})=1$, i.e., $s=4$. When $-3\leq n-s<-1$,
we have $1\leq 2N_{C_{4}}(C_{s})<3$. Then $N_{C_{4}}(C_{s})=1$, i.e., $s=4$.
In this case, $N_{C_{4}}(G)=1$. If $T_{n+4-s}$ is isomorphic to one of
$D_{n+2-n}$, $E_{6}$, $E_{7}$ and $E_{8}$, it yields that
$N_{P_{3}}(T_{n+4-s})=n+3-s$ from Figure 2. In this case, we get
$2N_{C_{4}}(C_{s})=1$. It contradicts the fact that $2N_{C_{4}}(C_{s})$ is
even. So we have $C_{s}=C_{4}$ and $T_{n+4-s}=P_{n}$.

We derive that $G=C_{4}+P_{n}$ from the above discussion, i.e., $G$ is a non-
isomorphic cospectral graphs of $\widetilde{D}_{n}$ if and only if
$G=C_{4}+P_{n}$. Similarly, we obtain that $G$ is a non-isomorphic cospectral
graphs of $\widetilde{E}_{6}$ if and only if $G=C_{6}+K_{1}$. ∎

Next, we prove that there is an eigenvalue $\lambda$ of the power hypergraph
$\widetilde{D}^{(k)}_{n}$ (resp. $\widetilde{E}^{(k)}_{6}$) but $\lambda$ is
not an eigenvalue of $(C_{4}+P_{n})^{(k)}$ (resp. $(C_{6}+K_{1})^{(k)}$).
Therefore, we obtain the main result in this section.

######  Theorem 3.8.

All Smith’s graphs are determined by the high-ordered spectrum.

###### Proof.

We will prove that $\widetilde{D}_{n}$ and $\widetilde{E}_{6}$ are DHS. By
Lemma 3.7, we know that graph $G$ is a non-isomorphic cospectral graphs of
$\widetilde{D}_{n}$ (resp. $\widetilde{E}_{6}$) if and only if $G=C_{4}+P_{n}$
(resp. $C_{6}+K_{1}$). Then we only need to prove that $\widetilde{D}_{n}$
(resp. $\widetilde{E}_{6}$ ) is not high-ordered cospectral with
$(C_{4}+P_{n})$ (resp. $C_{6}+K_{1}$ ).

Next, we prove that there is an eigenvalue $\lambda$ of the hypergraph
$\widetilde{D}^{(k)}_{n}$ such that $\lambda$ is not an eigenvalue of the
hypergraph $(C_{4}+P_{n})^{(k)}$. Since $D_{n}$ is an induced subgraph of
$\widetilde{D}_{n}$ and $2\cos{\frac{1}{2n+2}}\pi$ is an eigenvalue of $D_{n}$
(Page 77 in [9]), we know that $(2\cos{\frac{1}{2n+2}}\pi)^{\frac{2}{k}}$ is
an eigenvalue of $\widetilde{D}^{(k)}_{n}$ by Lemma 2.1. However, we know that
$(2\cos{\frac{1}{2n+2}}\pi)^{\frac{2}{k}}$ is not an eigenvalue of
$(C_{4}+P_{n-3})^{(k)}$ by Lemma 2.3. Then $\widetilde{D}_{n}$ is DHS.

Since $D_{3}$ is an induced subgraph of $\widetilde{E}_{6}$ and
$2\cos{\frac{\pi}{8}}$ is an eigenvalue of $D_{3}$, it yields that
$(2\cos{\frac{\pi}{8}})^{\frac{2}{k}}$ is an eigenvalue of
$\widetilde{E}^{(k)}_{6}$ by Lemma 2.1. However, we know that
$(2\cos{\frac{\pi}{8}})^{\frac{2}{k}}$ is not an eigenvalue of
$(C_{6}+K_{1})^{(k)}$ by Lemma 2.3. Then $\widetilde{E}_{6}$ is DHS. ∎

##  4 High-ordered cospectral invariants of trees

Studying the high-ordered cospectral invariants of graph $G$ is helpful to
study the high-ordered spectral characterizations of $G$. In this section, we
give some high-ordered cospectral invariants of trees by using the spectral
moments of power hypertrees. As their applications, we give infinitely many
pairs of trees with same spectrum but different high-ordered spectrum.

Let $T$ and $T^{*}$ be two high-ordered cospectral trees. Then
$\mathrm{S}_{d}(T^{(k)})=\mathrm{S}_{d}({T^{*}}^{(k)})$ for all positive
integer $k$ and $d$. Let $\mathfrak{T}({m})$ denote the set of trees with $m$
edges. From the spectral moment formula shown in Lemma 2.5, we obtain

$\displaystyle\sum\limits_{m=1}^{\frac{d}{k}}\frac{1}{2}{{(k-1)}^{\left({|E(T)|-m}\right)(k-1)}}{k^{m(k-2)+1}}{\sum\limits_{\widehat{T}\in{\mathfrak{{T}}({m})}}{c_{\frac{2d}{k}}(\widehat{T})}}\left(N_{T}(\widehat{T})-N_{T^{*}}(\widehat{T})\right)=0$
(4.1)

for all positive integer $k\geq 2$. From Equation (4.1), we obtain some high-
ordered cospectral invariants of trees.

######  Theorem 4.9.

If a tree $T$ is high-ordered cospectral with tree $T^{*}$, then

${\sum\limits_{\widehat{T}\in{\mathfrak{T}(m)}}{c_{d}(\widehat{T})}}N_{T}(\widehat{T})={\sum\limits_{\widehat{T}\in{\mathfrak{T}({m})}}{c_{d}(\widehat{T})}}N_{T^{*}}(\widehat{T})$

for all positive even integer $d$ and all $m\in[\frac{d}{2}]$.

###### Proof.

Since $T$ and $T^{*}$ are high-ordered cospectral, we have
$\mathrm{S}_{d}(T^{(k)})=\mathrm{S}_{d}({T^{*}}^{(k)})$ for all positive
integer $k\geq 2$. Since $\mathrm{S}_{d}(T^{(k)})=0$ if $k\nmid d$, we assume
$k\mid d$ in the following proof. Let $d=kz$. From Equation (4.1), we have

$\displaystyle\sum\limits_{m=1}^{z}\frac{1}{2}{{(k-1)}^{\left({|E(T)|-m}\right)(k-1)}}{k^{m(k-2)+1}}{\sum\limits_{\widehat{T}\in{\mathfrak{T}({m})}}{c_{2z}(\widehat{T})}}\left(N_{T}(\widehat{T})-N_{T^{*}}(\widehat{T})\right)=0$
(4.2)

for all positive integer $k$. Let
$f_{m}(k)=\frac{1}{2}{{(k-1)}^{\left({|E(T)|-m}\right)(k-1)}}{k^{m(k-2)+1}}$.
Let
$y_{m}={\sum\nolimits_{\widehat{T}\in{\mathfrak{T}({m})}}{c_{2z}(\widehat{T})}}\left(N_{T}(\widehat{T})-N_{T^{*}}(\widehat{T})\right)$
for all $m\in[z]$. From Equation (4.2), we have
$\sum\nolimits_{m=1}^{z}{f_{m}(k)y_{m}}=0$. Since
$\frac{{{f_{i}}(k)}}{{{f_{i-1}}(k)}}=\frac{{{k^{k-2}}}}{{{{(k-1)}^{k-1}}}}$
for $i=2,3,\cdots,z$, we get
${{f_{m}}(k)}={{f_{1}}(k)}(\frac{{{k^{k-2}}}}{{{{(k-1)}^{k-1}}}})^{m-1}$. Then
$\sum\nolimits_{m=1}^{z}{{{f_{1}}(k)}(\frac{{{k^{k-2}}}}{{{{(k-1)}^{k-1}}}})^{m-1}y_{m}}=0$.
For $z\geq 1$, take $k$ as $k_{1},k_{2},\ldots,k_{z}$, which are $z$ different
positive integers. It is following that

$\displaystyle\left[{\begin{array}[]{*{20}{c}}{{f_{1}}({k_{1}})}&{{f_{1}}(k_{1})}\frac{{{k_{1}^{k_{1}-2}}}}{{{{(k_{1}-1)}^{k_{1}-1}}}}&\cdots&{{f_{1}}(k_{1})}(\frac{{{k_{1}^{k_{1}-2}}}}{{{{(k_{1}-1)}^{k_{1}-1}}}})^{z-1}\\\
{{f_{1}}({k_{2}})}&{{f_{1}}(k_{2})}\frac{{{k_{2}^{k_{2}-2}}}}{{{{(k_{2}-1)}^{k_{2}-1}}}}&\cdots&{{f_{1}}(k_{2})}(\frac{{{k_{2}^{k_{2}-2}}}}{{{{(k_{2}-1)}^{k_{2}-1}}}})^{z-1}\\\
\vdots&\vdots&\vdots&\vdots\\\
{{f_{1}}({k_{z}})}&{{f_{1}}(k_{z})}\frac{{{k_{z}^{k_{z}-2}}}}{{{{(k_{z}-1)}^{k_{z}-1}}}}&\cdots&{{f_{1}}(k_{z})}(\frac{{{k_{z}^{k_{z}-2}}}}{{{{(k_{z}-1)}^{k_{z}-1}}}})^{z-1}\\\
\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{y_{1}}}\\\ {{y_{2}}}\\\
\vdots\\\ {{y_{z}}}\\\ \end{array}}\right]=0.$ (4.11)

Then the coefficient matrix of the Equation (4.11) is a Vandermonde matrix.
And since $k_{1},k_{2},\ldots,k_{z}$ are distinct, the determinant of the
coefficient matrix are not equal to zero. Then $y_{m}=0$ for all $m\in[z]$. ∎

Let $m=|E(T)|=|E(T^{*})|$ in Theorem 4.9, we directly get the following high-
ordered cospectral invariants.

######  Theorem 4.10.

Let $T$ and $T^{*}$ be two high-ordered cospectral trees. Then
$c_{d}(T)=c_{d}(T^{*})$ for any positive even number $d\geq{2|E(T)|}$.

Let the set of trees with $m$ edges
${{\mathfrak{T}}({m})}=\\{{\widehat{T}}_{1},{\widehat{T}}_{2},\ldots,{\widehat{T}}_{|{{\mathfrak{T}}({m})}|}\\}$.
From Theorem 4.9, we get the following equation after taking $d$ as
$d_{1},d_{2},\ldots,d_{|{\mathfrak{T}(m)}|}$.

$\displaystyle\left[{\begin{array}[]{*{20}{c}}{{c_{{d_{1}}}}({\widehat{T}_{1}})}&{{c_{{d_{1}}}}({\widehat{T}_{2}})}&\cdots&{{c_{{d_{1}}}}({\widehat{T}_{|{\mathfrak{T}(m)}|}})}\\\
{{c_{{d_{2}}}}({\widehat{T}_{1}})}&{{c_{{d_{2}}}}({\widehat{T}_{2}})}&\cdots&{{c_{{d_{2}}}}({\widehat{T}_{|{\mathfrak{T}(m)}|}})}\\\
\vdots&\vdots&\vdots&\vdots\\\
{{c_{{d_{|{\mathfrak{T}(m)}|}}}}({\widehat{T}_{1}})}&{{c_{{d_{|{\mathfrak{T}(m)}|}}}}({\widehat{T}_{2}})}&\cdots&{{c_{{d_{|{\mathfrak{T}(m)}|}}}}({\widehat{T}_{|{\mathfrak{T}(m)}|}})}\\\
\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{h_{1}}}\\\ {{h_{2}}}\\\
\vdots\\\ {{h_{|{\mathfrak{T}(m)}|}}}\\\ \end{array}}\right]=0,$ (4.20)

where $h_{i}=N_{T}(\widehat{T}_{i})-N_{T^{*}}(\widehat{T}_{i})$ for all
$i\in[|{\mathfrak{T}(m)}|]$. If there exist
$d_{1},d_{2},\ldots,d_{|{\mathfrak{T}(m)}|}$ such that the coefficient matrix
of Equation (4.20) is nonsingular, we get
$N_{T}(\widehat{T})=N_{T^{*}}(\widehat{T})$ for all
$\widehat{T}\in\mathfrak{T}(m)$. By the formula for the spectral moment
coefficients for trees, i.e. Equation (2.3), we can calculate the spectral
moment coefficients $c_{d}({\widehat{T}})$. We calculate the $d$-th order
spectral moment coefficients of trees with 3 edges for $d=6,8$, the $t$-th
order spectral moment coefficients of trees with 4 edges for $t=8,10,12$ and
the $l$-th order spectral moment coefficients of trees with 5 edges for
$l=10,12,14,16,18,20$ (see Table 1 and Table 2).

(a) $P_{2}$

(b) $P_{3}$

(c) $P_{4}$

(d) $S_{4}$

(e) $P_{5}$

(f) $Q_{5}$

(g) $S_{5}$

(h) $P_{6}$

(i) $Q_{6}$

(j) $R_{6}$

(k) $H_{6}$

(l) $J_{6}$

(m) $S_{6}$

Figure 3: Trees with edges at most $5$

| $c_{d}(P_{4})$ | $c_{d}(S_{4})$
---|---|---
$d=6$ | 6 | 12
$d=8$ | 32 | 72
(a) The spectral moment coefficients of trees with 3 edges

| $c_{d}(P_{5})$ | $c_{d}(Q_{5})$ | $c_{d}(S_{5})$
---|---|---|---
$d=8$ | 8 | 16 | 48
$d=10$ | 60 | 140 | 480
$d=12$ | 300 | 804 | 3120
(b) The spectral moment coefficients of trees with 4 edges

Table 1: The spectral moment coefficients of trees with 3 or 4 edges | $c_{d}(P_{6})$ | $c_{d}(Q_{6})$ | $c_{d}(R_{6})$ | $c_{d}(H_{6})$ | $c_{d}(J_{6})$ | $c_{d}(S_{6})$
---|---|---|---|---|---|---
$d=10$ | 10 | 20 | 20 | 40 | 60 | 240
$d=12$ | 96 | 216 | 228 | 504 | 792 | 3600
$d=14$ | 588 | 1484 | 1652 | 3976 | 6552 | 33600
$d=16$ | 2944 | 8304 | 9728 | 25216 | 43680 | 252000
$d=18$ | 13158 | 41328 | 50832 | 140832 | 257184 | 1668240
$d=20$ | 54730 | 190800 | 245880 | 724320 | 1398600 | 10206000
Table 2: The spectral moment coefficients of trees with 5 edges

We obtain the following higher-ordered cospectral invariants about the number
of subtrees.

######  Theorem 4.11.

If a tree $T$ is high-ordered cospectral with a tree $T^{*}$, then
$N_{T}(\widehat{T})=N_{T^{*}}(\widehat{T})$ for any tree $\widehat{T}$ within
5 edges.

###### Proof.

Since $T$ and $T^{*}$ are high-ordered cospectral, we know that $T$ and
$T^{*}$ are cospectral. Then $N_{T}(P_{2})=N_{T^{*}}(P_{2})$ and
$N_{T}(P_{3})=N_{T^{*}}(P_{3})$. From Table 1 and Table 2, we get the
coefficient matrix of Equation (4.20) when $m=3,4,5$. It is easy to check that
these coefficient matrices are nonsingular. Then
$N_{T}(\widehat{T})=N_{T^{*}}(\widehat{T})$ for all
$\widehat{T}\in\mathfrak{T}(m)$, $m=1,2,3,4,5$. ∎

By the above high-ordered cospectral invariants, we get infinitely many pairs
of cospectral trees with different high-ordered spectrum. In [19], Schwenk
gave a useful method to construct cospectral trees and proved his famous
conclusion:“Almost all trees are not DS”. Let $T_{u}$ and $T_{v}$ be the tree
$T$ rooted at vertices $u$ and $v$, respectively. For any rooted tree $F$, the
coalescences $F\cdot T_{u}$ and $F\cdot T_{v}$, as shown as in Figure 4, are
cospectral trees [22]. We have $N_{F\cdot T_{v}}(R_{6})-N_{F\cdot
T_{u}}(R_{6})=d$, where $d$ is the degree of the root of $F$. By Theorem 4.11,
we know that $F\cdot T_{u}$ is not high-ordered cospectral with $F\cdot
T_{v}$.

######  Proposition 4.12.

Cospectral trees $F\cdot T_{u}$ and $F\cdot T_{v}$ have different high-ordered
spectrum.

(a) $F\cdot T_{u}$

(b) $F\cdot T_{v}$

Figure 4: The coalescences $F\cdot T_{u}$ and $F\cdot T_{v}$

##  5 Conclusions and future work

In this paper, we show that the high order spectrum of graphs can determine
more graphs than the spectrum of graphs. Up to now, we have not found any non-
isomorphic high order cospectral graphs, that is to say, we have not found any
graphs which are not DHS. We are very interested in the existence and
construction of non-isomorphic high order cospectral graphs.

The Reconstruction Conjecture of graphs is a famous conjecture in graph theory
[13, 21, 12]. In [14], Kelly proved that “The Reconstruction Conjecture” is
true for trees. This result is often called “The Reconstruction Theorem for
Trees”. By the Theorem 4.11 in this paper, we know that the number of subtrees
are high order cospectral invariants of trees. From Kelly’s Lemma [14], we
know that the number of subgraphs are reconstructible parameter of graphs. We
believe that the high order spectral characterizations of trees are closely
related to the “Reconstruction Theorem for Trees”. Based on the above
introduction, we have the following conjecture.

######  Conjecture 5.13.

Two trees are isomorphic if and only if they have same high-ordered spectrum.

## References

## References

  * [1] Romain Boulet and Bertrand Jouve.  The lollipop graph is determined by its spectrum.  Electronic Journal of Combinatorics, 15(1):1948–1953, 2008.
  * [2] Kauê Cardoso, Carlos Hoppen, and Vilmar Trevisan.  The spectrum of a class of uniform hypergraphs.  Linear Algebra and its Applications, 590:243–257, 2020.
  * [3] Lixiang Chen and Changjiang Bu.  A reduction formula for the characteristic polynomial of hypergraph with pendant edges.  Linear Algebra and its Applications, 611:171–186, 2021.
  * [4] Lixiang Chen, Changjiang Bu, and Jiang Zhou.  Spectral moments of hypertrees and their applications.  Linear and Multilinear Algebra, 2021.  Doi:10.1080/03081087.2021.1953431.
  * [5] Gregory J. Clark and Joshua N. Cooper.  A Harary-Sachs theorem for hypergraphs.  Journal of Combinatorial Theory, Series B, 149:1–15, 2021.
  * [6] Joshua Cooper and Aaron Dutle.  Spectra of uniform hypergraphs.  Linear Algebra and its Applications, 436(9):3268–3292, 2012.
  * [7] Dragoš Cvetković and Peter Rowlinson.  Spectra of unicyclic graphs.  Graphs and Combinatorics, 3(1):7–23, 1987.
  * [8] Dragoš M. Cvetković.  Graphs and their spectra.  Publikacije Elektrotehničkog fakulteta. Serija Matematika i fizika, (354-356):1–50, 1971.
  * [9] Dragoš M. Cvetković, Horst Sachs, and Michael Doob.  Spectra of Graphs: Theory and Applications.  VEB Deutscher Verlag der Wissenschaften, 1980.
  * [10] Chris D. Godsil and Brendan D. McKay.  Constructing cospectral graphs.  Aequationes Mathematicae, 25(1):257–268, 1982.
  * [11] Hs H. Günthard and Hans Primas.  Zusammenhang von graphentheorie und mo-theorie von molekeln mit systemen konjugierter bindungen.  Helvetica Chimica Acta, 39(6):1645–1653, 1956.
  * [12] Frank Harary.  On the reconstruction of a graph from a collection of subgraphs.  In Theory of Graphs and its Applications (Proc. Sympos. Smolenice, 1963), pages 47–52. Publ. House Czechoslovak Acad. Sci. Prague, 1964\.
  * [13] Paul J. Kelly.  On isometric transformations.  PhD thesis, University of Wisconsin, 1942.
  * [14] Paul J. Kelly.  A congruence theorem for trees.  Pacific Journal of Mathematics, 7(1):961–968, 1957.
  * [15] Mirko Lepović and Ivan Gutman.  No starlike trees are cospectral.  Discrete Mathematics, 242(1-3):291–295, 2002.
  * [16] Lek-Heng Lim.  Singular values and eigenvalues of tensors: a variational approach.  In 1st IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, pages 129–132. IEEE, 2005.
  * [17] Russell Merris.  Laplacian matrices of graphs: a survey.  Linear Algebra and its Applications, 197:143–176, 1994.
  * [18] Liqun Qi.  Eigenvalues of a real supersymmetric tensor.  Journal of Symbolic Computation, 40(6):1302–1324, 2005.
  * [19] Allen J. Schwenk.  Almost all trees are cospectral.  In Frank Harry, editor, New Directions in the Theory of Graphs, pages 275–307, New York, 1973. Academic Press.
  * [20] John H. Smith.  Some properties of the spectrum of a graph.  Combinatorial Structures and their Applications, pages 403–406, 1970.
  * [21] Stanislaw M. Ulam.  A collection of mathematical problems, volume 8 of Interscience Tracts in Pure and Applied Mathematics.  Interscience, New York-London, 1960.
  * [22] Edwin R. van Dam and Willem H. Haemers.  Which graphs are determined by their spectrum?  Linear Algebra and its Applications, 373:241–272, 2003.
  * [23] Edwin R. van Dam and Willem H. Haemers.  Developments on spectral characterizations of graphs.  Discrete Mathematics, 309(3):576–586, 2009.
  * [24] Lothar von Collatz and Ulrich Sinogowitz.  Spektren endlicher grafen.  Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg, 21(1):63–77, 1957.
  * [25] Wang Wei.  A simple arithmetic criterion for graphs being determined by their generalized spectra.  Journal of Combinatorial Theory, Series B, 122:438–451, 2017.
  * [26] Jiang Zhou, Sun Lizhu, Wenzhe Wang, and Changjiang Bu.  Some spectral properties of uniform hypergraph.  Electronic Journal of Combinatorics, 21:4–24, 2014.