text:
\section*{Introduction}

Formation of robust gapless modes in the presence of defects is a defining property of the topological matter \cite{halperin1982quantized,teo2010topological}. Moreover, this property shows a great promise with regards to applications of topological materials, such as dissipationless lines \cite{buttiker1988absence}, spintronics, and quantum computation\cite{qi2011topological}. Therefore, the ability to control the presence and properties of defects is of both fundamental and practical interest.

 \begin{figure}[h!]
 \centering
	\includegraphics[width=0.8\linewidth]{"lcm_comp2.pdf"}
	\caption{\textbf{The dynamics of the LCM under quench isolating a single site in the sample's middle.} \textbf{a} The scheme of the "cutt-off" scenario. At the t=0 the hoppings from the central site are set to zero. Simultaneously the on-site parameters are changed.  \textbf{b} Teleportation of the marker value from the central site which was cut off. The evolution of the LCM distribution is described by Eq. ~\ref{eq:chern_m} as in Refs.\cite{d2015dynamical,privitera2016quantum, caio2019topological, toniolo2018time} . The insert shows the behaviour of $C(\mathbf{r},t)$ non-local currents close to the site. Shown here is one of the many possible non-local currents described by Eq. ~\ref{eq:chern_m}, for more details see Supplementary materials. \textbf{c} Distribution of $C(\mathbf{r},t)$ Eq. ~\ref{eq:nq_streda} at different times. The blue arrows on insert shows the absence of LCM currents from the isolated site.}
	\label{fig:comp_lcm}
\end{figure}
Samples with defects are inhomogeneous and require special care from a theoretical perspective. Global topological indices, e.g. Chern number \cite{thouless1982quantized}, are not applicable directly to the systems with topological defects, as they characterise the whole system. Recently, a family of quasitopological quantities, coined local topological markers, has been developed and studied \cite{kitaev2006anyons,bianco2011mapping, loring2011disordered}. They provide an estimation of the global topological index based on the parameters of the neighborhood of a given point in space and are not necessarily strictly quantized. The robust characteristic is the average of the markers over large areas of a system \cite{bianco2014chern}.

 Out of equilibrium systems reveal another striking difference between the global and local topological quantities. Global indices of infinite systems are conserved during unitary evolution unless the symmetry of the problem changes dynamically \cite{toniolo2018time,caio2015quantum,mcginley2019classification}. In finite samples, characterised by local markers, the situation is different. The topological index - an average of local markers over large areas - may change as was demonstrated by numerical simulations \cite{d2015dynamical,privitera2016quantum, caio2019topological} and analytically \cite{toniolo2018time}.

Neither local topological markers\cite{kitaev2006anyons,bianco2011mapping, loring2011disordered}, nor their equations of motions are unique. It would be highly desirable that the dynamics for the markers would be local. That is the marker would flow from a site to its neighborhood, not teleport to distant sites. In an insightful paper \cite{caio2019topological} it was conjectured that the dynamics of the marker satisfy this property and governed by local currents emanating from the system's boundaries, defined implicitly through the lattice continuity equation on the marker. We demonstrate that such currents could not be defined without the modification of equations of motions for the markers used before \cite{d2015dynamical,privitera2016quantum, caio2019topological, toniolo2018time} . Most clearly it can be seen in the numerical experiment shown in Fig.\ref{fig:comp_lcm} \textbf{b}. There the teleportation of the marker from the disconnected central site is shown.

In the present manuscript we suggest a new physical perspective on the dynamics of the local topological markers. We have used the connection of markers to the charge response to a magnetic field. New equations of motion have a clear physical meaning, do obey the continuity equation, and are much easier to measure in an experiment as require only density measurements. Furthermore, our approach does not rely on single-particle projectors onto filled states. Therefore it might be used in the interacting systems.

Equipped with the new approach, we numerically study the possibility to control the position of line defects in a finite 2-d Chern insulator. For a start, we studied a predictable setting depicted schematically in Fig.\ref{fig:setting}. The system is prepared so that the central region is in a topologically nontrivial phase and surrounded by the trivial insulator, as indicated by the initial distribution of the local topological marker.  We considered the dynamics of local markers and their currents, corresponding to the slow movement to the right of the topological domain. We demonstrate the possibility of quasiadiabatic control of topological marker distribution in the system. We show that the steady-state distribution of the local marker catches the location of gapless modes.

\section*{Local Marker Dynamics}

\textit{Local Chern marker.} Two dimensional non-interacting systems are classified topologically by Chern number \cite{thouless1982quantized} $C$ in the absence of symmetries \cite{altland1997nonstandard} to protect. Strictly speaking, topological phases with $C\neq 0$ realize in the thermodynamic limit on a torus. Real-world samples subjected to open boundary conditions necessarily contain topological defects; one can be certain that at least the boundary vacuum - Chern insulator is present in the system. Further, the sample may contain macroscopically large patches with parameters corresponding to different Chern numbers. One needs a substantially local quantity that allows to label different parts of the systems according to their Chern number to work in such settings. Local topological markers \cite{kitaev2006anyons,bianco2011mapping, loring2011disordered} do precisely the job.

Physically the markers may be connected to the charge response to applied magnetic \cite{bianco2013orbital,loring2011disordered,mitchell2018amorphous, kitaev2006anyons} or electric \cite{peru} field.  For concreteness, we will use the Local Chern Marker\cite{bianco2011mapping} (LCM). \footnote{The Bott index's\cite{loring2011disordered}  application to the systems with zero-modes is problematic. Much of the discussion to follow can be applied as well to the Kitaev's marker closely related to the LCM both formally \cite{bianco2011mapping} and physically \cite{mitchell2018amorphous}}. For this marker the local version of the Streda formula \cite{streda1982theory} holds\cite{bianco2013orbital} up to insignificant corrections appearing on the system borders or defects:

 \begin{equation}  \label{eq:chern_streda}
 \mathrm{C_S}(\mathbf{r})= \phi_0\frac{\delta n(\mathbf{r})}{\delta B_z},
\end{equation}

here $\phi_0=$ is the flux quanta. The variation of the average density $n(\mathbf{r})$ is taken with respect to a uniform magnetic field which is supposed to be turned on adiabatically slow. The Kitaev marker \cite{kitaev2006anyons} can be obtained from a very similar equation with the uniform magnetic field replaced by a local one \cite{mitchell2018amorphous}.

 A generic example of the marker distribution in an inhomogeneous sample is presented in Fig.~\ref{fig:setting}\textbf{b}. The marker is very close to quantized value in the ''bulk'' of a patch in a trivial or topological phase. In vicinity of defects LCM develops sharp negative peaks. A physical intuition for this behaviour is that defects serve as electron sources for the currents flowing towards the bulk when a small magnetic field is applied.

The definition of the markers is typically \cite{bianco2013orbital,loring2011disordered, kitaev2006anyons} given in terms of the  projectors $\hat{P}$ onto the single-particle filled states. The LCM in particular is defined as the spatially resolved commutator of  $\hat{X}$ and $\hat{Y}$ operators dressed in projectors:

\begin{equation}  \label{eq:chern_m}
 C(\mathbf{r})=  - 2\pi i  \sum_{\alpha}\left\langle\mathbf{r_\alpha}|\left[\hat{P} \hat{X} \hat{P}, \hat{P} \hat{Y} \hat{P}\right]| \mathbf{r_\alpha}\right\rangle.
 \end{equation}

This the formula highlights two more of the marker's properties. First, it is a quasi-local quantity due to exponential decay of $\hat{P}(\bm{r},\bm{r}')$ with the distance between $\bm{r}$ and $\bm{r}'$ \cite{lieb1972finite}. Second, it sums over a finite sample to zero, as a commutator is traceless.

\textit{LCM dynamics}

It is very tempting to use the markers to understand the topological properties out of equilibrium. How the time-dependent markers should be defined is not immediately clear, however. The markers are not strictly speaking operators, obeying the Heisenberg equation of motion. Rather they are operator-valued functions of the system's state. An appealing approach is to use the same function (e.g. Eq. \ref{eq:chern_m} for LCM) of the state and allow it to evolve \cite{privitera2016quantum,toniolo2018time,caio2019topological,d2015dynamical}. It was revealed that in finite-size systems the markers can change \cite{d2015dynamical,privitera2016quantum}, in contrast to global topological indices \cite{d2015dynamical,caio2015quantum,toniolo2018time}. Importantly, their evolution reflects the change in a topological phase \cite{d2015dynamical,privitera2016quantum, caio2019topological}.

A remarkably simple and appealing picture of the dynamics of the LCM was presented in Ref. \cite{caio2019topological}. It was conjectured that the marker's evolution is governed by local currents emanating from the edges of a sample. The currents were defined through a lattice continuity equation. The claim is highly non-trivial, as the LCM is not a fully local quantity. Let us stress that, the global conservation of the LCM - $\sum C(\bm{r})=0$ - does not guarantees the existence of a local continuity equation. In fact, such a continuity equation can not be satisfied by $C(\mathbf{r},t)$ in general.

Most clearly it can be observed in the following numerical experiment, shown in Fig.~\ref{fig:comp_lcm}\textbf{(a)}. Consider a homogeneous sample in the topological phase. At the time $t=0$ a central site is cut off from the rest of the system. That is all the hoppings from it are quenched to zero. Simultaneously, the on-site parameters are changed. Were $C(\mathbf{r},t)$ satisfy a lattice continuity equation, no change of the marker would be observed on the site.

The Fig. \ref{fig:comp_lcm}(c) demonstrates that the on-site value of the marker does change. The time-dependent Eq.~\ref{eq:chern_m} predicts the teleportation of the marker from the disconnected site. Therefore, the numerical evidences \cite{caio2019topological} of the current-like local behaviour of the LCM are seeking an explanation and elaboration. We demonstrate in the Supplementary how the non-local currents in this case may be explicitly defined.

We suggest a new perspective on the markers' evolution, generalizing the Eq.~\ref{eq:chern_streda}:

\begin{equation}  \label{eq:nq_streda}
\mathrm{C}_S(\mathbf{r}, t)= \phi_0\frac{\delta n(\mathbf{r}, t)}{\delta B_z},
\end{equation}
here the magnetic field is supposed to be turned on adiabatically at $t=-\infty$.

Serving as good topological marker as the defined by means of Eq.~\ref{eq:chern_m}, the new approach has important advantages. First, this formula has a clear physical interpretation as a spatial resolution of the Hall conductivity. Second, it is easier to measure.  The typical protocol of reconstructing the LCM requires the knowledge of off-diagonal elements of the projectors $P(\bm{r},\bm{r}')$ \cite{caio2019topological,irsigler2019microscopic}. This is a much harder experimental task \cite{ardila2018measuring}, then the density measurements \cite{bakr2009quantum}. Third, the evolution of the LCM as described by the Eq.~\ref{eq:nq_streda} satisfies the lattice continuity equation:

\begin{equation} \label{eq:cont_eq}
\begin{split}
 \partial_{t} \mathrm{C}_S(\mathbf{r}, t)=-\text{div} \mathbf{J}^c ,\\
J^c(\mathbf{r},\bm{r}', t) = \frac{\delta J^e(\mathbf{r},\bm{r}', t)}{\delta B_z}.
\end{split}
\end{equation}

where $J^e(\mathbf{r},\bm{r}', t)$ is the average current on the bond between sites $\bm{r}$ and $\bm{r}'$.

Therefore the Eq.~\ref{eq:nq_streda} forces the dynamics to be local, as we can see both numerically (Fig.~\ref{fig:comp_lcm}\textbf{(c)}) and analytically. Moreover, it provides a physical explanation why the current-like behaviour of the LCM was observed\cite{caio2019topological} in the first place. The reason is that the marker currents can be directly connected to the electric currents as Eq.~\ref{eq:cont_eq} demonstrates. See Supplementary for the further discussion of the quench dynamics as in Ref.~\cite{caio2019topological}.

\section*{Controlling an inhomogeneous topological lattice.}
	\begin{figure}[h!]
		\centering
		\includegraphics[width=0.9\linewidth]{"main2.pdf"}
		\caption{\textbf{The movement of the "topological" area inside the Chern insulator in the trivial phase.} \textbf{a} The evolution of LCM corresponding to one-site shift of the area in the slice indicated in \textbf{c}. \textbf{b} Corresponding evolution of $y$ component of electric current in units of $et_ha/\hbar$. The current was obtained after shifting of chemical potential by $0.1\Delta$. \textbf{d}-\textbf{g} The spatial distribution of the LCM and its currents on the lattice at different times. $\tau=100$ $\hbar/t_h$.  Lattice size is $22$ by $21$, topological domain size is $10$ by $10$.
		}
		\label{fig:flow}
	\end{figure}

The control over the spatial distribution of topological properties is a very important problem both fundamentally and practically. The tunable position of robust zero modes is particularly fruitful, as these modes are key to the most impressive possible applications of topological matter \cite{buttiker1988absence, qi2006topological}. From the local marker perspective the task might be stated as the possibility to dynamically control the marker distribution. As a starting point on the road we studied the slow dynamics of the topological markers.

We considered the slow movement of a domain in topological phase inside an otherwise trivial sample as presented in Fig.~\ref{fig:setting}. After a time $\tau$, the domain with ''topological''  parameters shifts by one site to the right under a linear time-protocol (See Methods for the details).  Does a shift in the parameter's distribution mean a real shift of the topological domain? To address the question we calculated the electric currents shifting the position of the chemical potential of the system by $\sim 0.1\Delta$. The steady-state amplitude of the electric currents in the y-direction in the middle of the sample and the distribution of currents on the bonds of the lattice is presented in Fig.~\ref{fig:flow} \textbf{b}-\textbf{c}. The shift of edge currents confirms that topological domain moved one site to the right.

As the system evolves, the distribution of the marker changes. The distribution of the local Chern marker follows the shift of the topological area, provided that the process of moving was sufficiently slow. The protocol period should be much larger than the inverse of bulk gap $\tau \gg 1/\Delta$. The Fig.~\ref{fig:flow} \textbf{a} demonstrates how the distribution evolves during the time $\tau=100$. In this regime, the steady-state distribution resembles that of the equilibrium ground state transferred as a whole one site to the right.

The local currents of the marker were observed near the borders of the topological domain. These currents are shown by the arrows in Fig.~\ref{fig:flow} \textbf{d}-\textbf{g}.  The defects played the role of the electron reservoirs for the bulk in the presence of the small magnetic field.

	\begin{figure}
		\centering
		\includegraphics[width=1\linewidth]{"oscill2.2.pdf"}
		\caption{\textbf{Topological marker behavior on the site in three different modes.} \textbf{a} Time-dependence of LCM on the site $x = x_r+1, \: y=y_m$ highlighted in blue in the insert. The dark blue curve corresponds to $C_S(r,t)$ from the Streda-based formula~\ref{eq:nq_streda}, the light blue $C(r,t)$ to the formula~\ref{eq:chern_m}. The topological domain here moves three sites to the right, during $\Delta t=300$. In the beginning, the site is to the right of the topological domain boundary. At $t=100$ it appears on the boundary, then $t=200$ in an intermediate mode, finally, at $t=300$ in the topological bulk. In the insert the initial position if highlighted. \textbf{b} Fouirer components of $C_S(r,t)$ oscillations in three modes defined on the previous figure. $\Delta$ is the energy gap. Lattice size is $22$ by $21$, topological domain size $10$ by $10$.}
		\label{fig:ch_afr}
	\end{figure}

\textit{Oscillations of LCM.}

We compared the time-dependency of $C_S(t)$ and $C(t)$ (Eq.~\ref{eq:chern_m},\ref{eq:nq_streda}). They agree qualitatively, though differ in details. As it can be seen from Fig.\ref{fig:ch_afr} \textbf{a} the amplitude of oscillations of the first is several times higher than the amplitude of the second. The nature of these oscillations can be explained by the fact that inhomogeneous samples always contain zero modes. The small energy of hamiltonian perturbation excites these levels and makes densities on the border oscillate. The rest of the system also responds to these fluctuations.

We considered the process of moving a topological domain three sites to the right. In Fig.\ref{fig:ch_afr} \textbf{a} the time-dependence of LCM on the site with coordinates $x = x_r+1, \: y=y_m$ is presented. When the site is on the topological edge at $t=100$, the LCM have maximum in both the value and the amplitude of oscillations. From $100$ to $200$ there is a drift of the marker towards its final value $C=1$ with a constant speed determined by $1/\tau$. The site continues to move into topological bulk from $200$ to $300$ and undergoes free-evolution without the change in Hamiltonian for the time remaining.

The site position determines the frequency of the oscillations and their magnitude. Amplitude-frequency response is presented in the Fig.\ref{fig:ch_afr} \textbf{b} for oscillations in different regimes. In the case when the site is moving from the edge at the time from $100$ to $200$, the dominant frequency coincides with the transition from the zero-energy state to the bottom of the band( Van Hove singularities In the thermodinamic limit). When the site moves into the bulk, only transitions between two energy bands survive. This is shown in the Fig.\ref{fig:ch_afr} \textbf{b}. The oscillations are less pronounced, the further from the border between different phases the site is - another indication of the locality of the dynamics of the marker.

\section*{Discussion}

 We provided a new clear physical picture of the local marker's evolution. We demonstrated that the evolution of the markers can be enforced to be fully local, despite the markers complex quasi-local definition. Local marker currents were directly connected to the local electrical currents. With the new tool at hand, we studied the possibility to control the line defects in a non-homogeneous finite sample. Our work paves the way to the development of tunable devices harnessing topologically protected modes on line defects.

The setting we studied may be realized in an experiment. Most notably, modern cold atomic platforms have the technologies to create and control topological interfaces \cite{goldman2016creating} in a single sample. In our approach, the dynamics of the marker can be traced by density measurements only. Such a technique is available at the moment \cite{bakr2009quantum}. The usual approach requires the knowledge of the single-particle density matrices \cite{irsigler2019microscopic} which is much harder in practice \cite{ardila2018measuring}.

Let us suggest possible extensions of our work. Interacting topological insulators and in particular fractional Chern insulators are very interesting to address with our method. Our equations of motion are ready can be applied to many-body systems as they do not rely on single-particle projectors as usual. In an interacting system the projectors to filled states are not defined, complicating the generalization of the local markers. On the other hand, recent equilibrium calculations indicate that the Streda-based formula Eq.\ref{eq:chern_streda} can be used as a local marker for fractional phases\cite{wang2021measurable}. Another important task is to find the optimal method to control the local marker's distribution. This requires further analytical and numerical studies of the nonhomogeneous topological systems out of equilibrium.

\section*{Methods}
\label{sec:methods}

\textit{Model.} We consider a model of Chern insulator ~\cite{qi2006topological} with Hamiltonian:

\begin{equation} \label{eq:Ch_eq}
H(t)=\sum_{n} t_h (c_{n,\sigma}^{\dagger} \frac{\hat{\sigma}_{z}-i \hat{\sigma}_{x}}{2}  c_{n+\hat{x},\sigma} + c_{n,\sigma}^{\dagger} \frac{\hat{\sigma}_{z}-i \hat{\sigma}_{y}}{2}  c_{n+\hat{y},,\sigma} + h. c.) + m(\mathbf{r}, t) \sum_{n} c_{n}^{\dagger} \hat{\sigma}_{z}  c_{n},
\end{equation}
for fermi-particles on a finite $N_x\times N_y$ square lattice with the lattice constant $a=1$. Here, Pauli matrices $\sigma_i$ act in the spin subspace. In the homogeneous case $m(\mathbf{r},t) = const$, the model's phase is determined by the ratio of parameters $m$ and $t_h$:

\begin{equation} \label{m_det}
\begin{array}{l}
-2t_h<m<0 \quad \text{topological;} \quad C=1 \\
0<m<2t_h \quad \text{topological;} \quad C=-1 \\
m<-2t_h, \; m>2t_h \quad \text{trivial;} \quad C=0\\
\end{array}
\end{equation}

The phases are separated by the Dirac semi-metal phases. The energy and time scale are fixed by setting $t_h$ and $\hbar$ to $1$.

 To calculate the LCM according to the Eq.\ref{eq:nq_streda} the density of electrons calculated numerically in the presence The magnetic field is coupled to the system by means of the Peierl's substitution, which introduces a site-dependent phase factor in the hopping matrix $t(\bm{r}_1,\bm{r}_2)$:

	\begin{equation}
	t(\bm{r}_1,\bm{r}_2) = t_h\cdot e^{i\frac{2\pi}{\phi_0}\int_{\bm{r}_1}^{\bm{r}_2}\bm{A}(\bm{r})\cdot \bm{dr}},
	\label{eq:hop}
	\end{equation}

The variation in Eq.~\ref{eq:nq_streda} is taken with respect to a uniform magnetic field with a flux small compared to the flux quantum.

\textit{Time protocol.}

 \begin{figure}[h!]
 \centering
	\includegraphics[width=0.8\linewidth]{"setting3.pdf"}
	\caption{\textbf{Scheme of the setting} \textbf{a} The domain with ''Chern number'' $C=1$ (blue in figure) is drifting inside a topologically trivial sample (grey $C=0$). The dark blue arrows indicate the drift direction of the topologically non-trivial domain. \textbf{b} The initial distribution of the local Chern  \cite{bianco2011mapping}. \textbf{c} The evolution of the on-site magnetic field $m$, controlling the phase in the model~\ref{eq:Ch_eq}. }
	\label{fig:setting}
\end{figure}

The space and time dependency of $m$ corresponds to the slow drift of a topological domain inside a system in a trivial phase (see Fig.~\ref{fig:setting}). Different phases are created by selecting the parameter $m$ as dictated by the Eq. (\ref{m_det}). In our case $m_1=-3$ and $m_2=-1$ were chosen for trivial and topological phase correspondingly. Initially the system is prepared in the ground state of $H(0)$. For the domain to move, we change the parameter $m$ at the boundaries of two different phases on the lattice from $m_1$ to $m_2$. The time protocol for such process can be expressed as follows:
\begin{align}  \label{eq:m_par}
m(x, t)= \left  \{ \begin{array}{lllll}
m_1, & x<x_l+[t / \tau] \\
m_2(1-t/\tau)+m_1\cdot t/\tau & x=x_l+[t/\tau] \\
m_2, & x_l+[t/\tau]< x< x_r+[t/\tau] \\
m_1(1-t/\tau)+m_2\cdot t/\tau & x=x_r+[t/\tau] \\
m_1, & x>x_r+[t/\tau]
\end{array} \right.
\end{align}
where $x_l$, $x_r$ denote two sites on the left and right border at the beginning of time protocol (See Fig.~\ref{fig:setting}.) and $\tau$ determine the period of protocol when the domain shifts on one site. \\

\textit{Unitary evolution.}
The evolution of the system is governed by the time-ordered exponent:
\begin{equation} \label{eq:unit_evol}
| \psi (t) \rangle = T\left\{\exp \left(-\frac{i}{\hbar} \int_{t_{0}}^{t} H\left(t^{\prime}\right) d t^{\prime}\right)\right\} | \psi (t_0) \rangle
\end{equation}
Numerically we approximated the state $| \psi (t) \rangle$ using the Trotter-Suzuki decomposition \cite{trotter1959product,suzuki1976generalized}.

\textit{Electric current.} For detecting edge modes we used the following operator for electrical currents
\begin{equation}
\hat{J}(\bm{r}_2, \bm{r}_2)= i \sum_{\sigma,\sigma'}\left( c_{\bm{r}_2, \sigma}^{\dagger} t(\bm{r}_1,\bm{r}_2) c_{\bm{r}_1, \sigma'}-\text { h.c. }\right) .
\end{equation}

\input{output.bbl}

\section*{Acknowledgements}

We thank Peru d'Ornellas for many insightful discussions and in particular for sharing the details of the ongoing work \cite{peru}. Useful comments from Oleg Dubinkin are gratefully acknowledged. This work was carried out in the framework of the Russian Quantum Technologies Roadmap. A.A.M. was also also supported by the "Basis" foundation under Grant No. 18-3-01

\section*{Author contributions statement}

D.B.G. and A.R.Y. performed the numerical simulations. A.A.M. initiated and directed the project. Eq.~\ref{eq:nq_streda} is due to A.N.R.  All the authors contributed to the results analysis and writing the manuscript.

\end{document}


text:
\part{\par \addvspace{4ex}\@afterindentfalse
 \secdef\@part\@spart}
\def\@part[#1]#2{\ifnum \c@secnumdepth >\m@ne
  \refstepcounter{part}
  \addcontentsline{toc}{part}{Part \arabic{part}: #1}
 \else \addcontentsline{toc}{part}{#1}
 \fi
 {\parindent 0pt \raggedright
  \ifnum \c@secnumdepth >\m@ne
   \reset@font\large\rm PART
   \ifcase\arabic{part} \or ONE \or TWO \or THREE \or FOUR \or FIVE
    \or SIX \or SEVEN \or EIGHT \or NINE \or TEN \else \fi
   \par \nobreak
  \fi
  \reset@font\LARGE \rm #2 \markboth{}{}\par }
 \nobreak \vskip 3ex \@afterheading}
\def\@spart#1{{\parindent 0pt \raggedright
  \reset@font\LARGE \rm #1\par}
 \nobreak \vskip 3ex \@afterheading}

\newif\if@levelone

\def\@startsection#1#2#3#4#5#6{\if@noskipsec \leavevmode \fi
   \par \@tempskipa #4\relax
   \@afterindenttrue
   \ifnum #2=\@ne
     \global\@levelonetrue
   \else
     \global\@levelonefalse
   \fi
   \ifdim \@tempskipa <\z@ \@tempskipa -\@tempskipa \@afterindentfalse\fi
   \if@nobreak
     \everypar{}\ifnum#2=2 \vskip 0pt plus1pt\fi
   \else
     \addpenalty{\@secpenalty}\addvspace{\@tempskipa}
   \fi \@ifstar
     {\@ssect{#3}{#4}{#5}{#6}}{\@dblarg{\@sect{#1}{#2}{#3}{#4}{#5}{#6}}}}

\def\SFB@hangraggedright{\rightskip=\@flushglue \let\\=\@centercr \parindent=0pt}

\newif\if@firstsection
\@firstsectiontrue
\def\section{%
 \if@firstsection \fixfootnotes\@firstsectionfalse \fi%
 \@startsection{section}{1}{\z@}
 {-24pt plus -12pt minus -1pt}{6pt}
 {\SFB@hangraggedright\reset@font\normalsize\bf}}
\def\@startsection{subsection}{2}{\z@{\@startsection{subsection}{2}{\z@}
 {-18pt plus -9pt minus -1pt}{6pt}
 {\SFB@hangraggedright\reset@font\normalsize\bf}}
\def\@startsection{subsubsection}{3}{\z@{\@startsection{subsubsection}{3}{\z@}
 {-18pt plus -9pt minus -1pt}{6pt}
 {\SFB@hangraggedright\reset@font\normalsize\it}}
\def\@startsection{paragraph}{4}{\z@{\@startsection{paragraph}{4}{\z@}
 {12pt plus 2.25pt minus 1pt}{-0.5em}{\reset@font\normalsize\bf}}
\def\@startsection{subparagraph}{5}{\parindent{\@startsection{subparagraph}{5}{\parindent}
 {12pt plus 2.25pt minus 1pt}{-0.5em}{\reset@font\normalsize\it}}
\setcounter{secnumdepth}{4}
\def\@sect#1#2#3#4#5#6[#7]#8{%
  \ifnum #2>\c@secnumdepth
    \def\@svsec{}%
  \else
    \refstepcounter{#1}%
    \if@levelone
      \ifSFB@appendix
        \edef\@svsec{}%
      \else
        \edef\@svsec{\csname the#1\endcsname\hskip 1em}%
      \fi
    \else
      \edef\@svsec{\csname the#1\endcsname\hskip 1em}%
    \fi
  \fi
  \@tempskipa #5\relax
  \ifdim \@tempskipa>\z@
    \begingroup #6\relax
    \if@levelone
      \ifSFB@appendix
        \@hangfrom{\hskip #3\relax\@svsec}{\interlinepenalty \@M
        APPENDIX \csname the#1\endcsname:\hskip 0.5em\uppercase{#8}\par}%
      \else
        \@hangfrom{\hskip #3\relax\@svsec}{\interlinepenalty \@M
        \uppercase{#8}\par}%
      \fi
    \else
      \@hangfrom{\hskip #3\relax\@svsec}{\interlinepenalty \@M #8\par}%
    \fi
  \endgroup
   \csname #1mark\endcsname{#7}%
   \addcontentsline{toc}{#1}{\ifnum #2>\c@secnumdepth \else
   \protect\numberline{\csname the#1\endcsname}\fi #7
 \else
  \def\@svsechd{#6\hskip #3\@svsec \if@levelone \uppercase{#8}\else #8\fi
  \csname #1mark\endcsname{#7}
  \addcontentsline{toc}{#1}{\ifnum #2>\c@secnumdepth \else
  \protect\numberline{\csname the#1\endcsname}\fi#7
  }\fi
 \@xsect{#5}}

\def\@ssect#1#2#3#4#5{%
  \@tempskipa #3\relax
    \ifdim \@tempskipa>\z@
      \begingroup #4\@hangfrom{\hskip #1}{%
      \interlinepenalty \@M
      \if@levelone
        \uppercase{#5}%
      \else {#5}%
      \fi\par}%
     \endgroup
    \else
      \def\@svsechd{#4\hskip #1\relax
      \if@levelone
        \uppercase{#5}%
      \else {#5}%
      \fi}
    \fi
  \@xsect{#3}%
}

\newif\ifSFB@appendix
\def
\section{Introduction}

The vast majority of stars do not form alone. They appear to form in a
hierarchy in structures of tens to tens of thousands of stars
(\citealp{testi00,gutermuth05,sanchez07,
andre07,goldsmith08,
andre10,
bressert10,difrancesco10,gutermuth09,
gouliermis10}).  Such structures are a natural consequence of
the gravo-turbulent model of star formation (e.g. \citealp{klessen00,
bonnell01, bonnell03,bate03, bonnell08, bate09, offner09}).

Hierarchical distributions are not in equilibrium, and will rapidly
dynamically evolve into dense star clusters or loose associations
(e.g. \citealp{aarseth72, Goodwin1998,bate98,boily99,kroupa03,goodwin04,
allison09, allison10, moeckel09, Fellhauer2009,gieles10}).
Such evolution is especially violent if the stars are initially dynamically cool
(see \citealp{allison09,allison10}) as many observations suggest they are (e.g. \citealp{Walsh04,
difrancesco04,peretto06,walsh07,andre07,kirk07,gutermuth08b}).

In a clustered phase many
interesting processes may occur such as rapid dynamical mass
segregation (\citealp{allison09,allison10}), binary disruption and
modification (\citealp{heggie74,kroupa95,parker09}), the
formation of high-order multiples like the Trapezium system (\citealp{aarseth74,
zinnecker08,allison10}), and star-disc interactions
affecting planetary system formation (\citealp{boffin98,watkins98,pfalzner05,thies05,thies10}).
Therefore an understanding of the collapse of hierarchical distributions is important
to understand the formation of star clusters and the possible
importance of these effects.

The evolution of initially substructured stellar distributions into
smooth star clusters has been studied by many authors (e.g. \citealp{aarseth72,
Goodwin1998,boily99,kroupa03,goodwin04,allison09,Fellhauer2009,
moeckel09,Smith2011}).  In particular, \def\citename##1{##1}\@internalcite{Fellhauer2009}
attempted to quantify how rapidly an initially clumpy
distribution could evolve into a smooth star cluster.

In this paper we particularly extend the work of \def\citename##1{##1}\@internalcite{Fellhauer2009}
and \def\citename##1{##1}\@internalcite{Smith2011} to investigate how a collapsing clumpy
distribution in a static gas potential is able to erase its
substructure and form a smooth cluster.  \def\citename##1{##1}\@internalcite{Fellhauer2009}
presented a semi-analytic model for the erasure of substructure, but
they did not properly account for two-body effects in their
simulations.  Here we revisit their analysis with an accurate $N$-body
code.

In Section~2 we present our initial conditions, in Section~3
we study the erasure of substructure before examining the rates at
which substructure is erased in Section~4.  Finally we discuss our
results in Section~5, and draw our conclusions in Section~6.

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[height=2.4in,width=3.2in]{rep1.eps} &
\includegraphics[height=2.4in,width=3.2in]{rep2.eps}
\\ \includegraphics[height=2.4in,width=3.2in]{rep3.eps} &
\includegraphics[height=2.4in,width=3.2in]{rep4.eps} \\
\end{array}$
\end{center}
\refstepcounter\@captype \@dblarg{\@caption\@captype}{An xyz-plot of evolution of substructure in a standard model
  simulation. The upper-left panel shows the initial stellar
  distribution (t=0.0 Myr). Stars are initially distributed in well
  defined clumps. The standard model star-forming region has a
  crossing-time $T_{\rm{cr}}^{\rm{sc}}=260$ kyr. As the simulation
  evolves, we show snapshots at 2 (upper-right panel), 4 (lower-left
  panel), and 10 (lower-right panel) crossing-times. Individual clumps
  `puff-up' significantly in less than 2 crossing-times. By 4
  crossing-times, only a few clearly defined clumps are visible, and by 10
  crossing-times almost all substructure has been erased.}
\label{repsim}
\end{figure*}

\section{Initial conditions}

We perform our $N$-body simulations using the direct $N$-body integration
code {\sc{nbody6}} (\citealp{aarseth03}).  The advantage of {\sc{nbody6}} is that
it is able to rapidly and accurately model stellar dynamics, and
two-body encounters in particular.

Our initial conditions are similar to those of \citealp{Fellhauer2009}
(hereafter F09).  Our young star forming regions have a total mass
of $1000 M_\odot$.  We assume they convert gas to stars with an
efficiency $\epsilon$ which ranges between 0.1 and 0.8. Therefore the
mass of stars is $\epsilon
\times 1000 M_\odot$, and the mass of gas $(1 - \epsilon)
\times 1000 M_\odot$.

We simulate the gas with a static background Plummer
potential with a Plummer scale radius of $R^{\rm{sc}}_{\rm{pl}}$ which
ranges between 0.02 and 1~pc.  We set a limiting cutoff radius for the
gas potential of $ 5 \times R^{\rm{sc}}_{\rm{pl}}$.  We
acknowledge that this is not ideal, as the gas will also dynamically
evolve.  However, it is presently impossible to simply model live
background gas and so we follow previous studies in including a static
background (e.g. \citealp{moeckel09}, F09 ,\citealp{Smith2011})

The stars are distributed within the gas potential in $N_0$ subclumps
which follow the underlying gas Plummer distribution.  $N_0$ ranges
from 4 to 32 resulting in a mass per clump of $M_{\rm pl} = 6$ to $80 M_\odot$ (where
the clump stellar mass is $M_{\rm pl} = (\epsilon \times 1000)/N_0
M_\odot$).

Subclumps are distributed within the Plummer sphere according to the
prescription of \def\citename##1{##1}\@internalcite{aarseth74}.  Their bulk
velocities are then scaled to a desired virial ratio $Q_{\rm{i}} = T/|\Omega|$
(where $T$ is the total kinetic energy and $\Omega$ the total
potential energy) where $Q_{\rm{i}}=0.5$ is virial equilibrium and our scaling
ranges between $Q_{\rm{i}}=0$ and 0.5.

Each clump is assumed to be a virialised Plummer sphere with a Plummer scale
radius of $R_{\rm pl} = 0.01$~pc and a cut-off radius beyond which no
stars are placed of $R_{\rm cut} = 5 R_{\rm pl} =  0.05$~pc.  We assume
that sub-clumps are virialised initially as their relaxation time is
so short that they will rapidly virialise (but this process is also
effective at destroying clumps as we shall see).

We take equal-mass stars of mass $0.5 M_\odot$ (roughly the average
mass of a star from a standard IMF).  This means that our clumps
contain from 12 to 160 stars depending on the values of $\epsilon$ and
$N_0$.  Again we acknowledge that equal-mass stars are not realistic.
In particular, differences in stellar masses will have a significant
effect on two-body encounters which we are attempting to examine in
particular detail.  However, introducing a range of stellar masses
would significantly increase stochasticity and add another free
parameter, so we choose to ignore it for now.

Therefore the important parameters are:
\begin{itemize}
\item star formation efficiency $\epsilon$;
\item gas Plummer scale radius $R^{\rm{sc}}_{\rm{pl}}$;
\item number of subclumps $N_0$;
\item virial ratio of the stellar distribution $Q_{\rm{i}}$.
\end{itemize}
From these parameters it is possible to calculate a number of very
useful quantities.

The filling factor, $\alpha$, is the fraction of the volume which
contains subclumps
\begin{equation}
\alpha=\frac{R_{\rm{pl}}}{R^{\rm{sc}}_{\rm{pl}}}\rm.
\end{equation}

The crossing time of the whole system, $T^{\rm sc}_{\rm cr}$, or of an
individual clump, $T_{\rm cr}$ is the typical time taken to cross the
whole system or individual clump (the typical size divided by the
typical speed).

The two-body relaxation time is a measure of how rapidly the internal
velocities of a clump will change by order their own magnitude and is
given by

\begin{equation}
t_{\rm{relax}}=0.1 \frac{N_{\rm{part}}}{{\rm{ln}}(N_{\rm{part}})}
t_{\rm{cr}}\rm,
\end{equation}
where $N_{\rm{part}}$ is the total number of particles in the system,
and $t_{\rm{cr}}$ is the crossing-time of the system.

A list of all the simulations and their parameters can be seen
in Table \ref{pars}.  Three different random realisations of each
parameter set are performed.

\begin{table*}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|c|}
\hline $\alpha$ & $\epsilon$ & $N_0$ & $R^{\rm{sc}}_{\rm{pl}}$&
$R^{\rm{sc}}_{\rm{cut}}$ & $M^{\rm{sc}}_{\rm{pl}}$ &
$T^{\rm{sc}}_{\rm{cr}}$ & $Q_{\rm{i}}$ & $\sigma^{\rm{sc}}_{\rm{3D}}$
& $M_{\rm{star}}$ & $M_{\rm{gas}}$ & $R_{\rm{pl}}$ & $R_{\rm{cut}}$ &
$M_{\rm{pl}}$ & $T_{\rm{cr}}$ & t$_{\rm{relax}}$\\ & & & (pc) & (pc) &
(M$_\odot$) & (kyr) & & (km s$^{-1}$) & (M$_\odot$) & (M$_\odot$) &
(pc) & (pc) & (M$_\odot$) & (kyr) & (kyr)\\ \hline 0.05 & 0.32 & 16 &
0.20 & 1.00 & 1000 &  260 & 0.5 & 2.5 & 320 & 680 & 0.01 & 0.05 & 20.0
& 20 & 21.7 \\  0.05 & 0.32 & 16 & 0.20 & 1.00 & 1000 &  260 &
0.$\dot{3}$ & 2.0 & 320 & 680 & 0.01 & 0.05 & 20.0 & 20 & 21.7\\ 0.05
& 0.32 & 16 & 0.20 & 1.00 & 1000 &  260 & 0.1 & 1.4 & 320 & 680 & 0.01
& 0.05 & 20.0 & 20 & 21.7\\ 0.05 & 0.32 & 16 & 0.20 & 1.00 & 1000 &
260 & 0.0 & 0.0 & 320 & 680 & 0.01 & 0.05 & 20.0 & 20 & 21.7\\ \hline
0.01 & 0.32 & 16 & 1.00 & 5.00 & 1000 & 2950 & 0.5 & 1.1 & 320 & 680 &
0.01 & 0.05 & 20.0 & 20 & 21.7\\ 0.02 & 0.32 & 16 & 0.50 & 2.50 & 1000
& 1043 & 0.5 & 1.6 & 320 & 680 & 0.01 & 0.05 & 20.0 & 20 & 21.7\\ 0.10
& 0.32 & 16 & 0.10 & 0.50 & 1000 &   93 & 0.5 & 3.6 & 320 & 680 & 0.01
& 0.05 & 20.0 & 20 & 21.7\\  0.20 & 0.32 & 16 & 0.05 & 0.25 & 1000 &
33 & 0.5 & 5.0 & 320 & 680 & 0.01 & 0.05 & 20.0 & 20 & 21.7\\ 0.50 &
0.32 & 16 & 0.02 & 0.10 & 1000 &    8 & 0.5 & 8.0 & 320 & 680 & 0.01 &
0.05 & 20.0 & 20 & 21.7\\ \hline 0.05 & 0.32 &  4 & 0.20 & 1.00 & 1000
&  260 & 0.5 & 2.5 & 320 & 680 & 0.01 & 0.05 & 80.0 & 10 & 31.5
\\ 0.05 & 0.32 &  8 & 0.20 & 1.00 & 1000 &  260 & 0.5 & 2.5 & 320 &
680 & 0.01 & 0.05 & 40.0 & 15 & 27.4 \\ 0.05 & 0.32 & 32 & 0.20 & 1.00
& 1000 &  260 & 0.5 & 2.5 & 320 & 680 & 0.01 & 0.05 & 10.0 & 29 & 19.4
\\ \hline 0.05 & 0.10 & 16 & 0.20 & 1.00 & 1000 &  260 & 0.5 & 2.5 &
96 & 904 & 0.01 & 0.05 &  6.0 & 37 & 17.9 \\ 0.05 & 0.20 & 16 & 0.20 &
1.00 & 1000 &  260 & 0.5 & 2.5 & 200 & 800 & 0.01 & 0.05 & 12.5 & 26 &
20.2 \\ 0.05 & 0.25 & 16 & 0.20 & 1.00 & 1000 &  260 & 0.5 & 2.5 & 248
& 752 & 0.01 & 0.05 & 15.5 & 24 & 21.7 \\ 0.05 & 0.50 & 16 & 0.20 &
1.00 & 1000 &  260 & 0.5 & 2.5 & 496 & 504 & 0.01 & 0.05 & 31.0 & 17 &
25.5 \\ 0.05 & 0.60 & 16 & 0.20 & 1.00 & 1000 &  260 & 0.5 & 2.5 & 600
& 400 & 0.01 & 0.05 & 37.5 & 15 & 26.1 \\ \hline  0.10 & 0.10 & 16 &
0.10 & 0.50 & 1000 &   93 & 0.5 & 3.6 & 96 & 904 & 0.01 & 0.05 &  6.0
&  37 & 17.9 \\ 0.10 & 0.20 & 16 & 0.10 & 0.50 & 1000 &   93 & 0.5 &
3.6 & 200 & 800 & 0.01 & 0.05 & 12.5 &  26 & 20.2 \\ 0.10 & 0.25 & 16
& 0.10 & 0.50 & 1000 &   93 & 0.5 & 3.6 & 248 & 752 & 0.01 & 0.05 &
15.5 &  24 & 21.7 \\ 0.10 & 0.50 & 16 & 0.10 & 0.50 & 1000 &   93 &
0.5 & 3.6 & 496 & 504 & 0.01 & 0.05 & 31.0 &  17 & 25.5 \\ 0.10 & 0.70
& 16 & 0.10 & 0.50 & 1000 &   93 & 0.5 & 3.6 & 696 & 304 & 0.01 & 0.05
& 43.5 &  14 & 26.1 \\ \hline 0.05 & 0.80 & 32 & 0.20 & 1.0 & 1000 &
260 & 0.5 & 2.5 & 800 & 200 & 0.01 & 0.05 & 25.0 &  18 & 22.9 \\ 0.05
& 0.40 & 16 & 0.20 & 1.0 & 1000 &   260 & 0.5 & 2.5 & 400 & 600 & 0.01
& 0.05 & 25.0 &  18 & 22.9 \\ 0.05 & 0.20 & 8 & 0.20 & 1.0 & 1000 &
260 & 0.5 & 2.5 & 200 & 800 & 0.01 & 0.05 & 25.0 &  18 & 22.9
\\     \hline
\end{tabular}
\refstepcounter\@captype \@dblarg{\@caption\@captype}{A complete list of the parameters of all the simulations in
  our parameter study. The table is split by horizontal lines into
  sets (set 1 to set 6 from top to bottom). Each set is chosen to test
  the influence of a specific parameter on the formation rate of the
  cluster (see text for further details). Columns give the filling
  factor $\alpha$, the SFE $\epsilon$, number of clumps $N_0$,
  followed by the Plummer radius $R^{\rm{sc}}_{\rm{pl}}$, the cut-off
  radius $R^{\rm{sc}}_{\rm{cut}}$, the total mass
  $M^{\rm{sc}}_{\rm{pl}}$, and the crossing time of the star-forming
  region $T^{\rm{sc}}_{\rm{cr}}$. The following two columns are
  initial Virial ratio $Q_{\rm{i}}$, and corresponding velocity
  dispersion of clumps with respect to their clumps within the region
  $\sigma^{\rm{sc}}_{\rm{3D}}$. The next two columns are the mass in
  stars $M_{\rm{star}}$ and mass in gas $M_{\rm{gas}}$ (modelled as an
  analytical background) within the star-forming region. Finally we
  show the Plummer radius $R_{\rm{pl}}$, the cut-off radius
  $R_{\rm{cut}}$, the mass $M_{\rm{pl}}$, the crossing time
  $T_{\rm{cr}}$, and the relaxation time of an individual clump
  $t_{\rm{relax}}$.}
\label{pars}
\end{table*}

\section{Results}
\label{results}

First we shall examine in detail our `standard model'.  This is a
virialised $Q_{\rm{i}} = 0.5$ cluster with
a star formation efficiency of $\epsilon = 0.32$ and $N_0 = 16$.  This
cluster has a stellar mass of $320 M_\odot$ in 16 $20 M_\odot$ clumps
(40 equal-mass stars per clump).  The cluster has a filling factor of
$\alpha = 0.05$ and a total crossing time of 260~kyr.  Each clump
has a crossing time of 20~kyr, and a relaxation time of 22~kyr.

Figure~\ref{repsim} shows the evolution of the standard model for 10
crossing times (2.6~Myr).  Initially (t=0 Myr) the stars have a highly clumpy and
sub-structured distribution. By 0.5~Myr the initial structure is
already less obvious, and by $>1$~Myr the cluster has a fairly smooth
appearance with very little evidence of the initial clumps.

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[height=3.2in,width=3.2in]{nbody6_30.eps} &
\includegraphics[height=3.2in,width=3.2in]{sbox_30.eps} \\
\end{array}$
\end{center}
\refstepcounter\@captype \@dblarg{\@caption\@captype}{The resulting stellar distribution after 0.7 Myr of evolution
  from the standard model initial conditions in; an {\sc{nbody6}} simulation
  (left panel), and a {\sc{superbox}} simulation (right panel). Substructure
  has been erased more rapidly in the {\sc{nbody6}} simulations as a result
  of a realistic handling of two-body relaxation.}
\label{nbody6vssbox}
\end{figure*}

It is expected that the erasure of substructure would be due to one,
some, or all of the following mechanisms.
\begin{enumerate}
\item Internal scattering and the ejection of stars from a clump by internal
  two-body interactions.
\item Tidal stripping of clumps by the gas potential.
\item Tidal encounters between clumps.
\item Collisions between clumps or stars.
\end{enumerate}

An examination of Figure~\ref{repsim} suggests that internal
scattering is a crucial factor in the erasure of substructure.  By
0.5~Myr, the appearance of the cluster is already quite smooth (we
shall return later to quantify the erasure of substructure, but for now we
will use a `by eye' examination).  In only two system crossing times
interactions between clumps cannot have been important.  The high
density of the initial clumps also suggests that tidal stripping by
other clumps or the gas potential can not have been responsible for the
smoothness.  The only process that works so rapidly on clumps is
internal relaxation as the clumps are 20 internal relaxation times old
by 0.5~Myr.

The effect of internal relaxation is to eject stars from the clumps
forming a smooth background of stars, and also to increase the size of
the clumps (`puffing-up').  The puffing-up of clumps is obvious in
Figure~\ref{repsim} where even isolated clumps at 0.5~Myr are clearly
larger than initially. Binary formation within clumps is observed. This likely plays a significant role in enhancing scattering and ejection of stars from clumps.

As well as introducing a background of ejected stars, puffing-up has
two effects which further increase the rate of clump mergers.  As
clumps are larger they are more susceptible to tidal stripping, and
the filling factor increases -- significantly increasing the rate of
clump collisions.

In the clump merger simulations of F09 these effects arising from
two-body encounters were missed as the {\sc{superbox}} code used for these simulations
damps two-body encounters entirely.  This is illustrated in Figure
\ref{nbody6vssbox} showing two simulations with the same initial
conditions\footnote{It should be noted that the {\sc{superbox}} runs have
  $10^5$ much lower-mass particles initially in each clump and here we
  display only 40 for a fair comparison.}  having evolved for 0.7
Myr (2.7 crossing-times). Therefore it is worth revisiting
the results of F09 when applied to star clusters in light of this
vitally important physical process.  We note that the F09 results
apply to situations where two-body encounters can be considered
negligible.

\begin{figure}
  \centering \epsfxsize=9cm \epsfysize=7cm \epsffile{set1comp.eps}
  \refstepcounter\@captype \@dblarg{\@caption\@captype}{A graph of evolution of number of clumps N with time $\tau$
    (in units of the star-forming region crossing time). Both curves
    are results for simulations with the standard model initial
    conditions. The solid (black) line is the results for the {\sc{nbody6}}
    run as measured using a minimum spanning tree based clump
    finder. The dashed (red) line is the results for the {\sc{superbox}}
    run. As can be qualitatively assessed by eye in Figure
    \ref{nbody6vssbox}, substructure is more quickly erased in the
    {\sc{nbody6}} simulations as a result of a realistic handling of two-body
    relaxation. Equivalently, the {\sc{nbody6}} star cluster forms more
    rapidly.}
  \label{set1comp}
\end{figure}

\@startsection{subsection}{2}{\z@{Clump counting with a Minimum Spanning Tree}

A minimum spanning tree (MST) is the shortest path linking a set of
$n$ points with no closed loops.  The MST is useful in that it always
has $n-1$ connections (edges) with a unique total length (see e.g. \def\citename##1{##1}\@internalcite{cartwright04}
 or \def\citename##1{##1}\@internalcite{allison09} for uses of the MST in
astronomy).  We use the algorithm described in \def\citename##1{##1}\@internalcite{allison09} to construct an
MST for a simulation at some point in time.

In order to identify clumps we introduce a cutting length
$l_{\rm{cut}}$.  If an edge has a length greater than $l_{\rm{cut}}$
then it is removed and we examine the subset of connections remaining
after cutting the longest edges.  If a subset contains more than
one-third of the stars initially within a subclump then we define it
as a clump.  (Note that \def\citename##1{##1}\@internalcite{gutermuth09} use a similar
method, as do many friends-of-friends clump finders.)  However, we
introduce an extra element as our clumps are not just
close in physical space, but in phase space.  Therefore we build our
MST in 6D phase space before applying the cut.

A significant problem with clump finding is that it can be very
sensitive to the cutting length used.  In our initial conditions we
have the luxury of knowing what clumps are present and where they
are.  This allows us to fine tune our cutting length to get the right
answer initially (a bad choice of cutting length can result in
garbage).  We also check a number of simulations by eye to see that
the structures selected as clumps are indeed clumps, and that no
structure have been missed.  This is not ideal, but the best we can do
at the moment with no way of properly selecting a cutting length
(\def\citename##1{##1}\@internalcite{gutermuth09} propose a way of selecting the cutting length,
but this only works well when structures are distinct).
\newline
\indent An illustration of the method is given in Figure \ref{set1comp} which
quantifies the evolution of the {\sc{nbody6}} and {\sc{superbox}} simulations
illustrated in Figure \ref{nbody6vssbox}.  The solid lines shows the rapid
decrease in the number of clumps in the {\sc{nbody6}} simulation compared to
those in the {\sc{superbox}} simulation in the dashed line.  This agrees with
the by eye assessment of the far more rapid erasure of substructure in
the former simulation (Figure \ref{repsim}).

\section{A parameter study of cluster formation rates}

With a quantitative measure of the evolution of the substructure and a
better understanding and modelling of the physical processes behind
this erasure, we can revisit the study of F09 to quantify the rate at
which structure is erased.

F09 suggested that the rate of erasure of substructure measured by the
number of clumps at time $\tau$, $N(\tau)$ compared to the initial
number of subclumps $N_0$ could be well-fitted by an equation of the
form
\begin{equation}
\label{expeqn}
N(\tau)=(N_0-1) {\rm{exp}}(-\eta \tau)+1 \rm,
\end{equation}
where $\eta$ is a free parameter that depends on the initial
conditions of the simulation.  A large value of $\eta$ corresponds to
a rapid loss of substructure and the rapid appearance of a smooth
cluster.  Therefore we refer to $\eta$ as the `cluster formation rate'.

As shown in Figure \ref{expfit} we agree with F09 that an exponential
decay is indeed a good fit to the rate at which substructure is lost.

\begin{figure}
  \centering \epsfxsize=9cm \epsfysize=7cm \epsffile{nclumpexpfit.eps}
  \refstepcounter\@captype \@dblarg{\@caption\@captype}{A graph of evolution of number of clumps N with time $\tau$
    (in units of the star-forming region crossing time). The solid
    (black) curve is the results of a Set 4 simulation as measured
    using the minimum spanning tree based clump finder. The dashed
    (red) line is an example of an exponential fit, using Equation
    \ref{expeqn}, to the simulation results. As presented in the key,
    this fit provides us with a value for the cluster formation rate
    parameter $\eta$. A high resulting value of $\eta$ indicates that
    substructure has been erased rapidly, and thus that a smooth
    cluster has formed rapidly.}
  \label{expfit}
\end{figure}

\@startsection{subsection}{2}{\z@{Results of the parameter study on cluster formation rate}
\label{parstudyresults}

In Figure \ref{varypars} we show the variation of the cluster
formation rate $\eta$ with different parameters in our study (see
Table \ref{pars}).  Please note the differing y-axis scales in the
panels of Figure \ref{varypars}.  Remember also that each clump has the
same size in each simulation.

\@startsection{subsection}{2}{\z@{The initial Virial ratio $Q_{\rm{i}}$}
The upper-left panel of Figure \ref{varypars} shows how cluster
formation rate depends on the initial virial ratio. As the initial
virial ratio becomes increasingly sub-virial (cool), the cluster formation rate
steadily increases. This occurs as a result of increased clump-clump
interactions - for sub-virial initial conditions clumps tend to fall
into a more compact configuration within the gas potential. Cluster
formation rate steeply rises as $Q_{\rm{i}}$ falls below 0.1. In this
case, clumps tend to fall towards the centre of the gas potential well
on time-scales $\sim1$ free-fall time, resulting in multiple
simultaneous clump-clump collisions. {\it{The initial Virial ratio
    $Q_{\rm{i}}$ is a key parameter controlling the cluster formation
    rate}}.

\@startsection{subsection}{2}{\z@{The filling factor $\alpha$}
The upper-right panel of Figure \ref{varypars} shows how cluster
formation rate depends on the filling factor $\alpha$. For
$\alpha<0.2$, the cluster formation rate is fairly constant as it is
dominated by the ejection of stars from clumps rather than clump-clump
interactions. However
for very high filling factors ($\alpha\sim0.5$) the cluster formation
rate becomes high. When filling factor is this high, clumps are
almost overlapping in the initial conditions, and this is enhanced by
clumps puffing-up, hence rapid merging
occurs.

\begin{figure*}
\begin{center}$
\begin{array}{cc}
\includegraphics[height=6.8cm,width=9cm]{varyQ.eps} &
\includegraphics[height=6.8cm,width=9cm]{varyAlpha.eps}
\\ \\  \includegraphics[height=6.8cm,width=9cm]{varyN0.eps} &
\includegraphics[height=6.8cm,width=9cm]{varySfe.eps}
\\ \\ \includegraphics[height=6.8cm,width=9cm]{fixtrelax.eps} &
\includegraphics[height=6.8cm,width=9cm]{varytrelax.eps} \\
\end{array}$
\end{center}
\refstepcounter\@captype \@dblarg{\@caption\@captype}{Plots of the dependency of the cluster formation rate
  parameter $\eta$ (y-axis) versus the parameters investigated in the
  parameter study (x-axis); initial Virial ratio $Q_{\rm{i}}$, filling
  factor $\alpha$ (upper-right), the initial number of clumps $N_0$
  (middle-left panel), region star formation efficiency $\epsilon$
  (middle-right panel), region star formation efficiency $\epsilon$
  for fixed clump relaxation time $t_{\rm{relax}}$ (bottom-left
  panel), and clump relaxation time $t_{\rm{relax}}$ (bottom-right
  panel). A description and discussion of these results is provided in
  the text (see Section \ref{parstudyresults}). Please note varying
  scale on y-axis. The two upper panels have a matching y-scale to
  each other. The four panels beneath them, all have matching y-scale
  to each other.}
\label{varypars}
\end{figure*}

\@startsection{subsection}{2}{\z@{The the initial number of clumps $N_0$ and star formation efficiency $\epsilon$}
The middle-left panel of Figure \ref{varypars} shows how the cluster
formation rate depends on the the initial number of clumps $N_0$, and the
middle-right panel shows how cluster formation rate depends on star
formation efficiency $\epsilon$.  At a first glance, the cluster formation
rate appears to depend on both $N_0$ and $\epsilon$. A trend towards
increasing cluster formation rate for increasing $N_0$ is visible in
the middle-left panel. Meanwhile, as star formation efficiency falls,
cluster formation rate increases. This trend is the reverse of
    what is seen in the F09 simulations.

However by varying $N_0$ and $\epsilon$ we vary the initial number of
stars in a clump. This results in a variation in the relaxation time
$t_{\rm{relax}}$ of any individual clump and hence the speed at which
clumps puff-up (see Set 3 and Set 4 of Table
\ref{pars}). As $N_0$ increases, $t_{\rm{relax}}$ falls. Meanwhile as
$\epsilon$ increases, $t_{\rm{relax}}$ grows. Therefore it is
difficult to separate the additional effects of varying
$t_{\rm{relax}}$ on cluster formation rate.

\@startsection{subsection}{2}{\z@{A fixed clump relaxation time $t_{\rm{relax}}$}
In the bottom left panel of Figure \ref{varypars} we show how the
cluster formation rate varies for clumps with a fixed internal
relaxation time (selecting values of $N_0$ and $\epsilon$ such that
the relaxation time is 22.9~Kyr). For a constant $t_{\rm{relax}}$, the
cluster formation rate is roughly constant.  This strongly suggests
that it is the internal relaxation of clumps -- the ejection of stars
to make a background and the puffing-up of clumps to enhance
clump-clump collisions and tidal interactions that are the crucial
physical parameters.

\@startsection{subsection}{2}{\z@{The relaxation time of individual clumps}
This conclusion is further supported by the
bottom-right panel of Figure \ref{varypars} where the cluster
formation rate can be seen to decrease with greater clump internal
relaxation time.  Note that these simulations are all for an initially
virialised clump distribution, and we also exclude the $\alpha=0.5$
filling factor simulation as the cluster formation times in these
simulations are dominated by other effects.  It appears that {\it{the
clump relaxation time is a key parameter controlling the cluster
formation rate}}.

\section{Discussion}
\label{discussion}
As discussed in the introduction, star formation is a messy and complex
process that does not initially produce a smooth, relaxed star
cluster.  To give the roughly spherical, smooth star clusters that we
often observe, sub-structure must be erased (see also \citealp{aarseth72,
Goodwin1998,boily99,kroupa03,goodwin04,allison09,Fellhauer2009}).

The speed at which a cluster can erase its substructure depends on its
initial virial ratio (see also \citealp{goodwin04}), but also
critically on the rate at which sub-structure evolves internally from
its initial state.  Dense, low-$N$ clumps have a short internal
relaxation time and will disperse rapidly (see also \citealp{kroupa03}).
Two-body encounters within a clump eject stars, forming a
general stellar background, as well as causing clumps to increase
significantly in size, making clump-clump interactions more likely, and
making them more susceptible to tidal stripping. However, these
effects are secondary to internal relaxation.

By considering the middle-right
panel of Figure \ref{varypars} we can gauge the importance of clump
merging in our simulations. The data points are from Set 3 and Set 4
of our parameter study. By increasing the star formation efficiency,
we additionally increase the mass of the clumps. For clumps to merge,
their relative impact velocities must be of order the velocity
dispersion of the clumps. Therefore more massive clumps should merge
more easily as seen in the corresponding F09 simulations. Instead we
see the opposite - a decreasing cluster formation rate with increasing
clump mass. This indicates that merging plays a minor role in the
cluster formation process. Instead it is dominated by the effects of
inter-clump two-body encounters. For higher star formation
efficiencies, the initial number of stars within a clump
increases. Thus two-body encounters occur less frequently and
consequently the cluster formation rate falls (as seen in the
middle-right panel of Figure \ref{varypars}).

The theory of clump merging developed in \def\citename##1{##1}\@internalcite{Fellhauer2002} and F09
is therefore not applicable in low-$N$ ($N \sim 1000$) systems such as in
the star-forming regions modelled in this study. However, in high-$N$
systems such as mergers of clusters within cluster complexes, the
effects of two-body encounters are far less important. In
these scenarios, the F09 theory should remain valid.

An important conclusion to draw from this analysis is that clusters
can change their appearance rapidly when they are young.  A cluster
that appears smooth and relaxed at an age of a few~Myr may well not
have formed that way.  A number of authors have recently emphasised
that clusters evolve rapidly and that current conditions are not
always a good indicator of past conditions (e.g. \citealp{bastian08,
allison09}).

\section{Summary $\&$ Conclusions}
Both observations and theory agree that stars form with a complex
clumpy distribution within a star forming
region. Small clumps containing
$\sim10$'s of stars form embedded within the envelope of molecular gas
from which they formed. However, observations of star clusters of $>$
a few~Myr in age often show smooth, relaxed distributions (e.g. the
Orion Nebula Cluster at $\sim 3$~Myr).  In order for a smooth
spherical star cluster to
form, the star-forming region must erase its initial substructure.

We investigate the mechanisms by which substructure is erased by
modelling star forming regions using the {\sc{nbody6}} code. Our stars are initially
distributed in clumps, and embedded in a static potential to mimic the
gravitational influence of the gas envelope on stellar dynamics. We
conduct a parameter study to investigate the key parameters effecting
the rate at which substructure is erased. Key parameters include the
initial Virial ratio of the clumps within the gas potential, the
filling factor of clumps within the gas, the initial number of clumps,
and the star formation efficiency of the total star-forming region.

We find a number of new and different results, than of those presented
in a similar study in \citealp{Fellhauer2009} (F09). These differences
arise predominantly due to the proper treatment of two-body
encounters in our simulations. Our key results may be summarised in
the following:

\begin{enumerate}
\item Clusters form predominantly from stars that are scattered out of
  clumpy substructure, and not by clump merging.
\item As a result, the rate at which a cluster forms is a strong
  function of the relaxation time within the clumps. Unlike in F09,
  the star formation efficiency of the region does not effect the
  cluster formation rate.
\item The initial virial ratio of the clumps is also a key parameter
  controlling the rate at which a cluster forms. The lower the initial
  virial ratio, the more rapidly substructure is erased and a cluster
  forms.
\end{enumerate}

As inter-clump scattering has been demonstrated to be of such
importance to cluster formation, it is vital that models considering
the stellar dynamics of star-forming regions do so correctly. The use
of softened gravity between stars will result in suppression of
two-body encounters, and as such a key channel by which clusters form
will be missed.

Furthermore, if star clusters form by scattering of stars from clumps,
there is an increased likelihood that a moving subclump can leave a
trail of stars which maintain a velocity signature of the clump from
which they originated. Similar trails are reported in simulations of
massive stars clusters within a dark matter halo
(\citealp{Assmann2010}). Such velocity structures may be observable in
young clusters with the advent of Gaia. If so, we anticipate that
these observations could provide strong constraints on the recent
formation history of young clusters. We defer a detailed study of this
topic to a latter paper (\citealp{Smith2011} in prep.).

\section*{Acknowledgements}
MF announces financial support through FONDECYT grant 1095092. RS is financed by GEMINI-CONICYT fund 32080008 and a COMITE MIXTO grant. PA is financed through a CONCICYT PhD Scholarship.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Proposed by Darwin and Fowler \cite{Darwin1922a,Darwin1922b} method of steepest descent is a typical approach for the derivation of canonical (Gibbs) ensembles (see \cite{Fowler1952,Huang1978,Pathria2011}) or other asymptotic problems on the border on combinatorics and statistical physics \cite{Flajolet2009}. The method bases on the properties of the complex integrals and allows one to efficiently calculate the desired limit of the typical combinatorial problems. Thought, typically Darwin-Fowler approach is concerned with the asymptotics, in this article, we focus mainly on the non-asymptotic case. With the introduced combinatorial approach based on Bell polynomials we solve exactly the Darwin-Fowler problem purely combinatorically without integral representation.

 Bell Polynomials have been applied recently to the wide range of problems of statistical physics e.g. in the description of gas of clusters \cite{AFronczak2013,Siudem2013}, partition function for ideal gases \cite{Zhou2018}, series expansion for quantum partition functions \cite{Hen2018}, general lattice models description \cite{AFronczak2014} and Ising model \cite{Siudem2014}. With this work, we complete this list with the application to the canonical ensemble, which additionally reveals unexpected combinatorial origins of this fundamental to statistical mechanics distribution. With the introduced approach we show that for the typical  degeneracies (i.e. constant one and harmonic oscillator's) the most probable configurations are given by the well-known combinatorial Lah and Stirlings numbers.

The paper is composed as follows: In the first section, we introduce the classic formulation of the Darwin-Fowler method. The second section focuses on the definition and basic properties of Bell polynomials and their relationship to Lah and Stirling  numbers. In the third section, we use the derived approach to solving the Darwin-Fowler procedure for the non-asymptotic case.
The work is concluded with the analysis of the few examples of the degeneracy of energy levels, which reveals their relationship with the   Lah and Stirling numbers.

\section{Darwin-Fowler method}
Let us consider an ensemble consisting of the total of $\mathcal{K}$ systems with possible energies $E=\{\varepsilon_1,\,\varepsilon_2,\,\dots\}$ with degeneracies $\Omega=\{\omega_1,\, \omega_2,\,\dots\}$, where the energy level $\varepsilon_\ell$ has degeneracy equal to $\omega_\ell$. Furthermore, we assume that there exists a~quantum of the energy $\varepsilon$, which without loss of generality means that $\varepsilon_\ell= \varepsilon\ell$. The above allows one to formulate the following condition for $\{n_\ell\}$ i.e. the numbers of systems with the energies from $E$
\begin{equation}\label{eq:cond}
    \begin{cases}
    \sum_{\ell=1}^{\infty}n_\ell=\mathcal{K},\\
	\sum_{\ell=1}^{\infty} \varepsilon_\ell n_\ell=\varepsilon\mathcal{N}.
    \end{cases}
    \Longrightarrow\quad\;
        \begin{cases}
    \sum_{\ell=1}^{\infty}n_\ell=\mathcal{K},\\
	\sum_{\ell=1}^{\infty}\ell n_\ell =\mathcal{N}.
    \end{cases}
\end{equation}
In the Darwin-Fowler method we examine the probability distribution over the energy levels
\begin{equation}\label{eq:P-def}
    \mathcal{P}_{\mathcal{N},\,\mathcal{K}}(\ell)=\frac{n^\star_{\mathcal{N},\,\mathcal{K}}(\ell)}{\mathcal{K}},
\end{equation}
where $n^\star_{\mathcal{N},\,\mathcal{K}}(\ell)$ is the number of the systems at $\ell$-th energy level for the average configuration of ensemble consisting of $\mathcal{K}$ systems and having total energy $\varepsilon \mathcal{N}$. We assume that every single configuration occurs with equal probability (with respect to the degeneracies) thus
\begin{equation}\label{eq:oczekiwana}
	n^\star_{\mathcal{N},\,\mathcal{K}} = \frac{\sum_{\{m_r\}_{\mathcal{N},\,\mathcal{K}}}m_\ell W_{\mathcal{N},\,\mathcal{K}}(\{m_r\})}{\sum_{\{m_r\}_{\mathcal{N},\,\mathcal{K}}} W_{\mathcal{N},\,\mathcal{K}}(\{m_r\})},\quad\mathrm{where}\quad	W_{\mathcal{N},\,\mathcal{K}}(\{m_\ell\})=\mathcal{K}!\prod_{\ell}\frac{\omega_\ell^{m_\ell}}{m_\ell!},
\end{equation}
and the summation is taken over all possible configurations satisfying Eq.  (\ref{eq:cond}) and factor $W_{\mathcal{N},\,\mathcal{K}}(\{m_\ell\})$ counts the number of possible realisations of the configuration $\{m_\ell\}$. With the introduced notion one can realize that for the function $\Gamma_{\mathcal{N},\,\mathcal{K}}$ defined as
\begin{equation}\label{eq:Gamma}
	\Gamma_{\mathcal{N},\,\mathcal{K}}(\omega_1,\,\omega_2,\,\dots)=\mathcal{K}!\sum_{\{m_r\}_{\mathcal{N},\,\mathcal{K}}}\left(\frac{\omega_1^{m_1}}{m_1!}\cdot \frac{\omega_2^{m_2}}{m_2!} \cdots\right),
\end{equation}
one can express Eq. (\ref{eq:oczekiwana}) for  the number $n^\star_{\mathcal{N},\,\mathcal{K}}$  in the following way
\begin{equation}
n^\star_{\mathcal{N},\,\mathcal{K}}=\omega_\ell\frac{\partial \ln\Gamma(u_1,\,u_2,\,\dots)}{\partial u_\ell}\Biggr|_{u_r=\omega_r,\,\,r=1,\,2,\,\dots}.\label{eq:wartoscocz}
\end{equation}

As we have noticed, the problem of determining the $P_{\mathcal{N},\,\mathcal{K}}$ distribution can be reduced to determining the value of the function $\Gamma_{\mathcal{N},\,\mathcal{K}}$ and its derivative. Typically, in the Darwin-Fowler method, this is done by using the steepest descend method (see \ref{sec:steepestdescend}). Instead of this approximate approach, we solve the problem exactly, using Bell polynomials, which we introduce in the next section.

\section{Bell Polynomials and their combinatorics}

Bell Polynomials (introduced in \cite{Bell1934}, see also sec. 3.3 in \cite{Comtet1974}) are inseparably linked with the famous Fa\'a di Bruno's formula i.e. the generalization of the chain rule (derivation of the composition of two functions) to higher derivatives. However, the formula has been published \cite{FaadiBruno1855}  30 years before Bell was born, so naturally,  it was stated without that notion (see Eq. (\ref{eq:Riordan})). We, however, use its version presented in Eq. (\ref{eq:Riordan-Bell}), where we introduced Bell polynomials. For the historical background of the Fa\'a di Bruno's formula  see \cite{Johnson2002,Craik2005}.  Fa\'a di Bruno's formula (in clasical formulation in Eq. (\ref{eq:Riordan}) and with the usage of Bell polynomials $B_{n,\,k}$ in Eq. (\ref{eq:Riordan-Bell})) states
\begin{align}\label{eq:Riordan}
\frac{d^n}{dx^n} F(G(x)) &=\sum_{\{m_r\}_n} {\frac {n!}{m_{1}!\,m_{2}!\,\cdots \,m_{n}!}}\cdot F^{(m_{1}+\cdots +m_{n})}(G(x))\cdot \prod _{j=1}^{n}\left[\frac{G^{(j)}(x)}{j!}\right]^{m_{j}}=\\
&= \sum_{k=1}^n F^{(k)}(G(x))\cdot B_{n,k}\left(G'(x),G''(x),\dots,G^{(n-k+1)}(x)\right),\label{eq:Riordan-Bell}
\end{align}
where $F$ and $G$ are analytical functions and summation in Eq. (\ref{eq:Riordan}) is taken over  integers $m_i$ such
\begin{equation*}
1\cdot m_{1}+2\cdot m_{2}+3\cdot m_{3}+\cdots +n\cdot m_{n}=n.
\end{equation*}
Bell Polynomials $B_{nk}$ from Eq.~(\ref{eq:Riordan-Bell}) can be defined in the two equivalent ways -- combinatorically and analytically. Let us recall those two approaches in the following subsections.

\subsection{Analytic definition of Bell Polynomials}

As in the previous section, let us consider two analytic functions $F$ and $G$ such that $F(0)=0$,  which means that
\begin{equation*}
F(x)=\sum _{n=1}^{\infty }f_{n} \frac{x^{n}}{n!},\;\;\; G(x)=\sum _{n=0}^{\infty }g_{n} \frac{x^{n}}{n!}.
\end{equation*}
Then we ask about the series expansion of the composition of both functions, similarly to the consideration in Fa\'a di Bruno's formula (cf. Eqs. (\ref{eq:Riordan}, \ref{eq:Riordan-Bell}))
\begin{equation*}
G(F(x))=\sum _{n=0}^{\infty }c_{n} \frac{x^{n}}{n!}.
\end{equation*}
Let us then define Bell polynomials $B_{n,k}$ and connect them with the coefficients $c_n,\,f_n,\,g_n$ as follows
\begin{equation}\label{eq:Belldef1}
\begin{cases}
c_0=g_0,\\
c_n=\sum _{k=1}^{n}g_{k}B_{n,\,k}(f_{1},\dots ,f_{n-k+1})\;\;\mathrm{for}\;\;n>0.
\end{cases}
\end{equation}
Combining Eqs. (\ref{eq:Riordan}) and (\ref{eq:Belldef1}) one obtains the following analytical definition of Bell polynomials
\begin{equation}\label{eq:Bell-definicja}
B_{n,k}(a_{1},a_{2},\dots ,a_{n-k+1})=\sum_{\{m_r\}_{n,k}} \frac{n!}{m_{1}!m_{2}!\cdots m_{n-k+1}!}\prod _{j=1}^{n}\left(\frac{a_j}{j!}\right)^{m_{j}},
\end{equation}
where the summation is taken over all non-negative integers $\{m_r\}$ which satisfy
\begin{align}\label{eq:m}
\begin{cases}1\cdot m_{1}+2\cdot m_{2}+3\cdot m_{3}+\cdots +n\cdot m_{n}&=n,\\
 m_{1}+ m_{2}+ m_{3}+\cdots + m_{n}&=k.
 \end{cases}
\end{align}
Let us note the similarity of conditions (\ref{eq:cond}) and (\ref{eq:m}), which makes Bell polynomials a natural tool for describing Darwin-Fowler's formalism. Let us also compare definition of function $\Gamma_{\mathcal{N},\,\mathcal{K}}$ given in the Eq. (\ref{eq:Gamma}) with Bell polynomial given by the Eq. (\ref{eq:Bell-definicja}).

\subsection{Combinatorial definition of Bell Polynomials}

Now let us ask seemingly totally different and non-connected to the previous one question: what is the number of possible decompositions of a set of $n$ elements into $k$ clusters (subsets)? Additionally, we assume that every cluster of size $l$ has $a_l\geq 0$ possible configurations. Firstly, let us fix one of the set divisions. Such decomposition is described by the sequence of the non-negative integers  $\{m_r\}$ which satisfy
\begin{align*}
\begin{cases}1\cdot m_{1}+2\cdot m_{2}+3\cdot m_{3}+\cdots +n\cdot m_{n}&=n,\\
 m_{1}+ m_{2}+ m_{3}+\cdots + m_{n}&=k,
 \end{cases}
\end{align*}
where $m_\ell$ describes the number of subsets of size $\ell$. In this situation, the number of possible implementations of such a division is equal to

\begin{equation*}
    \frac{n!}{m_{1}!m_{2}!\cdots m_{n-k+1}!}\prod _{j=1}^{n}\left(\frac{a_j}{j!}\right)^{m_{j}},
\end{equation*}
which, summed over all possible partitions $\{m_r\}_{\mathcal{N},\,\mathcal{K}}$ leads to Bell polynomials. The formal proof of the equivalence between analytic and combinatorial definition of Bell Polynomials can be found in \cite{Comtet1974}, however let us describe this fact in the following example.
\begin{example}
Let us consider two examples of Bell Polynomials with $n=6$, $k=3$ and $n=6$, $k=4$. In both cases the coefficients of the polynomials count the number of a possible partition of the set into clusters of sizes given by indices of $a$:
\begin{align*}
    B_{{6,3}}(a_{1},a_{2},a_{3},a_{4},a_5,a_6)&=\underbrace{15}_{=\binom{6}{4}}a_{1}^2a_{4}+\underbrace{60}_{=\binom{6}{3}\binom{3}{2}}a_{1}a_{2}a_{3}+\underbrace{15}_{=\frac{1}{3!}\binom{6}{2}\binom{4}{2}}a_{2}^{3},\\
B_{{6,4}}(a_{1},a_{2},a_{3},a_{4},a_5,a_6)&=\underbrace{20}_{=\binom{6}{3}} a_1^3 a_3+\underbrace{45}_{=\frac{1}{2!}\binom{6}{2}\binom{4}{2}=} a_1^2 a_2^2.
\end{align*}
\end{example}
In the following section we need the formula for the derivation of Bell polynomials
\begin{equation}\label{eq:Bellder}
    \frac{\partial B_{n,\,k}(a_1,\,a_2,\,\dots)}{\partial a_\ell}=\binom{n}{\ell}B_{n-\ell,\,k-1}(a_1,\,a_2,\,\dots),
\end{equation}
which we prove in \ref{sec:Bellder}. Let us also recall another important property of Bell Polynomials after \cite{Comtet1974}
\begin{equation}\label{eq:Bell-relation}
B_{n,k}(abx_1, ab^2x_2, ...)=a^kb^nB_{n,k}(x_1, x_2 ,...).
\end{equation}

\subsection{Ordinary Bell polynomials}
The Bell polynomials $B_{n,k}$ discussed in the previous sections are sometimes called exponential polynomials (see \cite{Comtet1974}). Let us consider a minor modification of them $\widehat{B}_{n,k}$, which for the sake of distinction, are called ordinary Bell polynomials (see Eq. [3o] in~\cite{Comtet1974}).
\begin{equation}\label{eq:OBP}
    \widehat{B}_{n,k}(a_{1},a_{2},\dots )={\frac {k!}{n!}}B_{n,k}(1!a_{1},2!a_{2},\dots)=\sum_{\{m_r\}_{n,k}} \frac{k!}{m_{1}!m_{2}!\cdots m_{n-k+1}!}\prod _{j=1}^{n}a_j^{m_{j}},
\end{equation}
As we will see in the following sections, ordinary polynomials are more natural for describing the Darwin-Fowler procedure, rather than exponential ones.

\subsection{Bell transformation and combinatorial numbers}
Bell polynomials also allow one to conveniently express known combinatorial numbers in a compact form. There are several versions of the so-called Bell transform (see \cite{Birmayer2019,Lushny}), but we will focus on expressing Lah and Stirling numbers by exponential Bell polynomials (see \cite{Comtet1974,Charalambides2002}). As we will see in the following sections, these numbers are closely related to the distributions in the Darwin-Fowler approach.

\begin{itemize}
    \item Unsigned Lah numbers (see A105278 in \cite{oeis} and Eq. [3h] in \cite{Comtet1974})
\begin{equation}\label{eq:Lah}
    L(n,k)=B_{n,k}(1!,2!,\dots)=\binom{n-1}{k-1}\frac{n!}{k!}.
\end{equation}
\item Unsigned Stirling numbers of the first kind (see A008275 in \cite{oeis} and Eq. [3i] in \cite{Comtet1974})
\begin{equation}\label{eq:S1}
    |S_1(n,k)|=B_{n,k}(0!,1!,2!,\dots).
\end{equation}
\item Stirling numbers of the second kind (see A008277 in \cite{oeis} and Eq. [3g] in  \cite{Comtet1974})
\begin{equation}\label{eq:S2}
    S_2(n,k)=B_{n,k}(1,1,1,\dots).
\end{equation}

\end{itemize}

\section{Bell Polynomial approach to  Darwin-Fowler model}
Using Bell polynomials (both exponential and ordinary) enables to  transform the expression for $\Gamma_{\mathcal{N},\,\mathcal{K}}$ given by Eq. (\ref{eq:Gamma}) as follows
\begin{equation}\label{eq:GammaB}
	\Gamma(\omega_1,\,\omega_2,\,\dots)=\frac{\mathcal{K}!}{\mathcal{N}!}B_{\mathcal{N},\,\mathcal{K}}(1!\omega_1,\,2!\omega_2,\,\dots)=\widehat{B}_{\mathcal{N},\,\mathcal{K}}(\omega_1,\,\omega_2,\,\dots).
\end{equation}
Thus, the expression for a $n^\star_{\mathcal{N},\,\mathcal{K}}$ (see Eqs. (\ref{eq:wartoscocz}) and (\ref{eq:OBP})) which we are looking for can be further simplified
\begin{equation}\label{n_sr_Belle}
n^\star_{\mathcal{N},\,\mathcal{K}}(\ell)=\frac{\mathcal{N}!}{(\mathcal{N}-\ell)!} \frac{\omega_\ell B_{\mathcal{N}-\ell, \mathcal{K}-1}(1!\omega_1, 2!\omega_2, \cdots)}{B_{\mathcal{N},\mathcal{K}}(1!\omega_1, 2!\omega_2, \cdots)}=\mathcal{K}\underbrace{\frac{\omega_\ell \widehat{B}_{\mathcal{N}-\ell, \mathcal{K}-1}(\omega_1, \omega_2, \cdots)}{\widehat{B}_{\mathcal{N},\mathcal{K}}(\omega_1, \omega_2,\cdots)}}_{=\mathcal{P}_{\mathcal{N},\,\mathcal{K}}(\ell)},
\end{equation}
which is the exact result for finite $\mathcal{K}$ and $\mathcal{N}$. As illustrations of this main result let us consider specific forms of degeneration in the following sections.

\section{Special cases of the degeneracy}

\subsection{Constant degeneracy}
Firstly let us assume that there is no degeneracy, i.e. $ \omega_\ell=\omega=\mathrm{constans}$, which changes Eq. (\ref{n_sr_Belle}) to the following form (see Eq. (\ref{eq:Bell-relation}))
\begin{equation*}
	n^\star_{\mathcal{N},\mathcal{K}}(\ell)=\frac{\mathcal{N}!}{\left(\mathcal{N}-\ell\right)!}\frac{B_{\mathcal{N}-\ell,\,\mathcal{K}-1}(1!,\,2!,\,3!,\,\dots)}{B_{\mathcal{N},\,\mathcal{K}}(1!,\,2!,\,3!,\,\dots)}=\frac{\mathcal{N}!}{\left(\mathcal{N}-\ell\right)!}\frac{L(\mathcal{N}-\ell,\,\mathcal{K}-1)}{L(\mathcal{N},\,\mathcal{K})},
\end{equation*}
where we can spot  Lah numbers  (see Eq. (\ref{eq:Lah})), which can be further simplified
\begin{align}\nonumber
	n^\star_{\mathcal{N},\mathcal{K}}(\ell)&=\frac{\mathcal{N}!}{\left(\mathcal{N}-\ell\right)!}\frac{\binom{\mathcal{N}-\ell-1}{\mathcal{K}-2}\frac{\left(\mathcal{N}-\ell\right)!}{\left(\mathcal{K}-1\right)!}}{\binom{\mathcal{N}-1}{\mathcal{K}-1}\frac{\mathcal{N}!}{\mathcal{K}!}}=\mathcal{K}\frac{\binom{\mathcal{N}-\ell-1}{\mathcal{K}-2}}{\binom{\mathcal{N}-1}{\mathcal{K}-1}}=\\&=\mathcal{K}\underbrace{\left(\mathcal{K}-1\right)\frac{\left(\mathcal{N}-\ell-1\right)!\left(\mathcal{N}-\mathcal{K}\right)!}{\left(\mathcal{N}-1\right)!\left(\mathcal{N}-\mathcal{K}-\ell+1\right)!}}_{\mathcal{P}_{\mathcal{N},\mathcal{K}}(\ell)}.\label{eq:kanonicznyLah}
\end{align}
Typically we would look for the asymptotic form of Eq. (\ref{eq:kanonicznyLah}) using the steepest descent method, which is described in \ref{sec:steepestdescend}. However, we will use the compact form of Lah numbers and apply the Stirling's approximation in  $\spadesuit$ and obtain
\begin{align*}
\mathcal{P}_{\mathcal{N},\mathcal{K}}=&\left(\mathcal{K}-1\right)\frac{\left(\mathcal{N}-\ell-1\right)!}{\left(\mathcal{N}-1\right)!}\frac{\left(\mathcal{N}-\mathcal{K}\right)!}{\left(\mathcal{N}-\mathcal{K}-\ell+1\right)!}\stackrel{\spadesuit}{\approx}\\
\stackrel{\spadesuit}{\approx}& \frac{\left(\mathcal{K}-1\right)\left(\mathcal{N}-\ell-1\right)^{\mathcal{N}-\ell-1}e^{-\left(\mathcal{N}-\ell-1\right)}}{\left(\mathcal{N}-1\right)^{\mathcal{N}-1}e^{-\left(\mathcal{N}-1\right)}}\frac{\left(\mathcal{N}-\mathcal{K}\right)^{\mathcal{N}-\mathcal{K}}e^{-\left(\mathcal{N}-\mathcal{K}\right)}}{\left(\mathcal{N}-\mathcal{K}-\ell+1\right)^{\mathcal{N}-\mathcal{K}-\ell+1}e^{-\left(\mathcal{N}-\mathcal{K}-\ell+1\right)}}=\\
=& \frac{\left(\mathcal{K}-1\right) e\left(\mathcal{N}-\mathcal{K}\right)^{\ell-1}}{\left(\mathcal{N}-\ell-1\right)^{\ell}}\left(1-\frac{\ell}{\mathcal{N}-1}\right)^{\mathcal{N}-1} \left(1+\frac{\ell-1}{\mathcal{N}-\mathcal{K}-\ell+1}\right)^{\mathcal{N}-\mathcal{K}-\ell+1} =\\
=& e\frac{\mathcal{K}-1}{\mathcal{N}-\mathcal{K}}\left( \frac{\mathcal{N}-\mathcal{K}}{\mathcal{N}-\ell-1}\right)^\ell \left(1-\frac{\ell}{\mathcal{N}-1}\right)^{\mathcal{N}-1} \left(1+\frac{\ell-1}{\mathcal{N}-\mathcal{K}-\ell+1}\right)^{\mathcal{N}-\mathcal{K}-\ell+1}.
\end{align*}
Reasonably assuming that the average energy $U$ of the system is constant even in the limit of large $\mathcal{K}$ i.e.
$\mathcal{N}=U \mathcal{K}$ one gets final form of the energy distribution as
\begin{align*}
\mathcal{P}_{\mathcal{N},\,\mathcal{K}}(\ell)\approx& e \underbrace{\frac{\mathcal{K}-1}{(U-1)\mathcal{K}}}_{\xrightarrow{\mathcal{K}\rightarrow \infty} (U-1)^{-1}}\underbrace{\left( \frac{(U-1)\mathcal{K}}{U\mathcal{K}-\ell-1}\right)^\ell}_{\xrightarrow{\mathcal{K}\rightarrow \infty} \left(\frac{U-1}{U}\right)^\ell}\underbrace{\left(1-\frac{\ell}{U\mathcal{K}-1}\right)^{U\mathcal{K}-1}}_{\xrightarrow{\mathcal{K}\rightarrow \infty} e^{-\ell}}\times\\ & \times\underbrace{\left(1+\frac{\ell-1}{(U-1)\mathcal{K}-\ell+1}\right)^{(U-1)\mathcal{K}-\ell+1}}_{\xrightarrow{\mathcal{K}\rightarrow \infty} e^{\ell-1}},
\end{align*}
which results in the expected formula for the distribution
\begin{equation}\label{eq:Pkanon}
\mathcal{P}_\infty(\ell)=\frac{e^{\beta \varepsilon \ell}}{Z},
\end{equation}
where $U=(1-e^{\beta \varepsilon})^{-1}$ and $Z=U-1$.

\subsection{Harmonic oscillators}
For the one-dimensional harmonic oscillator the weight of energy $\varepsilon_\ell$ is equal to
\begin{equation*}
\omega_\ell=(\ell+1),
\end{equation*}
which implies the following form of the energy distribution from Eq. (\ref{n_sr_Belle})
\begin{equation}\label{eq:HO1}
\mathcal{P}_{\mathcal{N},\,\mathcal{K}}(\ell)=(\ell+1)\frac{ \widehat{B}_{\mathcal{N}-\ell, \mathcal{K}-1}(2, 3, \dots)}{\widehat{B}_{\mathcal{N},\mathcal{K}}(2, 3,\dots)}.
\end{equation}
Eq. (\ref{eq:HO1}) can be further generalized for the $D$-dimensional harmonic oscillator, because degenerations are then given as
\begin{equation*}
\omega_\ell=(\ell+1)(\ell+2)\dots(\ell+D)=:(\ell+1)^{(D)},
\end{equation*}
where $(\ell+1)^{(D)}$ denotes the rising factorial for brevity. With such degeneracies one face the problem of  determination the values of  Bell polynomials $B_{n,k}(2^{(D)},3^{(D)},\dots)$, which can be done with the following formula
\begin{equation}\label{Bell_Stirling}
B_{n,k}(2^{(D)},3^{(D)},\dots)=\left[ (D-1)!\right]^k \sum_{r=k}^n |S_1(n,r)|S_2(r,k)D^r,
\end{equation}
where $S_1, S_2$ are the Stirling numbers of the first and second kind respectivelly, see Eqs. (\ref{eq:S1}) and (\ref{eq:S2}). Let us note that  Eq. (\ref{Bell_Stirling}) follows from Eq. (8.50) in \cite{Charalambides2002} and allows one to obtain the following final formula for the most probable configuration for $D$-dimensional harmonic oscillator
\begin{equation}\label{nsroscylatory}
n^\star_{\mathcal{N},\mathcal{K}}=\frac{\mathcal{N}!}{(\mathcal{N}-\ell)!}\frac{(\ell+1)^{(D)}}{(D-1)!}\frac{\sum_{r=\mathcal{K}-1}^{\mathcal{N}-\ell}{|S_1(\mathcal{N}-\ell,r)|S_2(r,\mathcal{K}-1)D^r}}{\sum_{r=\mathcal{K}}^\mathcal{N} |S_1(\mathcal{N}, r)| S_2(r,\mathcal{K})D^r}.
\end{equation}

As we can see, the above Eq. (\ref{nsroscylatory}) combines the energy distribution for a harmonic oscillator with the Stirling numbers.

\section{Acknowledgement}
We would like to thank Agata Fronczak for stimulating discussions (and pointing out our mistakes at the early stage of the work). GS work has been supported by the National Science Centre of Poland (Narodowe Centrum Nauki, NCN) under grant no. 2015/18/E/ST2/00560.