%~Mouliné par MaN_auto v.0.40.4 (2e4f8ce3) 2026-04-13 09:38:05
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_de-Sousa_20250850.bib}

\usepackage{siunitx}
\usepackage{subcaption}
\usepackage{algorithm}
\usepackage{listings}
\usepackage{booktabs}

\newcommand{\GA}{\mathrm{GA}}
\newcommand{\NLP}{\mathrm{NLP}}
\newcommand{\NOP}{\mathrm{NOP}}

\newcommand{\EG}{\mathit{EG}}

\newcommand{\void}{\,\cdot\,}

\newcommand{\tot}{\mathrm{tot}}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}
\newcommand*\Dif[1]{\mathop{}\!\mathrm{d}^{#1}\!}

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

%%%---------------------------------------------------------------------------

%%% BAR & OVERLINE

%%% Better widebar macro
\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
 }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
 }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

%%% ENVIRONEMENTS

%%% Env. {system} (mimicks {cases} but only supports a single column)
\makeatletter
\newenvironment{system}{
  \let\@ifnextchar\new@ifnextchar
  \left\lbrace
  \def\arraystretch{1.2}
  \array {@{}l@{}}
}{\endarray\right.}
\makeatother

%%%---------------------------------------------------------------------------

%%% VERTICAL BAR

