A Sansica
email: sansica.andrea@jaxa.jp

Y Ohmichi

J.-Ch Robinet

A Hashimoto

Affiliations

Laminar supersonic sphere wake unstable bifurcations

The laminar sphere unstable bifurcations are sought at a Mach number of M∞ = 1.2. Global stability performed on steady axisymmetric base flows determines the regular bifurcation critical Reynolds number at Re reg cr = 650, identifying a steady planar-symmetric mode to cause the loss of the wake axisymmetry. When global stability is performed on steady planar-symmetric base flows, a Hopf bifurcation is found at Re Hopf cr = 875 and an oscillatory planar-symmetric mode is temporally amplified. Despite some differences due to highly compressible effects, the supersonic unstable bifurcations present remarkably similar characteristics to their incompressible counterparts, indicating a robust laminar wake behavior over a large range of flow speeds. A new bifurcation for steady planar-symmetric base flow solutions is found above Re > 1000, caused by an anti-symmetric mode consisting of a 90 ○ rotation of the dominant mode. To investigate this reflectional symmetry breaking bifurcation in the nonlinear framework, unsteady nonlinear calculations are carried out up to Re = 1300 and dynamic mode decomposition (DMD) based on the combination of input data low-dimensionalization and compressive sensing is used. While the DMD analysis confirms dominance and correspondence in terms of modal spatial distribution with respect to the global stability mode responsible for the Hopf bifurcation, no reflectional symmetry breaking DMD modes were found, asserting that the reflectional symmetry breaking instability is not observable in the nonlinear dynamics. The increased complexity of the wake dynamics at Re = 1300 can be instead explained by nonlinear interactions that suggest the low-frequency unsteadiness to be linked to the destabilization of the hairpin vortex shedding limit cycle.

I. INTRODUCTION

For many engineering applications, the appearance of largescale unsteadiness represents a design operability limiting factor. The increase of drag and flow-induced vibrations consequent to the unsteadiness onset are not only detrimental for the overall performances but can eventually lead to structural failure due to fatigue. Fundamental studies on three-dimensional (3D) unsteady dynamics associated with axisymmetric bodies, like a sphere, have shown been to be physically relevant to understand more complex configurations, such as launcher after-bodies or re-entry capsules. The incompressible behavior of the space-time dynamics of a flow past a sphere has been well documented in the past both experimentally [START_REF] Taneda | Visual observations of the flow past a sphere at Reynolds numbers between 10 4 and 10 6[END_REF][START_REF] Wu | Sphere wakes in still surroundings at intermediate Reynolds numbers[END_REF][START_REF] Ormières | Transition to turbulence in the wake of a sphere[END_REF][START_REF] Schouveiler | Self-sustained oscillations in the wake of a sphere[END_REF][START_REF] Gumowski | Transition to a time-dependent state of fluid flow in the wake of a sphere[END_REF] and numerically. [START_REF] Natarajan | The instability of the steady flow past spheres and disks[END_REF][START_REF] Tomboulides | Direct and large-eddy simulation of axisymmetric wakes[END_REF][START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF][START_REF] Ghidersa | Breaking of axisymmetry and onset of unsteadiness in the wake of a sphere[END_REF][START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF][START_REF] Fabre | Bifurcations and symmetry breaking in the wake of axisymmetric bodies[END_REF][START_REF] Pier | Local and global instabilities in the wake of a sphere[END_REF][START_REF] Meliga | Unsteadiness in the wake of disks and spheres: Instability, receptivity and control using direct and adjoint global stability analyses[END_REF][START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF] While the interest of the scientific community for the incompressible regime is still very high and numerous studies have been recently done, [START_REF] Skarysz | Experimental investigation of the wake behind a rotating sphere[END_REF][START_REF] Choi | Flow around in-line sphere array at moderate Reynolds number[END_REF][START_REF] Rodriguez | Fluid dynamics and heat transfer in the wake of a sphere[END_REF][START_REF] Eshbal | Measurement of vortex shedding in the wake of a sphere at Re = 465[END_REF][START_REF] Bhowmick | Population distribution in the wake of a sphere[END_REF][START_REF] Yin | On the wake flow behind a sphere in a pipe flow at low Reynolds numbers[END_REF] very limited information exists regarding the effects of compressibility, especially at supersonic speeds.

The scientific community currently seems to agree on the stabilizing effect of the Mach number (M∞ = U * ∞ /c * ∞ with U * ∞ and c * ∞ being the dimensional free-stream velocity and speed of sound, respectively) both at subsonic and supersonic speeds for sphere wakes [START_REF] Meliga | Effect of compressibility on the global stability of axisymmetric wake flows[END_REF][START_REF] Nagata | Investigation on subsonic to supersonic flow around a sphere at low Reynolds number of between 50 and 300 by direct numerical simulation[END_REF][START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF][START_REF] Riahi | A pressure-corrected immersed boundary method for the numerical simulation of compressible flows[END_REF] and other flow configurations as thin wakes [START_REF] Turkyilmazoglu | The absolute instability of thin wakes in an incompressible/compressible fluid[END_REF] and boundary layers on rotating disks. [START_REF] Turkyilmazoglu | Absolute and convective instabilities in the compressible boundary layer on a rotating disk[END_REF] The first (pitchfork or regular) and secondary (supercritical Hopf) unstable bifurcations for an incompressible flow past a sphere are responsible for the loss of spatial axisymmetry and for causing the flow to become unsteady, respectively. Although they are seen to persist in the compressible regime, the corresponding critical Reynolds numbers move toward higher values for increasing Mach number, indicating its stabilizing effect on the wake dynamics. However, no unstable bifurcations at supersonic speeds were found in the Reynolds number (Re = U * ∞ D * s /ν * ∞ with ν * ∞ being the dimensional free-stream kinematic viscosity and D * s being the dimensional sphere diameter) range considered (Re < 600), and the flow system was shown to be globally stable and, more precisely, temporally steady and spatially axisymmetric. Evidence of the existence of unstable bifurcations in the supersonic regime was recently given in a set of direct numerical simulations at higher Reynolds numbers by Nagata et al., [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] who found the flow to first lose its spatial axisymmetry and then become unsteady via the shedding of hairpin vortices downstream of the sphere, like for its subsonic counterpart. The experiments of Nagata et al. [START_REF] Nagata | Experimental investigation of transonic and supersonic flow over a sphere for Reynolds numbers of 10 3 -10 5 by free-flight tests with schlieren visualization[END_REF] at freeflight conditions examined schlieren visualizations of a sphere at 0.9 < M∞ < 1.6 and 3900 < Re < 380 000 and also showed that an unsteadiness existed at Re = 8100 and M∞ = 1.4. However, neither the mechanisms driving the onset of instabilities nor a characterization of the unsteady behavior was given in detail. An interesting and less investigated flow regime that presents wake dynamic instabilities is the transonic one. Riahi et al. [START_REF] Riahi | A pressure-corrected immersed boundary method for the numerical simulation of compressible flows[END_REF] performed calculations of incompressible and compressible spheres with an improved new immersed boundary method and accurately captured different flow configurations. Among these, an unsteady (alternating hairpin) wake was found at M∞ = 0.95 and Re = 600, for which a shock wave forms on the sphere wake. Other relevant studies showing the effect of cooled/heated [START_REF] Nagata | Direct numerical simulation of flow around a heated/cooled isolated sphere up to a Reynolds number of 300 under subsonic to supersonic conditions[END_REF] or rotating [START_REF] Nagata | Direct numerical simulation of flow past a transversely rotating sphere up to a Reynolds number of 300 in compressible flow[END_REF] laminar sphere are worth mentioning for the Mach number range investigated (0.3 < M∞ < 2.0).

Despite some evidence in the literature, a precise identification and characterization of the regular and Hopf bifurcations of a supersonic sphere wake are still lacking. The present work aims at extending the analysis done in the work of Sansica et al., [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] for which no evidence of the unstable bifurcations was found at M∞ = 1.2 and Reynolds numbers up to Re = 380. The current objective is therefore to perform fully 3D global stability analysis on a sphere wake at the same Mach number but in a higher Reynolds number range between Re = 500 and 1300, with the intent to stretch the boundaries of the state-of-the-art knowledge toward an unknown and less documented region of the M-Re plane.

