Does Explosive Nuclear Burning occur in Tidal Disruption Events of White Dwarfs by Intermediate Mass Black Holes ?

TL;DR

This study uses high-resolution 3D and 1D simulations to explore whether nuclear detonations occur during white dwarf tidal disruptions by intermediate mass black holes, finding no shock-induced detonations at current resolutions.

Contribution

The paper provides unprecedented high-resolution simulations showing the absence of shock waves and detonations in white dwarf tidal disruption events, challenging previous assumptions.

Findings

Nuclear reactions decrease with higher resolution due to spurious heating.

No shock wave generation observed in 3D simulations at high resolution.

Shock waves and potential detonations found in 1D simulations at high resolution.

Abstract

We investigate nucleosynthesis in tidal disruption events (TDEs) of white dwarfs (WDs) by intermediate mass black holes (IMBHs). We consider various types of WDs with different masses and compositions by means of 3 dimensional (3D) smoothed particle hydrodynamics (SPH) simulations. We model these WDs with different numbers of SPH particles, $N$, from a few $10^4$ to a few $10^7$, in order to check mass resolution convergence, where SPH simulations with $N>10^7$ (or a space resolution of several $10^6$ cm) have unprecedentedly high resolution in this kind of simulations. We find that nuclear reactions become less active with increasing $N$, and that these nuclear reactions are excited by spurious heating due to low resolution. Moreover, we find no shock wave generation. In order to investigate the reason for the absence of a shock wave, we additionally perform 1 dimensional (1D) SPH and mesh-based simulations with a space resolution ranging from $10^4$ to $10^7$ cm, using characteristic flow structure extracted from the 3D SPH simulations. We find shock waves in these 1D high-resolution simulations. One of these shock waves triggers a detonation wave. However, we have to be careful of the fact that, if the shock wave emerged at a bit outer region, it could not trigger the detonation wave due to low density. Note that the 1D initial conditions lack accuracy to precisely determine where a shock wave emerges. We need to perform 3D simulations with $\lesssim 10^6$ cm space resolution in order to conclude that WD TDEs become optical transients powered by radioactive nuclei.

Tables

Table 1. Table 1: Summary of initial conditions.

Model	$M_{\text{wd}}$	$R_{\text{wd}}$	$M_{\text{bh}}$	$\beta$	$R_{\text{p}}$	Compositions	CC	IMBH	$N$	Comments
ONeMg	$1.2$	$0.35$	$100$	$3.0$	$0.52$	16O $60$% 20Ne $35$% 24Mg $5$%	w/	PW	up to $13$M
CO1	$0.6$	$0.75$	$500$	$5.0$	$1.4$	12C $50$% 16O $50$%	w/o	NP	up to $6.3$M	Rosswog’s run 8
CO2	$0.6$	$0.75$	$500$	$5.0$	$1.4$	12C $50$% 16O $50$%	w/	TR	up to $6.3$M
CO3	$0.6$	$0.75$	$500$	$5.0$	$1.4$	12C $50$% 16O $50$%	w/	PW	up to $25$M
He	$0.3$	$1.00$	$500$	$5.0$	$2.4$	4He $100$%	w/	PW	up to $25$M

Full text

hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: general — black hole physics — white dwarfs

1 Introduction

The number of candidates for tidal disruption events (TDEs), in which stars are tidally disrupted by a massive black hole, has been rapidly increasing (e.g. Komossa, 2015). Various stellar types can be considered, including a main sequence star, a red giant star, and a white dwarf (WD). So far, several high energy transients have been proposed to be TDEs of WDs (WD TDEs) (Krolik & Piran, 2011; Shcherbakov et al., 2013; Jonker et al., 2013).

Among various implications, finding WD TDEs is of special importance to address the existence of an intermediate mass black hole (IMBH) – a black hole (BH) disrupting a WD can be an IMBH. 111WDs can also be tidally destroyed by a stellar-mass BH. However, this is beyond the scope of this paper. We will investigate it elsewhere. This is because a supermassive BH (SMBH) with more than $10^{5}M_{\odot}$ just swallows a WD rather than disrupts it (Luminet & Pichon, 1989; Kobayashi et al., 2004). WD TDEs can thus be an important probe to IMBHs, since a large number of WD TDEs may be detected with the aid of current transient surveys (e.g., intermediate Palomar Transient Factory) and next-generation transient surveys (e.g. the Zwicky Transient Facility and the Large Synoptic Survey Telescope). In the future, the detections of WD TDEs could constrain the abundance of IMBHs in the universe, although just a few IMBH candidates have ever been discovered by X-ray observatories, e.g., M82 X-1 (Matsumoto et al., 2001) and HLX-1 (Farrell et al., 2009). The abundance of IMBHs will place constraints on the formation scenarios of SMBHs (e.g. Rees, 1984).

A WD TDE can probably be observed as various types of transients. We focus on the possibility of the WD TDE observed as an optical transient resulting from its thermonuclear explosion, although there are many studies for possible signatures of WD TDEs as other types of transients (Zalamea et al., 2010; Clausen & Eracleous, 2011; Haas et al., 2012; MacLeod et al., 2014; Cheng & Bogdanović, 2014; East, 2014; Shiokawa et al., 2015; Ioka et al., 2016). In a WD TDE, the WD is heated by compression in the direction perpendicular to the orbital plane. Hereafter, the direction perpendicular to the orbital plane is called the $z$-direction. Then, the WD could undergo explosive nuclear burning, yielding radioactive nuclei such as 56Ni. (Luminet & Pichon, 1989; Rosswog et al., 2008, 2009; Haas et al., 2012; Sell et al., 2015). Rosswog et al. (2008, 2009) studied WD TDEs and their nucleosynthesis signatures. Using Rosswog et al.’s data of $0.6M_{\odot}$ carbon-oxygen (CO) WD, MacLeod et al. (2016) predicted the light curve and spectrum of the WD TDE.

