Detecting Population III Stars through Tidal Disruption Events in the Era of JWST and Roman

We consider Pop III stars in the main-sequence stage. Massive Pop III stars are primarily radiation-pressure dominated, so their structures can be approximated using a polytropic model with an index of γ = 4/3 or $n=\tfrac{1}{\gamma -1}=3$ (Bromm et al. 2001b). Throughout this work, we use the mass–radius relationship of Pop III stars adopted from Bromm et al. (2001b):

$$\begin{eqnarray}&&{R}_{\star }\simeq 0.7{R}_{\odot }{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{0.45}{\left(\displaystyle \frac{Z}{{10}^{-9}}\right)}^{0.09},\end{eqnarray} \tag{ 1 }$$

where M_⋆ is the stellar mass, R_⋆ is the stellar radius, and Z is the metallicity of Pop III stars. When Z increases, cooling through metal lines is enhanced and the effective temperature of a Pop III star drops, which means the average stellar density also decreases given ρ_⋆ ∝ T³ for a polytropic star with n = 3 (Fowler & Hoyle 1964). Therefore, R_⋆ increases with increasing Z. The value of the metallicity Z considered has varied from zero to a very small nonzero number in previous treatments (Bromm et al. 2001b; Schaerer 2003; Murphy et al. 2021; Klessen & Glover 2023). Many studies show that there likely exists a critical metallicity at which the transition from metal-free Pop III stars to metal-poor Population II (Pop II) stars occurs, and that value is ∼10⁻³–10⁻⁵ Z_⊙ depending on the IMF of Pop III stars (Bromm et al. 2001a; Schneider et al. 2002; Yoshida et al. 2004; Wise et al. 2012; Jaacks et al. 2018). Furthermore, Jaacks et al. (2018) show that the mean gas metallicity increases with decreasing redshift z and reaches ∼10⁻⁵ Z_⊙ at z ∼ 10. Hereafter, in this work in order to probe the metallicity dependence of our results, we calculate fiducial properties for Pop III stars with two chosen values for the metallicity: Z = 10⁻⁵ and a lower metallicity of Z = 10⁻⁹.

Given the stellar structure outlined above, we can calculate the average stellar density ρ_⋆ of Pop III stars based on Equation (1). In Figure 1, we show the variation of ρ_⋆ for Pop III stars as a function of their mass in the typical mass range of 30–300M_⊙ and compare that with the case of Pop I main-sequence stars. For the latter, we assume that their masses can extend to few hundreds of solar masses (Crowther et al. 2010; Rickard & Pauli 2023), and we adopt the mass–radius relation from Kippenhahn & Weigert (1990):

$$\begin{eqnarray}{R}_{\star }\sim \left\{\begin{array}{ll}{R}_{\odot }{\left(\tfrac{{M}_{\star }}{{M}_{\odot }}\right)}^{0.8}; & {M}_{\star }\leqslant 1{M}_{\odot }\\ {R}_{\odot }{\left(\tfrac{{M}_{\star }}{{M}_{\odot }}\right)}^{0.57}; & {M}_{\star }\gt 1{M}_{\odot }\end{array}\right..\end{eqnarray} \tag{ 2 }$$

It can be seen that for both types of stars, ρ_⋆ drops as the M_⋆ increases. Furthermore, Pop III stars have ρ_⋆ ∼ ρ_⊙ for the lower-metallicity case Z ∼ 10⁻⁹ and ρ_⋆ ∼ 0.1ρ_⊙ as Z approaches the upper limit of 10⁻⁵.

Zoom In Zoom Out Reset image size

A star is tidally disrupted by a MBH when it approaches within its tidal disruption radius (R_T ). To first order, R_T can be calculated using

$$\begin{eqnarray}&&{R}_{T}={R}_{\star }{\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{{M}_{\star }}\right)}^{1/3},\end{eqnarray} \tag{ 3 }$$

where M_BH is the black hole (BH) mass. Therefore, for a fixed M_BH, R_T only depends on the average stellar density (${R}_{T}\propto {\rho }_{\star }^{-1/3}$). Figure 2 shows a comparison of R_T for Pop III stars (of different masses and metallicities) and with that for solar-type stars disrupted by MBHs with masses ranging from M_BH = 10⁵–10⁹ M_⊙. It can be seen that the R_T for Pop III stars at very low metallicity (Z ∼ 10⁻⁹) are similar to that of the Sun, while the R_T for Pop III stars at relatively high metallicities (Z ∼ 10⁻⁵) are a few times larger, which is consistent with the density comparison seen in Figure 1.

Zoom In Zoom Out Reset image size

Moreover, since ${R}_{T}\propto {M}_{\mathrm{BH}}^{1/3}$ while the gravitational radius of a BH R_g ≡ GM_BH/c² ∝ M_BH, there exists a upper limit of M_BH (called the Hill mass) beyond which an approaching star is swallowed by the MBH as a whole. One can see in Figure 2 that the high-metallicity (Z ∼ 10⁻⁵) Pop III stars can be disrupted by MBHs with masses up to 10⁹ M_⊙. This can be used as an important indicator for Pop III TDE detection, since Pop I main-sequence stars can only be disrupted by MBHs with M_BH ≲ 10⁸ M_⊙, except in the extreme case when the MBH has a close to maximal spin (Kesden 2012). Although the estimated mass of the BH in sources like UHZ1 (at z ≈ 10.1) is estimated to be ∼4 × 10⁷ M_⊙, it is expected that even more massive SMBHs could exist at these epochs, although they are expected to be extremely rare. Meanwhile, MBHs with masses less than that of UHZ1 could be significantly more numerous at these early epochs.

After the star is disrupted, about half of the stellar debris is ejected from the system and the rest remains bound to the MBH with a spread in its specific binding energy. The rate that the bound stellar debris orbits back to the pericenter is called the debris mass fallback rate (${\dot{M}}_{\mathrm{fb}}$), which can be calculated as (Evans & Kochanek 1989; Phinney 1989)

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{fb}}(t)\approx \displaystyle \frac{1}{3}\displaystyle \frac{{M}_{\star }}{{t}_{\mathrm{fb}}}{\left(\displaystyle \frac{t}{{t}_{\mathrm{fb}}}\right)}^{-5/3},\end{eqnarray} \tag{ 4 }$$