The paper is organized as follows: the problem formulation, comprehensive of the description of the governing equations and stability problem, is presented in Sec. II; code features and numerical simulation setup details are reported in Sec. III; the results obtained on base flow solutions and by global stability analysis around regular and Hopf bifurcations are described in Sec. IV; Sec. IV also includes global stability and dynamic mode decomposition (DMD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] analyses up to Re = 1300, for which a new bifurcation is found; some conclusions and future work opportunities are discussed in Sec. V.

II. PROBLEM FORMULATION

A. Governing equations

The compressible 3D Navier-Stokes (N-S) equations for a perfect gas can be written in the non-dimensional form as

∂q ∂t = N (q), (1)
where q = [ρ, ρu, ρE] T is the state vector in the conservative form (with ρ, u, and E being the fluid density, the velocity vector, and the total energy, respectively) and t is the time. The differential nonlinear N-S operator N can be explicitly expanded as

N (q) = -∇ ⎛ ⎜ ⎝ ρu ρu ⊗ u + pI -τ ρEu + pu -τu + q ⎞ ⎟ ⎠ (2) with p = (γ -1)ρE - 1 2 u ⋅ u, τ = μ[(∇ ⊗ u + ∇ ⊗ u T ) - 2 3 (∇ ⋅ u)I], q h = - μCp Pr ∇T, (3)
where p is the pressure, τ is the viscous stress tensor, Cp is the heat capacity at constant pressure, μ is the dynamic viscosity, Pr is the Prandtl number, T is the temperature, and q h is the heat flux. The Prandtl number is assumed to be constant and, being the fluid considered air, equal to Pr = 0.72. The dynamic viscosity is assumed to follow Sutherland's law as

μ = T 3/2 1 + Ts T + Ts , (4)
where Ts = 110.4 K/T * i,∞ with T * i,∞ being the dimensional freestream stagnation temperature (the superscript * indicates dimensional quantities). The array of the streamwise, vertical, and transverse directions is indicated by x = [x, y, z] T .

The time and length scales are made dimensionless using U * ∞ /D * s and D * s , respectively, where U * ∞ is the free-stream velocity and D * s is the sphere diameter. The dimensionless frequency or Strouhal number (non-dimensionalized by the sphere diameter and free-stream velocity) is defined as St = ω/2π. Similarly to the studies of Nagata et al., [START_REF] Nagata | Investigation on subsonic to supersonic flow around a sphere at low Reynolds number of between 50 and 300 by direct numerical simulation[END_REF][START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] Sansica et al., [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] and Riahi et al., [START_REF] Riahi | A pressure-corrected immersed boundary method for the numerical simulation of compressible flows[END_REF] for the high-Mach and low-Reynolds number regimes considered, the dimensional sphere diameters are in the range D * s = 30 μm -80 μm. A practical application regards gas-particle multi-phase flows in a solid rocket motor at take off, [START_REF] Shimada | Computational fluid dynamics of multiphase flows in solid rocket motors[END_REF] where the alumina particles of size 30 μm-200 μm released in the exhaust supersonic jet could be modeled as tiny spheres. However, the main objective here is to provide a better understanding of the wake dynamics of axisymmetric bodies at supersonic speeds to help unraveling the physics of more complex configurations, such as launcher after-bodies [START_REF] Meliga | Effect of compressibility on the global stability of axisymmetric wake flows[END_REF][START_REF] Rigas | Symmetry breaking in a 3D bluff-body wake[END_REF] or re-entry capsules, [START_REF] Sansica | Global stability analysis of jaxa h-ii transfer vehicle re-entry capsule[END_REF] for which some similarities exist.

B. Stability problem

The stability problem is based upon the use of the linearized N-S equations. The first step to obtain this linearized set of equations is to assume that the nonlinear system in Eq. ( 1) admits an equilibrium solution, q b , defined by N (q b ) = 0 and referred to as the fixed point or base flow. The standard small perturbation technique is used to decompose the instantaneous flow into base flow and small disturbances q(x, t) = q b (x) + εq ′ (x, t), with ε ≪ 1. By assuming that the perturbations are infinitesimal, all nonlinear fluctuating terms are ignored and the linearized N-S equations can be written as

∂q ′ ∂t = L (q ′ ), (5)
where

q ′ = [ρ ′ , ρ ′ u b + ρ b u ′ , ρ ′ E b + ρ b E ′ ]
T is the state vector of conservative perturbation variables and L = ∂N /∂q |q b is the Jacobian operator obtained by linearizing the N-S operator N around the base flow q b . By choosing the normal mode or wave solution q ′ (x, t) = q(x) exp(λt) + c.c (with c.c. indicating the complex conjugate) and being L the discrete form of L, the eigenproblem Lq = λq is obtained. The complex eigenvalue can be split into its real and imaginary parts λ = σ + iω, where σ is the temporal growth rate and ω is the angular frequency. While the angular frequency characterizes the oscillatory behavior, the temporal growth rate indicates whether the equilibrium state bifurcates to another solution. This bifurcation is expressed in a linear framework by the existence of eigenmodes with a corresponding positive growth rate. It is important to stress that no explicit forcing is imposed and the current approach can only detect instabilities that are self-sustained in the numerical domain under exam. Convective instabilities that might grow and/or decay within the domain considered will be "lost" among the modes in the continuous stable branch.

III. NUMERICAL METHOD

An in-house solver is used to compute the compressible N-S equations, both nonlinear and linearized formulations, on multiblock structured grids with a finite-volume approach. To compute the convective fluxes, the Roe scheme [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] is extended to the third order in combination with a MUSCL approach and all flux limiters are set to be inactive. A second order centered scheme is used to differentiate the viscous terms. A dual time stepping method is used to advance the N-S solution in time by an implicit pseudo-unsteady approach. [START_REF] Jameson | Time-dependent calculations using multigrid with applications to unsteady flows past airfoils and wings[END_REF] The steady base flow calculations are computed by setting the Courant-Friedrichs-Lewy (CFL) number equal to 10, allowing the filtering of possible unsteadiness and obtaining converged steady fixed points until the residuals of the state variables in the L 2 -norm are at least 10 -8 . The boundary conditions used are as follows: noslip velocity, adiabatic temperature, and pressure extrapolation on the walls; a uniform velocity is imposed at the numerical domain inflow; characteristic boundary conditions are set at the domain lateral boundaries and outflow to minimize wave reflections. A series of validation test cases for both incompressible and compressible (up to supersonic) flow regimes are performed to verify the correct functioning of the nonlinear solver and are presented in the Appendix.

Some preliminary tests (not shown here) on nonlinear unsteady calculations have been carried out to select the time integration parameters. For both nonlinear and linear unsteady calculations, the values that minimize computational costs and assure a solution independent of the time integration parameters are as follows: the CFL number for the pseudo-unsteady sub-iterations is set to 10, the maximum number of sub-iterations is 35, and the integration time step size is Δt = 9.31 × 10 -3 . The boundary conditions selected for the (nonlinear) steady base flow solutions are also used for the unsteady calculations in both nonlinear and linearized forms.

However, it is important to emphasize that, although the same discretization schemes and boundary conditions are used for both nonlinear and linear calculations, some adaptation to comply with the linearization procedure had to be taken into account. Concerning the spatial discretization, the Roe scheme is based on the Jacobian matrix of the new flux function associated with the linearized equations. [START_REF] Crouch | Predicting the onset of flow unsteadiness based on global instability[END_REF] The boundary conditions are also linearized and modified: zero-velocity perturbations are enforced on the sphere wall and domain inlet, and the characteristic boundary conditions are evaluated on the base flow solution.