To initiate explosive nuclear burning in a WD TDE, not only an adiabatic compression but also a shock heating are required in order to reach high temperature. As seen in Figure 1, a material consisting of helium (He), CO, or oxygen-neon-magnesium (ONeMg) needs to be adiabatically compressed by more than four orders of magnitude to rise its temperature from $10^{6}$ K to $3\times 10^{8}$ K (He), or $3\times 10^{9}$ K (CO and ONeMg), above which the explosive nuclear burning is triggered. However, it is difficult to compress a WD by more than four orders of magnitude. As seen from Figure 1 of Rosswog et al. (2009), $\beta$, which is the ratio of the tidal radius of the WD to a pericenter distance between the WD and IMBH, is permitted up to about $20$. The scale height of the WD in the $z$-direction, $z_{\textrm{min}}$, is estimated as $z_{\textrm{min}}/R_{\textrm{wd}}\sim\beta^{-3}$ at the pericenter (Luminet & Carter, 1986; Brassart & Luminet, 2008; Stone et al., 2013), where $R_{\textrm{wd}}$ is the original radius of the WD. Therefore, the WD is compressed by at most a factor of $8000$. Furthermore, we overestimate the compression of the WD in the above discussion, since the WD is not only compressed in the $z$-direction, but also elongated in the direction of the orbital plane.

An additional heating by the shock wave is therefore required for initiating the nuclear reactions, however it is not clear whether previous simulations have high resolution enough to detect such a shock wave. Rosswog et al. (2009) have performed smoothed particle hydrodynamics (SPH) simulations with up to $4$ million particles. Since more than $100$ particles are aligned in the $z$-direction, these simulations seem to resolve the shock wave. However, the WD is elongated in the direction of the orbital plane, and the number of particles in the $z$-direction should be fewer than $100$, where the shock front should be perpendicular to the $z$-direction (Kobayashi et al., 2004). Haas et al. (2012) have performed mesh-based simulations for a WD TDE of $0.6M_{\odot}$ CO WD encountered by a $1000M_{\odot}$ IMBH with $\beta=6$, and have suggested that explosive nuclear burning is triggered at the pericenter passage. If the original radius of the WD is $10^{9}$ cm, the scale height of the WD in the $z$-direction, $z_{\textrm{min}}$ should be $5\times 10^{6}$ cm. On the other hand, the finest mesh size of their simulation is about $10^{7}$ cm, larger than $z_{\textrm{min}}$.

In this paper, we perform high-resolution simulations of WD TDEs. We aim to investigate whether these simulations accurately follow a thermonuclear explosion in WD TDEs, and what is numerically required to follow the explosion.

This paper is structured as follows. In Section 2, we describe our simulation method. In Section 3, we present the simulation results. In Section 4, we discuss in detail the reason why the nuclear burning falsely occurs in 3D SPH simulations with low resolution. In Sections 5, we summarize this paper.

2 Method

In this section, we describe our simulation method. We overview our SPH code in section 2.1. We then describe setup of 3D SPH simulations in section 2.2, and of 1D simulations in section 2.3.

2.1 SPH code

Our SPH code solves the vanilla ice SPH equations. We adopt Wendland $C^{2}$ kernel for the SPH kernel interpolation (Wendland, 1995; Dehnen & Aly, 2012). Our SPH code is applicable to 3D and 1D planar geometry. The number of neighbor particles of a given particle is about $120$ (3D) and about $5$ (1D) unless otherwise specified. Neighbor particles are defined as particles which are inside a sphere centered at a given particle with its kernel-support radius, where the SPH kernel function reaches zero at the kernel-support radius. We choose artificial viscosity proposed by Monaghan (1997). The artificial viscosity is dependent on the strength of a shock wave (Morris & Monaghan, 1997). The viscosity from shear motion is suppressed by Balsara switch (Balsara, 1995). We calculate self gravity among SPH particles with adaptive gravitational softening (Price & Monaghan, 2007). We optimize the SPH and self-gravity calculations on distributed-memory systems, using FDPS (Iwasawa et al., 2015, 2016) and explicit AVX instructions (e.g. Tanikawa et al., 2012, 2013).

We use the Helmholtz equation of state (EoS) with (or without) Coulomb corrections (Timmes & Swesty, 2000), which considers partially relativistic and partially degenerate electrons including electron-positron pairs, ions, radiation, and Coulomb interactions. We include nuclear reactions with Aprox13 (Timmes et al., 2000) which solves $(\alpha,p)(p,\gamma)$ and $(\gamma,p)(p,\alpha)$ links as well as the $\alpha$-chain reaction network. The nuclear reactions are solved implicitly if a particle has high density ($>5\times 10^{7}$ g cm*-3*) and high temperature ($>3\times 10^{9}$ K), and otherwise solved explicitly. We can avoid overheating and overcooling with this implicit method, even if we use a large timestep, say $10^{-6}$ s222 Raskin et al. (2010) said that the timestep should be less than $10^{-12}$ s in the explicit method, in order to avoid the overcooling of photo disintegration.. We adopt the routines to calculate the Helmholtz EoS and Aprox13 developed by the Center for Astrophysical Thermonuclear Flashes at the University of Chicago.

