Large-scale regularity for the stationary Navier-Stokes equations over non-Lipschitz boundaries

In this paper we address the large-scale regularity theory for the stationary Navier-Stokes equations in highly oscillating bumpy John domains. These domains are very rough, possibly with fractals or cusps, at the microscopic scale, but are amenable to the mathematical analysis of the Navier-Stokes equations. We prove: (i) a large-scale Calder\'on-Zygmund estimate, (ii) a large-scale Lipschitz estimate, (iii) large-scale higher-order regularity estimates, namely, $C^{1,\gamma}$ and $C^{2,\gamma}$ estimates. These nice regularity results are inherited only at mesoscopic scales, and clearly fail in general at the microscopic scales. We emphasize that the large-scale $C^{1,\gamma}$ regularity is obtained by using first-order boundary layers constructed via a new argument. The large-scale $C^{2,\gamma}$ regularity relies on the construction of second-order boundary layers, which allows for certain boundary data with linear growth at spatial infinity. To the best of our knowledge, our work is the first to carry out such an analysis. In the wake of many works in quantitative homogenization, our results strongly advocate in favor of considering the boundary regularity of the solutions to fluid equations as a multiscale problem, with improved regularity at or above a certain scale.


INTRODUCTION
In this work we consider the large-scale boundary regularity for the stationary Navier-Stokes equations    −∆u ε + ∇p ε = −u ε · ∇u ε in B ε in a domain with a rough bumpy boundary. The no-slip boundary condition is prescribed only on the lower part Γ ε 1 of ∂B ε 1,+ . The boundary is rough in two aspects: (i) possible lack of regularity at the microscopic scale as the boundary may have fractals or inward cusps; (ii) bumpiness, i.e., the boundary is highly oscillating. The functions u ε = (u ε 1 (x), u ε 2 (x), u ε 3 (x)) ∈ R 3 and p ε = p ε (x) ∈ R denote respectively the velocity field and the pressure field of the fluid. The definitions of B ε r,+ and Γ ε r are given in Subsection 1.4. We will show large-scale regularity estimates, including Lipschitz estimate (see Theorem A in Subsection 1.1 below), C 1,γ estimate (see Theorem B) and C 2,γ estimate (see Theorem C). These improved regularity results at the large scales are generally false at the small scales due to the roughness of the boundary. The tools developed in this paper enable us to decouple the large-scale regularity from the small-scale properties of the boundary. Therefore, our results (i.e., Theorems A, B and C) show that stationary incompressible Newtonian fluids are regular above the microscopic scale, regardless of the irregularity of surfaces at the microscopic scale.
Before going into the details of our results and of the mathematical analysis, let us give some more general perspectives. The study of fluids over rough boundaries plays a prominent role in the field of hydrodynamics, at least for three reasons.
First, rough, bumpy or corrugated surfaces are ubiquitous in nature and engineering. They appear at any scales from geophysics (see for instance [69] for the fractal-like coremantle boundary in the Earth) to zoology [70] and microfluidics [80]. At the microstructure, the geometry may be anything from fractal to periodic and crenellated. No surface is perfectly smooth, and the lack of smoothness may actually enable us to resolve certain oddities, such as the no-collision paradox for a sphere dropped in a viscous fluid under the action of gravity [76,58,30,39,53]. Moreover, certain roughness patterns are either selected by biological processes and environmental pressure such as scales of sharks for their drag reduction properties, or designed for industrial applications especially in aeronautics, microfluidics and for the transport of fluids in pipes [70,31,65].
Second, the study of roughness is strongly tied to the derivation of boundary conditions in fluid mechanics. The question whether fluids slip or not over surfaces is still a matter of active debate. Experiments show that there is no universal answer and that the slip behavior depends a lot on the geometry and microstructure of the surface [14,64]. A widespread idea is that roughness favors slip. To give one specific example where finding the most accurate boundary condition is critical, let us cite the field of glaciology. The assessment of various friction laws for the flow of a glacier over a rough bedrock is crucial in order to understand the speed of glacier discharge and eventually estimate the sea level rise as a result of global warming [59,68].
Third, the study of the impact of roughness on the behavior of fluids accompanied the development of turbulence research, as underlined by Jiménez in [56]: Turbulent flows over rough walls have been studied since the early works of Hagen (1854) and Darcy (1857), who were concerned with pressure losses in water conduits. They have been important in the history of turbulence. Had those conduits not been fully rough, turbulence theory would probably have developed more slowly. The pressure loss in pipes only becomes independent of viscosity in the fully rough limit, and this independence was the original indication that something was amiss with laminar theory. Flows over smooth walls never become fully turbulent, and their theory is correspondingly harder.
Investigations of the effect of roughness on fluid flows span three distinct regimes. In the laminar regime, studies focus on the drag reducing properties of roughness elements [13,35]. As for the onset of turbulence [72,77], there are some indications that roughness lowers the critical Reynolds number for the transition from the laminar to turbulent regime [79]. In the fully turbulent regime, a similarity hypothesis for the flow over flat surfaces and for the flow over rough surfaces was put forward [78]. The extent to which such a universal law holds is still being disputed [56,23,34,71].
The three main directions raised above are reflected in the mathematical works. The literature is vast. Therefore we do not aim for exhaustivity.
First, there is an extensive body of works that deal with wall (or friction) laws, or in other words, effective or homogenized boundary conditions. One aims at replacing rough boundaries by fictitious, smooth or flat boundaries. In that line of research, it is well-known that Navier-slip boundary conditions provide refined approximations for fluids above bumpy boundaries. Under some quantitative ergodicity assumptions, one can get error estimates. Historically, periodic roughness profiles were first looked at [3,1,54,55]. Analysis of almost-periodic [41] or random stationary ergodic [37,12] boundary oscillations was done more recently. Let us also mention a few works that address non-stationary fluids [18,50], for which the analysis is less developed due to its inherent difficulties. We also point out that some authors attempted to justify boundary conditions arising in fluid mechanics starting from boundary conditions at the microscopic scale; see for instance [22,19,15] for the derivation of the no-slip boundary condition from a perfect slip boundary condition at the microscale, or [27] for the computation of the homogenized effect starting from Navier-slip boundary conditions at the microscale.
A second topic is the study of the effect of roughness on singular limits. The topics of rotating fluids and of the homogenized effect of bumpiness on Ekman pumping has been studied in numerous papers [36,38,29,28]. The paper [40] carries out an analysis of the vanishing viscosity limit in a specific scaling regime. There are also studies concerned with equations in singularly perturbed domains such as the Stokes equations in rough thin films [25] or water waves above a rough topography in the shallow regime [26].
Third, rough domains pose considerable numerical challenge. This aspect has certainly driven the development of wall laws in a model reduction perspective; see for instance [1,32]. Other approaches are being elaborated, such as Direct Numerical Simulations [21], Lattice Boltzmann Methods that are adapted to intricate geometries [79] and Large Eddy Simulations [4,16] that in this context cause important parametrization issues of the small scales.
In this work, we tackle these questions from the angle of regularity theory. The following two general objectives in regularity theory motivate our results: (i) identify building blocks describing the local behavior of solutions, (ii) estimate the decay of certain excess quantities at various scales. We prove that fluids above bumpy boundaries, that are very rough at the microscopic scale, have improved regularity at large scales. Our results are in the spirit of large-scale regularity estimates pioneered in [9] for periodic homogenization, and later extended to stochastic homogenization; see for instance [6,7,45,46] and [5] for the higherorder large-scale regularity theory. Our research program was started with the works [62,63] concerned with uniform regularity estimates above highly oscillating boundaries for elliptic equations. In the paper [51], the large-scale Lipschitz and C 1,γ estimates for the stationary Navier-Stokes equations were established above Lipschitz boundaries. A local Navier wall law was also obtained. Finally, let us also mention the paper [82] that deals with the large-scale regularity of elliptic equations above arbitrarily rough microstructures.
1.1. Outline of the main results of the paper. In this paper, we study the large-scale regularity for stationary incompressible viscous fluids modeled by the Stokes or Navier-Stokes equations, in domains that are very rough and bumpy at the microscale. Our results show that the large-scale regularity is completely independent of small-scale properties of the boundary.
Let us stress some novel aspects of our results. We refer to Subsection 1.2 for a further comparison with a few related works, and to Subsection 1.3 for an outline of the proofs.
First we consider John domains, whose boundaries allow for fractals and inward cusps. Hence, the boundaries considered in this paper get closer to the modeling of real boundaries found in nature, that in particular do not need to be graphs. John domains have in a broad sense the minimal properties for the analysis of incompressible fluids. Indeed, we rely on a Bogovskii operator in John domains to estimate the pressure. For precise definitions and a more complete discussion, we refer to Subsection 1.4 below.
Second, beyond the Lipschitz estimate, we prove higher-order C 1,γ and C 2,γ estimates for γ ∈ [0, 1), as stated in Theorems B and C. These require the construction of boundary layer correctors, which is at the heart of the paper in Section 4; see Subsection 4.2 for the first-order boundary layers and Subsection 4.3 for the second-order boundary layers. As far as we know, the present work is the first to construct the second-order boundary layers with a linear growth in the direction tangential to the boundary. To make the analysis more tractable, we assume that the boundary is periodic for the structure result of second-order boundary layers; see Theorem 4.3 and Theorem 4.4. We are aware of the papers [11,17], where a refined second-order approximation is constructed for the Stokes equations in a two-dimensional rough channel. However, the boundary layers considered in [11,17] only involve data spanned by linear and quadratic polynomials of the vertical variable, x 2 and x 2 2 in this two-dimensional case, which are bounded on the bumpy boundary. In our threedimensional situation, the class of "no-slip Stokes polynomials" (see Subsection 4.1) is much richer and involves boundary data with linear growth at spatial infinity.
Third, we provide explicit quantitative regularity estimates in the non-perturbative regime. Fourth, in the vein of the seminal works [9,10] and of [61,48], we provide pointwise estimates for the large-scale decay of the velocity and pressure parts of the Green function associated to the Stokes system in bumpy John half-spaces; see Subsection 1.3 and Appendix B. These estimates are pivotal to construct the first-order boundary layers in Subsection 4.2.
We now state the three main theorems of the paper.
Notice that Theorem A, as well as the subsequent results, holds in the non-perturbative regime for arbitrarily large M in (1.1). This is due to the energy subcritical nature of the stationary Navier-Stokes equations, which makes it an easier problem than the nonstationary Navier-Stokes system. Remark also that the powers of M in the right-hand side of (1.2) are explicit.
For higher-order C 1,γ and C 2,γ regularity results, we measure the oscillation of the solution with respect to modified polynomials that vanish on the bumpy boundary. These modified polynomials are polynomials of degree one and two that are corrected by the firstorder and second-order boundary layers.
Notice that by the Caccioppoli and Poincaré inequalities, the velocity estimate in (1.3) can also be replaced equivalently by a large-scale estimate of |∇u ε − ∇w(x/ε)|.
While Theorem B holds for arbitrary bumpy John half-spaces, for the next result, we work in periodic John domains. As we outlined above, the extra periodicity assumption makes the analysis of the second-order boundary layers more manageable.
We point out that the building blocks in Q 1 (Ω) and Q 2 (Ω) are defined through the firstorder and second-order boundary layers and play roles of correctors of Stokes system in the bumpy John domain Ω. It turns out that the above three regularity results, Theorems A, B and C, hold also for the linear Stokes equations, with a linear dependence on the size M of the solutions inḢ 1 (B ε 1,+ ). Therefore, these statements immediately imply the Liouville theorems for Stokes equations in bumpy John half-spaces with sublinear (see Corollary 3.1), subquadratic or subcubic growth (see Theorem 5.8).
1.2. Comparison to two closely related works. To further underline the novelty of our work, let us compare our results to the ones of two tightly linked papers.
In the paper [51], the first and second authors of the present work carried out the analysis of the large-scale Lipschitz and C 1,γ regularity for the stationary Navier-Stokes system. The results there, similar to Theorem A and Theorem B here, hold outside the perturbative regime, that is for arbitrarily large M in (1.1). The main differences between [51] and the present work are: (i) In [51] the bumpy boundary is given by a Lipschitz graph without structure, while here we work in bumpy John domains as defined in Definition 1.2 that are not necessarily graphs, without structure for the large-scale Lipschitz and C 1,γ regularity. (ii) In [51] the analysis relies on a compactness method originating from [9] and the first-order boundary layer correctors are needed to prove the large-scale Lipschitz estimate in Theorem A, while here we resort to a quantitative method, which enables us to by-pass the use of the first-order boundary layers for the large-scale Lipschitz regularity; see Subsection 1.3.
(iii) In [51] no analysis of the higher-order large-scale regularity is carried out, while here we build the second-order boundary layer correctors that make it possible to prove Theorem C. (iv) In [51], no pressure estimate is established, while in the present paper, we establish the pressure estimates in all cases, following the strategy developed recently in [49]. (v) In [51], the nonlinear estimates are not explicit, while here the dependence on M in (1.2), (1.3) and (1.4) is given as an explicit polynomial in M .
In the paper [82], the third author of the current paper carried out an analysis of the largescale Lipschitz regularity for linear elliptic equations in domains with arbitrary roughness at small scales and quantitative Reifenberg flatness at large scales. Hence, those domains are much rougher than the bumpy John domains considered here. We underline that the discrepancy in these assumptions on the domains comes from the fact that for incompressible Navier-Stokes equations, as opposed to elliptic equations, we have to estimate the pressure in terms of the velocity, which relies on a Bogovskii-type operator as in [51]; see Subsection 1.3 and Appendix A. To address this point we work in bumpy John domains defined by Definition 1.2.

