The imprint of the analogue Hawking effect in subcritical flows

We study the propagation of low frequency shallow water waves on a one dimensional flow of varying depth. When taking into account dispersive effects, the linear propagation of long wavelength modes on uneven bottoms excites new solutions of the dispersion relation which possess a much shorter wavelength. The peculiarity is that one of these new solutions has a negative energy. When the flow becomes supercritical, this mode has been shown to be responsible for the (classical) analog of the Hawking effect. For subcritical flows, the production of this mode has been observed numerically and experimentally, but the precise physics governing the scattering remained unclear. In this work, we provide an analytic treatment of this effect in subcritical flows. We analyze the scattering of low frequency waves using a new perturbative series, derived from a generalization of the Bremmer series. We show that the production of short wavelength modes is governed by a complex value of the position: a complex turning point. Using this method, we investigate various flow profiles, and derive the main characteristics of the induced spectrum.


I. INTRODUCTION
The idea to use fluid flows to mimic the Hawking effect of black holes [1], which allows them to spontaneously emit a thermal radiation, has been intensively studied at the theoretical level [2][3][4][5]. More recently, several experimental studies have been set up, in various media as diverse as Bose-Einstein condensates [6,7], optical fibers [8,9] or surface waves, either in water [10][11][12][13], or in superfluids (therein called "ripplons" [14]). To obtain such a setup, one needs a fluid flow whose velocity crosses the speed of waves. Moreover, because the Hawking effect necessarily involves short wavelength modes, it also necessary to take into account dispersive effects that arise at short distances. One promising possibility is to use surface waves on flowing water [15,16], since their propagation speed is much lower than that of sound waves. When neglecting capillarity and dissipation, surface waves propagate on a one-dimensional homogeneous flow with a frequency ω and a wave number k that obey the dispersion relation [17,18] where Ω = ω − vk is the comoving frequency. In this equation, g is the local gravitational field, v the flow velocity, and h B the depth of water. To characterize the flow, it is convenient to introduce the Froude number, defined as F = v/c, where c = √ gh B is the propagation speed of long wavelength waves. If F > 1 (resp. F < 1), the flow is called "supercritical" (resp. 'subcritical'). The transition from subcritical to supercritical, a "transcritical flow", is the analog of a black hole if the flow accelerates, and a white hole if the flow decelerates. Unfortunately, it is experimentally delicate to obtain controllable transcritical flows to study the analog Hawking radiation. One difficulty is caused by the appearance of an undular pattern, or "undulation" that deforms the free surface and whose amplitude raises with the Froude number [17,19,20]. Instead, experimental studies have so far focused on flows of high Froude numbers but that stay subcritical [10][11][12][13]. In these cases, due to dispersive effects one still observes the production of the negative norm mode responsible for the Hawking effect in transcritical flows. However, theoretical treatments on the analog Hawking effect have mostly focused on transcritical flows, and it is therefore presently unclear what governs the spectrum of this negative norm mode production for subcritical flows. Recent numerical works [21][22][23] have indicated that the spectrum is in general quite different when the flow is subcritical with respect to the transcritical case.
In this work, we provide an analytical characterization of the low frequency scattering when the flow is inhomogeneous, i.e. when c, v, and h B depend on the position x. For this, we developed a new mathematical approach, based on a generalization of the Bremmer series [24][25][26]. This allows us to obtain a perturbative expansion of the various scattering coefficients in gradients of the background. More precisely, the "small parameter" of the expansion will be played by the variation of the height of the flow, i.e. |h B |. The first order treatment shows that the values of the scattering coefficients are mainly governed by complex turning points. At very low frequencies, all the turning points reduce to a single complex horizon, which is the locus where the Froude number F (x) reaches 1 if it is analytically continued to complex positions x. Although we focus in this work on water waves, we believe that our conclusions are still valid for subcritical flows (or what replaces it) in other analog gravity systems where dispersion decreases the velocity at short wavelengths, e.g., optical fibers [8,9], or sound in a duct [27]. The paper is divided as follows. In the first section we present the setup and the wave equation.
In the second one, we derive the perturbative series, and show that the first order is given in terms of contour integrals involving the complex turning points. In the last section, we apply this general framework to specific flow examples, discuss the various regimes and the relevant physics for present experiments. In Appendices, we provide very general proofs, preparing our results for further extensions.

A. Surface wave equation
We shall consider the propagation of water waves in the so-called "weak dispersive regime". In this regime, the dispersion relation (1) is approximated by the first two terms of the low k expansion, i.e. k tanh(h B k) ∼ h B k 2 − h 3 B k 4 /3. As mentioned in the introduction, it is necessary to take into account dispersive effects since the scattering processes we are interested in involve short wavelength modes. When the maximum value of the Froude number F max is close to 1, a case referred to as 'near critical flows', the weak dispersive regime provides a good approximation of the scattering coefficient. However, we believe that even for flows that are not near critical, the qualitative features we describe will be very similar. In the weak dispersive regime, gravity waves are described by the action [20] In the following, v is assumed to be positive, so that water flows from left to right. The field φ encodes the fluctuations of the velocity potential of the flow at the surface. It is directly related to the change of height of the free surface. In the presence of a (linear) wave, the water depth becomes h B (x) + δh(t, x). The surface elevation δh is then given by where this follows from Bernouilli's equation [16,20]. Minimizing the action (2) gives us the equation of motion for the field We assume that the background flow is stationary, and therefore, we look for solutions of Eq. (4) at fixed frequency, i.e. of the form φ = Re (φ ω (x)e −iωt ), where φ ω (x) is a complex stationary mode. Since time-dependent solutions are obtained by taking the real part, it is enough to work with ω > 0. The modes φ ω (x) satisfy the equation Before trying to solve this equation, it is useful to analyze its main properties. As it is derived from an action, it possesses a canonically conserved norm, given by This norm plays a crucial role in the characterization of the scattering. For positive frequency modes, the sign of the norm coincides with that of the energy. As we shall see, due to dispersion, the system possesses negative energy modes, or equivalently, modes with a negative norm (6).
This is characteristic of unstable flows [28]. The generation of a negative norm mode by sending a positive norm one is referred to as "anomalous scattering". For transcritical flows, this scattering (in the smooth limit) is the classical analog of the Hawking effect 2 . For subcritical flows, it is still present, but it was so far unclear what governs the spectrum, i.e., the values of the scattering coefficients for various frequencies. As we shall deal exclusively with stationary modes, it is more convenient to work with the conserved current rather than the norm (6). Because of the dispersive term, this current is not the standard Klein-Gordon current, but has a more complicated form [30]. Starting from the action (2), it reads For any mode solution of (5), the current is (exactly) conserved, i.e. ∂ x J = 0. This current represents the amount of norm that is transported by a mode. Its conservation is of course equivalent to that of the norm (6). This can be directly seen from the identity ∂ t ρ + ∂ x J = 0 (see App. B 3), which follows from the application of the Noether theorem to Eq. (2). Notice also that for ω > 0, ωJ[φ ω ] is the energy current [20]. When the background flow is homogeneous, i.e. c, v, and h B are constant, the solutions are superpositions of plane waves e ikωx . Here, k ω is the wave number, or momentum, and satisfies the dispersion relation Below a certain threshold frequency ω crit , this equation possesses 4 distinct roots (see Fig. 1). Two of them have long wavelengths, while the two others have short wavelengths. The first two are the usual left-mover (noted k u , as it moves against the flow, i.e. "upstream") and right-mover (noted k d , for "downstream"). The two other roots, which are absent when the flow velocity vanishes, are due to both the nonzero flow and dispersion. One of them, k − , has a positive value but a negative norm. The other have a negative value and a positive norm, and is denoted k + . (The index refers to the sign of the norm.) The negative norm mode described by k − will play a crucial role in the following. The aim of this work is to study, when the flow becomes inhomogeneous, how these four modes mix. In particular, we shall see how the propagation of a long wavelength left-mover k u generates the short wavelength modes, as was experimentally observed in [10,13]. For flows that become critical, this generation is the classical analog of the Hawking effect. For subcritical flows, such a mode conversion still exist, but the law governing the scattering coefficients was so far not known analytically. This is what we aim at characterizing.