2.2 Setup of 3D SPH simulations

We follow the evolution of WD TDEs by means of our SPH code coupled with the nuclear reaction network. We adopt SPH modeling for the WDs, and use a fixed potential to model the IMBH gravity.

For the IMBH potential, we adopt three kinds of models: Newtonian potential (NP), the Paczyński-Wiita (PW) potential (Paczyńsky & Wiita, 1980) with the modification by Rosswog (2005), and a generalized Newtonian potential obtained by Tejeda & Rosswog (2013), hereafter called TR potential, for the treatment of a BH with no spin (see also Tejeda et al. (2017) for the treatment of a spinning BH). The IMBH is located on the coordinate origin. None of them considers the IMBH spin.

Our initial conditions are summarized in Table 1. We relax the configurations of particles for these WDs in the same way as Tanikawa et al. (2015) (see also Sato et al., 2015, 2016). These WDs have no spin. We change the number of particles for the WDs, $N$, from a few $10^{4}$ up to a few $10^{7}$. The parameter set of model CO1 is the same as run 8 of Rosswog et al. (2009). However, the WD radii might be different between theirs and ours due to different ways to make initial conditions. This might make a difference between the pericenter distance of WDs in Rosswog’s run 8 and in our model CO1, although $\beta=5$ in both models.

Additionally, we investigate three models. The first and second ones are the same as models ONeMg and CO1, respectively, except without solving a nuclear reaction network. We call these models “ONeMg w/o nuc” and “CO1 w/o nuc”. The third one is the same as model ONeMg, except that the number of neighbor particles is set to be proportional to $N^{1/3}$, which is called “ONeMg-H”. This proportionality means that the kernel-support radii of particles are equal among different $N$ models if the density of the particles is equal. Note the kernel-support radii are proportional to $N^{-1/3}$ if the number of neighbor particles is fixed.

The WD in each simulation approaches the IMBH on a parabolic orbit. The initial separation between the WD and IMBH is set at 8 times the tidal radius of the WD, where the tidal radii are $7.1\times 10^{9}$ cm for models CO, $1.2\times 10^{10}$ cm for model He, and $1.5\times 10^{9}$ cm for model ONeMg. Our simulations follow the evolution of these WD TDEs for $60$ s for models CO, $80$ s for model He, and $10$ s for model ONeMg. At the end of the simulations, nuclear reactions have already ceased.

2.3 Setup of 1D simulations

We construct 1D initial conditions, by extracting 1D flow structure in the $z$-direction from the results of model ONeMg w/o nuc in the 3D SPH simulations. Details are described in section 3.2. Using these initial conditions, we perform 1D SPH simulations, and FLASH simulations.

Since the 1D SPH numerical methods are the same as the 3D SPH simulations except for the dimension, we overview the FLASH simulations here. The FLASH code (Fryxell et al., 2000) is an Eulerian code. We adopt uniform mesh, although the FLASH code supports adaptive mesh refinement. We choose the piecewise parabolic method (Colella & Woodward, 1984) for the gas hydrodynamic solver. We use the Helmholtz EoS and nuclear reactions with Aprox13, which are also used for our SPH simulations. The timestep is chosen to be the minimum value of the hydrodynamics timestep and nuclear reaction timestep. The hydrodynamics timestep is 10 % of the Courant-Friedrichs-Lewy number, and the nuclear reaction timestep is 1 % of the ratio of the specific internal energy to the specific nuclear energy-generation rate.

We describe the 1D initial conditions in section 3.2. Note that we make the 1D initial conditions, based on the results of the 3D SPH simulations (which is described in section 3.1).

3 Results

In section 3.1, we present the results of our 3D SPH simulations. The results show that the amount of the materials experiencing explosive nuclear burning is decreased with increasing $N$, and that the explosive nuclear burning found in these simulations is a numerical artifact due to low resolution, not physical effects. Furthermore, we do not find shock waves in our 3D SPH simulations. The absence of the shock waves may be also due to low mass resolution, even if $N>10^{7}$. In order to fix this problem, we perform 1D simulations with high resolution, by extracting data from the 3D SPH simulations. We show the results in section 3.2.

3.1 3D simulations

In Figure 2, we summarize the masses of the original, Si group, and Fe group elements at the time just after nuclear reactions have ceased. Mass accreted by an IMBH at this time is negligible. The difference in the nuclear burning products does not come from mass accretion, but comes from the nuclear reactions.

In model ONeMg, the amount of original (unburned) elements increases, and the amounts of Si and Fe group elements decrease, with increasing $N$. In models CO, the amount of Si group elements decreases monotonically with increasing $N$. Although the amount of the original elements decreases first, it finally increases from some $N$. The dependence of Fe group on $N$ is opposite to that of the original elements. In model He, the dependence of each element on $N$ is similar to those of models CO, although the dependence is smaller. Overall, the amount of unburned materials increases at large $N$ (say $N>10^{6}$) regardless of the WD masses and compositions.

Since the nucleosynthesis of model ONeMg indicates the strongest dependence on $N$, we show thermodynamical quantities of model ONeMg in Figure 3. At the time, nuclear reactions are still in progress. The density of particles is not sensitive to $N$ (see the top panels). On the other hand, their temperature strongly depends on $N$ (see the second top panels). High-temperature (say $>2\times 10^{9}$ K) region becomes smaller as $N$ increases. As the high-temperature region becomes smaller, the amount of the unburned materials increases (see the second bottom panels), and the amount of 56Ni decreases (see the bottom panels), since the materials consisting of ONeMg experience explosive nuclear burning at a temperature above $2.5\times 10^{9}$ K.

