Acoustic black hole analogy to analyze nonlinear acoustic wave dynamics in accelerating flow fields Open Access

The nonlinear behavior of acoustic waves emitted or scattered by accelerating boundaries has recently attracted attention as nonlinear amplitude modulations and frequency side-bands can reveal distinct information about the motion of a wave emitter or scatterer. A comprehensive review of different (potential) applications is given in the work of Gasperini et al.,¹ who emphasize the difficulty of predicting nonlinear Doppler shifts and propose a computational method based on a suitable coordinate transformation that can account for the accelerating motion of wave scattering boundaries. Using a similar approach, we recently confirmed the analytical predictions by Christov and Christov² regarding the nonlinear Doppler waveform of an oscillating emitter by the means of numerical simulations.³ However, the background medium is assumed to be at rest in the aforementioned studies, whereas it is known that Doppler shifts are also caused by the motion of the background medium itself. This is harnessed, for instance, in medical Doppler ultrasonography to estimate the velocity of blood flows or moving tissue and to identify pathological patterns in the Doppler spectrum.^4,5 While considerable research efforts have been made to improve Doppler signal processing techniques for medical applications,⁶ a more fundamental understanding of nonlinear Doppler modulations in moving flow fields is still missing.

In the present work, we derive a physical model and present a computational method to simulate acoustic waves in moving background flow fields. We then employ the acoustic black hole analogy in order to investigate the dynamics of an acoustic wave propagating in an accelerating background flow field.

An acoustic black hole is a convergent fluid flow with sonic transition, where the point of sonic transition forms an acoustic event horizon from which small signal sound waves cannot escape.⁷ Unruh⁸ discovered that the propagation of sound waves in an irrotational and inviscid fluid flow, and in particular, the associated acoustic metric, corresponds exactly to the motion of a scalar field in a curved (3 + 1)-dimensional Lorentzian space–time geometry.⁹ Since then, the acoustic black hole analogy and other back hole analogue systems have gained considerable attention because the kinematic connection to curved space–time geometry allows for the investigation of certain kinematic aspects of gravitational black holes, such as Hawking radiation^8,10–13 or rotational super-radiance,^14,15 in laboratory experiments.^16,17 Apart from the investigation of astrophysical phenomena, the study of acoustic black holes may also contribute to an improved understanding of the acoustic wave behavior in accelerating background flow fields, which may, in turn, help to improve relevant technical applications such as medical Doppler ultrasonography.

In principle, the presence of a sonic transition is not essential to the wave dynamics discussed in this work, which merely depend on the flow acceleration. However, as it is further discussed in Sec. II A, a steady-state sonic horizon has the beneficial feature that it provides a fixed reference point of the propagating waveform that is subject to a constant background flow acceleration, even if the background flow field is non-uniform in space. This allows to parameterize the wave-flow configuration by combining a suitable Helmholtz number with the acceleration at the sonic horizon, which again depends on the size of the steady-state acoustic black hole. We will then make use of this parametrization to study the amplitude modulations along geometrically similar wave characteristics.

Starting from the work of Unruh,⁸ Visser,⁹ and Balbinot et al.,¹⁸ we first derive a quasi-one-dimensional wave equation for a moving and compressible background medium based on first-order perturbation theory (see Sec. II B). As we aim to study the dynamics of acoustic waves in a smooth background flow field, it is briefly motivated in Sec. II C how, under certain conditions, a smooth sonic transition is supported in de Laval nozzle flows¹⁹ and Bondi–Parker-type flows,^20,21 where the latter describe the mass accretion and plasma wind flows around compact stellar objects.

In Sec. III, we derive a suitable coordinate transformation in physical space, conveyed by a generic change of variables, which is instrumental in maintaining the well-posedness of the wave equation and used as basis for a numerical solution. While this approach is well-established in the context of wave emitters moving relative to a quiescent background medium,^1–3 the well-posedness of the governing equation is still a particular problem when the wave emitter approaches or even exceeds the sonic speed.² Despite the fact that we consider a wave emitter of zero velocity relative to the moving background medium, so that no sonic shock can form, the governing wave equation becomes ill-posed for supersonic flows if the wave equation is formulated in a fixed Eulerian reference frame, which is passed by the supersonic flow. Therefore, a space- and time-dependent coordinate transformation recovers the proper two-sided wave propagation behavior at any point in physical space. This idea evolves along similar lines as the utilization of Galilean and Lorentz invariants of the Navier–Stokes equations and the wave equation as presented by Gregory et al.²² and Wang.²³ Gregory et al.²² pointed out that the description of Doppler phenomena in moving fluids can greatly simplify under the appropriate coordinate transformation.

Finally, we consider the canonical case of a spherically symmetric steady-state acoustic black hole, for which quasi-planar and linear wave propagation is obtained in the limit of very large radii of the sonic horizon. The focus is on the geometrical similarity of the wave characteristics as described by the Helmholtz number, and the effect of the sonic horizon's surface gravity, that is, the background flow acceleration, on the acoustic wave dynamics.

