Non Uniqueness of Power-Law Flows

We apply the technique of convex integration to obtain non-uniqueness and existence results for power-law fluids, in dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\ge 3$$\end{document}d≥3. For the power index q below the compactness threshold, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \in (1, \frac{2d}{d+2})$$\end{document}q∈(1,2dd+2), we show ill-posedness of Leray–Hopf solutions. For a wider class of indices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \in (1, \frac{3d+2}{d+2})$$\end{document}q∈(1,3d+2d+2) we show ill-posedness of distributional (non-Leray–Hopf) solutions, extending the seminal paper of Buckmaster & Vicol [10]. In this wider class we also construct non-unique solutions for every datum in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}L2.


Introduction
This paper studies non-uniqueness and existence of solutions of the following model of non-Newtonian flows in d dimensions, d ≥ 3 where the velocity field v and and the pressureπ are the unknowns, Dv = 1 2 (∇v +∇ T v), and the non-Newtonian tensor A is given by the following power law for some ν 0 , ν 1 ≥ 0 and q ∈ (1, ∞). A natural energy associated with the system (1) is Let us consider a distributional solution v to (1), (2) with spatial mean zero, on a ddimensional flat torus. The formula (3) together with A(Q)Q ∼ |Q| q explains why v ∈ L ∞ (L 2 ) ∩ L q (W 1,q ) is called an energy solution. If such solution satisfies additionally the energy inequality e(t) ≤ e(0) (t-a.e.), then it is called a Leray-Hopf solution.
For the problem (1) we show two non-uniqueness and one existence result. In short: (A) In the regime 1 < q < 2d/(d + 2): There are non-unique Leray-Hopf solutions.
Our results are sharp concerning the power-law index q. The regime 1 < q < (3d + 2)/(d + 2) includes the case of the incompressible Navier-Stokes equation in d ≥ 3. The precise formulations can be found in Sect. 1.3.

Background of power-law flows. Model (1) with a slightly different choice of A(Q), namely
with q ≥ 2 was introduced to wide mathematical community by Ladyzhenskaya at her 1966 Moscow ICM speech; her formula (30) in [24] corresponds exactly to (1), (4). With q = 2, both models (1), (2) and (1), (4) reduce to the (incompressible) Navier-Stokes equations. The Ladyzhenskaya's choice: (4) with q ≥ 2 and our (2) with q ≥ 2 are analytically equivalent. In particular, the non-Newtonian tensor A(Q) is in both cases nonsingular at Q = 0, and distributional solutions are well-defined for velocity fields in the class v ∈ L 2 loc , Dv ∈ L q loc .
Ladyzhenskaya's rationale for analysing (1) was twofold: on the one hand, relaxation q ≥ 2 helps to avoid the traps of the Navier-Stokes case q = 2. At the same time, the choice of power-laws for the tensor A is both consistent with first principles of continuum mechanics and widely used in applications. Let us elaborate on each of these points.
The model (1) with power-law for A of type (2) or (4) agrees with the constitutive relations for incompressible, viscous fluids. Recall that in deriving the Navier-Stokes equation one restricts the admissible relations between the Cauchy stress tensor T and D (dictated by the material frame indifference) by the Ansatz of linear dependence between T and D (i.e. by the Stokes law), cf. [22]. The power law model relaxes this Ansatz, but remains well within the frame indifference principle.
Of course studying an arbitrary model that is merely consistent with the first principles may be applicationally void. This is not the case of (1) however. The power-laws have been proposed independently in 1920's by Norton [33] in metallurgy and by de Waele [19] and Ostwald [34] in polymer chemistry. The related timeline can be found in section 1 of [36]. For details, the interested reader may consult also the monographs [3,28,37,38] and the recent survey [4] with its references. Just in order to fix the hydrodynamical intuition, let us observe that q < 2 in (1) models the case when the fluid is more viscous (roughly, 'solid-like') for small shears ('external forces') and less viscous ('liquid-like') for large shears e.g. ice pack, ketchup, emulsion paints, hair gel, whereas q > 2 means reverse behavior e.g. cornstarch-water solution, silicone-based solutions.
Let us note that, despite the mathematical interest in q ≥ 2 in context of gaining regularity compared to Navier-Stokes equations, the 'shear-thinning' case q ≤ 2 appears to be more meaningful for applications, where models of type (1) with ν 0 > 0 appear as Bird-Carreau-Yasuda models (or called by a subset of those names). In particular, experimental fits for the threshold value 6/5 and above can be found on p. 174 of [3]. Furthermore, even parameter choices well-into our Leray-Hopf non-uniqueness regime are suggested, cf. p.18 of [39]. (In both [3] and [39] n = q − 1, d = 3. A discrepancy between appearing there a and our model is insignificant for our results.) From the applicational perspective, our result may be seen as invalidating certain choices of parameters and data.

Essential analytical results for power-law fluids.
Consider the system (1), (2). For q > 2d d+2 the space W 1,q of system's energy embeds compactly into L 2 loc of the convective term div (v ⊗ v). Hence one may expect an existence proof of Leray-Hopf solutions via compactness methods. Indeed, a relevant statement can be found in [20], which is itself the final step in a chain of attempts of many authors, including Frehse and Nečas with collaborators [21,29] to improve the lower bound on q. To be precise, the energy inequality e(t) ≤ e(0) is not stated explicitly in [20]; however it can be proven e.g. along the lines of proof of Theorem 3.3 of [4].
What is known about existence and uniqueness of solutions to (1) can be thus sketched as follows 1.3. Our contribution. The short version of our results presented at the very beginning of the paper, recast graphically to facilitate comparison with Fig. 1, reads Observe that Fig. 2 complements Fig. 1 sharply with respect to q.
Let us now present the detailed statements of our results. We always consider system (1) on the d-dimensional flat torus T d , with v having its spatial mean zero, and we define the notion of distributional solution in the usual way.  loc (T d × (0, 1)) with Dv ∈ L max{1,q−1} loc (T d × (0, 1)) is a distributional solution to (1), (2) on the space-time domain for any divergence-free ϕ ∈ C 1 (T d × [0, 1]) vanishing at t = 0 and t = 1, and any ψ ∈ C 1 (T d ).

Non-uniqueness in the Leray-Hopf class
Our first theorem and its corollary show that below the compactness exponent, i.e. for q < 2d d+2 , multiple Leray-Hopf solutions may emanate from the same L 2 initial data. In fact, we produce solutions v ∈ C(L 2 ) ∩ C(W 1,q ) with quite arbitrary pre-determined profile e of the (total) energy (3).
Analysing the proof of Theorem A one realises that choosing an infinite family of nonincreasing energy profiles {e α } α∈A with a common C 1 bound, one can produce infinitely many distinct Leray-Hopf solutions with the same initial datum.

