Compressible potential flows around round bodies: Janzen-Rayleigh expansion inferences

The subsonic, compressible, potential flow around a hypersphere can be derived using the Janzen-Rayleigh expansion (JRE) of the flow potential in even powers of the incident Mach number $\mathcal{M}_\infty$. JREs were carried out with terms polynomial in the inverse radius $r^{-1}$ to high orders in two dimensions (2D), but were limited to order $\mathcal{M}_\infty^4$ in three dimensions (3D). We derive general JRE formulae to arbitrary order, adiabatic index, and dimension. We find that powers of $\ln(r)$ can creep into the expansion, and are essential in 3D beyond order $\mathcal{M}_\infty^4$. Such terms are apparently absent in the 2D disk, as we confirm up to order $\mathcal{M}_\infty^{100}$, although they do show in other dimensions (e.g. at order $\mathcal{M}_\infty^2$ in 4D) and in non-circular 2D bodies. This suggests that the disk, which was extensively used to study basic flow properties, has additional symmetry. Our results are used to improve the hodograph-based approximation for the flow in front of a sphere. The symmetry-axis velocity profiles of axisymmetric flows around different prolate spheroids are approximately related to each other by a simple, Mach-independent scaling.


Introduction
While ideal flows (inviscid, with no heat conduction or additional energy dissipation effects; Landau & Lifshitz 1959) are an extreme limit, they play an important role in research, for example (i) as a basis for more realistic flows, with additional effects such as viscosity; (ii) for modeling the bulk of weakly-interacting Bose-Einstein condensate (BEC) superfluids, which can be approximated as an inviscid, compressible fluid with a polytropic index γ = 2; (iii) for modeling flow regimes which are not sensitive to the level of weak viscosity, such as in front of a round object; and (iv) for code validation and pedagogical reasons. Janzen-Rayleigh expansions (JREs) can be broadly identified as an expansion of the flow variables in terms of the Mach number. JREs were used by Janzen (1913) and Rayleigh (1916) as a method to study D'Alembert paradox with the addition of compressibility effects. They considered an inviscid, compressible, subsonic flow of a fluid with a polytropic equation of state (EoS), with no external forces or initial vorticity. Specifically, a steady flow was assumed around a disk in two dimensions (2D) or around a sphere in three dimensions (3D), with an incident uniform flow far from the body. Introducing the flow potential, a scalar non-linear partial differential equation (PDE) was obtained, and expanded in the incident Mach number squared, M 2 ∞ . Using the same setup, the JRE was used to solve the flow around various blunt objects (Goldstein & Lighthillm 1944;Hasimoto 1951;Hida 1953;Longhorn 1954;Imai 1957;Kaplan 1957;Van Dyke 1958). The JRE has been generalized to several areas of research, such as vortex flows (e.g. Moore & Pullin 1991, 1998Meiron, Moore & Pullin 2000;Leppington 2006;Crowdy & Krishnamurthy 2018), porous channel flows (Majdalani 2007;Maicke & Majdalani 2008;Maicke, Saad & Majdalani 2010;Cecil, Majdalani & Batterson 2015), and acoustics (Slimon, Soteriou & Davis 2000;Moon 2013), and could be used in a wide range of other applications, as we demonstrate below.
The 2D flow around a disk has been researched extensively (Janzen 1913;Rayleigh 1916;Van Dyke & Guttmann 1983;Guttmann & Thompson 1993), for example in search for a solution to the transonic controversy, namely, "the existence, or non-existence, of a continuous transonic flow, that is, without a shock wave, around a symmetrical wing profile, with zero incidence with respect to the undisturbed velocity" (Ferrari 1966). These and other problems that require a high-order expansion could not be explored as thoroughly in the 3D case, because previous JREs for the sphere were limited to second, i.e. M 4 ∞ , order (Kaplan 1940;Tamada 1940). Indeed, a power series in r −1 yields a nonphysical behaviour at the third, i.e. M 6 ∞ , expansion order. A higher order expansion in 3D is also needed, for example, to derive the flow in front of a sphere, in order to model axisymmetric bodies in various fields of physics (e.g. Keshet & Naor 2016, and references therein). Flows in higher dimensions, d > 3, are also important, mainly for theoretical and pedagogical purposes.
We study the steady, inviscid, compressible flow around a hypersphere in d spatial dimensions, and provide explicit results for the d = 2 (disk) and d = 3 (sphere) cases. We derive the JRE in 3D beyond the presently available second-order, to an arbitrary order. As an example of the usage of such an expansion, we compare the axisymmetric flow in front of a sphere to that in front of a spheroid, and show that a simple scaling approximates the flow well for the prolate case. Furthermore, three orders in the JRE are sufficient to adequately describe these flows, for Mach numbers ranging from the incompressible to the sonic regimes.
The paper is organized as follows. In §2, we derive the JRE equations from the hydrodynamical ones. In §3, we show that each term in the JRE of the flow potential around a hypersphere is a finite sum of a product of powers of the radial coordinate, powers of its logarithm, and a set of orthogonal functions (Jacobi polynomials) of the polar coordinate. In §4, we outline the semi-analytic algorithm to compute the JRE and a numerical pseudospectral method we use to solve the non-linear compressible flow. We present results from our numerical solver and from the semi-analytical JRE, and highlight the compressibility to compute the JRE, in §5. The JRE is used to improve a hodographic approximation for the flow in front of the sphere in §6. In §7, we show that the axisymmetric flow in front of prolate spheroids is well approximated by that of the scaled flow in front of a sphere. We summarize and discuss our results in §8. In appendix §A, we provide explicit values of the coefficients of the series representation of the JRE for low orders of the flows around a disk and a sphere. Appendix §B discusses the hodographic approximation and general results. In appendix §C we detail the numerical solver.

