Gravitational wave emission from unstable accretion discs in tidal disruption events

TL;DR

This paper investigates gravitational wave emissions from unstable accretion disks formed after tidal disruption events, estimating their detectability by LISA through analytical and numerical methods.

Contribution

It provides the first combined analytical and numerical analysis of gravitational waves from unstable accretion disks post-tidal disruption, highlighting potential detectability.

Findings

Numerical strain measurements are two orders of magnitude lower than analytical estimates.

Disks affected by the Papaloizou-Pringle instability could be detectable by LISA in certain scenarios.

Analytical estimates suggest possible gravitational wave signals from these disks.

Abstract

Gravitational waves can be emitted by accretion discs if they undergo instabilities that generate a time varying mass quadrupole. In this work we investigate the gravitational signal generated by a thick accretion disc of $1 M_{\odot}$ around a static super-massive black hole of $10^{6}M_{\odot}$, assumed to be formed after the tidal disruption of a solar type star. This torus has been shown to be unstable to a global non-axisymmetric hydrodynamic instability, the Papaloizou-Pringle instability, in the case where it is not already accreting and has a weak magnetic field. We start by deriving analytical estimates of the maximum amplitude of the gravitational wave signal, with the aim to establish its detectability by the Laser Interferometer Space Antenna (LISA). Then, we compare these estimates with those obtained through a numerical simulation of the torus, made with a 3D smoothed particle hydrodynamics code. Our numerical analysis shows that the measured strain is two orders of magnitude lower than the maximum value obtained analytically. However, accretion discs affected by the Papaloizou-Pringle instability may still be interesting sources for LISA, if we consider discs generated after deeply penetrating tidal disruptions of main sequence stars of higher mass.

Equations

$$g_{\mathrm{sg}}=\frac{G{M}_{\ast }}{R_{\ast }^2},$$

$$g_{\mathrm{t}}\approx \frac{\partial g_{\mathrm{h}}}{\partial r}{R}_{\ast }=-\frac{\partial }{\partial r}\left(\frac{G{M}_{\mathrm{h}}}{r^2}\right)R_{\ast }=\frac{2G{M}_{\mathrm{h}}}{r^3}{R}_{\ast },$$

$$r_{\mathrm{t}}=R_{\ast }{\left(\frac{M_{\mathrm{h}}}{M_{\ast }}\right)}^{1/3}\simeq 7\times 10^{12}{r}_{\ast }{\left(\frac{M_6}{m_{\ast }}\right)}^{1/3}\text{cm},$$

$$r_{\mathrm{p}}\le 2r_{\mathrm{g}}=\frac{2G{M}_{\mathrm{h}}}{c^2}=3M_6\times 10^{11}\text{cm},$$

$$\beta =\frac{r_{\mathrm{t}}}{r_{\mathrm{p}}}.$$

$$1\le \beta \lesssim 20\times r_{\ast }{M}_6^{-2/3}{m}_{\ast }^{-1/3}.$$

$$\Omega (R)={\Omega }_0{\left(\frac{R_0}{R}\right)}^2,$$

$$\frac{(\omega +m{\Omega }_{-}{)}^2+m{g}_{-}/R_{-}}{(\omega +m{\Omega }_{+}{)}^2+m{g}_{+}/R_{+}}={\left(\frac{R_{+}}{R_{-}}\right)}^{2m}\frac{(\omega +m{\Omega }_{-}{)}^2-m{g}_{-}/R_{-}}{(\omega +m{\Omega }_{+}{)}^2-m{g}_{+}/R_{+}},$$

$$g(R)=\frac{GM}{R_0^2}{\left[\frac{R_0}{R}\right]}^3\left[1-\frac{R}{R_0}\right].$$

$$\square {\bar{h}}_{\mu \nu }=-\frac{16\pi G}{c^4}{T}_{\mu \nu },$$

$${\partial }^{\mu }{\bar{h}}_{\mu \nu }=0,$$

$${\bar{h}}_{\mu \nu }=h_{\mu \nu }-\frac{1}{2}{\eta }_{\mu \nu }h,$$

$$h_{ij}^{\mathrm{TT}}(t,\mathbf{x})\simeq \frac{4G}{D{c}^4}{\Lambda }_{ij,kl}(\mathbf{n}){\int }_{|x^{\prime }|<s}{\text{d}}^3{x}^{\prime }{T}^{kl}\left(t-\frac{D}{c}+\frac{x^{\prime }\cdot \mathbf{n}}{c};x^{\prime }\right),$$

$$h_{ij}^{\mathrm{TT}}(t,\text{n})=\frac{1}{D}\frac{2G}{c^4}{\Lambda }_{ij,kl}(\text{n}){\ddot{M}}^{kl}\left(t-\frac{D}{c}\right).$$

$$M^{kl}=\frac{1}{c^2}\int {\text{d}}^{3}x{T}_{00}{x}^k{x}^l,$$

$$h_{+}$$

$$h_{\times }$$

$$h\approx (h_{+}^2+h_{\times }^2{)}^{1/2}.$$

$$h\simeq \frac{1}{D}\frac{4G}{c^2}\frac{E_{\mathrm{kin}}}{c^2},$$

