{'fraction_non_alphanumeric': 0.031382299576641355, 'fraction_numerical': 0.01397755527182313, 'mean_word_length': 4.51705388037568, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Triangulating Black Hole Forming Stellar Collapses through Neutrinos

In the event of a black hole (BH) forming stellar collapse, the neutrino signal should terminate abruptly at the moment of BH formation, after a phase of steady accretion. Since neutrinos are expected to reach Earth hours before the electromagnetic signal, the combined detection of the neutrino burst through multiple neutrino telescopes could allow to promptly determine the angular location of a nearby stellar collapse in the sky with high precision. In this paper, we contrast the triangulation pointing procedure that relies on the rise time of the neutrino curve, often considered in the literature, to the one that takes advantage of the cutoff of the neutrino curve at the moment of BH formation. By forecasting the neutrino signal expected in the IceCube Neutrino Observatory, Hyper-Kamiokande and DUNE, we devise a strategy to optimize the identification of the rise and cutoff time of the neutrino curve. We show that the triangulation method developed by employing the end tail of the neutrino curve allows to achieve at least one order of magnitude improvement in the pointing precision for a galactic burst, while being insensitive to the neutrino mixing scenario. The triangulation pointing method based on the cutoff of the neutrino curve will also guarantee a better performance for BH forming collapses occurring beyond our own Galaxy.


I. INTRODUCTION
Being copiously produced, neutrinos carry key information about and strongly affect the physics of the latest stages of life of stars with a zero-age main sequence mass between 10 and 150 M [1][2][3][4]. For core-collapse supernova explosions, the bounce of the electron-degenerate iron core induces a shock wave in the infalling stellar matter that is revived by neutrinos, after stalling [5][6][7][8][9]. However, if the explosion mechanism fails and matter continues to accrete onto the transiently stable protoneutron star, a stellar-mass black hole (BH) forms.
Recent theoretical work [10][11][12][13][14] as well as the observation of disappearing red giants [15][16][17][18][19][20] hint towards the possibility that BH formation may occur for 10-30% of all massive stars undergoing core collapse, and possibly more. Such collapses are expected to emit a faint electromagnetic signal due to the stripping of the hydrogen envelope [21][22][23], hence neutrinos (and probably gravitational waves) could be the only messengers of such transient events.
The early detection of neutrinos is especially important in aiding the prompt reconstruction of directional information about the stellar collapse. To this purpose, the neutrino signal could be exploited for pointing and triangulation [34][35][36][37][38][39][40][41][42][43][44]. In particular, in water Cherenkov neutrino detectors, such as Super-Kamiokande [45] and the 1 Neutrino curve is used in analogy to light curve in this work.
upcoming Hyper-Kamiokande [46], the electron-neutrino elastic scattering can be employed to determine the angular location of the stellar collapse, allowing to achieve a few-degree pointing precision [34,47]. The upgrade of Super-Kamiokande with gadolinium foresees an increased tagging efficiency that will improve the pointing precision by 50%, while a sub-degree precision could be potentially reached with Hyper-Kamiokande [35]. A similar performance is expected from DUNE [43,48]. Triangulation by computing the time delays of the neutrino signal from standard core collapse supernovae among various neutrino telescopes does not allow to obtain such a precision. However, for BH forming collapses, an improved pointing precision of triangulation is guaranteed because of the higher neutrino event statistics achievable with respect to standard core-collapse supernovae [37,39,40].
In this paper, we investigate the possibility of improving on the determination of the angular location of a BH forming collapse. We work under the assumption of a negligible or absent supersonic accretion flow in the surroundings of the collapsing proto-neutron star, leading to the formation of a naked BH. In fact, in such a scenario, one should expect a sharp drop of the neutrino signal [25][26][27][28], with no smearing of the tail determined by the coherent scattering of neutrinos off nuclei, possibly occurring in the event of a non-negligible accretion flow [31]. We show that, by focusing on the tail of the neutrino curve instead of on its bounce time to compute the time delays, one can achieve at least an order of magnitude better determination of the angular location of the stellar collapse because of the larger event rate expected at the end of the accretion phase, which is more easily distinguishable from the background signal.
This paper is organized as follows. In Sec. II, we introduce the neutrino emission properties for our BH forming collapse model. The main features of the neutrino tele-scopes adopted in this work are outlined in Sec. III. The method adopted to determine the timing of the neutrino signal is introduced in Sec. IV, a comparison between the determination of the time delay among pairs of detectors is discussed by employing the rise time as well as the cutoff time of the neutrino curve. In Sec. V, our findings on the triangulation of the neutrino signal are presented. The dependence our findings on the distance at which the BH forming collapse occurs can be found in Sec. VI, followed by a discussion on our main results and caveats in Sec. VII and concluding remarks in Sec. VIII.

II. NEUTRINO EMISSION PROPERTIES FOR BLACK HOLE FORMING COLLAPSES
The neutrino emission properties for our benchmark BH forming collapse are extracted from a onedimensional (1D) spherically symmetric hydrodynamical simulation of BH forming collapse from the Garching group [3,49]: the 40 M model with slow accretion rate (s40c) for which the gravitational instability occurs after 2100 ms. The employed high density nuclear equation of state is the Lattimer and Swesty one with nuclear incompressibility modulus K = 220 MeV [50]. We refer the interested reader to Sec. 2.5 of Ref. [3] for more details on this model.
It is worth to be noted that BH forming collapses have been simulated in multi-D, see e.g. Refs. [51][52][53][54][55], and 3D hydrodynamical simulations with sophisticate neutrino transport exhibit strong signatures of hydrodynamical instabilities in the neutrino curve [54]. However, since we are interested in the rise time of the neutrino curve as well as in its drop, and these features are not affected by the dimensionality of the simulation, a 1D model suits our purpose well.
Dedicated simulations following the dynamical implosion of the transiently stable neutron star into a naked BH foresee two different scenarios: either the naked or low accreting neutron star collapses into a BH by cooling or a softening of the nuclear equation of state leads to a phase transition where exotic states of matter form [25][26][27][28][29][30]. In the first case, the neutrino signal terminates abruptly at the moment of BH formation, whereas in the second scenario a gradual decrease of the neutrino signal is foreseen as more and more matter approaches the event horizon.
An alternative scenario with respect to the one of naked BH formation described above foresees a substantial accretion flow surrounding the collapsing protoneutron star. Under the assumption of spherical accretion, the arrival time of neutrinos may be affected by the coherent scattering of neutrinos in the infalling medium surrounding the BH [31]; this would result in a longer tail of the neutrino curve (so-called echo).
Triangulation through neutrinos could be carried out both in the case of the formation of a naked BH and when a non-negligible accretion affects the tail of the neutrino curve. However, in this work, we choose to focus on the scenario of naked BH formation since we are interested in quantitatively exploring the best circumstances under which triangulation could be carried out for BH forming collapses.
Our benchmark model has an instantaneous cutoff of the neutrino emission properties at the time of BH formation. Since we are interested in the tail of the neutrino curve, we assume that the signal predicted in Ref. [28] is representative of naked BH formation scenarios. In fact, Ref. [28] followed the hyperonization of the cooling core during the collapse of a protoneutron star to a BH.
We take into account a linear tail in the neutrino signal due to the increasing gravitational redshift lasting for approximately 0.5 ms [28] The non-radial neutrino emission may contribute to soften the cutoff of the neutrino curve [56], adding an exponential contribution at the end of the neutrino tail; we assume that this effect is also enclosed in the representative tail of the neutrino curve. In addition, the contribution from neutrinos on unstable circular orbits is expected to decay exponentially with a decay time of O(0.05) ms [28,57], however this subleading component is neglected. We will return on the impact of our modeling assumptions in Sec. VII.
By relying on the inputs from our benchmark BH stellar collapse model, for each flavor ν β (with β = e, µ or τ and ν µ = ν τ = ν x =ν x hereafter) and post-bounce time t, we compute the neutrino differential flux for a stellar collapse occurring at a distance d as follows: with E being the neutrino energy, L ν β the luminosity of the neutrino flavor ν β , and E ν β the neutrino mean energy. The neutrino energy distribution is [58,59]: (2) where ξ ν β (t) is a normalization factor [ dEφ ν β (E, t) = 1] and the pinching factor is defined such that with E 2 ν β being the second energy moment. Figure 1 shows the temporal evolution of the flux (Eq. 1) of neutrinos and antineutrinos of all flavors for our BH forming collapse occurring at 10 kpc from the observer. After the neutronization burst, we can see the long-lasting accretion phase characteristic of our BH forming collapse model together with the abrupt termination of the neutrino signal at the moment of BH formation.
While propagating in the stellar envelope, neutrinos undergo flavor mixing [3]. In addition to the well known Mikheyev-Smirnov-Wolfenstein (MSW) effect [60][61][62], relevant at large distances from the core, neutrino-neutrino interaction is important in the innermost region, and especially in the neutrino decoupling region [3,[63][64][65]. In addition, if the matter background should have significant fluctuations, flavor conversion physics would be further affected, as it happens in the presence of turbulence [66,67]. The modeling of the neutrino propagation and its flavor evolution in the dense stellar core is an area of active research due to the conceptual and numerical challenges that it entails, see e.g. Refs. [3,63,64] for dedicated discussions. In order to factor in any uncertainty in the neutrino signal linked to flavor mixing, in this paper, we rely on two extreme mixing scenarios to illustrate the maximum variability of the neutrino signal. The first scenario that we consider is the absence of neutrino mixing (referred to as "no-mix"), such that the neutrino signal is not affected by flavor mixing. The second extreme scenario is the one of full flavor mixing (referred to as "full-mix"), such that the flux of ν e and the one of ν x are swapped under the assumption of full flavor conversion; and a similar assumption holds forν e andν x .

III. NEUTRINO TELESCOPES
In order to explore the pointing precision for a BH forming collapse by relying on the final tail of the neutrino signal or on the bounce time, we consider three of the largest existing and upcoming neutrino telescopes: Hyper-Kamiokande (HK), the Deep Underground Neutrino Experiment (DUNE), and the IceCube Neutrino Observatory (IC). In addition to their size, the geographic distribution of these detectors provides a power-

A. Hyper-Kamiokande
Hyper-Kamiokande [46] is an upcoming water Cherenkov neutrino telescope, which will be located in Japan and will be the successor of Super-Kamiokande, well known for its sensitivity to supernova neutrinos due to its high photocoverage [45]. Following Ref. [41], we assume that the total mass of Hyper-Kamiokande will be 374 kt with detection efficiency HK = 100% [68]. The angular coordinates of the geographical location of Hyper-Kamiokande are reported in Table I.
The main channel of neutrino detection is inverse beta decay (IBD):ν e +p → e + +n. Other sub-dominant channels are neutrino-electron elastic scattering and neutrinooxygen charged current interactions, but these are not considered in this work. The neutrino event rate is computed by folding the neutrino differential flux defined in Eq. 1 with the IBD differential cross section σ IBD [70]. The kinematic threshold for IBD is E thr = m e + + ∆m np = 1.804 MeV, while Hyper-Kamiokande is designed to have a threshold of ∼ 3 MeV for the positron kinetic energy. However, because of the typical neutrino energy from BH forming collapses and the fact that the background is large at small energies, we use E e,thr = 7 MeV [47] as threshold for the positron kinetic energy. This corresponds to a threshold on the neutrino energy of 8.293 MeV. The resultant event rate for our BH forming collapse is where N p is the number of protons in Hyper-Kamiokande and E e is the positron energy. The neutrino event rate expected in Hyper-Kamiokande for our BH forming collapse benchmark model is shown in the top panel of Fig. 2 as a function of time. The "no-mix" case shows higher event numbers than the "full-mix" one due to the higher flux ofν e compared to ν x as seen in Fig. 1. This is a feature of the accretion phase, although the exact relative evolution as a function of time for the two mixing scenarios depends on the mass of the supernova model, see e.g. Refs. [54,[71][72][73][74].
Concerning the background, we base our estimate on Ref. [41] which quote 0.1 s −1 events per 20 kt. For Hyper-Kamiokande, this corresponds to a total background rate of 1.87 × 10 −3 ms −1 . The time resolution in Hyper-Kamiokande is O(ns) [46], which is several orders of magnitude below the requirements for triangulation.

B. DUNE
DUNE will be built at Sanford, South Dakota in the US (see Table I for its latitude and longitude). DUNE will be a liquid argon time projection chamber where scintillation light as well as drifted ionization signals will be used to precisely reconstruct timing and topology of the recorded events. Its fiducial mass is expected to be 40 kt. The main detection channel is the charged current interaction between electron neutrinos and argon: ν e + 40 Ar → e − + 40 K , whose cross section, σ Ar , is provided in Refs. [75,76]. The sub-leading neutral current channel, ν + 40 Ar → ν + 40 Ar , is not included due to the uncertainties in its cross section. Other sub-leading channels that are not considered in this work areν e charged current interactions and elastic neutrino-electron scattering.
The main background in DUNE is expected to be coming from fast neutrons at a rate of 4 s −1 for an electron energy threshold of 6 MeV, assuming no shielding and filtering [77]. In addition, solar neutrinos are expected to contribute with a rate of 110 per day with a 5 MeV threshold [78], and the spallation background for the same threshold is 4 × 10 −4 s −1 [79]. Hence, in total, the background rate is assumed to be 4.0017×10 −3 ms −1 .
The main charged current channel proceeds through a The event rate is plotted for Hyper-Kamiokande, DUNE and IceCube from top to bottom, respectively and for the two extreme mixing scenarios considered in this work. The detector background has been plotted for each detector, but it is only visible for IC in the form of an offset of the event rate. In each panel, the inset plots highlight the rise of the neutrino curve.
number of different excited states, where the most important has a threshold neutrino energy of 5.888 MeV [77]. Hence, in order to control the background, we adopt E thr = 12 MeV as neutrino energy threshold. While DUNE = 86% is used for the detection efficiency [80].
The resultant event rate for our BH forming collapse, R DUNE , is defined by substituting in Eq. (4) the number of argon targets as well as the efficiency of DUNE and the ν e + 40 Ar cross section and it is shown in the middle panel of Fig. 2 as a function of time. The time resolution for the time projection chamber is 0.6 ms, but this can be improved to sub µs resolution by using the scintillation light detected by the photon detector system [75,81].

C. IceCube Neutrino Observatory
The IceCube Neutrino Observatory consists of a cubic kilometer of instrumented ice at the South Pole; Table I reports its angular coordinates. This neutrino telescope is primarily intended for the detection of TeV-PeV neutrinos and, hence, is more sparsely instrumented than Hyper-Kamiokande and DUNE. Neutrinos in IceCube are primarily detected via IBD through the Cherenkov radiation emitted as the produced positron transverses the ice. Due to its sparse instrumentation, IceCube cannot detect individual neutrino events at MeV energies, instead the neutrino burst from a BH forming collapse could be detected as a general increase in the background noise rate [82].
Since there is no reconstructed energy, the IBD kinematic threshold of 1.804 MeV is used as the threshold for detection. The effective mass of IceCube is ∼ 3.5 Mton [82], and the time resolution is a few ns [83]. The IceCube event rate, R IC , is implemented by following the procedure outlined in Sec. V of Ref. [74]. While, by taking into account the detector dead time, the overall background rate is expected to be 1.48 × 10 3 ms −1 [82]. The resultant event rate for our BH forming collapse is shown in the bottom panel of Fig. 2 as a function of time.

IV. TIMING THE NEUTRINO SIGNAL
The angular location of the collapse of a massive star could be determined through triangulation, if the neutrino signal happens to be recorded by several neutrino detectors in different geographical locations [34,[37][38][39][40][41]. In fact, the time delay in the signal detected in different neutrino telescopes allows to determine the neutrino arrival direction. We wish to compare the neutrino triangulation outcome by employing the rise time 2 (i.e., the approach most commonly adopted in the literature) and the cutoff of the neutrino curve at the moment of BH formation. To this end, in this section we evaluate the accuracy and precision of timing by taking into account the neutrino event rate expected in the neutrino telescopes introduced in Sec. III. For the sake of simplicity, in this section, we assume that all detectors are at the same location, i.e. there is no delay between the detectors due to the neutrino time of flight. This calculation will be useful when correcting for the differences between the detectors (see Sec. V).
The expected event rate is the sum of the background rate and the event rate due to neutrinos from the BH forming collapse, as shown in Fig. 2. In order to account for random coincidences, the pure background rate is added before and after the neutrino signal from the BH forming collapse. The combined event rate is binned in 0.01 ms bins, and 10 4 realizations of the neutrino event rate are drawn for each detector by assuming Poisson statistics. For each realization, the rise time and the end time of the neutrino signal are computed together with the time differences. From the set of realizations, the mean and uncertainty on each number are determined as detailed in this section.

A. Determination of the rise time of the neutrino curve
The rise times of the neutrino event rates in Hyper-Kamiokande and DUNE are determined by relying on the first non-isolated detected event; i.e., the first event which is followed by another one in the time span dt. This method is similar to the one used in Ref. [41], where dt = 15 ms was adopted. For our BH forming collapse at 10 kpc, we assume dt = 1 ms for Hyper-Kamiokande and dt = 0.66 ms for DUNE; this choice minimizes the number of random coincidences and the effect of prebounce neutrinos. Even with our choice of dt, the uncertainty in the determination of the rise time of the neutrino curve tends to be dominated by outliers, hence the first non-isolated events before t reject = −10 ms for Hyper-Kamiokande and t reject = −20 ms for DUNE are rejected.
The approach described above is not applicable to Ice-Cube due to its large background rate. Instead, the following fit for the event rate is adopted [84]: where t 0 is the bounce time that we aim to determine and τ is a parameter fitted to the signal for each realization. This fitting function is applied between −50 ms (i.e. before the bounce, when the background dominates) and the first 100 ms post bounce. Then, the parameters R IC , τ and t 0 are fitted to the simulated signal for each event rate realization, as illustrated in the left panel of Fig. 3. The mean and standard deviation of the rise time obtained over 10 4 realizations are listed in the top of Table II for Hyper-Kamiokande, DUNE and IceCube, for both mixing scenarios considered in this paper. One can see that the rise time strongly depends on the mixing scenario and the error on its determination is larger  Fig. 2). This procedure allows to achieve an optimal determination of the time delays. Notice that the points in the left panel are offset on the y-axis by random numbers between -0.5 and 0.5 for readability.
in the "no-mix" scenario, because of the slower rise time of the neutrino curve, as shown in Fig. 2. The estimated rise time is lower in the "full-mix" scenario for Hyper-Kamiokande and IceCube, since theν x signal starts slightly earlier thanν e and rises faster. The relatively large uncertainty in DUNE is mainly due to the pre-bounce neutrinos that produce a tail of estimated bounce times down to t reject . The peak due to the neutronization burst is not observable in DUNE in the "fullmix" scenario and, hence, the estimated mean rise time is larger. Since the pre-bounce neutrinos are not visible in the "full-mix" scenario in DUNE, the uncertainty is improved. This behavior could be avoided by tailoring the timing method to the mixing scenario-by using a smaller value of dt in the "no-mix" scenario; however, this approach is not optimal given the current uncertainties on flavor mixing.
The time difference for pairs of detectors is calculated by estimating the rise time from one realization in each detector and then by computing the difference between the bounce times measured in each of the two detectors. The figures reported in Table II are the mean and standard deviation over 10 4 realizations. One can see that the best determination of the time difference can be obtained by combining IceCube and Hyper-Kamiokande because of their very large statistics.   For Hyper-Kamiokande and DUNE, the determination of the cutoff time of the neutrino curve follows a proce-dure similar to the one adopted for the rise time. The last non-isolated neutrino event is determined by selecting the last event, which is less than dt after the second to last event. For Hyper-Kamiokande, dt = 0.5 ms is used since a smaller value might discard events in the tail of the signal. For DUNE, dt = 1.5 ms is adopted to maintain a low number of discarded events that are part of the signal as well as a low number of outliers. If the last non-isolated event occurs at times larger than t reject = 10 ms after the end of the signal, it is counted as an outlier and rejected.
For IceCube, the end of the event rate is fitted through two plateaus connected by a linear function, as shown in the right panel of Fig. 3. The first plateau, the tail duration (expected to be around 0.5 ms), as well as the end of the tail are determined by the fit. The second plateau is taken to be the background rate determined after the tail. Note that if the actual position of the cutoff is used as the initial guess for the fit, the results are artificially pushed to this particular value. To avoid this, we draw the initial value from a Gaussian of width 0.05 ms centered on the actual cutoff.
The average cutoff times computed through this method are listed in the bottom part of Table II for each of the considered neutrino telescopes as well as for their pairs. Note that the error in the determination of the cutoff time of the signal is worse for the "full-mix" scenario with respect to the "no-mix" case, as expected, given the smaller event rate (see Fig. 2). Notably, since the number of expected events is larger at the end of the neutrino curve (due to the sustained accretion phase) than for its initial rise, the determination of the end time is less sensitive to the flavor mixing scenario and it is on average more accurate.

V. TRIANGULATION
The angular location of a BH forming collapse with respect to two neutrino telescopes can be constrained by measuring the time delay of the signal measured in the two neutrino detectors. Detectors with a larger separation are expected to be less affected by the uncertainty in the determination of the delay times. Hence, the best chances for locating a BH forming collapse with high accuracy and precision are expected for detectors homogeneously distributed across the Earth. In this section, we present our findings for the determination of the BH forming collapse angular location and its uncertainty by employing the rise time as well as the cutoff time of the neutrino curve and by taking into account the geographical displacement between the neutrino detectors as from Table I.

A. Pinpointing the angular location
The geographical location of the detectors adopted in this paper is reported in Table I; hence, the distance between each pair of detectors is: These large distances demonstrate that all three detectors are well separated, and the maximal time of flight for a relativistic neutrino traveling between two detectors is about O(30) ms for all three detector combinations.
Locations on the sky can be defined in the equatorial coordinate system 3 , which is fixed against background stars. We assume that the BH forming collapse is observed at noon UTC on vernal equinox. For a BH forming collapse at right ascension α and declination δ, the travel direction of the neutrinos can be calculated as n = −(cos α cos δ, sin α cos δ, sin δ) T .
The time delay between the detectors i and j is where r i is the position of the detector i in Cartesian coordinates, and c is the speed of light. The measure of the neutrino time delays is not only determined by the location of the BH forming collapse with respect to the neutrino detectors, but also by the method employed for the timing. In fact, the different detector materials, sizes, and technologies lead to time estimates that deviate beyond statistical fluctuations-this uncertainty is enclosed in a bias term, and the measured time delay is defined as [41]: The bias B i,j between the detectors i and j can be estimated as the mean time differences in Table II, where t true i,j = 0. The measured time delays defined as in Eq. (9) and assuming a BH forming collapse at the Galactic Center detected at vernal equinox at noon are listed in Table III. The trend is very similar to what was found in Table II, and the physical explanations are the same. By subtracting the biases, the expected time delays are recovered with good precision (results not shown here). The time delay between a pair of detectors allows to determine the angle θ spanned by the two lines connecting the detectors and the location of the BH forming collapse in the sky. This angle defines a cone in the sky. Using three detectors, three cones would intersect in two points, while a unique location can be determined with four detectors. In this paper, we focus on combining three detectors together in order to show the potential of the method employing the cutoff time of the neutrino curve.

B. Determination of the angular uncertainty
The uncertainty on the angle θ for each pair of detectors can be determined without full triangulation [34,40]. In fact, for a given pair of detectors, where D i,j is the distance between the two detectors (see Eqs. 6). The resulting uncertainty on the angle is [34,40]: which allows to determine the maximal and minimal uncertainties on the angular location of the BH forming collapse.
In order to derive the uncertainty on t true i,j , the bias uncertainty and the one related to t measured i,j are added in quadrature (see Eq. (9)): This approach leads to an increase of the uncertainty on t true i,j with respect to t measured i,j by approximately √ 2. Although the employment of the true bias might provide an optimistic pointing accuracy, our main focus is on pointing precision and, for that, Eq. (12) is quite conservative [41]. The minimal and maximal uncertainties, δ(θ i,j ), computed through this approach are reported in Table IV for the method employing the rise time of the neutrino curve and the one relying on the cutoff time. By using the cutoff time rather than the rise time, one can see that the minimal uncertainty is significantly below 1 deg.

C. Combining multiple detectors
When combining the time differences between n different detectors, correlations between the detector combinations must be accounted for. This is done by including the covariance matrix V . In analogy with Eq. (12), we account for the bias correction by using where each matrix is of dimension m × m, with m given by the binomial coefficient m = n 2 . The entries in the covariance matrices correspond to the time delay vector t = (t 1,2 , t 1,3 , . . . , t 1,n , t 2,3 , . . . , t n−1,n ) T .
In order to determine the preferred regions on the sky, we rely on the chi-square: where t true is computed using Eq. (9) with the values from Tables III and II and t(α, δ) is the time delay vector computed by assuming that the BH forming collapse has right ascension α and declination δ. For two detectors, there is only one entry in t, and the matrix product reduces to a simple scalar product. Assuming no correlations, Eq. (15) reduces to a sum over such scalar product [37,39,41]. However, when more detectors are included, the time delays are correlated by definition, and using a diagonal covariance matrix underestimates the uncertainty in determining the direction.
The sky regions constrained by combining two and three detectors using the rise times and cutoff times are shown in the left and right panels of Fig. 4, respectively. 14h 16h 18h 20h 22h 0h 2h 4h 6h 8h 10h The colored bands represent the regions of the sky constrained by using the indicated two-detector combinations, whereas the black contours represent the combinations of three detectors. As expected, the pointing precision of the cutoff times is superior to that of the rise times.

VI. DEPENDENCE OF THE POINTING PRECISION ON THE DISTANCE OF THE BLACK HOLE FORMING COLLAPSE FROM EARTH
The flux of neutrinos from a BH forming collapse scales with the inverse distance squared, hence affecting the determination of the angular location of a BH forming collapse. In this section, we investigate the dependence of the timing and pointing precision from the distance at which the BH forming collapse is located with respect to Earth.
A. Dependence of the pointing precision on the distance: rise-time method The results presented so far are based on the assumption of a BH forming collapse at d = 10 kpc from Earth. However, as shown in Eq. (1), the neutrino flux scales as d −2 . Consequently, we should expect that the typical delay between the first event and a following event increases with d 2 . When determining the rise time in Hyper-Kamiokande and DUNE, the time span dt and the rejection time t reject introduced in Sec. IV A should scale as functions of d. For Hyper-Kamiokande, the event rate is high enough that we do not need to increase t reject , and we keep it constant at −10 ms. Whereas, for DUNE, we use a rejection time t reject = − (d/10 kpc) 2 × 20 ms for distances above 10 kpc and t reject = −20 ms for smaller distances. The time span dt is simply scaled with d 2 to match its value at 10 kpc.
The rise time in IceCube can be determined by adapting the fitting function introduced in Eq. (5) to the signal (see also Fig. 3). This approach continues to give a good estimate of the timing up to 60 kpc. Above this distance, the uncertainty becomes so large that no directional information can be extracted through triangulation.
The minimal and maximal angular uncertainties as a function of the distance of the BH forming collapse derived from the rise time are shown in the left panel of Fig. 5. The uncertainty increases gradually until a distance of 25 kpc above which the timing in DUNE degrades rapidly. Notice that, although the bands for IC-DUNE and HK-DUNE appear to narrow in the proximity of 25 kpc, the vertical width does not shrink significantly while the horizontal one scales as a function of the distance. On the other hand, the combination of IceCube and Hyper-Kamiokande displays a gradual increase in uncertainty until 60 kpc where the minimum uncertainty on the angle is ∼ 50 • corresponding to a very large fraction of the sky.
B. Dependence of the pointing precision on the distance: cutoff-time method The determination of the cutoff time for a BH forming collapse at 10 kpc has been introduced in Sec. IV B. As d varies, we use dt = 0.5 ms for Hyper-Kamiokande and dt = 1.5 ms for DUNE. As for t reject , we use t reject = 10 ms for Hyper-Kamiokande independently of d, since the event rate is so large that it allows for an excellent discrimination of the signal from the background even for large d. On the other hand, for DUNE, we use t reject = (d/10 kpc) 2 × 10 ms and fix t reject = 125 ms for d providing an otherwise larger t reject .
As for the cutoff fit for IceCube (see Fig. 3), we cannot assume that we a priori know the time of the cutoff within ±0.05 ms for any d, hence we draw the initial guess from  11)) on the determination of the BH forming collapse location in the "no-mix" scenario as a function of the distance for the rise time method (left panel) and the cutoff time method (right panel). The minimum and maximum angular uncertainties are calculated using Eq. (11), and the uncertainty is corrected for the effect of the bias by multiplying with √ 2 (see Eq. 12). The pointing accuracy is significantly better when using the cutoff time method (right). In addition, due to their large statistics, the combination of IceCube and Hyper-Kamiokande allows to achieve a better performance. a wider distribution. An appropriate choice for the width is σ draw = 25R IC,bckgd /R 2 IC , where R IC,bckgd and R IC are the background and signal rates respectively. Due to the linear tail, we require σ draw > 0.05 ms. At d = 10 kpc, a period with pure background of 10 ms is added after the signal. The period with background is increased to 3 σ draw when 3 σ draw > 10 ms.
The minimal and maximal angular uncertainties from the cutoff time method are shown in the right panel of Fig. 5. At distances up to ∼ 15 kpc, the maximal uncertainty is below 10 • . At larger distances the precision decreases notably and, as for the rise time method, the timing in DUNE degrades fast. The combination of IceCube and Hyper-Kamiokande provides the highest precision at all distances. For BH forming collapses throughout the Milky Way (up to 23 kpc) the precision is well below 10 • . At larger distances, the noise in IceCube leads to larger uncertainties, and at ∼ 60 kpc the precision is comparable to the one obtained with the rise time method (see left panel of Fig. 5).
The uncertainty achieved by combining all three detectors cannot be directly extrapolated from the bands for two detector combinations in Fig. 5, but the general trends are similar for two reasons. First, the width of the combined region will never be larger than the smallest uncertainty of the two detector combinations. Second, in the region where DUNE has no sensitivity (rise time beyond ∼ 30 kpc) the only constraint will come from the combination of IceCube and Hyper-Kamiokande, and hence this gives a very accurate picture of the full three detector uncertainty. Consequently the width of the bands in Fig. 5 are a good proxy for the uncertainty in the three detector case, although precise numbers require additional calculations that go beyond the scope of this work.

VII. DISCUSSION
The employment of the cutoff time method is clearly promising for the determination of the angular location of a naked stellar-mass BH, in agreement with the findings of Refs. [37,39,40]. One of the reasons is the high event rate due to the long-lasting accretion phase, preceding the instant of BH formation. This, combined with the steepness of the tail of the neutrino curve, allows for a better discrimination of the signal with respect to the background than in the case of the rise time method. In the following, we discuss some of the unknowns potentially affecting our analysis.
• Astrophysical uncertainties. If the newly-formed BH is not naked, a neutrino echo could affect the tail of the neutrino curve because of the nonnegligible supersonic accretion flow surrounding the central compact object [31]. Hence, one should expect a smearing of the tail of the neutrino curve and a less precise determination of the BH forming collapse location. Still, we expect that the cutoff time method will perform better than the rise time method due to the sharper transition between the time bins where the event rate is dominated by the signal and the time bins dominated by the background.
The modeling of the tail of the neutrino curve has been inspired by Ref. [28] for our benchmark BH forming collapse. Of course, one should keep in mind eventual variations on the steepness of the tail due to current uncertainties in the theoretical modelling of BH forming collapses and the lack of multidimensional hydrodynamical simulations tracking the instant of BH formation.
In addition to the linear tail, an exponential component with a decay time of O(0.05) ms [28,56,57], due to the neutrino trajectories being bent by the gravitational field, should be expected. This time scale is comparable to the precision that we obtain with Hyper-Kamiokande and IceCube, and should not increase the overall uncertainty by more than a factor of two. Hence, we expect that the impact of this correction on our results should be minor.
• Neutrino mixing. Because of the current uncertainties involved in the modeling of the neutrino flavor evolution in neutrino-dense media [3,63,64], we have chosen two extreme scenarios: absence of flavor mixing or full flavor conversion. Importantly, the rise time method strongly depends on the mixing scenario because of the different rise of the electron and non-electron neutrino curves (Fig. 2); whereas, the cutoff method is less sensitive to it, see Table IV, due to the anyway large event rate. In this sense, despite the astrophysical uncertainties, the cutoff time method is more promising than the rise time method.
• Proposed methodology. Our results are in agreement with the findings of Refs. [37,39]. However, despite comparable theoretical input for the neutrino curve with respect to Ref. [37], our statistical method allows to achieve an accuracy below ms even for the rise time method (see Tables II  and IV), while Ref. [37] assumed a few ms precision. Our work differs from Ref. [39] since we focus on next generation neutrino detectors and optimize the statistical technique as well as the detection strategy.
• Bias. A word of caution should be added concerning the bias often neglected in the existing literature, see e.g. Refs. [37,39]. The true bias of the time differences is estimated by the mean values in Table II. For the rise time method, the bias strongly depends on the mixing scenario. As a consequence, it is difficult to provide reliable theoretical input for determining the bias. On the contrary, the bias for the cutoff time method is almost independent of mixing since the signal in all flavors is cut off simultaneously.
A data driven bias could be determined as well [41]. This does not rely on any model assumptions, but also cannot account for the full bias. Since for the cutoff time method, the bias is largely independent 14h 16h 18h 20h 22h  6. Impact of the bias on the determination of the location of the BH forming collapse in the "no-mix" scenario. The true location of the BH forming collapse is marked with a black dot. The 1, 2, and 3 σ confidence level contours are plotted in gray for the cases without bias correction, in red for the rise time method, and in blue for the cutoff time method; the latter is magnified in the plot inset. Correcting for the the bias affects the pointing accuracy and it should not be neglected. The red and blue 1σ contours correspond to the black lines in Fig. 4.
on the mixing scenario, the true bias determined in Table II is more robust to our purpose than a data driven one.
In order to highlight the impact of the bias on the determination of the angular location of a BH forming collapse, Fig. 6 shows the uncertainty contours for the rise and cutoff time methods for the cases with bias correction as in Fig. 4 (red and blue contours) and the ones without bias correction (gray contours), i.e. when t measured i,j is used to estimate the BH location in Eq. (15). One can clearly see that, without including the bias, the best fit is shifted with respect to the true location of the stellar collapse and the uncertainty bands are underestimated.
• Experimental uncertainties. Despite the fact that we have developed an improved strategy to determine the cutoff time of the neutrino curve, our method remains sensitive to our knowledge of the detector backgrounds. For example, the possibility of further reducing the backgrounds (e.g. through gadolinium enrichment in Hyper-Kamiokande) as well as a better knowledge of the backgrounds (e.g. in DUNE) could surely improve our findings concerning the determination of the angular accuracy on the location of the BH forming collapse.
In addition, the employment of the neutral current channel in DUNE could boost the event rate, allowing a better timing performance, especially as the distance of the BH forming collapse increases.

VIII. CONCLUSIONS
The collapse of a massive star could lead to BH formation. In this case, neutrinos could carry precious input to locate the BH forming collapse in the sky. In this work, we explored the possibility of improving on the determination of the angular location of a BH forming collapse through triangulation pointing, under the assumption that accretion in the surrounding of the central compact object is halted and naked BH formation occurs.
By forecasting the neutrino event rate from a BH forming collapse expected in the upcoming Hyper-Kamiokande and DUNE as well as the existing IceCube Neutrino Observatory, we show that focusing on the tail of the neutrino curve rather than on its bounce time, as often considered in the literature, allows to achieve at least an order of magnitude better determination of the angular location of the stellar collapse. This is due to the larger event rate expected during the accretion phase preceeding BH formation, which is more easily distinguishable from the detector backgrounds. Notably, despite the existing uncertainties in the modeling of neutrino mixing in the stellar core, we find that the employment of the cutoff time of the neutrino signal for triangulation pointing is less sensitive to flavor mixing.
This work also highlights the importance of the bias in the determination of the angular location of the stellar collapse. In fact, neglecting it could lead to a wrong esti-mation of the true location of the BH forming collapse as well as to underestimate the uncertainty bands. Given the similar behavior for all flavors at cutoff, the bias correction is also less sensitive to the mixing scenario when the cutoff time method is employed.
Our improved strategy for triangulation pointing of BH forming collapses contributes to the growing set of real-time strategies to improve the pointing capabilities through multi-messenger data. However, work remains to be done to pin down the theoretical uncertainties in the modeling of the neutrino signal at BH formation in multi-dimensions as well as to improve our understanding of the detector responses and backgrounds.