On infinite order differential operators in fractional viscoelasticity

TL;DR

This paper explores the properties of viscoelastic models characterized by infinite derivatives, focusing on Bessel models that mimic fractional Maxwell behavior at short times, advancing understanding of complex material responses.

Contribution

It introduces and analyzes properties of viscoelastic models with infinite derivatives, specifically examining Bessel models and their fractional Maxwell-like behavior.

Findings

Bessel models exhibit fractional Maxwell behavior at short times

Models involve both integer and fractional order derivatives

Provides insights into complex viscoelastic responses

Abstract

In this paper we discuss some general properties of viscoelastic models defined in terms of constitutive equations involving infinitely many derivatives (of integer and fractional order). In particular, we consider as a working example the recently developed Bessel models of linear viscoelasticiy that, for short times, behave like fractional Maxwell bodies of order $1/2$.

Equations

$$\epsilon (t)=J(0+)\sigma (t)+(\dot{J}\ast \sigma )(t),\sigma (t)=G(0+)\epsilon (t)+(\dot{G}\ast \epsilon )(t),$$

$$P(D)=\sum _{k=0}^N{p}_k{D}^k,Q(D)=\sum _{k=0}^N{q}_k{D}^k,$$

$$P(D)\sigma (t)=Q(D)\epsilon (t).$$

$$\sum _{r=k}^N{p}_r{\sigma }^{(r-k)}(0+)=\sum _{r=k}^N{q}_r{\epsilon }^{(r-k)}(0+),k=1,2\dots ,N,$$

$$P(s)\tilde{\sigma }(s)=Q(s)\tilde{\epsilon }(s),\forall s\in \mathbb{C}:\text{Re}(s)>\alpha ,$$

$$\Phi \left(D\right)f(t):=\sum _{n=0}^{\infty }{a}_n{D}^{n}f(t).$$

$$P_{\nu }(D_{\ast })\sigma (t)=Q_{\nu }(D_{\ast })\epsilon (t),\nu >-1,$$

$$P_{\nu }(z)=I_{\nu }(\sqrt{z}),Q_{\nu }(z)=I_{\nu +2}(\sqrt{z}),$$

$$\begin{array}{rl}& \Theta (\lceil \nu /2\rceil -k)\sum _{n=0}^{\infty }{p}_n{\sigma }^{(\lceil \nu /2\rceil +n-k)}(0+)+\\ & \Theta (k-\lceil \nu /2\rceil -1)\sum _{n=k-\lceil \nu /2\rceil }^{\infty }{p}_n{\sigma }^{(\lceil \nu /2\rceil +n-k)}(0+)=\\ & =\Theta (\lceil \nu /2\rceil +1-k)\sum _{n=0}^{\infty }{q}_n{\epsilon }^{(\lceil \nu /2\rceil +n+1-k)}(0+)+\\ & \Theta (k-\lceil \nu /2\rceil -2)\sum _{n=k-\lceil \nu /2\rceil -1}^{\infty }{q}_n{\epsilon }^{(\lceil \nu /2\rceil +n+1-k)}(0+),\end{array}$$

$$p_k=\frac{1}{k!\Gamma \left(\nu +k+1\right)4^{k+\frac{\nu }{2}}},q_k=\frac{1}{k!\Gamma \left(\nu +k+3\right)4^{k+1+\frac{\nu }{2}}},$$

$$I_{\nu }(\sqrt{s})\tilde{\sigma }(s)=I_{\nu +2}(\sqrt{s})\tilde{\epsilon }(s).$$

$$s\tilde{J}(s;\nu )=\frac{I_{\nu }(\sqrt{s})}{I_{\nu +2}(\sqrt{s})}.$$

$$\begin{array}{rl}J(t;\nu )=& 2\left(\frac{\nu +2}{\nu +3}\right)+4(\nu +1)(\nu +2)t\\ & -4(\nu +1)\sum _{n=1}^{\infty }\frac{1}{j_{\nu +2,n}^2}\exp \left(-j_{\nu +2,n}^{2}t\right),\end{array}$$

$${\overline{P}}_{\nu }(s)\tilde{\sigma }(s)={\overline{Q}}_{\nu }(s)\tilde{\epsilon }(s),$$