A matrix-free method [START_REF] Edwards | Krylov methods for the incompressible Navier-Stokes equations[END_REF][START_REF] Bagheri | Matrix-free methods for the stability and control of boundary layers[END_REF] is used to solve the eigenproblem Lq = λq. It is possible to introduce the exponential propagator M = exp(LΔT) that linearly advances the perturbation solution in time as q ′ (t n+1 ) = Mq ′ (t n ) with t n+1 = t n + ΔT. An Arnoldi algorithm [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF][START_REF] Lehoucq | ARPACK user's guide: Solution of large scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF][START_REF] Barkley | Direct optimal growth analysis for time steppers[END_REF] is coupled to the linear solver [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF][START_REF] Loiseau | Investigation of the roughness-induced transition: Global stability analyses and direct numerical simulations[END_REF][START_REF] Guiho | Instabilities in oblique shock wave/laminar boundary-layer interactions[END_REF][START_REF] Loiseau | Time-stepping and Krylov methods for large-scale instability problems[END_REF] to extract the leading eigenmodes of M. The time step between two consecutive linear snapshots is ΔT = 400 × Δt and ΔT = 200 × Δt for the study of the regular and Hopf bifurcations, respectively. The Krylov subspace dimension is set between 60 and 100 to assure the convergence of the least temporally damped/most temporally amplified eigenmode(s) to be lower than at least 10 -6 .

A. Simulation setup

The numerical domain is composed of six blocks and its topology is displayed in Figs. 1(a) and 1(b). The grid resolution is (nx, ny, nz) = (281, 141, 141) for each block, considering a fine clustering of grid points in the vicinity of the sphere walls and in the wake region. To avoid spurious reflections, the grid is stretched near all lateral boundaries and outflow. Inflow, outflow, and lateral boundaries are set to be Lin = 20 D * s , Lout = 60 D * s , and L lat = 30 D * s , respectively. The center of the sphere is located at (x, y, z) = (0, 0, 0) and its non-dimensional diameter is unitary. The location of boundary conditions and the main characteristic scales are presented in Figs. 1(c) and 1(d). The Reynolds number based on free-stream quantities and sphere diameter spans from Re = 500 to 1300, having evidence of the existence of both regular and Hopf bifurcations in this range. [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] The flow conditions are summarized in Table I. Grid resolution and domain size sensitivity studies for both base flow solutions and global stability results are presented in the Appendix, ensuring the suitability of the present numerical setup.

IV. RESULTS

As a first exploratory campaign, fully 3D unsteady nonlinear simulations have been carried out for Re = 600, 850, and 1000. A probe located just downstream of the sphere at (x, y, z) = (1.0, 0.0, 0.0) is used to record the streamwise velocity time histories for the three cases. Similarly to the work of Nagata et al., [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] while the cases at Re = 600 and 850 are steady, an unsteady periodic behavior at St = 0.134 appears for the Re = 1000 case (further details about the unsteady spatio-temporal dynamics are given in Sec. IV C).

The sphere wakes of these three cases fundamentally differ not only temporally but also spatially: while the wake is axisymmetric for Re = 600, the planar-symmetry characterizes the wakes at Re = 850 and 1000. The zero-streamwise velocity iso-line (black solid line) is plotted over the contours of the density gradient for the two perpendicular x-z (first row) and x-y (second row) planes for the solutions at Re = 600 [(a)-(d)], Re = 850 [(b)-(e)], and Re = 1000 [(c)-(f)] in Fig. 2. The projection of the 3D streamlines onto both planes is also represented (white solid lines). The unsteady solution at Re = 1000 has been time-averaged over five periods (eight samples per period) of the periodic behavior recorded by the streamwise velocity probe. As the selection of the plane where the axisymmetry is lost is equiprobable, for better clarity, the anti-symmetry plane of any solution presented throughout the whole paper has been rotated to match the x-y plane. The bow shock that forms in front of the sphere is rather Reynolds-insensitive and it is steady also for the Re = 1000 case. The separated region downstream of the body

Free-stream Mach number

M∞ = 1.2 Free-stream stagnation temperature T * i,∞ = 287 K Free-stream stagnation pressure p * i,∞ = 1.013 × 10 5 Pa

Reynolds number

Re ∈ [500; 1300] switches from an axisymmetry structure at Re = 600 to a planarsymmetric one at Re = 850 and 1000. In a time-averaged sense, the separated region for Re = 1000 is significantly reduced in length with respect to the steady Re = 850 case. The streamlines upstream of the bow shock are all parallel to the flow direction and spiral axisymmetrically downstream of the sphere for Re = 600. For the Re = 850 and 1000 cases, the streamlines are instead spiraling asymmetrically in the x-y plane and cause the sphere wake to bend upwards. The sphere wake flow at M∞ = 1.2 for these three representative cases, therefore, changes both its temporal and spatial characteristics in this Reynolds number range: the wake is first steady and axisymmetric at Re = 600, then steady and planar-symmetric at Re = 850, and finally unsteady and planar-symmetric at Re = 1000. Similarly to what happens at incompressible [START_REF] Natarajan | The instability of the steady flow past spheres and disks[END_REF][START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF][START_REF] Tezuka | Three-dimensional global linear stability analysis of flow around a spheroid[END_REF][START_REF] Fabre | Bifurcations and symmetry breaking in the wake of axisymmetric bodies[END_REF][START_REF] Meliga | Unsteadiness in the wake of disks and spheres: Instability, receptivity and control using direct and adjoint global stability analyses[END_REF][START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF] and compressible highsubsonic [START_REF] Meliga | Effect of compressibility on the global stability of axisymmetric wake flows[END_REF][START_REF] Nagata | Investigation on subsonic to supersonic flow around a sphere at low Reynolds number of between 50 and 300 by direct numerical simulation[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] speeds, the sphere wake first loses its axisymmetry and then its steady state, indicating the existence of a regular and a Hopf bifurcation, respectively.

To determine and characterize these two bifurcations from a global stability point of view, an appropriate choice of the base flow solutions needs to be done. Considering the type of bifurcations that the fully 3D unsteady nonlinear calculations indicate, it is possible to conclude that for increasing Reynolds numbers, (a) steady axisymmetric base flow solutions become unstable with respect to the regular bifurcation and (b) steady planar-symmetric base flow solutions become unstable with respect to the Hopf bifurcation. In Sec. IV A, different approaches are used to obtain the base flow solutions depending on the nature of the bifurcation to be investigated and their main features are described.

A. Baseflow solutions

Based on the flow features identified by the fully 3D unsteady nonlinear calculations, the two following different strategies for the calculation of the base flow solutions used in the global stability analyses are adopted:

(1) Base flow solutions for the regular bifurcation-steady forced symmetry: The objective is to identify the bifurcation that causes the loss of the wake axisymmetry. In this sense, axisymmetric (and steady) solutions represent the stable/unstable base flows below/above the regular bifurcation critical Reynolds number (Re reg cr ). Thus, a new set of steady axisymmetric base flow solutions is obtained by artificially forcing the flow symmetry taking into consideration only a quarter of sphere in the Reynolds number range Re = 500-750. The numerical domain is therefore longitudinally sliced via two perpendicular cuts passing through the sphere center and applying symmetry boundary conditions to the newly formed lateral boundaries. Although this set of base flow solutions has been obtained on a quarter of domain, it is important to stress that the stability calculations presented in Sec. IV B are fully 3D. Thus, once the steady solution in the quarter of domain has sufficiently converged, a 3D axisymmetric solution must be constructed by revolution.

(2) Base flow solutions for the Hopf bifurcation-steady fully 3D:

The objective is to identify the bifurcation that causes the onset of the unsteadiness in the planar-symmetric solutions.

In this regard, steady (and planar-symmetric) solutions represent stable/unstable base flows below/above the Hopf bifurcation critical Reynolds number (Re Hopf cr ). For this reason, the implicit pseudo-unsteady approach [START_REF] Jameson | Time-dependent calculations using multigrid with applications to unsteady flows past airfoils and wings[END_REF] at CFL = 10 is used to artificially filter any possible unsteadiness and obtain the fully 3D steady base flow solutions in the Reynolds number range Re = 600-1000. Since the steady calculations are fully 3D and no axisymmetry is artificially imposed, these base flow solutions can be either steady axisymmetric or steady planarsymmetric. While only the steady planar-symmetric solutions are used for the characterization of the Hopf bifurcation, this set of calculations can also serve to nonlinearly identify the regular bifurcation critical Reynolds number and verify the corresponding linear calculations. The wake axi-symmetry is, in fact, broken between Re = 645 and 650, as it will be confirmed by the global linear stability results presented in Sec. IV B. II summarize the Reynolds number evolution of the base flow separated region lengths for the two calculation strategies described above. The separation lengths for some selected Reynolds numbers obtained with unsteady fully 3D calculations are also reported. When the wake is planar-symmetric, the separated region length is calculated by operating an average in the azimuthal direction. If an unsteadiness onsets, the separation length is also time-averaged over five periods (eight samples per period) of the periodic behavior recorded by a streamwise velocity probe in the wake region.