Janzen-Rayleigh Expansion
Consider an isentropic flow in d-dimensions with no external forces of a perfect fluid with a polytopic, ideal gas EoS of adiabatic index γ, p ∝ ρ γ . (2.1) The equations governing the flow are the continuity equation, and the Euler equation, Here, ρ, u, and p are respectively the mass density, velocity, and pressure, c ≡ (dp/dρ) 1/2 = (γp/ρ) 1/2 is the sound velocity, and ∇ is the del operator in d-dimensions.
(2.4) For our inviscid, steady flow, (2.3) yields the Bernoulli principle, namely, the quantity wu 2 + c 2 is constant along streamlines, where we define w ≡ (γ − 1)/2. We consider the flow along a body for a uniform flow incident from infinity, where r is the radial coordinate, z is the coordinates along the flow, and the subscript ∞ indicates the far region, tending to spatial infinity. Using the boundary conditions (2.5) and the assumption that every streamline starts at infinity, the Bernoulli equation may be written as Considering the subsonic regime, we henceforth restrict the discussion to a potential flow, writing the velocity as the gradient of the flow potential, u ≡ ∇φ . (2.7) Isolating c 2 from (2.6), substituting it into (2.4), and using the potential (2.7), we obtain a single PDE for the potential φ (Rayleigh 1916), We normalize the variables to obtain them in a dimensionless form, by taking φ → (u ∞ R)φ, and r → Rr (which also normalize the velocity u → u ∞ u), where R is a characteristic length scale, for example the radius of a hypersphere. Defining M ∞ ≡ u ∞ /c ∞ as the Mach number at infinity, we then arrive at the dimensionless We supplement this equation for the potential with two boundary conditions (BCs): the uniform incident flow from infinity, and the slip, no-penetration condition on the surface of the body (n · u = 0), namely φ (r → ∞) = z = r cos θ andn · ∇φ = 0 . (2.10) Here,n is the normal to the body, and we use hyperspherical coordinates: a radial coordinate 0 r ∞, a polar angle 0 θ π measured with respect to the uniform flow at infinity, and (d − 2) additional angles. For axisymmetric flows such as in the case of a hypersphere, the flow is by symmetry independent of these additional angles.
To arrive at the JRE, we substitute the expansion into (2.9), and isolate the different powers of M ∞ . This leads to a recursive PDE for φ m (for m 1), with an initial equation along with BCs φ m (r → ∞) = δ m,0 r cos θ (2.14) at spatial infinity, andn · ∇φ m = 0 (2.15) on the body. Here, δ i,j is the Kronecker delta symbol.
The zeroth-order equation (2.13) is the Laplace equation for φ 0 , corresponding to the incompressible limit M ∞ → 0. For higher orders, equation (2.12) can be regarded as a set of Poisson equations for φ m at each order m, with a source term being a function of the lower-order solutions, {φ i } m−1 i=0 , and their derivatives. In general, the solution to (2.12) at any order m 1 under the BCs (2.14) and (2.15) is a sum of an inhomogeneous solution and a homogeneous solution, solves the Laplace equation. The JRE (2.11) is constructed by iteratively determining the functions φ m . One starts by deriving φ 0 , obtained as the solution to the zeroth-order Laplace equation (2.13) under the BCs. Increasingly higher-order functions φ m are then incrementally derived, by solving the corresponding Poisson equations (2.12) under the same BCs. At each order m 1, the inhomogeneous solution φ (in) m is fixed by the source term in the corresponding (2.12), which is explicitly written once all lower order functions φ m are determined. This solution is then combined with an expansion φ (ho) m of homogenous solutions, the coefficients of which are determined by applying the BCs to the resulting (2.16).

