Derivation of the Navier - Stokes - Poisson system with radiation for an accretion disk

TL;DR

This paper derives a 2-D Navier-Stokes-Poisson system with radiation from a 3-D compressible fluid model, demonstrating convergence of weak solutions to strong solutions under certain conditions, advancing understanding of accretion disk dynamics.

Contribution

It establishes the rigorous derivation and convergence of 3-D compressible radiation fluid solutions to 2-D models relevant for accretion disks, incorporating rotation, gravitation, and radiation effects.

Findings

Weak solutions converge to 2-D strong solutions under small Froude number.

Convergence holds for all times less than the maximal lifespan of the 2-D solutions.

The derived models incorporate radiation effects in accretion disk dynamics.

Abstract

We study the 3-D compressible barotropic radiation fluid dynamics system describing the motion of the compressible rotating viscous fluid with gravitation and radiation confined to a straight layer. We show that weak solutions in the 3-D domain converge to the strong solution of the rotating 2-D Navier-Stokes-Poisson system with radiation for all times less than the maximal life time of the strong solution of the 2-D system when the Froude number is small or to the strong solution of the rotating pure 2-D Navier- Stokes system with radiation.

Equations

$${\partial }_t\varrho +{\mathrm{div}}_x(\varrho \vec{u})=0,$$

$${\partial }_t(\varrho \vec{u})+{\mathrm{div}}_x(\varrho \vec{u}\otimes \vec{u})+{\nabla }_{x}p(\varrho )+\varrho \vec{\chi }\times \vec{u}={\mathrm{div}}_x\text{mbox}\text{F}S+\varrho {\nabla }_x\varphi +\varrho {\nabla }_x|\vec{\chi }\times \vec{x}{|}^2+{\vec{S}}_F.$$

$${\vec{S}}_F=({\sigma }_a+{\sigma }_s){\int }_0^{\infty }{\int }_{{\mathcal{S}}^2}\vec{\varsigma }I{\mathrm{d}}_{\vec{\varsigma }}\sigma \mathrm{d}\nu ,$$

$$-\Delta \varphi =4\pi G(\eta \varrho +(1-\eta )g)\text{mbox}in(0,T)\times {\Omega }_{ϵ}.$$

$$\varphi (t,\vec{x})=G{\int }_{\text{mbox}{\text{FF}R}^3}\frac{\eta \varrho (t,\vec{y})+(1-\eta )g(\vec{y})}{\left|\vec{x}-\vec{y}\right|}\mathrm{d}\vec{y}.$$

$$\int_{\mbox{\FF R}^{3}}|\nabla K({\vec{x}}-{\vec{y}})\big{|}\big{(}{\eta\varrho(t,\vec{y})+(1-\eta)g(\vec{y})}\big{)}\big{|}\ {\rm d}\vec{y}<\infty,$$

$$\frac{1}{c}{\partial }_{t}I+\vec{\varsigma }\cdot {\nabla }_{x}I=S,$$

$$S=\sigma_{a}(B-I)+\sigma_{s}\Big{(}\frac{1}{4\pi}\int_{S^{2}}I\ {\rm d}_{\vec{\varsigma}}\sigma-I\Big{)}.$$

$$0\le {\sigma }_s(\varrho ),{\sigma }_a(\varrho )\le c_1,$$

$$\sigma_{a}(\varrho)B(\nu,\varrho)\big{(}1+B(\nu,\varrho)\big{)}\leq h(\nu),\ h\in L^{1}(0,\infty)$$

$${\partial }_{\varrho }{\sigma }_a(\varrho ),{\partial }_{\varrho }{\sigma }_s(\varrho ),{\partial }_{\varrho }B(\varrho ,\nu ),B(\varrho ,\nu )\le c_2.$$

$$p\in C([0,\infty ))\cap C^1((0,\infty )),p(0)=0,p^{\prime }(\varrho )>0\text{mbox}forall\varrho >0,$$

$$\underset{\varrho \to \infty }{\lim }\frac{p^{\prime }(\varrho )}{{\varrho }^{\gamma -1}}=a>0$$

$$\mbox{\F S}=\mu\big{(}\nabla_{x}{\vec{u}}+\nabla_{x}^{t}{\vec{u}}-\frac{2}{3}{\rm div}_{x}{\vec{u}}\,\mbox{\F I}\big{)}+\xi\ {\rm div}_{x}{\vec{u}}\ \mbox{\F I},$$

$$\varrho (0,\vec{x})={\tilde{\varrho }}_{0,ϵ}(\vec{x}),\vec{u}(0,\vec{x})={\tilde{\vec{u}}}_{0,ϵ}(\vec{x}),I(0,\vec{x},\vec{\varsigma },\nu )={\tilde{I}}_{0,ϵ}(\vec{x},\vec{\varsigma },\nu ),\vec{x}\in {\Omega }_{ϵ},\varsigma \in {\mathcal{S}}^2,\nu \in \text{mbox}{\text{F}R}^{+}$$

$$\vec{u}{|}_{\partial \omega \times (0,ϵ)}=\vec{0}$$

$$\vec{u}\cdot \vec{n}{|}_{\omega \times \{0,ϵ\}}=0,[\text{mbox}\text{F}S({\nabla }_x\vec{u})\vec{n}]\times \vec{n}{|}_{\omega \times \{0,ϵ\}}=\vec{0}.$$

$$u_3=0\text{mbox}on\omega \times \{0,ϵ\}.$$

$$I(t,{\vec{x}},\vec{\varsigma},\nu)=0\ \mbox{for}\ ({\vec{x}},\vec{\varsigma})\in\Gamma^{1}_{-}\equiv\left\{({\vec{x}},\vec{\varsigma})\ \Big{|}\ (\vec{x},\vec{\varsigma})\in\partial\omega\times(0,\epsilon)\times{\mathcal{S}}^{2},\ \vec{\varsigma}\cdot\vec{n}\leq 0\right\}.$$