where the fallback time, t_fb, is the orbital time of the most tightly bound debris.

In this work, we adopt the results from Guillochon & Ramirez-Ruiz (2013, hereafter GR13) who performed high-resolution hydrodynamical simulations to study this disruption process. In particular, adopting polytropic stellar models, GR13 provide the following fitting formulae for the debris fallback time and the peak mass fallback rate:

$$\begin{eqnarray}&&{t}_{\mathrm{fb}}={B}_{\gamma }{M}_{6}^{1/2}{m}_{\star }^{-1}{r}_{\star }^{3/2}\,\mathrm{yr},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{fb},\mathrm{peak}}={A}_{\gamma }{M}_{6}^{-1/2}{m}_{\star }^{2}{r}_{\star }^{-3/2}\,{M}_{\odot }\ {\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 6 }$$

where M₆ = M_BH/10⁶ M_⊙, m_⋆ = M_⋆/M_⊙, and r_⋆ = R_⋆/R_⊙. Furthermore, for a polytropic star with γ = 4/3, the coefficients are as follows:

$$\begin{eqnarray}&&{B}_{4/3}=\displaystyle \frac{-0.38670+0.57291\sqrt{\beta }-0.31231\beta }{1-1.2744\sqrt{\beta }-0.90053\beta },\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{A}_{4/3}=\exp \left[\displaystyle \frac{27.261-27.516\beta +3.8716{\beta }^{2}}{1-3.2605\beta -1.3865{\beta }^{2}}\right].\end{eqnarray} \tag{ 8 }$$

Here β ≡ R_T /R_p is called the penetration factor, where R_p is the pericenter distance of the stellar orbit. A deep plunging orbit is commonly associated with β ≫ 1, whereas mild or partial disruption is denoted by β ∼ 1 and β ≲ 1, respectively. GR13 showed that stars with polytropic index γ = 4/3 are fully disrupted when β ≳ 1.85, and under this condition ${\dot{M}}_{\mathrm{fb},\mathrm{peak}}$ remains similar. Hence, we mark this as the critical penetration factor β_c and use β = β_c = 1.85 for the calculations throughout this Letter unless otherwise specified. For β = 1.85, B_4/3 ∼ 0.08 and A_4/3 ∼ 3.

For Pop III stars, one can then use the mass–radius relationship (Equation (1)) and rewrite Equations (5) and (6) as

$$\begin{eqnarray}&&{t}_{\mathrm{fb}}\simeq 9\left(\displaystyle \frac{{B}_{4/3}}{0.08}\right){\left(\displaystyle \frac{{m}_{\star }}{300}\right)}^{-0.3}{\left(\displaystyle \frac{Z}{{10}^{-5}}\right)}^{0.1}{M}_{6}^{0.5}\,\mathrm{day},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\begin{array}{l}{\dot{M}}_{\mathrm{fb},\mathrm{peak}}\simeq 2.7\times {10}^{3}\left(\displaystyle \frac{{A}_{4/3}}{3}\right){\left(\displaystyle \frac{{m}_{\star }}{300}\right)}^{1.3}{\left(\displaystyle \frac{Z}{{10}^{-5}}\right)}^{-0.1}\\ \qquad \qquad \qquad \times \,{M}_{6}^{-0.5}\,{M}_{\odot }\ {\mathrm{yr}}^{-1}.\end{array}\end{eqnarray} \tag{ 10 }$$

We plot t_fb and ${\dot{M}}_{\mathrm{fb},\mathrm{peak}}$ as functions of M_BH for various Pop III stars compared to the Sun in Figure 3. It can be seen that t_fb for Pop III stars are usually around a few days, which is shorter compared to that of a solar-type star wherein t_fb ∼a few tens of days for the case of M_BH = 10⁶ M_⊙. Furthermore, Pop III stars have extremely high ${\dot{M}}_{\mathrm{fb},\mathrm{peak}}$, which exceeds that of Pop I star TDEs by several orders of magnitude. For example, when M_BH = 10⁶ M_⊙, for Pop III star TDEs ${\dot{M}}_{\mathrm{fb},\mathrm{peak}}$ can reach ${10}^{4-6}{\dot{M}}_{\mathrm{Edd}}$, where ${\dot{M}}_{\mathrm{Edd}}\equiv {L}_{\mathrm{Edd}}/\eta {c}^{2}$ is the Eddington accretion rate of the black hole, with L_Edd being the Eddington luminosity, c the speed of light, and η the radiative efficiency with a nominal value of 0.1. This hyper-Eddington debris mass fallback rate is a consequence of the very large masses of Pop III stars coupled with their short debris mass fallback timescales.

Zoom In Zoom Out Reset image size

After the peak of the flare, the debris mass fallback rate drops with time following t^−5/3. At late times, the fallback rate should transit from super-Eddington to sub-Eddington. We can calculate the timescale over which the fallback rate stays super-Eddington using ${\dot{M}}_{\mathrm{fb}}={\dot{M}}_{\mathrm{Edd}}$:

$$\begin{eqnarray}&&{t}_{\mathrm{Edd}}\approx 38.5{\left(\displaystyle \frac{{m}_{\star }}{300}\right)}^{3/5}{\left(\displaystyle \frac{{t}_{\mathrm{fb}}}{10\,\mathrm{days}}\right)}^{2/5}{M}_{6}^{-3/5}\ \mathrm{yr}.\end{eqnarray} \tag{ 11 }$$