Non-uniqueness of distributional solutions
If we drop the ambition to control the energy and require only to pre-determine the profile of the kinetic part of the energy T d |v| 2 (t), then we produce non-unique solutions for exponents below the scalingcritical one, i.e. for q < 3d+2 d+2 . Moreover, they enjoy the regularity v ∈ C(L 2 ) ∩ C(W 1,r ) for any r < 2d d+2 . This is our second result.
Moreover, fix 0 ≤ T 1 < T ≤ 1 and two energy profiles e 1 , e 2 as above, such that e 1 (t) = e 2 (t) for t ∈ [0, T ]. There exists v 1 , v 2 satisfying 1), 2) and such that v In particular, choosing T 1 = 0, T = 1/2 and e 1 , e 2 to be as above, non-increasing and e 1 ≡ e 2 , the corresponding v 1 , v 2 are two distinct distributional solutions, which belong to C(L 2 ) ∩ C(W 1,r ), dissipate the kinetic energy, and share the same initial datum.

Existence of multiple solutions for any L 2 data
In Theorems A, B the initial data are attained strongly (in particular we can add initial values to the distributional formulas for solutions, extending test functions to non-vanishing ones at t = 0), but they are constructed in the convex integration scheme, thus possibly non-generic. This issue is addressed in our third theorem. It shows existence of energy solution emanating from any solenoidal vector field in L 2 , for power laws below the scaling exponent.
Remark 1. Theorems A, B, C hold for Ladyzhenskaya's choice (4) instead of our choice (2), in the case q ≥ 2. On the contrary, for q < 2, they do not hold for (4), since (4) contains also linear dissipation (see Sect. 1.1). However, Ladyzhenskaya's interest was limited to q ≥ 2 (as a vehicle to mitigate Navier-Stokes difficulties), and in this range the theorems hold both for (2) and for (4), whereas, for q < 2, our choice (2) is more common in applications than (4).

Differences between our non-uniqueness and existence results.
The non-uniqueness Theorems A, B, focus on possibly strongest notions of solutions: they allow, respectively, for full-or kinetic energy inequality and strong attainment of a (constructed) initial datum, but they do not produce non-unique solutions for any initial datum. Conversely, Theorem C provides existence of many weak solutions for an arbitrary solenoidal initial datum in L 2 . In particular, this is the first existence proof for the case of q ≤ 2d d+2 . The obtained solutions are, however, much weaker than that of Theorems A, B : they do not allow for any kind of energy inequality (in fact, even their kinetic energies are in a sense pathologically large) and the initial datum is attained merely in a weak sense.
For the Euler equations a result similar to our Theorem C (existence of distributional solutions for any initial datum, with unnaturally high energy) was first proven in [40]. The difference with our Theorem C is that, in the case of the Euler equations (no dissipation), the problem of controlling higher order derivatives is completely absent.

The 3d
Navier-Stokes case. Theorems B, C cover also the case of non-unique weak solutions of three-dimensional Navier-Stokes equations, first proven in [10]. Our Theorem B shows that ∇v ∈ L 6/5− . This probably holds for solutions constructed in [10] as well, though the best regularity claimed there is curl v ∈ L 1 . Theorem C produces infinitely many weak solutions for any divergence-free datum in L 2 (but with unnaturally high energies).
We stay close to the concentration-oscillation method developed for the transport equation in [31,32], and localised to avoid dimension loss in [30], see also [5].
The basic picture of the construction, as in any convex integration scheme applied to the equations of fluid dynamics, is the following. Given an exact flow (v, π ), i.e. a solution to (1), one tries to distinguish the good ('laminar', 'averaged') component of v, i.e. v and the remainder, thought to be responsible for turbulence (interestingly, the case q = 3 in (2), where scaling (6) fails, is the Smagorinsky model for turbulence). A typical averaging process · does not commute with nonlinear quantities, thus applying · to (1) yields for u = v Above, u is a well-behaved flow and the Reynolds stress R encodes the difference between u = v and the exact v itself. The rough idea behind producing non-unique solutions to (1) is to reverse-engineer the above picture. We can thus consider the following relaxation of (1) Assume we have identity (9) with certain (u 0 , π 0 , R 0 ). It is easy to find at least one smooth solution of (10), since R 0 is at our disposal. If one can produce another u 1 , q 1 such that (u 1 , q 1 , R 1 ) solves (9) and R 1 is strictly smaller than R 0 , there is a hope to iteratively diminish the Reynolds part R n to 0 with n → ∞. Consequently, in the limit one produces an exact solution v, π. Non uniqueness in the above procedure may be specified in at least two ways: • either by enforcing v to be equal at some times, say for t ∈ [0, 1/3], to a given regular solution v 1 and for t ∈ [2/3, 1] to another regular solution v 2 , as for instance in [6] or [31]. • or by specifying a kinetic energy profile, see e.g. [10,18], or the present work.

Organisation of proofs.
In Sect. 2, we state the main proposition of the paper, i.e. Proposition 1, which contains the inductive step described above, from (u 0 , π 0 , R 0 ) to (u 1 , π 1 , R 1 ), with R 1 "much smaller" than R 0 . Section 3 gathers preliminary material. In Sect. 4 we introduce a generalisation of Mikado flows that serves as a building block for u 1 given u 0 . Next, in Sect. 5, assuming a solution (u 0 , π 0 , R 0 ) to (9) is given, we define (u 1 , π 1 , R 1 ). Estimates for (u 1 − u 0 ) and R 1 occupy Sect. 6. Section 7 concludes the proof of the main Proposition 1. Having it in hand, we prove Theorem A in Sect. 8. The proofs of Theorems B-C follow similar lines and therefore are only sketched in Sects. 9-10.
1.8. Notation. We use mostly standard notation, e.g. T d denotes the d-dimensional torus [0, 1] d ,Ẇ 1,q is a homogenous Sobolev space, C ∞ 0 (T d ; B) are smooth functions with mean zero, domain T d and values in set B (the target set will be sometimes omitted). We take N = {1, 2, . . . }.
We suppress the variables and the spatial domain of integration, if no confusion arises. We use | · | . instead of · . for norms. For L p -norms on the torus T d , we will abbreviate | · | L p (T d ) to | · | L p or even to | · | p . In other cases, e.g. when taking the L p -norm on R d , we will explicitly write the underlying domain, where the norm is calculated, e.g. | · | L p (R d ) . The finite-dimensional norm is | · |. The projection onto null-mean functions is We will call d × d (symmetric) matrices (symmetric) tensor. For a tensor T , we denote its traceless part byT := T − 1 d tr (T )Id. The space of symmetric tensors will be denoted by S, its open subset of positive definite tensors by S + . If R is a symmetric tensor, div R is the usual row-wise divergence. We use two types of constants M's, which are uniform over iterations, and C's which are not (both possibly with subscripts), for details see Sect. 6.1. All constants may vary between lines.
Further notation is introduced locally when needed.

Main Proposition: An Iteration Step
Recall that S is the space of symmetric tensors.

Definition 2.
A solution to the Non-Newtonian-Reynolds system is a triple (u, π, R) where with spatial null-mean u, π, satisfying in the sense of distributions.

Remark 2.
Despite smoothness of u, we can not require that (10) is satisfied in the classical sense or π, R are smooth (in space), because of non-smoothness of A(Du).

