Galerkin force model for transient and post-transient dynamics of the fluidic pinball

Abstract We propose an aerodynamic force model associated with a Galerkin model for the unforced fluidic pinball, the two-dimensional flow around three equal cylinders with one radius distance to each other. The starting point is a Galerkin model of a bluff-body flow. The force on this body is derived as a constant-linear-quadratic function of the mode amplitudes from first principles following the pioneering work of Noca (On the evaluation of time-dependent fluid-dynamic forces on bluff bodies. PhD thesis, California Institute of Technology, 1997), Noca et al. (J. Fluids Struct., vol. 13, issue 5, 1999, pp. 551–578) and Liang & Dong (Virtual force measurement of POD modes for a flat plate in low Reynolds number flows. AIAA Paper 2014-0054). The force model is simplified for the mean-field model of the unforced fluidic pinball (Deng et al., J. Fluid Mech., vol. 884, 2020, p. A37) using symmetry properties and sparse calibration. The model is successfully applied to transient and post-transient dynamics in different Reynolds number regimes: the periodic vortex shedding after the Hopf bifurcation and the asymmetric vortex shedding after the pitchfork bifurcation comprising six different Navier–Stokes solutions. We foresee many applications of the Galerkin force model for other bluff bodies and flow control.


Introduction
The literature on aerodynamic forces on bodies associated with proper orthogonal decomposition (POD) or any other Galerkin model is surprisingly sparse. On the one hand, force computations are at the heart of engineering fluid mechanics. On the other hand, systematic investigations and interpretations of the aerodynamic force in the Galerkin framework are mostly missing. Considering POD as a linear decomposition of the flow field realizations, Brunton & Rowley (2009) observed that 'While POD modes and the low order model allow for accurate reconstruction of the flow field and preserve Lagrangian coherent structures, it is not clear that this model is directly useful for reconstructing body forces quickly and accurately, since lift and drag forces depend nonlinearly on the flow field, meaning that contributions from different POD modes cannot be added independently'.
The pioneering early work of Noca (1997) and Noca, Shiels & Jeon (1999) reveals that the instantaneous fluid dynamic forces on the body can be expressed with only the velocity fields and their derivatives. Liang & Dong (2014) applied it to the velocity-based POD modes, and derived a force expression in terms of the force of each POD mode and the force from the interaction between the POD modes. The Galerkin force model proposed in this work reveals that any force component is a constant-linear-quadratic function of the mode amplitudes.
The starting point of our investigation is a working Galerkin model based on a low-dimensional modal expansion of an incompressible viscous fluid flow around a stationary body. Intriguingly, mean-field theory (Stuart 1958(Stuart , 1971) was the first foundation of many Galerkin models, building on weakly nonlinear generalizations of stability analyses. Mean-field theory delivered the first derivation of the Landau model (see, e.g. Landau & Lifshitz 1987) for super-and subcritical Hopf bifurcations. The Landau model is experimentally supported for the onset of vortex shedding behind the cylinder wake (Schumm, Berger & Monkewitz 1994;Zielinska & Wesfreid 1995). Generalizations explain the cross-talk between different frequencies over the base flow (Luchtenburg et al. 2009;Shaabani-Ardali, Sipp & Lesshafft 2020), special cases of 'quasi-laminar' interactions foreshadowed by Reynolds & Hussain (1972).
A few decades later, the pioneering wall turbulence POD model by Aubry et al. (1988) allows employing snapshot data far a low-dimensional encapsulation of the Navier-Stokes dynamics. Since then, numerous empirical reduced-order models have been proposed (Rempfer 2000;Kunisch & Volkwein 2002;Rowley, Colonius & Murray 2004;Ilak & Rowley 2006;Bergmann, Bruneau & Iollo 2009;Taira et al. 2017). Control-oriented versions have been developed by , Barbagallo, Sipp & Schmid (2009), Bagheri, Brandt & Henningson (2009), Hinze & Volkwein (2005) and Gerhard et al. (2003). A working Galerkin model can predict the flow and thus the force. Theories for aerodynamic forces have a rich history documented in virtually every fluid mechanics textbook (see, e.g. Panton 1984). There are several force formulae for different cases. Potential flow theory for finite bodies can only explain the force due to accelerations of the body and predicts vanishing drag (d'Alambert paradox). The Zhukovsky formula derives the lift for the potential flow around streamlined cylinders, while the drag computation is still excluded by the d'Alambert paradox. The lifting line theory by Prandtl (1921) extends Zukovsky's formula for finite wings and adds a drag estimate from the created trailing edge vortices. Kirchhoff (1869) laid the first practical foundation for bluff-body drag by allowing for a separation with infinitely thin shear layer. Until today, the drag and lift forces of a body are inferred from the downstream velocity profile (Schlichting & Gersten 2016). These are arguably the most common force theories.
In the Galerkin modelling literature, unsteady forces have been formulated as functions of mode coefficients, as in Bergmann & Cordier (2008) and Luchtenburg et al. (2009). The force formulae are generally calibrated from the reconstructed flow field. Noca et al. (1999) offered an expression for the unsteady forces on an immersed body in an incompressible flow, which only requires knowledge of the velocity field and its time derivative. Based on this idea, Liang & Dong (2014) presented a velocity POD mode force survey method to measure the forces from POD modes on a flat plate. It has shown that the force superposition of each mode of a full POD model can accurately predict the instantaneous forces, and the leading six POD modes are enough to predict the drag force with 5 % error.
In this study, we focus on the unforced 'fluidic pinball', the flow around three equidistantly placed cylinders in cross-flow (Bansal & Yarusevych 2017). Following Chen et al. (2020), the gap distance between the cylinders is chosen to be one radius and the triangle formed by the centres of three cylinders points upstream. This distance allows for an interesting 'flip-flopping' dynamics. The advantage of the fluidic pinball is that already the two-dimensional laminar flow exhibits a surprisingly rich dynamics which has recently been accurately modelled (Deng et al. 2020). As the Reynolds number increases, the flow behaviour changes from a globally stable fixed dynamics to a periodic symmetric vortex shedding after a Hopf bifurcation, to asymmetric vortex shedding after a subsequent pitchfork bifurcation, followed by quasi-periodic and chaotic behaviour. Intriguingly, the post-pitchfork regime with three unstable steady solutions as well as two stable asymmetric limit cycles and one unstable symmetric limit cycle is adequately described by a single five-dimensional Galerkin model. Apparently, the force model for multiple transients of this pitchfork regime is already a challenge.
In the present work, we propose a Galerkin force model for the transient dynamics of the unforced fluidic pinball at different Reynolds numbers. We derive the unsteady forces from the Navier-Stokes equations yielding a constant-linear-quadratic expression of the mode amplitudes of the Galerkin expansion. The consistent form with Liang & Dong (2014) strengthens the theoretical basis of the force expression. Any known symmetric property of the modes is usually considered in the relative modal analysis (Rigas et al. 2014;Podvin et al. 2020), particularly advised for symmetry-breaking instabilities of flows around a symmetric configuration (Fabre, Auguste & Magnaudet 2008;Borońska & Tuckerman 2010). Since the fluidic pinball exhibits a mirror symmetry, we further investigate the force expression under the reflectional (Z2) symmetry. The drag and lift contributions must come from the specific subsets of the constant-linear-quadratic polynomial functions, which is consistent with the drag-and lift-producing modes identified in Liang & Dong (2015).
The manuscript is organized as follows. Section 2 derives the aerodynamic force from a Galerkin model. Section 3 describes the simulation and Galerkin model of the fluidic pinball. In § 4, the force models for the transition of a simple Hopf bifurcation and for the transition of a simple pitchfork bifurcation are discussed. Next, the force model with the elementary modes of two successive bifurcations for the multi-attractor case is investigated in § 5, together with a optimization based on the correction of the mean-field distortion. We summarize the results and outline future directions of research in § 6.

Galerkin force model
In this section, the derivation of a Galerkin force model is described and discussed. Based on the framework of a Galerkin expansion ( § 2.1), the drag and lift forces are expressed as constant-linear-quadratic functions of the mode amplitudes in § 2.2. Alternatively, the forces can consistently be derived from the momentum balance, as elaborated in appendix A. The force model can be further simplified under symmetry considerations in § 2.3.

The Galerkin framework
The fluid flow satisfies the non-dimensionalized incompressible Navier-Stokes equations where p and u are respectively the pressure and velocity flow fields, ν = 1/Re, with the Reynolds number Re. Here, ∂ t , ∇, , ⊗ and · respectively denote the partial derivative in time, the nabla and Laplace operators as well as the outer and inner tensor products. All the variables have been non-dimensionalized, with the cylinder diameter D, the oncoming velocity U, the time scale D/U and the density ρ of the fluid. It is assumed that there exists at least one steady solution (u s , p s ), satisfying the steady Navier-Stokes equations For the Galerkin framework, the space of the square-integrable vector fields L 2 (Ω) is introduced in the observation domain Ω. The associated inner product for two velocity fields u(x) and v(x) reads ( 2.3) The velocity field is decomposed in a basic mode u 0 and a fluctuating contribution. The basic mode may be the steady Navier-Stokes solution u s or the time-averaged flowū. The fluctuation is represented by a Galerkin approximation of N orthonormal space-dependent modes u i (x), i = 1, . . . , N, with time-dependent amplitudes a i (t) where the basic mode u 0 is associated with a 0 ≡ 1 following Rempfer & Fasel (1994). The orthonormality condition reads The Galerkin expansion (2.4) satisfies the incompressibility condition and the boundary conditions by construction. The evolution equation for the mode amplitudes a i is derived by a Galerkin projection of the Navier-Stokes equation (2.1) onto the modes u i ) Ω for the viscous, convective and pressure terms in the Navier-Stokes equations (2.1), respectively. Details are provided by Noack, Papas & Monkewitz (2005). Thus, a linear-quadratic Galerkin system (Fletcher 1984) can be derived, (2.6)

Drag and lift forces on a body
Let Γ be the boundary of the body in the flow domain Ω and n the unit normal pointing outward from the surface element dS. The α-component F ν α (α = x, y, z) of the viscous force vector F ν on the boundary is expressed by where e α is the unit vector in the α-direction and S α,β = (∂ α u β + ∂ β u α )/2 the strain rate tensor with indices α, β = x, y, z. Similarly, the α-component of the global pressure force, exerted on an immersed body, is defined as Without external forces, the viscous and pressure forces in Ω counter-balance the inertial terms provided by the left-hand side of (2.1). The drag force is defined as the projection on e x of the pressure and viscous forces exerted on the body The lift force is similarly defined as the projection on e y of the resulting pressure and viscous forces exerted on the body (2.10) The drag and lift coefficients read Employing the Galerkin approximation (2.4), the viscous force (2.7) can be re-written as  Deng, B.R. Noack, M. Morzyński and L.R. Pastur Similarly, from the pressure Poisson equation the expression for the pressure field is derived as (2.16) The boundary conditions for partial pressures p jk are discussed by Noack et al. (2005). Integrating (2.8) with (2.15) shows that the pressure force is a quadratic polynomial of the a j values Taking the steady solution as the basic mode u 0 = u s with a 0 ≡ 1 implies that a with a i = δ 0i is a fixed point of (2.6) and the total force can be expressed as a constant-linear-quadratic expression in terms of the mode coefficients  A crucial step relies on the choice of the u i modes for the decomposition of (2.4). These could be the POD modes, as usually considered in fluid flows. However, a better choice could be to decompose the flow field on a basis of modes that are becoming active when the system is undergoing a bifurcation. This choice of the so-called bifurcation modes will be investigated in § 3.3.

The Navier-Stokes equations under the Z 2 -symmetry
When the fluid flow configuration exhibits a mirror symmetry, the Navier-Stokes equations (2.1) possess at least one symmetric steady solution (u s , p s ), satisfying (2.2). The Z 2 -symmetry of the velocity and pressure fields, with respect to the (x, z)-plane defined by y = 0, implies where the symmetric components (u s , v s , p s ) ∈ U s and the antisymmetric components (u a , v a , p a ) ∈ U a , U s and U a being respectively the symmetric and antisymmetric subspaces of the system. Other steady solutions can exist, which break the symmetry of the system. We will consider the symmetric steady solution (u s , p s ) as the reference point of (2.1) in the Reynolds decomposition of the flow field as (2.25).
The dynamics under consideration can include transient and post-transient regimes. Here, we introduce the T-averaged flow fieldsū T (x, t) as where T is a time scale to be chosen. When the flow field is oscillating in time, an appropriate choice for T is the period of the local oscillation. The mean flow field is further defined asū and only focuses on the post-transient limit.
When two mirror-conjugated attractors co-exist, it is convenient to introduce the ensemble-averaged flow fieldū • T (x, t) as whereū ± T (x, t) are the T-averaged flow field on the way to each individual attractor. This definition could be readily extended to more than two conjugated attractors. As an ensemble average on mirror-conjugated attracting sets, the ensemble-averaged flow field u • T (x, t) belongs to the symmetric subspace U s . At this point, it is most convenient to introduce the Reynolds decomposition of the flow field, in the form where the mean-field deformation u Δ (x, t) accounts for the distortion of the flow field from the symmetric steady solution u s (x) to the ensemble-averaged flow fieldū • (2.26) The fluctuation flow field u (x, t) has a vanishing time average, meaning that u( subspace U s and u (x, t) to the antisymmetric subspace U a . Thus, a symmetry-based decomposition of (2.1) results in a symmetric and an antisymmetric part, yielding Integrating (2.27a) on the spatial domain Ω, both the left-and right-hand sides yield a time-evolving force vector aligned on e y , while integrating (2.27b) yields a time-evolving force vector aligned on e x . The former is the resulting lift force applying to the boundaries of the fluid domain, while the latter is the drag force. Thus, the Z 2 -symmetry applied to equations (2.20a) and (2.20b) yields where C • D is the drag coefficient of the symmetric steady solution. The vanishing terms in (2.28) can be easily derived from the definitions of q ν α;j and q p α;jk in § 2.2 as As a result, the drag contribution must come from the symmetric subsets of the constant-linear-quadratic polynomial functions, and from the antisymmetric subsets for the lift contribution.

Galerkin model of the fluidic pinball
The force model derived in § 2 is applied to a configuration of three equidistantly placed cylinders in a cross-flow, known as the 'fluidic pinball' configuration (Noack & Morzyński 2017). The flow configuration and the direct Navier-Stokes solver are described in § 3.1. As the Reynolds number is increased, the flow undergoes two subsequent supercritical Hopf and pitchfork bifurcations. The corresponding force dynamics at different Reynolds numbers is reported in § 3.2. The bifurcation modes, newly introduced by Deng et al. (2020), are defined in § 3.3. They provide the orthogonal basis for the Galerkin projection.  axis and aligned with the symmetry axis of the cylinder cluster. All quantities will be non-dimensionalized with the cylinder diameter D, the velocity U ∞ and the unit fluid density ρ. Considering the symmetry of this configuration, a Cartesian coordinate system will be used in the following discussion, with its origin in the middle of the rightmost two cylinders. In this study, no external force will be applied to these three cylinders. A no-slip condition is applied on the cylinders and the velocity in the far wake is assumed to be U ∞ .
Here, the Reynolds number is defined as Re = U ∞ D/ν, where ν is the kinematic viscosity of the fluid. A no-stress condition is applied at the output of the domain. The resolution of the Navier-Stokes equations (2.1) is based on a second-order finite-element discretization method of the Taylor-Hood type (Taylor & Hood 1973), on an unstructured grid of 4225 triangles and 8633 vertices and an implicit integration of the third-order in time. The instantaneous flow field is calculated with a Newton-Raphson iteration until the residual reaches a tiny tolerance prescribed. This approach is also used to calculate the steady solution, which is derived from the steady Navier-Stokes equations (2.2). The direct Navier-Stokes solver used herein has been validated in Noack et al. (2003) and Deng et al. (2020), with a detailed technical report (Noack & Morzyński 2017). The grid used for the simulations was shown to provide a consistent flow dynamics, compared to a refined grid, see Deng et al. (2020).

Flow features and the corresponding force dynamics
Different from Deng et al. (2020), where the viscous contribution to the forces has been ignored, the lift C L and drag C D coefficients are here calculated from the resulting force F of pressure and viscous components exerted on the three cylinders.
The flow characteristics depend on the Reynolds number Re. Following the literature on clusters of cylinders (Chen et al. 2020), the characteristic length scale is chosen to be the cylinder diameter D and not the transverse width 5D/2 of the configuration. This width loses its dynamic significance for the large distances considered in other studies.
For Reynolds numbers Re < Re H ≈ 18, the symmetric steady solution u s (x) was found to be stable and is the only attractor of the system. A supercritical Hopf bifurcation occurs at Re = Re H , associated with the cyclic release of counter-rotating vortices in the wake of the three cylinders from the shear layers that delimit the configuration, forming a von Kármán street of vortices. The corresponding Reynolds number based on the transverse width of the fluidic pinball is 45, i.e. is well aligned with typical onsets of vortex shedding behind bluff bodies. For the critical value Re = Re PF ≈ 68, the system undergoes a supercritical pitchfork bifurcation. As a result, two additional  (unstable) steady solutions occur, namely u + s (x) and u − s (x), which break the reflectional symmetry of the configuration, as shown with the lift coefficients of the steady solutions in figure 2. The mean field inherits the asymmetry of the steady solutions, with the jet between the two downstream cylinders being deflected upward or downward. As reported in Deng et al. (2020), at Re = Re PF , the statistically symmetric limit cycle, associated with the statistically symmetric vortex shedding, becomes unstable with respect to two mirror-conjugated statistically asymmetric limit cycles, associated with statistically asymmetric von Kármán streets of vortices. Figure 3 shows the time evolution of the lift and drag coefficients at Re = 80, when the initial condition is either the symmetric steady solution u s (figure 3a) or the asymmetric steady solution u + s (figure 3b). In both cases, the asymptotic regime is the same. However, when starting from the symmetric steady solution in figure 3(a), a long-living plateau of the drag coefficient is reached around t ≈ 775, which corresponds to the transient exploration of the unstable limit cycle, centred on the symmetric T-averaged flow field u 98 (x, 775). Note that, during the transient dynamics from the steady solution to the is the drag coefficient of the symmetric steady solution at the corresponding Reynolds number. unstable limit cycle, the drag coefficient is monotonically increasing, before reaching the transient plateau. The drag coefficient is further increased when leaving the unstable limit cycle towards the asymptotically stable limit cycle, the latter being centred on the asymmetric mean flow fieldū + . Figure 4 shows another representation of the transient dynamics for Re = 30, 80 and 100, starting from different initial conditions in the plane ( being the drag associated with the symmetric steady solution at the Reynolds number under consideration. The black cross (×) stands for the symmetric steady solution u s while the asymmetric u + s and u − s steady solutions are respectively represented by a red circle and a blue square, when they exist, at Re = 80 and 100. As can be observed in this figure, different from what happens at Re = 80, the transient dynamics from the symmetric steady solution at Re = 100 first reaches one of the two asymmetric steady solutions, before evolving toward the stable attracting limit cycle.

The bifurcation modes of the fluidic pinball
In the case of two subsequent supercritical Hopf and pitchfork bifurcations, Deng et al. (2020) have shown that the reduced-order model must comprise 5 modes (3.1) Hence, in the decomposition of (2.4), the number of modes is restricted to N = 5. For dynamic interpretability, the basic mode u 0 (x) is chosen to be the symmetric steady , of the velocity field associated with the five elementary degrees of freedom solution u s (x). The first three modes u 1,2,3 (x) are associated with the Hopf bifurcation, the last two modes u 4,5 (x) with the pitchfork bifurcation. We will refer to these modes as the irreducible bifurcation modes of the system. Modes u 3 (x) and u 5 (x) are symmetric. The instability-related modes u 1,2 (x) and u 4 (x) are antisymmetric. Modes u 1,2 (x) span the subspace associated with the limit cycle of the Hopf bifurcation, while u 4 (x) accounts for the symmetry breaking of the pitchfork bifurcation. In Deng et al. (2020), modes u 1,2 (x) are provided by the first two dominant POD modes, while mode u 4 (x) is defined as where u ± s (x) are the two additional (asymmetric) steady solutions arising from the supercritical pitchfork bifurcation. Mode Here,ū(x) will be restricted to the symmetric mean flow field, associated with the statistically symmetric limit cycle, whether this limit cycle is stable or unstable. Similarly to u 4 (x), mode u 5 (x) is defined as These two modes, together with modes u 1,2,3 , are shown in figure 5 after ortho-normalization by a Gram-Schmidt procedure, and the corresponding time-dependent amplitudes a i (t), i = 1, . . . , 5, in the full-flow dynamics are shown in figure 6 when starting from either the symmetric steady solution (figure 6a) or the asymmetric steady solution (figure 6b).

Galerkin force model associated with the supercritical Hopf and pitchfork bifurcation
As already mentioned, the fluidic pinball undergoes a supercritical Hopf bifurcation at Re = Re HP and a subsequent supercritical pitchfork bifurcation at Re = Re PF > Re HP . The Galerkin force models are derived for the supercritical Hopf bifurcation in § 4.1 and for the supercritical pitchfork bifurcation in § 4.2.

Force model associated with the supercritical Hopf bifurcation
The symmetric steady solution u s ∈ U s is stable at low Reynolds numbers. At Re Re HP , it undergoes a supercritical Hopf bifurcation. The resulting Galerkin expansion reads (4.1) and the corresponding mean-field Galerkin system with σ = σ 1 − βa 3 and ω = ω 1 + γ a 3 , where σ 1 and ω 1 are the initial growth rate and frequency depending on the Reynolds number. For a direct supercritical Hopf bifurcation, σ 1 , ω 1 , β > 0, σ 3 < 0 and β 3 > 0. We refer to Deng et al. (2020) for details. Introducing (4.1) in (2.7) and (2.8), the total force can be written as (2.18) with N = 3 degrees of freedom. From symmetry considerations, as u 1,2 ∈ U a and u 0,3 ∈ U s , the coefficients l x;1 , l x;2 , q x;13 , q x;23 , l y;0 , l y;3 , q y;11 , q y;12 , q y;22 , q y;33 are vanishing. Finally, the drag formulae (2.28) simplify to C D = C • D + l x;3 a 3 + q x;11 a 2 1 + q x;12 a 1 a 2 + q x;22 a 2 2 + q x;33 a 2 3 , (4.3a) C L = l y;1 a 1 + l y;2 a 2 + q y;13 a 1 a 3 + q y;23 a 2 a 3 . (4.3b) Here again, C • D is the drag coefficient associated with the symmetric steady solution. The unknown parameters in the force model can be identified by a least-squares approach, according to the known force dynamics and the relevant mode amplitudes. However, for the mean-field Galerkin system (4.2), the slaving relation between the degree of freedom a 3 to the oscillating degrees of freedom a 1 , a 2 imposes an additional sparsity in the force model. We employ the SINDy (sparse identification of nonlinear dynamics) algorithm (Brunton, Proctor & Kutz 2016) to arrive at simpler and more interpretable models. A L1-regularization can be introduced in the LASSO (least absolute shrinkage and selection operator) regression process, which includes the L1-norm of the vector of coefficients in calculating the loss function. Another option in the SINDy algorithm is the sequential thresholded least-squares regression, which iteratively applies the least-squares regression and eliminates terms with weight smaller than a given threshold. Both regression algorithms benefit from simplicity, only requiring one sparsity parameter λ. The optimal λ balances the accuracy and complexity of the identified model. To evaluate the performance of the identified model, the complexity is presented with the number of non-zero coefficients and the accuracy by the coefficient of determination,denoted as the r 2 score (Draper & Smith 1998). A detailed review of this sparsity parameter can be found in   physical constraints of energy-preserving quadratic nonlinearities successfully identifies the sparse model, benefiting from the Galerkin projection of the Navier-Stokes equations (Loiseau, Noack & Brunton 2018). The LASSO algorithm is applied to a scenario starting with the unstable symmetric steady solution at Re = 30. The training data used for the sparse regression are provided by the force coefficients and the mode amplitudes from the direct numerical simulation (DNS) starting with the symmetric steady solution to the final asymptotic regime. The resulting transient dynamics and the asymptotically attracting limit cycle are shown in the three-dimensional space of the time-delayed coordinates of C L and C D in figure 7.
The possible over-fitting terms, such as the slaving relation between a 3 and a 2 1 , a 2 2 , can be suppressed with a larger L1-penalty parameter for the LASSO algorithm. The choice of the L1-penalty parameter drives the sparsity of the identified model. A too small L1 will lead to a complex model with few eliminated terms; on the contrary, a too large L1 can jeopardize accuracy. Both cases weaken the robustness of the identified model, and the same is observed for the sequential thresholded least-squares regression. The influence of the sparsity parameter λ and the comparison of these two regression methods are presented in appendix B.
Gradually increasing the L1-penalty from 0 to nearly 1, the terms a 1 a 2 , a 3 , a 2 2 , a 2 1 are eliminated subsequently in the drag model, while a 2 3 is always retained. The sparsity parameter λ, here the L1-penalty, is chosen as the largest value without any known over-fitting term. Hence, according to the order of elimination, a 3 is the over-fitting term in the drag model due to the slaving relation between a 3 and a 2 1 , a 2 2 . The details of this choice can be found in appendix B. Finally, the identified force model reads The force model is highly accurate, as corrobororated by the r 2 scores of 0.9991 and 0.9942 for the drag and lift formulae, respectively. As shown in figure 8, the dynamics of the force model compares well with the real force transient dynamics, starting from the symmetric steady solution at Re = 30. In the drag model (4.4a), the coefficient of a 3 is vanishing. Mode u 3 actually contributes to the increase of the drag through a 2 3 , as evidenced by the positive coefficient of the a 2 3 term. This is an interesting result, since the effect of the bifurcation mode u 3 is to decrease the length of the recirculation bubble in the T-averaged flow fieldū , resulting in an increase of the drag through the quadratic term a 2 3 . This quadratic dependency is also reported in .
It is also worth noticing that a 2 3 contributes to the mean value of C D while a 2 1 , a 2 2 accounts for the instantaneous oscillations of C D , as C D oscillates at twice the vortex shedding frequency. For C L , the oscillatory pair (a 1 , a 2 ) fits well with the phase of the initial transient part, while the pair (a 1 a 3 , a 2 a 3 ) resolves the phase dependency of the post-transient part of the dynamics.

4.2.
Force model associated with the supercritical pitchfork bifurcation Next, we consider the supercritical pitchfork bifurcation, which breaks the symmetry of the symmetric steady solution u s at Re Re PF . In this case the antisymmetric mode u 4 describes the antisymmetric instability, which corresponds to an unstable eigenmode with a real eigenvalue. The resulting Galerkin expansion reads (4.5) and the corresponding mean-field Galerkin system where σ 4 is the positive initial growth rate, which depends on the Reynolds number. For a direct supercritical pitchfork bifurcation, σ 4 , β 4 > 0, σ 5 < 0 and β 5 > 0, see Deng et al. (2020) for details. Substituting (4.5) in (2.7) and (2.8), with N = 2 in (2.18), and with u 4 ∈ U a and u s , u 5 ∈ U s , the force model becomes C D = C • D + l x;5 a 5 + q x;44 a 2 4 + q x;55 a 2 5 , (4.7a) C L = l y;4 a 4 + q y;45 a 4 a 5 . (4.7b) Five parameters, namely l x;0 , l x;5 , q x;44 , q x;55 , l y;4 , q y;45 need to be identified. In the fluidic pinball, the pitchfork bifurcation occurs after the primary Hopf bifurcation as the Reynolds number is increased. However, the transient dynamics observed at Re = 100, when starting close to the symmetric steady solution, first exhibits the static symmetry breaking, which is typical of the pitchfork bifurcation, before developing the cyclic release of vortices, which is characteristic of the Hopf bifurcation. The early stage of the transient dynamics, starting from the symmetric steady solution and evolving toward one of the asymmetric steady solutions, is shown in figure 9. The time evolutions of the lift C L (t) and drag C D (t) coefficients are shown in figure 10.
Only the degrees of freedom associated with the pitchfork bifurcation are active in this early stage of the transient dynamics, as also shown in figure 17(a). The degrees of freedom associated with the Hopf bifurcation will only become active further in time during the transient dynamics, which will be further discussed in § 5.4. Accordingly, a force model is derived for the transition after a simple pitchfork bifurcation. The training data are the lift C L (t) and drag C D (t) coefficients and the relevant mode amplitudes in (4.7) from the early to final stage of the transient dynamics. The observed slaving of a 5 in a 2 4 may reduce the robustness of the identified model. Gradually increasing the L1-penalty parameter in the LASSO regression, the optimized force model reads C D = 3.58248992 + 0.04367604a 5 − 0.08525184a 2 5 , (4.8a) with r 2 = 0.9949 for the drag model and r 2 = 0.9992 for the lift model. The over-fitting term a 2 4 has been eliminated in the sparse formula of the drag force. Note that the mode u 5 contributes to the drag through a 5 , while a 2 5 acts in decreasing the drag, as indicated by the sign of their associated coefficients in (4.8a). Figure 10 compares the evolutions of the drag and lift coefficients in the full-flow dynamics (solid black line) to their prediction by the force model (4.8) (red dashed curve), during the early stage of the transient dynamics at Re = 100. The derived force model is well aligned with the real force dynamics using only two active degrees of freedom of the pitchfork bifurcation in the dynamics of the system.

Galerkin force model for multiple invariant sets
We focus on the regime after the pitchfork bifurcation Re Re PF = 68 and before the quasi-periodic behaviour Re Re QP = 104. This flow has 6 invariant sets: 3 unstable fixed points, 2 stable asymmetric mirror-conjugated periodic orbits and one meta-stable symmetric limit cycle. Section 5.1 investigates the dynamics of the fluidic pinball at Re = 80, when the degrees of freedom associated with the Hopf bifurcation are first activated before the degrees of freedom associated with the pitchfork bifurcation. The predictive power of the force model is assessed in § 5.2. Section 5.3 introduces two additional degrees of freedom in the force model, in order to take into account the distortion of the shift mode when the attractor is reached. The robustness of the force model is emphasized in § 5.4 by considering the flow dynamics at Re = 100, where the pitchfork degrees of freedom are activated before the Hopf degrees of freedom during the transient dynamics.

Force model at Re = 80
At Re = 80, the system has already undergone a supercritical Hopf bifurcation and a supercritical pitchfork bifurcation. The trajectories issued from u s and u ± s are shown in the time-delayed embedding state space (C L (t), C L (t − τ ), C D (t)) of figure 11. The force model will rely on five degrees of freedom at minimum, namely the three degrees of freedom associated with the Hopf bifurcation a i , i = 1, 2, 3 and the two degrees of freedom a i , i = 4, 5, associated with the pitchfork bifurcation. As a generalization of (4.3) and (4.7), the force model reads C D = C • D + l x;3 a 3 + q x;11 a 2 1 + q x;12 a 1 a 2 + q x;22 a 2 2 + q x;33 a 2 3 + l x;5 a 5 + q x;44 a 2 4 + q x;55 a 2 5 + q x;14 a 1 a 4 + q x;24 a 2 a 4 + q x;35 a 3 a 5 , (5.1a) C L = l y;1 a 1 + l y;2 a 2 + q y;13 a 1 a 3 + q y;23 a 2 a 3 + l y;4 a 4 + q y;45 a 4 a 5 + q y;15 a 1 a 5 + q y;25 a 2 a 5 + q y;34 a 3 a 4 . (5.1b) Due to symmetry reasons, only 2 linear terms (a 3 , a 5 ) and 9 quadratic terms (a 2 1 , a 1 a 2 , a 1 a 4 , a 2 2 , a 2 a 4 , a 2 3 , a 3 a 5 , a 2 4 , a 2 5 ) are left in (5.1a) for the drag coefficient. For the lift coefficient, only 3 linear terms, a 1 , a 2 , a 4 , and 6 quadratic terms, a 1 a 3 , a 1 a 5 , a 2 a 3 , a 2 a 5 ,  a 3 a 4 , a 4 a 5 , are left in (5.1b). The training data are taken from the DNS starting from the three steady solutions, with the real force dynamics, see the black curves in figure 12, and the relevant mode amplitudes, see figure 6. The coefficients of the force models are identified by the sequential thresholded least-squares regression with the optimal sparsity parameter λ. We note that the LASSO regression can also be used here. See appendix B for the comparison of these two methods. The resulting force model reads The good accuracy of the identified drag model can be determined from the high r 2 score of 0.9816. The drag model of (5.2a) preserves both the basic forms of the drag model for the Hopf and pitchfork bifurcations and the signs of the coefficients. This indicates that the identified model is robust. The only remaining cross-term a 3 a 5 provides the coupling relation between the degrees of freedom associated with both bifurcations. A robust sparse formula for the lift model is more difficult to derive, due to the oscillating dynamics of the lift and the fact that a 4 and a 5 also oscillate at the fundamental frequency. With respect to the basic lift model of two bifurcations, a balanced method is used here to solve the difficulty of the identification. Starting with a large L1-penalty, the derived under-fitted system can figure out the most elementary features of the dynamics, eliminating a 1 a 5 , a 2 a 5 , a 4 a 5 . This is reasonable as a 3 is approximately ten times larger than a 5 , which means that most of the mean-field distortion comes from u 3 . However, if the L1-penalty is too large, the term a 4 a 5 can disappear from the lift model, making the resulting model inconsistent with (5.2b). In order to balance sparsity and robustness, a 4 a 5 needs to be reintroduced into the library. The sparse formula of the lift model in (5.2b) is determined by least-squares regression, constraining the parameters of a 1 a 5 , a 2 a 5 to zero. The r 2 score of the identified lift model is 0.9673.
The identified force dynamics in (5.2) (dashed red line) is compared to the real force dynamics (solid black line) at Re = 80 in figure 12. The force model based on the least-order model can reproduce the main features of the real force dynamics. The drag model of (5.2a) shows how the degrees of freedom of the Hopf (a 2 1 , a 2 2 , a 2 3 ) and pitchfork 918 A4-18 (a 5 , a 2 5 ) bifurcations contribute to the drag force, as well as the coupling between these degrees of freedom (a 3 a 5 ). The lift model of (5.2b) shows that the lift oscillations occur through the coupling of the oscillating degrees of freedoms a 1 , a 2 to a 3 , while the coupling between the degrees of freedoms a 4 and a 5 contribute to the mean value of C L . Hence, the mean lift coefficient can be simplified with fewer terms, asC L = l y;4 a 4 + q y;45 a 4 a 5 + q y;34 a 3 a 4 , which is consistent with the Krylov-Bogoliubov assumption (Jordan & Smith 1999).

Assessing the predictive power of the force model
The time evolution of the drag and lift coefficients in the fluidic pinball are shown in figure 12 as solid black lines. The evolutions of the drag and lift coefficients in the model (5.2) are shown with dashed red lines. The model reproduces correctly the time scales of the force dynamics as well as the transient and asymptotic amplitudes of the forces. However, it is observed that the fine details of the transient dynamics, at the early stage of the linear instability, are not satisfactorily reproduced in the identification process (figure 12(a,b) at t ≈ 590 and 475 respectively). The ranges of time concerned, in both cases, are also associated with oscillations in a 4 , as observed during the initial stage at t ≈ 590 in figure 6(a) and t ≈ 475 in figure 6(b). This strongly suggests that the oscillations of a 4 be triggered by the degrees of freedom associated with the Hopf bifurcation. This means that the degrees of freedom of the pitchfork bifurcation are affected by the degrees of freedom of the Hopf bifurcation, at least when the distance from the bifurcation point is large enough, which is the case at Re = 80.
In addition, as recalled in § 3.3, at Re ≈ 68, both the steady symmetric solution and the symmetric-based limit cycle undergo a supercritical pitchfork bifurcation. We emphasize that this coincidence of two local pitchfork bifurcations might not occur by chance, as mentioned in Deng et al. (2020). As a result of these two simultaneous bifurcations, the degrees of freedom involved in the pitchfork bifurcation of the fixed point might not coincide with those involved in the pitchfork bifurcation of the limit cycle. For this reason, it is reasonable to introduce two distinct sets of degrees of freedom for each of them, namely a 4 , a 5 at the fixed point and a 6 , a 7 at the limit cycle. These two additional degrees of freedom will complete the mean-field model with more details and will take into account the mean-field distortion during the transition from the fixed point to the limit cycle. The new resulting mean-field Galerkin system, with seven degrees of freedom, is derived in appendix E, while the new resulting force model is discussed in the next subsection.

The need for additional modes
All our attempts to smooth out the kinks observed at the beginning of the exponential growth, in both C D and C L in the frame of the force model (5.2), failed, even when over-fitting the model without any sparsity. This strongly indicates that five degrees of freedom might not be sufficient to account for the force evolution on the full time range.
Digging into this idea, it becomes manifest that the way u 3 (x) is built, namely as the difference between the statistically symmetric mean flow field, associated with the unstable limit cycle, and the symmetric steady solution u s (x), 3) see figure 12(a), does not allow us to satisfactorily account for the complete dynamics of the lift and drag forces. This also indicates that u 3 (x) gets distorted when the system is evolving along the manifold, which connects the unstable limit cycle to one of the two conjugated stable limit cycles. In other words, the mean-field distortion on the attractors associated with the two asymmetric mean flow fieldsū ± , namely do not coincide exactly with u 3 (x). The asymmetric mean flow fieldsū ± only focus on the post-transient dynamics, as shown in figure 12(c), which can be expressed withū ± T (x, 700). The difference between u ± 3 (x) and u 3 (x) is asymmetric and can be decomposed into a symmetric and an antisymmetric part, respectively u 6 (x) and u 7 (x): As a result, the modes u 6 (x) and u 7 (x) can be defined as, After orthogonal normalization by a Gram-Schmidt procedure, the resulting modes are shown in figure 13, with their mode amplitudes in figure 14. When comparing the definitions of u 6 and u 7 in (5.6) and of u 4 and u 5 in (3.2) and (3.3), it is not surprising that the spatial structure of u 6 , respectively u 7 (see figure 13), be so similar to the spatial structure of u 4 , respectively u 5 (see figure 5). We also mention that, u 6 , u 7 as defined in (5.6), would be equivalent to the pitchfork modes u 4 , u 5 built on the periodic solutions instead of being built on the steady solutions. However, after the Gram-Schmidt procedure, the u 6 , u 7 modes of figure 13 have been transformed into corrective modes of u 4 , u 5 when departing from the steady solutions and approaching the asymptotic limit cycles. The corrective modes u 6 , u 7 should be slaved to u 4 , u 5 along the mean-field distortion of u 3 . The corresponding slaving relation will not be discussed in this paper. Hence, the combination of u i , i = 4, . . . , 7, works as a flexible pitchfork mode expansion, which adapts the whole phase space where all the invariant sets (steady/periodic) locate.
In figure 14, the transient dynamics of a 6 , a 7 is shown to be also similar to a 4 , a 5 in figure 6. Not surprisingly, the opposite initial bump of a 6 , a 7 helps to better fit the dynamics on the manifold. Besides, a 6 , a 7 show no contribution close to the steady solutions, as their role is to adapt the modes u 4 , u 5 when approaching the stable limit cycle. The force model identification is more challenging with these two additional modes. High robustness is required for our force model without losing the identified terms in § 5.1. Compared to the force formula (5.1) with five modes, 8 new terms are introduced in the drag formula, namely a 7 , a 1 a 6 , a 2 a 6 , a 4 a 6 , a 2 6 , a 3 a 7 , a 5 a 7 , a 2 7 , and 7 additional terms are considered in the lift formula, namely a 6 , a 3 a 6 , a 5 a 6 , a 1 a 7 , a 2 a 7 , a 4 a 7 , a 6 a 7 . Due to the similar transient dynamics of a 4 , a 5 and a 6 , a 7 , the corrective degrees of freedom a 6 , a 7 can easily replace a 4 , a 5 in the identified model. Hence, the original structure of the force model with five modes could be lost. To avoid possible over-fitting, we need to free the active terms gradually and constrain the parameters of a 4 , a 5 during the sparse regression to ensure the robustness of the result. In addition, the newly introduced terms should work as a corrective function to the original force model with five degrees of freedom. In other words, the new force model with seven degrees of freedom should inherit the original structure of (5.2).
Based on the structure of the drag model (5.2a), the terms a 7 , a 7 a 7 , a 3 a 7 , a 4 a 6 and a 5 a 7 are introduced in the extended model. The terms a 1 a 6 , a 2 a 6 , a 2 6 are firstly set to zero because their corresponding terms a 1 a 4 , a 2 a 4 , a 2 4 in (5.2a) are vanishing. In order to improve the robustness of the regression results, the terms a 5 and a 2 5 are constrained with the values from (5.2a). Increasing the L1-penalty of the LASSO regression, l x;7 , q x;46 and q x;77 vanish successively, and an obvious under-fitting starts when losing q x;35 . The introduced terms q x;37 , q x;57 are robust with little possibility of over-fitting. Eventually, the drag model reads C D = 3.77331204 + 0.05888312a 5 − 0.00169970a 2 1 − 0.00156775a 2 2 + 0.00513885a 2 3 + 0.00786294a 3 a 5 + 0.00950204a 3 a 7 − 0.25910470a 2 5 − 0.06264888a 5 a 7 .
(5.7) Equation (5.7) preserves the original form of (5.2a), with tiny changes to the coefficients. This extended model fits well the dynamics of the drag coefficient, with the r 2 score increasing to 0.9981, which also can be seen with the red dashed curve of figure 15(a,b). As already mentioned, the drag monotonically increases with the development of the vortex shedding. This is obvious, for instance, from figure 15, when the lift starts to oscillate and the drag to increase. The positive signs of q y;33 , q y;35 and q y;37 , in the drag model of (5.7), are responsible for this monotonic increase of the drag. Compared to the drag model with only a 5 in § 5.1, the contribution to the drag of a 5 and a 7 is more subtle. They contribute to an increase of the drag through a 5 , a 5 a 3 and a 7 a 3 , while they promote a decrease of the drag through a 2 5 and a 5 a 7 . As a non-trivial result, the statistically asymmetric (stable) limit cycles have a larger drag than the statistically symmetric (unstable) limit cycle, while the asymmetric steady solutions have a lower drag than the symmetric steady solution. This is obvious in figure 11 when considering the relative positions of the three steady solutions and three limit cycles along the C D axis. Note that the parameters q x;11 , q x;22 and q x;57 all own negative signs but are relatively small. The two parameters q x;11 , q x;22 solely contribute to the oscillating dynamics, as discussed in § 4.1, while q x;57 optimizes the fitting result when evolving toward the attracting limit cycles.
Analogously, for the lift model, the values of l y;4 and q y;45 are taken from the identified lift model in (5.2b), while q y;17 , q y;27 are set to zero for consistency with the structure of (5.2b), in which q y;15 , q y;25 are absent. Based on the structure of model (5.2b), the terms a 6 , a 6 a 7 , a 3 a 6 , a 5 a 6 and a 4 a 7 are introduced in the extended model. The final sparse form is identified by the LASSO regression on gradually increasing the L1-penalty. A sparse lift model, compatible with the structure of (5.2b), is derived as C L = 0.00762433a 1 + 0.01102097a 2 − 0.10179203a 4 − 0.03129798a 6 − 0.00141416a 1 a 3 − 0.00289952a 2 a 3 + 0.00656293a 3 a 4 − 0.01082375a 3 a 6 + 0.05914386a 4 a 5 + 0.02365784a 4 a 7 − 0.03348935a 5 a 6 . (5.8) In addition to the lift model of (5.2b), the lift model of (5.8) contains the terms a 6 , a 3 a 6 and a 5 a 6 , as well as the coupling between a 4 to a 7 . The r 2 score has increased to 0.9952. Both the oscillating dynamics in the early stage and the symmetry-breaking stage are better reproduced for the lift coefficient, as the red dashed curve of figure 15(c,d) proves.
With the two additional degrees of freedom a 6 , a 7 , the time evolution of the drag and lift coefficients are well reproduced, as shown in figure 15. Without notable changes of the original lift structure, the phase of the lift dynamics is now correctly caught along with the complete transient dynamics.

Force model at Re = 100
In § 4.2, we derived a basic force formula for the primary stage of the transient evolution at Re = 100, when only the degrees of freedom of the pitchfork bifurcation were involved. We now consider the complete force evolution at Re = 100. Figure 16 shows trajectories issued from the three different steady solutions in the three-dimensional time-delayed embedding space of C L and C D . The black trajectory, issued from the symmetric steady solution u s (black cross × in figure 16) first approaches the asymmetric steady solution u + s (red point) before escaping out of it and eventually reaching the stable (statistically asymmetric) limit cycles aroundū + .
The same mode decomposition strategy is proposed, resulting in a reduced-order model with 7 modes. The mode amplitudes from two DNS, starting from either the symmetric steady solution u s (a) or the asymmetric steady solution u + s (b), are shown in figure 17. As already observed in figure 10, the drag coefficient (solid black line) in figure 18(a) exhibits a minimal value for a transient state around t ≈ 700. This transient state is the asymmetric steady solution u + s (red circle of figure 16). In the frame of our modal decomposition (3.1), u + s is approximated as u + s ≈ u s + a 4 (700)u 4 + a 5 (700)u 5 , (5.9) with only a 4 and a 5 being active in the dynamics of the fluidic pinball, as can be seen in figure 17(a). From (4.8a), the drag coefficient only depends on a 5 and a 2 5 , which actually contribute to the transitory increase and an overall decrease on the drag. This is fully consistent with the transition of the drag coefficient observed in figures 10(a) and 18(a) from t = 300 to 700; a 5 is found to contribute to the initial rising of C D , around t ≈ 650, while a 2 5 contributes to the subsequent decrease of the drag coefficient, around t = 700. 918 A4-23 The degrees of freedom associated with the Hopf bifurcation become active later during the transient dynamics, when the state space orbit leaves the unstable asymmetric steady solution u + s toward the stable attracting limit cycle aroundū + s . The training data are the real force coefficients and the mode amplitudes taken from the DNS starting with the three different steady solutions to the final asymptotic regimes. Following the same calibration procedure as for Re = 80, we first apply the LASSO regression for the force model with the five leading degrees of freedom, and then introduce the two additional degrees of freedom a 6 , a 7 into the regression for optimization. Performing the sparse regression in this way can prevent the elimination of a 4 , a 5 and ensure the corrective effect of a 6 , a 7 , thereby improving the robustness of the identification. The force model at Re = 100 reads C D = 3.58248992 + 0.04367604a 5 − 0.00302817a 2 1 − 0.00354079a 2 2 + 0.00158873a 2 3 + 0.02169661a 3 a 5 + 0.02223079a 3 a 7 − 0.08525184a 2 5 − 0.04763643a 5 a 7 , (5.10a) C L = 0.00346208a 1 + 0.00269236a 2 − 0.13611053a 4 + 0.05962648a 6 + 0.00029274a 1 a 3 − 0.00045784a 2 a 3 + 0.00389912a 3 a 4 − 0.02102284a 3 a 6 + 0.09194312a 4 a 5 + 0.02056288a 4 a 7 − 0.10980990a 5 a 6 , (5.10b) with r 2 = 0.9984 for the drag model of (5.10a), and r 2 = 0.9901 for the lift model of (5.10b). As shown in figure 18, the force model fits well the time evolution of the drag and lift coefficients. Moreover, (5.7), (5.8) and (5.10) own the same active terms. Henceforth, the drag force model preserves the same structure with the same signs of the active terms as the Reynolds number is increased. In addition, although the transient dynamics at Re = 80 and 100 is qualitatively very different, with the seven degrees of freedom differently activated during the transient, the force model of (5.7) and (5.8) is still consistent at Re = 100, with the correctly identified mean-field model. For the lift model (5.10b), we notice the same structure with the sign changes for the terms a 1 a 3 and a 6 , compared to (5.8), which is acceptable for the oscillating dynamics. Compatible with the basic lift force model, the lift force model with seven degrees of freedom also correctly identifies the force transitions, as shown in figure 18(c,d).

Conclusions and outlook
We proposed aerodynamic force formulae complementing mean-field POD Galerkin models for the unforced fluidic pinball. The starting point is a general Galerkin method for unsteady incompressible viscous flow around a stationary body. First, the instantaneous force is derived as a constant-linear-quadratic function of the mode amplitudes from first principles. The viscous and pressure contributions to the force are directly obtained from the Galerkin expansion and lead to a constant-linear-quadratic force in terms of the mode amplitudes. These terms lead to corresponding changes in the flow from which the force can also be derived. One contribution from the convective term describes the momentum flux contribution. The additional contribution from the local acceleration requires the Galerkin system to replace the time derivatives of the mode amplitudes by a state function. In contrast to the pioneering work by Noca et al. (1999), the derivation is valid for arbitrary multiply connected domains.
The drag and lift formula is simplified for the fluidic pinball model exploiting the symmetry of the modes. Approximately half of the terms can be discarded on the grounds of symmetry. A second simplification is performed with a sparse calibration of the remaining coefficients. The sparsity parameter λ penalizes any non-vanishing term and yields sparse human-interpretable expressions. The challenges of the purely projection-based approach is discussed in appendix C, and the challenges of using standard POD modes is elaborated in appendix D.
The sparse force model methodology is applied to three transient dynamics: (i) the periodic regime of statistically symmetric vortex shedding at Re = 30, (ii) the periodic regime of statistically asymmetric vortex shedding at Re = 100 and (iii) the same regime at Re = 80 but with metastable statistically symmetric periodicity.
The transient dynamics at Re = 30 from the steady solution to the limit cycle is resolved by standard third-order mean-field Galerkin model with two oscillatory modes for vortex shedding and one shift mode for the mean-field distortion . The drag formula includes the squares of all mode amplitudes consistent with the second harmonic fluctuations. The drag monotonically increases during the transient. The lift formula includes the amplitudes of the von Kármán modes and their products with the shift mode, consistent with expectations. Its oscillation increases until the limit cycle is reached.
The dynamics at Re = 100 after the Hopf and pitchfork bifurcation has three unstable fixed points, one symmetric steady solution and a mirror-symmetric pair of asymmetric ones. The transients from these fixed points terminate in one of the asymmetric limit cycles corresponding to the asymmetric shedding states. This dynamics is described by a fifth-order Galerkin model (Deng et al. 2020), where the first three modes resolve the Hopf bifurcation and the next two modes the pitchfork bifurcation. The associated drag formula contains the terms of Re = 30. In addition, the drag is modified by linear and quadratic terms with the shift modes associated with the Hopf and pitchfork instabilities. These additional terms vanish without pitchfork bifurcation and do not introduce harmonics of vortex shedding. Similarly, the lift formula generalizes the expression at Re = 30.
The intermediate Reynolds number 80 leads to a more complex force model, as the transients may pass through a meta-stable symmetric limit cycle. The accuracy of the force model could significantly be increased by two additional Galerkin modes, which resolve variations between symmetric and asymmetric limit cycles. The drag and lift formulae were correspondingly longer and good agreement with computational data is achieved.
Summarizing, the sparse force model describes multi-attractor behaviour of the unforced fluidic pinball even for a complex dynamics with three steady and three periodic solutions. For this configuration, we have the advantage of a thorough understanding of the dynamics via a low-dimensional mean-field Galerkin model. We envision successful applications of sparse regression for aerodynamic forces for turbulent flows, e.g. for the bi-stable behaviour of the Ahmed body wake (Grandemange, Gohlke & Cadot 2013;Östh et al. 2014;Barros et al. 2017).
The force formula may be particularly instructive for drag reduction with active control (Choi, Jeon & Kim 2008). Given a Galerkin model, the force formula indicates beneficial regions of the state space. Thus, an upfront kinematical insight is gained in which direction control needs to 'push' the attractor. For instance, the third-order mean-field model and the force formula imply that stabilization is required for drag reduction, consistent with the earlier studies of Protas (2004) and Bergmann & Cordier (2008). Future generalizations may also profit from stochasticity (Sapsis & Lermusiaux 2009).

Appendix A. Forces from the momentum balance
The forces can be alternatively derived from the residual of the Navier-Stokes equations in the domain Ω. This domain is assumed to enclose the obstacle and extend sufficiently far away from the obstacle such that the free-stream condition u = e x can be applied on the left, top and bottom boundaries of the fluid domain Ω. The domain boundary ∂Ω contains the surface of the immersed body Γ and the outer surface S ∞ . It should be noted that the surface element dS on the body points inside the body, i.e. opposite to the direction in § 2.2. The force in direction e α is derived from the integrated momentum balance in that direction (e α , R(u, p)) Ω = 0. (A2) Four terms are obtained. The first contribution is the viscous term. This term can be converted into a skin friction integral over Γ and S ∞ . The contribution over the outer integral vanishes under free-stream conditions. The remaining contribution is the viscous force applied to the immersed body where e α · (∇u + (∇u) T ) · n = 2 α,β=x,y,z S α,β n β . The second contribution is the pressure term which can analogously reduce to the pressure force on the immersed body Not surprisingly, we arrive at the formula of § 2.2. The force exerted on the body is equal but opposite to the force exerted on the fluid. The third term is the local acceleration ⎛ where m t α;j = (e α , u j ) Ω . The fourth term arises from the convective acceleration ⎛ where q c α;jk = (e α , ∇ · [u j ⊗ u k ]) Ω . The volume integral over Ω can be converted into a momentum flux surface integral over the boundary. Making use of the momentum balance (A2), the third and fourth contributions from the acceleration terms equal the total force This force formula contains constant, linear and quadratic terms of the mode amplitudes as well as their time derivatives. The state-dependent formula (2.18) may be obtained from (A7) by replacing the time derivatives with (2.6). The total forces on the immersed body are here again represented by a constant-linear-quadratic expression. The above mentioned formulae dresses Newton's second law F = ma in a Galerkin framework for fluid flow. Equation (A7) corresponds to 'ma' and is purely based on the fluid motion. Equation (2.18) corresponds to 'F ' and allows distinguishing between the contribution of viscous and pressure stresses.

Appendix B. Influence of the sparsity parameter and regression methods
In the SINDy algorithms, the sparsity parameter is either the L1-penalty for the LASSO regression or the threshold for the sequential thresholded least-squares (STLS) regression. We denote the L1-penalty and the threshold as the sparsity parameter λ in both cases. These two methods can, however, lead to different results. We can choose the one with a better performance according to the actual needs.
In § 4.1, we derived the sparse drag model with three degrees of freedom at Re = 30. Benefiting from the low cost of computation for the regression test, we can iteratively run the algorithm changing the sparsity parameter λ and investigate the performance changes of the identified model. The performance of the identified drag model by these two regression methods when varying λ is illustrated in figures 19(a) and 20(a), together with a comparison with the real force dynamics for three typical values of λ.
The sparsity parameter starts with 0 (pure least-square regression) and increases up to nearly 1. The structures of the resulting models at λ = 0.95 for the LASSO regression, see figure 19(d), and λ = 0.45 for the STLS regression, see figure 20(c), are identical, where only a 2 3 remains. However, the STLS regression is more sensitive to the sparsity parameter, as shown in figure 20(a). The terms a 3 , a 2 1 and a 2 2 are eliminated at the same time. The remaining a 2 3 is replaced by a 3 with λ > 0.48, and the identified models are obviously under-fitted, as shown in figure 20(d). In contrast, the LASSO regression eliminates the terms gradually, first a 1 a 2 , then a 3 , and eventually a 2 1 together with a 2 2 . In figure 19(a), the elimination of a 2 1 and a 2 2 only reduces the r 2 score by 0.003. But the loss of the fluctuating drag dynamics indicates an under-fitting. Hence, the optimal λ is found for 0.85.
To figure out the reason for the failure of the identification with the STLS regression when λ > 0.48, we compare the evolution of the coefficients with increasing λ in figure 21; a 1 a 2 is the first eliminated term in both cases. The coefficients in the initial stage before the elimination of a 3 are almost the same. After the elimination of a 3 with the LASSO regression, as shown in figure 21(a), the coefficients of a 2 1 and a 2 2 become of order O(10 −3 ). Since the STLS regression algorithm thresholds the terms with smaller coefficients, the tiny coefficients of a 2 1 and a 2 2 will be set to zero simultaneously. When the STLS regression is used with a too large sparsity parameter λ, the term with larger coefficient can survive. As illustrated in figure 21(b), the remaining term a 2 3 is replaced by a 3 with a larger coefficient. This explains the reason why the STLS regression final converges to a 3 , which is obviously the wrong term for the real drag force dynamics in figure 20(d).
We apply the same analysis for the sparse drag model with five degrees of freedom at Re = 80, as described in § 5.1. The evolutions of the performances under the two regression methods are shown in figure 22. The STLS regression goes in the wrong direction for λ > 0.17. After checking the list of coefficients, the key term a 2 3 is deleted irretrievably, resulting in the inability of the model to fit correctly. However, the regression result right before the critical value provides the most simplified and relevant drag model of (5.2a) with r 2 = 0.9816. The LASSO regression is much safer in the elimination of terms; a 2 3 can survive during the regression in all of the range of λ from 0 to almost 1. This further indicates that the key terms have a better robustness in the LASSO regression. From figure 22(a), the optimal λ is chosen at 0.85, involving six terms and r 2 = 0.9791. Although there are only five terms remaining when λ = 0.9, the resulting model is not stable. It returns to six terms and r 2 = 0.9755 at λ = 0.95, with different active terms compared to the model at λ = 0.85. At the optimal value, the identified drag model consists of terms a 5 , a 2 1 , a 1 a 2 , a 2 2 , a 2 3 , a 3 a 5 , where a 2 5 is missing. Since a 1 a 2 is of order O(10 −4 ), we can directly apply the least-square regression to the updated library by deleting a 1 a 2 and adding a 2 5 . The regression result is the same as for the STLS regression. condition for the pressure, i.e. the normal derivative of p in the outward direction n must vanish on the whole domain boundary ∂Ω, ∂ n p = n · ∇p = 0. (C1) In this study, we apply a no-slip condition on the velocity without the above-mentioned Neumann boundary condition on pressure. Hence, the partial pressure fields p jk cannot be determined to a constant pressure field. Analogously, q p α;jk cannot be solved with an exact value. Even if we assume Neumann boundary conditions for the pressure field p, it is still a numerically challenging work since the pressure field are expanded to numerous partial pressure fields p jk , see (2.15).
Without considering the pressure force contribution, we can reconstruct the viscous force from the viscous force contribution of the bifurcation modes. The resulting viscous force model only contains linear terms and reads C ν D = 1.01814664 + 0.00159948a 3 − 0.0023798a 5 + 0.00601715a 7 , C ν L = 0.000267167a 1 + 0.00004522a 2 − 0.01409768a 4 − 0.0055717a 6 .
The viscous force contributions of each bifurcation mode is explicitly computed without any symmetry assumption, no sparsity can be expected in this model. Yet, after eliminating terms with a coefficient less than O(10 −5 ), the resulting force models (C2) only involve the terms associated with the bifurcations modes with the appropriate symmetry, indicated as the symmetric modes u 3 , u 5 , u 7 in C ν D and the symmetric modes u 1 , u 2 , u 4 , u 6 in C ν L . The performance of the force model using the real viscous force contribution of the seven bifurcation modes is illustrated in figure 23. The r 2 score for the viscous drag model is 0.9786 and 0.9183 for the viscous lift model. The accuracy and the predictive ability of the force model are acceptable for the drag model with only three items and the lift model with only four terms.
Due to the lack of boundary conditions for the pressure field contribution, we only focus on the reconstruction of the viscous force with the purely projection-based approach. The contribution to the viscous drag and lift forces, given by max |q ν α;j a j |, with α = x, y, are shown in figure 24. The main force contribution comes from the leading 50 POD modes. The viscous force reconstructed with the N leading POD mode amplitudes reads The viscous drag C ν D and lift C ν D coefficients reconstructed with different numbers of POD modes are compared to the real force dynamics in figure 25.
For a sequential N, the error of the reconstructed force coefficients with N leading POD modes can be also evaluated with the r 2 score. A higher r 2 score indicates less error in the reconstructed force. As expected, the error tends to decrease when the number of POD modes is increased in figure 26. To achieve r 2 > 0.999, N = 36 leading POD modes are required for the drag force, and N = 51 modes for r 2 > 0.9999. For the lift force, these two critical numbers are respectively N = 30 and N = 86. In actual situations, the model with r 2 > 0.999 has enough accuracy. Note that no sparsity is involved in the model because the force contribution of each POD mode is computed explicitly.
We now focus on the regression-based approach, we set a truncation of the model with N = 10, 20, 50 leading POD modes, and try to use the sparse regression to find a drag model with a balance between accuracy and complexity. To be noted, the drag force considered here involves both the pressure and viscous contributions to the force. To reach the same r 2 score, it requires more POD modes due to the additional quadratic complexity of the pressure force contribution.
The identified system coefficients are recorded in table 1, and the model performance is exemplified in figure 28.