Numerical-relativity simulation for tidal disruption of white dwarfs by a supermassive black hole

We study tidal disruption of white dwarfs in elliptic orbits with the eccenticity of $\sim 1/3$--$2/3$ by a non-spinning supermassive black hole of mass $M_{\rm BH}=10^5M_\odot$ in fully general relativistic simulations targeting the extreme mass-ratio inspiral leading eventually to tidal disruption. Numerical-relativity simulations are performed by employing a suitable formulation in which the weak self-gravity of white dwarfs is accurately solved. We reconfirm that tidal disruption occurs for white dwarfs of the typical mass of $\sim 0.6M_\odot$ and radius $\approx 1.2 \times 10^4$\,km near the marginally bound orbit around a non-spinning black hole with $M_{\rm BH}\alt 4\times 10^5M_\odot$.


I. INTRODUCTION
Tidal disruption of ordinary stars and/or white dwarfs by supermassive black holes has been revealed to be one of the major sources of bright electromagnetic transients (see, e.g., Refs. [1-3]), which have been actively observed in the last decade. In addition, gravitational waves emitted by tidal disruption of white dwarfs closely orbiting supermassive black holes could be observable by Laser Interferometer Space Antenna (LISA) [4]. Electromagnetic signals associated with tidal excitation (e.g., Ref. [5]) or mass stripping (e.g., Refs. [6][7][8][9] for related works) or tidal disruption (e.g., Refs. [10,11]) of white dwarfs can be an important electromagnetic counterpart of gravitational waves. Because the expected event rate is not so high [12] that the signal-to-noise ratio of gravitational waves for the LISA sensitivity is unlikely to be very high, the discovery of the possible electromagnetic counterparts will help extract gravitational waves from the noisy data in the LISA mission.
The condition for mass shedding and tidal disruption during the cross encounter of stars with supermassive black holes is often described by the so-called β -parameter defined by where r p is the periastron radius for the orbit and r t is the Hill's radius [13] defined by with R * the stellar radius, M * the stellar mass, and M BH the mass of the supermassive black hole, respectively. Since the early 1990s (see, e.g., Refs. [14][15][16]), a large number of numerical simulations have been performed in the last three decades (see, e.g., Refs. [17,18] for reviews of the latest works and Refs. [19][20][21][22][23][24] for some of the most advanced works). They have shown that mass stripping can take place at the close encounter if β is larger than about 0.5, and tidal disruption can take place if β 1 for stars in parabolic orbits (see, e.g., Refs. [17,25,26] for Newtonian simulation works, and also early semianalytical work [27,28]). It is also shown that for close orbits around a black hole, the general relativistic effect can significantly reduce the critical value of β for the tidal disruption [23]. Indeed, general relativistic works show that for circular orbits near the innermost stable circular orbit of black holes, the mass shedding can occur even for β ∼ 0.4 [29,30]. However, the previous analyses have been carried out in Newtonian gravity or in relativistic gravity of a black hole with Newtonian (or no) gravity for the companion star or in a tidal approximation with a relativistic tidal potential [29][30][31]. To date, no fully general relativistic (the so-called numericalrelativity) simulation, i.e., a simulation with no approximation except for the finite differencing, has been done for the tidal disruption problem with β 1 (but see Refs. [32][33][34] for a head-on and an off-axis collision).
Numerical-relativity simulation is suitable for the tidal disruption problem for the case that the orbit at the tidal disruption is highly general-relativistic. This is particularly the case for tidal disruption of white dwarfs by supermassive black holes because it can occur only for orbits very close to the black-hole horizon. Advantages of the numerical-relativity simulation are: (i) the redistribution of the energy and angular momentum of the star can be followed in a straightforward manner and (ii) we can directly follow the matter motion after the tidal disruption including the subsequent disk formation.
In this paper, we present a result of numerical-relativity simulations for tidal disruption of white dwarfs of typical mass (0.6-0.8M ) by a supermassive black hole with relatively low mass (M BH = 10 5 M ) for the first time. For simplicity, the white dwarfs are modeled by the Γ = 5/3 polytropic equation of state. As a first step toward more detailed and systematic studies, we focus on tidal disruption of white dwarfs in mildly elliptic orbits aiming at confirming that our numerical-relativity approach is suitable for reproducing the criteria of tidal disruption, which has been already investigated in many previous works referred to above.
The paper is organized as follows. In Sec. II, we describe our formulation for evolving gravitational fields, matter fields, and for providing initial data of a star in elliptic orbits around supermassive black holes. In Sec. III, numerical results are presented paying particular attention to the criterion for tidal disruption. Section IV is devoted to a summary. Throughout this paper, we use the geometrical units of c = 1 = G where c and G are the speed of light and gravitational constant, respectively. The Latin and Greek indices denote the space and spacetime components, respectively.

