On the optimal design of wall-to-wall heat transport
Charles R. Doering Ian Tobasco
math.AP physics.flu-dyn








We consider the problem of optimizing heat transport 
through an incompressible fluid layer. Modeling passive scalar transport
by advection-diffusion, we maximize the mean rate of
total transport by a divergence-free velocity field. Subject
to various boundary conditions and intensity constraints, we prove
that the maximal rate of transport scales linearly in the r.m.s.kinetic energy and, up to possible logarithmic corrections, as the
1/3rd power of the mean enstrophy in the advective regime. This
makes rigorous a previous prediction on the near optimality of convection
rolls for energy-constrained transport. Optimal designs for enstrophy-constrained
transport are significantly more difficult to describe:
we introduce a “branching” flow design with an unbounded number
of degrees of freedom and prove it achieves nearly optimal transport.
The main technical tool behind these results is a variational principle
for evaluating the transport of candidate designs. The principle
admits dual formulations for bounding transport from above and below.
While the upper bound is closely related to the “background method”,
the lower bound reveals a connection between the optimal design problems
considered herein and other apparently related model problems from
mathematical materials science. These connections serve to motivate
 designs.



On the optimal design of wall-to-wall heat transport
    Charles R. Doering and Ian Tobasco
    December 30, 2023
====================================================





§ INTRODUCTION




 §.§ The wall-to-wall optimal transport problem


This paper concerns a class of optimal design problems from fluid
dynamics that asks to maximize the overall transport of heat through
an incompressible fluid layer. Passive scalar transport by an incompressible
fluid is governed by the advection-diffusion equation

    ∂_tT+𝐮·∇ T=κΔ T

where T(𝐱,t) is the scalar field undergoing transport,
referred to as temperature throughout, 𝐮(𝐱,t)
is the velocity vector field of the fluid, and κ is the coefficient
of molecular diffusivity. In general, the velocity field 𝐮
and temperature T may depend on both space 𝐱=(x,y,z)
and time t. Due to incompressibility, 𝐮 must remain
divergence-free. Thinking of 𝐮 as being in our control,
we set ourselves the task of choosing it to maximize the overall transport
of heat determined by eq:advection-diffusion-kappa.

This is a rich class of optimal design problems and we are interested
in the dependence of any solutions, i.e., optimal designs,
on various constraints that may be imposed. We discuss specific constraints
for the velocities later on, but let us handle the temperature field
first. Supposing the fluid is contained between two impenetrable parallel
planar walls at a distance h, we fix the temperature at the walls
by imposing the constant Dirichlet boundary conditions

    T|_z=0=T_hot  and   T|_z=h=T_cold.

If the velocity field 𝐮 is regular enough —eventually
our constraints on it will ensure this —the advection-diffusion
equation eq:advection-diffusion-kappa admits a unique
solution T satisfying eq:Dirichlet-dimesional for
every essentially bounded initial temperature field T|_t=0=T_0(𝐱).
We see, therefore, that the overall heat transport specified by eq:advection-diffusion-kappa
should depend in general on 𝐮 and T_0. However, as
the partial differential equation (PDE) eq:advection-diffusion-kappa
is dissipative, any dependence on the initial temperature T_0
is eventually lost as t→∞, and the resulting heat transport
can be thought of as being set by 𝐮 alone.

In this paper, we study the optimal design of wall-to-wall heat transport
in the long-time limit, subject to various boundary conditions and
intensity constraints on the velocity field 𝐮. To simplify
matters, we consider all fields to be periodic in the wall-parallel
variables x and y with periods l_x and l_y. That is,
we take 𝐱 to belong to the domain Ω=𝕋_xy^2×[0,h]_z,
identifying 𝕋_xy^2 with [0,l_x]×[0,l_y]
in the usual way. We turn now to discuss the precise measure of overall
heat transport that will be optimized throughout. 



  §.§.§ Finite-time wall-to-wall optimal transport


According to the advection-diffusion equation eq:advection-diffusion-kappa
and the boundary conditions eq:Dirichlet-dimesional,
the vertically averaged rate of heat transport per unit area up to
time t=τ is given by

    J_τ=1/τ1/l_xl_yh∫_0^τ∫_Ωk̂·(𝐮T-κ∇ T) d𝐱dt=κ/h(T_hot-T_cold)+_0^τ_ΩwT d𝐱dt.

Here, 𝐮=uî+vĵ+wk̂ and  denotes
an average over the integration domain. We are interested in determining
those velocity fields which maximize the overall heat transport J_τ.
Of course, unless 𝐮 is suitably constrained the optimal
transport sup J_τ will be infinitely large. It is natural to prescribe the overall
magnitude of 𝐮, and to enforce suitable boundary
conditions at the walls ∂Ω. The resulting optimal
design problems take the form

    sup_𝐮(𝐱,t)
    
    ||𝐮||=U
    
    +b.c.
     J_τ

where the parameter U sets the advective intensity of the admissible
velocity fields.

We consider two classes of admissible velocity fields which we refer to as being “energy-”
or “enstrophy-constrained”. In the energy-constrained class, we take

    ||𝐮||^2=_0^τ_Ω|𝐮|^2 d𝐱dt

in eq:finite-timepblm so that the constraint
||𝐮||=U sets the average kinetic energy available for
advection. As for boundary conditions, the no-penetration ones

    w|_∂Ω=0

are well-suited to this class. We call the problem that results the 
finite-time energy-constrained wall-to-wall optimal transport
problem. The finite-time enstrophy-constrained[For various boundary conditions including the ones considered here,
the mean square rate of strain ||∇𝐮||_L^2(Ω)^2
and enstrophy ||∇×𝐮||_L^2(Ω)^2 are
the same.] problem arises from taking

    ||𝐮||^2=h^2_0^τ_Ω|∇𝐮|^2 d𝐱dt

and enforcing the no-slip boundary conditions

    𝐮|_∂Ω=0

in eq:finite-timepblm. The essential results of this
paper hold as well for the stress-free boundary conditions

    w|_∂Ω=0  and  ∂_zu|_∂Ω=∂_zv|_∂Ω=0,

although our focus is mostly on the no-slip ones. 









  §.§.§ Infinite-time wall-to-wall optimal transport


Having introduced the finite-time energy- and enstrophy-constrained
wall-to-wall optimal transport problems, we turn to discuss their
infinite-time analogs which are the focus of this paper. Let ⟨·⟩ denote
the (limit superior) space and long-time average

    ⟨ f⟩ =lim sup_τ→∞ _0^τ_Ωf(𝐱,t) d𝐱dt.

As an integration by parts shows, the space and long-time averaged
heat transport determined by eq:advection-diffusion-kappa
satisfies

    lim sup_τ→∞ J_τ=κ⟨ |∇ T|^2⟩ .

Note this depends on 𝐮 but not on the initial temperature
T_0 so long as it is bounded. 
In direct analogy with the finite-time optimal transport
problems, we define the infinite-time energy-constrained wall-to-wall
optimal transport problem by

    sup_𝐮(𝐱,t)
    ⟨ |𝐮|^2⟩ =U^2
    
    w|_∂Ω=0
     κ⟨ |∇ T|^2⟩

and the infinite-time enstrophy-constrained problem by

    sup_𝐮(𝐱,t)
    ⟨ |∇𝐮|^2⟩ =U^2/h^2
    𝐮|_∂Ω=0 κ⟨ |∇ T|^2⟩ .

It is these infinite-time optimal design problems that we study in
the remainder of this paper. As we never return to the finite-time
problems, we discontinue the use of the distinguishing phrases from
now on.

A word is in order regarding the sense in which we consider eq:energy-constrained-in-units
and eq:enstrophy-constrained-in-units to be solved. We
do not claim that there must exist maximizers for either problem. 
Although this certainly merits investigation, and is related
to questions of Γ-convergence <cit.> of
the finite-time problems to the infinite-time ones, we choose in this
paper to focus instead on the maximum value of transport which
is always well-defined. To the maximum value is associated maximizing
sequences, i.e., near optimizers which we may seek to describe. Even
in the steady versions of eq:energy-constrained-in-units
and eq:enstrophy-constrained-in-units —where
all fields are assumed to be independent of time and optimal designs
are guaranteed to exist —determining the maximal transport
is a non-trivial task. 

The energy- and enstrophy-constrained wall-to-wall optimal transport
problems eq:energy-constrained-in-units and eq:enstrophy-constrained-in-units
were introduced in <cit.> and studied further
in <cit.> by a combination of asymptotic and numerical methods.
Similar methods have since been applied to study other related optimal
transport problems <cit.>.
A key question left unresolved by these works is whether the local
maximizers constructed therein actually achieve heat transport comparable
to that of global optimizers. In this paper, we present a new mathematically
rigorous approach to answering this question. Our methods do not rely
on the use of Euler-Lagrange equations; as these are non-concave maximization
problems with many local maximizers, critical point conditions do
not suffice to identify global optimizers. Rather, our starting point
is a new variational formula for evaluating wall-to-wall heat transport,
which is useful both for proving a priori upper bounds on optimal
transport as well as lower bounds on the transport of candidate designs.
For the energy-constrained problem, we prove that the convection roll
designs from <cit.> achieve globally optimal
heat transport up to a universal prefactor in the advection-dominated
regime. For the enstrophy-constrained problem, we construct a new
class of “branching” designs featuring a large and potentially
unbounded number of degrees of freedom. A well-chosen branching design
achieves optimal transport up to possible logarithmic corrections. 

The wall-to-wall optimal transport problem is naturally related to
the study of transport in turbulent fluids. One consequence of our
results is a proof that any flows arising in Rayleigh's original two-dimensional
model of buoyancy-driven convection between stress-free walls
<cit.> must achieve significantly sub-optimal rates
of heat transport in the large Rayleigh number regime Ra≫1. Indeed,
while our results imply the existence of incompressible flows achieving transport 
consistent with the proposed “ultimate scaling” law Nu∼ Ra^1/2 (up to logarithmic corrections), 
such transport is impossible in Rayleigh's original model <cit.>.
In fact, our analysis leads us to wonder whether such logarithmic corrections to
scaling should always hold, independent of dimension or boundary
conditions. Behind these claims is a more or less explicit connection
between the fluid dynamical optimal design problems considered herein
and other apparently related model problems from the study of “energy-driven
pattern formation” in materials science <cit.>.
We discuss these considerations in detail at the end. A preliminary version of our methods and results was announced in <cit.>.




 §.§ Main results and methods




  §.§.§ Non-dimensionalization


We are concerned with the dependence of energy- and enstrophy-constrained
wall-to-wall optimal transport  eq:energy-constrained-in-units
and eq:enstrophy-constrained-in-units in their parameters. 
We make use of two standard non-dimensional quantities.
The Peclét number

    Pe=Uh/κ

is a dimensionless measure of the intensity of advection relative
to that of diffusion. Transport by eq:advection-diffusion-kappa
is dominated by advection when Pe≫1 and by diffusion when Pe≪1.
The Nusselt number Nu is a dimensionless measure of the enhancement
of heat transport by convection over that of pure conduction. In the
fluid layer geometry, 

    Nu(𝐮)=κ⟨ |∇ T|^2⟩/κ/h^2(T_hot-T_cold)^2=1+h/κ1/T_hot-T_cold⟨ wT⟩ .

Note this does not depend on the initial temperature T_0. 

Such non-dimensionalization reduces the number of free parameters
in eq:energy-constrained-in-units
and eq:enstrophy-constrained-in-units to three: the dimensionless group Pe and the aspect
ratios of the domain h/l_x and h/l_y. That
is, it suffices to take

    h=κ=T_hot=1  and   T_cold=0

and study the dependence of the resulting non-dimensionalized optimal
transport problems

    sup_𝐮(𝐱,t)
    ⟨ |𝐮|^2⟩ =Pe^2
    
    w|_∂Ω=0
     Nu(𝐮)    and    sup_𝐮(𝐱,t)
    ⟨ |∇𝐮|^2⟩ =Pe^2
    𝐮|_∂Ω=0 Nu(𝐮)

on Pe, l_x, and l_y. Henceforth, we understand the Nusselt number to be given by

    Nu(𝐮)=⟨ |∇ T|^2⟩ =1+⟨ wT⟩

where T is determined from 𝐮 by solving the advection-diffusion
equation

    ∂_tT+𝐮·∇ T=Δ T

with Dirichlet boundary conditions 

    T|_z=0=1  and   T|_z=1=0

and any essentially bounded initial data T|_t=0=T_0 (the choice
of which is immaterial to our results). The domain Ω=𝕋_xy^2× I_z
where 𝕋_xy^2 is identified with [0,l_x]×[0,l_y]
and I_z=[0,1]. As always, 𝐮 is understood to be divergence-free. 



  §.§.§ Summary of main results


Our results concern the asymptotic dependence of optimal transport
in the advective regime Pe≫1. Concerning energy-constrained transport,
we find that the maximal transport rate scales linearly in the
r.m.s. kinetic energy as Pe→∞. More precisely, we prove
the following result:

 There exist positive constants C
and C' so that

    CPe≤sup_𝐮(𝐱,t)
    ⟨ |𝐮|^2⟩ =Pe^2
    
    w|_∂Ω=0
     Nu(𝐮)≤1/2Pe

for all Pe≥ C'. The constant C is independent of all parameters
and C' depends only on the aspect ratios of the domain. 

As noted in <cit.>, the a priori upper bound Nu≲ Pe
can be proved by a quick application of the maximum principle and
the Cauchy-Schwarz inequality. On the other hand, to prove the lower
bound one must construct a certain family of admissible velocity fields
{𝐮_Pe} and prove that their Nusselt numbers scale
linearly in Pe in the advective regime. Such a construction was
described using methods of matched asymptotic analysis in <cit.>
(albeit with no attempt to control the errors in the ensuing estimates).
Our construction is inspired by that one: we consider a convection
roll system as in fig:rollsrough and choose the number
of rolls to scale optimally in Pe. Our approach to evaluating Nu allows to rigorously justify the predictions from <cit.>
regarding the (near) optimality of such flows.

The enstrophy-constrained problem turns out to be much more
difficult to resolve. We prove that the maximal enstrophy-constrained transport rate scales, up to possible
logarithmic corrections, as the 2/3rd power of the r.m.s. rate-of-strain
as Pe→∞. Furthermore, we obtain a bound on
the size of any corrections to this scaling:

 There exist positive constants C,
C', and C” so that


    C1/log^4/3PePe^2/3≤sup_𝐮(𝐱,t)
    ⟨ |∇𝐮|^2⟩ =Pe^2
    𝐮|_∂Ω=0
     Nu(𝐮)≤ C'Pe^2/3

for all Pe≥ C”. The constants C and C' are independent
of all parameters and C” depends only on the aspect ratios of
the domain.

 
The same bounds apply to enstrophy-constrained optimal transport between no-penetration or stress-free walls.
Indeed, by a simple inclusion argument, maximal transport between impenetrable walls is never less than for stress-free walls, and both are bounded below by maximal transport between no-slip walls.
Since the a priori upper bound Nu≲ Pe^2/3 applies so long as w|_∂Ω=0 (this is what is proved in sec:aprioribds), the result follows.

This result concerning the 2/3-scaling law of enstrophy-constrained
wall-to-wall optimal transport —modulo logarithms —was
first announced in our paper <cit.>. The present
paper provides all the mathematical details of the analysis outlined
there, as well as a much more complete discussion of our general approach
to the optimal design of heat transport. The bulk of it is
devoted to motivating and evaluating the branching flows depicted
in fig:branchingrough, which are the key to proving the
logarithmically corrected lower bound from thm:mainbounds_w2w_ens.

After this work was completed, a computational study of the Euler-Lagrange equations for the enstrophy-constrained wall-to-wall optimal transport problem reported convincing numerical evidence for velocity fields that produce Nu ∼ Pe^2/3 in three dimensions <cit.>. Interestingly, numerical studies of the two-dimensional problem have thus far failed to produce heat transport scaling of this sort <cit.>. Whereas the velocity fields produced in these two-dimensional studies feature near-wall “recirculation zones”, which serve to enhance heat transport at moderate Pe, they come nowhere near the complexity of our branching flows. The three-dimensional computations, however, do exhibit branching of a fully three-dimensional character. Whether or not such three-dimensional branching flows can be constructed so as to eliminate the logarithmic gap in thm:mainbounds_w2w_ens as Pe→∞ remains to be seen.






  §.§.§ Outline of the approach


thm:mainbounds_w2w_eng and thm:mainbounds_w2w_ens
contain two types of statements: a priori upper bounds on the
Nusselt number Nu that hold for all velocity fields, and matching
lower bounds on Nu for suitable designs. Methods to establish rigorous
upper bounds on convective transport go back at least to Howard in
the context of turbulent buoyancy-driven convection <cit.>,
and Constantin and one of the authors who developed the “background method”
to prove upper bounds on Nu (albeit absent Howard's hypothesis
of statistical stationarity) <cit.>.
Although a suitably adapted background method can be applied here
<cit.> we do not proceed in this way. Instead, we present a new
method for establishing upper bounds based on the fact that, for steady
velocity fields, there exists a variational principle for evaluating
heat transport. In the time-dependent case, this leads to new variational
bounds on Nu that imply the background method. The bound we obtain
is as follows:

    Nu(𝐮)≤inf_η ⟨ |∇η|^2+|∇Δ^-1[∂_tη+div(𝐮η)]|^2⟩

where η must satisfy

    η|_z=0=1  and  η|_z=1=0.

Here and throughout Δ^-1 denotes the inverse Laplacian operator
with vanishing Dirichlet boundary conditions. The bound eq:unsteady_ub
is sharp for steady flows; in that case eq:unsteady_ub
becomes an equality and η need not depend on time.

In contrast, methods to establish rigorous lower bounds on Nu are
far and few between. The righthand side of eq:unsteady_ub
is a convex minimization. Therefore, on general grounds, there should
exist a concave maximization that is its dual. We find that

    Nu(𝐮)-1≥sup_ξ ⟨ 2wξ-|∇ξ|^2-|∇Δ^-1[∂_tξ+div(𝐮ξ)]|^2⟩

where ξ must satisfy

    ξ|_z=0=0  and  ξ|_z=1=0.

As with eq:unsteady_ub, the bound eq:unsteady_lb
becomes sharp for 𝐮 that do not depend on time. 