This timescale is shown in Figure 4. It is clearly seen that Pop III TDEs have super-Eddington fallback rates for a much longer duration (up to hundreds of years) compared to a solar-type star TDE (a few years) considering M_BH = 10⁶ M_⊙.

Zoom In Zoom Out Reset image size

We also note that there exist alternative models for the disruption process of stars. For example, a recent semi-analytical work by Bandopadhyay et al. (2024) showed that the fallback time could be almost independent of the stellar mass (in Pop I star TDEs). Yet these different models likely lead to consistent first-order results for the peak fallback rates, which matters the most for the observed flux (see Appendix A for comparison between models and more discussion).

The extremely high debris mass fallback rate seen in Figure 3 means that even if just a small fraction of the debris reaches the vicinity of the MBH, a super-Eddington accretion flow could result, which will launch powerful winds due to the large radiation pressure (Dai et al. 2018; Bu et al. 2023). Furthermore, outflows are also expected to be powered by debris stream collisions (Shiokawa et al. 2015; Bonnerot & Rossi 2019; Lu & Bonnerot 2020).

It has been proposed that these outflows are responsible for producing the optical/UV emission observed from TDEs in the local Universe (Loeb & Ulmer 1997; Strubbe & Quataert 2009; Lodato & Rossi 2011; Metzger & Stone 2016; Roth et al. 2016, 2020; Dai et al. 2021). There have also been sophisticated numerical simulations studying super-Eddington black hole accretion, outflow, and emission (Sadowski & Narayan 2016; Dai et al. 2018; Jiang et al. 2019; Thomsen et al. 2022). However, no simulation has been done yet to study super-Eddington accretion flows with Eddington ratios as high as 10^5–7. Given lack of guidance from simulations, we adopt an analytical model proposed by Strubbe & Quataert (2009, hereafter SQ09) to calculate the emission properties of Pop III star TDEs, while making a few changes to adapt model parameters consistent with more recent studies.

We first provide a brief summary of the SQ09 model below. SQ09 assume that during the early, super-Eddington phase a sizable fraction of the fallback material constitutes an outflow. Furthermore, they also assume that a fraction of the wind kinetic energy is converted to thermal energy, which in turn produces blackbody radiation. This assumption appears to hold even when the outflows are powered by debris stream collisions, and it is supported by recent simulations (Dai et al. 2018; Zanazzi & Ogilvie 2020; Andalman et al. 2022; Thomsen et al. 2022). SQ09 assume a spherical wind geometry and a constant wind velocity; under these circumstances, the wind density profile can be approximated as

$$\begin{eqnarray}\rho (r)=\left\{\begin{array}{ll}{\dot{M}}_{\mathrm{wind}}/(4\pi {r}^{2}{v}_{\mathrm{wind}}), & r\lesssim {R}_{\mathrm{edge}}\\ 0, & \mathrm{outside}\end{array}\right..\end{eqnarray} \tag{ 12 }$$

Here the wind mass rate ${\dot{M}}_{\mathrm{wind}}$ is assumed to be a constant fraction of the debris mass fallback rate so that ${\dot{M}}_{\mathrm{out}}\equiv {f}_{\mathrm{out}}{\dot{M}}_{\mathrm{fb}}$. We use a fiducial value f_out = 0.5 (instead of f_out = 0.1 used by SQ09), inspired by recent numerical simulations of super-Eddington accretion flows (e.g., Dai et al. 2018; Jiang et al. 2019; Thomsen et al. 2022). We recognize that the exact value of f_out has a dependence on the Eddington ratio, but our results do not significantly vary when f_out > 0.5. R_edge ≡ v_wind t denotes the edge of the wind. The wind velocity v_wind is assumed to be the escape velocity ${v}_{\mathrm{esc},L}\equiv \sqrt{2{{GM}}_{\mathrm{BH}}/{R}_{L}}$ at the wind-launching radius R_L , which is set to be 2R_p (the circularization radius). Therefore, we have ${v}_{\mathrm{esc},L}=\sqrt{{{GM}}_{\mathrm{BH}}/{R}_{p}}$, which can be further simplified as

$$\begin{eqnarray}&&{v}_{\mathrm{esc},L}=0.1{\beta }^{1/2}{m}_{\star }^{1/6}{r}_{\star }^{-1/2}{M}_{6}^{1/3}c.\end{eqnarray} \tag{ 13 }$$

Furthermore, SQ09 calculate the temperature of the wind by assuming that the gas thermal energy density and kinetic energy density are similar at the wind-launching site:

$$\begin{eqnarray}&&{{aT}}_{L}^{4}=\displaystyle \frac{1}{2}{\rho }_{\mathrm{fb},L}{v}_{\mathrm{esc},L}^{2}.\end{eqnarray} \tag{ 14 }$$

Here a = σ/4, with σ being the Stefan–Boltzmann constant, T_L the gas temperature at R_L , and ρ_fb,L the gas density at R_L , which can be calculated as ${\rho }_{\mathrm{fb},{\rm{L}}}\equiv {\dot{M}}_{\mathrm{fb}}/(4\pi {R}_{L}^{2}{v}_{\mathrm{esc},L})$. It is further assumed that the wind expands adiabatically so that its temperature scales with gas density:

$$\begin{eqnarray}&&T\propto {\rho }^{1/3}.\end{eqnarray} \tag{ 15 }$$

The photosphere radius of the wind, R_ph, is located where

$$\begin{eqnarray}&&{R}_{\mathrm{ph}}{\kappa }_{s}\rho ({R}_{\mathrm{ph}})\sim 1,\end{eqnarray} \tag{ 16 }$$

where κ_s is the electron scattering opacity, which is taken to be 0.35 cm² g⁻¹ assuming the hydrogen abundance for Pop III stars is ∼75%, while the rest of the mass is mostly dominated by helium (Bromm et al. 2009).

Initially, the wind gas density is very high due to the large fallback rate, so R_edge κ_s ρ(R_edge) ≫ 1 and the photosphere almost coincides with the edge of the wind, which gives

$$\begin{eqnarray}&&{R}_{\mathrm{ph}}\approx {R}_{\mathrm{edge}}=1.3\times {10}^{3}{f}_{v}{\beta }^{1/2}{m}_{\star }^{1/6}{r}_{\star }^{-1/2}{M}_{6}^{-2/3}\left(\displaystyle \frac{t}{\mathrm{day}}\right)\,{R}_{{\rm{s}}},\end{eqnarray} \tag{ 17 }$$

where R_s is the Schwarzschild radius of a nonspinning black hole. The effective temperature of the photosphere at this stage is given by

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{ph}} & \approx & 1.5\times {10}^{4}{f}_{v}^{-1/3}{m}_{\star }^{-7/72}{r}_{\star }^{1/24}{M}_{6}^{5/36}\\ & & \times {\left(\displaystyle \frac{t}{\mathrm{day}}\right)}^{-7/36}{\left(\displaystyle \frac{{t}_{\mathrm{fb}}}{\mathrm{yr}}\right)}^{-1/18}{\rm{K}}.\end{array}\end{eqnarray} \tag{ 18 }$$