II. BASIC EQUATIONS FOR THE TIME EVOLUTION
A. Gravitational field First, we reformulate the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism [35,36] in numerical relativity to a form suitable for the simulation of high-mass ratio binaries, in particular for accurately computing a weak self-gravity of white dwarfs. Throughout this paper, high-mass ratio binaries imply those composed of a very massive black hole of mass M BH 10 5 M and a white dwarf (or an ordinary star) of mass of M * = O(M ) with the radius R * 10 3 km, for which the compactness defined by M * /R * is smaller than 10 −3 .
We consider the two-body problem with a compact orbit of the orbital separation r 30M BH . With such setting, the magnitude of the gravitational field generated by the black hole, which is defined by g µν − η µν is, of order M BH /r > 10 −2 .
Here g µν and η µν are the spacetime metric and Minkowski metric, respectively. On the other hand, the magnitude of the gravitational field generated by white dwarfs and ordinary stars is of order M * /R * < 10 −3 , which is much smaller than that by the black hole. To accurately preserve the nearly equilibrium state of such stars during their orbits, an accurate computation of the gravitational field by them is required. However, if we simply solve Einstein's equation, a numerical error for the computation of the black-hole gravitational field can significantly affect the gravitational field for the white dwarfs/ordinary stars. To avoid this numerical problem, we separate out the gravitational field into the black hole part and other part, although we still solve fully nonlinear equations. The idea employed here is similar to that of Ref. [33], but we develop a formalism based on the BSSN formalism.
In a version of the BSSN formalism [36], the basic equations are written in the form: where α is the lapse function, β j is the shift vector,γ i j is the conformal three metric defined from the three metric γ i j bỹ γ i j := γ −1/3 γ i j with γ = det(γ i j ), W := γ −1/6 ,Ã i j is the conformal trace-free extrinsic curvature defined from the extrinsic curvature K i j byÃ i j = W 2 (K i j − γ i j K/3) with K := K k k , andΓ i := −∂ jγ i j . ρ h , J i , and S i j are quantities defined from the energy-momentum tensor, T µν , by ρ h = T µν n µ n ν , J i = −T µν n µ γ νi , and S i j = T µν γ µi γ ν j with n µ the timelike unit vector normal to spatial hypersurfaces.
In this problem, we employ the so-called puncture gauge [37], in which the evolution equations for α and β i are written as where B i is an auxiliary three-component variable and η B is a constant of order M −1 BH . By introducing a static black-hole solution for the geometric variables, α 0 , β i 0 ,γ 0 i j , W 0 ,Ã 0 i j , and K 0 and by writing all the variables by we then write down the equations for α s , β i s ,γ s i j , W s ,Ã s i j , K s , andΓ i s (these are denoted by a representative variable Q s as follows). Specifically, the evolution equations (4)-(7) and (8)-(10) of the geometrical variables (denoted by a representative variable Q) are schematically written in the form Then, for the decomposition of Q = Q 0 + Q s with F(Q 0 ) = 0 (under the conditions of ∂ t Q 0 = 0), we write the equation for Q s as In numerical simulation, F(Q 0 ) obtained from finite difference is nonzero which contains the truncation error of evolving the stationary background metric numerically. Here, we added the second-term in the right-hand side to explicitly subtract the leading error of evolving the background metric so that the right-hand side of the evolution equation of Q s does not have the zeroth order terms in Q s Any static black-hole solutions can be used for α 0 , β i 0 , · · · , but in the BSSN formalism with the puncture gauge, the metric relaxes to a solution in the limit hypersurface with K 0 = 0. Using such a trumpet-puncture black hole also allows us to construct the initial data in the conformal-thin-sandwich (CTS) formalism [38] (see Sec. II C). Thus, in the present formalism, it is appropriate to employ such a solution. In the nonspinning black hole, the analytic solution is known and is written as [39] and K 0 = 0 where R is a function of r determined by [40] We note that r = 0 corresponds to R = 3M BH /2 and the event horizon is located at R = 2M BH (i.e., r ≈ 0.78M BH ) in this solution.