In the following, we discuss important parameters that govern the acoustic wave dynamics in a stationary spherical flow with sonic transition. Figure 1 (left) shows the schematic of a spherical sonic horizon at radial position $r=rh$⁠, where the magnitude of the inward directed radial flow velocity u coincides with the sonic speed c. In the context of the present work, R and $Ṙ$ are the instantaneous position and velocity, respectively, of a virtual spherical wave emitting boundary moving along with the flow as indicated in Fig. 1 (left). Thus, there is no relative motion between the background medium and the wave emitter, so that any nonlinear Doppler shifts are exclusively imparted by the non-uniform motion of the background medium. For the following derivations, it is also advantageous to write the radial velocity distribution as

The parameter $n(r,t)$ in Eq. (1) can formally be both time- and space-dependent, and it can be specified in such a way that the continuity equation for a compressible spherically symmetric flow is satisfied. In the present work, we perform numerical solutions for the particular case of a stationary spherically symmetric flow, for which the continuity equation $4πr2ρ(r)u(r)=4πR2ρ(R)Ṙ$ yields

which is equivalent to Eq. (1) for

Hence, we have n = 2 for an incompressible flow and $n≤2$ for a compressible flow with monotonically increasing density toward the center.

For a given stationary velocity field $u(r)$⁠, R may be obtained as solution of Eq. (1) subject to initial conditions for R and $Ṙ$⁠. In the context of the present work, this solution is required as a boundary condition for the governing wave equation and must, therefore, be obtained analytically. In order to facilitate the analytical solution of Eq. (1), we take n as spatial-temporally constant. This assumption is justified as we are only interested in a confined region around the sonic horizon and the principle effect of the magnitude of n. Furthermore, n may also be chosen in such a way that Eq. (1) represents a non-spherical flow geometry, as for instance the quasi-one-dimensional flow in a de Laval nozzle. This allows us to apply the following considerations to any quasi-one-dimensional flow with sonic transition.

With $n=const$⁠, it can be inferred from the work of Visser⁹ that the acoustic surface gravity $gh$ at the sonic horizon of a steady-state acoustic black hole associated with a quasi-one-dimensional flow follows as

The term “surface gravity” reflects the analogy to the acceleration experienced by a freely falling observer at the event horizon of a gravitational Schwarzschild black hole.²⁴ For a nearly constant speed of sound in the proximity of the steady-state sonic horizon, $gh$ represents the acceleration of the stationary background flow field. The change of the fluid acceleration with respect to r is proportional to the acceleration itself and determines how quickly neighboring wave characteristics diverge away from each other as indicated in Fig. 1 (right). This suggests that the acoustic surface gravity is an important parameter that governs the rate at which an outward propagating wave is red-shifted.

The variation of the background flow field over a waveform further depends on the wavelength itself. Therefore, we adopt the Helmholtz number²⁵ to characterize the size of the acoustic black hole relative to the emitted wavelength. With $fa$ being the emitted wave frequency and $λa=c/fa$ the emitted wavelength, the Helmholtz number is given by $He=2πrh/λa$⁠. Substituting $rh$ from Eq. (4) and dropping the dimensionless parameter n as well as the constant $2π$⁠, the Helmholtz number can be rewritten as

meaning that the ratio $c2/λa$ is uniquely defined by the values of $He, gh$⁠, and $c∂c/∂r|rh$⁠, where the latter contribution is negligible if c varies slowly over the acoustic wavelength. This relation suggests a geometrical similarity of the wave characteristics for a particular value of He and a given flow field as parameterized by n. A particular value of He can be obtained from different combinations of $gh$ and $λa$⁠. As the Helmholtz number depends, in this particular form, on the flow acceleration at a distinct reference point rather than the radius of a spherical wave emitter, its applicability is not restricted to spherical symmetry, but extends to any quasi-one-dimensional geometry, where the corresponding flow field can be parameterized by the scaling factor n in Eq. (1).

Because it is of particular relevance to the present work, the derivation of a wave equation describing the propagation of acoustic waves in a moving fluid is briefly repeated here, drawing on the works of Unruh,⁸ Visser,⁹ and Balbinot et al.¹⁸ We start from the three-dimensional representation, which is then simplified to the case of spherical symmetry at the end of this section. Neglecting viscous stresses, the dynamics of the fluid motion and the acoustic waves are governed by the Euler equations, that is, the conservation equations of mass and momentum, given by

In Eqs. (6) and (7), u, ρ, and p are the fluid velocity, density, and pressure, respectively, and ψ is the gravitational potential. Further assuming the flow to be irrotational (⁠$∇×u=0$⁠), a velocity potential $ϕ$ exists so that

Integration of Eq. (7) and using Eq. (8) gives the transient Bernoulli equation

