What causes the fragmentation of debris streams in TDEs?

The polytropic sphere is initially relaxed without the black hole gravitational potential. A velocity damping is added in this phase in order to remove possible noise in the initial random displacement of the particles. After this first relaxation, we checked that the density profile of the sphere was the expected polytropic sphere, as described in Lodato et al. (2009).

In the case of a rotating star, after this first phase a rigid rotation is imposed and the sphere is relaxed once more, this time without the velocity damping (that would slow down the stellar rotation), similarly to Sacchi & Lodato 2019. The amount of rotation is described by the dimensionless parameter $\alpha=\omega/\omega_{\textup{b}}$, where $\omega$ is the stellar angular velocity and $\omega_{\textup{b}}$ is the break-up velocity, defined as the velocity at which the centrifugal force equals the stellar self-gravity. The non-rotating case corresponds to $\alpha=0$, while $\alpha=1$ indicates the maximally rotating case. The star is relaxed until it reaches a new equilibrium state that we monitor through the central density and thermal energy. Usually the value of $\alpha$ varies during the relaxation process (Sacchi & Lodato, 2019). However, the value of $\alpha$ in our simulation is relatively small, $\alpha=0.2$. For this slow value of stellar rotation a negligible reduction is observed during this second relaxation phase.

Once the star is relaxed, whether rotating or not, it is injected into a parabolic orbit around a $10^{6}\,M_{\odot}$ black hole, which is modelled as a Keplerian potential. The star is injected from a distance equal to 3 $r_{\textup{t}}$, with penetration factor $\beta=r_{\textup{t}}/r_{\textup{p}}=1$, where $\beta$ is the ratio between the tidal radius and pericentre distance $r_{\textup{p}}$. The penetration factor indicates how close to the hole the star passes during its orbit. Values of $\beta$ greater than 1 mean that the star undergoes a deep plunge in the black hole potential well, while values smaller than 1 mean a far away passage. A passage with $\beta=1$ guarantees a complete stellar disruption, for $\beta<1$ instead one only gets partial disruption or the stripping of stellar material (Guillochon & Ramirez-Ruiz, 2013).

After the star gets past the pericenter, it is disrupted. When it is sufficiently far from the source of the Keplerian potential, the latter is replaced with a sink particle of mass $10^{6}\,M_{\odot}$ and a relatively large accretion radius of 5 tidal radii. This means that every SPH particle that gets closer than 5 tidal radii from the black hole is considered to be accreted and removed from the simulation. This of course does not allow us to observe the physics of stream-stream shocks, circularization and accretion, however this choice is due to the fact that we are focusing our interest on the stream evolution properties.

The first crucial result obtained through our simulation is that the stream of gas generated by the tidal disruption of a non-rotating star modelled as a polytropic sphere with adiabatic index $\gamma=5/3$ is affected by gravitational instabilities that bring it to fragment into smaller almost spherical blobs, as already shown by Coughlin & Nixon (2015).

The presence of fragmentation is assessed via visual analysis of the stream appearance and through the analysis of the mean density and the density fluctuations, evaluated as the standard deviation of the mean density.

Figure 1 shows how initially both the mean value of the density (solid black line) and the fluctuations (dashed black line) fall as a power-law with power index $n\approx-1.67$. After almost 30 $t_{\textup{min}}$ the density fluctuations reach their minimum and start growing, soon overcoming the mean density of the stream. The turning point of the fluctuations is interpreted, in analogy with Cosmology (Peebles, 1980), as the moment where fragmentation starts.

In order to confirm this technique we analyse also the projected density of the stream taking snapshots or our simulation. Figure 2 shows a section of the stream at four different moments: $t=$10, 25, 35 and 55 $t_{\textup{min}}$ (that is as 13 months, almost 3 years, 4 years and 6 years after the stellar disruption). These moments corresponds to an instant way before the fragmentation, right before, right after and far after the fragmentation occurred. From these snapshots it is clear that the fragmentation occurs indeed when the density fluctuations reach the turning point.