B. Hydrodynamics
In this paper we model white dwarfs simply by the polytropic equation of state, where P and ρ are the pressure and rest-mass density, respectively, κ the polytropic constant, and Γ adiabatic index for which we set to be 5/3. For the hydrodynamics, we solve the continuity and Euler equations, with ∇ µ the covariant derivative with respect to g µν and where ε and u µ are the specific internal energy and four velocity, respectively. In this work we do not solve the energy equation, and determine ε simply by ε = κρ Γ−1 /(Γ − 1) which is derived from the condition that the specific entropy is conserved for the fluid elements. The continuity and Euler equations are solved in the same scheme as that used in Refs. [41,42]. The motivation for using the polytropic equation of state comes from the fact that our primary purpose of this paper is to explore the tidal disruption condition for a relatively low value of β < 1 and the formation of shocks by the tidal compression does not play any role. We here focus only on the process of tidal disruption and subsequent short-term evolution of the tidally disrupted material. After the tidal disruption, the fluid is highly elongated and during the long-term evolution of the fluid elements with different specific energy and angular momentum, they collide and shocks are likely to be formed. For such a phase, the shock heating will play an important role. Our plan is to follow this phase by solving the energy equation with a more general equation of state.

C. Initial condition
First, we describe the formulation employed in this paper for computing the initial data in which white dwarfs are approximately in an equilibrium state in their comoving frame. From Eq. (28), we have where h is the specific enthalpy defined by h := 1 + ε + P/ρ. To derive Eq. (30), we used Eq. (27). For the isentropic fluid, the first law of thermodynamics is written as where dQ denotes the variation of a quantity Q in the fluid rest frame. In the polytropic equations of state employed in this work, we obtain the relation In this situation, Eq. (30) is rewritten to Then, we define k µ := u µ /u t . Using this quantity, Eq. (33) is written to where L k denotes the Lie derivative with respect to k µ . The second term of Eq. (34) is written as where we used u µ u µ = −1. Thus, Eq. (34) is written to We consider an initial condition for a system composed of a star of mass M * and radius R * , for which the center is located on the x-axis, around a massive black hole of mass M BH M * and M BH R * which is located at a coordinate origin. We assume that the star predominantly moves toward the y-direction with the identical specific momentum. Thus being a constant to be determined. Here the term of β i 0 is added to simplify the iteration process for computing quasiequilibrium states. Then, u t is calculated from In the present context, L k (hu i ) can be assumed to be zero for i = y and z, because the star has momentarily translation invariance for the motion toward the yand z-directions. By contrast, with respect to the x-direction, the star receives the force from the massive black hole. Since the radius of the star, R * , is much smaller than the orbital separation, x 0 , and x 0 is larger than the black-hole radius of where C is an integration constant. We note that Eq. (38) is not an exact first integral of the Euler equation but can be considered as an approximate one for obtaining an initial condition in which the star is in an approximate equilibrium state. For computing initial conditions, we assume the line elements of the form where ψ is the conformal factor. Using the Isenberg-Wilson-Mathews formalism [43,44], the basic equations are written as where and ∆ is the flat Laplacian.Â i j is defined from the extrinsic curvature,Ã i j , byÂ i j = ψ 6Ã i j and K is set to be zero. Using the CTS decomposition [45,46] with trumpet-puncturê Equation (42) is rewritten as Note that although there are some works in constructing binary black holes initial data with trumpet-puncture [47][48][49], this is, to our knowledge, the first attempt combining the CTS decomposition and puncture method with the limit (trumpet) hypersurface in constructing quasi-equilibrium initial data in nonvacuum spacetime. We assume that the contribution to the extrinsic curvature from the black hole is negligible because the orbital momentum of the black hole is negligible in this problem, and thus, we set the black hole at rest (however, it is straightforward to take into account the small black-hole motion [50] in our formalism.) For a solution of the initial data, we have to determine the free parameters, A, C, and V . In the polytropic equation of state, we can consider κ as well as the central density ρ c as free parameters. In the following, we first consider that V and rest mass of the star are input parameters and A, C, and κ are parameters to be determined during the iteration process in numerical computation. Our method to adjust κ to a desired value will be described later.
To determine these three parameters we need three conditions, for which we choose the following relations. First, we fix the location of the surface of white dwarfs along the x-axis as x = x 1 (referred to as point 1) and x = x 2 (point 2). Typically, we choose x 1 + x 2 = 2x 0 . At the surface, h = 1, and thus, Eq. (38) gives where u t 1 and u t 2 are the values of u t at points 1 and 2. In addition, we fix the rest mass of the star which is defined by where m * is approximately equal to the gravitational mass M * because the star is only weakly self-gravitating. Using the condition (48), the values of C and A are determined, and subsequently, h is determined from Eq. (38). In the polytropic equation of state, the rest-mass density is written as and thus, from Eq. (49), κ is determined for given values of m * and x 2 − x 1 . Once these free parameters are determined, the rest-mass density are obtained from Eq. (50). For realistic setting, we have to obtain the desired values of the mass of the star and the value of κ. The value of κ is controlled by varying the stellar diameter x 2 − x 1 for a given value of m * .
To take into account the effect of the black-hole gravity, we employ the puncture formulation by setting where ψ 0 , α 0 , β k 0 , andÂ 0 i j denote the solutions of vacuum Einstein's equation shown already in Sec. II A. Then we numerically solve the equations for φ , X, β k s , andÂ s i j from Eqs. (40), (41), (47), and (46). The initial data is prepared using the octree-mg code [51], an open source multigrid library with an octree adaptive-mesh refinement (AMR) grid, which we modified to support a fourth-order finite-difference elliptic solver.