$$h_{\mathrm{n}}(f{)}^2=f{S}_{\mathrm{n}}(f),$$

$$|h_{\mathrm{c}}(f){|}^2=4f^2|\tilde{h}(f){|}^2,$$

$$E_{\mathrm{kin}}=\frac{1}{2}{M}_{\ast }{v}_{\mathrm{kepl}}^2=\frac{1}{2}{M}_{\ast }\frac{G{M}_{\mathrm{h}}}{r_{\mathrm{p}}},$$

$$h=\frac{r_{\mathrm{s}\ast }{r}_{\mathrm{gh}}}{D{r}_{\mathrm{p}}},$$

$$f=\frac{1}{2\pi }{\left(\frac{G{M}_h}{r_{\mathrm{p}}^3}\right)}^{1/2}.$$

$$\frac{M_{\ast }}{M_{\text{sun}}}\approx \frac{R_{\ast }}{R_{\text{sun}}}\Rightarrow m_{\ast }\approx r_{\ast },$$

$$h\approx 10^{-22}\times \beta {\left(\frac{D}{16\text{Mpc}}\right)}^{-1}{m}_{\ast }^{1/3}{M}_6^{2/3},$$

$$f\approx 10^{-4}\text{Hz}\times {\beta }^{3/2}{m}_{\ast }^{-1}.$$

$$h_{\mathrm{c}}\approx h{\mathcal{N}}_{\mathrm{c}}^{1/2},$$

$${\mathcal{N}}_{\mathrm{c}}\approx f\tau ,$$

$$1\le \beta \lesssim 20.$$

$${\mathcal{N}}_{\mathrm{c}}\approx 57.$$

$$h_{\mathrm{c}}\approx 4\times 10^{-21},f=10^{-3}\text{Hz}.$$

Full text

††pubyear: 2019††pagerange: Gravitational wave emission \textcolorblackfrom unstable accretion discs in tidal disruption events–A

1 Introduction

Tidal disruption events (TDEs) occur when a star wanders too close to a super-massive black hole (SMBH), getting disrupted by the tidal field of the hole. After the disruption, roughly half of the star circularizes and forms a disc-shaped structure that accretes onto the SMBH, while the other half escapes on hyperbolic orbits with different orbital energies. These events are important electromagnetic sources, in UV-optical (e.g. Gezari et al. 2008, Komossa et al. 2008, Gezari et al. 2009, Gezari et al. 2017), in X-rays (e.g. Bade et al. 1996, Komossa & Greiner 1999, Komossa & Bade 1999, Greiner et al. 2000), in radio (e.g. Zauderer et al. 2011) and also in $\gamma$-rays (e.g. Bloom et al. 2011, Levan et al. 2011 and Cenko et al. 2012). The bolometric lightcurve of these events is expected to decrease with time as $t^{-5/3}$, like illustrated by Rees (1988) and Phinney (1989). However, recent studies (Lodato et al. 2009, Guillochon & Ramirez-Ruiz 2013) have shown that this trend is reached only at late times and that lightcurves in specific bands might show a different evolution (e.g. Lodato & Rossi 2011).

The mechanism of accretion onto the central object may affect the formation of the accretion disc. In particular there are three parameters that we need to consider in this process: the circularization time, $t_{\rm circ}$, the viscous accretion time $t_{\rm visc}$ and the radiative cooling time $t_{\rm cool}$ (Evans & Kochanek 1989, Bonnerot et al. 2016). If $t_{\rm visc}>t_{\rm circ}$, the accretion process starts only when the disc is formed. In this scenario, we have two possibilities: if $t_{\rm cool}<t_{\rm circ}$ the disc is thin, otherwise the disc is thick, so we have a torus (see, e.g., Lodato 2007).

If we consider a thick disc, the evolution of the system may undergo a global non-axisymmetric hydrodynamic instability, that is the Papaloizou-Pringle instability (PPI), first described by Papaloizou & Pringle (1984). For this instability to arise, the torus needs to have a shallow specific angular momentum profile and a well-defined inner and outer radii. Recently, Nealon et al. (2018), using initial conditions motivated by Bonnerot et al. (2016), have shown that a torus formed after a TDE may be susceptible to the PPI, with a frequency close to the Keplerian one and a radius approximately two times the pericenter of the initial stellar orbit. Additionally, for very low initial magnetic fields\textcolorblack, they also suggested that the PPI may drive super-Eddington accretion onto the SMBH faster than the magneto rotational instability (MRI). \textcolorblackAround the same time, magnetohydrodynamic simulations of PPI susceptible tori by Bugli et al. (2018) illustrated that the presence of magnetic fields is able to quickly suppress the PPI.

The PPI generates a time dependent density distribution, that leads to the emission of gravitational waves (GWs). These waves might be detected by the Laser Interferometer Space Antenna (LISA)111https://www.elisascience.org.To date, the GW signal from unstable accretion discs formed during TDEs has not been investigated yet. Anyway there are few studies on the GW emission associated with the phase of tidal disruption of the star. Some of these are based on full general relativistic hydrodynamics codes (see, e.g., Haas et al. 2012, Anninos et al. 2018), other on smoothed particle hydrodynamics (SPH) codes (Rosswog et al. 2009, Kobayashi et al. 2004).

