Black hole interference patterns in flavour oscillations

Motivated by neutrino astronomy, we consider a plane wave of coupled and massive flavours, scattered by a static black hole, and describe analytically and numerically the corresponding oscillation probability in the surrounding space. Both the interpretation as particles travelling along geodesics and as scattered waves are studied, and consistently show a non-trivial and potentially long range interference pattern, in contrast to the spatially uniform transition probability in a flat spacetime. We introduce a numerical method for studying the oscillations around black holes, which accounts for the full curved geometry and flavour wave mixing. Whilst limited to the region immediately around the black hole, this numerical approach has the potential to be used in more general contexts, revealing the complex interference patterns which defy analytic methods.


Introduction
The emerging field of neutrino astronomy at observatories such as the Ice Cube [1], Antares [2], KM3Net 2.0 [3] and the Baikal Neutrino Telescope [4] promises to improve our understanding of the fundamental properties of neutrino physics. If these are to be combined with multi messenger observations from gravitational wave events at LIGO/Virgo ( [5] and see [6] for a review), an understanding of the behaviour of neutrinos in curved spacetimes, in particular their flavour oscillations, will be required for proper interpretation of the results. Neutrino interferometry has, for example, been proposed as a means by which to discover asymmetry in neutrinos and antineutrinos, and thus distinguish Dirac and Majorana neutrinos [7,8].
The phenomenon of flavour oscillations is described by Quantum Mechanics and, as explained in more detail below, is based on the study of the wave function of neutrinos, more specifically of their phase. The phase shift induced by gravitational fields on particles' wave functions was originally observed with neutrons in the Earth's gravitational field [9], with the corresponding effect first studied theoretically in [10]. For neutrinos oscillations, the effect of gravitational fields on the wave function along geodesics should also have phenomenological implications [11], due for example to the gravitational redshift which increases the oscillation length [12].
General studies of gravitational effects on the phase of wave functions, for different spins, were done in [13], [14], [15], [16], [17], and more specific studies involve supernovae explosions [18], neutrino optics and the Lens-Thirring effect [19], or the study of three flavour oscillations in curved space time [20]. There is also a potential violation of the Equivalence Principle in relation to flavour oscillations in curved space time, as discussed in [21], [22], [23].
We study here the effect of a curved background on the flavour transition probability, motivated by the scenario where neutrinos emitted by an astrophysical event pass by a BH before being observed. Such a scenario has been studied in [24], where two neutrino beams are lensed by a massive object, and the oscillation probability is calculated at the focal point. We propose instead to assume an incident plane wave of two mixed flavours, which is scattered by a BH. In this stationary process, we calculate analytically and numerically the phase shift which is gravitationally induced at every point of space, and which allows us to determine the oscillation interference pattern. This study therefore takes into account delocalised effects arising from a wave description, and not only the interference of two specific geodesics.
Section 2 reviews features related to oscillation probability, in flat spacetime and also in the presence of a BH. In the latter case, the description is based on particles travelling along geodesics, and we approximate an incident plane wave as a set of parallel and coherent beams, with different impact factors. In this cylindrically symmetric configuration, the resulting scattering pattern exhibits paraboloids of minimum probability, aligned with the axis of symmetry, at large distances. Although highly localised, these paraboloids would in principle result in directions far from the BH along which no oscillation of neutrino flavours would be observed.
We then turn to the wave description in Section 3, where we give generic analytical features of a massive wave scattered by a BH. The resulting equations can be solved analytically only in very limited cases, and the aim of this section, together with Appendix A, is mainly to understand the scattering properties of the model. Section 4 covers the numerical implementation of the study, in the vicinity of the BH, following the dynamical evolution of plane wave flavour mass eigenstates on a fixed Schwarzschild background in horizon penetrating coordinates. The resulting interference pattern thus takes into account the self interference of the waves, the full effect of the curved geometry and the resulting backward scattering. Appendix B gives further details of the code used and convergence tests. Whilst we are far from being able to directly observe neutrinos in the vicinity of BHs, such interference patterns may affect the weak interactions around black holes or massive objects such as neutron stars, for example, in a merger scenario.

