Towards analogue black hole merger

We study the effects of the wavevector-dependent losses on polariton condensates. We demonstrate that because of these losses, a single vortex becomes a center of a convergent flow, which allows describing it by an analogue Kerr black hole metric with a dynamically evolving origin. For a pair of vortices, we find an analogue of the 3rd Kepler's law and estimate the emission rate of the gravitational waves. We simulate an analogue of the inspiral phase of a black hole merger. Our work therefore suggests that polariton condensates with quantum vortices represent a setting with a fully self-consistent dynamical metric for broad analogue studies.


Introduction
Analogue physics is based on the similarities between the mathematical models describing different systems, where the solutions found in one field can be used in the other fields.This cross-inspiration has already helped to discover the Higgs field and its excitation, the Higgs boson, by analogy with symmetry breaking in superconductors [1].It has also helped to explain the extraordinary properties of graphene by Klein tunneling [2][3][4], an effect predicted in highenergy physics [5], but never observed in its original form.There is currently an important activity in a particular branch of analogue physics: the analogue gravity [6,7].The topic of analogue black hole simulations was started by the seminal work of Unruh [8], that focused on an analogue of the Hawking emission [9], originally expected to occur at the astrophysical black hole horizons, but then predicted by analogy and observed experimentally [10,11] at subsonic-supersonic transitions in classical and quantum fluids.
The first studies on the analogue gravity were considering conservative fluids, both classical [8] and quantum [12].The presence of a steady fluid flow therefore required the presence of a drain (and a source).In 1D, the fluid is usually injected at one edge of the system and removed at the other edge [10,13].In 2D, the so-called "draining bathtub" configuration is used, with a local drain in the center [14,15].Quantum fluids, such as atomic BECs (Bose-Einstein condensates), were rather considered in 1D, without losses, in a dynamical configuration with a moving potential [11,16].The direct effect of the losses on the Hawking emission in analogue systems started to be considered more recently [17].
Exciton-polariton (polariton) condensates [18], formed from strongly coupled excitons and photons in microcavities [19] under sufficiently strong pumping, are a well-known example of a lossy interacting quantum fluid [20] showing superfluidity [21], topological defects [22], and allowing to generate a wide variety of flows.Polariton lifetime usually ranges from 0.2 to several hundreds of picoseconds [23][24][25].This lifetime may seem shorter than in other systems (e.g.evaporation times of cold atoms or liquid helium), but polaritons also exhibit much faster dynamics, both in their propagation (thanks to their small effective mass, m pol ∼ 5×10 −5 m 0 with m 0 being the free electron mass) and in their energy relaxation (thanks to efficient polaritonpolariton and polariton-phonon interactions), the latter allowing, in some cases, to achieve thermal equilibrium in spite of the short lifetime [23,26].In analogue gravity studies, these losses have mostly been considered as being wavevector-independent: they are usually described mathematically by terms like −i Γ(x, y)ψ(x, y, t )/2, where Γ(x, y) is the polariton decay rate (with a possible spatial profile, allowing to describe the draining bathtub).If a polariton flow is superfluid at some point, these losses necessarily lead to the formation of a horizon [27,28].This enabled the demonstration of the formation of a horizon in a 1D polariton wire experimentally [29].The possibilities to measure the dispersion of the excitations of the polariton condensates in order to extract the speed of sound, as well as to obtain the flow velocity from the phase of the condensate wavefunction, have inspired a lot of experimental activity [30][31][32][33].2D polariton superfluids may also enable the study of analogue gravity through the prism of Andreev-Hawking phenomena [34,35].Particularly important for the present work are the studies of vortex dynamics, which have been already carried out in polariton condensates [36][37][38][39].
So far, whatever the platform, in most of the studies of analogue gravity the metric was static, determined by the conditions of the experiment and unable to react dynamically [6].For example, in the studies of superradiance [40,41], the energy of the black hole was decreasing due to the accumulation of negative-energy Bogoliubov excitations, but the metric itself was not modified by this accumulation because that would represent a higher-order correction.The examples of dynamical effects include an analogue black hole ringdown, that is, a decay of a weak perturbation of the horizon, that has been demonstrated recently [42].In particular, in the studies of backreaction [43][44][45][46], all changes to the metric remain weak.In general, most of the settings allow simulating a "large" black hole and a weak excitation propagating in its curved spacetime, but not a strong perturbation or a modification of this spacetime.Analogue Penrose effect with quantum vortices [47] has allowed to strongly modify the angular momentum of a black hole and therefore its metric, but only partially because the black hole's attraction remained the same.Meanwhile, the ultimate goal of achieving emergent analogue gravity [48][49][50] obviously requires a configuration with a metric capable of responding dynamically and evolving with time.
The most well-known dynamical process already observed experimentally with real black holes is the so-called black hole merger [51].It involves two black holes on a close orbit, spiraling in due to the emission of gravitational waves [52,53].The inspiral process is followed by the merger of two horizons into a single one and then by a ringdown of this horizon to its final shape [54,55].The whole process requires a dynamical spacetime for its description: the motion of a black hole means that the overall spacetime changes; the merger and the ringdown of the horizons imply a significant transformation of the shape of each of the two horizons, with released energy which starts to become comparable with that of a single black hole [56].
In this work, we demonstrate that the wavevector-dependent losses inherent to polariton condensates make a single quantum vortex in such condensate a center of a convergent flow, and thus make it an analogue black hole (more precisely, a Kerr black hole, because of the angular momentum of the quantum vortex).Capable of motion, these quantum vortices thus represent a model system for studying the dynamics of a system of gravitating bodies (analogue black holes).We numerically consider a pair of same-sign vortices and describe their motion with an analogue "post-Newtonian" formalism.We derive modified analogues of the Newton's law of universal gravitation and of the 3rd Kepler's law, and find the solutions of the equations of motion for the inspiral stage.We compare these results with the well-known expressions for astrophysical black holes.The validity of our analytical solutions is confirmed by numerical simulations.