Figure 3 and Table

Below Re = 650, the fully 3D approach predicts an axisymmetric solution, suggesting that the regular bifurcation threshold is around this Reynolds number value. Similarly to the incompressible counterpart (as shown in the inset of Fig. 3, where the incompressible results from Sansica et al. [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] are added), the separation length calculated with the fully 3D approach above the regular bifurcation is FIG. 3. Reynolds number evolution of the separation region length for the steady forced symmetry (solid line with full circle symbols), steady fully 3D (dashed line with empty circle symbols), and unsteady fully 3D (dashed-dotted line with empty diamond symbols) approaches. All separation lengths are obtained by averaging in the azimuthal direction. The separation lengths for the unsteady fully 3D calculations are also time-averaged. The inset plot shows the results for the incompressible sphere wake flow, taken from the work of Sansica et al. [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] The same line styles are used to indicate the different wake configurations of the incompressible solutions. shorter (in an azimuthally averaged sense) than the forced axisymmetric solution at the same Reynolds number. Another similarity with the incompressible sphere wake is the fact that the separated region length initially decreases just above the regular bifurcation to then monotonically increase after reaching a minimum. Below Re = 875, the unsteady fully 3D calculations predict steady solutions and indicate that the Hopf bifurcation must be around this Reynolds number value. In agreement with the incompressible counterpart, the span and time-averaged separated region above Re = 875 are significantly shorter than the corresponding steady solutions.

Although the incompressible and supersonic results are qualitatively showing the same trends, some quantitative differences exist. While the non-dimensional separation length (or equivalently the size of the separation with respect to the sphere diameter) is much larger in the supersonic regime, the rate of growth of the unstable base flow separation lengths with the Reynolds number around the bifurcation is higher for the incompressible case. Another important quantitative discrepancy regards the Reynolds number difference between the two bifurcations, which is about four times larger for the supersonic case. These differences again confirm the stabilizing effect of the Mach number on the onset of instabilities.

B. Regular bifurcation

The fully 3D axisymmetric base flows, artificially obtained by the revolution of the steady solutions calculated on a quarter of sphere, are used to investigate the regular bifurcation for the Reynolds numbers Re = 600, 620, 640, 645, 648, 649, 650, 675, and 700. Global stability analysis is performed on these base flows, and the growth rate Reynolds number evolution of the least temporally damped/most temporally amplified non-oscillatory mode is plotted in Fig. 4(a). The shaded area for negative values of σ indicates temporally damped modes. The regular bifurcation threshold can be identified for a critical Reynolds number of Re reg cr = 650. As confirmed by the fully 3D base flow calculations presented in Sec. IV A and in agreement with the work of Nagata et al., [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] this critical Reynolds number value defines the transition from a steady axisymmetric solution to a steady planar-symmetric one. Interestingly, the slope of the growth rate Reynolds evolution around the incompressible regular bifurcation is lower than the supersonic one (see Ref. 24). The full eigenspectrum for the Re = 700 case is presented in Fig. 4(b), where it is possible to see that only the non-oscillatory eigenvalue at (σ, St) = (0.024, 0.0) is temporally amplified. Although not shown here, the modes on the continuous branch are spatially localized in the far wake. Generally, while modes below St = 0.05 preserve axisymmetry, those at higher frequencies have an m = 1 azimuthal character and resemble the mode shape responsible for the Hopf bifurcation (shown later in Sec. IV C).

The corresponding eigenmode is represented in Fig. 5 by plotting the real part of the streamwise, vertical, and transverse eigenvelocities. To better appreciate the azimuthal character of the mode, the streamwise, vertical, and transversal eigenvelocities (û, v, and ŵ, respectively) have been transformed from the Cartesian to cylindrical coordinates to obtain the streamwise, radial, and azimuthal (ûx, ûr, and ûθ , respectively) ones. All eigenvelocities in the present paper have been normalized with the maximum of the real part of the streamwise eigenvelocity. The mode is characterized by an azimuthal distribution corresponding to an m = 1 mode, and the perturbations in the streamwise eigenvelocity extend far downstream of the sphere body. Two "petal-like" structures appear due to the expansions generated on the sphere downstream of the bow shock. It must be emphasized that, except for the petal-like structures appearing on the sphere caused by high compressibility effects, the resemblance of the eigenmode with its incompressible counterpart is evident.

C. Hopf bifurcation

Before describing the global stability results for the Hopf bifurcation, nonlinear unsteady fully 3D simulations are performed at Re = 800, 850, 875, 900, 950, and 1000 to give a qualitative validation against that of Nagata et al. [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] and provide a verification of the corresponding linear calculations. The nonlinear unsteady fully 3D calculations are restarted from the corresponding steady planarsymmetric fully 3D solution and let evolve in time. No forcing has been explicitly added to the simulation, and the unsteadiness slowly builds up from numerical noise. Similarly to the work of Nagata et al., [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] an unsteadiness only appears above Re = 875. The time evolution of the streamwise velocity recorded by a probe in the sphere wake at (x, y, z) = (1.0, 0.0, 0.0) is reported in Fig. 6 calculated during the linear transient phase, the term "nonlinear" refers here to the type of calculation to avoid confusion with the "linear" growth rate calculated by global stability analysis. Figure 6(c) shows a zoom of the linear transient phase in the natural logarithmic scale for the absolute value of the base flow subtracted streamwise velocity time evolution. The calculated nonlinear growth rate is σNL ≈ 0.054. The spatial flow configuration is visualized by the planar and 3D views of the Q-criterion iso-surfaces in Fig. 7, where the instantaneous density gradient contours are also added to the planar views. The Q-criterion iso-surfaces are colored by the streamwise velocity (the light to dark colors correspond to the levels from -0.35 to 1.15). Despite the existence of an axisymmetric bow shock in front of the sphere, these visualizations show that the unsteady dynamics in the sphere wake is characterized by planar-symmetric hairpin vortex shedding. A first important remark must be made on the fact that a very similar behavior characterizes the unsteady incompressible sphere wake above the Hopf bifurcation threshold, both in terms of non-dimensional frequency value and spatial structures.

To describe the Hopf bifurcation and precisely identify the threshold value, global stability analysis is performed for Re = 850, 870, 873, 874, 875, 876, 880, 900, 950, and 1000 on the fully 3D steady planar-symmetric base flows presented in Sec. IV A. The growth rate Reynolds number evolution of the least temporally damped/most temporally amplified oscillatory mode is plotted in Fig. 8(a), where the shaded area for negative values of σ indicates stable base flow cases. The Hopf bifurcation threshold is identified for a critical Reynolds number of Re Hopf cr = 875, defining the transition from steady to unsteady wake dynamics. The full eigenspectrum for the Re = 1000 case is presented in Fig. 8(b), showing that the oscillatory eigenvalue at (σ, St) = (0.054, 0.124) is temporally amplified. First of all, it is useful to remark that the growth rate predicted by linear stability analysis and the one calculated in the transient phase of the nonlinear unsteady calculations at Re = 1000 (σNL ≈ 0.054, see Sec. IV C) are the same, verifying the correct functioning of the linear solver with respect to its nonlinear counterpart. As already seen for other flow configurations or regimes, [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF][START_REF] Barkley | Linear analysis of the cylinder wake mean flow[END_REF] the linear stability analysis predicted frequency is instead slightly underestimated with respect to the nonlinear one. The real part of the eigenvelocities corresponding to the (σ, St) = (0.054, 0.124) mode is visualized in Fig. 9. Similarly to the regular bifurcation, the temporally amplified mode above the Hopf bifurcation in the supersonic regime spatially resembles its incompressible counterpart. The velocity perturbations are packed into planar-symmetric structures whose wavelength is about 3 × Ds and extend in the wake far downstream of the sphere. The planar-symmetry of this mode is not broken, i.e., it has the same spatial symmetry of the base flow. A quantitative difference between this supersonic Mach number and the incompressible counterpart is that the slope of the growth rate Reynolds evolution around the critical instability onset value is lower (see Ref. 24), meaning that for the supersonic case, the planar-symmetric steady state is preserved over a larger range of Reynolds numbers. The appearance of the Hopf bifurcation with respect to the regular one is, in fact, retarded. This retardation can be confirmed by comparing the difference between the Hopf and regular critical Reynolds numbers, being about 60 in the incompressible regime and 225 at M∞ = 1.2.

