# Black hole quasibound states from a draining bathtub vortex flow

###### Abstract

Quasinormal modes are a set of damped resonances that describe how an excited open system is driven back to equilibrium. In gravitational physics these modes characterise the ringdown of a perturbed black hole, e.g. following a binary black hole merger. A careful analysis of the ringdown spectrum reveals the properties of the black hole, such as its angular momentum and mass. In more complex gravitational systems the spectrum might depend on more parameters, and hence allows us to search for new physics. In this letter we present a hydrodynamic analogue of a rotating black hole, that illustrates how the presence of extra structure affects the quasinormal mode spectrum. The analogy is obtained by considering wave scattering on a draining bathtub vortex flow. We show that due to vorticity of the background flow, the resulting field theory corresponds to a scalar field on an effective curved spacetime which acquires a local mass in the vortex core. The obtained quasinormal mode spectrum exhibits long-lived trapped modes, commonly known as quasibound states. Our findings can be tested in future experiments, building up on recent successful implementations of analogue rotating black holes.

Introduction. The linear response of a perturbed black hole can be well-understood in terms of damped resonances called quasinormal modes (QNMs) qnmreview; Konoplya:2011qq. These damped modes possess a complex frequency whose real part corresponds to the oscillation frequency and whose imaginary part gives the lifetime. The QNM spectrum of a black hole is completely characterized by the black hole parameters, and does not depend on the initial conditions of the perturbations. With the on-going development of gravitational wave astronomy, significant efforts are being dedicated to use QNM spectra to reveal information about the structure of black holes, thereby allowing us to test General Relativity and alternative theories of Gravity implications. In this letter, we shall illustrate how the QNM spectrum can be drastically altered by additional structures in the context of Analogue Gravity.

Analogue Gravity, pioneered by Unruh in 1981 unruh1981, explores the possibility of testing gravitational effects in a broad variety of systems AGreview. For instance, it was shown in surfacewave that surface waves propagating on an inviscid, irrotational, and shallow fluid flow are equivalent to a scalar field propagating on an effective spacetime. This spacetime is completely determined by its propagation speed and the background fluid flow . In particular, it is possible to model a rotating black hole using an axisymetric fluid flow . The fluid configuration will exhibit an ergosphere at if and an event horizon at if . These two features are sufficient to give rise to many interesting effects that occur around astrophysical black holes, including Hawking radiation and superradiance AGreview. In particular, this allows us to experimentally test the universality and robustness of these effects. The past decade has seen an increase of interest in experimental realizations of analogue black holes, resulting in the measurement of both classical vanc_expmnt; faccio_hwkrad; germain_watertank and quantum steinhauer analogue Hawking radiation in (1+1)-dimensional systems. Experimental research on rotating, (2+1)-dimensional systems began more recently and is already bearing fruit. For example, superradiant amplification of surface waves was observed in a draining bathtub (DBT) vortex flow SR_obsvn) and rotating black hole analogues are also being explored using photon fluids prain_exp.

In this work, we present a simplified model of a rotational DBT-type fluid flow motivated by realistic velocity profiles SR_obsvn; stepanyants; andersen seen in experiments. Our model is also the analogue of a black hole with additional internal structure owing to the non-vanishing vorticity of the background flow in the centre of the vortex. Within our approximation, the vorticity causes perturbations to acquire a local effective mass close to the horizon. This is similar to certain studies of gravity, where a massless field can acquire a local mass by coupling to another field Thomas1; Thomas2. By studying the QNM spectrum, we find that our model also admits long-lived trapped modes in addition to usual QNMs. These are known in the literature as quasi-bound states (QBSs) qnmreview and are predicted to occur in particular for massive fields around Kerr black holes QBS_kerr; dolandirac. Our findings are of relevance for both hydrodynamics and gravity. Our aim is to understand the consequences of internal structures on the QNM spectrum. The question arises in gravity when one introduces modifications to the usual Kerr metric environment as well as in fluid flows where the core structure of a vortex deviates from the irrotational case. With Analogue Gravity experiments gaining an increasing amount of momentum, it is likely that the challenge of measuring the QNM spectrum will be tackled in the near future. In order to perform such experiments, it is necessary to understand the dependence of the QNM spectrum on the specifics of the background flow.