The model
We will be dealing with polariton condensation under non-resonant pumping [18,26].As said in the introduction, the polariton losses due to the finite quality factor of the cavity and to the finite exciton lifetime have usually been described as being wavevector-independent.This is, however, just an approximation.The polariton decay rate reads: where |x k | 2 is the exciton fraction of the polariton, Γ x the exciton lifetime and Γ c the cavity photon lifetime, which are both assumed to be constant, which again is an approximation [57][58][59].
Near zero exciton-photon detuning, and for small k, the exciton fraction dependence on the wave vector can be approximated by: where ħΩ is the Rabi splitting and m c is the photon mass in the cavity.As discussed in [60], this gives: Depending on the sign of Γ x − Γ c , a convergent flow can be realised using different methods.If Γ x − Γ c > 0, the decay increases with k.In other words, the decay increases when getting closer to the core of the vortex, which provokes the convergent flow that can lead to the formation of an event horizon and to an attractive interaction between vortices.This condition can be realized in high-Q cavities or structures with fast excitonic decay.
A more standard situation corresponds to Γ x − Γ c < 0. In this case, a convergent flow can be realised by considering non-resonant pumping leading to the formation of a Bose-Einstein condensate.For the condensate and for its excited states, the total net decay or growth rate η k includes the rates of scattering in W i n,k and out W out ,k of the modes: As analysed in [60], W out ,k − W i n,k also scales as k 2 and the net scattering rate η k reads: where m pol is the polariton mass.The presence of a stationary condensate in the ground state means a zero net rate for the corresponding state: η 0 = 0.It also implies the absence of gain for higher-energy states, that is Λ < 0, which is precisely equivalent to having k 2 -losses.Such dependence of the losses on the wave vector has been first introduced phenomenologically as the only way to reproduce experimental results [61,62], providing a very good agreement with them.
Losses depending on the total particle energy were previously suggested to describe relaxation in superfluids [63] and Bose-Einstein condensates [64,65].However, the works on polaritons and polariton condensates have introduced a different shape of the loss term, depending explicitly on the kinetic energy only (k 2 ).Recent experiments have confirmed this k 2 -dependence of the polariton decay [66].
We therefore study the behavior of the polariton condensate described by the Gross-Pitaevskii equation with k 2 -losses and saturated gain: where m is the polariton mass, U is the confining potential (for example, a cylindrical mesa), α is the polariton-polariton interaction constant, Λ is the constant characterizing the k 2 -losses, and γ is the prefactor of the gain term, saturated by the total density n t ot .We do not describe local reservoir saturation [67][68][69], because we consider the case when the condensate dynamics is much faster than that of the reservoir.
Since the losses depend on the wave vector, and the largest wave vectors are present in the core region of the quantum vortex solution for a conservative condensate, it is natural to expect that the effect of the losses appearing in the modified equation ( 6) will also be concentrated in the core region.The relation between the wave vector variation of the excitonic and photonic fractions and the spatial variation of the properties for a quantum vortex was already understood in Ref. [70] (with consequences for the spatial density distribution) and used very recently in the description of the vortex kinematics via the Berry curvature in a recent work [71], but the excitonic core in that work had a longer lifetime, and so the vortices were repelling each other.
In the following sections, we solve Eq. ( 6) numerically to find the properties of a single vortex.We then look for an analytical description of the flows created by this vortex, taking the solution for a quantum vortex in the conservative case at large r as a starting point and treating the loss term as a perturbation.We then study the behavior of a pair of vortices both numerically and analytically.