Figure 4 shows the structure of model ONeMg on a plane perpendicular to the orbital plane. Although most of the materials are burned out for $N=0.098$M, the kernel-support radius is comparable to the scale height of the structure in the $z$-direction. For $N=6.3$M, only a small portion of the WD experiences nuclear reactions. However, the scale height of the small portion is comparable to the kernel-support radius. A large portion of the WD, where the scale height is much larger than the kernel-support radius, is not burned.

Overall, the kernel-support radius becomes sufficiently smaller than the scale height as $N$ increases. Instead, nuclear reactions become inactive. The nuclear reactions occur only when the SPH simulations do not resolve the scale height of the structure in the $z$-direction. Therefore, we conclude that nuclear reactions artificially occur due to numerical effects, not due to physical effects, even if $N$ is a few $10^{7}$.

We investigate model ONeMg w/o nuc in order to confirm that nuclear reactions falsely occur, not due to any errors in our nuclear reaction network, but due to the resolution effect of our SPH simulations. Figure 5 shows thermodynamical quantities of model ONeMg w/o nuc on the orbital plane. The density distributions are similar to those in model ONeMg. The temperature distributions are different from those in model ONeMg. Overall, the temperature in this model is lower than in model ONeMg, since the nuclear reaction network is turned off.

The distribution of the maximum temperature each particle experiences from $t=0$ s to $t=6$ s is a good indicator for where explosive nuclear burning would occur from $t=0$ s to $t=6$ s if the nuclear reaction network was turned on. Materials which reach a temperature of $2.5\times 10^{9}$ K experience explosive nuclear burning if their density is more than several $10^{5}$ g cm*-3*. Therefore, explosive nuclear burning would occur in red regions. These red regions are almost coincident with the regions where original components are burned out and a large amount of 56Ni is produced in model ONeMg (see Figure 3). Moreover, these red regions shrink as $N$ increases, which is consistent with the results in model ONeMg.

Figure 6 shows the temperature distributions on cross sections in model ONeMg w/o nuc. The locations of the cross sections are the same as those in model ONeMg. In $N=0.098$M, temperature is more than $2.5\times 10^{9}$ K almost in the region. Therefore, explosive nuclear burning would occur. The kernel-support radii are comparable to the scale heights over all the range of $x$-coordinate. As $N$ increases, the kernel-support radii becomes smaller, and the temperature becomes lower. In $N=6.3$M, the kernel-support radii at $x=0.5\times 10^{9}$ cm are smaller than the scale height by a factor of a few, and the temperature is lower than $10^{9}$ K. Explosive nuclear burning would not occur there. At $x=0.3\times 10^{9}$ cm, the kernel-support radii are comparable to the scale height, and the temperature is more than $2.5\times 10^{9}$ K. Explosive nuclear burning would occur.

These results are consistent with those in model ONeMg. Therefore, explosive nuclear burning in model ONeMg occurs when SPH simulations fail to resolve the scale heights of WD TDEs. In other words, nuclear reactions seen in the low-resolution runs are a numerical artifact due to spurious heating of SPH simulation.

We comment the reasons why a leading part of the WD remains high temperature with increasing $N$, and why a trailing part of the WD gets lower temperature with increasing $N$ (see the bottom panels of Figure 5). At the pericenter passage, the leading part is closer to an IMBH than the trailing part. The leading part has smaller scale height than the trailing part. Our SPH simulations almost resolve the scale height of the trailing part unless $N<10^{5}$ (see Figure 6). However, our SPH simulations fail to resolve the scale height of the leading part, even if $N>10^{7}$. Consequently, the leading part gets high temperature, even if $N>10^{7}$.

There is another evidence that the nuclear reactions occur only when the kernel-support radii of particles are comparable to the scale height of the structure. In model ONeMg-H, the kernel-support radii do not become smaller by design even if $N$ increases. Figure 2 shows that the nucleosynthesis of model ONeMg-H is much less sensitive to $N$ than that of model ONeMg.

Finally, we search for shock waves in model ONeMg w/o nuc, focusing on regions where the SPH kernel is smaller than the scale height. However we find no shock wave. If the shock waves were generated, temperature would increase discontinuously in the $z$-direction. However, we do not find such a temperature structure, for example in the right panel of Figure 6.

3.2 1D simulations

We do not find shock waves in model ONeMg w/o nuc, despite the fact that SPH kernel-support radii are smaller than the scale height (see the right panel of Figure 6). We suspect that the absence of the shock waves is due to the low resolution of our simulations, even with $N>10^{7}$. In order to clarify the issue, we perform 1D planar simulations with high resolution, by extracting local flow structures from model ONeMg w/o nuc with $N=6.3$M. For the 1D planar simulations, we use our 1D SPH code and FLASH code.

We make initial conditions for the 1D simulations as follows. We re-perform 3D SPH simulations in order to record physical quantities of particles at every timestep. The recorded particles are located at regions pointed by white crosses in the rightmost panel of Figure 5. From bottom to top, particles get higher temperature. So, we call these regions “low-T region”, “middle-T region”, and “high-T region” from bottom to top.