Haas et al. (2012) focus on ultra-close TDEs of white dwarfs (WDs) by a rotating intermediate-massive black hole (IMBH). They assume to have a $10^{3}\text{M}_{\sun}$ BH, with $1\text{M}_{\sun}$ WD of radius $6000\,\text{km}$ on a parabolic orbit around it, at an average distance from us of $30\,\text{kpc}$. They find that the signal is a burst, with an amplitude $\approx 10^{-18}$ and frequency of few Hz. Moreover, they illustrate that the BH spin does not affect the GW signal in a significant way. Anninos et al. (2018) perform simulations of $0.2-0.6\text{M}_{\sun}$ WD on a parabolic orbit around a static $10^{3}-10^{4}\text{M}_{\sun}$ BH, assuming that the source is at an average distance of $10\,\text{Mpc}$. They derive that the GW strain is $\approx 10^{-22}$ and the frequency of the GW signal is $10^{-2}-10^{-1}\,\text{Hz}$, so these sources might be possible interesting targets for LISA.

On the other side, Rosswog et al. (2009), with an SPH code, have simulated WDs tidally disrupted by IMBHs obtaining the same conclusions as Anninos et al. (2018). Instead Kobayashi et al. (2004), using an SPH code implemented with general relativity, have investigated the GW signal emitted during TDEs of main sequence (MS) and helium stars around a SMBH (both static and rotating) of mass $10^{6}\text{M}_{\sun}$, showing that it might be detectable by LISA, in the case of strong encounters and assuming that these events take place at a distance $\leq 20\,\text{Mpc}$. They have also shown that the signal is in most part insensitive to the particular equation of state and structure of the star. Moreover, they have calculated that, for a solar type star, the peak value for the GW strain is around $\sim 10^{-22}$, in the case of static SMBHs, while it is approximately one order of magnitude higher, $\sim 10^{-21}$, in the case of rotating holes. They also have studied the signal analytically, approximating the star with a point-mass particle, and they have found that the numerical results are well \textcolorblackdescribed by this approximation. \textcolorblackIn the field of GW astronomy, all these works are important since underline the basic features of the GW emission associated with TDEs. In particular they all show that the strain associated with the phase of disruption has a similar shape, independent on the particular structure of the star. In the context of this paper, where we focus on disruption of solar type stars from SMBHs, the most relevant work is that of Kobayashi et al. (2004), of which our work can be seen as an extension to describe what happens after the disruption itself.

\textcolorblackAs for the GW signal generated by the PPI instead, Kiuchi et al. (2011) have illustrated that the PPI in thick self-gravitating tori may produce a detectable signal for GW detectors, both at high frequencies ($100\sim 200\,\text{Hz}$) in the case of stellar BHs, and in the low frequencies range (mHz) in the case of SMBHs. For this latter case, which is the one that LISA might observe, they have in particular considered a system formed by a SMBH of $10^{6}$M*☉* and a thick accretion disc of $\sim 10^{5}$M*☉*, system that they assume it might have formed by the collapse of a super-massive star. The orbital radius of the torus is $\approx 10^{7}\,\text{km}$ and the system is at a distance of $10\,\text{Gpc}$. They have made both a numerical study, with a general relativistic grid code, and an analytical study, finding that, if they assume to have the PPI with one over-density, the GW signal peak could reach $10^{-19}-10^{-18}$, with a frequency of $\sim 10^{-3}\,\text{Hz}$. While Kobayashi et al. (2004) have found a good \textcolorblackagreement between the numerical and the analytical results for the GW signal from the stellar disruption phase, for the signal generated by an unstable accretion disc Kiuchi et al. (2011)’s analytical estimates are one order of magnitude higher than their numerical results.

In this paper we link the previous work of Kobayashi et al. (2004) and Kiuchi et al. (2011), since we investigate the GW signal from tori formed by TDEs \textcolorblackof MS stars. In particular, we start from $1\text{M}_{\sun}$ TDE remnant unstable to the PPI, as shown in Nealon et al. (2018), with no self-gravity, around a non-rotating $10^{6}\text{M}_{\sun}$ SMBH. We study this signal both through an analytical and a numerical analysis.

The structure of the paper is the following. In Section 2 we illustrate the theory behind our work, i.e. TDE physics, the main features of the PPI and the basics of GW emission. In Section 3 we describe our analytical estimates, while in Section 4 the numerical study. In Sections 5 and 6 we discuss our results and give our conclusions respectively.

2 Theory

2.1 Tidal disruption physics

Here we consider a standard scenario for stellar disruptions by a SMBH, where a solar type star of mass $M_{*}$ and radius $R_{*}$ is on a parabolic orbit around the hole of mass $M_{\rm h}$. Let us suppose to work under the impulse approximation, i.e. the interaction between the star and the SMBH takes place only at the pericenter of the stellar orbit, $r_{\rm p}$. To determine the minimum approach at which the star can be tidally disrupted, we have to equate the gravitational acceleration on the stellar surface due to the self-gravity, $g_{\rm sg}$,