Remark 3. (R vsR)
Use of the trace-free Reynolds stress simplifies computation, in particular proof of the energy iterate Proposition 14. The difference between R toR is facilitated by the ambiguity of pressure: (u, π, R) solves (10) ⇐⇒ (u, π − 1 d tr R,R) solves (10).
As observed in the introduction, the crucial point in the convex integration scheme is, given (u 0 , q 0 , R 0 ), to produce an appropriate correction (u 1 , q 1 , R 1 ) which decreases R i , improves the energy gap, and retains as much regularity as possible. This single iteration step is given by Proposition 1. Let ν 0 , ν 1 ≥ 0 and q < 2d d+2 be fixed. Fix an arbitrary e ∈ C ∞ ([0, 1]; [ 1 2 , 1]). There exist a constant M such that the following holds. Let (u 0 , π 0 , R 0 ) be a solution to the Non-Newtonian-Reynolds system (10), as in Definition 2. Let us choose any δ, η, ∈ (0, 1]. Assume that and Then, there is another solution (u 1 , π 1 , R 1 ) to (10) (as in Definition 2) such that

Control of A.
We collect the needed growth estimates for A(Q) and for A(Q)Q.

Lemma 1 (Growth estimates for A). Let
The proof is standard. For convenience of the reader, we added it in Appendix.

Nash-type decomposition.
Let us denote the set of positive-definite d ×d symmetric tensors by S + . We recall Lemma 2.4 in [17] Lemma 2. For any compact set N ⊂ S + there exists a finite set K ⊂ Z d and smooth functions k : N → [0, 1], such that any R ∈ N has the following representation: 3.3. The role of oscillations. The convex integration paradigm is to use fast oscillations of corrector functions (correcting u i to u i+1 in our case, roughly speaking) to inductively diminish error terms (in our case Reynolds stresses R i ). Thus for a function f and λ ∈ N let us define Observe that f λ has the same L p norms as f since we work on T d , and a factor λ appears for each derivative, i.e.

It holds
. Then for any r ∈ [1, ∞] Proof. The case r = ∞ follows the proof of Lemma 2.6 in [31]. For the case r < ∞, since v is null-mean, let us solve the Laplace equation div ∇h = v and define G := ∇h. It holds (div G) λ = λ −1 div (G λ ) and thus, integrating by parts and using Hölder The Sobolev embedding for the null-mean G yields |G| r ≤ C|∇G| L min(r ,d+1) = C|∇ 2 h| L min(r ,d+1) . This is controlled thanks to Calderón-Zygmund theory by |v| L min(r ,d+1) .
Even when the l.h.s. of (17) is replaced with T d |av λ |, the decorrelation between frequencies of a and v λ allows to improve the generic Hölder inequality to (for the proof cf. Lemma 2.1 of [31]): Proposition 3 (Improved Hölder). Let f, g be smooth maps on T d . Let r ∈ [1, ∞].

Antidivergence operators.
We provide now various inverse divergence operators, needed for construction of R 1 in Proposition 1, with appropriate estimates. The purpose of the bilinear inverse divergences below is to extract oscillations of one function, say g λ , out of the product f g λ . The last of them, R 2 N , is an operator with symmetric tensor values, such that div div R 2 N f = f for every null-mean real function f ; it facilitates construction of the R lin term of R 1 , cf. (60). (i) (div −1 : symmetric antidivergence) There exists div −1 : and for the fast oscillating u λ (ii) (R N : improved symmetric bilinear antidivergence) For any N ≥ 1 there exists a bilinear operator R N : (iii) (R N : improved symmetric bilinear antidivergence on tensors) For any N ≥ 1 there exists a bilinear operatorR N : (iv) (R 2 N : improved symmetric bilinear double antidivergence) For any N ≥ 1 there exists a bilinear operator The proof is standard, cf. [18,30,31] and can be found in Appendix.
For the operator defined in (iv), we use the notation R 2 N to denote that this operator acts as a double antidivergence. It does not coincide in general with R N • R N .
The above bilinear antidivergences may be thought of as approxima- where the gap between R N and R ∞ closes as N → ∞, similarly for R 2 N and R 2 ∞ .

Mikado Flows
In this section we introduce the building blocks of our construction, namely the concentrated localized traveling Mikado flows. This section is essentially a rearrangement of known material: Mikado flows were first introduced in [17] in the framework of the Euler equations; the concentrated Mikado flows (or fields) were defined in [31], inspired by the intermittent Beltrami waves of [10]; the traveling version of concentrated Mikados we present here is close in spirit to the intermittent jets introduced for the first time in [6], and resemble the construction in [30].
The original Mikado flows of [17] are fast oscillating pressureless stationary solutions to Euler equations having the form where k ∈ Z d is a direction. For a finite set of directions K (given by the decomposition Lemma 2) one can choose functions k ∈ C ∞ 0 (T d , R) so that the following holds for any k, k ∈ K and λ ∈ N Satisfying property (i) is equivalent to choosing k so that ∇ k · k ≡ 0, then also (ii) follows. Having (iii) is a normalisation of − T d ( k ) 2 . Disjointness of supports (iv) is ensured in d ≥ 3 via an appropriate choice of an anchor point ζ k for the cylinder B ρ (0) + {ζ k + sk} s∈R + Z d (which is the periodisation of the cylinder B ρ (ζ k ) + {sk} s∈R with radius ρ and axis being the line passing through ζ k with direction k). Such choice is possible in view of the first part of the following lemma.

Lemma 3. (Disjoint periodic tubes and blobs)
There exist {ζ k } k∈K ⊆ R d and ρ > 0 such that for all k, k ∈ K , k = k : for every s ∈ R.
Notice the difference between (27) and (28): in (27) there are two free parameters s, s , whereas in (28) there is only one free parameter s. Part (ii) of the statement is not strictly needed for our proof (because we are always assuming d ≥ 3). We think however that Part (ii) could be of interest in view of future applications to two dimensional problems, and thus we decided to include it in the statement. Proof of the Lemma can be found in Appendix.
The convex integration approach uses the properties (i)-(iv) to diminish a given Reynolds stress R 0 of a given solution (u 0 , π 0 , R 0 ) to the Non-Newtonian-Reynolds system (10) by correcting u 0 roughly as follows. Thanks to (iii), we can decompose R 0 via the Nash Lemma 2 into Let us add to u 0 the correctorũ = k (R 0 ) k λ k. Recall the notation Since the term P =0 k λ k ⊗ k λ k is λ-periodic and null mean, applying R ∞ of (24) to the r.h.s. above yields R 1 of order λ −1 , such that div R 1 = div (ũ ⊗ũ − R 0 ). So picking λ large, i.e. letting k λ oscillate fast, allows to deal with the error R 0 . The property (i) allows to control divũ.

Concentrated Mikado flows. Since
fast oscillations, in general, blow up derivatives of the correctorũ. Thus controlling Sobolev norms of velocity fields appearing over convex integration steps seems problematic. This issue may be circumvented by a concentration mechanism, introduced in [31] and critically inspired by [10]. Let us briefly explain it. For n ≤ d, take a compactly supported smooth function f : R n → R, rescale it to f μ (x) = μ a f (μx), μ ≥ 1, and periodize without renaming to f μ : T d → R. This is concentrating and results in This procedure yields the 'concentrated Mikado' k μ k satisfying Having now an interplay between λ and μ one can expect to control certain Sobolev norms by choosing a, p appropriately. However, to preserve the properties (i) and (ii) of (26), i.e. ∇ f μ · k ≡ 0 (or in other words: k μ k being the Euler flow), the underlying function f μ cannot depend on the direction k. It means that the underlying real function is not compactly supported in R d , but at best in R d−1 . Thus at best n = d − 1, but then (26). This leads to the choice a = d−1 2 above and consequently to Scaling (30) would force us to prove our results with d substituted by d − 1, so that e.g. q could vary only in the interval 1 < q < 2(d−1) (d−1)+2 (which in particular requires d ≥ 4, as in [27]).
Summing up, the concentrated Mikado k μ k satisfies properties (26) (i)-(iv) of the original Mikado, but has unsatisfactory scaling (v).

