No hair theorems for analogue black holes
Abstract
We show that transonic onedimensional flows which are analogous to black holes obey nohair theorems both at the level of linear perturbations and in nonlinear regimes. Considering solutions of the GrossPitaevskii (or Kortewegde Vries) equation, we show that stationary flows which are asymptotically uniform on both sides of the horizon are stable and act as attractors. Using Whitham’s modulation theory, we analytically characterize the emitted waves when starting from uniform perturbations. Numerical simulations confirm the validity of this approximation and extend the results to more general perturbations and to the (nonintegrable) cubicquintic GrossPitaevskii equation. When considering time reversed flows that correspond to white holes, the asymptotically uniform flows are unstable to sufficiently large perturbations and emit either a macroscopic undulation in the supersonic side, or a nonlinear superposition of soliton trains.
pacs:
I introduction
On a superficial level there are similarities between general relativity and hydrodynamics as they both rest on non linear partial differential equations depending on space and time. Interestingly, under some specific circumstances, this mere analogy is replaced by a precise correspondence. For instance, when working in the hydrodynamical approximation, i.e., for long wave lengths, the linearized wave equation in a moving fluid is identical to the d’Alembert equation of a scalar field propagating in a certain four dimensional curved spacetime Unruh (1981). This remark has been the starting point of investigations aimed at understanding the validity domain of this correspondence. It is clear that the strict equivalence is lost when including dispersive effects which affect short wave length modes propagating in media Jacobson (1991). However an interesting and non trivial question is to identify the consequences of these dispersive effects: for instance, how would they affect the asymptotic properties of the Hawking radiation emitted by an analogue black hole flow Unruh (1995)? After a thorough analysis, it was found that they do not significantly affect the spectrum when certain inequalities involving the short dispersive length and the spatial gradient of the flow are met, see Refs. Coutant et al. (2012); Robertson (2012) for recent updates.
The present work aims to address similar questions when including nonlinear effects. It is clear that the equations of general relativity and hydrodynamics differ at this level. However, as we explicitly show, both the linear and nonlinear stability of analogue black hole flows are in direct correspondence with those of black holes in general relativity Misner et al. (1973). Let us briefly remind the main properties of the latter. First, stationary black holes are characterized by very few parameters, namely their mass, angular momentum, and electric charge. Second, when perturbed during a finite time interval, the time evolution subsequently brings back the solution to one of those stationary solutions. These two properties are generally referred to as “no hair theorems”. The precise validity limits of these theorems is still the subject of investigations Chrusciel et al. (2012). Moreover, several hairy black hole solutions have been found when admitting exotic matter fields, see Refs. in Chrusciel et al. (2012).
When studying the behavior of onedimensional transonic flows, solutions of the GrossPitaevskii (GP) or Kortewegde Vries (KdV) equation, we observe similar properties. We work with timeindependent external potentials (shape of the obstacle when using the KdV equation) which vary only in a finite domain of . In this case, when considering the set of stationary flows which are asymptotically homogeneous (AH) as – a condition analogous to asymptotic flatness for black holes – a simple counting of integration constants reveals that there is at most a discrete set of solutions. Moreover, when the potential consists of a single step, we explicitly demonstrate that the solution is unique. In other words, one finds a single series of solutions which can be parameterized by the conserved current characterizing the flow. When the potential is smooth, we numerically found a series of solutions which smoothly connect to this series in the step like limit. In certain cases, we also found a disconnected series of solutions. These can be considered as hairy black holes as they contain a large fraction of a soliton attached to the sonic horizon.
The analogy with general relativity is significantly reinforced when considering the stability of the AH black hole solutions. At the linear level, numerical and analytical results establish that local perturbations decay in time so as to reach an AH stationary solution. This conclusion applies both to the GP and the KdV equations. Hence, at the linear level, the AH black hole solutions of these equations are attractors. It has to be pointed out that the time reversed flows, analogous to white holes, are not attractors: in these flows, the scattering of linear perturbations generates a macroscopic undulation Coutant and Parentani (2014). When adding non linearities, several behaviors are found at late time. Our findings are in agreement with, and add to, the results of Refs. Mayoral et al. (2011); Michel and Parentani (2013, 2015); de Nova et al. (2015).
When taking into account the nonlinearities of the GP and KdV equations, by a combination of analytical and numerical methods, we show that there is a large domain of initial black hole configurations which evolve towards an AH stationary solution. In particular, using Whitham modulation theory applied to the GP equation Kamchatnov (2000), we show that the “hair” present in the initial configuration is expelled away from the horizon by three nonlinear waves. At at small amplitudes, these three waves give back the three types of linear waves emitted by an analogue black hole flow. As a result, this emission can be considered as the nonlinear version of the stimulated Hawking effect Rousseaux et al. (2008); Macher and Parentani (2009a); Weinfurtner et al. (2011). The predictions of Whitham’s theory are confirmed and extended by numerical analysis. Interestingly it seems that the integrability of the equations plays an important role both in analogue models and in general relativity. In the first case, it provides the invariants needed to apply Whitham’s theory and eventually characterize the timeevolution from a perturbed configuration to an AH, stable stationary solution. In the second case, an integrable nonlinear sigma model is used to constrain the set of asymptotically flat, stationary and axisymmetric solutions of EinsteinMaxwell theory, leading to the uniqueness theorem Chrusciel et al. (2012).
The paper is organized as follows. In Section II, we first study the stationary solutions of the GP equation and show that the set of AH solutions is discrete. Then we study the linear stability of black hole flows and show that perturbations all decay polynomially in time. In Section III, we consider the nonlinear stability of black hole flows. We first use Whitham’s theory in subsection III.1, and then present numerical results for the full GP equation in subsection III.2. Some new results on whitehole flows are given in Section IV. We conclude in Section V. In appendixes we present various technical aspects concerning notions used in the main body of the text. Appendix A reminds the derivation of Whitham’s equations and derives the properties of the nonlinear waves used in the main text. Appendixes B and C extend the analysis to KdVlike equations, while Appendix D treats the case of a cubicquintic GP equation.
Ii Uniqueness and linear stability of blackhole flows
ii.1 Asymptotically homogeneous transonic flows
We consider a onedimensional flowing condensate whose wavefunction satisfies the GP equation. In a unit system where the atomic mass and reduced Planck constant are equal to , this equation reads
(1) 
The static external potential is denoted by and the effective onedimensional coupling by . In order to have stable homogeneous configurations for constant and , we shall consider only condensates with repulsive interactions between atoms, so that . The stationary solutions can be written as
(2) 
where , , and are real. is the density of the condensate and its local velocity. Plugging Eq. (2) into Eq. (1) gives
(3) 
where is the conserved current.
ii.1.1 Homogeneous potentials
To prepare the analysis of transonic flows associated with spacedependent and , it is useful to recall the main properties of the stationary solutions when and are independent of . In this case, spatially bounded solutions of Eq. (3) exist if and only if
(4) 
All bounded solutions are periodic and can be written using elliptic functions Kamchatnov (2000). Two of them are homogeneous, and given by the two positive roots of the righthand side of Eq. (3). They merge when . The smallest one, , describes a supersonic flow since the condensate velocity is larger than the speed of longwavelength perturbations . Instead, the flow described by is subsonic as . Interestingly, the stationary perturbations on top of these two homogeneous solutions are radically different. Stationary flows with a density close to the supersonic value contain a zero frequency modulation of the density (hereafter called an undulation). Instead, flows with mean densities close to contain a train of solitons, see Fig. 1.
Let us briefly explain how to get the main properties of these solutions. More details can be found in Michel and Parentani (2013). Multiplying Eq. (3) by and integrating the resulting equation gives
(5) 
where , are three constants. They are related to the current and frequency through
(6) 
When there is no double root, the range of is a connected component of the domain where the polynomial on the righthand side of Eq. (5) is positive. The density is therefore asymptotically bounded only if that domain is itself bounded. This requires that the three constants are all real. We order them as . Then, the two domains in which the righthand side of Eq. (5) is positive are and . The second one corresponds to a divergent solution. The only bounded one thus corresponds to the first domain, and oscillates between and . The bounded solution is also characterized by a fourth parameter giving the phase of the oscillations.
In the limit , becomes the homogeneous supersonic solution with density discussed after Eq. (4). When is close to , describes a smallamplitude undulation on top of this homogeneous solution, see Fig. 1. In the limit , one gets a homogeneous subsonic flow of density , and when is close to , contains a train of widelyseparated solitons.
Interestingly, when there exists another solution, called “shadow soliton”, which diverges at a finite value of and also goes to the subsonic flow asymptotically. This last solution must be discarded when working with homogeneous functions and because of its diverging character. However, it must be included in the forthcoming analysis of the steplike case as the divergence can be erased by the change of and .
ii.1.2 Steplike potentials, unicity theorem
To obtain stationary transonic flows, and/or must depend on . To have simple solutions, we assume that and are given by
(7)  
with and . Then, for any real value of , there exists a globally homogeneous solution of Eq. (3) given by
(8) 
Its frequency increases with so that Eq. (4) is always satisfied on each side. Explicitly, one finds
(9) 
This solution is transonic iff
(10) 
The flow is then subsonic for and supersonic for . When , it corresponds to a black hole flow. This can be deduced from the large reduction of the wave vector experienced by counterpropagating waves, see next subsection and Unruh (1981); Macher and Parentani (2009a); Barcelo et al. (2011). When and are negative, the transonic flow is analogous to a white hole. In that case, long wave length incoming modes are converted into short wave length modes.
Interestingly, having fixed , , and of Eq. (7), when obeys Eq. (10), the branch of solutions of Eq. (9) gives the unique AH transonic flow. Indeed, for this interval of , all the other stationary bounded transonic solutions contain an undulation in the supersonic region and/or a train of solitons in the subsonic one, see Fig. 1. The uniqueness can be understood when considering the set of integration constants which characterize the general stationary solution. They are all fixed by the requirement that the solution be AH on both sides.
Yet, there is another series of stationary AH transonic flows. They contain a halfsoliton for , and correspond to the “waterfall” configuration studied in Larré et al. (2012). The corresponding value of is given by
(11) 
where is the density for . Quite remarkably, these solutions are bounded only for values of larger than the upper bound of Eq. (10). When decreasing along the branch of waterfall solutions, the difference between the densities on the two sides of the horizon decreases. It vanishes precisely when , at which point the homogeneous and the waterfall solutions coincide. Decreasing further, only the homogeneous solution of Eq. (8) remains bounded. The series of Eq. (9) and Eq. (11) are represented in Fig. 2.
In brief, in the case of step like potentials, for all , there is a unique AH transonic solution. Hence in this case, stationary transonic flows obey an analogous version of the black hole uniqueness theorem Israel (1967).
These results may be understood as follows. The analysis of subsection II.1.1 applies on both sides of the horizon . We thus have six integration constants . Notice that these are not all independent when considering global solutions, since the current and frequency take the same values on both sides. A global solution is in fact determined by 4 parameters, including and . Given these two quantities, the space of solutions is thus twodimensional and parametrized by two integration constants. These fix the phase and amplitude of the undulation on the supersonic side, and the properties of the soliton train on the other side, as explained in Michel and Parentani (2015).
ii.1.3 Smooth potentials and hairy black holes
When and smoothly vary with , the above arguments are no longer sufficient to characterize the complete set of stationary bounded solutions of Eq. (1). Yet, the set of AH stationary transonic flows remains discrete at fixed . Indeed, imposing that the flow is AH and supersonic for completely specifies a solution of Eq. (3) once and are fixed. Then the condition that the solution be also AH for gives an extra condition, which constrains the frequency to belong to a (possibly empty) discrete set of branches . (Notice that the homogeneous solution of Eq. (8) still exists if the variations of and are tuned in such a way that be a constant.)
To generalize the step like potentials of Eq. (7), we considered the continuous profiles
(12)  
characterized by the domains and where and linearly interpolate between their asymptotic values. To obtain the stationary flows, we numerically solved Eq. (5). As expected, we found AH solutions which are smoothly connected to those of the former subsection in the limit . More interestingly, we also found a series of “hairy” AH solutions when is smaller than . These solutions are disconnected from those of the main series because they contain an almost complete soliton which is attached to the sonic horizon, and which mainly lives in the subsonic region, see Fig. 3. Our numerical study indicates that its center is rejected at in the steplike limit, or in the limit where . Our simulations also indicate that this solution is unstable under nonlinear perturbations, which may trigger the emission of the soliton at infinity. We hope to characterize this instability in a future work. On the right panel of Fig. 2, we show the values of for the AH solutions obtained with Eq. (12). The red dots describe solutions of the main series with . The green dots (crosses) correspond to hairy solutions obtained with . It is clear that these solutions belong to two different series of . In Fig. 3 we represent the density profile of three solutions of the main series and one hairy solution with .
In conclusion, when dealing with smooth profiles for and , one always finds the main series of AH solutions characterized by a smooth density profile across the horizon. In addition, a discrete number of hairy solutions can exist. When we found such solutions, we observed that they are disconnected from the solutions of the main series because they contain some solitonic configuration attached to the region where the potential varies.
ii.2 Linear stability of black hole flows
When considering the stability of transonic flows, it appears that those corresponding to black holes are much more stable than those describing white holes. This is non trivial because the Smatrix describing the scattering of linear perturbations on a white hole horizon is the inverse of that describing the scattering on the time reversed black hole flow Macher and Parentani (2009a). This relation implies that white hole flows possess the same degree of stability as black holes, namely, they are dynamically stable. (This means that the spectrum of linear perturbations contains no complex frequency modes.) Yet, the late time evolution of perturbations scattered on a sonic horizon strongly depends on whether the modes experience a redshift (black hole), or a blueshift (white hole). White hole flows are studied in Section IV where it is shown that they exhibit a linear infrared instability. In the present subsection we consider the linear stability properties of blackhole flows and show that they are stable under linear perturbations.
To study the behavior of perturbations in qualitative terms, we present the results of some numerical simulations. We choose a spatial domain containing the sonic horizon and look at the latetime evolution of the perturbations in . Numerical simulations indicate that perturbations decay as up to logarithmic factors. These observations hold for all AH solutions, both in a steplike potential of Eq. (8) and in smooth and . To represent the various typical cases, we consider the scattering of a gaussian perturbation which is initially localized in the subsonic region, on the sonic horizon, or in the supersonic region. On the left panel of Fig. 4, we work in the linear regime and study the evolution of a perturbation of relative amplitude of order and of width equal to (in units of the healing length). On the right panel, to show that similar results hold for much larger perturbations, we work with of order for the same spatial width. Then, we show the time dependence of the integral of over the segment , which contains the sonic horizon localized on . At early times, the density fluctuation increases as the perturbation enters the integration region, or remains roughly constant when the perturbation is already in the integration region at . At late times, in the three cases, one clearly sees that the average of goes to zero as . The fact that the temporal behavior is very similar on the left panel (linear regime) and on the right panel (nonlinear regime since ), reveals that the decay law of in is extremely robust. It implies that all perturbations are diluted away so as to give at late time an AH stationary solution.
This behavior is similar to that of linear perturbations propagating in a black hole metric Misner et al. (1973). We briefly present at the end of the Section the nature of the correspondence between the present case and relativistic fields. We now give the analytical elements needed to understand this behavior. The full calculation is done in Appendix C for the linearized KdV equation and a square perturbation.
We look for perturbed solutions of the form
(13) 
where , see Eq. (2). To first order in , Eq. (1) gives
(14) 
Asymptotically, on either side, the background quantities and are constant. Hence we can look for solutions of the form
(15) 
Plugging this form into Eq. (14), we obtain a system of two linear equations on and , which has nontrivial solutions if the dispersion relation
(16) 
(shown in Fig. 5) is satisfied. For , there are two roots in in the subsonic region. In the supersonic region, when is smaller than a critical frequency , see Eq. (25) in Macher and Parentani (2009a), there are two additional roots. For , there are thus six asymptotic plane waves (two in the subsonic region and four in the supersonic one), which give the asymptotic behavior of the globally defined bounded modes. These modes span a threedimensional vector space. The scattering process is characterized by the relationship between the incoming and outgoing modes. The in (out) basis is composed of the three modes which asymptotically contain one wave with group velocity oriented towards (away from) the sonic horizon. These two bases are related by a Bogoliubov transformation. The behavior of its coefficients and the link with the Hawking effect can be found in Macher and Parentani (2009a, b).
We consider some initial perturbation and look for its latetime behavior. We assume that decays fast enough at infinity so that the coefficients of its expansion into incoming modes are finite. In each of the asymptotic regions, has the form
(17) 
where the subscript labels the branches of the dispersion relation. The coefficients and are determined by the overlap between and incoming modes, and by the expansion of the latter into asymptotic plane waves. Possible divergences in the expansion into outgoing modes come only from those of the 9 Bogoliubov coefficients. These were computed analytically in the steplike regime in Finazzi and Parentani (2012) and numerically for a smooth flow in Macher and Parentani (2009a). In both cases, it was found that the only divergence occurs for , and only affects the coefficients relating longwavelength to shortwavelength modes. However, as explained in Appendix C, the leading terms in these coefficients do not contribute to the amplitude of as they come in pairs with opposite signs which cancel each other. At late times, can be computed using a saddlepoint approximation. The saddlepoints are located where
(18) 
In this expression, represents the phase of the prefactor ( or ). In the limit at fixed , we must thus look for points where , i.e., where the group velocity vanishes. (If the initial conditions are smooth, divergences in can arise only from terms in in the expression of , where is a smooth function, coming from the integral giving the overlap, analogous to Eq. (105) for the KdV equation. Hence the divergences of can only arise through divergences of .) Such divergence comes only at , for the two roots which merge there. A straightforward calculation shows that close to this point diverges as . Performing a Gaussian integration, we obtain that the dominant waves decrease as . This confirms that the dominant waves at late times have a frequency near , in accordance with the numerical results of Fig. 4. In principle, this linear analysis is valid provided at all times. However, numerical simulations indicate that the latetime behavior remains the same even for “large” initial perturbations, with close to (see Fig. 5, right panel). This is due to the dilution of the perturbation: as explained in the next section, homogeneous black hole solutions are stable, so that initial nonlinear perturbations will generally decay in time. Hence the perturbation close to the horizon effectively enters the linear domain at late times. In brief, numerical simulations and the above discussion establish that the AH black hole flows are linearly stable and that, in the vicinity of the sonic horizon, all linear perturbations decay as . This is very reminiscent to the no hair theorem of general relativity Misner et al. (1973).
To conclude this subsection, we briefly sketch the link between the decay we obtained and that exhibited by free scalar fields propagating in black hole metric. The comparison should be done at two different levels, namely when including or not the dispersive terms of Eq. (16). When ignoring dispersion effects and the quantum pressure, the simplest way to proceed consists in separating into its real and imaginary parts Barcelo et al. (2011). Straightforward algebra shows that the imaginary part of obeys
(19) 
where are , and where Einstein’s summation conventions are used. The matrix is given by
(20) 
where we introduced the local value of the sound speed . Ignoring some subtleties related to the conformal part of the metric in 1+1 dimensions, Eq. (19) is the d’Alembert equation in a stationary black hole with line element given by
(21) 
When the flow is transonic, it describes a black hole metric and the event horizon is located where . The analogy goes deeper as the “analog” surface gravity given by (evaluated on the horizon) governs the redshifting of the solutions of Eq. (19) which propagate against the flow. Namely the wave number of a localized wave packet leaving the near horizon region is redshifted according to Macher and Parentani (2009a), exactly as found in General Relativity. At the classical level, this guarantees that the decay of the solutions of Eq. (19) will be similar to that found in general relativity Misner et al. (1973). Similarly, at the quantum level, the scattering of vacuum fluctuations on the sonic horizon should produce a thermal flux with a temperature given by the standard expression Unruh (1981).
Yet, to reproduce the late time decay law in , the dispersive effects of Eq. (16) must be taken into account, as they are governed by the behavior near the critical frequency .^{1}^{1}1When including them, the separation into real and imaginary parts of is no longer useful as they now couple to each other, see Eq. (A6) in Macher and Parentani (2009a). In fact, this coupling complicates the relationship between the present settings and those characterizing dispersive fields in the context of alternative theories of gravity which break the local Lorentz invariance at short distances, see e.g. Jacobson (2010). The interested reader might consult the appendixes A and B of Macher and Parentani (2009a) for further discussions. To clarify the role of , it is interesting to consider the decay represented in Fig. 4 when sending while keeping fixed all quantities appearing in Eq. (19). In Fig. 6 we represent the decay of the averaged squared density fluctuations for four different values of , namely (blue, continuous), (purple, dashed), (red, dotted), and (orange, dotdashed). To compare the four cases, the time is given in units of . For the same reason, the vertical axis represents the squared relative density perturbation , averaged over a domain which scales as in the unit system used so far, see Eq. (1). We clearly see that decreasing leaves the early evolution unchanged. This means that the decay of linear density perturbations at early times is governed by the relativistic equation (19). At latetime we find the decay law in . We verified that the value of when it starts to decay in is proportional to the squared of the Fourier component of the initial value of the relative density perturbation with the critical wave vector . Hence, even when working at fixed , the evolution of smooth perturbations (i.e., with no Fourier component for ) is governed by Eq. (19). This is guaranteed by the fact that the propagation of the perturbation in the black hole flow redshifts the wave vectors . In brief, Fig. 6 establishes that the dispersive effects of Eq. (16) reduce the decay rate of perturbations near a black hole horizon when compared with the twodimensional relativistic result. Importantly, this new behavior confirms that the hair is lost at late time.
Iii Nonlinear stability of black hole flows
iii.1 Analytical results
We now turn to the nonlinear evolution of initial perturbations on transonic black hole flows. In this subsection, we show analytical results obtained using Whitham’s modulation theory Whitham (1965). Numerical results, which do not rely on the assumptions of this method, are shown afterwards. The reader interested in the derivation of Whitham’s equations for the problem at hands may consult Appendix A and the textbook Kamchatnov (2000).
Whitham’s modulation theory offers a general scheme for finding approximate quasiperiodic solutions of (quasi) integrable nonlinear partial differential equations. The solutions we consider oscillate on a short scale while their amplitudes, mean values and wavelengths vary on much larger ones. The key idea is to average some of the conservation laws of the equation over the fast scale to obtain coupled equations for the slow evolution of a set of effective parameters, called Riemann invariants, which describe the solution locally. This method is particularly useful when the equation is integrable by the inverse scattering method. As reminded in Appendix A, the AKNS scheme then provides a general way to find the set of Riemann invariants. To the best of our knowledge, Eq. (1) with dependent functions and is not integrable. However, since Whitham’s equations are local, when working with and of Eq. (7), one can use the integrability for uniform to determine the solutions on each side of the discontinuity at . The globally defined solutions are then obtained by using the matching conditions, i.e., continuity of and .
We look for the time dependence of solutions characterized by asymptotic densities and velocities as different from those of the solution given by Eq. (8). Our aim is to show that, for a wide range of values of , the solution locally converges to that particular homogeneous solution. Doing so we shall first extend the stability analysis of the previous section by including nonlinear effects, as well as perturbations extending to infinity. Second, we shall obtain the main features of the emission process which progressively replaces the initial configuration by the homogeneous and stationary blackhole flow of Eq. (8). We shall see that, in a large domain of parameter space, three macroscopic, nonlinear, scaleinvariant waves are emitted. These can be conceived as the result of a nonlinear stimulated Hawking radiation. As shown below, these features are well reproduced by numerical simulations and, along with the above linear analysis, they provide a precise description of the latetime behavior of the solution when approaching the homogeneous black hole flow. They confirm that the solution given by Eq. (8) acts as an attractor.
In this work we consider approximate solutions described by at most 4 Riemann invariants. Since the steplike dependence of and introduces no scale, one can further restrict our attention to scaleinvariant solutions, for which the Riemann invariants depend only on . Our goal is to find the domain of parameter space in which the timedependent solution interpolates between given asymptotic values of and as , and a homogeneous blackhole solution in some spatial interval which contains and grows linearly in time. That is, we look for the solutions of the Whitham equations (55) describing functions and such that