[TABLE]

where $G$ is the gravitational constant, to the tidal field exerted on the star, $g_{\rm t}$,

[TABLE]

where $g_{\rm h}$ is the gravitational acceleration due to the presence of the central object and $r$ is the distance between the centre of mass of the star and the SMBH. Comparing equation (1) and (2), we obtain that the two fields are comparable at a distance $r_{\rm t}$, called tidal radius, given by

[TABLE]

where we have introduced the dimensionless parameters $m_{*}$, $r_{*}$ and $M_{6}$, such that $M_{*}=m_{*}\,\text{M}_{\sun}$, $R_{*}=r_{*}\,\text{R}_{\sun}$ and $M_{\rm h}=M_{6}\times 10^{6}\,\text{M}_{\sun}$.

Thus, the disruption of the star occurs when $r_{\rm p}\leq r_{\rm t}$. However, for deeply penetrating events, the star may be directly swallowed by the hole. This occurs when

[TABLE]

where $r_{\rm g}$ is the gravitational radius of the SMBH and $c$ is the speed of light. To quantify how close the star is with respect to the SMBH and consequently the strength of the disruption, we define a dimensionless parameter, called the penetration factor $\beta$, as

[TABLE]

Given the two constraints above, that $r_{\rm p}\leq r_{\rm t}$ and that $r_{\rm p}>2r_{\rm g}$, $\beta$ can vary in the range

[TABLE]

2.2 The Papaloizou Pringle instability

The PPI was first studied by Papaloizou & Pringle (1984) and Blaes & Glatzel (1986), but a more simplified description of the PPI is covered in Pringle & King (2007). They consider a cylindrical flow of an incompressible fluid, with no $z$-dependence, that satisfies the Rayleigh criterion. In particular, they assume that this fluid has a rotational velocity, $\Omega$, given by

[TABLE]

where $R$ is the cylindrical radius, $R_{0}$ is a reference radius and $\Omega_{0}$ is the angular velocity at $R_{0}$. The solution of the perturbation equation for this flow is given by (e.g. Blaes & Glatzel 1986, Pringle & King 2007)

[TABLE]

where $\omega$ is the mode frequency, $m$ is the azimuthal wavenumber, $R_{-}$ and $R_{+}$ are the inner and outer radii respectively and $g$ is the effective gravity

[TABLE]

Equation (8) tells us that the growth rate of the unstable modes is independent of the mass of the torus, instead depending on the geometry of the torus as defined by the inner and outer boundaries.

One of the solutions of equation (8) is an unstable growing mode; this takes place when there is an exchange of angular momentum and energy between a \textcolorblacknegative-energy wave, that travels from the inner edge of the torus \textcolorblack(where $\omega<\Omega$), with a wave that travels from the outer edge \textcolorblack(where $\omega>\Omega$), in the opposite direction \textcolorblackwith respect to the local medium but with the same angular phase speed of the \textcolorblacknegative-energy wave. \textcolorblackThis redistribution of momentum and energy happens at the corotation radius $R_{\rm cor}$, defined as the radius where $\omega=\Omega$. When this happens, we have the PPI. In particular Nealon et al. (2018) have found that, for the thick torus they have studied, the fastest unstable growing mode corresponds to $m=1$ and manifests as one over\textcolorblack-dense region forming in the torus.

2.3 Gravitational wave emission

In general relativity (GR), Einstein’s equations, under the weak field approximation, become a wave equation plus a gauge condition222The Greek indices range from 0 to 3 (the four space-time coordinates), while the Latin indices from 1 to 3 (only spatial coordinates). (Einstein, 1918)

[TABLE]

where $T_{\mu\nu}$ is the stress-energy tensor and $\bar{h}_{\mu\nu}$ is defined as

[TABLE]

where $h_{\mu\nu}$ are small perturbations of the Minkowskian metric $\eta_{\mu\nu}$ and $h$ is trace of $h_{\mu\nu}$. If we assume that (e.g. Buonanno 2007)

the internal motion of the source is slow compared to the speed of light, 2. 2.

the self gravity of the source is negligible, 3. 3.

the signal is detected at a distance $D$ very far from the source, 4. 4.

the Transverse-Traceless gauge holds, where only two components of $h_{\mu\nu}$ are independent,333In the TT gauge we have $\bar{h}_{ij}$=$h_{ij}$.

the solution of the wave equation may be written as

[TABLE]

where $\Lambda_{ij,kl}$ is the TT-operator, $s$ is the typical size of the system and $\boldsymbol{n}$ is the direction of propagation of the wave. With the previous assumptions, it is possible to write a multipole expansion of $h^{\rm TT}_{ij}$; however, the monopole term is zero (due to the conservation of mass of the system), and the dipole term is also zero (due to the linear momentum conservation). For these reasons, the first non-vanishing term in the expansion is the quadrupole term and so, after some algebraic passages (see, e.g., Buonanno 2007), equation (13) reads

[TABLE]