$$\begin{array}[]{c}I(t,{\vec{x}},\vec{\varsigma},\nu)=I(t,\vec{x},\vec{\varsigma}-2(\vec{\varsigma}\cdot\vec{n})\vec{n},\nu)\\ \mbox{for}\ ({\vec{x}},\vec{\varsigma})\in\Gamma^{2}_{-}\equiv\left\{({\vec{x}},\vec{\varsigma})\ \Big{|}\ (\vec{x},\vec{\varsigma})\in\omega\times\{0,\epsilon\}\times{\mathcal{S}}^{2},\ \vec{\varsigma}\cdot\vec{n}\leq 0\right\}.\end{array}$$

$$(x_h,ϵx_3)\in {\Omega }_{ϵ}\mapsto (x_h,x_3)\in \Omega ,\text{mbox}where{x}_h=(x_1,x_2)\in \omega ,x_3\in (0,1),$$

$${\nabla }_{ϵ}=({\nabla }_h,\frac{1}{ϵ}{\partial }_{x_3}),{\mathrm{div}}_{ϵ}\vec{u}={\mathrm{div}}_h{\vec{u}}_h+\frac{1}{ϵ}{\partial }_{x_3}{u}_3,$$

$${\vec{x}}_h=(x_1,x_2),{\vec{u}}_h=(u_1,u_2),{\nabla }_h=({\partial }_{x_1},{\partial }_{x_2}),$$

$${\mathrm{div}}_h{\vec{u}}_h={\partial }_{x_1}{u}_1+{\partial }_{x_2}{u}_2.$$

$${\partial }_t\varrho +{\mathrm{div}}_{ϵ}(\varrho \vec{u})=0,$$

$${\partial }_t(\varrho \vec{u})+{\mathrm{div}}_{ϵ}(\varrho \vec{u}\otimes \vec{u})+{\nabla }_{ϵ}p(\varrho )+\varrho \vec{\chi }\times \vec{u}$$

$$={\mathrm{div}}_{ϵ}\text{mbox}\text{F}S({\nabla }_{ϵ}\vec{u})+\frac{1}{F{r}^2}\varrho {\nabla }_{ϵ}\varphi +\varrho {\nabla }_{ϵ}|\vec{\chi }\times \vec{x}{|}^2+{\vec{S}}_F,$$

$$\partial_{t}I+{\vec{\varsigma}}\cdot\nabla_{x}I=S=\sigma_{a}\left(B-I\right)+\sigma_{s}\Big{(}\frac{1}{4\pi}\int_{{\cal S}^{2}}I\ {\rm d}_{{\vec{\varsigma}}}\sigma-I\Big{)},$$

$${\nabla }_{ϵ}\varphi (t,\vec{x})=ϵ{\int }_{\Omega }\eta \varrho (t,\vec{y})\frac{(x_1-y_1,x_2-y_2,ϵ(x_3-y_3))}{(|{\vec{x}}_h-{\vec{y}}_h{|}^2+{ϵ}^2(x_3-y_3{)}^2{)}^{\frac{3}{2}}}\mathrm{d}\vec{y}$$

$$+{\int }_{\text{mbox}{\text{FF}R}^3}(1-\eta )g(\vec{y})\frac{(x_1-y_1,x_2-y_2,ϵx_3-y_3)}{(|{\vec{x}}_h-{\vec{y}}_h{|}^2+(ϵx_3-y_3{)}^2{)}^{\frac{3}{2}}}\mathrm{d}\vec{y}=:ϵ\eta {\vec{\Phi }}_1+(1-\eta ){\vec{\Phi }}_2=:\vec{\Phi },$$

Full text

We show that weak solutions in the 3-D domain converge to the strong solution of

— the rotating 2-D Navier–Stokes–Poisson system with radiation in $\omega$ as $\epsilon\to 0$ for all times less than the maximal life time of the strong solution of the 2-D system when the Froude number is small ($Fr={\mathcal{O}}(\sqrt{\epsilon}))$,

— the rotating pure 2-D Navier–Stokes system with radiation in $\omega$ as $\epsilon\to 0$ when $Fr={\mathcal{O}}(1)$.

1 CEA, DAM, DIF, F-91297 Arpajon, France

2 Institute of Mathematics of the Academy of Sciences of the Czech Republic

Žitná 25, 115 67 Praha 1, Czech Republic

3 Charles University, Faculty of Mathematics and Physics

Mathematical Inst. of Charles University, Sokolovská 83, 186 75 Prague 8, Czech Republic

4Dpto. Ecuaciones Diferenciales y Análisis Numérico and IMUS,

Facultad de Matemáticas, Universidad de Sevilla, C/Tarfia, s/n, 41012 Sevilla, Spain

Key words: Navier–Stokes–Poisson system, radiation, rotation, Froude number, accretion disk, weak solution, thin domain, dimension reduction.

1 Introduction

Our aim in this work is the rigorous derivation of the equations describing objects called “accretion disks” which are quasi planar structures observed in various places in the universe.

From a naive point of view, if a massive object attracts matter distributed around it through the Newtonian gravitation in presence of a high angular momentum, this matter is not accreted isotropically around the central object but forms a thin disk around it. As the three main ingredients claimed by astrophysicists for explaining the existence of such objects are gravitation, angular momentum and viscosity (see [21] [22] [24] for detailed presentations), a reasonable framework for their study seems to be a viscous self-gravitating rotating fluid system of equations.

These disks are indeed three-dimensional but their size in the “third” dimension is usually very small, therefore they are often modeled as two-dimensional structures. Our goal in this paper is to derive rigorously the fluid equations of the disk from the equations set in a “thin” cylinder of thickness $\epsilon$ by passing to the limit $\epsilon\to 0^{+}$ and applying recent techniques of dimensional reduction introduced and applied in various situations by P. Bella, E. Feireisl, D. Maltese, A. Novotný and R. Vodák (see [2], [17], [27] and [28]).

The mathematical model which we consider is the compressible barotropic Navier–Stokes–Poisson system with radiation ([8], [9], [10]) describing the motion of a viscous radiating fluid confined in a bounded straight layer $\Omega_{\epsilon}=\omega\times(0,\epsilon)$, where $\omega\subset\mbox{\F R}^{2}$ has smooth boundary. Moreover, as we suppose a global rotation of the system, some new terms appear due to the change of frame.