From the recorded data in each region, we extract physical quantities at the moment when $z$-direction relative velocities between the outermost particle and a particle on the orbital plane reaches a peak. Based on distribution of density and velocity of the extracted data along the $z$-direction, we map particles for 1D simulations, and assign physical quantities for 1D mesh. We fix the temperature to be $10^{7}$ K. For the FLASH simulations, we put atmosphere with a constant density of $1$ g cm*-3* if there is no WD material. The total mass of atmosphere is negligible, ($\lesssim 0.0001$ % of the total mass of WD materials).

We describe the recipe of the 1D SPH and FLASH simulations in the following. We consider only hydrodynamics for both 1D SPH and FLASH simulations. We do not solve nuclear reactions unless otherwise specified. Furthermore, neither do we consider self gravity among particles, nor the IMBH gravity. The number of particles we use is $N=10$, $10^{3}$, and $10^{4}$ for the 1D SPH simulations. The number of grids we use is $N_{\textrm{g}}=40$, $1600$ and $3200$. Note that the space resolution in the largest $N$ and $N_{\textrm{g}}$ is $\sim 10^{4}$ cm and $\sim 10^{5}$ cm for the 1D SPH and FLASH simulations, respectively. In some cases of $N_{\textrm{g}}=3200$, we perform additional simulations, solving hydrodynamics coupled with the nuclear reaction network.

In the cases of the low-T and high-T regions, we perform additional 1D SPH simulations with $N=10^{4}$ and FLASH simulations with $N_{\textrm{g}}=3200$. In these simulations, we change the initial $z$-direction velocity by increasing the $z$-direction velocity for the low-T region, and by decreasing the $z$-direction velocity for the high-T region, so that density at $z=0$ in 1D simulations is equal to that in 3D simulations when materials are most compressed. The reason for this additional setups is the following. We underestimate and overestimate the compression of materials in the low-T and high-T regions, respectively, without the correction of the $z$-direction velocity (hereafter, $v_{\textrm{z}}$ correction). These underestimate and overestimate come from the facts that we do not consider an IMBH and self gravity, and that materials are not elongated in the direction of the orbital plane in these 1D simulations, respectively.

Figure 7 shows the results of the 1D simulations. We set the initial time (i.e. $t=0$ s) to be the time when these 1D simulations are started. For the low-T region with $N=10$, the temperature on the orbital plane becomes high due to the spurious heating. Actually, the spurious heating is seen in the results of low-resolution simulations, i.e. the 1D SPH simulations with $N=10$ and the FLASH simulations with $N_{\textrm{g}}=40$ (see Figure 8), regardless of the regions. This is consistent with our 3D SPH simulations in section 3.1.

As seen in Figure 7, the results of the 1D SPH simulations are in good agreement with those of FLASH simulations, if the simulations have high resolution, i.e. the 1D SPH simulations with $N=10^{4}$ and the FLASH simulations with $N_{\textrm{g}}=3200$. For this comparison, atmosphere should be ignored. Furthermore, we investigate the convergence of these results against $N$ and $N_{\textrm{g}}$. Comparing the results of $N=10^{3}$ and $10^{4}$ in the 1D SPH simulations, and those of $N_{\textrm{g}}=1600$ and $3200$ in the FLASH simulations, we find these results are quite similar. Therefore, these results are converged against $N$ and $N_{\textrm{g}}$.

As seen in the high-resolution results of Figure 7, in all the regions, shock waves are observed. In order to investigate whether these shock waves trigger detonation waves, we perform FLASH simulations coupled with the nuclear reaction network for the middle-T region. We consider two cases where materials consist of CO and ONeMg.

Fluid motion for the CO and ONeMg cases is the same as fluid motion for the case without considering nuclear reactions until a shock wave emerges at $t\sim 0.0391$ s. The middle-T panels of Figure 7 show the fluid motion in the case without considering nuclear reactions. At $t=0$ s, all the materials shrink. At $t\sim 0.0273$ s, materials at $z\sim 0$ cm bounce back. At $t\sim 0.0391$ s, a shock wave is generated at $z\sim 1.7\times 10^{7}$ cm. At this time, materials at $z\lesssim 1.7\times 10^{7}$ cm are expanding, and the rest are shrinking. The fraction of the expanding materials is $99.96$ %. We emphasize that the shrinking materials (i.e. $0.04$ % of materials) are not atmosphere.

Figure 9 shows the time evolution of $z$-direction velocity, density, temperature, and mass fractions of nuclear elements for the two cases after the shock wave emerges at $t\sim 0.0391$ s. For both the cases, the shock waves which emerges at $t\sim 0.0391$ s proceed in the positive $z$-direction (see panels (a) and (b) of Figure 9). Hereafter, a shock wave proceeding in the positive $z$-direction is referred as “forward shock wave”. Only for the CO case, a reverse shock wave emerges at $z\sim 1.7\times 10^{7}$ cm, which can be seen as the discontinuous decrease of $z$-direction velocity (see Figure 10 and panel (a) of Figure 9). The reverse shock wave is formed, such that materials behind the reverse shock wave experience explosive nuclear burning (see panels (g) and (i) of Figure 9), get high temperature (see panel (e) of Figure 9), rapidly expand, and push back materials in front of the reverse shock wave. The forward shock wave for the CO case has higher velocity than for the ONeMg case, since the former is energized by more active nuclear reactions than the latter.

