# Explosive Ejections Generated by Gravitational Interactions

P. R. Rivera-Ortiz, A. Rodríguez-González, J. Cantó, Luis A. Zapata

► To cite this version:

P. R. Rivera-Ortiz, A. Rodríguez-González, J. Cantó, Luis A. Zapata. Explosive Ejections Generated by Gravitational Interactions. The Astrophysical Journal, 2021, 916, 10.3847/1538-4357/ac05bb. insu-03705295

HAL Id: insu-03705295

https://insu.hal.science/insu-03705295v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Explosive Ejections Generated by Gravitational Interactions

P. R. Rivera-Ortiz¹ ID, A. Rodríguez-González² ID, J. Cantó³, and Luis A. Zapata⁴ ID

¹ Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG), F-38000 Grenoble, France

² Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Ap. 70-543, 04510 D.F., Mexico

³ Instituto de Astronomía, Universidad Nacional Autónoma de México, Ap. 70-264, 04510 D.F., Mexico

⁴ Instituto de Radioastronomía y Astronomía, Universidad Nacional Autónoma de México, Apdo. Postal 3-72 (Xangari), 58089 Morelia, Michoacán, Mexico

Received 2020 September 24; revised 2021 May 12; accepted 2021 May 25; published 2021 July 27

## Abstract

During the fragmentation and collapse of a molecular cloud, it is expected that it will have close encounters with (proto)stellar objects that can lead to the ejection of a fraction of them as runaway objects. However, the duration and consequences of such encounters are perhaps small enough for there to be no direct evidence of their occurrence. As a first approximation, herein we analytically study the interaction of a massive object that moves at high velocity into a cluster of negligible-mass particles with an initial number density distribution $\propto R^{-\alpha}$. We have found that the runaway conditions of the distribution after the encounter are related to the mass and the velocity of the star and the impact parameter of each particle to the stellar object. Then, the cluster particles are gravitationally accelerated by the external approaching star, destroying the cluster, and the dispersion and velocities of the particles have explosive characteristics. We compare this analytical model with several numerical simulations and, finally, apply our results to the Orion fingers in the Orion BN/KL region, which show an explosive outflow that could be triggered by the gravitational interaction of several (proto)stellar objects.

*Unified Astronomy Thesaurus concepts: Interstellar dynamics (839); Interstellar medium (847); Star forming regions (1565)*

## 1. Introduction

During recent decades the star formation processes had been extensively studied, which has led to a general comprehension of the stages that drive a molecular cloud to collapse, fragment, and form new stars (McKee & Ostriker 2007; Krumholz 2014), arriving at the formation of stellar clusters in a star-forming region where the star-forming efficiency and the conditions to create and preserve a gravitational bound have to be considered (Kirk et al. 2014). This problem has been analyzed by Fall et al. (2010) and Kruijssen et al. (2012). In any case, and even when a forming cluster is bounded, the local conditions could eject some stars as a result of close gravitational interactions from a multiple-star system (Irrgang et al. 2018). These encounters are expected, at least in the densely populated regions, but there is not any observational direct evidence of them. Therefore, the duration and characteristics of such encounters, and their influence in the interstellar medium, are still unknown.

It has been suggested recently that the explosive outflows reviewed by Bally et al. (2017) may have produced as a consequence of these kind of close encounters. The closest of these explosive outflows is Orion BN/KL, which may have been produced by a close gravitational interaction of several protostellar objects (Becklin & Neugebauer 1967; Kleinmann & Low 1967). This outflow shows a characteristic filamentary structure emitting in H₂ known as the *Orion fingers* with a kinetic energy of around $10^{47}-10^{48}$ erg distributed almost isotropically around a common center, which is the same common origin for the runaway stellar objects BN, x, and the binary I, with masses around $10 M_{\odot}$ (Bally et al. 2020). However, motivated by the Orion BN/KL morphology, Allen & Burton (1993) proposed that this outflow could have an explosive origin where some dense clumps were ejected at very high velocities and interact with the surrounding environment, with the resulting wakes being the actual fingers. At the tip of these fingers there is highly excited emission, e.g., Fe II, produced by the high velocities of these bullets, which has a

proper motion as large as 300 km s⁻¹ (Allen & Burton 1993; Cunningham 2006; Nissen et al. 2007). Zapata et al. (2009) reported a set of streamers emitting in CO $J=2 \to 1$, which are related to the fingers, and that follow a Hubble law that is a signature of other explosive outflows like DR21 (Zapata et al. 2013) and G5.89 (see Zapata et al. 2019, 2020).

The problem of a star ejected from a stellar cluster or the evaporation of members from a stellar cluster has been previously investigated by several authors, including Blaauw (1961), King (1966), Perets & Šubr (2012), and Wang et al. (2019). They analyzed the close encounters of cluster members, which lead to the formation of close binaries and another star taking the excess energy as kinetic energy. Nevertheless, this mechanism does not take into account the formation of explosive outflows (reviewed by Bally et al. 2017). It has been suggested that a close encounter of an already runaway star disrupts a stellar cluster. In particular, Bally et al. (2015) proposed that an external runaway star disrupted a forming stellar cluster, and the debris of this disruption is the explosive outflow from Orion BN/KL.

In order to explain the explosive outflow of Orion BN/KL, Bally et al. (2015) proposed a qualitative model to describe the close interaction of a forming stellar cluster that is disturbed by an external high-velocity and massive object destroying the former cluster through unstable orbits, then producing smaller bounded systems and releasing the excess energy in the gaseous environment through several close interactions of the involved objects expelling gas clumps that could form the filamentary structure of the explosive outflow. Evidence exists in the literature of high-velocity (up to several hundreds of kilometers per second) stellar objects, such as runaway neutron stars (Irrgang et al. 2018). Dorigo Jones et al. (2020) analyzed a sample of runaway stars with masses from 10–60 $M_{\odot}$ and velocities from 20–200 km s⁻¹, which makes it plausible that a runaway star could impact and disrupt a cluster. Also, Bally et al. (2020) performed and reviewed
several hydrodynamical simulations that investigate the probability of a system of four stellar objects forming a system of runaway stars such as the one observed in Orion BN/KL and there is evidence of runaway stars moving away from a common point (Rodríguez et al. 2017). Nevertheless, the mechanism that triggered explosive outflows has not been deeply explored, since it cannot explain the acceleration and dispersion of the clumps that create the fingers compared with the observed velocities. Also, the repeated encounters seem to be against the idea of a single explosive event.

