Geometrical dependence of optimal bounds in Taylor–Couette ﬂow

This paper is concerned with the optimal upper bound on mean quantities (torque, dissipation and the Nusselt number) obtained in the framework of the background method for the Taylor–Couette ﬂow with a stationary outer cylinder. Along the way, we perform the energy stability analysis of the laminar ﬂow, and demonstrate that below radius ratio 0 . 0556, the marginally stable perturbations are not the axisymmetric Taylor vortices but rather a fully three-dimensional ﬂow. The main result of the paper is an analytical expression of the optimal bound as a function of the radius ratio. To obtain this bound, we begin by deriving a suboptimal analytical bound using analysis techniques. We use a deﬁnition of the background ﬂow with two boundary layers, whose relative thicknesses are optimized to obtain the bound. In the limit of high Reynolds number, the dependence of this suboptimal bound on the radius ratio (the geometrical scaling) turns out to be the same as that of numerically computed optimal bounds in three different cases: (1) the perturbed ﬂow only satisﬁes the homogeneous boundary conditions but need not be incompressible; (2) the perturbed ﬂow is three-dimensional and incompressible; (3) the perturbed ﬂow is two-dimensional and incompressible. We compare the geometrical scaling with the observations from the turbulent Taylor–Couette ﬂow, and ﬁnd that the analytical result indeed agrees well with the available direct numerical simulations data. In this paper, we also dismiss the applicability of the background method to certain ﬂow problems and therefore establish the limitation of this method.

A. Kumar transport, and mixing efficiency) on input parameters. The lack of analytical solutions of the Navier-Stokes equations in the fully turbulent regime has forced the scientific community to adopt a multi-faceted approach to this problem, in which simple physical theories and reduced models are proposed, and then corroborated by direct numerical simulations (DNS) and/or results from laboratory experiments. However, the inability to perform simulations and experiments in the extreme parameter regimes that often concern atmospheric, oceanic and astrophysical flows and engineering applications leaves these theories unsubstantiated.
In these extreme parameter regimes, an alternative approach that can provide meaningful information is to obtain rigorous bounds on the aforementioned global properties. The first method to obtain bounds was developed by Howard (1963) and Busse (1969), but it was not until the 1990s that bounding techniques gained general popularity, with the introduction of the so-called 'background method' by Doering and Constantin (Doering & Constantin 1992, 1996Constantin & Doering 1995). The background method is based on ideas from Hopf to produce a priori estimates for the solutions of the Navier-Stokes equations with inhomogeneous boundary conditions (Hopf 1955). So far, it has been applied to many different fluid mechanics problems (Doering & Constantin 1992, 1996Constantin & Doering 1995;Caulfield & Kerswell 2001;Tang, Caulfield & Young 2004;Whitehead & Doering 2011;Goluskin & Doering 2016;Fantuzzi 2018;Fantuzzi, Pershin & Wynn 2018;Kumar & Garaud 2020;Arslan et al. 2021a,b;Fan, Jolly & Pakzad 2021;Kumar et al. 2022). See Fantuzzi, Arslan & Wynn (2022) for a recent review.
In the background method, we write the total flow field as a sum of two flow fields: the background flow and the perturbed flow. To obtain a bound on the desired bulk quantity requires choosing a background field that satisfies a certain integral constraint (extracted from the governing equations of the perturbed flow). Generally, one takes one of the following two routes. The first route is to specify a functional form of the background flow and then use standard inequalities. This route leads to an analytical but suboptimal bound on the bulk quantity as a function of system parameters. The second route is to find the best possible bound (optimal bound) through a variational formulation of the background method in which one solves the corresponding Euler-Lagrange equations, usually numerically. Numerous studies pertaining to the background flow have concentrated on the scaling of optimal bounds as a function of the principal flow parameter, such as the Reynolds number and the Rayleigh number. However, only a handful of them studied the variation of these bounds with the shape of the domain. One such study is by Wen et al. (2013), where the authors were interested in determining the dependence on aspect ratio of the optimal bound on heat transfer in porous medium convection.
In this paper, we are concerned with the question of whether it is possible to obtain the analytical expression for the dependence of optimal bounds on the geometrical parameters of the system. Indeed, while the numerically obtained optimal bounds usually follow an easily identifiable simple power law in the principal flow parameter, the variation of the optimal bounds with geometrical parameters is not so readily apparent. Furthermore, we also aim to determine whether this analytical form bears any resemblance to the actual dependence of the corresponding bulk quantity on system geometry in fully turbulent flows. This question is motivated by engineering applications where the geometry plays an important role.
In a recent study, we attempted to provide bounds on the friction factor in the context of pressure-driven helical pipe flows (Kumar 2020). We focused in particular on the dependence of this bound on the geometrical parameters: the curvature and torsion of 948 A11-2

A. Kumar
Reynolds number tends to infinity, but also that this dependence is the same as the one obtained from the suboptimal analytical bound.
The arrangement of the paper is as follows. We begin by describing the problem configuration, the definitions of the relevant mean quantities, and the relations between those quantities in § 2. In § 3, we perform the energy stability analysis of the laminar flow. In § 4, we obtain analytical bounds on the mean quantities. Section 5 presents optimal bounds obtained in the three cases listed above and compares the results with the analytical bounds from § 4. In § 6, we show that the background method cannot be applied to certain flow problems past certain Reynolds numbers. Finally, § 7 presents a discussion, comparison with DNS results, the broad applicability of the present study and open problems.

Problem set-up
Consider the flow of an incompressible Newtonian fluid of density ρ and kinematic viscosity ν between two coaxial circular cylinders, where the inner cylinder rotates with constant angular velocity Ω and the outer cylinder is stationary. The radius of the inner cylinder is R i , and the radius of the outer cylinder is R o . The quantity η = R i /R o is referred to as the radius ratio hereafter, and d = R o − R i is the gap width. We non-dimensionalize the variables as follows: where p 0 is the reference pressure, and x, u, t and p denote the non-dimensional position, velocity, time and pressure, respectively. The starred variables are the corresponding dimensional quantities. In non-dimensional form, the governing equations are ∇ · u = 0, (2.2) ∂u ∂t + u · ∇u = −∇p + 1 Re is the Reynolds number, which, along with the radius ratio η, characterizes the flow field fully. Note that instead of the Reynolds number, one can also use the Taylor number (1 + η) 6 64η 4 Re 2 (2.5) to characterize the flow field. We use a cylindrical coordinate system (r, θ, z). The boundary conditions are (u r , u θ , u z ) = (0, 1, 0) at r = r i , (2.6) (u r , u θ , u z ) = (0, 0, 0) at r = r o , (2.7) where r i and r o are the non-dimensional inner and outer cylinder radii. In this paper, we will assume that the flow is periodic in the z direction with non-dimensional length L. The 948 A11-4 domain of interest, denoted by V, is therefore given by V = {(r, θ, z) | r i r r o , 0 θ < 2π, 0 z < L}. (2.8) At sufficiently small Reynolds numbers, or equivalently, at small Taylor numbers, the flow is laminar and can be expressed as (2.9) Before proceeding further, it is useful to introduce a few convenient notations. We use angle brackets for the volume integration, and overbar for the long-time average of a quantity: The L 2 -norm of a quantity is henceforth denoted as In what follows, the three quantities that we are interested in bounding are the energy dissipation rate, the torque and the equivalent of a Nusselt number (defined based on the transverse current of azimuthal velocity). These quantities are not independent, as we now demonstrate. We start by writing the dimensional expression of the time-averaged torque required to rotate the inner cylinder: where τ * rθ denotes the shear-stress. In non-dimensional form, the torque is given by In a statistically stationary state, the work done by the torque to rotate the inner cylinder eventually dissipates in the fluid, i.e.