where $\ddot{M}^{kl}$ is the second time derivative of the moment of inertia of the system, $M^{kl}$, defined as

[TABLE]

with $T_{00}/c^{2}=\rho$.

As shown in Buonanno (2007) and Maggiore (2007), from equation (13) we can derive the expression for a wave propagating in any direction. In particular, if we consider a wave along the $\boldsymbol{z}$ direction, we get the following expressions for \textcolorblack$h_{11}\doteq h_{+}$ and \textcolorblack$h_{12}\doteq h_{\times}$ (see Maggiore 2007)

[TABLE]

The GW strain can be calculated as

[TABLE]

To make a simple estimate of the GW amplitude associated to a source, we can approximate equation (18) using the expression given in Thorne (1987),

[TABLE]

where $E_{\rm kin}$ is the kinetic energy of the moving source. In order to see if this peak value may be detectable by LISA, we need to compare it with the sensitivity curve of the instrument. This curve (see Amaro-Seoane et al. 2017) is calculated in terms of the characteristic amplitude of the noise, $h_{\rm n}$, as a function of the frequency, $f$. In particular, $h_{\rm n}(f)$ is defined as (Moore et al. 2014, Maggiore 2018)

[TABLE]

where $S_{\rm n}(f)$ is the noise spectral density. The quantity related to the GW signal that we compare to $h_{\rm n}(f)$ is the characteristic amplitude of the signal, $h_{\rm c}$, defined as (Maggiore 2018)

[TABLE]

where $\tilde{h}(f)$ is the Fourier transform of the strain.

3 Method

In this work we study the GW signal associated to $1$M*☉* TDE remnant, simulated by Nealon et al. (2018), unstable to the PPI with $m=1$. They assume that the remnant is not magnetised and that the star disrupted by the SMBH is a solar type star, with $\beta=5$, as in Bonnerot et al. (2016). They study the evolution of the system for 20 orbits, although the growth of the unstable PPI modes is expected to be suppressed by the MRI before this.

Our work is divided into two parts: an analytical study, to have some first estimates of the expected signal, and a numerical study, to derive the GW strain associated to the torus.

3.1 Analytical study

For the analytical estimates, we proceed as follows. The PPI involves a displacement of a mass $\approx 1\text{M}_{\sun}$, moving roughly on a Keplerian orbit around the SMBH at a distance $\approx 2r_{\rm p}$ (see Bonnerot et al. 2016). A similar argument has been made by Kobayashi et al. (2004), when considering the GW signal of the disruption phase. Indeed, the GW signal associated to the stellar disruption can be estimated through equation (19), expressing $E_{\rm kin}$ as

[TABLE]

where $v_{\rm kepl}$ is the Keplerian velocity of the star. Thus, substituting equation (22) into equation (19), we find that

[TABLE]

where $r_{\rm s*}$ is the Schwarzschild radius of the star and $r_{\rm gh}$ is the gravitational radius of the hole. We expect the frequency of the signal to be, approximately, the Keplerian frequency

[TABLE]

In particular, if we consider a MS star, for which it is reasonable to assume that mass and radius scales in the same way, i.e.

[TABLE]

we have that equation (23) and equation (24) become

[TABLE]

To compare these estimates with the LISA sensitivity curve, we need to derive $h_{\rm c}$ through equation (21). It can be shown (Maggiore 2018) that this equation may be expressed in a more useful way as

[TABLE]

where $\mathcal{N}_{\rm c}$ is the number of cycles spent in the detector bandwidth $[f_{\rm min},f_{\rm max}]$ (Maggiore 2007). Since our source is monochromatic, we can approximate this quantity as

[TABLE]

where $\tau$ is the emission time (Colpi & Sesana 2017).

With our analytical estimates we are assuming the raw strain of the GWs to be the same in the case of TDEs and unstable accretion discs. However, we see that this is not true for $h_{\rm c}$, since in the case of TDEs (as explored by Kobayashi et al. 2004) we have a burst signal with a duration that is simply the inverse of the frequency. This means that, in this scenario, $\mathcal{N}_{\rm c}\approx 1$. Instead for the signal emitted during the PPI, that we focus on here, we have to consider a period of time of $\approx 20$ orbits, which implies $\mathcal{N}_{\rm c}$ is higher. As for the GW frequency, we expect this to be the Keplerian frequency, as in the case of TDEs, but multiplied by the PPI azimuthal number $m$. However, since for our system $m=1$, the expected frequency is just the Keplerian one.

3.2 Amplitude of the maximum strain

If we consider the disruption of a solar type star by a non-rotating $10^{6}$M*☉* hole, we know from equation (6) that $\beta$ varies in the range

[TABLE]

Since, for a given star, the GW strain and frequency depend only on $\beta$, we have that these quantities range from $h\textcolor{black}{\approx}10^{-22}$ and $f\textcolor{black}{\approx}10^{-4}\,\text{Hz}$, for $\beta=1$, to $h\textcolor{black}{\approx}2\times 10^{-21}$ and $f\textcolor{black}{\approx}9\times 10^{-3}\,\text{Hz}$, for $\beta=20$. From the GW strain it is possible to derive $h_{\rm c}$ through equation (28), where the number of cycles is given by (see equation 29)