Concentrated localized Mikado flows.
A natural idea to deal with the 'loss of dimension' in (30) is to localise the Ansatz (25). Let us thus take a smooth radial cutoff function φ and define : R d → R via (x) = φ(|x|). We want to retain gains stemming from concentrating, and since now : R d → R, while of (25) allowed merely for d − 1 concentrations, it is better to concentrate in , thus producing μ . We periodize this function without renaming it and allow to oscillate at an independent frequency λ 1 . Hence our new Ansatz reads Let us now state and prove a result gathering needed properties of the cutoff μ,λ 1 .
Let ρ > 0 and {ζ k } k∈K be given by Lemma 3. We now set for every k ∈ K and Property (32) now follows from the scaling properties of μ , whereas (33) from (27), and assumed μ ≥ ρ −1 .

Concentrated localized traveling Mikado flow.
Unsurprisingly, introducing ddimensional cutoff destroys the properties (i)-(iii) of standard Mikados. The most severe loss, due to its critical scaling, is not having (ii) anymore. A crucial idea how to handle this issue, introduced in [10], is to let the cutoff function travel in time along l k with speed ω. This leads to a corrector term Y k (see below), whose time derivative compensates lack of (ii). At the same time Y k is of order 1 ω , so it can be controlled by choosing ω large.
The concentrated localized traveling Mikado flow is our final Ansatz. It will be denoted by W k , but it is important to bear in mind that it is determined by the parameters The next proposition concerns our final Mikado flows W k and Mikado correctors Y k . Proposition 5. Let K ⊂ Z d be a fixed finite set of directions. Let k λ 2 be the function used to produce the standard Mikado (25) with its properties (i)-(iv). Let μ,λ 1 be the localisation provided by Lemma 4. Define the functions W k : the functions W k , Y k are spatially λ 1 -periodic are have the following properties: Proof. The spatial λ 1 -periodicity of W k , Y k follows from the assumption λ 2 /λ 1 ∈ N.
Since W k , Y k are obtained from stationary functions by means of a Galilean shift, for (36) and (37) it suffices to estimate the respective stationary functions.
with the second inequality due to (25) and (32). Since by assumption λ 1 μ < λ 2 , we obtain (36) estimate for W k . A similar computation yields the estimate for Y k . For (37) we compute with the equality valid because the normalisation of (32) holds. Since the normalisation (iii) of the standard Mikado k implies that ( k λ 2 ) 2 − 1 is null-mean, and it oscillates at the frequency λ 2 , by Proposition 2 we have (32). We reached (37). Disjointness of supports of follows from (33): Assume that for some ( thus contradicting (33). The Mikado functions have the form Remark 7. Let us compare our W k with the concentrated Mikado.
(a) W k is not divergence free, i.e. (i) does not hold. Furthermore W k is now time dependent. (b) W k does not satisfy (ii). There appears Mikado corrector Y k to compensates this deficiency, see (38). This means however that the new term Y k must be appropriately estimated. (c) Property (iii) holds approximately, see (37). (d) Supports are pairwise disjoint (now in space-time). (e) The scaling with a dimension loss (30) is now improved to (36), which is our main gain.
Proposition 5 yields the following estimates in relation to (b), (c), (e): In order to deal with (a), let us recall the heuristics of an ideal antidivergence operator R ∞ of Remark 6. A short computation involving (34), (36), and (32) yields Therefore, seeking smallness of will motivate the choice of the relations between the parameters μ, λ 1 , λ 2 , ω in Sect. 7.
Notice that we did not add the term λ 2 (40), as it is the product of the first and the third term in (40), and thus its smallness is implied by smallness of these terms.

Definition of
Let (u 0 , π 0 , R 0 ) be a solution to the Non-Newtonian-Reynolds system (10), and δ, η ∈ (0, 1] as in Proposition 1. We define where • u p is a perturbation based on the Mikado flows of Sect. 4, aimed at decreasing R 0 , • u c is a corrector restoring solenoidality of u 1 and compensating for our Mikado flows not solving Euler equations. Since (u 0 , π 0 , R 0 ) solves (10), it holds in the sense of distributions 5.1. Decomposition of R 0 and energy control. In general, R 0 is only continuous (recall Definition 2 and Remark 2). Since it is convenient to work with smooth objects, we regularize R 0 (extended for times outside [0, 1] by R 0 (x, 0) and R 0 (x, 1), respectively) with the standard mollifier φ in space and time. Thus Now R 0 is smooth and (12) implies Next, we decomposeR 0 into basic directions. In order to stay within S + of Lemma 2, we shift and normaliseR 0 via The role of √ 2 + .. is to avoid the degeneracy |R 0 (x, t)| = 0, whereas the role of γ 0 is to pump energy into the system, thus facilitating the step (11) → (14). Observe that γ 0 > 0 because of (11). The choice (44) yields in particular ≥ 2|R 0 (x, t)| and hence Remark 8. The set N ⊂ S + of Lemma 2 is fixed by (46) uniformly over the convex integration iterations. Define Thanks to Lemma 2 it holds The disjoint supports of W k (t), W k (t), k = k imply Recall the notation (48) and (50) to write We therefore have Remark 9. Observe that for the original Mikados, or their concentrated version, the second line of (52) vanishes. These additional terms will be taken care of by (37) and (38).
In order to avoid troublesome solenoidality correctors of the last in (52), let us (Helmholtz) project it onto divergence free vectors by P H = I d − ∇ −1 div and balance the identity by incorporating ∇ −1 div into the pressure, with the new pressurẽ Applying P =0 to both sides of the resulting identity, we arrive at

Choice of u c .
The corrector term u c has the following roles: (i) to cancel the highestorder bad term 'div (W k ⊗ W k )' of (53) via (38), (ii) to render the entire perturbation u p + u c solenoidal and (iii) null-mean. For (i), observe that (38) implies Thus taking will allow to cancel, with a part of the time derivative of u I c , the bad term 'div (W k ⊗W k )' of (53).
For (ii), observe that thanks to P H , u I c is already solenoidal. Therefore it suffices to compensate lack of solenoidality of u p . We will now define u I I c accordingly. By the definition (49) of u p , the definition (34) of W k , and since div k λ 2 k = 0 (cf. the property (i) of (26) Therefore we define where R 2 N is the double antidivergence given by Proposition 4, and N will be fixed later (see the discussion at the beginning of Sect. 6).
Since div div R 2 N = Id P =0 , div u I I c + div u p = 0. As u I c , u I I c are null-mean, to take into account the condition (iii), it suffices to define

Reynolds stresses.
Let us distribute ∂ t (u I c + u I I c ) andR 0 in (41) as follows We rewrite the r.h.s. of (58) further, recasting it into a divergence form.
(i). First line of r.h.s. of (58)) The definition (49) of u p , the definition (34) of W k , and 0 = k · ∇ k λ 2 (x) via property (i) of (26), give together Using the above formula and the definition (56) of u I I c , we define the antidivergence of the first line of r.h.s. of (58) (ii). Second line of r.h.s. of (58)) (iii). Third line of r.h.s. of (58)) Here we use an important idea of [10]. Via the definition (55) of u I c and the property (54) we have Adding the above identity to (53) that expresses div (u p ⊗u p −R 0 ), the 'div (W k ⊗ W k )' term cancel out and one has Let us thus define, leaving out ∇˜ since it will be accounted for by the pressure perturbation q p , so that div R quadr = ∂ t u I c + div (u p ⊗ u p −R 0 ) − ∇˜ . (iv. Fourth line on r.h.s. of (58)) Let (v. Last line of r.h.s. of (58)) Let