In the following three reasons, we conclude that the reverse shock wave accompanies a detonation wave standing at $z\sim 1.7\times 10^{7}$ cm. First, upstream materials flow across $z\sim 1.7\times 10^{7}$ cm at a speed of $\sim 10^{9}$ cm s*-1* (see panel (a) of Figure 9), which is larger than the sound velocity of the upstream materials ($\lesssim 5\times 10^{8}$ cm s*-1*). Second, CO elements are burned out, and Si and Fe group elements are synthesized just behind $z\sim 1.7\times 10^{7}$ cm (see panels (g) and (i) of Figure 9). Third, this fluid structure is long-lived. It does not decay during the simulation, although we only perform the simulation until $t\sim 0.0422$ s due to high calculation cost of the nuclear reaction network. From $t\sim 0.0391$ s to $t\sim 0.0422$ s, $7$ % of materials flow across $z\sim 1.7\times 10^{7}$ cm. After they experience explosive nuclear burning, they consist of $75$ % Fe group elements, $20$ % Si group elements, and $4$ % He.

If we continue the simulation, the detonation wave will be sustained until at least $t\sim 0.0625$ s in the following reason. Materials flow supersonically across $z\sim 1.7\times 10^{7}$ cm from $t\sim 0.0391$ s to $t\sim 0.0625$ s in the case without considering nuclear reactions. This should be also true in the CO case if the detonation wave stands at $z\sim 1.7\times 10^{7}$ cm during this time. This is because upstream materials of the detonation wave should not receive the influence of nuclear reactions due to their supersonic motions. Then, for the CO case, materials should flow supersonically into the detonation wave, experience explosive nuclear burning, and sustain the detonation wave as fuels. During this time, more than $50$ % of materials should flow across the detonation wave. We expect a large amount of radioactive nuclei is synthesized.

There is only the single detonation wave at $z\sim 1.7\times 10^{7}$ cm; the forward shock wave at $t=0.0410$ s and $z\sim 2.3\times 10^{7}$ cm does not accompany a detonation wave (see the left panels of Figure 9). We explain the reason for the formation of only the single detonation wave at $z\sim 1.7\times 10^{7}$ cm. The forward shock wave creates a hotspot at $t\sim 0.0391$ s and $t\sim 1.7\times 10^{7}$ cm. The hotspot may possibly initiates double detonation waves initially; one is the detonation wave standing at $\sim 1.7\times 10^{7}$ cm after $t\sim 0.0391$ s, and the other may be a detonation wave which follows the forward shock wave creating the hotspot. The former detonation wave is sustained by high-density fuels ($\sim 3\times 10^{7}$ g cm*-3*), while the latter detonation wave may soon decay due to the supply of only low-density fuels ($<10^{5}$ g cm*-3*), as seen in panel (c) of Figure 9.

For the convergence check of a space resolution, we also perform FLASH simulations with a space resolution of $\sim 10^{6}$ cm to $\sim 10^{7}$ cm for the middle-T region in the CO case. We find the same detonation wave as in a space resolution of $\sim 10^{5}$ cm only when the space resolution is $\lesssim 10^{6}$ cm. Therefore, we need simulation with a space resolution of $\lesssim 10^{6}$ cm in order to follow such a detonation wave.

For the ONeMg case, the forward shock wave raises the temperature of materials to $>2\times 10^{9}$ K (see panel (f) of Figure 9), and triggers nuclear reactions slightly (see panels (h) and (j) of Figure 9). However, the nuclear reactions cease soon. The forward shock wave does not excite explosive nuclear burning triggering a reverse shock wave and a detonation wave (see panel (b) of Figure 9).

A forward shock wave does not always excite explosive nuclear burning triggering a reverse shock wave and a detonation wave even when materials consist of CO. We also perform a FLASH simulation coupled with the nuclear reaction network for the low-T region with $v_{\textrm{z}}$ correction in which materials consist of CO. Figure 11 shows the time evolution. At $t\sim 0.0996$ s, a forward shock wave emerges at $z\sim 7.3\times 10^{7}$ cm (see the top panel of Figure 11). The forward shock wave raises the temperature of the materials (see the bottom panel of Figure 11). It triggers nuclear reactions, however does not initiate a reverse shock wave nor a detonation wave as seen in the top panel of Figure 11. This is quite similar to the ONeMg case in the middle-T region (see panels (b) and (f) of Figure 9).

The initiation of a detonation wave strongly depends on where a forward shock wave emerges. Figure 12 shows density and temperature in the low-T region with $v_{\textrm{z}}$ correction and in the middle-T region when forward shock waves emerge. The forward shock wave emerges at density $\sim 10^{6}$ and $\sim 10^{7}$ g cm*-3* for the former and latter cases, respectively. According to table 6 in Seitenzahl et al. (2009), spontaneous initiation of a detonation wave requires a hotspot with size of $\sim 10^{7}$ and $\sim 10^{5}$ cm for density $\sim 10^{6}$ and $\sim 10^{7}$ g cm*-3*, respectively. Since a mesh size is $\sim 10^{5}$ cm, the hotspot is too small to generate a detonation wave for the former case, and is sufficiently large for the latter case. Our results about the initiation of the detonation wave are consistent with the criteria in Seitenzahl et al. (2009).

As seen in table 8, 9, 10, and 11 of Seitenzahl et al. (2009), the hotspot size to generate a detonation wave becomes rapidly larger as carbon fraction becomes smaller. Even when the carbon fraction decreases from $50$ % to $30$ %, a hotspot size required to generate a detonation wave becomes larger by a factor of $\sim 10$. Therefore, the reason why no detonation wave occurs for the ONeMg case in the middle-T region (see the right panels of Figure 9) is that a hotspot is too small to generate a detonation wave in ONeMg compositions.