where $h(p)$ is the difference of the fluid's enthalpy between the total pressure p and some undisturbed reference pressure p₀. By introducing the first-order perturbations $ϕ=ϕ0+εϕ1, p=p0+εp1, ρ=ρ0+ερ1$⁠, and $u=u0+εu1$ about the background quantities,⁹ where the subscript 0 indicates the background quantities and the terms scaled by ε are small amplitude acoustic perturbations, substituting into Eqs. (6) and (7), and further linearizing the enthalpy difference according to $h(p0+εp)≈h0+εp/ρ0$⁠, we obtain

Solving Eq. (13) for p₁ gives

where

is the material or Lagrangian derivative operator. Inserting Eq. (14) into the linearized barotropic relation $dp/dρ=c2≈p1/ρ1$ of the acoustic perturbation and then substituting for ρ₁ in Eq. (11), one arrives at the wave equation⁹ 

Similar derivations based on first-order perturbations are found in the works of Ewert and Schröder²⁶ and Godin,²⁷ where Ewert and Schröder²⁶ derived an additional term that represents the effect of rotational velocity components in addition to the irrotational flow field. In spherical coordinates and under the assumption of spherical symmetry, Eq. (11) is rewritten as $∂ρ1/∂t+(1/r2)∂[r2(ρ0u1+ρ1u0)]/∂r=0$⁠, where r and u₀ are the radial coordinate and the radial background flow velocity, respectively. Repeating the above steps analogously, the wave equation becomes

where $D/Dt=∂/∂t+u0∂/∂r$ is the Lagrangian/material derivative operator analogously to Eq. (15). Expanding the derivative operators in Eq. (17) gives

As the focus of the present work is on the effect of the background flow acceleration on the acoustic wave dynamics in a continuous flow field, we make the assumption of a stationary and smooth (shock-free) sonic transition. In this section, we briefly indicate under which conditions a smooth sonic transition is physically possible. In principle, a smooth sonic transition can be established in de Laval nozzles¹⁹ or in spherically symmetric Bondi–Parker-type flows.^20,21 As pointed out by Williams and Dyson¹⁹ and briefly outlined in the following, the prototypical Bondi–Parker equations are very similar to the equations describing the quasi-one-dimensional steady-state flow in a de Laval nozzle.

We assume that the barotropic relation $dp/dρ=c2$ holds for the background fluid state, which is indeed justified under the assumption of a shock-free flow. With the surface area $A(r)$ passed by the flow in the normal direction, a stationary, inviscid, and quasi-one-dimensional flow with mass flow rate $ρ0u0A$ and gravitational potential ψ is governed by the continuity and momentum equations^28,29

Equations (19) and (20) can be combined to give the radial velocity gradient

A singular velocity gradient (shock) can only be avoided if the parenthesized term on the right-hand side of Eq. (21) vanishes. The gradient of the gravitational potential in a Bondi–Parker-type flow is given by²⁹ 

where G is the gravitational constant and M the mass of a compact stellar object. The critical/sonic point is at $rh=GM/(2c2)$⁠,^29–31 while $(∂A/∂r)/A=2/r$ for a spherical flow, so that a smooth sonic transition at $r=rh$ is indeed supported by the above equations.¹⁹ 

The stationary quasi-one-dimensional flow in a de Laval nozzle is governed by the same set of equations, however, with $∂ψ/∂r=0$⁠. Here, the parenthesized term vanishes at the point of sonic transition under the appropriate boundary conditions since $∂A/∂r=0$ at the nozzle throat, effectively resulting in non-convergent/non-divergent streamlines.¹⁹ Kovalenko and Eremin³² demonstrate the stability of the spherically symmetric transonic Bondi-type flow against spherically symmetric acoustic perturbations. They further suggest that the breakdown of spherical symmetry is responsible for the onset of instability and shock formation. By linear stability analysis, Ray and Bhattacharjee³¹ show that while the Bondi–Parker equations admit a stationary transonic solution associated with the state of lowest energy of the system,²⁰ this solution is physically not realizable. Ray and Bhattacharjee³¹ further show that a transonic solution is only supported by linear stability analysis if the flow is transient. However, while the radial velocity profile becomes time-dependent under the assumption of a transient flow, the instantaneous velocity profiles are still qualitatively similar to the solution of the stationary Bondi–Parker equations.³¹ Therefore, one may assume that the timescale of the variation of the background flow velocity u₀ is negligibly small as compared to the timescale of the acoustic wave problem. Maicke et al.³³ discuss various distinct three-dimensional effects leading to the departure from the shock-free condition in a supersonic de Laval nozzle flow. Moreover, the realization of nearly quasi-one-dimensional conditions, precluding the occurrence of perturbations of the background flow field, requires superfluidic conditions as can be achieved with Bose–Einstein condensates^34,35 or superfluidic helium II.³⁶ 

