A Note on the Lambert W Function: Bernstein and Stieltjes Properties for a Creep Model in Linear Viscoelasticity

Abstract

The purpose of this note is to propose an application of the Lambert W function in linear viscoelasticity based on the Bernstein and Stieltjes properties of this function. In particular, we recognize the role of its main branch, $W_0\left(t\right)$, in a peculiar model of creep with two spectral functions in frequency that completely characterize the creep model. In order to calculate these spectral functions, it turns out that the conjugate symmetry property of the Lambert W function along its branch cut on the negative real axis is essential. We supplement our analysis by computing the corresponding relaxation function and providing the plots of all computed functions.

1. Introduction

The Lambert function $W\left(z\right)$ is defined as the root of the transcendental equation

$$W\left(z\right)exp\left(W\left(z\right)\right)=z.$$

(1)

The mathematical history of the W function goes back to the 18th century as outlined in the seminal paper [1] where the interested reader can be informed about the analytical and numerical properties of this function. Further details can be found in several papers, see, e.g., [2,3,4,5,6], Section 4.13 of the handbook [7], and in the recent book [8] with references therein. For our purposes, we outline the paper [9]. The applications are found in many areas of applied science as outlined, e.g., in [8,10,11] and references therein. Let us outline also the applications in probability, see, e.g., [12,13,14], which could be considered and possibly revised in view of the present results in linear viscoelasticity.

In our analysis, we restrict our attention to the main branch $W_0\left(z\right)$ of the W function on the real semiaxis $z=t\ge 0$. We claim that the Lambert function $W_0\left(t\right)$ can be a candidate involved as a possible model for dimensionless creep compliance (i.e., a material function) in linear viscoelasticity. Indeed, according to a previous idea of [15], the property of being a Bernstein function (that is, with a completely monotonic derivative) for such material function is proved to be fundamental as stated in Chap. 2 of the book [16] and references therein. In addition, the property of also being a Stieltjes function provides further results for this viscoelastic model that, as far as we know, are new in the framework of the literature of linear viscoelasticity.

The plan of this paper is as follows. In Section 2, we summarize the essentials of linear viscoelasticity in order to point out the properties of the rate of creep of our model based on the Lambert function. As a consequence, the rate of creep, being a completely monotonic (CM) function with the additional requirement to be a Stieltjes function, turns out to be expressed in terms of two different spectral functions. In Section 3, we carry out a numerical and analytical analysis of our new model based on the properties of the Lambert function. We also illustrate with plots the related quantities which can better characterize the model itself. The consistence of our results are validated with MATHEMATICA. In the Conclusions section, we summarize our final remarks.

2. Essentials of Linear Viscoelasticity

According to the theory of linear viscoelasticity, viscoelastic bodies are characterized by two different interrelated functions, called material functions, which are causal in time (i.e., null for $t<0$). The first material function is the creep compliance $J\left(t\right)$, defined as the strain response to a unit step function of stress. The second material function is the relaxation modulus $G\left(t\right)$, which is defined as the stress response to a unit step function of strain. For more details, see, e.g., [17,18,19] and the recent book [16].

Taking $J\left(0^{+}\right)=J_0>0$ (thus $G\left(0^{+}\right)=G_0=1/J_0$), a viscoelastic body is assumed to show a nonvanishing instantaneous response both in creep and relaxation tests. Therefore, it is convenient to introduce the following dimensionless functions $\psi \left(t\right)$ and $\varphi \left(t\right)$:

$$J\left(t\right)=J_0[1+q\psi \left(t\right)],G\left(t\right)=G_0\varphi \left(t\right),$$

(2)

where $\psi \left(t\right)$ is a non-negative increasing function with

$$\psi \left(0\right)=0,$$

(3)

and $\varphi \left(t\right)$ is a non-negative decreasing function with

$$\varphi \left(0\right)=1.$$

(4)

The nondimensional quantity q takes into account a suitable scaling of the strain, according to convenience, in experimental rheology. However, in the following we assume for simplicity $q=1$ without loss of generality.

Viscoelastic bodies may be distinguished in solid like and fluid like, whether $J(+\infty )$ is finite or infinite, so that $G(+\infty )=1/J(+\infty )$ is nonzero or zero, respectively.

As pointed out in most treatises on linear viscoelasticity, (e.g., in [16,18,19]), the relaxation modulus $G\left(t\right)$ can be derived from the corresponding creep compliance $J\left(t\right)$ through the Volterra integral equation of the second kind

$$G\left(t\right)=\frac{1}{J_0}-\frac{1}{J_0}{\int }_0^t\frac{dJ}{d{t}^{\prime }}G(t-t^{\prime })d{t}^{\prime }.$$

(5)