In the region of temporally damped eigenmodes (σ < 0) presented in Fig. 8(b), one mode at (σ, St) = (-0.021, 0.142) seems to be moving away from the continuous branch and get closer to the temporally amplified region (σ > 0). It, therefore, seems natural to wonder what happens to this mode for a further increase of

FIG. 9.

Global mode at (σ, St) = (0.054, 0.124) for the Re = 1000 case: 2D iso-lines in the x-z and x-y planes (first and second rows, respectively), 3D iso-surfaces (third row), and 2D projections of the iso-lines in the y-z planes at x = 1, 4, 8, 12, 16, 20 (fourth row) of the real part of the streamwise (first column), radial (second column), and azimuthal (third column) eigenvelocities for the levels ûx = ±0.20, ûr = ±0.20, and ûθ = ±0.20, respectively. The white and black iso-lines and iso-surfaces indicate the positive and negative perturbation eigenvelocities, respectively. The base flow density gradient contours (levels between 0.0 and 2.0) are added to the planar visualizations with the zero-streamwise velocity iso-line indicated by the dashed black line.

the Reynolds number. The evolution and the characteristics of this mode are discussed in Sec. IV D.

D. Reflectional symmetry breaking bifurcation

The analysis of the Hopf bifurcation suggests the possibility to have multiple temporally amplified modes for increasing Reynolds number. While it is important to say that the bifurcation following the Hopf one should be analyzed by performing the stability of the limit cycle with a Floquet analysis (i.e., on the unsteady periodic base flow), the objective of the present section is to spatially and temporally characterize the second bifurcation of the steady planarsymmetric base flow solutions. The observability of this instability in the nonlinear dynamics will be addressed in Sec. IV E. For this reason, and with the same principle adopted to obtain the base flow solutions between Re = 600 and 1000 to study the Hopf bifurcation, steady fully 3D calculations are performed for only three other Reynolds numbers, Re = 1100, 1200, and 1300. The corresponding base flow solutions still present the same planar-symmetry found up to Re = 1000, and the azimuthally averaged base flow separation lengths are Lsep = 2.992, 3.083, and 3.171 for Re = 1100, 1200, and 1300, respectively.

These steady planar-symmetric base flow solutions are used to perform global stability analysis. The mode mentioned in Sec. IV C, which at Re = 1000 was moving away from the continuous branch, becomes temporally amplified between Re = 1000 and 1100, indicating the existence of a new bifurcation. In addition, a handful of other modes move away from the continuous branch and tend toward positive values of the growth rate. As well as the mode responsible for the Hopf bifurcation (named M1) and the second mode just found to become temporally amplified at Re = 1100 (named M2), the other modes moving away from the continuous branch are named M3, M4, M5, and M6 based on their vicinity to the temporally amplified region of the spectrum and indicated in Fig. 10(a), where the eigenspectrum for the Re = 1300 case is shown. The Reynolds number evolution of modes M1-M6 is tracked by using dark gray full circle symbols and arrows in Fig. 10(a) but also explicitly plotted in Fig. 10(b). Modes M3-M6 approach more or less rapidly σ = 0, but all remain stable for this Reynolds range, at least up to the highest Reynolds number investigated. While the identification of further bifurcations due to these modes is not intended, it is relevant to look at the corresponding eigenfunctions. Since mode M3 is just a longer wavelength version of M2 and modes M5 and M6 are practically identical to M1 and M2, respectively, only the M1, M2, and M4 modes are plotted in Figs. 11 and12. While mode M1 virtually presents no differences with respect to the mode responsible for the Hopf bifurcation reported in Fig. 9 except for being localized slightly closer to the sphere, it is interesting to see that the newly temporally amplified mode M2 corresponds to a 90 ○ rotation of mode M1 with respect to the x-axis. Following the definition proposed by Fabre, Auguste, and Magnaudet, 11 who refer to the dominant mode that retains the base flow symmetry as "reflectional symmetry preserving" and to the sub-dominant rotated mode as "reflectional symmetry breaking," the newly found bifurcation is named "reflectional symmetry breaking" bifurcation. Although at supersonic speeds, this reflectional symmetry breaking bifurcation confirms the intuition of Fabre, Auguste, and Magnaudet 11 and what

FIG. 13.

Reynolds number evolution of the separation region length for the steady fully 3D (dashed line with empty circle symbols) and unsteady fully 3D (dasheddotted line with empty diamond symbols). All separation lengths are obtained by averaging in the azimuthal direction. The separation lengths for the unsteady fully 3D calculations are also time-averaged.

found by Sansica et al. [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] for incompressible cases, for which an antisymmetric sub-dominant mode becomes temporally amplified if the Reynolds number is sufficiently large. Although M4 is just a longer wavelength version of M1, it was preferred to plot this mode rather than the less temporally damped M3 for comparisons that will be done in Sec. IV E. It should be noted that none of these modes breaks the planar-symmetry of the base flow. Finally, the reader should also keep in mind that M3 and M6 are also rotated of 90 ○ , in the same way as for mode M2.

E. Observability of the reflectional symmetry breaking bifurcation in the nonlinear wake dynamics

The reflectional symmetry breaking bifurcation described in Sec. IV D was found for steady planar-symmetric base flow solutions. However, the bifurcation successive to the Hopf one should be instead sought by performing a limit cycle stability study on the unsteady periodic base flow with a Floquet analysis (out of scope for the present work). The observability of the reflectional symmetry breaking bifurcation in the nonlinear wake dynamics is therefore addressed by performing nonlinear unsteady calculations at Re = 1100, 1200, and 1300. The nonlinear unsteady calculations have been restarted from the corresponding steady planarsymmetric base flow solutions, in the same way done for the nonlinear analysis shown in Sec. IV C for Re = 1000. As the flow features at Re = 1100 and 1200 are very similar to those at Re = 1300, the analysis will mainly focus on the comparisons of the latter with the Re = 1000 case.

At saturation, the nonlinear time-averaged solution is calculated by collecting eight samples per period of the frequency at St = 0.148, over five periods of the lowest peak at St = 0.024. The time and azimuthally averaged separation lengths (dashed-dotted line with empty diamond symbols) are plotted in Fig. 13 with the azimuthally averaged separation lengths of the corresponding base flow solutions (dashed line with empty circle symbols). For increasing Reynolds number, the nonlinear time-averaged separation length tends to reach a minimum, similarly to what found at incompressible regimes. [START_REF] Citro | Stability and sensitivity analysis of the secondary instability in the sphere wake[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF][START_REF] Bouchet | Hydrodynamic forces acting on a rigid fixed sphere in early transitional regimes[END_REF] The time evolution of the streamwise velocity recorded in the sphere wake at (x, y, z) = (1.0, 0.0, 0.0) is reported in Fig. 14(a) for the Re = 1300 case. The initial linear transient is removed up to t ≈ 350 [indicated by the vertical gray dashed line in Fig. 14(a)], and the corresponding PSD distribution, calculated over about 70 periods of the dominant hairpin shedding frequency, is shown in Fig. 14(b). It is possible to see that the spectrum is much richer with respect to the one at Re = 1000 (Fig. 6) and many peaks appear also at lower frequencies [St = O( 10 -2 )]. The mid frequency energy peaks [St = O( 10 -1 )] are in the same Strouhal number range of the linear dominant and reflectional symmetry breaking instabilities and represent possible candidates to be their nonlinear counterparts.

