Backreaction in an analogue black hole experiment

TL;DR

This paper demonstrates that a fluid analogue of a black hole can exhibit backreaction effects, where the system's global parameters evolve due to scattering processes, enabling the study of dynamic black hole spacetimes in laboratory settings.

Contribution

It introduces a rotating draining vortex flow model showing that fluid systems can exhibit backreaction effects, linking scattering processes to global parameter changes.

Findings

Memory effect encodes scattering in the total mass of the system.

Backreaction causes significant global changes in background parameters.

Dynamical metric captures wave-background interactions in evolving black hole analogues.

Abstract

In general relativity, the interaction between a black hole and the fields around it (a process known as backreaction) proceeds via the evolution of the black holes mass and angular momentum. Analogue models of gravity, particularly fluid mechanical analogues, have been very successful in mimicking the propagation of fields, and the effects they experience, around black holes. However, hydrodynamic black holes are externally driven systems whose effective mass and angular momentum are set by experimental parameters, and as such no significant internal backreaction processes are expected to take place. We show, using a rotating draining vortex flow, that a fluid system of finite size exhibits a memory that keeps track of scattering processes in the system. This memory is encoded in the total mass of the system and hence, the backreaction arises as a significant global change in the background parameters, as opposed to a small local correction. More importantly, this backreaction is encapsulated by a dynamical metric, raising the possibility of studying wave-background interaction around evolving black hole spacetimes.

Tables

Table 1. Table 1: Details of the different experiments performed to measure the height change Δ ​ H Δ 𝐻 \Delta H . Q 𝑄 Q is the average flow rate over the course of the experiment; H 0 subscript 𝐻 0 H_{0} is the initial water height; f 𝑓 f and a 𝑎 a are the incident wave frequency and amplitude respectively; Δ ​ t Δ 𝑡 \Delta t is the window of wave incidence; H ˙ 0 subscript ˙ 𝐻 0 \dot{H}_{0} is the gradient of the linear height change at early times; τ 𝜏 \tau is the decay of the exponential decrease at late times and t c subscript 𝑡 𝑐 t_{c} is approximately the time when the behaviour switches from linear to exponential. The models used to extract these parameters are stated in in Eqs. (4) and (7). The uncertainties on the parameters are computed using the residuals from the fitting.

Exp.	$Q\left[\ell /\min \right]$	$H_0$ [mm]	$f$ [Hz]	$a$ [mm]	$\Delta t$ [s]	${\dot{H}}_0$ [mm/s]	$\tau$ [s]	$t_c$ [s]
$1$	$14.0\pm 0.4$	$20\pm 1$	$4$	$1.6\pm 0.1$	$120$	$-0.0118\pm 0.0002$	-	-
$2$	$14.0\pm 0.4$	$19\pm 1$	$2,3,4$	$2.3\pm 0.5$	$120$	$-0.0098\pm 0.0002$	-	-
$3$	$29.4\pm 0.4$	$65\pm 1$	$4$	$2.1\pm 0.5$	$600$	$-0.0389\pm 0.0005$	$94\pm 4$	$383\pm 3$

Equations

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }(\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu }\varphi )=0,$$

$$\dot{M}=-{\int }_{\gamma }{\rho }_{f}h{\mathbf{v}}_{\perp }\cdot d\mathbf{l},$$

$${\dot{H}}_0\simeq -\frac{1}{\mathcal{A}}{\int }_{\gamma }h{\mathbf{v}}_{\perp }\cdot d\mathbf{l},$$

$$H(t)=H_0+{\dot{H}}_{0}t+\mathcal{O}(t^2).$$

$$f(t,\theta ,r,z)=\sum _m{f}_m(r,z)e^{im\theta -i\omega t},$$

$$\dot{H}_{0}=-\frac{2\pi}{\mathcal{A}}\sum_{m}\tfrac{1}{2}\mathrm{Re}[h_{m}^{*}\mathbf{v}_{m}]\cdot\mathbf{r}\big{|}^{r_{2}}_{r_{1}},$$

$$H(t)=H_{0}+\tfrac{1}{2}\dot{H}_{0}(t+t_{c})-\dot{H}_{0}\tau\log\big{(}2\cosh\big{[}\tfrac{t-t_{c}}{2\tau}\big{]}\big{)}.$$

$$g_{\mu \nu }=\left(\begin{array}{cc}-(gH-V_{\perp ,i}{V}_{\perp }^i)& -V_{\perp ,i}\\ -V_{\perp ,j}& {\delta }_{ij}\end{array}\right),$$

$${\int }_V{\partial }_t\rho dV=-{\iint }_{\mathcal{S}}\mathbf{J}\cdot d\mathbf{A}=-{\int }_V\nabla \cdot (\rho \mathbf{V})dV,$$

$${\int }_{\mathcal{A}}{\int }_0^H{\rho }_f\nabla \cdot \mathbf{V}dzdA=0.$$

$$\begin{split}\dot{M}=\rho_{f}\int_{\mathcal{A}}\partial_{t}H~{}dA=-\int_{\mathcal{A}}\Big{[}\mathbf{J}_{\perp}\cdot\bm{\nabla}_{\perp}H-J_{z}(z=0)\\ +\int^{H}_{0}\bm{\nabla}_{\perp}\cdot\mathbf{J}_{\perp}dz\Big{]}dA,\end{split}$$

$$\dot{M}=-{\int }_{\mathcal{A}}{\nabla }_{\perp }\cdot (H{\mathbf{J}}_{\perp })dA=-{\int }_{\gamma }H{\mathbf{J}}_{\perp }\cdot d\mathbf{l},$$

$$H\to H+ϵh,\mathbf{V}\to \mathbf{V}+ϵ\mathbf{v},$$