As ${R}_{\mathrm{edge}}{\kappa }_{s}\rho ({R}_{\mathrm{edge}})\propto {\dot{M}}_{\mathrm{wind}}/{R}_{\mathrm{edge}}$ decreases with t, after a certain time R_edge κ_s ρ(R_edge) drops to 1 and the photosphere starts to recede afterwards. This transition time is calculated to be

$$\begin{eqnarray}&&{t}_{\mathrm{edge}}\approx 9.6{f}_{\mathrm{out},0.5}^{3/8}{f}_{v}^{-3/4}{B}_{\gamma }^{1/4}{M}_{6}^{1/8}{R}_{p,3{R}_{{\rm{s}}}}^{3/8}{m}_{\star }^{1/8}{r}_{\star }^{3/8}\,\mathrm{day},\end{eqnarray} \tag{ 19 }$$

where f_out,0.5 ≡ f_out/0.5, and f_v = v_wind/v_esc,L (this fraction is set as 1 throughout this Letter unless otherwise specified). As we will show later, t_edge is also the time when the flare bolometric luminosity reaches the peak. We plot t_edge for Pop III stars and the Sun in Figure 5. It is clearly seen that the transition occurs on timescales of a few tens of days for Pop III stars as compared to a few days for the Sun in case of M_BH = 10⁶ M_⊙. This leads to larger photosphere radii and higher peak luminosities for Pop III TDEs.

Zoom In Zoom Out Reset image size

After the photosphere starts to recede, the photosphere radius and temperature at this later stage can be estimated following Equations (15) and (16), and given as

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{ph}} & = & {\kappa }_{s}{\dot{M}}_{\mathrm{out}}/(4\pi {v}_{\mathrm{wind}})\\ & \approx & 5.4{f}_{\mathrm{out},0.5}{f}_{v}^{-1}\left(\displaystyle \frac{{\dot{M}}_{\mathrm{fb}}}{{\dot{M}}_{\mathrm{Edd}}}\right){R}_{p,3{R}_{{\rm{s}}}}^{1/2}{R}_{{\rm{s}}},\end{array}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{ph}}\approx 2.5\times {10}^{5}{f}_{\mathrm{out},0.5}^{-1/3}{f}_{v}^{1/3}{\left(\displaystyle \frac{{\dot{M}}_{\mathrm{fb}}}{{\dot{M}}_{\mathrm{Edd}}}\right)}^{-5/12}{M}_{6}^{-1/4}{R}_{p,3{R}_{{\rm{s}}}}^{-7/24}\,{\rm{K}},\end{eqnarray} \tag{ 21 }$$

where ${R}_{p,3{R}_{{\rm{s}}}}={R}_{p}/3{R}_{{\rm{s}}}$.

Specifically, for Pop III stars, using their mass–radius relationship (Equation (1)), we have

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{edge}} & \simeq & 85{f}_{v}^{-3/4}{\left(\displaystyle \frac{{f}_{\mathrm{out}}}{0.5}\right)}^{3/8}{\left(\displaystyle \frac{{B}_{\gamma }}{0.08}\right)}^{1/4}{\left(\displaystyle \frac{\beta }{1.85}\right)}^{-3/8}\\ & & \times \,{\left(\displaystyle \frac{{m}_{\star }}{300}\right)}^{15/44}{\left(\displaystyle \frac{Z}{{10}^{-5}}\right)}^{3/44}{M}_{6}^{-1/8}\ \mathrm{day}.\end{array}\end{eqnarray} \tag{ 22 }$$

Before t_edge:

$$\begin{eqnarray}&&{R}_{\mathrm{ph}}\propto {M}_{6}^{1/3}{m}_{\star }^{-2/33}{z}_{\star }^{-1/22},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{ph}}\propto {M}_{6}^{1/9}{m}_{\star }^{-2/33}{z}_{\star }^{-1/264}.\end{eqnarray} \tag{ 24 }$$

After t_edge:

$$\begin{eqnarray}&&{R}_{\mathrm{ph}}\propto {m}_{\star }^{28/33}{z}_{\star }^{3/22},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{T}_{\mathrm{ph}}\propto {M}_{6}^{2/9}{m}_{\star }^{-4/11}{z}_{\star }^{-17/264}.\end{eqnarray} \tag{ 26 }$$