As is common in the context of acoustic black holes,^7,8,18,37 we take the speed of sound as a constant in our analysis. This simplifying assumption is associated with a first-order approximation of the background flow compressibility, as it is shown by expanding $p0(ρ0)$ around some reference density $ρ0,ref$ under the assumption of a barotropic flow. With $dp0/dρ0=c2$⁠, the Taylor series expansion yields

showing that the variation of c with respect to ρ₀ is a second-order effect. With $c=const$⁠, the second term on the right-hand side of Eq. (18) vanishes, $∂c/∂r=0$ in Eq. (4), and Eq. (5) reduces to

Again employing the critical point $rh=GM/(2c2)$ of a spherical Bondi–Parker-type flow,^29–31 as well as Eq. (4) for $c=const$⁠, Eq. (22) is conveniently expressed as

Equation (20) readily gives the density gradient

Substituting Eq. (26) into Eq. (18) yields

We will use Eqs. (24), (25), and (27) in the further course of this study. A general advantage of the above formulation is that the effect of a variable and non-uniform background density field is absorbed in the right-hand side term that exclusively depends on the background flow field. The appropriate velocity distribution can then be treated as an input to the coefficients of the wave equation, and it can be obtained either analytically, numerically, or experimentally. However, the coefficient $c2−u02$ of the Laplacian $∂2ϕ1/∂r2$ in Eq. (27) changes its sign at the point of sonic transition so that the problem is ill-posed for regions where $u02>c2$⁠. In Sec. III, we will introduce a time- and space-dependent coordinate transformation in physical space in order to extend the validity range of Eq. (27) to supersonic background flow fields.

Regarding the numerical procedure and the uniqueness of the solution, it is a particular problem that the coefficient of the Laplacian on the left-hand side of Eq. (27) switches its sign at the point of sonic transition, rendering the problem ill-posed. The same problem is pointed out by Christov and Christov,² however, for the wave equation describing a moving wave emitter approaching the sonic shock condition in a quiescent background medium. In the present work, the acoustic wave does not form a shock as the wave emitter motion is aligned with the background flow, and the problem merely originates in the circumstance that the wave equation is formulated in a fixed Eulerian reference frame, which the fluid passes at the speed of sound at the sonic horizon. Hence, at the point of sonic transition, an outward propagating wave “freezes” and appears to have a zero frequency (infinite wavelength). In a reference frame moving with the background flow, by contrast, proper two-sided wave propagation at finite wavelength is observed. This behavior can be recovered by a suitable Galilean-type coordinate transformation, which emphasizes the fact that the well-posedness is not primarily a property of the governing equation itself, but rather of the reference frame in which it is formulated.

In order to switch to a more suitable reference frame to describe the wave dynamics for transonic background flows, we introduce the moving physical domain $Ω(t)=[R(t),r∞]$ with coordinates (r, t), moving wave emitting boundary $R(t)$⁠, and fixed boundary $r∞$⁠, and the fixed computational domain $Θ=[XR,X∞]$ with coordinates $(ξ(r,t),t)$ and fixed boundaries $XR$ and $X∞$⁠. Next, we invoke the time-dependent coordinate transformation $r: [XR,X∞]→[R(t),r∞],(ξ,t)↦r(ξ,t)$ to implement the generic change of variables,

A similar approach for a wave emitter moving relative to a quiescent background medium is presented by Gasperini et al.¹ The derivative operators conveying the change of variables follow as

Applying the derivative operators to Eqs. (28) and (29) and substituting into Eq. (27), the transformed wave equation,

is obtained, where the spherical Laplacian gives rise to the additional contribution $(2/r)∂Φ1/∂ξ$⁠. Only the derivatives of $Φ1$ need to be discretized. The derivatives of ξ, ψ, and $U0$ depend on the specification of the coordinate transformation and the background flow field and can, therefore, be seen as space- and time-dependent coefficients.

In particular, the coordinate transformation can be chosen in such a way that a non-negative coefficient of the Laplacian $∂2Φ1/∂ξ2$ is obtained, thereby guaranteeing the well-posedness of the problem. It follows from Eq. (33) that the Laplacian coefficient is non-negative if the condition

is satisfied. The term $∂ξ/∂t$ is the velocity of a fixed point in Θ as viewed from a moving point in $Ω(t)$⁠.

We opt for a linear mapping between the moving physical domain $Ω(t)$ and the fixed computational domain Θ as described by the linear function

and choose $XR=0$ and $X∞=1$⁠. One is, in principle, free in the choice of $XR$ and $X∞$⁠. However, we will choose the sonic horizon radius on the order of unity and scale the emitted wavelength accordingly, so that the Jacobian of the coordinate transformation $∂ξ/∂r$ is on the order of unity with the above choices of $XR$ and $X∞$⁠. The derivatives of ξ follow as

From the above relations and Eqs. (1) and (25), the derivatives of ψ and $U$ follow as