Armed with these observations, we describe a new duality-based approach
to producing candidate designs. Consider the steady enstrophy-constrained
wall-to-wall optimal transport problem

    max_𝐮(𝐱)
    _Ω|∇𝐮|^2=Pe^2
    𝐮|_∂Ω=0
    Nu(𝐮),

whose optimal value bounds the unsteady maximum from below. Appealing
to the steady version of eq:unsteady_lb, we find that
eq:enstrophy-constrainedpblm-steady can be rewritten
as

    min_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    𝐮|_∂Ω=0,ξ|_∂Ω=0
     _Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ_Ω|∇𝐮|^2·_Ω|∇ξ|^2

where ϵ=Pe^-2. Indeed, the optimal values of eq:enstrophy-constrainedpblm-steady
and eq:optl_design_pblm are reciprocals and their optimizers
are in correspondence. Thus, solving the steady enstrophy-constrained
problem eq:enstrophy-constrainedpblm-steady for Pe≫1
is equivalent to solving eq:optl_design_pblm for ϵ≪1.
We refer to eq:optl_design_pblm as an “integral” formulation
of wall-to-wall optimal transport.

The family of variational problems eq:optl_design_pblm
is non-convex and singularly perturbed. The situation shares important
similarities with other model problems from the field of “energy-driven
pattern formation” in materials science <cit.>.
These include the study of branching patterns in micromagnetics <cit.>
and wrinkling cascades in thin elastic sheets <cit.>.
For such problems, it is known that certain patterns which, at a glance,
look like fig:branchingrough provide nearly optimal ways
of matching low energy states that are geometrically incompatible
but forced to coexist. We discuss such connections further in sec:energy-driven-pattern-formation.

Of course, eq:optl_design_pblm does not derive from materials
science but instead from fluid dynamics. We note the striking similarities
between it and Howard's variational problem, the latter of which gave
birth to the field of variational bounds on turbulent transport
<cit.>. It was recognized by Busse <cit.>
that Howard's problem should admit multiply-scaled optimizers. The
resulting construction is known as Busse's “multi-α” technique.
After suitable modifications (wall-to-wall optimal transport and Howard's
problem are in the end quite distinct) Busse's techniques can also be used to
study eq:optl_design_pblm. We consider these connections
further in sec:Implications-for-RBC.