A. Kumar
of the evolution equation of the total kinetic energy. The dissipation per unit mass non-dimensionalized by Ω 3 R 3 i /d is given by From the divergence-free condition (2.2), and the boundary conditions (2.6) and (2.7), along with the use of the divergence theorem, one finds that As a result, the non-dimensional dissipation can also be written as Using (2.13), (2.14) and (2.17), we finally obtain a relation between the non-dimensional torque and the non-dimensional dissipation as which is the non-dimensional version of (2.14). Another quantity of interest is the transverse current of azimuthal velocity as defined in Eckhardt, Grossmann & Lohse (2007), given by where ω * = u * θ /r * is the local angular velocity. As shown by Eckhardt et al. (2007), J ω * is independent of the radial direction. In an analogy with Rayleigh-Bénard convection, one defines the Nusselt number as the ratio of the transverse current of azimuthal velocity to its corresponding value in the laminar regime, i.e. (2.22) Substituting r * = R i in the right-hand-side of (2.21), one obtains the following relation between the torque and the transverse current of azimuthal velocity: implying that the Nusselt number can also be written as where G lam and ε lam are the values of the non-dimensional torque and dissipation in the laminar regime, respectively.

Energy stability analysis
We begin by discussing the energy stability of the laminar flow u lam . The importance of energy stability analysis in the context of bounding theories comes from the fact that bounds on mean quantities introduced in the previous section are by definition saturated by the laminar state below the energy stability threshold. The energy stability of the laminar Taylor-Couette flow has been studied before both theoretically and numerically, by e.g. Serrin (1959) and Joseph (1976). In these studies, the general conclusion was that at the energy stability threshold, the least stable perturbations are axisymmetric Taylor vortices. However, as we will demonstrate in this section, this commonly accepted result does not hold below a certain radius ratio (η < 0.0556). Instead, we find that the least stable perturbations at the energy stability threshold in that case are fully three-dimensional. We begin by defining the functional whereṽ is a perturbation over the laminar flow that satisfies the homogeneous boundary conditions at the inner and outer cylinders. From the governing equations, one can show that the laminar flow u lam is energy stable when H(ṽ) is non-negative (see e.g. Serrin 1959;Ding & Marensi 2019). We will consider three types of constraints on the perturbationsṽ: no constraints, other than the homogeneous boundary conditions (case 1); three-dimensional (3-D) incompressible perturbations (case 2); and two-dimensional (2-D) (z-invariant) incompressible perturbations (case 3). We perform an energy stability analysis for each of these cases, and present the results as a function of the radius ratio. We note that recently, Ding & Marensi (2019) also studied the energy stability of the laminar state in Taylor-Couette flow but only for the axisymmetric perturbations. The critical Taylor number Ta c defining the energy stability threshold is the largest Taylor number for which the functional H(ṽ) is non-negative. For clarity, we add superscripts and use the notation Ta nc c , Ta 3D c and Ta 2D c when referring to case 1, case 2 and case 3, respectively. The statement of the non-negativity of the functional H(ṽ) can be posed as a convex optimization problem, where we require the minimum value of H to be non-negative. Then it can be shown using the corresponding Euler-Lagrange equations that the non-negativity of the functional H(ṽ) is equivalent to the non-negativity of the smallest eigenvalue in the eigenvalue problem A. Kumar satisfies the homogeneous boundary condition at r = r i is given bỹ The critical Reynolds number for energy stability is the smallest value of Re for which this solution also satisfies the homogeneous boundary condition at r = r o . We then obtain the critical Taylor number using (2.5), which leads to (3.4) In case 2 (3-D incompressibleṽ) and case 3 (2-D (z-invariant) incompressibleṽ), we must turn to numerical computations to calculate the critical Taylor number. To find the eigenvalues of (3.2), we first transform the equations into a generalized eigenvalue problem using the spatial discretization described in § 5, and then use the DGGEV routine by LAPACK for the computation. Let us call the critical wavenumbers of the least stable perturbation at the energy stability threshold 2π/L c (where L c would then be known as the critical aspect ratio) in the z-direction, and m c in the θ-direction. We use the bisection algorithm in the Taylor number, and the ternary search algorithm in aspect ratio or azimuthal wavenumber (depending on the case at hand), to determine accurately Ta c , L c and m c . The dependence of the critical Taylor number for energy stability on the radius ratio η is shown in figure 1 for all three cases. The critical axial wavenumber (2π/L c ) and the critical azimuthal wavenumber (m c ) of the corresponding perturbations in case 2 and case 3 are shown in figure 2. From figure 1(a), we see that the critical Taylor number increases as we go from case 1 (green line) to case 3 (red line), which is not surprising since we increase correspondingly the number of constraints on the perturbations. In all three cases, the critical Taylor number increases monotonically with decreasing η, and tends to infinity as η → 0. By contrast, the critical Taylor number tends to a constant in the small gap width limit (η → 1): in case 1, Ta nc c → 4π 4 ≈ 389.6364, whereas in case 2 and case 3, Ta 3D c → 6831 and Ta 2D c → 31 641, which are, respectively, 17.5 and 81.2 times larger than in case 1. In this limit, the marginally stable perturbation in case 2 recovers the well-known axisymmetric Taylor vortices (Serrin 1959;Joseph 1976). In case 3, the marginally stable perturbation is composed of vortices whose axis is parallel to the cylinder axis (Harrison 1921). Figure 1(b) shows a zoomed-in version of figure 1(a) for small values of η. We also show, for case 2 (blue line), a separate curve that assumes that perturbations are axially symmetric (dashed blue line). For large radius ratio, the two are identical, confirming that the axisymmetric Taylor vortices are indeed the least stable perturbations. However, we note that below radius ratio η s = 0.0556, the marginally stable perturbation switches from the axisymmetric Taylor vortices to being fully 3-D. Figure 3 shows the marginally stable 3-D flow and Taylor vortices at η = η s . A distinctive feature of the marginally stable 3-D flow, compared to marginally stable axisymmetric Taylor vortices, is that one end of a typical vortex lies near the outer cylinder, but the other end lies at one of the two lines that are offset from the inner cylinder. Also, the critical aspect ratio corresponding to marginally stable 3-D flow is larger than the one corresponding to the Taylor vortices. In fact, with further decrease in the radius ratio, the axisymmetric critical aspect ratio corresponding to the marginally stable 3-D flow grows, whereas the one corresponding to Taylor vortices shrinks, as can 948 A11-8 been seen from figure 2(a). The decrease of the aspect ratio of the critical perturbations implies that the term ∇ṽ 2 2 increases rapidly as η → 0, which causes the corresponding critical Taylor number for axisymmetric flows to do the same. This explains why the axisymmetric perturbations are no longer preferred for very low η. At η = 0.0188, the In both cases, the radius ratio is η s = 0.0556, and the corresponding critical Taylor numbers are equal. The streamlines are coloured according to the magnitude of the velocity. In both the cases, the velocity field has been scaled such that the maximum magnitude is 1. A typical vortex is shown using relatively thicker lines in both cases. Note that only half the vortex is shown in the axial direction. critical Taylor number for the marginally stable Taylor vortices becomes even larger than the one corresponding to the two-dimensional flow (Ta 2D c ). Given that we were able to compute the critical Taylor number in case 1 analytically as a function of η, it is worth investigating whether the dependence of Ta c on η in cases 2 and 3 is similar to that of case 1. To do so, we look at figures 1(c) and 1(d), which show the ratios Ta  (3.5) However, the same is not valid for case 3, where Ta 2D c /Ta nc c varies substantially with η. The spikes in figure 1(d), which are not visible in figure 1(a), correspond to the discrete change in critical azimuthal wavenumber when η varies, shown in figure 2(b).