$${\overline{P}}_{\nu }(s)={\left(\frac{4}{s}\right)}^{\frac{\nu }{2}}{I}_{\nu }(\sqrt{s})=\sum _{k=0}^{\infty }\frac{1}{k!\Gamma \left(\nu +k+1\right)}{\left(\frac{s}{4}\right)}^k,$$

$${\overline{Q}}_{\nu }(s)={\left(\frac{4}{s}\right)}^{\frac{\nu }{2}}{I}_{\nu +2}(\sqrt{s})=\sum _{k=0}^{\infty }\frac{1}{k!\Gamma \left(\nu +k+3\right)}{\left(\frac{s}{4}\right)}^{k+1}.$$

$${\overline{P}}_{\nu }(D)\sigma (t)={\overline{Q}}_{\nu }(D)\epsilon (t),\nu >-1,$$

$$P_{\nu }^{(N)}(z)$$

$$Q_{\nu }^{(N)}$$

$$P_{\nu }(D_{\ast })\sigma (t)=\underset{N\to \infty }{\lim }{P}_{\nu }^{(N)}(D_{\ast })\sigma (t),Q_{\nu }(D_{\ast })\epsilon (t)=\underset{N\to \infty }{\lim }{Q}_{\nu }^{(N)}(D_{\ast })\epsilon (t).$$

$$\left|P_{\nu }^{(N)}(D_{\ast })\sigma (t)\right|\le \alpha (t),\left|Q_{\nu }^{(N)}(D_{\ast })\epsilon (t)\right|\le \beta (t),$$

$$\mathcal{L}\left\{P_{\nu }(D_{\ast })\sigma (t);s\right\}$$

$$\mathcal{L}\left\{Q_{\nu }(D_{\ast })\epsilon (t);s\right\}$$

$$\mathcal{L}\left\{D_{\ast }^{\alpha }f(t);s\right\}=s^{\alpha }\tilde{f}(s)-\sum _{k=0}^{\lceil \alpha \rceil -1}{f}^{(k)}(0+)s^{\alpha -k-1},\alpha \in {\mathbb{R}}^{+},$$

$$\mathcal{L}\left\{Q_{\nu }^{(N)}(D_{\ast })\epsilon (t);s\right\}=Q_{\nu }^{(N)}(s)\tilde{\epsilon }(s)-\sum _{n=0}^N{q}_n\sum _{k=0}^{\lceil \nu /2\rceil +n}{\epsilon }^{(k)}(0+)s^{(\nu /2)+n-k}.$$

$$\sum _{k=0}^{\lceil \alpha \rceil -1}{f}^{(k)}(0+)s^{\alpha -k-1}=s^{\{\alpha \}-{\delta }_{\alpha }}\sum _{r=1}^{\lceil \alpha \rceil }{f}^{(\lceil \alpha \rceil -r)}(0+)s^{r-1},$$