By either analogy, we are led to construct self-similar branching
flows as candidates for eq:optl_design_pblm. The streamlines
depicted in fig:branchingrough are symmetric about z=1/2;
each half of the domain is made up of n convection roll systems
coupled through n-1 transition layers. In the bulk there are large
anisotropic convection rolls at some horizontal length-scale l_bulk.
Streamlines refine away from the bulk until there results an isotropic
convection roll system at some much smaller length-scale l_bl.
The entire construction can be modeled by a single length-scale function
ℓ(z) that interpolates through the layers. In terms of ℓ,
we find the optimal branching construction to be picked out by the
solution of 

    min_ℓ(z)
    ℓ(z_bulk)=l_bulk
    ℓ(z_bl)=l_bl l_bl+∫_z_bulk^z_bl(ℓ')^2 dz+ϵ(1/l_bulk^2+∫_z_bulk^z_bl1/ℓ^2 dz+1/l_bl)^2,

which satisfies

    ℓ(z)    ∼ϵ^1/6log^1/6 1/ϵ √(1-z)   z∈[z_bulk,z_bl],
    
    l_bulk   ∼ϵ^1/6log^1/6 1/ϵ, l_bl∼ϵ^1/3log^1/3 1/ϵ.

Although this analysis does not prove that optimal designs must exhibit
fluctuations according to these rules, it does yield designs sufficient
to obtain the asserted lower bounds from thm:mainbounds_w2w_ens.
The lower bounds from thm:mainbounds_w2w_eng on energy-constrained
optimal transport are much simpler to obtain, and serve as a test
case for our approach.


We close this introduction with a heuristic argument for the optimality
of our designs. Our streamlines branch according to the power law

    ℓ(z)∼ c(ϵ)√(z∧(1-z))

with a prefactor depending on ϵ=Pe^-2. Consider a cloud
of “hot” or “cold” particles released at one of the walls.
We can model their ensuing motion by solving the stochastic differential
equation 

    dX_t=𝐮(X_t)dt+dW_t

where W_t is the usual Brownian motion. Diffusion is a much slower
means of transport than advection. In the amount of time that it takes
a typical particle to advect vertically a distance ∼ d from
the wall, it will wander by diffusion horizontally a distance ∼√(d).
Therefore, to keep typical “hot” and “cold” particles separate
along their transit from wall to wall, one should “fit” the streamlines
of 𝐮 to the underlying diffusion by taking ℓ∼√(d).





 §.§ Outline of the paper


sec:aprioribds proves the a priori upper bounds
from thm:mainbounds_w2w_eng and thm:mainbounds_w2w_ens
and establishes the variational principles and bounds on Nu alluded
to above. The proof of the lower bounds from thm:mainbounds_w2w_eng
and thm:mainbounds_w2w_ens is spread across sec:OptlDesignProb_w2w,
sec:energy-constraineddesign, and sec:enstrophy_constrained_design.
In sec:OptlDesignProb_w2w we describe our general approach
to the optimal design of heat transport. In sec:energy-constraineddesign
we test our methods on the steady energy-constrained problem and obtain
a proof of the lower bound part of thm:mainbounds_w2w_eng.
In sec:enstrophy_constrained_design we consider the steady
enstrophy-constrained problem and prove the lower bound part of thm:mainbounds_w2w_ens.
We conclude in sec:Implications-for-RBC and sec:energy-driven-pattern-formation
with a discussion of bounds on turbulent heat transport, and
a discussion of wall-to-wall optimal transport as a problem of energy-driven
pattern formation. 



 §.§ Notation


Having non-dimensionalized, we employ the domain Ω=𝕋_xy^2× I_z
where 𝕋_xy^2 is identified with [0,l_x]×[0,l_y]
and I_z=[0,1]. The spatial average of an integrable function
f on Ω is denoted by

    _Ωf=1/|Ω|∫_Ωf(𝐱) d𝐱

where |Ω|=|𝕋_xy^2|=l_xl_y. Generally speaking,
 indicates a well-defined average over the indicated domain
of integration. Some distinguished averages used in this paper include

    f=_𝕋_xy^2f=1/|𝕋_xy^2|∫_𝕋_xy^2f(x,y,·) dxdy

which averages over the periodic variables x and y, the (limit
superior) space and long-time average

    ⟨ f⟩ =lim sup_τ→∞ _0^τ_Ωf=lim sup_τ→∞ 1/τ1/|Ω|∫_0^τ∫_Ωf(𝐱,t) d𝐱dt

and the truncated space and time average

    ⟨ f⟩ _τ=_0^τ_Ωf=1/τ1/|Ω|∫_0^τ∫_Ωf(𝐱,t) d𝐱dt.


We use the standard L^2- and Ḣ^1-norms for functions
on Ω,

    ||f||_L^2(Ω)=√(∫_Ω|f|^2)  and  ||f||_Ḣ^1(Ω)=√(∫_Ω|∇ f|^2).

The set of smooth and compactly supported functions on Ω is
C_c^∞(Ω). The Sobolev space H_0^1(Ω)
is its completion in the norm ||·||_Ḣ^1(Ω).
We use (·,·) to denote the duality pairing of
H^-1 with H_0^1. We denote by Δ^-1 the inverse
Laplacian operator with vanishing Dirichlet boundary conditions, which
is well-defined from H^-1(Ω)→ H_0^1(Ω). 

The notation X≲ Y means that there exists a positive constant
C not depending on any parameters such that X≤ CY. We use
the notations X∧ Y=min{X,Y} and X∨ Y=max{X,Y}.



 §.§ Acknowledgements


We thank R. V. Kohn and A. N. Souza for helpful discussions.
This work was supported by NSF Awards DGE-0813964 and DMS-1812831 (IT), DMS-1515161 and DMS-1813003
(CRD), a Van Loo Postdoctoral Fellowship (IT) and a Guggenheim Foundation
Fellowship (CRD).



§ A PRIORI BOUNDS ON WALL-TO-WALL OPTIMAL TRANSPORT


We begin our analysis of wall-to-wall optimal transport by proving
the a priori upper bounds from thm:mainbounds_w2w_eng
and thm:mainbounds_w2w_ens. Unless otherwise explicitly
stated, we consider throughout that ⟨ |𝐮|^2⟩ <∞
so that eq:advection-diffusion is well-posed.

The upper bound from thm:mainbounds_w2w_eng on energy-constrained
transport is straightforward to prove, and we dispatch with it first.

We have that

    Nu(𝐮)≤1+1/2⟨ |w|^2⟩ ^1/2

whenever w|_∂Ω=0.


Let us recall the argument from <cit.>. First,
note that Nu does not depend on the initial temperature T_0.
Thus, we can take T_0=1-z and conclude by the maximum principle
that the associated solution of eq:advection-diffusion
satisfies

    0≤ T≤1  a.e.


Note also that because w vanishes at ∂Ω,

    ⟨ w⟩ =0.

Combining this with Jensen's inequality and the definition of the
Nusselt number eq:Nusseltdefn, we have that

    Nu-1     =⟨ wT⟩ =⟨ w(T-1/2)⟩
       ≤⟨ |w|⟩ ||T-1/2||_L_tx^∞≤1/2⟨ |w|^2⟩ ^1/2.


The remainder of this section is on upper bounds for enstrophy-constrained
transport. We prove the following bound:

 There exists a positive constant C such
that

    Nu(𝐮)≤1+C⟨ |∇ w|^2⟩∧⟨ |∇ w|^2⟩ ^1/3.

whenever w|_∂Ω=0. This constant is independent of
all parameters.


To the authors' knowledge, there are at least three proofs of prop:ensbound_w2w.
For one, it can be obtained via an application of the background method
<cit.>. It can also be seen as a consequence of Seis' arguments
from <cit.>. Our proof of prop:ensbound_w2w
is different from either of these: we obtain it via a new approach
using a Dirichlet-type variational principle for the functional Nu(𝐮). 

It should be mentioned that we are not the first to notice the variational
structure of the advection-diffusion equation. The existence of a
variational principle for advection-diffusion in bounded domains appears
to have been first reported in <cit.>, where
it was used to systematically derive “best approximation” finite
element schemes. Around the same time, as described in <cit.>,
variational principles for computing effective complex conductivities
in periodic homogenization were discovered by Gibiansky and Cherkaev
(the relevant corrector equation is again divergence-form but not
self-adjoint). We learned about the existence of such principles from
the papers <cit.>,
whose formulas for computing effective diffusivities in periodic homogenization
inspired the formulas for Nu obtained below. Let us also mention
the related work <cit.> which discusses non-standard
variational principles for PDEs at large. It was a pleasant surprise
to learn that the seemingly ad hoc change of variables introduced
in <cit.> for handling the Euler-Lagrange equations
of wall-to-wall optimal transport turn out to be similar to those
employed in previous works, and that behind it all is a variational
principle for Nu. 

The remainder of this section is organized as follows. First we establish
a variational principle for Nu in the steady case where the reasoning
is most transparent. We then extend the arguments to general unsteady
flows, where the variational principle turns into a variational bound
as anticipated in eq:unsteady_ub. prop:ensbound_w2w
follows immediately thereafter. Later on in sec:OptlDesignProb_w2w,
we obtain the dual formula to bound Nu from below. In order to
highlight the key step in the proof —a certain symmetrizing
change of variables for the advection-diffusion equation —we
refer to this as the “symmetrization method”.



 §.§ The symmetrization method for steady velocity fields


We start with the case where 𝐮 is an arbitrary divergence-free
vector field belonging to L^2(Ω;ℝ^3). In this
case,

    Nu(𝐮)=1+_Ωwθ=1+_Ω|∇θ|^2

where θ=T-(1-z) is the deviation of the temperature field
from the conductive state. That is, θ is the unique (essentially
bounded) weak solution of 

    𝐮·∇θ=Δθ+w

with zero Dirichlet boundary data θ|_∂Ω=0. To
change variables, we let θ_± be the unique weak solutions
of the pair of formally adjoint PDEs

    ±𝐮·∇θ_±=Δθ_±+w

with θ_±|_∂Ω=0, and observe that θ=θ_+.
Then, we define η,ξ∈ H_0^1(Ω) by

    η=1/2(θ_+-θ_-)  and  ξ=1/2(θ_++θ_-)

and observe they satisfy the equivalent system of PDEs

    𝐮·∇η=Δξ+w
    𝐮·∇ξ=Δη.

We claim the change of variables (θ_+,θ_-)↔(η,ξ)
yields a variational formula for Nu.

Testing the second equation in eq:steadyPDEsystem against
ξ and integrating by parts shows that ∇ξ⊥∇η
in L^2(Ω), since

    ∫_Ω∇η·∇ξ=-∫_Ω𝐮·∇ξξ=0.

Therefore, 

    Nu-1=_Ω|∇θ_+|^2=_Ω|∇η|^2+|∇ξ|^2

or, using the first PDE in eq:steadyPDEsystem,

    Nu-1=_Ω|∇η|^2+|∇Δ^-1[𝐮·∇η-w]|^2.

Consider the righthand side of eq:Nu-non-local-formula
as it depends on η. Since 𝐮∈ L^2 and is divergence-free,
the righthand side is well-defined for η∈ H_0^1(Ω)∩ L^∞(Ω).
Now consider the variational problem

    inf_η∈ H_0^1(Ω)∩ L^∞(Ω) _Ω|∇η|^2+|∇Δ^-1[𝐮·∇η-w]|^2,

which is strictly convex so that any minimizer must be unique. By
a first variation argument, we see that η is a minimizer of
eq:steady_min-pblm_for_Nu if and only if it satisfies
the Euler-Lagrange equation

    Δη=𝐮·∇Δ^-1[𝐮·∇η-w].

This is a rewrite of the system eq:steadyPDEsystem with
ξ defined by 

    ξ=Δ^-1[𝐮·∇η-w].

Since that system possesses solutions in the class H_0^1∩ L^∞,
it immediately follows that eq:steady_min-pblm_for_Nu
has a minimizer. It also follows that the minimal value in eq:steady_min-pblm_for_Nu
is Nu-1. 

Relabeling η as η+1-z and using that 𝐮 is
divergence-free yields the following variational principle for heat
transport:

 Let 𝐮(𝐱) be
a divergence-free vector field in L^2(Ω;ℝ^3).
Then,

    Nu(𝐮)=min_η∈ H^1(Ω)∩ L^∞(Ω)
    η|_z=0=1, η|_z=1=0
     _Ω|∇η|^2+|∇Δ^-1div(𝐮η)|^2.





 §.§ Upper bounds on unsteady transport by symmetrization


With some care, the previous argument can be adapted to the time-dependent
case. Here, the useful change of variables arises from the pair of
PDEs

    ±(∂_t+𝐮·∇)θ_±=Δθ_±+w

which are formally adjoint in space and time. These are in obvious
analogy with eq:adjointsteadyPDEs from the steady case.
However, since parabolic PDEs are generically only well-posed forward
in time, making sense of the “-” equation presents an added
difficulty. To deal with this, we will reverse the sense of time between
the equations, performing the change of variables t→τ-t for
appropriately chosen τ≫1. In the limit τ→∞, we
recover the unsteady variational bound.

Define the admissible set of test functions

    𝒜={η:η∈ L^∞([0,∞);L^2(Ω))}
    ∩{η(t)∈ H^1(Ω)∩ L^∞(Ω),∂_tη(t)∈ H^-1(Ω) a.e. t}.



 Let 𝐮(𝐱,t) be
a divergence-free vector field with bounded mean energy ⟨ |𝐮|^2⟩ <∞.
Then,

    Nu(𝐮)≤inf_η∈𝒜
    η|_z=1=0, η|_z=0=1
     ⟨|∇η|^2+|∇Δ^-1[∂_tη+div(𝐮η)]|^2⟩ .



We begin by introducing the (approximately) symmetrized variables.
Let θ_+=T-(1-z) and note it solves

    ∂_tθ_++𝐮·∇θ_+=Δθ_++w

on [0,∞)×Ω and vanishes at ∂Ω. Let
θ_-^τ be the unique essentially bounded weak solution
of

    -∂_tθ_--𝐮·∇θ_-=Δθ_-+w

on Ω_τ=[0,τ]×Ω that vanishes at ∂Ω
and has final time data θ_-^τ(τ)=0. 



Next, define the symmetrized variables η^τ
and ξ^τ by

    η^τ=1/2(θ_+-θ_-^τ)  and  ξ^τ=1/2(θ_++θ_-^τ).

Observe η^τ and ξ^τ solve

    ∂_tη^τ+𝐮·∇η^τ=Δξ^τ+w
    ∂_tξ^τ+𝐮·∇ξ^τ=Δη^τ

on Ω_τ, vanish at ∂Ω, and remain bounded
in L_tx^∞ uniformly in time. In particular, by the maximum
principle,

    ||η^τ||_L_tx^∞∨||ξ^τ||_L_tx^∞≲||θ_+||_L_tx^∞∨||θ_-^τ||_L_tx^∞≤ C(T_0).


We proceed as in the steady case, accumulating errors that vanish
as τ→∞. From the second PDE in eq:symmetrized_timedept
we find that

    ⟨∇η^τ·∇ξ^τ⟩ _τ    =_0^τ1/|Ω|(-Δη^τ,ξ^τ) dt=_0^τ1/|Ω|(-∂_tξ^τ-div(𝐮ξ^τ),ξ^τ) dt
        =-_0^τ1/|Ω|d/dt1/2||ξ^τ||_L^2(Ω)^2+_0^τ_Ω𝐮ξ^τ·∇ξ^τ
        =1/2τ1/|Ω|(||ξ^τ(0)||_L^2(Ω)^2-||ξ^τ(τ)||_L^2(Ω)^2).

Therefore,

    |⟨∇η^τ·∇ξ^τ⟩ _τ|≤ C(T_0)1/τ.

Since 

    Nu-1=⟨ |∇θ|^2⟩ =⟨ |∇θ_+|^2⟩

and

    ⟨ |∇θ_+|^2⟩ _τ=⟨ |∇η^τ|^2+|∇ξ^τ|^2⟩ _τ+2⟨∇η^τ·∇ξ^τ⟩ _τ,

we conclude that

    Nu-1=⟨ |∇η^τ|^2+|∇Δ^-1[∂_tη^τ+𝐮·∇η^τ-w]|^2⟩ _τ+O(1/τ).


Next, we prove that η^τ is approximately minimal, with
an error that vanishes as τ→∞. Let η∈𝒜
vanish at ∂Ω but be otherwise arbitrary, and consider
the difference

    A_τ=⟨ |∇η|^2+|∇Δ^-1[∂_tη+𝐮·∇η-w]|^2⟩ _τ -⟨ |∇η^τ|^2+|∇Δ^-1[∂_tη^τ+𝐮·∇η^τ-w]|^2⟩ _τ.

Using the convexity of |·|^2, we can expand around η^τ
and use eq:symmetrized_timedept to arrive at the lower
bound

    A_τ/2   ≥⟨∇(η-η^τ)·∇η^τ+∇Δ^-1(∂_t+𝐮·∇)(η-η^τ)·∇Δ^-1[(∂_t+𝐮·∇)η^τ-w]⟩ _τ
        =-1/|Ω|_0^τ((∂_t+𝐮·∇)ξ^τ,η-η^τ)+((∂_t+𝐮·∇)(η-η^τ),ξ^τ) dt 
        =-1/|Ω|_0^τd/dt[∫_Ω(η-η^τ)ξ^τ] dt
       ≳-1/τ1/|Ω|(||η||_L_t^∞L_𝐱^2∨||η^τ||_L_t^∞L_𝐱^2)·||ξ^τ||_L_t^∞L_𝐱^2≥-C(T_0,η,Ω)1/τ.

Combining this with eq:Nu_unsteady_uptoerror, we find
that

    Nu-1≤⟨ |∇η|^2+|∇Δ^-1[∂_tη+𝐮·∇η-w]|^2⟩ _τ+C(T_0,η,Ω)1/τ.

Taking τ→∞ yields the inequality

    Nu-1≤⟨ |∇η|^2+|∇Δ^-1[∂_tη+𝐮·∇η-w]|^2⟩ .

This holds for all η∈𝒜 that vanish at ∂Ω.
Changing variables by η→η+1-z and optimizing yields the
result.

Even if 𝐮 depends on time, η can be taken to be
independent of time and still used to bound Nu. The simplified
version of thm:primal_Nubound_w2w that results is analogous
to thm:steady-primal-varprin, but for unsteady heat transport.

 Let 𝐮(𝐱,t) be
a divergence-free vector field with bounded mean energy ⟨ |𝐮|^2⟩ <∞.
Then,

    Nu(𝐮)≤inf_η∈ H^1(Ω)∩ L^∞(Ω)
    η|_z=1=0, η|_z=0=1
     _Ω|∇η|^2+⟨|∇Δ^-1div(𝐮η)|^2⟩ .





 §.§ Proof of prop:ensbound_w2w


We prove prop:ensbound_w2w by deducing it from cor:symmsteadybound_w2w
with a good choice of test function η. Evidently, we must find
a convenient way to bound the non-local term appearing there. Since
∇Δ^-1div is an L^2-orthogonal projection,
one has the bound

    ⟨|∇Δ^-1div(𝐮η)|^2⟩≤⟨|𝐮(η-c)|^2⟩

for an arbitrary constant c. In the case that 𝐮 satisfies
no-slip boundary conditions, one can deduce prop:ensbound_w2w
by choosing η≈ c thereby localizing the righthand side
to a small neighborhood of ∂Ω. Then, a straightforward
application of Poincaré's inequality yields the result.

The final step in the preceding argument requires all components of
𝐮 to vanish at ∂Ω. This is not useful for
dealing with no-penetration boundary conditions. Nevertheless, prop:ensbound_w2w
holds in this more general case. The key is to approach the non-local
term from cor:symmsteadybound_w2w by duality. Observe
that

    ∫_Ω|∇Δ^-1div 𝐦|^2=sup_θ∈ H_0^1(Ω) ∫_Ω2∇θ·𝐦-|∇θ|^2

whenever 𝐦∈ L^2(Ω). Thus, the inequality

    ∫_Ω|∇Δ^-1div 𝐦|^2≤ C

and the statement that

    ∫_Ω2∇θ·𝐦≤ C+∫_Ω|∇θ|^2  ∀ θ∈ H_0^1

are one and the same. Taking 𝐦=𝐮η where η
depends only on z, we conclude it will be useful to have bounds
of the form

    ∫_Ω2wθη'(z)≤ C+∫_Ω|∇θ|^2  ∀ θ∈ H_0^1.

The following preliminary result allows us to establish bounds of
this type.

Let w,θ∈ H_0^1(Ω).
Then,

    |wθ(z)|≲|z∧(1-z)|/|𝕋_xy^2|||∂_zθ||_L^2(Ω)||∂_zw||_L^2(Ω)  a.e.

 


The reader familiar with the background method may recognize that
this inequality also plays a key role in carrying out that approach
to a priori bounds. In particular, it is useful for verifying
the spectral constraint. See sec:onbackgroundmethod for
more on the connection between the symmetrization method and the background
method.


By the usual approximation arguments, we can take 𝐮 and
θ to be smooth. Differentiating and applying the Cauchy-Schwarz
inequality, we find that

    |d/dzw^2|≤2/|𝕋_xy^2|||∂_zw||_L^2(𝕋_xy^2)||w||_L^2(𝕋_xy^2).

Integrating this from 0 to z yields 

    ||w||_L^∞([0,z];L_xy^2)^2   ≤2/|𝕋_xy^2|∫_0^z||∂_zw(z')||_L_xy^2||w(z')||_L_xy^2 dz' 
       ≤2/|𝕋_xy^2|||∂_zw||_L^2(Ω)||w||_L^2([0,z];L_xy^2)
       ≤2/|𝕋_xy^2|||∂_zw||_L^2(Ω)||w||_L^∞([0,z];L_xy^2)|z|^1/2.

Therefore,

    ||w||_L^∞([0,z];L_xy^2)≤2/|𝕋_xy^2||z|^1/2||∂_zw||_L^2(Ω).

Similarly, we find that

    ||θ||_L^∞([0,z];L_xy^2)≤2/|𝕋_xy^2||z|^1/2||∂_zθ||_L^2(Ω).


Now consider the product wθ. We have that

    |d/dzwθ|≤|∂_zwθ|+|w∂_zθ|≤1/|𝕋_xy^2|(||∂_zw||_L_xy^2||θ||_L_xy^2+||w||_L_xy^2||∂_zθ||_L_xy^2)

for all z. Therefore,

    ||wθ||_L^∞([0,z])   ≤1/|𝕋_xy^2|∫_0^z||∂_zw(z')||_L_xy^2||θ(z')||_L_xy^2+||w(z')||_L_xy^2||∂_zθ(z')||_L_xy^2 dz'
       ≤1/|𝕋_xy^2|(||∂_zw||_L^2(Ω)||θ||_L^2([0,z];L_xy^2)+||∂_zθ||_L^2(Ω)||w||_L^2([0,z];L_xy^2))
       ≤|z|^1/2/|𝕋_xy^2|(||∂_zw||_L^2(Ω)||θ||_L^∞([0,z];L_xy^2)+||∂_zθ||_L^2(Ω)||w||_L^∞([0,z];L_xy^2)).

Applying eq:rootgrowth_w and eq:rootgrowth_theta
we find that 

    ||wθ||_L^∞([0,z])≤4|z|/|𝕋_xy^2|||∂_zw||_L^2(Ω)||∂_zθ||_L^2(Ω).

The argument above is symmetric under z→1-z, so we immediately
obtain the inequality

    ||wθ||_L^∞([1-z,1])≤4|1-z|/|𝕋_xy^2|||∂_zw||_L^2(Ω)||∂_zθ||_L^2(Ω).

These two combine to prove the result.

Now we are ready to prove prop:ensbound_w2w, which immediately
implies the upper bound part of thm:mainbounds_w2w_ens.
We follow the plan laid out above.



We apply cor:symmsteadybound_w2w with an appropriate
class of test functions {η_δ}. Given δ∈(0,1/2],
we define η_δ by

    η_δ(z)=
    1-1/2δz     0≤ z≤δ
    1/2   δ≤ z≤1-δ
    1/2δ(1-z)     1-δ≤ z≤1
    .

Note these are admissible in eq:symmsteadybound_w2w.
Thus, 

    Nu(𝐮)≤inf_δ∈(0,1/2]{_Ω|∇η_δ|^2+⟨|∇Δ^-1div(𝐮η_δ)|^2⟩} .

The first integral appearing above is simple to estimate and it satisfies

    _Ω|∇η_δ|^2∼1/δ.

So,

    Nu≲inf_δ∈(0,1/2]{1/δ+⟨|∇Δ^-1div(𝐮η_δ)|^2⟩} .


Now we must estimate the non-local term. We claim that

    ⟨|∇Δ^-1div(𝐮η_δ)|^2⟩≲δ^2⟨ |∇ w|^2⟩

for all δ∈(0,1/2]. We argue at a.e. time t.
By duality, 

    ∫_Ω|∇Δ^-1div(𝐮η_δ)|^2=sup_θ∈ H_0^1(Ω) 2(div(𝐮η_δ),θ)-∫_Ω|∇θ|^2.

Since  η_δ only depends on z, div(𝐮η_δ)=wη_δ'
and we have by Fubini that

    _Ω|∇Δ^-1div(𝐮η_δ)|^2=sup_θ∈ H_0^1(Ω) _I_z2θ wη_δ'-_Ω|∇θ|^2.

Therefore to show eq:claim_for_nopenetration it suffices
to prove the inequality

    _I_z2θ wη_δ'≤ Cδ^2_Ω|∇ w|^2+_Ω|∇θ|^2  ∀ 𝐮∈ H_0^1(Ω;ℝ^d),θ∈ H_0^1(Ω)

Recalling the formula for η_δ from eq:testfunctions,
we see we must prove that

    _0^δ|θ w(z)| dz≲δ(_Ω|∇ w|^2)^1/2(_Ω|∇θ|^2)^1/2  ∀ 𝐮∈ H_0^1(Ω;ℝ^d),θ∈ H_0^1(Ω).

Applying lem:lingrowth_wtheta proves this result and
hence the desired estimate eq:claim_for_nopenetration. 

Assembling eq:intermediate_Nu_fromabove and eq:claim_for_nopenetration,
we conclude that 

    Nu≲inf_δ∈(0,1/2]{1/δ+δ^2⟨ |∇ w|^2⟩} .

If ⟨ |∇ w|^2⟩≥1 we may choose
δ∼⟨ |∇ w|^2⟩ ^-1/3 to conclude
that


    Nu≲⟨ |∇ w|^2⟩ ^1/3.

To handle the case ⟨ |∇ w|^2⟩≤1
we treat the choice δ=1/2 more carefully in the above.
Since η_1/2=1-z, eq:intermediate_Nu_fromabove_0
and eq:claim_for_nopenetration combine to prove that

    Nu≤1+⟨|∇Δ^-1div(𝐮η_1/2)|^2⟩≤1+C⟨ |∇ w|^2⟩ .

From eq:Nu-bd-1 and eq:Nu-bd-2 we conclude
the result.





§ OPTIMAL DESIGN OF STEADY WALL-TO-WALL TRANSPORT 


This section begins the proof of the lower bounds from thm:mainbounds_w2w_eng
and thm:mainbounds_w2w_ens. While every admissible velocity
field yields a lower bound on the maximal rate of heat transport,
it is not at all clear what sorts of features are required for velocity
fields to achieve maximal (or nearly maximal) transport. It is natural
to wonder how the overall character of optimal designs depends on
the intensity budget. One possibility is captured by the convection
rolls pictured in fig:rollsrough. This design is a relatively
simple one, as the number of lengthscales required to describe it
remains independent of the intensity budget. A second, much more complicated
possibility is captured by the branching designs from fig:branchingrough.
There, the total number of lengthscales is allowed to depend on the
intensity budget and can be unbounded as Pe→∞. In any case,
one requires a general method by which to evaluate Nu to allow
for comparison between candidate designs. 

The best scenario would be to develop an ansatz-free approach to evaluating
heat transport that, by its functional form, suggests optimal designs.
In this section, we achieve this for a general class of steady (i.e.,
time-independent) wall-to-wall optimal transport problems, including
the energy- and enstrophy-constrained ones as special cases. The class
of problems we have in mind are of the form

    sup_𝐮(𝐱)
    
    ||𝐮||=Pe
    
    +b.c. Nu(𝐮)

where ||·|| denotes any norm in which the advective intensity
of 𝐮 may be measured. As described in the introduction,
the steady energy-constrained problem arises from employing the (volume-averaged)
L^2-norm

    ||𝐮||=(_Ω|𝐮|^2)^1/2

to measure advective intensity, while the steady enstrophy-constrained
one arises from the (volume-averaged) Ḣ^1-norm

    ||𝐮||=(_Ω|∇𝐮|^2)^1/2

In any case, we require that 𝐮∈ L^2 in order that
its heat transport be well-defined. As our aim in this section is
to present a general approach to intensity-constrained optimal transport,
we leave the boundary conditions unspecified. Of course, we do not
claim that there exist optimizers at this level of generality.

The principal result of this section is that the general wall-to-wall
optimal transport problem eq:optimalw2w_generalintensity-1
can be reformulated as the double minimization 

    inf_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    
    +b.c. _Ω|∇Δ^-1div(𝐮ξ)|^2+1/Pe^2||𝐮||^2·_Ω|∇ξ|^2

in the velocity field 𝐮 and a new variable ξ. The
boundary conditions for 𝐮 remain the same as for eq:optimalw2w_generalintensity-1,
while ξ is required to vanish at ∂Ω. As will become
clear, ξ plays a role in the analysis of Nu similar to that
of η from the a priori bounds of sec:aprioribds —in
fact, these variables are dual. The optimal values in eq:optimalw2w_generalintensity-1
and eq:optimaldesgnpblm_general-1 are reciprocals, and
their optimizers are related through a certain change of variables.
We refer the reader forward to sec:energy-constraineddesign
and sec:enstrophy_constrained_design for the application
of these observations to energy- and enstrophy-constrained optimal
transport. Presently, our goal is to establish the connection between
eq:optimalw2w_generalintensity-1 and eq:optimaldesgnpblm_general-1,
and to illustrate how the latter suggests optimal designs. As in sec:aprioribds,
our approach centers on the existence of a variational principle for
Nu(𝐮); it is dual to the one appearing there. After achieving
this duality and using it to obtain eq:optimaldesgnpblm_general-1,
we proceed to make some general remarks on the construction of near
optimal designs. 


What distinguishes our approach from others in the literature (e.g.,
<cit.>) is that it does not rely on first order
optimality conditions or Euler-Lagrange equations to describe designs.
If one wants to discuss globally optimal designs this is an important
distinction to make, as the wall-to-wall optimal transport problem
is a non-concave variational problem and checking a first order optimality
condition is not guaranteed to produce a global optima. Instead, our
approach centers on the existence of a variational principle for the
heat transport functional Nu(𝐮), dual to the one described
in sec:aprioribds. Using this observation, we are able
to recast the steady wall-to-wall optimal transport problem as a double
optimization in 𝐮 and a new auxiliary variable ξ
(dual to the variable η from sec:aprioribds). The
solution of this latter problem is guaranteed to produce globally
optimal steady designs. 





 §.§ Dual variational formulations for transport


Recall from the analysis of a priori bounds on transport in
sec:aprioribds that there is a variational principle
for heat transport in the steady case, and that Nu can be written
as the optimal value of a certain convex minimization problem:

    Nu(𝐮)-1=min_η|_∂Ω=0 _Ω|∇η|^2+|∇Δ^-1div(𝐮η)|^2.

For the precise statement, see thm:steady-primal-varprin.
As this is convex it should, in principle, admit a dual formulation. 

 Let 𝐮 be a divergence-free
vector field in L^2(Ω;ℝ^3). Then,

    Nu(𝐮)-1=max_ξ∈ H_0^1(Ω)∩ L^∞(Ω) _Ω2wξ-|∇ξ|^2-|∇Δ^-1div(𝐮ξ)|^2.



Recall from the proof of thm:steady-primal-varprin the
PDE system

    𝐮·∇η=Δξ+w,
    𝐮·∇ξ=Δη

and the formula

    Nu-1=_Ω|∇ξ|^2+|∇η|^2.

Testing the first equation in eq:steadyPDEsystem-1 with
ξ, the second with η, and integrating by parts yields the
string of equalities

    ∫_Ωwξ    =∫_Ω|∇ξ|^2+∫_Ωξ𝐮·∇η
        =∫_Ω|∇ξ|^2-∫_Ωη𝐮·∇ξ=∫_Ω|∇ξ|^2+|∇η|^2.

Thus,

    Nu-1=_Ω2wξ-|∇ξ|^2-|∇η|^2=_Ω2wξ-|∇ξ|^2-|∇Δ^-1div(𝐮ξ)|^2.

Note in the last step we used the PDE system again.

Now consider the maximization

    sup_ξ∈ H_0^1(Ω)∩ L^∞(Ω) _Ω2wξ-|∇ξ|^2-|∇Δ^-1div(𝐮ξ)|^2.

Reasoning with its Euler-Lagrange equation just as in the proof of
thm:steady-primal-varprin, we deduce that this maximization
problem is well-posed in the given admissible class, with optimal
value Nu-1. 

As claimed in the introduction, there is a corresponding result holding
for unsteady velocities which allows to bound Nu from below, but
not necessarily to evaluate it. This result was described in eq:unsteady_lb,
and its proof is similar to that of thm:primal_Nubound_w2w.
We remark that although the variational formulas for steady heat transport
from thm:steady-primal-varprin and thm:dual_Nubound_w2w
are strongly dual —they provide convex and concave alternatives
for evaluating Nu whose optimal values agree —such strong
duality need not hold for unsteady flows. More precisely, we note
that for general velocity fields the bounds eq:unsteady_ub
and eq:unsteady_lb need not coincide. In particular,
there will be a duality gap for any velocity field which satisfies

    lim inf_τ→∞ ⟨ wT⟩ _τ<lim sup_τ→∞ ⟨ wT⟩ _τ.

These fields have the peculiar property that the lim sup and lim inf
alternatives for defining the space and long-time average heat transport
Nu do not coincide.



 §.§ An integral formulation of wall-to-wall optimal transport


Having written Nu for steady flows as the optimal value of the
concave maximization problem eq:dual_steady_w2w, we can
now give a useful reformulation of the entire class of steady wall-to-wall
optimal transport problems from eq:optimalw2w_generalintensity-1.
This “integral” formulation of steady optimal transport will be
used to design and evaluate nearly optimal flows in the subsequent. 

Let F(Pe) denote the optimal value of the steady optimal transport
problem eq:optimalw2w_generalintensity-1,

    F(Pe)=sup_𝐮(𝐱)
    
    ||𝐮||=Pe
    
    +b.c. Nu(𝐮).

 Applying thm:dual_Nubound_w2w, we find that

    F(Pe)=sup_𝐮(𝐱),ξ(𝐱)
    
    ||𝐮||=Pe
    
    +b.c. _Ω2wξ-|∇ξ|^2-|∇Δ^-1div(𝐮ξ)|^2.

The boundary conditions for 𝐮 remain unchanged, while
according to thm:dual_Nubound_w2w we must require that
ξ|_∂Ω=0. Performing the substitution

    ξ→λξ,  λ∈ℝ

and maximizing over λ yields the equivalent variational problem

    1/F(Pe)=inf_𝐮(𝐱),ξ(𝐱)
    
    ||𝐮||=Pe
    
    +b.c. _Ω|∇Δ^-1div(𝐮ξ)|^2+|∇ξ|^2/(_Ωwξ)^2.

Changing variables via the substitutions

    𝐮→ Pe𝐮/||𝐮||  and  ξ→||𝐮||/Peξ

allows to eliminate the intensity constraint on 𝐮 altogether,
so that

    1/F(Pe)=inf_𝐮(𝐱),ξ(𝐱)
    
    +b.c. _Ω|∇Δ^-1div(𝐮ξ)|^2+1/Pe^2||𝐮||^2_Ω|∇ξ|^2/(_Ωwξ)^2.

Given the scaling symmetries of the above, we may impose the constraint

    _Ωwξ=1

on the minimization without altering the result. This yields the promised
integral reformulation of wall-to-wall optimal transport


    1/F(Pe)=inf_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    
    +b.c. _Ω|∇Δ^-1div(𝐮ξ)|^2+1/Pe^2||𝐮||^2·_Ω|∇ξ|^2,

and proves the equivalence between eq:optimalw2w_generalintensity-1
and eq:optimaldesgnpblm_general-1. 

We turn to discuss the integral formulation of optimal transport just derived.
Observe eq:optimaldesgnpblm_general consists of two types
of terms, each of which prefers a different kind of design. The first
term

    _Ω|∇Δ^-1div(𝐮ξ)|^2

prefers 𝐮ξ to be divergence-free and we refer to it
as the “advection term” throughout. This preference is strong
in the advective regime Pe≫1, as it appears at leading order
in Pe^-1 in the functional above. The remaining terms

    ||𝐮||^2  and  _Ω|∇ξ|^2

contribute at higher order in Pe^-1, and act to regularize designs.
Any admissible design must satisfy the “net-flux” constraint

    _Ωwξ=1

as well as boundary conditions. While patterns such as the convection
roll and branching ones depicted in fig:designs can be
easily made to satisfy such constraints —see sec:energy-constraineddesign
and sec:enstrophy_constrained_design for details —determining
the optimal lengthscales for such designs requires performing an optimization
as in eq:optimaldesgnpblm_general.

Evidently, the most difficult term to evaluate is the advection one
eq:advectionterm. Before turning to discuss its analysis
in detail, and what it implies for near optimal designs, we make two
general remarks. In order to get a hint as to what designs eq:optimaldesgnpblm_general
prefers in the advective limit Pe→∞, one can entertain
the “limiting” wall-to-wall optimal transport problem

    inf_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    
    +b.c. _Ω|∇Δ^-1div(𝐮ξ)|^2

obtained by formally taking Pe = ∞ in (<ref>). This, however, is an ill-posed variational problem. Its optimal value
is zero as there exist admissible sequences {(𝐮_k,ξ_k)}
satisfying the net-flux constraint eq:netflux and achieving

    lim_k→∞ _Ω|∇Δ^-1div(𝐮_kξ_k)|^2=0.

Yet, no smooth enough admissible pair (𝐮,ξ) can satisfy
the net-flux constraint and simultaneously achieve 

    div(𝐮ξ)=0.

Indeed, if 𝐮ξ were divergence-free then by averaging
eq:div-free in the periodic variables x and y we
would find that the flux of ξ by 𝐮 through each slice
{z=const.} is independent of the slice, i.e.,

    wξ(z)=_𝕋_x,y^2wξ

is constant in z. Applying the boundary conditions which require
at least that ξ|_∂Ω=0, we conclude that wξ
must vanish throughout the entire domain. This contradicts the net-flux
constraint eq:netflux. Therefore, wall-to-wall optimal
transport is a singularly perturbed variational problem: the regularizing
terms from eq:optimaldesgnpblm_general which at first
glance appear to contribute at higher order in Pe^-1 are crucial
for determining the character of optimal designs. 

Our second observation is more straightforward: it is regarding the
disappearance of the intensity constraint in the passage from eq:optimalw2w_generalintensity-1
to its integral formulation eq:optimaldesgnpblm_general.
Since eq:optimaldesgnpblm_general is invariant under
the rescaling

    𝐮→λ𝐮  and  ξ→1/λξ,  λ≠0,

the magnitudes of any of its minimizers are not uniquely determined.
Still, if (𝐮_Pe,ξ_Pe) achieves optimality in eq:optimaldesgnpblm_general,

    𝐮=Pe/||𝐮_Pe||𝐮_Pe

solves the wall-to-wall problem eq:optimalw2w_generalintensity-1. 



 §.§ Analysis of the advection term


For a class of designs {(𝐮_α,ξ_α)}_α∈ℐ
to compete in the minimization eq:optimaldesgnpblm_general,
it must at least achieve

    inf_α∈ℐ _Ω|∇Δ^-1div(𝐮_αξ_α)|^2=0.

How difficult is it for an admissible pair (𝐮,ξ) to
make this advection term nearly zero? First, note that in such
a situation, the vertical flux of ξ by 𝐮 through
each slice {z=const.} must be nearly independent of the
slice, 

    d/dzwξ≈0.

By the net-flux constraint eq:netflux, it follows that

    wξ≈1

in nearly all of the domain. This is an example of a “design principle”
for wall-to-wall optimal transport: any nearly optimal design must
achieve eq:design principle with equality in the limit
Pe→∞. Although eq:design principle does not
completely characterize optimal designs, it does give a necessary
condition for constructing competitive ones. This will be particularly
useful later on in sec:enstrophy_constrained_design,
where we devise a functional form for the branching depicted in fig:branchingrough. 

The advection term eq:advectionterm contains a wealth
of information for evaluating designs beyond eq:design principle,
but to use it in practice one must confront its non-locality. In
the wall-to-wall domain Ω=𝕋_xy^2×[0,1]_z≅[0,l_x]×[0,l_y]×[0,1],
we may take advantage of periodicity to represent a function f∈ L^2(Ω)
via its Fourier series 

    f(𝐱)=∑_𝐤∈ℤ_l_x,l_y^2f̂_𝐤(z)e^i𝐤·𝐱'

where 𝐱'=(x,y) and 

    ℤ_l_x,l_y^2={ (k_x,k_y):l_x/2πk_x,l_y/2πk_y∈ℤ} .

We employ the Fourier transform 

    f̂_𝐤(z)=_𝕋_xy^2e^i𝐤·𝐱'f(𝐱',z) dxdy