[TABLE]

We obtain this number considering that the torus simulated by Nealon et al. (2018) shows the PPI at a frequency which is the Keplerian frequency of the star tidally disrupted (so calculated at $r_{\rm p}$), while the torus is located around $2r_{\rm p}$. So $\tau$ is 20 times the period at $2r_{\rm p}$. The peak signal expected from this source with respect to the LISA sensitivity curve is illustrated in figure 1, where we have that the black curve is the sensitivity curve of the instrument, while the blue dots stand for the signal from a $1$M*☉* \textcolorblackTDE remnant. The parameter $\beta$ increases from left to right. We show with a red triangle the case simulated by Nealon et al. (2018), which has a characteristic strain and a frequency given by

[TABLE]

The signal goes above the sensitivity curve of LISA \textcolorblackfor $\beta\gtrsim 6$. \textcolorblackTo find this critical value of $\beta$, $\beta_{\rm cr}$, we infer from the plot the frequency at which the signal intersects the LISA sensitivity curve, and then we derive $\beta_{\rm cr}$ using equation  (27), considering $m_{*}=1$. Therefore, the case simulated by Nealon et al. (2018) appears not to be visible by the interferometer.

It is interesting to extend these calculations for a larger range of stellar masses. So, we investigate the signal from MS stars in the interval $1\leq m_{*}\leq 100$. However we have to consider, that, increasing $m_{*}$ also increases the minimum $\beta$ that is still in the frequency range visible to LISA\textcolorblack, since, as we can see from equation (27), we have

[TABLE]

\textcolor

blackAlso in figure 1, we have plotted the signal for \textcolorblack$m_{*}=10,100$ (grey and green respectively), where $\beta$ increases from left to right. For $m_{*}=10$, we have considered $5\leq\beta\leq 92$, while for $m_{*}=100$ we have $22\leq\beta\leq 430$. \textcolorblackIn particular, we have that the signal overtakes the instrument sensitivity curve for $\beta_{\rm cr}\approx 15$ if $m_{*}=10$, and $\beta_{\rm cr}\approx 45$ if $m_{*}=100$.

It is interesting to note that for each mass the signal which corresponds to the highest value of $\beta$ is at a frequency of $\approx 9\times 10^{-3}\,\text{Hz}$, which is indeed approximately the frequency of the innermost stable orbit for a hole of $10^{6}M_{\sun}$. \textcolorblackThis can be also seen from equation (27), from which we derive that the maximum frequency, $f_{\rm max}$, is independent on $m_{*}$ and is given by

[TABLE]

\textcolor

blackwhere $\beta_{\rm max}=r_{\rm t}/2r_{\rm gh}$.

From these simple, order of magnitude estimates of the process, it seems that the TDE of a solar type star by a static $10^{6}$M*☉* hole may be detectable by LISA, when we consider high values of $\beta$. However, we have to consider \textcolorblackthree important simplifications that probably make our calculations overestimate the real value:

not all of the stellar mass is involved in the PPI, i.e. the amplitude of the PPI perturbations is $<1$; 2. 2.

the torus is not located at $2r_{\rm p}$ but it is spreading out; 3. 3.

\textcolor

blackother effects, such as the MRI, would quench the PPI after a few orbits (see section 6), reducing the number of cycles $\mathcal{N}_{\rm c}$.

Thus, it is reasonable to expect that our analytical analysis overestimates the GW strain associated with the torus. So it makes sense to assume that the strain associated with the PPI, $h^{\rm PPI}$, can be written as

[TABLE]

where $\xi$ is a factor that varies in the range (0,1]. To quantify the effect of the assumptions listed above and find $\xi$ we have thus performed a numerical study described in Section 4.

4 Numerical study

We have repeated the simulation of Nealon et al. (2018): using the PHANTOM code, we simulate a $1$M*☉* torus of $N=5\times 10^{6}$ particles and we require the cross section parameter, $d$, to be

[TABLE]

where $r_{0}=0.5\,\text{AU}\textcolor{black}{\approx}2r_{\rm p}$, $r_{\rm p}$ being the pericenter of the Bonnerot et al. (2016)’s simulation. $r_{0}$ is the distance of the maximum density of the torus. First of all\textcolorblack, we relax the particles in order to reduce numerical artefacts due to the fact that they are set on a grid, then we set the torus in a fixed Keplerian potential and we add $m=1$ density perturbation to have the PPI. As in Nealon et al. (2018) we also only implement viscosity in order to capture shocks. The simulation runs for $\sim 20$ orbits, and each period $T$ is given by

[TABLE]

As in Nealon et al. (2018), during the first 3 orbits an over-density develops in the torus that orbits with the Keplerian frequency. The PPI continues to grow until 5th-6th orbits, when the over-density reaches its peak and a shock has developed, which spreads from the inner to the outer radii. Then the overdensity remains beyond this and the shock decreases.

Since we are interested in deriving $h_{+}$, $h_{\times}$ and $h$ associated with the system, we need to discretize equations (16) and  (17). We start by discretizing the momentum of inertia of the system, introduced with equation (15), as