Concerning gravitation a modelization difficulty appears as we consider the restriction to $\Omega_{\epsilon}$ of the solution of the Poisson equation in $\mbox{\F R}^{3}$: when the thickness of the cylinder tends to zero, a simple argument shows that the gravitational potential given by the Poisson equation in the whole space goes to zero. So if we want to recover the presence of gravitation at the limit, and then keep track of the physical situation, we will have to impose some scaling conditions. In fact as the limit problem will not depend on $x_{3}$, the flow is stratified and we expect that the scaling involves naturally the Froude number; see also [7].

More precisely, the system of equations giving the evolution of the mass density $\varrho=\varrho(t,{\vec{x}})$ and the velocity field ${\vec{u}}={\vec{u}}(t,{\vec{x}})=(u_{1},u_{2},u_{3})$, as functions of the time $t\in(0,T)$ and the spatial coordinate ${\vec{x}}=(x_{1},x_{2},x_{3})\in\Omega_{\epsilon}\subset\mbox{\F R}^{3}$, reads as follows:

[TABLE]

[TABLE]

On the right-hand side of (1.2) the radiative momentum $\vec{S}_{F}$ appears, given by

[TABLE]

where the unknown function $I=I(t,{\vec{x}},{\vec{\varsigma}},\nu)$ is the radiative intensity; see below for more details concerning the quantities describing the radiative effects.

The gravitational body forces are represented by the force term $\varrho\nabla_{x}\phi$, where the potential $\phi$ obeys Poisson’s equation

[TABLE]

Above, $G$ is the Newton constant and $g$ is a given function, modelling the external gravitational effect. Solving (1.4) in the whole space and supposing that $\varrho$ is extended by 0 outside $\Omega_{\epsilon}$, we have

[TABLE]

The parameter $\eta$ may take the values 0 or 1: for $\eta=1$ self-gravitation is present and for $\eta=0$ gravitation acts only as an external field (some astrophysicists consider self-gravitation of accretion disks as small compared to the external attraction by a given massive central object modelled by $g$, see [24]). Note that for the simplicity reasons we assume the external gravitation to be time independent.

We suppose that $g$ belongs to the regularity class such that integral (1.5) converges. Moreover, since in the momentum equation the term $\nabla_{x}\phi$ appears, we also need that

[TABLE]

where $K({\vec{x}}-{\vec{y}})=\frac{1}{|{\vec{x}}-{\vec{y}}|}$.

The effect of radiation is incorporated into the system through the radiative intensity $I=I(t,{\vec{x}},\vec{\varsigma},\nu)$, depending, besides the variables $t,{\vec{x}}$, on the direction vector $\vec{\varsigma}\in\mathcal{S}^{2}$, where $\mathcal{S}^{2}$ denotes the unit sphere in $\mbox{\F R}^{3}$, and the frequency $\nu>0$. The action of radiation is then expressed in term of integral average $\vec{S}_{F}$ with respect to the variables ${\vec{\varsigma}}$ and $\nu$.

The evolution of the compressible viscous barotropic flow is coupled to radiation through radiative transfer equation [4] which reads

[TABLE]

where $c$ is the speed of light. The radiative source $S:=S_{a}+S_{s}$ is the sum of an emission–absorption term $S_{a}:=\sigma_{a}(B(\nu,\varrho)-I)$ and a scattering contribution $S_{s}:=\sigma_{s}(\widetilde{I}-I)$, where $\widetilde{I}:=\frac{1}{4\pi}\int_{{\mathcal{S}}^{2}}I\ {\rm d}_{{\vec{\varsigma}}}\sigma$. The radiation source $S$ then reads

[TABLE]

We further assume:

• Isotropy: The coefficients $\sigma_{a}$, $\sigma_{s}$ are independent of $\vec{\varsigma}$.
• Grey hypothesis: The coefficients $\sigma_{a}$, $\sigma_{s}$ are independent of $\nu$.

The function $B=B(\nu,\varrho)$ measures the distance from equilibrium and is a barotropic equivalent of the Planck function.

Furthermore, we take

[TABLE]

[TABLE]

for any $\varrho\geq 0$. Note that relations (1.8–1.9) represent “cut-off” hypotheses at large density.

We need one more assumption on the radiative quantities,

[TABLE]

Assumption (1.9) is needed in the a priori estimate to get existence of a weak solution, assumption (1.10) will be important later in order to get estimates of the remainder in the relative entropy inequality.

Our system is globally rotating at uniform velocity $\chi$ around the vertical direction $\vec{e}_{3}$ and we denote $\vec{\chi}=\chi\vec{e}_{3}$. The Coriolis acceleration $\varrho\vec{\chi}\times{\vec{u}}$ and the centrifugal force term $\varrho\nabla_{x}|\vec{\chi}\times\vec{x}|^{2}$ is therefore present (see [5]).

The pressure is a given function of density satisfying hypotheses

[TABLE]

[TABLE]

for a certain $\gamma>3/2$.

The viscous stress tensor 𝕊 fulfils Newton’s rheological law determined by

[TABLE]

where $\mu>0$ is the shear viscosity coefficient and $\xi\geq 0$ is the bulk viscosity coefficient.

Finally, the system is supplemented with the initial conditions

[TABLE]

and with the boundary conditions. Here, the situation is more complex. For the velocity, we consider the no slip boundary conditions on the boundary part $\partial\omega\times(0,\epsilon)$ (the lateral part of the domain)

[TABLE]

and slip boundary condition on the boundary part $\omega\times\{0,\epsilon\}$ (the top and bottom part of the layer)

[TABLE]

Let us remark that we have $\vec{n}=\pm\vec{e}_{3}$ on $\omega\times\{0,\epsilon\}$, hence the first condition in (1.15) can be rewritten as

[TABLE]

We imposed the slip condition on the boundary $\omega\times\{0,\epsilon\}$ in order to avoid difficulties in passing to the “infinitely thin” limit; using the no slip boundary condition on the top and bottom part of the layer would imply that the velocity converges to zero when we let $\epsilon\to 0^{+}$.

Similar problem we meet with the radiative intensity. We consider at the lateral part of the boundary the condition

[TABLE]

Considering the same condition also on the top and bottom part of the layer (i.e., for $\vec{x}\in\omega\times\{0,\epsilon\}$) would lead to a situation we try to avoid: in the limit, the radiation disappears. We therefore consider

[TABLE]

This boundary condition is called specular reflection. More details needed for our paper will be given later, see also [1] for further comments and different possibilities.