and proceed to decompose the advection term mode-by-mode.

 Let 𝐮∈ L^2(Ω;ℝ^3)
and ξ∈ L^∞(Ω). The advection term satisfies 

    _Ω|∇Δ^-1div(𝐮ξ)|^2=_I_z|wξ-_Ωwξ|^2+𝒬(div(𝐮ξ))

where 𝒬=∑_𝐤≠0Q_𝐤
and 𝒬_𝐤 is the positive semi-definite quadratic
form given by

    𝒬_𝐤(f)=_I_z× I_z'G_𝐤(z,z')f̂_𝐤(z)f̂_𝐤^*(z') dzdz'

with kernel

    G_𝐤(z,z')=csch(|𝐤|)/|𝐤|×sinh(|𝐤|z)sinh(|𝐤|(1-z'))     z≤ z'
    sinh(|𝐤|z')sinh(|𝐤|(1-z))     z≥ z'
        𝐤≠0.



This follows from the Green's function representation for -Δ^-1
on Ω with vanishing Dirichlet boundary conditions. Calling
J=div(𝐮ξ) and applying Parseval's identity, we
see that the advection term can be written as

    _Ω|∇Δ^-1div(𝐮ξ)|^2    =1/|Ω|(J,-Δ^-1J)=∑_𝐤_I_z(-d^2/dz^2+|𝐤|^2)^-1Ĵ_𝐤Ĵ_𝐤^*
        =∑_𝐤_I_z× I_z'G_𝐤(z,z')Ĵ_𝐤(z)Ĵ_𝐤^*(z') dzdz'

where for 𝐤≠0 the functions G_𝐤
are as defined above. For 𝐤=0, we have 

    G_0(z,z')=
    z(1-z')     z≤ z'
    
    z'(1-z)     z≥ z'
    .

We must only address the form of the 0th term now.

We recognize that

    _I_z× I_z'G_0(z,z')Ĵ_0(z)Ĵ_0(z') dzdz'=_Ω|∇Δ^-1J(z)|^2.

By periodicity,

    J=div(𝐮ξ)=d/dzwξ.

Thus,

    _Ω|∇Δ^-1J(z)|^2    =_Ω|d/dz(d/dz)^-2d/dzwξ(z)|^2
        =_I_z|wξ-_I_zwξ|^2=_I_z|wξ-_Ωwξ|^2

as required.

As the quadratic form 𝒬 from lem:advection_Fourierests_w2w
is non-negative, the advection term satisfies the lower bound

    _Ω|∇Δ^-1div(𝐮ξ)|^2≥_I_z|wξ-_Ωwξ|^2.

Note this quantifies the design principle eq:design principle.
The appearance of 𝒬 in lem:advection_Fourierests_w2w
also makes clear why this principle alone does not suffice to characterize
optimal designs.

We end this section by recording some useful estimates on the kernels
{G_𝐤} from the definition of 𝒬. These
will be used later in sec:enstrophy_constrained_design.

 Let A⊂ I_z be Borel measurable.
Then, 

    ||G_𝐤||_L^1(A× A)≲|A|/|𝐤|(|A|∧1/|𝐤|)  𝐤≠0

and 

    ||G_0||_L^1(A× A)≲|A|^2_Az∧(1-z) dz.



To see the first estimate, observe that G_𝐤 satisfies
the pointwise estimate

    |G_𝐤(z,z')|≤1/|𝐤|e^-|𝐤||z-z'|  𝐤≠0.

Now,

    ∫_I_z× I_ze^-|𝐤||z-z'|A(z)A(z') dzdz'≤∫_I_z[∫_e^-|𝐤||z-z'|dz']A(z) dz=2/|𝐤||A|.

Therefore,

    ||G_𝐤||_L^1(A× A)≲|A|/|𝐤|^2.

On the other hand, we have that

    ||G_𝐤||_L^1(A× A)≤||G_𝐤||_L^∞||1_A× A||_L^1≤|A|^2/|𝐤|.

Combining these two bounds gives the first result.

Now we prove the second estimate. We need to show that

    ∫ G_0(z,z')A(z)A(z') dzdz'≲|A|∫_Az∧(1-z) dz.

By symmetry, 

    ∫ G_0(z,z'))A(z)A(z') dzdz'     =2∫_0^1∫_0^z'z(1-z')1_A(z)1_A(z') dzdz'
       ≤2∫_0^1∫_0^z'z'(1-z')1_A(z)1_A(z') dzdz' 
       ≤2|A|∫_0^1z(1-z) dz.

Since 

    z(1-z)≤ z∧(1-z)  ∀ z∈[0,1]

we conclude the desired result.




§ ENERGY-CONSTRAINED TRANSPORT AND CONVECTION ROLL DESIGNS


In the previous section, we considered the general class of steady
wall-to-wall optimal transport problems eq:optimalw2w_generalintensity-1
and produced their equivalent integral formulations eq:optimaldesgnpblm_general-1.
In this section and the next, we use these formulations to study the
steady energy- and enstrophy-constrained problems. The subsequent
analyses are largely independent. Nevertheless, the reader may find
it helpful to study the energy-constrained problem first as its proof
is much shorter and its technical details much less burdensome.

Here we discuss energy-constrained transport. The main result of this
section is a proof of the lower bound from thm:mainbounds_w2w_eng.
Recall from sec:OptlDesignProb_w2w that the steady energy-constrained
optimal transport problem

    max_𝐮(𝐱)
    _Ω|𝐮|^2=Pe^2
    
    w|_∂Ω=0
     Nu(𝐮)

admits the integral formulation 

    min_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    
    w|_∂Ω=ξ|_∂Ω=0
     _Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ_Ω|𝐮|^2·_Ω|∇ξ|^2

with ϵ=Pe^-2. The optimal values of eq:steady_w2w_eng
and eq:optimaldesign_w2weng are reciprocals and their
optimizers are related through symmetrization. 

Our goal now is to identify the scaling law of eq:optimaldesign_w2weng
in the advective regime ϵ≪1. Combined with the results
of sec:OptlDesignProb_w2w, this completes the proof of
the lower bound half of thm:mainbounds_w2w_eng. 

 Let 𝔈(ϵ)
denote the optimal value of eq:optimaldesign_w2weng.
Then, 

    𝔈(ϵ)∼ϵ^1/2

whenever ϵ≲1∧ l_x^4∧ l_y^4. 

The a priori lower bound 𝔈(ϵ)≳ϵ^1/2
is implied by the corresponding bound Nu≲ Pe from thm:mainbounds_w2w_eng.
The remainder of this section is regarding the upper bound 𝔈(ϵ)≲ϵ^1/2.
We prove it by constructing the convection roll designs depicted in
fig:rollsrough. Such designs can be parameterized using
two variables: the number of rolls and their wall-to-wall extent.
Carrying out the optimization from eq:optimaldesign_w2weng
with respect to these variables yields the desired upper bound. The
condition that ϵ be small enough in the statement above
is required to fit what would be, in the absence of an overall horizontal
period, an optimal number of rolls inside the domain. Given the symmetry
between x and y, we may suppose that l_x≤ l_y in what
follows. 



 §.§ Convection roll designs


The integral formulation eq:optimaldesign_w2weng requires
designing a velocity field 𝐮 and a test function ξ.
For the velocity field, we introduce a family of streamfunctions of
the form

    ψ(x,z)=χ(z)Ψ(x)  χ∈ C_c^∞(I_z), Ψ∈ C^∞(𝕋_x).

Each such ψ gives rise to a divergence-free velocity field by

    𝐮=(-∂_zψ,0,∂_xψ).

These are two-dimensional flows as their ĵ-component vanishes
identically. Although we do not claim that optimizers must be of this
form, we will prove that such a construction suffices to capture the
optimal scaling law of eq:optimaldesign_w2weng.

Next we must describe test functions ξ well-suited to the velocity
fields. Recall the design principle eq:design principle,
which states that for a design to be competitive it must satisfy 

    wξ≈_Ωwξ=1.

This rules out taking, for instance, ξ=ψ as it would result
in zero flux thru each slice {z=const.}. We can, however,
choose ξ to depend only on x and z as does ψ. Then
by Parseval's identity we can rewrite the flux as 

    wξ=∑_𝐤ŵ_𝐤(z)ξ̂_𝐤^*(z)=∑_𝐤(ik_x^1/2ψ̂_𝐤)(k_x^1/2ξ̂_𝐤)^*.

Taking

    |∂_x^1/2ξ|∼|∂_x^1/2ψ|

allows to satisfy eq:designprinciple-recalled.

Now we make the convection roll construction concrete. Given δ∈(0,1/2)
and l such that 

    1/l∈2π/l_xℕ,

we define

    ψ_δ,l(x,z)     =χ_δ(z)· l^1/2Ψ(x/l)
    ξ_δ,l(x,z)     =χ_δ(z)· l^1/2Ψ'(x/l)

where 

    Ψ(x)=c_0cos x

and c_0 is chosen so that (Ψ')^2=1. Here,
the cut-off functions {χ_δ}⊂ C_c^∞(I_z)
are required to satisfy 


  * χ_δ=1 on [δ,1-δ],

  * |χ_δ|≲1 and |χ'_δ|≲1/δ,

  * _I_zχ_δ^2=1.

The constants in these assumptions are independent of all parameters.
In what follows, we often neglect to record the subscripts δ
and l as the meaning is clear.

First, we check admissibility.

The convection roll construction described above is admissible for
eq:optimaldesign_w2weng. 


All conditions in admissibility are clear, except for the net-flux
constraint which we verify now. Given the above, we find that

    wξ=∂_xψξ=χ_δ^2(Ψ')^2=χ_δ^2

for all z, so that

    _Ωwξ=_I_zwξ=_I_zχ_δ^2=1

as required.

Next, we estimate the advection term from eq:optimaldesign_w2weng.

The convection roll construction satisfies 

    _Ω|∇Δ^-1div(𝐮ξ)|^2=∫_I_z|wξ-1|^2≲δ.

In particular, the quadratic form 𝒬 from lem:advection_Fourierests_w2w
vanishes on it.


We apply lem:advection_Fourierests_w2w. Note that since
all functions entering into the construction are independent of y,
we may work with k in place of 𝐤=(k_x,0). We start
with the k=0 term from lem:advection_Fourierests_w2w:
it satisfies the estimate

    ∫_I_z|wξ-1|^2=∫_I_z|χ^2-1|^2≲δ.

Now we address the k≠0 terms. We must compute the quadratic
form 𝒬 from lem:advection_Fourierests_w2w,
and to do so we must compute Ĵ_k for k≠0 where J=𝐮·∇ξ.
Note that by the form of the convection roll construction,

    J     =∇^⊥ψ·∇ξ=-∂_zψ∂_xξ+∂_xψ∂_zξ
        =-(χ_δ'(z)· l^1/2Ψ(x/l))(χ_δ(z)·1/l^1/2Ψ”(x/l)) 
            +(χ_δ(z)1/l^1/2Ψ'(x/l))(χ_δ'(z)· l^1/2Ψ'(x/l))
        =1/2(χ_δ^2)'(z)·Θ(x/l)

where Θ=(Ψ')^2-ΨΨ”. Since Θ=c_0^2sin^2+c_0^2cos^2=c_0^2,

    J=1/2(χ_δ^2)'(z)

which is entirely a function of z. This shows that Ĵ_k=0
for k≠0. Hence, 𝒬 vanishes on the convection roll
construction. 

Continuing, we estimate the higher order terms from eq:optimaldesign_w2weng. 

The convection roll construction satisfies 

    _Ω|∇ξ|^2∨_Ω|𝐮|^2≲l/δ+1/l.



Clearly, 

    _Ω|∇ξ|^2∼_Ω|∇ψ|^2

since ||Ψ^(k)||_L^2(𝕋_x)∼1 for all k.
Noting that

    _Ω|𝐮|^2=_Ω|∇ψ|^2=_Ω(χ'_δ)^2lΨ^2+χ_δ1/l(Ψ')^2≲l/δ+1/l

we conclude the result.


Combining the above yields the following estimate on 

    E(ϵ;δ,l)=_Ω|∇Δ^-1div(𝐮_δ,lξ_δ,l)|^2+ϵ_Ω|𝐮_δ,l|^2·_Ω|∇ξ_δ,l|^2.



The convection roll construction satisfies

    E(ϵ;δ,l)≲δ+ϵ(l/δ+1/l)^2.





 §.§ Proof of prop:lowerbound_engconstr_w2w


Finally we choose δ and l to prove the desired bound.


cor:convection_roll_est holds for all admissible δ
and l, i.e., so long as δ∈(0,1/2) and l^-1∈2π l_x^-1ℕ.
To optimize the result, we take 

    l∼√(δ)  and  δ∼ϵ^1/2

which we can do so long as ϵ≲1∧ l_x^4. Under
this condition, 

    inf_δ,l E(ϵ;δ,l)≲ϵ^1/2.






§ ENSTROPHY-CONSTRAINED TRANSPORT AND BRANCHED FLOW DESIGNS 


We now consider the steady enstrophy-constrained wall-to-wall
optimal transport problem in the framework of sec:OptlDesignProb_w2w.
The main result of this section is a proof of the lower bound from
thm:mainbounds_w2w_ens. As in the previous section on
energy-constrained transport, we exploit the fact that
the enstrophy-constrained problem

    max_𝐮(𝐱)
    _Ω|∇𝐮|^2=Pe^2
    𝐮|_∂Ω=0 Nu(𝐮)

can be written in integral form as

    min_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    𝐮|_∂Ω=0,ξ|_∂Ω=0
     _Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ_Ω|∇𝐮|^2·_Ω|∇ξ|^2

where ϵ=Pe^-2. This form of the problem suggests the possibility
of analyzing a certain multiple scales ansatz for 𝐮 and
ξ (we explain the intuition behind this further in sec:Implications-for-RBC
and sec:energy-driven-pattern-formation). As proved below,
such an ansatz turns out to capture the scaling of the optimal value
of eq:optimaldesign_w2wenstr in ϵ up to possible
logarithmic corrections. The precise statement is as follows:

 Let 𝔈(ϵ)
denote the optimal value of eq:optimaldesign_w2wenstr.
Then, 

    ϵ^1/3≲𝔈(ϵ)≲ϵ^1/3log^4/31/ϵ

whenever ϵ≲1 and ϵlog1/ϵ≲ l_x^6∧ l_y^6.

The a priori lower bound 𝔈(ϵ)≳ϵ^1/3
is implied by the upper bound Nu≲ Pe^2/3 from thm:mainbounds_w2w_ens.
(For a proof which is more self-contained, see the discussion surrounding
eq:Howard_LB.) To prove the upper bound 𝔈(ϵ)≲ϵ^1/3log^4/31/ϵ
we must construct a suitable class of designs and estimate their heat
transport. The successful ones are as depicted in fig:branchingrough.
In contrast with the convection roll designs considered previously,
such “branching” designs are evidently more complicated to analyze.
The main challenge of course lies with estimating the advection term.
Here, we make use of lem:advection_Fourierests_w2w and
lem:Gestimates. Note the requirement that ϵ
be small enough is to ensure that our construction fits into the given
domain. As in the previous section, we need only consider the case
l_x≤ l_y by symmetry. 

Combined with the results of sec:OptlDesignProb_w2w,
prop:lowerbound_enstrconstr_w2w completes the proof of
the lower bound from thm:mainbounds_w2w_eng. 


Let us state a refinement of the result advertised above, that we
prove in this section.

There exist steady 2D divergence-free velocity vector fields {𝐮_Pe}
satisfying no-slip boundary conditions and such that

    _Ω|∇𝐮_Pe|^2=Pe^2  and   Nu(𝐮_Pe)∼Pe^2/3/log^4/3Pe.







 §.§ The branching construction


The integral formulation eq:optimaldesign_w2wenstr requires
the construction of a divergence-free velocity field 𝐮
and a test function ξ (the latter of which plays the role of
temperature in this approach). For the velocity fields, we will use
a streamfunction ψ(x,z) whose streamlines are as in fig:branchingrough.
That figure can be thought of as consisting of many individual convection
roll systems which have been carefully fit together. In the bulk,
there are large anisotropic convection rolls at some horizontal lengthscale
l_bulk. At the walls, there are much smaller isotropic convection
rolls at some other lengthscale l_bl≪ l_bulk. Between the
bulk and the walls, streamlines branch and refine through several
transition layers, a single one of which is shown in fig:transitionlayer.
As the construction is symmetric about z=1/2, we only need describe
it for z∈[1/2,1].

Counting upwards from the bulk, we understand by the jth transition
layer that part of the domain where z∈[z_j,z_j+1]. The points
{z_j}_j=1^n marking the edges of the layers satisfy

    1/2<z_bulk=z_1<z_2<⋯<z_n=z_bl<1.

At the horizontal slice {z=z_j} the velocity components fluctuate
at lengthscale l_j. These decrease monotonically according as

    l_bulk=l_1>l_2>⋯>l_n=l_bl.

In what follows, we think of the parameters {z_j}_j=1^n
and {l_j}_j=1^n as playing a distinguished role in specifying
the branching design.

Given such a streamfunction ψ and its corresponding two-dimensional
velocity field

    𝐮=(-∂_zψ,0,∂_xψ),

we must choose a “temperature” field ξ well-suited to the
minimization eq:optimaldesign_w2wenstr. Recall the design
principle eq:design principle discussed in sec:OptlDesignProb_w2w,
which requires that 

    wξ≈_Ωwξ=1

throughout the domain. For our purposes, it will suffice to set

    ξ=w

and enforce that w^2≈1. Such considerations
significantly constrain the way that streamlines may branch. We note
that while eq:pointwise_slaving may not necessarily hold
for optimal designs, it greatly simplifies the ensuing analysis. And,
as claimed in prop:lowerbound_enstrconstr_w2w and proved
below, such a choice introduces at most a logarithmic error in our
estimates of enstrophy-constrained optimal transport.











We are now ready to give the precise functional form of our branching
construction. Let points {z_j}_j=1^n satisfying eq:zk_inequalities
and lengthscales {l_j}_j=1^n satisfying 

    1/l_j∈2π/l_xℕ

and eq:lk_inequalities be given. Let

    Ψ(x)=c_0cos x

where c_0 is chosen so that (Ψ')^2=1. We define
{ψ_j}_j=1^n by

    ψ_j(x)=l_jΨ(x/l_j),

and set

    ψ(x,z)=ψ_1(x)  for z∈[1/2,z_1]

which corresponds to the bulk. In the boundary layer, we set

    ψ(x,z)=g(z-z_n/1-z_n)ψ_n(x)  for z∈[z_n,1].

Here, g∈ C^∞([0,1]) is a fixed cutoff function satisfying
the matching conditions g(0)=1, g(1)=0, g'(0)=g'(1)=0, as
well as the integral condition

    ∫_0^1(g(t))^2 dt=1.

In the jth transition layer we take

    ψ(x,z)=f(z-z_j/z_j+1-z_j)ψ_j(x)+f(z_j+1-z/z_j+1-z_j)ψ_j+1(x)  for z∈[z_j,z_j+1]

where f∈ C^∞([0,1]) is a second fixed cutoff function.
We require it to satisfy the Pythagorean condition

    (f(t))^2+(f(1-t))^2=1  for t∈[0,1]

as well as the matching conditions f(0)=1, f(1)=0, and f'(0)=f'(1)=0.
One might choose, for instance, 

    f(t)=√(1/2-1/2tanh(τ·t-1/2/t^2(1-t)^2))  for 0<t<1

for some τ∈(0,∞). fig:cutofffunction shows
such a cutoff function. This completes our construction of a general
class of branching designs; we turn now to select an optimal one for
use in eq:optimaldesign_w2wenstr.













 §.§ Admissibility 


Our first task is to check the admissibility of this construction.
The Fourier series of ψ can be written in the form 

    ψ=∑_j=1^nχ_j(z)ψ_j(x),

where the amplitude functions {χ_j}_j=1^n belong to
H^2(I_z)∩ L^∞(I_z). (Such regularity is guaranteed
by the definition of the cutoff functions f and g.) The support
of the jth amplitude function satisfies

    supp χ_j⊂[z_j-1,z_j+1]∪[1-z_j+1,1-z_j-1]

and the expansion is nearly diagonal in the sense that 

    χ_jχ_j'≢0|j-j'|≤1.

Thus,

    supp χ_jχ_j+1⊂[z_j,z_j+1]∪[1-z_j+1,1-z_j]  ∀ j

and the same is true for products of derivatives thereof. Note also
that by eq:pythagorean_conditions,

    ∑_j=1^nχ_j^2=1   z∈[1-z_bl,z_bl]

and by eq:squareintegralisone_for_g,

    _z_bl^1χ_n^2 dz=1.

 

The branching construction defined above is admissible for eq:optimaldesign_w2wenstr. 


That the boundary conditions for 𝐮 and ξ are met
is clear. Here, we check that the net-flux constraint _Ωwξ=1
is satisfied. Since

    ∂_xψ=∑_j=1^nχ_jψ_j'

and {ψ_j'}_j=1^n forms an L^2-orthonormal set,
we can write that

    _Ωwξ=_Ω(∂_xψ)^2=∫_I_z∑_j=1^nχ_j^2 dz.

By eq:chisquaresum, the integrand is equal to one for
z∈[1-z_bl,z_bl]. Since 

    ∫_z_bl^1χ_n^2 dz=1-z_bl

by eq:blsquareintegral_is_one, we conclude that

    _Ωwξ=2(∫_1/2^z_bl+∫_z_bl^1)∑_j=1^nχ_j^2 dz=2(z_bl-1/2+1-z_bl)=1.


Next, we record some technical requirements that will greatly simplify
the identification of an optimal branching construction. These requirements
are compatible with the upper bound from prop:lowerbound_enstrconstr_w2w,
and we have not been able to improve upon the scaling of this result
by removing them. First, we require the transition layer thicknesses
{δ_k}_k=1^n, which are defined by

    δ_k=z_k+1-z_k  1≤ k≤ n-1  and  δ_bl=δ_n=1-z_n,

and the horizontal lengthscales {l_k}_k=1^n to satisfy
the relations

    δ_1≳δ_2≳⋯≳δ_n

and

    l_k≲δ_k  1≤ k≤ n-1    and    δ_n∼ l_n.

The latter guarantees that a certain anisotropy is present throughout
the construction which will simplify, amongst other things, the estimation
of the higher order terms from eq:optimaldesign_w2wenstr.
Second, we require that 

    δ_bulk=z_bulk-1/2∼1.

Note the constants implicit in eq:assumption2-0-eq:assumption3
are not allowed to depend on any parameters. Third, we require that
the refinement of lengthscale through each transition layer occur
by period doubling, i.e., 

    l_k+1=1/2l_k  1≤ k≤ n-1.

This last requirement will serve to simplify the Fourier analysis
involved in estimating the non-local advection term.



 §.§ Estimating the efficiency of branching


In this section we estimate each of the terms from eq:optimaldesign_w2wenstr
for the branching construction. The requirements laid out in sec:branchingconstruction
and sec:Admissibility are understood to hold. The constants
appearing below are only allowed to depend on those implicit in eq:assumption2-0-eq:assumption3,
and so do not depend on any parameters.

First, we deal with the advection term. By lem:advection_Fourierests_w2w,

    _Ω|∇Δ^-1div(𝐮ξ)|^2=_I_z|wξ-1|^2 dz+𝒬(𝐮·∇ξ)

where 𝒬=∑_𝐤≠0𝒬_𝐤
and 

    𝒬_𝐤(f)=_I_z× I_z'G_𝐤(z,z')f̂_𝐤(z)f̂_𝐤^*(z') dzdz'.

As the construction is two-dimensional, k_y does not play a role.
For ease of reading, we denote 𝐤=(k_x,0) simply by k
in what follows. 

 The branching construction satisfies

    𝒬(𝐮·∇ξ)=∑_i=1^n-1𝒬_k_i^sum(𝐮·∇ξ)+𝒬_k_i^diff(𝐮·∇ξ)

where

    k_i^sum=1/l_i+1+1/l_i  and   k_i^diff=1/l_i+1-1/l_i  for 1≤ i≤ n-1.



Using the formula for the branching construction given in eq:branchingstreamfunction_w2w,

    J     =∇^⊥ψ·∇∂_xψ=∑_j,j'∇^⊥(χ_iψ_i)·∇∂_x(χ_jψ_j)=∑_|j-j'|≤1∇^⊥(χ_jψ_j)·∇∂_x(χ_j'ψ_j')
        =∑_j=1^n∇^⊥(χ_jψ_j)·∇∂_x(χ_jψ_j) 
            +∑_j=1^n-1∇^⊥(χ_jψ_j)·∇∂_x(χ_j+1ψ_j+1)+∇^⊥(χ_j+1ψ_j+1)·∇∂_x(χ_jψ_j)
        =J_self+J_nbr.

Given our choice of fundamental streamfunction eq:masterstreamfunction,
these expressions can be made explicit and we do so now.

The general term in the first sum, J_self, satisfies 

    ∇^⊥(χ_jψ_j)·∇∂_x(χ_jψ_j)     =(-χ_j'ψ_j,χ_jψ_j')·(χ_jψ_j”,χ_j'ψ_j')
        =-χ_j'χ_jψ_jψ_j”+χ_jχ_j'ψ_j'ψ_j'=1/2(χ_j^2)'(z)·Θ(x/l_j)

where Θ=(Ψ')^2-ΨΨ”. Using eq:masterstreamfunction,
we see that Θ=c_0^2(cos^2+sin^2)=c_0^2
so that

    J_self=c_0^2∑_i=1^n(χ_i^2)'.

In particular, we find that J_self is constant in the periodic
variable x, so that the Fourier coefficient [J_self]_k^∧
vanishes identically except for when k=0. For k=0, note that

    [J_self]_0^∧=J_self=J=1/2d/dz(∂_xψ)^2.


Continuing, we see that the general term in the second sum, J_nbr,
satisfies

    ∇^⊥(χ_jψ_j)·∇∂_x(χ_j+1ψ_j+1)+∇^⊥(χ_j+1ψ_j+1)·∇∂_x(χ_jψ_j)
        =-(χ_j'ψ_j,χ_jψ_j')·(χ_j+1ψ_j+1”,χ_j+1'ψ_j+1')-(χ_j+1'ψ_j+1,χ_j+1ψ_j+1')·(χ_jψ_j”,χ_j'ψ_j')
        =-χ_j'ψ_jχ_j+1ψ_j+1”+χ_jψ_j'χ_j+1'ψ_j+1'-χ_j+1'ψ_j+1χ_jψ_j”+χ_j+1ψ_j+1'χ_j'ψ_j'
        =χ_jχ_j+1'Θ_j,j+1+χ_j'χ_j+1Θ_j+1,j

where

    Θ_j,j+1    =ψ_j'ψ_j+1'-ψ_j”ψ_j+1
    Θ_j+1,j    =ψ_j+1'ψ_j'-ψ_j+1”ψ_j.

Given eq:masterstreamfunction, we find that

    Θ_j,j+1    =sin(x/l_j)sin(x/l_j+1)+l_j+1/l_jcos(x/l_j)cos(x/l_j+1)
    Θ_j+1,j    =sin(x/l_j)sin(x/l_j+1)+l_j/l_j+1cos(x/l_j)cos(x/l_j+1).

Applying standard trigonometric identities, 

    Θ_j,j+1    =1/2(1+l_j+1/l_j)cos(k_j^diffx)+1/2(l_j+1/l_j-1)cos(k_j^sumx)
    Θ_j+1,j    =1/2(1+l_j/l_j+1)cos(k_j^diffx)+1/2(l_j/l_j+1-1)cos(k_j^sumx)

where k_j^diff and k_j^sum are as in the statement of
the result. In sum, 

    J_nbr=∑_j=1^nχ_jχ_j+1'Θ_j,j+1+∑_j=1^nχ_j'χ_j+1Θ_j+1,j.


From eq:self-interactions, eq:trans-interactions,
and the decomposition J=J_self+J_nbr it is clear which wave-numbers
are present in 𝒬. We see that Ĵ_k is not identically
zero if and only if 

    k∈{0}∪{k_j^diff:1≤ j≤ n-1}∪{k_j^sum:1≤ j≤ n-1}.

Since l_j+1≠ l_j, we see that 

    0∉{k_j^diff:1≤ j≤ n-1}∪{k_j^sum:1≤ j≤ n-1}.

For general choices of lengthscales {l_j} these two sets of
wavenumbers may intersect; however, given our special choices of lengthscales
in eq:l_relations, we find that

    {k_j^diff:1≤ j≤ n-1}∩{k_j^sum:1≤ j≤ n-1}=∅.

(If k_i^sum=k_j^diff then 2^i+1+2^i=2^j+1-2^j
from which the contradiction 3·2^i=2^j follows.) Therefore,

    𝒬(𝐮·∇ξ)=∑_k≠0𝒬_k(𝐮·∇ξ)=∑_j=1^n-1𝒬_k_j^sum(𝐮·∇ξ)+𝒬_k_j^diff(𝐮·∇ξ).


Now we estimate each of the non-zero contributions to the advection
term picked out by lem:diagonal_interactions.

The branching construction satisfies 

    _I_z|wξ-1|^2 dz≲δ_n

and

    𝒬_k_j^sum(𝐮·∇ξ)∨𝒬_k_j^diff(𝐮·∇ξ)≲l_j^2/δ_j  1≤ j≤ n-1.



We begin with the k=0 term. Since 

    wξ=∑_j=1^nχ_j^2

we find that

    ∫_I_z|wξ-1|^2 dz≲δ_n.


Next we wish to estimate 𝒬_k_j^sum and 𝒬_k_j^diff.
Recall that

    𝒬_k_j^sum    =∫_I_z× I_z'G_k_j^sum(z,z')Ĵ_k_j^sum(z)Ĵ_k_j^sum^*(z') dzdz'
    𝒬_k_j^diff    =∫_I_z× I_z'G_k_j^diff(z,z')Ĵ_k_j^diff(z)Ĵ_k_j^diff^*(z') dzdz'.

By eq:self-interactions and eq:trans-interactions,

    Ĵ_k=[J_nbr]_k^∧=∑_j=1^nχ_jχ_j+1^'Θ_j,j+1(k)+χ_j^'χ_j+1Θ_j+1,j(k)

for k≠0. It follows from eq:nointerference that

    Ĵ_k_j^sum    =1/2(l_j+1/l_j-1)χ_jχ_j+1'+1/2(l_j/l_j+1-1)χ_j'χ_j+1
    Ĵ_k_j^diff    =1/2(1+l_j+1/l_j)χ_jχ_j+1'+1/2(1+l_j/l_j+1)χ_j'χ_j+1.

In particular, 

    supp Ĵ_k_j^sum∪supp Ĵ_k_j^sum⊂supp χ_jχ_j+1'∪supp χ_j'χ_j+1⊂[z_j,z_j+1]∪[1-z_j+1,1-z_j]=I_j

by eq:cross_support. Thus,

    𝒬_k_j^sum    =∫_I_j× I_jG_k_j^sum(z,z')Ĵ_k_j^sum(z)Ĵ_k_j^sum^*(z') dzdz'
    𝒬_k_i^diff    =∫_I_j× I_jG_k_j^diff(z,z')Ĵ_k_j^diff(z)Ĵ_k_j^diff(z') dzdz'.


Now we estimate these quadratic forms. By Hölder's inequality,

    𝒬_k_j^sum≤||G_k_j^sum||_L^1(I_j× I_j)||Ĵ_k_j^sum||_L^∞^2  and  𝒬_k_j^diff≤||G_k_j^diff||_L^1(I_j× I_j)||Ĵ_k_j^diff||_L^∞^2.

Observe that 

    ||χ_jχ_j+1'||_L^∞∨||χ_j'χ_j+1||_L^∞≲1/δ_j

as a result of eq:assumption2-0. Combining this with
the first part of lem:Gestimates applied with A=I_j,
we find that

    𝒬_k_j^sum≲δ_j/k_j^sum(δ_j∧1/k_j^sum)·(l_j/l_j+11/δ_j)^2  and  𝒬_k_j^diff≲δ_j/k_j^diff(δ_j∧1/k_j^diff)·(l_j/l_j+11/δ_j)^2.

It follows from eq:l_relations that 

    1/l_j≲ k_j^sum∧ k_j^diff,

so we can simplify these estimates to 

    𝒬_k_j^sum∨𝒬_k_j^diff≲δ_jl_j(δ_j∧ l_j)·(l_j/l_j+11/δ_j)^2.

Using the first part of eq:assumption2-1 followed by
eq:l_relations, we conclude that

    𝒬_k_j^sum∨𝒬_k_j^diff≲l_j^4/δ_jl_j+1^2≲l_j^2/δ_j.

This completes the proof.


We turn to estimate the higher order terms from eq:optimaldesign_w2wenstr.

The branching construction satisfies

    _Ω|∇𝐮|^2∨_Ω|∇ξ|^2≲1/l_1^2+∑_j=1^n-1δ_j/l_j+1^2+l_n^2/δ_n^3.



Note that 

    ∂_xxψ=∑_j=1^nχ_jψ_j”,  ∂_xzψ=∑_j=1^nχ_j'ψ_j', and  ∂_zzψ=∑_j=1^nχ_j”ψ_j.

Therefore, by orthogonality,

    _Ω|∇∇ψ|^2=∑_j=1^n_Ω|χ_jψ_j”|^2+2|χ_j'ψ_j'|^2+|χ_j”ψ_j|^2.

For j≠1, we see from an application of eq:assumption2-0
that

    _Ω|χ_jψ_j”|^2=||χ_j||_L^2(I_z)^2(ψ_j”)^2≲δ_j+δ_j-1/l_j^2≲δ_j-1/l_j^2.

For j=1, we have instead that

    _Ω|χ_1ψ_1”|^2=||χ_1||_L^2(I_z)^2(ψ_1”)^2≲1/l_1^2.


Similarly, we find that

    _Ω|χ_j'ψ_j'|^2    =||χ_j'||_L^2(I_z)^2(ψ_j')^2≲1/δ_j-1+1/δ_j≲1/δ_j,
    _Ω|χ_j”ψ_j|^2    =||χ_j”||_L^2(I_z)^2(ψ_j)^2≲(1/δ_j-1^3+1/δ_j^3)l_j^2≲l_j^2/δ_j^3

for all j. Therefore,

    _Ω|∇∇ψ|^2≲1/l_1^2+∑_j=1^n-1(δ_j/l_j+1^2+1/δ_j+l_j^2/δ_j^3)+1/δ_n+l_n^2/δ_n^3.

By eq:assumption2-1, 

    δ_j/l_j+1^2≳1/δ_j∨l_j^2/δ_j^3  and  l_n^2/δ_n^3≳1/δ_n.

 The result follows.


We now assemble the previous estimates. Let

    E(ϵ;{z_k},{l_k})=_Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ_Ω|∇𝐮|^2·_Ω|∇ξ|^2

where (𝐮,ξ) are constructed from {z_k}_k=1^n
and {l_k}_k=1^n as described in sec:branchingconstruction.
It will be convenient in what follows to think of estimating E
in terms of some smoothly interpolated version of these parameters. 

Let ℓ(z) be any smooth,
monotonic function defined on [z_bulk,z_bl] that satisfies

    ℓ(z_k)=l_k  1≤ k≤ n.

Then, the branching construction corresponding to {z_k}_k=1^n
and {l_k}_k=1^n satisfies 

    E(ϵ;{z_k},{l_k})≲ l_bl+∫_z_bulk^z_bl(ℓ'(z))^2 dz+ϵ(1/l_bulk^2+∫_z_bulk^z_bl1/(ℓ(z))^2 dz+1/l_bl)^2.



Collecting the results above and using eq:l_relations,
we conclude that

    _Ω|∇Δ^-1div(𝐮ξ)|^2≲ l_bl+∑_j=1^n-1l_j^2/δ_j∼ l_bl+∑_j=1^n-1|l_j+1-l_j/δ_j|^2δ_j

and


    _Ω|∇𝐮|^2∨_Ω|∇ξ|^2≲1/l_1^2+∑_j=1^n-1δ_j/l_j+1^2+l_n^2/δ_n^3∼δ_bulk/l_bulk^2+∑_j=1^n-11/l_j^2δ_j+1/l_bl.

By Jensen's inequality and the definition of ℓ(z), 

    |l_j+1-l_j/δ_j|^2=|_z_j^z_j+1ℓ'(z) dz|^2≤_z_j^z_j+1|ℓ'|^2 dz

so that

    _Ω|∇Δ^-1div(𝐮ξ)|^2≲ l_bl+∑_j=1^n-1∫_z_j^z_j+1|ℓ'|^2 dz=l_bl+∫_z_bulk^z_bl(ℓ')^2 dz.

Also, as l_j+1∼ l_j by eq:l_relations, 

    1/l_j^2δ_j=∫_z_j^z_j+11/l_j^2 dz∼∫_z_j^z_j+11/ℓ^2 dz.

Therefore, 

    _Ω|∇𝐮|^2∨_Ω|∇ξ|^2≲δ_bulk/l_bulk^2+∑_j=1^n-1∫_z_j^z_j+11/ℓ^2 dz+1/l_bl=δ_bulk/l_bulk^2+∫_z_bulk^z_bl1/ℓ^2 dz+1/l_bl.





 §.§ Proof of prop:lowerbound_enstrconstr_w2w


The result of the previous analysis is that the branching construction
from sec:branchingconstruction satisfies the efficiency
estimate

    E(ϵ;{z_k},{l_k})≲ l_bl+∫_z_bulk^z_bl(ℓ')^2 dz+ϵ(1/l_bulk^2+∫_z_bulk^z_bl1/(ℓ)^2 dz+1/l_bl)^2

where ℓ(z) is obtained from {l_k}_k=1^n by smooth
and monotonic interpolation. Now to prove prop:lowerbound_enstrconstr_w2w,
we will optimize the righthand side in the free parameters ℓ(z),
l_bulk, and l_bl, and then back out admissible choices of
{z_k}_k=1^n and {l_k}_k=1^n from the result.
To ensure that the requirements from sec:branchingconstruction-sec:Admissibility
hold, we must carry out this optimization under the constraint that

    0≤ℓ'(z)≲1   z∈[z_bulk,z_bl].

That the minimizer of 

    min_ℓ(z)
    
    l(z_bulk)=l_bulk
    ℓ(z_bl)=l_bl l_bl+∫_z_bulk^z_bl(ℓ')^2 dz+ϵ(1/l_bulk^2+∫_z_bulk^z_bl1/(ℓ)^2 dz+1/l_bl)^2

satisfies eq:monotonicity_of_l will be verified later
on. 

First, let us determine the optimal form of ℓ(z). We consider
that ϵ≪1 throughout this preliminary discussion, which
should serve to motivate the choices made in the formal proof that
follows. Consider the contributions to eq:1d_lengthscale_minpblm
coming from the transition layers where z∈[z_bulk,z_bl].
We can identify the scaling of their minimum value by balancing the
corresponding integrands. This yields

    ℓ'(z)∼ϵ^1/2(∫_z_bulk^z_bl1/ℓ^2 dz)^1/21/ℓ(z).

It is natural to impose the boundary condition ℓ(1)=0 to determine
ℓ. We find that 

    ℓ(z)∼ c(ϵ)(1-z)^1/2

where c(ϵ) must be determined by substitution into eq:balance.
Thus,

    c(ϵ)    ∼ϵ^1/6(∫_z_bulk^z_bl1/1-z dz)^1/6=ϵ^1/6log^1/6(1-z_bulk/1-z_bl).

Since

    l_bl∼ c(ϵ)(1-z_bl)^1/2  and   l_bulk∼ c(ϵ)(1-z_bulk)^1/2

we conclude that

    1-z_bulk/1-z_bl∼l_bulk/l_bl.

Anticipating that l_bl≪ l_bulk for ϵ≪1, we conclude
that the optimal form of the smooth lengthscale function ℓ(z)
is given by

    ℓ(z)∼ϵ^1/6log^1/6(l_bulk/l_bl) (1-z)^1/2,   z∈[z_bulk,z_bl].

Such an ℓ yields the estimates

    ∫_z_bulk^z_bl(ℓ')^2 dz    ∼ϵ(∫_z_bulk^z_bl1/ℓ^2 dz)^2∼ϵ^1/3log^1/3(l_bulk/l_bl)∫_z_bulk^z_bl1/1-z dz
        =ϵ^1/3log^1/3(l_bulk/l_bl)log(1-z_bulk/1-z_bl)∼ϵ^1/3log^4/3(l_bulk/l_bl)

for ϵ≪1. 

Next, we determine the optimal choices for l_bulk and l_bl
in this asymptotic regime. Plugging eq:powerlaw_for_l
back into eq:1d_lengthscale_minpblm yields the resulting
minimization

    min_l(z_bulk)=l_bulk
    ℓ(z_bl)=l_bl l_bl+ϵ^1/3log^4/3(l_bulk/l_bl)+ϵ(1/l_bulk^2+1/l_bl)^2.

Critical point tests yield the optimal scalings

    l_bulk∼ϵ^1/6log^1/61/ϵ  and   l_bl∼ϵ^1/3log^1/31/ϵ

for ϵ≪1. Note this is consistent with the hypothesis that
l_bl≪ l_bulk in this regime. To summarize, the smooth lengthscale
function ℓ(z) picked out by our analysis of eq:1d_lengthscale_minpblm
scales as

    ℓ(z)∼ϵ^1/6log^1/61/ϵ (1-z)^1/2,   z∈[z_bulk,z_bl],

where

    1-z_bulk∼1  and  1-z_bl∼ϵ^1/3log^1/31/ϵ.


We are now ready to prove the upper bound from prop:lowerbound_enstrconstr_w2w.


Our plan is to verify the existence of a branching construction, as
described in sec:branchingconstruction, whose parameters
{z_k}_k=1^n and {l_k}_k=1^n are consistent with
the optimal smooth lengthscale function ℓ(z) from eq:powerlaw_for_l-1.
Once we verify the requirements of sec:branchingconstruction-sec:Admissibility
hold, the desired bound E(ϵ;{z_k},{l_k})≲ϵ^1/3log^4/31/ϵ
follows as above. For the reader's convenience, we recall the requirements
that must be checked: these are eq:zk_inequalities, eq:lk_inequalities,
and eq:assumption2-0-eq:l_relations.

We start by defining 

    ℓ(z)=ϵ^1/6log^1/61/ϵ (1-z)^1/2   z∈[1/2,1]

in obvious analogy with eq:powerlaw_for_l-1. To choose
the horizontal lengthscales {l_k}_k=1^n, we set 

    l_bulk=l_x/2π1/k_bulk

where k_bulk∈ℕ satisfies

    k_bulk-1<l_x/π1/ϵ^1/6log^1/61/ϵ≤ k_bulk,

and take

    l_k=l_bulk/2^k-1,   k=1,…,n.

Note to ensure k_bulk≥1 we must require that ϵ^1/6log^1/61/ϵ≲ l_x.
This condition is given in the statement of prop:lowerbound_enstrconstr_w2w.
Note also that eq:lk_inequalities and eq:l_relations
hold. 

Now as eq:eqn_for_cts_l is strictly decreasing, we may
define the points {z_k}_k=1^n by 

    ℓ(z_k)=l_k,   k=1,…,n.

This gives

    z_k=1-c_11/2^2(k-1),   k=1,…,n

where

    c_1=(l_bulk/ϵ^1/6log^1/61/ϵ)^2=1/4π^21/ϵ^1/3log^1/31/ϵl_x^2/k_bulk^2.

By eq:bulk_wavenumber, c_1≤1/4 so that
z_1=1-c_1≥3/4 as required by eq:assumption3.
Note eq:zk_inequalities and eq:assumption2-0
are satisfied as well.

Finally, we fix n∈ℕ by enforcing eq:bl_bulk_conditions_for_l,
which states here that

    l_bulk/2^n-1∼ϵ^1/3log^1/31/ϵ.

To achieve this, let us define n∈ℕ via the inequalities

    n-1<log_2(2π/ϵ^1/3log^1/31/ϵ1/k_bulk)≤ n.


Having chosen {z_k}_k=1^n and {l_k}_k=1^n, we
may invoke the definitions from sec:branchingconstruction
to produce a branching construction (𝐮,ξ). Note we have
checked each requirement from sec:branchingconstruction-sec:Admissibility
except for eq:assumption2-1. That l_n∼δ_n
follows from eq:bl_bulk_conditions_for_l and eq:powerlaw_for_l-1.
Now we show that l_k≲δ_k for all k. Since δ_k=z_k+1-z_k
and l_k∼|l_k+1-l_k|, this requires showing that

    1≲|z_k+1-z_k/l_k+1-l_k|

for all k. Noting z'(ℓ)<0, we only need to show that

    1≲|z'(ℓ)|  ℓ∈[l_bl,l_bulk].

Differentiating eq:eqn_for_cts_l implicitly, we find
that

    |z'(ℓ)|=2(1-z)^1/2/ϵ^1/6log^1/61/ϵ≳δ_bl^1/2/ϵ^1/6log^1/61/ϵ∼l_bl^1/2/ϵ^1/6log^1/61/ϵ∼1

as required.

In sum, we have produced a branching construction (𝐮,ξ)
consistent with the requirements of sec:branchingconstruction-sec:Admissibility
whose parameters {z_k}_k=1^n and {l_k}_k=1^n
interpolate the desired smooth lengthscale function ℓ(z) from
eq:powerlaw_for_l-1. The estimates proved in sec:Energy-estimates
apply, and we may immediately conclude from cor:branchingestimate-final
and the discussion surrounding eq:bl_bulk_conditions_for_l
and eq:powerlaw_for_l-1 that 

    E(ϵ;{z_k},   {l_k})  =_Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ_Ω|∇𝐮|^2·_Ω|∇ξ|^2
       ≲ l_bl+∫_z_bulk^z_bl(ℓ')^2 dz+ϵ(1/l_bulk^2+∫_z_bulk^z_bl1/(ℓ)^2 dz+1/l_bl)^2
       ≲ϵ^1/3log^1/31/ϵ+ϵ^1/3log^4/3(ϵ^1/6log^1/61/ϵ/ϵ^1/3log^1/31/ϵ)+ϵ(2/ϵ^1/3log^1/31/ϵ)^2
       ≲ϵ^1/3log^1/31/ϵ+ϵ^1/3log^4/31/ϵ+ϵ^1/3/log^2/31/ϵ≲ϵ^1/3log^4/31/ϵ

for ϵ≲1. Thus, prop:lowerbound_enstrconstr_w2w
is proved. 



§ IMPLICATIONS FOR THE ANALYSIS OF TURBULENT HEAT TRANSPORT


There is a long history, originating in the works of Malkus <cit.>
and Howard <cit.>, of variational methods for the
analysis of turbulent heat transport, the primary focus of which has been on absolute or a priori upper bounds. 
Consider the usual setup
of Rayleigh-Bénard convection (RBC), wherein an incompressible
fluid layer is heated from below and cooled from above, and is subjected
to a constant downwards-pointing gravitational force. The temperature
field T(𝐱,t) undergoes transport by means of advection-diffusion,

    ∂_tT+𝐮·∇ T=Δ T.

The advecting velocity 𝐮(𝐱,t) is coupled back
to temperature field T through a suitable momentum equation. This
could be, for instance, Darcy's law as it is for convection in a fluid
saturated porous layer. Here, we are concerned with convection in
a fluid layer for which, in the Bousinessq approximation, eq:adv-diff-eqn
is supplemented with the buoyancy forced incompressible Navier-Stokes
equations

    ∂_t𝐮+𝐮·∇𝐮+∇ p=PrΔ𝐮+PrRa𝐤̂T

and

    div 𝐮=0.

The two non-dimensional parameters are the Prandtl number Pr, the
ratio of the fluid's kinematic viscosity to its thermal diffusivity,
and the Rayleigh number Ra, a ratio of the intensities of driving
to damping forces which is proportional here to the bulk buoyancy
force across the layer. Altogether, eq:adv-diff-eqn-eq:div-free-1
constitute the equations of Rayleigh-Bénard convection in a fluid
layer <cit.>. For boundary conditions we continue to
assume that the temperature field is imposed at the top and bottom
of the layer by

    T|_z=1=0  and   T|_z=0=1,

while the velocity field is taken to satisfy either the no-slip boundary conditions

    𝐮|_∂Ω=0

or the stress-free boundary conditions

    w|_∂Ω=0  and  ∂_zu|_∂Ω=∂_zv|_∂Ω=0.

All fields are assumed to be periodic in the xy-plane.

The rate of heat transport in RBC can be measured by the Nusselt number
Nu, which evidently depends on Pr and Ra in some unknown
and complicated way. (It can also depend on the initial data, as well
as on the aspect ratios of the fluid layer.) Determining this relationship
and/or establishing absolute bounds on it continues to be the subject
of numerous works across the physical and mathematical literatures.
To date, the best known upper bound holding uniformly in Pr and for no-slip
velocity boundary conditions states that

    Nu≲ Ra^1/2

for Ra≫1 <cit.>.
This bound also holds for stress-free velocity boundary conditions
in the three-dimensional layer Ω=𝕋_xy^2× I_z,
but more is known in the two-dimensional case where Ω=𝕋_x× I_z:
in two dimensions with stress-free boundary conditions, one has that
Nu≲ Ra^5/12 uniformly in Pr for Ra≫1 <cit.>.
(In the formal limit where Pr=∞ and eq:momentum-eqn
is replaced with Stoke's equation, the situation is quite different
<cit.>.)
There is little to no evidence, however, that any of these finite
Pr bounds are in fact sharp, i.e., that there exist solutions of the equations of motion
eq:adv-diff-eqn-eq:div-free-1 satisfying
Nu∼ Ra^1/2 as Ra→∞ (or Nu∼ Ra^5/12 for stress-free
boundaries in two dimensions). 
















































In light of all this, we note that the main fluid dynamical contribution of this paper is a proof that when the momentum equation eq:momentum-eqn is replaced by the enstrophy-constraint

    ⟨ |∇𝐮|^2⟩ =Ra(Nu-1)

which it implies, the upper bound Nu≲  Ra^1/2 becomes asymptotically sharp up to logarithmic corrections.  That is, for all large enough Ra there exist velocity and temperature fields satisfying eq:adv-diff-eqn, eq:div-free-1, and eq:momentum-balance along with the requisite boundary conditions such that

    Ra^1/2/log^2Ra≲ Nu ≲ Ra^1/2.

This follows from thm:mainbounds_w2w_ens upon taking Pe^2 = Ra(Nu-1). Therefore, either the well-known bound eq:Ra1/2bd on RBC is asymptotically sharp as Ra→∞ and Pr is fixed, or details from the momentum equation eq:momentum-eqn beyond the balance eq:momentum-balance are essential for determining the scaling law of maximal turbulent heat transport. 



































































































































The remainder of this section places our analysis of wall-to-wall optimal transport into its proper fluid dynamical context. 
To keep the discussion at a reasonable length, we do not attempt to summarize the vast literature on the
subject but instead focus on two of the most well-known methods for
proving a priori bounds on transport: the variational approach
of Howard, and the background method of Constantin-Doering. Our plan is to recall just enough about these methods to allow for comparison with the techniques developed in this paper. 
For Howard's approach see sec:Howardapproach, while for the background method see sec:onbackgroundmethod. sec:momentumeqn concerns the role of the momentum equation.

Before we proceed, let us mention the existence of the recently developed “auxiliary functional” method for producing bounds on time-averaged quantities <cit.>.
While the background method may ultimately be derived by a particular
choice of auxiliary functional —the same is true for the
recently proposed method of Seis <cit.> —it
is not yet clear if there exists any auxiliary functional that yields
an improvement to scaling beyond Nu≲ Ra^1/2. Although
for ordinary differential equations the auxiliary functional method
always yields sharp bounds on long-time averages <cit.>,
it remains to be seen if such a situation holds for general PDEs.




 §.§ On the variational approach of Howard




  §.§.§ Howard's variational problem


If RBC is to be taken as a predictive model for turbulent convection,
one naturally asks: which of its solutions are actually realizable by experiment? Setting aside dynamical stability
as a possible selection principle, Malkus introduced in <cit.> the idea that perhaps amongst all possible solutions of the equations of motion, those that are realized maximize their heat transport overall. An operational approach to establishing upper bounds inspired by Malkus' idea is to search for a larger admissible set of velocity and temperature fields, which contains all solutions of RBC, amongst which the maximal transport can analytically
be determined. This is Howard's variational approach.

Following Howard <cit.>, we observe that if 𝐮 and T arise in RBC, they must satisfy two identities known as the “power integrals”. To derive the first of these, dot the momentum equation
eq:momentum-eqn into 𝐮, integrate by parts
and average in space and time. Changing variables by θ=T-(1-z) yields
 the first of Howard's identities

    Ra⟨ wθ⟩ =⟨ |∇𝐮|^2⟩ .

(Note this is simply a restatement of eq:momentum-balance
from above.) A similar manipulation involving the
temperature equation eq:adv-diff-eqn yields the second
identity

    ⟨ wθ⟩ +⟨ wθ⟩ ^2-⟨ |wθ|^2⟩ =⟨ |∇θ|^2⟩ .

Consider now the problem of maximizing Nu amongst
all divergence-free vector fields 𝐮 and scalar fields
θ that vanish at the walls and furthermore satisfy eq:powerintegral1
and eq:powerintegral2. Since the equations of motion
of RBC imply these constraints, the resulting maximum sets
an upper bound on Nu for RBC. 


Setting aside matters of statistical stationarity <cit.>,
one can give an equivalent formulation of the variational problem
described above which makes it tractable for analysis. Under certain
further assumptions on the solutions of RBC (the “requirements of
homogeneity” from <cit.>), Howard deduced that the minimization

    min_𝐮(𝐱),θ(𝐱)
    𝐮|_∂Ω=0,θ|_∂Ω=0
    _Ωwθ=1
     _Ω|wθ-1|^2+ϵ_Ω|∇𝐮|^2_Ω|∇θ|^2

is equivalent to the maximization sup Nu described above, and
that its optimal value can be used to produce an a priori bound
on RBC (the algebraic manipulations in the proof of this are like
those performed in sec:OptlDesignProb_w2w in the derivation
of the integral formulation of steady wall-to-wall optimal transport). 

The minimization eq:Howard's_pblm is known as Howard's
problem. It bears striking resemblance to our integral formulation
of steady enstrophy-constrained wall-to-wall transport

    min_𝐮(𝐱),ξ(𝐱)
    𝐮|_∂Ω=0,ξ|_∂Ω=0
    _Ωwξ=1
     _Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ_Ω|∇𝐮|^2_Ω|∇ξ|^2,

obtained in sec:OptlDesignProb_w2w. To find the relationship
between eq:Howard's_pblm and eq:steady_optl_design_enstrophy,
we apply lem:advection_Fourierests_w2w along with the
net-flux constraint _Ωwξ=1 and decompose the advection
term as

    _Ω|∇Δ^-1div(𝐮ξ)|^2=_Ω|wξ-1|^2+𝒬(div(𝐮ξ))

where 𝒬 is the positive semi-definite quadratic form
defined in lem:advection_Fourierests_w2w. 

This last equation reveals the precise distinction between Howard's problem eq:Howard's_pblm and our integral formulation in eq:steady_optl_design_enstrophy. Because 𝒬 is positive semi-definite, it is evident that the minimum in eq:Howard's_pblm is not smaller than the minimum in eq:steady_optl_design_enstrophy. As a result, Howard's upper bound on heat transport is not lower than ours. Though the improvement in scaling in our approach is limited by eq:log_Nu_bds to at most a logarithmic correction, it remains to be seen whether such a correction holds as an absolute upper bound. We turn now to consider the difference between the optimizers of eq:Howard's_pblm and eq:steady_optl_design_enstrophy.  




  §.§.§ Busse's multi-α technique


As shown by Howard and Busse <cit.>,
the optimal value of Howard's problem eq:Howard's_pblm
scales as ϵ^1/3 for ϵ≪1. Thus, Howard's approach
to bounds on RBC yields Nu≲ Ra^1/2 and no better. The
a priori lower bound implicit in this result is due to Howard;
the upper bound was obtained by Busse as an application of his “multi-α”
technique, which seeks to produce asymptotically valid solutions of
the Euler-Lagrange equations of eq:Howard's_pblm involving
multiple horizontal wave numbers. Busse's multi-α analysis turns out to share parallels with our construction of branching flows, which we would like to discuss now. 

We start by recalling Howard's lower bound:


    min_𝐮(𝐱),θ(𝐱)
    𝐮|_∂Ω=0,θ|_∂Ω=0
    _Ωwθ=1
     _Ω|wθ-1|^2+ϵ_Ω|∇𝐮|^2_Ω|∇θ|^2≳ϵ^1/3

for ϵ≪1. Let (𝐮,θ) be admissible, which we can take to be smooth. Let
δ∈(0,1/2) be such that 

    0≤|wθ|≤1/2  for z∈[0,δ],  and  wθ(δ)=1/2.

(If there does not exist such a δ, then _Ω|wθ-1|^2≳1≫ϵ^2/3.)
By its definition,

    _Ω|wθ-1|^2≥1/|I_z|∫_0^δ|wθ-1|^2≳δ.

lem:lingrowth_wtheta states that

    1/|𝕋_xy^2|||∂_zθ||_L^2(Ω)||∂_zw||_L^2(Ω)≳|wθ(z)|/|z∧(1-z)|  ∀ z.

Taking z=δ and squaring, we conclude that

    _Ω|∇𝐮|^2_Ω|∇θ|^2≳1/δ^2.

Therefore, the optimal value in the lefthand side of eq:Howard_LB
is bounded below by 

    inf_δ∈(0,1/2){δ+ϵ1/δ^2}∼ϵ^1/3

for ϵ≪1, and eq:Howard_LB is proved. 

Now we discuss Busse's upper bound: it asserts the existence
of admissible pairs {(𝐮_ϵ,θ_ϵ)}
satisfying 

    _Ω|w_ϵθ_ϵ-1|^2+ϵ_Ω|∇𝐮_ϵ|^2_Ω|∇θ_ϵ|^2≲ϵ^1/3

for ϵ≪1. Busse's multi-α technique is analogous to our branching construction from sec:enstrophy_constrained_design.
Arguing as in that section, we find that our branching construction
with lengthscale ℓ(z) satisfies the estimates

    _Ω|wθ-1|^2≲ l_bl  and  _Ω|∇𝐮|^2_Ω|∇θ|^2≲(1/l_bulk^2+∫_z_bulk^z_bl1/ℓ^2 dz+1/l_bl)^2

so long as 0≤ℓ'(z)≲1. Since branching is admissible
for Howard's problem, we find its optimal value is bounded above by

    min_ℓ(z)
    ℓ(z_bulk)=l_bulk
    ℓ(z_bl)=l_bl
    
    0≤ℓ'(z)≲1
     l_bl+ϵ(1/l_bulk^2+∫_z_bulk^z_bl1/ℓ^2 dz+1/l_bl)^2.

Choosing

    ℓ(z)    ∼1-z   z∈[z_bulk,z_bl],
    
    l_bulk   ∼1, l_bl∼ϵ^1/3

yields eq:Busse_UB. Although Busse's construction
is usually described in terms of discrete wavenumbers {α_k}_k=1^n
and points {z_k}_k=1^n, for ϵ≪1 these can be
seen to arise from interpolation of the continuous lengthscale ℓ(z)∼1-z,
similar to the presentation in sec:enstrophy_constrained_design.

Coming back to wall-to-wall optimal transport, we can now discuss the difference between the optimizers of Howard's problem eq:Howard's_pblm and our integral formulation in eq:steady_optl_design_enstrophy. As the analysis in sec:enstrophy_constrained_design indicates, adding 𝒬 to Howard's problem eq:Howard's_pblm should change the preferred lengthscale for branching from Busse's linear law ℓ∼1-z to our square root one ℓ∼ c(ϵ)√(1-z). The estimates obtained there show that

    𝒬∼∫_z_bulk^z_bl(ℓ')^2dz.

Thus, the 1D problem eq:1dpblm_Busse for selecting the lengthscale function ℓ turns into eq:1d_w2w for wall-to-wall optimal transport. It remains to be seen whether the true optimizers of eq:steady_optl_design_enstrophy exhibit branching with these preferred lengthscales. Presumably, developing such fine detailed knowledge of the minimizers would help resolve the question of logarithmic corrections to scaling.  



 §.§ On the background method




  §.§.§ Background method for RBC


In <cit.>, Constantin and one of the authors introduced
an alternate method to Howard's for establishing a priori bounds
on RBC, which can be applied without any assumptions of statistical
stationarity or homogeneity. We recall the argument now, with the
goal of connecting it to the symmetrization method from
sec:aprioribds. We follow the presentation in <cit.>.

Let 𝐮 and T arise from RBC and decompose the temperature
field into the sum of stationary “background” and fluctuating
parts,

    T(𝐱,t)=τ(z)+θ(𝐱,t)

where τ(0)=1 and τ(1)=0. Then, 

    1/2d/dt(_Ω|θ|^2+1/PrRa_Ω|𝐮|^2)+1/2_Ω|∇ T|^2=1/2∫_0^1|τ'|^2-H_τ(𝐮,θ)

where H_τ is the quadratic form

    H_τ(𝐮,θ)=_Ω1/Ra|∇𝐮|^2+1/2|∇θ|^2+wθ(τ'-1).

Provided that H_τ≥0 for all divergence-free vector fields
𝐮(𝐱) and scalar fields θ(𝐱)
vanishing at ∂Ω, we can drop the last term from the
dissipation equation and take a long-time average to find the inequality

    ⟨ |∇ T|^2⟩≤∫_0^1|τ'|^2.

This proves the following variational bound:

    Nu≤inf_τ(z)
    τ(0)=1,τ(1)=0
    
    H_τ≥0
     ∫_0^1|τ'|^2.

Those background fields τ which satisfy H_τ≥0 are
known as spectrally stable. 

As proved in <cit.>, there exist spectrally stable background
fields {τ_δ} satisfying 

    ∫_0^1|τ_δ'|^2∼1/δ

for all δ≤ Ra^-1/2. Minimizing the resulting bound Nu≲1/δ
over this range of δ proves that Nu≲ Ra^1/2. We
note the remarkable similarity between the background fields constructed
in <cit.> and those constructed for the symmetrization
method in eq:testfunctions, which we turn to discuss now.



  §.§.§ Background method for optimal transport


As observed in <cit.>, one can obtain a priori bounds on
optimal transport via a suitable modification of the background method.
Here, our goal is to show that the symmetrization method from sec:aprioribds,
when properly abstracted and optimized, yields an a priori bound
on transport whose value is exactly the same as that obtained
in <cit.>. This begs the question of whether better background fields might
be constructed to improve upon the scaling Nu≲ Ra^1/2 (albeit by at most a logarithmic amount). Numerical evidence looks to point in the opposite direction, as the optimal bounds found in <cit.> scale ∼ Ra^1/2. We are not aware of a proof demonstrating this at the present time.


The modified background method from <cit.> is as follows. Let
T solve the advection-diffusion equation eq:adv-diff-eqn.
Performing the background decomposition

    T(𝐱,t)=τ(z)+θ(𝐱,t)

with τ(0)=1 and τ(1)=0 and introducing a Lagrange multiplier
λ∈ℝ, we find that

    1/2d/dt_Ω|θ|^2+1/2_Ω|∇ T|^2=1/2∫_0^1|τ'|^2+λ/2Pe^2-H_τ,λ(𝐮,θ)

where H_τ,λ is the quadratic form

    H_τ,λ(𝐮,θ)=_Ωλ/2|∇𝐮|^2+1/2|∇θ|^2+wθτ'.

If H_τ,λ≥0 for all divergence-free vector fields
𝐮(𝐱) and scalar fields θ(𝐱)
vanishing at ∂Ω, the dissipation equation eq:energy-dissapation-w2w
implies that

    ⟨ |∇ T|^2⟩≤∫_0^1|τ'|^2+λ Pe^2.

Thus,

    Nu≤inf_τ(z),λ
    τ(0)=1,τ(1)=0
    
    H_τ,λ≥0
    {∫_0^1|τ'|^2+λ Pe^2}

In parallel with the background method discussed above, we 
refer to background fields τ satisfying H_τ,λ≥0
as being spectrally stable at Lagrange multiplier λ. 

On the other hand, the symmetrization method from sec:aprioribds
yields the bound

    sup_𝐮(𝐱,t)
    ⟨ |∇𝐮|^2⟩ =Pe^2
    𝐮|_∂Ω=0
     Nu(𝐮)≤inf_η(𝐱)
    η|_z=1=0,η|_z=1=1
    {_Ω|∇η|^2+Pe^2sup_𝐮(𝐱)
    𝐮|_∂Ω=0
    _Ω|∇𝐮|^2=1
     _Ω|∇Δ^-1div(𝐮η)|^2}

when carried out optimally. As it turns out, these bounds are one and the same. 

Let U_bm(Pe) and U_symm(Pe) denote the optimal values appearing
on the righthand sides of eq:backgroundmethod_w2w and
eq:symmmethod_w2w, respectively. We have that U_bm=U_symm.


As the following proof shows, the minimization in eq:symmmethod_w2w
can be performed over η depending z alone without changing the resulting value.


We prove this in two steps: first we show that U_symm≤ U_bm
and then we prove the reverse inequality. In both cases, we will use
the fact that

    ∫_Ω|∇Δ^-1div 𝐦|^2=sup_θ(𝐱)
    θ|_∂Ω=0
     ∫_Ω2𝐦·∇θ-|∇θ|^2

for all 𝐦∈ L^2(Ω;ℝ^3). 

We begin by showing that U_symm≤ U_bm. Taking 𝐦=𝐮τ
in eq:duality_equality, we see that a background field
τ(z) satisfies H_τ,λ≥0 if and only if

    ∫_Ω|∇Δ^-1div 𝐮τ|^2≤λ∫_Ω|∇𝐮|^2

for all divergence-free 𝐮 that vanish at ∂Ω.
Therefore, 

    U_bm≥inf_τ(z),λ
    τ(0)=1,τ(1)=0
    
    H_τ,λ≥0
    {∫_0^1|τ'|^2+Pe^2sup_𝐮(𝐱)
    𝐮|_∂Ω=0
    _Ω|∇𝐮|^2=1
     _Ω|∇Δ^-1div(𝐮τ)|^2}≥ U_symm

since enlarging the admissible set only decreases the resulting minimal
value.

Now we prove that U_symm≥ U_bm. Parameterizing the admissible
set from eq:symmmethod_w2w via the level sets of 

    M(η)=sup_𝐮(𝐱)
    𝐮|_∂Ω=0
    _Ω|∇𝐮|^2=1
     _Ω|∇Δ^-1div(𝐮η)|^2,

we can write that 

    U_symm=inf_λinf_η(𝐱)
    η|_z=1=0,η|_z=1=1
    
    M(η)=λ{_Ω|∇η|^2+Pe^2λ} .

Extend the definition of H_τ,λ to functions of 𝐱
by taking

    H_η,λ(𝐮,θ)=_Ωλ/2|∇𝐮|^2+1/2|∇θ|^2+𝐮θ·∇η.

By eq:duality_equality,

    M(η)=λinf_𝐮(𝐱),θ(𝐱)
    𝐮|_∂Ω=0,θ|_∂Ω=0
     H_η,λ(𝐮,θ)=0

and the latter happens if and only if H_η ,λ≥ 0.














Using that
_Ω|∇η|^2 is convex in η and that {η:H_η,λ≥0}
is also convex, we can replace η with its periodic average
τ=η to deduce that

    U_symm≥inf_λinf_τ(z)
    τ|_z=1=0,τ|_z=1=1
    
    H_τ,λ≥0
    {_Ω|τ'|^2+Pe^2λ} =U_bm

as desired.





 §.§ On the realizability of optimal heat transport by buoyancy-driven convection


We return to the full system eq:adv-diff-eqn-eq:div-free-1 now. One may wonder if buoyancy forces are capable of producing flows,
time-dependent or steady, that realize near-optimal heat transport.
The answer depends upon the way in which flow intensity is constrained.

First, we note that the energy-constrained wall-to-wall optimal transport problem corresponds
to RBC in a fluid saturated porous layer where the Navier-Stokes momentum
equation eq:momentum-eqn is replaced by Darcy's law.
This implies the balance law ⟨ |𝐮|^2⟩ =Ra(Nu-1)
which, when combined with the result of thm:mainbounds_w2w_eng,
yields the optimal scaling Nu∼ Ra in this setting. Direct numerical simulations
of time-dependent high-Ra porous medium convection <cit.>
are consistent with this scaling, indicating that buoyancy forces
can produce flows realizing optimal heat transport insofar
as scaling is concerned. On the other hand, asymptotic and numerical
investigations indicate that the best possible transport by steady
flows satisfies Nu∼ Ra^0.6 <cit.>.

Second, we observe that the enstrophy-constrained optimal transport problem corresponds to
Rayleigh's original model of buoyancy-driven convection in a fluid
layer <cit.>. There, steady convection also appears
to be strongly sub-optimal with the highest computationally observed
scaling being Nu∼ Ra^0.31 <cit.>.
To date, there are no turbulent high-Ra direct numerical simulations
indicating heat transport scaling much higher than Nu∼ Ra^1/3.

We close our discussion of fluid dynamical implications by commenting
on the certain sub-optimality of heat transport in Rayleigh's original model.
Rayleigh imposed eq:adv-diff-eqn-eq:div-free-1
in two-dimensions with stress-free velocity boundary conditions and
the usual Dirichlet temperature ones. Although RBC in a fluid layer
must obey the bound Nu≲ Ra^1/2 in any dimension and for
any boundary conditions, the result of <cit.>
is that in two-dimensions and with stress-free boundaries
Nu≲ Ra^5/12. Nevertheless, by combining the relevant balance law ⟨ |∇𝐮|^2⟩ =Ra(Nu-1) implied by the Navier-Stokes momentum equation eq:momentum-eqn with the result of thm:mainbounds_w2w_ens and the remark immediately thereafter, we conclude that optimal heat transport in the setting of Rayleigh's model must satisfy Nu∼ Ra^1/2 (up to logarithmic corrections). 
Our analysis is consistent with all the requirements of Rayleigh's model except for the Navier-Stokes momentum equation eq:momentum-eqn.
Thus, buoyancy-driven convection in two-dimensions between stress-free
boundaries must yield strongly sub-optimal rates of heat transport
as compared with what happens if eq:momentum-eqn is not
imposed. This underscores the importance of using the momentum equation —rather
than only a balance law it implies —for determining the
asymptotic heat transport of turbulent RBC.


















§ OPTIMAL TRANSPORT AS ENERGY-DRIVEN PATTERN FORMATION 


There is a second scientific context, other than the fluid dynamical one, in which the methods behind our analysis of wall-to-wall optimal heat transport have played a fundamental role. This is the subject of “energy-driven pattern formation” in mathematical materials science  <cit.>. 

Perhaps the key methodological contribution of this paper is the reformulation
of the general steady wall-to-wall optimal transport problem

    sup_𝐮(𝐱)
    
    ||𝐮||=Pe
    
    +b.c.
     Nu(𝐮)

in its integral form

    inf_𝐮(𝐱),ξ(𝐱)
    _Ωwξ=1
    
    +b.c.
     _Ω|∇Δ^-1div(𝐮ξ)|^2+ϵ||𝐮||^2_Ω|∇ξ|^2.

This change of viewpoint, accomplished in sec:OptlDesignProb_w2w,
hinges on the fact that the Nusselt number of a steady velocity field
𝐮 can be written as the maximal value of a certain non-local
functional in ξ. The resulting problem eq:integral_formulation
is equivalent to the original one eq:steady_w2w_pblm,
and optimizers correspond. In the examples of energy- and enstrophy-constrained
optimal transport considered in sec:energy-constraineddesign
and sec:enstrophy_constrained_design, where ||·||
is the (volume-averaged) L^2- or Ḣ^1-norm, the integral
formulation eq:integral_formulation plays a key role
in the construction of divergence-free velocity fields that achieve
nearly optimal transport. As that analysis shows, the complexity of
the successful construction —whether it can be described
using few lengthscales or many —depends strongly on the
choice of norm. 

Besides its practical use for the estimation of optimal transport,
eq:integral_formulation shares striking similarities
with other non-convex and singularly perturbed variational problems
from mathematical materials science. The study of patterns selected
by energy minimization principles in this field is known as energy-driven pattern
formation. It is important to note that the wall-to-wall optimal transport problem is variational by definition.
Thus, our observation is not that there exists some variational formulation
for it, but rather that the specific formulation eq:integral_formulation
reminds of various model problems from energy-driven pattern formation.
From this point of view, it is no surprise that the (nearly) optimal
patterns constructed in this paper for wall-to-wall transport —convection
rolls and branching flows —bear similarities with other
well-appreciated patterns from materials science including domain
branching in micromagnetics <cit.>
and wrinkling cascades in thin elastic sheets <cit.>.
What eq:integral_formulation offers is a functional analytic
framework in which to make such connections precise. 

We discuss below two model problems from energy-driven pattern formation
and their connections to wall-to-wall optimal transport. We leave
their general scientific introduction to the references therein, focusing
instead on the salient features of their analysis. This discussion
provides an alternate viewpoint on the role of branching patterns
in the variational analysis of transport, which complements the older
purely fluid dynamical arguments of Busse <cit.>.
We hope these remarks prove useful to the reader interested in our
approach.



 §.§ Magnetic domain branching in a uniaxial ferromagnet


Our first example comes from micromagnetics and concerns the patterns
formed by magnetic domains in a uniaxial ferromagnet. The energetic
description is as follows. We take as the magnet the domain Ω=(-L,L)_x×[0,1]_y,z^2
where x is the the preferred direction of magnetization and L
is the magnet's (non-dimensionalized) length. On Ω we define
a magnetization vector field 𝐦(𝐱)=m_1î+m_2ĵ+m_3k̂
which is required to be of unit size |𝐦|=1, and is extended
by zero to the rest of space ℝ^3\Ω. The
micromagnetic energy that results is

    ∫_all space|∇Δ^-1div 𝐦|^2+∫_magnetQ(1-m_1^2)+ϵ|∇𝐦|

where the divergence is understood in the distributional sense. Strictly
speaking, this is a “sharp interface” model in which the total
variation norm

    ∫_Ω|∇𝐦|=∑_i=1^3∫_Ω|∇ m_i|=∑_i=1^3sup_𝐯∈ C_c^1(Ω;ℝ^n)
    
    ||𝐯||_L^∞(Ω)≤1
     ∫_Ωm_idiv 𝐯

features instead of the Ḣ^1-norm (for more on this reduction
see <cit.>). The first term appearing in eq:micromag_energy
is called the magnetostatic energy; it accounts for the cost of the
magnetic field induced by 𝐦 in the ambient space. The
second term is the anisotropy energy and it arises from an underlying
crystalline anisotropy which prefers 𝐦 to be ±î.
The third term is the interfacial energy. It permits 𝐦
to be discontinuous, but limits the total area of any interfaces across
which 𝐦 jumps. The parameters Q and ϵ set
the relative strengths of these effects. The magnetostatic and interfacial
energies have direct analogs in the wall-to-wall problem eq:integral_formulation;
the anisotropy term does not. Note that, due to the constraint |𝐦|=1,
this functional is non-convex.

There are various designs for 𝐦 one can entertain in minimizing
eq:micromag_energy. One is the so-called Kittel structure,
in which 𝐦 is independent of x and ±î-valued
throughout the magnet, alternating between these at some to be determined
lengthscale l in the yz-plane. This design costs no anisotropic
energy and the optimal l is selected by minimizing its magnetostatic
and interfacial costs. Another important design is the Landau-Lifshitz
structure, in which 𝐦 is independent of x and ±î-valued
except for in a thin boundary layer near x=± L. There, it is
taken to be perpendicular to î in such a way as to eliminate
the magnetostatic energy completely, thus coupling the thickness of
the boundary layer to the lengthscale l of oscillations in the
bulk. This is a sharp-interface version of the convection roll design
described in sec:energy-constraineddesign. Finally, there
is the Privorotskiĭ construction, which plays the role of the branching
flows from sec:enstrophy_constrained_design. It too involves
a very large number of distinct lengthscales which interpolate between
a preferred lengthscale in the bulk l_bulk and a significantly
smaller one at the boundary l_bl. We refer the reader for more
details to <cit.> including a description
of the relevant regimes.

What can be proved regarding this non-convex, non-local minimization
problem? Following the reference <cit.> we assume
that 𝐦(x,y,z) is periodic in (y,z) and identify [0,1]_y,z^2
with 𝕋_y,z^2. Then there exist positive constants
C and C' such that the minimum micromagnetic energy satisfies

    CQ^1/3ϵ^2/3L^1/3≤minimum micromagnetic energy≤ C'Q^1/3ϵ^2/3L^1/3

for all sufficiently large Q and sufficiently small ϵ/L.
The proof of this result requires two kinds of arguments. The upper
bound comes from estimating the cost of an optimal Privorotskii construction
(the conditions on Q, ϵ, and L ensure that the result
is significantly less than those obtained by the Kittel and Landau-Lifshitz
structures). The lower bound asserts that the Privorotskii construction
cannot be beat as far as scaling is concerned. The original proof
of it can be found in <cit.>, but we note the existence
of a second more recent proof in <cit.> which
utilizes the end-point Gagliardo-Nirenberg interpolation inequality

    ||f||_L^4/3(𝕋^2)≲||∇ f||_L^1(𝕋^2)^1/2||f||_H^-1(𝕋^2)^1/2

holding for all mean-zero and periodic functions f. 



 §.§ Blistering patterns in thin elastic sheets


Our second example comes from elasticity theory. Consider a thin elastic
sheet of (non-dimensional) thickness h which is strongly bonded
to the top of a large rubber block, except for on some known sub-domain
Ω⊂ℝ^2. Applying biaxial compression to the
block causes the sheet to blister in the unbonded domain. The result
is a complex pattern of wrinkles and folds whose details can be modeled
through the minimization of a certain non-convex and singularly perturbed
variational problem. As in <cit.>,
we consider minimization of the internal elastic energy under clamped
boundary conditions. In the Föppl-von Karman model, the elastic
energy (per unit thickness) is given by

    ∫_blistered
    region|e(𝐯)+1/2∇ϕ⊗∇ϕ|^2+h^2|∇∇ϕ|^2

where the “in-plane” displacement parallel to the top of the block
is 𝐯(𝐱) and the “out-of-plane” displacement
perpendicular to it is ϕ(𝐱). Here, e(𝐯)
denotes the symmetric part of the in-plane displacement gradient ∇𝐯.
Taken together, the in- and out-of-plane displacements yield the map
(𝐱,0)↦(𝐱+𝐯(𝐱),ϕ(𝐱))
which describes the deformation of the blister. At the edge of the
blister ∂Ω we impose the clamped boundary conditions

    𝐯|_∂Ω=-λ𝐱, ϕ|_∂Ω=0, and ∂_νϕ|_∂Ω=0.

The parameter λ is positive and sets the amount of overall
compressive strain. The first term in the energy is called the membrane
term. It prefers the in-plane strain e(𝐯)+1/2∇ϕ⊗∇ϕ
to vanish. The second one is called the bending term, and it prefers
the out-of-plane displacement to vary on longer lengthscales or not
at all. The relative strength of these effects is determined by the
parameter h, which is understood to be small.

There are significant parallels between the elastic energy functional
eq:FvK and the integral formulation of wall-to-wall transport
eq:integral_formulation. Of course, the bending term
from eq:FvK and the higher order terms from eq:integral_formulation
act to regularize designs. More interestingly, we observe a similarity
between the membrane term from eq:FvK and the advection
term and net-flux constraint from eq:integral_formulation.
Let us introduce a streamfunction ψ for the divergence-free
velocity field 𝐮 (we work with a two-dimensional fluid
layer now) and rewrite the advection term as

    _fluid layer|∇Δ^-1div(𝐮ξ)|^2=_fluid layer|∇Δ^-1J(ψ,ξ)|^2

where J(ψ,ξ)=∇^⊥ψ·∇ξ. Recall also
that the net-flux constraint requires 

    _fluid layerwξ=1.

As pointed out in sec:An-integral-formulation —see
the discussion surrounding eq:div-free —for
smooth enough designs (𝐮,ξ) the advection term cannot
vanish while the net-flux constraint and boundary conditions w|_∂Ω=ξ|_∂Ω=0
hold. As eq:integral_formulation makes clear, wall-to-wall
optimal transport is precisely about balancing these competing effects.
Regarding elasticity, we ask: what does it take for the membrane term
to nearly vanish? This can be answered with the aid of the lower bound

    _blistered
    region|e(𝐯)+1/2∇ϕ⊗∇ϕ|^2≳|_blistered
    region1/2∇ϕ⊗∇ϕ-λId_2×2|^2 +_blistered
    region|∇∇(ΔΔ)^-1∇∇ϕ|^2,

which is sharp for certain domains and boundary conditions.[Equality can be achieved if a suitable Airy stress function solving
ΔΔΨ=∇∇ϕ with appropriate boundary
conditions exists.] For the in-plane strain to nearly vanish, the bulk average of 1/2∇ϕ⊗∇ϕ
must be nearly constant and equal to a known multiple of the identity.
At the same time, ϕ must nearly satisfy the degenerate Monge-Ampère
equation ∇∇ϕ=0. It follows from the results of
<cit.> that these are incompatible constraints,
i.e., the membrane term cannot vanish while the bending term remains
finite. The situation is remarkably similar to that of wall-to-wall
optimal transport. 

The scaling law of the minimum energy for blistering is known. As
proved in <cit.>, there exist
constants C and C' depending only on Ω so that 

    Ch≤minimum elastic energy≤ C'h

for small enough h. The upper bound comes from a branching construction
involving finer and finer oscillations in ∇ϕ at a certain
lengthscale depending on the distance from the blister edge ∂Ω.
As opposed to the corresponding result for the wall-to-wall problem,
there is no logarithmic correction to scaling in eq:min-energy-est_FvK.
This can be explained with the help of eq:-1scaling and
eq:-2scaling: whereas the advection term has a -1
scaling in its quadratic nonlinearity J(ψ,ξ), the membrane
term has a -2 scaling in ∇∇ϕ and therefore
permits much stronger oscillations. As a result, branching can be
more easily accommodated in blistering than in optimal transport.
The lower bound from eq:min-energy-est_FvK asserts that
branching indeed achieves the minimum energy up to a prefactor depending
only on the domain. Its proof reminds of the proof of Howard's lower
bound given after eq:Howard_LB. For details we refer
the reader to <cit.> for the case where Ω is
a square with periodic boundary conditions at opposite sides, and
to <cit.> for the more general case of an arbitrary
domain Ω with suitably smooth boundary.





















































































































































































































































































































plain