We note that this first-order calculation in SQ09 assumes that the outflows start to be launched around the peak of the debris mass fallback rate and they ignore the outflows produced prior. The luminosity increases in the initial phase as the wind builds up, and therefore the characteristic rising time for the flare will be t_edge instead of t_fb. While this assumption can cause issues for a solar-type star TDE which has t_edge ≲ t_fb, it works well for Pop III star TDEs in which t_edge ≫ t_fb. Therefore, Pop III star TDE flares should rise on timescales of t_edge, i.e., a few weeks to months in the intrinsic TDE frame.

Furthermore, SQ09 also consider disk emission after the fallback rate drops to the sub-Eddington level. At this phase, the gas can radiate efficiently and a thin accretion disk is expected to form, which produces emission mainly in the X-ray and EUV bands. However, for Pop III star TDEs, as seen from Figure 4, the fallback rate stays super-Eddington for a very long time (tens to hundreds of years). Therefore, in this work we focus only on the emission produced from the winds launched in the super-Eddington phase and ignore the disk emission at very late times.

We first calculate the TDE luminosity and spectral energy distribution (SED) in the rest frame of the host galaxy based on the properties of the photosphere described in Section 2.3. Assuming blackbody emission, the energy spectrum and bolometric luminosity can be computed using

$$\begin{eqnarray}&&\nu {L}_{\nu }=4{\pi }^{2}{R}_{\mathrm{ph}}^{2}\nu {B}_{\nu }({T}_{\mathrm{ph}}),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&L=4\pi \sigma {R}_{\mathrm{ph}}^{2}{T}_{\mathrm{ph}}^{4}.\end{eqnarray} \tag{ 28 }$$

We plot the emission SEDs in Figure 6 and explore the dependence on various model parameters. The default set of parameters for a fiducial Pop III star TDE model is M_BH = 10⁶ M_⊙, M_⋆ = 300M_⊙, Z = 10⁻⁵, f_out = 0.5, f_v = 1, and β = β_c = 1.85. We stick to these parameters for calculations unless specified otherwise.

Zoom In Zoom Out Reset image size

We illustrate the evolution of the emission from a single TDE in Figure 6(a), which shows the SEDs at different epochs for the fiducial case. It is seen that initially the luminosity increases and the peak of the SED shifts slightly toward lower frequency as a result of increasing R_ph and decreasing T_ph during this phase. After t_edge, R_ph recedes while T_ph increases, so the SED evolves in reverse, i.e., the luminosity decreases and the SED peak shifts toward higher frequency.

Figure 6(b) shows the SEDs at t = t_edge (the peak of the flare) from the TDEs of Pop III stars of different masses and metallicities in comparison to that of a solar-type star, all around a 10⁶ M_⊙ MBH. It can be immediately noticed that a significant fraction of the emission energy resides in the UV/optical wavelength regimes for all cases. The SEDs of Pop III star TDEs, compared to that of the solar-type TDE, have lower peak temperatures (with the SED peaking in the optical band instead of the UV) and slightly larger luminosities. However, the mass or metallicity of a Pop III star does not make a significant difference to its tidal flare emission. Increasing the metallicity leads to a slightly smaller luminosity and increasing the stellar mass shifts the SED toward slightly longer wavelengths.

Figure 6(c) shows the impact of M_BH on the SED. It can be seen that more massive black holes produce more luminous flares. Interestingly, the recent detection of sources such as UHZ1 (Bogdán et al. 2023) and GNz-11 (Maiolino et al. 2024) suggest that MBHs of M_BH ≳ 10⁶⁻⁷ M_⊙ at z ≥ 10 could be much more abundant than previously estimated. Therefore, it is likely that some Pop III TDEs can produce very luminous flares with intrinsic bolometric luminosity larger than 10⁴⁵ erg s⁻¹.

Finally, Figure 6(d) shows how various other parameters, namely, f_out (outflow fraction), f_v (ratio between wind velocity and escape velocity), and β (stellar orbital penetration parameter), affect the SED. One can see that overall the choice of these parameters only mildly impacts the SED. We note that increasing f_out or decreasing f_v generally both lead to a larger peak luminosity. A lower β means the star is only partially disrupted, which will reduce the outflow mass, so the effect induced is similar to having a lower f_out. Increasing β beyond β_c only slightly increases the peak fallback rate so it barely also affects the peak luminosity of the flare.

Next, we consider the evolution of the flare flux in specific wavelength bands. Figure 7 shows the optical (430–750 THz) and UV (750–3 × 10⁴ THz) light curves of TDEs of Pop III and solar-type stars by a MBH of M_BH = 10⁶ M_⊙. For the fiducial Pop III star model, we also mark the epochs (t = t_edge/2, t_edge, 2t_edge, 4t_edge) on the light curves. One can see that both the optical and UV luminosities increase in the initial phase when R_ph traces the edge of the wind. However, the behavior of the UV and optical light curves after t_edge are different. For a solar-type star, since the peak of the flare SED stays in the UV regime, both the UV and the optical light curve reach the peak around t_edge and decay afterwards. However, for a Pop III star, the TDE flare SED shifts from UV to optical bands around t_edge and shifts back to the UV band afterwards. Therefore, while the tidal flare optical light curve still peaks at t_edge, the UV light curve continues to rise for tens to few hundreds of days after t_edge. Furthermore, at late times, while the optical luminosity decays rather closely following the debris mass fallback rate (∝t^−5/3), the UV light curve has a much shallower slope (∝t^−0.57) due to the temperature evolution in this phase. This also means that the UV light curve has an evolution timescale longer than that of the optical light curve.

Zoom In Zoom Out Reset image size

