TEMPORAL REGULARITY OF SYMMETRIC STOCHASTIC p -STOKES SYSTEMS

. We study the symmetric stochastic p -Stokes system, p ∈ (1 , ∞ ) , in a bounded domain. The results are two-folded. First, we show that in the context of analytically weak solutions, the stochastic pressure – related to non-divergence free stochastic forces – enjoys almost − 1 / 2 temporal derivatives on a Besov scale. Second, we verify that the velocity component u of strong solutions obeys 1 / 2 temporal derivatives in an exponential Nikolskii space. Moreover, we prove that the non-linear symmetric gradient V ( εu ) = ( κ + | εu | ) ( p − 2) / 2 εu , κ ≥ 0 , has 1 / 2 temporal derivatives in a Nikolskii space.

Key words and phrases.SPDEs, Non-linear Laplace-type systems, Strong solutions, Regularity, Stochastic p-Stokes system, power-law fluids, generalized fluids.
The system (1.1) is called stochastic symmetric p-Stokes system and describes the evolution of power-law fluids in the regime of laminar flow.Ultimately, one wants to substitute (1.1a) by du + (u • ∇)u − div S(εu) + ∇π dt = G(u) dW (t) (1.2) to include convective effects.This might lead to turbulent flow.Maybe most prominent is the case p = 2. Then the system reduces to the infamous stochastically forced Navier-Stokes equations.A derivation of this model and how meaningful choices of G might look like can be found in e.g.[MR04] (see also [BFH17] and the references therein).
While the system (1.1) does not bother with the mathematical challenges of turbulence, it is still difficult to obtain improved regularity.There are several reasons: (a) Solutions to non-linear equations lack smoothness even for smooth data.(b) Limited regularity of the stochastic data G(•) dW leads to limited regularity of the solution.(c) Stochastic integration in general Banach spaces causes technical restrictions.(d) Delicate interplay between the divergence free constrain (1.1b) and the pressure.In this article we focus on the temporal regularity for solutions to the system (1.1).Before we comment on the available literature on regularity and finally state the main contributions of this paper, we briefly discuss some results related to the existence of non-Newtonian fluids.
1.1.Existence.A mathematical theory for the existence of weak solutions to the deterministic counterpart of (1.1) enriched by the convective term was initiated by Ladyzhenskaya [Lad67;Lad69] and Lions [Lio69].Motivated by the works, many authors contributed to the construction of solutions, e.g.[MNR93; MNRR96; MNR01; Wol07; DRW10].Ultimately leading to the construction of a weak solution for p > 2n/(n + 2).Below the compactness threshold W 1,p x → → L 2 x wild solutions may emerge, see e.g.[BV19].
The stochastic case is less known.An answer for solution independent noise coefficient was given in [TY11;Yos12].The authors require the condition p > 9/5 in the three dimensional setting.It was extended in [Bre15] to cover more general noise coefficients and a unified condition p > (2n + 2)/(n + 2).
A construction for stochastic electro-rheological fluids, here p also depends on Ω × (0, T ) × O, is done in [BG19].
1.2.Regularity.Regularity is needed in order to control discretizations.In particular, regularity moderates the speed of convergence of numerical schemes, see e.g.[BDSW21].Therefore, and out of mathematical curiosity, many authors have investigated the regularity of solutions to (1.1) and related models.
1.2.1.Non-linear.In general, there is a barrier that does not allow for arbitrary high regularity even though the data is smooth, cf.[IM89].Still partial results can be obtained.
C 1,α -regularity has been obtained in e.g.[CM99; AM02; DER07; Fan07].Hölder regularity is sub-optimal from a numerical point of view.Instead, regularity of the expression V (εu) = (κ + |εu|) (p−2)/2 εu is more suitable for the investigation of numerical algorithms as observed by Barrett and Liu in [BL93].Indeed, the operator V measures the monotonicity of S in a natural way, i.e., for all A, B ∈ R n×n it holds A first result for non-singular shear thinning fluids, κ = 1 and p ∈ (1, 2), has been obtained in [Nau88].The author verifies local regularity V (εu) ∈ W 1,2 loc .A similar result has been obtained for shear thickening fluids, p ≥ 2, in [Bre12].Additionally, local regularity for the non-linear tensor S and the pressure π was found in [JKL14].They show ∇S(εu), ∇π ∈ L p loc .In [DKS14] regularity transfer of the data to the solution in terms of Campanato and weighted BMO (bounded mean oscillation) spaces is shown.
Spatial regularity results for the parabolic p-Stokes system has been obtained in [AMS04].The authors show that ∇u ∈ L ∞ t L 2 loc and ∇V (εu) ∈ L 2 loc (O T ).It was extended to the ϕ-Stokes problem in [BK17].
As far as we know, improved temporal regularity for the parabolic symmetric p-Stokes system has not been analyzed in the literature.Only a few results that neglect the appearance of the pressure but keep the symmetric gradient are available.For example, the parabolic symmetric p-Laplace system has been studied by Frehse and Schwarzacher in [FS15] and by Burczak and Kaplický in [BK16].In [FS15] it is shown that x .Here B denotes a Besov space (see Section 2.1 for more details).1.2.2.Stochastic -spatial regularity.In general, the tools from deterministic theory are not available for stochastic equations.In particular, one cannot test the equation by a test function.Instead, one needs to perform an expansion for a suitable functional.This method is called Itô's formula.While some test functions can be recovered in this way others cannot.Still, the deterministic regularity theory acts as a guiding example on what regularity can be expected.
It has been observed by Gess [Ges12] that stochastic partial differential equations, arising as a stochastically perturbed gradient flow of an energy, can lead to improved regularity results.He proposes conditions that ensure the existence of a strong solution.
Fortunately, the Helmholtz-Leray projection Π div (see (2.3)) allows to reformulate (1.1a) and (1.1b) into ) In other words, the Helmholtz-Leray projection factorizes the equations for the velocity field and the pressure.It enables to seek for the velocity independently of the pressure.Now it is possible to interpret (2.20a) as a stochastically perturbed gradient flow of the energy induced by the potential ϕ κ (t) := t 0 (κ + s) p−2 s ds on the space W 1,p 0,div .Indeed, (2.20a) can be written as where DJ (u) = Π div div S(εu).Therefore, it is natural to look for regularity estimates that match the concept of strong solutions.
A similar approach is used by Breit and Gmeineder in [BG19].They prove spatial regularity as in the deterministic case also for stochastically forced electrorheological fluids and obtain ∇u x .An analogous line of argumentation can be used to establish spatial regularity for the p-Stokes system.1.2.3.Stochastic -temporal regularity.In sharp contrast to spatial regularity, where similar results as in the deterministic theory can be expected, is temporal regularity for stochastic partial differential equations.In general, temporal regularity is substantially worse compared to deterministic equations.This is due to the action of the random data G(•) dW on the solution u.
Even in the simplest case G ≡ 1, we expect u to not exceed the regularity threshold induced by the Wiener process W .A precise regularity result was derived by Hytönen and Veraar in [HV08].They show that P-a.s.W ∈ B 1/2 Φ2,∞ , where Φ 2 (t) = exp(t 2 ) − 1, and W ∈ B 1/2 p,q for any p ∈ [1, ∞] and q ∈ [1, ∞).Therefore, a natural space for temporal regularity is the Nikolskii-space of functions with exponential second moments and 1/2 derivatives.
Only recently, precise regularity results for stochastic integrals in 2-smooth Banach spaces have been found by Ondreját and Veraar in [OV20].Maybe most important, stochastic integration behaves as deterministic integration in the sense that improved regularity of the integrand leads to improved regularity of the stochastic integral itself.For example, if the integrand is bounded in time, then the stochastic integral is as regular as a Wiener process, i.e., Exactly this fact is used in [Wic21] to derive temporal regularity for the stochastic p-Laplace system.In particular, we manage to recover the limiting temporal In other words, the regularity of W fully transfers to regularity of u.Additionally, exploiting the V -coercivity, it is possible to derive non-linear gradient regularity The stochastic p-Stokes system is even worse.A major obstacle in the derivation of higher regularity for the system (1.1) is the pressure.If G(•) dW is not divergencefree, then the pressure is very rough in time (see e.g.[LRS03;CP12]).This can be seen by formally decomposing the pressure into a deterministic and a stochastic component π = π det + π sto , where While the deterministic pressure can be treated by classical arguments, we observe that the stochastic pressure looks like white noise in time.This leads to severe difficulties in numerical discretizations as soon as we approximate the velocity field by discretely divergence-free fields rather than exactly divergence-free fields.More details can be found in [FQ21; FPV21] and the references therein.
1.3.Main results.The main contributions of this paper are two-folded.1.3.1.Temporal regularity of stochastic pressure for weak solutions.We apply the general theory of monotone stochastic partial differential equations developed by Liu and Röckner [LR10], see also the book [LR15], to the p-Stokes system.More recently, the theory has been refined to cover more general equations, see e.g.[RSZ22;AV22].In this way, existence of weak solutions is almost immediate.The main novelty is a careful consideration of the regularity of the stochastic pressure.
1.3.2.Temporal regularity of velocity field for strong solutions.We exploit the gradient flow structure (1.5) on the space of solenoidal vector fields in order to use the existence theory of strong solutions initiated by Gess [Ges12].In doing so, we provide a sufficient condition for the existence of strong solutions to the p-Stokes system.
Theorem 2. Let the assumptions of Theorem 1 be satisfied.Additionally, let q ∈ [1, ∞), p ≥ 2, J (u 0 ) ∈ L q ω be F 0 -measurable and G satisfy Assumption 33.Then there exists a unique strong solution (u, π) to (1.1).Moreover, it holds We extend the techniques developed in [Wic21] to obtain improved temporal regularity for the velocity field and its non-linear symmetric gradient.
Exponential integrability requires slightly stronger moment estimates.In a similar way to the a priori bound (1.6a) for weak solutions, one can establish an estimate of the form In other words, if the initial condition has N 2 moments, so does the solution at later times and the assumption of part (b) is satisfied.
Theorem 5. Let p ≥ 2, Assumption 33 be satisfied and (u, π) be a strong solution to (1.1).Additionally, assume In Section 2 we introduce the mathematical framework.In particular, we discuss function spaces, the Helmholtz decomposition, the Bogovskii operator, potentials, mapping properties of stochastic integrals and gradient flows.
Section 3 deals with the construction of weak solutions.Moreover, regularity properties of the pressure are discussed.
Strong solutions are constructed in Section 4.
In Section 5 we present the temporal regularity of the velocity and the non-linear gradient.
More details on Besov norms are presented in the Appendix A.