A. Kumar
For small radius ratio, it is possible to predict the asymptotic behaviour of Ta 3D c and Ta 2D c . We find that both Ta 3D c /Ta nc c and Ta 2D c /Ta nc c decrease as η → 0, as can be seen in figures 1(c) and 1(d). By construction, the asymptotic values of the ratios have to be larger than 1. Therefore, in the small radius ratio limit, we can obtain the asymptotic behaviour of Ta 3D c and Ta 2D c as

An analytical bound
In this section, we obtain a simple, suboptimal, analytical bound on the torque, the rate of energy dissipation and the Nusselt number defined in § 2. We use the well-known background method (Doering & Constantin 1992) whose exact formulation in the context of the present problem is given in Appendix A. As usual, we define U to be the background flow and v to be the perturbed flow, such that the total flow is u = U + v.

Optimal bounds in Taylor-Couette flow
The background flow U is divergence-free and satisfies the same boundary conditions as u, so the perturbed flow v satisfies the homogeneous version of the boundary conditions. For mathematical convenience (see Appendix A), we further define the so-called 'shifted perturbation'ṽ = v − φ (see (A17)), and we simply refer toṽ as the perturbation from here onward. As shown in Appendix A, a bound on the rate of energy dissipation, can be obtained for any choice of the background flow for which the functional (see (A22)) is positive semi-definite. In (4.1), the constant a is a balance parameter that takes values between 0 and 2. This bound is identical to the one obtained by Ding & Marensi (2019), after noting that they used a different non-dimensionalization. While showing that H(ṽ) is non-negative, we do not impose the incompressibility constraint on the perturbationsṽ and assume only thatṽ satisfies the homogeneous boundary conditions. We make a choice of the background flow U for which is non-zero only in boundary layers, which are assumed to have thicknesses δ i and δ o near the inner and outer cylinders, respectively. In particular, the selected background flow U is then where Λ is an O(1) constant, i.e. independent of Re. The decision to allow for different boundary layer thicknesses is inspired from the work of Kumar (2020), who speculated in the context of helical pipe flows that by doing so, it is possible to capture important geometrical aspects of problems that would otherwise not appear. As we are interested primarily in deriving bounds at asymptotically high Reynolds numbers, for convenience, we define rescaled boundary layer thicknesses as where, by construction, h i , h o are greater than 0 and are O(1). Our goal in this section is to adjust the relative size of the boundary layers (h i /h o ) to optimize the bound (4.1) simultaneously for different values of η in the limit of high Reynolds number.

A. Kumar
We start by obtaining a simple estimate for the quantity In deriving the result, we have used the fundamental theorem of calculus in the first line and Hölder's inequality in the second line, followed by an integration in r to obtain the third line. Finally, we used Young's inequality to obtain the last line. In a similar manner, one can also show that Next, we note that Using estimates (4.6)-(4.9) along with the expression of the background flow (4.4), we finally obtain a simple bound on term II in (4.2) as This shows that the functional H is positive semi-definite as long as (4.12) Using (2.9) and (4.4) in (4.1), we then obtain an upper bound on the dissipation as follows: The upper bound obtained is called ε a b ; we use 'b' in the subscript to signify that it is a bound, and use 'a' in the superscript to signify that it is obtained analytically. In the final 948 A11-12 step of the procedure, we adjust the values of the unknown parameters h i , h o , Λ and a to optimize the bound (4.13) while satisfying the constraint (4.12). The optimal values of the parameters, in the limit of high Reynolds number, are Using (4.5a,b) and (2.5), we can now write the optimal choice of boundary layer thicknesses δ i and δ o in the limit of Re → ∞ (or equivalently Ta → ∞) as The corresponding bound on the dissipation in the limit of Re → ∞ is then given by (4.16) Here, we added '∞' in the subscript to indicate that it is the main term of the bound in the limit Re → ∞. Using the relationship (2.24), we obtain an equivalent upper bound on the Nusselt number in the high-Reynolds-number limit as This expression contains a dependence on both the Taylor number (the principal flow parameter) and the radius ratio (the geometrical parameter). To separate out the geometrical dependence in (4.17), we define and call it the geometrical scaling of the bound on Nu. This geometrical scaling is defined in such a way that χ(1) = 1 (the relevance of η = 1 being that it corresponds to the plane Couette flow case). Finally, by combining (4.16) with the relation (2.20), we obtain an upper bound on the torque as a function of the Reynolds number:  had previously obtained a bound on the torque in Taylor-Couette flows by considering a background flow with a single boundary layer. The bound obtained by Constantin is also proportional to Re 2 , as in (4.19). But the coefficient in front has a different dependence on the radius ratio η. The reason for this difference is that we chose a background flow with two boundary layers and adjusted their relative thicknesses to optimize the bound. We will see later that this optimization procedure enables us to capture the actual dependence of optimal bounds on the radius ratio.

Optimal bounds
In this section, we now proceed to obtain optimal bounds on the bulk quantities, i.e. the best possible bounds within the framework of the background method. As described

A. Kumar
in § 1, we consider three scenarios, case 1, case 2 and case 3, in which we impose constraints incrementally on the perturbed flow field and obtain numerically the optimal bounds in each case, which allow us to examine systematically the hypothesis stated in the Introduction.
The general development of the background method for Taylor-Couette flow is presented in Appendix A. In what follows, we first describe our numerical algorithm, then proceed to present the results.