[TABLE]

where $\rm a$ is the index that runs over the number of particles and $m_{\rm a}$ is the mass of the $\rm a$-th particle. Here we are only using the motion of the star material to calculate the GW signal (this is justified in appendix A). Since in the simulation we use same mass particles, we can write $m_{\rm a}=m$ and so we have

[TABLE]

Then we estimate the derivative of equation (39) numerically using central differencing

[TABLE]

where $j$ is the index related to time. Finally we substitute $\ddot{M}_{11}$, $\ddot{M}_{22}$ and $\ddot{M}_{12}$ into equations (16) and (17).

4.1 Numerical results

We numerically derive $\textcolor{black}{h^{\rm PPI}_{+}}$, $\textcolor{black}{h^{\rm PPI}_{\times}}$ from the SPH simulations for different directions, \textcolorblackwhere the superscript ‘PPI’ stands for the strain measured from the simulation. In particular in figure 2 we show $\textcolor{black}{h^{\rm PPI}_{+}}$ (left) and $\textcolor{black}{h^{\rm PPI}_{\times}}$ (right), calculated for the wave propagating in the $\boldsymbol{z}$ direction, that is the direction perpendicular to the stellar orbit. Both waveforms reach the peak between $\sim 1\times 10^{-3}\,\text{yr}$ and $\sim 2\times 10^{-3}\,\text{yr}$, that is around the fifth and sixth orbit, when the overdensity is stronger (see Nealon et al. 2018), and in particular the peak value for the $+$ polarization is lower than the one of the $\times$ polarization by a factor $\sim 2$.

In figure 3 we plot the strain of the wave $\textcolor{black}{h^{\rm PPI}}$ (see equation 19) with respect to the time $t$ in years, along the direction of propagation $\boldsymbol{z}$. It is interesting to note that the peak value is $\sim 4\times 10^{-24}$, which is two orders of magnitude lower than the expected analytical estimate of $\sim 5\times 10^{-22}$ for the raw strain, shown in equation (26). So for the torus formed after a TDE of a solar mass star, disrupted by a static $10^{6}\text{M}_{\sun}$ hole, with $\beta=5$, we have that $\xi\sim 10^{-2}$ and so

[TABLE]

If we assume that the same scaling factor holds also for different $\beta$ and for different stellar masses, we can extrapolate our results by shifting down by $10^{-2}$ the signals shown in figure 1. Thus, the expected signal would be visible only for $M_{*}\gtrapprox 10$M*☉* \textcolorblackand for values of $\beta$ higher than $\beta_{\rm cr}$ defined above.

4.2 Resolution test

We have run the same simulation with different resolutions. In particular we have \textcolorblackused $10^{5}$, $10^{6}$, $5\times 10^{6}$ and $9\times 10^{6}$ particles. For each simulation we plot the strain of the wave with respect to the time, as shown in figure 4. From this plot, it seems that, already with $5\times 10^{6}$ particles, there is convergence. In particular the peak of the strain is in similar position for the two highest resolutions. The curve corresponding to $9\times 10^{6}$ particles has only been run for long enough to demonstrate convergence. We thus regard our estimates from the simulation to be robust.

5 Discussion

According to our analytical study, the maximum peak value associated to 1M*☉* remnant disrupted by a static $10^{6}$M*☉* hole, with $\beta=5$, is just under the sensitivity curve of LISA. However, if we increase $\beta$, the signal becomes visible to the interferometer ($\beta\gtrsim 6$). Similar results can be obtained also for all the stars in the range $1\leq m_{*}\leq 100$, but, while increasing $m_{*}$, we have to check that the GW frequency remains in the range accessible to LISA.

The estimates made for the raw strain are similar to the ones in Kobayashi et al. (2004). However we need to consider that $h_{\rm c}$ is increased by a factor $\mathcal{N}_{\rm c}^{1/2}\approx 7.5$. In fact, since the GW signal of a TDE is a burst signal, we have that $\mathcal{N}_{\rm c}\approx 1$, which means $h_{\rm c}\approx h$. On the contrary, for the GW signal associated to the PPI in thick accretion discs, which is emitted on a longer time interval (in our system, 20 orbits), we have $h_{\rm c}\approx 7.5h$.

As for the numerical study, we show that the GW \textcolorblackstrain derived numerically from a 3D SPH simulation, performed in a classical frame, is two orders of magnitude lower than the analytical expectation

[TABLE]

This difference is due to the assumptions that we have made for our analytical analysis. First of all, the derivation of the strain is made for a tidally disrupted star, which is approximated by a point mass particle. Instead, the real system is made of a torus, that is spreading out due to the PPI. For this reason, considering the disc as a point where all the mass is located is not precise enough. Moreover, with our analytical study, we assume that all the total gas mass of the disc, since it is enclosed in the particle, undergoes the PPI, but in the real physical system in general only a part of the mass is involved. This is another feature that we are neglecting in the first part of our study.