Mathematical setup
Let O ⊂ R n , n ∈ N, be a bounded domain (further assumptions on O will be needed for the stability of the Helmholtz-Leray projection and the Bogovskii operator).For some given T > 0 we denote by I := (0, T ) the time interval and write O T := I×O for the time space cylinder.Moreover let (Ω, F, (F t ) t∈I , P) denote a stochastic basis, i.e., a probability space with a complete and right continuous filtration (F t ) t∈I .We write f g for two non-negative quantities f and g if f is bounded by g up to a multiplicative constant.Accordingly we define and .We denote by c a generic constant which can change its value from line to line.Inner products and duality pairings are denoted by (•, •) and •, • , respectively.
2.1.Function spaces.As usual, for 1 ≤ q < ∞ we denote by L q (O) the Lebesgue space and W 1,q (O) the Sobolev space.Moreover, W 1,q 0 (O) denotes the Sobolev spaces with zero boundary values.It is the closure of C ∞ 0 (O) (smooth functions with compact support) in the W 1,q (O)-norm.We denote by W −1,q (O) the dual of W 1,q 0 (O).The space of mean-value free Lebesgue functions is denoted by L q 0 (O).The space of smooth, compactly supported and divergence-free vector fields is called C ∞ 0,div (O) and its closure within the W 1,p -norm is abbreviated by W 1,p 0,div (O).We do not distinguish in the notation between scalar-, vector-and matrix-valued functions.
For a Banach space (X, • X ) let L q (I; X) be the Bochner space of Bochnermeasurable functions u : I → X satisfying t → u(t) X ∈ L q (I).Moreover, C 0 (I; X) is the space of continuous functions with respect to the norm-topology.We also use C 0,α (I; X) for the space of α-Hölder continuous functions.Given an Orliczfunction Φ : [0, ∞] → [0, ∞], i.e. a convex function satisfying lim t→0 Φ(t)/t = 0 and lim t→∞ Φ(t)/t = ∞ we define the Luxemburg-norm The Orlicz space L Φ (I; X) is the space of all Bochner-measurable functions with finite Luxemburg-norm.For more details on Orlicz-spaces we refer to [DHHR11].
Given h ∈ I and u : I → X we define the difference operator τ h : {u : The Besov-Orlicz space B α Φ,r (I; X) with differentiability α ∈ (0, 1), integrability Φ and fine index r ∈ [1, ∞] is defined as the space of Bochner-measurable functions with finite Besov-Orlicz norm (2.1) The case r = ∞ is commonly called Nikolskii-Orlicz space and abbreviated by N α,Φ = B α Φ,∞ .When Φ(t) = t p we write B α p,r (I; X) and call it Besov space.If X is reflexive, we denote by B −α p ,q (I; X ) = B α p,q (I; X) .Similarly, given a Banach space (Y, • Y ), we define L q (Ω; Y ) as the Bochner space of Bochner-measurable functions u : Ω → Y satisfying ω → u(ω) Y ∈ L q (Ω).The space L q F (Ω × I; X) denotes the subspace of X-valued progressively measurable processes.We abbreviate the notation L q ω L q t L q x := L q (Ω; L q (I; L q (O))) and L q− := r<q L r .2.2.Helmholtz-decomposition.The Helmholtz decomposition is a powerful tool in the analysis of fluids.It allows for a complete decoupling of the governing equations of the velocity field and the pressure.For a recent survey we refer to [BNPB12].More details can be found e.g. in [FMRT01, Chapter 2 Section 3].
The general idea is to decompose a vector field into a divergence-free vector field and a gradient of a potential.Let v ∈ L 2 x and G v be a solution to ∀ξ ∈ W 1,2 x : Note that G v is defined up to a constant.We will prescribe the constant by the condition G v O = 0. Now, we define the Helmholtz-Leray projection Π div and its orthogonal complement Π ⊥ div by In order to understand the range of the operator Π div and Π ⊥ div we introduce the spaces The next lemma can be found e.g. in [Tem79, Theorem 1.4].
Remark 7. Since Π div and Π ⊥ div are L 2 -orthogonal projections, an equivalent definition to (2.3) is given by Stability of the Helmholtz-decomposition is particularly important, when it comes to the reconstruction of the pressure.Clearly, as an L 2 -projection it trivially holds Higher order stability is closely linked to regularity of solutions to the Laplace equation.Indeed, if div v ∈ L 2 x , then (2.2) is equivalent to the Neumann-Laplacian Therefore, it is not surprising that improved stability properties of the Laplacian transfer to stability of the Helmholtz-Leray projection, cf.[Tem79, Remark 1.6].
Theorem 8. Let O be a C k+1 -domain and q ∈ (1, ∞).Then Proof.The claims immediately follow by Theorem 8 and the density of smooth functions.
2.3.Bogovskii operator.The right inverse of the divergence is nowadays called Bogovskii operator.It traces back to his work [Bog80] and is a useful tool for the regularity analysis of the pressure.

Monotonicity and potentials. A continuous, convex and strictly increasing function
is called an N -function.
Definition 11.Let ϕ be an N-function.We say that ϕ is uniformly convex, if ϕ is C 1 on [0, ∞) and C 2 on (0, ∞) and assume that uniformly in t > 0. The constants hidden in are called the characteristics of ϕ.
Associated to an uniformly convex N -function ϕ we define the tensors Sometimes it is convenient to change the growth of an N -function near 0. One possibility is to introduce the shifted N -function ϕ κ for κ ≥ 0 by We define S κ and V κ analogously to (2.12).However, we neglect in the notation the dependence of S κ and V κ on κ, since the characteristic of ϕ κ is uniform in κ ≥ 0.
Lemma 14 ([DE08, Lemma 32]).Let ϕ be an uniformly convex N -function.Then for all δ > 0 there exists C δ > 0 such that for all t, u ≥ 0 The following inequality follows immediately from Lemma 13 and Lemma 14.
Lemma 15 (Young type inequality).Let ϕ be an uniformly convex N-function.Then for each δ > 0 there exists C δ ≥ 1 (only depending on δ and the characteristics of ϕ) such that Lemma 16 (Change of Shift [DFTW20, Lemma 42]).Let ϕ be an uniformly convex N-function.Then for each δ > 0 there exists C δ ≥ 1 (only depending on δ and the characteristics of ϕ) such that for all P, Q ∈ R n×n and t ≥ 0.
2.5.Structural assumptions on the stochastic data.Let U be some separable Hilbert space and {u j } j∈N be a complete orthonormal system of U .We render the stochastic input via a cylindrical noise on the abstract space U .
Assumption 17 (Cylindrical Wiener process).We assume that W is an U -valued cylindrical Wiener process with respect to (F t ).Formally W can be represented as where {β j } j∈N are independent 1-dimensional standard Brownian motions.
We construct the noise coefficient G as a two parameter map.In its first parameter it acts as a Nemytskii operator whereas the dependence on the Wiener process is linear.
G is uniquely determined by the sequence of functions {g j }.In particular, regularity assumptions and summation properties need to be imposed on {g j } to obtain an operator G that is stable on specific function spaces.
The construction of the stochastic integral in the framework of Hilbert spaces can be done by the Itô isometry and requires the integrand to be an Hilbert-Schmidt operator, i.e., if . Moreover, we have the equivalence The construction of the stochastic integral in general Banach spaces is delicate and one needs to look for a suitable generalization of the Itô isometry.It turns out that γ-radonifying operators are natural in the context of stochastic integration.The following result was obtained by Ondreját and Veraar and contains optimal stability of the stochastic integral driven by a cylindrical Wiener process for γradonifying operators with values in a separable 2-smooth Banach space.For more details on γ-radonifying operators we refer to the survey [Nee10].
From now on we assume that W is given by (2.15) and G is of the form (2.16).
2.6.Gradient flow.The Helmholtz-Leray projection Π div allows to reformulate (1.1a) and (1.1b) into ) Now it is possible to interpret (2.20a) as a stochastically perturbed gradient flow of the energy induced by the potential ϕ κ (t) := t 0 (κ + s) p−2 s ds on the space W 1,p 0,div .It has been observed by Diening and Kreuzer [DK08] that the energy gap of In fact, the result can be generalized to energy gaps between arbitrary functions u, v ∈ W 1,p x .Lemma 20.Let u, v ∈ W 1,p x .Then Proof.Due to a Taylor expansion it holds where the first and second order Gateaux derivatives are given by This establishes the claim.
Clearly, if we insert a minimizer v = v * in (2.23), then we recover (2.22).However, the estimate (2.23) ensures that J is strongly convex.
In order to identify (2.20a) as a perturbed gradient flow, we need to find a pointwise representation of the gradient DJ (u).This establishes (2.27).
The second part follows by Riesz's representation theorem.
Remark 22.We want to emphasize that in general it is not possible to decouple the projection and the divergence, i.e., Π div div should be understood as one operator (defined by the duality relation (2.9)) rather than the composition of two.However, for smooth functions it is exactly the composition of the divergence and the Helmholtz-Leray projection.