Background flow and wave equation. The DBT vortex is a particularly simple hydrodynamic model based on the assumption of a highly symmetric flow. It is described by the irrotational and incompressible velocity profiles and where (circulation) and (drain rate) are positive constants, and the polar coordinates. The associated effective metric mimics a rotating black hole as it possesses a horizon at , and an ergosphere at . Its QNM spectrum has been studied in the literature, see e.g. QNM_DBT_stab; QNM_DBT; DBT_resonances.

Even though the DBT model successfully describes the behaviour of the flow sufficiently far away from the centre, a vortex which forms under experimental conditions contains a core in which the flow is no longer irrotational. A more realistic formula for is given by the well known Rankine vortex lautrup,

(1) |

where is the radius of the vortex core and is the Heaviside step function. An analytically amenable interpolation of this formula, which is a smooth function of , was proposed by Rosenhead rosenhead and later studied by Mih mih2; vatistas; hydraulic,

(2) |

Notice that there exists many more complicated models, offering a more accurate description of the vortex flow depending on the precise initial and boundary conditions of the flow lautrup. In this letter, we are interested in the main deviations introduced by a rotational core. Hence we will work with Eq. (2) as its analytic simplicity will lend itself to our frequency domain simulations. Since we are dealing with a two-dimensional axisymmetric model, the radial component is constrained by the incompressibility condition, which leads to the same radial velocity as in the DBT vortex, i.e. . Using these velocity profiles, we investigate the QNM spectrum.

The equations governing an effective (2+1)-dimensional ideal fluid flow in the shallow water regime (single layer approximation buhble) are given by,

(3) |

These equations are valid in the regime , where is the height of the free surface and is the scale of variation in the -plane. Perturbations u to the background velocity can be expressed using a Helmholtz decomposition, , where the cograd operator is defined as and is the unit vector in the direction perpendicular to the (2+1) fluid ( can also be seen as the curl of the three dimensional vector field ). In the regime of short wavelengths (but still shallow water ), which amounts to a WKB approximation, the curl-free component obeys the wave equation,

(4) |

and the other component is obtained through . Moreover, the surface elevation can be obtained from

(5) |

In particular, the frequency content of and will be identical. Eq. (4) describes the propagation of a scalar field with a mass proportional to the background vorticity . We note that wave equation (4) becomes exact in two particular cases. The first is a solid body rotation with , in which case the perturbations are called inertia gravity waves buhble. The second is an irrotational flow, in which case and the wave equation reduces to its standard form gravwaves. The problem of waves scattering on a Rankine-type vortex has been addressed in the literature, usually in the regime kopiev; coste. In our case, the effects we are interested in arise in the regime where is of the order of .

Since the velocity profile we assume is axisymmetric, Eq. (4) is separable and a generic perturbation can be written as . To simplify, one can perform a Boyer-Lindquist type transformation (see e.g. (DBT_resonances)) and introduce a radial tortoise coordinate through . Using this coordinate, the horizon is located at , while at large , we have . The wave equation for a single frequency and azimuthal number reduces to,

(6) |

where the effective potential is given by,

(7) |

This differs from the usual potential since we are now using Eq. (2) for and the non-vanishing vorticity contributes the term , where . Since the potential is symmetric under the transformation {, }, we restrict ourselves to in the frequency domain. Hence, co- (counter-) rotating waves are defined by ().

The QNM boundary conditions are that the wave is purely in-going on the horizon and purely out-going at spatial infinity. By solving Eq. (6) in the corresponding limits and , the boundary conditions can be expressed as,

(8) |

where and are constants, and is the angular frequency at the horizon.

Computing quasinormal modes. Solving the wave equation subject to the boundary conditions (8) selects a discrete set of complex frequencies . Several methods are known in the literature to accomplish this objective qnmreview. We study the problem using three distinct methods. We start with a WKB method which allows us to make the distinction between the QNMs and QBSs and approximate their frequencies. Next, we implement a frequency domain simulation using a continued fraction method which allows us to accurately compute the quasinormal frequencies to 6-digit precision. Our last approach is a direct time domain simulation of Eq. (4), giving the time evolution of an initial perturbation of the vortex flow. Such a time evolution could be directly tested in future experiments. As we show all three methods give consistent results. (In addition, by setting , we have checked that all three methods reproduce the irrotational QNM spectrum for the DBT profile Cardoso04; DBT_resonances.)