Our proof will be based on the relative entropy inequality, developed by Feireisl, Novotný and coworkers in [13] and [12]. Recall, however, that the relative entropy inequality was first introduced in the context of hyperbolic equations in the work of C. Dafermos [6], then developed by A. Mellet and A. Vasseur [19], L. Saint-Raymond [25] and finally extended to the compressible barotropic case by P. Germain [15].

Remark 1.1

The relativistic version of system (1.1–1.7) has been introduced by Pomraning [23] and Mihalas and Weibel–Mihalas [20] and investigated more recently in astrophysics and laser applications (in the inviscid case) by Lowrie, Morel and Hittinger [16] and Buet and Desprès [3], with a special attention to asymptotic regimes.

In the remaining part of this section we suitably rescale our system of equations and formulate the primitive and the target system. Section 2 contains definition of the weak solution to our system. Section 3 deals with the existence of solutions to the target system. In Section 4 we present the relative entropy inequality and state the convergence result for our thin disk model. Last Section 5 contains the proof of the convergence result.

1.1 Formal scaling analysis, primitive system and target system

We rescale our problem to a fixed domain. To this aim, we introduce

[TABLE]

however, keep the notation $\varrho$ for the density, ${\vec{u}}$ for the velocity and $I$ for the radiative intensity. We further denote

[TABLE]

[TABLE]

[TABLE]

Moreover, in order to identify the appropriate limit regime, we perform a general scaling. Since we are only interested in the behaviour of the Froude number, we set all other non-dimensional numbers immediately equal to one.

The continuity equation reads now

[TABLE]

the momentum equation is

[TABLE]

[TABLE]

and the transport equation has the form

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and recall

[TABLE]

We denote (cf. (1.13))

[TABLE]

We now distinguish two cases with respect to the behaviour of the Froude number, namely $Fr\sim 1$ and $Fr\sim\sqrt{\epsilon}$. In order to avoid technicalities, we directly consider either $Fr=\sqrt{\epsilon}$ or $Fr=1$. Furthermore, according to the choice of the Froude number, we have to consider the correct form of the gravitational potential, namely in the former the self-gravitation and in the latter the external gravitation force. In the latter, we could also include the self-gravitation, it would, however, disappear after the limit passage $\epsilon\to 0^{+}$.

Supposing $Fr=\sqrt{\epsilon}$ and $\eta=1$, we get the primitive system

[TABLE]

[TABLE]

[TABLE]

Next, taking $Fr=1$ and $\eta=0$, the primitive system reads

[TABLE]

[TABLE]

[TABLE]

Our goal is to investigate the limit process $\epsilon\to 0^{+}$ in the systems of equations (1.26–1.28) and (1.29–1.31), respectively, under the assumptions that initial data $[\varrho_{0,\epsilon},{\vec{u}}_{0,\epsilon},I_{0,\epsilon}]$ converge in a certain sense to $[r_{0},{\vec{V}}_{0},J_{0}]=[r_{0,h},({\vec{w}}_{0,h},0),J_{0,h}]$.

Let us return back to the former, i.e. $Fr=\sqrt{\epsilon}$ and $\eta=1$. As the target system does not depend on the vertical variable $x_{3}$, we expect that the sequence $[\varrho_{\epsilon},{\vec{u}}_{\epsilon},I_{\epsilon}]$ of weak solutions to (1.26–1.28) will converge to $[r,\vec{V},J]$ for $\vec{V}=[\vec{w},0]$, where $\vec{w}=(w_{1},w_{2})$ and the triple $[r(t,{\vec{x}}_{h}),{{\vec{w}}}(t,{\vec{x}}_{h}),J(t,{\vec{x}}_{h},{\vec{\varsigma}},\nu)]$ solves the following 2-D rotating Navier–Stokes–Poisson system with radiation in the domain $(0,T)\times\omega$

[TABLE]

[TABLE]

[TABLE]

with the formula

[TABLE]

where

[TABLE]

Above, $\mbox{\F I}_{h}$ is the unit tensor in $\mbox{\F R}^{2\times 2}$,

[TABLE]

and

[TABLE]

When $Fr=1$, we also expect that the sequence $[\varrho_{\epsilon}(t,{\vec{x}}),{\vec{u}}_{\epsilon}(t,{\vec{x}}),I_{\epsilon}(t,{\vec{x}},{\vec{\varsigma}},\nu)]$ of weak solutions to (1.29–1.31) will converge to $[r,{\vec{V}},J]$, where the velocity vector $\vec{V}$ is as above, $[r(t,{\vec{x}}_{h}),{{\vec{w}}}(t,{\vec{x}}_{h}),J(t,{\vec{x}}_{h},{\vec{\varsigma}},\nu)]$ solves now the 2-D rotating Navier–Stokes system with radiation and external gravitational force

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

Observe that, through formula (1.35), the gravitational contribution in the target momentum equation for $Fr=\sqrt{\epsilon}$ is the tangential gradient of a single layer potential which actually is different from the analogous quantity deriving from the solution of the 2-D Poisson equation $-\Delta_{h}\phi=Gr$, which would lead to the well-known logarithmic expression.

Finally we check, as stressed by Maltese and Novotný [17], that the bulk viscosity coefficient is modified in the limit (compare (1.36) with (1.12)).

Our aim is now to prove that solutions of (1.26–1.28) and (1.29–1.31) converge in a certain sense (to be precised) to the unique solution of (1.32–1.34) and (1.38–1.40), respectively.

Note also that considering the boundary conditions of the type (1.17) on the whole boundary of $\Omega$ we would get in the limit that $J\equiv 0$. Our method would yield that the solutions to (1.26–1.28) and (1.29–1.31), respectively, would converge to the same system as above, however, without the radiation.

2 Weak solutions of the primitive system

We consider the rescaled problems (1.26–1.28) and (1.29–1.31), respectively, with boundary conditions

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

We define the adapted functional space

[TABLE]

In the weak formulation of the Navier–Stokes–Poisson system, equation of continuity (1.26) is replaced by its weak version

[TABLE]

satisfied for all $\tau\in(0,T]$ and any test function $\varphi\in C^{\infty}_{c}([0,T)\times\overline{\Omega})$.