Finally, R, $Ṙ$⁠, and $R¨$ need to be specified. For a stationary spherical flow as described by Eq. (1), these quantities follow by solving Eq. (1) for R under the boundary condition $u(rh)=const=−c$ (steady-state flow) and the initial condition $R(t=0)=R0$⁠, which gives

where it is assumed that $n=const$⁠.

An explicit finite difference time domain (FDTD) method based on central differences in space is used to solve the transformed wave equation [Eq. (33)] in the fixed computational domain Θ. Without special treatment, the explicit FDTD solution may be contaminated by dispersive numerical noise,^38,39 which can be amplified by the frequency side-bands that are generated by the accelerating background flow field. In order to counteract the development and amplification of dispersive numerical noise, we opt for sixth-order accurate central differences in space, where this relatively high order helps to delay the onset of dispersive numerical noise.⁴⁰ The ξ-discrete solution is explicitly advanced in time using the predictor–corrector method by Dey and Dey.⁴¹ Originally developed for ordinary differential equations, the predictor–corrector method was later applied to the linear wave equation in the context of seismic wave propagation by Nascimento and Pestana,⁴² who showed that the method effectively mimics an attenuation term weighted by a factor γ. This term predominantly acts on the higher wave frequencies so that the diffusive effect on the numerical solution vanishes with increasing grid resolution while preserving the stabilizing effect. In the present work, we use the recommended value $γ=0.5$⁠, for which the predictor–corrector method evolves into a second-order Runge–Kutta scheme when applied to an ordinary differential equation.⁴³ However, on application to Westervelt's wave equation,⁴⁴ we previously found a convergence rate between first and second orders, which is mainly attributed to the treatment of the mixed derivatives as source terms as described in our previous work,³ in which a more detailed description of the numerical method and rigorous convergence studies are provided. For a detailed discussion of the stability properties of the predictor–corrector method, we refer to the work of Dey.⁴³ 

We consider a spherically symmetric steady-state acoustic black hole for a Bondi-type (inward directed) flow as described by Eq. (27), with the gradient of the gravitational potential given by Eq. (25). The input flow field is given by Eq. (1), where the boundary motion is specified by Eqs. (44) and (45). We choose n = 2 (incompressible flow field) for our base configuration, but we also present results for n < 2.

The acoustic waves are emitted from a virtual spherical boundary with instantaneous radius $R(t)$ collapsing along with the spherically symmetric fluid flow (see Fig. 1). We arbitrary choose a sound speed of $c=1500 m/s$⁠. Given the dispersion relation $c=λafa$⁠, the values of He [see Eq. (24)] and $gh$ [see Eq. (4) for $∂c/∂r=0$] uniquely define the emitted wave frequency $fa$ and the background flow velocity distribution. The simulations comprise a dimensionless time of $fat=27.5$ emitted wave periods. The sonic horizon position $rh$ is implicitly specified by Eq. (1), where $u0(rh)=−c$⁠. The initial wave emitter position $R0>rh$ (the wave emitter is “falling” into the acoustic black hole) is chosen in such a way that the sonic horizon approximately intersects the 14th wave period for a slowly varying background velocity u₀. The position $r∞$ of the far field boundary is chosen in such a way that the waves do not reach the boundary. Apart from the grid convergence study in Sec. IV B, the grid resolution is $Nppw=400$⁠, where $Nppw$ is the number of grid points per emitted wavelength $λa$⁠. The Courant–Friedrich–Lewy number is always $CFL=cΔt/Δr=0.1$⁠, where $Δt$ is the time step size and $Δr$ the spatial step size as specified by $Nppw$ and $λa$⁠.

Figure 2 depicts an instantaneous acoustic waveform intersected by the sonic horizon at $r/rh=1$ for five different grid resolutions. The acoustic perturbation potential $ϕ1$ is normalized by the excitation amplitude $Δϕ1,a$⁠. The acoustic surface gravity is $gh=7.5×106 m/s$⁠, and the Helmholtz number as given by Eq. (24) is $He=20$⁠. This configuration involves the largest surface gravity and the smallest Helmholtz number considered in this study. Hence, it involves the largest gradients of the background flow velocity relative to the acoustic wavelength and constitutes, therefore, the most challenging case for the numerical method. The waveform in Fig. 2 is recorded at $fat=27.5$⁠. The figure on the right shows a detail of the wave profile close to the outer domain boundary $r∞$⁠, which corresponds to the longest travel distance of the wave. It can be seen that for $Nppw=400$⁠, a sufficiently well-converged solution is obtained. This required resolution agrees with the findings of our previous work,³ in which a more rigorous analysis of the accuracy of the numerical method is presented.