Pressure.
In order to balance for ∇˜ and ensure null-trace ofR 1 , we choose π p such that having freedom in choosing c(t), we set it so that π p (x, t) is null-mean. with Notice that tensor R 1 is symmetric, since its components are symmetric (cf. respective definitions and Proposition 4).

Estimates
We continue the proof of the main Proposition 1. In the previous section, given (u 0 , π 0 , R 0 ) solving the non-Newtonian-Reynolds system, we defined u 1 = u 0 +u p +u c , π 1 , and the new Reynolds stress R 1 , required by Proposition 1. The perturbation u p , the corrector u c and the error R 1 depend on the six parameters which satisfy the condition (35). The mollification parameter helps to avoid degeneracies or singularities of A. Let us immediately fix it so that where η > 0 is the parameter appearing in the statement of the main Proposition 1.
In this section we estimate u p , u c , R 1 and the energy gap of the new solution u 1 in terms of the remaining five parameters μ, λ 1 , λ 2 , ω, N . They will be appropriately chosen in Sect. 7 so that (13a)- (14) hold, thus concluding the proof of the main Proposition 1.
Remark 10. (Silent assumptions) For results of this section, we assume without writing it explicitly at each occasion: (u 1 , π 1 , R 1 ) is the triple constructed in the previous section, the set of directions K is fixed by Remark 8, the relations (35) hold, the assumptions of Proposition 1 hold, the choice (69) holds.
More precisely, we denote by M any constant depending only on the following parameters On the other side, we will denote by C (possibly with subscripts) any constant depending not only on the universal quantities (70), but also on (u 0 , π 0 , R 0 ), δ, η given in the assumptions of Proposition 1, i, r if we are estimating the W i,r norm N the (not yet fixed) parameter in (67), Constants C will be controlled within each iteration step (i.e. in the proof of Proposition 1) by appropriate choices μ, λ 1 , λ 2 , ω.

Estimates for velocity increments.
We will use now the improved Hölder inequality (18) and the preliminary estimates to control u p , u c .

Proposition 7 (Estimates for the principal increment u p ). For every r ∈ [1, ∞]
Moreover, there is a universal constant M > 0 such that Proof. The definition (49) of u p yields Using (74) to control |a k | ∞ and (36) to control |W k (t)| r , we obtain (75). An analogous computation gives (76). Notice that for r = 2 the power of μ in (75) is 0: for this reason, to reach (77), one needs more care. Recall from Proposition 5 that W k is λ 1 -periodic. Therefore we may apply improved Hölder inequality (Proposition 3) to the r.h.s. of (78) and use (36) for Let us now use (73), (74) to control a k , reaching (77). Now we deal with the corrector u c . In order to shorten the related formulas, let us introduce Observe that r → L(r ) is a non-decreasing map.

Proposition 8 (Estimates for the corrector u c ). For every r ∈ (1, ∞) it holds
Proof. Recall that u c = (u I c + u I I c ) − − T d u p by its definition (57), with u I c defined by (55) and u I I c defined by (56). The Calderón-Zygmund estimate, (74), and (36) give For the estimate of u I I c inẆ i,r we use the inequality (23) for j = i + 1, s = ∞, which yields Now we estimate a k by (74) and k μ,λ 1 by (32). Using the assumption (35), we arrive at Recall u c = (u I c + u I I c ) − − T d u p by its definition (57). Therefore, putting together (81), (82), and (17) with r = 1 and v = k , yields (80).

Estimates on the Reynolds stress. Recall (67)
In this section we estimate each term of R 1 . For our further purposes L 1 estimates suffice, but due to using Calderón-Zygmund theory, some estimates are phrased as L r ones.

Proposition 9 (Estimates on the principal Reynolds R quadr ). For every r
Proof. Recall the definition (62) of R quadr . Let us estimate its three terms in order of their appearance. (i) The first term of R quadr is the sum over k ofR 1 (∇a 2 k , P =0 (W k ⊗ W k )). The term P =0 (W k ⊗ W k )) is null mean and it oscillates at the frequency λ 1 , since W k does. Therefore (21) with (74) and (36) give The second term of R quadr is the sum over k of a 2 k (− W k ⊗ W k − k ⊗ k). We use estimate (74) to control a k terms and (37) to control the Mikado terms (iii) The last term is the sum over k of div −1 P H P =0 (∂ t a 2 k )Y k . We deal with div −1 via (19) and with P H via Calderón-Zygmund, control ∂ t a k using (74), and use (36) to estimate Y k . Hence (86) Together, (84), (85), (86) yield (83).

Proposition 10 (Estimate on R lin ). It holds
Proof. Recall the definition (60) of R lin . It involves three terms, which we estimate in order of their appearance. (i) The first term of R lin is the sum over k of Using (21) with |u| s = | k | ∞ , the assumption (35), and disposing of a k as usual, one has (88) (ii) The second term of R lin is the sum over k of We observe that Using the computation (89) in (23) with j = 0 and |u| s = | k | ∞ , we get (iii) The third term of R lin equals u 0 ⊗ u p + u p ⊗ u 0 , so we write using (75) Putting together (88), (90), (91), and observing that the right-hand sides of both (88) and (90) are estimated by C d,N ωμ − d (35), one has (87).

Proposition 12 (Estimates on the dissipative Reynolds R A ). For q ∈ (1, ∞) being the growth parameter of A it holds
for ν 0 > 0, q≤2, any r > 1, Proof. By definition (65), we have Therefore the inequality (15) gives the pointwise estimate Using Jensen inequality and 1 q−1 ≥1 in the first case, and Hölder inequality with q − 1, q−1 q−2 in the last case, one has For any s ∈ (1, ∞) the estimate (80) controls |Du c | s via λ 2 L(s), whereas λ 2 μ d/2−d/s controls |Du p | s thanks to (76). This closes the case q < 2 of (93). Recalling that u 1 = u 0 + u p + u c and that C may contain norms of u 0 , we obtain the case q ≥ 2.
Remark 11. For the current purpose of proving Proposition 1 and thus Theorem A, the case q ≤ 2 of (93) suffices. We included already the case q ≥ 2, because it is needed to prove Theorem B.
Immediately from the definition of R moll in (64) and the choice of in (69) we have