JRE for a hypersphere
Simple solutions exist for the flow around a d-dimensional hypersphere, for which the Laplace equation has known analytic solutions, and the Poisson equation is easily solved. Choosing the characteristic length R as the hypersphere radius, the normalized slip BC (2.15) here simplifies to ∂ r φ (r = 1, θ) = 0 .
(3.1) Considering the hyperspherical and axial symmetries, it is useful to write the general solution to the Laplace equation (2.13) as the infinite sum of positive and negative radial powers (Feng et al. 2011 where we define µ ≡ cos θ, the normalized Jacobi polynomials and numerical coefficients A n and B n . Here, J (α,β) n (µ) and C (α) n (µ) are respectively the standard Jacobi and Gegenbauer polynomials.
The source term in (2.12) is a sum of multiple terms, each composed of even derivatives in θ and odd multiples of functions φ. Therefore, the polar dependence will always exhibit the same symmetry as the zero-order solution φ 0 . As the hypersphere is isotropic, the incompressible flow also shows a backward-forward symmetry. We conclude that the functions J (d) n appearing in φ m on all orders m have only odd n. Next, we construct the hypersphere JRE by considering increasingly larger orders m.

Zeroth order: m = 0
For the m = 0 order, the infinity BC (2.14) allows only the radially linear and decaying terms in the solution (3.2). The slip BC (3.1) restricts the solution further, allowing for only the decaying term proportional to J Henceforth, we omit the dimension superscript (d) unless necessary.
3.2. First order: m = 1 The first-order potential, φ 1 , is determined by which, by the zeroth-order (2.13), simplifies to Plugging the m = 0 solution (3.4) into (3.6) gives the Laplace equation with the explicit source term, The orthogonal decomposition of Jacobi polynomials (Chaggara & Koepf 2010) then yields

I. S. Wallerstein and U. Keshet
This result may be compactly written in the form m,k,n are the expansion coefficients of the order-m Poisson source term, for radial order k and angular order n.
As the Poisson equation is linear, suffice to solve, for arbitrary k and n, the equation This equation has the particular, inhomogeneous solution , (3.11) provided that k / ∈ {−d − n, n − 2}. These two exceptional values of k do not occur at order m = 1, as seen from (3.8), but they may appear at higher orders, as discussed below in §3.3. We may now expand the inhomogeneous solution as where the numerical coefficients s 1,k,n are determined by equating (3.8) and (3.9). From (3.8) and (3.11) we see that the largest (i.e. least negative) power of r for m = 1 is includes only radially declining terms. For higher orders, (2.12) combines the derivatives of multiple φ m functions, but the highest radial power k max in the source term remains no larger than −d − 1. Consequently, the inhomogeneous solution for all m 1 orders includes only radially declining terms. This conclusion, combined with the BCs, indicate that the homogeneous solution may be expanded with only negative powers of r for any m > 0, m,n are numerical coefficients, to be determined below. The full solution for m = 1 now becomes (3.14) where we introduced the numerical coefficients These coefficients may now be determined from the slip BC (3.1), implying that . (3.17) As the coefficients s (in) 1,k,n are known from (3.8) and (3.9), the solution (3.14) is completely specified by (3.15) and (3.17).

Higher orders: m 2
Substituting φ 0 and φ 1 in the Poisson equation (2.12) for the next, m = 2 order, results again in an equation of the form (3.9), but now for φ 2 and with different coefficients. Equations of the same form persist for increasingly higher orders m, as long as the source term in (2.12) is free of non-zero terms s m,k,n which satisfy one of the aforementioned special conditions, k = −d−n or k = n−2. Indeed, the source term consists of differential operators of the form ∇ 2 f (r, µ) and ∇f (r, µ) · ∇g(r, µ), where f and g are constructed from lower order φ terms. As long as f and g are sums of terms, each of which is a product of powers of r and Jacobi polynomials in µ, the source term would remain of the form (3.9).
For the special cases k ∈ {−d − n, n − 2}, the solution to (3.10) is no longer given by (3.11), which formally diverges. It is sufficient to consider the former case, k = −d − n, as the latter, k = n − 2, is then obtained by the transformation n → −n − d + 2, under which J (d) n remains invariant. Here, the inhomogeneous solution to (3.10) becomes Once a ln(r) term of (3.18) appears in the inhomogeneous solution for φ m , as a result of a special k term in the source-term expansion (3.9), a power-series in r as in (3.14) is no longer sufficient. Indeed, when x = ℓ + 1 and y = p ln (r) for integers ℓ and p, contains logarithmic, and not only polynomial, radial terms. Since ln(r) cannot be expanded as a power series of r −1 valid over the full 1 r < ∞ range, one must then take into account logarithmic terms in the expansion of φ m and in the resulting source terms of higher order functions.
Generalizing the prototypical source term in (3.9) to include a logarithm of r to some power l, we must therefore also consider the generalized equation (3.21) As long as k, d, and n do not satisfy one of the special cases the inhomogeneous solution to this equation is (3.24) The case k = n − 2 is again obtained using the transformation n → −n − d + 2.
( 3.25) The overall solution for φ m may now be expanded as the finite sum where the numerical coefficients s m,k,n,ℓ are determined using the BCs in analogy with the above m = 2 discussion. The expansion (3.26) is complete, as the inhomogeneous solution yields only powers of r and ln (r), and the homogeneous solution yields only powers of r, so the resulting higher-order source terms are again of the form (3.21). One can prove, by induction, the following rules for the m 1 indices in (3.26), and 1 n 1 + 2m ; (3.28) the summation limits for the logarithmic term are discussed in §5.

Semi-analytical and numerical solvers for disk and sphere flows
We explicitly solve equation (2.12) for the hypersphere in the d = 2 and d = 3 cases, namely, we derive the flows around a 2D disk and around a 3D sphere. The normalized Jacobi functions reduce to the Chebyshev polynomials T n in 2D, and to the Legendre polynomials P n in 3D, As discussed in §3, we consider the generalized expansion (2.11) and (3.26) of the flow potential in each case, where a and b are the expansion coefficients in 2D and 3D, respectively and the summation index n takes only odd values. We compute these JREs analytically, following the steps outlined in §3. The m = 0 term is obtained from (3.4). For each order m > 0, we compute the source terms on the RHS of (2.12), and decompose them into a sum of terms proportional to r k J n (µ) ln ℓ (r) as in (3.21). The inhomogeneous solution φ (in) m is then obtained as the corresponding sum of solutions of the form (3.23)-(3.25), and added to the homogeneous solution φ (ho) m of (3.13). The coefficients s (ho) m,n in the latter sum are determined from the slip BC, as demonstrated for m = 1 in (3.17). This process is repeated until the necessary accuracy of the JRE is reached.
The JRE results are compared to the flow obtained from a numerical solution to the non-linear equation (2.9). We follow Pham et al. (2005) and use the pseudospectral collocation method (Boyd 2001) for a fast convergence. In the pseudospectral method, the solution to a PDE is expanded in terms of base functions, each being a product of individual (and typically orthogonal) functions for each coordinate. In collocation methods, the PDE is then evaluated and solved at a set of points, usually taken as the roots of the basis functions. This leads to a set of linear equations for the basis function coefficients, which are solved to give an approximate PDE solution. In general for pseudospectral methods, if the solution is smooth (all derivatives of all orders are continuous) then the convergence is exponential in the number N of collocation points (Boyd 2001), O e −N . Pham et al. (2005) use the Chebyshev-Fourier set of basis functions, invoking cos(nθ) and sin(nθ) with integer n as angular basis functions. As discussed in §3, the analytical solution for φ is a sum of cos(nθ) with odd positive n, so we expand the flow potential in odd cosine functions. To achieve the same accuracy as Pham et al. (2005), we need only a fourth of the angular basis functions.
For the radial part, we apply an inversion map ̺ ≡ r −1 and expand in Chebyshev Polynomials T k (̺), to resolve both small and large r behavior. The pseudospectral expansion thus becomes (4.5) Here, k max and n max are the radial and angular resolution, i.e. the number of collocation points, and c k,n are real coefficients. The BCs are guarantees if they are satisfied by the potential φ 0 , chosen as the incompressible solution (3.4). Since the effective computational domain and symmetry are the same in our 2D and 3D frameworks, we expand both flows, around a disk and around a sphere, in the same basis functions (4.5).
Appendix §C provides more details on the the calculations of c k,n and demonstrates the numerical convergence.

Results
We calculate the JRE (4.3) and (4.4) up to order m = 30 for the disk and m = 18 for the sphere. Appendix §A provides explicit expressions for the coefficients a and b up to order m = 3 (i.e. M 6 ∞ ) for a general γ. For the pseudospectral code we use a resolution up to 32 Chebyshev and 64 odd cosine functions, but as shown below, much lower orders and resolutions are sufficient to capture the behaviour of the flow (see figure 2 and 6). The results are presented in the following figures mainly for γ = 7/5. Different choices of γ give qualitatively similar results, as demonstrated for γ = 5/3 in figure 3b. Figure 1 shows the flow around a disk (left panel) and a sphere (right panel) for incident Mach numbers at infinity approaching the respective critical/sonic Mach numbers. The latter are tuned to yield a sonic flow at the equator of the hypersphere, M(r = 1, θ = π/2) = 1, leading to M disk flow is symmetric under both x and z reflections, i.e. θ → −θ and θ → π − θ, so it is sufficient to plot only one quadrant of the x − z plane. We thus utilize all four quadrants to show the differences between the JRE obtained up to different, m = 0 to m = 5, orders (see labels). To show the JRE corrections to the flow, we plot the streamlines (arrows) that pass through a set of equidistant points at x = ±2, so the differences between quadrant arrows are meaningful. Comparing the different quadrants shows that the compressible effects captured by higher orders m raise the velocity and lower the density in the equatorial plane; most but not all of the change already transpires as m = 0 is raised to m = 1. Figure 2 shows the compressible contribution to the flow around a sphere; qualitatively similar results are obtained for the disk. The contribution is computed based on JREs of different orders, as well as on the numerical solution, and shown for the polar velocity on the sphere (column a), the polar velocity at the equatorial plane (column b), and the radial velocity along the symmetry axis (column c). We plot these profiles for two different flows, with incident Mach numbers M ∞ = 0.1 to demonstrate a subsonic case (row 1), and M ∞ = M c ≃ 0.5619 for the critical case (row 2). As seen from row 1, at low Mach numbers, the JRE converges rapidly; the m = 1 JRE is sufficient for accurately (within ∼ 0.005% in u for M ∞ = 0.1) capturing the compressible effects. Convergence is slower for larger M ∞ , but manageable (compressible effects captured within ∼ 5% in u for m = 1) even in the critical limit.
For the flow around a disk in 2D, we find no logarithmic terms in the flow potential, for any γ. Namely, all computed JRE (4.3) coefficients a with ℓ = 0 vanish, as illustrated in table 2 up to m = 3 order. We confirm this behavior for arbitrary γ up to order m = 30. Using a prescribed γ speeds up the JRE computations, allowing us to reach higher orders. We thus compute the JRE up to order m = 50 for the specific cases γ ∈ {1, 7/5, 5/3, 2}, corresponding respectively to isothermal, ideal diatomic, ideal monatomic, and weaklyinteracting bose, gasses. In all of these cases, we find no logarithmic terms in the flow potential.
In contrast, logarithmic terms are unavoidable in the flow around a sphere in 3D, for any γ. Indeed, the JRE (4.4) shows the first logarithmic term at order m = 3, for k = −8, n = 7 and ℓ = 1, as indicated in table 4. This coefficient is proportional to 5 + 7γ + 2γ 2 , which vanishes only for negative, non-physical γ values. In addition to such ℓ = 1 terms for order m 3, we find ℓ = 2 terms for order m 6, ℓ = 3 terms for order m 9, and so on, with the highest logarithmic term order increasing by one every three orders in m. Furthermore, the term with the highest logarithm power is also proportional to ∝ r −2m−2 P 2m+1 (µ). This behavior is verified up to order m = 18 for arbitrary γ > −1, and up to order m = 30 for the specific values of γ chosen above.
While proving that these effects persist as m increases to infinity is beyond the scope of the present work, we hypothesize that they do: Conjecture 1. Logarithmic JRE terms never appear in the flow of a polytropic γ 1 fluid around a disk.

Example: axial hodographic approximation
The flow in front of a blunt body has important implications in diverse applications, and is well approximated in both subsonic and supersonic regimes by expanding the perpendicular gradients q ≡ ∂ θ u θ | θ=0 in terms of the parallel (normalized) velocity u = |u r | (Keshet & Naor 2016) In the subsonic regime of this hodographic approximation, q = q(u) is defined in the 0 < u < 1 region between stagnation and spatial infinity, and is expanded around stagnation in the form q (u) ≈ q 0 + q 1 u + q 2 u 2 + q 3 u 3 , which we designate as a 3rd order hodographic approximation, Hodo(3). It can be shown that q 0 does not vary much with respect to the incompressible case, q 1 = −1/2, and q 2 is small with respect to q 3 (Keshet & Naor 2016). Using the JRE for the sphere, here we calculate the coefficients q i analytically, for given order m. The resulting expressions for the coefficients are provided for 1 m 4, as a function M ∞ and γ, in appendix §B. For illustration, table 1 provides the numerical values of these coefficients for the critical Mach number with γ = 7/5. The coefficients converge rapidly, reaching three-digit accuracy for m = 4. The effect of compressibility can be seen to be small, as q 2 and q 3 are smaller than q 0 and q 1 by two orders of magnitude. We confirm the small deviation of q 0 from its incompressible value 3/2, the precise result q 1 = −1/2 for all orders, and that q 2 ≈ q 3 /3 is small albeit non-negligible. Figure 3 shows the compressible contribution to the radial velocity along the symmetry axis of a sphere, with and without the JRE corrections, for a flow at the critical Mach number. Results are shown both for γ = 7/5 (left panel) and γ = 5/3 (right panel), based on the pseudospectral code at resolution of (8, 8) (dotted-dashed curve), the full JRE of order m = 3 (solid), on the hodographic approximation of Keshet & Naor (2016, dashed), and on our improved hodographic approximation (dotted). The corrected hodographic approximation provides a much better fit to the JRE and to the actual flow. The two panels slightly differ because they pertain to different M ∞ values; for the same M ∞ , the axial flow for different 1 < γ < 2 values are nearly indistinguishable.

Example: flow in front of spheroids
It has been suggested (Keshet & Naor 2016) that the flow in front of an axisymmetric body is well approximated by the flow in front of a sphere, rescaled to give the same nose curvature. Consider general prolate or oblate spheroids of the form x 2 + y 2 /α+z 2 /α 2 = 1, chosen to be axisymmetric along theẑ, flow axis, and with unit curvature at the nose. We shift the z coordinate such that the resulting spheroid overlaps with the unit sphere at the nose z → z + 1 − α 2 . The pseudospectral code (details in §C) is modified to solve the flow around such spheroids. Figure 4 demonstrates the shapes of a sphere (α = 1) and of four shifted prolate spheroids (α ∈ {2, 3, 4, 5}), along with the incompressible flow along the latter, most prolate case. Figure 5 shows both the incompressible (left panels) and the critical (right panels) flows in front of these bodies. The normalized velocity u (upper panels) varies with α, M ∞ (and slightly with γ, see figure 5d and 3). However, a nearly universal result is obtained for the scaled velocity approximately insensitive to the Mach number, prolate body profile, and value of γ. For oblate spheroids (α < 1), the velocity does not scale to the universal curve.

Summary and discussion
We generalize the Janzen-Rayleigh expansion (JRE) of a hypersphere to arbitrary dimension, providing an exhaustive solution that supplements previous approaches with additional, usually necessary, logarithmic radial terms. Such a generalization is found to be essential in 3D, required for obtaining the correct solution for the sphere even at third (m = 3, i.e. M 6 ∞ ) order, although it is apparently not needed for the special case of a disk in 2D. The resulting, arbitrary-order JRE is a useful tool for studying various problems, as we demonstrate by generalizing and extending previous solutions for the axial flow of a sphere, and by presenting a simple approximate scaling that generalizes the axial flow for prolate spheroids.
The JRE of a potential, compressible flow is derived in general (in §2) and around a hypersphere of arbitrary dimension d (in §3). The expansion is based on combining the continuity equation (2.4) and the Bernoulli equation (2.6) into a single equation (2.9) for the flow potential, φ, which depends on the incident Mach number M ∞ far away from the hypersphere. An expansion of φ in powers of the M ∞ results in a set of recursive equations (2.12), with the zero order being the incompressible solution (3.4). The equations for higher orders, φ m , are derived as a sum of particular Poisson equations (3.21), with an exhaustive solution that is a finite sum (3.26) of terms combining powers of r, Jacobi polynomials J (d) n (cos θ), and, in addition, powers of ln (r) that emerge when the Poisson source terms satisfy one of the special conditions (3.22).
Previous JRE calculations were typically carried out with only part of the general solution, each order being considered as a sum of powers of r and Jacobi polynomials alone. This approach works well in 2D, where one does not encounter any divergence when thus avoiding logarithmic terms, at least up to order m = 50. However, this approach severely limits the JRE in 3D, as without logarithmic terms, the m = 3 order diverges. The inclusion of logarithmic terms allows us to calculate terms much higher than available for any previous work in 3D (Kaplan 1940;Tamada 1940;Lighthill 1960;Fuhs & Fuhs 1976), and with higher accuracy (Frolov 2003). Furthermore, we are able to compute the JRE to arbitrary order in any dimension, providing interesting insights on the flow. For instance, the JRE of the 4D hypersphere shows logarithmic terms at order m = 2, indicating that the 3D sphere is not unique in requiring such terms.
After deriving the general solutions to all possible Poisson equations that can emerge in the problem (3.23, 3.24, and 3.25), we calculate the JRE analytically for general γ to order m = 30 in 2D and m = 18 in 3D. For select values of γ, we proceed to compute 50 JRE orders in 2D and 30 orders in 3D.
The m = 3 order JRE includes a γ-dependent logarithmic term proportional to 5 + 7γ +2γ 2 (table 4), which does not vanish for γ = 1 (i.e. w = 0). This shows that both nonlinear terms on the RHS of (2.12) contribute to the logarithm. We find that the highest order of any logarithmic term in the expansion increases by one every three orders in m, and when this occurs, only one logarithmic term appears, multiplied by a single power of r and a single Legendre polynomial ∝ r −2m−2 P 2m+1 (µ) (see table 4 for m = 3). We speculate that these findings are true for higher orders in the JRE and all physical values of γ (conjecture 2).
We expect the absence of logarithmic terms in the 2D disk JRE to persist for all orders m and for all γ (conjecture 1). The absence of logarithmic terms is the disk JRE is not, however, a property of 2D flows in general. For instance, consider the incompressible flow of potential φ 0 = (r+1/r) cos θ+B 3 r −3 cos(3θ) around an algebraic body, which is defined by the no-slip BCs and this potential, unique in leading to a finite polynomial equation (Wallerstein & Keshet, in preparation). Here, B 3 is a coefficient in the general solution of the Laplace equation (3.2). The non-circular shape of the body leads to an infinite number of homogeneous terms in the first-order JRE potential φ 1 , and thus to an infinite number of both inhomogeneous and homogeneous terms in φ 2 . It is nevertheless possible to isolate the finite number of terms that contribute to the nonzero term ∝ r −5 cos(5θ) ln(r) in φ 2 .
Regardless of logarithmic terms, the angular part the solution includes Jacobi polynomials of odd n for all JRE orders, so the flow is symmetric fore and aft of the hypersphere. Hence, as long as the JRE converges and no additional effects are included, there is no net force on the hypersphere, and the flow is drag-free for all dimensions. The d'Alembert paradox thus persists in all dimensions.
Even though we calculate the JRE to high orders, many physical features can be adequately captured with only a few low orders (figure 2), even for a sonic flows. For example, we find an accuracy better than ∼ 1% in u for m = 5. The low JRE orders (m 5) show that compressibility effects are more prominent in the vicinity of the hypersphere, and in particular at the equatorial plane (figures 1 and 2).
With the ability to calculate the JRE to any desired accuracy, we compare it to other approximation methods, such as the hodographic approximation for the flow on the symmetry axis in front of a sphere, which has important physical implications even in the inviscid regime (e.g. Keshet & Naor 2016). The hodographic approximation of Keshet & Naor (2016) performs better than the JRE of order m = 0, but worse than the JRE of order m = 1 for any γ (figure 3a and 3b). We use the JRE to improve the hodographic approximation (appendix §B), such that it performs better than JRE of order m = 2, but worse than the JRE of order m = 3 again for any γ, but can be continued to the supersonic regime (Keshet & Naor 2016).
It has been speculated that the flow in front of the sphere can be applied to axisymmetric bodies with a similarly scaled nose curvature (Keshet & Naor 2016, and references therein), such as prolate spheroids (figure 4). While the flows in front of different spheroids are not identical (figure 5a and 5b), they can be approximately mapped onto each other with a universal scaling, independent of γ and M ∞ , for flows ranging from the incompressible (figure 5c) to the sonic (5d) regimes. The velocity u (α) r in front of prolate spheroids with semi-axes α > 0 and α 1/2 can be well approximated by a scaled JRE for a sphere It is natural to ask if such scalings can be identified using the generalized JRE in front of other bodies and in other dimensions, but this is beyond the scope of the present work.
We thank Y. Lyubarsky, D. Kogan and E. Grosfeld for helpful discussions. This research has received funding from the GIF (Grant No. I-1362-303.7 / 2016), and was supported by the Ministry of Science, Technology & Space, Israel, by the IAEC-UPBC joint research foundation (Grants No. 257/14 and 300/18), and by the Israel Science Foundation (Grant No. 1769/15).

Appendix A. Explicit JRE coefficients for low orders
Expressions for the coefficients of the JRE up to order M 6 ∞ , valid for any γ, are provided in table 2 for a disk, and in tables 3 and 4 for a sphere. Each table is broken to parts, each corresponding to a different m index of M 2m ∞ . Columns correspond to the n appearing in cos (nθ) or P n (cos θ) basis functions, and the rows correspond to k of the r k terms. All coefficients in these tables pertain to non-logarithmic terms, i.e. to ∝ ln ℓ (r) terms with ℓ = 0, except for the very last term in table 4, which has ℓ = 1.

Appendix B. Hodographic approximation of the solution of a radial flow
In §6 we discuss the hodographic approximation for the flow in front of a sphere, and use the JRE to obtain a more accurate expansion. To relate between the radial velocity and the radius we use Eq. (7) from Keshet & Naor (2016), with W 2 = 2/ (γ + 1) , S 2 = 2/ (γ − 1) (notice that S 2 = w −1 ) and which is the Mach number with respect to the stagnation point. To complete the integral we calculate the coefficients q i for i ∈ {0, 1, 2, 3} using the JRE. We provide explicit expressions for general γ in table 5 which are valid for 0 M ∞ M c (with q 1 ≡ 1/2 identically (Keshet & Naor 2016)). These values give a good approximation for the radial velocity near the stagnation point, but deviates far from the sphere. To ensure the correct BC far from the body we add a fifth term to q (u) such that the denominator in the integral in (B 1) vanishes at r → ∞ i.e. q 4 = 1 − q 0 − q 1 − q 2 − q 3 .

Appendix C. Details and convergence of the pseudospectral solver
In order to capture the behaviour of the flow far and near the sphere (we review the implementation for the 3D case. The 2D case is similar with a difference only in the differential operators) we map the radial domain from r ∈ [1, ∞] to ̺ = r −1 ∈ [0, 1]. The differential operators with respect to ̺ are (C 1) The pseudospectral method requires the expansion of the solution as a series of functions (Boyd 2001) (not required to be orthogonal) and evaluating the equation at collocation points. To ensure the polar BCs (no polar velocity on the poles i.e. Neumann BC for 3D and periodic boundary in 2D) we expand in cosine functions only, furthermore, the back and front reflection symmetry of the sphere limits the cosine functions to only odd cosines (as analytically shown for the JRE in §3), cos [(2n + 1) θ]. In the radial direction we use the Chebyshev polynomials of the first kind (4.5). We take the collocation points to be a m,k,n m = 0 m = 1 m = 2 m = 3 n = 1 n = 1 n = 3 n = 1 n = 3 n = 5 n = 1 n = 3 n = 5 n = 7 k = +1 1 k = −1 1     Table 5. Coefficients of the first four terms in the Taylor expansion of q (u), for different JRE orders m, at general Mach number and γ. ̺ i = cos πi 2k max , 1 i k max − 1; θ j = πj 2 (n max + 2) , 1 j n max + 1. (C 2) We do not consider the collocation points θ ∈ {0, 1} because the expansion in odd cosines fulfills the BCs identically. On the points ̺ ∈ {0, 1} we solve the BCs at infinity and at the sphere respectively and not the non-linear equation (2.9). To solve the nonlinear equation we use a simple Newton-Raphson algorithm where at each step we solve a linear set of equations. In all computation we take the error to be the maximal deviation of the equation at the collocation points and take a tolerance of 10 −11 . Most of the computations converge rapidly and do not require more than seven Newton-Raphson steps. Figure 6 shows the compressible contribution to the polar velocity on the sphere for various resolutions of the pseudospectral code. The code is quite converged, for most angles, except close to the equator and poles. A resolution of four by four shows the largest relative error (with respect to resolution of 16 by 32), which is smaller than 2% at the equator. To test the robustness of the computation, we expand the solution with different functions, e.g. odd and even cosines as well as in sine functions, this resulted in the vanishing coefficient of the even cosine and sine functions. For the radial functions we expand in a linear combination of Chebyshev polynomials which fulfills the BC identically and in Legendre polynomials. All combinations gave the same results with respect to the tolerance.
(C 3) (we omit the azimutal coordinate ϕ from symmetry considerations). In this set of coordinates the body is described by the equation cosh(µ b ) = α/(α − 1). We normalize this coordinate and inverse it such that the computational domain is [0, 1] × [0, π/2] just as in the case of the sphere. This changes the differential operators to (after variable change µ = µ b ̺ −1 ): We take φ 0 to be a function that fulfills both BCs