It is interesting to compare the results we obtain numerically with the ones of Kiuchi et al. (2011). \textcolorblackIn contrast to our work, they show that their numerical results are one order of magnitude lower than the analytical estimates (figure 4 of Kiuchi et al. 2011), not two like in our case. However their study present some differences with respect to ours \textcolorblackthat may justify this discrepancy. First of all, they consider a disc much more massive, formed in a completely different scenario. In particular, the torus-SMBH mass ratio is bigger than in our case, which implies that it is reasonable to expect a significant contribution to the GW signal also from the BH. Then, they consider a self-gravitating torus. Self gravity is not a requirement of the PPI, since this instability can occur in any torus that is not accreting. However, self-gravity \textcolorblackplays an important role in the formation of clumps of matter inside the torus, that are the sources of the GWs. Moreover, we also need to consider that for sources with a non negligible self-gravity, it is not clear that the derivation made in paragraph 2.3 still holds (see Buonanno 2007). For detailed discussions about the role of self-gravity in relativistic discs see Korobkin et al. (2011) and Mewes et al. (2016). Finally they use a grid code with GR. Kobayashi et al. (2004) have shown that for TDEs with low $\beta$, GR does not change the estimates of GW emission in a significant way. However, this could be different in the case of accretion discs, which are very massive (as in Kiuchi et al. 2011) and for this reason have a non-negligible interaction with the SMBH.

6 Conclusions

In this work we have studied the GW signal associated with the PPI in a thick disc, with a shallow specific angular momentum profile and an inner and outer radii well defined. This disc has just resulted from the TDE of a $1$M*⊙* star around a non rotating SMBH of $10^{6}$M*⊙*, with $\beta=5$.

First of all we have made some analytical estimates of the maximum amplitude, to see if the signal could be above the LISA sensitivity curve. From this study we have found that the signal might be visible for discs with masses in the range $[1$M${}_{\sun},100$M${}_{\sun}]$, for high values of $\beta$, \textcolorblackin particular $\beta\gtrsim 6\,\text{for}\,1\text{M}_{\sun}$, $\beta\gtrsim 15\,\text{for}\,10\text{M}_{\sun}$, and $\beta\gtrsim 45\,\text{for}\,100\text{M}_{\sun}$.

Then we have performed a numerical 3D SPH simulation, using the \textcolorblackPHANTOM code, and compared the numerical results with the previous estimates. We have found that the numerical maximum strain is two orders of magnitude lower than the analytical one. This lowering can be justified considering that, during our analytical analysis, we neglect the disc is spreading out due to the PPI, and that not all the disc mass is involved in the instability.

\textcolorblackIt is important to remember that in our simulation we do not consider magnetic fields. The inclusion of magnetic fields from the beginning could suppress the GW emission or even avoid it, according to the particular structure of the disc after circularization. In this regard, Bugli et al. (2018) have illustrated that the MRI can allow an initial formation of the PPI, but the disc needs to already show a dominant $m=1$ mode. The final $m=1$ mode present in Bonnerot et al. (2016), that justifies the perturbation in our model as in Nealon et al. (2018), may not occur if magnetic fields are included a priori.

In conclusion, our numerical study suggests that this source could still be detectable by LISA, if we consider tori with masses in the range $[10$M${}_{\sun},100$M${}_{\sun}]$, resulting after \textcolorblackdeeply penetrating TDEs.

\textcolor

black

Acknowledgements

MT and GL have received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement NO 823823 (RISE DUSTBUSTERS project). RN has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 681601).

Appendix A The contribution of the black hole to the gravitational wave signal

In equation (38), we calculate the GW signal of the system only taking into account the contribution of the torus and not of the hole. This can be explained in the following ways.

First of all we can approximate the system SMBH-torus as a binary system where we have that the total mass of the system, $M=M_{\rm h}+M_{\rm d}$, where $M_{\rm d}$ is the mass of the disc, can be approximated with the mass of the black hole, since we have that $M_{\rm h}$ is between $10^{6}$ and $10^{4}$ times the mass of the torus, if we consider tori with masses in the range $[1\text{M}_{\sun},100\text{M}_{\sun}]$. The reduced mass of the system, defined as

[TABLE]

is approximately equal to the mass of the disc, for the same reasons above. So it is reasonable to say that the position of the centre of mass of the system corresponds to the position of the black hole, since the black hole is much more massive than the torus. These considerations suggest that the mass of the black hole is not a significant factor in the amplitude of the GWs, and as such we do not consider it in our derivation of equation (38).

A more formal way to show that the SMBH does not contribute to the GW signal emission in our study, is the the following. In the centre of mass frame, we have that

[TABLE]

where $r_{1}$ and $r_{2}$ are the displacements from the centre of mass of the hole and the disc, respectively. If we want to derive the moment of inertia of the black hole we have

[TABLE]

In particular, if we assume $M_{\rm h}=10^{6}\text{M}_{\sun}$ and $M_{\rm d}=1\text{M}_{\sun}$, we have that the moment of inertia of the hole becomes

[TABLE]

while if $M_{\rm d}=100$M*☉* we have

[TABLE]

The above illustrates that the GW signal emitted by the SMBH is negligible with respect to that emitted by the material in the disc. Hence when we calculate the GW signal, we only consider the mass in the disc.