B. Characterization of the background flow
We assume that the fluid flows over a smooth obstacle. The height of fluid h B (x) varies monotonically from an asymptotic value on the left side, to a minimum value h min and increases again to a constant value on the right side. When the obstacle is smooth enough, the (unperturbed) free surface stays approximately flat, and in this case, the other background quantities are directly deduced from the height by the relations where g is the local gravitational acceleration, and q the (conserved) flow rate (water flux per unit width, expressed in m 2 · s −1 ). As we shall see, the most relevant quantity to describe the flow is the local value of the Froude number where the second equality is satisfied when (9) is. In realistic flows, curvature effects (but also dissipation) will deform the free surface, and the relation between these quantities becomes more intricate 3 . Unless otherwise specified, we will treat the three functions v, c, and h B as independent, thereby leaving the possibility to include corrections to Eq. (9). However, to keep control on the various approximations, we assume that the gradients, in units of the dispersive scale, are essentially of the same order, i.e., that |h B v /v| and |h B c /c| are of the same order as |h B |, and we refer to the "smooth limit" as |h B | 1. This is automatically the case if the three functions are related by Eq. (9). Far from the obstacle, we assume that the background quantities v, c, and h B are constant. Over this obstacle, the flow velocity v increases to a maximum v max while the wave speed decreases to a minimum c min . At the top of the obstacle, the Froude number reaches its maximum F max (see Fig. 2). The assumptions we make are in practice nontrivial. First it assumes that no turbulence is formed by the flow close to the free surface. While this is reasonable for subcritical flows, we also assumed that no undulation appears at the free surface. When increasing the Froude number, even below 1, such an undulation is more likely to form, as observed in [10,13]. We believe that the presence of such an undulation could be treated by our framework (see the remark of footnote 4), but the computations will be more involved. We feel that such an analysis goes beyond the scope of the present paper.

C. The WKB approximation
When v(x), c(x), and h B (x) vary, plane waves are no longer solutions of Eq. (5). In the limit of a very smooth background |h B | → 0, solutions of Eq. (5) are given by WKB modes, i.e., locally plane waves characterized by a local momentum k ω (x). This local momentum is a solution of the Hamilton-Jacobi equation, which is nothing else than the dispersion relation (8) in an inhomogeneous background Each solution of k j of this equation depends on both ω and x. To lighten the notations, we shall drop this dependance when unnecessary. Throughout this paper, we also assume that the 4 roots are real and distinct for all x. Since the flow stays subcritical all along, i.e. F max < 1, this is realized below a threshold frequency ω min = min x (ω crit (x)). For near critical flows, the value of this frequency reads For ω < ω min , a WKB mode is then given by where the subscript j indicates the corresponding Hamilton-Jacobi root, i.e. j ∈ {u, +, −, d}. A careful analysis (see e.g. Appendix A of [31] or Appendix A of this work) shows that in the limit of smooth backgrounds, the amplitude simply reads where v g is the group velocity of the corresponding mode, and Ω its co-moving frequency, defined after Eq. (1). Moreover, at the level of the WKB approximation, the current J is easy to compute and one sees that Hence, the WKB amplitude (14) normalizes the current of a WKB mode to ±1 [30]. This property of WKB modes comes from the fact that the current J is an adiabatic invariant of the problem [32]. Unfortunately, the WKB approximation precisely consists in neglecting the mode mixing, which is what we are after. To overcome this problem, we first notice that for ω < ω min , no crossing occurs, i.e. the 4 roots of (11) are distinct for all x. Therefore, the 4 WKB modes in (13) are perfectly well-defined functions of x. Instead of using them as approximate solutions of the wave equation, we shall use them as a new basis to represent exact solutions of the wave equation (5). This allows us to recast the wave equation in an equivalent form, adapted to a perturbative expansion of the scattering coefficients in the background gradients.