In summary, the bolometric luminosity of a Pop III star TDE in the rest frame of the host galaxy increases with a larger M_BH, a smaller metallicity, a higher outflow fraction or a slower wind velocity, which all likely promote the potential observability of these TDEs. Furthermore, a higher Pop III mass will additionally serve to lengthen the evolution timescale of these light curves. We calculate the properties of the observed and redshifted TDE emission in the next section.

So far, we have established that the emission from Pop III TDEs mostly resides in UV/optical wave bands in the rest frame of the TDE host galaxy. In the subsequent calculations, we adopt a canonical value of redshift z = 10 for Pop III star TDEs and calculate the redshifted TDE emissions observed at z = 0.

With ν_e and ν_o denoting the rest-frame frequency and the observed frequency, respectively, the two are related as

$$\begin{eqnarray}&&{\nu }_{{\rm{e}}}=(1+z){\nu }_{{\rm{o}}}.\end{eqnarray} \tag{ 29 }$$

Therefore, the observed Pop III star TDE SED should be redshifted following ν_o = ν_e /(1 + z) and peaks around or slightly below 10¹⁴ Hz, which corresponds to λ ∼ 10³ nm. This means that a large fraction of the Pop III star TDE emission should be redshifted to the near-infrared (NIR) regime, and such events can be potentially detected by the NIRCam on JWST covering the wavelength range of 600–5000 nm as well as the WFI on Roman with a wave band of 480–2300 nm.

The observed specific flux, F_ν (ν_o), is related to the rest-frame specific luminosity, L_ν (ν_e), following (Hogg et al. 2002)

$$\begin{eqnarray}&&{F}_{\nu }({\nu }_{{\rm{o}}})=\displaystyle \frac{1+z}{4\pi {D}_{L}^{2}}{L}_{\nu }({\nu }_{{\rm{e}}}).\end{eqnarray} \tag{ 30 }$$

Here D_L denotes the luminosity distance, which can be expressed as

$$\begin{eqnarray}&&{D}_{L}=(1+z){D}_{c},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray*}&&{D}_{c}=\displaystyle \frac{c}{{H}_{0}}{\int }_{0}^{z}\displaystyle \frac{{dz}}{\sqrt{{{\rm{\Omega }}}_{M}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}}\end{eqnarray*}$$

is the comoving distance, and Ω_M , Ω_Λ, H₀ is the matter, energy density parameter, and Hubble constant at the current time, respectively.

We calculate the observed fluxes of Pop III TDEs at different epochs throughout the event for different stellar models using Equation (30). The results are plotted in Figure 8. It can be seen that for all models the fluxes in the NIR band stay above the detection limit of the JWST NIRCam (∼11 nJy, corresponding to the F150W filter) and Roman WFI (∼20.9 nJy, corresponding to F106 filter) throughout these epochs. Also, no significant differences in the flux level are observed due to the different masses or metallicities of Pop III stars. Moreover, we note that the observed flux from the TDE of a disrupted solar-type star is also detectable by NIRcam and WFI.

Zoom In Zoom Out Reset image size

Additionally, we have examined the dependence of observed flux on other model parameters, and the results are shown in Figure 9. Consistent with the results illustrated in Figures 6(c) and (d), one can see that the NIR flux increases significantly with increasing M_BH and moderately with decreasing f_v , but is barely affected by f_out and β. In particular, the NIR flux becomes undetectable even around peak luminosity when M_BH < 10⁵ M_⊙, unless the wind is exceptionally slow (f_v ≲ 0.1). In summary, in all models, the observed Pop III star TDE NIR flux around flare peak stays significantly above the detection limit of NIRCam and WFI, as long as M_BH ≳ 10⁶ M_⊙. Therefore, our results demonstrate that the detection of Pop III star TDEs at z ∼ 10 with both JWST or Roman is feasible.

Zoom In Zoom Out Reset image size

Next, we calculate the NIR light curves of Pop III star TDEs as would be observed in the JWST NIRcam band. The NIR light curves are shown in Figure 10, where one can see that the NIR flux can exceed both the JWST NIRCam and Roman WFI detection limits even when the TDE originates at z = 15. Moreover, at such high redshifts the time-dilation effect due to cosmological redshift is very strong, since t_o = (1 + z)t_e. Taking a Pop III star TDE at z = 10 as an example, the observed NIR light curve (above the detection limit) rises on timescales ∼100–1000 days and decays over 10³–10⁴ days. At late time, the NIR flux still decays following a pattern relatively close to the classical mass fallback rate (∝t^−5/3).

Zoom In Zoom Out Reset image size

Furthermore, inspired by the recent discovery of UHZ1, we include the case of a Pop III star disrupted by a MBH with M_BH = 4 × 10⁷ M_⊙ hosted by galaxy at z = 10.1 (orange line). As expected, the observed NIR flux is significantly higher in this case, which increases the chance of observing such events in current and upcoming NIR surveys.

Interestingly, given the observed Pop III star TDE flares rise over a few years, there is still a prospect of identifying such TDEs detected during the rising phase as transients, if multiple detections are made within the typical operation time of surveys. However, during the flare-decay phase, as the evolution timescale is very long these events will have almost near-constant brightness during the limited operation time, and therefore they are more likely to be miscategorized as active galactic nuclei (AGNs) based on their photometry. One promising method for distinguishing Pop III star TDEs from the TDEs of Pop I stars or AGNs is through spectroscopic follow-ups. Metal lines should be present in the spectra of either Pop I star TDEs or AGNs even at z ∼ 10 given that gas is likely no longer metal-free (Yang et al. 2023). However, if a Pop III star TDE occurs in a previously quiescent galaxy, a metal-free spectrum should be produced, since all gas supplied to the MBH likely originates from the tidally disrupted Pop III star.