The search for a criterion able to identify effectively the time at which the stream fragments is of particular interest as the one presented in the Introduction, suggested by Coughlin et al. (2016b), does not give not an overly accurate estimate. To be more precise, fig. 3 shows the mass distribution ${\rm d}M/{\rm d}(\rho r^{3})$ as a function of $\rho r^{3}$ for two stellar structures: the so far studied $\gamma=5/3$ and a more compact case $\gamma=2$, more susceptible to fragmentation. The distribution is shown at two times: when the stream is still far before the fragmentation point ($t=0.3\,t_{\textup{frag}}$, dashed line) and far after ($t=2\,t_{\textup{frag}}$, solid line). Fig. 3 shows how the criterion based on the tides, $\rho r^{3}$ be greater than the black hole mass $M_{\textup{h}}=10^{6}\,M_{\odot}$ (indicated by the vertical dotted line in the figure), is satisfied by almost all of the stream, eccept for a small fraction of it ($\lesssim 10\%$) already at the earliest time. The fact that the more compact scenario shows a distinct plateau at high densities after disruption might imply that a sharper criterion should have its critical point further up in the $\rho r^{3}$ ladder, instead of setting it at the black hole mass.

As suggested by Coughlin et al. (2016b), based on simulations described in Coughlin & Nixon (2015) and Coughlin et al. (2016a), the noise due to the numerical methods affects the time at which fragmentation occurs. The aforementioned convergence test has been performed in order to prove that the fragmentation is indeed physical and not only a result of numerical features of the code. The convergence has been considered satisfied when increasing the number of SPH particles employed in the simulation, the time at which fragmentation occurs would not increase appreciably.

Figure 4 shows the time at which fragmentation occurs (hereafter $t_{\textup{frag}}$) as a function of the number of SPH particle of our simulations. When the simulation is performed using a number of particles smaller than $\approx 10^{6}$, $t_{\textup{frag}}$ is strongly dependent on the resolution of the simulation, while above one million SPH particles it can be considered roughly constant. Note that the time of fragmentation occurs when the luminosity has already dropped by $\sim(30)^{-5/3}\sim 3\times 10^{-3}$ from the peak. We consider our choice of using, as mentioned in the previous section ($1.5\times 10^{6}$ particles), satisfactory.

We now consider the mechanisms responsible for the fragmentation. We compare the various contributions of the forces playing an active role in our picture by comparing their characteristic time scales.

Representing the stream as a cylindrical fluid is a good assumption for the dynamics of these events, as also pointed out by Coughlin & Nixon 2019. In their work, the authors wrote the Lagrangian describing the motion of the core of the stream considering the effects of the black hole and the core mass gravitational forces, focusing on the fallback rate temporal behaviour and the bound core fate.

In this work, we are interested in studying the fluid behaviour of the debris stream, its equilibrium state and its (eventual) fragmentation. This allow us to derive the temporal evolution of the stream quantities using simple physical considerations, introducing time scales related to the main forces at play and discussing their evolution with time. Four are the forces that are actively involved in determining the dynamics of a fluid stream. Apart from the self-gravity of the stream and the tidal effects of the black hole, the internal forces of the gas (i.e. its pressure) and the background stretching of the gas are also at play. Indeed, the stream is affected by a stretching in the volume along the radial direction of the black hole, due to the differential acceleration of the particles composing the stream. This phenomenon is very important because it affects all the quantities, introducing a time-dependency that becomes important as time passes by.

Each of these forces corresponds to a typical time scale: the stream self-gravity is linked to the free-fall time $t_{\textup{ff}}$, the tidal force generated by the black hole corresponds to a typical dynamical time $t_{\textup{d}}$, the gas pressure of the stream si associated to the sound crossing time $t_{\textup{sc}}$ and finally we have a stretching time scale $t_{\textup{s}}$. As pointed out, we are interested in the debris stream: in the following analysis all these quantities are always referred to and computed for the part of the stream with positive energy.

The first force we consider is the stream self-gravity. This force will obviously favor the fragmentation and gravitational collapse and thus the fragmentation of the debris. The free-fall time is given by (Jeans, 1902):