A. The Bremmer representation
The idea to use the WKB modes as a basis has been widely studied and used for second order differential equations, where it is called the Bremmer series. It has a wide range of applications, from scattering theory of the Schrödinger equation [26], or wave propagation in inhomogeneous media [33], to particle production in early cosmology [34,35]. Here however, we must use an extended version of this method, as the problem is intrinsically higher order. As explained in the previous section, the S-matrix is 4 × 4. In Appendix A, we present detailed proofs of how to extend the Bremmer series for higher order equations. In this section, we present the method without technical calculation, in order to focus on the physics and the significance of this new representation. The key idea is to write general, exact solutions of Eq. (5) as superpositions of WKB waves, where the various amplitudes are x-dependent, i.e.
At this level, the function A j (x) are unspecified functions, and are not given by (14). Since this introduces 4 unknown functions, instead of 1, we impose 3 extra conditions. The idea is to decompose also the first, second and third derivatives of φ(x) on the WKB basis, and the forth derivative will then be given by the equation of motion (5). Explicitly, we assume, in addition to (16), (17c) Using these 3 conditions and the main ansatz (16), we show that the knowledge of φ(x) is equivalent to the knowledge of A u (x), A + (x), A − (x) and A d (x). The 4 equations combine to give the single matrix equation where V is the Vandermonde matrix of the 4 roots k u , k + , k − , and k d , i.e.
Because the three roots are distinct, det(V ) = 0, and hence, the relation between φ and its derivatives and (A u , A + , A − , A d ) is one-to-one. Therefore, the wave equation (5) can now be entirely recast in an equivalent equation for the local amplitudes A j (x). To obtain the equation satisfied by the local amplitudes A j (x), we plug the ansatz (18) in the wave equation (5). Since the first three derivatives of φ are given by (18) The off-diagonal elements of M are easy to interpret: they give the coupling between the different WKB branches, due to the varying background. Those are responsible for the nontrivial scattering. On the other hand, the diagonal terms of Eq. (20) represent the adiabatic evolution of the amplitudes A j (x). To further simplify the equation, we can integrate these diagonal terms by working with normalized amplitudes. For this, we define where N j is chosen so that the first term of (20) disappears. This gives a first order equation on N j , which directly integrate as N j = exp x M jj (x )dx . As we show in App. A 2 and B, this leads to We recognize here nothing else than the WKB amplitude given in Eq. (14). This is not a surprise, as N j gives the adiabatic evolution of the amplitudes. We shall refer to these new coefficients a j (x) as the local scattering coefficients. At the level of the WKB approximation, they are constant. When the background varies, these coefficients becomes non constant, meaning that the propagation of one mode excites the other ones, leading to nontrivial asymptotic scattering coefficients. These coefficients are governed by a first order equation, directly obtained from (20), and which reads This equation possesses several key features, that we now wish to underline. First, this equation is strictly equivalent to the original equation (5). No approximation have been used so far, but this rewriting is very adapted to a perturbative resolution. Second, the coupling coefficients M j are proportional to derivatives of the background. In the limit |h B | 1, they are small and have a slowly varying phase (see App. A 2 for their exact expressions). Because of this, the coefficients a j mainly couple through the change of their WKB phases e i (k (x )−k j (x ))dx . This structure implies that the scattering will become significant when this phase difference satisfies a resonance condition (see next section). The last key property of Eq. (23) is obtained when computing the conserved current (7) in terms of the local scattering coefficients. Since J involves only the first three derivatives of φ, the ansatz (18) guarantees that the computation of J is identical as in the case of plane waves. After some effort (shown in App. B 3), we show that Once again, this equation is exact. It guarantees that the scattering governed by Eq. (23) conserves the norm of Eq. (6). Also, from Eq. (23), the conservation of the current implies that the matrix M has some symmetric/antisymmetric properties, something that is not transparent from their explicit expressions (given in App. A 2).

B. The complex turning points
We now turn to the evaluation of the scattering coefficients in the limit of smooth backgrounds |h B | 1. Since 4 modes exist on both sides, there are 4 incoming legs and 4 outgoing ones, and the complete S-matrix is 4 × 4. In a quantum mechanical language, the scattering coefficients can be seen as transition amplitudes for a mode transition k → k j . These transitions can then be estimated in perturbation theory of Eq. (23), i.e. in an expansion in the matrix elements M j , which are small in smooth backgrounds |h B | 1. For the present purpose, we shall consider only one specific scattering mode, but our results easily extends to the others. We consider a long wavelength mode coming in from the right, which means that a u (+∞) = 1, and Fig. 3). This fixes half of the asymptotic values of the local scattering coefficients. The other half gives the scattering coefficients (see Fig. 3). T and R are the transmission and reflection coefficients between the two long wavelength modes, while α and β describe the generation of the short wavelength modes. Using (24), the conservation of the current imposes the following relation between the scattering coefficients We see that the coefficient β contributes with the unusual sign.
T α β Obstacle At zeroth order in M, a u (x) ∼ 1, while the other a j (x) vanish. When inserting this on the right hand side of (23), and integrating from −∞ to +∞, we obtain the first order expression in M of the scattering coefficients. At this order, T is 1, while the other three are given by an integral expression. To start, we focus on the computation of α. By solving Eq. (23) at leading order, only the coefficient M +u contributes and we have This gives the first order expression for the coefficient α. It is possible, starting from Eq. (23) to derive an expression for α at any order in M. It has been shown in various cases that the obtained series is generally convergent [36,37], and that the convergence is usually quite fast [26,38]. In Fig. 4, we give a diagrammatic representation of the perturbative series. In a regime where the various scattering coefficients are small, |α| 1, |β| 1, and |R| 1, it is legitimate to truncate this series at first order, since higher orders will be essentially given by higher products in these quantities. This is true in the smooth limit |h B | 1, but we also need ω to be sufficiently far from ω min , otherwise we would have |α| = O(1) (although ω min − ω might in practice be quite small and |α| 1 still valid, see e.g. Fig. 5(a)). We now assume that this is the case, and study the consequence of the first order result (26). The main contribution of the integral governing Eq. (26) comes from the saddle point of the exponential. This saddle point satisfies the equation By assumption, this equation is not verified by any real x. However, when the background functions are analytic, there exist complex solutions x * ∈ C. If x * were real, it would correspond to a turning point, and, hence, in our case we call x * a complex turning point. Decomposing it in real and imaginary parts, It follows from (26) that the α coefficient is given by the contribution of the saddle point as where C is a constant prefactor (discussed below). In Eq. (29), x 0 is a real reference point, which can be chosen anywhere. If several turning points are present, α is given by a sum of the contribution (29) for each of them. Usually, the ones that are the closest to the real axis give the dominant contributions, while the others produce only exponentially small corrections 4 . As we shall see in Sec. IV B 2, very symmetric flows typically have two main interfering contributions to the scattering coefficients. Taking the modulus of (29), we have As explained in App. A 3, the complex turning point must be chosen such that the contour integral in (30) has a positive imaginary part. It follows that α is generally exponentially small, which is a common feature of low gradients or adiabatic limits [26,39]. The prefactor C in (30) is rather delicate to obtain. In the smooth limit |h B | → 0, we show that it tends to 1 (see App. A 3). However, this limit fails at reproducing the ultra low frequency behavior of the coefficients. The reason is that the limits |h B | → 0 and ω → 0 do not commute. Indeed, when the gradients are nonzero but small, in the limit ω → 0, the prefactor vanishes as as we show in App. B 2. The characteristic frequency ω s is estimated in Eq. (B14). The key point is that ω s is proportional to |h B |, and, hence, becomes very small in the smooth limit, and in particular, ω s ω min (ω min defined in Eq. (12)). To summarize, the prefactor C is characterized by two regimes. When ω s ω ω min it is 1, but for ω ω s it is given by Eq. (31). A similar computation for β leads to a similar expression, The complex turning point for β is a priori different from that of α, since it obeys a different resonance condition k u (x β * ) − k − (x β * ) = 0. In App. B 2 we show that when ω s ω min , the prefactor is essentially the same as for α. As we see from (30) and (32), in the smooth limit, the scattering coefficients are exponentially small.
The last coefficient R, giving the mode mixing between the two long wavelength modes (see Fig. 3), can also be evaluated perturbatively, and possesses a contour integral expression as (30) and (32) involving a different complex turning point x R * . The first order expression of R is however more delicate. Indeed, for very low frequencies, this first order expression becomes of order 1, meaning thatx the perturbative treatment breaks down. On the contrary, the expressions for α and β stay small in the limit ω → 0. The reason for this discrepancy can be seen in the expressions of the matrix elements of M. While M du is proportional to the background gradients |h B |, M +u and M −u are further suppressed by a factor O(ω 1/2 ). Heuristically, we explain this by the fact that R governs a transition involving only long wavelength modes, while α and β involve a short wavelength one, which improves the accuracy of the perturbative treatment even at low frequencies ω ω min . Fortunately, it has been numerically obtained in various works [22,23,40] that the reflection coefficient R stays small for all frequencies.
Since we are mainly interested in the production of short wavelength modes by k u , we shall ignore the mode k d in the sequel. Note that to obtain second order estimates of α and β, k d can no longer be ignored since it will appear as an intermediate state in the transitions k u → k + or k u → k − (see Fig. 4).