Similarly, the momentum equation (1.27) is replaced by

[TABLE]

[TABLE]

[TABLE]

for any $\vec{\varphi}\in C^{\infty}_{c}([0,T)\times\overline{\Omega};\mbox{\F R}^{3})$ such that ${\vec{\varphi}}|_{[0,T]\times\partial\omega\times\{0,1\}}={\vec{0}}$ and $\varphi_{3}|_{[0,T]\times\omega\times\{0,1\}}=0$. Above, $j=1$ if $\eta=1$ (i.e. $Fr=\sqrt{\epsilon}$) and $j=2$ if $\eta=0$ (i.e. $Fr=1$). The radiative transport equation is satisfied in the following sense

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for all $\varphi\in C^{\infty}([0,T]\times\overline{\Omega})$, where {\vec{N}}=\big{(}n_{1},n_{2},\frac{1}{\epsilon}n_{3}\big{)} with ${\vec{n}}$ the external normal to $\Omega$. Note that ${\vec{\varsigma}}\cdot{\vec{n}}=\pm\varsigma_{3}$ on $\omega\times\{0,1\}$.

Moreover, denoting

[TABLE]

and

[TABLE]

the energy inequality

[TABLE]

[TABLE]

[TABLE]

holds for a.e. $\tau\in(0,T)$, where $j=1,2$ as above. Its validity is closely connected to the following result.

Lemma 2.1

[Darrozes–Guiraud] Under our assumptions, we have for a.a. $(t,\vec{x})\in(0,T)\times\partial\Omega$

[TABLE]

The proof of the lemma can be found in [1].

We are now in position to define weak solutions of our primitive system.

Definition 2.1

We say that $\varrho,{\vec{u}},I$ is a weak solution of problem (1.26–1.28) and (1.29–1.31), respectively, if

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and if $\varrho$, ${\vec{u}}$, $I$ satisfy the integral identities (2.5), (2.6), (2.7) together with the total energy inequality (2.10) and the integral representation of the gravitational force (1.22).

We have the following existence result for the primitive system

Proposition 2.1

Assume that $\omega\subset\mbox{\F R}^{2}$ is a domain with compact boundary of class $C^{2+\nu}$, $\nu>0$. Suppose that the stress tensor is given by (1.12) and $p$ verifies (1.11), the boundary conditions are given by (2.1–2.4) and the initial data satisfy the conditions

[TABLE]

[TABLE]

[TABLE]

Let $\gamma>3/2$ if $\eta=0$ or $\gamma>\frac{12}{7}$ if $\eta=1$ and let the external force $g\in L^{p}(\mbox{\F R}^{3})$ for $p=1$ if $\gamma>6$ and $p=\frac{6\gamma}{7\gamma-6}$ for $\frac{3}{2}<\gamma\leq 6$.

Then problems (1.26–1.28) and (1.29–1.31), respectively, admit at least one finite energy weak solution according to Definition 2.1.

More details can be found in [11]. Note that the different boundary conditions for the radiation intensity do not cause any troubles due to Lemma 2.1. Using this result, in fact, the existence of the solution can be shown using the approach given in [9] when $\chi=0$ (non rotating case) and for no slip condition on $\partial\Omega$. It is first easy to see that the centrifugal term can be treated in the same way as the gravitational term in [9] and that the Coriolis term may be absorbed in the energy by a Gronwall argument. Finally the slip conditions on top and bottom of the domain may be accommodated using the argument of Vodák [28].

3 Strong solution of the target system

We consider our target system (1.32–1.35) and (1.38–1.41), respectively, with the boundary conditions

[TABLE]

and

[TABLE]

Let $(\overline{r},\vec{0},\overline{J})$ be a given constant state with $\overline{r}>0$, and $\overline{J}=B(\nu,\overline{r})$. We denote

[TABLE]

where $\mbox{\F V}_{0}$ is the initial vorticity (recall that $V^{0}_{ij}=\partial_{j}w_{i}^{0}-\partial_{i}w_{j}^{0}$),

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

for an arbitrary fixed $\alpha$ such that $3<\alpha<6$.

The following result holds

Proposition 3.1

Let $p\in C^{2}(0,\infty)$.

Let $(r^{0},{\vec{w}}^{0},J^{0})\in H^{3}(\omega;\mbox{\F R}^{4})$, $\inf_{\omega}J^{0}>0$, $\inf_{\omega}r^{0}>0$ and assume the following compatibility condition

[TABLE]

holds, where $\phi_{h}=\phi$ (see (1.35)) for $\eta=1$ and $\phi_{h}=\widetilde{\phi}$ (see (1.41)) for $\eta=0$ and $V^{0}=({\vec{w}}^{0},0)$.

There exist positive constants $\delta\leq 1$ and $\Gamma>0$ depending on the data such that if $E_{0}\leq\Gamma\delta$, the triple $(r,\vec{w},J)$ is the unique classical solution to the Navier–Stokes–Poisson system with radiation (1.32–1.35) and (1.38–1.41) in $(0,T)\times\omega$ for any $T>0$ such that

[TABLE]

[TABLE]

[TABLE]

Moreover, there exists $\delta_{1}>0$ such that if $e_{0}\leq\delta_{1}$, then

[TABLE]

and

[TABLE]

[TABLE]

The proof of the Proposition 3.1 follows from [10] and [18].

Remark 3.1

In fact this solution $(r,\vec{w},J)$ can be defined in the whole domain $\Omega=\omega\times(0,1)$ by the triple $(r,\vec{V},J)$, where $\vec{V}=(\vec{w},0)$ and all quantities are constant in $x_{3}$.

Another possible strong solution can be constructed on short time intervals when no restriction on the size of the initial data is imposed. The result reads

Proposition 3.2

Let $p\in C^{2}(0,\infty)$. Let $(r^{0},{\vec{w}}^{0},J^{0})\in H^{3}(\omega;\mbox{\F R}^{4})$, $\inf_{\omega}J^{0}>0$, $\inf_{\omega}r^{0}>0$ and assume the following compatibility condition

[TABLE]

holds, where $\phi_{h}=\phi$ (see (1.35)) for $\eta=1$ and $\phi_{h}=\widetilde{\phi}$ (see (1.41)) for $\eta=0$ and $V^{0}=({\vec{w}}^{0},0)$.