Figure 3 shows the distribution of the acoustic perturbation potential $ϕ1$ in the space–time plane. The spatial axis is normalized by the sonic horizon radius $rh$ and the time axis by the emitted wave frequency $fa$⁠. The distribution is depicted for different values of the surface gravity $gh$ [see Eq. (4) for $∂c/∂r=0$] and the Helmholtz number He [see Eq. (24)]. Starting from subfigure (a), which represents the smallest acoustic black hole relative to the fixed emitted wavelength, the distributions in (b)–(d) are obtained by increasing the sonic horizon radius $rh$ while leaving the other parameters unchanged, so that $gh$ decreases and He increases accordingly. Initially, the wave emitting boundary with instantaneous radius $R(t)$ is located outside the sonic horizon as indicated by the black solid line. Subsequently, the wave emitting boundary collapses and eventually passes the sonic horizon. It can be seen that the wave characteristics diverge away from the sonic horizon over time, which corresponds to a gradual red-shift of the outgoing wave train.

Figure 4 shows the instantaneous wave profiles corresponding to the configurations in Fig. 3, recorded at $fat=27.5$⁠. For reference, the black dashed line indicates a $1/r$-decay of the waveform, starting from the sonic horizon. The reference value for the wave amplitude at the sonic horizon is obtained by linear interpolation between the two neighboring wave peaks. The general observation is that the entire wave is increasingly amplified with increasing magnitude of $gh$⁠. In all cases, the rate of amplification is most pronounced in the innermost region close to the wave emitting boundary so that the maximum wave amplitude is observed at some distance away from the location of wave emission, in contrast to the $1/r$-decay dictated by the spherically symmetric geometry. It appears that the dimensionless distance $r/rh$ of the maximum wave amplitude increases with a decreasing magnitude of $gh$⁠.

Figure 5 depicts the wave amplitude modulation at the sonic horizon over time, where the amplitude modulation factor is defined as $AM(rh,t)=ϕ1(rh,t−th)/ϕ1(rh,th)$ for $t≥th$⁠. The time $th$ marks the time instant at which the wave emitting boundary with instantaneous position $R(t)$ passes the sonic horizon at $rh$⁠. It can be seen that the growth rate of AM increases with $gh$⁠. Again, this illustrates the overall amplification of the waveform imparted by the acceleration of the background flow. Even though we are lacking an analytical expression for AM at this time, it is interesting to note that the curvatures of the AM lines appear to be qualitatively related to the curvatures of the corresponding wave emitter trajectories in Fig. 3.

For larger values of $gh$ in Fig. 3, the wave characteristics possess a distinct curvature, whereas they flatten out as $gh$ decreases. At the same time, the characteristics tend to evolve more and more parallel to the sonic horizon line, which means that for very small $gh$ and very large He, the entire wave train degenerates into a standing wave as viewed from the reference frame of the sonic horizon. This is demonstrated by the space–time diagram of $ϕ1$ in Fig. 6(a), in which the sonic horizon $rh$ is further increased while again leaving the other parameters unchanged. Figure 6(b) depicts the corresponding instantaneous wave profile at $fat=27.5$⁠, indicating that the waveform becomes planar and exhibits increasingly linear behavior as $rh$ increases. The increasingly linear behavior manifests in a reduction of the frequency side-bands and amplitude modulations.

As discussed in Sec. II A, the same Helmholtz number He [see Eq. (24)] can be obtained for different combinations of the emitted wavelength $λa$ and the surface gravity $gh$⁠. Figure 7 shows such a variation for the smallest (a) and the largest (d) acoustic black hole in Fig. 3, which are represented by subfigures (a) and (b) in Fig. 7. Figures 7(c) and 7(d) are based on the same Helmholtz numbers as (a) and (b), respectively. However, the surface gravity in (c) is equal to the surface gravity in (b), and the surface gravity in (d) is equal to the surface gravity in (a), which means that in both (c) and (d), the emitted wavelength is adjusted to obtain the corresponding values for He. It follows from Eq. (24) that the required scaling factor of $λa$ is the inverse of the surface gravity ratios, which is 4 in this particular case. The emitted wave frequency changes according to $fa=c/λa$⁠. The comparison between Figs. 7(a) and 7(c) and Figs. 7(b) and 7(d), respectively, shows that for identical values of He, the wave characteristics evolve along virtually identical trajectories, however in a space–time plane normalized by $rh$ and $fa$⁠. This means that the characteristics corresponding to the same Helmholtz number He for different combinations of $gh$ and $λa$ are not identical in physical space–time, albeit geometrically similar.

However, as shown by the instantaneous wave profiles in Fig. 8, the radial distribution of the acoustic perturbation $ϕ1$ is different for identical values of He but different values of $gh$⁠. While the difference is negligible for $He=80$⁠, it becomes more pronounced for $He=20$⁠. Furthermore, Fig. 8 echoes our previous observations; a larger surface gravity induces a stronger amplification of the waveform and is accompanied by a shift toward closer radii $r/rh$ of the wave amplitude maximum.