IV. APPLICATION TO NEAR CRITICAL FLOWS
A. The simplest example: Case of a short obstacle We shall start by analyzing a simple example. This will allow us to present the techniques to explicitly evaluate Eq. (30), in a case where the computations stay relatively simple. For this we assume that the mode mixing essentially takes place in a close vicinity of the top of the obstacle, i.e., where F max is reached. In Sec. IV B 2 we will give a more precise meaning to this "short obstacle limit". Under this assumption, the evolution of the Froude number is well approximated by a second order Taylor expansion near its maximum The parameter d characterizes the length of variation of the Froude number near its maximum value. As we shall see, this quantity directly affects the scattering coefficient. In addition, to simplify the discussion, we present the results in two steps depending on the ratio ω/ω min (but without assuming anything concerning the ratio ω/ω s ). We first study the limit ω/ω min 1, and in a second part, study the corrections in ω/ω min .

Low frequency limit
In the limit ω ω min , the resonance conditions for α and β become the same, and reduce to Since this condition gives the location of the horizon when the flow is transcritical, for low frequencies, the complex turning point can be interpreted as a complex horizon. Using the profile of Eq. (33), it is given by To obtain the correct sign of the integral in Eq. (30), we must choose Im(x * ) > 0. Moreover, when ω → 0, the roots of the Hamilton-Jacobi equation (11) become simpler, and we find We notice here that the root difference (36) scales like (1 − F ) 1/2 . Therefore, in the limit 1 − F max 1, it is only necessary to consider the variations of the function 1 − F (x). The other background quantities can be approximated by their value near F max , since taking into account extra terms will produce subleading corrections in 1 − F max . In this limit, the profile of Eq. (33) gives We now use this expression to compute the complex integral governing the scattering coefficient α through Eq. (30). For convenience, we chose the reference point x 0 = 0, then We deduce the amplitude of the α coefficient where α 0 is short for α ω ω min . We see that in the regime ω ω min , all the frequency dependence is in the prefactor C ω (we added the index ω with respect to Eq. (30) to emphasize this point).
Since we assume nothing concerning the ratio ω/ω s , C ω varies from 1 to ω/ω s when ω decreases. By a similar computation, we show that β 0 has the same amplitude for very low frequencies, i.e. |β ω→0 | 2 ∼ |α ω→0 | 2 . This can be seen by direct computation, but comes in fact from a more general property of the mode equation. Indeed, the change φ ω → (φ −ω ) * leaves the mode equation (5) invariant, and as can be seen by looking at the roots of (11), exchanges the role of α and β. This leads to the relation The property above has been widely used in Hawking radiation studies [3,31,41]. Here also, it significantly simplifies the computations of the scattering coefficients.

Frequency dependence
The corrections to Eq. (39) in ω/ω min are more delicate to obtain. These corrections have two origins. The first is the shift of the value of the complex turning point, and the second is the exact expression of the roots k u , k + , and k − . For the latter, one needs to solve the Hamilton-Jacobi equation (11). Unfortunately, one cannot simply compute the corrections to Eq. (36) for ω ω min because such corrections will not be accurate close to the turning point. We can still simplify Eq. (11) by discarding the last root k d , which plays essentially no role at first order in perturbation theory. The Hamilton-Jacobi equation (11) is then reduced to a third order equation in k. To obtain it, we carefully take the square root of (11) so as to select the relevant branch (see Fig. 1), and expand the result up to O(k 3 ). This gives A direct comparison of this equation with Eq. (11) shows that the three roots k u , k + , and k − are approximated by the roots of (41) up to small corrections in 1 − F 1 5 . We now obtain 5 To see this, we first notice that since we performed an expansion in k, the highest error made is on the value of k − (ω = 0), which has the highest value (see Fig. 1). From Eq. (11), it is given by h the roots by solving this equation using the Cardan-Tartaglia method. To start, the associated discriminant gives the equation for all the complex turning points, i.e. it gives the condition for two roots to merge, Similarly to Eq. (36), at leading order in 1 − F max , it is only necessary to consider the variations of 1 − F (x), while h B and c are well approximated by c min and h min . Doing so, the complex turning point for α is given by The other complex turning point x β * is another solution of Eq. (42). For ω = 0, both emerge from the complex horizon x 0 * , the difference scaling like O((ω/ω min ) 2/3 ). To simply express the roots, we introduce the auxiliary functions The Cardan-Tartaglia method then gives the three roots as combinations of U + ω and U − ω . In particular, We now evaluate this near the top of the obstacle, that is, using Eq. (33). At leading order in 1 − F max , we have We are now ready to compute the complex integral governing the coefficient α in Eq. (30). We start by writing where we defined the functions I ± by Combining the preceding results, and applying Eq. (29), we finally obtain This gives the expression of α ω in the flow profile of Eq. (33). Eq. (49) should be valid up to ω ω min , as long as |α ω | 1. By using the same method, Eq. (32) leads to a similar expression for the coefficient β ω . To obtain it, one can either redo the calculation of the complex integral, or more quickly, carefully apply Eq. (40). The results are presented on Fig. 5(a). Since the full expression of the I ± functions is rather complicated, it is instructive to look at the limit ω → 0, and see how α ω (resp. β ω ) deviates from Eq. (39). Interestingly, the functions I ± are not differentiable for → 0 and therefore, small frequency corrections display non-analytic terms. Indeed, after some efforts, one can show that This gives approximate expressions for the scattering coefficients As we observe on Fig. 5(a), the low frequency expressions Eq. (51) are quite accurate up to ω ω min (where the perturbative expression (30) can no longer be trusted). At this level we would like to emphasize several qualitative features displayed by Eq. (51) that are maintained for more general profiles. First, when ω → 0, |α ω | 2 ∼ |β ω | 2 and both vanish as O(ω) due to the prefactor (see Eq. (31)). Second, when ω/ω min increases, |α ω | 2 becomes larger than |β ω | 2 . Third, the corrections in ω/ω min display non-analytic terms, in O(ω ln(ω)). It is also quite instructive to analyse the behavior of the ratio r ω = |β ω /α ω | 2 . Indeed, the linearity of the logarithm of this ratio in ω has been used in the literature as a sign for the thermality of the emitted spectrum. Moreover, this ratio is also independent of the prefactor C ω . For low frequencies, Eq. (51) gives On Fig. 5(b), we plotted the evolution of r ω , using both Eq. (49) and the low frequency expression (52). As we see, despite the presence of non-analytic corrections, r ω looks fairly linear in ω. However, this ratio alone misses several features of the scattering that differs from the Hawking regime, and in particular the low frequency ω ω min behavior of α and β.