where $\rho$ is the density of the stream and the $2\pi$ factor accounts for the cylindrical geometry we are working with.

The tidal forces of the black hole tend to prevent any gravitational collapse. The dynamical time is computed as

where $\rho_{\textup{h}}=3M_{\textup{h}}/4\pi r^{3}$ is the “density" of the black hole over a sphere of radius $r$.

Pressure tends to stabilize the stream, the relative thermal time is just the sound crossing time across the transverse direction

where $H$ is the transverse size of the stream and $c_{\textup{s}}$ is the sound speed of the gas. The sound speed is calculated from the density, knowing the polytropic equation of state that describes the gas :

The stream width $H$ is obtained by averaging the distance of the particles composing the stream, considered again as a cylinder, from its axis: operatively we “sliced" the stream along its length and for each slice we computed the mean distance of the particles within the slice from the cylinder axis, finally we mediated all the radii obtained along the stream. Figure 5 shows $H$ as a function of the distance from the black hole normalized to the tidal radius of the original star (solid black line). The dashed black line shows the $r^{1/2}$ behaviour found by Coughlin et al. (2016b) in disagreement with the $r^{1/4}$ prediction of Kochanek (1994).

The last time scale is associated to the stream stretching along its axis. This apparent force is due to the fact that right after the stellar disruption each gas particle lies on a different orbit determined by its energy. The spread in energy of the debris is given by the change of the black hole gravitational potential across the star at the time of disruption, that is at the first pericenter passage, in the impulse approximation (Lacy et al., 1982; Rees, 1988). We suppose that the length $l$ of the stream changes with time according to the prescription $l=l_{0}a(t)$, where $l_{0}$ is the length of the stream at the initial time and $a(t)$ represents a dimensionless parametre of the stretching. This stretching introduces a proper timescale, defined in analogy with the Hubble time in Cosmology (Peebles, 1980) as

where the dot indicates the time derivative. The idea to consider the expansion and the contraction of the fluid volume in analogy with Cosmology using dimensionless parameters and an appropriate coordinate system has been already used to study the Solar wind fluctuations (Grappin & Velli, 1996; Landi et al., 2014; Del Zanna et al., 2015) as well as the collapse of pre-stellar cores (Toci et al., 2018).

Figure 6 shows the behaviour of the four described time scales as a function of time. The solid black line represents the tidal time, the dashed one the free-fall time, the dotted line is the thermal time and dash-dotted one is the stretching time. They are all normalized to the minimum return time (the time at which the first of the stellar debris comes back to its orbital pericenter, in our case $\approx 41$ days) while the time on the $x$-axis is normalized to the fragmentation time $t_{\textup{frag}}\approx 30\,t_{\textup{min}}$. The vertical dotted red line has been drawn to facilitate the comparison of the time scales at the fragmentation time.

The first consideration one can draw from figure 6 is that all the time scales have comparable magnitudes (within a factor $1.1$), crossing several times during the stream evolution. The only exception is the tidal time that is a factor $2.5-12$ longer than the others, although at the beginning of the simulation, this time scale should have been the smallest one in order for the star to be tidally disrupted. This immediately shows that tidal forces are not responsible for fragmentation, where it happens.

Apart from the tidal time scale, all of the other time scales share a very similar evolution. The thermal time is the smallest one along almost the entirety, hence pressure is the most relevant force of the stream evolution, and it prevents its fragmentation. The time scales that disentangles from the others is the stretching time, that accelerates its growth after roughly $4\,t_{\textup{min}}$. When fragmentation occurs the stretching time is almost a factor 2 larger then both the free-fall and thermal time, although it is still a 10 smaller then the tidal time.

The second important outcome of figure 6 is that the fragmentation occurs right after the free-fall time becomes shorter than the thermal time. This suggests that the fragmentation is a consequence of the pressure failing in balancing the stream self-gravity rather then an overcome of the black hole tidal forces or of the stream stretching by the gravitational collapse. This is better shown in figure 7: the solid black line represents the ratio between the free-fall time and the thermal time while the dashed black line the ratio between free-fall time and stretching time. The two red dotted lines highlight the time at which fragmentation occur (the vertical one) and the line on which the two considered time scales are equals (horizontal one).