$${\delta }_{\alpha }=\{\begin{array}{rl}& 0,\alpha \in \mathbb{Z},\\ & 1,\alpha \notin \mathbb{Z}.\end{array}$$

$$\sum _{k=0}^{\lceil \alpha \rceil -1}{f}^{(k)}(0+)s^{\alpha -k-1}=s^{\{\alpha \}-{\delta }_{\alpha }}\sum _{k=0}^{\lceil \alpha \rceil -1}{f}^{(k)}(0+)s^{\lceil \alpha \rceil -k-1},$$

$$\sum _{k=0}^{\lceil \alpha \rceil -1}{f}^{(k)}(0+)s^{\lceil \alpha \rceil -k-1}=\sum _{r=1}^{\lceil \alpha \rceil }{f}^{(\lceil \alpha \rceil -r)}(0+)s^{r-1},$$

$$\begin{array}{rl}\mathcal{L}\left\{Q_{\nu }^{(N)}(D_{\ast })\epsilon (t);s\right\}& =Q_{\nu }^{(N)}(s)\tilde{\epsilon }(s)-\frac{s^{\{\nu /2\}-{\delta }_{\nu /2}}}{s}\times \\ & \times \sum _{n=0}^N{q}_n\sum _{r=1}^{\lceil \nu /2\rceil +n+1}{s}^r{\epsilon }^{(\lceil \nu /2\rceil +n+1-r)}(0+).\end{array}$$

$$\displaystyle\begin{split}\sum_{n=0}^{N}q_{n}\sum_{r=1}^{M+n}s^{r}\,b_{M+n-r}=\sum_{k=1}^{N+M}s^{k}\,\Big{[}&\Theta(M-k)\sum_{h=0}^{N}q_{h}\,b_{M+h-k}+\\ &\Theta(k-M-1)\sum_{h=k-M}^{N}q_{h}\,b_{M+h-k}\Big{]}\,,\end{split}$$

Full text

In Fract. Calc. Appl. Anal., Vol 20, No 4 (2017), pp. 854–867, DOI: 10.1515/fca-2017-0045.

1. Introduction

In recent years, the quest for potential biomedical applications of ordinary and fractional viscoelastic models has attracted much attention, see e.g. [1; 4; 8; 9; 14; 17; 18; 19; 20; 21]. In particular, a generalization of the viscoelastic model introduced by Giusti and Mainardi in [9] has lead to the Bessel models [4].

First of all, it is worth recalling some generalities on linear viscoelasticity. As widely discussed in the literature (see e.g. [3; 7; 11; 15; 16]), a linear viscoelastic body is defined as a linear system in which the stress plays the role of excitation function for a certain material, while the strain act as the response function, or vice versa. The stress-strain relation, also known as constitutive relation, for a given material can then be expressed mathematically in two different forms, either as an integral equation or as a differential equation.

First, we denote with $\mathcal{H}^{N}$ the Heaviside class of functions (i.e. set of causal functions $f(t)$ such that $f\in C^{N}([0,+\infty[)$, with $N\in\mathbb{N}$, and with $\sigma,\varepsilon\in\mathcal{H}^{N}$ the uniaxial stress and strain functions for a certain system, respectively.

Given these general assumptions, let us start from the integral form. Then, under the hypothesis of sufficiently well behaved causal histories, the constitutive equations can be written in the following forms

[TABLE]

where the $\ast$ represents the convolution product. In this context we call $J(t)$ the creep compliance and $G(t)$ the relaxation modulus. These function are also often referred to as material functions of the system. It is also worth remarking that both these functions are non-negative on $t\geq 0$, $J(t)$ is non-decreasing and $G(t)$ is non-increasing.

Let us now turn our attention to the differential form of the constitutive equation for a given linear system. Indeed, these kind of stress-strain relation are usually encountered while attempting a bottom-up description of a physical system. In fact, it is way easier to write down a differential equation describing the stress-strain relation for a material starting from the fundamental principles of classical mechanics. Moreover, this formalism admits the familiar description of linear viscoelastic systems in terms of networks of springs and dash-pots.

In the following we will use $D$ to designate the time derivative operator and we will use $P(D)$ and $Q(D)$ to denote two linear differential operators defined as

[TABLE]

where $p_{k},q_{k}\in\mathbb{R}$, $\forall k\in\mathbb{N}\cup\{0\}$. Clearly, the order of these operators does not exceed $N$.

Then, the stress-strain relation for a given linear system takes the form

[TABLE]

As shown in [11], if the initial conditions are such that

[TABLE]

together with some regularity conditions on the relaxation function, then the two forms of the constitutive equation are equivalent.

Moreover, if $\sigma(t),\varepsilon(t)=O\left(\exp\left(\alpha\,t\right)\right)$ as $t\to\infty$, for some $\alpha\in\mathbb{R}$, and if these functions meet the conditions (1.4), then, in the Laplace domain, we have that

[TABLE]

where we have used the notation $\mathcal{L}\{f(t)\,;\,s\}:=\widetilde{f}(s)$ with $\mathcal{L}$ denoting the Laplace transform operator.

It is worth remarking that $P(s)$ and $Q(s)$ are just polynomials in $s$ whose order does not exceed $N$. Moreover, Eq. (1.5) states that, if both the stress and the strain meet the condition (1.4), then the contribution of the initial conditions, while passing from (1.3) to (1.5) can be completely neglected.

Now, everything works quite well as long as we work within the framework of a finite linear system, i.e. if $N<\infty$. If $N=\infty$ we have to deal with constitutive equation involving infinite order differential operators. Therefore, it is worth recalling some general information about this kind of operators.

Let $\Phi(z)=\sum_{n=0}^{\infty}a_{n}\,z^{n}$ and let $f(t)$ be an entire function of $t$, then

[TABLE]

In the following $\Phi(D)$ will be referred to as differential operator and the series $\sum_{n}a_{n}\,z^{n}$ will be referred to as the generating power series of $\Phi(D)$. If the set $\{a_{n}\neq 0\,|\,n\in\mathbb{N}_{0}\}$ is infinite, then $\Phi(D)$ is said to be a infinite order differential operator. In particular, if the series in (1.6) is convergent $\forall t$, then the differential operator $\Phi(D)$ is said to be applicable to the (entire) function $f(t)$. For further details on these operators, we invite the interested reader to refer to the monograph by Sikkema [22].

The aim of this paper is to understand what happens if we consider a linear viscoelastic model whose constitutive equation is given by an equation involving infinite order differential operators. Without loss of generality, we will focus our attention on the specific case of a model of the Bessel class [5]. Nonetheless, the arguments presented in the following sections can be readily extended to any constitutive equation of order $N+\alpha$, with $\alpha\in\mathbb{R}$ and $N$ either finite (integer) or infinite.

2. The stress-strain equation for a Bessel body

Let $\sigma,\varepsilon\in\mathcal{H}^{\infty}$, and let us consider a stress-strain relation given by

[TABLE]

where $D_{\ast}$ denotes the Caputo’s fractional derivative and $P_{\nu}$ and $Q_{\nu}$ are some infinite order fractional differential operators defined in terms of the functions

[TABLE]

where $I_{\alpha}(z)$ stands for the series representation of the modified Bessel function of the first kind of order $\alpha\in\mathbb{R}$.

Now, if we assume that the initial data for Eq. (2.1) meet the condition

[TABLE]

for all $k\in\mathbb{N}$, where

[TABLE]

$\Theta$ is the Heaviside generalized function and $\lceil x\rceil$ is the celling function of $x$, then the corresponding constitutive equation, in the Laplace domain, reduces to

[TABLE]

Basically, the condition (LABEL:initial-condition) guarantees that the contribution these initial data must not appear in the Laplace domain, namely they have to be vanishing or cancel in pair-balance. Moreover, this guarantees the mutual consistency of both the integral representation of the stress-strain relation and the differential form.

Now, from Eq. (2.5) one can easily deduce the creep compliance for the model (in the Laplace domain), that turns out to be exactly the one that defines the Bessel models [5]. Indeed, we have that

[TABLE]

The Bessel creep compliance, in the time domain, is a Bernstein function given by (see [5; 10])

[TABLE]

where $j_{\nu+2,\,n}$ represents the $n$th positive real root of the Bessel functions of the corresponding order. Moreover, it is worth remarking (see [10]) this function can be thought of as the generating function for an infinite network of ordinary springs and dash-pots.

Furthermore, this last remark is also supported by the following argument. If the stress and the strain histories meet the condition (LABEL:initial-condition), then we have that (2.1) reduces to an ordinary infinite order differential equation. Indeed, we can always recast (2.5) as

[TABLE]

where

[TABLE]

[TABLE]

Then, inverting back to the time domain, we get the new constitutive equation

[TABLE]

where $D$ is, again, the ordinary derivative with respect to time.

To prove all these results we just need to focus our attention on the behaviour of the initial conditions. In the following discussion we will generalize the arguments presented in [11] concerning the treatment of initial conditions in linear viscoelasticity to the case of infinite order differential operators.

First, let us assume $P_{\nu}(D_{\ast})\,\sigma(t)$ and $Q_{\nu}(D_{\ast})\,\varepsilon(t)$ to be convergent $\forall t\geq 0$. If this is the case, then we can introduce some cut-off operators $P^{(N)}_{\nu}(D_{\ast})$, $Q^{(N)}_{\nu}(D_{\ast})$ defined as

[TABLE]

for which we have that

[TABLE]

We can now perform our analysis on the cut-off operators, for which the background formalism is well established, and then extend the results taking the limit for $N\to\infty$.

Given the previous assumptions, we have that $P^{(N)}_{\nu}(D_{\ast})\,\sigma(t)$ and $Q^{(N)}_{\nu}(D_{\ast})\,\varepsilon(t)$ are two sequences of functions that converge almost everywhere on $t\geq 0$. If we further assume that $\exists\alpha,\beta\in L_{1}([0,\infty[)$ such that

[TABLE]

for all $N\geq 1$ and almost everywhere on $t\geq 0$, then, taking profit of the dominated convergence theorem, we have that

[TABLE]

This remark then allows us to deal with finite order operators without losing the connection with the infinite order one.

For sake of brevity, let us focus only on the strain part. The corresponding analysis for the stress part of the equation can be carried out in the same way.

First, let us compute the Laplace transform of $Q^{(N)}_{\nu}(D_{\ast})\,\varepsilon(t)$. Recalling the Laplace transform of the Caputo’s fractional derivative

[TABLE]

then one can easily infer that

[TABLE]

Now, one can easily show that

Lemma 1**.**

Let $\alpha\in\mathbb{R}^{+}$ then

[TABLE]

where $\{\alpha\}$ is the fractional part of $\alpha$ and $\delta_{\alpha}$ is such that

[TABLE]

Remark*.*

Notice that for $\alpha\in\mathbb{R}^{+}$ we have that $\alpha=\{\alpha\}+\lceil\alpha\rceil-\delta_{\alpha}$.

Proof.

Starting from the RHS of (2.17), taking profit of the previous remark, one just needs to extract the integer part of $\alpha$ form $s^{\alpha}$. Then, the LHS of (2.17) is obtained by changing the summation index from $k$ to $r=\lceil\alpha\rceil-k$. Indeed, taking profit of the last remark,

[TABLE]

then, defining the index $r=\lceil\alpha\rceil-k$, we get that

[TABLE]

that concludes our proof. ∎

Then, by means of the result in Lemma 1 we can rewrite Eq. (2.16) as

[TABLE]

Following a similar argument as the one by Gurtin and Sternberg in [11], one can infer that

Lemma 2**.**

Let $a,M\in\mathbb{N}$ and let $q_{n},\,b_{n},\,s\in\mathbb{R},\,\forall n\in\mathbb{N}_{0}$, then

[TABLE]

for all $N\in\mathbb{N}$.

Proof.

This can be proved by induction.

First, for $N=1$ the LHS of Eq. (2.22) is given by

[TABLE]

Beside, the RHS of Eq. (2.22) is given by,

[TABLE]

Thus, it is clear that for $N=1$ both sides of are equal Eq. (2.22) and, therefore, Eq. (2.22) is verified for $N=1$.

Now, let us suppose that Eq. (2.22) is true for $N=\mathcal{N}$, for $\mathcal{N}>1$, then for $N=\mathcal{N}+1$ we have that the LHS is given by

[TABLE]

Beside, the RHS of Eq. (2.22) is given by,

[TABLE]

Thus, Eq. (2.22) holds for $N=\mathcal{N}+1$ and the proof of the induction step is complete. In conclusion, we have that Eq. (2.22) holds $\forall N\in\mathbb{N}$. ∎

We can now use the result in Lemma 2 to show that Eq. (2.21) can be rewritten as

[TABLE]

The proof is quite straightforward, indeed we just need to apply the result in Lemma 2 to the second term in Eq. (2.21), noticing that $M=\lceil\nu/2\rceil+1$ and $b_{n}=\varepsilon^{(n)}(0+)$.

If we perform the same analysis for $P^{(N)}_{\nu}(D_{\ast})\,\sigma(t)$ we get

[TABLE]

Then, by taking the limit for $N\to\infty$ we get $\mathcal{L}\left\{Q_{\nu}(D_{\ast})\,\varepsilon(t)\,;\,s\right\}$ and $\mathcal{L}\left\{P_{\nu}(D_{\ast})\,\sigma(t)\,;\,s\right\}$, respectively.

Hence, if we now assume that the initial data meet the condition in Eq. (LABEL:initial-condition), then

[TABLE]

reduces to (2.5). Then, all the subsequent results discussed in the first part of this section follow immediately from (LABEL:initial-condition) and (2.5).

2.1. Bessel $\&$ Fractional Maxwell of order $1/2$

From the previous discussion we have that the constitutive equation for a Bessel model of order $\nu>-1$ can be expressed as in Eq. (2.10), provided that we assume the validity of the condition in (LABEL:initial-condition).

Now, following the discussion in [6], it is interesting to see that, even if we are dealing with ordinary infinite order differential equation, as $t\to 0^{+}$ we have that the stress-strain relation reduces to

[TABLE]

which is actually Maxwell-like fractional model of order $1/2$ (see e.g. [15; 16]).

2.2. Entire functions and applicability

Now, let us focus our attention on the regularized operators $\overline{P}_{\nu}$ and $\overline{Q}_{\nu}$ defined in terms of the power series (2.8) and (2.9), respectively. In particular, it is easy to prove the following statement

Theorem 1**.**

Let $\overline{P}_{\nu}(z)$ and $\overline{Q}_{\nu}(z)$ be the functions defined as in (2.8) and (2.9). Then we can note that $\overline{P}_{\nu}(z)$ and $\overline{Q}_{\nu}(z)$ are entire functions of order $\rho=1/2$ and type $\sigma=1$.

Proof.

Let us focus our attention on the function $\overline{P}_{\nu}(z)$, then the proof for $\overline{Q}_{\nu}(z)$ follows by the very same procedure.

First, given the definition of $\overline{P}_{\nu}(z)$, we can immediately infer that it is an entire function because it is defined in terms of a modified Bessel function.

Now, given an entire function represented by the power series $\sum_{n=0}^{\infty}a_{n}\,z^{n}$, then the order $\rho$ and the type $\sigma$ of such a function are given by (see [12])

[TABLE]

Thus, for $\overline{P}_{\nu}(z)$ we have that

[TABLE]

from which it follows that

[TABLE]

Recalling Stirling’s approximation, i.e.

[TABLE]

then,

[TABLE]

Hence, for $\overline{P}_{\nu}(z)$ we can conclude that

[TABLE]

Analogously, it is easy to see that

[TABLE]

from which we can infer that

[TABLE]

result that concludes our proof. ∎

Remark*.*

We could have come to the same results in Theorem 1, alternatively, by recalling that the functions (2.8) and (2.9) are basically a rescaled version of the modified Bessel functions. Now, the latter can be seen as particular realizations of the Mittag-Leffler functions, for which a complete discussion of the order and type has been carried out by Kiryakova in [13].

This result is particularly interesting for two distinct reasons.

First, the case in which an infinite order differential operator is defined in terms of a generating power series corresponding to an entire function has been studied by several authors, see e.g. [2; 22]. In particular, according to these papers, in this case it seems to be easier to prove the applicability of such operators to certain subsets of the space of entire functions. For example, in [2] Cha et al have been able to provide a sufficient condition for two entire functions, $f(z)$ and $g(z)$, to be such that the combination $\sum_{n}f^{(n)}(0)\,D^{n}g(z)/n!$ represents an entire function. Moreover, they have also deduced the growth of the resulting function, induced by the types of both $f(z)$ and $g(z)$.

Secondly, it is interesting to note that the (fractional) order of asymptotic constitutive equation (2.10), as $t\to 0^{+}$, matches exactly with the order of the entire functions defining the differential operators. Whether there is a formal connection between these two results is still unknown, but it is surely worth of further studies.

3. Conclusions

After a brief review of the basic formalism of linear viscoelasticity, we have discussed some properties of viscoelastic models defined in terms of constitutive equations involving infinitely many derivatives. In particular, we have considered, as a working example, the recently developed Bessel models of linear viscoelasticiy.

In Section 2, we have been able to connect the constitutive equation for a Bessel model with the theory of ordinary infinite order differential equations. We have also determined the form of the initial data that guarantees the mutual consistency of both the integral and the differential form of the stress-strain relation. This condition also allows for a constitutive equation, in the Laplace domain, that does not involve the explicit contribution of initial data. Finally, in Section 2.1 and 2.2 we have observed that, even thought the stress-strain relation for a Bessel body involves infinitely many ordinary derivatives, the asymptotic behaviour as $t\to 0^{+}$ is effectively described by a fractional Maxwell model of order $1/2$. We have also pointed out that the order of the asymptotic (fractional) model is the same as the order of the entire functions that define the regularized operators (2.8) and (2.9).

Acknowledgements

The work of the author has been carried out in the framework of the activities of the National Group of Mathematical Physics (GNFM, INdAM).

I am particularly grateful to Prof. Francesco Mainardi for all of the support and mentorship he has shown me over the past eight years and to Prof. Roberto Casadio for helpful comments and discussions, and also for his patience and constant support. I would also like to extend my thanks to I. Colombaro, V. A. Diaz, Prof. M. Fabrizio, R. Garra, Prof. A. Kamenchtchik, Prof. V. Kiryakova, Prof. S. Rogosin and Prof. T. Ruggeri for many helpful discussions.