Monotonic profiles
We shall start by analyzing the case of a profile whose Froude number, increase monotonically from a minimum to a maximum value. This case is very useful to better understand the more general profiles of Sec. IV B 2, but it also has its interests in its own right. To model such a profile, we assume that the Froude number is given by The flow starts from a low Froude number F min = F 0 − D on the left side, and smoothly rises to reach F max = F 0 + D. The parameter γ gives the slope of the profile. As in the preceding section, it is simpler to first look at the low frequency limit ω ω min and in a second time discuss the corrections in ω/ω min . At low frequencies, the physics is dictated by the complex The prefactor C ω is given by Eq. (B16). We clearly observe three distinct regimes: for ω ω s , |α ω | 2 ∼ |β ω | 2 ∝ ω, for ω s ω ω min , |α ω | 2 ∼ |β ω | 2 almost constant in ω, and for ω s ω ω min , |α ω | 2 increases, while |β ω | 2 decreases. Right panel (b): Ratio r ω as a function of ω. In both plots, we have chosen the flow parameters such that d(1 − F max )/h min = 1, and ln(ω s /ω min ) −4.4. The solid lines are obtained using the full functions I ± , while the dashed ones are the approximations (51), and (52). Note that the present treatment cannot be trusted too close to ln(ω/ω min ) ≈ 0. horizon, i.e. the location satisfying Eq. (34), which governs the common value of α 0 and β 0 . From equation (53), we find the complex horizon 6 to be We then compute the low frequency value of α and β (see App. B 4), and we find We see that the value of the coefficient depends not only on the slope γ, but also on the height of the step, i.e. the parameter D. Moreover, we notice that it depends on 1 − F max with a different power law than in the short obstacle case (compare (55) to (39)). The reason for this is that unlike in the short obstacle case, the imaginary part of the complex horizon of (54) is independent of 1 − F max , hence α 0 depends on it only through the roots k u − k + . To obtain the corrections for ω = 0, we follow the same procedure as in Sec. IV A 2, and compute only the leading order correction. A rather tedious computation shows that for small ω/ω min , and 1 − F max 1, We notice that unlike in (51), the above equation shows no non-analytic terms. This turns out to be an accident of the profile of Eq. (53), where the leading order non-analytic terms (in O(ω ln(ω))) contributes only to the phase of α ω . This is no longer true for the next-to-leading corrections in ω/ω min or if the profile slightly differs from (53).

Non-monotonic profiles
We are now ready to analyze a more general class of flow profile, which have a similar shape than the ones studied numerically and experimentally [10,13,22,23] (undulation excluded). For this we assume that the Froude number is given by The maximum value of the Froude number F max is by assumption smaller than 1. The flow starts from a low Froude number F min = F 0 − D on the left side, rises to reach F max and then decreases again to F min . On the left side (resp. right side), the slope is controlled by the parameter γ l (resp. γ r ). Again, we first consider the low frequency limit ω ω min and then discuss the corrections in ω/ω min = 0. Depending on the parameters of the profile, we distinguish three different regimes: • long obstacles γ l,r L 1, asymmetric γ l = γ r , • long obstacles γ l,r L 1, symmetric γ l = γ r , • short obstacles γ l,r L 1.
When the obstacle is long (L is large), there exist two complex horizons, located around the top of each slopes. In the limit γ l,r L 1, they are given by the monotonic profile expression (54) centered at ±L/2, i.e.
The scattering coefficients are then given by a sum of two interfering contributions |α ω | 2 ∼ |α l | 2 + |α r | 2 + 2|α l ||α r | cos Re where α l,r are given by the single turning point expression (55) with γ = γ l,r . This equation allows us to draw several conclusions concerning the behavior of the scattering coefficient. If the profile is asymmetric, γ l = γ r , the one with the biggest slope dominates in the expression for |α ω | 2 . Indeed, the imaginary part of the corresponding turning point lies closer to the real axis, as it is inversely proportional to γ l,r . On the contrary, if the flow is symmetric, γ l = γ r , the two contributions have the same weight and interfere through the phase of Eq. (59). This phase shift is accumulated not only along the real line, between x l R and x r R , but also in the complex plane, from x l,r R to x l,r * . For a long obstacle (γ l,r L 1) and ω ω min , the phase shift of (59) is given by where ζ l,r are defined after Eq. (B35). When decreasing L, the two complex horizons keep the same imaginary part, but their real part get closer. At a certain critical value L = L c , they merge into a single solution 7 . For lower values L < L c , one of the solutions migrates closer to the real axis, while the other moves afar. If the profile is not perfectly symmetric, one observes something similar, but instead of merging together, the roots first get closer, and around the critical value of L, repel each other so that the one with the smallest imaginary part approaches the real axis, while the other moves afar (see Fig. 6). This mechanism is very similar to the "avoided crossing", well-known in quantum mechanics [42]. When L < L c , we enter in the regime of a short obstacle. In this case, one of the complex horizons dominates in the expression (59) for α. This solution has a real part close (equal if γ l = γ r ) to zero, i.e., it lies close to the top of the obstacle. This case becomes very similar to the one studied in Sec. IV A. At L = 0, the complex horizon closest to the real axis is found to be If we additionally have 1 − F max 1 − F min , (61) exactly reduces to the short obstacle case of Sec. IV A, meaning that Eq. (33) becomes a good approximation to describe the scattering.
Re(x $ ) When relaxing the assumption ω ω min , α ω and β ω differ but are still given by an interfering sum as in Eq. (59). As ω increases, the two turning points x α * and x β * emerge from x 0 * and migrate in different directions in the complex plane. The effect of this migration is twofold. First, for long obstacles, the relative location of the left and right turning point changes, and therefore, the phase (60) between the two interfering contributions in Eq. (59) varies. For some values of the frequency, this phase will be a multiple of 2π, and the coefficients show a dip, as the first order estimate in Eq. (59) vanishes. Such dips have been numerically observed in [22]. Second, the amplitudes of the single turning point contributions, i.e. |α l | and |α r | (resp. |β l | and |β r | for β ω ) are altered as in Eq. (56). On Fig. 7, we represented the evaluation of α ω and β ω for a symmetric and an asymmetric profile. On Fig. 8, we show how the coefficients oscillate in ω due to interferences between the two turning points as in Eq. (59), and how this affects the ratio r ω .