f 0 -3f 1 f 1 f 0 -2f 1 2f 1 f 0 -f 1 3f 1 f 0 f 0 + f 1 2f 0 + 2f 1 f 0 + 2f 1 2f 0 -f 1 2f 0 2f 0 + f 1
The instantaneous Q-criterion iso-surfaces for the Re = 1000 (first column) and Re = 1300 (second column) cases are plotted in the planar views reported in Fig. 15 and colored by the streamwise velocity (the light to dark colors correspond to the levels from -0.35 to 1.15). No large differences can be discerned, except for a slightly shorter wavelength of the shed hairpin vortices for the Re = 1300 case, in agreement with a slightly higher shedding frequency (St = 0.134 for the Re = 1000 case vs St = 0.148 for the Re = 1300 case). It is also important to say that for the Reynolds number range considered, no wake modulation, oscillation, or helical motions were found like at lower Mach numbers. [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] One last comment regards the Reynolds evolution of the problem nonlinearity. While from the spectrum in Fig. 6 for the Re = 1000 case, it is easily recognizable the presence of the fundamental frequency at St = 0.134 and its first (St = 0.268) and second (St = 0.402) harmonics, for the Re = 1300 case, the analysis is relatively more complex. The latter is strongly nonlinear, and it is necessary to distinguish the fundamental frequencies, here f 0 and f 1 , and what results from nonlinear interactions, i.e., the harmonics nf 0 and mf 1 and quadratic or triadic interactions: nf 0 ± mf 1 with n, m being the real positive integers. By taking f 0 = 0.148 (peak "P7" in Fig. 14) and f 1 = 0.041 (peak "P2"), we can, for example, obtain 2f 0 = 0.296 (peak "P12"), f 0 + f 1 = 0.189 (peak "P8"), f 0 -f 1 = 0.107 (peak "P5"), and so on. The nonlinear interactions for all peaks in Fig. 14 for PSD values above 10 -4 are summarized in Table III. By subtracting the mean value and normalizing the oscillation amplitude of the lift-coefficient time-histories, the phase plots [ ĈL(t) vs ĊL(t)] reported in Fig. 16 reinforce the perception of different nonlinear dynamics complexity at Re = 1000 and Re = 1300. The phase space passes from a limit cycle forming a single loop (Re = 1000) to a new limit cycle with three or four loops (Re = 1300). This indicates that the low frequency peak found at St = 0.041 for the Re = 1300 case may be linked to the destabilization of the limit cycle at St = 0.148.

F. Dynamic mode decomposition analysis

Since the Floquet analysis is out of scope for the present work, an alternative methodology to shed some light on the observability of the reflectional symmetry breaking bifurcation can be given by a DMD analysis [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] performed on the unsteady nonlinear calculations. FIG. 18. DMD mode at St = 0.134 for the Re = 1000 case: 2D iso-lines in the x-z and y-z planes (first and second rows, respectively), 3D iso-surfaces (third row), and 2D projections of the iso-lines in the x-y planes at x = 1, 4, 8, 12, 16, 20 (fourth row) of the real part of the streamwise (first column), radial (second column), and azimuthal (third column) eigenvelocities for the levels ûx = ±0.20, ûr = ±0.20, and ûθ = ±0.20, respectively. The white and black iso-lines and iso-surfaces indicate the positive and negative perturbation eigenvelocities, respectively. The time-averaged density gradient contours (levels between 0.0 and 2.0) are added to the planar visualizations with the zero-streamwise velocity iso-line indicated by the dashed black line.

For the Re = 1000 and 1300 cases, 350 snapshots with a temporal interval of 1.057 are collected at saturation. The snapshots are 3D flow fields composed of density, velocities, and pressure for each domain cell and therefore very large. The conventional DMD approach [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] is limited by the large memory required. For this reason, low-dimensionalization is applied to the input dataset. [START_REF] Ohmichi | Preconditioned dynamic mode decomposition and mode selection algorithms for large datasets using incremental proper orthogonal decomposition[END_REF] Eighty one proper orthogonal decomposition (POD) modes are obtained by incremental POD [START_REF] Arora | Stochastic optimization for PCA and PLS[END_REF] and used for the dimensionality reduction. Total-least-squares DMD [START_REF] Hemati | De-biasing the dynamic mode decomposition for applied Koopman spectral analysis of noisy datasets[END_REF] was then applied to the low dimensional datasets. The dominant DMD modes are identified by means of the compressive sensing technique based on the greedy approach. [START_REF] Ohmichi | Numerical investigation of wake structures of an atmospheric entry capsule by modal analysis[END_REF] Further details of the modal analysis framework can be found in the work of Ohmichi, Kobayashi, and Kanazaki, [START_REF] Ohmichi | Numerical investigation of wake structures of an atmospheric entry capsule by modal analysis[END_REF] where the "FBasis" modal analysis tool developed at the Japan Aerospace Exploration Agency (JAXA) is described in detail.

The DMD growth rates are reported in Fig. 17 For the Re = 1000 case, the compressive sensing technique identifies as dominant DMD modes those at frequencies (in order of dominance) St = 0.134, 0.268, and 0.402. These modes correspond to the three peaks found in the PSD distribution. The most dominant mode is the one at St = 0.134, and the other two only represent its first and second higher harmonics, which will therefore not be discussed. The eigenvelocities corresponding to the dominant fundamental DMD mode are reported in Fig. 18. A strong link with the nonlinear dynamics exists, and the footprint of the shed hairpins seems clearer in the DMD eigenfunctions, as visible from the azimuthal eigenvelocity component. One of the main discrepancies between the global stability and the DMD modes consists of the different development of the eigenwake in the streamwise direction. As shown in Fig. 19, both eigenwakes have similar wavelengths and bend upwards following the deflection of the separated region (indicated by the dashed line). While the wake of the global stability mode is confined in a narrow region behind the sphere, the DMD mode diverges while moving downstream of the body. The distortion of the basic flow by nonlinear effects and larger diffusion is at the origin of this spreading effect. The latter is taken into account in the time-averaged flow, around which the DMD modes are calculated. Despite these differences, the similarities with the temporally amplified global mode presented in Fig. 9 are evident. Both the global and the DMD modes are arranged spatially in packets of similar wavelength shape, have the same azimuthal character, and extend very far downstream of the body in the sphere wake region. To further highlight the similarities between the dominant global stability and DMD modes, two movies ("Re1000-GSA-Mode.mp4" and "Re1000-DMD-Mode.mp4" for the global stability and DMD modes, respectively) of the iso-surfaces of the eigenvelocities for the same levels presented in Figs. 9 and 18 (top: ûx, middle: ûr, and bottom: ûθ ) during six periods of their corresponding frequency are provided as the supplementary material.

The compressive sensing technique for the Re = 1300 case indicates as most dominant modes those at (in order of dominance) St = 0.148, 0.041, 0.296, 0.190, and 0.107. While the St = 0.296 mode is just the first harmonic of St = 0.148, the eigenvelocities corresponding to the remaining four DMD modes are reported in Figs. 20 and21(d)] DMD modes present some similarities with the dominant mode but at a higher level of spatial complexity due to nonlinear interactions. All the dominant modes have retained planar-symmetry, and no reflectional symmetry breaking ones are present. To make sure the reflectional symmetry breaking mode is not present among the sub-dominant modes, all DMD modes were carefully checked, but no evidence was found. It can therefore be concluded that the global stability mode responsible for the reflectional symmetry breaking bifurcation is not observable in the nonlinear wake dynamics.

V. CONCLUSIONS