As explained section 1, a stream of gas debris obeying to a polytropic equation of state with $\gamma\gtrsim 5/3$ is prone to fragmentation fueled by gravitational instabilities. The results of the previous section prove that also the $\gamma=5/3$ case is susceptible of fragmentation.

An interesting factor that can be taken into account and could potentially slow down the fragmentation process is the stellar rotation. In this paper, we will focus on a configuration where the stellar initial rotation shares direction and sense with the orbital angular momentum of the star itself.

We will not treat the case where the rotation is able to prevent the stellar tidal disruption. This occurs when the stellar rotation is sufficiently fast and its axis is parallel and opposite to the stellar orbital angular momentum, or within some tens degrees from this configuration. In this case the stream does not form or it is extremely faint (Sacchi & Lodato, 2019) and therefore is not relevant to the problem at hand.

If at least one component of the initial stellar rotation lies in the orbital plane there will be another force acting upon the debris stream: the centrifugal force inherited from the star. This force will add support against the collapse, counter-acting to some extent the self-gravity of the stream and making it harder for it to fragment.

In the case considered here, where the rotation axis is perpendicular to the orbital plane, no additional force is added to the picture described above, the only thing that gets modified is the stellar internal structure and mass-energy distribution, that gets wider. This is due to a swelling of the rotating star with respect to the non-rotating case. A wider mass-energy distribution would cause a faster stretching of the stream, thus slowing down the fragmentation process.

To better understand this configuration and test our prediction, we performed a numerical simulation, adding an initial stellar rotation ($\alpha=0.2$).

The simulations confirm our predictions about the role of stellar rotation aligned with the orbital angular momentum of the star: the time at which fragmentation occurs gets delayed by $\approx 15\%$. The delay is caused by the bigger stretching that a wider initial mass-energy distribution provokes.

Figure 8 shows the stretching parameter for the non rotating (solid black line) and rotating case (dashed black line) as a function of time (at the fragmentation time specific of each configuration). One can notice how, in the rotating case, the stretching is bigger, but only by $\approx 5\%$, this is however sufficient to produce the aforementioned amount of delay in the fragmentation time. This is due to the fact that the time scales have the same magnitude, thus a small changing in the value of one of the forces can lead to a significant difference in the evolutionary pattern of the stream.

Here, we consider the dynamics of the debris stream after disruption as composed of test particles freely falling in the radial direction in the gravitational field of the black hole. This approach is the same as the one of Coughlin et al. (2016b), who first considered the evolution of the structure of the debris stream in a semi-analytical way. They adopt the same simplification of radial freely falling particles and find that this is a very good approximation for the post-disruption hydrodynamical stream. Coughlin et al. (2016b) find approximate solutions for the position and velocity of the stream elements in the limit where the particles are close to the marginally bound orbit, that corresponds to the stream center of mass. Here, as described below, we generalize the solution of Coughlin et al. (2016b) and find an exact analytical solution for the whole of the stream. We then argue that it is the deviation from the “close to marginally bound” orbits that causes fragmentation.

The equation of motion of the debris is:

where bound, marginally bound, and unbound orbits are defined by $E<0$, $E=0$ and $E>0$, respectively. The solutions to these equations are well know since they are the standard Friedmann equations used in Cosmology to describe a matter dominated universe in the closed, flat and open case, respectively (Friedman, 1922). For the marginally bound orbit (corresponding to a flat Universe) we have (see also eq. 8 in Coughlin et al. 2016b)

for the unbound debris (corresponding to an open Universe) we have a parametric solution:

while for the bound debris (corresponding to a closed Universe) we have:

where $r_{0}=GM/|2E|$ and $t_{0}=GM/|2E|^{3/2}$ are the definition of a scale radius and a scale time, that depend on the energy of the debris $E$. These expressions give parametrically the exact positions of the debris elements across the whole stream (cf. the approximate positions given in Eq. 19 in Coughlin et al. 2016b)

It is instructive to approximate the above solutions for small $\eta$, which corresponds to orbits close to the marginally bound case. This is obtained by expanding the hyperbolic and the trigonometric sine and cosine as:

In this limit, for both bound and unbound orbits, we obtain:

The relative position with respect to the center of mass of the stream is then

that implies that $a\propto t^{4/3}$ (cf. again Coughlin et al. 2016b, their eq. 16). Actually, our exact solution allows us to compute the evolution of $a(t)$ for the whole stream. Figure 11 (left panel) shows the logarithmic derivative of the stretching parameter $a(t)$ for the bound (dashed line) and the unbound (solid line) portion of the stream. One can see that initially both are close to the expected value of 4/3, while at late times (which occurs progressively earlier in physical time the farther we move away from the marginally bound orbit) the stretching of the unbound debris slows down, eventually reaching freely streaming orbits, where $a(t)\propto t$, while the bound debris are stretched faster, as the tide of the black hole increases. Similarly, one can compute the stretching timescale $t_{\rm s}=a/\dot{a}$, shown in Fig. 11 (right panel), which grows faster (slower) than $3t/4$ for the unbound (bound) portion of the stream.

Coughlin et al. 2016b obtain their evolution by solving in a relatively complicated way the equation of motion of the debris, by introducing a parameter $\xi$, defined as:

which is a measure of how far from the marginally bound orbit a stream element is. They further make the reasonable but in principle not justified assumption that the stream velocity follows a self-similar solution:

and finally solve numerically for $f$ from an ordinary differential equation. They also provide an approximate solution for $f(\xi)$ in the form $f=1/\xi-1$, which they then use to compute the evolution of the stream far from the marginally bound orbit.

Our approach allows us to bypass all this and obtain an analytical and exact solution for $f(\xi)$, at the same time allowing us to demonstrate that the radial velocity is indeed self-similar.

Let us first consider the behaviour close to the marginally bound orbit, Eq. (18). The radial velocity is readily obtained:

and, recalling the definition of $\xi$, we obtain

which has the properties $f(2/3)=1$ and $f^{\prime}(2/3)=-5/2$ (cf. Coughlin et al. 2016b). We thus see that, indeed, the ansatz by Coughlin et al. (2016b) that the solution would be self-similar with respect to the variable $\xi$ is indeed correct, at least for particles close to the marginally bound orbit, as it only depends on how far does the particle lie with respect to the stream center of mass. We shall see in a moment that the self-similarity extends to the whole stream exactly.

The function in Eq. (23) approximates very well the full solution in the vicinity of the marginally bound orbit but fails farther from it. Actually, an exact and closed form analytical expression for $f(\xi)$ can be obtained by differentiating Eqs. 14 and 15:

After some simple algebra, we can thus obtain $f(\xi)$ parametrically:

for the bound portion of the stream ($\xi>2/3$) and

