claimClaim remarkRemark hypothesisHypothesis assumptionAssumption

    Analysis of a dilute polymer modelMarvin Fritz, Endre Süli, Barbara Wohlmuth

    # Analysis of a dilute polymer model with a time-fractional derivative††thanks: Submitted to the editors .
    MF is supported by the State of
    Upper Austria.

    Marvin Fritz
    Computational Methods for PDEs, Johann Radon Institute for Computational and Applied Mathematics, Altenberger Str. 69, 4040 Linz, Austria ().
<EMAIL_ADDRESS>

    Endre Süli
    Mathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Road, Oxford OX2 6GG, United Kingdom ().
<EMAIL_ADDRESS>

    Barbara Wohlmuth
    Department of Mathematics, Technical University of Munich, Boltzmannstr. 3,
    85748 Garching bei München, Germany ()
<EMAIL_ADDRESS>

    ###### Abstract

    We investigate the well-posedness of a coupled Navier–Stokes–Fokker–Planck system with a time-fractional derivative. Such systems arise in the kinetic theory of dilute solutions of polymeric liquids, where the motion of noninteracting polymer chains in a Newtonian solvent is modelled by a stochastic process exhibiting
    power-law waiting time, in order to capture subdiffusive processes associated with non-Fickian diffusion.
    We outline the derivation of the model from a subordinated Langevin equation. The elastic properties of the polymer molecules immersed in the solvent are modelled by a finitely extensible nonlinear elastic (FENE) dumbbell model, and the drag term in the Fokker–Planck equation is assumed to be corotational.
    We prove the global-in-time existence of large-data weak solutions to this time-fractional model of order $\alpha\in(\tfrac{1}{2},1)$, and derive an energy inequality satisfied by weak solutions.

    ###### keywords:

    Existence of weak solutions; Navier–Stokes–Fokker–Planck system; dilute polymer model; Hookean and FENE-type bead-spring-chain model; Riemann–Liouville fractional derivative; time-fractional PDE

    35Q30, 35Q84, 35R11, 60G22, 82C31, 82D60

    ##
    1 Introduction

    This paper is concerned with the existence of weak solutions to a system of nonlinear partial differential equations that arises in the kinetic theory of dilute solutions of polymeric fluids. Within this class of models we focus on finitely-extensible nonlinear elastic, FENE-type, dumbbell models with a corotational drag term. In contrast to previous literature on the analysis of these models we assume power law waiting times in the derivation of the system, which results in the appearance of a time-fractional derivative in the Fokker–Planck equation describing the evolution of the probability density function. This raises new questions about the study of well-posedness, and we provide rigorous results concerning the existence of global-in-time weak solutions to the system of partial differential equations featuring in the model.

    Dilute polymer models are derived and extensively described in the monograph [16] and in the book by Öttinger [58]; see also [67] for a mathematically rigorous derivation of the Hookean bead-spring-chain model from Brownian dynamics. We also refer to the papers [42, 31] for a comparison of several FENE-type dumbbell models. Such systems are of microscopic-macroscopic type since they involve a coupling of the (macroscopic) Navier–Stokes equations for the description of incompressible fluid flow and the Fokker–Planck equation for the microscopic processes associated with
    the statistical properties of polymer molecules immersed in the fluid.
    Concerning the weak and strong well-posedness of FENE-type models, we refer to the works [33, 35, 51, 73, 63]. More general dilute polymer models are analyzed in [6, 7, 8, 11, 10]. Further, we mention the papers [45, 46, 64, 50, 6, 9, 21, 44, 18], which, similarly to the discussion herein, are concerned with dumbbell models that assume a corotational drag term in the Fokker–Planck equation. In such models it is supposed that polymer molecules are not stretched by the surrounding solvent, although they are allowed to rotate without stretching; see, for example, [36].

    Time-fractional differential equations have been the focus of considerable attention in the mathematical and engineering literature in recent years. Such equations are nonlocal in time and have an innate history effect. They are of relevance in applications where memory effects are present and hereditary properties of materials are studied; see, for example, the textbooks on viscoelasticity [49, 71], hydrology [66], financial economics [23], and mechanical processes [2, 61]. The time-fractional Fokker–Planck system, in particular, allows subdiffusive behaviour and has been previously studied in [54, 55, 56, 5, 4, 28, 29, 30, 37] with regards to its derivation and applicability. The articles [62, 38, 39, 40, 41] have investigated the numerical analysis and the simulation of solutions to the time-fractional Fokker–Planck equation. The time-fractional model considered herein has been explored computationally in [15], albeit in the simpler setting of a linear (Hookean) elastic spring force instead of the FENE spring model that we study here.

    We employ a spatial Galerkin approximation in conjunction with a compactness argument to prove the existence of weak solutions to the time-fractional Navier–Stokes–Fokker–Planck system under consideration. More specifically, we discretize the system in space and derive appropriate energy bounds, which then enable us to pass to the limit in the discretized system. Spatial discretizations of dilute polymer models were previously considered in [9, 13, 14]. In addition, weak solutions to time-fractional PDEs have been investigated using the Galerkin approach in the publications [26, 27, 25]. There have also been initial steps in the analysis of a decoupled time-fractional Fokker–Planck equation with time-dependent forces; see the papers [24, 53, 40, 41, 52]. However, the coupling of the time-fractional Fokker–Planck equation to the Navier–Stokes system gives rise to new technical complications, which have not been addressed previously.

    In Section 2 we derive the model from the Langevin equation assuming power-law waiting time. In this way time-fractional derivatives in the sense of Riemann–Liouville appear in the associated Fokker–Planck equation. By mimicking the technique for the derivation of the standard dumbbell model, a time-fractional Navier–Stokes–Fokker–Planck system is obtained. In Section 3 we introduce several function spaces of Sobolev-type and recall some important results from the theory of fractional derivatives, including chain inequalities and embedding theorems. In Section 4 we transform the model in order to make it amenable to the subsequent analysis. We then equip the model with suitable initial and boundary conditions and we make use of the associated Maxwellian to rescale the Navier–Stokes–Fokker–Planck system. In Section 5 we finally state and prove a theorem asserting the existence of large-data global-in-time weak solutions to the model with a time-fractional derivative of order $\alpha\in(\tfrac{1}{2},1)$.

    ##
    2 Derivation of the time-fractional FENE-type system

    In this section we derive the time-fractional Navier–Stokes–Fokker–Planck model, admitting both Hookean and FENE-type bead-spring-chains. The classical Hookean bead-spring-chain model is derived from a system of stochastic differential equations; see, for example, the articles [67, 7] and the thesis [72, Section 1.3]. The time-fractional Fokker–Planck equation is derived from a Langevin equation in [47]. We shall emphasize the differences in the derivation of the time-fractional Navier–Stokes–Fokker–Planck model and examine the steps where the time-fractional derivative is introduced.

    In this work, the time-fractional derivative of order $\alpha\in(0,1)$ is understood in the sense of Riemann–Liouville and is given by

    (1)

    $\partial_{t}^{\alpha}w(t):=\partial_{t}\int_{0}^{t}\frac{(t-s)^{-\alpha}}{\Gamma(1-\alpha)}w(s)\,\textup{ds},$

    where $\Gamma(\alpha):=\int_{0}^{\infty}t^{\alpha-1}e^{-t}\,\textup{d}t$ is Euler’s Gamma function.

    ###
    2.1 Derivation

    We shall idealize each polymer molecule as a pair of massless beads connected by a massless elastic spring. It is assumed that the resulting, so called, dumbbell is suspended in a Newtonian solvent, whose motion is governed by the incompressible Navier–Stokes equations for the velocity $u$ and the pressure $p$ of the fluid. Let us denote the position vectors (with respect to an arbitrary, but fixed, reference point in $\mathbb{R}^{3}$) of the centers of mass of the two beads at time $t$ by $x_{i}(t)\in\mathbb{R}^{3}$ for $i\in\{1,2\}$. At time $t$, the center of mass of the dumbbell is at $x_{c}(t):=\frac{1}{2}\big{(}x_{1}(t)+x_{2}(t)\big{)}$ and the elongation (or conformation) vector pointing from $x_{1}(t)$ to $x_{2}(t)$ is $q_{1}(t):=x_{2}(t)-x_{1}(t)$. We denote the vector pointing in the opposite direction by $q_{2}(t):=-q_{1}(t)$ and assume that $q_{1}(t)$ and $q_{2}(t)$ are contained, for all $t\geq 0$, within a given convex open set $D\subset\mathbb{R}^{3}$ that satisfies $0\in D$, and $-q\in D$ whenever $q\in D$.

    Three kinds of force act on the $i$-th bead in the dumbbell suspended in the fluid: a drag force (Stokes drag) arising from the motion of the bead through the solvent, an elastic spring force, and a random force, which is assumed to be Brownian, modelling the random collisions of the bead with the molecules of the surrounding Newtonian solvent. As each of the two beads is assumed to be massless, Newton’s second law implies that

    (2)

    $\text{Drag Force}_{i}+\text{Spring Force}_{i}+\text{Brownian Force}_{i}=0,\quad i=1,2.$

    To define the drag force
    acting on the $i$-th bead of the dumbbell suspended in the fluid, we apply Stokes’ law and get

    $\displaystyle-\zeta\Big{(}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}\mathrm{t}}x_{i}(t)-u(x_{i},t)\Big{)},\quad i\in\{1,2\}.$

    Here $\zeta$ is the friction coefficient and $u(\cdot,\cdot)$ stands for the fluid velocity. We note that $\zeta$ carries the SI unit $[\text{kg}/\text{s}]$ and is linear in the dynamic viscosity $\eta$. For a sphere of radius $R$ it reads $\zeta=6\pi R\eta$.

    The elastic spring force $F:D\to\mathbb{R}^{3}$ of the spring connecting the two beads is

    $F(q):=HU^{\prime}(\tfrac{1}{2}|q|^{2})q,$

    where $U$ is a nonnegative continuously differentiable potential and $H$ is the spring constant having the SI unit $[\text{kg}/\text{s}^{2}]$.
    In the case of the Hookean dumbbell model the spring force is linear and is given by $F(q)=Hq$ with $q\in D=\mathbb{R}^{3}$ and the corresponding potential is $U(s)=s$ for $s\in[0,\infty)$, while in the case of a classical FENE model one has, instead,

    $D=B_{|q_{\text{max}}|}(0),\quad F(q)=\frac{Hq}{1-|q|^{2}/|q_{\text{max}}|^{2}},\quad U(s)=-\frac{|q_{\text{max}}|^{2}}{2}\ln\!\bigg{(}1-\frac{2s}{|q_{\text{max}}|^{2}}\bigg{)}$

    for $q\in D$ and $s\in[0,\tfrac{1}{2}|q_{\text{max}}|^{2})$, where $B_{|q_{\text{max}}|}(0)$ is an open ball in $\mathbb{R}^{3}$ with radius $|q_{\text{max}}|$ centered at the origin and $|q_{\text{max}}|>0$ is a strict upper bound on the maximal extension to which a dumbbell can be stretched.

    The Brownian force acting on the $i$-th bead at time $t$ is denoted by $B_{i}(t)$ and is formally defined by

    $B_{i}(t):=\sqrt{2k_{B}\mu_{T}\zeta}\,\frac{\mathop{}\!\mathrm{d}W_{i}(t)}{\text{d}t},\quad i\in\{1,2\},$

    where $k_{B}$ is the Boltzmann constant in $[\text{kg m}^{2}/(\text{s}^{2}\text{K})]$, $\mu_{T}$ denotes the temperature in $[\text{K}]$, and $W(t)=\big{(}W_{1}(t),W_{2}(t)\big{)}^{\mathrm{T}}$ is a vector of two independent Wiener processes.
    As each Wiener process is distributed according to $\mathcal{N}(0,t)$, we find that the SI unit for $\frac{\mathop{}\!\mathrm{d}W_{i}(t)}{\text{d}t}$ is $[1/\sqrt{\text{s }}]$.

    By introducing the notation

    (3)

    $X(t):=\begin{pmatrix}x_{1}(t)\\
    x_{2}(t)\end{pmatrix},\quad b(X(t),t):=\begin{pmatrix}u(x_{1}(t),t)+\zeta^{-1}F(q_{1}(t))\\
    u(x_{2}(t),t)+\zeta^{-1}F(q_{2}(t))\end{pmatrix},$

    and dividing Eq. 2 by $\zeta$, we get as impulse balance
    the so-called Langevin equation, an It$\widehat{\rm o}$ stochastic differential equation of the form

    $\mathop{}\!\mathrm{d}X(t)=b(X(t),t)\,\textup{d}t+\sqrt{\frac{2k_{B}\mu_{T}}{\zeta}}\mathop{}\!\mathrm{d}W(t),\quad t\geq 0.$

    The partial differential equation describing the evolution of the probability density function
    $\tilde{\psi}:=\tilde{\psi}(x_{1},x_{2},t)$
    of the random variable $X(t)$ is the standard Fokker–Planck equation [60]:

    (4)

    $\partial_{t}\tilde{\psi}=\sum_{i=1}^{2}\left(-\textup{div}_{x_{i}}\!\big{(}b_{i}((x_{1}^{\rm{T}},x_{2}^{\rm{T}})^{\rm T},t)\,\tilde{\psi}\big{)}+\frac{k_{B}\mu_{T}}{\zeta}\Delta_{x_{i}}\tilde{\psi}\right),$

    with $x_{1}$ and $x_{2}$ considered to be column-vectors in $\mathbb{R}^{3}$.

    At this point, we deviate from the usual derivation of the Fokker–Planck equation for the evolution of the probability density function of the stochastic process $X$
    and introduce a subordination of the Langevin equation. This allows us to model trapping events to the motion of the particles. For $\alpha\in(0,1)$, let $U_{\alpha}(\cdot)$ be an $\alpha$-dependent subordinator with expectation $\mathbb{E}[e^{-\lambda U_{\alpha}(\tau)}]=\text{exp}(-\tau\Phi_{\alpha}(\lambda))$, where

    $\Phi_{\alpha}(\lambda):=\tau_{0}^{\alpha-1}\lambda^{\alpha}$

    is the so-called Laplace exponent and $\tau_{0}$ is a characteristic time-scale (to be fixed). We note that the limiting value of $\alpha=1$ results in the standard integer-order case. The inverse subordinator $S_{\alpha}^{t}$, defined as the first-passage time of $U_{\alpha}$, is then given by

    $S_{\alpha}^{t}:=\inf_{\tau>0}\{\tau\,:\,U_{\alpha}(\tau)>t\}.$

    Suppose that $Y_{\alpha}(t)$ is a solution of the stochastic differential equation

    $\mathrm{d}Y_{\alpha}(t)=b(Y_{\alpha}(t),U_{\alpha}(t))\mathrm{d}t+\sqrt{\frac{2k_{B}\mu_{T}}{\zeta}}\mathrm{d}W(t),\quad t\geq 0.$

    Define $X(t):=Y_{\alpha}(S_{\alpha}^{t})$. It then follows from Theorem 1 in [48] that if $b$ is twice continuously differentiable with respect to its
    variables and satisfies the usual Lipschitz condition, then the probability density function of the process $X$ is a solution of the time-fractional
    Fokker–Planck equation
    resulting from the replacement of $\tilde{\psi}$ on the right-hand side of (4) by $\tau_{0}^{1-\alpha}\partial_{t}^{1-\alpha}\tilde{\psi}$, i.e.,

    $\partial_{t}\tilde{\psi}=\sum_{i=1}^{2}\left(-\textup{div}_{x_{i}}\!\big{(}b_{i}((x_{1}^{\rm{T}},x_{2}^{\rm{T}})^{\rm T},t)\,\tau_{0}^{1-\alpha}\partial_{t}^{1-\alpha}\tilde{\psi}\big{)}+\frac{k_{B}\mu_{T}}{\zeta}\Delta_{x_{i}}\tau_{0}^{1-\alpha}\partial_{t}^{1-\alpha}\tilde{\psi}\right).$

    Next, we perform the linear change of variables $(x_{1},x_{2})\mapsto(\frac{1}{2}(x_{1}+x_{2}),x_{2}-x_{1})=:(x,q)$, whereby we have identified a point $x\in\Omega\subset\mathbb{R}^{3}$ in the (macroscopic) flow domain $\Omega$ with the center of mass of the dumbbell, and have, as before, denoted by $q$ the vector pointing from $x_{1}$ to $x_{2}$.
    By recalling the definition Eq. 3 of $b(\cdot,\cdot)$ and setting $\psi(x,q,t):=\tilde{\psi}(x-\tfrac{1}{2}q,x+\tfrac{1}{2}q,t)$, we find that

    (5)

    $\displaystyle\partial_{t}\psi\,+$
    $\displaystyle\tau_{0}^{1-\alpha}\textup{div}_{x}\bigg{(}\frac{u(x-\tfrac{1}{2}q,t)+u(x+\tfrac{1}{2}q,t)}{2}\partial_{t}^{1-\alpha}\psi\bigg{)}$

    $\displaystyle\quad+\tau_{0}^{1-\alpha}\textup{div}_{q}\bigg{(}\big{(}u(x+\tfrac{1}{2}q,t)-u(x-\tfrac{1}{2}q,t)\big{)}\partial_{t}^{1-\alpha}\psi-\frac{2F(q)}{\zeta}\partial_{t}^{1-\alpha}\psi\bigg{)}$

    $\displaystyle=\frac{k_{B}\mu_{T}\tau_{0}^{1-\alpha}}{2\zeta}\Delta_{x}\partial_{t}^{1-\alpha}\psi+\frac{2k_{B}\mu_{T}\tau_{0}^{1-\alpha}}{\zeta}\Delta_{q}\partial_{t}^{1-\alpha}\psi.$

    To proceed, we assume ‘local homogeneity’, i.e., that the spatial variation of the velocity field over the microscopic length-scale of a single dumbbell is negligibly small. Consequently, the arithmetic mean $\big{(}u(x-\frac{1}{2}q,t)+u(x+\frac{1}{2}q,t)\big{)}/2$ can be approximated by $u(x,t)$ in the second term on the left-hand side of Eq. 5. In the case of the third term, we use Taylor expansion of $u$ about the point $x$ to obtain

    (6)

    $\displaystyle u(x+\tfrac{1}{2}q,t)-u(x-\tfrac{1}{2}q,t)$
    $\displaystyle=\nabla_{\!x}u(x,t)q+\mathcal{O}(|q|^{3})$

    $\displaystyle=\Big{(}\sigma\big{(}u(x,t)\big{)}+\omega\big{(}u(x,t)\big{)}\Big{)}q+\mathcal{O}(|q|^{3}),$

    where we have further split the gradient of $u$ into its symmetric and antisymmetic parts, which are, respectively, defined as follows:

    (7)

    $\sigma(u):=\frac{\nabla_{\!x}u+(\nabla_{\!x}u)^{\mathrm{T}}}{2},\qquad\omega(u):=\frac{\nabla_{\!x}u-(\nabla_{\!x}u)^{\mathrm{T}}}{2}.$

    In the approximation Eq. 6, we omit the $\mathcal{O}(|q|^{3})$ term and further, we also omit the symmetric part of the gradient of $u$, i.e., we consider the, so called, corotational model.
    While the omission of the $\mathcal{O}(|q|^{3})$ term from Eq. 6 can be justified on the grounds that $|q|\ll 1$, our omission of the term $\sigma(u(x,t))$ from the additive decomposition $\nabla_{x}u(x,t)=\sigma(u(x,t))+\omega(u(x,t))$ is for purely technical reasons.
    In the time-fractional corotational model considered here, polymer molecules are therefore allowed to rotate, but they are forced to do so without stretching. This modelling assumption weakens the coupling between the Navier–Stokes equation and the (time-fractional) Fokker–Planck equation; for example, if the initial datum for $\psi$ happens to be spherically symmetric with respect to $q$ and independent of the spatial variable $x$, then the Fokker–Planck equation is decoupled from the Navier–Stokes equation. In this respect, the model that we study here is no
    different from the corotational model considered (in the case of $\alpha=1$) in the works [45, 46, 64, 50, 6, 9, 21, 44, 18].
    Thus, we consider the following corotational Fokker–Planck equation with a time-fractional derivative:

    (8)

    $\displaystyle\partial_{t}\psi+\tau_{0}^{1-\alpha}\textup{div}_{x}(u\partial_{t}^{1-\alpha}\psi)+\tau_{0}^{1-\alpha}\textup{div}_{q}\big{(}\omega(u)q\partial_{t}^{1-\alpha}\psi\big{)}$

    $\displaystyle\quad=\frac{k_{B}\mu_{T}\tau_{0}^{1-\alpha}}{2\zeta}\Delta_{x}\partial_{t}^{1-\alpha}\psi+\frac{2k_{B}\mu_{T}\tau_{0}^{1-\alpha}}{\zeta}\Delta_{q}\partial_{t}^{1-\alpha}\psi+\tau_{0}^{1-\alpha}\textup{div}_{q}\Big{(}\frac{2F(q)}{\zeta}\partial_{t}^{1-\alpha}\psi\Big{)}.$

    We shall rely on the fact that $q^{\mathrm{T}}\omega(u)q\equiv 0$ thanks to the skew-symmetry of $\omega(u)$. In the general noncorotational case $q^{\mathrm{T}}(\nabla_{x}u)q$ is of course not identically nonzero in the flow domain. For $\alpha=1$ at least, the proof of the existence of large-data global-in-time weak solutions to the general noncorotational Navier–Stokes–Fokker–Planck system relies, instead, on an entropy estimate (cf. [12]). It is this entropy estimate that needs to be replicated in the time-fractional case, for $\alpha\in(0,1)$, as a key
    ingredient of the proof of the existence of global-in-time large-data weak solutions. The analysis of the general noncorotational time-fractional model is deferred to future work.

    Let $\psi(x,q,t)$ denote from now on the probability density function for a collection of $N\gg 1$ dumbbells. As we are dealing with a dilute polymeric fluid, the polymer molecules suspended in the fluid are assumed not to interact with each other and they move without self-interaction. The function $\psi$ therefore satisfies the same partial differential equation, (8), as in the case of a single dumbbell. The only difference is in the choice of the initial datum $\psi^{0}\geq 0$ for the Fokker–Planck equation.
    A typical choice of $\psi^{0}$ in the present context is

    $\psi^{0}(x,q)=\frac{1}{N}\sum_{j=1}^{N}\alpha_{j}(q)\Psi_{j}(x),$

    where $\alpha_{j}(q)\geq 0$ for all $q\in D$ and $\int_{D}\alpha_{j}(q)\mathop{}\!\mathrm{d}q=1$, $j=1,\ldots,N$, and $\Psi_{j}\geq 0$, for all $x\in\Omega$ and $\int_{\Omega}\Psi_{j}(x)\mathop{}\!\mathrm{d}x=1$, $j=1,\ldots,N$. For example, one may choose $\Psi_{j}$ as a mollifier (i.e. a nonnegative $C^{\infty}_{0}$ approximation to the Dirac measure) concentrated at a point $z_{j}\in\Omega$, $j=1,\ldots,N$, with $\{\Psi_{j}\}_{j=1}^{N}$ forming a scaled partition of unity; here $z_{j}$ can be thought of as the initial location of the center of mass of the $j$-th dumbbell.

    For the sake of simplicity we shall confine our attention to the case when the boundary $\partial\Omega$ of the macroscopic flow domain $\Omega\subset\mathbb{R}^{3}$ has no inflow or outflow parts. The polymeric fluid under consideration is therefore confined to $\Omega$, and its velocity will be supposed to satisfy the no-slip boundary condition $u(x,t)=0$ for all $(x,t)\in\partial\Omega\times(0,T)$.
    Thus the
    total number $N\gg 1$ of polymer molecules contained in $\Omega$ remains constant in time. The number density
    $\varrho_{P}:=\varrho_{P}(x,t)$, in $[m^{-3}]$, of the polymer molecules contained in $\Omega$ is called the polymer number density and it is related to $N$ by $N=\int_{\Omega}\varrho_{P}(x,t)\mathop{}\!\mathrm{d}x$. The polymer number density is further related to the probability density function $\psi$ satisfying the Fokker–Planck equation (8) by

    (9)

    $\displaystyle\varrho_{P}(x,t)=N\int_{D}\psi(x,q,t)\mathop{}\!\mathrm{d}q,$

    with the understanding that $\psi$ has the SI unit of $[m^{-6}]$. We shall supplement the Fokker–Planck equation (8) with no-flux (homogeneous Neumann) boundary conditions on $\partial\Omega\times D\times(0,T)$ and on $\Omega\times\partial D\times(0,T)$, which will ensure that $\int_{\Omega}\int_{D}\psi(x,q,t)\mathop{}\!\mathrm{d}q\mathop{}\!\mathrm{d}x$ is constant in time, and is therefore equal to

    $\int_{\Omega}\int_{D}\psi(x,q,0)\mathop{}\!\mathrm{d}q\mathop{}\!\mathrm{d}x=\int_{\Omega}\int_{D}\psi^{0}(x,q)\mathop{}\!\mathrm{d}q\mathop{}\!\mathrm{d}x=1.$

    Having established the “microscopic” equation that describes the statistical properties of polymer molecules in the continuum, we turn our attention to the continuum mechanical “macroscopic” equations of motion of the incompressible fluid in which the polymer molecules are suspended. We note that polymeric fluids are non-Newtonian fluids and the presence of the $N\gg 1$ polymer molecules contributes an additional term, $\tau=\tau(\psi)$, to the stress tensor appearing in the balance of linear momentum equation in the Navier–Stokes system: a symmetric polymeric extra stress tensor, which we shall define below. We assume that the evolution of the velocity field $u$ and the pressure $p$ is governed by the incompressible Navier–Stokes system

    (10)

    $\displaystyle\varrho(\partial_{t}u+(u\cdot\nabla_{\!x})u)-\eta\Delta_{x}u+\nabla_{\!x}p$
    $\displaystyle=\textup{div}_{x}\tau(\psi)$

    $\displaystyle\quad\text{in }\Omega\times(0,T),$

    $\displaystyle\textup{div}_{x}u$
    $\displaystyle=0$

    $\displaystyle\quad\text{in }\Omega\times(0,T),$

    supplemented with the homogeneous Dirichlet boundary condition $u=0$ on $\partial\Omega\times(0,T)$ and the initial condition $u(0)=u^{0}$ in $\Omega$
    at $t=0$.
    As before, $\eta$ denotes the dynamic viscosity and $\varrho$ stands for the macroscopic density, which is assumed to be constant in space and time. The polymeric extra-stress tensor $\tau=\tau(\psi)$ is defined by the so-called Kramers expression, see, e.g., [16],

    (11)

    $\tau(\psi):=\varrho_{P}\big{(}\mathcal{C}(\psi)-k_{B}\mu_{T}I_{3}\big{)},$

    where, as before, $\varrho_{P}$ is the polymer number density, $k_{B}$ and $\mu_{T}$ are the Boltzmann constant and the absolute temperature, respectively, and $I_{3}$ is the $3\times 3$ identity matrix.
    Finally, the $3\times 3$ symmetric tensor $\mathcal{C}$ appearing in the expression for $\tau$ is defined by

    (12)

    $\mathcal{C}(\psi):=\frac{\int_{D}F(q)q^{\mathrm{T}}\psi\mathop{}\!\mathrm{d}q}{\int_{D}\psi\mathop{}\!\mathrm{d}q},$

    see, e.g., [67]. Because the polymer number density $\varrho_{P}$ is related to $\psi$ by (9), the polymeric extra stress tensor is simplified to

    $\tau(\psi)=\tau^{1}(\psi)+\tau^{2}(\psi),$

    where

    (13)

    $\displaystyle\tau^{1}(\psi)$
    $\displaystyle:=N\int_{D}F(q)q^{\mathrm{T}}\psi\mathop{}\!\mathrm{d}q,$

    (14)

    $\displaystyle\tau^{2}(\psi)$
    $\displaystyle:=-Nk_{B}\mu_{T}\int_{D}\psi\mathop{}\!\mathrm{d}q\ I_{3}.$

    Recall that $\psi$ is assumed to have the SI unit $[m^{-6}]$.

    We note that the definition of $\mathcal{C}(\psi)$ results in the SI unit $[kg/(s^{2}m^{2})]$ for the polymeric extra stress tensor $\tau=\tau(\psi)$.
    The tensor $\tau$ is responsible for coupling the velocity $u$ and the pressure $p$ to the probability density function $\psi$.
    Dividing the Navier–Stokes momentum equation by the macroscopic density $\varrho$ and introducing the kinematic viscosity as $\nu:=\eta/\varrho$, we arrive at the Navier–Stokes system in its dimensional form:

    (15)

    $\displaystyle\partial_{t}u+(u\cdot\nabla_{\!x})u-\nu\Delta_{x}u+\frac{1}{\varrho}\nabla_{\!x}p$
    $\displaystyle=\frac{1}{\varrho}\textup{div}_{x}\tau(\psi)$

    $\displaystyle\quad\text{in }\Omega\times(0,T),$

    $\displaystyle\textup{div}_{x}u$
    $\displaystyle=0$

    $\displaystyle\quad\text{in }\Omega\times(0,T).$

    In the next section we shall perform a nondimensionalization of the Navier–Stokes–Fokker–Planck system. In order to distinguish a nondimensionalized
    quantity from its dimensional form, we shall use the subscript dim for dimensional quantities in cases where confusion might arise; so, for example, we shall write $\tau_{\text{dim}}$ to indicate the original dimensional form of the polymeric extra stress tensor, while $\tau$ will henceforth signify its nondimensionalized form. Similarly $\mathcal{C}_{\mathrm{dim}}$ will denote the original dimensional form of $\mathcal{C}$ (cf. (12)), while $\mathcal{C}$ will signify its form following nondimensionalization.

    ###
    2.2 Nondimensionalization

    Next, we transform the Navier–Stokes–Fokker–Planck system, see (15) and (8), into its nondimensionalized form. To this end, we define the quantities $\widehat{x}$, $\widehat{t}$ and $\widehat{u}$ by setting

    $x=L_{0}\widehat{x},\quad t=T_{0}\widehat{t},\quad u(x,t)=U_{0}\widehat{u}(\widehat{x},\widehat{t}),$

    where $L_{0}$ and $T_{0}$ stand for the characteristic macroscopic length-scale and time-scale, respectively, and $U_{0}$ denotes the characteristic velocity of the macroscopic flow.
    In a similar manner, we introduce $\widehat{q}$ by letting $q=l_{0}\widehat{q}$, where $l_{0}$ is a characteristic microscopic length-scale, recall that $\psi$ was assumed to have the SI unit $[m^{-6}]$, and we define the dimensionless quantity $\widehat{\psi}$ by

    (16)

    $\displaystyle\widehat{\psi}(\widehat{x},\widehat{q},\widehat{t}):=(L_{0}l_{0})^{3}\,\psi(x,q,t).$

    We point out that scaling $\psi$ by a constant does not affect the definition of the $\mathcal{C}_{\text{dim}}(\psi)$ given by (12); i.e., for any positive constant $a$ we have $\mathcal{C}_{\text{dim}}(a\psi)=\mathcal{C}_{\text{dim}}(\psi)$.
    Scaling $\psi$ by a positive constant, however, affects its relationship to the polymer number density $\varrho_{P}(x,t)$. To avoid the nondimensionalization of $\varrho_{P}(x,t)$, we do not consider the definition (11) of $\tau_{\text{dim}}(\psi)$ but of its equivalent form given by (13) and (14).

    Consequently, we obtain from the time-fractional Fokker–Planck equation (8) the following partial differential equation

    $\displaystyle\frac{1}{T_{0}}\partial_{\widehat{t}}\widehat{\psi}+\frac{\tau_{0}^{1-\alpha}U_{0}}{T_{0}^{1-\alpha}L_{0}}\textup{div}_{\widehat{x}}(\widehat{u}(\psi)\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi})+\frac{\tau_{0}^{1-\alpha}U_{0}}{T_{0}^{1-\alpha}L_{0}}\textup{div}_{\widehat{q}}\big{(}\omega(\widehat{u}(\psi))\widehat{q}\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}\big{)}$

    (17)

    $\displaystyle\quad=\frac{k_{B}\mu_{T}\tau_{0}^{1-\alpha}}{2\zeta T_{0}^{1-\alpha}L_{0}^{2}}\Delta_{\widehat{x}}\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}+\frac{2k_{B}\mu_{T}\tau_{0}^{1-\alpha}}{\zeta l_{0}^{2}T_{0}^{1-\alpha}}\Delta_{\widehat{q}}\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}+\frac{2H\tau_{0}^{1-\alpha}}{\zeta T_{0}^{1-\alpha}}\textup{div}_{\widehat{q}}(\widehat{F}(\widehat{q})\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi})$

    defined on the nondimensionalized domain $\widehat{\Omega}\times\widehat{D}\times(0,\widehat{T})$. Here we have used the notation $\widehat{u}(\psi)$ to indicate that the velocity field is understood to depend on the original (dimensional) probability density function $\psi$, and will only be expressed as a function depending on the nondimensionalized probability density function $\widehat{\psi}$ in the next step.
    We set the characteristic macroscopic and microscopic time-scales to, respectively, $T_{0}:=L_{0}/U_{0}$, $\tau_{0}:=T_{0}$ and the nondimensionalized force in case of the standard FENE type model to

    (18)

    $\displaystyle\widehat{F}(\widehat{q});=\frac{\widehat{q}}{1-|\widehat{q}|^{2}/|\widehat{q}_{\max}|^{2}},\qquad\widehat{q}_{\max}:=q_{\max}/l_{0}.$

    The prefactors of the last two terms on the right-hand side of (2.2) are equal if the microscopic length-scale $\ell_{0}$ is defined as follows:

    (19)

    $\displaystyle\ell_{0}^{2}:=\frac{k_{B}\mu_{T}}{H}.$

    Next we introduce the Deborah number $\lambda$, defined as the ratio of the time it takes for the material to adjust to applied stresses/deformations and the characteristic time scale, $T_{0}$, by

    $\lambda:=\frac{\zeta/(4H)}{T_{0}}=\frac{\zeta}{4HT_{0}}=\frac{\zeta\ell_{0}^{2}}{4k_{B}\mu_{T}T_{0}}.$

    Finally, we define the nondimensional parameter

    $\varepsilon:=\frac{k_{B}\mu_{T}}{2\zeta U_{0}L_{0}},$

    whose numerator and denominator both carry the SI unit of $\mathrm{J=[kg\,m^{2}/s^{2}]}$ corresponding to energy. Thus, we obtain from (2.2) by multiplying it with $T_{0}$ the following nondimensional equation:

    $\displaystyle\partial_{\widehat{t}}\widehat{\psi}+(\widehat{u}(\psi)\cdot\nabla_{\widehat{x}})\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}+\textup{div}_{\widehat{q}}\,\big{(}\omega(\widehat{u}(\psi)){\widehat{q}}\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}\big{)}\,\qquad\qquad\qquad$

    (20)

    $\displaystyle=\varepsilon\Delta_{\widehat{x}}\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}+\frac{1}{2\lambda}\textup{div}_{\widehat{q}}\,(\nabla_{\widehat{q}}\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}+\widehat{F}(\widehat{q})\,\partial_{\widehat{t}}^{1-\alpha}\widehat{\psi}).$

    We note in passing that the ratio of the diffusion coefficients $\varepsilon$ and $1/(2\lambda)$ featuring in the equation (2.2) is equal to $(\ell_{0}/(2L_{0}))^{2}$, which is $\ll 1$, so the first-term on the right-hand side of (2.2), called the center-of-mass diffusion term, is frequently neglected in practical considerations. Crucially, we shall retain this term in the equation and will continue to work with a strictly positive center-of-mass diffusion coefficient $\varepsilon$, as is implied by the derivation of the Fokker–Planck equation (2.2) performed above.

    In the same manner, we multiply the Navier–Stokes system (15) by $T_{0}/U_{0}$ and arrive at the following nondimensionalized system posed on $\widehat{\Omega}\times(0,\widehat{T})$:

    (21)

    $\displaystyle\partial_{\widehat{t}}\widehat{u}+(\widehat{u}\cdot\nabla_{\widehat{x}})\widehat{u}-\frac{1}{\text{Re}}\Delta_{\widehat{x}}\widehat{u}+\nabla_{\widehat{x}}\widehat{p}$
    $\displaystyle=\textup{div}_{\widehat{x}}\widehat{\tau}(\widehat{\psi}),$

    $\displaystyle\textup{div}_{\widehat{x}}\widehat{u}$
    $\displaystyle=0,$

    where Re stands for the Reynolds number and $\widehat{p}$ for a scaled pressure. Defining $\tau(\widehat{\psi})$ in a proper way allows us to replace $\widehat{u}(\psi)$ in our nondimensionalized Fokker–Planck equation by $\widehat{u}(\widehat{\psi})$ or, in compact notation, by $\widehat{u}$.
    Thanks to the reformulation of the stress, we reconsider the two terms $\widehat{\tau}^{1}(\widehat{\psi})$ and $\widehat{\tau}^{2}(\widehat{\psi})$, appearing in $\widehat{\tau}(\widehat{\psi}):=\widehat{\tau}^{1}(\widehat{\psi})+\widehat{\tau}^{2}(\widehat{\psi})$, separately. The definition (13) of $\tau^{1}(\psi)$ in combination with the nondimensionalization steps and the definition of $\widehat{F}$ and $\widehat{\psi}$, see (18) and (16), respectively yields

    (22)

    $\displaystyle\frac{1}{\varrho}\frac{T_{0}}{U_{0}}\frac{1}{L_{0}}Hl_{0}^{2}l_{0}^{3}$
    $\displaystyle\int_{\widehat{D}}\widehat{F}(\widehat{q})\widehat{q}^{\mathrm{T}}\widehat{\psi}\mathop{}\!\mathrm{d}\widehat{q}\,(L_{0}l_{0})^{-3}N$

    $\displaystyle=\ \frac{N}{\varrho L_{0}^{3}}\frac{k_{B}\mu_{T}}{U_{0}^{2}}\int_{\widehat{D}}\widehat{F}(\widehat{q})\widehat{q}^{\mathrm{T}}\widehat{\psi}(\widehat{q})\mathop{}\!\mathrm{d}\widehat{q}=\gamma\,\widehat{\mathcal{C}}(\widehat{\psi})=:\widehat{\tau}^{1}(\widehat{\psi}),$

    with the dimensionless quantities $\gamma$ and $\widehat{\mathcal{C}}(\widehat{\psi})$ being defined by

    $\gamma:=\frac{k_{B}\mu_{T}N}{\varrho U_{0}^{2}L_{0}^{3}}\quad\mbox{and}\quad\widehat{\mathcal{C}}(\widehat{\psi}):=\int_{\widehat{D}}\widehat{F}(\widehat{q})\widehat{q}^{\mathrm{T}}\widehat{\psi}\mathop{}\!\mathrm{d}\widehat{q}.$

    To define $\widehat{\tau}^{2}(\widehat{\psi})$, we start from the scaled macroscopic momentum balance equation and the definition (14):

    (23)

    $\displaystyle\frac{1}{\varrho}\frac{T_{0}}{U_{0}}\frac{1}{L_{0}}(-k_{b}\mu_{T})l_{0}^{3}$
    $\displaystyle\int_{\widehat{D}}\widehat{\psi}\mathop{}\!\mathrm{d}\widehat{q}\,(L_{0}l_{0})^{-3}N\ I_{3}$

    $\displaystyle=\ -\frac{N}{\varrho L_{0}^{3}}\frac{k_{B}\mu_{T}}{U_{0}^{2}}\int_{\widehat{D}}\widehat{\psi}\mathop{}\!\mathrm{d}\widehat{q}\ I_{3}=\gamma\int_{\widehat{D}}\widehat{\psi}\mathop{}\!\mathrm{d}\widehat{q}\ I_{3}=:\widehat{\tau}^{2}(\widehat{\psi}).$

    For the sake of notational simplicity we drop from now on all $\,\widehat{\cdot}\,$ symbols.

    We close this section by recalling from (21) and (2.2) our nondimensionalized model problem:

    (24)

    $\displaystyle\partial_{t}\psi+(u\cdot\nabla_{\!x})\partial_{t}^{1-\alpha}\psi+\textup{div}_{q}\big{(}\omega(u)$
    $\displaystyle q\partial_{t}^{1-\alpha}\psi\big{)}$

    $\displaystyle-\tfrac{1}{2\lambda}\textup{div}_{q}(\nabla_{\!q}\partial_{t}^{1-\alpha}\psi+\partial_{t}^{1-\alpha}\psi)$
    $\displaystyle=\varepsilon\partial_{t}^{1-\alpha}\Delta_{x}\psi$

    $\displaystyle\quad\text{in }\Omega\times D\times(0,T),$

    $\displaystyle\partial_{t}u+(u\cdot\nabla_{\!x})u-\frac{1}{\text{Re}}\Delta_{x}u+\nabla_{\!x}p$
    $\displaystyle=\textup{div}_{x}\tau(\psi)$

    $\displaystyle\quad\text{in }\Omega\times(0,T),$

    $\displaystyle\textup{div}_{x}u$
    $\displaystyle=0$

    $\displaystyle\quad\text{in }\Omega\times(0,T),$

    and we supplement this system of equations with the following boundary conditions:

    $\displaystyle\left(\tfrac{1}{2\lambda}(\nabla_{\!q}\partial_{t}^{1-\alpha}\psi+U^{\prime}q\partial_{t}^{1-\alpha}\psi)-\omega(u)q\partial_{t}^{1-\alpha}\psi\right)\cdot n_{\partial D}$
    $\displaystyle=0\qquad\text{on }\Omega\times\partial D\times(0,T),$

    $\displaystyle\varepsilon\nabla_{\!x}\partial_{t}^{1-\alpha}\psi\cdot n_{\partial\Omega}$
    $\displaystyle=0\qquad\text{on }\partial\Omega\times D\times(0,T),$

    $\displaystyle u$
    $\displaystyle=0\qquad\text{on }\partial\Omega\times(0,T).$

    The nondimensional initial data $u^{0}(x)$ and $\psi^{0}(x,q)$ for the velocity field and the probability density function, respectively, are obtained from their dimensional counterparts. We note that the scaling (16) ensures
    that $\int_{\Omega}\int_{D}\psi^{0}(x,q)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}q=1$; this then guarantees, thanks to the homogeneous Neumann boundary conditions on $\psi$, that $\int_{\Omega}\int_{D}\psi(x,q,t)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}q=1$ for all $t\geq 0$.

    ##
    3 Mathematical preliminaries

    In this section, we introduce some useful definitions and results regarding the fractional derivative in the sense of Riemann–Liouville and recall the Aubin–Lions lemma, which is a key result featuring in proofs of existence of weak solutions to nonlinear PDEs based on compactness arguments.

    For a Hilbert space $H$ with inner product $(\cdot,\cdot)_{H}$ and norm $\|\cdot\|_{H}$, we shall denote the duality pairing between $H$ and its dual space $H^{\prime}$ by $\langle\cdot,\cdot\rangle_{H}$. We shall denote the inner product on the Bochner space $L^{2}(0,T;H)$ by $(\cdot,\cdot)_{L^{2}H}$, and we shall write $(\cdot,\cdot)_{L^{2}_{t}H}$ when in this inner product the temporal interval of integration is $(0,t)$ for some $t\in(0,T)$ rather than the complete interval $(0,T)$, i.e.,

    $(u,v)_{L^{2}_{t}H}:=\int_{0}^{t}(u(s),v(s))_{H}\,\text{d}s\qquad\forall\,u,v\in L^{2}(0,T;H).$

    The norm induced by this inner product will be denoted by $\|\cdot\|_{L^{2}_{t}H}$.

    ###
    3.1 Riemann–Liouville kernels

    The Riemann–Liouville kernel function $g_{\alpha}$ of order $\alpha$ is defined by $g_{\alpha}(t):=t^{\alpha-1}/\Gamma(\alpha)$, $t\in(0,T)$, for $\alpha>0$ and $g_{0}(t):=\delta_{0}(t)$ (the Dirac distribution concentrated at $0$) for $\alpha=0$. We observe that $g_{\alpha}\in L^{p}(0,T)$ for any $\alpha\in(1-1/p,1)$ and $p\in[1,\infty)$, and the kernel function satisfies the following semigroup property; see [22, Theorem 2.4]:

    (25)

    $g_{\alpha}*g_{\beta}=g_{\alpha+\beta}\qquad\forall\,\alpha,\beta\geq 0.$

    We note that when $\alpha\in(0,1)$, one can bound the $L^{p}(0,t)$-norm of a function $u:(0,T)\to\mathbb{R}$ by its convolution with $g_{\alpha}$ as follows: for any $t\in(0,T]$, we have that

    (26)

    $\displaystyle\|u\|_{L^{p}(0,t)}^{p}:=\int_{0}^{t}|u(s)|^{p}\,\textup{ds}$
    $\displaystyle\leq t^{1-\alpha}\int_{0}^{t}(t-s)^{\alpha-1}|u(s)|^{p}\,\textup{ds}$

    $\displaystyle\leq T^{1-\alpha}\Gamma(\alpha)\big{(}g_{\alpha}*|u|^{p}\big{)}(t).$

    This implies that the space

    $L^{p}_{\alpha}(0,T):=\big{\{}u:(0,T)\to\mathbb{R}:\sup_{t\in(0,T)}(g_{\alpha}*|u|^{p})(t)<\infty\big{\}},$

    is indeed a subspace of $L^{p}(0,T)$.
    If the order $\alpha$ of the kernel function $g_{\alpha}$ is larger than 1, then one can exploit the semigroup property of the kernel and apply Young’s convolution inequality (cf. Lemma 3.2 in [57]) as follows:

    $(g_{1+\alpha}*u)(t)=(g_{1}*g_{\alpha}*u)(t)=\int_{0}^{t}(g_{\alpha}*u)(s)\,\textup{ds}\leq\|g_{\alpha}\|_{L^{1}(0,t)}\|u\|_{L^{1}(0,t)},$

    for any $u\in L^{1}(0,T)$ and any $t\in(0,T]$.

    ###
    3.2 Time-fractional derivative

    We can rewrite the definition of the Riemann–Liouville derivative stated in Eq. 1 in a compact form by using the convolution operator $*$ as $\partial_{t}^{\alpha}w=\partial_{t}(g_{1-\alpha}*w)$.
    We refer to the classical textbooks [22, 3] and the newer monographs [32, 20] regarding fractional calculus and fractional differential equations.

    We define the fractional Riemann–Liouville–Bochner space for $\alpha\in(0,1)$ and $p\in[1,\infty)$ on $(0,T)$ with values in $H$ by

    $\mathcal{W}^{\alpha,p}(0,T;H):=\big{\{}u\in L^{p}(0,T;H):g_{1-\alpha}*u\in W^{1,p}(0,T;H)\big{\}}.$

    Here, the convolution $\ast$ is of course understood to be with respect to the temporal variable $t\in(0,T)$. In the limit,
    when $\alpha=1$, we have that $g_{1-\alpha}=g_{0}=\delta$, and then

    $\mathcal{W}^{1,p}(0,T;H):=W^{1,p}(0,T;H):=\big{\{}u\in L^{p}(0,T;H):\partial_{t}u\in L^{p}(0,T;H)\big{\}}.$

    However for $0<\alpha<1$,
    the Riemann–Liouville space $\mathcal{W}^{\alpha,p}(0,T;H)$ differs from the fractional-order Sobolev–Bochner space

    $W^{\alpha,p}(0,T;H):=\Big{\{}u\in L^{p}(0,T;H):(s,t)\mapsto\tfrac{\|u(t)-u(s)\|_{H}}{|t-s|^{\alpha+1/p}}\in L^{p}((0,T)\times(0,T))\Big{\}},$

    which can be confirmed by noting that the function $g_{\alpha}$ is an element of $\mathcal{W}^{\alpha,p}(0,T):=\mathcal{W}^{\alpha,p}(0,T;\mathbb{R})$ for $\alpha\in(1-\tfrac{1}{p},1)$ but not of $W^{\alpha,p}(0,T)$; see [19, Proposition 3.13].

    ######
    Remark 3.1.

    Even though the space $\mathcal{W}^{\alpha,p}(0,T)$ is not a subspace of the Sobolev–Slobodeckiĭ space $W^{\alpha,p}(0,T)$, it is nevertheless continuously embedded into $C([0,T])$, the space of uniformly continuous functions defined on $[0,T]$, for $\alpha\in(1-\frac{1}{p},1]$ and $p\in[1,\infty)$; see, [19, Remark 6.2].

    We also introduce the following Riemann–Liouville space incorporating a homogeneous initial condition at $t=0$, albeit in a somewhat nonstandard manner:

    $\displaystyle\mathcal{W}^{\alpha,p}_{0}(0,T;H)$
    $\displaystyle:=\big{\{}u\in\mathcal{W}^{\alpha,p}(0,T;H):(g_{1-\alpha}*u)(0)=0\big{\}}.$

    We note that the function $g_{1-\alpha}*u:[0,T]\to H$ has a well-defined trace at $t=0$ (even when the function $u$ itself might not have one) thanks to the continuous embedding

    $g_{1-\alpha}*u\in W^{1,p}(0,T;H)\hookrightarrow AC([0,T];H).$

    For a given element $z\in H$, the convolution $g_{1-\alpha}*z$ should be understood to mean the function $t\mapsto(g_{1-\alpha}*g_{1})(t)z\in H$; recall that $g_{1}(t)\equiv 1$ for all $t\geq 0$. Thus, $z\in H$ is in this context now thought of as the mapping $t\mapsto g_{1}(t)z\in\mathcal{W}^{\alpha,p}(0,T;H)$, for $\alpha\in(0,1)$, $p\in[1,\infty)$ and $0<\alpha p<1$, or if $\alpha=1$ and $p\in[1,\infty)$.
    Thanks to the semigroup property (25) we then have that

    $t\mapsto(g_{1-\alpha}*z)(t)=z\,g_{2-\alpha}(t)=\frac{z}{\Gamma(2-\alpha)}t^{1-\alpha}\in C([0,T];H)$

    for any $\alpha\in[0,1]$. Thus, for $\alpha\in(0,1]$, $p\in[1,\infty)$ and $z\in H$ we define the following ‘translated’ Riemann–Liouville space:

    (27)

    $\mathcal{W}^{\alpha,p}_{z}(0,T;H):=\big{\{}u\in L^{p}(0,T;H):u-z\in\mathcal{W}^{\alpha,p}(0,T;H),~{}(g_{1-\alpha}*u)(0)=0\big{\}}.$

    Note that if $\alpha\in(0,1)$, $p\in[1,\infty)$ and $0<\alpha p<1$, or if $\alpha=1$ and $p\in[1,\infty)$, then $u-z\in\mathcal{W}^{\alpha,p}(0,T;H)$ if, and only if $u\in\mathcal{W}^{\alpha,p}_{0}(0,T;H)$, and therefore, for such $\alpha$ and $p$ we have that $\mathcal{W}^{\alpha,p}_{z}(0,T;H)=\mathcal{W}^{\alpha,p}_{0}(0,T;H)$ irrespective of the choice of $z\in H$.

    Next, we state an inverse convolution (or deconvolution) property. Its name stems from the fact that convolution with the kernel $g_{\alpha}$ acts as an inverse mapping on the operator of taking $\alpha$-th fractional derivative, up to a term that involves the initial value at $t=0$.

    ######
    Lemma 3.2 (Inverse convolution).

    Let $\alpha\in(0,1]$ and $p\in[1,\infty)$. Suppose further that $H$ is a Hilbert space and $z\in H$. Then, for any $t\in(0,T)$, we have the following equalities:

    (28)

    $\displaystyle(g_{\alpha}*\partial_{t}^{\alpha}u)(t)$
    $\displaystyle=u(t)-(g_{1-\alpha}*u)(0)g_{\alpha}(t)\quad$
    $\displaystyle\forall\,u\in\mathcal{W}^{\alpha,p}(0,T;H),$

    (29)

    $\displaystyle(g_{\alpha}*\partial_{t}^{\alpha}u)(t)$
    $\displaystyle=u(t)$
    $\displaystyle\forall\,u\in\mathcal{W}_{z}^{\alpha,p}(0,T;H).$

    ######
    Proof 3.3.

    We start with the proof of the equality Eq. 28.
    Recall that for any function $u\in\mathcal{W}^{\alpha,p}(0,T;H)$ we have $g_{1-\alpha}*u\in AC([0,T];H)$, and the fundamental theorem of calculus for absolutely continuous functions therefore yields, for any $t\in[0,T]$,

    $(g_{1-\alpha}*u)(t)-(g_{1-\alpha}*u)(0)=\int_{0}^{t}\partial_{s}(g_{1-\alpha}*u)(s)\,\textup{ds}=(g_{1}*\partial_{t}^{\alpha}u)(t).$

    We convolve this equality with the kernel $g_{\alpha}$ and make use of the semigroup property Eq. 25 to obtain

    $(g_{1}*u)(t)-(g_{1-\alpha}*u)(0)g_{1+\alpha}(t)=g_{1+\alpha}*\partial_{t}^{\alpha}u,$

    where we have used that $g_{\alpha}*1=g_{\alpha}*g_{1}=g_{1+\alpha}$, because $\alpha\Gamma(\alpha)=\Gamma(1+\alpha)$.
    Next, we differentiate this equality in $t$ and observe that $\partial_{t}(g_{1}*u)=u$, $\partial_{t}g_{1+\alpha}=g_{\alpha}$, and $\partial_{t}(g_{1+\alpha}*v)=g_{\alpha}*v$, which yields Eq. 28.
    We finally note that Eq. 29 follows trivially from Eq. 28 and Eq. 27.

    The following result is a direct consequence of the interaction between fractional derivatives and kernel functions.

    ######
    Corollary 3.4.

    The following identities hold:

    (30)

    $\displaystyle\partial_{t}^{\alpha}(g_{\alpha}*u)$
    $\displaystyle=\partial_{t}(g_{1-\alpha}*g_{\alpha}*u)=\partial_{t}(1*u)=u$

    $\displaystyle\forall\,u\in L^{1}(0,T;H),$

    $\displaystyle\partial_{t}^{1-\alpha}\partial_{t}^{\alpha}u$
    $\displaystyle=\partial_{t}(g_{\alpha}*\partial_{t}^{\alpha}u)=\partial_{t}u$

    $\displaystyle\forall\,u\in W_{0}^{1,1}(0,T;H).$

    We shall require the following special case of the classical Aubin–Lions lemma; see [65]. Suppose that the Hilbert spaces $V,H,Z$ form a Gelfand triple $V\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}H\hookrightarrow Z$. then, the following classical compact embeddings hold:

    (31)

    $\displaystyle W^{1,1}(0,T;Z)\cap L^{p}(0,T;V)$
    $\displaystyle\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}L^{p}(0,T;H),$

    $\displaystyle p\in[1,\infty),$

    $\displaystyle W^{1,r}(0,T;Z)\cap L^{\infty}(0,T;V)$
    $\displaystyle\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}C([0,T];H),$

    $\displaystyle r\in(1,\infty);$

    see [65]. Several fractional counterparts of the Aubin–Lions lemma have been proposed; see [59, 70, 43]. We make use of the following result; see [59, Corollary 3.2]:

    $\displaystyle\mathcal{W}^{\alpha,1}(0,T;Z)\cap L^{p}(0,T;V)$
    $\displaystyle\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}L^{r}(0,T;H),$

    $\displaystyle p\in(1,\infty),\quad r\in[1,p),\quad\alpha\in(0,1).$

    The proof can be easily adapted to the limit case $r=p$ if the $\alpha$-th fractional derivative is in a better space than $L^{1}(0,T;Z)$. This is done for Caputo derivatives in [43]. In fact, we obtain

    (32)

    $\displaystyle\mathcal{W}^{\alpha,r}(0,T;Z)\cap L^{p}(0,T;V)$
    $\displaystyle\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}L^{p}(0,T;H),$

    $\displaystyle r\in(1,\infty),\quad\alpha\in(0,1).$

    The classical chain rule does not hold for fractional derivatives, but one can use the following inequality as a remedy; see [68, Theorem 2.1]:

    (33)

    $\frac{1}{2}\partial_{t}^{\alpha}\|u\|^{2}_{H}+\frac{1}{2}g_{1-\alpha}(t)\|u\|_{H}^{2}\leq(u,\partial_{t}^{\alpha}u)_{H}\quad\forall\,u\in\mathcal{W}_{z}^{\alpha,2}(0,T;H),$

    for $z\in H$ and almost all $t\in(0,T)$.

    ##
    4 Model revisited

    Having summarised the required results from fractional calculus, we revisit the mathematical model that we have derived in Section 2.
    Let us assume for the moment that the solution $\psi$ to the Fokker–Planck equation belongs to $\mathcal{W}^{1-\alpha,p}(0,T;H)\cap C([0,T];H)$ for some $\alpha\in(0,1)$ and a suitable Hilbert space $H$, to be chosen. As $\psi\in C([0,T];H)$, it follows that $\|(g_{\alpha}\ast\psi)(t)\|_{H}\leq\frac{t^{\alpha}}{\Gamma(1+\alpha)}\|\psi\|_{C([0,T];H)}$, and therefore $(g_{\alpha}\ast\psi)(0)=0$. Hence, $\psi\in\mathcal{W}^{1-\alpha,p}_{0}(0,T;H)$. It then follows from (28), with $\alpha$ replaced by $1-\alpha$ and $u=\psi$ there, that $(g_{1-\alpha}\ast\partial_{t}^{1-\alpha}\psi)(t)=\psi(t)$ for $t\in(0,T)$.
    Motivated by these properties, we introduce the auxiliary function $\phi$ by

    (34)

    $\phi:=\partial_{t}^{1-\alpha}\psi=\partial_{t}(g_{\alpha}*\psi),$

    whereby $\psi=g_{1-\alpha}*\phi$. We then have that $\partial_{t}\psi=\partial_{t}(g_{1-\alpha}\ast\phi)=\partial_{t}^{\alpha}\phi$; and, thanks to the assumed continuity of $\psi$ (i.e. $\psi\in C([0,T];H)$) it
    makes sense to require attainment of the initial condition $\psi(0)=\psi^{0}$, i.e. $(g_{1-\alpha}\ast\phi)(0)=\psi^{0}$. We shall therefore introduce the substitution $\phi:=\partial_{t}^{1-\alpha}\psi$ in (24), which results in the following system of equations:

    (35)

    $\displaystyle\partial_{t}u+(u\cdot\nabla_{\!x})u-\nu\Delta_{x}u+\nabla_{\!x}p-\textup{div}_{x}\tau(g_{1-\alpha}*\phi)$
    $\displaystyle=0,$

    $\displaystyle\textup{div}_{x}u$
    $\displaystyle=0,$

    $\displaystyle\partial_{t}^{\alpha}\phi+(u\cdot\nabla_{\!x})\phi+\textup{div}_{q}(\omega(u)q\phi)-\tfrac{1}{2\lambda}\textup{div}_{q}(\nabla_{\!q}\phi+U^{\prime}q\phi)-\varepsilon\Delta_{x}\phi$
    $\displaystyle=0,$

    subject to the initial conditions $u(0)=u^{0}$ and $(g_{1-\alpha}*\phi)(0)=\psi^{0}$ for a given nonnegative $\psi^{0}$ that fulfils $\int_{D}\psi^{0}\,\text{d}q=1$. Furthermore, we equip the system with the following boundary conditions:

    (36)

    $\displaystyle\begin{aligned} u&=0\qquad\text{on }\partial\Omega\times(0,T),\\
    \left(\tfrac{1}{2\lambda}(\nabla_{\!q}\phi+U^{\prime}q\phi)-\omega(u)q\phi\right)\cdot n_{\partial D}&=0\qquad\text{on }\Omega\times\partial D\times(0,T),\\
    \varepsilon\nabla_{\!x}\phi\cdot n_{\partial\Omega}&=0\qquad\text{on }\partial\Omega\times D\times(0,T).\end{aligned}$

    ###
    4.1 The Maxwellian and Maxwellian-weighted function spaces

    We introduce the normalized Maxwellian by

    (37)

    $M(q)=\frac{e^{-U(\tfrac{1}{2}|q|^{2})}}{\int_{D}e^{-U(\tfrac{1}{2}|s|^{2})}\mathop{}\!\mathrm{d}s}.$

    Moreover, we define the (Maxwellian-weighted) Hilbert spaces

    $\displaystyle\mathcal{H}=\{h\in L^{2}(\Omega;\mathbb{R}^{d}):\textup{div}\,h=0\},\qquad~{}\mathcal{H}_{0}$

    $\displaystyle=\{h\in\mathcal{H}:h\cdot n_{\partial\Omega}=0\text{ on }\partial\Omega\},$

    $\displaystyle\mathcal{V}=\{v\in H^{1}(\Omega;\mathbb{R}^{d}):\textup{div}\,v=0\},\qquad~{}\mathcal{V}_{0}$

    $\displaystyle=\{v\in\mathcal{V}:v|_{\partial\Omega}=0\text{ on }\partial\Omega\},$

    $\displaystyle\mathcal{Y}=L^{2}(\Omega\times D),\qquad\quad\,\,\,\widehat{\mathcal{Y}}=L^{2}_{M}(\Omega\times D)$

    $\displaystyle=\{y\in\mathcal{Y}:\|y\|_{\widehat{\mathcal{Y}}}:=\|M^{1/2}y\|_{\mathcal{Y}}<\infty\},$

    $\displaystyle\mathcal{X}=H^{1}(\Omega\times D),\qquad\quad\widehat{\mathcal{X}}=H_{M}^{1}(\Omega\times D)$

    $\displaystyle=\{\phi\in\mathcal{X}:\|\phi\|_{\widehat{\mathcal{X}}}<\infty\},$

    $\displaystyle\mathcal{Z}=H^{1}(D;H^{1}(\Omega)),~{}\widehat{\mathcal{Z}}=H_{M}^{1}(D;H^{1}(\Omega))$

    $\displaystyle=\{\zeta\in\mathcal{Z}:\|\zeta\|_{\widehat{\mathcal{Z}}}<\infty\},$

    where the norms on ${\widehat{\mathcal{X}}}$ and ${\widehat{\mathcal{Z}}}$ are defined by $\|\phi\|_{\widehat{\mathcal{X}}}^{2}:=\|\phi\|_{\widehat{\mathcal{Y}}}^{2}+\|\nabla_{\!q}\phi\|_{\widehat{\mathcal{Y}}}^{2}+\|\nabla_{\!x}\phi\|_{\widehat{\mathcal{Y}}}^{2}$ and $\|\zeta\|_{\widehat{\mathcal{Z}}}^{2}:=\|\zeta\|_{\widehat{\mathcal{X}}}^{2}+\|\nabla_{\!x}\nabla_{\!q}\zeta\|_{\widehat{\mathcal{Y}}}^{2}$. Obviously, $H_{M}^{2}(\Omega\times D)\subseteq{\widehat{\mathcal{Z}}}$, where $H^{2}_{M}(\Omega\times D)$ is the subspace of $H^{1}_{M}(\Omega\times D)$ consisting of all functions defined on $\Omega\times D$ whose second (weak) partial derivatives belong to ${\widehat{\mathcal{Y}}}=L^{2}_{M}(\Omega\times D)$.
    We refer to [6] regarding theoretical results on these weighted Hilbert spaces. In particular, we have the Gelfand triples

    $\displaystyle\mathcal{V}\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}\mathcal{H}\hookrightarrow\mathcal{V}^{\prime},\quad$

    $\displaystyle\mathcal{V}_{0}\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}\mathcal{H}_{0}\hookrightarrow\mathcal{V}_{0}^{\prime},$

    $\displaystyle\mathcal{X}\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}\mathcal{Y}\hookrightarrow\mathcal{X}^{\prime},\quad$

    $\displaystyle{\widehat{\mathcal{X}}}\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}{\widehat{\mathcal{Y}}}\hookrightarrow{\widehat{\mathcal{X}}}^{\prime},$

    where $\mathcal{V}^{\prime}$, $\mathcal{V}^{\prime}_{0}$, $\mathcal{X}^{\prime}$ and ${\widehat{\mathcal{X}}}^{\prime}$ denote the dual space of, respectively,
    $\mathcal{V}$, $\mathcal{V}_{0}$, $\mathcal{X}$ and ${\widehat{\mathcal{X}}}$.

    Using the definition of the normalized Maxwellian $M$, see Eq. 37, we have that

    $M(q)\nabla_{\!q}M(q)^{-1}=-M(q)^{-1}\nabla_{\!q}M(q)=\nabla_{\!q}U(\tfrac{1}{2}|q|^{2})=U^{\prime}(\tfrac{1}{2}|q|^{2})q.$

    We introduce the scaled variable ${\widehat{\phi}}=\phi/M$ and with the formula

    $M\nabla_{\!q}{\widehat{\phi}}=\nabla_{\!q}\phi+M\nabla_{\!q}M^{-1}\phi=\nabla_{\!q}\phi+U^{\prime}q\phi$

    we can rewrite the fractional Fokker–Planck equation in Eq. 35 as

    $\partial_{t}^{\alpha}\phi+(u\cdot\nabla_{\!x})\phi+\textup{div}_{q}\big{(}\omega(u)q\phi\big{)}=\tfrac{1}{2\lambda}\textup{div}_{q}(M\nabla_{\!q}\widehat{\phi})+\varepsilon\Delta_{x}\phi.$

    As was indicated earlier, we shall confine ourselves here to considering the corotational model, i.e.,
    $\omega(v)=-\omega(v)^{\mathrm{T}}$, $q^{\mathrm{T}}\omega(v)q=0$; if $\textup{div}\,v=0$ it then follows that

    (38)

    $\big{(}M{\widehat{\phi}}\,\omega(v)q,\nabla_{\!q}{\widehat{\phi}}\big{)}_{\mathcal{Y}}=\frac{1}{2}\big{(}M\omega(v)q,\nabla_{\!q}{\widehat{\phi}}^{2}\big{)}_{\mathcal{Y}}=-\frac{1}{2}\big{(}\textup{div}_{q}(M\omega(v)q),{\widehat{\phi}}^{2}\big{)}_{\mathcal{Y}}=0;$

    see [9, 6].
    We note in passing that partial integration yields the following equalities:

    (39)

    $\displaystyle-2\big{(}M\omega(u)q\widehat{\varphi},\nabla_{\!q}{\widehat{\phi}}\big{)}_{\mathcal{Y}}$
    $\displaystyle=\big{(}\nabla_{\!x}(M\widehat{\varphi}\nabla_{\!q}{\widehat{\phi}})q,u\big{)}_{\mathcal{Y}}+\big{(}u\cdot q,\textup{div}_{x}(M\widehat{\varphi}\nabla_{\!q}{\widehat{\phi}})\big{)}_{\mathcal{Y}}$

    $\displaystyle=\big{(}M\nabla_{\!x}\widehat{\varphi}(\nabla_{\!q}{\widehat{\phi}})^{\mathrm{T}}q,u\big{)}_{\mathcal{Y}}+\big{(}M\widehat{\varphi}\nabla_{\!x}\nabla_{\!q}{\widehat{\phi}}\,q,u\big{)}_{\mathcal{Y}}$

    $\displaystyle\quad+\big{(}u\cdot q,M\nabla_{\!x}\widehat{\varphi}\cdot\nabla_{\!q}{\widehat{\phi}}\big{)}_{\mathcal{Y}}+\big{(}u\cdot q,M\widehat{\varphi}\,\textup{div}_{x}\nabla_{\!q}{\widehat{\phi}}\big{)}_{\mathcal{Y}}.$

    We recall that the stress tensor $\tau(\psi)=\tau_{1}(\psi)+\tau_{2}(\psi)$ is of the form given by (22) and (23); i.e.,

    (40)

    $\displaystyle\tau^{1}(\psi)=\gamma\,\mathcal{C}(\psi),\quad\tau^{2}(\psi)=\gamma\int_{D}\psi\mathop{}\!\mathrm{d}q\ I_{3},\quad\mathcal{C}(\psi):=\int_{D}F(q)q^{\mathrm{T}}\psi\mathop{}\!\mathrm{d}q,$

    where $\gamma>0$ is a dimensionless constant.

    For $\mathcal{C}(M{\widehat{\psi}})$, we are in a setting that allows us to deduce the following bound; see also [9, Eq. (3.7)]:

    (41)

    $\displaystyle\int_{\Omega}|\mathcal{C}(M{\widehat{\psi}})|^{2}\mathop{}\!\mathrm{d}x$
    $\displaystyle=\int_{\Omega}\left|\,\int_{D}F(q)q^{\mathrm{T}}M{\widehat{\psi}}\mathop{}\!\mathrm{d}q\,\right|^{2}\mathop{}\!\mathrm{d}x$

    $\displaystyle\leq\int_{D}M|F(q)q^{\mathrm{T}}|^{2}\,\text{d}q\,\int_{\Omega\times D}M|{\widehat{\psi}}|^{2}\mathop{}\!\mathrm{d}(x,q)$

    $\displaystyle\leq C\|{\widehat{\psi}}\|_{{\widehat{\mathcal{Y}}}}^{2}\quad\forall\,{\widehat{\psi}}\in{\widehat{\mathcal{Y}}}.$

    ##
    5 Existence of weak solutions

    In this section, we prove the existence of a weak solution to the time-fractional Navier–Stokes–Fokker–Planck system. We use a Galerkin procedure and discretize the partial differential equations in space and derive suitable energy estimates. We will emphasize the places where the time-fractional derivative comes into play. We shall then pass to the limit in the sequence of Galerkin approximations to deduce the existence of a weak solution. We shall proceed step by step and prove this result through several lemmas. While, for physical reasons, the derivation of the system was in the previous sections discussed in the case of $d=3$ space dimensions, the analysis below applies to both two and three space dimensions. We begin by introducing the concept of a weak solution to the time-fractional corotational Navier–Stokes–Fokker–Planck system under consideration.

    ######
    Definition 5.1.

    Suppose that $d\in\{2,3\}$.
    We call the pair $(u,{\widehat{\phi}})$ a weak solution to the system Eq. 35, Eq. 36 provided that

    $\displaystyle u$
    $\displaystyle\in L^{\infty}(0,T;\mathcal{H}_{0})\cap L^{2}(0,T;\mathcal{V}_{0})\cap W^{1,\tfrac{8}{4+d}}(0,T;\mathcal{V}_{0}^{\prime}),$

    $\displaystyle{\widehat{\phi}}$
    $\displaystyle\in L^{2}(0,T;{\widehat{\mathcal{X}}}),\quad\partial_{t}^{\alpha}{\widehat{\phi}}\in L^{\tfrac{8}{4+d}}\big{(}0,T;{\widehat{\mathcal{Z}}}^{\prime}\big{)},$

    satisfies the initial conditions $u(0)=u^{0}$, $(g_{1-\alpha}*{\widehat{\phi}})(0)={\widehat{\psi}}^{0}:=\psi^{0}/M$ and the variational problems

    (42)

    $\displaystyle\langle\partial_{t}u,v\rangle_{L^{8/(4-d)}\mathcal{V}}+\big{(}(u\cdot\nabla_{\!x})u,v\big{)}_{L^{2}\mathcal{H}}+\nu(\nabla_{\!x}u,\nabla_{\!x}v)_{L^{2}\mathcal{H}}$

    $\displaystyle+k_{B}\mu_{T}\big{(}\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}),\nabla_{\!x}v\big{)}_{L^{2}\mathcal{H}}=0,$

    (43)

    $\displaystyle\langle\partial_{t}^{\alpha}{\widehat{\phi}},\widehat{\zeta}\rangle_{L^{8/(4-d)}{\widehat{\mathcal{Z}}}}-(u{\widehat{\phi}},\nabla_{\!x}\widehat{\zeta})_{L^{2}{\widehat{\mathcal{Y}}}}+\frac{1}{2\lambda}(\nabla_{\!q}{\widehat{\phi}},\nabla_{\!q}\widehat{\zeta})_{L^{2}{\widehat{\mathcal{Y}}}}+\varepsilon(\nabla_{\!x}{\widehat{\phi}},\nabla_{\!x}\widehat{\zeta})_{L^{2}{\widehat{\mathcal{Y}}}}$

    $\displaystyle+\frac{1}{2}(\nabla_{\!x}({\widehat{\phi}}\nabla_{\!q}\widehat{\zeta})q,u)_{L^{2}{\widehat{\mathcal{Y}}}}-\frac{1}{2}\big{(}u\cdot q,\textup{div}_{x}({\widehat{\phi}}\nabla_{\!q}\widehat{\zeta})\big{)}_{L^{2}{\widehat{\mathcal{Y}}}}=0.$

    for all $v\in L^{8/(4-d)}(0,T;\mathcal{V}_{0})$ and $\widehat{\zeta}\in L^{8/(4-d)}(0,T;{\widehat{\mathcal{Z}}})$. In these variational problems and hereafter, for a Hilbert space $H$ and $p\in[1,\infty]$, subscripts of the form $L^{p}H$ and $L^{p}_{t}H$ appearing in the various inner products, norms, and duality pairings, signify $L^{p}(0,T;H)$ and $L^{p}(0,t;H)$, respectively.

    We summarize the assumptions that we require for proving the existence of a weak solution in the sense of Definition 5.1 in Assumption 5 below.

    Let the following assumptions hold:

      *
    •

    $D\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, is a bounded open ball centered at the origin, $\Omega\subseteq\mathbb{R}^{d}$ is a Lipschitz domain (i.e., bounded, open, connected set in $\mathbb{R}^{d}$, with a Lipschitz-continuous boundary $\partial\Omega$), and $T<\infty$ is a fixed final time;

      *
    •

    $u^{0}\in\mathcal{H}_{0}$, ${\widehat{\psi}}^{0}\in{\widehat{\mathcal{X}}}$ with ${\widehat{\psi}}^{0}\in H^{1}_{M}\big{(}D;H^{1+d/2+\delta}(\Omega)\big{)}$ for $\delta>0$ arbitrarily small;

      *
    •

    $\tau(\psi)=\tau_{1}(\psi)+\tau_{2}(\psi)$ is given by (40), with the identity matrix $I_{3}\in\mathbb{R}^{3\times 3}$ replaced by the identity matrix $I_{d}\in\mathbb{R}^{d\times d}$ in the definition of $\tau_{2}(\psi)$, and $\mathcal{C}$ satisfies Eq. 41;

      *
    •

    $\alpha\in(1/2,1)$;

      *
    •

    $k_{B},\mu_{T},\nu,\lambda,\varepsilon>0$.

    The main result of this paper is the following theorem asserting the existence of global-in-time large-data weak solutions to the time-fractional Navier–Stokes–Fokker–Planck system under consideration.

    ######
    Theorem 5.2.

    Let Section 5 hold;
    then, there exists a weak solution $(u,{\widehat{\phi}})$ to the system Eq. 35, Eq. 36 in the sense of Definition 5.1.

    In order to prove this theorem, we state several lemmas, which will eventually imply Theorem 5.2. We begin by constructing a sequence of Galerkin approximations $\{(u_{k},{\widehat{\phi}}_{k})\}_{k=1}^{\infty}$ to the system of partial differential equations under consideration, resulting in a system of fractional-order ordinary differential equations, which admits a local-in-time solution $(u_{k},{\widehat{\phi}}_{k})$
    for each $k\geq 1$ thanks to standard theory. We then derive an energy estimate for the sequence of Galerkin approximations, which is uniform with respect to $k$; this then implies that, for each $k\geq 1$, the local-in-time solution of the fractional-order system of ordinary differential equations can be extended to the entire time-interval $[0,T]$; it also implies the existence of a weakly/weakly-$*$ convergent subsequence $(u_{k_{j}},{\widehat{\phi}}_{k_{j}})$. Finally, we pass to the limit $j\to\infty$ and apply a compactness argument to deduce that the limiting pair of functions
    $(u,{\widehat{\phi}})$ is in fact a weak solution to the system of partial differential equations in the sense of Definition 5.1. The Galerkin method has been applied to various time-fractional PDEs; see, e.g., [27, 26, 25, 69]; it has also been applied to Navier–Stokes–Fokker–Planck systems in [34, 9, 13, 14], with an integer-order Fokker–Planck equation.

    ###
    5.1 Galerkin discretization

    We follow the construction of [17, Section 2.1] and conclude by the Hilbert–Schmidt theorem [17, Lemma A.4] the existence of a countable set $\{h_{j}\}_{j=1}^{\infty}$ of eigenfunctions in $\mathcal{V}_{0}\cap H^{1+\frac{d}{2}+\delta}(\Omega)^{d}$, with $\delta>0$ arbitrarily small, whose linear span is dense in $\mathcal{H}_{0}$ such that the $h_{j}$, $j\in\{1,2,\dots\}$, are orthonormal in $\mathcal{H}$ and orthogonal in $H^{1+\frac{d}{2}+\delta}(\Omega)^{d}$ in the sense that $(h_{j},h_{i})_{H^{1+\frac{d}{2}+\delta}(\Omega)}=\lambda_{j}\delta_{i,j}$ for any $i,j\in\{1,2,\ldots\}$ and $\lambda_{j}>0$ for all $j=1,2,\ldots$. Similarly, we fix a countable set $\{y_{j}\}_{j=1}^{\infty}$ in $H^{2}_{M}(\Omega\times D)$ that forms an orthogonal system in ${\widehat{\mathcal{X}}}$ and an orthonormal system in ${\widehat{\mathcal{Y}}}$.
    We then define the $k$-dimensional linear spaces

    $\displaystyle\mathcal{H}_{k}:=\text{span}\{h_{1},\dots,h_{k}\},\quad\widehat{\mathcal{Y}}_{k}:=\text{span}\{y_{1},\dots,y_{k}\},$

    and we consider the Galerkin approximations

    (44)

    $\begin{gathered}u_{k}(t)=\sum_{j=1}^{k}u_{k}^{j}(t)h_{j},\quad{\widehat{\phi}}_{k}(t)=\sum_{j=1}^{k}{\widehat{\phi}}^{j}_{k}(t)y_{j},\end{gathered}$

    where
    $u^{j}_{k}$ and ${\widehat{\phi}}^{j}_{k}$ are real-valued coefficient functions for all $j\in\{1,\dots,k\}$.
    The canonical orthogonal projection onto the finite-dimensional space $\mathcal{H}_{k}$ is defined by $\Pi_{\mathcal{H}_{k}}:\mathcal{H}\to\mathcal{H}_{k}$, $h\mapsto\sum_{j=1}^{k}(h,h_{j})_{\mathcal{H}}h_{j}$, and in the same way for $\Pi_{{\widehat{\mathcal{Y}}}_{k}}:{\widehat{\mathcal{Y}}}\to{\widehat{\mathcal{Y}}}_{k}$.
    For $h=\sum_{j=1}^{\infty}(h,h_{j})_{\mathcal{H}}h_{j}$ we have that

    $\|h\|^{2}_{H^{1+\frac{d}{2}+\delta}(\Omega)}=\sum_{j=1}^{\infty}\lambda_{j}|(h,h_{j})_{\mathcal{H}}|^{2},\quad\|\Pi_{\mathcal{H}_{k}}h\|^{2}_{H^{1+\frac{d}{2}+\delta}(\Omega)}=\sum_{j=1}^{k}\lambda_{j}|(h,h_{j})_{\mathcal{H}}|^{2},$

    from which we conclude via the Sobolev embedding theorem that, for each $k\geq 1$,

    $\|\Pi_{\mathcal{H}_{k}}h\|_{W^{1,\infty}(\Omega)}\leq C\|\Pi_{\mathcal{H}_{k}}h\|_{H^{1+\frac{d}{2}+\delta}(\Omega)}\leq C\|h\|_{H^{1+\frac{d}{2}+\delta}(\Omega)}.$

    The Galerkin equations read as follows: we wish to find a tuple $(u_{k},{\widehat{\phi}}_{k})\in\mathcal{H}_{k}\times{\widehat{\mathcal{Y}}}_{k}$ for each $k\geq 1$ such that $u_{k}(0)=u_{k}^{0}:=\Pi_{\mathcal{H}_{k}}u^{0}$, $(g_{1-\alpha}*{\widehat{\phi}}_{k})(0)={\widehat{\psi}}_{k}^{0}:=\Pi_{{\widehat{\mathcal{Y}}}_{k}}{\widehat{\psi}}^{0}$, and

    (45)

    $\displaystyle(\partial_{t}u_{k},v_{k})_{\mathcal{H}}$
    $\displaystyle+((u_{k}\cdot\nabla_{\!x})u_{k},v_{k})_{\mathcal{H}}+\nu(\nabla_{\!x}u_{k},\nabla_{\!x}v_{k})_{\mathcal{H}}$

    $\displaystyle+\,k_{B}\mu_{T}(\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k}),\nabla_{\!x}v_{k})_{\mathcal{H}}=0,$

    (46)

    $\displaystyle(\partial_{t}(g_{1-\alpha}*{\widehat{\phi}}_{k}),\widehat{\zeta}_{k})_{\widehat{\mathcal{Y}}}$
    $\displaystyle+((u_{k}\cdot\nabla_{\!x}){\widehat{\phi}}_{k},\widehat{\zeta}_{k})_{\widehat{\mathcal{Y}}}+\frac{1}{2\lambda}(\nabla_{\!q}{\widehat{\phi}}_{k},\nabla_{\!q}\widehat{\zeta}_{k})_{\widehat{\mathcal{Y}}}$

    $\displaystyle+\varepsilon(\nabla_{\!x}{\widehat{\phi}}_{k},\nabla_{x}\widehat{\zeta}_{k})_{\widehat{\mathcal{Y}}}-(\omega(u_{k})q{\widehat{\phi}}_{k},\nabla_{\!q}\widehat{\zeta}_{k})_{\widehat{\mathcal{Y}}}=0,$

    for all $v_{k}\in\mathcal{H}_{k}$ and $\widehat{\zeta}_{k}\in{\widehat{\mathcal{Y}}}_{k}$.

    ######
    Lemma 5.3.

    Suppose that Section 5 holds; then, for each $k\geq 1$,
    there exists a local-in-time solution $(u_{k},{\widehat{\phi}}_{k})$ to the Galerkin system Eq. 45, Eq. 46.

    ######
    Proof 5.4.

    Let $U_{k}(t):=(u_{k}^{1}(t),\ldots,u_{k}^{k}(t))^{\mathrm{T}}$ and $\widehat{\Phi}_{k}(t)=(\widehat{\phi}_{1}^{k},\ldots,\widehat{\phi}_{k}^{k})^{\rm T}$. With this notation the Galerkin subsystem (45) becomes an initial-value problem for a system of ordinary differential equations of the form $\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}\mathrm{t}}U_{k}=F(t,U_{k},\widehat{\Phi}_{k})$, while, by noting that $\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}\mathrm{t}}(g_{1-\alpha}*{\widehat{\phi}}_{k})=\partial_{t}^{\alpha}{\widehat{\phi}}_{k}$, the Galerkin subsystem (46) takes the form of an initial-value problem for a system of fractional-order ordinary differential equations $(\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}\mathrm{t}})^{\alpha}\widehat{\Phi}_{k}=G(t,U_{k},\widehat{\Phi}_{k})$. As the functions $F$ and $G$ are continuous with respect to their arguments and locally Lipschitz continuous with respect to their second and third arguments, we can appeal to the generalization of the Cauchy–Lipschitz theorem stated in Theorem 5.1 of [22] to deduce the existence of a unique continuous solution, defined on a time interval $[0,T_{k}]$ where $0<T_{k}\leq T$, where $U_{k}$ is, in fact, a continuously differentiable function of $t$ by the classical Cauchy–Lipschitz theorem.

    ####
    5.1.1 Energy estimates

    Next, we derive a $k$-uniform energy estimate, which will allow us to extend, for each $k\geq 1$, the corresponding local-in-time Galerkin solution, whose existence is guaranteed by Lemma 5.3, to the entire time interval $[0,T]$; it will also enable us to extract weakly converging subsequences of Galerkin approximations. We begin by deriving a bound on the solution to the Galerkin approximation of the Navier–Stokes equation; we shall then derive a bound on the solution to the Galerkin approximation of the Fokker–Planck equation. At the end, we will add the two bounds and apply Gronwall’s lemma to obtain a $k$-uniform energy estimate.

    ######
    Lemma 5.5.

    Let Section 5 hold; then the following bound on the Galerkin solution $u_{k}$, in terms of $u_{k}^{0}$ and ${\widehat{\phi}}_{k}$, holds for all $t\in(0,T_{k})$:

    (47)

    $\frac{1}{2}\|u_{k}(t)\|_{\mathcal{H}}^{2}+\frac{\nu}{2}\|\nabla_{\!x}u_{k}\|_{L^{2}_{t}{\mathcal{H}}}^{2}\leq\frac{1}{2}\|u^{0}_{h}\|_{\mathcal{H}}^{2}+\frac{T^{2-2\alpha}k_{B}^{2}\mu_{T}^{2}}{2\nu(1-\alpha)^{2}}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}.$

    ######
    Proof 5.6.

    We take the test function $v_{k}=u_{k}(t)$ in the equation Eq. 42, which gives

    $\displaystyle\frac{1}{2}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t}\|u_{k}\|_{\mathcal{H}}^{2}+\nu\|\nabla_{\!x}u_{k}\|_{\mathcal{H}}^{2}$
    $\displaystyle=-k_{B}\mu_{T}\big{(}\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}),\nabla_{\!x}u_{k}\big{)}_{\mathcal{H}},$

    and we can further bound the right-hand side from above to deduce that

    $\frac{1}{2}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t}\|u_{k}\|_{\mathcal{H}}^{2}+\nu\|\nabla_{\!x}u_{k}\|_{\mathcal{H}}^{2}\leq\frac{\nu}{2}\|\nabla_{\!x}u_{k}\|_{\mathcal{H}}^{2}+\frac{k_{B}^{2}\mu_{T}^{2}}{2\nu}\|\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k})\|_{\mathcal{H}}^{2}.$

    By bounding the term $\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k})$ as in Eq. 41 we arrive at the inequality

    (48)

    $\displaystyle\frac{1}{2}\frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t}\|u_{k}\|_{\mathcal{H}}^{2}+\frac{\nu}{2}\|\nabla_{\!x}u_{k}\|_{\mathcal{H}}^{2}$
    $\displaystyle\leq\frac{k_{B}^{2}\mu_{T}^{2}}{2\nu}\|g_{1-\alpha}*{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}.$

    Next, we note that, for any $t\in(0,T_{k})$,

    (49)

    $\displaystyle\int_{0}^{t}\|(g_{1-\alpha}*{\widehat{\phi}}_{k})(s)\|_{\widehat{\mathcal{Y}}}^{2}\,\textup{ds}\leq\int_{0}^{t}\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}\big{)}^{2}(s)\,\textup{ds}\leq\|g_{1-\alpha}\|_{L^{1}_{t}}^{2}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}.$

    This follows by observing that $g_{1-\alpha}$ only depends on the scalar variable $s$, which permits pulling the ${\widehat{\mathcal{Y}}}$-norm inside of the convolution, followed by applying Young’s convolution inequality (cf. Lemma 3.2 in [57]) in the resulting integrand.

    We then integrate (48) over the interval $[0,t]$ where $t\in(0,T_{k})$ and use (49) to bound the right-hand side of the resulting inequality. Finally we
    note that $g_{1-\alpha}$ is integrable on $(0,t)$ and its integral is bounded by $T^{1-\alpha}/(1-\alpha)$.
    This gives Eq. 47.

    Having derived a bound on $u_{k}$, we move on to the derivation of a bound on ${\widehat{\phi}}^{k}$ by testing the Galerkin system (46).

    ######
    Lemma 5.7.

    Let Section 5 hold and let $\gamma>0$ be arbitrary but fixed; then, the following bound on the sequence of Galerkin solutions $\{(u_{k},{\widehat{\phi}}_{k})\}_{k=1}^{\infty}$ holds:

    (50)

    $\displaystyle\frac{\gamma}{2}\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)+\frac{\gamma\,T^{-\alpha}}{16\,\Gamma(1-\alpha)}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}+\frac{\gamma}{2\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}+\gamma\varepsilon\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}$

    $\displaystyle\quad\leq C(\alpha,\gamma)\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{H^{1}(D;W^{1,\infty}(\Omega))}^{2}\int_{0}^{t}g_{2\alpha-1}(s)\|u_{k}(s)\|_{\mathcal{H}}^{2}\,\textup{ds}+C(\alpha,\gamma,T)\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{X}}}^{2}.$

    ######
    Proof 5.8.

    We note again that, thanks to the inverse convolution property, see Eq. 29, ${\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}=g_{\alpha}*\partial_{t}^{\alpha}{\widehat{\phi}}_{k}$.
    We take this function as the test function in the variational Fokker–Planck equation Eq. 43, i.e., $\widehat{\zeta}={\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}=g_{\alpha}*\partial_{t}^{\alpha}{\widehat{\phi}}_{k}$, which gives

    (51)

    $\displaystyle(\partial_{t}^{\alpha}{\widehat{\phi}}_{k},g_{\alpha}*\partial_{t}^{\alpha}{\widehat{\phi}}_{k})_{\widehat{\mathcal{Y}}}+\frac{1}{2\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\varepsilon\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle\quad=g_{\alpha}(t)\cdot\Big{(}\big{(}(u_{k}\cdot\nabla_{\!x}){\widehat{\phi}}_{k},{\widehat{\psi}}_{k}^{0}\big{)}_{\widehat{\mathcal{Y}}}+\frac{1}{2\lambda}\big{(}\nabla_{\!q}{\widehat{\phi}}_{k},\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\big{)}_{\widehat{\mathcal{Y}}}$

    $\displaystyle\qquad+\varepsilon\big{(}\nabla_{\!x}{\widehat{\phi}}_{k},\nabla{\widehat{\psi}}_{k}^{0}\big{)}_{\widehat{\mathcal{Y}}}-\big{(}\omega(u_{k})q{\widehat{\phi}}_{k},\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\big{)}_{\widehat{\mathcal{Y}}}\Big{)}=:R.$

    We then use the fractional chain inequality Eq. 33 to bound the left-hand side of Eq. 51 from below, which yields

    $\frac{1}{2}\partial_{t}^{\alpha}\|{\widehat{\phi}}_{k}-g_{\alpha}{\widehat{\psi}}_{k}^{0}\|^{2}_{\widehat{\mathcal{Y}}}\leq(\partial_{t}^{\alpha}{\widehat{\phi}}_{k},{\widehat{\phi}}_{k}-g_{\alpha}{\widehat{\psi}}_{k}^{0})_{\widehat{\mathcal{Y}}}=(\partial_{t}^{\alpha}{\widehat{\phi}}_{k},g_{\alpha}*\partial_{t}^{\alpha}{\widehat{\phi}}_{k})_{\widehat{\mathcal{Y}}}.$

    Regarding the right-hand side of Eq. 51, we integrate the last term containing $\omega(u_{k})$ by parts, see
    Eq. 39, and get

    $\displaystyle-\big{(}M\omega(u_{k})q{\widehat{\phi}},\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\big{)}_{\mathcal{Y}}$
    $\displaystyle=\big{(}M\nabla_{\!x}{\widehat{\phi}}_{k}(\nabla_{\!q}{\widehat{\psi}}_{k}^{0})^{\mathrm{T}}q,u_{k}\big{)}_{\mathcal{Y}}+\big{(}M{\widehat{\phi}}_{k}\nabla_{\!x}\nabla_{\!q}{\widehat{\psi}}_{k}^{0}q,u_{k}\big{)}_{\mathcal{Y}}$

    $\displaystyle\quad+\big{(}u_{k}\cdot q,M\nabla_{\!x}{\widehat{\phi}}_{k}\cdot\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\big{)}_{\mathcal{Y}}+\big{(}u_{k}\cdot q,M{\widehat{\phi}}_{k}\textup{div}_{x}\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\big{)}_{\mathcal{Y}}.$

    We apply Hölder’s inequality to obtain the following bound on the right-hand side, $R$, of the equality Eq. 51:

    $\displaystyle R\leq{}$
    $\displaystyle g_{\alpha}(t)\cdot\!\Big{(}\|u_{k}\|_{\mathcal{H}}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{L^{2}(D;L^{\infty}(\Omega))}$

    $\displaystyle+\frac{1}{2\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}\|\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}+\varepsilon\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}\|\nabla_{\!x}{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}$

    $\displaystyle+C\|u_{k}\|_{\mathcal{H}}\|q\|_{L^{\infty}(D)}\|M^{1/2}\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\|_{L^{2}(D;W^{1,\infty}(\Omega))}\big{(}\|{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}+\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}\big{)}\Big{)}.$

    Hence, thanks to Young’s inequality, we arrive at the following bound on $R$:

    $\displaystyle R\leq$
    $\displaystyle~{}{}\frac{1}{\varepsilon}g_{\alpha}(t)^{2}\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{L^{2}(D;W^{1,\infty}(\Omega))}^{2}\|u_{k}\|_{\mathcal{H}}^{2}+\frac{\varepsilon}{4}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle{}+\frac{1}{4\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\frac{g_{\alpha}(t)^{2}}{4\lambda}\|\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}^{2}+\frac{\varepsilon}{4}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\varepsilon g_{\alpha}(t)^{2}\|\nabla_{\!x}{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle{}+\delta\|{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}+\frac{\varepsilon}{4}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+C(\varepsilon,\delta)g_{\alpha}(t)^{2}\|M^{1/2}\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\|_{L^{2}(D;W^{1,\infty}(\Omega))}^{2}\|u_{k}\|_{\mathcal{H}}^{2},$

    where $\delta>0$ is sufficiently small, to be chosen appropriately later on. After combining the lower bound on the left-hand side of Eq. 51 and the upper bound on the right-hand side we have that

    $\displaystyle\frac{1}{2}\partial_{t}^{\alpha}\|{\widehat{\phi}}_{k}-g_{\alpha}{\widehat{\psi}}_{k}^{0}\|^{2}_{\widehat{\mathcal{Y}}}+\frac{1}{2\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\varepsilon\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle\quad\leq\frac{1}{\varepsilon}g_{\alpha}(t)^{2}\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{L^{2}(D;W^{1,\infty}(\Omega))}^{2}\|u_{k}\|_{\mathcal{H}}^{2}+\frac{\varepsilon}{4}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle{}\qquad+\frac{1}{4\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\frac{g_{\alpha}(t)^{2}}{4\lambda}\|\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}^{2}+\frac{\varepsilon}{4}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\varepsilon g_{\alpha}(t)^{2}\|\nabla_{\!x}{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle{}\qquad+\delta\|{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\frac{\varepsilon}{4}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+C(\varepsilon,\delta)g_{\alpha}(t)^{2}\|M^{1/2}\nabla_{\!q}{\widehat{\psi}}_{k}^{0}\|_{L^{2}(D;W^{1,\infty}(\Omega))}^{2}\|u_{k}\|_{\mathcal{H}}^{2},$

    and absorbing terms on the right-hand side into the left-hand side gives

    $\displaystyle\frac{1}{2}\partial_{t}^{\alpha}\|{\widehat{\phi}}_{k}-g_{\alpha}{\widehat{\psi}}_{k}^{0}\|^{2}_{\widehat{\mathcal{Y}}}+\frac{1}{4\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+\frac{\varepsilon}{2}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}$

    $\displaystyle\quad\leq C(\varepsilon,\delta)g_{\alpha}(t)^{2}\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{H^{1}(D;W^{1,\infty}(\Omega))}\|u_{k}\|_{\mathcal{H}}^{2}+\delta\|{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}^{2}+C(\varepsilon,\delta)g_{\alpha}(t)^{2}\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{X}}}^{2}.$

    We note that $g_{\alpha}^{2}=\frac{\Gamma(2\alpha-1)}{\Gamma(\alpha)^{2}}g_{2\alpha-1}$ is integrable for $\alpha\in(\tfrac{1}{2},1)$ and $g_{1}*g_{2\alpha-1}=g_{2\alpha}$, which is continuous, bounded, and monotonically increasing on $[0,T]$ for $\alpha\in(\tfrac{1}{2},1)$. We integrate the inequality over $(0,t)$ and exploit the representation $\partial_{t}^{\alpha}v=\partial_{t}(g_{1-\alpha}*v)$ of the Riemann–Liouville derivative, which gives

    (52)

    $\displaystyle\frac{1}{2}\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)+\frac{1}{4\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}+\frac{\varepsilon}{2}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}-\delta\|{\widehat{\phi}}_{k}\|^{2}_{L^{2}_{t}{\widehat{\mathcal{Y}}}}$

    $\displaystyle\quad\leq C(\delta,\alpha)\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{H^{1}(D;W^{1,\infty}(\Omega))}^{2}\int_{0}^{t}g_{2\alpha-1}(s)\|u_{k}(s)\|_{\mathcal{H}}^{2}\,\textup{ds}+C(\delta)\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{X}}}^{2}g_{2\alpha}(T).$

    Further, we derive a lower bound on the first term of the left-hand side by noting that $(g_{1}*v)(t)\leq T^{\alpha}\Gamma(1-\alpha)(g_{1-\alpha}*v)(t)$, see Eq. 26, and therefore we have that

    $\displaystyle\frac{1}{2}\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)$

    $\displaystyle\quad\geq\frac{T^{-\alpha}}{2\Gamma(1-\alpha)}\int_{0}^{t}\|{\widehat{\phi}}_{k}(s)-{\widehat{\psi}}_{k}^{0}g_{\alpha}(s)\|_{{\widehat{\mathcal{Y}}}}^{2}\,\textup{ds}$

    $\displaystyle\quad\geq\frac{T^{-\alpha}}{2\Gamma(1-\alpha)}\int_{0}^{t}\big{|}\,\|{\widehat{\phi}}_{k}(s)\|_{\widehat{\mathcal{Y}}}-g_{\alpha}(s)\|{\widehat{\psi}}_{k}^{0}\|_{{\widehat{\mathcal{Y}}}}\big{|}^{2}\,\textup{ds}$

    $\displaystyle\quad=\frac{T^{-\alpha}}{2\Gamma(1-\alpha)}\int_{0}^{t}\|{\widehat{\phi}}_{k}(s)\|_{\widehat{\mathcal{Y}}}^{2}-2g_{\alpha}(s)\|{\widehat{\phi}}_{k}(s)\|_{\widehat{\mathcal{Y}}}\|{\widehat{\psi}}_{k}^{0}\|_{{\widehat{\mathcal{Y}}}}+g_{\alpha}^{2}(s)\|{\widehat{\psi}}_{k}^{0}\|_{{\widehat{\mathcal{Y}}}}^{2}\,\textup{ds},$

    where we applied the reverse triangle inequality in the second estimate.
    The function $g_{\alpha}$ belongs to $L^{2}(0,t)$ for any $\alpha\in(\tfrac{1}{2},1)$ and the integral of $g_{\alpha}^{2}$ is positive. We apply Hölder’s inequality in the second term of the integrand and note that the $L^{2}(0,t)$-norm of $g_{\alpha}$ has the upper bound $C(\alpha)T^{\alpha-1/2}$. We thus have that

    $\displaystyle\frac{1}{2}\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)$
    $\displaystyle\geq\frac{T^{-\alpha}}{2\Gamma(1-\alpha)}\left[\frac{1}{2}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}-2C(\alpha)T^{\alpha-1/2}\|{\widehat{\psi}}_{k}^{0}\|_{{\widehat{\mathcal{Y}}}}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}\right]$

    $\displaystyle\geq\frac{T^{-\alpha}}{2\Gamma(1-\alpha)}\left[\frac{1}{4}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}-C(T,\alpha)\|{\widehat{\psi}}_{k}^{0}\|_{{\widehat{\mathcal{Y}}}}^{2}\right],$

    where we have applied Young’s inequality in the last step. We multiply the energy estimate Eq. 52 by $\gamma>0$ and obtain for $\delta=\frac{T^{-\alpha}}{16\Gamma(1-\alpha)}$ the estimate Eq. 50.

    Next, we combine the estimates on $u_{k}$ and ${\widehat{\phi}}_{k}$, see Lemma 5.5 and Lemma 5.7, to obtain a $k$-uniform bound.

    ######
    Lemma 5.9.

    Let Section 5 hold; then, the following $k$-uniform estimate on the Galerkin solution $(u_{k},{\widehat{\phi}}_{k})$ holds:

    (53)

    $\displaystyle\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)+\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{X}}}}^{2}+\|u_{k}(t)\|_{\mathcal{H}}^{2}+\|\nabla u_{k}\|_{L^{2}_{t}{\mathcal{H}}}^{2}$

    $\displaystyle\quad\leq C\big{(}\alpha,T,\|u^{0}\|_{\mathcal{H}}^{2},\|M^{1/2}{\widehat{\psi}}^{0}\|^{2}_{H^{1}(D;H^{1+d/2+\delta}(\Omega))}\big{)}.$

    ######
    Proof 5.10.

    We add the integrated velocity inequality Eq. 47 to the bound Eq. 50 and obtain the following combined bound:

    $\displaystyle\frac{\gamma}{2}\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)+\frac{\gamma\,T^{-\alpha}}{16\,\Gamma(1-\alpha)}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}+\frac{\gamma}{2\lambda}\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}+\gamma\varepsilon\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}$

    $\displaystyle\quad+\frac{1}{2}\|u_{k}(t)\|_{\mathcal{H}}^{2}+\frac{\nu}{2}\|\nabla u_{k}\|_{L^{2}_{t}{\mathcal{H}}}^{2}$

    $\displaystyle\leq C(\alpha,\gamma)\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{H^{1}(D;W^{1,\infty}(\Omega))}^{2}\int_{0}^{t}g_{2\alpha-1}(s)\|u_{k}(s)\|_{\mathcal{H}}^{2}\,\textup{ds}+C(\alpha,\gamma,T)\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{X}}}^{2}$

    $\displaystyle\quad+\frac{1}{2}\|u_{k}^{0}\|_{\mathcal{H}}^{2}+\frac{T^{2-2\alpha}k_{B}^{2}\mu_{T}^{2}}{2\nu(1-\alpha)^{2}}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{Y}}}}^{2}.$

    We now choose $\gamma$ such that

    $\frac{\gamma\,T^{-\alpha}}{16\,\Gamma(1-\alpha)}\geq\frac{T^{2-2\alpha}k_{B}^{2}\mu_{T}^{2}}{\nu(1-\alpha)^{2}}.$

    Hence, we can absorb the last term on the right-hand side into the second term on the left-hand side, and
    the combined energy inequality thus becomes

    (54)

    $\displaystyle\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)+\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{X}}}}^{2}+\|u_{k}(t)\|_{\mathcal{H}}^{2}+\|\nabla u_{k}\|_{L^{2}_{t}{\mathcal{H}}}^{2}$

    $\displaystyle\quad\leq C(\alpha,T)\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{H^{1}(D;W^{1,\infty}(\Omega))}^{2}\int_{0}^{t}g_{2\alpha-1}(s)\|u_{k}(s)\|_{\mathcal{H}}^{2}\,\textup{ds}$

    $\displaystyle\qquad+C(\alpha,T)\big{(}\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{X}}}^{2}+\|u_{k}^{0}\|_{\mathcal{H}}^{2}\big{)},$

    where we took the minimum of each prefactor of the norms on the left-hand side and divided the inequality by this value. Gronwall’s lemma then implies that

    (55)

    $\displaystyle\big{(}g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2}\big{)}(t)+\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{X}}}}^{2}+\|u_{k}(t)\|_{\mathcal{H}}^{2}+\|\nabla u_{k}\|_{L^{2}_{t}{\mathcal{H}}}^{2}$

    $\displaystyle\quad\leq C(\alpha,T)\cdot\big{(}\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{X}}}^{2}+\|u_{k}^{0}\|_{\mathcal{H}}^{2}\big{)}\cdot\textup{exp}\bigg{(}\frac{T^{2\alpha-1}}{2\alpha-1}\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|^{2}_{H^{1}(D;W^{1,\infty}(\Omega))}\bigg{)}.$

    We note that the initial conditions of the Galerkin system are defined by $u_{k}^{0}=\Pi_{\mathcal{H}_{k}}u^{0}$ and ${\widehat{\psi}}_{k}^{0}=\Pi_{{\widehat{\mathcal{Y}}}_{k}}{\widehat{\psi}}^{0}$. Therefore, we have that $\|u_{k}^{0}\|_{\mathcal{H}}^{2}\leq\|u^{0}\|_{\mathcal{H}}^{2}$ and

    $\|M^{1/2}{\widehat{\psi}}_{k}^{0}\|_{H^{1}(D;W^{1,\infty}(\Omega))}^{2}\leq C\|M^{1/2}{\widehat{\psi}}^{0}\|_{H^{1}(D;H^{1+d/2+\delta}(\Omega))}^{2}.$

    We insert these bounds into the right-hand side of the inequality LABEL:Eq:EnergyBack1 and we thus arrive at the desired $k$-uniform energy estimate Eq. 53.

    ###
    5.2 Convergence of subsequences

    Having derived the $k$-uniform energy estimate Eq. 53 stated in Lemma 5.9,
    we shall extract weakly/weakly-$*$ converging subsequences of Galerkin solutions $(u_{k},{\widehat{\phi}}_{k})$. We shall also prove strong convergence of a subsequence $u_{k_{j}}$ in $L^{2}(0,T;\mathcal{H}_{0})$ in order to pass to the limit $j\to\infty$ in the nonlinear terms in the variational Navier–Stokes–Fokker–Planck system.

    ######
    Lemma 5.11.

    Let Section 5 hold and assume that $r\in[1,\infty)$ for $d=2$ and $r\in[1,6)$ for $d=3$; then, the sequence of Galerkin solutions $(u_{k},{\widehat{\phi}}_{k})$ from Lemma 5.3 contains a subsequence $(u_{k_{j}},{\widehat{\phi}}_{k_{j}})$ that admits the following convergences as $j\to\infty$:

    (56)

    $\displaystyle u_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup u$

    $\displaystyle\text{weakly-$*$ in }L^{\infty}(0,T;\mathcal{H}_{0}),$

    $\displaystyle u_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup u$

    $\displaystyle\text{weakly\phantom{-*} in }L^{2}(0,T;\mathcal{V}_{0})\cap L^{8/d}(0,T;L^{4}(\Omega)^{d}),$

    $\displaystyle{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup{\widehat{\phi}}$

    $\displaystyle\text{weakly\phantom{-*} in }L^{2}(0,T;{\widehat{\mathcal{X}}}),$

    $\displaystyle\partial_{t}u_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup\partial_{t}u$

    $\displaystyle\text{weakly\phantom{-*} in }L^{8/(4+d)}(0,T;\mathcal{V}_{0}^{\prime}),$

    $\displaystyle u_{k_{j}}$
    $\displaystyle\longrightarrow u$

    $\displaystyle\text{strongly\hskip 2.84526pt in }L^{2}\big{(}0,T;L^{r}(\Omega;\mathbb{R}^{d})\big{)},$

    $\displaystyle u_{k_{j}}$
    $\displaystyle\longrightarrow u$

    $\displaystyle\text{strongly\hskip 2.84526pt in }C([0,T];\mathcal{V}_{0}^{\prime}),$

    $\displaystyle\partial_{t}^{\alpha}{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup\partial_{t}^{\alpha}{\widehat{\phi}}$

    $\displaystyle\text{weakly\phantom{-*} in }L^{8/(4+d)}(0,T;{\widehat{\mathcal{Z}}}^{\prime}),$

    $\displaystyle{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\longrightarrow{\widehat{\phi}}$

    $\displaystyle\text{strongly\hskip 2.84526pt in }L^{2}(0,T;{\widehat{\mathcal{Y}}}),$

    $\displaystyle\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k_{j}})$
    $\displaystyle\longrightarrow\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}})$

    $\displaystyle\text{strongly\hskip 2.84526pt in }L^{2}\big{(}0,T;L^{2}(\Omega;\mathbb{R}^{d\times d})\big{)},$

    $\displaystyle g_{1-\alpha}*{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup g_{1-\alpha}*{\widehat{\phi}}$

    $\displaystyle\text{weakly-$*$ in }L^{\infty}(0,T;{\widehat{\mathcal{Y}}})\cap L^{2}(0,T;{\widehat{\mathcal{X}}}),$

    $\displaystyle g_{1-\alpha}*{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\longrightarrow g_{1-\alpha}*{\widehat{\phi}}$

    $\displaystyle\text{strongly\hskip 2.84526pt in }C([0,T];{\widehat{\mathcal{X}}}^{\prime})\cap L^{2}(0,T;{\widehat{\mathcal{Y}}}).$

    ######
    Proof 5.12.

    In Lemma 5.9 we stated various $k$-uniform bounds on $u_{k}$ and ${\widehat{\phi}}_{k}$. Thanks to the Banach–Alaoglu and Eberlein–Šmulian theorems, see [1, Theorem 8.10], there are weakly/weakly-$*$ converging subsequences $u_{k_{j}}$ and ${\widehat{\phi}}_{k_{j}}$. In particular, we obtain the convergences

    (57)

    $\displaystyle u_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup u$

    $\displaystyle\text{weakly-$*$ in }L^{\infty}(0,T;\mathcal{H}_{0}),$

    $\displaystyle u_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup u$

    $\displaystyle\text{weakly\phantom{-*} in }L^{2}(0,T;\mathcal{V}_{0}),$

    $\displaystyle{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup{\widehat{\phi}}$

    $\displaystyle\text{weakly\phantom{-*} in }L^{2}(0,T;{\widehat{\mathcal{X}}}).$

    We shall establish the strong convergence of $u_{k_{j}}$ in $L^{2}(0,T;\mathcal{H}_{0})$ by applying the Aubin–Lions compactness lemma; see Eq. 31. To this end, we need to bound the time derivative of $u_{k}$ in a suitable dual space. Let us therefore consider an arbitrary element $v\in L^{8/(4-d)}(0,T;\mathcal{V}_{0})$ and bound each of the terms appearing on the right-hand side of Eq. 45 below by means of Hölder’s inequality:

    $\displaystyle\int_{0}^{T}|\langle\partial_{t}u_{k},v\rangle_{\mathcal{V}_{0}}|\,\textup{d}t$
    $\displaystyle=\int_{0}^{T}\Big{|}-((u_{k}\cdot\nabla_{\!x})u_{k},\Pi_{\mathcal{H}_{k}}v)_{\mathcal{H}}$

    $\displaystyle\quad-\nu(\nabla_{\!x}u_{k},\nabla_{\!x}\Pi_{\mathcal{H}_{k}}v)_{\mathcal{H}}-k_{B}\mu_{T}\big{(}\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k}),\nabla_{\!x}\Pi_{\mathcal{H}_{k}}v\big{)}_{\mathcal{H}}\Big{|}\,\textup{d}t$

    $\displaystyle\leq C\int_{0}^{T}\Big{(}\|u_{k}\|_{L^{4}(\Omega)}\|u_{k}\|_{\mathcal{V}}\|\Pi_{\mathcal{H}_{k}}v\|_{L^{4}(\Omega)}$

    $\displaystyle\quad+\|u_{k}\|_{\mathcal{V}}\|\Pi_{\mathcal{H}_{k}}v\|_{\mathcal{V}}+\|\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k})\|_{\mathcal{H}}\|\Pi_{\mathcal{H}_{k}}v\|_{\mathcal{V}}\Big{)}\,\textup{d}t.$

    Hence, using Ladyzhenskaya’s inequality, we have that

    $\displaystyle\int_{0}^{T}|\langle\partial_{t}u_{k},v\rangle_{\mathcal{V}_{0}}|\,\textup{d}t$

    $\displaystyle\leq C\int_{0}^{T}\Big{(}\|u_{k}\|_{\mathcal{H}}^{1-d/4}\|u_{k}\|_{\mathcal{V}}^{1+d/4}\|\Pi_{\mathcal{H}_{k}}v\|_{L^{4}(\Omega)}+\|u_{k}\|_{\mathcal{V}}\|\Pi_{\mathcal{H}_{k}}v\|_{\mathcal{V}}+\|{\widehat{\phi}}_{k}\|_{\widehat{\mathcal{Y}}}\|\Pi_{\mathcal{H}_{k}}v\|_{\mathcal{V}}\Big{)}\,\textup{d}t$

    $\displaystyle\leq C\Big{(}\|u_{k}\|_{L^{\infty}{\mathcal{H}}}^{1-d/4}\|u_{k}\|_{L^{2}{\mathcal{V}}}^{1+d/4}\|v\|_{L^{8/(4-d)}{\mathcal{V}}}+\|u_{k}\|_{L^{2}{\mathcal{V}}}\|v\|_{L^{2}{\mathcal{V}}}+\|{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{Y}}}}\|v\|_{L^{2}{\mathcal{V}}}\Big{)}$

    $\displaystyle\leq C\|v\|_{L^{8/(4-d)}{\mathcal{V}}}.$

    This then implies that $\partial_{t}u_{k}$ is bounded in $L^{8/(4+d)}(0,T;\mathcal{V}_{0}^{\prime})$.
    It follows by the Aubin–Lions lemma Eq. 31 that

    (58)

    $\displaystyle\partial_{t}u_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup\partial_{t}u$

    $\displaystyle\text{weakly\phantom{-*} in }L^{8/(4+d)}(0,T;\mathcal{V}_{0}^{\prime}),$

    $\displaystyle u_{k_{j}}$
    $\displaystyle\longrightarrow u$

    $\displaystyle\text{strongly\hskip 2.84526pt in }L^{2}\big{(}0,T;L^{r}(\Omega;\mathbb{R}^{d})\big{)},$

    $\displaystyle u_{k_{j}}$
    $\displaystyle\longrightarrow u$

    $\displaystyle\text{strongly\hskip 2.84526pt in }C([0,T];\mathcal{V}_{0}^{\prime}),$

    where $r\in[1,\infty)$ for $d=2$ and $r\in[1,6)$ for $d=3$.

    Similarly, we consider an arbitrary element $\widehat{\zeta}\in L^{\frac{8}{4-d}}(0,T;{\widehat{\mathcal{Z}}})$ and we recall that ${\widehat{\mathcal{Z}}}$ was defined in the beginning of Section 4.1. We test the Galerkin equation of ${\widehat{\phi}}_{k}$ with $\Pi_{\mathcal{H}_{k}}\widehat{\zeta}$ giving

    (59)

    $\displaystyle\int_{0}^{T}\!\!|\langle\partial_{t}^{\alpha}{\widehat{\phi}}_{k},\Pi_{\mathcal{H}_{k}}\widehat{\zeta}\rangle_{{\widehat{\mathcal{X}}}}|\,\textup{d}t$
    $\displaystyle=\!\int_{0}^{T}\!\!\Big{|}-((u_{k}\cdot\nabla_{\!x}){\widehat{\phi}}_{k},\Pi_{\mathcal{H}_{k}}\widehat{\zeta})_{\widehat{\mathcal{Y}}}-\frac{1}{2\lambda}(\nabla_{\!q}{\widehat{\phi}}_{k},\nabla_{\!q}\Pi_{\mathcal{H}_{k}}\widehat{\zeta})_{\widehat{\mathcal{Y}}}$

    $\displaystyle\quad-\varepsilon(\nabla_{\!x}{\widehat{\phi}}_{k},\nabla_{x}\Pi_{\mathcal{H}_{k}}\widehat{\zeta})_{\widehat{\mathcal{Y}}}+(\omega(u_{k})q{\widehat{\phi}}_{k},\nabla_{\!q}\Pi_{\mathcal{H}_{k}}\widehat{\zeta})_{\widehat{\mathcal{Y}}}\Big{|}\mathop{}\!\mathrm{d}t.$

    We note that $u_{k}$ is bounded in $L^{8/d}(0,T;L^{4}(\Omega)^{d})$ by the following interpolation result

    $\int_{0}^{T}\|u_{k}\|_{L^{4}}^{8/d}\,\textup{d}t\leq\int_{0}^{T}\|u_{k}\|_{\mathcal{H}}^{8/d-2}\|u_{k}\|^{2}_{\mathcal{V}}\,\textup{d}t\leq\|u_{k}\|_{L^{\infty}\mathcal{H}}^{8/d-2}\|u_{k}\|_{L^{2}\mathcal{V}}^{2}.$

    Regarding the last term in (59), we integrate by parts and estimate by Hölder’s inequality

    $\displaystyle-$
    $\displaystyle\big{(}\omega(u_{k})q{\widehat{\phi}},\nabla_{\!q}\widehat{\zeta}\big{)}_{\widehat{\mathcal{Y}}}$

    $\displaystyle=\big{(}\nabla_{\!x}{\widehat{\phi}}_{k}(\nabla_{\!q}\widehat{\zeta})^{\mathrm{T}}q,u_{k}\big{)}_{\widehat{\mathcal{Y}}}+\big{(}{\widehat{\phi}}_{k}\nabla_{\!x}\nabla_{\!q}\widehat{\zeta}q,u_{k}\big{)}_{\widehat{\mathcal{Y}}}+\big{(}u_{k}\cdot q,\nabla_{\!x}{\widehat{\phi}}_{k}\cdot\nabla_{\!q}\widehat{\zeta}\big{)}_{\widehat{\mathcal{Y}}}$

    $\displaystyle\quad+\big{(}u_{k}\cdot q,{\widehat{\phi}}_{k}\textup{div}_{x}\nabla_{\!q}\widehat{\zeta}\big{)}_{\widehat{\mathcal{Y}}}$

    $\displaystyle\leq C\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{Y}}}}\big{(}\|\nabla_{\!q}\widehat{\zeta}\|_{L^{8/(4-d)}{\widehat{\mathcal{Y}}}}+\|\nabla_{\!x}\nabla_{\!q}\widehat{\zeta}\|_{L^{8/(4-d)}{\widehat{\mathcal{Y}}}}\big{)}\|q\|_{L^{\infty}}\|u_{k}\|_{L^{8/d}L^{4}}$

    $\displaystyle\quad+\|{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{X}}}}\|\nabla_{\!x}\nabla_{\!q}\widehat{\zeta}\|_{{L^{8/(4-d)}{\widehat{\mathcal{Y}}}}}\|q\|_{L^{\infty}}\|u_{k}\|_{L^{8/d}L^{4}}$

    $\displaystyle\quad+C\|u_{k}\|_{L^{8/d}L^{4}}\|q\|_{L^{\infty}}\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{Y}}}}\big{(}\|\nabla_{\!q}\widehat{\zeta}\|_{L^{8/(4-d)}{\widehat{\mathcal{Y}}}}+\|\nabla_{\!x}\nabla_{\!q}\widehat{\zeta}\|_{L^{8/(4-d)}{\widehat{\mathcal{Y}}}}\big{)}$

    $\displaystyle\quad+\|u_{k}\|_{L^{8/d}L^{4}}\|q\|_{L^{\infty}}\|{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{X}}}}\|\textup{div}_{x}\nabla_{\!q}\widehat{\zeta}\|_{L^{8/(4-d)}{\widehat{\mathcal{Y}}}}$

    $\displaystyle\leq C\|u_{k}\|_{L^{8/d}L^{4}}\|q\|_{L^{\infty}}\|{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{X}}}}\|\widehat{\zeta}\|_{L^{8/(4-d)}(0,T;{\widehat{\mathcal{Z}}})}.$

    Using this estimate, we can bound (59) as follows:

    (60)

    $\displaystyle\int_{0}^{T}|\langle\partial_{t}^{\alpha}{\widehat{\phi}}_{k},\Pi_{\mathcal{H}_{k}}\widehat{\zeta}\rangle_{{\widehat{\mathcal{X}}}}|\,\textup{d}t$

    $\displaystyle\leq C\Big{(}\|u_{k}\|_{L^{\infty}\mathcal{H}_{0}}\|{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{X}}}}\|\widehat{\zeta}\|_{L^{2}{\widehat{\mathcal{X}}}}+\|\nabla_{\!q}{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{Y}}}}\|\nabla_{\!q}\widehat{\zeta}\|_{L^{2}{\widehat{\mathcal{Y}}}}$

    $\displaystyle\quad+\|\nabla_{\!x}{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{Y}}}}\|\nabla_{x}\widehat{\zeta}\|_{L^{2}{\widehat{\mathcal{Y}}}}+\|u_{k}\|_{L^{8/d}L^{4}}\|q\|_{L^{\infty}}\|{\widehat{\phi}}_{k}\|_{L^{2}{\widehat{\mathcal{X}}}}\|\widehat{\zeta}\|_{L^{8/(4-d)}(0,T;{\widehat{\mathcal{Z}}})}\Big{)}$

    $\displaystyle\leq C\|\widehat{\zeta}\|_{L^{8/(4-d)}(0,T;{\widehat{\mathcal{Z}}})}.$

    Hence, we obtain the $k$-uniform boundedness of $\partial_{t}(g_{1-\alpha}*{\widehat{\phi}}_{k})=\partial_{t}^{\alpha}{\widehat{\phi}}_{k}$ in the space $L^{8/(4+d)}(0,T;{\widehat{\mathcal{Z}}}^{\prime})$, which is continuously embedded in $L^{8/(4+d)}(0,T;(H^{2}_{M}(\Omega\times D))^{\prime})$.
    Therefore, we are in the setting of the Gelfand triple

    ${\widehat{\mathcal{X}}}\mathrel{\mathrlap{{\mspace{4.0mu}\lhook}}{\hookrightarrow}}{\widehat{\mathcal{Y}}}\hookrightarrow\big{(}H_{M}^{2}(\Omega\times D)\big{)}^{\prime}.$

    We thus obtain from the fractional Aubin–Lions lemma, see (32), that

    (61)

    $\displaystyle\partial_{t}^{\alpha}{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\relbar\joinrel\rightharpoonup\partial_{t}^{\alpha}{\widehat{\phi}}$

    $\displaystyle\text{weakly\phantom{-*} in }L^{8/(4+d)}(0,T;{\widehat{\mathcal{Z}}}^{\prime}),$

    $\displaystyle{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\longrightarrow{\widehat{\phi}}$

    $\displaystyle\text{strongly\hskip 2.84526pt in }L^{2}(0,T;{\widehat{\mathcal{Y}}}).$

    The convolution $g_{1-\alpha}*{\widehat{\phi}}_{k}$ is bounded in $L^{2}(0,T;{\widehat{\mathcal{X}}})$ thanks to Young’s convolution inequality

    $\|g_{1-\alpha}*{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{X}}}}\leq\|g_{1-\alpha}\|_{L^{1}_{t}}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{X}}}}\leq CT^{1-\alpha}\|{\widehat{\phi}}_{k}\|_{L^{2}_{t}{\widehat{\mathcal{X}}}}.$

    Moreover, $g_{1-\alpha}*{\widehat{\phi}}_{k}$ is bounded in $L^{\infty}(0,T;{\widehat{\mathcal{Y}}})$ by the following chain of estimates:

    $\displaystyle\|g_{1-\alpha}*{\widehat{\phi}}_{k}\|_{L^{\infty}{\widehat{\mathcal{Y}}}}$

    $\displaystyle\quad\leq\sup_{t\in(0,T)}\int_{0}^{t}g_{1-\alpha}(t-s)\|{\widehat{\phi}}_{k}(s)\|_{\widehat{\mathcal{Y}}}\,\textup{ds}$

    $\displaystyle\quad\leq\sup_{t\in(0,T)}\int_{0}^{t}g_{1-\alpha}(t-s)\|{\widehat{\phi}}_{k}(s)-{\widehat{\psi}}_{k}^{0}g_{\alpha}(s)\|_{\widehat{\mathcal{Y}}}\,\textup{ds}+(g_{1-\alpha}*g_{\alpha})(t)\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}$

    $\displaystyle\quad\leq\sup_{t\in(0,T)}\int_{0}^{t}g_{1-\alpha}(t-s)\|{\widehat{\phi}}_{k}(s)-{\widehat{\psi}}_{k}^{0}g_{\alpha}(s)\|^{2}_{\widehat{\mathcal{Y}}}\,\textup{ds}+\frac{1}{4}\int_{0}^{t}g_{1-\alpha}(t-s)\,\textup{ds}+\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}}$

    $\displaystyle\quad=\sup_{t\in(0,T)}(g_{1-\alpha}*\|{\widehat{\phi}}_{k}-{\widehat{\psi}}_{k}^{0}g_{\alpha}\|_{\widehat{\mathcal{Y}}}^{2})(t)+\frac{1}{4}g_{2-\alpha}(T)+\|{\widehat{\psi}}_{k}^{0}\|_{\widehat{\mathcal{Y}}},$

    and the first term on the right-hand side is bounded by Lemma 5.9. Since we have already proved a bound on $\partial_{t}(g_{1-\alpha}*{\widehat{\phi}}_{k})=\partial_{t}^{\alpha}{\widehat{\phi}}_{k}$, see (60), we may use the Aubin–Lions lemma, see Eq. 31, to obtain the following strong convergence results:

    (62)

    $\displaystyle g_{1-\alpha}*{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\longrightarrow g_{1-\alpha}*{\widehat{\phi}}$

    $\displaystyle\text{strongly\hskip 2.84526pt in }L^{2}(0,T;{\widehat{\mathcal{Y}}}),$

    $\displaystyle g_{1-\alpha}*{\widehat{\phi}}_{k_{j}}$
    $\displaystyle\longrightarrow g_{1-\alpha}*{\widehat{\phi}}$

    $\displaystyle\text{strongly\hskip 2.84526pt in }C([0,T];{\widehat{\mathcal{X}}}^{\prime}).$

    Lastly, we note that the mapping $Mg_{1-\alpha}*\varphi\mapsto\mathcal{C}(Mg_{1-\alpha}*\varphi)$ is linear and continuous thanks to (41), and therefore we have from (62)1 that

    (63)

    $\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}}_{k_{j}})\longrightarrow\mathcal{C}(Mg_{1-\alpha}*{\widehat{\phi}})\quad\text{ strongly in }L^{2}\big{(}0,T;L^{2}(\Omega;\mathbb{R}^{d\times d})\big{)}.$

    ###
    5.3 Passage to the limit

    Next, we pass to the limit $j\to\infty$ in the time-integrated $k_{j}$-th Galerkin system Eq. 45, Eq. 46. Specifically, we shall use the convergence results stated in the preceding lemma to show that the weak limits, $u$ and ${\widehat{\phi}}$, satisfy the variational Navier–Stokes–Fokker–Planck system in the sense of Definition 5.1.

    ######
    Proof 5.13 (Proof of Theorem 5.2).

    We consider the time-integrated Galerkin system

    (64)

    $\displaystyle\int_{0}^{T}\Big{(}\langle\partial_{t}u_{k_{j}},v\rangle_{V}+((u_{k_{j}}\cdot\nabla_{\!x})u_{k_{j}},v)_{\mathcal{H}}$

    $\displaystyle\quad+\nu(\nabla_{\!x}u_{k_{j}},\nabla_{\!x}v)_{\mathcal{H}}+k\mu(C(Mg_{1-\alpha}*{\widehat{\phi}}_{k_{j}}),\nabla_{\!x}v)_{\mathcal{H}}\Big{)}\eta(t)\,\textup{d}t=0$

    (65)

    $\displaystyle\int_{0}^{T}-(g_{1-\alpha}*{\widehat{\phi}}_{k_{j}},\widehat{\zeta})_{\widehat{\mathcal{Y}}}\eta^{\prime}(t)+\Big{(}((u_{k_{j}}\cdot\nabla_{\!x}){\widehat{\phi}}_{k_{j}},\widehat{\zeta})_{\widehat{\mathcal{Y}}}+\frac{1}{2\lambda}(\nabla_{\!q}{\widehat{\phi}}_{k_{j}},\nabla_{\!q}\widehat{\zeta})_{\widehat{\mathcal{Y}}}$

    $\displaystyle\quad+\varepsilon(\nabla_{\!x}{\widehat{\phi}}_{k_{j}},\nabla\widehat{\zeta})_{\widehat{\mathcal{Y}}}-(\omega(u_{k_{j}})q{\widehat{\phi}}_{k_{j}},\nabla_{\!q}\widehat{\zeta})_{\widehat{\mathcal{Y}}}\Big{)}\eta(t)\,\textup{d}t=0,$

    for all $v\in\mathcal{H}_{k_{j}}$, $\eta\in C_{0}^{\infty}(0,T)$ and $\widehat{\zeta}\in{\widehat{\mathcal{Y}}}_{k_{j}}$.
    Passing to the limit $j\rightarrow\infty$ in (64) using Eq. 57–Eq. 63
    is standard, and results in (42). It therefore remains to pass to the limit $j\rightarrow\infty$ in (65). In particular, the convergence of the linear terms follow immediately by weak convergence, and we only consider the two nonlinear terms.
    We note that ${\widehat{\phi}}_{k_{j}}\to{\widehat{\phi}}$ strongly in $L^{2}\big{(}0,T;{\widehat{\mathcal{Y}}})$ and $\omega(u_{k_{j}})\to\omega(u)$ weakly in $L^{2}(0,T;\mathcal{H}_{0})$, from which we deduce that

    $\int_{0}^{T}(\omega(u_{k_{j}})q{\widehat{\phi}}_{k_{j}},\nabla_{\!q}\widehat{\zeta})_{\widehat{\mathcal{Y}}}\eta(t)\,\textup{d}t\longrightarrow\int_{0}^{T}(\omega(u)q{\widehat{\phi}},\nabla_{\!q}\widehat{\zeta})_{\widehat{\mathcal{Y}}}\eta(t)\,\textup{d}t,$

    as $j\to\infty$. With the same reasoning, we are able to show that

    $\int_{0}^{T}((u_{k_{j}}\cdot\nabla_{\!x}){\widehat{\phi}}_{k_{j}},\widehat{\zeta})_{{\widehat{\mathcal{Y}}}}\eta(t)\,\textup{d}t\longrightarrow\int_{0}^{T}((u\cdot\nabla_{\!x}){\widehat{\phi}},\widehat{\zeta})_{{\widehat{\mathcal{Y}}}}\eta(t)\,\textup{d}t\quad\mbox{as $j\to\infty$}.$

    We use the density of $\cup_{k=1}^{\infty}\mathcal{H}_{k}$ in $\mathcal{V}$ and of $\cup_{k=1}^{\infty}{\widehat{\mathcal{Y}}}_{k}$ in $H^{2}_{M}(\Omega\times D)$, which completes the proof by observing that the tuple $(u,{\widehat{\phi}})$ satisfies the variational form of the time-fractional Navier–Stokes–Fokker–Planck system as stated in Definition 5.1.

    It remains to check that the initial conditions are satisfied. First, we obtain the convergence $u_{k_{j}}(0)\to u(0)$ in $\mathcal{V}_{0}^{\prime}$ as $j\rightarrow\infty$; see again Eq. 56. However, by definition, $u_{k_{j}}(0)=\Pi_{\mathcal{H}_{k_{j}}}u^{0}$, which converges to $u^{0}$ in $\mathcal{H}_{0}$
    as $j\rightarrow\infty$. By the uniqueness of the limit it follows that $u(0)=u^{0}$. Regarding the solution of the Fokker–Planck equation, we use again the strong convergence Eq. 56 to conclude $(g_{1-\alpha}*{\widehat{\phi}})(0)={\widehat{\psi}}^{0}$.

    Having proved that a weak solution tuple $(u,{\widehat{\phi}})$ to the time-fractional system in the sense of Definition 5.1 exists, we return to the original variable $\psi:=g_{1-\alpha}*(M{\widehat{\phi}})$, whose evolution is governed by the time-fractional Fokker–Planck equation Section 2.2. Indeed, $g_{1-\alpha}*\psi=g_{2-2\alpha}*\phi$, which is continuous for $\alpha>1/2$ and $(g_{1-\alpha}*\psi)(0)=0$ as we have originally assumed in the model transformation. In this sense, we have also shown the existence of a variational solution tuple $(u,\psi)$ to the original time-fractional model.

    ## Conclusions and outlook

    In this paper, we investigated the well-posedness of a coupled Navier–Stokes–Fokker–Planck system with a time-fractional derivative. Such systems arise in the kinetic theory of polymeric liquid solutions with noninteracting polymer chains. We outlined the derivation of the model from a subordinated Langevin equation and considered the case of a finitely extensible nonlinear elastic (FENE-type) dumbbell model with a corotational drag term. We proved the existence of large-data global-in-time weak solutions to the corotational time-fractional model of order $\alpha\in(\tfrac{1}{2},1)$ and derived a uniform energy inequality by considering a nonstandard and novel testing procedure.
    The existence of weak solutions to the general noncorotational time-fractional FENE model is an open problem, which will be studied in a forthcoming paper by using a different testing procedure; see [12] for the integer-order setting (corresponding to $\alpha=1$). Concerning the numerical approximation of the time-fractional system considered here we refer the reader to the recent paper
    [15].

    ## References

      *
    [1]

    H. W. Alt, Linear Functional Analysis: An Application-Oriented
    Introduction, Springer, 2016.

      *
    [2]

    T. M. Atanackovic, S. Pilipovic, B. Stankovic, and D. Zorica, Fractional Calculus with Applications in Mechanics: Wave Propagation, Impact
    and Variational Principles, John Wiley & Sons, 2014.

      *
    [3]

    D. Baleanu, K. Diethelm, E. Scalas, and J. J. Trujillo, Fractional
    Calculus: Models and Numerical Methods, World Scientific, 2012.

      *
    [4]

    E. Barkai, Fractional Fokker–Planck equation, solution, and
    application, Physical Review E, 63 (2001), p. 046118.

      *
    [5]

    E. Barkai, R. Metzler, and J. Klafter, From continuous time random
    walks to the fractional Fokker–Planck equation, Physical Review E, 61
    (2000), p. 132.

      *
    [6]

    J. W. Barrett, C. Schwab, and E. Süli, Existence of global weak
    solutions for some polymeric flow models, Mathematical Models and Methods
    in Applied Sciences, 15 (2005), pp. 939–983.

      *
    [7]

    J. W. Barrett and E. Süli, Existence of global weak solutions
    to some regularized kinetic models for dilute polymers, Multiscale Modeling
    & Simulation, 6 (2007), pp. 506–546.

      *
    [8]

    J. W. Barrett and E. Süli, Existence of global weak solutions
    to dumbbell models for dilute polymers with microscopic cut-off,
    Mathematical Models and Methods in Applied Sciences, 18 (2008), pp. 935–971.

      *
    [9]

    J. W. Barrett and E. Süli, Numerical approximation of
    corotational dumbbell models for dilute polymers, IMA Journal of Numerical
    Analysis, 29 (2009), pp. 937–959.

      *
    [10]

    J. W. Barrett and E. Süli, Existence and equilibration of
    global weak solutions to Hookean-type bead-spring chain models for dilute
    polymers, Mathematical Models and Methods in Applied Sciences, Vol. 22, No.
    5, 2012, 22 (2010).

      *
    [11]

    J. W. Barrett and E. Süli, Existence of global weak solutions
    to Fokker–Planck and Navier–Stokes–Fokker–Planck equations in
    kinetic models of dilute polymers, Discrete & Continuous Dynamical Systems
    – S, 3 (2010), p. 371.

      *
    [12]

    J. W. Barrett and E. Süli, Existence and equilibration of global
    weak solutions to kinetic models for dilute polymers I: Finitely extensible
    nonlinear bead-spring chains, Mathematical Models and Methods in Applied
    Sciences, 21 (2011), pp. 1211–1289.

      *
    [13]

    J. W. Barrett and E. Süli, Finite element approximation of
    kinetic dilute polymer models with microscopic cut-off, ESAIM: Mathematical
    Modelling and Numerical Analysis, 45 (2011), pp. 39–89.

      *
    [14]

    J. W. Barrett and E. Süli, Finite element approximation of
    finitely extensible nonlinear elastic dumbbell models for dilute polymers,
    ESAIM: Mathematical Modelling and Numerical Analysis, 46 (2012),
    pp. 949–978.

      *
    [15]

    J. Beddrich, E. Süli, and B. Wohlmuth, Numerical simulation of the
    time-fractional Navier–Stokes–Fokker–Planck equation and
    applications to polymeric fluids.

    submitted to Journal of Computational Physics (2023).

      *
    [16]

    R. B. Bird, C. F. Curtiss, R. C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, Volume 2: Kinetic Theory, Wiley, 1987.

      *
    [17]

    M. Bulicek, J. Málek, and E. Süli, Existence of global weak
    solutions to implicitly constituted kinetic models of incompressible
    homogeneous dilute polymers, Communications in Partial Differential
    Equations, 38 (2013), pp. 882–924.

      *
    [18]

    A. V. Busuioc, I. S. Ciuperca, D. Iftimie, and L. I. Palade, The
    FENE dumbbell polymer model: Existence and uniqueness of solutions for
    the momentum balance equation, Journal of Dynamics and Differential
    Equations, 26 (2014), pp. 217–241.

      *
    [19]

    A. Carbotti and G. E. Comi, A note on Riemann–Liouville
    fractional Sobolev spaces, Communications on Pure & Applied Analysis, 20
    (2021), p. 17.

      *
    [20]

    W. Chen, H. Sun, and X. Li, Fractional Derivative Modeling in
    Mechanics and Engineering, Springer, 2022.

      *
    [21]

    T. Debiec and E. Süli, Corotational Hookean models of dilute
    polymeric fluids: Existence of global weak solutions, weak-strong
    uniqueness, equilibration, and macroscopic closure, SIAM Journal on
    Mathematical Analysis, 55 (2023), pp. 310–346.

      *
    [22]

    K. Diethelm, The Analysis of Fractional Differential Equations: An
    Application-Oriented Exposition using Differential Operators of Caputo
    type, Springer, 2010.

      *
    [23]

    H. Fallahgoul, S. Focardi, and F. Fabozzi, Fractional Calculus and
    Fractional Processes with Applications to Financial Economics: Theory and
    Application, Academic Press, 2016.

      *
    [24]

    M. Fritz, Well-posedness and simulation of weak solutions to the
    time-fractional Fokker–Planck equation with general forcing.

    submitted to Discrete and Continuous Dynamical Systems – B
    (2023).

      *
    [25]

    M. Fritz, U. Khristenko, and B. Wohlmuth, Equivalence between a
    time-fractional and an integer-order gradient flow: The memory effect
    reflected in the energy, Advances in Nonlinear Analysis, 12 (2023).

      *
    [26]

    M. Fritz, C. Kuttler, M. L. Rajendran, B. Wohlmuth, and L. Scarabosio,
    On a subdiffusive tumour growth model with fractional time
    derivative, IMA Journal of Applied Mathematics, 86 (2021), pp. 688–729.

      *
    [27]

    M. Fritz, M. L. Rajendran, and B. Wohlmuth, Time-fractional
    Cahn–Hilliard equation: Well-posedness, regularity, degeneracy, and
    numerical solutions, Computers & Mathematics with Applications, 108 (2022),
    pp. 66–87.

      *
    [28]

    B. Henry, T. Langlands, and S. Wearne, Anomalous diffusion with
    linear reaction dynamics: From continuous time random walks to fractional
    reaction-diffusion equations, Physical Review E, 74 (2006), p. 031116.

      *
    [29]

    B. I. Henry, T. Langlands, and P. Straka, Fractional
    Fokker–Planck equations for subdiffusion with space-and time-dependent