$$\begin{split}\int_{\mathcal{A}}\Big{[}\partial_{t}H+&\bm{\nabla}_{\perp}\cdot(H\mathbf{V}_{\perp}+\epsilon^{2}h\mathbf{v}_{\perp})\Big{]}dA\\ &=\epsilon\int_{\mathcal{A}}\Big{[}\partial_{t}h+\bm{\nabla}_{\perp}\cdot(H\mathbf{v}_{\perp}+h\mathbf{V}_{\perp})\Big{]}dA,\end{split}$$

$${\dot{H}}_0\simeq -\frac{1}{\mathcal{A}}{\int }_{\gamma }h{\mathbf{v}}_{\perp }\cdot d\mathbf{l},$$

$$\frac{1}{2}{\mathbf{V}}_{\perp }^2+gH=\mathrm{const},{\nabla }_{\perp }\cdot (H{\mathbf{V}}_{\perp })=0.$$

$${\mathbf{V}}_{\perp }=-D/r{\mathbf{e}}_r+C/r{\mathbf{e}}_{\theta },$$

$$({\partial }_t+{\mathbf{V}}_{\perp }\cdot {\nabla }_{\perp }{)}^2\varphi -c^2{\nabla }_{\perp }^2\varphi =0,$$

$$r_h=D/c,r_e=\sqrt{C^2+D^2}/c.$$

$$\begin{array}{rl}& (c^2{r}^2-D^2)r^2{\partial }_r^2{\varphi }_m\\ +& \left[D^2+c^2{r}^2-2iD\left(\omega r^2-mC\right)\right]r{\partial }_r{\varphi }_m\\ +& \left[(\omega r^2-mC{)}^2-2imCD-m^2{c}^2{r}^2\right]{\varphi }_m=0.\end{array}$$

$$\begin{array}{r}{\varphi }_m(r)=A_m^h\left[1-\frac{{\tilde{\omega }}^2-\frac{2imCD}{r_h^4}-\frac{m^2{c}^2}{r_h^2}}{2c\left(\frac{c}{r_h}-i\tilde{\omega }\right)}(r-r_h)+...\right],\end{array}$$

$$\begin{array}{rl}{h}_m(r_h)=& i\tilde{\omega }{g}^{-1}{\varphi }_m(r_h)+c{g}^{-1}{\partial }_r{\varphi }_m{|}_{r=r_h},\\ {\mathbf{v}}_m(r_h)=& {\partial }_r{\varphi }_m{|}_{r=r_h}{\mathbf{e}}_r+im{r}_h^{-1}{\varphi }_m(r_h){\mathbf{e}}_{\theta },\end{array}$$

$${\dot{H}}_0=-\frac{\pi }{4gc\mathcal{A}}\sum _{m}F|A_m^h{|}^2,$$

$$F=\frac{{\tilde{\omega }}^4-\frac{c^8}{D^4}\left(4\omega m\frac{C}{c^2}+m^4\right)}{{\tilde{\omega }}^2+\frac{c^4}{D^2}}.$$

$$\begin{array}{r}\omega -{\mathbf{V}}_{\perp }\cdot \mathbf{k}=\pm ck(1-k^2{H}^2/6+...),\end{array}$$

$$\begin{array}{rl}\mathcal{L}=& -{\tilde{g}}^{\mu \nu }{\partial }_{\mu }\varphi {\partial }_{\nu }{\varphi }^{\ast }\\ =& |D_t\varphi {|}^2-c^2|{\nabla }_{\perp }\varphi {|}^2,\end{array}$$

$$\begin{array}{rl}{\tilde{g}}^{\mu \nu }=& \left(\begin{array}{cc}-1& -V_{\perp }^i\\ -V_{\perp }^j& c^2{\delta }^{ij}-V_{\perp }^i{V}_{\perp }^j\end{array}\right)\\ =& -{\delta }_t^{\mu }{\delta }_t^{\nu }+(c^2{\delta }^{ij}-V_{\perp }^i{V}_{\perp }^j){\delta }_i^{\mu }{\delta }_j^{\nu }-2V_{\perp }^i{\delta }_{(i}^{\mu }{\delta }_{t)}^{\nu }.\end{array}$$

$${\partial }_{\mu }\left(\frac{\partial \mathcal{L}}{\partial ({\partial }_{\mu }\varphi )}\right)-\frac{\partial \mathcal{L}}{\partial \varphi }=0.$$

$${\partial }_{\mu }{j}^{\mu }[\varphi ]=0\Rightarrow {\partial }_t\rho [\varphi ]+{\nabla }_{\perp }\cdot \mathbf{j}[\varphi ]=0,$$

$${\varphi }^a\to {\varphi }^a+\delta {\varphi }^a\Rightarrow \mathcal{L}\to \mathcal{L}+\delta \mathcal{L}.$$

Full text

Introduction. Analogue gravity, pioneered by Unruh in 1981 Unruh (1981), is a research programme which studies gravitational phenomena from general relativity (GR) using a wide variety of non-gravitational systems (see Barcelo et al. (2011) for a review). Unruh originally considered the propagation of sound waves through a fluid, and argued that if the fluid becomes supersonic in some region, the system exhibits a dumb hole horizon - the analogue of a black hole horizon. More generally, he showed that the wave equation describing the propagation of linear fluctuations $\phi$ through an ideal fluid is equivalent to the Klein-Gordon (KG) equation,

[TABLE]

describing the propagation of a massless scalar field on an effective curved spacetime. The effective metric $g_{\mu\nu}$ describing such a spacetime is completely determined by the properties of the fluid flow under consideration, raising the intriguing possibility of studying general relativistic wave phenomena in the laboratory.

Many different analogues have since been investigated in a variety of condensed matter systems Rousseaux et al. (2008); Steinhauer (2014, 2016); Belgiorno et al. (2010); Vocke et al. (2018); Fifer et al. (2019), including surface waves on top of a shallow fluid Schützhold and Unruh (2002). Although the analogy was originally conceived to investigate the trans-Planckian problem associated with Hawking radiation Unruh and Schützhold (2005), analogue gravity has enjoyed a number of other successes: notably surface wave experiments have been used to measure Hawking radiation Euvé et al. (2016); Weinfurtner et al. (2011, 2013), superradiance Torres et al. (2017), and quasi-normal ringing Torres et al. (2020).