V. CONCLUSION
In this paper, we studied the scattering of low frequency waves on a subcritical fluid flow, that is, whose Froude number stays below 1. We developed a new method, based on a generalization of the Bremmer series, where exact solutions of the wave equation are written as a local superposition of WKB modes (see Eq. (16)). The coefficients of this superposition, which we called local scattering coefficients, are position dependent and possess several useful properties. First, they are by construction slowly varying. At some locations along the flow, they can transit from one constant value to another. This can be interpreted as the creation of a new mode. Second, their asymptotic values directly give the scattering coefficients. Third, the local scattering coefficients are governed by a first order differential equation, Eq. : plot of the ratio r ω = |β ω /α ω | 2 as a function of ω, for three long obstacles less and less symmetric: γ r /γ l = 1 (solid), γ r /γ l = 1.2 (dot-dashed), γ r /γ l = 1.5 (dotted), all with γ l = 0.3D. We see that when the obstacle becomes asymmetric, the ratio becomes fairly linear, as in the case exposed in Fig. 5 is equivalent to the original wave equation and is adapted to a perturbative treatment at low gradients, i.e. |h B | 1. In Sec. III B, we expose the first order perturbative results of this series. We show that the coefficients are mainly governed by complex turning points, corresponding to the locations where two roots of the dispersion relation (11) merge when these are analytically continued in the complex plane. In general, there exist many turning points in the complex plane. Importantly, the ones closest to the real axis dominate while the others contribute as exponentially small corrections. Hence the scattering coefficients are governed by a few dominating contributions, taking the form of complex exponentials of contour integrals from the real line to the complex turning points, see Eqs. (30), and (32).
We then applied these results to a large class of flow profiles, so as to extract the generic features of the scattering coefficients α and β. By studying the behavior of the scattering coefficients as a function of the frequency ω, we distinguish three main regimes. For ultra-low frequencies, ω ω s , |α| 2 and |β| 2 both vanish linearly in ω, see Eq. (31). For intermediate frequencies ω s ω ω min , |α| 2 and |β| 2 share a constant value as in Eq. (39), and when ω/ω min becomes significant, they start drifting apart as shown by Eq. (51). When they do, |β| 2 is generically smaller than |α| 2 . Moreover, we show in Sec. IV B 2 that long obstacles generally produce two dominating complex turning points. If the obstacle is symmetric enough, these two contributions give rise to oscillations in |α| 2 and |β| 2 due to interferences, as illustrated in Figs. 7, and 8. All these features are in perfect agreement with what have been previously observed numerically in [21][22][23] and are complemented by analytic predictions for the parameters governing the various regimes.
In all, this analysis describes in detail what is the "imprint" of Hawking radiation when the flow accelerates but stays subcritical. The physics of the Hawking effect is dictated by horizons, and we have shown here that its imprint in subcritical flows is governed by complex turning points. In this regime, the spectrum becomes more complicated, as it is governed by nonlocal quantities. The study of complex turning points allowed us to provide a simple characterization of this spectrum. When increasing the Froude number, these turning points get closer to the real axis, until they reach it. Before they do, the present treatment breaks down, but it is expected that for increasing F the spectrum will smoothly change from the subcritical one to the Hawking one when F is sufficiently larger than 1. In the Hawking regime, the characteristic length of non-locality becomes smaller than the characteristic length of the gradients [5], and as a result, the spectrum becomes entirely governed by the surface gravity, i.e. the gradient of the Froude number at the horizon. The analytical study of this transition will be the aim of future investigations.
In this appendix, we derive the equation satisfied by the local scattering coefficients defined in Sec. III A. The method we present is a generalization of the Bremmer series [24,25,35]. Whereas the Bremmer series deals with second order differential equations, such as the Schrödinger equation, we consider higher order differential equations [36,37,43]. This is essential to describe dispersive effects of wave propagation, as in our case 8 . Higher order differential operators are also useful to study Schrödinger types of equation in momentum representation, see e.g. [44][45][46]. Our method also bears many similarities with the adiabatic series used in a wide variety of contexts, such as electronic transitions in molecular collisions [39] or particle creation in cosmology [34]. However, here the corresponding operator is not self-adjoint. Under certain conditions, there exists a quadratic conserved quantity, but it has no reason to be positive definite. This is the case for the wave equation of surface waves, see Eq. (24).
To understand the general structure behind this method, we first present it for a general differential equation of degree N , and then apply it to the surface wave equation (5) (App. B). We consider the differential equation Here, the functions f n (x) are assumed to be real while the g n (x) can be complex. To this equation, we associate the corresponding Hamilton-Jacobi equation where we introduced the Hamilton-Jacobi polynomial P HJ . We see that the g n 's do not appear in the Hamilton-Jacobi equation. The reason is that while the f n 's represent the background as perceived by the field φ(x), the g n 's represent the features of the wave equation that are absent of the Hamilton-Jacobi equation. In other words, they encode the possible orderings one can choose when promoting k in (A2) as the operator −i∂ x to obtain (A1). Hence, in the method we shall present, we treat the g n 's as small quantities, as the same order as the gradients of the background, i.e. the f n 's. Going back to (A2), since P HJ is a polynomial of degree N , it has N different roots. The key assumption of the following derivation, is that for all x, the N roots are real and distinct. In particular, no crossing, where one would have k j (x) = k =j (x) for some x, occurs. A common procedure with higher order ODEs is to trade the scalar equation of degree N (A1) for a vectorial equation of degree 1. For this, we gather φ and its derivatives in a column vector . . . 8 Since we adopt here the point of view of dispersive wave equations, our approach has several technical differences with respect to other generalizations of Bremmer series [36,37,43], such as the distinction between the f n and g n , or the adiabatic invariant of λ-canonical systems (see below).
where the −i's are here for future convenience. Eq. (A1) then takes the simple matricial form where (A5) C(x) is the N × N companion matrix associated with the polynomial P HJ . The key idea of the Bremmer approach is to "locally diagonalize C(x)", i.e. at fixed x, and then use the eigen-basis to rewrite the original equation (A1). The characteristic polynomial of C is simply P HJ , and therefore the roots of the Hamilton-Jacobi equation are the eigen-values of C. Hence, we write our field, solution of (A1), as where the phases are the primitive integral of the Hamilton-Jacobi roots, i.e.
Since this introduces N new unknown functions (A j ) j=1..N instead of one, we further impose a similar relation between all derivatives of φ in Φ and the local scattering coefficients A j . We define the A j such that where V is the Vandermonde matrix of the N roots k j , i.e.
Since all the roots are distinct, det(V ) = j<i (k i − k j ) = 0, and the correspondance (A8) between Φ and the A j 's is one-to-one. Physically, Eq. (A8) means that the N −1 first derivatives of φ act as if the A j (x) were constant. We cannot impose this to the N -th derivative, since the field must be a solution of the wave equation (A1). This last condition will instead give us a differential equation satisfied by the A j (x). To obtain it, we first notice that the Vandermonde matrix V is the diagonalizing matrix of the companion matrix C. Indeed, it is rather easy to check that We now have enough material to rewrite the mode equation (A1) in a simple manner. For this, we start by deriving the definition of the local scattering coefficients (A8) This form is still not well suited for a perturbative resolution. What we want is to partially integrate this equation by normalizing the local scattering coefficients. We define the normalized coefficients A j = N j a j . The prefactors N j are chosen in order to get rid of the diagonal elements in the matrix of Eq. (A11). This will recast the equation governing the local scattering coefficients into ∂ x a j e iS j = M · a j e iS j , where the diagonal elements of M are 0. This guarantees that at 0-th order, one recovers the WKB approximation, i.e. the local scattering coefficients are constant. Inserting A j = N j a j in Eq. (A11), we get the condition for the prefactors The matrix element [V −1 ∂ x V + V −1 DV ] jj is real and nonsingular, which guarantees that N j stays real and positive. After integrating this equation, we finally deduce the equation governing the local scattering coefficients