Numerical algorithm
Here, we first describe the general numerical framework used to compute the optimal bounds, and then provide further details of the algorithm in each of the specific cases. Finding the optimal bound begins with the same background method applied to the Taylor-Couette flow as in § 4, which is described in Appendix A. However, instead of using functional inequalities, we now follow the standard route toward optimal bounds, and derive a set of Euler-Lagrange equations that optimal solutions satisfy, given specific constraints in each case. The derivation is presented in Appendix A, and the equations are given in (A28a-d). In general, the Euler-Lagrange equations can have multiple solutions. However, we are interested in finding the unique solution that also satisfies the spectral constraint (A23). To find this particular solution, we use the two-step algorithm first introduced by Wen et al. (2013) in the context of porous medium convection. A remarkable property of this algorithm is that it eliminates the requirement of numerical continuation (Plasting & Kerswell 2003). As the two-step algorithm can be implemented at any value of the flow parameter, this flexibility has led to wider usage in several other studies of the background method to obtain the optimal bound numerically (Wen et al. 2015 Souza, Tobasco & Doering 2020). The first step of the algorithm uses a pseudo-time stepping scheme in which the Euler-Lagrange equations (A28a-d) are converted into a time-dependent system of partial differential equations (PDEs) as follows: where the index i ranges over the r, θ and z components ofṽ. Steady-state solutions of (5.1a-d) are equivalently solutions of the Euler-Lagrange equations (A28a-d). Note that we multiply the Fréchet derivatives with certain coefficients before introducing the time derivatives on the left-hand side. This makes the coefficient of the linear term (the Laplacian) a constant in the resultant time-dependent PDEs. Also, note that the coefficient in front of the Fréchet derivative with respect toṽ is positive, while the coefficients in front of the Fréchet derivatives with respect to U θ and a are negative. The reason is that we are maximizing the bound with respect toṽ while minimizing it with respect to U θ and a. Ding & Marensi (2019) proved that if the pseudo-time stepping scheme leads to a steady-state solution, then that solution must be the globally optimal solution of the Euler-Lagrange equations (A28a-d), i.e. the one that leads to optimal bounds. Conveniently, the same proof extends to the case where the perturbed flow satisfies only the homogeneous boundary conditions, and to the case where the perturbed flow is two-dimensional and incompressible. The proof of Ding & Marensi (2019) does not guarantee the existence of a steady-state solution to (5.1a-d). But in all the cases that we investigated, the pseudo-time stepping scheme did relax to a steady-state solution.
The second step of the two-step algorithm is a Newton iteration (see Wen et al. 2015) that has a faster convergence rate than the pseudo-time stepping scheme but requires a good initial guess. Naturally, we use the solution obtained at the end of the pseudo-time stepping scheme as the initial guess.
Solving the Euler-Lagrange equations in case 1 comes with two major simplifications. First, the pressure gradient term in (A28a) disappears, as we do not impose the incompressibility constraint on the perturbation. Second, it can be shown that the optimal perturbation depends only on the radial direction (see Appendix B). With these simplifications, the convergence of the pseudo-time stepping scheme is so rapid that the subsequent Newton iteration is not needed. Therefore, we use only the first step of the two-step algorithm described above. Furthermore, we found that it is also possible to solve the simplified Euler-Lagrange equations analytically in the limit Re → ∞ using the method of matched asymptotics (solutions are presented in Appendix B).
In case 2, it is also possible to make a simplification. Indeed, Ding & Marensi (2019) presented numerical evidence that the optimal solution does not depend on θ when the aspect ratio L (i.e. the height of the cylinder) is large enough. Therefore, we choose L = 20, which is sufficiently large to guarantee that the optimal flow is axisymmetric. To solve the system of time-dependent PDEs (5.1a-d), we consider the following Fourier decomposition in the z-direction: The radial direction is further discretized using the Chebyshev collocation method. We use a semi-implicit Crank-Nicolson scheme for the time integration, where we treat the linear terms implicitly and use the second-order Adams-Bashforth extrapolation for the nonlinear terms. We use an influence matrix method to solve for the pressure at each time step (see Peyret 2013, p. 236). The code is parallelized using MPI. Note that the pressurep in (5.2a,b), as compared to the one in Appendix A, has been multiplied with an appropriate factor such that it is precisely the gradient ofp that appears in the time-evolving PDEs (5.1a-d). Depending on the radius ratio and Taylor number considered, we vary the number modes in the z-direction from N = 200 to N = 6000, and the number collocation points in the r-direction from 120 to 320. The numerical strategy for solving the Euler-Lagrange equations in case 3 is similar to case 2 described above. The only difference is that for the 2-D incompressible perturbations, the flow quantities depend on the θ direction but are independent of z. Therefore, we consider the following decomposition instead: (5.3a,b) In this case, depending on the radius ratio and Taylor number considered, we vary the number modes in the θ-direction from M = 40 to M = 3000, and the number collocation points in the r-direction from 120 to 320.

Optimal bound results
In this subsection, we present the optimal bounds obtained using the numerical schemes described above for each of the three different sets of constraints on the perturbations. We begin by showing a typical optimal background flow profile at η = 0.6 and Ta = 10 6 in each case in figure 4. For comparison, we have also included the background flow profile constructed in (4.4) to derive the original analytical bound. As can be seen in figure 4, all four background flow profiles vary as cr, for some constant c, in the bulk region. This is expected intuitively as this type of background profile makes the sign-indefinite term (which is, in a loose sense, the hardest to control in the bulk region) in the spectral constraint (A23) zero. Near the cylinders, the background flows consist of two thin boundary layers. In order to meet the prescribed boundary conditions, the gradients in these thin layers are large, which makes the sign-indefinite term non-zero. However, as the perturbation has to satisfy the homogeneous boundary conditions, the net contribution from this term will still be smaller than the positive term in (A23) as long as the boundary layer thickness is small enough. In the optimal state, the boundary layers are of just the right size so that the positive term and the sign-indefinite term balance each other and the spectral constraint is marginally satisfied. When moving from case 1 to case 3, the restrictions on the perturbations increase, and this decreases the possibilities in which the sign-indefinite term can be negative. Therefore, the boundary layers become thicker, protruding more into the bulk region. Figure 5 shows the optimal bounds on the Nusselt number Nu b , as a function of the Taylor number Ta. We denote the bounds as Nu nc b for case 1 (figures 5a,b), Nu 3D b for case 2 (figures 5c,d), and Nu 2D b for case 3 (figures 5e,f ). We cover a wide range of parameters both in radius ratio (from η = 0.1 to η = 0.99) and in Taylor number. In figures 5(a), 5(c) and 5(e), the bound Nu b has been scaled with its expected asymptotic dependence on Ta, namely Ta 1/2 . The colour and shape of the symbols each correspond to a different radius ratio, as shown in the legend. The symbols in the plots in figure 5 correspond to data points computed using the numerical algorithm from the previous subsection, whereas the solid lines connecting the data points are calculated using interpolation, providing a guide to the eye. For every radius ratio value, the solid line is extended up to the highest Taylor number for which the computation is performed. Beyond this point, we extrapolate using a best fit of the form applied to the data of Nu b /Ta 1/2 computed from the last two decades in Ta. For each value of η, we thus define A(η) as the asymptotic limit of Nu b /Ta 1/2 as Ta → ∞. Table 1 summarizes the values of A(η) obtained from this fitting procedure for different radius ratios. We have added appropriate abbreviations in the superscript of A(η) to signify the case at hand. We remark that these extrapolations were necessary, especially for the small radius ratios, where the bound on the Nusselt number Nu b converges slowly to its asymptotic scaling in the Taylor number Ta. In figures 5(b), 5(d) and 5( f ), the bound Nu b has been scaled by Ta 1/2 as well as the geometrical scaling χ(η) obtained in (4.18). Note the striking collapse of the different radius ratio curves at high Taylor numbers in all three cases. Correspondingly, we also see from table 1 that the ratio A(η)/χ (η) is nearly independent of η, with less than 1.1 % variation in the average between the largest and smallest values. This suggests that the geometrical dependence of the bound on the Nusselt number at high Taylor number is χ(η) irrespective of the case considered. In case 1, the value of A nc (η)/χ (η) is close to 9/128, which is the exact asymptotic result that we obtained from the method of matched asymptotics in Appendix C. We also observe from table 1 that the value of A/χ (η) in cases 2 and 3 is very close to a constant for η > 0.5, but varies a little more for η < 0.5. This is likely due to the fact that the extrapolation is less accurate at small radius ratio because the computed data are further from being in the asymptotic regime compared with the case when the radius ratio is not small. For this reason, we assume that the average of A(η) calculated for η 0.5 is the correct asymptotic limit of Nu b /Ta 1/2 as Ta → ∞, and obtain Nu nc b,∞ = 9 8 Here, we have added '∞' in the subscript to point out that these are the main terms of the optimal bounds in the limit Ta → ∞.
In summary, we have shown that for case 1, case 2 and case 3, the optimal bounds are respectively a factor of 1.5, 12.46 and 27.66 better than the suboptimal bound (4.17) in the high Taylor number limit. Crucially, this improvement is uniform in the radius ratio η. We had obtained the analytical expression for the geometrical scaling χ(η) from a fairly simple suboptimal analytical bound calculated using a choice of background flow with two boundary layers whose thicknesses were adjusted to optimize the bound. During this procedure, we had not applied any constraint on the perturbed flowṽ (other than the homogeneous boundary conditions), and further, used standard calculus inequalities that are known to overestimate the bound on Nu b . Consequently, it is not at all self-evident why the optimal bounds should have the same geometrical scaling. The fact that the optimal bounds (5.5a-c), which are up to an order of magnitude better than the suboptimal bound (4.17), preserve the same geometrical dependence on radius ratio is therefore a simple yet remarkable result.