1.3.
Outline of the strategy for the proofs. We now point to some essential ingredients and ideas for the proofs. In particular, the fact that we cannot rely on any smoothness at the microscopic scale requires several technical innovations.
Analysis in John domains. We perform the analysis in bumpy John domains, as defined in Definition 1.2. This type of domains is a good compromise between: • on the one hand a high level of arbitrariness of the boundary, which is not a graph, includes certain fractals or cusps, does not oscillate with any structure, and hence approaches better the properties genuinely rough physical surfaces found in real fluids, • and on the other hand the possibility of being amenable to mathematical analysis, considering the fact already underlined above that we work with incompressible fluid models that involve estimating the pressure, rather than elliptic equations which can be studied in even rougher domains.
In John domains, we can rely on the Bogovskii operator of [2], whose properties are summarized in Theorem A.1. This operator is required from the beginning of our analysis in Subsection 2.1 in order to prove the Caccioppoli inequality for the Stokes system, which then implies the reverse Hölder inequality (2.2), as a starting point of the large-scale regularity theory.
All the boundary estimates of this work are mesoscopic estimates in the sense that they involve averaged quantities smoothing out the possibly rough microscales. Although it is a direct consequence of the Caccioppoli inequality, notice that the reverse Hölder inequality (2.2) is a large-scale estimate. Indeed, going from the Caccioppoli inequality (A.3) to (2.2) uses the Poincaré inequality that holds in balls large enough, typically at a scale greater than ε. At scales smaller than ε, inward cusps of highly oscillating bumpy John domains may be seen, preventing Poincaré's inequality to hold.
In a nutshell: in the works [63,51], tools were developed, particularly for the analysis of the first-order boundary layer correctors, to handle bumpy domains with a boundary given by the graph of a Lipschitz function without structure. Here, the analysis in bumpy John domains requires to push the techniques even further, to the limit, as it seems, of what is technically possible. There is one particular point, where we are completely unable to transfer the techniques used above Lipschitz graphs to the present context. Indeed, in [63,51] we used a domain decomposition method pioneered in [41] to study the wellposedness of the Stokes system for the first-order boundary layer correctors. We do not manage to adapt this strategy, in particular the technique of local energy estimates in the bumpy channel, to our current situation. In this paper, we develop a different argument to construct the first-order boundary layers, based on the large-scale Lipschitz estimate proved as an a priori estimate. We discuss this intricate point in more details at the beginning of Subsection 4.2.
Quantitative method for the large-scale regularity. We rely on a quantitative method for large-scale regularity, inspired by the Schauder's theory pioneered by [6,8,73], the Calderón-Zygmund theory motivated by [20] and [74,Chapter 4] and the pressure estimate developed in [47,48,49]. This method is based on a perturbation argument. The principle of this method is the following: (1) approximate the original rough problem, by a smooth problem at any mesoscopic scales and obtain suboptimal quantitative estimates; (2) use the improved regularity of the approximate problem, to get the scale-by-scale decay of excess quantities (measuring for instance, Hölder continuity, Lipschitz, C 1,γ , C 2,γ , or higher regularity) for the original rough problem, up to a small error; (3) conclude by a real variable argument such as Theorem 2.6 or an iteration lemma such as Lemma 3.10, which are in some sense black boxes oblivious to the equations.
In the context of homogenization, the homogenized limit problem with constant coefficients is the approximate problem. Here, the approximate problem is a Stokes problem in a domain with a flat boundary. Both problems have improved regularity, in the sense that the solutions are basically as smooth as one wishes. We remark that from a high-level point of view all the regularity estimates in this paper follow the above scheme. For the large-scale W 1,p regularity stated in Theorem 2.5, item (1) above corresponds to Lemma 2.7, item (2) corresponds to the estimate (2.17) and item (3) corresponds to Theorem 2.6. For the proof of Lipschitz estimate in Theorem A, item (1) corresponds to Lemma 3.2, item (2) corresponds to Lemma 3.5 and item (3) corresponds to Lemma 3.10. The proofs of higher-order regularity estimates in Theorems B and C follow a similar scheme.
Relations between the results and structure of the proofs. We conclude this part by dissecting the logical connections between the main results and the intermediate results stated in the paper. For the Stokes system, the Caccioppoli inequality stated in Lemma A.3 and the Poincaré-Sobolev inequality imply the reverse Hölder inequality of Lemma 2.2. The Gehring lemma then implies an improvement of integrability stated in Lemma 2.4. This large-scale Meyerstype estimate is crucial to get the quantitative approximation result, Lemma 2.7, of the system in bumpy John domains by systems in flat domains, at any mesoscopic scales. The large-scale Calderón-Zygmund estimate for Stokes system stated in Theorem 2.5 is then obtained as a combination of Lemma 2.7 and a real-variable argument in Theorem 2.6.
For the Navier-Stokes equations, we use the large-scale Calderón-Zygmund estimate of Theorem 2.5 in combination with a large-scale Sobolev embedding stated in Theorem 2.8 to bootstrap the integrability of the nonlinear term. In the end, we are able to prove the decay of an excess quantity associated to the nonlinear term (see Theorem 2.9), which reads a largescale Hölder estimate for the nonlinearity. The large-scale Calderón-Zygmund estimate also enables us to get an improved approximation result; see Lemma 3.2. These estimates, together with the smoothness of the approximate problems in flat domains (see Lemma 3.5), imply an excess decay estimate in Lemma 3.7, which leads to Theorem A by an iteration lemma (i.e., Lemma 3.10).
The large-scale Lipschitz regularity in Theorem A then makes it possible to construct the velocity and pressure parts of the Green function in bumpy John domains, and to estimate its decay at large scales. This is the purpose of Appendix B, where we prove estimates for the velocity part of the Green function (see Proposition B.3), its derivatives (see Proposition B.4), and the pressure part of the Green function (see Proposition B.5). These estimates are the key for our new proof of the existence of the first-order boundary layer correctors; see Theorem 4.1. In this way we are able to by-pass the difficulties posed by the method used in [41,29,28,63,51].
To the best of our knowledge, the present work is the first to carry out a thorough analysis of the second-order boundary layer correctors, allowing for linear growth of the boundary data in the tangential direction. Our idea is to use the first-order boundary layer correctors in an Ansatz for the second-order boundary layers. In Theorem 4.3 and Theorem 4.4, we give the detailed structure of the second-order boundary layers. For our analysis to go through, we need some good quantitative convergence/decay of the first-order boundary layers away from the boundary. Hence we work in a periodic framework, according to Definition 1.3; but this is by no means an optimal assumption. Other structures, such as almost-periodic structures with a non-resonance condition, or random ergodic with quantitative decorrelation properties at large scales, would certainly be manageable.
The key outcome of Section 4 handling the construction of boundary layers is summarized in Proposition 4.6 and Proposition 4.7. They are then used in Section 5 to run the excess decay method for the higher-order regularity. Theorem B and Theorem C are proved there.