Computation of the M-matrix elements and prefactors N j
To analyze the various coefficients of the M-matrix, we will first explicitly integrate the equation for the prefactors N k in (A13). For this we first need to compute the matrix V −1 . To do so, we introduce N reduced polynomials P j such that P HJ (k) = (k − k j )P j (k), i.e.
Using them, the coefficients of V −1 are easily expressed as The coefficients α j are symmetric polynomials of the N − 1 roots k =j . Their exact expression is rather involved, however, for our present purpose, it is not necessary to write them explicitly. Indeed, using (A17), a direct computation gives where the denotes derivative with respect to k. Moreover, since P HJ (k) = (k − k j )P j (k), the derivatives of P j can be expressed as derivatives of P HJ . In particular, the diagonal term of At this level, it is tempting to identify the right-hand side of this equation as the logarithmic derivative of P HJ (k j ). However, one should not forget that the expression P HJ (k j ) depends on x through both the root k j and the coefficients of P HJ (k), i.e. the f n (x) of Eq. (A2). To go further, we must take into account the contribution of the D-matrix of Eq. (A5). To easily express the coefficients of V −1 DV , we introduce the new polynomial Using it, the combination of (A17) and (A5) shows that Therefore, the prefactor equation (A13) rewrites If the polynomial Q has the good form, this equation directly integrates. Moreover, we recall that unlike P HJ , Q might be complex. While its real part contributes to the amplitude of the prefactor, its imaginary part generates a phase. We call the equation (A1) "λ-canonical", if there exists a function λ(x) > 0 such that Under this assumption, the prefactor equation (A13) directly integrates, and one finds By construction, this generalizes the adiabatic invariant of Sec. II C. The phase shift that appears in (A24) might directly affect the scattering coefficients, e.g., by altering the resonance condition (27). However, in our case this term stays negligible. We now turn to the computation of the matrix elements M j (x). For this, we first express Eq. (A18) in terms of the Hamilton-Jacobi polynomial, and obtain Combining this with the D-matrix, the prefactor N , i.e. Eqs. (A21) and (A24), the expression (A15) becomes where = sign(P HJ (k )). We have also dropped the phase proportional to Im(Q) in (A24) for more clarity, but it is straightforward to add it up.

Scattering coefficients in the smooth limit
At this level, we point out that the equation governing the local scattering coefficients, i.e. (A14), is fully equivalent to the initial wave equation (A1). No approximation has been made so far. However, (A14) gives a simple perturbative method to compute the scattering coefficients in the limit of slowly varying backgrounds. At first order, (A14) shows that the scattering coefficient describing the transition from the mode j to is given by Interestingly, in the smooth limit, that is when f n → 0 and g n → 0, α j→ has a universal behavior. To compute the integral (A27), we make a change of variable S = (k j − k )dx . Since no crossing occurs on the real line (by assumption), ∂ x S = k j (x) − k (x) = 0, and our change of variable is licit. Eq. (A27) becomes This integral can now be evaluated by a residue theorem. All we need to do is to pick up the contribution of all the singularities of the integrand such that Im(S) > 0. α j→ is then given by The main type of singularity would be a zero of k j − k , which is exactly a saddle point S * = S(x * ). Another type of singularity could arise if one of the background functions f n (x) has a pole in the complex plane. This second type of singularity usually gives subdominant contribution, as is the case in all the profiles considered in this paper. We shall thus consider the contributions of saddle points only. For this, one should obtain the corresponding residue. Close to the saddle point x ∼ x * , we have k ∼ k * − δk(x), and k j ∼ k * + δk(x), and hence To obtain this equation, we have neglected the contribution of the Q-term in Eq. (A26), but a (rather nontrivial) computation shows that it produces subdominant corrections in O(|h B |). Then, using the fact that δk ∝ (x − x * ) 1/2 near the saddle point x * , a little algebra shows that Using this, the residue theorem gives us the integral of Eq. (A27), As discussed in the text, in this sum, only the terms with the smallest Im(S * ) contribute significantly, while the other gives exponentially small corrections. Hence the scattering coefficients are usually given by a few contributions, coming from the complex turning points the closest to the real line. We point out here the close similarity between our derivation and the first order result of the Bremmer series [26,38] (in particular Eq. (A30)). In the latter case, π/3 is only the first term of a series, and in the adiabatic limit, it is possible to show that the entire sum is 1. Therefore, the adiabatic limit leads to α j→ ∼ e iS * rather than Eq. (A32). This is also what happens in the adiabatic limit of a time-dependent two-level system [39]. For this reason, it is reasonable to conjecture that this will also be the case here. If this conjecture holds, the replacement π/3 → 1 amounts to a partial resummation of some diagrams of the perturbative resolution of Eq. (A14), as shown in Fig. 9.
. . . Figure 9: We conjecture that the first order expression of the scattering coefficients (A32) can be improved by the replacement π/3 → 1. The conjecture is that this replacement consists in resumming the leading contribution in the smooth limit of the diagrams that involve only the roots k and k j . Once this is done, the higher order perturbative expressions are obtained by summing only over diagrams that involve at least an intermediate state that differs from the initial and final ones.
Appendix B: Application to the surface wave equation (N = 4)