Wavenumber spectrum of perturbation
In this subsection, we investigate the wavenumber spectrum of the perturbed flowṽ with a particular focus on the small-scale structures present inṽ. In the optimal state,ṽ contains only a finite number of modes, called the critical modes, in either the z-or θ -direction, depending on the case considered, i.e. where K 2 ⊂ N and K 3 ⊂ Z \ {0} are finite sets. Moreover, as we will demonstrate below, the smallest scales in the perturbation are present only near the boundaries. It is thus reasonable to hypothesize that the smallest length scale in the perturbationṽ is similar to the boundary layer thickness of the background flow U. To pursue this idea further, we divide the critical modes present in the perturbation into four different categories. If, for a given critical mode, more than 90 % of the contribution to its L 1 (dr) norm comes from the region then we say that mode is active only near the inner cylinder. Similarly, if it comes from the region then we say that it is active only near the outer cylinder. Finally, if more than 90 % of the contribution comes from regions S in and S out together, then we say that the mode is active near both the cylinders; otherwise, we say that the mode is active in the bulk. This way of categorizing the modes may seem somewhat arbitrary at first, but looking at the shape of different critical modes, it becomes readily apparent that any other appropriate definition would have led to the same conclusion. We use the following colour scheme to differentiate modes according to our classification: blue for the modes that are active near the inner cylinder, red for the modes that are active near the outer cylinder, green for the modes that are active near both the cylinders, and black for the modes that are active in bulk. Figures 6(b,d, f ,h) show the plots ofṽ θ,n c (r) for critical modes at Ta = 10 8 and for radius ratios η = 0.2, 0.4, 0.6 and 0.8. We now see that the plots ofṽ θ,n c (r) provide an unambiguous visual justification of our earlier classification of critical modes into four categories, therefore our classification is robust. We first apply this categorization to the optimal perturbations found in case 2, denoting the wavenumber of the critical mode with smallest length scale that is active near the inner cylinder as k s in , and the one that is active near the outer cylinder as k s out . Assuming that our hypothesis about the similarity of the boundary layer thickness in the background flow and

A. Kumar
ṽ θ,n c Figure 6. (a,c,e,g) Wavenumbers of the critical modes in the optimal perturbation as a function of Ta at η = 0.2, 0.4, 0.6 and 0.8, respectively. The colour indicates if the critical mode is active near the inner cylinder (blue), outer cylinder (red) or both cylinders (green), and in the bulk (black), according to the classification given in the main text. The blue and red solid lines are the theoretical predictions for the critical mode with the largest wavenumber active near the inner and outer cylinders (see (5.11a,b)), respectively. (b,d, f,h) Corresponding azimuthal componentsṽ θ,nc (r) of critical modes at the same radius ratios, at Ta = 10 8 .
the smallest length scale in the perturbation made above is correct, we may expect

Optimal bounds in Taylor-Couette flow
where δ i and δ o are given by (4.15a,b). Substituting their expressions in the above relation leads to From these relations, we obtain the dependence of k s in and k s out not only on Ta, but also on the radius ratio η. In particular, we predict that the smallest length scales in the perturbation should become larger as η → 0. Furthermore, at a given η, the small-scale structures near the outer cylinder are predicted to be 1/η times larger than the ones near the inner cylinder. Figures 6(a,c,e,g) show the wavenumbers of the critical modes in the optimal perturbations as a function of the Taylor number for four different radius ratio values η = 0.2, 0.4, 0.6 and 0.8, respectively. By fitting these plots, we find that the constant of proportionality in (5.10a,b) that best fits the data at high Taylor numbers is C = 0.244, therefore we expect These two relations are plotted in figures 6(a), 6(c), 6(e) and 6(g) with solid blue and red lines, respectively. We see that the smallest length scales in the critical perturbation near the inner and outer boundaries, respectively, indeed follow the relations (5.11a,b). Furthermore, these smallest scales achieve their asymptotic scaling in Taylor number quicker than the corresponding optimal bounds on Nusselt number shown in figure 5, without any need for extrapolation of the data. We therefore argue that (5.11a,b) and figure 6 together provide a strong validation of the analytical predictions from § 4.
We can use similar ideas to predict the scaling of the smallest length scales in optimal perturbations in case 3. Using (4.14a-d), one would anticipate m s i ∝ r i /δ i and m s o ∝ r o /δ o , where m s i is the largest wavenumber of a critical mode active near the inner cylinder, and m s o is the largest wavenumber of a critical mode active near the outer cylinder. The plots in figures 7(a,c,e,g) show the wavenumbers of critical modes in the 2-D optimal perturbations as a function of the Taylor number at radius ratios 0.2, 0.4, 0.6 and 0.8. We apply the same mode identification method, and use the same colour scheme to differentiate the critical modes, as before. From these plots, we can fit the data at high Taylor numbers, to measure the constant of proportionality in the expressions for m s i and m s o , leading to We see that the wavenumbers of the critical mode with smallest length scale that is active near the inner and outer cylinders are equal. The relation (5.12), shown as a solid green line in figures 7(a,c,e,g), does seem to predict the largest wavenumbers at high Taylor numbers correctly. As in figure 6, figures 7(b,d, f ,h) Figure 7. (a,c,e,g) Wavenumbers m c of the critical modes that constitute the 2-D optimal perturbation as a function of the Taylor number for radius ratios η = 0.2, 0.4, 0.6 and 0.8, respectively. We use the same colour scheme as in figure 6 to distinguish different critical modes. The solid green line is the relation (5.12), which predicts the largest critical wavenumber. (b,d, f,h) Plots ofṽ c θ,mc , the coefficient of cos m c θ in the azimuthal component of the velocity, at Ta = 10 8 and the same radius ratios as in the (a,c,e,g) plots.
i.e. (5.14) dimensions, modes that are active solely near the cylinders oscillate in the boundary layer to ensure that (5.14) is satisfied.