One particularly simple model of a rotating black hole is provided by surface waves propagating on a rotating, draining fluid flow - the so-called draining bathtub vortex (DBT). Much work has been devoted to understanding features of this model (e.g. Basak and Majumdar (2003a, b); Cardoso et al. (2004); Berti et al. (2004); Anacleto et al. (2011); Dolan et al. (2012); Dolan and Oliveira (2013); Richartz et al. (2015); Dempsey and Dolan (2016); Churilov and Stepanyants (2019); Banerjee et al. (2018)), most of which relies on the assumption of an inviscid, incompressible, irrotational fluid in shallow water. Modifications resulting from the violation of these last two assumptions have been considered Patrick et al. (2018); Torres et al. (2018), and black hole effects (superradiance Torres et al. (2017) and quasinormal ringing Torres et al. (2020)) have been shown to persist under such conditions.

Another key assumption, which all analogue systems rely upon, is that the waves propagate on a fixed background. This assumption is necessarily violated in both hydrodynamical as well as graviational systems, since the fluctuations drive the evolution of the background through non-linear terms in the equations of motion. This process in known as the backreaction. The usual justification for neglecting backreaction is that the non-linear terms appear at quadratic order in perturbation theory and have little influence on the fluctuations, which are studied at linear order. However, since these terms can grow in time, they will eventually become important in determining the dynamics of the background. In atmospheric physics, this has widely been studied under the name wave-mean interaction theory Bühler (2014) to predict large scale changes in atmospheric currents resulting from small perturbations. In studies of surface fluctuations, it has long been recognised that waves produce a second order mass flux in the direction of wave propagation Phillips (1966), and that this mass flux induces a drift velocity called the Stokes drift Kenyon (1969); Hasselmann (1963). Despite recognition from the fluid dynamics community, these effects have yet to be incorporated into the analogue gravity formalism.

In this letter, we study the backreaction in an analogue black hole simulator by scattering surface waves with a DBT vortex. Whilst the backreaction is to be expected in any non-linear system, it normally proceeds via a small local change in the background parameters. By contrast, we show that our system exhibits a significant global change in one of the background parameters: the total fluid mass. We demonstrate theoretically that this can be explained by way of the wave induced mass flux across the open boundary of a finite sized system. Due to the analogy with black hole physics, the evolving background induces changes in the effective metric, as energy and angular momentum are exchanged between the waves and the analogue black hole. Therefore, one could in principle use a system of this type to investigate backreaction due to superradiance and Hawking radiation.

Theory. The system we are considering is a stationary draining fluid flow, where water enters a tank via an inlet and exits at a drain in a continuous cycle. Such a setup has previously been employed as a simulator for black hole superradiance Torres et al. (2017) and ringdown Torres et al. (2020) processes. For simplicity, we assume cylindrical symmetry about the drain and adopt polar coordinates $(r,\theta,z)$ centred on the drain. Water occupies a region $z\in[0,H]$ in the vertical direction. The area of the tank $\mathcal{A}$ (in the direction perpendicular to the vertical) is assumed constant and is bounded by the surface $\gamma$. For an incompressible fluid with density $\rho_{f}=\mathrm{const}$, the system is fully characterised by the water height $H$ and velocity field $\mathbf{V}$. Let $\gamma$ be comprised of the surfaces $r=r_{1}$, encircling (and nearby) the drain, and $r=r_{2}$, which forms the outer wall of the tank. We assume the inlet condition is specified on a small section of $r=r_{2}$, but any $\theta$ dependence introduced is confined to a small layer at the edge of the tank which we neglect.

Now consider fluctuations (i.e. surface waves) described by $(h,\mathbf{v})$, which are switched on at $t=0$ and are created inside $\mathcal{A}$. It is well-known Phillips (1966) that linear surface waves propagating on a background flow produce a mass flux $\mathbf{j}$ in the direction of wave propagation (which in our case is in the plane $\perp$ to the vertical). If the background $(H,\mathbf{V})$ is stationary for $t<0$, then shortly after the onset of waves, the amount of mass $M$ contained within the system will be altered if there is a net mass flux $\mathbf{j}_{\perp}=\rho_{f}\mathbf{v}_{\perp}$ over $\gamma$, according to

[TABLE]

where the overdot denotes the time derivative. The assumptions involved in the derivation of this formula are detailed in Appendix A. Furthermore, since $M=\rho_{f}\int_{\mathcal{A}}HdA$, a change in total mass results in a change in the water height. If we assume that $H$ is approximately level over $\mathcal{A}$ (i.e. spatially uniform), then the water height adjusts according to,

[TABLE]

where subscript [math] indicates that strictly a quantity is evaluated at $t=0$. Expanding $H(t)$ around $t=0$ gives,

[TABLE]

The $\mathcal{O}(t^{2})$ corrections become significant once changes to the background become large, which means that terms quadratic in the background variables also contribute to Eq. (3). Once this happens, the $t$ dependence of $H$ and $\mathbf{V}_{\perp}$ becomes interlinked and to solve the coupled system, one must use a second equation (the equation of momentum conservation) in addition to the equation of mass conservation which is used in the derivation of Eq. (2). Although we do not attempt this here, we expect this to produce exponential behaviour of $H(t)$ at late times, since $\dot{H}$ will depend on the value of $H$.

Near $t=0$, fluctuations of frequency $\omega$ perceive a quasi-stationary background and may be written,

[TABLE]

where $f$ is a placeholder for $(h,\mathbf{v})$ and we have also used our assumption of cylindrical symmetry to make a decomposition into modes of azimuthal number $m\in(-\infty,\infty)$. Using this ansatz, Eq. (3) becomes