The Bremmer series for surface waves
To apply the preceding results to the problem at hand, we must recast Eq. (5) under the form of Eq. (A1).
From this we directly extract the functions f n (x) and g n (x).
Using Eq. (A23), we verify that our equation is λ-canonical with the function Moreover, the quantity P HJ (k j ) is directly related to the group velocity v j g = (∂ ω k j ) −1 . Indeed, deriving the equation P HJ (k j ) = 0 with respect to ω gives us P HJ (k j )∂ ω k j + ∂ ω P HJ = 0. (B4) We also check that the phase shift due to Im(Q) = 0 is negligible in the regime of interest, i.e. when |h B | 1 is valid. Hence, gathering the results of (B3) and (B4), we finally obtain where Ω is defined after Eq. (1). The expression for the M-matrix is with = sign(Ω(k )v g (k )). From Eq. (B2), we also obtain the expression for Q defined in (A20), which reads To have an estimation of the magnitude of these matrix elements, we compute them in the limit ω ω min . The matrix elements governing the transition k u → k + or k u → k − are given by On the other hand, the matrix element driving the transition k u → k d reads As we see, both are proportional to derivatives of the background, and are small in the smooth limit. But the matrix elements that govern α and β are further suppressed by a factor O(ω 1/2 ) at low frequencies, as we mentioned at the end of Sec. III B.

Scattering coefficients at ultra low frequencies
To obtain the ultra low frequency behavior of the scattering coefficients, we take the limit ω → 0 before the smooth limit |h B | → 0. Indeed, as mentioned in Sec. III B, these two limits do not commute, and the computation of App. A 3 is valid for |h B | → 0 at fixed ω. To correctly obtain the ω → 0 limit, we directly use the low frequency expressions of the matrix elements (B8) in the integral representation for α of Eq. (26) (the computation for β gives the same answer in that limit). The coefficient α is then given by As argued in Sec. IV A 1, in near critical flows 1−F max 1, the variations of 1−F (x) dominate with respect to the other background quantities. Therefore, at leading order in 1 − F max , we have Since we look at the zero frequency limit, in particular we have ω ω min , and hence, the phase of the integrant is evaluated using the zero frequency expressions of the roots, i.e. Eq. (36). Hence the saddle point of Eq. (B11) is the complex horizon defined by Eq. (34). However, one cannot use a residue theorem to compute the integral as was done in Sec. A 3. The reason is that if we make the same change of variable, the point S * = x * (k u (x ) − k + (x ))dx is a branch point and not a pole. Instead, we compute (B11) using a saddle point theorem [47]. This gives From this we deduce We see that in the ultra low frequency limit, α is still governed by the same exponential involving the complex turning point (complex horizon in this regime) as in Eq. (A32), but the prefactor is not 1. Instead, the prefactor vanishes as |C| 2 ∼ ω/ω s . The characteristic frequency ω s is defined from Eq. (B13) and reads We see that ω s /ω min is proportional to a derivative of the background, and therefore is small in the smooth limit. To obtain the simple expression (B14), we applied a saddle point approximation, which is valid if The above calculations and that of App. A 3 show that the prefactor is mainly characterized by two behaviors. When ω ω s , it vanishes as ω/ω s , but for ω s ω ω min , it approaches 1. To reproduce this behavior, we conjectured an explicit form of the prefactor as a function of ω and used it for the numerical plots of Figs. 5(a), 5(b), and 7. This form is The absolute value on ω is here to remind us that the prefactor is the same for α ω and β ω , and does not vary when one applies the relation (40). To conclude this subsection, we wish to underline the fact that this prefactor is an effective way of describing the ultra low frequency behavior of the scattering coefficients. Indeed, were we able to sum over all the singularities of Eq. (A29), we would presumably obtain an expression valid for all values of ω ω min , including the limit ω ω s . It is when we restrict the sum to its dominant contribution that we lose the possibility of taking the limit ω → 0. The "effective" prefactor described in Eqs. (B13) and (B16) allows us to recover the correct ultra low frequency behavior.

The conserved current
In this subsection, we briefly sketch how the conserved current J of Eq. (7) is obtained and how it applies to a superposition of plane waves. This conserved current can be directly obtained by applying the Noether theorem to the action (2) with the symmetry φ → e iλ φ. A slightly quicker way is to start from the conserved norm of Eq. (6). Indeed, the Noether theorem says that the norm density ρ = −Im φ * (∂ t + v∂ x )φ and the current J are related by the conservation law ∂ t ρ + ∂ x J = 0. (B17) By computing ∂ t ρ and using the wave equation (4), we deduce the current. The calculation is as follows: When we apply ∂ t ρ + ∂ x J = 0 to stationary solutions φ(t, x) = Re(φ ω (x)e −iωt ), we see that the current (7) is x-independent. We now want to apply this current to a local superposition of WKB waves as in Eq. (16) so as to obtain Eq. (24). Because J involves only derivatives of φ up to third order, the ansatz of Eqs. (16), (17a), (17b), and (17c) shows that the computation for a local WKB superposition is the same as for exact plane waves, i.e. when v, c, and h B are constant. Moreover, since J is a quadratic quantity in the field, it is enough to show this for a superposition of 2 plane waves φ = A 1 e ik 1 x + A 2 e ik 2 x , where k 1 and k 2 are solutions of the dispersion relation (8). Injecting the above form in (7), we see that where The factor Ω(k)v g (k) in the diagonal terms directly follows from Eq. (B4). However, it is much more delicate to show that the cross term J × is exactly 0. Since J is by construction independent of x, it has to vanish. One could presumably stop here, invoking the latter argument, but it would be instructive to understand why the first factor of Eq. (B21) is always 0. It follows from the fact that k 1 and k 2 are distinct solutions of the dispersion relation (8). We have proven this identity, but its derivation is rather involved and it is unclear what its physical interpretation is, or how one could generalize it. We present it in the form of the following lemma: Lemma: Let P be the polynomial P (k) = k 4 + ak 2 + bk + c.

Contour integral for a monotonic profile
In this section, we present the computation of the contour integral necessary to obtain the value α ω→0 in the profile (53), that is, Eq. (55). Interestingly, this computation is very similar to what happens for the Schrödinger equation in a tanh potential [26]. The integral leading to Eq. (55) is defined as where we used 1 − F max 1, and we have chosen x 0 = x R for convenience. Before computing the above integral, we rewrite the single step profile of Eq. (53) as We compute the contour integral above with the parametrization x (t) = x R + it∆ 0 , with x R and ∆ 0 given by Eq. (54). Using the identity we get We now interpret this integral as another contour integral, clockwise along the lower half unit circle (see Fig. 10). For this we define z = e −iπt , which implies dz/z = −iπdt, and rewrite (B32) as We deform this contour into the segment ]-1;1[, right below the pole at z = 0, which gives Finally, using the identity (z − i ) −1 = Pz −1 + iπδ(z), we get The second term of this equation being real, we have