There exist positive constant $T_{*}$ depending on the data such that on $(0,T_{*})$, there exists triple $(r,\vec{w},J)$, the unique classical solution to the Navier–Stokes–Poisson system with radiation (1.32–1.35) and (1.38–1.41) such that

[TABLE]

[TABLE]

Proof of Proposition 3.2 can be deduced from [10]. See [14] or [26] for a similar type of results.

4 Relative entropy inequality

Let us introduce, in the spirit of [17], a relative entropy inequality which is satisfied by any weak solution $(\varrho,{\vec{u}},I)$ of the rotating Navier–Stokes–Poisson system (1.1–1.6).

We define the relative entropy functional

[TABLE]

with

[TABLE]

where $(r,\vec{V},J)$ is a triple of ”arbitrary” smooth enough functions where only $r,\vec{V}$ are arbitrary and $J$ satisfies the transport equation for $J$ with the boundary condition (2.4) on $\omega\times\{0,1\}$ and (2.3) on $\partial\omega\times(0,1)$, or $J$ fulfills (2.3) on $\partial\omega\times(0,1)$ and is independent of $x_{3}$. Note that the latter case is exactly that we will need later.

Then we have

Lemma 4.1

Let all assumptions of Proposition 2.1 be satisfied and $I_{0}\in L^{2}(\Omega\times S^{2}\times(0,\infty))$. Let $Fr=\sqrt{\epsilon}$ and let $(\varrho,{\vec{u}},I)$ be a finite energy weak solution of system (1.26–1.28) (then $j=1$) or $Fr=1$ and let $(\varrho,{\vec{u}},I)$ be a finite energy weak solution of system (1.29–1.31) (then $j=2$) in the sense of Definition 2.1.

Then $(\varrho,{\vec{u}},I)$ satisfies the relative entropy inequality

[TABLE]

[TABLE]

where the remainder ${\mathcal{R}}$ is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for any triple $(r,\vec{V},J)$ of test functions such that

[TABLE]

and either

[TABLE]

*or ${\vec{V}}$ is independent of $x_{3}$, and $J$ satisfies the transport equation for $J$ with the boundary condition (2.4) at $\omega\times\{0,1\}$ and (2.3) at $\partial\omega\times(0,1)$, or $J$ fulfills (2.3) at $\partial\omega\times(0,1)$ and is independent of $x_{3}$. *

Proof:  Using $\varphi=\frac{1}{2}\ |\vec{V}|^{2}$ as test function in (2.5) we get

[TABLE]

Using $\varphi=-\vec{V}$ as test function in (2.6) yields

[TABLE]

[TABLE]

[TABLE]

Above, $j=1$ for $\eta=1$ and $j=2$ for $\eta=0$. Using $\varphi=-H^{\prime}(r)$ as test function in (2.5) leads to

[TABLE]

[TABLE]

Note that $rH^{\prime}(r)=H(r)+p(r)$, therefore $r\partial_{t}H^{\prime}(r)=\partial_{t}p(r)$ and $r\nabla_{\epsilon}H^{\prime}(r)=\nabla_{\epsilon}p(r)$. Employing these identities gives

[TABLE]

[TABLE]

Taking difference between the weak formulation of the transport equation for $I$ and $J$, using as test function $\varphi=I-J$ yields

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\vec{N}=(n_{1},n_{2},\frac{1}{\epsilon}n_{3})$. Adding (4.4–4.7) and (2.10) (without the part connected to the radiative transfer equation) and recalling that the boundary integrals in (4.7) are non-negative, we end up with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which yields (4.3).

4.1 Convergence result

We aim at proving the following result.

Theorem 4.1

Suppose that the pressure $p$ satisfy hypothesis (1.11), and that the stress tensor is given by (1.12).

Let $r_{0},{\vec{w}}_{0},J_{0}$ satisfy assumptions of Proposition 3.1 or 3.2 and let $T_{*}>0$ be the time interval of existence of the strong solution to the problem (1.32–1.35) or (1.38–1.41), respectively, corresponding to $r_{0},{\vec{w}}_{0},J_{0}$.

In addition to hypotheses of Proposition 2.1, we suppose that $I_{0}\in L^{2}(\Omega\times S^{2}\times(0,\infty))$, (1.10) and

• •

either $Fr=1$, $\eta=0$, $\gamma>\frac{3}{2}$ and $g\in L^{p}(\mbox{\F R}^{3})$ with $p=1$ for $\gamma>6$ and $p=\frac{6\gamma}{7\gamma-6}$ for $\gamma\in(\frac{3}{2},6]$, and

[TABLE]

for all ${\vec{x}}_{h}\in\omega$

• •

or $Fr=\sqrt{\epsilon}$, $\eta=1$ and $\gamma\geq\frac{12}{5}$.

Let $(\varrho_{\epsilon},{\vec{u}}_{\epsilon},I_{\epsilon})$ be a sequence of weak solutions to the 3D compressible Navier-Stokes-Poisson system with radiation (1.26–1.28) or (1.29–1.31)) with (2.1–2.4)) emanating from the initial data $\varrho_{0},{\vec{u}}_{0},I_{0}$.

Suppose that

[TABLE]

where ${\vec{V}}_{0}=[{\vec{w}}_{0},0]$ and all quantities are extended constantly in $x_{3}$ to $\Omega$.

Then

[TABLE]

[TABLE]

and the triple $(r,{\vec{w}},J)$ restricted to $\omega$ satisfies the 2D rotating Navier–Stokes–Poisson system with radiation (1.32–1.35) or (1.38–1.41), respectively, with the boundary condition (3.1–3.2) on the time interval $[0,T_{*}]$.

Remark 4.1

From (4.9) it follows in addition to (4.10)

[TABLE]

and

[TABLE]

Corollary 4.1

Suppose that the pressure $p$ satisfy hypothesis (1.11), and that the stress tensor is given by (1.12).

Assume that $[\varrho_{\epsilon,0},{\vec{u}}_{\epsilon,0},I_{\epsilon,0}]$, $\varrho_{\epsilon,0}\geq 0$ satisfy

[TABLE]

[TABLE]

[TABLE]