A. Setup
The simulation is performed using an AMR algorithm with the equatorial symmetry imposed on the z = 0 (equatorial) plane using the SACRA-TD code (for SACRA see Refs. [41,42]). We prepare two sets of finer domains, one of which comoves with a white dwarf and the other of which is located around the center and covers the massive black hole. Because the radius of the white dwarf, R * , is smaller than the black-hole horizon radius ∼ M BH , we need to prepare more domains for resolving the white dwarf. In addition to these domains, we prepare coarser domains that contain both the finer domains in their inside. All the domains are covered by (2N + 1, 2N + 1, N + 1) grid points for (x, y, z) with N being an even number.
Specifically, each domain is labeled by i which runs as 0, 1, 2, · · · , i fix , · · · , i BH , · · · , i max . The grid resolution for the domains with i fix ≤ i ≤ i BH is identical with that with i BH +1 ≤ i ≤ 2i BH − i fix + 1(< i max ), respectively. For 0 ≤ i ≤ i BH , the center of the domain is located at the origin, at which a black hole is present. Strictly speaking, the black hole moves due to the backreaction against the motion of the companion star, but this motion is tiny because of the condition M BH M * . For these domains, the ith level covers a half cubic region of is the grid spacing for the ith level, and the grid spacing for each level is determined by ∆x i+1 = ∆x i /2 (i = 0, 1, 2, · · · , i BH − 1 and i = i BH + 1, · · · , i max − 1) with ∆x i BH +1 = ∆x i fix and L i BH ∼ 0.8M BH .
For the moving domains that cover the white dwarf, the center is chosen to approximately agree with the location of the density maximum. In the present context, the local density maximum is approximately located along a geodesic around the supermassive black hole. The size of the finest domain with i = i max , L max , is chosen so that it is 1.3-1.5R * . We check the convergence of two different models with three grid resolutions as illustrated in Fig. 1. Higher resolution is used for model M8V17 to measure the spin up of the white dwarf more accurately (see Sec. III B). We obtain good convergence for TABLE I. Models considered in this paper and the fate (last column). M7V16a and M7V16b correspond to the models with R * = 8.5 × 10 3 and 7.0 × 10 3 km, respectively. For other models, R * ≈ 10 4 (M * /0.7M ) −1/3 km. r p and r p,A are periastron radius in the present coordinates and the Schwarzschild coordinates, respectively. TD and OC denote tidal disruption and appreciable oscillation of white dwarfs, and NN denotes that no appreciable tidal effect is found. both models, and thus, we employ N = 60 as the standard resolution in this paper.