[TABLE]

where $*$ denotes the complex conjugate and we have dropped the subscript $\perp$ on the velocity perturbation for conciseness. At $r=r_{2}$, $\mathbf{v}_{m}\cdot\mathrm{\mathbf{e}}_{r}$ vanishes everywhere except at the inlet, where we assume a rapid influx of water (exceeding the propagation speed of the fluctuations) required to drive the high flow velocities. Therefore, since the fluctuations are generated inside $\mathcal{A}$, the total mass flux receives no contribution at $r=r_{2}$ and is determined solely by the form of the fluctuations at $r=r_{1}$. In Appendix B, we show how Eq. (6) can be evaluated in the shallow water regime for a DBT vortex. This system exhibits an effective horizon at $r=r_{h}$ which we take as our inner boundary.

Experiments. To test the prediction of Eqs. (4) and (6), as well as to observe the time-scales over which they are valid, we scatter monochromatic surface waves with a DBT vortex in a controlled experiment. The vortex is generated by pumping water at flow rate $Q$ into a rectangular ($2.65~{}\mathrm{m}\times 1.38~{}\mathrm{m}$) tank and allowing it to drain through a hole of radius $d=2~{}\mathrm{cm}$ located in the centre. Once the vortex is in equilibrium (determined by the constancy of $Q$ and $H_{0}$) monochromatic surface waves are generated using a series of electrically controlled pistons. The change in water height $\Delta H=H-H_{0}$ is determined by illuminating the free surface with a laser sheet and tracking it’s average position with a high-speed camera. More details of the method can be found in Appendix C.

In Table I of Appendix C, we summarise the parameters used in each experiment. We display the measured $\Delta H(t)$ profiles in Fig. 1. In experiment 1 (panel A) we verify the linear decrease of $H$ at early times predicted by Eq. (4). In experiment 2 (panel B) we test the frequency dependence of Eq. (6) by varying the frequency from 2 to 4 Hz, finding no significant variation of $\dot{H}_{0}$ within this range. In experiment 3 (panel C) we demonstrate that at late times the behaviour of $H(t)$ deviates from a linear decline.

For the results of the first two experiments, we fit $\Delta H(t)$ with the linear decline predicted in Eq. (4) to provide a value for $\dot{H}_{0}$. Since we do not have a prediction for $H(t)$ over the full evolution, we devise a phenomenological model which has the expected behaviour at early and late times, i.e. linear then exponential,

[TABLE]

The extra parameters are the time $t_{c}$ at which the evolution switches from linear to exponential, and the exponential decay time $\tau$. Strictly speaking, this model has linear asymptotics in the limit $t\to-\infty$. However, the behaviour at $t=0$ is well approximated as linear provided $t_{c}>\tau$. The parameters obtained from the fits are also contained in Table I.

Discussion. Across all experiments, we find linear behaviour at early times, lending support to our prediction in Eq. (4). This behaviour persists whilst $\Delta H\ll H_{0}$ as expected, which can be seen by comparing $\Delta H$ in Fig. 1 with the corresponding value of $H_{0}$ in Table I. The observed lack of frequency dependence in the second experiment is supported by computation of Eq. (6) assuming a shallow water flow (see Appendix B), which predicts that the value $\dot{H}_{0}$ only varies by $\sim 1\%$ over the range $f\in[2,4]~{}\mathrm{Hz}$. Furthermore, by comparing the gradients across the different experiments we see that $\dot{H}_{0}$ depends on $H_{0}$, supporting our claim that the late time behaviour should be exponential. This is further evidenced by the late time tail in experiment 3.

In Fig. 2, we compare a simulation of the shallow water wave equation (see Eq. (B5) in Appendix B) for $H_{0}=2~{}\mathrm{cm}$ and $f=2~{}\mathrm{Hz}$ with the scattered wave from experiment 2 for the same frequency. A mathematical justification of our inviscid, shallow water treatment for this particular experiment can be found in Appendix C. This treatment is corroborated by the clear similarities between experiment and simulation in Fig. 2. The scattered wave in our set-up was measured using the air-water interface sensor described in Torres et al. (2017) and for the simulation, we used the flow parameters $C=0.013~{}\mathrm{m^{2}/s}$ and $D=0.001~{}\mathrm{m^{2}/s}$ with the velocity profile in Eq. (B2) of Appendix B. These values were chosen to be similar to that found from the flow measurements of Torres et al. (2017, 2020), which involved the same experimental apparatus, and the precise values were tuned to match the number of wavefronts on the left and right sides of the images 111Changing the values of $C$ and $D$ by 10% resulted in notable visual discrepancies between the two scattering patterns in Fig. 2. The predicted value of $\dot{H}_{0}$ changed by about 60% when using these values.. Using the form of the perturbations on the horizon, which requires knowledge of the transmission coefficient (see Appendix B), Eq. (3) predicts $\dot{H}_{0}=-0.011~{}\mathrm{mm/s}$. Comparing with the measured value $\dot{H}_{0}=-9.8\pm 0.2\times 10^{-3}~{}\mathrm{mm/s}$ from in Table I, we find good agreement, despite the violation of $H\neq H(r)$ and $\bm{\nabla}\times\mathbf{V}=0$ (also assumed by our prediction) in the vortex core where the horizon is located (see Appendix C for further discussion).

To relate the effect we have observed to the backreaction in GR, consider the following. For slowly evolving background, one can adopt different timescales for the background and the fluctuations (i.e. a Born-Oppenheimer approximation Born and Oppenheimer (1927)). In this approximation, which is also employed to study the backreaction in GR Balbinot et al. (2005a), the fluctuations of a shallow water, irrotational fluid obey Eq. (1), with the components of the evolving effective metric given by,

[TABLE]