where $r_{0},{\vec{w}}_{0},J_{0}$ belong to the regularity class of Propositions 3.1 and 3.2, and

[TABLE]

Let $[\varrho_{\epsilon},{\vec{u}}_{\epsilon},I_{\epsilon}]$ be a sequence of weak solutions to the 3-D compressible Navier–Stokes–Poisson system with radiation (1.1–1.6) emanating from the initial data $[\varrho_{\epsilon,0},{\vec{u}}_{\epsilon,0},I_{\epsilon,0}]$.

Then properties (4.9–4.10) hold.

5 Proof of Theorem 4.1

5.1 Preliminaries

We can easily verify that

[TABLE]

for any ${\vec{u}}\in W^{1,2}(\Omega;\mbox{\F R}^{3})$. However, for any ${\vec{u}}\in W^{1,2}_{0,\vec{n}}(\Omega;\mbox{\F R}^{3})$ we have

[TABLE]

Thus for any ${\vec{u}}\in W^{1,2}_{0,\vec{n}}(\Omega;\mbox{\F R}^{3})$

[TABLE]

provided $\mu>0$ and $\xi\geq 0$. Moreover, under same assumptions

[TABLE]

for any ${\vec{w}}\in W^{1,2}_{0}(\omega;\mbox{\F R}^{2})$.

Moreover, note that we also have the Poincaré inequality in the form

[TABLE]

for any ${\vec{w}}\in W^{1,2}_{0}(\omega;\mbox{\F R}^{2})$.

Due to the energy equality (2.10) and Korn’s inequality (5.2) above, we have the following bounds for the sequence $(\varrho_{\epsilon},{\vec{u}}_{\epsilon},I_{\epsilon})$

[TABLE]

[TABLE]

with the constant $C$ independent of $\epsilon$. These estimates hold if $\gamma\geq\frac{12}{5}$ (if $\eta=1$) or under the assumptions on $g$ from Theorem 4.1 (if $\eta=0$), for any $\gamma>\frac{3}{2}$. Note that the limit on $\gamma$ comes from the gravitational potential, as

[TABLE]

for $1<p<\infty$, with $C$ independent of $\epsilon$. Thus

[TABLE]

[TABLE]

if $\gamma\geq\frac{12}{5}$. On the other hand,

[TABLE]

[TABLE]

with $p$ from Theorem 4.1, as

[TABLE]

where we used the embedding $W^{1,p}\hookrightarrow L^{\frac{6\gamma}{5\gamma-6}}$.

Moreover, we may deduce the following estimate for the radiative intensity. Assuming that $I_{0}$ belongs to $L^{2}(\Omega\times{\mathcal{S}}^{2}\times\mbox{\F R}^{+})$, multiplying (1.6) by $I$ we get

[TABLE]

Consequently, denoting

[TABLE]

we deduce, after integrating integrating the above expression and using Lemma 2.1, that

[TABLE]

[TABLE]

[TABLE]

We now recall the necessary definitions of essential and residual sets.

5.2 Essential and residual sets

For two numbers $0<\underline{\varrho}\leq\overline{\varrho}<\infty$, the essential and residual subsets of $\Omega$ are defined for a.e. $t\in(0,T)$ as follows:

[TABLE]

For any function $h$ defined for a.e. $(t,x)\in(0,T)\times\Omega$, we write

[TABLE]

In the sequel we will choose $\underline{\varrho}=\inf_{(0,T)\times\Omega}r$ and $\overline{\varrho}=\sup_{(0,T)\times\Omega}r$.

From [17] we have

Lemma 5.1

Let $0<a<b<\infty$. There exists a constant $C=C(a,b)>0$ such that for all $\varrho\in[0,\infty)$ and $r\in[a,b]$

[TABLE]

where

[TABLE]

and $\underline{\varrho}=a$, $\overline{\varrho}=b$ in the definition of the essential set.

A consequence of this result is the lower bound

[TABLE]

[TABLE]

5.3 Estimates of the remainder

In what follows, we plan to use in the relative entropy inequality (4.2) as “smooth test functions” the solution to the 2-D rotating Navier–Stokes–Poisson system with rotation constructed in Section 3. To this aim, we slightly rearrange the terms in the remainder (4.3) in order to be able to use the validity of the 2-D system. However, we keep writing all the integrals over $\Omega$ and assume for a moment that all functions which are independent of $x_{3}$ are constant in this variable, and the third velocity component is zero. Denoting by $(\varrho_{\epsilon},{\vec{u}}_{\epsilon},I_{\epsilon})$ the solution of the primitive system we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In what follows, we will estimates the terms $R_{j}$ for $j=1,2,\dots,8$.

5.3.1 Estimate of $R_{1}$

We have

[TABLE]

[TABLE]

Recall that we used estimate (5.9), that fact that ${\vec{V}}=({\vec{w}},0)$ and note that due to Section 3 we know that

[TABLE]

5.3.2 Estimate of $R_{2}$

We first consider the part of the integral over the essential set and use again estimate (5.9) from Lemma 5.1.

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

For the residual part we consider separately the regions $\{\varrho<\underline{\varrho}/2\}$ and $\{\varrho>2\overline{\varrho}\}$. Then

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that we used $\int_{\Omega}\varrho_{\epsilon}\ {\rm d}{\vec{x}}=\int_{\Omega}\varrho_{0,\epsilon}\ {\rm d}{\vec{x}}\leq C$ independently of $\epsilon$.

5.3.3 Estimate of $R_{3}$

We use the fact that $(r,{\vec{w}},J)$ solves the target system. Therefore we have

[TABLE]

in the case when $Fr=\sqrt{\epsilon}$, and

[TABLE]

in the case when $Fr=1$. We will use this fact in the treatment of the term $R_{6}$.

5.3.4 Estimate of $R_{4}$

Since

[TABLE]

we have

[TABLE]

Therefore

[TABLE]

In order to estimate the first term, we use a similar approach as in the estimate of $R_{2}$. We divide the integral into three parts: over the essential set, the set where $\varrho_{\epsilon}<\underline{\varrho}/2$ and where $\varrho_{\epsilon}>2\overline{\varrho}$. Then

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and, finally,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Similarly, we will deal with $R_{4}^{2}$. Using Taylor formula and the regularity of the pressure, and dividing the integral over the essential and residual sets, we have

