1 Introduction

Let \({\mathcal {O}} \subset {\mathbb {R}}^n\) be a bounded domain, \(n \in {\mathbb {N}}\) and \(T > 0\) be finite. Given an initial datum \(u_0\) and a stochastic forcing \(G(\cdot ) \,\textrm{d}W\), we are seeking for a velocity field \(u: \Omega \times (0,T) \times {\mathcal {O}} \rightarrow {\mathbb {R}}^n\) and a pressure \(\pi : \Omega \times (0,T) \times {\mathcal {O}} \rightarrow {\mathbb {R}}\) that satisfy the relations

$$\begin{aligned} \,\textrm{d}u + \big ( - \text {div}S(\mathbb {\epsilon } u) + \nabla \pi \big )\,\textrm{d}t&= G(u) \,\textrm{d}W(t) \quad{} & {} \text { in } \Omega \times (0,T) \times {\mathcal {O}}, \end{aligned}$$
(1.1a)
$$\begin{aligned} \text {div}u&= 0 \quad{} & {} \text { in } \Omega \times (0,T) \times {\mathcal {O}}, \end{aligned}$$
(1.1b)
$$\begin{aligned} u&= 0 \quad{} & {} \text { on } \Omega \times (0,T) \times \partial {\mathcal {O}}, \end{aligned}$$
(1.1c)
$$\begin{aligned} u(0)&= u_0{} & {} \text { on } \Omega \times {\mathcal {O}}, \end{aligned}$$
(1.1d)
$$\begin{aligned} \langle \pi \rangle _{{\mathcal {O}}}&= 0{} & {} \text { in } \Omega \times (0,T), \end{aligned}$$
(1.1e)

where \(S(\xi ):= \left( \kappa + \left| \xi \right| \right) ^{p-2} \xi \), \(\xi \in {\mathbb {R}}^{ n \times n}\), \(p \in (1,\infty )\), \(\kappa \ge 0\), \(\mathbb {\epsilon } u:= 2^{-1}\big ( \nabla u + (\nabla u)^T \big )\) and . The stochastic input \(G(\cdot ) \,\textrm{d}W\) is rendered via a cylindrical Wiener process W on an abstract separable Hilbert space U and a suitable noise coefficient G.

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