Notations and definitions.
John domains. We first define John domains. These domains were introduced by John in [57] and named after John in [66]. Our analysis takes advantage of a key property of John domains, namely the existence of a right inverse of the divergence operator. Such an operator is usually called a Bogovskii operator; see Appendix A where we state the result of [2].
Examples of John domains are: Lipschitz domains, NTA domains, domains with inward cusps or certain fractals such as Koch's snowflake. Notice that domains with outward cusps are not John domains. For our work, we generalize the above definition from bounded domains to a class of unbounded domains.
Let Ω be a domain containing the upper half-space of R 3 and assume that ∂Ω ⊂ {−1 < x 3 < 0}. We say that Ω is a bumpy John domain (or a bumpy John halfspace) with constants (L, K), if for any x ∈ {x 3 = 0} and any R ≥ 1, there exists a bounded John domain Ω R (x) with respect to x R = x + Re d and with constant L ∈ (0, ∞) according to Definition 1.1 such that where B R,+ (x) = Q R (x) ∩ Ω. Here Q R (x), defined later, is a cube centered at x with side length 2R.
The above definition guarantees that the constants of John domains are rescaling-and translation-invariant. This is a natural requirement as we are considering unbounded domains. Definition 1.3. We say that Ω is a periodic bumpy John domain if the following holds: (i) Ω is a John domain with constant (L, K), (ii) and Ω is (2πZ) 2 -translation invariant, namely 2πz + Ω = Ω for any z ∈ Z 2 × {0}.
For simplicity, we assume K = 2 in the whole paper. Otherwise, the constant in our main theorems will also depend on K.
Throughout the paper, we assume that Ω is a bumpy John domain satisfying Definition 1.2, or a periodic bumpy John domain satisfying Definition 1.3. We will always specify in case periodicity is needed. In fact, periodicity is used to construct the second-order boundary layer correctors in Subsection 4.3 and hence is also an assumption of Theorem C, Proposition 4.7, Subsection 5.2 and Theorem 5.8 (ii). Let We refer to Ω ε as a highly oscillating bumpy John domain. Note that (1.6) A key fact about Ω ε is that Ω ε is still a John domain with the same constants as in Definition 1.2, as these constants are scale-invariant.
Throughout the paper, we use the following notations: Since the boundary could be very rough at small scales, B ε r,+ and Γ ε r may have disconnected components. Fortunately, this will not cause any issue since the solutions will be extended naturally by zero across the boundary. We also define Q r = Q r (0) = (−r, r) 3 , Q r (y) = y + Q r (0), Q ε r = Q r ∩ {x ∈ R 3 |x 3 > −ε}, Q ε r (y) = Q r (y) ∩ {x ∈ R 3 |x 3 > −ε} and Q r,+ (y) = Q r (y) ∩ {x ∈ R 3 |x 3 > 0}. From the definition of B ε r,+ , one has B ε r,+ ⊂ Q ε r and |Q ε r \ B ε r | ≤ 4εr 2 .
Weak solutions. We work in the framework of weak solutions of (NS ε ). A velocity/pressure pair (u ε , p ε ) ∈ H 1 (B ε 1,+ ) 3 × L 2 (B ε 1,+ ) is said to be a weak solution to (NS ε ) if u ε satisfies: (i) ∇ · u ε = 0 in the sense of distributions, (ii) ψu ε ∈ H 1 0 (Q 1 ) 3 for any cut-off function ψ ∈ C ∞ 0 (Q 1 ), and (iii) the weak formulation for any ϕ ∈ C ∞ 0 (B ε 1,+ ) 3 . The Poincaré inequality is a fundamental tool in our paper. Since the weak solution vanishes on the lower boundary Γ ε 1 , we extend it to Q ε 1 by zero across Γ ε 1 . This enables us to use for instance [42,Proposition 3.15], to get that: for all fixed bumpy John domain Ω with constant L ∈ (0, ∞) according to Definition 1.2, for all fixed r ≥ ε, and for all u ∈ H 1 (B ε r,+ ) such that u = 0 on Γ ε r , where C is an absolute constant independent of ε and r. Notice that this estimate is only valid at scales r ≥ ε. Indeed, below that scale the constant in (1.8) may degenerate because in particular of inward cusps at small scales.
Other frequently used notations. The notation C denotes a positive constant that varies from line to line, and may or may not be universal. Whenever needed, we make precise what the constant depends on. The notation x · y stands for the inner product of vectors x, y ∈ R N . The notation a b (resp. a b) means that there exists a universal constant C such that a ≤ Cb (resp. Ca ≥ b). The notation a ≈ b stands for a b and a b.
1.5. Outline of the paper. Section 2 is devoted to the proof of the large-scale Calderón-Zygmund estimate stated in Theorem 2.5. We then use this result to bootstrap the regularity and obtain a large-scale Hölder estimate for the nonlinear term in the Navier-Stokes equations; see Theorem 2.9. In Section 3, we prove Theorem A. In Section 4 we construct the first-order and second-order boundary layer correctors. Theorem B and Theorem C are proved in Section 5. There are three appendices. Appendix A is devoted to the results related to Bogovskii's operator in John domains. Appendix B handles the construction and estimates for the Green function associated to the Stokes system in bumpy John domains. Appendix C provides a proof for the iteration Lemma 3.10.