As Bally et al. (2015) proposed, the Orion BN/KL could have suffered an interaction between a forming stellar cluster and an external stellar object, disrupting the cluster members, but the origin of the high-velocity features known as the H₂ Orion fingers and the CO streamers, which have velocities the order of several hundreds of kilometers per second, is still unknown. Explaining the ejection of the residual clumps from a disrupting cluster can then lead to a better understanding of the explosive outflows. The nature of this cluster is actually unknown and its geometry, size, or mass are a matter of debate. Nissen et al. (2007) found an inferior limit for the leading clumps of the fingers of $10^{-5} M_{\odot}$ and according to Rivera-Ortiz et al. (2019a, 2019b; hereafter RO19a and RO19b, respectively) and Dempsey et al. (2020), the initial mass of each clump is in the planetary mass scale, so it could be reasonable to think that they come from the young protoplanetary disk of one of the stellar objects. Also, according to Morbidelli & Nesvorny (2012) and Birnstiel et al. (2012), the number, mass, and distribution of the planetary embryos in the protoplanetary disks are related to the time evolution of a star. These clumps would suffer a similar interaction as the one described in Bally et al. (2015), such as a Rutherford dispersion originated by a gravitational force, affecting them drastically and changing their dynamical properties via one gravitational impulse (Sastry & Alladin 1970) produced by the passage of a massive object. As a first approximation, we describe the interaction of a fast massive star with a cluster of negligible-mass objects.

In this paper we analyze the gravitational effect of a massive star on a cloud of smaller objects as the possible origin of the explosive outflows described before. In Section 2 we propose the analytical tools to describe the interaction, and in Section 3 we describe the N-body simulations used to calibrate the analytical model. We propose some interpretations applied to Orion BN/KL in Section 4 and finally we conclude in Section 5.

## 2. Model

In this paper we address the problem of the gravitational collision between a cluster of cloudlets of negligible mass and a massive object with mass M, moving at a high velocity $v_0$.

### 2.1. Static Massive Particle

We first consider the motion of one of these cloudlets in the frame of reference of the colliding object (the system xy, see Figure 1). In this frame of reference, the cloudlet approaches the massive object with a velocity $v_0$ (at infinity), along an initial straight trajectory with an impact parameter $\xi$. Considering that initial gravitational interaction is small enough to be neglected, the mechanical energy of the incident particle is positive (equal to $v_0^2/2$ per unit mass) and the trajectory is a hyperbole.

<figure><img src="image_3.png" /><figcaption>Figure 1. Gravitational interaction between a mass M moving with velocity $v_0$ and a massless cloudlet at rest. (a) In the frame of reference of the mass M, the cloudlet initially moves with velocity $v_0$ at infinity, along a straight trajectory with an impact parameter $\xi$. After the collision, the cloudlet deviates from its initial trajectory by an angle $\alpha_m$ and regains a velocity $v_0$ equal to its initial velocity. (b) In the frame of reference of the cluster, the cloudlet is initially at rest while the star is moving with velocity $v_0$. The cloudlet is accelerated and is thrown away at an angle $\alpha'_m$ with velocity $v'_m$ (see text).</figcaption></figure>

We define the gravitational radius $\xi_0$ as

$$ \xi_0 = \frac{GM}{v_0^2}, \qquad (1) $$

and the characteristic time $t_0$,

$$ t_0 = \frac{\xi_0}{v_0}. \qquad (2) $$

Using these quantities, we adimensionalize distances with $\xi_0$, velocities with $v_0$, and times with $t_0$.

Then, the hyperbolic trajectory of a cloudlet, in polar coordinates, is given by (Calkin 1989; see also Cantó et al. 2011),