Proposition 13 (Estimate on R moll ). It holds
where η is the parameter appearing in the assumptions of Proposition 1.

Estimates on the energy increment.
We intend to approach the desired energy profile e(t), i.e. perform the step (11) → (14). Let us thence define δ E as follows Recall quantities L of (79). We will show Proposition 14 (energy iterate). For q ∈ (1, ∞) being the growth parameter of A it holds Proof. Recall (51). Taking its trace and recalling thatR 0 is traceless we have By the definition (44) it holds d = 2d 2 + |R 0 | 2 + dγ 0 , therefore Integrating and using 2 + |x| 2 ≤ + |x|, we have (98) We estimate the first two terms of the r.h.s. of (98) using (69) and (43) as follows where in the second inequality we used the assumption e(t) ≥ 1 2 . This in (98) yields The first integral of r.h.s. of (99) involves a λ 1 -oscillating function P =0 |W k | 2 , recall Proposition 5. Therefore, using (17), then (74) to control a k and (36) for W k , we have For the integral following the second sum in (99), we use (37) and (74) to get We plug (100) and (101) to (99) and obtain Use u 1 = u 0 + u p + u c in the definition (95) of E to write for the time instant t where a cancellation occurs, thanks to the definition (45) of γ 0 . Inequality (102) allows to control the first term of the r.h.s. above. For the last term we use (16), next Hölder inequality with q, q q−1 , and finally u 1 = u 0 + u p + u c to get Thus, integrating in time over Consequently The terms in the second line above are estimated, using (80) for u c and (75) for u p , by C N (L(2) 2 + L(2) + μ −d/2 ). Observe that L(2) of this term can absorb λ −1 2 λ 1 μ of the first line. The terms in the last line are estimated by λ 2 μ (λ 2 L(q)) q , using (76) for Du p , (80) for Du c . We thus arrived at (96).

Proof of the Main Proposition 1
Having at hand the estimates of the previous section, we are ready to show that (u 1 , q 1 , R 1 ) constructed in Sect. 5 satisfy the inequalities (13a)- (14).
The estimates of the previous section have at their right-hand sides two type of terms: ones where the parameters λ 2 , λ 1 , μ, ω are intertwined, and the remaining ones. These remaining ones can be made small simply by choosing the relevant parameters large. The terms with λ 2 , λ 1 , μ, ω interrelated need more care, so let us focus on them. They contain two little technical nuisances: (i) appearance of N and (ii) estimates for some parts of R not holding in L 1 . Let us ignore these nuisances for a moment, which is easily acceptable after recalling (i) R ∞ of Remark 6 (which heuristically cancels the terms involving N ) and that (ii) estimates for R hold in L r for any r > 1, whereas an of room is assured by the assumed sharp inequality q < 2d d+2 . So for a moment let us consider estimates of Sect. 6 allowing for |R| 1 and disregarding the terms with N . After inspection, we see that smallness of their right-hand sides where λ 2 , λ 1 , μ, ω are intertwined, needed for Proposition 1 is precisely the smallness of (40), Remark 7. Therefore we will proceed as follows.
Firstly, guided by (40), we will choose relation between magnitudes of λ 2 , λ 1 , μ, ω. To this end we postulate and choosing relation between magnitudes means picking a, b, c so that (40) are strictly decreasing in λ.
Secondly, we will need to make sure that when r > 1 and N appear, the relations between magnitudes do not change. This will be achieved by choosing N large and r small in relation to a, b, c.
Finally, we will send λ → ∞ to reach (13a)- (14). Picking magnitudes a, b, c. The requirement that powers in (40) rewritten in terms of (103) are negative reads These conditions on a, b, c can be simultaneously achieved as follows.
(1) The conditions not involving b amount to the requirement From the assumption q < 2d d+2 of Proposition 1 it follows that d/q − d/2 > 1. Therefore satisfying (105) is possible with a large. More precisely, let us pick Then between 1 + a and (d/q − d/2)a there are at least two natural numbers. We then fix c ∈ N as the largest natural number satisfying (105). Notice that there is still at least one natural number between 1 + a and c.
This is possible, because, as observed in point (1), there is at least one natural number between 1+a and c and thus also between ad/2 and (d/2−1)a +c −1. The condition (107) automatically verifies the two conditions concerning b.
Let us denote by −ζ < 0 the largest power of those appearing in (104). We have just showed

Fixing N and r
This choice yields Using the definition (79) of L with (108) and (110), one has Importantly, fixing the gauge N freezes all C N 's in estimates to C. Let us fix also an exponent r 0 ∈ (1, ∞) (close to 1) such that This is possible because the l.h.s. above vanishes as r 0 → 1.
Recall that R 1 = −(R lin + R corr + R quadr + R moll + R A ) by its definition (67). By (87) we have, with C N now fixed to C by the choice (109) of N where for the second inequality we invoked μ = λ a by (103), (108), and (110). Similarly, using (92) and (111) For the L 1 -estimate of R quadr we need to switch to the L r 0 estimate, where r 0 was fixed in (112). We have, using (83) Thanks to the choice of r 0 in (112), we hence have Similarly, for the estimate of R A we use (93) with q ∈ (1, 2), obtaining Therefore by (108), (111) and q − 1 ∈ (0, 1) Together, the terms R lin , R corr , R quadr , R A are bounded in view of, respectively, (115), (116), (117), and (118) by Cλ −ζ with certain ζ > 0: Therefore, using for the remaining R moll the estimate (94), we have thus showing (13c) by taking λ large. Let us show the last remaining inequality (14). By (96), with C N fixed to C by the choice (109) of N , we have in view of (108) and (111) (120) The proof of Proposition 1 is concluded.