ESTIMATES FOR THE NONLINEARITY
The goal of this section is to obtain some regularity estimates for the nonlinearity −u ε ⊗ u ε for the Navier-Stokes equations. As usual, this follows from a bootstrap argument for the stationary Navier-Stokes equations. However, since there is no smoothness up to the boundary, we have to carry out a delicate large-scale bootstrap argument. (2.1) We extend u ε and F ε by zero to the whole of Q 1 = Q 1 (0); they are denoted again by u ε ∈ H 1 (Q 1 ) 3 and F ε ∈ L 2 (Q 1 ) 3×3 respectively. Note that we also have ∇u ε = 0 in Q 1 (0) \ B ε 1,+ . For any r ≥ ε and Q 16r (y) ⊂ Q 1 (0), Lemma A.3 and the Sobolev-Poincaré inequality imply − Qr(y) Here the constant C depends only on L.
We refer to [82, Lemma 2.2] for a proof in the case of elliptic equations. The proof is exactly the same here, since it relies directly on the use of the Caccioppoli and Sobolev-Poincaré inequalities. Notice that this estimate holds only at large-scales, namely, r ≥ ε, because the Sobolev-Poincaré inequality fails for r ≪ ε as inward cusps are allowed in John domains and these cusps can be seen at a scale less then ε. Then we can show a large-scale Gehring's inequality (self-improving property); see the proof of Lemma 2.4 below.
For p ∈ [1, ∞), define the averaging operator The important exponents for us are p = 6 5 and p = 2. For convenience, sometimes we write M 2 t as M t in Subsection 2.2. The following lemma collects useful properties of M t .
Here the constant C depends on p and p ′ , but not on s, t, t 1 or t 2 .
Now, we rewrite the reverse Hölder inequality (2.2) by using the averaging operator.
and Ω be a bumpy John domain with constant L according to Definition 1.2. Let (u ε , p ε ) be a weak solution of (2.1). For any ε ≤ t ≤ r ≤ 1 and where C depends only on L.
Proof. By Lemma 2.1 (i) and (iii), we have − Qr(y) where we have used (2.2) in the second inequality.
To conclude the proof, we use a covering argument to adjust the size of cubes. By covering the cube Q 32r (y) by a finite number of cubes Q r (y i ), we get the estimate for ε ≤ t ≤ r ≤ 1 such as Q 96r (y) ⊂ Q 1 , at the price of a larger constant C than in (2.9). Now, if we replace 32r by r, then the desire estimate holds for 32ε ≤ 32t ≤ r and B 3r (y) ⊂ Q 1 . Finally, notice that for r 32 ≤ t ≤ r, the desire estimate (2.8) is trivial by using the properties in Lemma 2.1 only. This concludes the proof. Remark 2.3 (Covering argument). The covering argument above to adjust the size of cubes should be a standard technique in analysis. Similar arguments may be used later in this paper.
As a consequence, we have the following improved integrability, which can be read as a large-scale Meyers-type estimate.
Lemma 2.4. Let L ∈ (0, ∞) and Ω be a bumpy John domain with constant L according to Definition 1.2. There exists some p 0 ∈ (2, ∞) so that for any ε ≤ t ≤ r ≤ 1 and Proof. First of all, Lemma 2.2 and Gehring's inequality [43,Proposition 5.1] imply that there exists some p 0 ∈ (2, ∞) so that for any ε ≤ t ≤ r and Q 3r ( Note that Hölder's inequality and (2.2) imply where we also used Lemma 2.1 (iv) and (i) in the last two inequalities. This proves the desired estimate for r ≥ 16t ≥ 16ε. Finally, for r 16 ≤ t ≤ r, the desired estimate follows from the covering argument (see the proof of Lemma 2.2).
The following theorem is a large-scale boundary Calderón-Zygmund estimate, or in other words, a large-scale boundary W 1,p estimate, for the linear Stokes system. Theorem 2.5. For all ε ∈ (0, 1 2 ), L ∈ (0, ∞) and p ∈ (2, ∞) the following statement holds. Let Ω be a bumpy John domain with constant L according to Definition 1.2 where the constant C depends only on L and p.
The proof of Theorem 2.5 relies on a combination of: (i) a real variable argument (see Theorem 2.6), and (ii) the quantitative approximation at sufficiently large scales s of the solution u ε to the Stokes system in the bumpy domain by a solution to a Stokes problem in a flat domain (see Lemma 2.7).
We first state the real variable result. The following theorem is taken from [ There exists κ 0 > 0, depending on λ, p, q, c 1 and N , with the property that if 0 < κ < κ 0 , then F ∈ L p (Q 0 ) and where C depends on λ, p, q, c 1 and N .
We now turn to the approximation. Fix t ≥ ε. To apply Theorem 2.6, we introduce an approximation of u ε at all scales Let s ≥ t be fixed. By the co-area formula and the fact that ∇u ε ≡ 0 below the bottom boundary, there exists some (2.14) Note that t 0 depends particularly on the specific solution u ε . But this is harmless as t 0 is bounded uniformly in ε in [1,2] and C is an absolute constant. Now, we construct an approximation of u ε in Q t 0 s (y) by considering the following Stokes system we may extend the solution w s naturally across this boundary. For our purpose, we need some regularity estimates for w s . First of all, the energy estimate implies (2.16) Second, by the classical regularity theory for the Stokes system over a flat boundary, we have This yields, see [81,Lemma 3.3] and [60, Remark 9.3]. The above higher integrability of w s plays an important role in the following lemma.
Lemma 2.7. Let L ∈ (0, ∞) and Ω be a bumpy John domain with constant L according to Definition 1.2. Let (w s , q s ) be given as above. Then there exists σ ∈ (0, 1 12 ] such that for where C depends only on L, and C θ depends on L, σ and θ. Proof. We rely on the variational definition of the weak solutions of (2.15). First of all, by (2.15), we see that u ε − w s ∈ H 1 0 (Q ε t 0 s (y)) 3 and ∇ · (u ε − w s ) = 0, since u ε has been extended by zero. Thus we can test (2.15) against u ε − w s to obtain (2.21) for any P ∈ R (to be determined later). Combining (2.20) and (2.21) and using the fact (2.22) Now, we are going to estimate the integrals on the right-hand side of the above equation. Note that 1 − η 2 ε,+ and ∇η ε, To estimate the first integral, we use the Poincaré inequality applied in R ε s to obtain (2.23) The last integral of ∇(u ε − w s ) in the above estimate will eventually be absorbed by the left-hand side of (2.22). The main difficulty to proceed is to obtain a certain estimate of smallness for ∇u ε over the thin strip R ε s . This can be done by using Lemma 2.4. In fact, if θ ∈ (0, 1) and s ≥ ε/θ, Lemma 2.4 yields Since for any z ∈ Q 6s (y), It is important to notice that in the last inequality, C is independent of θ and C θ depends on θ. By a similar argument as in [82], we can now estimate the right-hand side of (2.23) as with some σ ∈ (0, 1 2 ). This is the desired estimate of ∇u ε in R ε s . Later on we will insert it into (2.23) and then (2.22) to reach a conclusion.
Let us turn to the estimate of the second integral on the right-hand side of (2.22). Using Hölder's inequality and the Poincaré inequality, we have Unlike the previous argument, we want to gain the smallness for (2.27) below from The estimate for ∇u ε over R ε s is given (2.24). On the other hand, by (2.18) and the Hölder's inequality, we have ˆR ε s
Proof of Theorem 2.5. We will first prove a slightly weaker version of (2.13) when Q 57r (x) ⊂ Q 1 (0). Then, (2.13) can be recovered thanks to a covering argument at the price of enlarging the constant by a numerical factor; see Remark 2.3 for more details. When Q 57r (x) is far away from the boundary Γ ε 1 , the estimate (2.13) is a consequence of interior regularity. Hence it suffices to prove (2.13) when Q 57r (x) ∩ Γ ε 1 = ∅. Note that this case can be reduced to the case when x ∈ {z 3 = 0} by a covering argument as well as interior regularity. To apply Theorem 2.6 to Q 0 := Q r (x), λ := 56 and F := then the well-known interior estimate for the Stokes system applies. If Q s (y) is contained entirely in {z 3 < −ε}, then trivially u ε ≡ 0 in Q s (y). Hence, it suffices to focus on the typical boundary case Q s (y) with y ∈ {z 3 = 0}. Moreover, we assume s < r/2 so that Now, for each Q s (y) with y ∈ {z 3 = 0}, we will discuss two cases. Case 1: s≥4t. By (2.17) and (2.19), there exists w s solving (2.15) and satisfying Note that the above estimate only holds for s ≥ ε/θ. Therefore, we will use Lemma 2.1 and replace ∇u ε and ∇w s by M 2 t [∇u ε ] and M 2 t [∇w s ], respectively. Precisely, the above two inequalities imply for s ≥ 4t ≥ ε/θ, This ends the study of the two cases. We now apply Theorem 2.6 with λ := 56, For any given p > 2, we may choose θ sufficiently small with Cθ σ < κ 0 so that the requirement of Theorem 2.6 is satisfied. Consequently, we arrive at 32) for all ε ∈ (0, θ(p)) and ε/(4θ) ≤ t ≤ r and Q 57r (x) ⊂ Q 1 (0). Estimate (2.13) now follows by a covering argument (see the proof of Lemma 2.2) and Lemma 2.1 (in order to adjust the size of balls and relax the condition t ≥ (4θ) to t ≥ ε). To remove the smallness condition ε ∈ (0, θ(p)), we observe that the case θ(p) ≤ ε ≤ t ≤ r ≤ 1 2 is trivial as the constant C is allowed to depend on p.