where the indices $i$ and $j$ run over spatial dimensions, $\delta_{ij}$ is the Kronecker delta function and the $t$ dependence enters via $H(t)$ and $\mathbf{V}_{\perp}(t)=V_{\perp}^{i}(t)\mathbf{e}_{i}$. In our analysis, we have estimated only $H(t)$ a short time after the beginning of wave incidence and the full dynamics of $H(t)$ and $\mathbf{V}_{\perp}(t)$ warrant further investigation. Furthermore, since the energy and angular momentum density of the shallow water background flow are given by $\varepsilon=\tfrac{1}{2}H\mathbf{V}_{\perp}^{2}$ and $L=H\mathbf{r}\times\mathbf{V}_{\perp}$ respectively, a change in $H$ corresponds to changes in both $\varepsilon$ and $L$. Therefore, the height change mediates the exchange of energy and angular momentum between the waves and background.

Although the analogy persists at the linear level, the backreaction equation, describing how the effective metric evolves, is specific to the equations of motion of the system under consideration, which are the Euler equations for ideal fluids and the Einstein equations in GR (different backreaction equations have been studied, for example, in Balbinot et al. (2005b, a); Schützhold et al. (2005); Schützhold (2007)). Despite this, it may still be possible to learn about generic features of slowly backreacting spacetime geometries using analogue systems. In addition, it is possible extract scattering amplitudes from changes in global parameters (here the water height) for controlled scattering processes, e.g. scattering of waves with a single azimuthal wave number.

Conclusion. In this work, we have studied the backreaction of surface waves on a draining vortex flow. Our results demonstrate that surface waves interacting with an initially stationary vortex will trigger the evolution of the background out of equilibrium. Due to the flow being externally driven, it was previously unclear whether the background had the freedom to adjust to the presence of waves in analogue gravity simulators. Our findings show that the backreaction is indeed observable, and that the system does in fact have freedom to re-distribute energy and angular momentum between the incident waves and the background flow. In the shallow water regime, we have argued that this evolution is encapsulated by a dynamical effective metric. Although this metric does not evolve according to the Einstein equation, further similarities between slowly evolving gravitational and analogue spacetimes have yet to be investigated.

This realisation is important for a number of reasons. Firstly, one must ensure that any wave effects (e.g. stimulated Hawking radiation, superradiance and quasinormal ringing) are measured on a timescale much shorter than the time it takes for the background to change, so that the assumption of a stationary background is not violated. This may restrict the frequency range one can probe in an analogue gravity experiment, as non-linear effects will influence the low frequency behaviour which takes places over longer timescales. Secondly, the effect we have described is a global (as opposed to local) phenomenon. Thus, one can use the asymptotic value of the water height to obtain insight into scattering processes, similar to the role that the black hole mass plays in GR.

Based on previous experience, a similar behaviour is expected to occur in suitable quantum systems, and thus our findings suggest that analogue gravity experiments can be used to cross-validate backreaction models in a relativistic setting. This is an area of research where lack of experimental input is stalling theoretical development.

Appendix A

Starting from first principles, we derive a general formula for the rate of change of total mass contained in a system due to a non-zero mass flux over it’s boundary. Consider a volume $V$ bounded by a surface $\mathcal{S}$. In the absence of sources and sinks, a change in the density inside $V$ is compensated by the flow of mass over $\mathcal{S}$,

[TABLE]

where in the second equality we have used the definition of the mass flux $\mathbf{J}=\rho\mathbf{V}$, as well as the divergence theorem to convert the surface integral into one over the volume. The water tank we consider is bounded from below and open from above, $z\in[0,\infty)$, with water occupying the region $z\in[0,H]$ where $z=H$ is the free surface (or height) of the water. Thus the density may be represented by $\rho=\rho_{f}\Theta(H-z)$, where $\rho_{f}=\mathrm{const}$ is the fluid’s density and $\Theta$ is the Heaviside step function. Inserting into Eq. (A.1) yields,

[TABLE]

Here we have exploited the kinematic boundary condition at the free surface $V_{z}(z=H)=(\partial_{t}+\mathbf{V}_{\perp}\cdot\bm{\nabla}_{\perp})H$, where $\perp$ denotes the components perpendicular to the vertical coordinate $z$. The area $\mathcal{A}$ of the system in the $\perp$ plane is assumed constant. Splitting the divergence into $\perp$ and $z$ components, performing the $z$ integral and using the kinematic boundary condition, we arrive at,

[TABLE]

where we have used the fact that the total mass of the system is given by $M=\rho_{f}\int_{\mathcal{A}}HdA$. If the system is in equilibrium, all of the terms on the right hand side cancel leading to $\dot{M}=0$.

We can gain further intuition by considering the shallow water regime, where the expression for $\dot{M}$ simplifies significantly. The shallow water approximation consists in assuming $L\ll H$, where $L$ is a length scale characterising variations of the fluid in the $\perp$ plane. In this approximation, the fluid is treated as a single two dimensional layer, in which the $z$ dependence in $\mathbf{V}_{\perp}$ is neglected. Eq. (A.3) then reduces to,

[TABLE]

where $\gamma$ is the boundary of $\mathcal{A}$ and we have also assumed that $J_{z}(z=0)=0$ inside $\mathcal{A}$ (this is motivated by the set-up we consider in the main text). Thus, the system is in equilibrium if total mass flux over $\gamma$ is zero.

Now, let $(h,\mathbf{v})$ be linear fluctuations about the background, which are switched on at the time $t=0$. The presence of fluctuations contributes additional terms to the mass flux. Making the replacements,

[TABLE]

where $\epsilon$ is an order tracking parameter, Eq. (A.4) can be written,

[TABLE]