$$\begin{aligned} \,\textrm{d}u + \big ( (u\cdot \nabla )u - \text {div}S(\mathbb {\epsilon } u) + \nabla \pi \big )\,\textrm{d}t= G(u) \,\textrm{d}W(t) \end{aligned}$$
(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. [50] (see also [12] 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:

  1. (a)

    Solutions to non-linear equations lack smoothness even for smooth data.

  2. (b)

    Limited regularity of the stochastic data \(G(\cdot ) \,\textrm{d}W\) leads to limited regularity of the solution.

  3. (c)

    Stochastic integration in general Banach spaces causes technical restrictions.

  4. (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 [39, 40] and Lions [43]. Motivated by the works, many authors contributed to the construction of solutions, e.g. [28, 46,47,48, 62]. Ultimately leading to the construction of a weak solution for \(p > 2n/(n+2)\). Below the compactness threshold \(W^{1,p}({\mathcal {O}}) \hookrightarrow \hookrightarrow L^2({\mathcal {O}})\) wild solutions may emerge, see e.g. [15].

The stochastic case is less known. An answer for solution independent noise coefficient was given in [57, 63]. The authors require the condition \(p > 9/5\) in the three dimensional setting. It was extended in [10] 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 \(\Omega \times (0,T) \times {\mathcal {O}}\), is done in [13].

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. [11]. Therefore, and out of mathematical curiosity, many authors have investigated the regularity of solutions to (1.1) and related models.

1.2.1 Deterministic

In general, there is a barrier that does not allow for arbitrary high regularity even though the data is smooth, cf. [37]. Still partial results can be obtained.

\(C^{1,\alpha }\)-regularity has been obtained in e.g. [2, 19, 20, 29]. Hölder regularity is sub-optimal from a numerical point of view. Instead, regularity of the expression \(V(\mathbb {\epsilon } u) = (\kappa + \left| \mathbb {\epsilon } u\right| )^{(p-2)/2} \mathbb {\epsilon } u\) is more suitable for the investigation of numerical algorithms as observed by Barrett and Liu in [4]. Indeed, the operator V measures the monotonicity of S in a natural way, i.e., for all \(A,B \in {\mathbb {R}}^{n\times n}\) it holds

$$\begin{aligned} \big ( S(A) - S(B) \big ) : (A- B) \eqsim \left| V(A) - V(B)\right| ^2. \end{aligned}$$

A first result for non-singular shear thinning fluids, \(\kappa = 1\) and \(p \in (1,2)\), has been obtained in [51]. The author verifies local regularity \(V(\mathbb {\epsilon } u) \in W^{1,2}_{\textrm{loc}}({\mathcal {O}})\).

A similar result has been obtained for shear thickening fluids, \(p \ge 2\), in [9]. Additionally, local regularity for the non-linear tensor S and the pressure \(\pi \) was found in [38]. They show \(\nabla S(\mathbb {\epsilon } u), \nabla \pi \in L^{p'}_{\textrm{loc}}({\mathcal {O}})\).

In [21] 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 [1]. The authors show that \(\nabla u \in L^\infty (0,T; L^2_{\textrm{loc}}({\mathcal {O}}))\) and \(\nabla V(\mathbb {\epsilon } u) \in L^2_{\textrm{loc}}((0,T)\times {\mathcal {O}})\). It was extended to the \(\phi \)-Stokes problem in [17].

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 [33] and by Burczak and Kaplický in [16]. In [33] it is shown that \(\partial _t u \in B^{1/2}_{2,\infty }(0,T; L^2({\mathcal {O}}))\). Here B denotes a Besov space (see Sect. 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 [35] 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 \(\Pi _{\text {div}}\) (see (A. 6)) allows to reformulate (1.1a) and (1.1b) into

$$\begin{aligned} \,\textrm{d}u - \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \,\textrm{d}t&= \Pi _{\text {div}} G(u) \,\textrm{d}W(t) \quad{} & {} \text { in } \Omega \times (0,T) \times {\mathcal {O}}, \end{aligned}$$
(1.3a)
$$\begin{aligned} \nabla \pi \,\textrm{d}t&= \Pi _{\text {div}}^\perp \big ( \text {div}S(\mathbb {\epsilon } u) \,\textrm{d}t+ G(u) \,\textrm{d}W(t) \big )\quad{} & {} \text { in } \Omega \times (0,T) \times {\mathcal {O}}. \end{aligned}$$
(1.3b)

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 (A. 19a) as a stochastically perturbed gradient flow of the energy

$$\begin{aligned} {\mathcal {J}}(u) := \int _{\mathcal {O}} \varphi _{\kappa } ( \left| \mathbb {\epsilon } u\right| ) \,\textrm{d}x \end{aligned}$$
(1.4)

induced by the potential \(\varphi _{\kappa }(t):= \int _0^t (\kappa + s)^{p-2} s \,\textrm{d}s\) on the space \(W^{1,p}_{0,\text {div}}({\mathcal {O}})\). Indeed, (A. 19a) can be written as

$$\begin{aligned} \,\textrm{d}u = - D {\mathcal {J}}(u) \,\textrm{d}t + \Pi _{\text {div}} G(u) \,\textrm{d}W(t) \quad \text { in } \Omega \times (0,T) \times {\mathcal {O}}, \end{aligned}$$
(1.5)

where \(D {\mathcal {J}}(u) = \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } 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 [13]. They prove spatial regularity as in the deterministic case also for stochastically forced electro-rheological fluids and obtain \(\nabla u \in L^2(\Omega ; L^\infty (0,T; L^2({\mathcal {O}})))\) as well as \(\nabla V(\mathbb {\epsilon } u) \in L^2(\Omega \times (0,T) \times {\mathcal {O}})\). 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(\cdot ) \,\textrm{d}W\) on the solution u.

Even in the simplest case \(G\equiv 1\), we don’t expect u to exceed the regularity threshold induced by the Wiener process W. A precise regularity result was derived by Hytönen and Veraar in [36]. They show that \({\mathbb {P}}\)-a.s. \(W \in B^{1/2}_{\Phi _2,\infty }(0,T)\), where \(\Phi _2(t) = \exp (t^2)-1\), and \(W \not \in B^{1/2}_{p,q}(0,T)\) for any \(p \in [1,\infty ]\) and \(q \in [1,\infty )\). 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 [53]. 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.,

$$\begin{aligned} G \in L^\infty (0,T) \quad \Rightarrow \quad \int _0^{\cdot } G \,\textrm{d}W(t) \in B^{1/2}_{\Phi _2,\infty }(0,T). \end{aligned}$$

Exactly this fact is used in [61] to derive temporal regularity for the stochastic p-Laplace system. In particular, we manage to recover the limiting temporal regularity \(u \in L^2(\Omega ; B^{1/2}_{\Phi _2,\infty }(0,T; L^2({\mathcal {O}})))\). 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 \(V(\nabla u) \in L^2(\Omega ; B^{1/2}_{2,\infty }(0,T; L^2({\mathcal {O}})))\).

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(\cdot ) \,\textrm{d}W\) is not divergence-free, then the pressure is very rough in time (see e.g. [18, 41]). This can be seen by formally decomposing the pressure into a deterministic and a stochastic component \(\pi = \pi _{\textrm{det}} + \pi _{\textrm{sto}}\), where

$$\begin{aligned} \nabla \pi _{\textrm{det}}&= \Pi _{\textrm{div}}^{\perp } \text {div}S(\mathbb {\epsilon } u), \\ \nabla \pi _{\textrm{sto}}&= \Pi _{\textrm{div}}^{\perp } G(u) {\dot{W}}. \end{aligned}$$

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 [30, 31, 42] and the references therein.

1.3 Outline

In Sect. 2 we introduce the setup and present our main results.

In Sect. 3 we construct weak solutions and discuss regularity of pressure.

In Sect. 4 we prove the existence of strong solutions and verify temporal regularity of velocity and non-linear gradient.

More details on Besov norms, Helmholtz decomposition, Bogovskii operator, potentials, mapping properties of stochastic integrals, and gradient flows are presented in Appendix A.

2 Setup and Main Results

Let \({\mathcal {O}} \subset {\mathbb {R}}^n\), \(n \in {\mathbb {N}}\), be a bounded Lipschitz domain. For some given \(T>0\) we denote by \(I:= (0,T)\) the time interval and write \({\mathcal {O}}_T:= I \times {\mathcal {O}}\) for the time space cylinder. Moreover let \(\left( \Omega ,{\mathcal {F}}, ({\mathcal {F}}_t)_{t\in I}, {\mathbb {P}} \right) \) denote a stochastic basis, i.e., a probability space with a complete and right continuous filtration \(({\mathcal {F}}_t)_{t\in I}\). We write \(f \lesssim g\) for two non-negative quantities f and g if f is bounded by g up to a multiplicative constant. Accordingly we define \( > rsim \) and \(\eqsim \). We denote by c a generic constant which can change its value from line to line. Inner products and duality pairings are denoted by \(\left( \cdot , \cdot \right) \) and \(\langle \cdot , \cdot \rangle \), respectively.

2.1 Function Spaces

As usual, for \(1\le q < \infty \) we denote by \(L^q({\mathcal {O}})\) the Lebesgue space and \(W^{1,q}({\mathcal {O}})\) the Sobolev space. Moreover, \(W^{1,q}_0({\mathcal {O}})\) denotes the Sobolev spaces with zero boundary values. It is the closure of \(C^\infty _0({\mathcal {O}})\) (smooth functions with compact support) in the \(W^{1,q}({\mathcal {O}})\)-norm. We denote by \(W^{-1,q'}({\mathcal {O}})\) the dual of \(W^{1,q}_0({\mathcal {O}})\). The space of mean-value free Lebesgue functions is denoted by \(L^q_0({\mathcal {O}})\). The space of smooth, compactly supported and divergence-free vector fields is called \(C^\infty _{0,\text {div}}({\mathcal {O}})\) and its closure within the \(W^{1,p}\)-norm is abbreviated by \(W^{1,p}_{0,\text {div}}({\mathcal {O}})\). We do not distinguish in the notation between scalar-, vector- and matrix-valued functions.

For a Banach space \(\left( X, \left\| \cdot \right\| _X \right) \) let \(L^q(I;X)\) be the Bochner space of Bochner-measurable functions \(u: I \rightarrow X\) satisfying \(t \mapsto \left\| u(t)\right\| _X \in L^q(I)\). Moreover, \(C^0({\overline{I}};X)\) is the space of continuous functions with respect to the norm-topology. We also use \(C^{0,\alpha }({\overline{I}};X)\) for the space of \(\alpha \)-Hölder continuous functions. Given an Orlicz-function \(\Phi : [0,\infty ] \rightarrow [0,\infty ]\), i.e. a convex function satisfying \( \lim _{t \rightarrow 0} \Phi (t)/t = 0\) and \(\lim _{t \rightarrow \infty } \Phi (t)/t = \infty \) we define the Luxemburg-norm

$$\begin{aligned} \left\| u\right\| _{L^\Phi (I;X)} := \inf \left\{ \lambda > 0 : \int _I \Phi \left( \frac{\left\| u\right\| _X}{\lambda } \right) \,\textrm{d}s\le 1 \right\} . \end{aligned}$$

The Orlicz space \(L^\Phi (I;X)\) is the space of all Bochner-measurable functions with finite Luxemburg-norm. For more details on Orlicz-spaces we refer to [25]. Given \(h \in I\) and \(u:I \rightarrow X\) we define the difference operator \(\tau _h: {\{{u: I \rightarrow X}\}} \rightarrow {\{{u: I\cap I - {\{{h}\}} \rightarrow X}\}} \) via \(\tau _h(u) (s):= u(s+h) - u(s)\). The Besov-Orlicz space \(B^\alpha _{\Phi ,r}(I;X)\) with differentiability \(\alpha \in (0,1)\), integrability \(\Phi \) and fine index \(r \in [1,\infty ]\) is defined as the space of Bochner-measurable functions with finite Besov-Orlicz norm \(\left\| \cdot \right\| _{B^\alpha _{\Phi ,r}(I;X)}\), where

$$\begin{aligned} \begin{aligned} \left\| u\right\| _{B^\alpha _{\Phi ,r}(I;X)}&:= \left\| u\right\| _{L^{\Phi }(I ;X)} + {[{u}]}_{B^{\alpha }_{\Phi ,r}(I;X)}, \\ {[{u}]}_{B^\alpha _{\Phi ,r}(I;X)}&:= {\left\{ \begin{array}{ll} \left( \int _{I} h^{-r\alpha } \left\| \tau _h u\right\| _{L^\Phi (I\cap I - {\{{h}\}};X)}^r \frac{\,\textrm{d}h}{h} \right) ^\frac{1}{r} &{} \text { if } r \in [1,\infty ), \\ \mathop {\mathrm {ess\,sup}}_{h \in I} h^{-\alpha } \left\| \tau _h u\right\| _{L^\Phi (I\cap I - {\{{h}\}};X)} &{} \text { if } r = \infty . \end{array}\right. } \end{aligned} \end{aligned}$$
(2.1)

The case \(r = \infty \) is commonly called Nikolskii-Orlicz space and abbreviated by \(N^{\alpha ,\Phi } = B^{\alpha }_{\Phi ,\infty }\). When \(\Phi (t) = t^p\) we write \(B^\alpha _{p,r}(I;X)\) and call it Besov space. If X is reflexive, we denote by \(B^{-\alpha }_{p',q'}(I;X') = \big ( B^{\alpha }_{p,q}(I;X) \big )'\).

Similarly, given a Banach space \(\left( Y, \left\| \cdot \right\| _Y \right) \), we define \(L^q(\Omega ;Y)\) as the Bochner space of Bochner-measurable functions \(u: \Omega \rightarrow Y\) satisfying \(\omega \mapsto \left\| u(\omega )\right\| _Y \in L^q(\Omega )\). The space \(L^q_{{\mathcal {F}}}(\Omega \times I;X)\) denotes the subspace of X-valued progressively measurable processes. We abbreviate the notation \(L^q_\omega L^q_t L^q_x:= L^q(\Omega ;L^q(I;L^q({\mathcal {O}}))) \) and \(L^{q-}:= \bigcap _{r< q} L^r\).

2.2 Structural Assumptions on the Stochastic Data

Let U be some separable Hilbert space and \({\{{u_j}\}}_{j\in {\mathbb {N}}}\) be a complete orthonormal system of U. We render the stochastic input via a cylindrical noise on the abstract space U.

Assumption 1

(Cylindrical Wiener process). We assume that W is an U-valued cylindrical Wiener process with respect to \(({\mathcal {F}}_t)\). Formally W can be represented as

$$\begin{aligned} W = \sum _{j \in {\mathbb {N}}} u_j \beta ^j, \end{aligned}$$
(2.2)

where \({\{{\beta ^j}\}}_{j\in {\mathbb {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.

Definition 2

(Noise coefficient). Let \({\{{g_j}\}}_{j\in {\mathbb {N}}}: {\mathcal {O}} \times {\mathbb {R}}^N \rightarrow {\mathbb {R}}^N\). We define G by

$$\begin{aligned} \begin{aligned} G: L^2_{{\mathcal {F}}}(\Omega \times I; L^2_x) \times U&\rightarrow L^2_{{\mathcal {F}}}(\Omega \times I; L^2_x) , \\ (v,u)&\mapsto G(v)u := \sum _{j \in {\mathbb {N}}} g_j\big (\cdot ,v(\cdot )\big )( u_j, u)_U. \end{aligned} \end{aligned}$$
(2.3)

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 Itôs isometry and requires the integrand to be an Hilbert-Schmidt operator, i.e., if \(G \in L^2_{{\mathcal {F}}}\big (\Omega \times I; L_2(U; L^2_x) \big )\), then

$$\begin{aligned} {\mathcal {I}}_W(G) := \int _0^{\cdot } G \,\textrm{d}W(t) := \sum _{j\in {\mathbb {N}}} \int _0^{\cdot } G(u_j) \,\textrm{d}\beta ^j(t) \end{aligned}$$
(2.4)

converges in \(L^2\big (\Omega ; C(I; L^2_x) \big )\). Moreover, we have the equivalence

$$\begin{aligned} {\mathbb {E}} \left[ \sup _{t \in I} \left\| {\mathcal {I}}_W(G)\right\| _{L^2_x}^2 \right] \eqsim {\mathbb {E}} \left[ \int _I \left\| G\right\| _{L_2(U;L^2_x)}^2 \,\textrm{d}t\right] . \end{aligned}$$
(2.5)

2.3 Main Result I: Temporal Regularity of Pressure for Weak Solutions

We apply the general theory of monotone stochastic partial differential equations developed by Liu and Röckner [44], see also the book [45], to the p-Stokes system. More recently, the theory has been refined to cover more general equations, see e.g. [3, 55]. In this way, existence of weak solutions is almost immediate. The main novelty is a careful consideration of the regularity of the stochastic pressure.

Definition 3

Let \(u_0 \in L^2_\omega L^2_{\text {div}}\) be \({\mathcal {F}}_0\)-measurable. The tupel \((u,\pi )\) is called weak solution to (1.1) if

  1. (a)

    \(u \in L^2_\omega C_t^0 L^2_{\text {div}} \cap L^p_\omega L^p_t W^{1,p}_{0,x}\) is \(({\mathcal {F}}_t)\)-adapted,

  2. (b)

    for all \(t \in I\), \(\xi \in C^\infty _{0,\text {div}}\) and \({\mathbb {P}}\)-a.s. it holds

    $$\begin{aligned} \int _{{\mathcal {O}}} (u(t) - u_0) \cdot \xi \,\textrm{d}x+ \int _0^t \int _{{\mathcal {O}}} S( \mathbb {\epsilon } u) :\nabla \xi \,\textrm{d}x\,\textrm{d}s= \sum _{j\in {\mathbb {N}}} \int _0^t \int _{{\mathcal {O}}} g_j(\cdot ,u)\cdot \xi \,\textrm{d}x\,\textrm{d}\beta ^j(s), \end{aligned}$$
    (2.6)
  3. (c)

    \(\pi = \pi _{\textrm{det}} + \pi _{\textrm{sto}}\) with \(\pi _{\textrm{det}} \in L^{p'}_\omega L^{p'}_t L^{p'}_{0,x}\), \(\pi _{\textrm{sto}} \in L^{2}_\omega B^{-1/2}_{\Phi _2,\infty } \big ( W^{1,2}_x \cap L^2_{0,x} \big )\)

  4. (d)

    and for all \(\xi _{t,x} \in C^\infty _0({\mathcal {O}}_T)\) and \({\mathbb {P}}\)-a.s.

    $$\begin{aligned} \langle \pi , \text {div}\xi _{t,x} \rangle = \int _{{\mathcal {O}}_T} S(\mathbb {\epsilon } u): \nabla \Pi _{\text {div}}^\perp \xi _{t,x} \,\textrm{d}x \,\textrm{d}t +\int _{{\mathcal {O}}_T} {\mathcal {I}}_W\big (G(u)\big ) \partial _t \Pi _{\text {div}}^\perp \xi _{t,x} \,\textrm{d}x \,\textrm{d}t. \end{aligned}$$
    (2.7)

Here \(\Pi _{\text {div}}^\perp \) is the complement of the Helmholtz–Leray projection and is defined in A. 6.

The general theory of monotone stochastic partial differential equations requires further assumptions on the noise coefficient.

Assumption 4

We assume that \({\{{g_j}\}}_{j\in {\mathbb {N}}} \in C^0({\mathcal {O}} \times {\mathbb {R}}^n; {\mathbb {R}}^n)\) with

  1. (a)

    (sublinear growth) for all \(v \in L^2_{\text {div}}\)

    $$\begin{aligned} \left\| G(v)\right\| _{L_2(U;L^2_x)}^2 = \sum _{j\in {\mathbb {N}}} \left\| g_j\big (\cdot ,v(\cdot )\big )\right\| _{L^2_x}^2 \le c_{\text {growth}}(1+\left\| v\right\| _{L^2_x}^2), \end{aligned}$$
    (2.8)
  2. (b)

    (Lipschitz continuity) for all \(v_1, \, v_2 \in L^2_{\text {div}}\) it holds

    $$\begin{aligned} \left\| G(v_1) - G(v_2)\right\| _{L_2(U;L^2_x)}^2 \le c_{\text {lip}} \left\| u_1 - u_2\right\| _{L^2_x}^2. \end{aligned}$$
    (2.9)

We are ready to formulate our first main result.

Theorem 5

Let \(p \in (1,\infty )\) and \({\mathcal {O}}\) be a bounded \(C^2\)-domain. Additionally, let \(q \in [1,\infty )\), \(u_0 \in L^{q}_\omega L^2_{\text {div}}\) be \({\mathcal {F}}_0\)-measurable and Assumption 4 be satisfied. Then there exists a unique weak solution \((u,\pi )\) to (1.1) such that

$$\begin{aligned} {\mathbb {E}} \left[ \sup _{t\in I} \left\| u(t)\right\| _{L^2_x}^{q} + \left( \int _I {\mathcal {J}}\big (u(t) \big ) \,\textrm{d}t \right) ^{q/2} \right]&\lesssim {\mathbb {E}} \left[ \left\| u_0\right\| _{L^2_x}^{q} \right] + 1, \end{aligned}$$
(2.10a)

where \({\mathcal {J}}\) is given by (1.4).

Moreover, it holds \(\pi = \pi _{\textrm{det}} + \pi _{\textrm{sto}} \) with

$$\begin{aligned} \left\| \pi _{\textrm{det}}\right\| _{L^{p'q /2}_\omega L^{p'}_t L^{p'}_{x}}^{p' /2}&\lesssim \left\| u_0\right\| _{L^q_\omega L^2_x}+ 1, \end{aligned}$$
(2.10b)
$$\begin{aligned} \left\| \pi _{\textrm{sto}}\right\| _{L^{q}_\omega B^{-1/2 - \varepsilon }_{s,r} W^{1,2}_x}&\lesssim \left\| u_0\right\| _{L^q_\omega L^2_x} + 1, \end{aligned}$$
(2.10c)

for any \(s, r \in (1,\infty )\) and \(\varepsilon >0\).

2.4 Main result II: 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 [35]. In doing so, we provide a sufficient condition for the existence of strong solutions to the p-Stokes system.

Definition 6

Let \(u_0 \in L^2_\omega L^2_{\text {div}}\) be \({\mathcal {F}}_0\)-measurable. The tupel \((u,\pi )\) is called strong solution if it is a weak solution and additionally satisfies

  1. (a)

    \(\Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \in L^2_\omega L^2_t L^2_{\text {div}}\),

  2. (b)

    for all \(t \in I\) and \({\mathbb {P}}-a.s.\) it holds

    $$\begin{aligned} u(t) - u_0 - \int _0^t \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \,\textrm{d}s= \sum _{j\in {\mathbb {N}}} \int _0^t \Pi _{\text {div}} g_j(\cdot ,u) \,\textrm{d}\beta ^j(s) \end{aligned}$$
    (2.11)

    as an equation in \(L^2_{\text {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. Additional assumptions need to be imposed on the noise coefficient to ensure existence of strong solutions.

Assumption 7

We assume that \({\{{g_j}\}}_{j\in {\mathbb {N}}} \in C^1({\mathcal {O}} \times {\mathbb {R}}^n; {\mathbb {R}}^n)\) such that for all \(v \in W^{1,p}_{0,x} \cap L^2_{\text {div}}\) it holds

  1. (a)

    (stability)

    $$\begin{aligned} \left\| \mathbb {\epsilon } \Pi _{\text {div}} G(v) \right\| _{L_2(U;L^p_x)} \lesssim 1+\left\| \mathbb {\epsilon } v\right\| _{L^p_x}, \end{aligned}$$
    (2.12)
  2. (b)

    (boundary condition) and \(\Pi _{\text {div}} G(v) \in W^{1,p}_{0,x}\).

Clearly, \(G(v) = v\) satisfies Assumption 7. Generally, the stability can be verified using the Sobolev stability of \(\Pi _{\text {div}}\), cf. Theorem 32, and a suitable assumption on the data \({\{{g_j}\}}\), see Example 8. But it is non-trivial to verify the boundary condition; the Helmoltz–Leray projection is a non-local operator. So even if G(v) vanishes on the boundary, it is unclear if \(\Pi _{\text {div}} G(v)\) also vanishes.

Example 8

Let \(h \in C^\infty _{0,x}\). Define \(g_1(x,v):= \sin (x\cdot v) h(x)\) and \(g_j = 0\) for \(j \ge 2\). Set \(G(v)(x):= \sum _{j\in {\mathbb {N}}} g_j(x,v(x)) u_j\) for \(v \in W^{1,p}_{0,x} \cap L^2_{\text {div}}\), \(\{u_j\}_{j \in {\mathbb {N}}}\) an ONB of U.

Notice that the i-th derivative of the k-th component is given by

$$\begin{aligned} \partial _{x_i} (g_1)_k(v)(x) = \sin (x \cdot v(x)) \partial _{x_i} h_k(x) + \cos (x \cdot v(x))h_k(x)( v_i(x) + \sum _{l} x_l \partial _{x_i} v_l(x) ). \end{aligned}$$

This implies \(G(v) = 0\) and \(\nabla G(v) = 0\) on \(\partial {\mathcal {O}}\). Now, using the Sobolev stability of \(\Pi _{\text {div}}\), Theorem 32,

$$\begin{aligned} \left\| \epsilon \Pi _{\text {div}} g_1(\cdot ,v)\right\| _{L^p_x} \le \left\| \Pi _{\text {div}} g_1(\cdot ,v)\right\| _{W^{1,p}_x} \lesssim \left\| g_1(\cdot ,v)\right\| _{W^{1,p}_x}. \end{aligned}$$

Moreover, due to Poincaré’s inequality, boundedness of \({\mathcal {O}}\) and Korn’s inequality,

$$\begin{aligned} \left\| g_1(\cdot ,v)\right\| _{W^{1,p}_x} \lesssim \left\| \nabla g_1(\cdot ,v)\right\| _{L^p_x} \lesssim \left\| \nabla h\right\| _{L^p_x} + \left\| h\right\| _{L^\infty _x} \left\| \epsilon v\right\| _{L^p_x}. \end{aligned}$$

This allows us to establish (2.12).

Theorem 9

Let \(p \ge 2\), \(q \in [1,\infty )\) and \( {\mathcal {J}}(u_0) \in L^{q}_\omega \) be \({\mathcal {F}}_0\)-measurable. Moreover, let G satisfy Assumptions 4 and 7. Then there exists a unique strong solution \((u,\pi )\) to (1.1). Moreover, it holds

$$\begin{aligned} {\mathbb {E}} \left[ \left( \sup _{t \in I} {\mathcal {J}}\big ( u(t) \big ) + \int _0^T \int _{{\mathcal {O}}} \left| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right| ^2 \,\textrm{d}x \,\textrm{d}t \right) ^{q }\right] \lesssim {\mathbb {E}} \left[ {\mathcal {J}}\big ( u_0 \big )^{q} \right] + 1. \end{aligned}$$
(2.13)

We extend the techniques developed in [61] to obtain improved temporal regularity for the velocity field and its non-linear symmetric gradient.

Theorem 10

Let \(p \in (1,\infty )\), Assumption 4 be satisfied and \((u,\pi )\) be a strong solution to (1.1).

  1. (a)

    (Polynomial integrability in time) Let \(q \in [1,\infty )\). Then \(u \in L^{2}_\omega B^{1/2}_{q,\infty } L^2_x\) and

    $$\begin{aligned} \left\| u\right\| _{L^{2}_\omega B^{1/2}_{q,\infty } L^2_x} \lesssim \left\| u\right\| _{L^{2}_\omega L^\infty _t L^2_x} + \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_\omega L^2_t L^2_x} + 1. \end{aligned}$$
    (2.14)
  2. (b)

    (Exponential integrability in time) Let \(N_2(t):= t^2 \ln (t+1)\), \(\Phi _2(t):= e^{t^2} - 1\) and \(u \in L^{N_2}_\omega L^\infty _t L^2_x\). Then \(u \in L^{2}_\omega B^{1/2}_{\Phi _2,\infty } L^2_x\) and

    $$\begin{aligned} \left\| u\right\| _{L^{2}_\omega B^{1/2}_{\Phi _2,\infty } L^2_x} \lesssim \left\| u\right\| _{L^{N_2}_\omega L^\infty _t L^2_x} + \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_\omega L^2_t L^2_x} + 1. \end{aligned}$$
    (2.15)

Remark 11

Theorem 10 works for the full range \(p \in (1,\infty )\). But the existence of strong solutions is still open for \(p \in (1,2)\). Existence of strong solutions for deterministic power-law fluids including convection has been established in [6].

Exponential integrability requires slightly stronger moment estimates. In a similar way to the a priori bound (2.10a) for weak solutions, one can establish an estimate of the form

$$\begin{aligned} \left\| u\right\| _{L^{N_2}_{\omega } L^\infty _t L^2_x} \lesssim \left\| u_0\right\| _{L^{N_2}_{\omega } L^2_x} + 1. \end{aligned}$$

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 12

Let \(p \ge 2\), Assumption 7 be satisfied and \((u,\pi )\) be a strong solution to (1.1). Additionally, assume \({\mathcal {J}}(u) \in L^1_\omega L^\infty _t\). Then \(V(\mathbb {\epsilon } u) \in L^2_\omega B^{1/2}_{2,\infty } L^2_x\) and

$$\begin{aligned} \left\| V(\mathbb {\epsilon } u)\right\| _{L^2_\omega B^{1/2}_{2,\infty } L^2_x}^2 \lesssim {\mathbb {E}} \left[ \int _{I} \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 \,\textrm{d}t\right] + {\mathbb {E}} \left[ \sup _{t\in I} {\mathcal {J}}\big (u(t)\big ) \right] + 1. \end{aligned}$$
(2.16)

3 Temporal Regularity of Stochastic Pressure for Weak Solutions

The aim of this section is the discuss weak solutions to (1.1). In particular, we split the pressure into two terms; one is related to the diffusion operator \(S(\mathbb {\epsilon } u)\), while the other corresponds to stochastic data.

3.1 Existence of Weak Solutions

The general theory of monotone SPDEs of Liu and Röckner [44] covers the construction of the velocity variable. A similar construction for weak solutions to power-law fluids has been done in [10]. For the sake of completeness, we state the result and comment on the main ingredients.

Theorem 13

Let \(p \in (1,\infty )\), \(q \in [1,\infty )\), \({\mathcal {O}}\) be open, bounded, Assumption 4 be satisfied and \(u_0 \in L^q_\omega L^2_{\text {div}}\) be \({\mathcal {F}}_0\)-measurable.

Then there exists \(u \in L^q_\omega C^0_t L^2_{\text {div}} \cap L^{p q/2}_\omega L^{p q/2}_t W^{1,p}_{0,x} \) adapted to \(({\mathcal {F}}_t)\) such that u satisfies (2.6) and

$$\begin{aligned} {\mathbb {E}}\left[ \sup _{t\in I} \left\| u(t)\right\| _{L^2_x}^q + \left( \int _I {\mathcal {J}}\big (u(t) \big ) \,\textrm{d}t \right) ^{q/2} \right] \lesssim \left\| u_0\right\| _{L^q_\omega L^2_x}^q + 1. \end{aligned}$$
(3.1)

The abstract conditions postulated by Liu and Röckner, that guarantee the existence of a weak solution, can be verified along the lines of [45, Example 4.1.9]. In fact, [45, Example 4.1.9] only deals with the case \(p \ge 2\). However, the general case \(p \in (1,\infty )\) can be recovered, if we adjust the underlying Gelfand triple to

$$\begin{aligned} W^{1,p}_{0,x} \cap L^2_{\text {div}} \hookrightarrow L^2_{\text {div}} \hookrightarrow \big ( W^{1,p}_{0,x} \cap L^2_{\text {div}} \big )'. \end{aligned}$$

The key tools for the a priori bound (3.1) are Itô’s formula for \(F(u) = \left\| u\right\| _{L^2_x}^q\) 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 [60, Theorem 4.3.1] for more details.

Remark 14

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 \in [1,\infty )\)

$$\begin{aligned} u_0 \in L^q_\omega L^2_x \quad \Rightarrow \quad u \in L^q_\omega C^0_t L^2_x. \end{aligned}$$

The limit case \(q = \infty \) needs to be excluded since the moment transfer is hindered by the Wiener process W. The most one can hope for is \(u \in L^{\Phi _2}_\omega L^2_x\) since \(W \in L_\omega ^{\Phi _2} \backslash L^\infty _\omega \). More details about integrability of Brownian motions in Banach spaces can be found in [36].

3.2 Reconstruction of the Pressure

Thanks to the Helmholtz-Leray projection \(\Pi _{\text {div}}\) it is possible to target the construction of the velocity u and the pressure \(\pi \) individually. The pressure corresponds to the residual error when testing (2.6) by smooth gradients rather than smooth divergence free fields, i.e., we use (2.7) as a definition of the pressure \(\pi \). Additionally, we artificially decompose \(\pi = \pi _{\textrm{det}} + \pi _{\textrm{sto}} \) into the components

$$\begin{aligned} \nabla \pi _{\textrm{det}}&= \Pi _{\textrm{div}}^{\perp } \text {div}S(\mathbb {\epsilon } u), \end{aligned}$$
(3.2a)
$$\begin{aligned} \nabla \pi _{\textrm{sto}}&= \Pi _{\textrm{div}}^{\perp } G(u) {\dot{W}}, \end{aligned}$$
(3.2b)

and discuss them separately.

First, we have a look at the deterministic component.

Lemma 15

Let \({\mathcal {O}}\) be a bounded \(C^2\) domain and \( u \in W^{1,p}_{0,x}\). Then there exists \(\pi _{\textrm{det}} \in L^{p'}_{0,x}\) such that for all \(\xi \in C^\infty _{0,x}\)

$$\begin{aligned} \int _{{\mathcal {O}}} \pi _\textrm{det} \text {div}\xi \,\textrm{d}x = \int _{{\mathcal {O}}} S(\mathbb {\epsilon } u): \nabla \Pi _{\text {div}}^\perp \xi \,\textrm{d}x \end{aligned}$$
(3.3)

and

$$\begin{aligned} \left\| \pi _{\textrm{det}}\right\| _{L^{p'}_x} \lesssim \left\| S(\mathbb {\epsilon } u)\right\| _{L^{p'}_x}. \end{aligned}$$
(3.4)

Proof

We interpret (3.2a) distributionally, i.e., for \(\xi \in C^\infty _{0,x}\)

$$\begin{aligned} \langle \pi _{\textrm{det}}, \text {div}\xi \rangle := \int _{{\mathcal {O}}} S(\mathbb {\epsilon } u): \nabla \Pi _{\text {div}}^\perp \xi \,\textrm{d}x . \end{aligned}$$
(3.5)

Without loss of generality we assume \(\langle \pi _{\textrm{det}}\rangle _{{\mathcal {O}}} = 0\).

Since we want to identify \(\pi _{\textrm{det}}\) as a proper function, we substitute \(\xi = {\mathcal {B}} \xi _0\), where \({\mathcal {B}}\) denotes the Bogovskii operator (see Theorem 34) and \(\xi _0 \in C^\infty _{0,x}\) with \(\langle \xi _0\rangle _{{\mathcal {O}}} = 0\), to obtain

$$\begin{aligned} \langle \pi _{\textrm{det}}, \xi _0 \rangle = \int _{{\mathcal {O}}} S(\mathbb {\epsilon } u): \nabla \Pi _{\text {div}}^\perp {\mathcal {B}} \xi _0 \,\textrm{d}x . \end{aligned}$$
(3.6)

Thus, using the symmetry of \(\Pi _{\text {div}}^\perp \) and integration by parts, (3.6) is equivalent to

$$\begin{aligned} \pi _{\textrm{det}} = -{\mathcal {B}}^* \Pi _{\text {div}}^\perp \text {div}S(\mathbb {\epsilon } u), \end{aligned}$$
(3.7)

where \({\mathcal {B}}^*\) is the adjoint operator of \({\mathcal {B}}\).

It remains to verify the boundedness of

$$\begin{aligned} A:= -{\mathcal {B}}^* \Pi _{\text {div}}^\perp \text {div}: L^{p'}_{0,x} \rightarrow L^{p'}_{0,x} . \end{aligned}$$

We instead address the boundedness of the adjoint operator and use a duality argument.

The right hand side of (3.6) is estimated by Hölder’s inequality

$$\begin{aligned} \int _{{\mathcal {O}}} S(\mathbb {\epsilon } u): \nabla \Pi _{\text {div}}^\perp {\mathcal {B}} \xi _0 \,\textrm{d}x \le \left\| S(\mathbb {\epsilon } u)\right\| _{L_x^{p'}} \,{\big \Vert { \nabla \Pi _{\text {div}}^\perp {\mathcal {B}} \xi _0 }\big \Vert }_{L^p_x}. \end{aligned}$$

We further estimate, due to the Sobolev stability of the Helmholtz-Leray projection (Theorem 32) and the Bogovskii operator (Theorem 34 with \(s = 0\), \(q = p\)),

$$\begin{aligned} {\big \Vert { \nabla \Pi _{\text {div}}^\perp {\mathcal {B}} \xi _0 }\big \Vert }_{L^p_x} \lesssim {\big \Vert { \nabla {\mathcal {B}} \xi _0 }\big \Vert }_{L^p_x} \lesssim {\big \Vert {\xi _0 }\big \Vert }_{L^p_x}. \end{aligned}$$

Therefore, the linear operator \(A: L^p_{0,x} \rightarrow L^p_{0,x}\) is bounded.

The claim (3.3) follows by our construction. The inequality (3.4) is verified by a density argument

$$\begin{aligned} \left\| \pi _{\textrm{det}}\right\| _{L^{p'}_{x}} = \sup _{ \xi \in L^p_{0,x}} \frac{\langle \pi _{\textrm{det}}, \xi \rangle }{\left\| \xi \right\| _{L^p_x}} = \sup _{ \xi \in C^\infty _{0,x}, \langle \xi \rangle _{{\mathcal {O}}} = 0 } \frac{\langle \pi _{\textrm{det}}, \xi \rangle }{\left\| \xi \right\| _{L^p_x}} \lesssim \left\| S(\mathbb {\epsilon } u)\right\| _{L_x^{p'}}. \end{aligned}$$

\(\square \)

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 16

Let \({\mathcal {O}}\) be a bounded domain with locally Lipschitz boundary and \({\mathcal {I}}_W\big (G(u) \big ) \in L^2_\omega B^{1/2}_{\Phi _2,\infty } L^2_x\). Then there exists \({\mathcal {K}}_{\textrm{sto}} \in L^2_\omega B^{1/2}_{\Phi _2,\infty } \big ( W^{1,2}_{0,x} \cap L^2_{0,x}\big ) \) such that for all \(t \in I\) and \(\xi \in C^\infty _{0,x}\) it holds

$$\begin{aligned} \int _{{\mathcal {O}}} {\mathcal {K}}_{\textrm{sto}}(t) \text {div}\xi \,\textrm{d}x&= - \int _{{\mathcal {O}}} {\mathcal {I}}_W \big ( G(u) \big )(t) \Pi _{\text {div}}^\perp \xi \,\textrm{d}x. \end{aligned}$$
(3.8)

Moreover,

$$\begin{aligned} \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{L^2_\omega B^{1/2}_{\Phi _2,\infty } W^{1,2}_x} \lesssim \left\| {\mathcal {I}}_W \big ( G(u) \big )\right\| _{L^2_\omega B^{1/2}_{\Phi _2,\infty } L^2_x}. \end{aligned}$$
(3.9)

Proof

The proof proceeds similar to the one of Lemma 15.

First, note (3.8) is equivalent to

$$\begin{aligned} {\mathcal {K}}_{\textrm{sto}}(t) := -{\mathcal {B}}^* \Pi _{\text {div}}^\perp {\mathcal {I}}_W\big ( G(u) \big )(t), \end{aligned}$$
(3.10)

where \({\mathcal {B}}^*\) is the adjoint of the Bogovskii operator. In particular, \({\mathcal {K}}_{\textrm{sto}}\) is mean-value free.

Due to Theorem 34, one can check that \( {\mathcal {B}}^*: L^2_{0,x} \rightarrow W^{1,2}_{0,x} \cap L^2_{0,x}\) is bounded. Indeed, let \(\xi , \zeta \in C^\infty _{0,x}\) such that \(\langle \xi \rangle _{{\mathcal {O}}} = \langle \zeta \rangle _{{\mathcal {O}}} = 0\), then by Hölder’s inequality and Theorem 34 with \(s = -1\) and \(q = 2\),

$$\begin{aligned} \int _{{\mathcal {O}}} {\mathcal {B}}^* \xi \zeta \,\textrm{d}x = \int _{{\mathcal {O}}} \xi \cdot {\mathcal {B}} \zeta \,\textrm{d}x \le \left\| \xi \right\| _{L^2_x} \left\| {\mathcal {B}} \zeta \right\| _{L^2_x} \lesssim \left\| \xi \right\| _{L^2_x} \left\| \zeta \right\| _{W^{-1,2}_{x}}. \end{aligned}$$

Thus, using the density of smooth mean-value free functions within \(W^{1,2}_{0,x} \cap L^2_{0,x}\),

$$\begin{aligned} \left\| {\mathcal {B}}^* \xi \right\| _{W^{1,2}_{x}} \lesssim \left\| \xi \right\| _{L^2_x}. \end{aligned}$$
(3.11)

Finally, (3.11) together with \(L^2_x\)-stability of the Helmholtz-Leray projection establish

$$\begin{aligned} \left\| {\mathcal {K}}_{\textrm{sto}}(t)\right\| _{W^{1,2}_x} \lesssim \left\| \Pi _{\text {div}}^\perp {\mathcal {I}}_W \big ( G(u) \big )(t) \right\| _{L^2_x} \le \left\| {\mathcal {I}}_W \big ( G(u) \big )(t) \right\| _{L^2_x}. \end{aligned}$$

An application of the \(L^2_\omega B^{1/2}_{\Phi _2,\infty }\)-norm verifies (3.9). \(\square \)

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 \(\xi \in C^\infty _0\big ( [0,T); C^\infty _{0,x} \big )\) we define the distribution

$$\begin{aligned} \langle \pi _{\textrm{sto}}, \xi \rangle := -\int _I \langle {\mathcal {K}}_{\textrm{sto}}, \partial _t \xi \rangle \,\textrm{d}t. \end{aligned}$$
(3.12)

To shorten the notation we abbreviate \(X = W^{1,2}_{0,x} \cap L^2_{0,x}\).

Lemma 17

Let \(\alpha \in (0,1)\) and \(p,q \in (1,\infty )\). Moreover, assume that \({\mathcal {K}}_{\textrm{sto}} \in L^2_\omega B^{\alpha }_{p,q} X\). Then \(\pi _{\textrm{sto}} \in L^2_\omega B^{\alpha -1}_{p,q} X\) and

$$\begin{aligned} \left\| \pi _{\textrm{sto}}\right\| _{L^2_\omega B^{\alpha -1}_{p,q} X} \lesssim \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{L^2_\omega B^{\alpha }_{p,q} X}. \end{aligned}$$
(3.13)

Proof

Let \(\xi \in C^\infty _0\big ( I; C^\infty _{0,x} \big )\) be spatially mean-value free. Due to (3.12) and duality

$$\begin{aligned} \langle \pi _{\textrm{sto}},\xi \rangle \le \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{B^{\alpha }_{p,q} X} \left\| \partial _t \xi \right\| _{(B^{\alpha }_{p,q} X)'}. \end{aligned}$$
(3.14)

Note that \(\xi \) 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 27,

$$\begin{aligned} \left\| \partial _t \xi \right\| _{\big (B^{\alpha }_{p,q}({\mathbb {R}}; X) \big )'} \eqsim \left\| \partial _t \xi \right\| _{\big ({\tilde{B}}^{\alpha }_{p,q}({\mathbb {R}}; X)\big )'}. \end{aligned}$$

Since X is reflexive and \(p, q < \infty \) we can identify the dual using Lemma 28

$$\begin{aligned} \left\| \partial _t \xi \right\| _{\big ({\tilde{B}}^{\alpha }_{p,q}({\mathbb {R}}; X \big )'} \eqsim \left\| \partial _t \xi \right\| _{{\tilde{B}}^{-\alpha }_{p',q'}({\mathbb {R}}; X')}. \end{aligned}$$

An application of Lemma 29 yields

$$\begin{aligned} \left\| \partial _t \xi \right\| _{{\tilde{B}}^{-\alpha }_{p',q'}({\mathbb {R}}; X')} \lesssim \left\| \xi \right\| _{{\tilde{B}}^{1-\alpha }_{p',q'}({\mathbb {R}}; X')}. \end{aligned}$$
(3.15)

Going back to the Besov norm in terms of integrability of differences, cf. Theorem 26, and using that \(\xi \) is extended by zero

$$\begin{aligned} \left\| \xi \right\| _{{\tilde{B}}^{1-\alpha }_{p',q'}({\mathbb {R}}; X')} \eqsim \left\| \xi \right\| _{B^{1-\alpha }_{p',q'}({\mathbb {R}}; X')} = \left\| \xi \right\| _{B^{1-\alpha }_{p',q'} X'}. \end{aligned}$$

This concludes

$$\begin{aligned} \left\| \partial _t \xi \right\| _{\big (B^{\alpha }_{p,q}({\mathbb {R}}; X) \big )'} \lesssim \left\| \xi \right\| _{B^{1-\alpha }_{p',q'} X'}. \end{aligned}$$
(3.16)

Overall, revisiting (3.14) and using the density of smooth mean-value free functions in X together with (3.16),

$$\begin{aligned} \left\| \pi _{\textrm{sto}}\right\| _{(B^{1-\alpha }_{p',q'} X')'} = \sup _{\xi \in C^\infty _0( I; C^\infty _{0,x} ),\, \forall t:\langle \xi (t)\rangle = 0} \frac{\langle \pi _{\textrm{sto}}, \xi \rangle }{\left\| \xi \right\| _{B^{1-\alpha }_{p',q'}X}} \lesssim \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{B^{\alpha }_{p,q} X}. \end{aligned}$$

Ultimately, due to the equivalent norm on the dual space, cf. Lemma 28,

$$\begin{aligned} \left\| \pi _{\textrm{sto}}\right\| _{(B^{1-\alpha }_{p',q'} X')'} \eqsim \left\| \pi _{\textrm{sto}}\right\| _{B^{\alpha -1}_{p,q} X}. \end{aligned}$$

Taking the square and expectation verifies (3.13). \(\square \)

Remark 18

So far the stochastic pressure was always estimated on suboptimal spaces. For example in [18, 41] the authors infer \(\pi _{\textrm{sto}} \in L^1_\omega W^{-1,\infty }_t L^2_{0,x}\). They neglect all additional information on the temporal differentiability of \({\mathcal {K}}_{\textrm{sto}}\) and only use the boundedness.

Whether a corresponding result of Lemma 17 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.15), i.e.,

$$\begin{aligned} \left\| \partial _t \xi \right\| _{B^{-\alpha }_{\Phi _2^*,1} W^{-1,2}_x} \lesssim \left\| \xi \right\| _{B^{1-\alpha }_{\Phi _2^*,1} W^{-1,2}_x}. \end{aligned}$$
(3.17)

In order to establish (3.15), 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 \in \{ 1, \infty \}\). Particularly, \(\Phi _2(t) = e^{t^2}-1\) is to close to \(\infty \) and whence \(\Phi _2^*\) is to close to 1.

Still, in regard of the regularity result of the integrated pressure, cf. Lemma 16, it is natural to define the set

$$\begin{aligned} {\mathbb {B}}^{-\alpha }_{\Phi ,r}(0,T; X)&:= \left\{ u \in \big ( C^\infty _0( (0,T);X') \big )' \big | \, (3.18\hbox {b}) \right\} , \end{aligned}$$
(3.18a)
$$\begin{aligned} \exists g \in B^{1-\alpha }_{\Phi ,r}(0,T;X) \,&\forall \xi \in C^\infty _0((0,T);X') :\, \langle u, \xi \rangle = -\int _0^T \langle g, \partial _t \xi \rangle \,\textrm{d}t. \end{aligned}$$
(3.18b)

Corollary 19

In the setting of Lemma 16 we have \(\pi _{\textrm{sto}} \in {\mathbb {B}}^{-1/2}_{\Phi _2,\infty }(0,T; X)\) \({\mathbb {P}}\)-a.s.

Proof

We only need to verify that the distribution \(\pi _{\textrm{sto}}\) defined by (3.12) is well-defined. The claim that \(\pi _{\textrm{sto}}\) belongs to \( {\mathbb {B}}^{-1/2}_{\Phi _2,\infty }(0,T; X)\) follows by the definition (3.18) and the assumption on \({\mathcal {K}}_{\textrm{sto}}\).

Let \(\xi \in C^\infty _0\big ((0,T); X' \big )\). Due to duality and Hölder’s inequality

$$\begin{aligned} \int _0^T \langle {\mathcal {K}}_{\textrm{sto}}, \partial _t \xi \rangle \,\textrm{d}t&\le \int _0^T \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{X}\left\| \partial _t \xi \right\| _{X'} \,\textrm{d}t \\&\lesssim \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{L^{\Phi _2}_t X} \left\| \partial _t \xi \right\| _{L^{\Phi _2^*}_t X'}. \end{aligned}$$

Next, we take the square and expectation

$$\begin{aligned} {\mathbb {E}} \left| \langle \pi _{\textrm{sto}}, \xi \rangle \right| ^2 \lesssim \left\| {\mathcal {K}}_{\textrm{sto}}\right\| _{L^2_\omega L^{\Phi _2}_t X}^2 \left\| \partial _t \xi \right\| _{L^{\Phi _2^*}_t X'}^2. \end{aligned}$$

Thus, we have shown that \(\pi _{\textrm{sto}}\) is well-defined. \(\square \)

3.3 Proof of Theorem 5

The previous results lead to a simple proof of Theorem 5.

Proof of Theorem 5

Due to the factorization (1.3) it is possible to construct the velocity field u independently of the pressure \(\pi \). Afterwards the pressure is reconstructed.

The existence of the velocity field and verification of the a priori bound (2.10a) is done in Theorem 13.

Section 3.2 deals with the pressure reconstruction; Lemma 15 verifies the regularity of the deterministic pressure; the stochastic pressure is handled in Lemma 16 and 17. In particular, Lemma 16 ensures that \({\mathcal {K}}_{\textrm{sto}} \in L^2_\omega B^{1/2}_{\Phi _2,\infty } X \hookrightarrow L^2_\omega B^{\alpha }_{s,r} X\) for any \(\alpha \in (0,1/2)\) and \(s,r \in (1,\infty )\). Therefore Lemma 17 implies (2.10c). \(\square \)

4 Temporal Regularity 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.

4.1 Existence of Strong Solutions

Strong solutions for a related model have already been constructed in [13] 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 (A. 19a), which – at least formally – corresponds to a test with \(-\Pi _{\text {div}}\text {div}S(\mathbb {\epsilon } u)\).

An application of Itô’s formula to \({\mathcal {J}}(u)^{q}\), \(q\in [1,\infty )\), provides

$$\begin{aligned} \,\textrm{d}{\mathcal {J}}(u)^q&= q {\mathcal {J}}(u)^{q-1} D {\mathcal {J}}(u)[\,\textrm{d}u] \\&+ \frac{q(q-1)}{2} {\mathcal {J}}(u)^{q-2} \big ( D {\mathcal {J}}(u) [ \,\textrm{d}u] \big )^2 + \frac{q}{2} {\mathcal {J}}(u)^{q-1} D^2 {\mathcal {J}}(u)[\,\textrm{d}u, \,\textrm{d}u], \end{aligned}$$

where \(D{\mathcal {J}}(u)\) and \(D^2{\mathcal {J}}(u)\) are given by (A. 23) and (A. 24), respectively. Most importantly, the deterministic integral in the first term has a sign, i.e.

$$\begin{aligned} D {\mathcal {J}}(u)[\,\textrm{d}u]&= - \int _{{\mathcal {O}}} \left| \Pi _{\text {div}} \text {div}S(\epsilon u)\right| ^2 \,\textrm{d}x \,\textrm{d}t \\&\hspace{2em} + \sum _{j \in {\mathbb {N}}} \int _{{\mathcal {O}}} (\kappa + \left| \epsilon u\right| )^{p-2} \epsilon u : \epsilon \Pi _{\text {div}} g_j(\cdot , u) \,\textrm{d}x \,\textrm{d}\beta ^j. \end{aligned}$$

The Itô corrections need to be interpreted as

$$\begin{aligned} \big ( D {\mathcal {J}}(u) [ \,\textrm{d}u] \big )^2&= \sum _{j\in {\mathbb {N}}} \left( \int _{{\mathcal {O}}} (\kappa + \left| \epsilon u\right| )^{p-2} \epsilon u : \epsilon \Pi _{\text {div}} g_j(\cdot , u) \,\textrm{d}x \right) ^2 \,\textrm{d}t, \end{aligned}$$

and, recalling \(\varphi _{\kappa }(t) = \int _0^t (\kappa + s)^{p-2} s \,\textrm{d}s\),

$$\begin{aligned} D^2 {\mathcal {J}}(u)[\,\textrm{d}u, \,\textrm{d}u]&= \sum _{j\in {\mathbb {N}}} \int _{{\mathcal {O}}} \left( \varphi _{\kappa }''\big (\left| \mathbb {\epsilon } u\right| ) - \frac{\varphi _{\kappa }'\big ( \left| \mathbb {\epsilon } u\right| \big )}{\left| \mathbb {\epsilon } u\right| } \right) \left( \frac{\mathbb {\epsilon } u}{\left| \mathbb {\epsilon } u\right| }: \mathbb {\epsilon }\Pi _{\text {div}} g_j(\cdot , u) \right) ^2 \,\textrm{d}x \,\textrm{d}t \\&\quad + \sum _{j\in {\mathbb {N}}} \int _{{\mathcal {O}}} \frac{\varphi _{\kappa }'\big ( \left| \mathbb {\epsilon } u\right| \big ) }{\left| \mathbb {\epsilon } u\right| } \left| \mathbb {\epsilon }\Pi _{\text {div}} g_j(\cdot , u)\right| ^2 \,\textrm{d}x \,\textrm{d}t. \end{aligned}$$

A priori bounds for the energy can be established whenever the correction terms can be bounded with respect to the energy \({\mathcal {J}}(u)\).

The above computation have been made rigorous in the abstract theory developed by Gess in [35].

Proof of Theorem 9

We will only address the case \(q = 1\). More details on the general case can be found in [60, Theorem 4.3.1]. Our aim is to apply the result [35, Theorem 1.4], i.e., we need to check the assumptions (A1\()-(\)A6).

First of all, we relate our framework to the one used in [35]. We choose \(H = L^2_{\text {div}}\), \(S = W^{1,p}_{0,x}\) and \(V = H \cap S\). Additionally, we let \(\varphi (v) = {\tilde{\varphi }}(v) = {\mathcal {J}}(v)\) and \(B_t(v) \equiv \Pi _{\text {div}} G(v)\). Keep in mind that we use \(\varphi \) as the potential that defines the energy whereas Gess uses \(\varphi \) to denote the energy itself.

Ad (A1): First we show that \({\mathcal {J}}: W^{1,p}_{0,x} \rightarrow {\mathbb {R}}\) is continuous. Indeed, let \(u,v \in W^{1,p}_{0,x}\). Then, due to the fundamental theorem,

$$\begin{aligned} \left| {\mathcal {J}}(u) - {\mathcal {J}}(v)\right|&= \left| \int _{{\mathcal {O}}} \varphi _{\kappa }(\left| \mathbb {\epsilon } u\right| ) - \varphi _{\kappa }(\left| \mathbb {\epsilon } u\right| ) \,\textrm{d}x\right| \\&= \left| \int _{{\mathcal {O}}} \int _0^1 \varphi _{\kappa }'( \left| \mathbb {\epsilon } w_\theta \right| ) \frac{\mathbb {\epsilon } w_\theta }{\left| \mathbb {\epsilon } w_\theta \right| }: \mathbb {\epsilon }(u-v) \,\textrm{d}\theta \,\textrm{d}x\right| \\&= \left| \int _{{\mathcal {O}}} \int _0^1 (\kappa + \left| \mathbb {\epsilon } w_\theta \right| )^{p-2} \mathbb {\epsilon } w_\theta : \mathbb {\epsilon }(u-v) \,\textrm{d}\theta \,\textrm{d}x\right| \\&\le \int _{{\mathcal {O}}} \int _0^1 (\kappa + \left| \mathbb {\epsilon } w_\theta \right| )^{p-1} \,\textrm{d}\theta \left| \mathbb {\epsilon }(u-v)\right| \,\textrm{d}x, \end{aligned}$$

where \(w_\theta = v + \theta (u-v)\). Note that \(\left| \mathbb {\epsilon } w_\theta \right| \le \max \{ \left| \mathbb {\epsilon } u\right| , \left| \mathbb {\epsilon } v\right| \}\) for all \(\theta \in [0,1]\). This and Hölder’s inequality imply

$$\begin{aligned}&\int _{{\mathcal {O}}} \int _0^1 (\kappa + \left| \mathbb {\epsilon } w_\theta \right| )^{p-1} \,\textrm{d}\theta \left| \mathbb {\epsilon }(u-v)\right| \,\textrm{d}x \\&\quad \le \left( \int _{{\mathcal {O}}} (\kappa + \max \{ \left| \mathbb {\epsilon } u\right| , \left| \mathbb {\epsilon } v\right| \})^{p} \,\textrm{d}x \right) ^{(p-1)/p} \left( \int _{{\mathcal {O}}} \left| \mathbb {\epsilon }(u-v)\right| ^p \,\textrm{d}x \right) ^{1/p}. \end{aligned}$$

Finally, we arrive at

$$\begin{aligned} \left| {\mathcal {J}}(u) - {\mathcal {J}}(v)\right| \le \left( \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^p \,\textrm{d}x + \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } v\right| )^p \,\textrm{d}x \right) ^{(p-1)/p} \left\| \nabla (u-v)\right\| _{L^p_x}. \end{aligned}$$

The continuity is verified.

Moreover, using that \(\varphi '(t) = t^{p-1}\) is homogeneous of degree \(p-1\) and a substitution,

$$\begin{aligned} \varphi _{\kappa }(2t)&= \int _0^{2t} \frac{\varphi '(\kappa + s)}{\kappa + s } s \,\textrm{d}s = \int _0^{2t} 2^{p-1} \frac{\varphi '(\frac{\kappa }{2} + \frac{s}{2})}{\frac{\kappa }{2} + \frac{s}{2} } \frac{s}{2} \,\textrm{d}s \\&=2^p \int _0^{t} \frac{\varphi '(\frac{\kappa }{2} + s)}{\frac{\kappa }{2} + s } s \,\textrm{d}s = 2^p \varphi _{\frac{\kappa }{2}}(t). \end{aligned}$$

Since \(p \ge 2\) we find \(\varphi _{\kappa }(t)\) is increasing in \(\kappa \) for fixed t and whence \({\mathcal {J}}(2u) \le 2^p {\mathcal {J}}(u)\).

Ad (A2): Recall

$$\begin{aligned} \varphi _{\kappa }(t)&= \frac{1}{p}(\kappa + t)^p - \frac{\kappa }{p-1}(\kappa + t)^{p-1} + \frac{\kappa ^p}{p(p-1)}, \\ \varphi _{\kappa }'(t)&=(\kappa + t)^{p-2} t, \\ \varphi _{\kappa }''(t)&= (\kappa + t)^{p-2}\left( 1 + (p-2) \frac{t}{\kappa + t}\right) . \end{aligned}$$

Applying a weighted Young’s inequality one finds

$$\begin{aligned} \frac{1}{2p}(\kappa + t)^p - \frac{2^{p-1} - 1}{p(p-1)}\kappa ^p \le \varphi _{\kappa }(t) \le \frac{1}{p}(\kappa + t)^p + \frac{1}{p(p-1)}\kappa ^p. \end{aligned}$$
(4.1)

In other words \(\varphi _{\kappa }\) behaves like the pth-power shifted along \(\kappa \). It follows using the monotonicity of powers and (4.1)

$$\begin{aligned}&\int _{{\mathcal {O}}} \frac{1}{p}\left| \mathbb {\epsilon } u\right| ^p \,\textrm{d}x \le \int _{{\mathcal {O}}}\frac{1}{p} ( \kappa + \left| \mathbb {\epsilon } u\right| )^p \,\textrm{d}x \lesssim \int _{{\mathcal {O}}} \varphi _{\kappa }(\left| \mathbb {\epsilon } u\right| ) + \kappa ^p \left| {\mathcal {O}}\right| \eqsim {\mathcal {J}}(u) + 1. \end{aligned}$$
(4.2)

Let \({\{{w_k}\}}_{k \in {\mathbb {N}}} \subset W^{1,p}_{0,x}\). Now, using (A. 24), Hölder’s and Young’s inequalities and (4.2)

$$\begin{aligned} \sum _{k=1}^\infty D^2 {\mathcal {J}}(v)[w_k,w_k]&\le (p-1) \sum _{k=1}^\infty \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } v\right| )^{p-2} \left| \mathbb {\epsilon } w_k\right| ^2 \,\textrm{d}x \\&\le (p-1) \left( \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } v\right| )^{p} \,\textrm{d}x \right) ^{(p-2)/p} \sum _{k=1}^\infty \left\| \mathbb {\epsilon } w_k\right\| _{L^p_x}^2 \\&\le \frac{(p-1)p}{p-2} \int (\kappa + \left| \mathbb {\epsilon } v\right| )^p \,\textrm{d}x + \frac{(p-1)2}{p} \left( \sum _{k=1}^\infty \left\| \mathbb {\epsilon } w_k\right\| _{L^p_x}^2 \right) ^{p/2} \\&\lesssim {\mathcal {J}}(v) + 1 + \left( \sum _{k=1}^\infty \left\| \mathbb {\epsilon } w_k\right\| _{L^p_x}^2 \right) ^{p/2}. \end{aligned}$$

Ad (A3): Lemma 42 shows that \({\mathcal {J}}\) is strongly convex.

Ad (A4): Recall (A. 23). Therefore,

$$\begin{aligned} -D{\mathcal {J}}(v)[v] = -\int _{{\mathcal {O}}} \varphi _{\kappa }'(\left| \mathbb {\epsilon } v\right| ) \left| \mathbb {\epsilon } v\right| \,\textrm{d}x \eqsim - \int _{{\mathcal {O}}} \varphi _{\kappa }(\left| \mathbb {\epsilon } v\right| ) \,\textrm{d}x = - {\mathcal {J}}(v) \le 0. \end{aligned}$$

Ad (A5): Follows by the assumption on the noise coefficient, cf. (2.9).

Ad (A6): The major ingredient that enables strong solutions is the regularity of the noise. Recall (2.3) and (A2) with \(w_k = \Pi _{\text {div}}G(v)u_k = \Pi _{\text {div}} g_k(\cdot , v )\). It remains to estimate

$$\begin{aligned} \left( \sum _{k=1}^\infty \left\| \mathbb {\epsilon } w_k\right\| _{L^p_x}^2 \right) ^{p/2} = \left( \sum _{k=1}^\infty \left\| \mathbb {\epsilon }\Pi _{\text {div}} g_k(\cdot , v )\right\| _{L^p_x}^2 \right) ^{p/2} =\left\| \mathbb {\epsilon } \Pi _{\text {div}} G(v)\right\| _{L_2(U;L^p_x)}^p. \end{aligned}$$

An application of (2.12) and (4.2) show

$$\begin{aligned} \left\| \mathbb {\epsilon } \Pi _{\text {div}} G(v)\right\| _{L_2(U;L^p_x)}^p \lesssim \left\| \mathbb {\epsilon } v\right\| _{L^p_x}^p + 1 \lesssim {\mathcal {J}}(v) + 1. \end{aligned}$$

Overall, we have verified (A1\()-(\)A6) and therefore we may apply [35, Theorem 1.4] which provides us with the a priori bound

$$\begin{aligned} \sup _{t \in I} {\mathbb {E}} \left[ {\mathcal {J}}(u) \right] + {\mathbb {E}} \left[ \int _I \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 \,\textrm{d}s \right] \lesssim {\mathbb {E}} \left[ {\mathcal {J}}(u_0) \right] + 1. \end{aligned}$$

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.,

$$\begin{aligned}&{\mathbb {E}} \left[ \sup _{t\in I} \int _0^t \sum _{k} D {\mathcal {J}}(u)[\Pi _{\text {div}} g_k(\cdot , v)] \,\textrm{d}\beta (s) \right] \\&\quad = {\mathbb {E}} \left[ \sup _{t\in I} \int _0^t \sum _{k} \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^{p-2} \mathbb {\epsilon } u: \mathbb {\epsilon } \Pi _{\text {div}} g_k(\cdot , v)] \,\textrm{d}\beta ^k(s) \right] . \end{aligned}$$

Due to the Burkholder-Davis-Gundy inequality

$$\begin{aligned}&{\mathbb {E}} \left[ \sup _{t\in I} \int _0^t \sum _{k} \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^{p-2} \mathbb {\epsilon } u: \mathbb {\epsilon } \Pi _{\text {div}} g_k(\cdot , v)] \,\textrm{d}x \,\textrm{d}\beta ^k(s) \right] \\&\quad \eqsim {\mathbb {E}} \left[ \left( \int _I \sum _{k} \left( \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^{p-2} \mathbb {\epsilon } u: \mathbb {\epsilon } \Pi _{\text {div}} g_k(\cdot , v)] \,\textrm{d}x \right) ^2 \,\textrm{d}s \right) ^{1/2} \right] . \end{aligned}$$

Invoking Hölder’s inequality

$$\begin{aligned}&{\mathbb {E}} \left[ \left( \int _I \sum _{k} \left( \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^{p-2} \mathbb {\epsilon } u: \mathbb {\epsilon } \Pi _{\text {div}} g_k(\cdot , v)] \,\textrm{d}x \right) ^2 \,\textrm{d}s \right) ^{1/2} \right] \\&\quad \le {\mathbb {E}} \left[ \left( \int _I \left( \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^{p} \,\textrm{d}x \right) ^{2(p-1)/p} \left\| \mathbb {\epsilon } \Pi _{\text {div}} G(v)\right\| _{L_2(U;L^p_x)}^2 \,\textrm{d}s \right) ^{1/2} \right] . \end{aligned}$$

The Assumption 7 together with (4.2) show

$$\begin{aligned}&{\mathbb {E}} \left[ \left( \int _I \left( \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| )^{p} \,\textrm{d}x \right) ^{2(p-1)/p} \left\| \mathbb {\epsilon } \Pi _{\text {div}} G(v)\right\| _{L_2(U;L^p_x)}^2 \,\textrm{d}s \right) ^{1/2} \right] \\&\quad \lesssim {\mathbb {E}} \left[ \left( \int _I {\mathcal {J}}(u)^2+1 \,\textrm{d}s \right) ^{1/2} \right] \lesssim {\mathbb {E}} \left[ \left( \int _I {\mathcal {J}}(u)^2 \,\textrm{d}s \right) ^{1/2} \right] + 1. \end{aligned}$$

Finally, Young’s inequality verifies

$$\begin{aligned}&{\mathbb {E}} \left[ \sup _{t\in I} \int _0^t \sum _{k} D {\mathcal {J}}(u)[\Pi _{\text {div}} g_k(\cdot , v)] \,\textrm{d}\beta (s) \right] \\&\quad \le \delta {\mathbb {E}} \left[ \sup _{t \in I} {\mathcal {J}}(u) \right] + c_\delta \left( {\mathbb {E}} \left[ \int _I {\mathcal {J}}(u) \,\textrm{d}s \right] + 1 \right) . \end{aligned}$$

This allows for a Gronwall type argument that establishes the improved a priori estimate

$$\begin{aligned} {\mathbb {E}} \left[ \sup _{t \in I} {\mathcal {J}}(u) \right] +{\mathbb {E}} \left[ \int _I \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 \,\textrm{d}s \right] \lesssim {\mathbb {E}}\left[ {\mathcal {J}}(u_0) \right] + 1. \end{aligned}$$

\(\square \)

Remark 20

On a first view we could ask, if one could expand the q-energy (applying Itô’s formula to \(\left\| \epsilon u\right\| _{L^q_x}^q\)) 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

$$\begin{aligned} \int _{{\mathcal {O}}} \Pi _{\text {div}} \text {div}\big (\left| \mathbb {\epsilon } u\right| ^{p-2} \mathbb {\epsilon } u \big ) \cdot \Pi _{\text {div}} \text {div}\big ( \left| \mathbb {\epsilon } u\right| ^{q-2} \mathbb {\epsilon } u \big ) \,\textrm{d}x \ge 0, \end{aligned}$$

then it is possible to obtain improved gradient regularity.

Regularity questions about \(\text {div}S(\mathbb {\epsilon } u)\) – and ultimately on \(\nabla S(\mathbb {\epsilon } u)\) – are out of reach, since

$$\begin{aligned} \Pi _{\text {div}}^\perp \text {div}S(\mathbb {\epsilon } u) = \nabla \pi _{\textrm{det}} \in L^{p'}_\omega L^{p'}_t W^{-1,p'}_x \end{aligned}$$

is only a distribution. The derivation of higher regularity for the pressure is a non-trivial task.

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 [61].

4.2 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 [36]. They show that \({\mathbb {P}}\)-almost surely Wiener processes belong to the Besov space \(B^{1/2}_{\Phi _2,\infty }\), \(\Phi _2(t) = e^{t^2}-1\), and do not belong to the Besov space \(B^{1/2}_{p,q}\) for \(p,q < \infty \).

We show that strong solutions inherit similar temporal regularity properties from the driving Wiener process.

Proof of Theorem 10

We will only present part (b) of Theorem 10 in detail and comment on the necessary changes for part (a).

We start with the estimate for the Besov-Orlicz semi-norm

$$\begin{aligned} {\mathbb {E}} \left[ {[{u}]}_{B^{1/2}_{\Phi _2,\infty }(I;L^2_x)}^2 \right] = {\mathbb {E}} \left[ \left( \sup _h h^{-1/2} \left\| \tau _h u\right\| _{L^{\Phi _2}(I \cap I-{\{{h}\}}; L^2_x)} \right) ^2 \right] . \end{aligned}$$

Using the strong formulation (2.11)

$$\begin{aligned} {\mathbb {E}} \left[ {[{u}]}_{B^{1/2}_{\Phi _2,\infty }(I;L^2_x)}^2 \right]&\lesssim {\mathbb {E}} \left[ \left( \sup _h h^{-1/2} \left\| \int _{\cdot }^{\cdot + h} \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \,\textrm{d}s\right\| _{L^{\Phi _2}(I \cap I-{\{{h}\}}; L^2_x)} \right) ^2 \right] \\&\quad + {\mathbb {E}} \left[ \left( \sup _h h^{-1/2} \left\| \tau _h {\mathcal {I}}(G(u))\right\| _{L^{\Phi _2}(I \cap I-{\{{h}\}}; L^2_x)} \right) ^2 \right] \\&=: \textrm{I} + \textrm{II}. \end{aligned}$$

Invoking \(L^{\infty } \hookrightarrow L^{\Phi _2}\) and Hölder’s inequality

$$\begin{aligned} \textrm{I}&\lesssim {\mathbb {E}} \left[ \left( \sup _h h^{-1/2} \sup _{t \in I \cap I-{\{{h}\}}} \int _{t}^{t + h} \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \right\| _{L^2_x} \,\textrm{d}s\right) ^2 \right] \\&\le {\mathbb {E}} \left[ \sup _h \sup _{t \in I \cap I-{\{{h}\}}} \int _{t}^{t + h} \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \right\| _{L^2_x}^2 \,\textrm{d}s\right] \\&\le \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_\omega L^2_t L^2_x}^2 . \end{aligned}$$

The second term is bounded by Theorem 35 (e),  (A. 14) and (2.8)

$$\begin{aligned} \textrm{II} \le \left\| {\mathcal {I}}(G(u))\right\| _{L^{2}_\omega B_{\Phi _2,\infty }^{1/2}(0,T;L^2_x)}^2&\lesssim \left\| G(u)\right\| _{L^{N_2}_\omega L^\infty _t \gamma (U;L^2_x)}^2 \\&= \left\| G(u)\right\| _{L^{N_2}_\omega L^\infty _t L_2(U;L^2_x)}^2 \lesssim \left\| u\right\| _{L^{N_2}_\omega L^\infty _t L^2_x}^2 + 1, \end{aligned}$$

where \(N_2(t) = t^2 \ln (t+1)\). This establishes the semi-norm estimate.

The full norm follows by the embedding \(L^\infty \hookrightarrow L^{\Phi _2}\)

$$\begin{aligned} \left\| u\right\| _{L^2_\omega B^{1/2}_{\Phi _2,\infty } L^2_x}^2&\lesssim \left\| u\right\| _{L^2_\omega L^{\Phi _2}_t L^2_x}^2 + {\mathbb {E}} \left[ {[{u}]}_{B^{1/2}_{\Phi _2,\infty }(I;L^2_x)}^2 \right] \\&\lesssim \left\| u\right\| _{L^{N_2}_\omega L^\infty _t L^2_x}^2 + \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_\omega L^2_t L^2_x}^2 + 1. \end{aligned}$$

We have verified (2.15).

Part (a) follows similarly to part (b), but we use the embedding \(L^\infty _t \hookrightarrow L^q_t\) and Theorem 35 (b) together with an extrapolation result, cf. [53, Remark 3.4]. \(\square \)

Remark 21

Theorem 10 immediately implies Hölder regularity up to almost 1/2 by an application of a Sobolev embedding (cf. [58, Sect. 2.3.2]), i.e., for any \(\alpha \in (0,1/2)\)

$$\begin{aligned} u \in L^2_{\omega } B^{1/2}_{\Phi _2,\infty } L^2_x \hookrightarrow L^2_{\omega } C^{0,\alpha }_t L^2_x. \end{aligned}$$

4.3 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 22

Let \(p \ge 2\) and Assumption 7 be satisfied. Then there exists a constant \(C > 0\) such that for all progressively measurable \(u \in L^{p}_\omega L^\infty _t \big ( L^2_{\text {div}} \cap W^{1,p}_{0,x}\big )\) it holds

$$\begin{aligned} \left\| \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}\big (G(u)\big )\right\| _{L^p_\omega B^{1/2}_{p,\infty } L^p_x}^p \le C\big ( 1 + \left\| \mathbb {\epsilon } u\right\| _{L^{p}_\omega L^\infty _t L^p_x}^p \big ). \end{aligned}$$
(4.3)

Proof

Due to the linearity of \(\mathbb {\epsilon } \Pi _{\text {div}}\) and \({\mathcal {I}}\) we find

$$\begin{aligned} \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}\big (G(u)\big ) = {\mathcal {I}}\big (\mathbb {\epsilon } \Pi _{\text {div}}G(u)\big ). \end{aligned}$$

Therefore, using the time regularity of stochastic integrals, cf. Theorem 35 (b),

$$\begin{aligned} \left\| \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}(G(u))\right\| _{L^p_\omega B^{1/2}_{p,\infty } L^p_x}^p&= \left\| {\mathcal {I}}( \mathbb {\epsilon } \Pi _{\text {div}} G(u))\right\| _{L^p_\omega B^{1/2}_{p,\infty } L^p_x}^p \\&\lesssim \left\| \mathbb {\epsilon } \Pi _{\text {div}} G(u)\right\| _{L^p_\omega L^\infty _t \gamma (U; L^p_x)}^p. \end{aligned}$$

Since \(L^p_x\) is of 2-smooth, we can use (A. 14). This and the growth assumption (2.12) establish

$$\begin{aligned} \left\| \mathbb {\epsilon } \Pi _{\text {div}} G(u)\right\| _{L^{p}_\omega L^{\infty }_t \gamma (U;L^p_x)}^p \lesssim \left\| \mathbb {\epsilon } u\right\| _{L^{N_p}_\omega L^{\infty }_t L^p_x}^p + 1. \end{aligned}$$

\(\square \)

Now we are ready to prove the result on temporal regularity of the non-linear gradient.

Proof of Theorem 12

The V-coercivity, cf. Lemma 38, allows to rewrite

$$\begin{aligned} \left\| \tau _h V(\mathbb {\epsilon } u)\right\| _{L^2_x}^2&\eqsim \left( \tau _h S(\mathbb {\epsilon } u), \tau _h \mathbb {\epsilon } u \right) . \end{aligned}$$
(4.4)

It follows, using the symmetry and integration by parts,

$$\begin{aligned} \left( \tau _h S(\mathbb {\epsilon } u), \tau _h \mathbb {\epsilon } u \right) = \left( \tau _h S(\mathbb {\epsilon } u), \tau _h \nabla u \right) = -\left( \tau _h \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u), \tau _h u \right) . \end{aligned}$$

Applying the strong formulation (2.11)

$$\begin{aligned}&\left( \tau _h \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u), \tau _h u \right) \\&\hspace{2em }= \left( \tau _h \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u), \int _{\cdot - h}^\cdot \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u) \,\textrm{d}\tau \right) \\&\hspace{4em} + \left( \tau _h \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u), \tau _h \Pi _{\text {div}} {\mathcal {I}}(G(u)) \right) \\&\hspace{2em } =: \textrm{I} + \textrm{II}. \end{aligned}$$

The first term is estimated by Hölder’s and Young’s inequalities

(4.5)

The stochastic term needs a more refined analysis. Due to the symmetry

$$\begin{aligned} \textrm{II}&= -\left( \tau _h S(\mathbb {\epsilon } u), \tau _h \nabla \Pi _{\text {div}} {\mathcal {I}}(G(u)) \right) = - \left( \tau _h S(\mathbb {\epsilon } u), \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}(G(u)) \right) . \end{aligned}$$

Using Lemma 38 and Young’s inequality

$$\begin{aligned} \textrm{II}&\le \int _{{\mathcal {O}}} \left| \tau _h S(\mathbb {\epsilon } u)\right| \left| \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}(G(u))\right| \,\textrm{d}x \\&\eqsim \int _{{\mathcal {O}}} \left( \kappa + \left| \mathbb {\epsilon } u\right| + \left| \tau _h \mathbb {\epsilon } u\right| \right) ^{p-2} \left| \tau _h \mathbb {\epsilon } u\right| \left| \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}(G(u))\right| \,\textrm{d}x \\&\le \delta \left\| \tau _h V(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 + c_\delta \int _{{\mathcal {O}}} \left( \kappa + \left| \mathbb {\epsilon } u\right| + \left| \tau _h \mathbb {\epsilon } u\right| \right) ^{p-2} \left| \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}} (G(u))\right| ^2 \,\textrm{d}x\\&=: \textrm{II}_1 + \textrm{II}_2. \end{aligned}$$

The first term can be absorbed. The remaining term is estimated by Hölder’s and Young’s inequalities

$$\begin{aligned} \begin{aligned} \textrm{II}_2&\le h c_\delta \left( \frac{p-2}{p} \int _{{\mathcal {O}}} (\kappa + \left| \mathbb {\epsilon } u\right| + \left| \tau _h \mathbb {\epsilon } u\right| )^p \,\textrm{d}x+ \frac{2}{p} h^{-p/2} \left\| \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}} (G(u))\right\| _{L^p_x}^p \right) \\&\lesssim h c_\delta \left( \frac{p-2}{p} \left( \sup _{t \in I} \left\| \mathbb {\epsilon } u\right\| _{L^p_x}^p + 1 \right) + \frac{2}{p} h^{-p/2} \left\| \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}} (G(u))\right\| _{L^p_x}^p \right) . \end{aligned} \end{aligned}$$
(4.6)

We are in position to estimate the non-linear expression \(V(\mathbb {\epsilon } u)\). Integrate (4.4) in time and rescale by \(h^{-1}\). Moreover, take the supremum over h and expectation. Finally, apply (4.5) and (4.6) to find

Fubini’s theorem and the a priori estimate (2.13) imply

The time regularity of the stochastic integral, cf. Lemma 22, allows to conclude

$$\begin{aligned}&{\mathbb {E}} \left[ \sup _{h\in I} h^{-p/2} \int _{I_h} \left\| \tau _h \mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}\big ( G(u)\big )\right\| _{L^p_x}^p \,\textrm{d}s \right] \\&\hspace{2em} = {\mathbb {E}} \left[ {[{\mathbb {\epsilon } \Pi _{\text {div}} {\mathcal {I}}\big ( G(u) \big )}]}_{B^{1/2}_{2,\infty } L^p_x}^p \right] \lesssim \left\| \mathbb {\epsilon } u\right\| _{L^{p}_\omega L^\infty _t L^p_x}^p + 1. \end{aligned}$$

Choosing \(\delta > 0\) sufficiently small and (4.2) lead to

$$\begin{aligned} {\mathbb {E}}\left[ {[{V(\mathbb {\epsilon } u)}]}_{ B^{1/2}_{2,\infty } L^2_x}^2 \right]&\lesssim {\mathbb {E}} \left[ \int _{I} \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 \,\textrm{d}s\right] + \left\| \mathbb {\epsilon } u\right\| _{L^{p}_\omega L^\infty _t L^p_x}^p + 1 \\&\lesssim {\mathbb {E}} \left[ \int _{I} \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 \,\textrm{d}s\right] + {\mathbb {E}} \left[ \sup _{t\in I} {\mathcal {J}}(u) \right] + 1. \end{aligned}$$

Lastly, using the equivalence \(\varphi '(t) t \eqsim \varphi (t)\),

$$\begin{aligned} {\mathbb {E}} \left[ \left\| V(\mathbb {\epsilon } u)\right\| _{L^2_t L^2_x}^2 \right]&= {\mathbb {E}} \left[ \int _I \int _{\mathcal {O}} \varphi _{\kappa }'(\left| \mathbb {\epsilon } u\right| ) \left| \mathbb {\epsilon } u\right| \,\textrm{d}x \,\textrm{d}s \right] \\&\eqsim {\mathbb {E}} \left[ \int _I \int _{\mathcal {O}} \varphi _{\kappa }(\left| \mathbb {\epsilon } u\right| ) \,\textrm{d}x \,\textrm{d}s \right] = {\mathbb {E}} \left[ \int _I {\mathcal {J}}(u) \,\textrm{d}s \right] . \end{aligned}$$

All together, we proved

$$\begin{aligned} \left\| V(\mathbb {\epsilon } u)\right\| _{L^2_\omega B^{1/2}_{2,\infty } L^2_x}^2 \lesssim {\mathbb {E}} \left[ \int _{I} \left\| \Pi _{\text {div}} \text {div}S(\mathbb {\epsilon } u)\right\| _{L^2_x}^2 \,\textrm{d}s\right] + {\mathbb {E}} \left[ \sup _{t\in I} {\mathcal {J}}(u) \right] + 1. \end{aligned}$$

\(\square \)

Remark 23

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. [61].

In particular, if one could prove increased integrability, e.g. \( \Pi _{\text {div}}\text {div}S(\mathbb {\epsilon } u) \in L^2_{\omega } L^{q}_t L^2_x\) for some \(q > 2\), then it would also be possible to show \(V(\mathbb {\epsilon } u) \in L^2_\omega B^{1/2}_{q,\infty } L^2_x\).

The temporal regularity enables a similar error analysis for time discretizations of (1.1) as done in e.g. [14, 26] for the p-Laplace system.