In order to investigate the effect of the flow field scaling parameter n in Eq. (1), we start from the case in Fig. 3(a), where n = 2, and reduce the scaling parameter to n = 0.5 while leaving the remaining parameters unchanged. The magnitude of the radial flow velocity then increases toward the center with $r−1/2$ [see Eq. (1)], which corresponds to the free-fall condition in a Bondi-type flow.⁴⁵ Figure 9 shows both the space–time distributions of $ϕ1$ and the instantaneous wave profiles at $fat=27.5$ for both configurations. A comparison of (a) and (c) shows that a reduction of n reduces the curvature of the wave characteristics, which is due to the reduction of the background flow velocity gradients. Equation (4) implies that the surface gravity $gh$ reduces linearly with decreasing n. A comparison of (b) and (d) shows that the smaller magnitude of $gh$ for n = 0.5 manifests in a less pronounced amplification of the waveform, which is in agreement with the previous observations.

Finally, it is noted that we do not observe geometrical similarity for identical Helmholtz numbers and different values of the flow field parameter n. This can be seen by comparing Fig. 9(c) with Fig. 3(d). Both figures are based on the same Helmholtz number and also the same surface gravity, but the latter exhibits a significantly smaller expansion of the wave characteristics around the sonic horizon region. Hence, it appears that geometrical similarity of the wave characteristics strongly relies on a given background flow field.

The acoustic black hole analogy employed to analyze acoustic waves propagating in non-uniform background flow fields reveals several interesting features. Most importantly, the acoustic black hole analogy allows for a parametrization of the wave-flow configuration in terms of a Helmholtz number He that is expressed as a function of the acoustic black hole's surface gravity [see Eq. (24)]. The numerical simulations demonstrate that He describes geometrically similar wave characteristics, which can, however, be subject to different amplitude modulations depending on the magnitude of the surface gravity. In particular, an increase in the surface gravity and, hence, the overall magnitude of the background flow acceleration cause an increasing amplification of the entire waveform, which again becomes increasingly noticeable over time as the wave propagates through the accelerating flow field. Consequently, the maximum wave amplitude occurs at some distance away from the wave emitting boundary. A qualitatively similar effect is reported by Christov and Christov² for a wave emitter accelerating through a quiescent background medium in a one-dimensional Cartesian case, where the nonlinear Doppler amplification becomes more noticeable with increasing distance from the source. In principle, the above observations are expected to apply to subsonic flows in a very similar fashion as the amplitude modulations merely depend on the acceleration of the flow field and not on the magnitude of the flow velocity. In fact, Fig. 6 shows that an acoustic wave exhibits quasi-linear behavior across the sonic horizon if the surface gravity is small enough and He is large. Loosely speaking, this scenario is associated with very large acoustic black holes, having in mind that the acoustic black hole's size should be assessed in conjunction with the acoustic wavelength.

In addition to the merits regarding the study of nonlinear wave dynamics in accelerating flow fields, the acoustic black hole analogy may also prove useful for the verification and benchmarking of numerical methods for acoustic waves in moving fluids. On the one hand, the system is straightforward to describe geometrically, with the sonic horizon as a distinct reference point, and it is well defined in terms of a small number of parameters. On the other hand, it is rich in phenomena and complexity, which can be increased arbitrary by considering a non-stationary sonic horizon, circumferential flow components, and other three-dimensional effects, nonlinear wave steepening, or discontinuities in the background flow field.

A general advantage of the presented model formulation is that the effects of a compressible background medium are absorbed in terms that exclusively depend on the stationary background flow velocity [see Eq. (27)]. Hence, the problem reduces to the knowledge of the appropriate background flow field that satisfies the continuity equation for a compressible flow of a barotropic fluid. In particular, we can study the principle effect of a space-dependent distribution of ρ₀ by varying the background velocity field using the parameter n in Eq. (1). While n = 2 corresponds to the incompressible flow assumption, n <2 is expected in the region of sonic transition as it can be inferred from Eq. (3) and the works of Treves et al.⁴⁵ and Ray and Bhattacharjee.³¹ As shown by Fig. 9(a), a decrease in n reduces the curvature of the wave characteristics as a result of the smaller velocity gradients. Equation (1) implies that a reduction of n moves the sonic horizon toward larger radii for given values of R and $Ṙ$⁠. The reduced curvature of the wave characteristics manifests in a less pronounced amplification of the waveform [see Fig. 9(b)]. We may further anticipate the qualitative effect of a variable speed of sound, which is a second-order effect as discussed in Sec. II C. Due to the compression of the fluid toward the center point, the speed of sound should increase toward the center accordingly. Taking the speed of sound at the sonic horizon radius $rh$ as a reference, c increases toward the inner region and decreases toward the outer region. For an outward propagating wave, this has the qualitative effect that the characteristics in the inner region evolve slower toward the center point, whereas the characteristics in the outer region escape at a lower speed. Therefore, the variable speed of sound is expected to straighten the wave characteristics, which is qualitatively similar to the effect of the variable density field. It is further noted that a variable speed of sound also alters the surface gravity [see Eq. (4)].