where we have gathered the terms at linear order in $\epsilon$ on the right hand side. Assuming that for $t<0$ background is in equilibrium, the presence of fluctuations will trigger the evolution of the background through the $\mathcal{O}(\epsilon^{2})$ term in Eq. (A.6). To separate out the background and fluctuation dynamics, once can introduce two new time variables $T_{1}=t$ and $T_{2}=\epsilon^{2}t$, such that $h=h(T_{1})$, $H=H(T_{2})$ and similarly for the velocity field. Furthermore, since $\bm{\nabla}_{\perp}\cdot(H\mathbf{V}_{\perp})=0$ just before $t=0$, the leading order contribution to this term will also be $\mathcal{O}(\epsilon^{2})$. Requiring that the governing equations be satisfied to each order in $\epsilon$, the left and right hand sides of Eq. (A.6) must both vanish respectively.

To compute change in $H$ close to $t=0$, we can drop the term $\bm{\nabla}_{\perp}\cdot(H\mathbf{V}_{\perp})$ since the background is still approximately in equilibrium. Assuming that $H$ is independent of spatial coordinates, the height change at early times is dictated by,

[TABLE]

which is the result stated in Eq. (3) of the main text. Note that the quantity under the integral is related to the Stokes drift velocity $\mathbf{v}_{S}$ for small waves in shallow water by $H\mathbf{v}_{S}=h\mathbf{v}_{\perp}$ Phillips (1966). Thus, we may interpret the mechanism behind the height change as the Stokes drift integrated over the boundary of the system. At later times, the $\bm{\nabla}_{\perp}\cdot(H\mathbf{V}_{\perp})$ term in Eq. (A.6) (which is replaced by the first and third term on the right hand side of Eq. (A.3) for generic water depth) must also be taken into account, meaning one has to solve the full backreaction equation.

The result in Eq. (3) is not specific to shallow water and only relies on the assumption of small waves, since the perturbed $z$ integral in Eq. (A.3) from $H$ to $H+\epsilon h$ reduces to a multiplication, resulting in the same $\mathcal{O}(\epsilon^{2})$ term as in Eq. (A.6).

Appendix B

In this appendix, we deal with a $t$ independent background to derive $\dot{H}_{0}$ in Eq. (3). We can make considerable progress analytically by assuming a shallow fluid which is also irrotational, i.e. $\bm{\nabla}\times\mathbf{V}=0$. This second assumption means that the velocity field represents only a single degree of freedom $\Phi$ satisfying $\mathbf{V}=\bm{\nabla}\Phi$. Similarly, the fluctuations satisfy $\mathbf{v}=\bm{\nabla}\phi$. The equations for the stationary background are,

[TABLE]

The first of these is Bernoulli’s equation resulting from momentum conservation, and the second comes from Eq. (A.4). An approximate solution to Eq. (B.1) far from the vortex core where $H\simeq\mathrm{const}$ is,

[TABLE]

which is known in the literature as the draining bathtub vortex (DBT), and $C>0$ and $D>0$ are the circulation and drain constants. The perturbed shallow water equations are then combined into a wave equation for the fluctuations,

[TABLE]

where $c=\sqrt{gH}$ is the wave speed. The system exhibits an acoustic horizon $r_{h}$ and ergosphere $r_{e}$ (defined respectively as the boundaries of the regions satisfying $|\mathrm{\mathbf{e}}_{r}\cdot\mathbf{V}|>c$ and $|\mathbf{V}|>c$) which are given by,

[TABLE]

Since the horizon is a one-way membrane for the perturbations, $r_{1}=r_{h}$ is a natural choice of inner boundary for the system. The solution for the fluctuations on the horizon then determines the form of $\dot{H}_{0}$. Using the separation ansatz in Eq. (5) for the field $\phi$, Eq. (B.3) becomes,

[TABLE]

This has a regular singular point at $r=r_{h}$ and thus a solution may be sought in terms of a Frobenius expansion about $r=r_{h}$ Churilov and Stepanyants (2019). Imposing that the solution remain finite, the first two terms of the series are,

[TABLE]

where $A_{m}^{h}$ and $\tilde{\omega}=\omega-mC/r_{h}^{2}$ are the mode amplitude and intrinsic angular frequency on the horizon respectively. To find $h_{m}$ and $\mathbf{v}_{m}$, we use the relations,

[TABLE]

which upon insertion into Eq. (6) yields,

[TABLE]

where $F$ (which controls the sign of $\dot{H}_{0}$) is a function of the flow and wave parameters given by,

[TABLE]

Note that for $|\tilde{\omega}|\gg c^{2}/D$ we have $F\simeq\tilde{\omega}^{2}$, which is the result one would naively expect from the eikonal approximation. Unlike the energy current, $\dot{H}_{0}$ does not change sign with $\tilde{\omega}$. This is a result of the mass flux only being one of several terms contained in the energy current when $\mathbf{V}\neq 0$ (see Appendix D for discussion).

Appendix C

Here we expand on the details of our experiments. To observe the height change resulting from impinging waves, we illuminate the free surface from above using a Yb-doped laser, mean wavelength $457~{}\mathrm{nm}$, which is converted into a thin laser sheet (thickness $\sim 2~{}\mathrm{mm}$). This appears on the free surface as a line spanning nearly the full length of the tank, thereby allowing us to see the free surface in our data. We filmed this line from the side with a high-speed Phantom Miro Lab camera at $24~{}\mathrm{fps}$. The line of sight of the camera was at an angle of $\bar{\theta}$ to the free surface such that $\cos\bar{\theta}=0.91\pm 0.01$, which was necessary to avoid shadows of waves passing in between the laser-sheet and the camera, which obscured our vision of the free surface. This induces an uncertainty on the measured height change which we include in our error estimates of $\dot{H}_{0}$. We recorded the height before sending any waves to confirm that the background was steady, and then monitored the water height whilst exciting waves of frequency $f=\omega/2\pi$ and amplitude $a$ over a time window $\Delta t$. The free surface was identified by finding the pixel of maximum intensity in each image and interpolating using adjacent points to determine the maximum to a sub pixel accuracy.