B. Numerical results
In the present paper we focus on the case that the black-hole mass is M BH = 10 5 M , the white-dwarf mass is M * = 0.6, 0. cases where κ is chosen such that R * = 8.5×10 3 km and R * = 7.0 × 10 3 km.
The initial separation is set to be x 0 = 20M BH (it is ≈ 21.01M BH in the Schwarzschild coordinates), and V is chosen to be 0.160, 0.165, 0.170, 0.175, and 0.180 (see Table I With these settings, the white dwarf has an elliptic orbit around the black hole with the periastron at r p ≈ (4.4-10)M BH , and thus, the eccentricity is approximately defined by e = (x 0 − r p )/(x 0 + r p ) is ≈ 1/3-2/3. Here, x 0 (= 20M BH ) and r p are defined in the radial coordinates of the metric of g 0 µν , and thus, the values of e slightly change if we define it in the areal coordinate (Schwarzschild radial coordinate).
For the models mentioned above, the value of β is in the range between 0.33 and 0.72 and estimated by we find 0.50 ≤ β ≤ 0.7, and thus, the white dwarf is expected to be strongly perturbed by the black-hole tidal field for M * = 0.6-0.8M . By contrast, for V = 0.180, β < 0.35 with M * = 0.7M , and thus, the tidal force of the black hole is likely to be too weak to perturb the white dwarf.
For V = 0.170, β ≈ 0.49, 0.44, and 0.40 with M * = 0.6, 0.7, and 0.8M respectively. In these cases, tidal disruption is not very likely to take place but the tidal force from the black hole should induce the stellar oscillation on the white dwarf. Because for Γ = 5/3, the stellar radius depends only weakly on the stellar mass, the presence or absence of the tidal disruption is likely to depend primarily on the value of V (or the specific angular momentum of the white dwarfs) in the present setting. In the following, we will show that our code can reproduce all these expected phenomena. Figure 3 plots the evolution of the maximum density for V = 0.160, 0.165, 0.170, and 0.180 with M * = 0.7M . We note that for M BH = 10 5 M , the orbital period for these parameters are in the range from ≈ 220 s for V = 0.160 to ≈ 250 s for V = 0.180. The figure shows the results expected in the previous paragraphs: For M7V16 and M7V165, the white dwarfs are tidally disrupted while approaching the black hole irrespective of the white-dwarf mass. For M7V17 (β = 0.44), the white dwarf is perturbed by the black hole near the periastron but it is not tidally disrupted. After the close encounter, the white dwarf is in an oscillating state due to the instantaneous tidal force received from the black hole. By contrast, for M7V18, the maximum density is approximately preserved to be constant, suggesting no disruption occurs and the tidal effect is negligible. Note that such tidal field may still perturb the white dwarf and produce detectable electromagnetic or gravitational-wave signal if a sufficient amplitude of oscillation is induced.
In Fig. 3, the results of M7V16 (β = 0.65), M7V16a (β = 0.55) and M7V16b (β = 0.45) are also compared. As expected, for the first two models, the white dwarfs are tidally disrupted, while for the most compact white dwarf, the tidal disruption does not occur although it is perturbed significantly by the black-hole tidal force. This illustrates that the β parameter is a good indicator for assessing whether tidal disruption takes place or not irrespective of the white-dwarf radius. Figure 4 shows the evolution of the maximum density when stellar oscillation is induced. For M6V17 (β = 0.49), the white dwarf is significantly elongated by the tidal force from the black hole; the central density is decreased to less than 50% of the original value after passing through the periastron. Associated with the tidal effect, the mass is lost from the white dwarf. However, with the increase of the orbital radius, the central density increases again, resulting in a less massive white dwarf. This is also the case for M8V165 (β = 0.47) and M7V16b (β = 0.45). These results indicate that the critical value of β for the tidal disruption is ∼ 0.50 and the threshold value for exciting a high-amplitude oscillation is β ∼ 0.45. Figure 4 also shows that even for 0.40 β 0.45 an appreciable oscillation is excited by the tidal force.
In Fig. 5 and Table I, we summarize the fates of white dwarfs as a result of the tidal interaction. It is found that for β 0.5 tidal disruption takes place and for β 0.4, the white dwarfs are perturbed appreciably by the black-hole tidal field. All these results agree approximately with the expectation from the previous studies.
For M7V16, tidal disruption takes place but only a small fraction of the white dwarf matter falls into the black hole because the fluid elements have specific angular momentum large enough to escape capturing by the black hole. Most of the tidally disrupted matter approximately maintains the original elliptic orbit (see Fig. 6) although the matter has an elongated profile. To clarify the eventual matter distribution around the black hole, we will need to follow the matter motion for more than 10 orbits. This topic is one of our major research targets in the future.
For 0.4 β 0.5, the white dwarf will be continuously perturbed by the black-hole tidal force whenever it passes through the periastron. In addition the angular momentum is transported during the tidal interaction, and it will lead to the transport of the orbital angular momentum to the white dwarf resulting in a spin-up of it. According to a perturbation study for the stellar encounter, the energy deposition during the tidal interaction in one orbit is written approximately as [52] ∆E where f tid is a factor of O(0.1), which depends on β and the equation of state. Associated with the energy deposition near the periastron, the angular momentum deposition is also deposited. In one orbit it is approximately estimated by ∆J spin ≈ ∆E tid /Ω p [53] where Ω p = M BH /r 3 p,A , and thus, where f spin is a coefficient of the same order of the magnitude of f tid . Because the maximum spin angular momentum of the star is approximately written as M * √ M * R * , we find that ∆J spin can be more than 0.1% of the maximum spin if a white dwarf passes through a close orbit with β 0.4. We approximate the orbit angular momentum J orbit and the spin angular momentum J spin of the white dwarf as where M h := d 3 xψ 6 ρh and the volume average of quantity q is defined as q := 1 M h d 3 xψ 6 ρhq. In such decomposition, the sum of orbital and spin angular momentum equals to the total angular momentum of the white dwarfs. We analyzed the spin angular momentum gain of the white dwarfs for M8V17, and we indeed find ∆J spin /(M * √ M * R * β 9/2 ) ≈ 0.1-0.3 as shown in Fig. 7. Note that the spin up of white dwarf ∆J spin is about 10 −6 of the total angular momentum, and hence, it is not easy to determine ∆J spin accurately. Although we cannot achieve a good convergence in ∆J spin , we are able to obtain a noticeable rise in J spin during the close encounter, which suggests f spin ∼ 0.1-0.3, which is consistent with the above analytic result. For close orbits, the tidal angular-momentum transport can dominate over the orbital angular momentum loss by gravitational-radiation reaction. Assuming that gravitational waves are most efficiently emitted near the periastron at which we may approximate the orbit to be circular, the angular momentum dissipation by gravitational waves in one orbit can be written as [54] ∆J GW ≈ 64π 5 where e denotes the eccentricity. Thus, the ratio of ∆J tid to ∆J GW is written as Thus it is larger than unity for r p,A 7M BH /c 2 , R * ≈ 10 4 km, M BH = 10 5 M , M * = 0.7M , and f spin = 0.2. This is also the case for the ratio of ∆E tid /∆E GW where ∆E GW is the energy dissipated by gravitational waves in one orbit. Thus, near the tidal disruption orbit, the orbital evolution would be primarily determined not by the gravitational-wave emission but by the tidal effect. To clarify the eventual fate of such a white dwarf, we obviously need a long-term accurate simulation. Such a topic is one of our future targets. We note that both ∆J spin and ∆J GW are much smaller than the orbital angular momentum of order M * M BH r p,A . Thus, the cumulative effect of the tidal angular momentum transport plays an important role just prior to the tidal disruption. By repeated tidal interaction, the spin angular velocity of the white dwarfs is likely to be enhanced up to ∼ M 1/2 In addition, the stellar oscillation for which the oscillation energy is comparable to or larger than the rotational kinetic energy should be excited. As a result, mass loss could be induced, resulting in the increase of the stellar radius and enhancing the importance of tidal interaction. In this type of the system, the tidal disruption is unlikely to take place by one strong impact by the black-hole tidal force but likely to do as a result of a secular increase of the stellar radius (see, e.g., Refs. [9,55] for related studies).

IV. SUMMARY
We reported a new numerical-relativity code which enables us to explore tidal disruption of white dwarfs by a relatively low-mass supermassive black hole. As a first step toward more detailed future studies, we paid attention to the condition for tidal disruption of white dwarfs with typical mass range in elliptic orbits by a nonspinning supermassive black hole. We showed that our code is capable of determining the condition for the tidal disruption. As expected from previous general relativistic works (e.g., Refs. [23,30]), the tidal disruption takes place for β 0.5 and an appreciable oscillation of the white dwarfs are induced by the black-hole tidal effect for β 0.4 for orbits close to the black hole in the Γ = 5/3 polytropic equation of state. The critical value for the onset of the tidal disruption is smaller than that obtained by Newtonian analysis. For white dwarfs with M * = 0.6M and R * = 1.2 × 10 4 km, β can be larger than 0.4 even for M BH ≈ 4 × 10 5 M if the periastron radius is r p,A = 4M BH . Our result indicates that in such systems with a relatively lowmass (but not intermediate-mass) supermassive black hole for which gravitational waves in the late inspiral phase can be detected by LISA [4], tidal disruption can occur for typical white dwarfs. For spinning black holes with the dimensionless spin parameter of 0.9, r p,A can be smaller than ∼ 1.7M BH [56]. For such black holes, tidal disruption of typical-mass white dwarfs may occur even for M BH ≈ 10 6 M . Investigation of this possibility is a future issue.
There are several issues to be explored. The first one is to extend our implementation for spinning black holes. Since no analytic solution is known for the spacetime of spinning black holes on the limit hypersurface, we need to develop a method to provide g 0 µν for employing the formulation introduced in this paper. One straightforward way to prepare such data is just to numerically perform a simulation for a spinning black hole (in vacuum) until the hypersurface reaches the limit hypersurface as a first step, and then, the obtained data are saved and used in the subsequent simulations with white dwarfs. A more subtle issue along this line is to prepare the initial condition. For nonspinning black holes, we can assume that the conformal flatness of the three metric, and as a result, the initial-value equations are composed only of elliptic-type equations with the flat Laplacian. For the spinning black holes, the basic equations are composed of elliptictype equations of complicated Laplacian, and hence, the numerical computation could be more demanding, although in principle it would be still possible to obtain an initial condition. We plan to explore this strategy in the subsequent work.
For modeling realistic white dwarfs it is necessary to implement a realistic equation of state. If we assume that the temperature of the white dwarfs is sufficiently low and the pressure is dominated by that of degenerate electrons, it is straightforward to implement this.
More challenging issue is to follow the hydrodynamics of tidally disrupted white-dwarf matter for a long term. After the tidal disruption, the matter of the white dwarf is likely to move around the black hole for many orbits. During such orbits, the matter collides each other, and eventually, a compact disk will be formed after the circularization. Such disks are likely to be hot due to the shock heating, and thus, it can be a source of electromagnetic counterparts of the tidal disruption. In the presence of magnetic fields, magnetorotational instability [57] occurs in the disk, and the magnetic fields will be amplified. If the amplified magnetic field eventually penetrates the black hole and if the black hole is appreciably spinning, a jet may be launched through the Blandford-Znajek effect [58]. After the amplification of the magnetic fields, a turbulent state will be developed in the disk and mass ejection could occur by the effective viscosity or magneto-centrifugal force [59]. The ejecta may be a source of electromagnetic signals. One long-term issue is to investigate such scenarios by general relativistic magnetohydrodynamics.