The transformation of the governing wave equation for a moving background medium conveyed by a generic change of variables is shown to be instrumental to maintain the well-posedness of the governing wave equation in the presence of a transonic flow field. In the present work, we opt for a time-dependent linear transformation in physical space. The linearity of the transformation facilitates the specification of the involved derivatives as a function of the boundary motion and the motion of the background flow field [see Eqs. (36)–(43)]. However, the linearity of the transformation also imposes a constraint on the possible combinations of the acoustic black hole size and the size of the moving/deforming physical domain $Ω(t)$⁠. This is because the speed of any point in $Ω(t)$ relative to the moving background flow must remain below the sonic speed in order to preserve the well-posedness of the problem as discussed in Sec. III A. Drawing on the method of characteristics, this constraint can be relaxed by specifying the transformation in such a way that any point of the moving physical domain $Ω(t)$ moves exactly with the background fluid. For a uniform background flow field, the convective wave equation [Eq. (27) with zero right-hand side] would then reduce to the standard wave equation under this particular Galilean transformation as it readily follows from the transformed wave equation [Eq. (33)]. Any potential non-uniformity of the flow field generates additional terms, which are then exclusively related to nonlinear Doppler modulations, that is, the generation of frequency side-bands and amplitude modulations. Thus, the appropriate coordinate transformation is expected to greatly facilitate the study of Doppler shifts and nonlinear amplitude modulations in complicated flow situations, as has also been concluded by Gregory et al.²² For instance, the method of characteristics may permit full analytical solutions for certain flow situations, in which separations of length and/or time scales can be identified. An example is found in the work by Christov and Christov,² where an asymptotic analytical solution for the wave amplitude and frequency modulations caused by an oscillating wave emitting boundary in a quiescent background medium is given, under the assumption that the timescale of the boundary motion is large as compared to the emitted wave period. While Fig. 6 demonstrates linear wave behavior in the limit of a very small surface gravity, a more detailed analysis is needed to understand the asymptotic behavior for very large magnitudes of the surface gravity.

Furthermore, the change of variables may be extended to the time variable, which would admit Lorentz transformations as described by Wang.²³ While the transformation in physical space is sufficient as long as the wave emitter moves along with the background fluid flow as it is the case in the present work, Lorentz transformations may indeed help to preserve the well-posedness of the governing wave equation in the case of relative motion between the wave emitter and the background fluid. Finally, the coordinate transformation may also be harnessed to redistribute grid points in physical space such that sharp features of the waveforms are accurately resolved. In such a moving grid method, a fixed, uniform grid is implemented in transformed space and the coordinate transformation governing the grid point motion is linked to a measure of targeted accuracy.^46,47

Within the scope of an acoustic black hole analogy, the acoustic wave-flow configuration can be parameterized in terms of the Helmholtz number $He=c2/(ghλa)$⁠, where c, $gh$⁠, and $λa$ are the speed of sound, the acoustic black hole's surface gravity,⁹ and the emitted wavelength, respectively. The Helmholtz number in this particular form is shown to describe geometrically similar wave characteristics for different combinations of $gh$ and $λa$⁠. Its applicability is not restricted to spherically symmetric flows, but extends to any quasi one-dimensional geometry. However, the wave amplitude modulations along geometrically similar wave characteristics depend on the magnitude of $gh$⁠. While the acoustic black hole analogy is a useful means to characterize the acoustic wavelength relative to the acceleration of a stationary but non-uniform background flow field, the observations regarding the nonlinear wave amplitude modulations are expected to apply to subsonic flows in a very similar fashion as they merely depend on the flow acceleration. For small $gh$⁠, the linear wave behavior is recovered in the proximity of the sonic horizon. Therefore, the present work promotes the acoustic black hole analogy for the study of nonlinear acoustic wave dynamics and Doppler modulations in accelerating flow fields. Regarding the computational method, it is shown that a suitable Galilean-type coordinate transformation from a moving/deforming physical domain to a fixed computational domain can maintain the well-posedness of the governing wave equation in the presence of a transonic background flow. Following the suggestions by Gregory et al.²² and Wang,²³ the flexibility afforded by the coordinate transformation can be further exploited in order to isolate terms that are exclusively related to nonlinear Doppler effects from those terms that constitute the invariants of the wave equation. This is expected to greatly simplify the analysis of nonlinear Doppler modulations in moving flow fields.²² 

This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Grant Nos. 441063377 and 443546539.

The authors have no conflicts to disclose.

Soeren Schenke: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Fabian Sewerin: Formal analysis (supporting); Funding acquisition (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Berend van Wachem: Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (supporting); Writing – review & editing (equal). Fabian Denner: Funding acquisition (equal); Methodology (supporting); Project administration (equal); Resources (equal); Software (supporting); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.