for the unbound portion of the stream ($\xi<2/3$). In Figure 12 we plot the analytical function $f(\xi)$ as described above (dashed black line), the numerically computed function based on solving eq. (4) in Coughlin et al. (2016b) (solid red line), and the two approximations: (i) the function $f=1/\xi-1$ proposed by Coughlin et al. (2016b) (dotted line) and (ii) the function we propose in Eq. (23) above (dash-dotted line). As we can see, Eq. (23) provides a better approximation (with respect to $f=1/\xi-1$) to the actual solution close to the marginally bound orbits, but fails to reproduce the asymptotic behaviour, especially for the bound portion of the stream ($\xi\gg 1$). On the contrary, the exact solution, Eqs. (25) and (26), does recover the numerically computed one. We also plot in Fig. 13 the resulting radial velocity at a given time $t\approx 25\,t_{\rm min}$, where in red we show the distribution obtained from our SPH simulation and the various lines indicate the distribution obtained using our approximate solution close to the marginally bound debris (Eq. (23), dash-dotted line) and the one from Coughlin et al. (2016b) approximated function described above (dotted line). We do not plot our exact solution because it coincides exactly with the numerical simulation. We also indicate, with a dashed line, the marginally bound orbit velocity profile. The intersection of this line with the simulation data marks the location of the marginally bound debris. Again, the general conclusion is that our exact solution is an excellent representation of the simulation data, that our approximate expression (Eq. (23) is a good approximation close to the marginally bound debris but fails away from it, and that the Coughlin et al. (2016b) proposed function, while not approximating the actual solution at any point, gives a fair first order description of the asymptotic regimes.

We are now in a better position to evaluate the stability of the stream. We have seen that the stretching of the orbits initially proceeds as $t^{4/3}$, but then slows down in the unbound portion of the stream and accelerates in the bound portion.

Now that we know the time dependency of the stream density and the scaling factors, it is straightforward to derive all of the other quantities and time-scales behaviours.

In Section 3 we introduced four time scales, each linked to one of the force acting on the stream. Here we will neglect the tidal time as this is much larger than the others and thus is not expected to play a role in the fragmentation process.

We have then the stretching time, the thermal one and finally the free-fall time. The first two stabilise the stream while the latter is ultimately responsible for fragmentation.

The stretching time dependence is the easiest to derive as this time scale is simply defined as $t_{\textup{s}}=a/\dot{a}$ and therefore

at least initially, bearing no dependence on the polytropic index.

We recall that the thermal time is defined as $t_{\textup{p}}=H/c_{\textup{s}}\propto 1/\sqrt{\rho}$, where we used the condition of hydrostatic balance in the transverse direction. Since $\rho\propto a^{1}/(1-\gamma)$ (Eq. 10)

Finally, also the free-fall time scale, defined as $t_{\textup{ff}}=1/\sqrt{2\pi G\rho}$ is proportional to $1/\sqrt{\rho}$ and thus scales with time in the same way as the thermal time. Actually, the assumption of hydrostatic balance in the transverse direction implies that the stream is always close to marginal stability, according to the Ostriker (1964) criterion, $\Lambda\sim c_{\rm s}^{2}/2G$ and thus fragmentation would ensue relatively easily. However, as long as the stream stretching occurs on the same timescale, this can act to stabilise it.

It is therefore clear why the $\gamma=5/3$ case is the critical one: for this particular value of the polytropic index all of the relevant time scales share the very same time dependence, at least initially: $t_{\textup{ff,sc,s}}\propto t$, the stream stays therefore marginally stable for a long time. After some time however the time evolution of the thermal and free fall time scales slows down and they disentangle from the stretching time. This can be seen in figure 6, and is due to the fact that at late times (when the stream starts to freely stream) they tend to be proportional to $t^{3/4}$ while the stretching time remains proportional to $t$.

In the case of more compact stars ($\gamma>5/3$) the stream stretching along its axis is less efficient than all of the others time scales and therefore the streams fragments much faster than in the $\gamma=5/3$ case. In all of the cases where $\gamma<5/3$ the stretching grows faster than the other forces and is therefore able to sustain the stream and prevent its fragmentation.

Figure 14 shows with a solid black line the quantity described in (7) while the two dotted red lines highlight again the fragmentation time and the level at which the condition is satisfied. We can see that the condition is almost satisfied during all of the stream evolution. It rapidly grows above unity right before the stream collapses, hence confirming that the stream can be considered in hydrostatic equilibrium for the part of its evolution we are interested in.

All of the arguments presented above consider only the unbound part of the stream. Employing our analytical results we can however infer that in the bound part the stretching time falls rapidly below all of the other time-scale and thus should prevent the fragmentation of that portion of the stream. This is indeed what we observe in our simulations, too. Figure 15 shows the density fluctuations in the unbound (solid black line) and bound (dashed black line) part of the stream as a function of time. The absence of a turning point in the stream of bound material signals that little to no fragmentation is found in the bound debris before they enter the region in which they are considered to be accreted.