Once the free surface in each image is determined as a function of spatial coordinate, we extract the zero wavenumber contribution to the signal using a spatial Fourier transform. This gave a measurement of the average height across the observed region at each time step, which we correct for the camera angle $\bar{\theta}$, obtaining the variation of the background height over time. In all experiments, the initial height was determined to be sufficiently steady (less than $10\%$ of the total change over the experiment) and any slight variation was due to the difficulty in maintaining constant $Q$ over the experiment. We corrected for this (as well as the slightly different initial heights) by fitting the curve prior to the start of the waves with a straight line and subtracting this from the entire data set to give the height change $\Delta H$. This ensured that the height change measured was the result of the perturbations.

In all experiments, we observed a small amplitude oscillation about the linear behaviour at the frequency $f$ (confirmed by a peak in the Fourier transform of $\Delta H$) corresponding the oscillatory (linear) terms in Eq. (A.6). Furthermore, the amplitude of this oscillation decreased with $f$ which is expected, since integrating $\dot{H}_{0}$ in time brings out a factor $1/f$. We also observe a sharp increase in height immediately before the height starts to decrease. This is because our wave generator panel is initially fully retracted and must move forward to its equilibrium position, thereby reducing the area of the system and increasing mean height. In accordance with this explanation, a sharper increase was observed for larger piston amplitudes.

To justify a shallow water treatment, we make the following argument. The dispersion relation is obtained from the wave equation by making the replacements $i\partial_{t}\to\omega$ and $-i\bm{\nabla}_{\perp}\to\mathbf{k}$, where $\mathbf{k}$ is the wavevector in the $\perp$ plane. For $kH<1$, where $k=||\mathbf{k}||$, the dispersion relation for gravity waves Landau and Lifshitz (1987) may be expanded in powers of $kH$,

[TABLE]

where the first term on the right hand side corresponds to the wave equation in Eq. (B.3), and the second term is the leading order dispersive correction for $kH\neq 0$. To estimate the size of this term far from the vortex core where $||\mathbf{V}_{\perp}||\ll c$, take $k\simeq\omega/c$, which is valid up to corrections of $\mathcal{O}((kH)^{2})$ and $\mathcal{O}(||\mathbf{V}_{\perp}||/c)$. Using this estimate, only experiment 2 for $f=2,3~{}\mathrm{Hz}$ satisfies $kH<1$. For these frequencies, the dispersive correction in Eq. (C.1) is $0.05$ and $0.11$ respectively and thus, in particular for $f=2~{}\mathrm{Hz}$, our shallow water treatment seems reasonable.

A similar argument may be invoked to justify the assumption of an inviscid fluid. Indeed, corrections to the dispersion relation due to viscosity $\nu$ appear in powers of $\nu k/c$ Visser (1998); Ganguly (2017), which can be seen purely on dimensional grounds. For water at room temperature with $\nu=1~{}\mathrm{mm^{2}/s}$, $\nu k/c\sim\mathcal{O}(10^{-5}f)$ and can thus be safely neglected for all frequencies considered.

Finally, although the assumptions $H\neq H(r)$ and $\bm{\nabla}\times\mathbf{V}=0$ are satisfied far from the centre, they are violated in the vortex core. The effect of free surface gradients was studied in Unruh (2013); Coutant and Parentani (2014); Richartz et al. (2015) and it was shown that whilst the KG equation is preserved, $H=H(r)$ alters the form of the effective metric. Flows with vorticity, however, do not preserve the KG equation since the usual scalar mode $\phi$ becomes coupled to additional degrees of freedom Bergliaffa et al. (2004); Liberati et al. (2019). Since both effects become important in the vortex core where the drain is located, their inclusion into into the description of the fluctuations would be necessary to improve the prediction for $\dot{H}_{0}$. However, since both are satisfied far from the vortex core, the KG equation still holds in this region and we can still interpret wave propagation there using the analogy to GR.

Appendix D

Here we provide a discussion of the different notions of energy current and their relation to the mass flux. We begin with an action for the perturbations $\mathcal{S}=\int\mathcal{L}d^{2}xdt$, where the Lagrangian density is given by,

[TABLE]

where $\tilde{g}^{\mu\nu}=\sqrt{-g}g^{\mu\nu}$ is the inverse metric and the squared absolute value is obtained by taking the product of a quantity with its complex conjugate. We have also defined the material derivative $D_{t}=\partial_{t}+\mathbf{V}_{\perp}\cdot\bm{\nabla}_{\perp}$ for conciseness. The explicit form of $\tilde{g}^{\mu\nu}$ is,

[TABLE]

The wave equation in Eq. (B.3) is then obtained using the field theoretic version of the Euler-Lagrange equations,

[TABLE]

Note that applying Eq. (D.3) to the first line of Eq. (D.1) results in the KG equation of Eq. (1), whereas the second line yields directly the shallow water wave equation in Eq. (B.3).

A conserved current $j[\phi]$ is a quantity with components $j^{\mu}[\phi]$ satisfying,

[TABLE]

where the second form in terms of $(\rho,\mathbf{j})$ splits $j$ into the temporal and spatial parts which are the charge $\rho$ and current $\mathbf{j}$ respectively (not to be confused with the density and mass flux in the main text). Note, the two together are collectively called the 4-current. Square brackets are use to indicate that a quantity is a functional of $\phi$.

Conserved currents are derived as follows. Consider an infinitesimal transformation of the field which induces a shift in the Lagrangian,

[TABLE]

Noether’s theorem Srednicki (2007); Schwartz (2014) states that $j$ is a conserved current if the Lagrangian changes by a total derivative $\delta\mathcal{L}=\partial_{\mu}F^{\mu}$. The components of the current are given by,

[TABLE]

where $a=1,2$ with $\phi_{1}=\phi$ and $\phi_{2}=\phi^{*}$. We now consider two different conserved quantities: the norm and the energy.