From the above, we find two implications. First, the resolution of 3D SPH simulations is too low to resolve a forward shock wave even if $N>10^{7}$. Second, the forward shock wave does not always initiate a detonation wave even if materials consist of CO.

Finally, we discuss the reliability of the initiation of a detonation wave for the CO case in the middle-T region. In this case, the forward shock wave occurs, and create the hotspot generating the detonation wave where the density gradient is steep. As seen in Figure 12, the hotspot has density of $10^{7}$ g cm*-3*, however a mesh near the hotspot has density of $\sim 10^{6}$ g cm*-3*. If the forward shock wave emerges at a bit outer region, it can not trigger a detonation wave, since the forward shock wave propagates outward (see panel (a) of Figure 9). Our 1D initial conditions lack accuracy to precisely determine where a forward shock wave occurs. Therefore, the initiation of a detonation wave may be unreliable. In order to follow the initiation of a detonation wave accurately, we should perform 3D simulations with a space resolution of $\lesssim 10^{6}$ cm. Note that the detonation wave in the middle-T region occurs only when a space resolution is $\lesssim 10^{6}$ cm.

4 Discussion

In this section, we discuss why materials are falsely heated in 3D SPH simulation with low resolution, using model CO1 which is similar to model of Rosswog’s run 8.

Model CO1 provides a solid basis for the comparison, since this parameter set is the same as that of Rosswog’s run 8, and MacLeod et al. (2016) have investigated the observational features of this model. In our model, the WDs yield $0.13M_{\odot}$ and $0.15M_{\odot}$ of Fe group elements for $N=0.39$M and $N=0.79$M, respectively. On the other hand, Rosswog’s run 8 has $0.18M_{\odot}$ for $N=0.5$M. Our results are consistent with those of Rosswog’s run 8 from the view point of the nucleosynthesis.

Figure 13 shows thermodynamical quantities of model CO1 in the same way as in Figure 3. As seen in the bottom panels, the amount of 56Ni increases from $N=0.098$M to $N=0.39$M, and decreases from $N=1.6$M to $N=6.3$M, which is consistent with the amount of 56Ni drawn in Figure 2. The increase from $N=0.098$M to $N=0.39$M comes from the increase of the density of particles with increasing $N$ (see the top panels). Note that 56Ni is more synthesized (while artificially) at higher density if materials experience explosive nuclear burning. The low density in $N=0.098$M results in the lower temperature and larger unburned mass in $N=0.098$M than those in $N=0.39$M. Since the density is low, the nuclear reactions are not active so much.

The dependence of the density on $N$ can be explained as follows. As $N$ becomes smaller, a particle on the orbital plane has a larger kernel-support radius. SPH simulation calculates the density of the particle, considering more distant particles from the orbital plane. Particles are distributed more sparsely with increasing distance from the orbital plane. Therefore, the density of the particle becomes smaller with a larger kernel-support radius (smaller $N$).

The decrease of 56Ni from $N=1.6$M to $N=6.3$M is also due to the same reason as the decrease in model ONeMg. Figure 14 shows the structure of model CO1 on a plane perpendicular to the orbital plane. In $N=6.3$M, the kernel-support radii of the particles are smaller than the scale height, except around $x=10^{9}$ cm where explosive nuclear burning occurs. In $N=0.098$M, $0.39$M, and $1.6$M, the kernel-support radii seem smaller than the scale height, even where explosive nuclear burning occurs. However, the kernel-support radius is comparable to the scale height at the time earlier than $t=26$ s.

In the following, we clarify the reason why the nuclear reactions occur if the kernel-support radius is comparable to the scale height in the $z$-direction. For this purpose, we follow the time evolution of particles located at the crosses in Figure 13. In order to focus on whether the nuclear reactions start or not, we adopt model CO1 w/o nuc.

Figure 15 shows the density and temperature of these particles in $N=0.39$M and $6.3$M. Since the 12C + 12C timescale is shorter than the local dynamical timescale in a part of particles in both $N$, they would experience the explosive nuclear burning if the nuclear reaction network was turned on.

We can see the trajectories of these particles in Figure 16. The trajectories are independent of $N$ from $x=-4.5\times 10^{9}$ cm to $x=0$ cm. The particles are in the range of $z=\pm 0.4\times 10^{9}$ cm at $x=-4.5\times 10^{9}$ cm, and in the range of $z=\pm 0.1\times 10^{9}$ cm at $x=0$ cm. At $x=0.8\times 10^{9}$ cm, the particles in $N=6.3$M are distributed slightly more densely than those in $N=0.39$M. Since these particles approach the orbital plane, the scale height of the structure shrinks in the $z$-direction. On the other hand, the kernel-support radii of these particles keep nearly constant; the density of the structure does not grow. This is because the structure is compressed in the $z$-direction, but is extended in the direction of the orbital plane.

The entropy of these particles grows from some point. Although their entropy increases, no shock wave is found out. If there is a shock wave, the trajectories of these particles should be changed discontinuously. Moreover, their entropy in $N=0.39$M grows at smaller $x$ (i.e. earlier) than those in $N=6.3$M. These behaviours of their entropy are clearly not converged against $N$.