We now estimate the intrinsic rates of Pop III TDEs in the early Universe. We adopt the approach in Pfister et al. (2020, 2022) to calculate the TDE rate through two-body scattering. The differential TDE rate is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{d}^{2}{\rm{\Gamma }}}{d\mathrm{log}{m}_{\star }d\mathrm{log}\beta } & = & 8{\pi }^{2}{{GM}}_{\mathrm{BH}}\displaystyle \frac{{R}_{T}}{\beta }{\rm{\Phi }}({m}_{\star }){m}_{\star }\\ & & \times \,{\displaystyle \int }_{0}^{{E}_{m}}\zeta (q,f,\beta ,{m}_{\star })\,{dE}.\end{array}\end{eqnarray} \tag{ 32 }$$

Here Φ(m_⋆) is the stellar IMF, E is the specific orbital energy of the star with the maximum value at E_m = GM_BH/R_T corresponding to the orbit with R_T as the closest approach, and ζ(q, f, β, m_⋆) is a function of the loss-cone filling factor q, the stellar distribution function f, the stellar mass m_⋆, and the penetration parameter β. The TDE of a Pop III star can be induced by the scattering between a Pop III star and a normal star or that between two Pop III stars. However, at the redshift range that we are considering, Pop III stars contribute to only a few percent of total stellar mass (Magg et al. 2022). Therefore, we only include the scattering of Pop III stars by normal stars into our calculation. Further details of the terms in Equation (32) and the rate calculation are given in Appendix B.

We summarize the assumptions used for the TDE rate calculation as below:

• 1.  
  Ranges of parameters: we adopt a MBH mass range of M_BH = [10⁵ M_⊙, 10⁸ M_⊙], a Pop III stellar mass range of m_⋆ = [30M_⊙, 300M_⊙], a normal star in the range of m_⋆ = [0.1M_⊙, 10M_⊙], and a β in the range [β_c = 1.85, R_T /R_s] to calculate the total TDE rate by integrating Equation (32).
• 2.  
  IMF of the stellar populations: we adopt the general Kroupa IMF (Kroupa 2001) for the normal main-sequence stars. However, the IMF for Pop III stars is model dependent and not constrained very well yet. Different IMFs have been proposed in the literature, and a few examples are shown in Table 1 and Figure 11. However, using different Pop III star IMFs brings negligible effects to TDE rates, as explained in Appendix B.
• 3.  
  Mass fraction of Pop III stars: Magg et al. (2022) showed that at z ∼ 10 Pop III stars have a mass fraction between 1% and 15% among the entire stellar population depending on the timescale at which Pop III stars transit to Pop II stars. We adopt a Pop III star mass fraction of 7% throughout the galaxy for our calculation.
• 4.  
  Stellar density distribution: as there is a lack of studies of the structures of galaxies for different stellar populations at these extremely high redshifts, we follow classical papers such as those by Wang & Merritt (2004), Stone & Metzger (2016), and Pfister et al. (2020) and assume isothermal stellar density distributions in the Keplerian potential for both Pop III and normal stars, which gives

  $$\begin{eqnarray}&&\rho (r)={\rho }_{0}{\left(r/{r}_{\inf }\right)}^{-\alpha },\end{eqnarray} \tag{ 33 }$$

  where ρ(r) is the stellar density, α = 2, and ${r}_{\inf }\,={{GM}}_{\mathrm{BH}}/2{\sigma }^{2}$, with σ being the velocity dispersion of the host galaxy.

Zoom In Zoom Out Reset image size

Table 1. Different IMFs of Pop III Stars Suggested by Different Works in the Literature

IMF	Φ(m⋆)	α	mc	Γ	${\dot{N}}_{\max }$
				(gal−1 yr−1)	(Mpc−3 yr−1)
Greif et al. (2011)	$\propto {m}_{\star }^{-\alpha }$	−0.17	⋯	1.8 × 10−8	1.1 × 10−9
de Bennassuti et al. (2014)	$\propto {m}_{\star }^{\alpha -1}\exp (-{m}_{c}/{m}_{\star })$	−1.35	20M⊙	4.1 × 10−8	2.5 × 10−9
Jaacks et al. (2018)	$\propto {m}_{\star }^{-\alpha }\exp (-{m}_{c}^{2}/{m}_{\star }^{2})$	0.17	4.47M⊙	2.0 × 10−8	1.2 × 10−9

Notes. Φ(m_⋆) describes the IMF equation, with α and m_c denoting the slope and cutoff mass, respectively. Γ represents the rate of Pop III star TDEs from a single galaxy with M_BH = 10⁶ M_⊙. ${\dot{N}}_{\max }$ represents the Pop III star TDE volumetric rate computed using the BHMF from the the merger-driven model.

Download table as:  ASCIITypeset image

Under these assumptions, the TDE rate of Pop III stars from a single galaxy hosting a 10⁶ M_⊙ MBH is calculated and listed in Table 1. This intrinsic TDE rate is Γ ∼ 10⁻⁸ gal⁻¹ yr⁻¹, which is barely affected by the choice of Pop III star IMF.

One can further compute the volumetric rate of Pop III star TDEs using

$$\begin{eqnarray}\begin{array}{rcl}\dot{N} & = & \iiint \displaystyle \frac{{dn}}{d\mathrm{log}{M}_{\mathrm{BH}}}\times \displaystyle \frac{{d}^{2}{\rm{\Gamma }}}{d\mathrm{log}{m}_{\star }d\mathrm{log}\beta }\\ & & \times \,d\mathrm{log}{m}_{\star }\,d\mathrm{log}\beta \,\,d\mathrm{log}{M}_{\mathrm{BH}},\end{array}\end{eqnarray} \tag{ 34 }$$