Therefore, according to (2), the nondimensional relaxation function $\varphi \left(t\right)$ obeys the Volterra integral equation

$$\varphi \left(t\right)=1-q{\int }_0^t\frac{d\psi }{d{t}^{\prime }}\varphi (t-t^{\prime })d{t}^{\prime },q=1.$$

(6)

In the Appendix B, we calculate the solution of (6) by using the Laplace transform. Hereafter, we use the following notation for the Laplace transform pair:

$$f\left(t\right)÷\tilde{f}\left(s\right)=\mathcal{L}[f\left(t\right);s]:={\int }_0^{\infty }{e}^{-st}f\left(t\right)dt,$$

(7)

where the sign ÷ denotes the juxtaposition of the function $f\left(t\right)$ with its Laplace transform $\tilde{f}\left(s\right)$. In order to distinguish the Laplace transform, we use the same notation f for the original function, but overlined with a tilde and with the proper argument s.

According to this notation, the solution of (6) derived in Appendix A is written as

$$\varphi \left(t\right)={\mathcal{L}}^{-1}\left[\frac{1}{s\left(1+s\tilde{\psi }\left(s\right)\right)};t\right].$$

(8)

It is quite usual to require the existence of positive retardation and relaxation spectra for the material functions $J\left(t\right)$ and $G\left(t\right)$, as pointed out in the monograph [15]. This implies (as formerly proved in [20] and then revisited by [21], as well as in [16]) that $J\left(t\right)$ and $G\left(t\right)$ and consequently, the dimensionless functions $\psi \left(t\right)$ and $\varphi \left(t\right)$, turn out to be Bernstein and completely monotonic functions, respectively.

We recall that a completely monotonic (CM) function is a non-negative, infinitely derivable function whose derivatives alternate in sign for $t>0$. Moreover, a Bernstein function is a non-negative function whose derivative is CM. According to Bernstein’s theorem, $f\left(t\right)$ is the Laplace transform of a non-negative function, if and only if $f\left(t\right)$ is a CM function (the interested reader is referred to the excellent monograph [22] for more details on the mathematical properties of CM functions). Apply Bernstein’s theorem to write the rate of creep as follows:

$${\psi }^{\prime }\left(t\right)={\int }_0^{\infty }{e}^{-rt}K\left(r\right)dr={\int }_0^{\infty }{e}^{-t/\tau }H\left(\tau \right)d\tau ,$$

(9)

where $K\left(r\right)$ and $H\left(\tau \right)$ are both non-negative functions and denote the required spectra in frequency (r) and in time ($\tau =1/r$), respectively.

Then, from (9), the time spectrum can be determined using the transformation $\tau =1/r$, so that

$$H\left(\tau \right)=\frac{K(1/\tau )}{{\tau }^2}.$$

(10)

We recognize that the Laplace transform of the rate of creep is the iterated Laplace transform of the frequency spectrum, that is, the Stieltjes transform of $K\left(r\right)$. Therefore, the Titchmarsh formula provides the inversion of the Stieltjes transform, see, e.g., [23]. Indeed, since ${\psi }^{\prime }\left(t\right)$ is CM, taking into account (9) and (A10) with (3), we have for its Laplace transform:

$$\begin{array}{ccc}\hfill \mathcal{L}[{\psi }^{\prime }\left(t\right);s]& =& s\tilde{\psi }\left(s\right)={\displaystyle {\int }_0^{\infty }{e}^{-st}{\psi }^{\prime }\left(t\right)dt}\hfill \\ & =& {\displaystyle {\int }_0^{\infty }{e}^{-st}\left({\int }_0^{\infty }{e}^{-rt}K\left(r\right)dr\right)dt},\hfill \end{array}$$

and exchanging the order of integration, we finally obtain:

$$s\tilde{\psi }\left(s\right)={\displaystyle {\int }_0^{\infty }\left({\int }_0^{\infty }{e}^{-t(r+s)}dt\right)K\left(r\right)dr}={\displaystyle {\int }_0^{\infty }\frac{K\left(r\right)}{r+s}dr}.$$

(11)

We thus recognize that $K\left(r\right)$ is the inverse of the Stieltjes transform of $s\tilde{\psi }\left(s\right)$. Under suitable conditions, the inversion can be obtained by the Titchmarsh formula. (This formula is found without names in [23]. However, in some papers and books, it is cited with several names, just Titchmarsh in [16], Stieltjes–Perron in [24], Gross–Levi in [25], Bobylev–Cercignani in [26]. In any case, it results as a simple exercise in complex analysis as stated in [27] for the Mittag-Leffler function.) Thus,

$$K\left(r\right)=\pm \frac{1}{\pi }{\left.\Im \left[s\tilde{\psi }\left(s\right)\right]\right|}_{{\displaystyle s=r}{e}^{\mp i\pi }}.$$