A note on the applicability of the background method
In one of our previous studies (Kumar 2020), we presented a sufficient criterion to determine when the background method can be applied, for a given flow geometry and boundary conditions. We demonstrated that it can be used with any flow problem (tangential-velocity-driven or pressure-driven) with impermeable boundaries, provided that the boundaries have the shape of streamtubes of the flow where A is a constant skew-symmetric tensor, V 0 is a constant vector, and x is the position vector. For these types of problems, one can further show that the upper bound on the dissipation becomes independent of viscosity at high Reynolds numbers. In this section, we explore the complementary question of whether there exist flow configurations for which the background method cannot be applied. Indeed, the applicability of the background method depends on the existence of an incompressible background flow (which also satisfies the inhomogeneous boundary conditions) such that the following functional is positive semi-definite: for any perturbationsṽ that satisfy the homogeneous boundary conditions. Consequently, proving that the background method cannot be applied reduces to the problem of finding a perturbation or a family of perturbations such that there is no background flow U for which H is positive semi-definite. We start by giving a few examples where the applicability of the background flow can be dismissed rigorously. We first consider the case of Taylor-Couette flow with suction at the inner cylinder. The energy stability analysis of this problem was considered by Gallet, Doering & Spiegel (2010). The boundary conditions for this problem are where the Reynolds number Re is defined such that u · e r = −1 at the inner cylinder. The non-dimensional angular velocities of the inner and outer cylinders are ω i and ω o , respectively. In this problem, the flow is constricted to a narrow area as it moves from the outer cylinder (inlet) to the inner cylinder (outlet). We restrict ourselves to two dimensions, but the arguments given below are valid in three dimensions as well. The domain of interest We consider a perturbationṽ of the form v = (ṽ r ,ṽ θ ) = (0, v 0 r(r − r i )(r − r o )) , (6.4) whose amplitude is v 0 . Note thatṽ satisfies the homogeneous boundary conditions and is incompressible. We now demonstrate that for this perturbation, the spectral constraint can never be satisfied above a certain Reynolds number regardless of the choice of background flow U.
To show that the spectral constraint (6.2) is not satisfied, we have to show that the second term is negative and that its absolute value is larger than the first term. Being linear inṽ,

A. Kumar
the last term can be made arbitrarily small compared with the first two terms by choosing v 0 in (6.4) to be large enough for any given background flow U. As such, it does not play any role in the following argument.
The calculation of the first term is straightforward: In the calculation of the second term, we take advantage of the fact that the chosen perturbation (6.4) is independent of θ , so Now, using periodicity as well as the incompressibility condition satisfied by the background flow U, the following hold: Using (6.4) and (6.7a,b) in (6.6) gives From (6.5) and (6.8), we deduce that the spectral constraint (6.2) will not be satisfied if Re > 5r 2 i + 5r 2 o + 2 2r i , (6.9) a condition that is, remarkably, independent of the choice of U. Note that in the limit of r i /r o → 1, the Reynolds number beyond which the method fails goes to infinity. This limit recovers the case of a plane Couette flow with suction and injection at the walls (Doering, Spiegel & Worthing 2000), where the background method can indeed be applied, so (6.9) is consistent with these results. A similar type of condition on the Reynolds number can be derived in the problem of the Taylor-Couette flow with injection at the inner cylinder, i.e.
In this problem, the flow overall expands into a larger area as it moves from the inner cylinder (inlet) to the outer cylinder (outlet). For this case, we can use similar arguments but with the new perturbed flow v = (ṽ r ,ṽ θ ) = (v 0 r(r − r i )(r − r o ), 0) . (6.11) The perturbation this time is not incompressible but can be shown to yield a negative H(ṽ) regardless of the background flow U, for sufficiently large Reynolds number. However, noting that the perturbation (6.11) is radial, one may then expect that an incompressible 948 A11-24 perturbation, which is composed of vortices stretched in the radial direction, will also yield a negative H(ṽ). This observation led us to consider the streamfunction We define the corresponding velocity fieldṽ = (ṽ r ,ṽ θ ) as This velocity field is divergence-free and satisfies the homogeneous boundary conditions at the surface of the cylinders. The streamlines ofṽ are depicted in figure 8. Next, define a family of rotation operators Q ϕ : D σ (V) → D σ (V), indexed with ϕ, on the space of divergence-free vector fields that satisfies the homogeneous boundary conditions at ∂V as (6.14) A tedious calculation, first involving an integration in ϕ, and then using the arguments similar to the suction problem above, shows that This calculation implies that if (6.16) where · is the ceiling function, then for the integral (6.15) is negative, which implies that there is at least one ϕ ∈ [0, 2π] such that H(Q ϕ (ṽ)) < 0. More generally, there exists a set S ⊂ [0, 2π], depending on the background flow U, of positive measure (μ(S) > 0) such that H(Q ϕ (ṽ)) < 0 for any ϕ ∈ S, i.e. the spectral constraint is not satisfied. Note that condition (6.16) is basically saying that the vortices in the incompressible perturbed flow field (6.13a,b) should be stretched in the radial direction, which we expected from the example of the compressible perturbed flow (6.11).
The key message from these two problems is that if there is a converging flow, then one can rule out the applicability of the background method by creating a perturbation whose streamlines are perpendicular to the direction of the mean flow, while in the case of a diverging flow, one can use a perturbation whose streamlines are parallel to the direction of the mean flow instead. Of course, in both cases, we need to make sure that the perturbation satisfies the homogeneous boundary conditions. Combining these ideas suggests that one cannot apply the background method to flows in a converging-diverging nozzle, because one can choose the perturbation to be composed of vortices that stretch either in the perpendicular direction to the flow in the converging section or parallel to the flow in the diverging section. Using the same arguments, one would then also conjecture that the background method can in general not be applied to flows between rough walls. Indeed, in this case, one can always find vertical sections where the flow expands or compresses, and then one could use the same strategy to choose perturbations for which H(ṽ) < 0. However, note that in this case, the compression or expansion is small, i.e. the gap width on averages decreases by only a factor (1 − ) in the converging part, or increases by a factor (1 + ) in the diverging part, where is the non-dimensional roughness scale. This problem is analogous to the converging-diverging nozzle if the Reynolds number is based on the surface roughness . Therefore, for the Reynolds number based on the average gap width, we expect that the spectral constraint (6.2) will not be satisfied if Re −1 . This still leaves the problem open for the flow systems that do not have a converging or diverging section, for example, flow in tortuous channels. We believe that even for these problems, the spectral constraint will fail to hold past a certain Reynolds number for any background flow. Therefore, we conjecture that the sufficient condition for the applicability of the background method mentioned in the beginning of this section is also a necessary condition.