Temporal regularity of stochastic pressure for weak solutions
The aim of this section is the discuss weak solutions to (1.1).In particular, we define the concept of weak solutions and present a sufficient condition on the stochastic data for the existence of a unique weak solution (u, π).Moreover, we split the pressure into two terms.One is related to the diffusion operator S(εu), while the second corresponds to the stochastic data.
3.1.Weak solutions.The concept of weak solutions is defined as follows.
and for all ξ t,x ∈ C ∞ 0 (O T ) and P-a.s.
3.2.A sufficient condition for weak solutions.
Assumption 24.We assume that Assumption 24 is standard for the derivation of weak solutions.Sometimes one couples the condition on the noise coefficient and the dissipation of the monotone operator S. Since we are only interested in regularity of π, we do not proceed this way.

3.3.
Existence of weak solutions.The general theory of monotone SPDEs of Liu and Röckner [LR10] covers the construction of the velocity variable.A similar construction for weak solutions to power-law fluids has been done in [Bre15].For the sake of completeness, we state the result and comment on the main ingredients.), O be open, bounded, Assumption 24 be satisfied and u 0 ∈ L q ω L 2 div be F 0 -measurable.Then there exists u The abstract conditions postulated by Liu and Röckner, that guarantee the existence of a weak solution, can be verified along the lines of [LR15, Example 4.1.9].In fact, [LR15, Example 4.1.9]only deals with the case p ≥ 2. However, the general case p ∈ (1, ∞) can be recovered, if we adjust the underlying Gelfand triple to The key tools for the a priori bound (3.5) are Itô's formula for F (u) = u q L 2 x together with a Gronwall argument.While the case q = 2 is standard, the general case needs some minor changes.Since this is not the main focus of the article, we only refer to [Wic22, Theorem 4.3.1]for more details.
Remark 26.The integrability in probability of the solution u is purely determined by the integrability in probability of the initial condition u 0 , i.e., for all q ∈ [1, ∞) x .The limit case q = ∞ needs to be excluded, since the moment transfer is hindered by the Wiener process W .The most one can hope for is More details about integrability of Brownian motions in Banach spaces can be found in [HV08].
3.4.Reconstruction of the pressure.Thanks to the Helmholtz-Leray projection Π div it is possible to target the construction of the velocity u and the pressure π individually.The pressure corresponds to the residual error when testing (3.1) by smooth gradients rather than smooth divergence free fields, i.e., we use (3.2) as a definition of the pressure π.Additionally, we artificially decompose π = π det + π sto into the components and discuss them separately.
First, we have a look at the deterministic component.
Lemma 27.Let O be a bounded C 2 domain and u ∈ W 1,p 0,x .Then there exists Since we want to identify π det as a proper function, we substitute ξ = Bξ 0 , where B denotes the Bogovskii operator (see Theorem 10) and ξ 0 ∈ C ∞ 0,x with ξ 0 O = 0, to obtain Thus, using the symmetry of Π ⊥ div and integration by parts, (3.10) is equivalent to where B * is the adjoint operator of B.
It remains to verify the boundedness of x .We instead address the boundedness of the adjoint operator and use a duality argument.
The right hand side of (3.10) is estimated by Hölder's inequality We further estimate, due to the Sobolev stability of the Helmholtz-Leray projection (Theorem 8) and the Bogovskii operator (Theorem 10 with s = 0, q = p), Therefore, the linear operator A : L p 0,x → L p 0,x is bounded.The claim (3.7) follows by our construction.The inequality (3.8) is verified by a density argument Second, we investigate the stochastic pressure.The main difficulty is the limited time regularity.Therefore, we initially take a look at the time integrated pressure.
Lemma 28.Let O be a bounded domain with locally Lipschitz boundary and Proof.The proof proceeds similar to the one of Lemma 27.
First, note (3.12) is equivalent to where B * is the adjoint of the Bogovskii operator.In particular, K sto is mean-value free.
Due to Theorem 10, one can check that B * : then by Hölder's inequality and Theorem 10 with s = −1 and q = 2, Thus, using the density of smooth mean-value free functions within W 1,2 0,x ∩ L 2 0,x , Finally, (3.15) together with L 2 x -stability of the Helmholtz-Leray projection establish An application of the Φ2,∞ -norm verifies (3.13).Now, we have identified the regularity class for the integrated pressure.The next step is the transfer of the results to the pressure.For a smooth test function To shorten the notation we abbreviate X = W 1,2 0,x ∩ L 2 0,x .
Lemma 29.Let α ∈ (0, 1) and p, q ∈ (1, ∞).Moreover, assume that x be spatially mean-value free.Due to (3.16) and duality Note that ξ trivially extends by zero to the full space in a smooth way.Now, we can use the equivalent spectral characterization of Besov norms, cf.Corollary 41, .
Since X is reflexive and p, q < ∞ we can identify the dual using Lemma 42 ) .An application of Lemma 43 yields Going back to the Besov norm in terms of integrability of differences, cf.Theorem 40, and using that ξ is extended by zero (3.20) Overall, revisiting (3.18) and using the density of smooth mean-value free functions in X together with (3.20), Ultimately, due to the equivalent norm on the dual space, cf.Lemma 42, Taking the square and expectation verifies (3.17).
Remark 30.So far the stochastic pressure was always estimated on suboptimal spaces.For example in [LRS03; CP12] the authors infer They neglect all additional information on the temporal differentiability of K sto and only use the boundedness.
Whether a corresponding result of Lemma 29 remains valid on the scale of exponentially integrable Besov spaces boils done to the understanding of the interaction between derivatives and norms in the spirit of (3.19), i.e., In order to establish (3.19), we heavily relied on the spectral representation of Besov norms.However, the equivalence of the spectral representation and the definition in terms of integrated differences fails for p ∈ {1, ∞}.Particularly, Φ 2 (t) = e t 2 − 1 is to close to ∞ and whence Φ * 2 is to close to 1.
Still, in regard of the regularity result of the integrated pressure, cf.Lemma 28, it is natural to define the set Corollary 31.In the setting of Lemma 28 we have Φ2,∞ X P-a.s.
Proof.We only need to verify that the distribution π sto defined by (3.16) is welldefined.The claim that π sto belongs to Φ2,∞ X follows by the definition (3.22) and the assumption on K sto .
Next, we take the square and expectation .
Thus, we have shown that π sto is well-defined.
3.5.Proof of Theorem 1.The previous results lead to a simple proof of Theorem 1.
Proof of Theorem 1. Due to the factorization (1.3) it is possible to construct the velocity field u independently of the pressure π.Afterwards the pressure is reconstructed.
The existence of the velocity field and verification of the a priori bound (1.6a) is done in Theorem 25.
In Section 3.4 we deal with the pressure reconstruction.Lemma 27 verifies the regularity of the deterministic pressure.The stochastic pressure is handled in Lemma 28 and 29.In particular, Lemma 28 ensures that X for any α ∈ (0, 1/2) and s, r ∈ (1, ∞).Therefore, Lemma 29 implies (1.6c).