Bootstrap argument.
In this subsection, we apply the large-scale Calderón-Zygmund estimate proved previously to study the regularity of the stationary Navier-Stokes equations (NS ε ). Note that in Theorem 2.5, F ε is a general function. We will take advantage of the nonlinearity F ε = −u ε ⊗ u ε . As usual, the proof relies on a bootstrap argument.
Throughout this subsection, we set F ε = −u ε ⊗ u ε . To begin with, note that the Sobolev embedding theorem implies F ε ∈ L 3 , which yields M t [F ε ] ∈ L 3 . Hence, (2.13) holds with p = 3. To further improve the large-scale regularity, we need to lift the regularity of F ε from that of ∇u ε .
For any 0 ≤ a < b ≤ ∞, define a new maximal function |g|.
Note that M 1 (0,∞) is the usual Hardy-Littlewood maximal function. Clearly, by the L p boundedness of the Hardy-Littlewood maximal function, The following estimate is a sort of the large-scale Sobolev embedding theorem.
and Ω be a bumpy John domain with constant L according to Definition Then for any p > 3 and any q satisfying where the constant C depends only L, p and q.
Proof. Let p > 3 and q satisfy (2.33). Without loss of generality, we assume in addition We consider the cases x 3 ≥ t and x 3 < t separately. Assume first x 3 ≥ t and let N be the natural number so that 2 N −1 t < x 3 ≤ 2 N t. Note that u ε vanishes in a large portion of Q 2 N+1 t (x). By the triangle inequality and the Poincaré inequality, we have Using the definition of K q and M 1 (2t,5r) , we obtain On the other hand, if x 3 < t, then B 2t (x) has a relatively large portion not contained in Ω ε . Thus, the Sobolev-Poincaré inequality implies Using the same argument as (2.35), we see that M t [F ε ](x) has the same bound as (2.36) for x 3 < t.
Since by assumption, 1 2p < 1 q < 1 2p + 1 3 and 0 < α < min{ 1 3 , 1 q } is arbitrary, we may choose α so that 1 q > 1 2p + α. This implies p( 2 q − 2α) > 1. Thus, using the L boundedness of the Hardy-Littlewood maximal function, we obtain Now, observe that we may choose α < 1 q − 1 2p but sufficiently close to 1 for anyq > q, where we also used the fact that K m (r) ≤ K n (r) for any 1 ≤ m ≤ n. Finally, to recover the case with the exact exponent q, we may start with aq < q still satisfying 1 q < 1 2p + 1 3 . Then (2.37) holds for anyq >q, which includes the caseq = q. This proves the desired estimate. Now, a bootstrap argument between (2.13) and (2.34) In the following, we use this to prove a large-scale Hölder's estimate for F ε , which plays an important role in the Lipschitz estimate in the next section.
For every l > 3 and δ > 0 satisfying lδ < 6, we have where the constant C depends only on L, l and δ.
Proof. Note that, by a similar argument as at the end of the proof of Theorem 2.5, we only have to prove (2.38) when ε/N 0 ≤ t ≤ r ≤ 1/N 1 for some N 0 , N 1 ≥ 2. Let l > 3 and δ > 0 with lδ < 6 be given and fixed for the proof. First of all, by the Sobolev embedding theorem, Then, applying Theorem 2.8, we obtain that for any 3 ≤ p < ∞, Now, using Theorem 2.5 again combined with a covering argument, we derive from the last inequality that Now, let p > l. By the interpolation, we have For the given δ ∈ (0, 1), we want 4 − 2θ = 4 − 6/l + δ. This implies θ = 3/l − δ/2 and thus we may choose One can easily verify that θ ∈ (0, 1) by the assumption on l and δ. Consequently, we derive from (2.39) that Finally, we apply Theorem 2.8 to obtain for r ≤ 1/400, This ends the proof.
Note that if ∇u ε itself is in L p for p > 3, then Morrey's inequality implies that u ε is C 0,1−3/p , which implies, since u ε vanishes on the boundary, where C depends on M , L and p. Hence, (2.38) is consistent with the usual Morrey estimate.

LARGE-SCALE LIPSCHITZ ESTIMATE
In this section, we will establish the large-scale Lipschitz estimate of u ε and the oscillation estimate of p ε . We remark, for later use in Subsection 4.2, that Theorem A implies the following Liouville theorem for the Stokes system where Ω is a John domain in Definition 1.2. The proof of the following statement is standard.
3.1. Sep-up and approximation. First of all, we may write (NS ε ) as a linear Stokes system As in the classical regularity theory for Stokes system, we will use the large-scale C 0,α estimate of F ε in Theorem 2.9 to prove the large-scale Lipschitz estimate. The proof is based on the excess decay method.
Similarly to the large-scale Calderón-Zygmund estimate of Theorem 2.5, we also need to approximate the Stokes system (S ε ) at all scales greater than ε. Fix r ∈ [ε, 1 2 ] and let (v r , q r ) be the weak solution of the following Stokes system where we have automatically extended u ε across the bottom boundary by zero-extension and t 0 is a constant in the interval [1,2] chosen analogously as in (2.14) and (2.15). Note that (S r ) is a special case of (2.15) with s = r and y = 0, which means the estimates (2.16)-(2.18) hold also for (v r , q r ), in place of (w s , q s ). The following lemma is an analog of Lemma 2.7.
By examining the proof of Lemma 2.7, we obtain From Lemma 2.1 (iii) and Theorem 2.5 with p = 3, Next, we estimate the pressure by using the Bogovskii lemma. The issue is that, in general, B ε r,+ is not a John domain. By Definition 1.
Remark 3.3. The pressure estimate in John domains in the proof of Lemma 3.2 is a standard technique that we will frequently use throughout this paper. It allows us to transfer the pressure estimate to the estimates of ∇u ε and F ε . and (3.8) The quantity H can be dubbed as a zeroth-order excess quantity. In Section 5 we will consider higher-order excess quantities H 1st and H 2nd to address the large-scale C 1,γ and C 2,γ regularity. Moreover, for a pair of function (3.9) The following lemma states the comparability between H(v r , q r ; θr) and H(v r , q r ; θr).
(ii) Let P * ∈ P 1 be such that By the triangle inequality, we see that By the definition of P * and the Poincaré inequality, we have Consequently, Hence, we obtain (3.11) from (3.15) combined with (3.14) and (3.16). This completes the proof of Lemma 3.4. where C depends only on L and α.
Proof. By the regularity of the Stokes equations in flat domains, v r ∈ C 1,α (Q ε r/2 ) , q r ∈ C 0,α (Q ε r/2 ). Let e 3 = (0, 0, 1). The boundary C 1,α estimate of v r on {x 3 = −ε} implies The C 0,α estimate of q r implies Then by the Bogovskii lemma in a Lipschitz domain Q ε r/2 , we have Hence on the one hand, by (3.18), (3.19) and the Caccioppoli inequality in rectangular region Q ε r , we see that On the other hand, since v r (x) − P (x + εe 3 ), for any P ∈ P 1 , is a weak solution to (S r ) with the same pressure q r , we may apply (3.20) to this solution and obtain By the triangle inequality and the definition of H, this leads to, for any P ∈ P 1 , H(v r , q r ; θr) ≤ H(v r − P, q r ; θr) where we also used the observation P (x + εe 3 ) = P (x) + ε(∇P )e 3 . In particular, we may choose P = P * that minimizes Then, it is clear that Then the estimate (3.17) follows from (3.21) with P = P * , (3.22) and the Poincaré inequality. The proof is complete.

23)
where C depends only on L and α.

24)
where C depends only on L and α.
Proof. The triangle inequality and Lemma 3.6 imply where in the last line the energy estimate of (S r ) is applied. By the definition of H, we find (3.26) The Poincaré inequality and Lemma 3.2 imply , we obtain the desired estimate (3.24) by the definition of Φ in (3.8). This completes the proof.
3.3. Iteration. In the following two lemmas, we prove some properties of H and Φ needed when iterating (3.24).
Lemma 3.8. Let L ∈ (0, ∞) and Ω be a bumpy John domain with constant L according to Definition 1.2. Let (u ε , p ε ) be a weak solution of (S ε ). There exists a function h(r) defined on [ε, 1 2 ] such that (3.30) Here C depends only on L. Notice that the function h depends on u ε .
Proof. The proof is similar to [49, Lemma 6.1] and hence we provide the outline of the proof. Let P r ∈ P 1 be such that Then the inequality (3.28) follows from and (3.29) is trivial by definition. For (3.30), we observe that for any r 1 , r 2 ∈ [r, 2r] This completes the proof of Lemma 3.8.

31)
where C depends only on L.

35)
where the constant C depends only on C 0 , α, β and θ.
Proof of Theorem A: In the following proof, we actually only need to show (1.2) for the case N 0 ε ≤ r ≤ 1/N 1 for some N 0 , N 1 ≥ 2. The case 1/2 ≥ r ≥ 1/N 1 follows trivially by enlarging the size of the cube and a standard pressure estimate (see Remark 3.3); the case ε ≤ r ≤ N 0 ε follows from the case r = N 0 ε. From the previous lemmas, we can choose N 0 = 4 and N 1 = 16. Hence, we may assume without loss of generality that r ∈ [4ε, 1 16 ]. We apply Lemma 3.10 to H(r) = H(u ε , p ε ; r) and Φ(r) = Φ(u ε , p ε ; r). Choose θ sufficiently small so that we have Cθ α ≤ 1/2 in (3.24) in Lemma 3.7. We need to verify the conditions in Lemma 3.10. Note that (3.34b) is obvious and (3.34d)-(3.34f) follow from Lemma 3.8. To verify (3.34a) from Lemma 3.7 and verify (3.34c) from Lemma 3.9 (with ε replaced by 4ε), it suffices to note that Theorem 2.9 implies for any β ∈ (0, 2), δ ∈ (0, 1) with β + δ < 2 and r ∈ [ε. 1 2 ]. Hence, we may apply Lemma 3.10 with B 0 = C(M + M 4+2β+2δ ) to obtain where in the last inequality, we have used the Poincaré inequality and a standard pressure estimate (see Remark 3.3) to bound Φ u ε , p ε ; 1 2 ) by C(M + M 2 ). By the Caccioppoli inequality, (3.36), (3.37) and Lemma 2.1 (iii), for r ∈ [4ε, 1 16 ], which proves the desired estimate of the velocity u ε . Next, we give an estimate for the pressure. For r ∈ [ε, 1 4 ], we observe that Using the technique as in (3.6) and by the Bogovskii lemma, the desired estimate of ∇u ε just proved and (3.36), we have On the other hand, let N ∈ N be such that 2 N r ∈ [ 1 32 , 1 16 ]. Then Now, observe that for each j = 0, 1, · · · , N − 1, Finally, by the same trick as in (3.6), we obtain Summarizing up the above estimates, we obtain the desired estimate for the pressure p ε . This completes the proof of Theorem A. ✷
Let P (x) = (P 1 (x), P 2 (x), P 3 (x)) be the no-slip Stokes polynomials of degree 2 of u at 0. First of all, since u = 0 on ∂R 3 + , then we must have