Summary and implications
In this paper, we computed optimal bounds on mean quantities in the Taylor-Couette flow problem with a stationary outer cylinder, with particular focus on the dependence of these bounds on the system geometry. Along the way, we studied the energy stability of the laminar flow in § 3. The main finding of this section was that for a value of radius ratio r i /r o below 0.0556, the marginally stable flow at the energy stability threshold is not composed of the well-known axisymmetric Taylor vortices but is instead a fully 3-D flow field.
To uncover the functional dependence of the optimal bounds on the radius ratio at large Taylor number, we began by deriving a suboptimal but analytical bound with the use of standard inequalities and a choice of background flow with two boundary layers (one near the inner cylinder and one near the outer cylinder) whose thicknesses were then adjusted to optimize the bound. We then argued that the dependence on the radius ratio captured by this analytical bound should also be the same for the optimal bounds at large Taylor numbers. We verified this statement systematically by obtaining distinct optimal bounds in three circumstances. In the first case, we imposed no constraints on the perturbation other than the homogeneous boundary conditions (case 1). Next, we allowed for 3-D incompressible perturbations (case 2). And finally, we considered 2-D incompressible perturbations (case 3). In the high Taylor number limit, we see an improvement of 1.5, 12.46 and 27.66, respectively, over the analytical bound as we move from case 1 to case 3, and that improvement is the same for all radius ratios. This result is striking and non-trivial because there is no known transformation of variables that makes the Euler-Lagrange equations (A28a-d) of the optimal bounds independent of the radius ratio.
In § 6, we rigorously dismissed the applicability of the background method for two flow problems. The limitation of the background method is previously known in the context of Rayleigh-Bénard convection at infinite Prandtl number (Nobili & Otto 2017), where it was shown that using a different method, a tighter bound can be obtained as compared to the background method. Here, we have shown that past a certain Reynolds number, no bound can be obtained using the background method applied to Taylor-Couette flow with suction or injection at the inner cylinder, i.e. there is no background flow that satisfies the spectral constraint even when the incompressibility condition on the perturbation is imposed. Generalizing these results then suggests that the spectral condition may not be satisfied for flow problems that contain converging or diverging sections, such as flow in a converging-diverging nozzle or flow between the rough walls. Possibly, the auxiliary functional method (Chernyshenko et al. 2014) could be a way forward to obtain bounds for such problems. The current best-known implementation of this auxiliary functional method to obtain a bound on energy dissipation in flow problems uses only quadratic functionals, in which case it has been shown to be equivalent to the background method (Chernyshenko 2022). However, Chernyshenko (2022) also proposed potential ways to implement non-quadratic functionals, which might be able to produce a finite bound on the energy dissipation in the problems stated above. Another possible approach, this time numerical, is to consider finite-dimensional truncated models of the actual flow problem, where there is a systematic way to use higher-than-quadratic auxiliary polynomials with the help of the sum-of-squares method (Goulart & Chernyshenko 2012;Huang et al. 2015;Fantuzzi et al. 2016;Goluskin 2018;Kumar 2019;Olson et al. 2021).
Our study brings to light the significance (or lack of significance, to be more precise) of the incompressibility constraint on the perturbation while calculating optimal bounds, especially in the limit of high Reynolds number, which is generally of interest in turbulent flows. As we showed in the present study, dropping the incompressibility constraint on the perturbations altogether (case 1) still recovers the correct dependence of the bounds on both the principal flow parameter (Ta, or equivalently Re) and the domain geometry (through the radius ratio). One cannot help but wonder whether the same holds true for other flow problems, including, for instance, the case of convection. It is a fundamental question of concern, as not imposing the incompressibility constraint tremendously decreases the computational cost of the optimal bound calculation. In the particular example studied here, in fact, not imposing the incompressibility condition allowed us to solve the Euler-Lagrange equations analytically using the method of matched asymptotics. This could also be potentially helpful in other studies involving the background method where it is relatively difficult to establish the scaling of the optimal bound even numerically, perhaps because the bounds involve logarithms Fantuzzi, Nobili & Wynn 2020) or a scaling other than a simple power law A. Kumar . In such situations, from an analysis point of view, it is usually challenging to decide what combination of the background field and calculus inequalities might be required to obtain the correct asymptotic scaling of the optimal bound. When the bound is obtained numerically instead, it can be difficult to establish its actual functional dependence on the principal flow parameter. Therefore, in these situations, considering case 1 can be very useful as one can try solving the Euler-Lagrange equations analytically using the method of matched asymptotics. These ideas can also be of relevance to other variational approaches, such as the wall-to-wall transport problem (Hassanzadeh, Chini & Doering 2014;Tobasco & Doering 2017;Motoki, Kawahara & Shimizu 2018a,b;Doering & Tobasco 2019;Souza et al. 2020;Tobasco 2022;Kumar 2022), which asks the question of what is the maximum heat transfer for a fixed energy or enstrophy budget.
Finally, assuming that the conclusions of this study apply more broadly, case 3 is also relevant to flow problems that are frequently investigated not just in three dimensions but in two dimensions as well, such as Rayleigh-Bénard convection or internally heated convection. For example, it could be interesting to determine how the optimal bounds depend on the shape of the roughness of the wall in 2-D Rayleigh-Bénard convection with rough boundaries, a problem investigated previously by Goluskin & Doering (2016) using the background method.
Before proceeding further, we note that one can also use the direct method of Seis (2015) to derive an upper bound on the Nusselt number with the same Taylor number dependence as in this paper. However, a question of interest could be if a bound with the same geometrical scaling can also be derived. We would expect that if one makes estimates near both the cylinders in Seis's approach, and uses an analogous optimization procedure (similar to this paper), then one may be able to derive a suboptimal bound with the same geometrical scaling as in this paper.