Flavour oscillations
We consider here two scalar flavours, since the essence of the flavour oscillation process is independent of the spin. Also, for fermions, the spin flip effect induced by gravity actually does not occur in Schwarzschild spacetime [25]. The present study is based on plane waves, although this assumption is associated with several ambiguities in the context of flavour oscillations. As explained in the clear and thorough review [26], these ambiguities can be avoided with wave packets. Nevertheless, the plane wave assumption is enough to construct the essential features of the interference pattern we are interested in.

Oscillation probability in the absence of gravity
Consider two scalar flavours φ a , φ b , which satisfy the equation of motion ( The eigenmasses and the corresponding mass eigenstates are given by where α is the mixing angle and

The equation of motion satisfied by the mass eigenstates is
with equal-momentum plane wave solutions in flat space, given by Rotating back to the flavour eigenstates leads to flavour configurations To calculate the oscillation probability, one introduces at every point of space the normalised kets for the mass eigenstates, and for the flavour eigenstates. For the field configurations (7), the phases are ϕ ± = k · r − ω ± t and it is easy to check that the flavour oscillation probability has the known expression and is independent of the spatial coordinates. If one considers a beam of ultra-relativistic particles, the propagation time can be identified with the propagation length, in a system of units in which c = 1. The oscillation length is then and corresponds to the experimental observable.

Gravity-induced phase shift
In the presence of gravity, the phase of mass eigenstates depends on the spatial coordinates in a non-trivial way. We assume then that the mass eigenstates have the same amplitude, which is a valid approximation if one considers the situation where the masses are almost-degenerate [26]. For the stationary scattering process we consider, the time dependence of mass eigenstates is trivial, and a configuration can be written where r denotes the set of space coordinates. Assuming A + A − , the amplitude can be factorised in the linear combinations representing the flavour eigenstates. This amplitude therefore cancels out in the normalisation of states for the calculation of oscillation probability, which thus depends on the phases Φ ± − ω ± t only. As a consequence, the same argument as the one leading to eq.(10) gives here P a→b = sin 2 (2α) sin 2 where the phase shift we are interested in is and is the main subject of our numerical analysis presented in section 4. Time cannot be replaced by the distance travelled in the probability (13), since one doesn't consider specific geodesics, but one can nevertheless characterise the motion of surfaces of constant probability (13). For this, at a given time t one denotes by D a surface corresponding to a given probability, far enough from the BH for spacetime to be assumed flat. As time increases ∆ must compensate the change in (ω + − ω − )t, for the probability to stay constant on D, which must therefore be shifted and deformed in such a way that where v is the velocity of points on D. The gradient ∇∆ is perpendicular to D at each point, and its scalar product with the velocity of D at this point is therefore a constant characterising the system. For a stationary process as we consider here, the surfaces D are periodically generated and move away from the scattering centre.

Neutrino interferometry in Schwarzschild geometry
Before turning to the study of a massive wave scattered by a BH, we describe here the point of view of particles traveling along specific trajectories. The first calculation on neutrino interferometry was done in [24], and related comments can be found in [27], [28], [29] and [30]. We consider the motion of a test particle of mass m in the Schwarzschild metric with where R s = 2GM and we restrict the study to the equatorial plane θ = π/2, where the conserved quantities are p t = mf (r) dt ds and p ϕ = −mr 2 dϕ ds .
Consider a neutrino source at the position r A where the beam is split into two beams, which are later lensed by a BH, to meet again at the position r B . The flavour oscillation probability at the intersection point has been calculated in [24], after evaluating the phase along null geodesics. Assuming that the impact parameters b 1 , b 2 of the two intersecting beams are large compared to R s , and up to terms of order (m ± /k) 2 , the result is where R ≡ r A + r B and Motivated by a plane wave scattered by a BH, we consider here a source at minus infinity in the xdirection, which produces asymptotically parallel and coherent beams with different impact parameters −∞ < b < +∞. Far enough from the BH, we can use the asymptotically Cartesian coordinates (x, y) in the equatorial plane, centred on the BH, with the impact parameters measured along the y axis, as shown in Fig. 1.
The limit r A → ∞ in the probability (19) can be taken in the following way which leads to the oscillation probability at the intersection of two deflected beams with impact parameters b 1 , b 2 , given by where r is the radial coordinate of this intersection point. It is interesting to consider the lines of constant probability in the equatorial plane, assuming that |b i | >> R s . The deflection angle δ for massless particles with impact parameter b along the y axis is The equation for the asymptotically straight trajectory, after deflection, is then such that the two impact parameters b 1 , b 2 corresponding to beams intersecting at the point (x > 0, y) are The probability in eq. (22) can then be sketched, as in Fig.2, and shows a non-trivial interference pattern for x > 0. Note that only forward scattering is taken into account. As can be seen on Fig.2(a) (and Fig.2(c) for smaller initial momentum), the probability is described by a spherical pattern, with a radial oscillation length consistent with the expression (11), and which arises from the sine squared in the first line of eq. (22). However, one can see in the magnified and stretched views (b) and (d) that there is another oscillation in the orthoradial direction where the probability vanishes along lines which originate from the BH, and extend in almost radial directions. These lines are described by the vanishing of the product of cosines in the second line of eq. (22), which happens when where m m + m − . We focus on the lines which are close to the axis of symmetry, therefore where x >> y. Together with the expressions (25), the condition (26) gives the parabolas where A n ≡ a n 1 + 1 with We note that, since the present derivations are valid close to the axis of symmetry, we have ∆m 2 ∆b 2 << 16kr and thus the 2nd line of eq. (22) is negligible. The parabolas (27) corresponding to different integers n can be seen on the magnified views in Figs.2(b) and (d). Given the cylindrical symmetry of the problem, the BH therefore generates a family of paraboloids of vanishing oscillation probability, along the direction of the original momentum. This additional oscillation length in the orthoradial direction is much smaller than the radial one, since it depends on the masses and not the mass difference.
As can be seen in the non magnified views in Fig. 2(a) and (c), there is a third oscillation feature, again in the orthoradial direction, where the regions of maximum probability are "shifted" out of phase. These out of phase maxima arise from the product of sines squared in the 2nd line of eq. (22), which is neglected to derive the parabolas (27). This term becomes of the same order as the 1st line of eq.(22) in these sectors. The typical oscillation length of these features depends on the mass difference and thus is of the same order as the radial oscillation length.
As a consequence of these new oscillation features, the flavour detected depends not only on the radial distance to the BH, but also on the distance between the observer and the axis of symmetry of the problem.
Finally, note that several independent length scales appear in the problem, and different combinations could lead to similar interference patterns to those shown in Fig.2. These length scales arise from: 1. the eigenmasses m + m − ; 2. the difference ∆m 2 ; 3. the initial momentum k;

the BH mass M .
We therefore have 3 independent dimensionless parameters on which the interference pattern depends, and these are chosen here, for phenomenological purposes, consistently with ultra-relativistic particles with a small mass difference (compared to their masses), and such that the radial oscillation length is of the order of R s .

Massive wave scattered by a static black hole
In this section we describe some analytical features of a massive plane wave scattered by a BH. Studies of massless scalar wave scattering from a BH can be found in [31], [32], [33], [34], [35], and a pedagogical review is given in [36]. As in the situation of light scattering, both forward and backwards glory effects [37] do occur for spin zero particles [38], which is confirmed in this section. Unlike in the previous section, where the description was done in the equatorial plane θ = π/2, we follow convention for these studies and consider here a cylindrically symmetric system where the wave vector of the original plane wave is along the z-axis. In coordinates (16), the fields therefore depends on r and θ, but not ϕ.

Decomposition in partial modes
We consider here a mass eigenstate φ with mass m, satisfying the equation where the metric is given by eq. (16). The so-called tortoise coordinate r can then be introduced as with which the metric is conformally equivalent to and involves the Eddington-Finkelstein coordinates t ± r . In terms of the radial coordinate r , the wave equation for a cylindrically symmetric field reads whereφ ≡ rφ does not depend on ϕ. The next step is to decompose the cylindrically symmetric fieldφ on the Legendre polynomials basisφ where the trivial time dependence describes a stationary process. Given the identity the expansion (34) plugged in the evolution equation, eq. (33) gives, for every l, where the effective potential seen by each mode u l is, Although the coordinate r leads to an elegant formulation of the problem, it is not appropriate for the study at large distance r >> R s in the massive case, as explained in Appendix A.

Semiclassical approach
In the semiclassical approximation discussed in detail in [39], the variable l is replaced by a continuous parameter, with a dominant contributionl obtained from the stationary phase method for the calculation of the phase shift of the mode u l .l satisfies δ(l) = θ , where δ(l) is the classical deflection angle and θ is the observation angle, and the differential cross section of the scattering process is then shown to be dσ dΩ semiclassical ω 2 (l + 1/2) sin θ(dδ/dl)l .
If one compares this expression with the classical differential cross section for scattered particles one can make the identification |b|ω l + 1 2 .
The latter equation gives an important insight into the scattering process: the dominant mode ul in the scattering process corresponds to a beam with typical impact factor given by eq.(41).

Vicinity of the horizon
Near the BH horizon, the radial coordinate can be written and one would expect the potential (37) to be negligible compared to ω 2 . One should check though, for fixed r, that large values of l do not invalidate this approximation. Given that 3 √ 3R s /2 is the critical impact parameter for a null geodesic, one can assume that the dominant beams contributing to the scattering pattern in the vicinity of the BH are those corresponding to an impact parameter satisfying |b| ≤ 3 √ 3R s /2. A detailed discussion of critical impact parameters can be found in [40], which includes the case of massive particles. Since we assume ultra-relativist particles though, the null geodesic approximation is enough for the present discussion.
According to the semiclassical argument above, the dominant modes are then characterised bŷ and satisfy f (r)l (l + 1) As a consequence, the solution to eq.(36) is then approximately independent of l in the relevant range (43). If we assume an ingoing wave only (since no signal can escape this region), we obtain u l ∝ exp(−iωr ), and the full solution is then approximately where We note that [41] and [42] discuss the possibility of an outgoing massive wave in this regime, in addition to the ingoing wave, with a different amplitude which takes into account a reflection from the Horizon, related to Hawking radiation. This effect is dominant at energies which are comparable to the Hawking temperature, and since the present study is motivated by high-energy neutrinos, we neglect here the outgoing wave.
Assuming that mass eigenstates have the same amplitude, A + (θ) A − (θ), the oscillation probability is given by P vicinity a→b sin 2 (2α) sin 2 and depends on the Eddington-Finkelstein coordinate t + r only. This is consistent with the numerical analysis in Section 4, where the Kerr-Schild time coordinate t satisfies t +r = t+r , such that the probability (47) for fixed t is almost uniform in the vicinity of the BH, since we have there (ω + − ω − )r << 1 (the oscillation timescale is very long). Note that the expression (47) predicts an infinite number of oscillations, for fixed t, as one approaches the horizon, where r → ∞. However, this is an artifact of the choice of observer at infinity, for whom t diverges at the horizon. In terms of t , the number of oscillations is finite.

Asymptotic solution
Eq. (36) can be solved asymptotically, by analogy with the scattering problem in a Coulomb potential. The different steps are explained in Appendix A, and the solution can be written in terms of the momenta and the radial coordinate ρ ≡ r + pR s 2k ln(2kr) .
If A is a constant amplitude, the solution can then be expressed in terms of the non-scattered plane wave φ plane , and we find where and the complex phase δ l is expressed in terms of the momentum p, as explained in Appendix A. Note that the plane wave appearing in eq.(50) is "distorted" in the sense that it involves both radial coordinates, r and ρ, as can be seen in Appendix A. Eq.(66) in this appendix shows that the expression (50) is valid up to terms of order 1/r 2 , and the unavoidable mixture of radial coordinates r and ρ can be understood as a consequence of the long range of the gravitational interaction. Also, the expression (50) is formal only, since the sum appearing in the definition of h(θ) is actually divergent [36]. But if one considers a "thick beam" with maximum impact parameter B, instead of an infinite incident plane wave, then the semiclassical argument given above allows to set a cut off L for l, defined by L Bω. Finally, the result (50) is valid for each mass eigenstate, for which one can extract the phase as but this is left to next section, where the original equations of motion, eq. (30), are solved numerically for each mass eigenstate.

Numerical calculation of the oscillation probability
In this section we describe the results of an evolution of the flavour fields on a fixed Schwarzschild BH background. The purpose is to study the interference patterns produced in the oscillation probability as the plane waves pass the BH. We expect, and indeed find, that a stationary pattern is built up dynamically, radiating from the BH as the evolution progresses. Full details of the numerical methods and convergence tests are provided in Appendix B.

Initial conditions
The two mass eigenstates of the flavour fields are set up as plane waves in the x-direction, as described by eq.(6) with t = 0, superimposed on a fixed background geometry which describes a Schwarzschild BH in Kerr-Schild coordinates (t , r, θ, ϕ) (see Appendix B).
Numerical values are chosen such that m − m + , in order to have approximately the same amplitude for φ − and φ + . The values in units set by the BH radius R s are m − = 1.0/R s and m + = 0.99/R s . The wavenumber is k = π/R s , such that the corresponding wavelength of the mass eigenstates is λ = 2R s . The resulting oscillation length in flat space corresponding to these values would be of the order of L ∼ 2000R s .
These values are not phenomenologically realistic, but are set by numerical constraints. To model neutrinos around solar mass BHs, one should set m ∼ 10 10 /R s and k ∼ 10 17 /R s as in the previous sections, but this would mean resolving two very disparate timescales, which is computationally intractable using the current method.
There are thus two possible interpretations of the numerical results. Firstly, assuming a solar mass BH, the mass values chosen would correspond to m ∼ 10 −10 eV , so very light neutrino-like particles. Alternatively, for neutrinos themselves, the simulated BH would correspond to primordial BHs of order 10 −10 M , although the momentum k should also be made larger in this case.
Despite these issues, we could expect that some features of the numerical solution are still relevant to the neutrino scattering case. We also show results for a higher momentum case k = 2π/R s and observe the differences.

Evolution and results
We evolve the mass eigenstates from the initial conditions above according to the Klein Gordon equation eq.(1), and calculate, at each point of space and for one given time t the phase difference We thus expect that the system will settle into a stationary regime, with the oscillation probability Since the oscillation timescale is of order 2000R s which is much greater than the simulation time, the term ∆(r, θ) should dominate the probability pattern. We show the real part of one of the complex fields in Fig. 3 for two different initial momenta k. The other fields show a similar behaviour, and we see that the amplitude of the two fields is approximately the same, justifying our use of eq. (13) to calculate the probability. Fig. 3 also shows the spatial patterns of the oscillation probability per eq. (54). The results can be summarised as follows: • The scale in the probability varies by an amount of order 0.001, thus the region is overall very uniform, as expected from the analytic calculations in Section 3; • In the interference pattern that builds up around the BH, there appear to be two main contributions: 1. Dark lines build up away from the BH, resulting in complex but roughly paraboloidal structures of minimal probability after the BH.

Finer, approximately radial lines of maximum probability cross the darker ones.
Whilst it is difficult to draw comparisons between the plots here and in Fig. 2, since the latter uses the particle viewpoint in which time variation is converted to spatial variations along the particle path, one could imagine that the fine structure seen in the magnified plots (on orders much smaller than the oscillation length) is related to the structures seen here numerically. According to eq.(15), as time evolves (on the order of the oscillation period) paraboloids of constant probability move along the x-axis and gradually close up.
• The forward glory in the field (the enhanced region in the fields shown in the upper two plots of Fig. 3), subtends an angle at the BH of 16 o in the case of the lower momentum, and this angle halves when the momentum is doubled. This region corresponds to the band of low oscillation probability, between two (c) Oscillation probability k = π/Rs (d) Oscillation probability k = 2π/Rs of the finer radial lines of maximum probability. There is clearly a correspondence between the spatial profile of the forward glory and the oscillation probability, such that the lines of equal probability also become closer together -the structure becomes more fine -with increased momentum.
• The amplitude of the forward glory is 1.4 times as large for the higher momentum case, whereas the overall variation in the probability is roughly halved.
• Unlike the analytic approximations, the numerical simulation shows the full spatial dependence and includes the backscattering effects of the BH, resulting in a non-trivial pattern of enhanced/suppressed oscillation probability on the left of the BH.

Conclusions
In this work we have modelled the interference patterns in flavour oscillations which are generated as plane waves of neutrinos pass by a BH, using both analytic and numerical techniques. In all cases we find highly non trivial interference patterns which potentially extend to large distances from the BH. In general, these are paraboloidal shaped regions of enhanced or suppressed probability. These regions, whilst localised, could result in unexpected neutrino detection patterns, in the situation where a BH is roughly in the line of sight between the Earth and an astrophysical event such as a supernova explosion. Together with the information obtained from electromagnetic or gravitational signals, the unexpected neutrino detection pattern could then give essential features on either the BH or the astrophysical event at the origin of the neutrino flux.
The interference patterns we describe assume coupled massive scalars. A full analysis, using fermions, would give the same result in the non spinning case. However, in the case of a spinning BH, one would need to account for the spin flip effect in addition to the flavour oscillation phenomenon. As in flat space time and in the presence of a magnetic field [46], the resulting flavour oscillation probability would be further modulated, due to a new length scale introduced in the system by the BH spin.
Another simplification is that we consider plane waves of coupled flavours, whereas a more phenomenologically motivated model would consider wave packets. This would be more complicated analytically, but potentially simpler to study numerically, since it is localised and thus would not require a high resolution to be maintained across the whole grid. However, the plane wave situation is a first step in understanding the whole interference structure, which would partially be reconstructed by a wave packet passing through a BH.
We also assume a small difference in the eigenmasses of 1%, such that we can consider only the phase difference between the mass eigenstates and neglect the amplitude contributions to the probability. We will consider the impact of relaxing this assumption in future work.
Finally, due to computational constraints, the numerical study necessarily focused on the regions close to the BH, with the BH horizon scale corresponding to the initial momenta of the plane waves, rather than the oscillation length as in the analytic work. This is phenomenologically less well motivated, either corresponding to neutrinos around low mass, primordial BHs, or ultra-light neutrino-like particles interacting with solar mass BHs. We hope to develop the numerical methods to study more phenomenologically motivated scenarios in future work. However, the present method of study already provides a means by which to study more complex geometries in the future, in particular spinning BHs, for which analytic derivations are intractable. and A l , B l are constants of integration. In order to recover usual notations [36], the constants of integration A l , B l are traded for a global factor A and the complex-valued phase shift δ l such that where Re(δ l ) = η l .

B Appendix: Numerical Implementation
This appendix summarises the key features of the numerical code used in the simulations, which is based on the publicly available numerical relativity (NR) code GRChombo, itself built on top of the open source Chombo framework [44]. For a more full discussion of GRChombo see [45], or the website www.grchombo.org.
Here we describe the key features of the code as used in this work, in which the matter fields were evolved on a fixed BH background, neglecting any backreaction of the fields on the metric (thus note that the NR capabilities of the code were not utilised).

B.1 Background metric
We use the Kerr-Schild form of the Schwarzschild solution as the background on which the fields evolve, thus neglecting any backreaction of the flavour fields on the metric. This metric is horizon penetrating and thus simulates the full exterior solution, but necessitates excision at the singularity -in the static Schwarzschild case this can be done by simply setting the flavour fields within the horizon (in practice we do this for r < R s /2) to decay to zero. Since the curvature of the metric prevents signals from escaping, errors due to the excision do not (in principle) propagate outwards. The Kerr-Schild time coordinate t is related to the Schwarzschild coordinate t by adding the difference between the tortoise coordinate r and the standard Schwarzschild radial coordinate r: whilst the radial Kerr Schild coordinate is equal to the Schwarzschild radius (note this is for the non spinning case). The form of the metric in the standard 3 + 1 ADM decomposition is then: where α = (1 + 2M/r) −1/2 (70) where x i = x i are the cartesian coordinates on the grid and r 2 = x 2 + y 2 + z 2 . The trace of the extrinsic curvature, the only other component required for the evolution of the fields, is given by

B.2 Matter field evolution
We have two complex scalar fields, which represent the mass eigenstates of the flavour fields. For each field the real part φ r evolves as: where the second order Klein Gordon equation has been decomposed into two first order equations as is usual for numerical evolutions. In the mass eigenstate basis, the mass matrix is diagonal and so the potential gradient is simply given by The imaginary parts follow the same evolution, with φ r → φ i .

B.3 Numerical details and convergence
GRChombo is a multi-purpose numerical relativity code, which uses adaptive mesh refinement to increase resolution in areas of interest, and is fully parallelised using the MPI framework. Whilst the full NR capabilities were not required for this work, its flexibility, scaling and built in GR tools made it easy to adapt for the fixed background case.

B.3.1 Discretisation in space and time
The metric values and their derivatives are calculated exactly at each point using the analytic expressions above, whereas the flavour field derivatives use 4th order finite difference stencils of the grid values and a 4th order Runge-Kutta time integration for the evolution equations. We use symmetric stencils for spatial derivatives, except for the advection derivatives (of the form β i ∂ i F ) for which we use one-sided/upwinded stencils. The length of the domain is L = 64R S , and we use four (2:1) refinement levels with the coarsest having 768 3 grid points. This coarsest resolution must be sufficient to resolve the wavelengths of the plane wave flavour fields (λ 2R S ) far from the BH. Around the BH additional resolution is required to properly resolve the horizon area. Kreiss-Oliger dissipation is used to control numerical errors in the evolution of the matter fields, in particular that arising at the grid boundaries.

B.3.2 Boundary conditions
We use periodic boundary conditions in all three directions to be able to simulate plane waves travelling in the x-direction. Ideally one would allow the waves to enter and exit at the ±x boundary by using appropriate ingoing and outgoing boundary conditions, thus minimising the boundary reflections which will inevitably disturb the results. A "lazier" alternative is simply to put the boundaries further away from the region of interest and then only "trust" the region which is causally disconnected from the boundaries at each point in the simulation. Thus the simulation is stopped at t = L/4 and only the r < L/4 part is considered valid.
For future simulations we plan to develop the appropriate boundary conditions, as this will potentially double the size of the region which we are able to accurately evolve. Figure 4: Here we show that the spatial profile of the probability P is converging at approximately 4th order, by comparing the difference between the results at 3 successive resolutions, which correspond to base resolutions of 384 3 , 512 3 and 768 3 , each with 4 levels of 2:1 refinement on top.

B.3.3 Convergence and testing
We have checked convergence of the simulations, as illustrated in Fig. 4.
We also checked that in the flat space case, the probability remained spatially uniform and the oscillation period was that given by eq. (10) (∼ 2000R s for the values used), as expected.
Finally we checked that the solution was unchanged when the boundaries were put twice further out, with L = 128 and a coarsest resolution of 1024 3 , to confirm that the reflections resulting from the periodicity were not affecting the results within the central area for which we show results.