[TABLE]

[TABLE]

with

[TABLE]

Using the bound

[TABLE]

we can estimate

[TABLE]

[TABLE]

5.3.5 Estimate of $R_{5}$

We write

[TABLE]

[TABLE]

[TABLE]

Easily, as in the estimate of $R_{1}$ and $R_{2}$, we have

[TABLE]

with

[TABLE]

Due to (1.8) we have

[TABLE]

[TABLE]

Similarly, using also (1.10), we get

[TABLE]

with

[TABLE]

5.3.6 Estimate of $R_{6}$

For the gravitational potential, we have to consider both cases separately. We start with the simpler one. i.e. with the case $Fr=1$. Here, only the gravitational effect of other objects than the fluid itself is considered. Recall that

[TABLE]

We combine the term $R_{3}$ with $R_{6}$. Therefore we have to verify

[TABLE]

[TABLE]

First note that due to our assumption on the integrability of $g$ and proceeding similarly as in the estimate of $R_{2}$ (replacing the $L^{2}$ estimate of ${\vec{V}}-{\vec{u}}_{\epsilon}$ by the $L^{6}$ estimate) is is enough to verify that

[TABLE]

[TABLE]

for a.a. $\tau\in(0,T_{*})$. Moreover, it is not difficult to verify that (note that to get estimates independent of $\epsilon$ of the integral over $\mbox{\F R}^{3}$ is easy) it remains to verify

[TABLE]

for all ${\vec{x}}_{h}\in\omega$, $x_{3}\in(0,1)$, $t\in(0,T_{*})$ and $\overline{g}\in C^{\infty}_{c}(\mbox{\F R}^{3})$. As

[TABLE]

for a.a. $(x_{h},x_{3})\in\Omega$, $(y_{h},y_{3})\in\mbox{\F R}^{3}$, $\tau\in(0,T_{*})$, and

[TABLE]

the Lebesgue dominated convergence theorem yields the required identity (5.11).

The case of the self-gravitation ($Fr=\sqrt{\epsilon}$) is more complex. Here, we have to show that

[TABLE]

[TABLE]

where $H_{\epsilon}=o(\epsilon)$ as $\epsilon\to 0^{+}$. The derivative of the integral over $\omega$ with respect to $x_{3}$ is indeed zero. First of all, for $\gamma\geq\frac{12}{5}$, as in (5.5), using the decomposition to the essential and the residual set and proceeding as in the estimates of the remainder above, we can show that it is enough to verify that

[TABLE]

for a.a. $\tau\in(0,T_{*})$.

Using the change of the variables to integrate over $\Omega_{\epsilon}$ it is enough to show

[TABLE]

Note that

[TABLE]

where v.p. means the integral in the principal value sense. Thus it remains to show

[TABLE]

and

[TABLE]

We fix ${\vec{x}}_{0}\in\omega$, $\Delta>0$, sufficiently small, and denote $B_{\Delta}({\vec{x}}_{0})=\{{\vec{x}}\in\omega;|{\vec{x}}-{\vec{x}}_{0}|<\Delta\}$ and $C_{\Delta}({\vec{x}}_{0})=\{{\vec{x}}\in\Omega;|{\vec{x}}_{h}-{\vec{x}}_{0}|<\Delta,0<x_{3}<1\}$.

We first consider (5.13). Fix $\delta>0$. Using the change of variables (from $\Omega$ back to $\Omega_{\epsilon}$) it is not difficult to see that there exists $\Delta>0$ such that for any $0<\epsilon\leq 1$, $0<x_{3}<1$ it holds

[TABLE]

and for this $\Delta>0$ there exists $\epsilon_{0}>0$ such that we have for any $0<\epsilon\leq\epsilon_{0}$

[TABLE]

which yields (5.13).

In order to verify (5.14), we proceed similarly. Since $\frac{{\vec{x}}_{h}-{\vec{y}}_{h}}{|{\vec{x}}_{h}-{\vec{y}}_{h}|^{3}}$ is a singular integral kernel in the sense of Calderón–Zygmund, for a fixed ${\vec{x}}_{0}\in\omega$, $0<x_{3}<1$ and $\delta>0$ there exists $\Delta>0$ such that

[TABLE]

and

[TABLE]

We fix such $\Delta>0$. Using that

[TABLE]

for any ${\vec{y}}_{h}\in\omega$, $0<x_{3},y_{3}<1$, except ${\vec{x}}_{0}={\vec{y}}_{h}$, we see that for the above fixed $\Delta>0$ there exists $\epsilon_{0}>0$ such that for any $0<\epsilon\leq\epsilon_{0}$

[TABLE]

hence we get (5.14).

5.3.7 Estimate of $R_{7}$ and $R_{8}$

Repeating the arguments from the estimate of $R_{5}^{4}$ and $R^{5}_{5}$, using (1.8–1.10) (in particular, the Lipschitz continuity of $B$, $\sigma_{a}$ and $\sigma_{s}$ in the density) we easily verify that

[TABLE]

5.3.8 Conclusion

Collecting all of the previous estimates, plugging them into the relative entropy inequality and taking $\delta$ small enough, we end with the inequality

[TABLE]

where $K\in L^{1}(0,T)$ and

[TABLE]

where $H_{\epsilon}(\tau)\to 0$ for $\epsilon\to 0$ in $L^{1}(0,T_{*})$. Hence, it implies by virtue of Gronwall’s lemma

[TABLE]

for a.a. $\tau\in[0,T]$, which establishes (4.9). Returning back to the relative entropy inequality (4.2), we verify (4.10) which finishes the proof of Theorem 4.1.

Acknowledgements: B. D. is partially supported by the ANR project INFAMIE (ANR-15-CE40-0011) Š. N. is supported by the Czech Science Foundation, grant No. 201-16-03230S and by RVO 67985840. Part of this paper was written during her stay in CEA and she would like to thank Prof. Ducomet for his hospitality during her stay. M. P. was supported by the Czech Science Foundation, grant No. 201-16-03230S. M.A.R.B. was partially supported by MINECO grant MTM2015-69875-P (Ministerio de Economía y Competitividad, Spain) with the participation of FEDER. She would also like to thank Prof. Nečasová for her hospitality during the stay in Prague.