WKB method. In the WKB regime (which is valid in particular for ) solutions are accurately described by a wave , where is a slowly varying amplitude and is the local wavenumber. Then obeys the Hamilton-Jacobi equation,

(9) |

with

(10) |

The two functions conveniently represent the potential (7) for varying (see Fig. 1). Indeed, we see that, at the level of the WKB approximation, if lies outside the range , hence the solution is oscillatory and hence the wave propagates. On the contrary, if is inside this range, where , the solution is evanescent. Where , there is a turning point. If in addition, the point is a local extremum of , it is an equilibrium point, i.e. the corresponding trajectory satisfies and . Since we consider , it is sufficient to consider the extrema of .

Near a local maximum, waves can hover around the vortex analogue of a light ring DBT_resonances; Awesome_paper. The real part of the QNM frequencies can be approximated through the conditions at the location such that where prime denotes the derivative with respect to . In the eikonal limit (), this condition is fulfilled at a single ( independent) radius , which is the light ring. Its orbital frequency governs the QNM spectrum according to the formula CardosoLR,

(11) |

where is called the overtone number and the term under the square root is evaluated at , . In the potential of Eq. (6), there are two maxima, and hence two associated families of QNMs.

In addition to the usual QNMs, the existence of two peaks in the scattering potential means that quasibound states (QBSs) can also exist in the system. These QBSs can be understood as trapped modes in the potential well that must tunnel across the peaks to decay. Hence, the characteristic life-times of these modes are significantly larger than that of QNMs.

To estimate the frequency of these QBSs, we perform a scattering amplitude calculation using WKB modes everywhere except at the turning points. We then construct a global solution using connection formulas at the turning points. When a QBS exists, the scattering potential contains four turning points (see Fig. 1). We define , and as the ranges spanned by the dip and the two peaks respectively; for example, . The real part of the QBS frequencies are given by the Bohr-Sommerfeld condition,

(12) |

where is the WKB action evaluated over the range and indexes the different energy levels in the dip. In contrast to QNMs, this relation is satisfied only by a finite number of ’s.

Once the real part is determined, the imaginary part of a long-lived mode with is given by,

(13) |

where the transmission coefficients across the potential barriers ( for the inner/outer) are given by
in the WKB regime and are the actions evaluated over the ranges . In Eq. (13) we take the sign if corresponding to a superradiant amplification on the inner barrier, and the sign otherwise. By comparing Eqs. (11) and (13), we see that the life-time of QBSs are typically exponentially large compared to that of QNMs. A full derivation of Eqs. (12-13) is outlined in the Supplemental Material.

Numerical simulations.
Our first numerical method is a direct time domain simulation of Eq. (4) using a gaussian pulse parallel to one of the cartesian axes as initial data. The wave equation is numerically integrated using the method of lines (MOL) and a 4th-order Runge-Kutta (RK4) algorithm. The quasinormal frequencies are extracted by performing a time Fourier transform of the signal once the initial pulse has passed. The second method is a frequency domain simulation implemented through a continued fraction method (CFM) leaverkerr; nollert. The two methods are detailed in the Supplemental Material.

Results. In Fig. 2 we present our main results. We use the MOL simulations to reconstruct the full 2-dimensional pattern (panel A) and its decomposition in azimuthal modes (panels B1-B5) at fixed time after the initial pulse has passed the vortex. We also present the real (panel C1) and imaginary (panel C2) components of complex frequency spectrum, obtained by the three different methods. We observe an excellent agreement between the results from the MOL (red crosses) and the CFM (blue diamonds) simulations. The WKB approximation then allows us to identify the underlying structure of the spectrum, i.e. whether the mode is a QNM (solid black dots) or a QBS (solid green dots). We show that most of the excited eigenmodes modes are QBSs, and as a result the ringdown has a significantly prolonged lifetime.