Proof of Theorem A
We will iterate Proposition 1. Let us start at the trivial solution (u 0 , π 0 , R 0 ) ≡ 0 with δ 0 = 1. At the nth step we take δ n := 2 −n and η n := δ n+1 which is the assumption (12) of the step n + 1. Similarly, for any t ∈ [0, 1] at the n-th step we get, by (14) which is the assumption (11) of the step n + 1, since Consequently we obtain iteratively, as η n ≤ 2 −n , Inequalities (121) mean that {u n } ∞ n=0 is a Cauchy sequence in C(L 2 ) ∩ C(W 1,q ). Denote its limit by v ∈ C(L 2 ) ∩ C(W 1,q ). Send n → ∞ in the distributional formulation of (10). In particular, in order to pass to the limit in the dissipative term, take a test function ϕ and use (15) for q < 2d d+2 < 2 and the Hölder inequality to obtain The right-hand sides tend to 0 as n → ∞ thanks to (121). Consequently we see that v satisfies the distributional formulation of (1). For the 2 which via Hölder inequality and (121) allows to pass with n → ∞. This and lim n→∞ |u n − v| C(L 2 ) = 0 provided by (121) yields (7). Let us now focus on proving the last part of Theorem A, i.e. the non-uniqeness statement. Let us take the two energy profiles e 1 , e 2 and the respective triples (u 1 n , π 1 n , R 1 n ) and (u 2 n , π 2 n , R 2 n ) of our convex integration scheme (in what follows, superscripts denote the cases of e 1 , e 2 , respectively). At each iteration step n → n + 1 one picks value of λ i n (= λ of Sect. 7.3) that works for (u i n , π i n , R i n ). Observe that choosingλ n = max (λ 1 n , λ 2 n ) works simultaneously for both triples. Thus, without renaming the triples, let us make the choiceλ n for both (u i n , π i n , R i n ), i = 1, 2. It results in using identical Mikado flows W k for both iterations. Now we want to inductively argue that, thanks to the assumed e 1 (t) = e 2 (t) for t ∈ [0, T ], it holds u 1 n (t) = u 2 n (t) for every n and t ∈ [0, T − 1 2 7 d ]. Let us assume thence that u 1 n (t) = u 2 n (t) andR 1 n (t) =R 2 n (t) for times t ∈ [0, T − n i=0 2 −i 2 7 d ] (This holds for n = 0, since we begin with the zero triple). Formula (47), with (44) and (45) shows that a i k,n+1 (t) (i.e. every a i k (t), k ∈ K at the step n → n + 1) depends on e(t), R i, n (t) and u i n| [0,t] , with ≤ 2 −(n+1) 2 7 d being the mollification parameter, cf. (69) with the choice δ n+1 := 2 −(n+1) , and the u i n -dependence being nonlocal due to the dissipative term in (45). So by our inductive assumption we see that a 1 k, Consequently, via the definition (49), the principal perturbations u i p (t), i = 1, 2 at the step n → n + 1 are identical for times t ∈ [0, T − n+1 , since the correctors and the new errors are defined pointwisely in time.
Under the assumption that e 1 , e 2 are identical on [0, T ], we produced iteratively u 1 n (t), u 2 n (t) that agree for t ∈ [0, T − 1 2 7 d ] thus also their limits satisfy v Replacing T − 1 2 7 d with any fixed number T 1 strictly smaller than T requires only mollifying at the scales below T − T 1 instead of 1 2 7 d .

Sketch of the Proof of Theorem B
Let us indicate changes needed in proofs of Proposition 1 and Theorem A to reach Theorem . Now, we extend the allowed range of growths of A to q ∈ (1, 3d+2 d+2 ) at the cost of abandoning the control over the dissipative term 2 t 0 A(Dv)Dv of the energy. Recall that r ∈ (max{1, q − 1}, 2d d+2 ) is an additional exponent, fixed in the assumptions of Theorem B. The main observation is that is subcritical in the sense of choices made in Sect. 7.1.
Let us first consider modifications in proof Proposition 1. Replacing in Sects. 7.1, 7.2 q of Proposition 1 with r implies the following analogue of (111): for some positive ζ > 0. Consequently, (114) holds now also with r in place of q. Next, since q now may exceed 2, to control |R A (t)| 1 we use the entire (93) to write In any of these cases, right-hand sides are controlled by powers of λ 2 μ d 2 − d r and λ 2 L(r ), therefore we can reach (119). Finally, since we abandon the control over the dissipative term in the energy inequality, (45) and δ E of (95) (let us rename it to δẼ) loose their dissipative terms. The consequence of the latter is that (96) simplifies to Inequality (122) allows to prove (120) for δẼ as in Sect. 7.
The above modifications allow to prove Theorem B along Sect. 8 with the difference that now {u n } ∞ n=0 forms a Cauchy sequence in C(L 2 ) ∩ C(W 1,r ) and, in the case q > 2, we pass to the limit in the dissipative term via (15) for q ≥ 2 and the Hölder inequality that give

Sketch of the Proof of Theorem C
Let us first introduce the following modification of Definition 2 with spatial null-mean u, solving the Cauchy problem in the sense of distributions, where the data are attained in the weak L 2 -sense.
The drop of regularity between objects of Definition 3 and objects of Definition 2 stems from a different starting point for our iterations. To prove Theorems A, B, we started the iteration at the smooth triple (u 0 , π 0 , R 0 ) = (0, 0, 0) and added smooth perturbations in each iteration. To prove Theorem C we will start iterations with (v a ,π a , −v a ⊗v a ), where v a ,π a solves the Cauchy problem of a non-Newtonian-Stokes system: The smoothness of such v a ,π s is in general false, so even though at each step we again add smooth perturbations, the regularity at each step cannot be better than that of v a ,π a solving (124).
The starting point of our iterations is given by

Proposition 15 (Leray-Hopf solutions for non-Newtonian Stokes
solving the Cauchy problem (124) in the sense of distributions, so that lim t→0 |v a (t) − a| 2 = 0 and The proof uses monotonicity of A and for strong attainment of the initial datum, the energy inequality.
The main ingredient of proof of Theorem C is a version of Proposition 1 tailored to deal with the Cauchy problem.
Roughly speaking, given a solution to the Non-Newtonian-Reynolds system (123), which assume the given initial datum, we construct another solution to (123) with the same initial datum, and with a well-controlled Reynolds stress. The price we pay to keep the initial datum intact is growth of energy. More precisely, energy of the ultimately produced solution to the Cauchy problem with datum a for (1) is much above an energy of a non-Newtonian Stokes emanating from the same a. Hence we cannot reach energy inequality, even for merely the kinetic energy in the range q < 2d d+2 . This is why we do not distinguish the subcompact range q < 2d d+2 in Theorem C. We are ready to state Define u p , u c by (49), (57) respectively (with the new ). Let us introduce a smooth cutoff and define the perturbationsũ p ,ũ c as follows Due to (128), the new version of (77) reads Since (125) Our assumption now does not controlR 0 (t) for t ≤ 2σ , but we can always write, choosing σ which, together with the smallness of |ũ c (t)| L 2 , cf. (80), gives the case [σ/2, 4σ ] of (126a). On [0, σ/2], it holds u 1 = u 0 , as there χ ≡ 0, thence the respective part of (126a). The estimate (126b) holds on the whole time interval, because the Mikado flows are small in W 1,r for r < 2d d+2 by construction. The estimate on the new Reynolds stress (126c) on [σ, 1] is analogue to the corresponding estimate (13c) of Proposition 1, as on [σ, 1] the time cutoff χ ≡ 1. The estimate on [0, σ/2] is trivially satisfied, as on this time interval the cutoff χ ≡ 0 (here as in the proof of Proposition 1, thus giving the η in the second line of (126c), whereas the term R 0 (1 − χ 2 ) is responsible for |R 0 (t)| L 1 in the second line of (126c). There is, however, in the new Reynolds stress R 1 an additional term coming from the time derivative of the cutoff: div For the first term in (130) we just use (75) with r = 1: where we used the choice μ = λ a , as in Sect. 7. For the second term in (130) we have, invoking (19) Since by (111) L(2) ≤ Cλ −ζ , both terms of (130) are estimated by negative powers of λ. Thus they can be made as small as we wish by picking λ big enough.
Let us now justify (127) along the proof of Proposition 14. In (98) γ 0 changes to γ . Now we do not control |R 0 (t)| 1 on the entire time interval [0, 1], only on [4σ, 1] via (125). On this interval χ = 1 and hence for any t ∈ [4σ, 1] one has the following counterpart of (102) We use now u 1 = u 0 +ũ p +ũ c to write for any t ∈ [4σ, 1] The r.h.s. above can be made arbitrarily small in view of (131) and arguments analogous to that leading to Proposition 1, which yields (127).
Iterating Proposition 16, we can now complete the proof of Theorem C. Let us choose σ n := 2 −n along the iteration. We will choose δ n = 2 −n and η n = 2 −(n+9) d −1 as in proofs of Theorems A, B. There are two main differences between the current iterations and the iterations leading to Theorems A, B. Firstly, we initiate the iterations with the triple (v a ,π a , −v a ⊗ v a ), where v a ,π a are given by Proposition 15. Secondly, we will choose the now additional free parameter γ just to distinguish between different solutions. Namely, let us choose γ n = d −1 δ n except for γ 3 , which we require it to be a large constant, say K .
The condition (125) for the initial triple is void (empty interval where it shall hold) and over iterations it is satisfied thanks to the first case of (126c) and our choices for η n , δ n , σ n . The third iteration produces u 3 out of u 2 such that At this step the energies of the iteratively produced solutions branch: choosing two K 's that considerably differ, we will see that the kinetic energies on t ∈ [1/2, 1] of the finally produced solutions differ considerably. From step n = 4 onwards γ n = δ n , thus the first line of (126a) is analogous to (13a). Iterating Proposition 16 we thus obtain convergence of the sequence Note the open side of an interval above. Taking into account the regularity class of v a and r < q, we thus have which allows to pass to the limit in the distributional formulation of (124), since by choice r > max (1, q − 1).
Concerning the attainment of the initial datum a, for any q 0 < 2 the estimate (126b) yields |v ∞ | q 0 (t) → 0 as t → 0. Therefore v = v ∞ + v a satisfies |v − a| q 0 (t) → 0 as t → 0. From this and the fact that the L 2 norm of v(t) is uniformly bounded in time on [0, 1], it follows that v(t) a weakly in L 2 as t → 0 + . Let us finally argue for multiplicity of solutions. At the step n = 3 let us choose two different K , K . Let us distinguish the resulting u n 's and their limits v by, respectively, u n and u n , and v and v . On t ∈ [1/2, 1] (127) yields for n ≥ 4 whereas The same inequalities hold for the strong limits v, v . Therefore, for d|K − K | > 1 they must differ.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Proof of Lemma 1. Let us first consider a scalar function
It is Lipschitz for ν 1 = 0. In the range q ∈ (1, 2] f is (q − 1)-Hölder continuous for ν 0 = 0, ν 1 > 0; and locally Lipschitz for ν 0 > 0, ν 1 > 0. By the last statement we mean that (133) It is proven by writing with the second inequality given by splitting q − 2 = (q − 1 − δ) + (δ − 1) so that q − 1 − δ > 0 and thus the function t → t q−1−δ is increasing. For the integral in the second line of (134) we use that for γ ∈ (−1, 0] and |t| + |s| > 0 it holds Computation for (135) is contained in proof of Lemma 2.1 in [1]. Using (135) in (134) we obtain (133). In the range q ≥ 2 for f , the function t → t q−2 is increasing and integrable at 0, so (133) holds for any ν 0 ≥ 0, ν 1 ≥ 0. Altogether Replacing scalars with tensors in (136) will give (15). Details follow. Let us consider first the (global Hölder) case ν 0 = 0, q < 2, i.e. the first line of (136). We show the respective first line of (15) by considering two cases: (i) Q, P lie along a line passing through the origin and (ii) Q, P lie on a sphere centered at the origin. The case (i) is Q = t Q 0 , P = s Q 0 for some s, t ∈ R and |Q 0 | = 1, thus (15) here follows immediately from (136). The case (ii) is |Q| = |P| = l for some l > 0. Here and the latter is a l-independent constant. Both cases (i), (ii) yield for every nonzero Q, P arguments lie: on the sphere ∂ B |Q| (0) on the same line through origin We have just proven the first ν 0 = 0, q < 2 case of (15). The second (global Lipschitz) case ν 0 > 0, q < 2, i.e. the tensorial version of the second line of (136), can be proven analogously. But in fact the computation leading to (133) works well when applied immediately to tensor mappings both in the case ν 0 > 0, q < 2 and q ≥ 2, giving the remainder of (15). Estimate (16) follows from an argument that gave the case q ≥ 2 in (15), applied tof (t) = (ν 0 + ν 1 |t|) q−2 t 2 .
(ii). R N : C ∞ (T d ; R) × C ∞ 0 (T d ; R d ) → C ∞ 0 (T d ; S)) (Compare [18], Proposition 5.2 and Corollary 5.3.) We construct the two-argument improved symmetric antidivergence iteratively upon div −1 . Let us commence with Our aim is an antidivergence that extracts oscillations of u λ out of the product f u λ . Therefore (137) suggests to apply div −1 to u λ , and correct the remainder. Let us thence compute and define R 1 ( f, u) := f div −1 u − d k=1 R 0 ∂ k f, (div −1 u)e k . The choice (138) of R 0 and (139) yield div R 1 ( f, u) = f u − − f u.
Define inductively R N ( f, u) := f div −1 u − d k=1 R N −1 ∂ k f, (div −1 u)e k , per analogiam with the construction R 0 → R 1 . Using induction over N , one proves now that R N is bilinear, symmetric, div R N ( f, u) = f u − − T d f u, satisfies Leibniz rule: and estimates (21), similarly to proof of Lemma 3.5 of [30]. For instance, to show div R N ( f, u) = f u − − T d f u we have (140) as the initial step. Assuming div with the second equality due to (139). The estimate (21) holds for N = 1 since (19) and Jensen imply |R 0 ( f, u)| p ≤ C| f | r |u| s and thus with the last inequality due to (20). For the inductive step N − 1 → N , let us compute, using (137) with the second inequality valid via the inductive assumption, i.e. (21) for N − 1. The estimate (19) and interpolation yield the inductive thesis.
and since it is a linear combination of R N 's, it retains all its properties.

Proof of
Since d ≥ 3, the countable union of planes {ζ k + sk + s k } s,s ∈R + Z d has zero measure and, since k, k ∈ Z d , it is closed in R d . Therefore, for any ρ ∈ (0, ρ ), also k ∈K B ρ (0) + ζ k + sk + s k } s,s ∈R + Z d is closed in R d and, if ρ small enough, it is strictly contained in R d . We can thus find ζ k and ρ ∈ (0, ρ ) such that with the superset being open, thus concluding the proof of part (i). For part (ii), consider the set of directions in Z 3 given bỹ By part (i), we can find {ζ k } k∈K ⊆ R 3 and ρ > 0 such that (in R 3 ), for all k, k ∈ K and s ∈ R (the balls above are in R 3 ). Observe that, for every (k, 1) ∈K , the axis of the cylinder B ρ (ζ k ) + s(k, 1) + Z 3 ⊆ R 3 is never lying on the horizontal plane x 1 , x 2 . Therefore we can assume w.l.o.g. that eachζ k has the form ζ k = (ζ k , 0).
Let us now show (28), with such choice of {ζ k } k∈K . Assume by contradiction that there are two distinct directions k, k ∈ K and s ∈ R such that Then there are x ∈ B ρ (ζ k ) ⊆ R 2 , y ∈ B ρ (ζ k ) ⊆ R 2 , l, l ∈ Z 2 , such that x + sk + l = y + sk + l ∈ R 2 .
This implies that