Norm conservation. We consider first performing a phase rotation on $\phi$. Since the wave equation is linear in $\phi$, there is an internal symmetry $\phi\to\exp(i\alpha)\phi$, where $\alpha$ is a phase rotation, which leaves the equations of motion invariant. For infinitesimal $\alpha$, $\phi$ changes by $\delta\phi=i\alpha\phi$ and $\delta\mathcal{L}=0$. The conserved quantities associated with this transformation are the norm $\rho_{n}[\phi]$ and the norm current $\mathbf{j}_{n}[\phi]$ given by,

[TABLE]

Using the ansatz in Eq. (5) for an axisymmetric, quasistationary system, we have $\partial_{t}\rho_{n}=0$. Applying the divergence theorem to $\bm{\nabla}_{\perp}\cdot\mathbf{j}_{n}[\phi]=0$ over the region $r=[r_{h},\infty)$, we deduce $\int^{2\pi}_{0}[\mathbf{r}\cdot\mathbf{j}_{n}]^{\infty}_{r_{h}}d\theta=0$. Using asymptotic solutions to Eq. (B.5) (see e.g. Berti et al. (2004)), this integral yields,

[TABLE]

where the left (right) of the equation is the norm current expressed at infinity (the horizon) and $+/-$ denote out/in-going respectively.

Energy conservation. Time translation $t\to t-\delta t$ induces a change in the field $\delta\phi^{a}=\delta t\partial_{t}\phi^{a}$ and the Lagrangian $\delta\mathcal{L}=\delta t\partial_{t}\mathcal{L}$. This gives rise to conservation of the energy current, which has components,

[TABLE]

Defining $\partial_{t}\phi=\dot{\phi}$ and $p=\partial\mathcal{L}/\partial\dot{\phi}=D_{t}\phi^{*}$ (and similarly for the complex conjugate), one can see that $\rho_{e}=\mathcal{H}$ is the Hamiltonian density.

To make contact with fluid dynamics, we define $(E,\mathbf{I})=(\rho_{e},\mathbf{j}_{e})/2g$ and also $\phi\in\mathbb{R}$. Using $D_{t}\phi=-gh$ and $\bm{\nabla}_{\perp}\phi=\mathbf{v}_{\perp}$, we obtain the usual equation for wave energy conservation in shallow water,

[TABLE]

where,

[TABLE]

In fluid dynamics, these equations are obtained by contracting the linear shallow water equations with $(gh,H\mathbf{v}_{\perp}^{T})$ where superscript $T$ indicates the transpose. Following the same procedure outlined above, we evaluate the expression for $E$ to show $\partial_{t}E=0$, and applying the divergence theorem we obtain $\int^{2\pi}_{0}[\mathbf{r}\cdot\mathbf{I}]^{\infty}_{r_{h}}d\theta=0$. Evaluating this expression in both limits results again in Eq. (D.8) multiplied by a factor of $\omega c/g$, thus making contact between the norm current and the energy current.

Restricting our analysis to $\omega>0$, notice that if $|A_{m}^{+}|^{2}<|A_{m}^{-}|^{2}$ the energy current is negative (points towards $r=r_{h}$). However if $|A_{m}^{+}|^{2}>|A_{m}^{-}|^{2}$ (which occurs for $\tilde{\omega}<0$) the energy current becomes positive and points outward, allowing for the extraction of energy from the system to infinity (this is the phenomenon of superradiant scattering). This is understood as the mode carrying a negative energy across the horizon, $\rho_{e}(r_{h})<0$, thereby lowering the energy of the system.

As a brief aside, Eq. (D.9) applied to harmonic field modes (of the form in Eq. (5)) becomes,

[TABLE]

One can then observe that this energy $4$-current is $\omega$ multiplied by the norm $4$-current in (D.7). Moreover, Eq. (D.6) applied to angular coordinate translations leads to an angular momentum $4$-current $j^{\mu}_{l}$ satisfying the well-known property $j^{\mu}_{e}/j^{\mu}_{l}=\omega/m$.

Now to make contact with the mass flux, we notice that when $\mathbf{V}=0$ we have $\mathbf{I}=gHh\mathbf{v}_{\perp}$. Once integrated, this is precisely the term appearing in Eq. (3) multiplied by a factor $c^{2}$. Therefore, in the absence of a flow, there is a direct correspondence between the energy current in Eq. (D.11b) and the shallow mass flux $\mathbf{j}_{\mathrm{SW}}=h\mathbf{v}_{\perp}$ given by $\mathbf{I}=\mathbf{j}_{\mathrm{SW}}c^{2}$ (reminiscent of the celebrated result $E=mc^{2}$). However, when $\mathbf{V}\neq 0$, we see clearly from Eq. (D.11b) that the energy current and mass flux no longer coincide. This explains why Eq. (6) does not have the symmetry properties of $\tilde{\omega}$.

For the reader who is more familiar with the language of electromagnetism, we draw some comparisons here. The combination $\tfrac{1}{2}\mathrm{Re}[h_{m}^{*}\mathbf{v}_{m}]$ appearing in Eq. (6) is reminiscent of the complex Poynting vector, $\mathbf{S}$, whose real part describes the flux of power in electromagnetic waves, in terms of the phasor decomposition of the electric and magnetic fields ($\mathbf{E}$ and $\mathbf{H}=\mathbf{B}/\mu_{0}$ respectively), Meystre and Sargent (2007),

[TABLE]

The correspondence between $\mathbf{I}$ and $\mathbf{S}$ becomes exact upon defining,

[TABLE]

where $E_{0}=cB_{0}$ are reference electric/magnetic fields respectively and the correspondence requires $H\to\varepsilon_{0}$ and $gH\to 1/\mu_{0}\varepsilon$. Indeed, one can show using this correspondence that Maxwell’s equations in the vacuum are equivalent to the shallow water system of equations in an irrotational, quiescent ($\mathbf{V}=0$) fluid. Hence, in standing water, the mass flux is analogous to the electromagnetic power flux.