We investigate why the entropy increases despite the absence of a shock wave. Figure 17 shows the velocities in the $z$-direction and sound velocities of these particles. Particles distant from the orbital plane approach the orbital plane supersonically from $x=-3\times 10^{9}$ cm to $x=0.8\times 10^{9}$ cm, while particles near the orbital plane have small velocities in the $z$-direction. Therefore, the relative velocities between the outermost and innermost particles from the orbital plane are supersonic. Since all these particles have the kernel-support radii comparable to the scale height, the innermost particles interact supersonically with the outermost particles. Then, the particles obtain entropy.

The reason why particles in smaller $N$ obtain entropy earlier is as follows. The kernel-support radii of particles are larger in smaller $N$. Therefore, in smaller $N$, the innermost particles start interacting with the outermost particles when the scale height is still larger.

As $N$ becomes larger, the kernel-support radii become smaller. Each of the innermost and outermost particles interacts only with their adjacent particles. Their relative velocities are subsonic. Consequently, these particles do not gain entropy in large $N$.

In summary, particles on the orbital plane are heated and gain entropy, interacting with the outermost particles from the orbital plane due to their kernel-support radii comparable to the scale height. Such interactions should not happen in reality, and the heating is spurious. The spurious heating triggers explosive nuclear burning falsely. Eventually, our 3D SPH simulations fail to follow the nucleosynthesis in WD TDEs, even if our simulations have unprecedentedly high resolution. This should also be the case for previous simulations with lower resolution.

5 Summary

We perform 3D SPH simulations with $N>10^{7}$ to follow the TDEs of He WDs, CO WDs, and ONeMg WDs by IMBHs. We observe that the explosive nuclear burning occurs in all the WDs. However, the final compositions are strongly dependent on $N$. We find that the amount of unburned materials increases with increasing $N$. The nucleosynthesis of these WDs is not converged against $N$, although our simulations contain unprecedentedly large $N$, up to a few $10^{7}$ particles.

The reason why the explosive nuclear burning occurs in small-$N$ simulations is as follows. In such small-$N$ simulations, the scale height of a WD in the $z$-direction becomes comparable to the kernel-support radii of particles at the pericenter passage around an IMBH. On the other hand, a particle at the surface of the WD approaches the orbital plane supersonically. Then, the particle interacts with other particles on the orbital plane. Their interactions generate spurious heating, and trigger explosive nuclear burning falsely. If $N$ is sufficiently large, these particles do not interact with each other, since the kernel-support radius is sufficiently small.

By means of mesh-based simulations, Haas et al. (2012) have suggested that a WD TDE can experience explosive nuclear burning. However, we conjecture that their mesh-based simulations have not either resolved the scale height of the WD in the $z$-direction at the pericenter passage. They have investigated TDEs with $0.6M_{\odot}$ CO WDs which orbit around a $10^{3}M_{\odot}$ IMBH with $\beta=6$. This model is similar to our model CO. In our model, the scale height in the $z$-direction is about $2\times 10^{7}$ cm at the pericenter passage when $\beta=5$. If the scale height is almost equal to those in Haas et al. (2012), their simulations can not resolve the scale height, since their finest mesh size is about $10^{7}$ cm.

We find no shock wave generation in 3D SPH simulations with $N>10^{7}$. So, we perform 1D SPH and FLASH simulations with higher resolution than the 3D SPH simulations in order to examine whether the absence of shock waves results from the low resolution of the 3D SPH simulations. Eventually, we find the following two insights. First, a shock wave emerges in both 1D SPH and FLASH simulations. In other words, the resolution of the current 3D SPH simulations is not enough to resolve the shock wave, even if $N>10^{7}$. Second, the shock wave can trigger a detonation wave, only when it makes a hotspot large enough to generate a detonation wave. For the case of materials consisting of ONeMg, no shock wave triggers a detonation. For the case of materials consisting of CO, if the shock wave occurs in a low-density region ($\lesssim 10^{6}$ g cm*-3*), a hotspot formed by the shock wave is too small to generate a detonation wave.

Although a detonation wave occurs in a part of our 1D simulations, the initiation of the detonation wave may be unreliable for the following reason. The shock wave triggering the detonation wave occurs where the density gradient is steep. If the shock wave occurs at a bit outer region, it can not trigger the detonation wave, since it occurs in much lower density region by a factor of $10$. Our 1D modeling lacks accuracy to precisely determine where the shock wave occurs.

In order to conclude that WD TDEs become optical transients powered by radioactive nuclei, we need to perform 3D simulations with a space resolution of $\lesssim 10^{6}$ cm. SPH simulation with $N\sim 10^{9}$ seems to satisfy this requirement. However, this may not be true in the following reason. A shock wave emerges at the surface of a WD. Generally, SPH simulation has lower space resolution at the surface of an object than inside. Therefore, SPH simulation may need $N\gg 10^{9}$ particles.

A. Tanikawa thanks K. Kawana, N. Yoshida, Y. Guo, and K. Sugiura for fruitful discussions. Numerical computations were carried out on Cray XC30 at Center for Computational Astrophysics, National Astronomical Observatory of Japan, and on Cray XC40 at Yukawa Institute for Theoretical Physics, Kyoto University. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. This research has been supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, by MEXT program for the Development and Improvement for the Next Generation Ultra High-Speed Computer System under its Subsidies for Operating the Specific Advanced Large Research Facilities, and by Grants-in-Aid for Scientific Research (24540227, 26400222, 26800100, 16H02168, 16K17656) from the Japan Society for the Promotion of Science. Finally, we thank the anonymous referee for helpful advice.