and for ,

and for ,

of Eq. (8) in an open interval of containing .
Such global solutions can be built using two types of nonlinear waves found when are constant: dispersive shock waves (DSW) and simple waves (SW), see Appendix A.3. Two examples are shown in Fig. 7. Along a SW, and are monotonic functions of the sole variable . They are related through
(22) 
where the sign is positive for a rightmoving wave (in the fluid frame) and negative for a leftmoving one. On the other hand, a DSW interpolates between smallamplitude oscillations (which vanish at the edge of the wave) and a soliton.
When the variations of the density and velocity along a SW, or a DSW, are small, the wave has a small amplitude and propagates with a velocity (in the lab frame) close to the group velocity of long wavelength linear perturbations, i.e., with , in agreement with Eq. (22). For a SW, this result remains true whatever the amplitude because it can be described locally as a superposition of nondispersive waves. For a DSW instead, dispersion plays an important role when the amplitude becomes large, hence the more complicated expressions for the velocity, see Eqs. (57,58). As a result, in a black hole flow, provided their amplitudes are not too large, these waves propagate away from the horizon. They may thus be seen as a nonlinear version of outgoing wavepackets produced by the scattering of the initial configurations on the horizon. It has to be noticed that the three types of outgoing waves emitted by a black hole flow are governed by long wave length roots of the dispersion relation, see Fig. 5.
Global solutions can be obtained by matching exact solutions of the Whitham equation on each side of , imposing continuity of and . Let us first study the case of symmetric asymptotic conditions, i.e., and . We consider solutions with three waves (DSW or SW) in total, two propagating in and one in . This is motivated by the behavior of linear modes emitted by a black hole. At the linear level, two waves are emitted in the supersonic region and one in the subsonic region Macher and Parentani (2009a). Since the DSW and SW we are looking for are nonlinear versions of outgoing wavepackets, it is natural to assume they follow the same behavior. The validity of this hypothesis will come a posteriori from the existence of solutions when is sufficiently close to the density of Eq. (8).
Using the properties of SW and DSW outlined in Appendix A.3, we find two types of solutions depending on the sign of . If , the solution has one DSW for and two SW for , separated by two homogeneous regions. A schematic plot is shown on the left panel of Fig. 8. On the left of the DSW and on the right of the two SW, the solution is (by construction) homogeneous with density and velocity . Between the DSW and the leftmost SW, the density is precisely equal to of Eq. (8). This equality can be understood from the fact that is the only value of allowing to match two homogeneous solutions at . The conservation law Eq. (59) fixes the final value of the velocity around in terms of the asymptotic conditions on the left side:
(23) 
In general, the latetime value of the current close to the horizon differs from the initial one. The difference is carried by the nonlinear waves.
We now consider the validity domain of this solution. We find that it exists if and only if the two following conditions are satisfied:
(24) 
This domain is shown in blue in Fig. 9 for . (We remind the reader that, for steplike potentials, and only intervene in fixing the value of .) The important point is that for , these inequalities reduce to . This condition is equivalent to saying that the left asymptotic region is subcritical and the right one is supercritical, see Eq. (10). We also notice that the blue domain pinches off when saturates its lowest bound. For lower values of , the flow is globally subsonic. In brief, if the asymptotic conditions are such that the flow is transonic and of the black hole type, then one solution with one DSW and two SW exists provided is not too large.
For , one finds essentially the same results with one SW emitted to the left and two DSW emitted to the right. A schematic plot of the density profile is shown on the right panel of Fig. 8. The density and the velocity around are still given by and of Eq. (23). The conditions of existence of this solution are
(25) 
The argument below Eq. (III.1) can also be applied here, giving that for this solution exists provided is not too large. The corresponding domain is shown in orange in Fig. 9.
The same analysis can be carried out with different asymptotic conditions , . We must then consider 6 additional types of solutions. 4 of them are obtained from the above ones by replacing one of the two DSW by a SW or conversely. One solution has three SW. The last one has three DSW. The domains of existence of each of these solutions are given in Appendix. A.4. There, it is also shown that one of these 8 solutions always exists in a neighbourhood of any asymptotic conditions compatible with a homogeneous blackhole flow, i.e., and .
To summarize, when working with the steplike and of Eq. (7), we have shown that sufficiently small initial perturbations are expelled at infinity, leaving at late time the homogeneous flow . The latter therefore acts as a local attractor, in the sense that the solution uniformly converges to over any bounded interval. In addition, the velocity profile uniformly converges to of Eq. (23). To make this claim more precise, we propose the following
Conjecture.
There exist two sets of positive real numbers and such that for every initial data satisfying the three conditions:

and are homogeneous outside of some bounded interval, with the asymptotic values such that ;

and ;

and for every ;
we have, for every bounded interval ,
We have used here the standard notation for any function and any interval . Although a natural analogue of the asymptotic flatness condition for the initial data, condition 1. may turn out to be too restrictive. We expect that an exponential convergence, or even maybe polynomial ones, should be enough. Similarly, the sufficient level of regularity of the initial data – i.e. the actual value of – is not entirely clear to us at this stage, with and as natural candidates. Conditions 2. and 3. are, despite their dependence on the choice of , the important point, providing a sense in which the initial perturbation is small enough. We hope to be able to sort these questions out in a future work. Let us nonetheless emphasize that all the numerical simulations presented in the next section support the above conjecture. They also indicate that the conjecture should hold when replacing the steplike functions of Eq. (7) by smooth ones, such as those of Eq. (III.2). Finally, under the above three conditions, we conjecture that the latetime properties of the three nonlinear waves moving away from the horizon should only depend on the asymptotic initial conditions and . In other words, the Fourier components of the smooth profiles and are diluted away at very late time, as was rigorously shown in the linear treatment in the former section, and as was also found in our simulations, see below.
To conclude this section, we notice that the present analysis does not apply to whitehole flows. Indeed, a crucial point in our calculations is that the three non linear waves move away from the discontinuity at . For black hole flows, this is realized for the domain of initial conditions represented in Fig. 9. In whitehole flows instead, the nonlinear waves move towards in the supersonic region, see Section IV. When reaching this point, the Whitham theory breaks down. We then expect that whitehole flows will show a more complex behavior than blackhole ones, in accordance with Michel and Parentani (2015). New results concerning whitehole flows are presented in Section IV.
iii.2 Numerical results
We numerically solved the GP equation for several reasons. First, by solving the GP equation directly, the results of the previous section can be checked without relying on the approximations of Whitham’s modulation theory. Second, the differences in the emission process induced by smooth functions and can be studied. Finally, we wanted to study what happens outside the domain of existence of solutions with three waves. In all simulations, we used the code of Michel and Parentani (2015).
In Fig. 10, we show two solutions obtained numerically when starting at from a configuration with homogeneous velocity and density . On the left panel, and satisfy Eq. (III.1), while on the right panel they satisfy Eq. (III.1).
At late times, we observe that these solutions are superpositions of the three waves of Section III.1 plus perturbations whose amplitude decays in time as . We checked that the properties of the SW and DSW, in particular their amplitudes, the position of their edges, and their domains of existence, agree with those given by Whitham’s equations in the limit . We performed additional simulations by replacing the homogeneous initial conditions by smooth varying ones in a bounded domain, and we observed that this agreement was preserved. It would be interesting to identify sufficient conditions on local variations of and ensuring that the late time properties of the three nonlinear waves only depend on the asymptotic values of and . Moreover, the fact that perturbations still decay as shows that the three macroscopic waves do not change the latetime behavior of perturbations obtained in Section II.2 close to the horizon.
To complete our analysis, we replaced steplike and by smooth functions. We first chose their variations so that a solution with homogeneous density still exists, i.e., such that is a constant. We worked with functions of the form
(26) 
where . When the asymptotic conditions are inside the domain described by Eqs. (III.1,III.1), we observed that the latetime properties of the solution are the same as in the steplike case: the homogeneous blackhole solution is reached for through the emission of three nonlinear waves with the same properties. However, the typical formation time of the SW and DSW now depends on . It is linear in for . We verified that these results qualitatively extend to the case where takes different values for and . In that case, there exists no stationary solution with a homogeneous density. However, as discussed in subsection II.1, there is a oneparameter family of solutions with angular frequencies and asymptotically constant values of . Our numerical results indicate that, if the initial conditions are sufficiently close to this series, one of its solutions is reached at through the emission of three waves. In brief, all simulations confirm that initial configurations within the domain specified by Eqs. (III.1,III.1) all evolve, at late times, towards an AH solution.
We finally performed numerical simulations starting from initial conditions with asymptotic behaviors outside the domains described by Eqs. (III.1,III.1). These conditions can be violated in two different ways. In the first case, corresponding to crossing one of the two vertical boundaries in Fig. 9, the two waves in the supersonic region overlap each other, producing a complicate interference pattern. Yet, our simulations indicate that the solution still converges to a homogeneous blackhole flow at late times, as the overlapping waves still escape to infinity. In the second case, one of the waves in the region (respectively ) has its left (respectively, right) edge moving to the left (respectively, right). [The wave with one edge moving towards the horizon lies in the region when crossing the lower boundary, or in the region when crossing the upper one. The sign of then gives its type (DSW or SW).] When the wave is a SW, the latetime solution contains part of a soliton or shadow soliton Michel and Parentani (2013), and asymptotes to an asymptotically homogeneous density different from on the corresponding side.
When the corresponding wave is a DSW, the situation is more complicated. In the simplest case, the solution still becomes stationary at late times over any finite interval of . It then contains a stationary density modulation attached to . We also found cases where the solution apparently never reaches a stationary profile around , instead emitting a density modulation with a nonvanishing phase velocity, which at the nonlinear level corresponds to a soliton train. Our numerical investigation was not extensive enough to determine with confidence the conditions in which one or the other behavior occurs, although it seems that a stationary solution is reached when the DSW with an edge moving towards the horizon is on the left , while the emission of soliton trains occurs when it is on the right .
Iv Whitehole flows
In this section we briefly study the timeevolution of perturbations on whitehole flows. As discussed at the end of Section III.1, there is an important difference between the evolutions of black and whitehole flows: while the former expel perturbations at infinity, the latter have a tendency to accumulate them close to the horizon, as can be understood from the fact that whiteholes behave as the timereversed of black holes. One can thus expect that whitehole flows will act as “repellents” rather than attractors. We here show that it is indeed the case. When starting from the homogeneous solution (which is an attractor for black hole flows), depending on whether the perturbation gives rise to a (sufficiently large) decrease or an increase of the near horizon density, the flow is destabilized by nonlinear effects and develops either only a macroscopic undulation, or sends a train of solitons accompanied by a macroscopic undulation. When working to linear order, one finds that small perturbations generally leave at late times a stationary undulation with a large amplitude, thereby signaling an infrared instability of the background flow.
iv.1 Linear perturbations
Let us first consider whitehole solutions of the KdV equation, as they are technically simpler to characterize. We work with functions and given by Eq. (73), with , , and . The trivial solution then corresponds to a white hole flow. We consider some initial perturbation . The corresponding initial conditions on are
(27) 
The timeevolution of can be determined by expanding it into outmodes. Due to the symmetry between inmodes of blackhole flows and outmodes of whitehole flows Macher and Parentani (2009b), the calculation is similar to the one sketched in Appendix C. In particular, the structure of divergences is the same, except that those which multiplied incoming waves now multiply outgoing ones, and conversely. This introduces an important difference, as the divergence in of the coefficients of dispersive waves for an extended perturbation, which did not contribute for a blackhole as it multiplied incoming waves, now multiplies outgoing waves Mayoral et al. (2011). As explained in Coutant and Parentani (2014), it adds a saddlepoint contribution for , which generates a stationary undulation with a large amplitude.^{2}^{2}2In Coutant and Parentani (2014) only perturbations localized in were considered. In that case, corresponding to perturbations in with a vanishing integral, the amplitude of the undulation vanishes as . When the perturbation satisfies instead, we verified that its amplitude goes to a finite, nonvanishing constant for . Numerical simulations using the linearized KdV equation confirm that this extends to smooth white hole configurations, see Fig. 11. Since the scattering coefficients on a whitehole flow described by the GP equation have the same divergences for as those obtained using the KdV equation, see Macher and Parentani (2009b), the same argument tells us that a stationary undulation shall also be produced in a condensate. In all cases, the macroscopic character of the undulation amplitude indicates that nonlinear effects will play a crucial role, which implies that linear equations are unable to predict the late time evolution.
iv.2 Nonlinear evolution
The nonlinear evolution of perturbations on a whitehole flow was studied numerically in Mayoral et al. (2011); Michel and Parentani (2015). In the present appendix we report some new numerical results. To make contact with the nonlinear analysis of blackhole flows of the previous section, we focus on the case of the GP equation. We obtained qualitatively similar results for the KdV equation, with the signs of the perturbations to the boundary conditions reversed. A rule of thumb, which works for the GP and KdV equations as well as for a superluminal KdV equation obtained by changing the sign of the dispersive term, is that analogue white hole flows are most unstable to perturbations with the sign of the difference between the stationary soliton and the corresponding homogeneous solution. For the GP equation, the soliton is a local underdensity, so the strongest instabilities come from perturbations of the density with a negative sign. For the KdV equation instead, the soliton is a local surelevation of the free surface. Correspondingly, a positive perturbation on leads to generally wilder behaviors than a negative one.
To start the analysis, we work with steplike functions and , given by Eq. (7) with and . We first consider an initially homogeneous configuration with close to and . Two typical solutions at intermediate times are shown in Fig. 12.
When , a finiteamplitude undulation develops in the supersonic region. Its amplitude is linear in and is constant at late time. When instead, the horizon emits superposed soliton trains in the subsonic region (three of them can be seen in the figure), and a perturbed undulation in the supersonic one. A largeamplitude perturbation is produced periodically close to the horizon, with a frequency linear in , which then separates into several solitons with different velocities in the subsonic region, plus a perturbation propagating on top of the undulation in the supersonic one. It should be emphasized that for homogeneous initial configurations there is no threshold on the value of . When it is positive (negative), one obtains a stationary undulation (nonstationary soliton train). For both signs of , the nonlinear solution is characterized by the nonanalytic parameter . In other words, irrespectively of the perturbation amplitude, the Bogoliubovde Gennes equation (14) cannot be used to determine the late time configuration.
We now study the evolution of localized perturbations. To this end, we choose initial conditions of the form
(28) 
For , we observe only the emission of an undulation in the supersonic region. For , we observe the emission of a finite number of solitons in the subsonic region, as well as a perturbed undulation in the supersonic one. In both cases, the amplitude of the undulation is roughly linear in at early times provided . Interestingly, we observe that it slowly decreases in time due to nonlinear effects, apparently going to zero for . As in the case of homogeneous initial configurations, irrespectively of value of , the linear equation (14) cannot determine the late time solution. This is due to the accumulation of the low frequency configurations on the sonic horizon Coutant and Parentani (2014).
Finally, we performed numerical simulations with functions and of the form Eq. (III.2) with , so that a solution with a homogeneous density exists. The main difference with the above results is that, when working with localized perturbations, solitons are emitted only if is below a negative threshold value . For not too close to , we observed that the time needed to produce the first soliton scales as . This can be understood as the condition for obtaining a sufficiently large underdensity so as to allow for the emission of a soliton. The simulations we performed were not precise enough to accurately determine the scaling of in and , although we found that increases with and , going to 0 for or . It would be interesting to further investigate these questions along with the validity of the linear equation (14) for .
V Conclusions
In this paper we studied onedimensional transonic solutions of the GP and KdV equations. We showed that they exhibit behaviors which are in close analogy with those of black hole solutions of general relativity. At the level of stationary solutions, we showed that the set of solutions which are asymptotically homogeneous (AH) on both sides is discrete at fixed value of the current . When considering steplike potentials, we demonstrated that the series of solutions parameterized by is unique. For smooth potentials, we numerically found a series of AH flows which is smoothly connected to the above series, as the solutions coincide in the high gradient limit. However, under some conditions, we also found a second series of AH solutions which is disconnected from the first one. These hairy solutions possess a large fraction of a soliton attached to the sonic horizon. Our preliminary investigations indicate that they are less stable than those of the first series because the soliton can be sent away from the horizon by a sufficiently large perturbation. In the rest of the paper, we focused on the stability properties of the first series of AH solutions.
At the level of linear perturbations, we found, both analytically and numerically, that the near horizon amplitude of all localized perturbations decays in time as a power law. This establishes that AH black hole flows are linearly stable. It should be noticed that the matrix which governs this linear scattering is the same as that encoding the Hawking effect in the present settings. One clearly sees here the close link between the (stimulated) Hawking process, i.e., the wave amplification upon scattering on the horizon, and the expulsion of all incident perturbations away from the horizon. Moreover, in the limit where the dispersive momentum scale is sent to infinity this expulsion follows the relativistic prediction.
To study the stability when including nonlinearities of the GP equation, we proceeded along two different approaches. We first used the approximate scheme of Whitham’s modulation theory to characterize in analytic terms the late time evolution of the solutions for steplike potentials. We showed that some field configurations are expelled from the horizon region to infinity by three nonlinear waves, known in the literature as dispersive shock waves and simple waves. Importantly, in the vicinity of the sonic horizon, the solution tends to one of the AH that we formerly characterized. It should be pointed out that Whitham’s theory also provides a characterization of the domain of (homogeneous) initial conditions which evolve at late time to an AH transonic flow. Finally, when taking the limit of small amplitude, it can be verified that these results generalize those of the linear analysis.
We then performed numerical simulations. We first showed that the late time behavior of a much wider class of initial configurations agrees with that predicted by Whitham’s theory, namely, an AH transonic flow is obtained by the emission of three nonlinear waves plus perturbations that decay in time. We verified that the properties of the three waves are in agreement with those obtained using Whitham’s theory. We also showed that the time dependent perturbations, which are not accounted for in our nonlinear analytical approach, decay in time with the same power law as that found in the linear analysis. All these results indicate that AH transonic flows are local attractors for neighboring flows. In a future work, using the integrability of the GP equation, we hope to be able to demonstrate this property. Since AH black hole solutions are local attractors, they can be produced without finetuning the initial conditions nor the potential . This should help observing both the spontaneous and stimulated analogue Hawking emission.
We also numerically studied the behavior of solutions when the initial conditions are outside the validity domain of the solutions obtained with Whitham’s theory. In this regime of large deviations with respect to the attractor solutions, several behaviors have been observed. In some cases, we found that the emitted nonlinear waves can leave behind them an undulation which propagates backwards towards the black hole horizon. In other cases, we observed the emission of soliton trains. This variety of behaviors is similar to that observed in Section IV when studying the time behavior of white hole flows. For these flows, the late time properties are rather complicated even when the perturbations have a small amplitude. Yet, the observed behaviors can be separated into two types. In this respect, our analysis indicates that the AH solutions are a kind of “separators”, rather than attractors, as the type of the solution is chosen according to the sign of the density fluctuation in the near horizon region. When there is a sufficiently large increase of the density, the solution displays a macroscopic undulation in the supersonic domain, whereas it gives rise to an emission of soliton trains when there is a sufficiently large density decrease. These two types of behaviors have been already found in the context of the dynamical instability (called the black hole laser) which is found when a stationary flow crosses twice the sound speed de Nova et al. (2015); Michel and Parentani (2015). Our study clearly shows that it is the white hole horizon which is responsible for the wide variety of temporal evolutions that was observed.
The present work leaves several open questions which deserve further study. First, it would be interesting to prove rigorously our conjecture using inverse scattering techniques. In the same vein, investigating deformed GP or KdV equations including nonintegrable terms could shed light on the relations between the mathematical properties of the equation and the “nohair” results. As a first example, we study numerically the case of the cubicquintic GP equation in Appendix D. Another possible extension concerns higherdimensional systems. In general relativity, nohair and uniqueness results crucially depend on the dimensionality of space Chrusciel et al. (2012). It would certainly be enlightening to see how the dimensionality affects black hole stability in systems described by the GP or hydrodynamic equations.
Acknowledgement
We thank A. M. Kamchatnov, N. Pavloff, and G. Shlyapnikov for interesting discussions. This work was supported by the French National Research Agency under the Program Investing in the Future Grant No. ANR11IDEX000302 associated with the project QEAGE (Quantum Effects in Analogue Gravity Experiments).
Appendix A Whitham equations
In this appendix we give the main steps in obtaining the Whitham modulation equations for Eq. (1) and their solutions used in the main text. The interested reader will find in Kamchatnov (2000) and references therein a full derivation. When possible we use the same notations and conventions as in this reference.
The Whitham modulation theory Whitham (1965) was developed to study oscillating solutions of partial differential equations with slowly varying parameters. It rests on the two following ideas. First, if there is a clear separation between the fast scale of oscillations and the slow scale of variation of the parameters, an averaging procedure can decouple them. Second, averaging conservation laws generically gives the most accurate and bestcontrolled results. This theory is then particularly useful when one has enough conservation laws to characterize the space of solutions one is interested in. As such, it is no surprise that deep links exist with integrability. A generic way to obtain the Whitham equations for an integrable system is to use the AKNS scheme, developped in Ablowitz et al. (1974) and applied to the Whitham theory in Kamchatnov (1994). Here we briefly review this procedure, following the presentation of Kamchatnov (2000).
a.1 The AKNS scheme
The basic idea of the AKNS scheme is to reformulate a partial differential equation one wishes to solve as the compatibility condition of a linear system of the form
(29) 
where is a twocomponent complex vector and , are two matrices of the form
(30) 
where , , , , , and are differential operators, analytic in . Here is a complex spectral parameter, independent on and . To simplify the notations, from now on the dependence in will not be written explicitly when no confusion is possible. The compatibility condition of Eq. (29) is , i.e.,
(31) 
Importantly, this system must be equivalent to our original partial differential equation for all values of . Then will generate an infinite number of conserved quantities which can be used to solve the problem, either exactly using the inverse scattering method when it applies, or approximately using the Whitham equations.
To see this, it is convenient to define two linearly independent solutions and of the linear problem and the three scalar quantities
(32) 
Partial derivatives of , , and can be computed straightforwardly. We find
(33) 
Conversely, the compatibility conditions of the system Eq. (33) give back Eq. (31) provided . Eq. (33) can thus be seen as a rewriting of the original problem. This formulation is particularly useful for deriving conserved or slowlyvarying quantities, as we now explain.
We first notice that is directly related to the generalized wronskian of :