(12)

In addition, if ${\psi }^{\prime }\left(t\right)$ is a Stieltjes function (i.e., the Laplace transform of a CM function), then $K\left(r\right)$ turns out to be the Laplace transform of a non-negative function that we denote by $\rho \left(u\right)$. Indeed,

$$\begin{array}{cc}{\psi }^{\prime }\left(t\right)=\hfill & {\displaystyle {\int }_0^{\infty }{e}^{-rt}K\left(r\right)dr}={\displaystyle {\int }_0^{\infty }{e}^{-rt}\left({\int }_0^{\infty }{e}^{-ru}\rho \left(u\right)du\right)dr}\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }{e}^{-ru}\left({\int }_0^{\infty }{e}^{-rt}dr\right)\rho \left(u\right)du}\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }\left({\int }_0^{\infty }{e}^{-r(t+u)}d\tau \right)\rho \left(u\right)du}={\displaystyle {\int }_0^{\infty }{\displaystyle \frac{\rho \left(u\right)}{t+u}}du},\hfill \end{array}$$

(13)

where we have applied

$$K\left(r\right)={\int }_0^{\infty }{e}^{-ru}\rho \left(u\right)du.$$

(14)

Therefore, from (13), and applying the Titchmarsh formula, we obtain

$$\rho \left(u\right)=\pm \frac{1}{\pi }{\left.\Im \left[{\psi }^{\prime }\left(t\right)\right]\right|}_{{\displaystyle t=u}{e}^{\mp i\pi }}.$$

(15)

Hence, in the case that the rate of creep ${\psi }^{\prime }\left(t\right)$ is a Stieltjes function, we have another spectral function $\rho \left(u\right)$ related to the previous spectral function $K\left(r\right)$ by a Laplace transformation.

3. Application of the Lambert Function in Linear Viscoelasticity

After recalling the properties of the Lambert W function, it is worth to see how its plot looks like compared with its asymptotic representation for large arguments. According to Figure 1, the function $W_0\left(t\right)$ appears as a positive one, increasing from 0 at $t=0$ up to ∞ as $t\to \infty$. The two-term asymptotic representation,

$$W_0\left(t\right)\sim log\left(t\right)-loglog\left(t\right),t\to \infty ,$$

(16)

fits the function slowly as $t\to \infty$.

We can see that the $W_0\left(t\right)$ function is positive increasing up to ∞ according to its property of being a Bernstein function, that is, a non-negative function with a CM first derivative in $t\ge 0$, see, e.g., [4], p. 116. Furthermore, its derivative is known to be a Stieltjes function, that is, according to the definition stated in the previous section, a CM function represented by a real Laplace transform of a CM function, where both functions exhibit non-negative spectral functions, see [9].

Differentiating in (1) and solving for $W^{\prime }$, we obtain the following expressions for the derivative of $W_0$:

$$W_0^{\prime }\left(t\right)={\displaystyle \frac{1}{e^{W_0\left(t\right)}+t}}={\displaystyle \frac{W_0\left(t\right)}{t\left[1+W_0\left(t\right)\right]}},t>0,$$

(17)

and its two-term asymptotic representation is

$$W_0^{\prime }\left(t\right)\sim {\displaystyle \frac{1}{t}}-{\displaystyle \frac{1}{t+log\left(t\right)}},t\to \infty .$$

(18)

The plots of the derivative of $W_0\left(t\right)$ and its asymptotic representation are depicted in Figure 2.

In the proposed viscoelastic model, we assume for the creep that

$$\psi \left(t\right)=W_0\left(t\right).$$

(19)

From this assumption, and also taking into account that $W_0^{\prime }\left(t\right)$ is a Stieltjes function, we derive the corresponding spectral functions $K\left(r\right)$ and $\rho \left(u\right)$ presented in the previous section. Indeed, the spectral function $\rho \left(u\right)$ can be determined from (15) and (19), i.e.,

$$\rho \left(u\right)=-{\displaystyle \frac{1}{\pi }}\Im \left[W_0^{\prime }(-u)\right],u>0,$$

(20)

where $W_0^{\prime }\left(z\right)$ is computed over or down the negative semiaxis that would be its branch cut in the complex plane.

According to the theory of the Titchmarsh formula, the function $\rho \left(u\right)$ would be non-negative and continuous. However, our application of the Titchmarsh formula provides a non-negative function, but discontinuous, see Figure 3, where we note that $\rho \left(u\right)=0$ for $0<u<1/e$ and non-negative for $u>1/e$. In particular, $\rho \left(u\right)$ is decreasing from ∞ in the limit $u\to {(1/e)}^{+}$ to zero in the limit $u\to +\infty$.