Comparison with the DNS
We now analyse briefly our results from a more practical point of view, and ask the question of whether the dependence of the Nusselt number on the radius ratio obtained in this paper bears any relationship with that of the actual turbulent flow. Note that the asymptotic dependence of the optimal bound on the Taylor number is known to overestimate the actual Nusselt number in turbulent Taylor-Couette flows by a logarithmic factor in Ta (Grossmann et al. 2016). As such, we cannot compare directly our results to the data, but instead merely ask the question of whether the geometric prefactor g(η) in the expression Nu(η, Ta) = g(η) f (Ta) measured in turbulent Taylor-Couette flows bears any resemblance to the prefactor χ(η) obtained in our optimal bound calculation; see (4.18).
We first test this idea on the direct numerical simulations (DNS) data from Ostilla-Mónico et al. (2014) and Froitzheim et al. (2019). In figure 9(a), we have plotted Nu versus Ta from these DNS, and in figure 9(b), we show the same data divided by χ(η). We see that the rescaled data do become more compact and appear to fall on a single curve. This observation gives us confidence that the geometrical dependence of the bound χ(η) obtained in this paper is a good approximation to that of the actual Nusselt number Nu measured in turbulent Taylor-Couette flows. However, we note that the data have not yet reached the asymptotic scaling corresponding to the high Ta regime, so the comparison at this point remains tentative. We also note that a different prediction for Nu(η, Ta) has been obtained recently by Berghout et al. (2020) using the idea of Monin-Obukhov theory for thermally stratified turbulent boundary layers. Their scaling of the asymptotic limit of high Ta is given as where κ = 0.39 is the von Kármán constant. The geometrical dependence in (7.1) differs from χ(η) by a factor (1 + η 2 ) 2 . However, it is reassuring to see that both expressions are proportional to η 3 in the limit of small radius ratio. A definitive answer to the question of whether the geometrical scaling χ(η) given by our bound is exact or just an approximation would require a precise comparison with the turbulent data at very high Taylor numbers collected for a range of radius ratios spanning the entire interval (0, 1), which is at present a challenge for the numerical computations.
7.3. Further generalizations We end this paper by discussing a few important consequences and generalizations of our study as well as future outlooks. The first of these consequences concerns the bound on dissipation. The optimal bound on the Nusselt number for case 2 (3-D incompressible perturbations) combined with the relations (2.24) and (2.5) gives us the optimal bound on the dissipation This bound tends to 0.00846 in the limit η → 1, which is within 1 % of the optimal bound obtained by Plasting & Kerswell (2003) for the plane Couette flow, namely 0.008553. The consistency between the two results shows that our work can, in retrospect, be viewed as a generalization of the result of Plasting & Kerswell (2003) to Taylor-Couette flow for an arbitrary radius ratio. The second item is related to our previous work (Kumar 2020) on the dependence of the bound on the friction factor λ on the radius of curvature κ and torsion τ for a pressure-driven flow in a helical pipe. We were able to employ a boundary layer optimization technique together with standard inequalities similar to that used here to A. Kumar obtain the following analytical bound on the friction factor in high Re limit: However, the complexity of the helical pipe geometry makes it impossible in practice to compute the corresponding optimal bound. Nevertheless, in the light of results from the present study, and assuming that we captured the geometrical dependence correctly, one can in principle compute the prefactor in a limit where the optimal bound can be computed, namely the case of a straight pipe, for which κ = τ = 0. This bound was computed by Plasting & Kerswell (2005) to be λ 3D b,∞ (0, 0) = 0.27, and using this result, we then expect that the optimal bound for helical pipes in the limit of high Reynolds number is λ 3D b,∞ = 0.27I(κ, τ ). (7.5) Finally, the results presented in this paper potentially open the door to solving many important outstanding problems in engineering. Indeed, within that context, we are often interested in finding the optimal geometry of the system or the object involved that minimizes or maximizes a certain flow quantity subject to some physical constraint. These types of problems therefore demand a careful study of the effect of the domain shape on a flow quantity. From this perspective, our study has broader implications. Even though we ruled out the applicability of the background method to a large class of problems (see § 6), this still leaves a number of interesting problems open for analysis. For example, two problems that have been investigated using DNS before, but where an application of the background method can provide further insights, are the Taylor-Couette flow with axisymmetric grooved walls (Zhu et al. 2016), and pressure-driven flow in a pipe with an elliptic cross-section (Nikitin & Yakhot 2005). Another problem where the background method has been used previously but capturing the exact domain shape dependence in the bounds was not the primary focus is the flow of fluid in an arbitrary domain driven by moving boundaries (Wang 1997). Our study suggests an interesting avenue towards solving these problems, by using the background method together with perturbations that are not assumed to be incompressible, which, as we demonstrated here, can greatly simplify the calculation.
Funding. This paper is dedicated to Charlie Doering, whose work has been instrumental in motivating the author's research. A.K. thanks P. Garaud for a careful read of the paper and for providing comments that improved the quality of the paper. A.K. acknowledges use of the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315.
Declaration of interests. The author reports no conflict of interest.

Appendix A. The background method
In this appendix, we formulate the background method to obtain an upper bound on the quantity 1 Re ∇u 2 2 (A1) 948 A11-30 in Taylor-Couette flow. It is clear from (2.19) that an upper bound on this quantity immediately provides an upper bound on the dissipation ε. We begin by writing the total flow field u as a sum of two divergence-free flow fields: We call U the background flow, and require that it satisfies the same boundary conditions as u, and is a function of only space. We call v the perturbation, or perturbed flow, which satisfies homogeneous boundary conditions. The governing equation for the perturbation, obtained by substituting (A2) in (2.3), is given by ∂v ∂t + U · ∇U + U · ∇v + v · ∇U + v · ∇v = −∇p + 1 Re We then obtain the evolution equation of the energy in the perturbed flow by taking the dot product of (A3) with v and integrating over the volume: (A4) Now using integration by parts, we can write where, in the index notation, At the same time, one also has the identity ∇U : ∇v = |∇u| 2 − |∇U| 2 − |∇v| 2 2 .
Using (A5) and (A7) The introduction of a balance parameter a in the background formulation goes back to Nicodemus, Grossmann & Holthaus (1997). Now it can be shown within the framework of the background method that the quantity v 2 2 is uniformly bounded in time (see Doering & Constantin 1992, for example). As a result, the long-time average of the time derivative

A. Kumar
Now, to optimize the bound (A21) under the incompressibility constraint onṽ, we write the following Lagrangian: Letting the first variation (the Fréchet derivative) of this functional with respect toṽ,p, U and a tend to zero leads to In general, these equations do not have a unique solution. However, the solution to these equations for which the background flow also satisfies the spectral constraint (A23), or equivalently, all the eigenvalues of the eigenvalue problem (A26a,b) are non-negative, is unique.

Appendix B. A useful lemma
Here, we prove that the marginally stable perturbations in the energy stability analysis of § 3 or optimal perturbations in § 5 depend on radius only when they are not required to be incompressible.
LEMMA B.1. Let D(V) be the set of smooth velocity fields (not necessarily incompressible) that satisfy the homogeneous boundary conditions. For a given choice of the balance parameter 0 < a < 2 and of the unidirectional background flow U = U θ (r) e θ , the functional H(ṽ) (given by (A22)) achieves a minimum whenṽ is a function of the radial direction only. Furthermore, if the background flow satisfies dU θ /dr − U θ /r 0, then the optimal perturbed flow corresponds toṽ r =ṽ θ .
REMARK B.2. Although we do not prove that the optimal background flow satisfies dU θ /dr − U θ /r 0, this condition was found to hold in every numerical computation of optimal bounds in all the three cases considered in our paper, as well as for the choice of the background flow in analytical construction presented in § 4. Therefore, it is natural to make the assumption that dU θ /dr − U θ /r 0.
In the second part, for every perturbationṽ = (ṽ r ,ṽ θ ), we define a modified perturbation