Single vortex
The metric of the weak excitations of vortex flow with non zero radial velocity v r is formally similar to the one of a Kerr Black hole [47,[72][73][74].Indeed, a cylindrically-symmetric configuration with radial and azimuthal flows (v φ = ±ħ/(mr )), with an appropriate change of coordinates [75], allows writing the metric of the condensate as: with the local speed of sound c s = α|ψ| 2 /m, and v tot the total velocity.The event horizon is given by the condition v r = c s , whereas the static limit occurs when v t ot = c s .In a conservative case, v r = 0, so there is a static limit, but no event horizon.We first perform a numerical solution of the modified Gross-Pitaevskii equation ( 6) to find the stationary configuration with a single vortex and analyze its properties.We use the parameters typical for a GaAs cavity: m = 5 × 10 −5 m 0 , α|ψ| 2 ≈ 1 meV (for a stationary homogeneous solution) giving the healing length of the order of 1 µm, and the k 2 -decay coefficient Λ = 0.1.
Figure 1 shows v t ot and v r .The numerical solution is plotted with solid lines: speed of sound c s (black), radial velocity v r (red), and total velocity v t ot (blue).One clearly see the presence of an horizon and of a static limit.Additional numerical simulations (not shown) indicate that the radius of the horizon, as well as the ratio of the horizon to the static limit, increase with Λ.In practice, Λ can be controlled [60] by the temperature (independently from the polariton mass) or by the detuning (which also affects the polariton mass and the interaction constant).
We now turn to an analytical description of the system under study.Analytical solution for the whole vortex wavefunction is, unfortunately, unknown even in the conservative case.The description of the core in the non-conservative case turns out to be more complicated, so we focus on the opposite limit.The results are still useful, because the long distance interaction between vortices is provided by velocity profiles far from the core.
We therefore consider the limiting case r ≫ ξ.At large distances, by applying series expansion to the equation for the radial part of the wavefunction [76], one can obtain an approximate solution ψ ≈ n 1 − ξ 2 /r 2 e i θ .We then insert this solution into Eq.( 6).The damping term gives two contributions, proportional to the Laplacian: The first contribution is due to the spatial variation of the density profile, while the second contribution is due to the centrifugal term, which is dominant for r ≫ ξ.We obtain therefore: where r is measured from the vortex center.We therefore find that the losses exhibit a spatial pattern determined by the position of the vortex: if the vortex moves, the loss pattern moves as well.This is like having a draining bathtub configuration with a mobile drain accompanying each vortex.
To find analytically the expression for the radial velocity of the flow for a single vortex with k 2losses, we use the continuity equation: since the solution is stationary, all losses inside a certain circle with a radius R are compensated by an inward flow through this circle.The total losses can be found as Here, again, we do not consider the losses taking place inside the vortex core (r < ξ), where our approximation is not applicable.The compensation of these losses by the flow allows to find an asymptotic expression for the radial velocity, again valid for r ≫ ξ: which gives As we will show below, in the next section, the strength of the attraction between vortices is determined by the radial velocity itself proportional to the damping coefficient Λ.If Λ is sufficiently strong, it can overcome the well-known long-range repulsion [76] of same-sign vortices, allowing to observe their bound states.Figure 1 compares the numerical solution of the equation ( 6) with the analytical estimates that we plot only in its validity range (for r > 2ξ).In this limit, the numerical solution is in a good agreement with the behavior suggested by the analytics: it presents the same scaling log r /r for the radial velocity for r ≫ ξ.The position of the static limit is predicted correctly.However, the analytical formula is not valid in the region r ≈ ξ where the event horizon is present and does not give any insight about its existence.

Two bodies: zero angular momentum
In the Newtonian picture, the material points cannot have their proper rotation.The only angular momentum possible is their angular momentum with respect to the center of mass.If it is zero (zero impact parameter), then the Earth falls directly on the Sun following a straight line.For a non-zero angular momentum, the bodies follow elliptical trajectories (in the bound configuration).
In relativity, each of the bodies can have its own angular momentum which cannot be neglected.The frame dragging associated with it leads to the Lense-Thirring effect [77]: the spacetime is "dragged" in the direction of rotation of a rotating body, and everything inside this spacetime is also "dragged" in the same direction.The straight infalling trajectories for zero "external" angular momentum thus become curved [78] because of the "internal" momentum of the attracting body: the Earth does not follow a straight line anymore, but still falls on the Sun for zero impact parameter.
In the most simple case, each of the vortices is moving with the flow surrounding its core, which is defined by the velocity pattern created by the other vortex.In this case, the "external" angular momentum is zero, and all the rotation that can be observed is due to the frame-dragging ("internal" angular momentum of each vortex).In this case, the radial velocity given by Eq. ( 12) completely determines the radial motion of each of them via the following differential equation: where a is an adjustable parameter (a proportionality constant, which includes Λ from Eq. ( 12)).The solution of this differential equation can be explicitly written for t (r ): This analytical solution is compared with numerical simulations for a pair of vortices in Fig. 2. The black line shows the distance between the centers of the two vortices from numerical simulations, and the red dashed line shows a fit based on the analytical model ( 14).This configuration corresponds to a collision of two black holes with the same proper rotation but with zero external angular momentum (zero impact parameter).No emission of gravitational waves is required for the vortices to attract at large distances because there is no angular momentum and associated rotational energy to dissipate.However, mutual rotation of the black holes can still be observed due to frame dragging, that is, the rotation of the spacetime itself in the direction of the rotation of each black hole.It is indeed observed in numerical simulations.

Non-zero angular momentum: Kepler's law
We now turn to the case of non-zero external angular momentum.In this case, in the framework of Newtonian celestial mechanics, if the emission of gravitational waves is neglected, each of the two bodies is expected to follow an elliptic orbit forever, with the period of rotation given by the 3rd Kepler's law.To be able to describe the effect of the gravitational waves in our analogue system, we need to find an equivalent of this law first.For r ≫ ξ, log r /ξ is a slowly varying function with respect to r , and therefore we can assume that α = log r /ξ is approximately constant for a single rotation period.We can thus consider a simplified differential equation for the radial coordinate whose solution is r (t ) = −2at .This allows us to find the radial acceleration which allows us to conclude that the attractive force (which is responsible for this acceleration) behaves as and the corresponding potential energy is This potential energy is decreasing faster than the Newton's law of gravity in 3D.Then, using the fact that the centripetal acceleration for stationary circular motion (considering that the orbit is approximately circular) should be d 2 r /d t 2 = v 2 θ /r , we can find the azimuthal velocity which gives the rotation period This is not a simple power law but has two asymptotics.For r ≫ ξ, T ∼ r 2 , but for r ∼ ξ the logarithm cannot be neglected and T ∼ r , so the power law changes between r 1 and r 2 .This  differs from the 3rd Kepler's law: T ∼ r 3/2 , but all three power law exponents are of the same order.
These analytical results are compared with the results of numerical simulations for a pair of vortices in Fig. 3: the blue solid line shows the modified Kepler's law, while the black dots are the results of numerical simulations.The numerical results are quite well described by the approximate analytical solution.

Gravitational waves
We now study the energy loss suffered by an accelerated vortex via the emission of waves.The nature of these waves is dual: the density waves in the condensate can at the same time be considered as gravitational waves (because they modify the metric of the condensate, which is density-dependent), but also as electromagnetic waves (because the vortices behave as charges [73] interacting via an electromagnetic-like field).Since the dipole mechanism of wave emission is suppressed by the momentum conservation law, the strongest emission corresponds to the quadrupole moment.The expression for the intensity of the quadrupole emission is universal, and does not depend on the nature of the waves (electromagnetic vector waves or gravitational tensor waves) [79] : where the quadrupole moment of a pair of bodies on a circular orbit is (other components are similar), and its 3rd-order time derivative therefore behaves as ...
where ω is the rotation frequency.
Using the centripetal acceleration determined by the azimuthal velocity (20), we obtain Combining the energy loss (25) with the potential energy (19), we find the following differential equation where we neglect the time derivative of α because of the slow variation of log r /ξ, as before.This allows us to simplify and rewrite the equation as: The solution for the trajectory can finally be written as an inverse function t (r ): For large distances, the variation of the log r terms can be neglected, and the solution scales as r (t ) ∼ (t 0 − t ) 1/6 : it is a power law behavior, as for the inspiral of the astrophysical black holes [51], but with a different power (1/6 vs 1/4).For shorter distances, the logarithms and the exponential integral functions start to play an important role, and the resulting trend becomes qualitatively different: at large t , r − ξ ∼ 1/t .This is a consequence of the analytical approximation, which is not valid for r ∼ ξ.
This solution is compared with the results of a numerical experiment in Fig. 4. The analytical solution is plotted as a solid line, while the experimental points are shown as black dots, with a good agreement between the two.The average decrease of the distance is much slower than in the case with zero external angular momentum shown in Fig. 2. It also seems that the vortices cannot approach each other closer than a certain distance, due to the modification of the spatial profiles of their cores, as compared to the case of a single vortex seen in Fig. 1.This leads to repulsion at short distances, which overcomes the attraction.The short range interaction and the means to overcome the repulsion will be studied in future works.
Finally, in Fig. 5 we show the snapshots of the polariton density distributions at two different moments of time corresponding to early (panel a) and late (panel b) stages of the inspiral.The horizon and the static limit are shown with black lines.So far, our simulations did not allow us to reach the merger of the horizons, but we hope that it could be possible under different conditions (stronger losses or a different k-dependence of these losses).On the contrary, a merger of the static limits does occur already in the simulation shown in Fig. 5.
We note that from symmetry considerations, the velocity in the center of the system (between the two vortices) is necessarily zero, and thus this point cannot be "inside" an analogue Kerr black hole, whose interior is defined by the condition v r > c.On the other hand, real black holes, as for example in the case of Kerr-Newman metric, can have a complicated internal structure with several horizons [80][81][82], and similar effects occur in analogue black holes [11,83] (so far only in 1D).We can therefore conclude that if a merger is achieved, the resulting configuration could have two concentric horizons.

Conclusions
To conclude, we have demonstrated that the presence of a k 2 -loss term in quantum fluid creates a convergent current for quantum vortices.This type of term is naturally present in excitonpolariton quantum fluid, either freely decaying, if the exciton loss exceeds the photon loss, or non-resonantly pumped, when an equilibrium Bose-Einstein condensate forms.The metric associated with quantum vortex is similar to that of a Kerr black hole, showing both a horizon and a static limit.We show that the convergent "draining bathtub" flow attached to each vortex provokes an attraction between them and leads to the emergence of an analogue of attractive gravitational interaction for vortices.For a pair of vortices, we demonstrate how the emission of analogue gravitational waves leads to the inspiral of these analogue black holes, with increasing rotation frequency.

Figure 1 .
Figure 1.The speed of sound c s , radial velocity v r and total velocity v t ot as functions of r (normalized by the sound velocity at infinity c i n f and the healing length ξ).The positions of the horizon and the static limit (defined by v r = c s and v t ot = c s , respectively) are marked by vertical dashed lines.

Figure 2 .
Figure 2. Zero angular momentum merger: time evolution of the distance between the analogue black holes.Numerical simulation (black line) and analytical solution (14) (red dashed line).

Figure 3 .
Figure 3. Rotation period as a function of distance between the analogue black holes: black dots -numerical simulations, blue solid line -analytical solution (21).

Figure 4 .
Figure 4. Non-zero angular momentum merger: time evolution of the distance between the analogue black holes.Black dots -numerical simulation, red solid line -analytical solution (28).

Figure 5 .
Figure 5. Two snapshots of the analogue black holes at the early (a) and late (b) stages of the inspiral.False color shows the polariton density |ψ| 2 , black lines mark the position of the static limit and the horizon.