In order to validate the Stieltjes properties with the related Titchmarsh formula, we need to compute the iterated Laplace transform of the (discontinuous) spectral function $\rho \left(u\right)$ and verify that it is the original derivative $W_0^{\prime }\left(t\right)$ (see Equation (A15)) except for suitable small numerical errors. Indeed, by using MATHEMATICA, this result was verified because the relative error was small enough, as shown in Figure 4. Nonetheless, we obtained an alternative derivation of the $\rho \left(u\right)$ function in Appendix A. It is worth noting that the conjugate symmetry property of the Lambert W function along its branch cut on the negative real axis, i.e., (A7), turns out to be essential for this derivation.

We note that a simple reasoning (validated by MATHEMATICA) provides the value of one for the Riemann generalized integral in the positive real semiaxis, that is

$${\int }_0^{\infty }\rho \left(u\right)du=1.$$

(21)

Indeed, taking into account (20), performing the change of variables $v=-u$, and knowing that $W_0\left(0\right)=0$, we have

$$\begin{array}{ccc}\hfill {\int }_0^{\infty }\rho \left(u\right)du& =& -\frac{1}{\pi }\Im {\int }_0^{\infty }{W}_0^{\prime }\left(-u\right)du\hfill \\ & =& -\frac{1}{\pi }\Im {\int }_{-\infty }^0{W}_0^{\prime }\left(v\right)dv\hfill \\ & =& -\frac{1}{\pi }\underset{t\to -\infty }{lim}\Im \left[W_0\left(0\right)-W_0\left(t\right)\right]\hfill \\ & =& \frac{1}{\pi }\underset{t\to -\infty }{lim}\Im \left[W_0\left(t\right)\right].\hfill \end{array}$$

According to the property, see [9],

$$\underset{t\to -\infty }{lim}\Im \left[W_0\left(t\right)\right]=\pi ,$$

we conclude that Equation (21) is true.

Similar properties are shared by $\rho \left(u\right)/u$ for $u\in (0,\infty )$, as validate by MATHEMATICA, but of course with different plots as it is shown in Figure 5.

Indeed, the value of the infinite integral of $\rho \left(u\right)/u$ in the positive semiaxis is also one,

$${\int }_0^{\infty }\frac{\rho \left(u\right)}{u}du=1=W_0^{\prime }\left(0\right),$$

(22)

where this result is based on (13) for $t=0$, (19), as well as the fact that $W_0^{\prime }\left(0\right)=1$, according to Figure 2.

According to (12) and (19), we need the expression of the Laplace transform of $W_0\left(t\right)$ in order to calculate the spectral function $K\left(r\right)$. However, the latter is not available in the literature. Nevertheless, we can numerically evaluate the spectral function $K\left(r\right)$ from (14) with the expression given in (20) for the $\rho \left(u\right)$ function, i.e.,

$$K\left(r\right)=-\frac{1}{\pi }{\int }_0^{\infty }\Im \left[W_0^{\prime }\left(-u\right)\right]du,$$

(23)

as it is shown in Figure 6. Nonetheless, since the $\rho \left(u\right)$ function given in (20) is subjected to a numerical validation as discussed above, we provide an analytical derivation of (23) in Appendix A.

Figure 7 shows the plot of the spectral function $H\left(\tau \right)$ in time.

We recognize that $H\left(\tau \right)$ is not a CM function because its relation with $K\left(r\right)$ in the transformation $r\to \tau =1/r$ cannot preserve this relevant property.

Now, let us calculate the dimensionless relaxation function $\varphi \left(t\right)$ solving the Volterra integral Equation (6) for our model, i.e., (19). According to (8), we have to compute

$$\varphi \left(t\right)={\mathcal{L}}^{-1}\left[\frac{1}{s\left(1+s\widetilde{W_0}\left(s\right)\right)};t\right].$$

(24)

Figure 8 shows the numerical evaluation of (24). Note that the graph of $\varphi \left(t\right)$ confirms that it is a non-negative decreasing function with $\varphi \left(0\right)=1$, in agreement with (4).

4. Conclusions

In this paper, we pointed out some properties of the Lambert function taking advantage of the large existing literature on that function. In particular, we used these properties to propose a model in linear viscoelasticity that appears to be novel, as far as we know. This model was well characterized, as far as the creep was concerned, with two spectral functions. These spectral functions were numerically validated with MATHEMATICA, as well as analytically derived in Appendix A.

It is worth noting that this model provides an instructive slow-varying creep function, slower than a logarithmic law exhibited by the Lomintz and Becker models known in the geophysical literature, see [16], Section 2.9.1.

To complete the analysis of our model we provided the calculation of the corresponding relaxation function.