$$ r = \frac{\xi^{\prime 2}}{1 - \cos \theta + \xi' \sin \theta}, \qquad (3) $$

where $r$ and $\xi'$ are the dimensionless distance to the particle and impact parameter, respectively.

In Equation (3), the angle $\theta$ is measured with respect to the $+\hat{x}$ direction (see Figure 1). From this equation, we can observe that for $r \to \infty$, the cloudlet moves along incoming and outcoming asymptotic straight lines, which are given by the angles $\theta$ equal to 0 and $\theta_\infty$. This last sentence can be written as

$$ \cos \theta_{\infty} = \frac{1 - \xi'^{2}}{1 + \xi'^{2}}. \qquad (4) $$

We set the scattering angle, i.e., the angle between the incoming and outcoming directions of the movement, as $\alpha_m = \theta_\infty - \pi$ (see
Figure 1). From Equation (4), we obtain:

$$ \cos \alpha_m = \frac{\xi'^2 - 1}{\xi'^2 + 1}. \qquad (5) $$

Our next step is to obtain the x- and y-components of the cloudlet's velocity as a function of $\theta$ and $\xi'$. The dimensionless velocity can be written in polar coordinates as

$$ \mathbf{u} = u_r \hat{r} + u_\theta \hat{\theta} = \frac{dr}{d\tau} \hat{r} + r \frac{d\theta}{d\tau} \hat{\theta} = \frac{d\theta}{d\tau} \left( \frac{dr}{d\theta} \hat{r} + r \hat{\theta} \right), \quad (6) $$

where $\tau$ is the dimensionless time. In this problem, the angular momentum is conserved, which is given by

$$ \xi' = ru_{\theta} = r^2 \frac{d\theta}{d\tau}. \qquad (7) $$

Then, we finally obtain $u_x$ and $u_y$, putting Equations (3) and (7) in Equation (6), and writing the directions $\hat{r}$ and $\hat{\theta}$ as functions of $\hat{x}$, $\hat{y}$, and $\theta$:

$$ \mathbf{u} = u_x \hat{x} + u_y \hat{y} = -\frac{1}{\xi'} (\xi' + \sin \theta) \hat{x} + \frac{1}{\xi'} (1 - \cos \theta) \hat{y}. \quad (8) $$

Equation (8) can be used to obtain the velocity module $u = \sqrt{u_x^2 + u_y^2}$ and, after some algebra, we obtain the mechanical energy conservation equation,

$$ \frac{1}{2} u^2 - \frac{1}{r} = \frac{1}{2}. \qquad (9) $$

Besides, we find a relationship between $d\tau$ and $d\theta$, combining Equations (3) and (7), which is given by

$$ d\tau = \frac{\xi'^3}{(1 - \cos \theta + \xi' \sin \theta)^2} d\theta, \qquad (10) $$

which can be integrated to give

$$ \tau + c = \ln \left[ 1 + \xi' \cot \frac{\theta}{2} \right] + \frac{1 + \xi'^2}{2 \left[ 1 + \xi' \cot \frac{\theta}{2} \right]} - \frac{\xi'}{2} \cot \frac{\theta}{2}, \qquad (11) $$

with $c$ as an integration constant. By choosing $\tau=0$ for an angle $\theta=\pi$ in the limit $\xi'=0$ we find that $c = \frac{1}{2}$.

Defining $\mu = [1 + \xi' \cot \frac{\theta}{2}]$, the relation between the position angle, e.g., the angle $\theta$, and the time $\tau$ is given by

$$ \tau = \ln[\mu] + \frac{1 + \xi'^2}{2\mu} - \frac{\mu}{2}. \qquad (12) $$

For an angle $\theta = \pi$, according to Equations (3) and (12), the particle crosses the x-axis ($\theta = \pi$) at a time and position that depend on the impact parameter, such that

$$ \tau_* = \frac{\xi'^2}{2} \qquad (13) $$

$$ x_* = -\frac{\xi'^2}{2}. \qquad (14) $$

Note that, correctly, for $\xi=0$, $\tau_* = x_* = 0$.

## 2.2. Moving Massive Particle

In this subsection, we adopt the frame of reference where both the observer and the cloudlets are at rest with respect to each other, and the massive particle moves with velocity $v_0$

(see Figure 1). In this frame of reference, the cloudlet has coordinates

$$ x' = x + \tau, \quad y' = y, \qquad (15) $$

and velocities

$$ u_x' = u_x + 1, \quad u_y' = u_y. \qquad (16) $$

From this last equation and considering Equation (8) we find that the velocity is given by

$$ u' = \frac{[2(1 - \cos \theta)]^{1/2}}{\xi'}. \qquad (17) $$

Then, the direction in which the cloudlet moves after interacting with the massive particle and far away from it, i.e., for $\theta = \theta_\infty$ (see Equation (4)), forms an angle $\pi - \alpha'_m$ with respect to $\hat{x}$ (see Figure 1(b)). This angle is obtained by doing

$$ \cos(\pi - \alpha'_m) = -\cos \alpha'_m = \frac{u_x'}{u'} = \frac{-\sin \theta_\infty}{\sqrt{2(1 - \cos \theta_\infty)}}. \qquad (18) $$

Considering Equation (4) into Equation (18), we get

$$ \alpha'_m = \arccos \left( - \frac{1}{(1 + \xi'^2)^{1/2}} \right). \qquad (19) $$

Finally, the position $x'_*$ when the particle crosses the x-axis is

$$ x'_* = x_* + \tau_* = 0, \qquad (20) $$

which is independent of $\xi$ (see Equations (14) and (13)), and therefore every particle crosses the x-axis at the same point although at a different time.

Combining Equations (17) and (19) we obtain:

$$ u = \frac{2}{(1 + \xi'^2)^{1/2}}. \qquad (21) $$

Using the dimensional expressions we have:

$$ v = \frac{2v_0}{(1 + \xi'^2)^{1/2}}. \qquad (22) $$

## 2.3. The Velocity Distribution

Consider a cluster of radius $\xi_c$ containing $N_T$ clumps. The density distribution is

$$ n(r) = Ar^\alpha. \qquad (23) $$

Then,

$$ N_T = \int_0^{\xi_c} 4\pi r^2 n(r) dr = \frac{4\pi A}{3 + \alpha} \xi_c^{3+\alpha} \qquad (24) $$

for $\alpha \neq -3$ and thus,

$$ A = \frac{(3 + \alpha)N_T}{4\pi\xi_c^{3+\alpha}}. \qquad (25) $$

Let $dN(\xi)$ be the number of clumps with an impact parameter between $\xi$ and $\xi+d\xi$. We define the distribution $g(\xi)$ of impact parameters such that

$$ dN(\xi) = N_T g(\xi) d\xi. \qquad (26) $$

Clearly,

$$ \int_0^{\xi_c} g(\xi) d\xi = 1 \qquad (27) $$

where $\xi_c$ is the maximum impact parameter.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Diagram of the cloudlet cluster with size $\xi_c$. The x-axis is parallel to the direction of motion of the massive particle at infinity and, for a given impact parameter $\xi$, there is a distribution function $g(\xi)$ that accounts for the fraction of cloudlets that are between $\xi$ and $\xi + d\xi$ in the limits $-x_m$ and $x_m$.</figcaption></figure>

From Figure 2 we find

$$dN(\xi) = 2\pi\xi d\xi \left( 2 \int_{0}^{x_m} n(x) dx \right) \quad (28)$$

where $x_m = (\xi_c^2 - \xi^2)^{1/2}$ and $n(x) = Ar^\alpha = A(\xi^2 + x^2)^{\alpha/2}$.

Let

$$\int_{0}^{x_m} n \, dx = A \int_{0}^{x_m} (\xi^2 + x^2)^{\alpha/2} \, dx = AI_{\alpha} \quad (29)$$

where

$$I_{\alpha} = \int_{0}^{x_m} (\xi^2 + x^2)^{\alpha/2} dx. \quad (30)$$

Then, from Equations (28) and (30), we can express $g(\xi)$ as

$$g(\xi) = 4\pi A \xi I_{\alpha} / N_T = \frac{(3 + \alpha)\xi}{\xi^{3+\alpha}} I_{\alpha}, \quad (31)$$

where we have used Equation (26).

Let $f(v)$ be the velocity distribution function, $N_T f(v) dv = N_T g(\xi) d\xi$, where the left side is the number of clumps with velocities between $v$ and $v + dv$. Then,

$$f(v) = g(\xi) \left| \frac{d\xi}{dv} \right|. \quad (32)$$

Clearly,

$$\int_{v_{\min}}^{v_{\max}} f(v) dv = 1 \quad (33)$$

From Equation (22) we obtain

$$\left| \frac{dv}{d\xi} \right| = \frac{2\xi\xi_0 v_0}{(\xi^2 + \xi_0^2)^{1/2}} = \frac{4v_0^2 \xi_0^2}{\xi v^3} \quad (34)$$

to finally express $f(v)$ as

$$f(v) = \frac{4(3 + \alpha) v_0^2 \xi_0^2}{\xi_c^{3+\alpha} v^3} I_{\alpha} \quad (35)$$

where we have used Equation (33).

According to Equation (22), the distribution function for the velocity has a lower limit:

$$v_{\min} = \frac{2v_0}{(1 + \xi_c^{\prime 2})^{1/2}}, \quad (36)$$

and an upper limit:

$$v_{\max} = 2v_0. \quad (37)$$

Additionally, the velocity distribution given by Equation (35) has a maximum and/or a minimum given by the condition

$$\left[ \frac{dI_{\alpha}}{dv} \right]_{v_{*}} = \frac{3I_{a}(v_{*})}{v_{*}} \quad (38)$$

for any of maximum or minimum at $v_*$.

We consider the following particular cases:

$$\alpha = 0, \quad I = x_m \quad (39)$$

$$\alpha = -1, \quad I = \arcsin h(x_m / \xi) \quad (40)$$

$$\alpha = -2, \quad I = \frac{1}{\xi} \arctan(x_m / \xi). \quad (41)$$

Then, a cluster with uniform density, i.e., $\alpha = 0$, implies three different velocities that characterize the velocity distributions $v_{\min}$, $v_{\max}$, and $v_{*,max}$ where

$$v_{*,max} = \frac{4v_0}{\sqrt{3(1 + \xi_c^{\prime 2})}}. \quad (42)$$

So,

$$\frac{v_{*,max}}{v_{\min}} = \frac{2}{\sqrt{3}} \approx 1.1547, \quad (43)$$

$$f_* = \frac{9\sqrt{3}(1 + \xi_c^{\prime 2})}{32\xi_c^{\prime 3}v_0}, \quad (44)$$

and

$$f(v_{\max}) = \frac{3}{2v_0 \xi_c^{\prime 2}}.$$

This velocity distribution is almost flattened in magnitude. Then, to explore the spatial distribution, numerical simulations are performed and discussed in the next section.

# 3. N-body Simulation Results

## 3.1. Initial Conditions

In order to test the analytical study, we have performed several N-body simulations of a particle of $M_* = 10 M_\odot$ “colliding” with a cluster of $N_T$ small mass particles to reproduce the velocity distributions derived in the last section. The free parameters of the simulations are $\alpha$ and $v_0$, and the impact parameter of the massive particle respect to the cluster center is $y^*$.

We selected $\xi_c = 0.67$ au as the cluster radius, which is a small region considering the size of a protostellar envelope. However, a larger size would represent a weaker interaction where the internal forces should be considered. The number of total particles is $N_T = 200$, which, approximately, have a typical separation of around 0.11 au. The cluster is centered at the origin in the simulations and we consider a spatial
<table><caption>**Table 1**
Numerical Simulation Models of a 10 $M_{\odot}$ Particle Moving at $v_0$ toward a Cluster with a Number Density Parameter $\alpha$ from the Point ($-10$ au, $y^*, 0$)</caption><thead><tr><th>Models</th><th>$v_0$ (km s<sup>-1</sup>)</th><th>$y_*$</th><th>$\alpha$</th></tr></thead><tbody><tr><td>v200R0</td><td>200</td><td>0</td><td>0</td></tr><tr><td>v200R1</td><td>200</td><td>0</td><td>-1</td></tr><tr><td>v200R2</td><td>200</td><td>0</td><td>-2</td></tr><tr><td>v100R1</td><td>100</td><td>0</td><td>-1</td></tr><tr><td>v300R1</td><td>300</td><td>0</td><td>-1</td></tr><tr><td>v400R1</td><td>400</td><td>0</td><td>-1</td></tr><tr><td>v200R1s05</td><td>200</td><td><sup>1</sup>&frasl;<sub>2</sub> $\xi_c$</td><td>-1</td></tr><tr><td>v200R1s1</td><td>200</td><td>$\xi_c$</td><td>-1</td></tr><tr><td>v200R0m</td><td>200</td><td>0</td><td>0</td></tr></tbody></table>

distribution (the number of clumps per unit volume) of the form

$$n(r) = \frac{(3 + \alpha)N_T}{4\pi\xi_c^{3+\alpha}} r^{\alpha}, \qquad (45)$$

where we analyze the case $\alpha = 0$, and for an homogeneous distribution $\alpha = -1$ and $\alpha = -2$.

This cluster must be in dynamical equilibrium, then, assuming that every particle has a circular orbit around the distribution center of mass, the orbital velocity of a particle at a radius $r$ is

$$v(r) = \sqrt{\frac{N_T M_i G}{\xi_c}} \left(\frac{r}{\xi_c}\right)^{\frac{\alpha+2}{2}}, \qquad (46)$$

where $M_i$ is their individual mass and $G$ is the gravitational constant. The analytical model considers the cluster particles at rest, so we choose $M_i = 10^{-10} M_{\odot}$ to have a very small orbital velocity $\sim 10^{-2}$ km s$^{-1}$ that, effectively, allows us to use zero velocity for every cluster particle. Also, the force between two cluster particles is negligible. We use a random number generator that chooses a number $\eta$ uniformly distributed in the interval [0,1]. The value is related to the radial distance as

$$r = \xi_c \eta^{1/(3+\alpha)}, \qquad (47)$$

from which we can sample $\xi$ as a function of the random number $\eta$ (for more detail see Rodríguez-González et al. 2007). We assigned random directions to the position vector of each particle.

We also included a single massive particle with a mass of $10 M_{\odot}$, moving toward the particle distributions with $v_0 = 100, 200, 300$, and $400$ km s$^{-1}$. At $t = 0$, the massive particle starts moving from the Cartesian point ($-10$ au, $y^*, 0$), in a direction parallel to the x-axis toward the clump distribution. Table 1 shows the initial velocity $v_0$ and the initial position $y^*$ over the y-axis of the massive particle and the exponent $\alpha$ of the distribution of particles at rest for the simulations presented in this paper.

In order to obtain a statistically significant result to compare with the velocity distribution, we have developed sets of 10 random distributions in each of our models. The stability of our numerical solver is presented in the Appendix.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The graph represents the velocity $v_{\infty}$ (which is constant after the interaction) of each particle in the v200R1 model (open circles) related to the distance from the ejection point at an evolutionary time of 500 yr. The solid line is a model of an homologous expansion (Hubble law) at time = 500 yr, and maximum distance of $4 \times 10^4$ au.</figcaption></figure>

## 3.2. Comparison with Analytical Results

First of all, the consideration of negligible-mass particles allowed us to ignore the potential energy of the cluster, since the kinetic energy of the massive particle is greater by several orders of magnitude and in this approximation the massive particle is not affected by the interaction. For the low-mass particles, there is a very fast interaction that ejects them in almost every direction; Figure 3 shows the velocity $v_{\infty}$, constant after the interaction, related to the distance, in astronomical units, at an evolutionary time of 500 yr for the v200R1 model (open circles) and the homologous expansion model (solid line) with final distance of $4 \times 10^4$ au at 500 yr. The average difference between the dynamical age of each point in the simulation and the model is $\sim 0.65$ yr, corresponding to the time when the low-mass particles reach $v_{\infty}$. As we can see, the numerical model follows a Hubble-type law, however, the small difference is due to: (a) the time it takes ($\sim 0.25$ yr) for the massive particle to get to the center of the particle cluster and leave it, considering $v_0 = 200$ km s$^{-1}$, and (b) the time for the low-mass particles and the massive particle to get their terminal velocities ($\sim 0.4$ yr), after their dynamic interaction.

Figure 4's upper panels show the velocity as a function of the impact parameter $\xi$; the solid line is the solution presented in Equation (22) and the open circles are each of the low-mass particles from our numerical models. The lower panels show the histograms obtained by integrating the distribution function (Equation (35)) for bins with a width $\Delta v = 10$ km s$^{-1}$, as the solid line, and the numerical solution is showed in bins. The models v200R0, v200R1, and v200R2, corresponding to models with a initial massive particle velocity of 200 km s$^{-1}$ and particle distribution with $\alpha = 0, -1$, and $-2$ (left, center, and right columns, respectively), are presented. Figure 4 shows a very good
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Analytical and numerical results for the models v200R0, v200R1, and v200R2 in the left, center, and right columns, respectively. The velocity vs. the initial impact parameters are plotted in the upper panels using open circles for the numerical simulation particles and a solid line for the analytical results, and the final velocity distributions are plotted in the lower panels, using histograms for the numerical model results. The analytical solution (Section 2) is represented as a superimposed curve.</figcaption></figure>

agreement between the numerical and analytical results; the minimum velocity is given by Equation (36). For the models with $v_0 = 200 \text{ km s}^{-1}$ we use $\xi_c = 3\xi_0$ (three times the gravitational radius), therefore $v_{\min} = 126.4 \text{ km s}^{-1}$, shown in these three models (v200R0, v200R1 and v200R2) independently of $\alpha$. The maximum velocity is given by Equation (37), in this case it is $400 \text{ km s}^{-1}$ for particles with $\xi = 0$. The velocity of our numerical models tends to this value and it is more evident in the model with $\alpha = -2$, where the particles are more concentrated at the center of the cluster; in this case, the velocity distribution is vertically asymptotic in $v_{\max}$. In the constant-density case, $\alpha = 0$, $v_{*,\max} = 145.47 \text{ km s}^{-1}$ and for this $\alpha$ the distribution does not have a minimum (local or global). For the cases where $\alpha = -1$ there are no analytical expressions for the position of the maximum or minimum in the velocity distribution, but this can be obtained semianalytically, for the case of $\alpha = -1$ $v_{*,\max} = 151.2 \text{ km s}^{-1}$ and $v_{*,\min} = 375.1 \text{ km s}^{-1}$, and for $\alpha = -2$ $v_{*,\max} = 163.5 \text{ km s}^{-1}$ and $v_{*,\min} = 308 \text{ km s}^{-1}$. These analytical values are in good agreement with the numerical ones shown in Figure 4.

Moreover, we ran models with different velocities, fixing $M_*$ and $\alpha = -1$, therefore the gravitational radius is different in each of our models: $\xi_0 = 0.85, 0.22, 0.095$, and $0.053$ au for v100R1, v200R1, v300R1, and v400R1, respectively. These models, with the exception of v200R1, are shown in Figure 5. The description of the rows (and the plots in them) is the same as Figure 4 except that the columns correspond to the models (V100R1, V300R1, and V400R1, are in the left, center, and right columns, respectively), whereas v200R1 is in the second column of Figure 4. Like in Figure 4, the analytical solution is in accordance with the numerical solution. For all these models we ran the numerical simulations using $\xi_c = 0.66$ au,

which is three times the gravitational radius of models with $v_0 = 200 \text{ km s}^{-1}$. As one expects, the particles of the V100R1 model are distributed in a cluster with a radius smaller than its gravitational one. Therefore the velocity distribution has values around the maximum speed, which is $v_{\max} = 200 \text{ km s}^{-1}$ in this model. In our models with higher initial velocities, which means a lower gravitational radius, the cluster is distributed in a larger volume, because we have fixed the radius of the cluster and the initial impact parameters of the particles are larger, so the velocity distribution has a similar form but one distributed in a larger range of velocities (corresponding at a maximum velocity on each model).

### 3.3. Nonsymmetrical Collision

It is also important to analyze the distribution of the tangential velocities when the massive particle does not pass through the center of the particles’ distribution. In the models v200R1s05 and v200R1s1, the massive particle goes through the cluster distribution at half the radius of the cluster on the y-axis and on the edge of the cluster, on the y-axis as well. Figure 6 shows the v200R1s05 model in the upper row and v200R1s1 in the lower row, where the right column shows the velocity as a function of the impact parameter and the left column shows the histogram of the tangential velocities. As one can see, the velocities of the particles of v200R1s1 are minors, because the impact parameters of the particles are larger than for v200R1s05’s particles, or for the case where the massive particle crosses the center of the distribution where the particle density is greater, because it increases toward the center of the cluster. The plots show a maximum plane-of-the-sky velocity of about 100 and $90 \text{ km s}^{-1}$ for v200R1s05 and v200R1s1, respectively (see Table 1).
<figure><img src="image_7.png" /><figcaption>Figure 5. Analytical and numerical results for the models V100R1, V300R1, and V400R1 (left, center, and right columns, respectively). The velocity vs. the initial impact parameters are plotted in the upper panels using open circles for the numerical simulation particles and a solid line for the analytical results, and the final velocity distributions are plotted in the bottom panels, using histograms for the numerical model results. The analytical solution (Section 2) is represented as a superimposed curve.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Velocity distributions of the v200R1s05 and v200R1s1 models, upper and lower panels, respectively. The velocity as a function of the initial impact parameter (left column) and the final velocity distribution (right panels) is shown, with rotation angles $\theta = 30^{\circ}$ and $\phi = 30^{\circ}$.</figcaption></figure>

## 4. The Case of Orion BN/KL

Several examples of explosions generated by dynamical interaction are presented in the astrophysics literature, i.e., globular cluster disintegration and destruction of planetary systems in dynamical encounters (Spurzem et al. 2009) or in an explosive star formation region such as Orion BN/KL. In Orion

BN/KL there are around 200 clumps moving into the interstellar medium with actual velocities between 100–300 km s⁻¹. They are moving away from a common origin with a set of protostars that apparently interacted in the past, following a homologous expansion law, better known as the Hubble Law. RO19a and RO19b calculated the initial velocity distribution of the clumps, considering their deceleration due to their interaction with the interstellar medium. In Figure 7 from RO19b, they showed an initial velocity distribution of the clumps in Orion BN/KL with two maximum values of ~200 and 400 km s⁻¹ for the global and local maximum, respectively. In this section, we use our dynamical interaction model to propose an ejection mechanism to account for the explosive properties of the Orion BN/KL outflow. We will analyze the importance of the projection angles and the mass of the clumps in the final velocity distributions and also compare this with the observational data.

### 4.1. Projection Angles

Using the results of the v200R1 model we obtain the projection, on the plane of the sky, of the position and velocity for each particle in our models. To achieve this, we must take into account the rotation projections on the plane of the sky to get the projected positions and velocities for each one of the particles. We use the rotation matrices $R_x$ and $R_z$, which are rotations with an angle $\theta$ and $\phi$ in the x- and z-axes.

Using the projected velocity, one can calculate the tangential velocities (the plane-of-the-sky velocities):

$$v_t = \sqrt{v_{\text{tot}}^2 - v_r^2}, \qquad (48)$$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Velocity distributions of the v200R1 model. The upper panels are the velocity as a function of the initial impact parameter and the final velocity distribution, with the description being the same as Figure 6. The bottom panels are the radial velocity (bottom-left panel) using $\theta = 30^{\circ}$ and $\phi = 30^{\circ}$, and the plane-of-the-sky velocity (bottom-right panel).</figcaption></figure>

with $v_{\text{tot}} = \sqrt{v_x^2 + v_y^2 + v_z^2}$; the radial velocities of each of the particles is given by the projected z-velocity,

$$v_r = v_z'. \qquad (49)$$

For this analysis, we present the result using the rotation angles $\theta = 30^{\circ}$ and $\phi = 30^{\circ}$. These angles are selected in order to obtain a wide tangential velocity distribution with a maximum value of around 200 km s$^{-1}$, like the case of Figure 7 in RO19b. It is important to note that it is not the aim of this paper to explore in detail the combination of angles and/or the precise initial clump distribution that reproduce the initial velocities proposed for this object. In that case, one must consider other effects, such as the gas dynamics, but we are interested in proposing this kind of explosion as a possible mechanism for the formation of this type of object.

Figure 7 shows the velocity as a function of distance (upper-left panel), the total velocity, radial velocity, and tangential velocity distributions (upper-right, lower-left, and lower-right panels, respectively).

Using Equation (49), we have calculated the radial velocity distribution in a range of velocities between about -110 and 270 km s$^{-1}$, with maximum values of about -100 and 150 km s$^{-1}$. Using these radial velocities in Equation (48), we obtained the velocity on the sky plane for each of the particles in the v200R1 numerical simulation set. The tangential velocity distribution of this model has a maximum value of about 150 km s$^{-1}$, and the shape of this tangential velocity distribution is similar to the initial velocity distribution presented in Figure 7 of RO19b, even when the distribution is spread over a large velocity range.

## 4.2. The Massive Clumps

However, the mass of clumps in the region of Orion BN/KL can be estimated by using the total mass of the moving gas in the region and, for simplicity, dividing it by the total number of current observed clumps. Additionally, RO19b predicted the

initial mass of the clumps (see Figure 6 of that paper). The initial mass of each individual clump is around $10^{-2} M_{\odot}$, and the total mass of these particles is comparable with the mass of the more massive particle (the star mass, i.e., $10 M_{\odot}$). The effects of the interaction between the low-mass particles and their contribution in the global motion of this event are not considered in the first models presented in this work. Nevertheless, the high velocity of the massive particle, and its momentum, play a more important role than the mass of the particles. In order to prove this, we ran a final model, v200R0m, where each of the low-mass particles have a mass of 0.01 $M_{\odot}$. We have assigned a random direction for the orbital velocity, assuming a circular motion according to Equation (46). The particles are in a quasi-equilibrium with each other, while the massive particle collides with this cluster and interacts with it. We have calculated the random positions and orbit directions of each of the particles following this procedure:

1. For each particle, assign a random radius using

$$r = \xi_c \eta^{1/(\alpha+3)} \qquad (50)$$

where $\eta$ is a uniform random number between 0 → 1.

2. We can calculate the Cartesian coordinates, x, y, and z, using

$$\begin{aligned} x &= r[\sin(\theta)\cos(\phi)], \\ y &= r[\sin(\theta)\sin(\phi)], \\ z &= r[\cos(\theta)] \end{aligned} \qquad (51)$$

where $\theta = \arccos(2\eta_t - 1)$ and $\phi = 2\pi\eta_p$, and $\eta_t$ and $\eta_p$ are uniform random numbers between 0 → 1.

3. We assigned an orbital velocity, $v_r$, using Equation (46).

4. We calculated $v_x$, $v_y$, and $v_z$ using

$$\begin{aligned} v_x &= v_r[\cos(\phi)\cos(\theta)\sin(\chi) + \sin(\phi)\cos(\chi)] \\ v_y &= v_r[\cos(\theta)\sin(\phi)\sin(\chi) - \cos(\phi)\cos(\chi)] \\ v_z &= -v_r\sin(\theta)\sin(\chi) \end{aligned} \qquad (52)$$

where $\chi = 2\pi\eta_x$ and $\eta_x$ is a uniform random number between 0 → 1.

Similar to previous models, we have run a set of 10 simulations using different random distributions. In this model we used a uniform distribution where the massive particle has an initial velocity of 200 km s$^{-1}$. Figure 8 shows the histogram of the total velocity for the v200R0m model as the solid line, the v200R0 model (with clumps of negligible mass) is shown by the dashed line, and the dashed-dotted line is the analytical solution (Section 2). The shape of these histograms is similar, but the v200R0m model presents a maximum of the distribution in a small velocity, about 130 km s$^{-1}$, 20 km s$^{-1}$ lower than the v200R0 model. The minimum velocity of any particle is about 80 km s$^{-1}$, this is also 30 km s$^{-1}$ lower than the v200R0 model; also, the particles with higher velocities are larger than in the model where there are negligible low-mass particles. The maximum velocity obtained in our numerical simulation is very similar to the maximum velocity predicted by RO19b (see Figure 7 in that paper), but the shape of the velocity distribution does not fully agree with that predicted in RO19b. However, as shown in the v200R1c0.5 and v200R1c1 models (Figure 6), the projected velocity distribution is also a function of the position through which the massive particle cross the initial particle distribution (Section 3.3) and the projection angle where the observed
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Velocity distribution of the v200R0m model. The solid and dashed lines are the results of the v200R0m and v200R0 models, respectively, and the dashed-dotted line is the analytical result obtained in Section 2.</figcaption></figure>

velocities are calculated (Section 4.1). However, a study of these parameters, for this particular object, is outside the scope of this paper and will be addressed in subsequent works.

Finally, Figure 9 shows the position and velocities of a single simulation of the v200R0m model, with projection angles of $\theta=0$ and $\phi=30, 45, \text{ and } 60$, for the upper, middle, and bottom panels, respectively. In this figure we plotted the tangential velocity of all the particles using red or blue arrows for positive and negative radial velocities, respectively. We also present, in the upper panel, the position of each particle when they crossed the y-axis (where they were blown away by dynamic interaction, according to Equation (14)). As one can see, the particles are ejected from a very small volume (about 0.3 au), which is insignificant with the size of the event after 500 years (about $4 \times 10^4$ au) and is in accordance with the observations that suggest a single ejection point, at least with the current resolution. Thus, these types of explosions seem to be ejected in a singular place in space, as well as the Orion BN/KL event.

However, the number of particles approaching or moving away from us is dependent on the projected angle, but the explosion is not at all isotropic in the x-axis, which is the axis of movement of the massive particle. But the morphology of the explosion is similar to the one observed in Orion BN/KL. It is important to note that the massive particle in the Orion BN/KL explosion should be a runaway massive star with a velocity of at least 150 km s$^{-1}$, but the dynamical interaction of the massive star with clumps with a total mass comparable to the star's mass could substantially decrease the velocity of the star at the end of the interaction. However, it is not the goal of this work to study the dynamic interaction of low-mass particles.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Position and velocities for each of the low-mass particles for the v200R1 model. The upper, middle, and bottom panels are the projected position of the each particles for $\phi = 30^\circ, 45^\circ, \text{ and } 60^\circ$, respectively. The blue and red arrows show the radial velocity direction, moving toward or away from us, respectively. The size of the arrows indicates the tangential velocity.</figcaption></figure>
## 5. Conclusion

We presented analytical and numerical solutions of the dispersion of particles due to dynamical interaction with a single massive particle. We have considered that the particles gather into a cluster and that they have a negligible mass with respect to the massive particle. The dynamical interaction with the massive particle produces a quasi-isotropic ejection of the particles.

Then, we carried out a set of numerical simulations of spherical distribution of massless particles (N-body simulations) for verifying our analytical solutions and obtaining an observational result. We have obtained a very good agreement between the numerical and analytical results.

The main conclusions are that:

1. The gravitational ejection mechanism is able to accelerate small clumps or cloudlets (i.e., low-mass gas fragments) to jetlike velocities and, therefore, it should be thoroughly explored in future work.

2. The terminal velocity $v_∞$ of each particle is function of its own impact parameter.

3. The maximum terminal velocity is given by the limit when the particle has an impact parameter equal to zero and it is twice the velocity of the massive particle ($v_{max} = 2v_0$). The minimal terminal velocity is related to the cluster radius.

4. The ejection angle of each of the particles is linearly related to the terminal velocity, and therefore is related to the impact parameter. Compact distributions of particles that are dynamically disturbed by a massive particle, which passes through the center of the distribution, produce more collimated ejections than in the case of more scattered clusters.

5. The resulting dispersion has an explosive signature, such as: (a) a small-scale common origin, (b) an isotropic distribution, and (c) velocities proportional to the distance from the common origin. This is the result of a short time interaction, which could be the mechanism that produced the explosive outflows.

6. The distribution of the ejection velocities is a function of the exponent of the initial distribution of particles ($α$), the gravitational and cluster radius, the massive particle velocity, and the terminal velocity of each of the particles. The minima and maxima, local or global, in these distributions can be obtained analytically.

7. The off-center dynamical interaction produces a wider velocity distribution, and with smaller velocities.

In addition, a dynamical interaction between a massive object with a cluster of less-massive particles is able to increase the individual energy of the cluster members producing an explosive event, instead of forming a new cluster with a massive particle or letting the massive particle cross the cluster with a minor perturbation.

Finally, we considered that the Orion BN/KL ejection was generated by a dynamical interaction. To demonstrate this, we have run a set of simulations where the true mass of gas in the fingers in Orion has been considered. Our models show that the

interaction of a massive particle with a distribution of particles with the same mass as that observed in the Orion fingers BN/KL region produces an ejection of material in all directions with a velocity distribution comparable to those observed in this region. Certainly, a study of parameters, mass distribution of the clumps, impact parameters of the massive star, morphology of the cluster of clumps, etc., in addition to the projection angles in which they are observed, are parameters that must be explored in detail in a future study.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program, for the Project The Dawn of Organic Chemistry (DOC), grant agreement No. 741002. We acknowledge support from PAPIIT-UNAM grants IN-109518 and IG-100218. The authors acknowledge Bertrand Lefloch for his useful comments to improve this work.

# Appendix
## N-body's Solver and Stability

The numerical method is a symmetrized leapfrog integrator with a variable time step formalism, which is second-order accurate and is able to preserve energy. For the N-body solution, we have considered N particles with masses $m_i$ and a position given by $x_i$, $y_i$, and $z_i$. The force between a pair of particles produces an acceleration, and the new position of each of the particles is strongly dependent on the time step $\Delta t$. A very large time step would solve incorrect trajectories and a small time step reproduces the real trajectory of each particle, dramatically increasing the computation time. In order to have an appropriate time step, we used the time step

$$ \Delta t = A^* \sqrt{\frac{R_{\min}}{a_{\max}}} \quad (A1) $$

where $A$ is a constant between $0 \to 1$, $R_{\min}$ is the mean distance between a pair of particles, and $a_{\max}$ is the maximum acceleration of a single particle. In order to prove the solutions of our N-body solver, during the simulation time we used a single-particle distribution of the v200R0 model, and we carried out the N-body simulation, using different values of $A$ (0.005, 0.05, 0.5). Figure 10 shows the relative position, $\Delta X$, $\Delta Y$, and $\Delta Z$ (for the left, center, and right panels, respectively), at an evolutionary time = 100 yr, with respect to the position of the model with a smaller $A$ value ($A = 0.005$), where $\Delta X = (X_A - X_{0.005})/X_{0.005}$, as for the other coordinates ($Y$ and $Z$). The plus symbols are used for the model with $A = 0.5$ and the diamond symbols are used for the results of the model with $A = 0.05$. The plot range (in the vertical axis) is $-1 \times 10^{-6} \to 1 \times 10^{-6}$, giving, then, a 0.001% difference between the models with $A = 0.5$ and 0.005. The small difference between models with different time steps added to the convergence with the theoretical results, guaranteeing that the code adequately solves the system of equations for the models used in this work.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Relative differences between models with different A constants, 0.5, 0.05, and 0.005 (see Equation (A1)). We have calculated the difference of each of the particles with respect to the smaller time-step model (A = 0.005). The left, center, and right panels are the differences in the X, Y, and Z coordinates. The plus and diamond symbols are the models with A = 0.5 and A = 0.05, respectively.</figcaption></figure>

ORCID iDs

P. R. Rivera-Ortiz https://orcid.org/0000-0001-6867-4210

A. Rodríguez-González https://orcid.org/0000-0002-2606-8527

Luis A. Zapata https://orcid.org/0000-0003-2343-7937

References

Allen, D. A., & Burton, M. G. 1993, Natur, 363, 54

Bally, J., Ginsburg, A., Arce, H., et al. 2017, ApJ, 837, 60

Bally, J., Ginsburg, A., Forbrich, J., et al. 2020, ApJ, 889, 178

Bally, J., Ginsburg, A., Silvia, D., & Youngblood, A. 2015, A&A, 579, 130

Becklin, E. E., & Neugebauer, G. 1967, ApJ, 149, 799

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Blaauw, A. 1961, BAN, 15, 265

Calkin, M. G. 1989, AmJPh, 57, 374

Cantó, J., Raga, A. C., Esquivel, A., et al. 2011, MNRAS, 418, 1238

Cunningham, N. 2006, PhD thesis, Univ. Colorado at Boulder

Dempsey, R., Zakamska, N. L., & Owen, J. E. 2020, MNRAS, 495, 1172

Dorigo Jones, J., Oey, M. S., Paggeot, K., et al. 2020, ApJ, 903, 43

Fall, S. M., Krumholz, M. R., & Matzner, C. D. 2010, ApJL, 710, L142

Irrgang, A., Kreuzer, S., & Heber, U. 2018, A&A, 620, A48

King, I. R. 1966, AJ, 71, 64

Kirk, H., Offner, S. S. R., & Redmond, K. J. 2014, MNRAS, 439, 1765

Kleinmann, D. E., & Low, F. J. 1967, ApJ, 149, L1

Kruijssen, J. M. D., Pelupessy, F. I., Lamer, H. J. G. L. M., et al. 2012, MNRAS, 421, 1927

Krumholz, M. R. 2014, PhR, 539, 49

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Morbidelli, A., & Nesvorny, D. 2012, A&A, 546, 18

Nissen, H. D., Gustafsson, M., Lemaire, J. L., et al. 2007, A&A, 466, 949

Perets, H. B., & Šubr, L. 2012, ApJ, 751, 133

Rivera-Ortiz, P., Rodríguez-González, A., Hernández-Martínez, L., & Cantó, J. 2019a, ApJ, 874, 38

Rivera-Ortiz, P. R., Rodríguez-González, A., Hernández-Martínez, L., et al. 2019b, ApJ, 885, 104

Rodríguez, L. F., Dzib, S. A., Loinard, L., Zapata, L., Gómez, L., Menten, K., & Lizano, S. 2017, ApJ, 834, 140

Rodríguez-González, A., Cantó, J., Esquivel, A., et al. 2007, MNRAS, 380, 1198

Sastry, K. S., & Alladin, S. M. 1970, Ap&SS, 7, 261

Spurzem, R., Giersz, M., Heggie, D. C., et al. 2009, ApJ, 697, 458

Wang, L., Kroupa, P., & Jerabkova, T. 2019, MNRAS, 484, 1843

Zapata, L. A., Ho, P. T. P., Fernández-López, M., et al. 2020, ApJL, 902, L47

Zapata, L. A., Ho, P. T. P., Guzmán Ccolque, E., et al. 2019, MNRAS, 486, L15

Zapata, L. A., Schmid-Burgk, J., Ho, P. T. P., Rodríguez, L. F., & Menten, K. M. 2009, ApJ, 704, 45

Zapata, L. A., Schmid-Burgk, J., Pérez-Goytia, N., et al. 2013, ApJL, 765 L29