(4.2)
The linear part is familiar. So let us concentrate on the quadratic part. Note that there are no terms x 2 1 , x 2 2 , or x 1 x 2 , because u = 0 on the boundary. If there is no further restriction on u, then there are 9 free variables b ij , 1 ≤ i, j ≤ 3, as shown in (4.2). If ∇ · u = 0 in Q 1/2,+ (0), then we claim that P is also divergence-free. If this claim is true, then we must have b 11 + b 22 + 2b 33 = 0, b 31 = b 32 = 0.
Because of this restriction on the coefficients, the dimension for the homogeneous no-slip Stokes polynomials of degree 2 becomes 6. We can find basis polynomials as follows: Note that these polynomials are solutions to the stationary Stokes system with associated pressure L (2j) given by and L 2j (x) = 2x 3 , for j = 5, 6. Now, let us show the claim that P is divergence-free. Since u = P + O(|x| 3 ), we have that ∇ · u = ∇ · P + O(|x| 2 ) = 0 in {x 3 ≥ 0} ∩ B 1/2 (0). Because of ∇ · P = C 0 + C 1 · x for some C 0 ∈ R and C 1 ∈ R 3 , we see that C 0 + C 1 · x = O(|x| 2 ). Hence we must have C 0 = 0 and C 1 = 0; otherwise, it is easy to find a contradiction by taking x = δC 1 or −δC 1 for sufficiently small δ.
Similarly to the linear solution pairs (P (1j) , 0), the fundamental fact about the polynomial pairs constructed above is that (P (2j) , L (2j) ) are quadratic solutions of Stokes equations in the upper half-space To study the C 1,γ and C 2,γ regularity of (NS ε ), the linear and quadratic solutions of Stokes equations in R 3 + are not enough. We need to construct linear and quadratic solutions in Ω which vanish on ∂Ω, where Ω is a bumpy John half-space in the sense of Definition 1.2. These solutions will be constructed based on (P (1j) , 0) and (P (2j) , L (2j) ). Observe that P (ij) does not vanish on ∂Ω. Therefore we have to introduce new correctors, called boundary layers, in order to correct the boundary discrepancy on ∂Ω. Precisely, we will show the existence of weak solutions (with corresponding sublinear or subquadratic growth) of the following boundary layer equations Here a couple (v, q) ∈ H 1 loc (Ω) 3 ×L 2 loc (Ω) is said to be a weak solution of (BL , and (iii) the weak formulation:

First-order boundary layers. We consider the first-order boundary layer equations
for j ∈ {1, 2}. The solvability of (BL (j) 1st ) follows from the next statement.
where the constant C depends only on L.
In the work [51] the well-posedness of the system (BL (j) 1st ) was proved over Lipschitz graphs by a domain decomposition method: coupling of the Stokes problem in a bumpy channel Ω∩{x 3 < 0} with the Stokes problem in the flat half-space {x 3 > 0} via a nonlocal Dirichlet to Neumann boundary condition at the interface {x 3 = 0}. We face considerable technical difficulties when trying to adapt this strategy to the case of bumpy John domains. Indeed, the local energy estimates in the bumpy channel require to estimate the pressure, or to work with divergence-free test functions. In either case, we need to construct a Bogovskii operator for a sequence of exhausting domains containing Ω ∩ {|x ′ | ≤ k, x 3 < 0} with a constant uniform in k. The construction of the Bogovskii operator of Theorem A.1 by [2] relies on connecting any point in the bumpy John domain to a fixed neighborhood of a reference pointx. Such a procedure gives, for a slim domain such as Ω∩{|x ′ | ≤ k, x 3 < 0}, a constant in the estimate (A.1) that scales proportionately to the horizontal size k of the domain. We are unable to take advantage of the small vertical extent of the domain to provide a modified construction of the Bogovskii operator. This would be needed to carry out the downward iteration on the local energy estimates, also called Saint-Venant estimates, in [51].
Here we take advantage of the fact that we already proved large-scale Lipschitz estimates by the quantitative method, without relying on boundary layers as in [51]. Therefore, we develop a new strategy using the large-scale Lipschitz estimate to prove the existence of solutions to (BL We also define Ω <N and Ω >N in a similar manner. Proof of Theorem 4.1. We denote (v (1j) , q (1j) ) by (v, q), not to burden the notation.   Notice that F is a bounded function supported in a slim channel S := {x ∈ R 3 | 3 ≤ x 3 ≤ 4}. Thus, the problem is reduced to finding a weak solution of (4.9) satisfying sup ξ∈Z 2ˆΩ ∩(ξ+(0,1) 2 )×R We rely on the representation of w and q by the Green kernel Thanks to the properties of the Green function (G, Π), it suffices to prove that ∇w and q are well-defined and satisfy the estimate (4.10).
In the following proof, we take the zero-extension of (G, Π) as is done in Appendix B.

4.3.
Second-order periodic boundary layers. Let P (2j) be a no-slip Stokes polynomial of degree 2. We consider the second-order boundary layer equations Constructing solutions to (BL  Throughout this subsection, we assume Ω is a periodic bumpy John domain according to Definition 1.3. Consider the fundamental periodic domain We regard Ω p as a submanifold of T 2 ×R, where T = R/2πZ is the flat torus. By definition, Ω p is open and connected in T 2 ×R. Moreover, Ω p ∩{x 3 < 2} is diffeomorphic to a bounded John domain in R 3 . We thus have the Bogovskii operator on Ω p ∩ {x 3 < 2}. It is important to notice that there is a one-to-one correspondence between the functions in Ω p and the (2πZ) 2 -periodic functions in Ω. We say a function F defined in Ω is (2πZ) 2 -periodic if F (x) = F (x + z) for any x ∈ Ω and z ∈ (2πZ) 2 × {0}. In other words, if f ∈ L 2 loc (Ω p ), then there exists a locally L 2 function F defined in Ω so that F (x) = f (x), wherex is the representation in Ω p so that x −x ∈ (2πZ) 2 × {0}. In this sense and for convenience, we do not distinguish between F and f . Denote by L 2 (Ω p ) and H 1 0 (Ω p ) the closure of C ∞ 0 (Ω p ) under the norms Let H 1 0,σ (Ω p ) be the subspace of H 1 0 (Ω p ) 3 that consists of all the divergence-free functions, namely, . We now recall the Fourier series representation for the solutions of (BL (j) 1st ) on the flat half-space {x 3 > 0}. The same formulas are obtained in [51, Proposition 3] based on the periodic Poisson kernel. Note that [51] uses the fact that the equations are imposed on a domain whose boundary is given by the graph, but a similar proof is valid if we utilize the zero extension of the functions.
and moreover,v (4.25) Here C is a universal constant.
Construction of v (2j) for j ∈ {1, 3, 5, 6}. We construct the second-order boundary layers v (2j) corresponding to P (2j) = P (2j) (x) for j ∈ {1, 3, 5, 6}. These boundary layers are solutions to (BL (j) 2nd ) with subquadratic growth; see Theorem 4.3. We begin with the case j = 1, where P (21) (x) = x 2 x 3 e 3 . We recall that (v (21) , q (21) The difficulty in the analysis of (BL (1) 2nd ) is that the boundary value is not periodic and has linear growth as x 2 → ∞. We aim at eliminating the growth factor x 2 and at recovering the periodic structure. The key finding is the connection between the first-order and secondorder boundary layers on the boundary, namely v (21) − x 2 v (11) = 0, x ∈ ∂Ω. (4.26) This observation is the basis of the Ansatz for v (21) . Recall that v (11) converges exponentially fast to the constant α (11) ∈ R 3 , when x 3 → ∞ by the spectral gap near frequency 0 yielded by the periodicity; see (4.25). Hence the non-decaying divergence ∇ · (x 2 v (11) (x)) = v (11) 2 (x) can be corrected by adding a corrector −α Here η + (·) is a function on R satisfying η + (t) is smooth and non-negative, Below, we also need the cut-off η − defined in (4.8).
The following statement gives the existence and the structure of second-order boundary layers with subquadratic growth.  (21) , q (21) where (R (21) , Q (21) where the constant C depends only on L.
Proof. We aim at proving the existence of (R (21) , Q (21) ) and estimating it so that (v (21) , q (21) ) defined by the right-hand sides in (4.28) gives a weak solution of (BL (1) 2nd ). (Existence) By the previous discussion, we begin with a formal examination of x 2 v (11) for all x ∈ Ω p . Notice that the expression above simplifies for x 3 > 1: Then, by Proposition 4.2, we get v (11) This means that x 2 v (11) (x) − α (11) 2 x 3 η + (x 3 )e 3 is not divergence-free. Thus, our next goal is to construct a function to correct the divergence for x 3 > 1. Define which is an element of H 1 (Ω p,>0 ) 3 where Ω p,>0 := Ω p ∩ {x 3 > 0}. Of course, there is no unique way to construct a right-inverse of the divergence such as d. We may extend d(x) to the whole domain Ω p by multiplying it by η + (x 3 ) and still correct the divergence To check the divergence condition, we calculate Obviously, D is supported in Ω p,<2 := Ω p ∩ {x 3 < 2}, in which we can rely on the Bogovskii operator to find a right-inverse of the divergence. Let A :=´Ω p,<2 D. Let χ + (x 3 ) be a smooth cut-off function such that (4.31)

Hence, by Appendix A, there is a Bogovskii corrector
and B H 1 (Ω p,<2 ) ≤ C, where C depends only on the John constant L of Ω. We extend B by zero to the whole domain Ω p and denote it again by B ∈ H 1 0 (Ω p ) 3 . Let us combine the above correctors and define Note that C ∈ H 1 0 (Ω p ) 3 . In particular, ∇C L 2 (Ωp) ≤ C, where C depends only on the John constant L of Ω. By (4.25), the function C converges exponentially fast to −(2π) −2 A as x 3 → ∞, and its derivatives decay exponentially fast to 0 as x 3 → ∞.
Construction of v (22) and v (24) . The boundary layers corresponding to P (22) and P (24) can be constructed by using the Green function. The fact that P (22) and P (24) only depend on the vertical variable x 3 and that there is no growth in the tangential variable x ′ makes the analysis much easier than for P (2j) , j ∈ {1, 3, 5, 6} studied above. The proof of the following proposition is almost identical to the one of Theorem 4.1. Notice that here we state Theorem 4.5 in the periodic case only for convenience. Indeed we use these correctors in combination with (v (2j) , q (2j) ) for j ∈ {1, 3, 5, 6} whose existence is stated in Theorems 4.3 and 4.4 in periodic bumpy John domains. However, the existence of (v (2j) , q (2j) ) for j ∈ {2, 4} can be proved in general bumpy John domains according to Definition 1.2.
where the constant C depends only on L.

4.4.
Estimates of boundary layers. Before closing this section, we summarize the estimates of the boundary layers. The following propositions can be proved in a similar manner as in [51,Lemma 4] combined with a direct computation. The details are omitted here.
where C depends only on L.
where C depends only on L.
Hence, Q 1 (Ω) (resp. Q 2 (Ω)) is the vector space that contains all the "linear" (resp. "quadratic") solutions of the Stokes system in Ω vanishing on the boundary ∂Ω; see the Liouvilletype results at the end of this section stated in Theorem 5.8.
Remark 5.1. Note that the pressure part in estimate (1.3) of Theorem B is different from the Lipschitz estimate in whichP is − B ε 1/2,+ p ε . Actually, in (1.3),P 1 is the average of the corrected pressure over a small ball, i.e., for some (w, π) ∈ Q 1 (Ω); see (5.15). This is reasonable since we are concerned with the C 0,γ estimate of the pressure andP 1 plays a role similar to the zeroth-order term in the Taylor expansion of the pressure, if the boundary is flat. We emphasize thatP 1 depends on ε. The point here is thatP 1 is independent of r.
The critical fact we are going to use is that any (w, π) ∈ Q 1 (Ω) is a solution of the Stokes system in Ω that vanishes on ∂Ω. Hence, by rescaling, (u ε , p ε ) − (εw(x/ε), π(x/ε)) is still a weak solution with a no-slip boundary condition. This observation allows us to capture the regularity beyond the Lipschitz estimate. To this end, we define the first-order excess by Recall that (w, π) ∈ Q 1 (Ω) means that for some ℓ 1 , ℓ 2 ∈ R. We will also use the quantity Φ(u ε , p ε ; ρ) defined in (3.8).

4)
where C depends only on L.
Indeed, since r 0 is comparable to 1, the above estimate follows directly from the Poincaré inequality and Bogovskii's lemma. The proof is complete.
The above theorem directly implies the C 1,γ estimate for the velocity. To handle the pressure estimate in Theorem B, we need the following lemma.

5.2.
Large-scale C 2,γ estimate over periodic boundaries. The goal of this subsection is to prove the large-scale C 2,γ regularity stated in Theorem C. In this subsection, we assume Ω is a periodic bumpy John domain defined in Definition 1.3. The argument for C 2,γ estimate is similar to the C 1,γ estimate. Throughout, we assume (w 1 , π 1 ) ∈ Q 1 (Ω) and (w 2 , q 2 ) ∈ Q 2 (Ω). In other words, for some ℓ 1k , ℓ 2j ∈ R.
Proof. The proof follows from the strategy developed in Section 3, in particular from Lemma 3.2 to Lemma 3.7. Let (v r , q r ) be the solution of the approximate problem (S r ). We will first use the C 2,1 estimate of v r = (v r,1 , v r,2 , v r,3 ) at the lower boundary x 3 = −ε.
Precisely, in view of no-slip Stokes polynomials defined in Section 4.1, the It follows from (5.18), (5.19) and the triangle inequality that Now, for any θ ∈ (0, 1 4 ] and r > ε, integrating (5.20) in Q ε θr , we obtain Next, to see the oscillation estimate for the pressure, observe that is a solution of the Stokes system in Q ε r with a no-slip condition on x 3 = −ε. Applying Bogovskii's lemma to q * r and the Caccioppoli inequality to v * r (combined with (5.21)) in Lipschitz domains, we have Notice that L (2j) are linear functions. Thus, an application of (5.19) This, combined with (5.21), gives . This is the key second-order excess estimate we need for (v r , q r ) in Q ε r . To proceed, we follow the similar argument developed in Section 3. Precisely, using an analogue of Lemma 3.4, we may replace Q ε r by B ε r,+ , with an error bounded by Then, by taking the approximation estimate in Lemma 3.2, we can further replace (v r , q r ) by (u ε , p ε ) with the error . Therefore, combined with the energy estimate for (S r ), we now have Next, we insert the boundary layers into the above inequality. By (5.8) and which follows from Proposition 4.7, we obtain from (5.24) and (5.19) along with the energy estimate for (S r ) that In view of the definition of H 2nd , we arrive at H 2nd (u ε , p ε ; θr) Finally, by the Caccioppoli inequality in B ε 5r,+ , we obtain the desired estimate.
The following lemma is parallel to Lemma 5.4.
Proof of Theorem C: The estimate for the velocity is contained in (5.26). The estimate for pressure can be derived similarly as Theorem B. The details are left to the reader. ✷ 5.3. Liouville-type results. As an application of the construction of boundary layers and uniform regularity, a Liouville-type theorem for Stokes systems can be shown by the largescale Lipschitz, C 1,γ and C 2,γ estimates. We point out that our large-scale regularity results hold also for the linear Stokes system, although with linear dependence on M in the righthand sides of (1.2), (1.3) and (1.4). The proofs are simpler, using that the source term F ε = 0. To describe the Liouville-type theorem, consider the Stokes system in the entire Ω   where Ω is a bumpy John domain according to Definition 1.2. Let B R = B R (0). We state the Liouville-type theorem as follows. Its proof follows from a routine rescaling of the large-scale regularity estimates. Notice that this result complements Corollary 3.1 already stated above. If for some σ ∈ (0, 1), up to a constant for p).
then there exists a function p ∈ L 2 (Ω) unique up to a constant for which we havê

Namely, the pair (u, p) is a weak solution of the Stokes equations.
A direct application of Bogovskii's operator is the Caccioppoli inequality for the Stokes equations in a John domain: and Ω be a bumpy John domain with constant L according to Definition 1.2. Let ε ∈ (0, 1 2 ] and F ε ∈ L 2 (B ε 1,+ ) 3×3 , and let u ε ∈ H 1 (B ε 1,+ ) 3 be a weak solution to (A.2). Then we have for all r ∈ [ε, 1 2 ], ∇u ε 2 where the constant C depends only on L. In particular it is independent of ε and r.
where the constant C depends only on L. In particular it is independent of ε and r.
Remark A.4. Remark that in the right-hand side of (A.3) the norms of u ε and F ε are taken on the enlarged ball B ε 2r,+ , rather than on B ε r,+ . This is in contrast with the Caccioppoli inequality in a bumpy Lipschitz domain, see for instance [51,Lemma 18]. The reason for this modification is that by Definition 1.2 of the bumpy John domain Ω, in order to estimate the pressure, we use the Bogovskii operator in a John domain Ω ε r satisfying B ε r,+ ⊂ Ω ε r ⊂ B ε 2r,+ when r ≥ ε.

APPENDIX B. LARGE-SCALE ESTIMATES FOR THE GREEN FUNCTION
This appendix is devoted to the study of the Green function for the Stokes equations in a bumpy John half-space according to Definition 1.2. The large-scale estimates proved in Section 3 will be applied. The basic scheme is to derive estimates for the velocity part of the Green function directly from the interior and large-scale boundary Lipschitz estimates. For this we follow the strategy pioneered in [9,10]. Then, we deduce the estimates for the pressure part of the Green function from Bogovskii's lemma and the estimates for the velocity part.
We use B R (x) = Q R (x) to denote the cube centered at x with side length 2R. If the center is not important in the context, it is abbreviated as B R . Throughout this appendix, Ω ≤N , Ω ≥N , Ω <N , and Ω >N defined around (4.7) will be used. Moreover, letx denote the projection of x ∈ R 3 on ∂R is denoted by Y 1,2 0,σ (D). Note that, when the Lebesgue measure of D is finite, we have Y 1,2 0 (D) = H 1,2 0 (D) by the Sobolev inequality u L 6 (D) ≤ C ∇u L 2 (D) if u ∈ C ∞ 0 (D). Moreover, we see that Y 1,2 0,σ (D) as well as Y 1,2 0 (D) 3 is a Hilbert space with an inner product u, v =´D ∇u · ∇v.

(B.6)
Then, by using Lemma A.2 on each bounded John subdomain, one sees that there is a pressure p ∈ L 2 loc (Ω) for which we have (B.4), uniquely determined under the condition p(x) → 0 as x 3 → ∞.
When constructing the pressure component Π(x, y) in (B.5), we need a careful analysis since the domain is unbounded unlike in [24]. Here the oscillation estimate of p will play a crucial role. For an open set E, define the oscillation of p in E by (B.7) The following lemma shows a fundamental oscillation estimate for the pressure. where C is a universal constant. (ii) Let z ∈ ∂R 3 + and let R > 2. If −∆u + ∇p = 0 and ∇ · u = 0 in Ω ∩ B R (z) and u = 0 on ∂Ω ∩ B R (z), then where C depends on L and is independent of z and R.
Proof. The interior case (i) is classical and the proof is omitted. Let us prove the boundary case (ii). Since only the case where R is sufficiently large is nontrivial, we assume that R > 32. For any x ∈ Ω >2 ∩ B R/2 (z), the mean value property of harmonic functions yields Here we assume r = x 3 2 ≤ R 16 , hence r > 1. Note that if x 3 > R 8 , the oscillation can be handled by the interior estimate (B.8).
Recall thatx is the projection of x on ∂R 3 + . Let Ω 4r (x) be a John domain given by Definition 1.2 so that Ω ∩ B 4r (x) ⊂ Ω 4r (x) ⊂ Ω ∩ B 8r (x). Clearly B r (x) ⊂ Ω 4r (x), B R/4 (x) ⊂ B 3R/4 (z) and B R/2 (x) ⊂ B R (z). By (B.10) and the Bogovskii lemma, where we also used the Lipschitz estimate of u in the fourth inequality. Similarly, we have On the other hand, by the pressure estimate for the Stokes system (an analogue of Theorem A with linear dependence on M ), Finally, combining the above estimates, we arrive at (B.12) Since x ∈ Ω >2 ∩ B R/2 (z) is arbitrary, this implies the estimate (B.9) with 2R in the righthand side instead of R. Then a covering argument using (B.8) yields the desired estimate (B.9). The proof is complete.
Remark B.2. The interior oscillation estimate holds also for −∆u + ∇p = f and ∇ · u = 0 in B R provided f ∈ L q (B R ) 3 for some q > 3. Precisely, by the classical Schauder theory, where C is a universal constant. Now, we are ready to construct Π(x, y) and prove Properties (i)-(iii) of (G, Π). For a given f ∈ C ∞ 0 (B R (0) ∩ Ω) 3 with R > 32, we consider the Stokes equations (B.4) with u given by (B.5) and the associated pressure p ∈ L 2 loc (Ω). For x ∈ Ω >2 such that |x| ≥ 4R,

This further implies
osc (B.14) This shows that p(x) converges to a constant as x → ∞. By the assumption that p(x) → 0 as x 3 → ∞, we know the limiting constant is zero. Hence, in view of (B.14), we derive for all x ∈ Ω >2 satisfying |x| ≥ 4R. Moreover, by arguing in a similar manner as in Step 3 in the proof of Theorem 4.1 and using (B.15) instead of (4.15), we find that for sufficiently large R ′ ≥ R, with a constant C(R ′ ) depending on R ′ . On the other hand, for x with either x ∈ Ω ≤2 or |x| ≥ 4R, we can connect x to another pointx ∈ Ω >2 with |x| ≥ 4R by a chain of a finite number of cubes {B r i (z i ) | i = 1, 2, · · · , N } such that B 2r i (z i ) ⊂ Ω. Using Remark B.2 on each B r i (z i ), as well as (B.15) applied tox, we see that for any x, R is a constant depending only on q, x, and R. From (B.15) and (B.17), for each fixed x ∈ Ω, the map f → p(x) is a bounded linear functional on L q (B R (0) ∩ Ω) 3 . By the Riesz representation theorem, there is a unique function Π(x, ·) ∈ L q ′ (B R (0) ∩ Ω) 3 with q ′ ∈ [1, 3 2 ), so that Note that the above Π(x, ·) is only defined in B R (0) ∩ Ω for a fixed x. As x and R vary, we can obtain a family of such functions, which can be glued together by the uniqueness of p. Thus we have constructed a function Π(x, y) defined in the entire Ω × Ω satisfying Π(x, ·) ∈ L q ′ loc (Ω) 3 . To investigate the local integrability of Π(·, ·), let us fix R > 1 and define a functional S(f, g) for smooth f, g supported in B R (0) ∩ Ω by From (B.16), by taking a sufficiently large R ′ ≥ R, we see that Now we can prove that (G, Π) satisfies Properties (i)-(iii). Property (iii) is obvious from the arguments so far. Property (ii) follows from Property (iii) combined with the Lebesgue differentiation theorem. Here we use the fact that, for all φ ∈ C ∞ 0 (Ω) 3 , the function of ŷ belongs to L 2 loc (Ω) 3 because of (B.18). The integrability of Π(·, y) in Property (i) follows from the weak form (B.3). Consequently, we have constructed the Green function (G, Π) meeting Properties (i)-(iii).
We should point out that in the above argument for existence, the estimate, for example of Π(x, ·), is very rough, especially when x is close to the boundary ∂Ω. This is because the large-scale regularity of Π(x, ·) is not taken into consideration. In the following, we obtain some more careful estimates of (G, Π) by studying the equation (B.3).

B.2.
Large-scale estimates of the velocity component. For convenience, let G(x, y) and Π(x, y) be zero-extended for both x and y. Recall the symmetry G(x, y) = G t (y, x), where G t is the transpose of G. Thus by definition, G(x, y) = 0 if either x ∈ ∂Ω or y ∈ ∂Ω and x = y. Denote by δ(x) the distance from x to ∂Ω.
Notice that ∇ x G denotes the derivative of G with respect to the first variable, i.e.
Similarly, ∇ y G denotes the derivative of G with respect to the second variable.

(B.21)
Here C depends on L.
Notice that G and Π are zero-extended outside Ω. Therefore, the integrals above make sense even in the case when B 1 (x) or B 1 (y) intersect Ω c . For the estimates concerned with the oscillation of the pressure, on the contrary, we make precise when the balls intersect the boundary; see for instance Lemma B.1.
Case 1: We first consider the typical case y / ∈ B 10r (x). Since we are working on cubes, it is more convenient to define R := |x − y| ∞ = max 1≤i≤3 |x i − y i |, which is comparable to the usual distance |x − y|. In this case, we only need to prove where we have used the Poincaré inequality, the large-scale Lipschitz estimate in the last inequality, and the property that for y / ∈ B 10r (x), B R/4 (x) ∩ B R/2 (y) = ∅. Using the energy estimate for (B.23), we have Our goal is to estimate the left-hand side of (B.22). We would like to use the large-scale Lipschitz estimate applied to G(x, ·) in B R/2 (y). Recall that we assume δ(y) ≤ δ(x) = r ≤ R 10 . If δ(y) ≤ 2, the Poincaré inequality implies − B 1 (y) |G(x, z)| 2 dz |G(x, z)| 2 dz |G(x, z)| 2 dz 1/2 ≤ C r(1 + δ(y)) R 3 ≤ C δ(x)(1 + δ(y)) |x − y| 3 , as desired. Case 2: y ∈ B 10r (x) and δ(y) < r 20 . In this case, we need to slightly modify the previous argument. Note that r = δ(x) is comparable to R = |x − y| ∞ since R 10 ≤ r ≤ 2R. As before, let f ∈ L 2 (B R/2 (y) ∩ Ω) 3 and (u, p) solve (B.23). Since B R/2 (y) ∩ B r/4 (x) = ∅, by the interior estimate, the Poincaré inequality and the energy estimate, |u(x)| ≤ C − |G(x, z)| 2 dz This proves Case 2. Case 3: y ∈ B 10r (x) and δ(y) ≥ r 20 . This is the interior case, and we actually have the pointwise estimate as desired.
Analogously, we can also show the estimates for the gradient of G. |Π(x, z)| 2 dz (B.37) Here C depends on L.
Proof. We will carry out a delicate oscillation estimate of the pressure originating from [48]. We first consider the estimate (i), i.e., x 3 > 2 and y 3 > 2. Consider a point w ∈ Ω with w = y. Let t = |w − y| ∞ . We claim with C independent of t, w, and y. The operator osc is defined in (B.7). We prove the above claim by considering different situations. If w ∈ B y 3 /2 (y), then t < Next, if w / ∈ B y 3 /2 (y), we consider two subcases: (a) |w 3 | < t 4 ; (b) |w 3 | ≥ t 4 . Without loss of generality, we assume t > y 3 2 > 20. For the case (a), letŵ be the projection of w on ∂R 3 + . Using the interior and boundary pressure estimates in John domains from Lemma B.1 combined with a covering argument, osc B t/4 (w)∩Ω >2 Π(·, y) ≤ C − where we have also used (B.29) and (B.33) in the second inequality. Now, for the case (b), B t/4 (w) may be decomposed as a union of a finite number of cubes B t/16 (w i ) with i = 1, 2, · · · , K 0 where K 0 is an absolute constant, so that B t/8 (w i ) is contained in Ω >2 . Thus, Π(x, y) = Π(y).
This convergence is uniform on any compact set in {y 3 > 2}. We show that Π(y) ≡ 0. In fact, if f ∈ C ∞ 0 (Ω) 3 , the pressure of the Stokes equations with the source f is given by p(x) =ˆΩ Π(x, y) · f (y) dy.