Existence of strong solutions
Within this section we discuss the existence of strong solutions to (1.1).The main difference between weak and strong solutions lies in the fact that the latter can be interpreted in a point-wise manner, i.e., there is no need to interpret (1.1a) in a distributional sense.
Strong solutions for a related model have already been constructed in [BG19] by means of improved spatial regularity.They formally test (1.1) by the Laplacian.
We on the other hand rely on the gradient flow structure (2.20a), which -at least formally -corresponds to a test with −Π div div S(εu).

Strong solutions.
Definition 32.Let u 0 ∈ L 2 ω L 2 div be F 0 -measurable.The tupel (u, π) is called strong solution if it is a weak solution and additionally satisfies (a) Π div div S(εu) ∈ L 2 ω L 2 t L 2 div , (b) for all t ∈ I and P − a.s. it holds as an equation in L 2 div .In general, strong solutions only upgrade the regularity of the divergence-free component of the equation.In particular, the pressure does not enjoy higher regularity for strong solutions compared to weak solutions.

4.2.
A sufficient condition for strong solutions.
Assumption 33.We assume that Assumption 33 can be verified by using the Sobolev stability of Π div , cf.Theorem 8, and a suitable assumption on the data {g j }.

4.3.
Proof of Theorem 2. Our main tool in the derivation of energy bounds for strong solutions is the gradient flow structure of the equation.
Proof of Theorem 2. We will only address the case q = 1.General moments can be obtained by expanding the energy J (u) q .For more details we refer to [Wic22, Theorem 4.3.1].Our aim is to apply the result [Ges12,Theorem 1.4].Therefore, we need to check the assumptions (A1) − (A6).
First of all we relate our framework to the one used in [Ges12].We choose H = L 2 div , S = W 1,p 0,x and V = H ∩ S. Additionally, we let ϕ(v) = φ(v) = J (v) and B t (v) ≡ Π div G(v).Keep in mind that we use ϕ as the potential that defines the energy whereas Gess uses ϕ to denote the energy itself.
Ad (A1): First we show that J : W 1,p 0,x → R is continuous.Indeed, let u, v ∈ W 1,p 0,x .Then, due to the fundamental theorem, where This and Hölder's inequality imply Finally, we arrive at The continuity is verified.Moreover, using that ϕ (t) = t p−1 is homogeneous of degree p − 1 and a substitution, Since p ≥ 2 we find ϕ κ (t) is increasing in κ for fixed t and whence J (2u) ≤ 2 p J (u).Ad (A2): Recall Applying a weighted Young's inequality one finds In other words ϕ κ behaves like the pth-power shifted along κ.It follows using the monotonicity of powers and (4.3) Let {w k } k∈N ⊂ W 1,p 0,x .Now, using (2.25), Hölder's and Young's inequalities and (4.4) . Ad (A3): Lemma 20 shows that J is strongly convex.Ad (A4): Recall (2.24).Therefore, Ad (A5): Follows by the assumption on the noise coefficient, cf.(3.4).Ad (A6): The major ingredient that enables strong solutions is the regularity of the noise.Recall (2.16) and (A2) with An application of the Assumption 33 and (4.4) show x + 1 J (v) + 1.Overall, we have verified (A1) − (A6) and therefore we may apply [Ges12, Theorem 1.4] which provides us with the a priori bound x ds E [J (u 0 )] + 1.
A careful consideration of the proof shows that one can also obtain an estimate where the supremum in time is taken before the expectation.One additional term, due to the stochastic integral, needs to be considered, i.e., Due to the Burkholder-Davis-Gundy inequality Invoking Hölder's inequality The Assumption 33 together with (4.4) show Finally, Young's inequality verifies This allows for a Gronwall type argument that establishes the improved a priori estimate Remark 34.On a first view we could ask, if one could expand the q-energy of the symmetric p-Stokes system.However, additional difficulties arise to to the presence of the symmetric gradient.In particular, one needs to understand the interaction between projected symmetric p-and q-Laplacian.If one could establish a sign then it is possible to obtain improved gradient regularity.Regularity questions about div S(εu) -and ultimately on ∇S(εu) -are out of reach, since x is only a distribution.The derivation of higher regularity for the pressure is a non-trivial task.