Conclusion. Our results show that the internal structure of a vortex can significantly affect its QNM spectrum. This provides a fluid analogue of the problem of spectral stability in black hole physics environment. In addition, the conclusions drawn in this paper should be observable experimentally by analyzing the ringdown of a vortex flow. If the core of the vortex lies behind the horizon (), then the observed spectrum will be that of an ideal DBT. On the contrary, if the core is large enough () our results show that the presence of QBSs will significantly alter the time response of the system, by allowing long-lived modes to hover around the vortex.

The ringdown we described above would manifest itself as rotating spirals imprinted on the free surface. In fact, such spirals are known to appear frequently around vortex flows (see Fig. 3 in Supplementary Material), for which a theoretical description seems to be lacking in the literature. We believe that the study of resonance frequencies could provide a fair description of this phenomenon.

###### Acknowledgements.

Acknowledgements. The authors would like to thank Emanuele Berti, Thomas Sotiriou and Théo Torres for illuminating discussions. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 655524. M.R. acknowledges partial support from the São Paulo Research Foundation (FAPESP), Grant No. 2013/09357-9, and from the Fulbright Visiting Scholars Program. M.R. is also grateful to the University of Mississippi for hospitality while part of this research was being conducted. S.W. acknowledges financial support provided under the Royal Society University Research Fellow (UF120112), the Nottingham Advanced Research Fellow (A2RHS2), the Royal Society Project (RG130377) grants, the Royal Society Enhancement Grant (RGF/EA/180286) and the EPSRC Project Grant (EP/P00637X/1). SW acknowledges partial support from STFC consolidated grant No. ST/P000703/.## References

Supplemental Material

## Appendix A Example of vortex spirals

## Appendix B WKB formula for QBS frequencies

The complex frequency of a QBS can be approximated by performing a tunnelling amplitude calculation using semi-classical wave mechanics berrymount. A solution to the wave equation can be approximated as a WKB mode everywhere except at the turning points of the scattering potential where (see Fig. 1). To find solutions near the turning point, one can make the expansion . Equation (6) becomes the Airy equation whose solutions can be expressed in terms of the orthogonal functions and . These take the asymptotic forms,

(14) |

where , with defined as the classically allowed region and the classically forbidden region. A global solution of Eq. (4) is constructed by relating WKB modes where is approximately linear, to Airy functions of large arguments through Eq. (14). Let () be the right (left) moving mode in the region and () the growing (decaying) mode in the region . Close to the turning point where , these modes take the form,

(15) |

Now, let be the globally defined solution which satisfies,

(16) |

Relating the two solutions to one another through Eq. (14) we see that,

(17) |

where is the transfer matrix relating exponential modes to the oscillatory modes. Eq. (17) is the well-known single-turning-point connection formula berrymount; Froman. Notice that there is an ambiguity in choosing the second solution of the Airy equation (this is the standard irreversibility of connection formulae issue Froman). Indeed, any combination of and with a non-zero coefficient for will have the same asymptotics as in Eq. (14), and that would lead to a different connection formula. That ambiguity is lifted by requiring that the second solution has a vanishing Wronskian (defined as ). This defines , and ensures that the connection formula (16) preserves the Wronskian conservation CoutantCF. In other words, . The complex conjugate relates exponential modes to oscillatory modes in the mirror situation. The inverses and relate oscillatory to exponential modes.

This can now be applied to the scattering problem for the potential in Eq. (7). To do so we combine the (single turning point) connection formula (16) with WKB propagation in between. For this, in each region where WKB is valid, we decompose the solution on a WKB basis as and relate the amplitudes in each regions using the connection formula (17). Then, the amplitudes of plane waves on both asymptotics are related by,

(18) |

where the matrices are the propagation matrices of the exponential modes underneath the peaks and the oscillatory modes in the dip, i.e.

(19) |

where . Note that is anti-diagonal owing to the fact that a growing mode in one direction is a decaying mode in the other direction. However by looking at the boundary conditions in Eq. (8), we see that if the mode will appear to propagate in the opposite direction and therefore will swap places with the on the left of Eq. (18). QBSs are defined by the condition . Solving Eq. (18) for this condition, one finds that,

(20) |

where the ’s are the reflection coefficients across the potential barriers, given by

(21) |

where the sign acts only on and is if the coefficient is superradiant, and otherwise. To extract from Eq. (20), we assume the QBS modes have a long lifetime, i.e. . We can then expand the action as . The real part then gives the Bohr-Sommerfeld condition, i.e. Eq. (12). The imaginary part gives the life-time as,