The increase of Mach number has a stabilizing effect on the sphere wake dynamics. [START_REF] Meliga | Effect of compressibility on the global stability of axisymmetric wake flows[END_REF][START_REF] Nagata | Investigation on subsonic to supersonic flow around a sphere at low Reynolds number of between 50 and 300 by direct numerical simulation[END_REF][START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF][START_REF] Riahi | A pressure-corrected immersed boundary method for the numerical simulation of compressible flows[END_REF] At supersonic speeds, no evidence of the regular and Hopf bifurcations, respectively, responsible for the loss of axisymmetry and the appearance of unsteadiness in incompressible sphere wakes, was found in both nonlinear [START_REF] Nagata | Investigation on subsonic to supersonic flow around a sphere at low Reynolds number of between 50 and 300 by direct numerical simulation[END_REF][START_REF] Riahi | A pressure-corrected immersed boundary method for the numerical simulation of compressible flows[END_REF] and linear [START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] frameworks. The Reynolds number must be significantly increased in order for these unstable bifurcations to reappear. [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] A global stability analysis is here carried out for sphere wakes at M∞ = 1.2 in the Reynolds number range Re ∈ [500; 1300]. Steady axisymmetric and steady planar-symmetric base flow solutions are selected to investigate the regular and Hopf bifurcations, respectively. The steady axisymmetric base flows, obtained by forcing the axisymmetry on a quarter of sphere, are found to become unstable at a critical Reynolds number of Re reg cr = 650, identifying a steady planar-symmetric mode to be responsible for the occurrence of a regular bifurcation. The steady planar-symmetric base flows, instead, show the appearance of a Hopf bifurcation at Re Hopf cr = 875, where an unsteady planar-symmetric eigenmode is the cause for the periodic hairpin vortex shedding in the sphere wake. Some differences intrinsic of the highly compressible flow conditions exist with respect to the incompressible counterparts. The supersonic base flows are characterized by the formation of a bow shock in front of the sphere and the presence of expansions on the sphere body. The critical Reynolds numbers at M∞ = 1.2 move toward higher values, and the difference between the two critical thresholds is more than four times bigger than in the incompressible regime. However, remarkably similar features exist between incompressible and compressible (up to supersonic) configurations. The same types of unstable bifurcations persist, and the wake dynamics is analogous in both spatial and temporal senses. This work, therefore, closely links the laminar sphere wake dynamics over a large range of flow velocities, from incompressible to (at least low) supersonic regimes.

At Re = 1000, the global stability analysis shows that a subdominant mode is approaching the temporally amplified region of the spectrum. For this reason, the global stability analysis performed on steady axisymmetric base flows is extended up to Re = 1300. The sub-dominant mode becomes temporally amplified between Re = 1000 and Re = 1100 and its spatial distribution corresponds to a 90 ○ rotation of the dominant mode, confirming some speculations done on the sphere wake instabilities in the incompressible flow regime about the existence of a "reflectional symmetry breaking" bifurcation. [START_REF] Fabre | Bifurcations and symmetry breaking in the wake of axisymmetric bodies[END_REF][START_REF] Sansica | Three-dimensional instability of a flow past a sphere: Mach evolution of the regular and Hopf bifurcations[END_REF] The latter was, however, found by global stability analysis performed around steady planar-symmetric base flow solutions, while to study the bifurcation successive to the Hopf one, a periodic base flow solution should instead be used for a Floquet analysis. Floquet analysis being out of scope for the present work, the observability of this reflectional symmetry breaking bifurcation in the nonlinear wake dynamics is checked by carrying out unsteady nonlinear calculations in the range Re = 1000-1300. The wake spectral energy content interests a wide range of frequencies as the Reynolds number increases, and several PSD peaks, representing potential nonlinear counterparts of the linear instabilities, appear. For both Re = 1000 and Re = 1300, DMD analysis indicates the dominant mode to be related to the nonlinear hairpin vortex shedding and its spatial distribution closely resembles the corresponding global stability mode responsible for the Hopf bifurcation. At Re = 1300, many sub-dominant modes are also detected, but evidence of nonlinear reflectional symmetry breaking modes was not found. The wake dynamics seems rather to be governed by nonlinear interactions that suggest the low frequency unsteadiness to be related to the destabilization of the hairpin vortex shedding limit cycle.

The absence of rotated DMD modes shows the need to perform global stability Floquet analysis on periodic base flows to characterize physical bifurcations successive to the Hopf one. Evidence of wake modulations and wake helical motions superimposed to the hairpin vortex shedding exists at lower speeds, [START_REF] Nagata | Direct numerical simulation of flow past a sphere at a Reynolds number between 500 and 1000 in compressible flows[END_REF] suggesting the possibility to determine subsequent bifurcations and offering future investigation opportunities. It is in the authors' opinions that due to the stabilizing effect of the Mach number (as seen in Ref. 23), the same mechanisms presented here might be preserved at higher speeds, most likely just accompanied by a shift of the critical Reynolds numbers toward higher values. Only for much higher Mach numbers in the hypersonic range, for which the bow-shock forming in front of the body will "wrap" closely around the sphere, new physics might appear. Finally, another regime of interest could be the transonic one, for which shocks will form on the sphere surface or on the sphere wake and may be generating some buffet-like phenomenon.

SUPPLEMENTARY MATERIAL

See the supplementary material for the temporal animations of the global stability ("Re1000-GSA-Mode.mp4") and DMD ("Re1000-DMD-Mode.mp4") modes at Re = 1000 during six periods of their corresponding frequency. The iso-surfaces of the streamwise (top plot), radial (middle plot), and azimuthal (bottom plot) eigenvelocities for the levels ûx = ±0.20, ûr = ±0.20, and ûθ = ±0.20, respectively. The white and black iso-surfaces indicate the positive and negative perturbation eigenvelocities, respectively.

APPENDIX: SOLVER VALIDATION AND SENSITIVITY STUDIES 1. Validation

The same numerical setup and grid described in Sec. III are used for the validation of the nonlinear solver. The results obtained with the present solver are compared over a range of different speeds, from the incompressible regime to the supersonic one. For the incompressible regime, the experimental results of Taneda [START_REF] Taneda | Experimental investigation of the wake behind a sphere at low Reynolds numbers[END_REF] and the numerical ones by Tomboulides, Orszag, and Karniadakis, [START_REF] Tomboulides | Direct and large-eddy simulation of axisymmetric wakes[END_REF] Magnaudet, Rivero, and Fabre, [START_REF] Magnaudet | Accelerated flows past a rigid sphere or a spherical bubble. Part 1. Steady straining flow[END_REF]

Grid resolution and domain sensitivity studies

The sensitivity of both base flow solutions and global stability analysis to grid refinement and domain size is here presented. The reader should note that the six-block elliptical domain topology is maintained for all cases. The domain and grid resolution used to produce the results presented in Sec. IV are indicated as G1-D1.

For the grid refinement study, the domain size is kept fixed and a coarser grid with (nx, ny, nz) = (225, 113, 113) per block (named G0-D1) and a finer grid with (nx, ny, nz) = (337, 169, 169) per block (named G2-D1) are considered. For the domain size sensitivity study, a smaller domain with (Lin, Lout, L lat ) = (10, 30, 15) (named G1-D0) and a larger domain with (Lin, Lout, L lat ) = (40, 120, 60) (named G1-D2) are used. The grid distribution is kept fixed, and elliptical layers of grid points are simply removed or added to account for a smaller or larger domain. For this reason, the grid resolution per block only changes only in the radial direction and is (nx, ny, nz) = (225, 141, 141) and (nx, ny, nz) = (316, 141, 141) for the G1-D0 and G1-D2 cases, respectively. These details of the grid resolution and domain size are summarized in Table V. To evaluate the effect of grid refinement and domain size, the separation length and the least temporally damped/most temporally amplified eigenvalue are chosen as sensitivity indicators for the base flow solution and global stability analysis, respectively. Two Reynolds numbers are selected for each bifurcation in order to consider a stable and an unstable case: Re = 600 and Re = 650 for the regular and Re = 850 and Re = 900 for the Hopf bifurcation. Table VI reports the values of the sensitivity indicators for the examined cases along with the relative percentage error with respect to the values of the G1-D1 case. For all cases, the base flow separation length sensitivity is below 2.5% and the numerical grid/domain G1D1 used as reference throughout the paper can therefore be considered suitable for the base flow

FIG. 1 .

 1 FIG. 1. Assembled (a) and exploded (b) six-block grid topology. Location of boundary conditions (c) and domain characteristic length scales (d) with respect to the sphere diameter.

FIG. 2 .

 2 FIG. 2. Three-dimensional streamlines projected onto x-z [(a)-(c)] and x-y [(d)-(f)] planes superimposed to the density gradient contours for the Re = 600 [(a) and (d)], Re = 850 [(b) and (e)], and Re = 1000 [(c) and (f)] cases. The separated region is indicated by the solid black line.

FIG. 4 .

 4 FIG. 4. Regular bifurcation global stability analysis: (a) Reynolds number evolution of the growth rate of the least temporally damped/most temporally amplified non-oscillatory global mode; (b) eigenspectrum in the St-σ plane for the Re = 700 case. The stable part of the plots has been shaded in light gray.