Temporal regularity of strong solutions
Fortunately, the strong formulation of equation (1.1) is sufficient to derive improved temporal regularity.We already established a similar result for the p-Laplace system in [Wic21].
5.1.Exponential Besov-regularity.In general, it is delicate to derive temporal regularity on the limited threshold of order 1/2 for stochastic partial differential equations.A key example was derived by Hytönen and Veraar in [HV08].They show that P-almost surely Wiener processes belong to the Besov space B 1/2 Φ2,∞ and do not belong to the Besov space B 1/2 p,q for p, q < ∞.We show that strong solutions inherit similar temporal regularity properties from the driving Wiener process.
Proof of Theorem 3. We will only present part (b) of Theorem 3 in detail and comment on the necessary changes for part (a).
We start with the estimate for the Besov-Orlicz semi-norm Using the strong formulation (4.1) Invoking L ∞ → L Φ2 and Hölder's inequality The second term is bounded by Theorem 19 (e), (2.19) and (3.3) This establishes the semi-norm estimate.
The full norm follows by the embedding x + 1.We have verified (1.9).
Part (a) follows similarly to part (b), but we use the embedding L ∞ t → L q t and Theorem 19 (b) together with an extrapolation result, cf.[OV20, Remark 3.4].
Remark 35.Theorem 3 immediately implies Hölder regularity up to almost 1/2 by an application of a Sobolev embedding (cf.[Tri83, Section 2.3.2]),i.e., for any α ∈ (0, 1/2) x .5.2.Nikolskii-regularity of non-linear gradient.Before we can turn our attention to investigate the regularity properties of the non-linear gradient, we need to understand how the stochastic integral operator interacts with linear projections and gradients.
Lemma 36.Let p ≥ 2 and Assumption 33 be satisfied.Then there exists a constant C > 0 such that for all progressively measurable u ∈ L p ω L ∞ t L 2 div ∩ W 1,p 0,x it holds  εΠ div G(u) Now we are ready to prove the result on temporal regularity of the non-linear gradient.
Remark 37. Strong solutions to the symmetric stochastic p-Stokes system (1.1) are comparable in terms of temporal regularity to the stochastic p-Laplace system, cf.[Wic21].
In particular, if one could prove increased integrability, e.g.Π div div S(εu) ∈ L 2 ω L q t L 2 x for some q > 2, then it would also be possible to show V (εu) ∈ L 2 ω B 1/2 q,∞ L 2 x .The temporal regularity enables a similar error analysis for time discretizations of (1.1) as done in e.g.[BHL21; DHW22] for the p-Laplace system.

Appendix A. Fourier representation of Besov norms
There are many different characterizations of Besov spaces.A nice overview can be found e.g. in the books of Triebel [Tri83;Tri92].Some results derived from a specific characterization are not obvious for others.In this appendix we introduce the Fourier representation of Besov spaces and state some properties.
8) Theorem 8 can be proven using the fundamental solution of the Neumannproblem (2.7) and Calderon-Zygmund estimates.For example the result for Besselpotential spaces is presented in [Mik02, Lemma 3.6 & Lemma 3.8].Corollary 9. Let O be a bounded C 2 -domain and q

=
Due to the linearity of εΠ div and I we findεΠ div I G(u) = I εΠ div G(u) .Therefore, using the time regularity of stochastic integrals, cf.Theorem 19 (b), εΠ div I(G(u)) I(εΠ div G(u))
Since L p x is of 2-smooth, we can use (2.19).This and the growth assumption (4.2) establish εΠ div G(u) x ) .