here ${\rm{\Phi }}({M}_{\mathrm{BH}})={dn}/d\mathrm{log}{M}_{\mathrm{BH}}$ denotes the black hole mass function (BHMF), which is defined as the number of MBHs with masses between $\mathrm{log}{M}_{\mathrm{BH}}$ and $\mathrm{log}{M}_{\mathrm{BH}}+d\mathrm{log}{M}_{\mathrm{BH}}$ in unit comoving volume. We note that the BHMF at high redshifts sensitively depends on the seeding scenarios and growth channels of MBHs, much of which is currently uncertain. In this work, we adopt three different BHMFs from Trinca et al. (2022) that consider not only light, intermediate, and heavy MBH seeds but also different processes responsible for their growth. Figure 12 shows the BHMF at z = 10 for these three growth models, namely, Eddington-limited Bondi–Hoyle–Lyttleton accretion, super-Eddington accretion, and merger-driven growth (Hoyle & Lyttleton 1941; Bondi 1952). We focus on the MBHs with M_BH ≥ 10⁵ M_⊙ for this work, since the observed Pop III star TDE flux goes below the threshold of NIRCam and WFI when M_BH < 10⁵ M_⊙ (Figure 9). It is seen that the BHMF from the merger-driven growth channel generally has higher values in this M_BH range than the other two BHMF choices. The volumetric Pop III star TDE rates using these three different BHMFs are calculated and shown in Figure 13, which are around 10⁻¹⁰–10⁻⁹ Mpc⁻³ yr⁻¹. Out of the three, not surprisingly, the BHMF based on the merger-driven growth model produces the highest integrated volumetric TDE rate ${\dot{N}}_{\max }\approx {10}^{-9}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also include the BHMF inferred from calibrating models with the detection of UHZ1 in Figure 12, and note that the BHMF taken from Trinca et al. (2022; as well as other current BHMF models) underestimates the number density of MBHs in all the growth models. Hence, it is possible that the actual Pop III TDE rate is higher than our estimate. An updated BHMF at high redshift from additional observations of a population of z ∼ 10 black holes will be crucial for more accurately constraining the Pop III star TDE rate.

In this section, we estimate the total number of Pop III star TDEs that stand to be detected by the JWST and Roman telescopes. For a first-order estimate of the upper limit of the detection numbers, we ignore factors such as survey strategies and limitations. The total number can then be approximated using

$$\begin{eqnarray}&&N=\dot{N}\times V\times T.\end{eqnarray} \tag{ 35 }$$

Here $\dot{N}$ is the volumetric TDE rate obtained in the previous section, T is the duration of the observational survey, and V is the comoving volume from which the event can be detected. V can be computed using the following equation:

$$\begin{eqnarray}&&V=\displaystyle \frac{4\pi [d{\left({z}_{\max }\right)}^{3}-d{\left({z}_{\min }\right)}^{3}]}{3}\times {f}_{\mathrm{sky}},\end{eqnarray} \tag{ 36 }$$

where d(z) is the comoving distance dependent on redshift z, and f_sky is the sky coverage fraction of the telescope of survey.

For the case of Pop III star TDEs, we adopt ${z}_{\min }=10$ when the Pop III star population reaches peak, and ${z}_{\max }=15$ when both the Pop III stars and MBHs should have started to form. More importantly, we have predicted in Figure 10 that Pop III star TDEs at z = 15 are still bright enough to be detected by JWST and Roman. Moreover, for simplicity, we stick to a constant $\dot{N}=1.2\times {10}^{-9}\ {\mathrm{Mpc}}^{-3}\ {\mathrm{yr}}^{-1}$ (based on the Pop III star IMF by Jaacks et al. 2018 and the BHMF from the merger-driven growth model), and note that this likely gives an upper limit to the detected number.

JWST Cycle 1 has various deep surveys, such as COSMOS-Web with a total area of 0.54 deg² and duration of 255 hr (Casey et al. 2023), and the JWST Advance Deep Extragalactic Survey (JADES) with its sky coverage of $236\,{\mathrm{arcmin}}^{2}$ and a total duration of 426 hr considering both the deep and medium modes (Eisenstein et al. 2023). Using both the COSMOS-Web survey and the deep mode of the JADES survey, which can both probe galaxies at z > 10, the total detected number of Pop III TDEs in a year is about 5 × 10⁻⁴. Throughout the 10 yr expected lifespan of JWST, even if we assume COSMOS-Web–like surveys are continuously conducted, the total expected detection number is still only (10 yr/255 hr) × 5 × 10⁻⁴ ∼ 0.2 events in total. Therefore, the chance of detecting Pop III star TDEs using JWST is slim. If such TDEs were detected, it would indicate the BHMF at z > 10 is much larger than currently estimated, which would pose an interesting challenge in our understanding of the formation efficiency and evolution of MBHs in the early Universe. With the recent discovery of UHZ1 and its interpretation as an OBG arising from a heavy initial seed by Natarajan et al. (2024), it appears that black hole seed formation in the early Universe can occur via multiple seeding pathways, hence rendering the process significantly more efficient than previously believed.

A much more promising telescope for detecting Pop III star TDEs is the upcoming Roman, which is designed to have a sensitivity similar to JWST (Figure 10) and to conduct wide-field surveys over a very large FOV of 0.281 deg² (Mosby et al. 2020). Hence, from the High Latitude Wide Area Roman survey, which aims to cover a total sky area of ∼2000 deg² (Wang et al. 2022b), the number of detected Pop III TDEs can reach ≲60 in a year. Furthermore, there is another proposed survey, the Next-generation All-sky Near-infrared Community surveY (or NANCY; Han et al. 2023), which plans to perform an all-sky scan. This survery strategy will lead to the detection of ∼1000 Pop III TDEs in a single year. Once a large sample of Pop III star TDEs has been observed, the number can be used to put a strong constraints on the properties and mass fraction of Pop III stars as well as the BHMF at z ≳ 10, which will provide the much needed insights into the efficiency of early black hole formation.