FIG. 5 .

 5 FIG.5. Global mode at (σ, St) = (0.024, 0.0) for the Re = 700 case: 2D iso-lines in the x-z and x-y planes (first and second rows, respectively), 3D iso-surfaces (third row), and 2D projections of the iso-lines in the y-z planes at x = 1, 4, 8, 12, 16, 20 (fourth row) of the real part of the streamwise (first column), radial (second column), and azimuthal (third column) eigenvelocities for the levels ûx = ±0.05, ûr = ±0.05, and ûθ = ±0.05, respectively. The white and black iso-lines and iso-surfaces indicate the positive and negative perturbation eigenvelocities, respectively. The base flow density gradient contours (levels between 0.0 and 2.0) are added to the planar visualizations, with the zero-streamwise velocity iso-line indicated by the dashed black line.

FIG. 6 .FIG. 7 .

 67 FIG. 6.Unsteady nonlinear calculation at Re = 1000: (a) time evolution of the streamwise velocity in the sphere wake at (x, y, z) = (1.0, 0.0, 0.0) and (b) corresponding power spectral density. A zoom on the linear transient of the natural logarithm of the absolute value of the base flow subtracted streamwise velocity time evolution (c) shows the nonlinear growth rate σ NL .

FIG. 8 .

 8 FIG. 8. Hopf bifurcation global stability analysis: (a) Reynolds number evolution of the growth rate of the least temporally damped/most temporally amplified unsteady global mode; (b) eigenspectrum in the St-σ plane for the Re = 1000 case. The stable part of the plots has been shaded in light gray.

FIG. 10 .

 10 FIG. 10. (a) Eigenspectrum in the Stσ plane for the Re = 1300 case. (b) Reynolds number evolution of the growth rate of the modes leaving the continuous branch and indicated as M1-M6. The stable part of the plots has been shaded in light gray.

FIG. 11 .

 11 FIG. 11.Global stability modes at (first row-M1) (σ, St) = (0.137, 0.128), (second row-M2) (σ, St) = (0.062, 0.137), and (third row-M4) (σ, St) = (-0.016, 0.035) for the Re = 1300 case: 3D iso-surfaces of the real part of the streamwise (first column), radial (second column), and azimuthal (third column) eigenvelocities for the levels ûx = ±0.20, ûr = ±0.20, and ûθ = ±0.20, respectively. The white and black iso-lines and iso-surfaces indicate the positive and negative perturbation eigenvelocities, respectively.

FIG. 12 .

 12 FIG. 12. Global stability modes at (first column-M1) (σ, St) = (0.137, 0.128), (second column-M2) (σ, St) = (0.062, 0.137), and (third column-M4) (σ, St) = (-0.016, 0.035) for the Re = 1300 case: 2D iso-lines projected onto (first row) the x-y planes and (second row) several y-z planes at x = 1, 4, 8, 12, 16, 20 of the real part of the streamwise eigenvelocities for the levels ûx = ±0.20. The white and black iso-lines indicate the positive and negative perturbation eigenvelocities, respectively. The base flow density gradient contours (levels between 0.0 and 2.0) are added to the planar visualizations with the zero-streamwise velocity iso-line indicated by the dashed black line.

FIG. 14 .FIG. 15 .

 1415 FIG. 14. Unsteady nonlinear calculation at Re = 1300: (a) time evolution of the streamwise velocity fluctuations in the sphere wake at (x, y, z) = (1.0, 0.0, 0.0) and (b) corresponding power spectral density.

FIG. 16 .

 16 FIG. 16. Mean value subtracted and amplitude normalized lift-coefficient phase plots for the Re = 1000 (first row) and Re = 1300 (second row) cases. The darker color of the trajectory indicates the time advancement.

FIG. 17 .

 17 FIG. 17. DMD growth rates for the Re = 1000 (a) and Re = 1300 (b) cases represented by full black circles. The corresponding PSDs are also reported as a solid gray line. The modes selected using the compressive sensing technique are indicated by the dark dashedline circles.

  for the Re = 1000 (a) and Re = 1300 (b) cases. The PSD distributions presented in Figs. 6(b) and 14(b) have been added (solid dashed gray line) to facilitate the correspondence of each peak with the DMD modes.

FIG. 19 .

 19 FIG. 19. Re = 1000 case: comparison between the streamwise eigenvelocities of the global stability (top) and DMD (bottom) modes. The base flow (top) and time-averaged (bottom) density gradient contours (levels between 0.0 and 2.0) are added to the planar visualizations with the zero-streamwise velocity iso-lines indicated by the dashed black line.

FIG. 20 .FIG. 21 .

 2021 FIG. 20. DMD modes for the Re = 1300 case at (first row) St = 0.148, (second row) St = 0.041, (third row) St = 0.190, and (fourth row) St = 0.107: 3D iso-surfaces of the real part of the streamwise (first column), radial (second column), and azimuthal (third column) eigenvelocities for the levels ûx = ±0.20, ûr = ±0.20, and ûθ = ±0.20, respectively. The white and black iso-surfaces indicate the positive and negative perturbation eigenvelocities, respectively.

  20 and 21. The St = 0.148 mode [first row in Figs. 20 and 21(a)] is the same found as dominant in the DMD analysis at Re = 1000 and represents the nonlinear counterpart of the dominant linear global instability. The St = 0.041 [second row in Figs. 20 and 21(b)], St = 0.190 [third row in Figs. 20 and 21(c)], and St = 0.107 [fourth row in Figs.

  FIG. 22.Validation for incompressible sphere wake test cases.

  Case nx × ny × nz Lin × Lout × L lat G1-D1 281 × 141 × 141 20 × 60 × 30 G0-D1 225 × 113 × 113 20 × 60 × 30 G2-D1 337 × 169 × 169 20 × 60 × 30 G1-D0 225 × 141 × 141 10 × 30 × 15 G2-D2 316 × 141 × 141 40 × 120 × 60

TABLE I .

 I Inflow parameters.

TABLE II .

 II Separation region lengths and symmetry characteristics for the steady forced symmetry, steady fully 3D, and unsteady fully 3D calculations. AS and PS stand for axisymmetric and planar-symmetric solutions, respectively.

		Steady	Steady	Unsteady
		Forced symmetry	Fully 3D	Fully 3D
	Re	Lsep	Lsep	Lsep
	500	2.620 (AS)	×	×
	550	2.730 (AS)	×	×
	600	2.830 (AS)	2.830 (AS)	2.830 -(AS)
	620	2.865 (AS)	2.865 (AS)	×
	640	2.896 (AS)	2.896 (AS)	×
	645	2.904 (AS)	2.904 (AS)	×
	650	2.910 (AS)	2.900 (PS)	×
	675	2.948 (AS)	2.812 (PS)	×
	700	2.983 (AS)	2.774 (PS)	×
	750	3.045 (AS)	2.743 (PS)	×
	800	×	2.752 (PS)	2.752 (PS)
	850	×	2.779 (PS)	2.779 (PS)
	875	×	2.797 (PS)	2.797 (PS)
	900	×	2.815 (PS)	2.553 (PS)
	950	×	2.858 (PS)	2.330 (PS)
	1000	×	2.901 (PS)	2.220 (PS)

TABLE III .

 III Nonlinear interaction analysis for the Re = 1300 case. The peaks are numbered and indicated in Fig.14.

	Peak no.	1	2	3	4	5	6	7	8	9	10	11	12	13
	Strouhal	0.024	0.041	0.065	0.083	0.107	0.125 0.148	0.188	0.213	0.232	0.256	0.297	0.340

TABLE IV .

 IV Validation for compressible (up to supersonic) sphere wake separation lengths. The results are compared against those of Nagata et al.[START_REF] Nagata | Investigation on subsonic to supersonic flow around a sphere at low Reynolds number of between 50 and 300 by direct numerical simulation[END_REF] (NNT&F) and Johnson and Patel 54 (J and P).Mach Reynolds Present NNT&F Δerr (%) J and P Δerr (%) around 5% or below. A general satisfactory agreement is achieved, and the correct functioning of the nonlinear solver is confirmed.

	0.30	150	1.25	1.14	9.65	1.22	2.46
	0.80	150	1.78	1.74	2.30	×	×
	1.20	150	0.90	0.95	5.26	×	×

TABLE V .

 V Sensitivity studies on grid resolution and domain size. Details of the numerical setup.

DATA AVAILABILITY

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