(22) |

Notice that is the semiclassical time for a wave packet to go from one turning point to the other in the potential well. Moreover, in the semiclassical limit, one has , and hence the log in the above equation can be expanded, leading to the expression used in the core of the text, see Eq. (13).

## Appendix C Numerical simulations in the time domain

The wave equation (4) is a partial differential equation (PDE) in for the field . However, we can exploit the symmetry of the system to write in which case the wave equation becomes a second order PDE in for . To solve this equation we use the Method of Lines (MOL). The first step to discretise equation (4) by approximating the -derivatives using 5-point finite difference (FD) stencils, hence converting the differential equation into a matrix equation. Boundary conditions (BCs) in are implemented at the edge of the FD matrix. The first BC is placed just inside the horizon and is left free which is achieved using a one-sided stencil. This is justified since the horizon acts as a one-way membrane, hence the value of the field inside the horizon cannot affect the exterior. The second BC is a hard wall placed far from the vortex centre. The wall is placed sufficiently far away that reflections do not enter the region of interest within the timescale of the simulations. To solve in time we use a fourth order Runge-Kutta method (RK4) with a gaussian pulse as initial condition, i.e.

(23) |

which is centred on , where gives the spread, and is an arbitrary amplitude. The first time derivative is imposed such that the pulse propagates toward the vortex, i.e.

(24) |

The corresponding conditions on the -components are calculated through,

(25) |

and similarly for . The wave equation is then used to compute at the next time step. Once the pulse has passed the vortex, is purely outgoing and the response of the system can be observed.

The frequency of the response is extracted using the time Fourier transform of , i.e. , over a range . If is of the form , then is given by,

(26) |

We fit this model to our data, using non-linear regression to obtain the values of and .

## Appendix D Continued fraction method

The continued fraction method, in the context of black hole perturbations, was developed by Leaver leaverkerr; leaverrn and improved by Nollert nollert. The procedure, in the present context, consists in writing the solution of (4) as a Frobenius expansion

(27) |

For (27) to be a solution of the wave equation (4), the coefficients have to satisfy a -term recurrence relation. (For the DBT vortex, the corresponding recurrence relation has only four terms Cardoso04). By successive Gaussian eliminations, we can reduce the -term recurrence relation to a -term recurrence relation of the form

(28) |

where , , and are coefficients that depend on and .

Using the change of variable , expression (27) above is compatible with the boundary conditions given in (8) as long as the infinite sum is everywhere convergent. The ratio test guarantees that the sum converges if

(29) |

From the recurrence relation, one can show that and, therefore, the sum converges for all , independently of the frequency . To assess convergence at spatial infinity, however, we need a theorem due to Pincherle gautschi; Cook14. This theorem guarantees that convergence at spatial infinity only occurs when the recurrence coefficients satisfy the following equation,

(30) |

where the expression on the left-hand-side is a continued fraction, written in standard notation (this explains the naming of the method). For a given value, the equation above picks a discrete set of frequencies that guarantee convergence of (27) and, therefore, correspond to quasinormal modes or quasibound states of the problem under analysis.

The numerical implementation of the method follows closely the algorithm described in Secs. III. H and III. I of Ref. Konoplya:2011qq. To solve (30), we use a damped Newton’s method with initial guess given by the WKB method estimates. In order to guarantee at least 6-digit precision for the all the calculated frequencies, we truncate the continued fraction in (30) at 25000 terms.

## Appendix E Comparison of results obtained from WKB, CFM and MOL methods

In Fig. 4, we show the complex frequency spectrum for , but displaying all QBSs and QNMs (no overtones) even those that weren’t excited in the time domain simulation. Moreover, we consider both signs of and . This provides a better understanding of the structure of the spectrum. As we see, the QBSs branch off from the QNM line associated with the outer light ring (see Fig. 1). In addition, the symmetry is made manifest. In Table 1, we show the numerical values of modes observed in Fig. 2, in order to compare the three methods more precisely.

WKB | CFM | MOL | |
---|---|---|---|

+2 | |||

-1 | |||

-3 | |||

-4 | |||

-5 | |||

-6 |