%%% ''Evaluated at'' macro for CRAS
\newcommand\disprestr[3][]{{				
 \left.\kern-\nulldelimiterspace 
 #2 								
 \vphantom{\big\vert} 				
 \right\vert_{#3}^{#1}
}}
\newcommand{\restr}[3][]{
	\mathchoice{\disprestr[#1]{#2}{#3}}
			   {\disprestr[#1]{#2}{#3}}
			   {#2\rvert_{#3}^{#1}}
			   {#2\rvert_{#3}^{#1}}
	}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}

\title{Mathematical modeling of Poisson’s ratio and time-frequency interconversion of viscoelastic models based on the fractional Zener formulation}
\alttitle{Modélisation mathématique du coefficient de Poisson et de l'interconversion temps-fréquence des modèles viscoélastiques basée sur la formulation fractionnaire de Zener}

\author{\firstname{Tiago} \middlename{Lima} \lastname{de Sousa}\CDRorcid{0000-0002-8992-6293}\IsCorresp}
\address{Federal University of Technology -- Paraná (UTFPR), 330 Dr. Washington Subtil Chueire St., Ponta Grossa, Paraná, Brazil}
\email{sousa@professores.utfpr.edu.br}

\author{\firstname{Jucélio} \middlename{Tomás} \lastname{Pereira}\CDRorcid{0000-0002-2483-4339}}
\address{Federal University of Paraná (UFPR), 100 Cel. Francisco H. dos Santos Ave., Curitiba, Paraná, Brazil}
\email{jucelio.tomas@ufpr.br}

\author{\firstname{Jéderson} \lastname{da Silva}\CDRorcid{0000-0002-2714-3330}}
\address{Federal University of Technology -- Paraná (UTFPR), 731 João Miguel Caram Ave., Londrina, Paraná, Brazil}
\email{jedersonsilva@utfpr.edu.br}

\thanks{T.~L.~de~Sousa gratefully acknowledges the Foundation for Research Support of the State of Amazonas (FAPEAM, Brazil) for the financial support provided through a doctoral scholarship, without which the present work would not have been possible. The author also thanks the Postgraduate Program in Mechanical Engineering (PG-Mec) of the Federal University of Paraná (UFPR) and the Federal University of Technology -- Paraná (UTFPR) for their institutional support. Prof.\ J.~T.~Pereira acknowledges the financial support received from the National Council for Scientific and Technological Development (CNPq, Brazil).}
\CDRGrant[UKRC]{2019-155900}

\keywords{\kwd{Viscoelasticity} \kwd{Poisson's Ratio} \kwd{time domain modeling} \kwd{correspondence principle}}
\altkeywords{\kwd{Viscoélasticité} \kwd{coefficient de Poisson} \kwd{modélisation dans le domaine temporel} \kwd{principe de correspondance}}

\begin{abstract}
The proposed formulation, based on the four-parameter fractional Zener model, provides a versatile constitutive framework for describing the mechanical behavior of viscoelastic materials (VEMs). The central objective of this work is to theoretically and practically demonstrate that viscoelastic models identified in the time domain can be consistently transformed into their frequency-domain counterparts, and vice versa. The principal contribution lies in the development of a new mathematical model for the time-dependent Poisson’s ratio, formulated in the time domain and derived directly from the constitutive relations of the fractional Zener model. Artificial experimental datasets are employed to validate the effectiveness and internal consistency of the proposed interconversion methodology. Once the fractional parameters are identified through an optimization approach, the corresponding complex viscoelastic functions ---~namely, the complex Young’s modulus, complex shear modulus, and complex Poisson’s ratio~--- are obtained through analytical interconversion into the frequency domain. Overall, the proposed framework reinforces the theoretical foundation connecting time- and frequency-domain representations of viscoelastic behavior and advances the modeling and characterization of viscoelastic materials.
\end{abstract}

\begin{altabstract}
La formulation proposée, fondée sur le modèle de Zener fractionnaire à quatre paramètres, offre un cadre constitutif polyvalent pour décrire le comportement mécanique des matériaux viscoélastiques (VEM). L'objectif principal de ce travail est de démontrer, tant sur le plan théorique que pratique, que les modèles viscoélastiques identifiés dans le domaine temporel peuvent être transformés de manière cohérente en leurs équivalents dans le domaine fréquentiel, et inversement. La principale contribution réside dans le développement d'un nouveau modèle mathématique pour le coefficient de Poisson dépendant du temps, formulé dans le domaine temporel et dérivé directement des relations constitutives du modèle fractionnaire de Zener. Des jeux de données expérimentales artificielles sont utilisés pour valider l'efficacité et la cohérence interne de la méthodologie d'interconversion proposée. Une fois les paramètres fractionnaires identifiés par une approche d'optimisation, les fonctions viscoélastiques complexes correspondantes ---~à savoir le module de Young complexe, le module de cisaillement complexe et le coefficient de Poisson complexe~--- sont obtenues par interconversion analytique dans le domaine fréquentiel. Dans l'ensemble, le cadre proposé renforce les fondements théoriques reliant les représentations du comportement viscoélastique dans les domaines temporel et fréquentiel et fait progresser la modélisation et la caractérisation des matériaux viscoélastiques
\end{altabstract}

\COI{The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.}

\dateposted{2026-04-24}
\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\section{Introduction}

In structural engineering applications, viscoelastic materials (VEMs) are widely employed both as structural components and as integral elements of vibration and noise control systems~\cite{DIG, BFLB, GHGHH, BGM}. In such contexts, two- and three-dimensional stress-strain analyses require a thorough understanding of the characteristic viscoelastic functions that govern the mechanical behavior of these materials, which may be defined either in the time or frequency domains.

In the time domain, the mechanical behavior of isotropic, linear, and thermorheologically simple solid VEMs is typically described by the Young’s relaxation modulus, $E(t)$, the shear relaxation modulus, $G(t)$, and the time-dependent Poisson’s ratio, $\upsilon(t)$. Analogous functions are also commonly introduced based on creep tests.

Conversely, in the frequency domain, VEMs are usually characterized by their complex viscoelastic functions: the complex Young’s modulus, $E^{*}(\Omega)$, the complex shear modulus, $G^{*}(\Omega)$, and the complex Poisson’s ratio, $\upsilon^{*}(\Omega)$~\cite{Arikoglu, Benedetto, Takarli, Zhou, Pawlak, Sousa2024, Sousa2025}.

In recent years, several studies~\cite{Park, Emri, Chae, Gergesova, Barruetabeña, Lin, Pereira, Saprunov, Rouleau, Gergesova2016, Ciniello, Álvarez} have addressed the identification of relaxation and creep compliance (Young’s and shear). In contrast, only a limited number of works have focused on the indirect identification of Poisson’s ratio~\cite{Theocaris, Tsou, Ashrafi, Grassia, Chen}.

For perfectly incompressible materials, the Poisson’s ratio takes the value of $0.5$. Elastomers are generally regarded as incompressible materials, and thus their Poisson’s ratio is often assumed to be constant and equal to $0.5$~\cite{Sim, Espíndola, Hecht}. Such assumptions, however, neglect the effects of time and frequency. In practice, this parameter exhibits variability in both the frequency and time domains~\cite{Chen2017, Pritz, Pritz2000, Tschoegl, Pritz2007, Sousa2018, Cui2022, Cabrera}.

One of the pioneering studies on the indirect identification of Poisson’s ratio was presented by~\cite{Theocaris}, which introduced a methodology to determine the time-dependent Poisson’s ratio from the Young’s relaxation modulus, $E(t)$, and the shear creep compliance, $J(t)$. In this framework, the time-dependent Poisson’s ratio, $\upsilon(t)$, is obtained through the following relation:
\begin{equation}\label{v_t}
\upsilon (t)=\frac{1}{2}\bracks*{E_{0}J(t)-2+ \int_{0}^{t}\frac{\dif E(\xi)}{\dif \xi}J(t-\xi)\dif \xi }.
\end{equation}
Here, $E_{0}$ denotes the instantaneous modulus, defined as the limit $\lim_{t \to 0} E(t)$, and $\xi$ represents an auxiliary integration variable in the time domain. To evaluate the convolution integral, \cite{Theocaris} employed an approximate method. Subsequently, \cite{Tsou} determined Poisson’s ratio using an indirect procedure based on the tensile relaxation modulus and the time-dependent bending modulus. In that study, two materials were analyzed: polycarbonate and cellulose acetate. For polycarbonate, Poisson’s ratio was found to remain constant, whereas for cellulose acetate, it followed a parabolic function with a distinct maximum. The authors of~\cite{Tsou} noted that these results warrant further investigation, particularly since the experiments were conducted without accounting for temperature effects.

Further, another methodology for determining Poisson’s ratio from tensile relaxation and creep tests was presented by~\cite{Ashrafi}. In this work, a model is proposed in which Poisson’s ratio depends on a time-dependent bulk modulus and a time-dependent relaxation (or creep) modulus. To eliminate the convolution integral and reduce computational effort, the bulk modulus is assumed constant, and a viscoelastic function for Poisson’s ratio is identified using a Prony series representation. Reference~\cite{Ashrafi} concludes that there is essentially no significant difference between Poisson’s ratio obtained indirectly from relaxation tests and that obtained from creep tests.

Along the same research line, \cite{Grassia} investigated Poisson’s ratio in the time domain. Using a one-term Prony series, the authors indirectly determined a function for Poisson’s ratio based on the tensile relaxation modulus and the bulk modulus. Their results suggest that this coefficient may exhibit either a monotonically increasing trend or a non-monotonic trend with a minimum point. However, the study~\cite{Grassia} did not extend the analysis to a larger number of Prony terms, which, according to~\cite{ParkShapery}, could lead to a more accurate identification of viscoelastic functions.

One of the most relevant and recent contributions to the determination of Poisson’s ratio in the time domain was presented by~\cite{Chen}, which introduced a methodology that first identifies the Young’s and shear relaxation moduli in the time domain, expressed in terms of Prony series. Subsequently, the study establishes relations among the relaxation and creep moduli and Poisson’s ratio, thereby enabling the determination of the parameters associated with both the creep modulus and Poisson’s ratio. In this approach, the properties related to Poisson’s ratio are obtained under the assumption of a constant bulk modulus.

Another important study, \cite{Charpin}, considers a composite material consisting of elastic inclusions embedded in a Maxwell non-aging viscoelastic matrix. By incorporating the volume fraction and applying the method proposed in~\cite{Mori}, the authors derive Poisson’s ratio for the composite system. The results demonstrate that viscoelastic Poisson’s ratios may exhibit monotonic behavior (either increasing or decreasing) or nonmonotonic behavior.

 From this overview, it becomes evident that a clear gap remains in the literature regarding the analytical formulation of the viscoelastic Poisson’s ratio within a unified constitutive framework. Most available methodologies determine this parameter indirectly through numerical procedures or rely on simplifying assumptions-most notably the hypothesis of a constant bulk modulus. Such assumptions may compromise physical consistency, since experimental and theoretical studies~\cite{Pritz, Tschoegl1989, EmriProdan} have demonstrated that the bulk modulus is, in general, time-dependent. Moreover, certain reported results contradict fundamental thermodynamic and mechanical constraints. For instance, \cite{Tschoegl, Tschoegl1989} showed that, for linear isotropic viscoelastic solids, Poisson’s ratio must exhibit a monotonically increasing behavior in the time domain.

 In addition, constitutive modeling of linear isotropic viscoelasticity is commonly formulated in terms of shear and bulk moduli (or equivalently Lamé parameters), while Poisson’s ratio is treated as a secondary quantity derived a posteriori. This practice does not explicitly ensure analytical coherence between time- and frequency-domain representations, particularly when fractional constitutive models are employed. Considering that viscoelastic materials exhibit coupled volumetric and deviatoric responses that are inherently time- and frequency-dependent, a consistent formulation of Poisson’s ratio as a primary viscoelastic function becomes essential for reliable numerical modeling and structural analyses.

 The present work is developed under the assumptions of isotropy, linear viscoelasticity, and thermorheological simplicity, which define the scope of applicability of the proposed identification procedure. Within this controlled constitutive framework, the main contribution of this study consists in the analytical derivation of the time-dependent Poisson’s ratio based on the four-parameter fractional Zener model, together with its explicit interconversion into the frequency domain. Unlike previous approaches, the proposed formulation provides closed-form expressions in both domains and enables the indirect identification of Poisson’s ratio directly from experimentally accessible Young’s and shear relaxation moduli, without imposing a constant bulk modulus.

 It is important to acknowledge that, in practical experimental identification procedures, deviations from linearity, isotropy, or thermorheological simplicity may introduce structural model bias and additional metrological challenges. The present study intentionally isolates the analytical consistency and identifiability of the proposed formulation under controlled constitutive assumptions. While this approach does not eliminate the potential influence of model bias in real-material applications, it allows the intrinsic theoretical properties of the formulation to be rigorously assessed. In practical implementations, systematic residual patterns or persistent discrepancies in the identification process may indicate the need for model enrichment beyond the assumptions adopted herein.

Given the versatility of viscoelastic formulations in both time and frequency domains, this work centers on the analysis and demonstration of viscoelastic models in the time domain, followed by their interconversion into the frequency domain. This process emphasizes the theoretical and practical coherence between the two representations of viscoelastic behavior.

\section{Theoretical concepts}

\subsection{Introduction}

In general, viscoelastic materials (VEMs) are characterized by exhibiting both elastic and viscous behaviors simultaneously. Their mechanical response can be represented through rheological models composed of springs and dampers, where the springs account for the elastic component and the dampers for the viscous component. Several authors~\cite{Chae, Barruetabeña, Mead, Welch, Hu, Jrad} employ these models in various configurations, such as series, parallel, and series-parallel arrangements. According to \cite{Tschoegl1989, Mainardi}, several constitutive equations based on models with fractional-order derivatives arise from these associations and can be expressed in the following generalized form:
\begin{equation}\label{const_eq}
\bracks[\Bigg]{1+\sum_{r=1}^{n} a_{r}\frac{\Dif{\beta_r}}{\dif t^{\beta_{r}}}}\sigma (t)=\bracks[\Bigg]{k+\sum_{r=1}^{n} b_{r}\frac{\Dif{\beta_r}}{\dif t^{\beta_{r}}}}\varepsilon (t),
\end{equation}
where $n$, $k$, $a_{r}$, and $b_{r}$ are parameters associated with the adopted rheological model; $\sigma(t)$ denotes the stress history; $\varepsilon(t)$ represents the strain history; and $\frac{\Dif{\beta_r}(\void)}{\dif t^{\beta_{r}}}$ is a fractional differential operator of order $\beta_{r}$, defined as $\beta_{r}=r+\beta -1$, with $0< \beta < 1$~\cite{Mainardi}. In this work, the fractional derivative is defined according to the Riemann--Liouville formulation presented in~\cite{Mainardi}, considering only data for $t>0$.

According to~\cite{Kassem}, the presence of the integral operator in the viscoelastic constitutive equations, such as Eq.~\eqref{v_t} or Eq.~\eqref{Boltzmann}, implies that solutions to boundary-value problems derived for linear elasticity cannot be directly applied to VEMs. As a consequence of the dependence of the material response on the stress or strain history, the stress, strain, and compliance (or modulus) functions cannot simply be obtained by substituting their elastic counterparts. In this context, a method known as the Correspondence Principle, developed in~\cite{Read}, was introduced to enable the use of elastic relationships in linear viscoelastic materials.

\subsection{Correspondence principle}

The Correspondence Principle for Elastic-Viscoelastic (CPEV) materials involves the application of the Laplace transform (LT), which allows the elimination of the differential (or integral) operators, thereby yielding a constitutive equation in an algebraic form analogous to those of elastic materials. As a result, the transformed equations can be employed to model the mechanical behavior of the material and to determine the stress-strain fields in VEMs.

To illustrate the CPEV, the Laplace transform is first applied to Eq.~\eqref{const_eq}, yielding
\begin{equation}\label{constitive_equation}
a(s)\,\sigma(s) = b(s)\,\varepsilon(s)
\end{equation}
where $a(s) = 1 + \sum_{r=1}^{n}(s^{\beta_r} a_r)$ and $b(s) = k + \sum_{r=1}^{n}(s^{\beta_r} b_r)$.

From Eq.~\eqref{constitive_equation}, Hooke’s law in the Laplace domain can be expressed as
\begin{equation}\label{HookLaw_LaplaceDOmain}
\sigma(s) = \hat{E}(s)\,\varepsilon(s),
\qquad
\varepsilon(s) = \hat{D}(s)\,\sigma(s)
\end{equation}
in which $\hat{E}(s) = \frac{b(s)}{a(s)}$ and $\hat{D}(s) = \frac{a(s)}{b(s)}$ are referred to as the operational Young’s modulus and the operational creep modulus in tension, respectively~\cite{ParkShapery, EmriProdan, Rahman}.

Observe that Eq.~\eqref{HookLaw_LaplaceDOmain} is expressed for the case of uniaxial tension testing. An analogous formulation can be established for pure shear testing, where the relationship involves stress, $\tau(\void)$, and shear strain, $\gamma(\void)$. Accordingly, in the Laplace domain, the constitutive equations can be represented as
\begin{equation}\label{Tau_and_G}
\tau (s) = \hat{G}(s)\gamma (s),
\qquad
\gamma(s) = \hat{J}(s) \tau(s)
\end{equation}
where $\hat{G}(s)$ and $\hat{J}(s)$ denote the shear operational modulus and the shear creep operational modulus, respectively.

From Eqs.~\eqref{HookLaw_LaplaceDOmain} and~\eqref{Tau_and_G}, two fundamental relationships in linear viscoelasticity can be derived, namely
\begin{equation}\label{Gs_Times_Js}
\hat{E}(s)\hat{D}(s)=1,
\qquad
\hat{G}(s)\hat{J}(s)=1.
\end{equation}

As stated in~\cite{Tschoegl1989}, Eqs.~\eqref{HookLaw_LaplaceDOmain}, \eqref{Tau_and_G}, and~\eqref{Gs_Times_Js} share the same algebraic structure as the constitutive equations for elastic materials. These relations establish the foundation of the theory of linear viscoelastic behavior.

Retransforming Eqs.~\eqref{HookLaw_LaplaceDOmain} and~\eqref{Tau_and_G} into the time domain leads to Boltzmann’s superposition principle, Eq.~\eqref{Boltzmann}, where the dependence of stress (or strain) on the history of strain (or stress) becomes evident. Accordingly, the general form of the constitutive equations for homogeneous, isotropic, and linear VEMs can be expressed as
\begin{equation}\label{Boltzmann}
\begin{aligned}
	& \sigma (t)=\int_{0}^{t}E(t-\xi) \frac{\dif \varepsilon (\xi)}{\dif \xi}\dif \xi,
\\	& \varepsilon (t)=\int_{0}^{t}D(t-\xi) \frac{\dif \sigma (\xi)}{\dif \xi}\dif \xi,
\\	& \tau (t)=\int_{0}^{t}G(t-\xi) \frac{\dif \gamma (\xi)}{\dif \xi}\dif \xi,
\\	& \gamma (t)=\int_{0}^{t}J(t-\xi) \frac{\dif \tau (\xi)}{\dif \xi}\dif \xi.
\end{aligned}
\end{equation}
In this case, $E(t)$ and $D(t)$ represent the Young’s relaxation modulus and the creep compliance, respectively, for a uniaxial tension test. In the Laplace domain, these moduli are related to their corresponding Young’s and creep operational moduli as follows~\cite{ParkShapery, Tschoegl1989}:
\begin{equation}\label{Model_COnstitutive1}
\begin{aligned}
	& \hat{E} (s)=sE(s),
\\	& \hat{D} (s)=sD(s).
\end{aligned}
\end{equation}
Furthermore, $G(t)$ and $J(t)$ are the shear relaxation and the shear creep compliance, respectively, defined in the time domain, and are related to their respective operational moduli, in Laplace domain, in the following form~\cite{ParkShapery, Tschoegl1989}:
\begin{equation}\label{Model_COnstitutive2}
\begin{aligned}
	& \hat{G} (s)=sG(s),
\\	& \hat{J} (s)=sJ(s).
\end{aligned}
\end{equation}

\subsection{Poisson's ratio in Laplace domain}

References~\cite{DIG, Tschoegl1989} show that, by relating $\hat{E}(s)$ and $\hat{G}(s)$, it is possible to define a function, referred to here as the operational Poisson’s ratio, which is expressed as
\begin{equation}\label{Model_COnstitutive}
\hat{\upsilon}(s) = \frac{\hat{E}(s)}{2\hat{G}(s)} - 1
\qquad \text{or} \qquad
\hat{\upsilon}(s) = \frac{\hat{E}(s)\hat{J}(s)}{2} - 1
\end{equation}
where $\hat{\upsilon}(s) = s\upsilon (s)$, with $\upsilon (s)$ being the Laplace transform of the time-dependent Poisson’s ratio, $\upsilon (t)$. Accordingly, from Eqs.~\eqref{Model_COnstitutive1}, \eqref{Model_COnstitutive2}, and~\eqref{Model_COnstitutive}, one obtains
\begin{equation}\label{v_s_sEsJs}
\upsilon (s) = \frac{sE(s)J(s)}{2} - \frac{1}{s}.
\end{equation}
Applying the inverse Laplace transform (ILT) to Eq.~\eqref{v_s_sEsJs} yields the time-dependent Poisson’s ratio, as expressed in Eq.~\eqref{v_t}. Nevertheless, when Eq.~\eqref{v_t} is employed to indirectly determine the Poisson’s ratio using complex constitutive models, numerical methods are required. This approach can become computationally expensive and may lead to inaccurate results~\cite{Chen,Tschoegl}. In this context, the following section introduces a methodology to obtain the Poisson’s ratio in a simpler and more efficient manner.

\subsection{Fractional viscoelasticity}

\subsubsection{Young and shear constitutive models in time domain}

According to~\cite{Mainardi, Bagley, Galucio}, several constitutive models involving fractional derivatives have been proposed, such as the Maxwell, Kelvin--Voigt, and Zener models. Among these, various authors~\cite{Ciniello, Pritz2003, Balbino, Sousa2017, Sousa2024} have shown that the four-parameter fractional Zener model is effective in predicting the behavior of VEMs in both time and frequency domains. The differential equation of this rheological model, which relates the stress and strain histories, can be written as
\begin{equation}\label{eq_contitutive}
\sigma (t)+\tau_{E}^{\beta_{E}}\frac{\Dif{\beta_E}}{\dif t^{\beta_{E}}}\sigma (t)=E_{\infty}\varepsilon (t)+E_{0}\tau_{E}^{\beta_{E}}\frac{\Dif{\beta_E}}{\dif t^{\beta_{E}}}\varepsilon (t),
\end{equation}
where $\tau_{E}$ is the relaxation time of the Young’s relaxation modulus, $E_{\infty}$ is the equilibrium modulus, and $E_{0}$ is the instantaneous modulus. In the International System of Units (SI), the variables $E_{\infty}$ and $E_{0}$ have units Pa, $\tau_{E}$ has unit seconds.

Eq.~\eqref{eq_contitutive} is written in terms of uniaxial tension tests. This equation has a corresponding form considering a pure shear test, which can be presented as
\begin{equation}\label{eq_contitutive2}
\tau (t) + \tau_{G}^{\beta_{G}} \frac{\Dif{\beta_G}}{\dif t^{\beta_{G}}} \tau (t) = G_{\infty}\,\gamma (t) + G_{0}\,\tau_{G}^{\beta_{G}} \frac{\Dif{\beta_G}}{\dif t^{\beta_{G}}} \gamma (t)
\end{equation}
where $\tau_{G}$ denotes the relaxation time of the shear modulus, $G_{\infty}$ is the equilibrium modulus, and $G_{0}$ is the instantaneous modulus. Furthermore, in the SI system, $G_{\infty}$ and $G_{0}$ are expressed in Pa, while $\tau_{G}$ is expressed in seconds.

\subsubsection{Young's and shear operational moduli and operational Poisson's ratio}

Considering $t>0$, the application of the Laplace transform (LT) to Eq.~\eqref{eq_contitutive} provides the steady-state response, which can be expressed as
\begin{equation}\label{Ehat_s}
\hat{E}(s)=\frac{\sigma (s)}{\varepsilon (s)}=\frac{E_{\infty}+E_{0}\tau_{E}^{\beta_{E}}s^{\beta_{E}}}{1+\tau_{E}^{\beta {E}}s^{\beta{E}}}.
\end{equation}
In this expression, $\hat{E}(s)$ denotes the Young’s operational modulus in tension~\cite{Tschoegl1989, Mainardi}. Analogously, the shear operational modulus, $\hat{G}(s)$, is obtained by applying the LT to Eq.~\eqref{eq_contitutive2}, yielding
\begin{equation}\label{Ghat_s}
\hat{G}(s)=\frac{\tau (s)}{\gamma (s)}=\frac{G_{\infty}+G_{0}\tau_{G}^{\beta_{G}}s^{\beta_{G}}}{1+\tau_{G}^{\beta {G}}s^{\beta{G}}}.
\end{equation}
Furthermore, the operational Poisson’s ratio can be obtained by combining Eqs.~\eqref{Model_COnstitutive}, \eqref{Ehat_s}, and~\eqref{Ghat_s}, which results in
\begin{equation}\label{vHATs_2}
\hat{\upsilon}(s)=\frac{E_{\infty}}{2G_{\infty}}\bracks*{\frac{(r_{E}\tau_{E}^{\beta_{E}} \tau_{G}^{\beta_{G}})s^{\beta_{E}+\beta_{G}}+(r_{E}\tau_{E}^{\beta_{E}})s^{\beta_{E}}+(\tau_{G}^{\beta_{G}})s^{\beta_{G}}+1}{(r_{G}\tau_{E}^{\beta_{E}} \tau_{G}^{\beta_{G}})s^{\beta_{E}+\beta_{G}}+(\tau_{E}^{\beta_{E}})s^{\beta_{E}}+(r_{G} \tau_{G}^{\beta_{G}})s^{\beta_{G}}+1} }-1
\end{equation}
where $r_{E}=\frac{E_0}{E_\infty}$ and $r_{G}=\frac{G_0}{G_\infty}$. With this last expression, it becomes possible to establish a framework for modeling the time-dependent Poisson’s ratio.

\subsubsection{Proposed model for the time-dependent Poisson's ratio}

One way to obtain the time-dependent Poisson's ratio is by applying the ILT to Eq.~\eqref{vHATs_2}. However, obtaining the ILT of the term in brackets is not a trivial task. Therefore, such a term is expanded into partial fractions as follows:
\begin{equation}\label{vHATsC}
\hat{\upsilon}(s)=\frac{E_{\infty}}{2G_{\infty}}\bracks*{A + \frac{Bs^{\beta_{E}}}{\tau_{E}^{\beta_{E}}s^{\beta_{E}}+1}+\frac{Cs^{\beta_{G}}}{{r_{G}\tau_{G}^{\beta_{G}}s^{\beta_{G}}}+1} +\frac{Ds^{\beta_{E}+\beta_{G}}}{(\tau_{E}^{\beta_{E}}s^{\beta_{E}}+1)(r_{G}\tau_{G}^{\beta_{G}}s^{\beta_{G}}+1)}}-1.
\end{equation}

From Eqs.~\eqref{vHATs_2} and~\eqref{vHATsC}, one can obtain $A=1$, $B=\tau_{E}^{\beta_{E}}(r_{E}-1)$, $C=\tau_{G}^{\beta_{G}}(1-r_{G})$, and $D = \tau_{E}^{\beta_{E}} \tau_{G}^{\beta_{G}} (r_{G} - 1)(r_{E} - 1)$. Since $\hat{\upsilon}(s)=s\upsilon (s)$, Eq.~\eqref{vHATsC} can be rewritten as
\begin{multline}\label{vHATsD}
\upsilon (s) = \frac{E_\infty}{2\,G_\infty} \left[A\parens[\Bigg]{\frac{1}{s}} + \frac{B}{\tau_E^{\beta_E}}\parens[\Bigg]{\frac{s^{{\beta_E} - 1}}{(s^{\beta_E} + 1 / \tau_E^{\beta_E})}} + \frac{C}{r_G\tau_G^{\beta_G}}\parens[\Bigg]{\frac{s^{{\beta_G} - 1}}{\parens[\big]{s^{\beta_G} + 1 / (r_G\tau_G^{\beta_G})}}} + \dotsb \right.
\\	\left. \dotsb + \frac{D}{\tau_E^{\beta_E} r_G \tau_G^{\beta_G}} \parens[\Bigg]{\frac{s^{{\beta_E} + {\beta_G} - 1}}{(s^{\beta_E} + 1 / \tau_E^{\beta_E})\parens[\big]{s^{\beta_G} + 1 / {(r_G\tau_G^{\beta_G})}}}}\right] - \parens[\Bigg]{\frac{1}{s}}.
\end{multline}

Following the theory presented by~\cite{Das, Atanacković}, and applying the ILT to Eq.~\eqref{vHATsD}, one obtains the time-dependent Poisson's ratio:
\begin{multline}\label{vHATsC2}
\upsilon (t) = \frac{E_\infty}{2\,{G_\infty}}\left[A + \frac{B}{\tau_E^{\beta_E}}{E_{\beta_E}}\parens[\Bigg]{- \frac{t^{\beta_E}}{\tau_E^{\beta_E}}} + \frac{C}{{r_G}\tau_G^{\beta_G}}{E_{\beta_G}} \parens[\Bigg]{- \frac{t^{\beta_G}}{{r_G}\tau_G^{\beta_G}}} + \dotsb \right.
\\	\left. \dotsb + \frac{D}{\tau_E^{\beta_E}{r_G}\,\,\tau_G^{\beta_G}} {\int_0^t \parens[\Bigg]{E_{\beta_E} \parens[\Bigg]{ - \frac{(t - \xi)^{\beta_E}}{\tau_E^{\beta_E}}} {E_{\beta_G}} \parens[\Bigg]{- \frac{\xi ^{\beta_G}}{{r_G}\tau_G^{\beta_G}}}} \dif \xi}\right] - 1,
\end{multline}
where $E_{\beta(\void)}(\void)$ is the Mittag-Leffler function of order $\beta(\void)$. It is observed that Eq.~\eqref{vHATsC2} requires integrating the product of two Mittag-Leffler functions with different differentiation orders, whose analytical solution is difficult to obtain. A possible approach is numerical integration, although this is computationally costly. In practice, however, fractional-order viscoelastic functions, considering uniaxial and shear tests, often exhibit very close differentiation orders, i.e., $\beta =\beta_{E}=\beta_{G}$~\cite{Chen2017, Sousa2018}. Thus, Eq.~\eqref{vHATs_2} can be rewritten as
\begin{equation}\label{vHAT4}
\hat \upsilon (s) = \frac{E_\infty}{2\,{G_\infty}} \bracks*{\frac{({{r_E}\tau_E^\beta \tau_G^\beta}){s^{2\beta}} + ({{r_E}\tau_E^\beta + \tau_G^\beta})\,{s^\beta} + 1}{{({\tau_E^\beta \,{s^\beta} + 1})({{r_G}\tau_G^\beta {s^\beta} + 1})}}} - 1.
\end{equation}
Expanding Eq.~\eqref{vHAT4} into partial fractions, one obtains
\begin{equation}\label{vHAT5}
\hat \upsilon (s) = \frac{E_\infty}{{2\,{G_\infty}}}\braces*{ {1 + A'\frac{{{s^\beta}}}{{({\tau_E^\beta {s^\beta} + 1})}} + B'\frac{{{s^\beta}}}{{({{r_G}\tau_G^\beta {s^\beta}\, + 1})}}} } - 1,
\end{equation}
where the coefficients $A'$ and $B'$ are given by
\begin{equation}
\begin{aligned}
	& A' = \frac{({\tau_G^\beta \tau_E^\beta {r_E} - \tau_E^{2\beta}{r_E} - \tau_G^\beta \tau_E^\beta + \tau_E^{2\beta}})}{({\tau_G^\beta {r_G} - \tau_E^\beta})},
\\	& B' = - \tau_G^\beta \,\frac{({\tau_G^\beta r_G^2 - \tau_E^\beta {r_E}{r_G} - \tau_G^\beta {r_G} + \tau_E^\beta {r_E}})}{({\tau_G^\beta {r_G} - \tau_E^\beta})}.
\end{aligned}
\end{equation}
As $\hat{\upsilon}(s)=s\upsilon (s)$, Eq.~\eqref{vHAT5} can be rearranged and expressed as
\begin{equation}\label{vHAT6}
\upsilon (s) = \frac{E_\infty}{{2\,{G_\infty}}}\braces*{{\frac{1}{s} + \parens[\Bigg]{\frac{{A'}}{{\tau_E^\beta}}} \frac{{{s^{\beta - 1}}}}{{({{s^\beta} + 1/\tau_E^\beta})}} + \parens[\Bigg]{\frac{{B'}}{{r_G}\tau_G^\beta}} \frac{{{s^{\beta - 1}}}}{{\parens[\Big]{{s^\beta} + \frac{1}{{r_G}\tau_G^\beta}}}}} } - 1.
\end{equation}
Following the approach presented in~\cite{Das, Atanacković}, and applying the ILT to Eq.~\eqref{vHAT6}, the time-dependent Poisson’s ratio can be written as
\begin{equation}\label{v_time}
\upsilon (t) = \frac{{E_\infty}}{{2\,G_\infty}}
\bracks*{1 + \parens[\Bigg]{\frac{A'}{\tau_E^\beta}} E_\beta \parens[\Bigg]{- \frac{t^\beta}{\tau_E^\beta}} + \parens[\Bigg]{\frac{B'}{r_G \tau_G^\beta}} E_\beta \parens[\Bigg]{- \frac{t^\beta}{r_G \tau_G^\beta}}
} - 1.
\end{equation}
Equation~\eqref{v_time} represents the proposed model to describe the time-dependent Poisson’s ratio. This viscoelastic function is identified and discussed in the results and discussion section. Using the Mittag-Leffler relations~\cite{Mainardi}, when $t \rightarrow \infty$, the equilibrium Poisson’s ratio, $\upsilon_{\infty}$, is obtained as
\begin{equation}\label{v_oo}
\upsilon_\infty = \lim_{t \to \infty} \upsilon (t) = \frac{E_\infty}{{2\,{G_\infty}}} - 1.
\end{equation}
On the other hand, for $t \rightarrow 0$, the limit of Eq.~\eqref{v_time} becomes
\begin{equation}\label{Lim_v_o}
\lim_{t \to 0} \upsilon (t) = \frac{{{E_0}}}{{2\,{G_0}}} - 1.
\end{equation}
This relation allows the definition of the instantaneous Poisson’s ratio, $\upsilon_{0}$, as
\begin{equation}\label{v_o}
\upsilon_0= \frac{{{E_0}}}{{2\,{G_0}}} - 1,
\end{equation}
where $E_{0}$ and $G_{0}$ are the instantaneous moduli values for $t \rightarrow 0$.

\subsubsection{Complex Young’s and shear moduli and the complex Poisson's ratio}

According to~\cite{ParkShapery, Tschoegl1989}, the complex viscoelastic functions in the frequency domain are related to the operational functions as follows:
\begin{equation}\label{E_G_v_omega}
E^*(\Omega) = \restr{\hat E(s)}{s \to i\Omega},
\qquad G^*(\Omega) = \restr{\hat G(s)}{s \to i\Omega},
\qquad \upsilon^*(\Omega) = \restr{\hat \upsilon(s)}{s \to i\Omega},
\end{equation}
where $E^*(\Omega)$, $G^*(\Omega)$, and $\upsilon^*(\Omega)$ are defined as the complex Young’s modulus, the complex shear modulus, and the complex Poisson's ratio, respectively. From Eq.~\eqref{Ehat_s}, the complex Young’s modulus can be written as
\begin{equation}\label{E_omega}
E^*(\Omega) = \frac{\bar E_{0} + \bar E_{\infty} (i\,\Omega \tau_E)^\beta}{1 + (i\,\Omega \tau_E)^\beta},
\end{equation}
in which $\bar E_{0} = E_{\infty}$ and $\bar E_{\infty} = E_{0}$. Analogously, the complex shear modulus is obtained from Eq.~\eqref{Ghat_s} as
\begin{equation}\label{G_omega}
G^*(\Omega) = \frac{\bar G_{0} + \bar G_{\infty} (i\,\Omega \tau_G)^\beta}{1 + (i\,\Omega \tau_G)^\beta},
\end{equation}
where $\bar G_{0} = G_{\infty}$ and $\bar G_{\infty} = G_{0}$.

These moduli have real and imaginary components that represent storage and loss of energy, respectively, given by
\begin{equation}\label{G_E_omega}
E^*(\Omega) = E'(\Omega) + iE''(\Omega),
\qquad G^*(\Omega) = G'(\Omega) + iG''(\Omega),
\end{equation}
where $E'(\Omega)$, $G'(\Omega)$, $E''(\Omega)$, and $G''(\Omega)$ are defined as the storage Young’s modulus, the storage shear modulus, the loss Young’s modulus, and the loss shear modulus, respectively.

Furthermore, the complex Poisson's ratio, following Eq.~\eqref{vHAT4}, can be expressed as
\begin{equation}\label{v_omega}
\upsilon^*(\Omega) = \frac{\bar E_{0}}{2\,\bar G_{0}}
\bracks*{
\frac{
(r_E \tau_E^\beta \tau_G^\beta) (i\Omega)^{2\beta} + (r_E \tau_E^\beta + \tau_G^\beta) (i\Omega)^\beta + 1}{
(r_G \tau_E^\beta \tau_G^\beta) (i\Omega)^{2\beta} + (r_G \tau_G^\beta + \tau_E^\beta) (i\Omega)^\beta + 1}
} - 1,
\end{equation}
where $\upsilon^*(\Omega) = \upsilon'(\Omega) - i\upsilon''(\Omega)$, in which $\upsilon'(\Omega)$ represents the storage Poisson's ratio and $\upsilon''(\Omega)$ the loss component. Relating these components, one can define the Poisson's loss factor, $\eta_\upsilon(\Omega)$, as
\begin{equation}\label{eta_v_omega}
\eta_\upsilon(\Omega) = \frac{\upsilon''(\Omega)}{\upsilon'(\Omega)}.
\end{equation}
Graphical representations of these functions, namely the Poisson’s loss factor and the corresponding storage modulus, are presented in the section of results and discussions.

\subsubsection{Temperature effects on viscoelastic behavior across time and frequency domains}

Temperature is a factor that strongly influences polymer behavior. According to~\cite{Lakes}, the viscoelastic behavior can be predicted using the Time-Temperature Superposition Principle (TTSP), which establishes that the viscoelastic response can be shifted along the time (or frequency) scale. In practice, this means that each record of time $t$ or frequency $\Omega$ at a given temperature $T$ can be shifted to a reference temperature $T_{0}$, as follows:
\begin{equation}\label{t_r}
t_{r} = \frac{t}{\alpha_T(T, T_0)}
\qquad \text{or} \qquad
\Omega_{r} = \alpha_T(T, T_0)\,\Omega,
\end{equation}
where $t_{r}$ is the reduced time (in seconds) and $\Omega_{r}$ is the reduced frequency (in Hertz) \cite{DIG, Lakes, BrinsonBrinson, Capodagli}. In this formulation, $\alpha_T(T, T_0)$ is the temperature shift factor.

In their study, \cite{Landel} proposed a mathematical model for this factor, classically known as the Williams--Landel--Ferry (WLF) equation, expressed as
\begin{equation}\label{log_alpha_T}
\log {\alpha_T}(T, T_0) = \frac{-C_{1}^{T}(T - T_0)}{C_{2}^{T} + (T - T_0)},
\end{equation}
where $C_{1}^{T}$ and $C_{2}^{T}$ are material-dependent constants characteristic of each polymer.

\subsubsection{Viscoelastic functions considering temperature effects}

As discussed by~\cite{Ciniello, Mainardi2011}, considering the fractional Zener model, Eq.~\eqref{eq_contitutive}, the time-dependent Young's relaxation modulus, accounting for temperature effects, can be expressed as
\begin{equation}\label{E_tr}
E(t_r) = E_\infty {\bracks*{1 + r_E \, E_\beta \parens[\Bigg]{-\parens[\Bigg]{\frac{t}{\tau_E \, \alpha_T(T, T_0)}}^{\beta}}}},
\end{equation}
where $E_\beta(\void)$ is the Mittag-Leffler function of order $\beta$.

Similarly, from Eq.~\eqref{eq_contitutive2}, the time-dependent shear relaxation modulus, including temperature effects, can be written as
\begin{equation}\label{G_tr}
G(t_r) = G_\infty \bracks*{1 + r_G \, E_\beta \parens[\Bigg]{-\parens[\Bigg]{\frac{t}{\tau_G \, \alpha_T(T, T_0)}}^{\beta}}}.
\end{equation}
Analogously, the proposed model for the time-dependent Poisson’s ratio, Eq.~\eqref{v_time}, considering temperature effects, becomes
\begin{equation}\label{v_tr}
\upsilon(t_r) = \frac{E_\infty}{2\,G_\infty} \bracks*{1 + \parens[\Bigg]{\frac{A'}{\tau_E^\beta}} E_\beta \parens[\Bigg]{-\parens[\Bigg]{\frac{t}{\tau_E \, \alpha_T(T, T_0)}}^\beta} + \parens[\Bigg]{\frac{B'}{r_G \, \tau_G^\beta}} E_\beta \parens[\Bigg]{-\frac{1}{r_G} \parens[\Bigg]{\frac{t}{\tau_G \, \alpha_T(T, T_0)}}^\beta}} - 1.
\end{equation}
Furthermore, the complex Young's modulus, considering temperature effects, is
\begin{equation}\label{Eomega_r}
E^*(\Omega_r) =
\frac{\bar{E}_0 + \bar{E}_\infty (i \, \Omega_r \tau_E)^\beta} {1 + (i \, \Omega_r \tau_E)^\beta}.
\end{equation}
Similarly, the complex shear modulus is
\begin{equation}\label{Gomega_r}
G^*(\Omega_r) =
\frac{\bar{G}_0 + \bar{G}_\infty (i \, \Omega_r \tau_G)^\beta} {1 + (i \, \Omega_r \tau_G)^\beta}.
\end{equation}
Analogously, the complex Poisson's ratio, including temperature effects, can be expressed as
\begin{equation}\label{v_omega_r}
\upsilon^*(\Omega_r) =
\frac{\bar{E}_0}{2 \bar{G}_0}
\bracks*{
\frac{
\parens{r_E \tau_E^\beta \tau_G^\beta}\,(i \Omega_r)^{2\beta} + \parens{r_E \tau_E^\beta + \tau_G^\beta}\,(i \Omega_r)^\beta + 1}{
\parens{r_G \tau_E^\beta \tau_G^\beta}\,(i \Omega_r)^{2\beta} + \parens{r_G \tau_G^\beta + \tau_E^\beta}\,(i \Omega_r)^\beta + 1}
} - 1.
\end{equation}
The viscoelastic functions presented in this section form a basis for predicting the mechanical behavior of isotropic, linear, and thermorheologically simple viscoelastic materials in both the time and frequency domains. By interconversion, the corresponding viscoelastic functions can be obtained for either domain.

\section{Materials and methods}

\subsection{Artificial experimental data in time domain}

Figure~\ref{fig1} presents the time-dependent relaxation moduli for Young’s modulus and shear modulus at different temperatures. The curves of $E(T, t)$ and $G(T, t)$ reveal the progressive reduction of both moduli with time, characterizing the viscoelastic nature of the material and highlighting the pronounced influence of temperature on its mechanical response.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{Fig1_left.jpg}
\caption{}\label{fig1_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{Fig1_right.jpg}
\caption{}\label{fig1_b}
\end{subfigure}
\caption{Time-dependent relaxation moduli of \eqref{fig1_a}~Young’s modulus and \eqref{fig1_b}~shear modulus at various temperatures, illustrating the material’s viscoelastic response and its pronounced dependence on temperature.}\label{fig1}
\end{figure}

Starting from the viscoelastic parameters reported by~\cite{Sousa2018} for the EAR\textsuperscript{\textregistered}-C1002 material, sets of synthetic experiments were generated for the time-dependent Young's and shear relaxation moduli. Following the approach adopted by Saprunov, Gergesova, and Emri~\cite{Saprunov2013}, artificial datasets were created from analytical viscoelastic models and subsequently perturbed with random noise to emulate experimental uncertainty.

To make the artificial experiments more realistic, uncorrelated multiplicative noise was introduced at each time-temperature point. The noise was generated using MATLAB’s \texttt{rand} function, which produces uniformly distributed random numbers in the interval $(0,1)$. In order to obtain a zero-mean symmetric perturbation, the transformation $2\,\texttt{rand} - 1$ was applied, expanding the interval to $(0,2)$ and subsequently shifting it to $(-1,1)$. The resulting variable was then multiplied by $10\%$, yielding a bounded uniform perturbation in the interval $[-10\%, +10\%]$ of the computed relaxation modulus value.

Thus, for each time and temperature point, the pseudo-experimental relaxation modulus was defined as
\begin{equation}
M_{\mathrm{noisy}}(t,T) = M_{\mathrm{model}}(t,T)(1 + \varepsilon),
\end{equation}
where \(\varepsilon \sim \mathcal{U}(-0.10,\,0.10)\), and \(M_{\mathrm{model}}(t,T)\) is the theoretical relaxation modulus without noise.

The injected noise is therefore uniformly distributed, multiplicative, statistically independent between sampling points, and non-correlated (white noise), with no prescribed spectral bandwidth.

The synthetic experimental data were generated over a time interval consistent with typical physical tests, ranging from $1$ to $3600$~seconds. Similarly, the artificial experimental temperatures were selected within the elastic limits of the material (from \qty{-35}{\degreeCelsius} to \qty{60}{\degreeCelsius}), ensuring proper coverage of the transition regions.

The procedure of validating computational methods through artificial experimental data has also been widely employed in the literature. For instance, Shanbhag~\cite{Shanbhag2023} used synthetic data generated from known Prony-series representations of creep compliance and stress relaxation functions to test and verify the stability of the interconversion algorithm implemented in the PSI (\emph{Prony Series Interconversion}) program. In that study, artificial datasets spanning more than 20 decades in time were generated, enabling quantitative comparison between the original analytical model and the reconstructed viscoelastic functions. This approach ensured that numerical accuracy, stability, and the sensitivity of the interconversion procedure could be assessed independently of experimental uncertainties, following a methodology similar in spirit to that of Saprunov, Gergesova, and Emri~\cite{Saprunov2013}.

\subsection{Methodology for material property identification}

Each experimental point of the Young's relaxation modulus can be represented by $E^{\exp}(t_k, T_j)$ at the $j$-th temperature, where $1 \le j \le n_T$, and at the $k$-th time, with $1 \le k \le n_P(T_j)$. Here, $n_P(T_j)$ denotes the number of experimental time points at the $j$-th temperature, and $n_T$ is the total number of temperatures. Given the constitutive material model $E(t_k, T_j)$, the $kj$-th error component in the time domain, denoted as $\Delta E_{t_{kj}}$, is defined as
\begin{equation}\label{Error_E}
(\Delta E_{t_{kj}})^2 =
\bracks*{
\frac{E^{\exp}(t_k, T_j) - E(t_k, T_j)} {E^{\exp}(t_k, T_j)}
}^2,
\end{equation}
where the superscript ``exp'' indicates experimental data. The overall quadratic error for all curves of the Young's relaxation modulus, $\Delta E_{t_R}$, is then given by
\begin{equation}\label{Error_E_square}
(\Delta E_{t_R})^2 =
\frac{1}{n_T} \, \frac{1}{n_P(T_j)}
\sum_{j=1}^{n_T} \sum_{k=1}^{n_P(T_j)}
(\Delta E_{t_{kj}})^2.
\end{equation}

Similarly, in shear relaxation tests, each experimental point is represented by $G^{\exp}(t_k, T_j)$. Given the constitutive model $G(t_k, T_j)$, the $kj$-th error component in the time domain, denoted as $\Delta G_{t_{kj}}$, is defined as
\begin{equation}\label{Error_G}
(\Delta G_{t_{kj}})^2 =
\bracks*{
\frac{G^{\exp}(t_k, T_j) - G(t_k, T_j)} {G^{\exp}(t_k, T_j)}
}^2.
\end{equation}
The resulting quadratic error for all shear relaxation curves in the time domain, $\Delta G_{t_R}$, is then
\begin{equation}\label{Error_G_square}
(\Delta G_{t_R})^2 =
\frac{1}{n_T} \, \frac{1}{n_P(T_j)}
\sum_{j=1}^{n_T} \sum_{k=1}^{n_P(T_j)}
(\Delta G_{t_{kj}})^2.
\end{equation}
Consequently, the global quadratic error in the time domain, considering both Young's and shear relaxation tests, $(\Delta \EG_{t_{\tot}})^2$, is defined as
\begin{equation}\label{Error_Total_E_G_square}
(\Delta \EG_{t_{\tot}})^2 =
\frac{
(\Delta E_{t_R})^2 +
(\Delta G_{t_R})^2}{2}.
\end{equation}

The viscoelastic parameters of Young's and shear moduli are identified through a single integrated hybrid optimization process, in which common parameters are obtained for both moduli. Following~\cite{Chen2017, Sousa2018, Mainardi2011, Lakes2006, Waterman}, the influence of temperature and the fractional differentiation orders are assumed to be the same for both Young’s and shear relaxation moduli. Moreover, as noted by~\cite{Pritz, Tschoegl1989, Waterman}, the time-dependent Poisson's ratio of a viscoelastic material is physically meaningful only when it lies between 0 and 0.5 and increases monotonically over time.

Accordingly, the standard constrained optimization problem, denoted as $P_{t}$, for the four-parameter fractional Zener model, is formulated as
\begin{equation}\label{P_t}
P_t \colon
\begin{system}
	\text{minimize} \quad (\Delta \EG_{t_{\tot}})^2(\mathbf{x}) \colon \mathbb{R}^9 \to \mathbb{R},
\\	\text{where} \quad \mathbf{x} = \parens[\big]{E_\infty, E_0, \tau_E, G_\infty, G_0, \tau_G, \beta, C_1^T, C_2^T },
\\	\text{subject to:}
\\
\begin{system}
	\mathbf{x}^{\inf} \le \mathbf{x} \le \mathbf{x}^{\sup},
\\	0 < \upsilon(t) < 0.5,
\\	\frac{\dif \upsilon(t)}{\dif t} > 0,
\end{system}
\end{system}
\end{equation}
where $\mathbf{x}^{\inf}$ and $\mathbf{x}^{\sup}$ denote the lower and upper bounds of the design variables, respectively. It is assumed that the same fractional differentiation order, $\beta$, applies to both moduli, and that temperature effects influence Young's and shear moduli equivalently.

In order to provide a quantitative measure of the agreement between the reference analytical model and the identified curves, the relative root mean square error (RRMSE) is also evaluated. This metric is defined as
\begin{equation}
\mathrm{RRMSE} = \sqrt{ \frac{1}{N} \sum_{k=1}^{N} \parens[\Bigg]{ \frac{y_{\mathrm{model}}(k) - y_{\mathrm{ref}}(k)}{y_{\mathrm{ref}}(k)}}^2 }
\end{equation}
where $y_{\mathrm{model}}$ and $y_{\mathrm{ref}}$ denote the identified and reference responses, respectively, and $N$ is the number of discrete points. This dimensionless metric allows for a normalized comparison independent of the magnitude of the analyzed quantity.

\subsection{Computational structure}

The VEM characterization process employs a hybrid optimization technique~\cite{Sousa2018}. The pseudocode algorithm for identifying the viscoelastic material under study is presented in Algorithm~\ref{tab:pseudocode_vem}. The implementation is carried out in the Matlab computational environment.

\begin{algorithm}[h!]
\caption{Pseudocode algorithm for identifying the VEM}\label{tab:pseudocode_vem}
\begin{enumerate}[1.]
\item Acquire the experimental viscoelastic material data.
\item Define the lower and upper bounds of the design variables.
\item \label{tab:pseudocode_vem_3} For each hybrid optimization procedure $i$ ($i = 1, \dots, \NOP$, where NOP is the number of optimization runs):
\begin{enumerate}
\renewcommand{\theenumii}{\ref{tab:pseudocode_vem_3}\arabic{enumii}.}
\renewcommand{\labelenumii}{\theenumii}
\item Approximate the global minimum design vector using GA, $X_{\GA}(i)$:
\begin{enumerate}
\renewcommand{\theenumiii}{\alph{enumiii})}
\renewcommand{\labelenumiii}{\theenumiii}
\item Define the GA parameters for Matlab\textsuperscript{\textregistered}.
\item Perform optimization using the GA.
\item Obtain a vector close to the global optimum, $X_{\GA}$.
\end{enumerate}
\item Refine the global minimum design using NLP:
\begin{enumerate}
\renewcommand{\theenumiii}{\alph{enumiii})}
\renewcommand{\labelenumiii}{\theenumiii}
\item Define the NLP parameters for Matlab\textsuperscript{\textregistered}.
\item Obtain the refined vector $X_{\NLP}(i)$ (starting from $X_{\GA}(i)$).
\end{enumerate}
\end{enumerate}
\item Determine the ``global optimal point'' as the best vector among $X_{\NLP}$ results.
\end{enumerate}
\end{algorithm}

Initially, a global search of the admissible parameter space is performed using a Genetic Algorithm (GA), aiming to mitigate the effects of multiple local minima typically associated with fractional viscoelastic models. Due to the strong nonlinearity and parameter coupling inherent to such formulations, purely gradient-based methods may exhibit sensitivity to the initial guess.

The parameter vector obtained from the GA stage is subsequently used as the initial point in a Nonlinear Programming (NLP) procedure, which provides local refinement and high-accuracy convergence. In the NLP stage, a strict stopping criterion of $1\times10^{-11}$ is adopted to ensure numerical precision.

To enhance robustness and reduce the influence of the stochastic nature of GA, ten independent optimization runs are performed. Each GA execution employs 500 individuals, a mutation rate of 5\%, and 1000 generations. All runs converge, after the NLP refinement stage, to the same optimal solution within numerical tolerance, indicating stability and practical identifiability of the identified parameters within the prescribed bounds.

It is worth emphasizing that, in contrast to analytical or asymptotic identification approaches frequently adopted in the literature ---~such as methods based on closed-form $L_2$ norm minimization or asymptotic analyses, which may provide theoretical guarantees of uniqueness under restrictive assumptions~--- the present contribution focuses on a numerical identification framework tailored to highly nonlinear fractional viscoelastic models.

Furthermore, the use of synthetic datasets is a deliberate methodological choice aimed at isolating the intrinsic stability, robustness, and practical identifiability of the proposed algorithm under controlled constitutive conditions. The objective of this study is therefore to assess the computational performance and consistency of the identification strategy itself, rather than to provide experimental validation.

Table~\ref{tab:material_limits} presents the ranges of design variables for the fractional Zener model, considering temperature effects.

\begin{table}[h!]
\caption{Ranges of material property values employed in the optimization process.}\label{tab:material_limits}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}lcc@{}}
\toprule
\textbf{Variable} & \textbf{Nomenclature} & \textbf{Interval Limits} \\
\midrule
WLF material constant 1 & $C_1$ & $[0, 1000]$ \\
WLF material constant 2 & $C_2 \ [\unit{\degreeCelsius}]$ & $[0, 2000]$ \\
Equilibrium modulus & $E_\infty, G_\infty \ [\unit{\pascal}]$ & $[10^1, 10^7]$ \\
Instantaneous modulus & $E_0, G_0 \ [\unit{\pascal}]$ & $[10^1, 10^{10}]$ \\
Relaxation times & $\tau_E, \tau_G \ [\unit{\second}]$ & $[10^{-4}, 10^0]$ \\
Fractional derivative order & $\beta$ & $(0,1)$ \\
\bottomrule
\end{tabular}
\end{table}

\section{Results and discussions}

Based on the sets of artificial experimental data previously presented and by applying the proposed methodology, the Young's and shear relaxation moduli were identified. The numerical results are summarized in Table~\ref{tab:identificados}, under the column of identified properties.

\begin{table}[h!]
\caption{Identified numerical properties.}\label{tab:identificados}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}cc@{}}
\toprule
\textbf{Constants} & {\textbf{Identified properties}} \\
\midrule
$E_\infty$ [Pa]					& 2.1286675E+06 \\
$E_0$ [Pa]						& 1.7237572E+09 \\
$\tau_E$ [s]						& 1.7322984E-01 \\
$G_\infty$ [Pa]					& 7.1190339E+05 \\
$G_0$ [Pa]						& 6.2021921E+08 \\
$\tau_G$ [s]						& 1.7322984E-01 \\
$\beta$							& 4.5241761E-01 \\
$C_1$							& 1.4331363E+01 \\
$C_2$ [\unit{\degreeCelsius}]	& 1.0311494E+02 \\
\midrule
\textbf{NLP error} & \textbf{1.1500000E-02} \\
\bottomrule
\end{tabular}
\end{table}

A close inspection of Table~\ref{tab:identificados} reveals a clear separation between the instantaneous ($E_0$, $G_0$) and equilibrium moduli ($E_\infty$, $G_\infty$) elastic responses, which is characteristic of fractional Zener-type viscoelastic behavior. The significant ratio between the instantaneous and equilibrium moduli confirms the pronounced relaxation effect embedded in the artificial dataset. Moreover, the equality of the characteristic relaxation times $\tau_E$ and $\tau_G$ demonstrates that the identification procedure consistently recovered the imposed constitutive structure, indicating internal coherence of the proposed framework.

\begin{figure}[!h]
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{Fig2_left.jpg}
\caption{}\label{fig2_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
\includegraphics[width=\textwidth]{Fig2_right.jpeg}
\caption{}\label{fig2_b}
\end{subfigure}
\caption{Artificial experimental data and corresponding fitted models for the relaxation moduli of Young’s modulus and shear modulus~\eqref{fig2_a}, and for the time-dependent Poisson’s ratio~\eqref{fig2_b}.}\label{fig2}
\end{figure}

Regarding the time-dependent relaxation moduli (Young's and shear), the fitting results of the theoretical model to the experimental data are shown in Figure~\ref{fig2}\eqref{fig2_a}. The identification results reveal excellent agreement, with nearly perfect fits for both the Young's and shear relaxation moduli.

It is important to clarify that the fitted curves presented in Figure~\ref{fig2} were obtained exclusively using the identified parameters listed in Table~\ref{tab:identificados}, without any use of reference or input parameters from the artificial dataset. The quantitative agreement between fitted and experimental curves is measured through the nonlinear programming objective function (NLP error), whose final value is reported in Table~\ref{tab:identificados}. The obtained value ($1.15 \times 10^{-2}$) indicates a very low global fitting error, confirming the accuracy of the identification procedure.
No systematic deviation is observed along the time axis, indicating that the fractional formulation adequately captures the memory effects governing the relaxation process over the entire analyzed interval.

Furthermore, using the proposed model for the time-dependent Poisson’s ratio (Eq.~\eqref{v_tr}), the corresponding viscoelastic function in the time domain was constructed, as illustrated in Figure~\ref{fig2}\eqref{fig2_b}. The resulting curve exhibits a monotonically increasing behavior, in accordance with the theoretical framework discussed by~\cite{Tschoegl, Tschoegl1989}, representing a physically consistent response.

Figure~\ref{fig2}\eqref{fig2_b} compares the reference Poisson’s ratio curve (Sousa et~al.~\cite{Sousa2018}), used for synthetic data generation, with the curve obtained from the identification procedure considering 10\% noise. The close agreement between both curves demonstrates the robustness of the proposed identification strategy, indicating that the identified fractional parameters remain stable even under perturbed data conditions.

For the time-domain Poisson’s ratio shown in Figure~\ref{fig2}\eqref{fig2_b}, the relative root mean square error (RRMSE) between the reference model and the identified curve is $3.02 \times 10^{-2}$, indicating very good agreement between the two responses, even in the presence of 10\% noise.

It is worth emphasizing that the smooth transition between the instantaneous and equilibrium values of Poisson’s ratio is directly governed by the fractional exponent $\beta$. The identified value ($\beta \approx 0.45$) indicates an intermediate memory regime between purely elastic and classical exponential relaxation, highlighting the ability of the fractional model to describe gradual and non-exponential lateral deformation effects.

By applying interconversion to the frequency domain, the complex Poisson's ratio was obtained, as illustrated in Figure~\ref{fig3}. The results show that the storage Poisson's ratio follows a monotonically decreasing trend, while its loss factor exhibits a distinct maximum along the frequency axis. These outcomes are consistent with the classical theory of linear viscoelasticity~\cite{DIG, Pritz, Tschoegl, Tschoegl1989, TKE}.

In the frequency domain (Figure~\ref{fig3}), the RRMSE associated with the storage Poisson’s ratio is $3.01 \times 10^{-2}$, demonstrating consistent agreement with the time-domain results. For the loss factor, the root mean square error (RMS) was adopted due to near-zero values in the reference data. The obtained value of $2.43 \times 10^{-3}$ indicates excellent agreement between the identified and reference curves, confirming the accuracy of the proposed interconversion procedure in capturing dissipative effects.

Figure~\ref{fig3} also presents a comparison between the reference frequency-domain response and the response reconstructed from the parameters identified under noisy conditions. The small deviation observed near the transition region confirms that the interconversion procedure preserves the physical consistency of the viscoelastic response, even when the time-domain identification is performed with noisy input data.

The presence of a well-defined peak in the loss factor identifies the frequency band in which volumetric-deviatoric coupling leads to maximum energy dissipation. This behavior further confirms the thermodynamic admissibility and physical consistency of the identified parameter set.

\begin{figure}[t]
\includegraphics[width=.48\textwidth]{Fig3_left.jpeg}
\hfill
\includegraphics[width=.48\textwidth]{Fig3_right.jpeg}
\caption{Poisson's ratio in the frequency domain.}\label{fig3}
\end{figure}

Furthermore, by performing the interconversion of the viscoelastic model from the time domain to the frequency domain, the complex Young’s and shear moduli are obtained, as presented in Figure~\ref{fig:4}. This procedure enables the evaluation of the material’s viscoelastic behavior in terms of its frequency-dependent properties. The storage modulus $E'(\Omega)$ and the loss factor $\eta(\Omega)$ were computed based on the parameters identified from the time-domain response, thereby establishing a consistent link between the two representations of the viscoelastic material model.

\begin{figure}[t]
\includegraphics[width=.48\textwidth]{Fig4_left.jpg}
\hfill
\includegraphics[width=.48\textwidth]{Fig4_right.jpg}
\caption{Graphs of complex Young’s modulus in frequency: storage and loss factors.}\label{fig:4}
\end{figure}

As expected for viscoelastic materials, the storage modulus increases progressively with frequency, reflecting the stiffening effect observed at higher excitation rates~\cite{BrinsonBrinson}. Simultaneously, the loss factor exhibits a characteristic maximum corresponding to the transition region between predominantly elastic and predominantly viscous responses. The ability to reproduce both domains from a single identified parameter set demonstrates the internal coherence and robustness of the proposed interconversion framework.

\begin{figure}[t]
\includegraphics[width=0.5\linewidth]{Fig5.jpeg}
\caption{Wicket plot showing the relationship between storage modulus and loss factor for Young’s and shear responses over the analyzed frequency range. The curve provides a compact representation of the stiffness-damping interplay and highlights the frequency band associated with maximum energy dissipation.}\label{fig:5}
\end{figure}

In addition, the wicket plot shown in Figure~\ref{fig:5} provides a compact and insightful representation of the relationship between stiffness and damping across the frequency spectrum. This plot simultaneously displays the evolution of the storage modulus and the corresponding loss factor, allowing for a clear visualization of how energy storage and energy dissipation mechanisms interact. As observed, the storage modulus increases progressively with frequency, reflecting the stiffening behavior typical of viscoelastic materials at higher excitation rates. In contrast, the loss factor exhibits a pronounced peak that identifies the frequency range in which dissipative effects are most significant. This characteristic behavior is consistent with fractional-order viscoelastic models and further validates the physical consistency of the identified parameters. Overall, the wicket plot highlights the frequency band in which the material demonstrates optimal damping performance, thereby complementing the traditional frequency-domain curves.

Therefore, it can be concluded that the proposed method for the interconversion of viscoelastic models enables the determination of material properties according to the fractional Zener constitutive model. This is achieved while simultaneously considering both tensile and shear tests in the time domain and also employing a viscoelastic model for the time-dependent Poisson's ratio, followed by the interconversion to obtain the complex Poisson's modulus, the complex Young's modulus, and the complex shear modulus.

Additionally, repeated independent optimization runs converged to the same optimal solution within numerical tolerance, providing evidence of practical identifiability and numerical stability of the adopted hybrid strategy.

\section{Conclusion}

This study establishes an analytically consistent formulation for the viscoelastic Poisson’s ratio within the framework of the four-parameter fractional Zener model. Rather than introducing a new constitutive law, the main contribution lies in deriving an explicit representation of the time-dependent Poisson’s ratio directly from the fractional relaxation moduli (Young’s and shear), under the assumptions of isotropy, linear viscoelasticity, and thermorheological simplicity. By treating Poisson’s ratio as a primary viscoelastic function, the proposed formulation ensures theoretical coherence between constitutive parameters and avoids simplifying assumptions such as a constant bulk modulus.

A unified interconversion strategy was further implemented, enabling consistent transitions between time- and frequency-domain representations of the relaxation moduli and the Poisson’s ratio. This unified treatment reinforces the analytical equivalence between domains and facilitates the incorporation of the proposed formulation into numerical simulations and structural analyses involving viscoelastic materials.

The identification procedure was evaluated using synthetic experimental datasets for the Young’s and shear relaxation moduli of the EAR\textsuperscript{\textregistered}-C1002 material. A hybrid optimization strategy combining a Genetic Algorithm and Nonlinear Programming was employed to ensure global exploration and local refinement of the parameter space. The repeated convergence of independent optimization runs to the same solution indicates numerical stability and practical identifiability within the prescribed parameter bounds.

The obtained results demonstrate physically consistent behavior: in the time domain, the viscoelastic Poisson’s ratio exhibits a monotonically increasing trend, whereas in the frequency domain the storage Poisson’s ratio decreases monotonically and its associated loss factor presents a distinct maximum. These characteristics are consistent with established theoretical expectations for isotropic linear viscoelastic materials.

Overall, the proposed analytical framework contributes to the consistent modeling and indirect identification of the viscoelastic Poisson’s ratio within a fractional constitutive setting, offering a coherent basis for future applications in parameter identification, numerical modeling, and structural vibration analysis involving viscoelastic materials.

\printCOI

\section*{Supplementary material}

Supporting information for this article is available on the journal’s website under \printDOI\ or from the authors.

\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{Supplementary-material-ExpArt-Mod-E-time.xlsx}}
\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{Supplementary-material-ExpArt-Mod-G-time.xlsx}}

\printbibliography

\end{document}