Modeling elastoviscoplastic materials using physics-informed neural networks

Abstract

Elastoviscoplastic (EVP) materials pose significant challenges for predictive modeling, especially under large amplitude oscillatory shear (LAOS) conditions due to their nonlinear rheological behavior. While various models have been proposed to describe yielding and viscoelastic transitions, their application to experimental data remains limited due to computational complexity, poor generalizability, and difficulty in fitting noisy data. This work re-evaluates rheological modeling of EVP materials by using a physics-informed neural network (PINN) approach that embeds a modified Saramito model into the training framework. By directly fitting time-dependent stress data on typical EVP materials, this method circumvents the need for gradient estimation from noisy measurements and offers a data-efficient, differentiable, and generalizable alternative to conventional fitting methods. A shear-thinning formulation is introduced to reflect the decreasing viscosity at higher strain amplitudes. The proposed PINN approach is validated using both synthetic and experimental data, demonstrating stable recovery of physical parameters, enhanced interpretability, and improved predictive capability across a range of strain amplitudes. This framework bridges microstructural deformation modes with rheological modeling, offering a powerful tool for understanding and predicting nonlinear viscoelastic behavior in soft matter.

---

Introduction

EVP materials are encountered in many applications, including coating, droplet-based inkjet printing, and 3D printing,^1,2 where they are subjected to various stress-loading conditions. In medical studies and diagnostics, there is a growing interest in analyzing and modeling EVP behavior of synovial fluid and human blood.^3,4 They exhibit a complex mechanical response characterized by the interplay of elastic, plastic, and viscous behaviors.^3,5,6 The yield stress is a key property of EVP materials, marking the point at which they transition from solid-like deformation to fluid-like flow.⁷ For EVP materials in this paper, we specifically use the term ‘yield’ to refer the elastic-to-plastic transition observed in oscillatory strain measurements. This definition may seem different from the classic concept of Bingham yield, which implies a rigid body below a critical stress threshold.⁸

Fig. 1 shows a typical amplitude sweep curve of type III nonlinear behavior to illustrate the elastic, plastic, and viscous-dominant regions depending on the applied strain amplitude. The elastic-dominant region is located at low strain amplitudes where storage modulus (G′) is significantly larger than loss modulus (G″), and both are relatively strain-independent.⁹ In this region, material behaves like a solid, storing energy elastically, implying that the microstructure remains mainly intact.^10–12 Beyond the linear, elastic-dominant regime, G′ declines while G″ increases, as the material begins to yield due to structural reformation.^13,14 This intermediate regime is known as the plastic or transition region, and often extends to the crossover point where G′ = G″. Above this point, where G″ overtakes G′, the post-yield liquid-like behavior reflects dominance of viscous dissipation as the irreversible flow results from structural breakdown leading to bulk fluidization.^11,15 Concentrated nanoemulsions, as one of the examples of EVP materials with type III nonlinear behavior,¹⁶ have rich structures, properties and applications,^8,17,18 highlighting the importance of studying their rheology. Generally, dilute emulsions show Newtonian behavior, while emulsions with sufficiently high volume fraction of dispersed phase become viscoelastic and EVP.¹⁹

Fig. 1 Dynamic moduli of type III nonlinear behavior for a typical elastoviscoplastic material from strain amplitude sweep test in oscillatory mode.

Several studies have tried to develop a proper model for mimicking the viscoelasticity,^20–22 such as the Oldroyd model (Fig. 2a),²³ and the Jeffrey's model (Fig. 2b).²⁴ These models are not able to account for the nonlinearity due to yielding in EVP materials. Many of these efforts involve modified versions of the Oldroyd model to account for the yield stress.^20,25,26

Fig. 2 Mechanical analogs for different models, with components arrangement for: (a) Oldroyd model, (b) Jeffrey's model, (c) Saramito model, and (d) strain-separated Saramito model.

The complexity of numerical fitting is the primary challenge in the fitting of yielding in conditional discrete EVP models. Model fitting typically relies on gradient-based optimization methods, which begin with an initial set of parameter guesses and iteratively update them by evaluating whether the model satisfies the governing equations under the current parameter set. This approach also requires finite difference calculations of gradients of the data to minimize the residual of differential equations, which can become misleading when dealing with noisy data. Even if we circumvent the disturbed gradients of noisy data by directly solving the equation, a numerical solver must go through all the data points, which means the solver's iteration loop should be inside the optimization iteration loop. This can make the convergence extremely time-consuming, while the discontinuity in the discrete model may also result in diverging during solving the equation in each loop.

To address all the associated challenges with numerical fitting, we adopt a physics-informed neural network (PINN) approach for fitting constitutive equations.^27,28 In recent years, PINN has emerged as an efficient mathematical tool for solving and fitting differential equations^27,28 in many engineering and science fields, including material science,^29,30 chemistry,^31–33 fluid mechanics,^34–36 and specifically rheology.^37–40

Unlike classic neural network approaches, which act like a black-box predictor with no physical intuition, PINN embeds a constitutive model to set a meaningful theoretical base.²⁷ PINN fitting algorithm uses gradient descent logic or more sophisticated variants of the gradient descent, like the Adam optimizer.⁴¹ The only difference in comparison to conventional numerical approaches is that neural networks' output can serve as an intermediate between the experimental data and the model, to alleviate the challenges related to a lack of data points or noise. The detailed instructions of this method for rheology models with sets of examples are provided by Mahmoudabadbozchelou et al.³⁹

PINN optimizes the neural network to align its output predictions with experimental data by minimizing a corresponding loss function known as data loss. Simultaneously, it optimizes the coefficients of a user-suggested governing equation in a way that its solution matches the expected values by minimizing a corresponding loss function known as physics loss. The simultaneous optimization of these two loss functions is based on the minimization of the weighted summation of them. Unlike traditional numerical methods, which can become unstable due to differentiation of noisy data, PINNs can still converge to physically consistent solutions as regularization can prevent overfitting to the noise of data.⁴² Moreover, the gradient calculation for a neural network is more efficient. This is due to the fact that the neural network itself is a continuous, differentiable function, even though we only evaluate it at discrete points, and the automatic differentiation can provide all the gradient values simultaneously. In contrast, conventional finite difference approximations require step-by-step calculations, which are much slower.

Aside from the advantages of PINN, some limitations of the conventional gradient descent approach persist. First of all, an improper learning rate is prone to time-consuming convergence and skipping the global minimum of loss function, or getting trapped in local minima instead of locating the global minima.⁴³ Another significant limitation arises from the quality of experimental data. Sparsity and noise in data can still affect the accuracy of PINN.⁴⁴ Thus, PINNs require careful weighting between data-driven loss and physics-based residual loss. This balancing prevents overfitting of the neural network to noise and forces the output to follow the physics of the model, which introduces additional hyperparameters (α). Depending on the complexity of the model, feeding one dataset from one rheological measurement protocol may not be enough for the PINN to recover correct parameters for a complex model.⁴⁵

Few studies used PINN or any other fitting method to directly fit rheological models to LAOS data.^39,46 To the best of our knowledge, none of them ended up with a model that accurately covers various strain amplitudes with a single set of fitting parameters. In the present work, therefore, we address two challenges. First, we suggest a model that can capture physical behavior at various strain amplitudes from linear to nonlinear regions. Second, we develop a new 2-step PINN framework to fit the suggested model simultaneously to several stress datasets of different strain amplitudes. Instead of directly penalizing residuals by gradient calculation, we enforce the governing equation through a cumulative reconstruction of stress with integration, which eliminates the need for explicit initial-condition constraints commonly used in PINN formulations, leading to faster convergence. In the following sections, we present a brief review of the EVP models and the associated challenges with conceptual interpretations and numerical complexities in fitting them to measured experimental data. Then, we address these challenges by suggesting an updated version of the model and the corresponding fitting approach to overcome the limitations.

Modeling EVP behavior

One of the most popular EVP models was proposed by Pierre Saramito,²⁰ which is simple enough even for numerical simulations.^47,48 In this model, a controlling slit (also known as a solid friction element) was introduced to represent yielding behavior (Fig. 2c). We use τ notation to present the stress attributed to the lower section of the model marked with a blue dashed box in Fig. 2c. Considering the stress contribution of the recoverable dashpot (η) added to τ, the total stress (σ) becomes:

σ = η[small gamma, Greek, dot above] + τ (1)

In the section of the model attributed to τ, the spring (G) is in series with the unrecoverable dashpot (illustrated in red) with η_p notation. This unrecoverable dashpot itself is in parallel with a slit element with activation stress of τ₀. When τ < τ₀, the slit remains stationary and do not let the parallel η_p to engage. Therefore, the model effectively simplifies to a linear Kelvin–Voigt model based on active G and η. However, when the applied stress in this section exceeds τ₀, the slit begins to move, allowing the red dashpot in Fig. 2c to activate, which introduces unrecoverable plastic deformation and nonlinearity to the model. It should be noted that this recovery refers to the recovery of strain. The equation for τ is a conditioned Maxwell model as follows:where only represents the relaxation of the blue box section in Fig. 2c.

Considering yield as a point where unrecoverable strain becomes effective, we can relate the yield stress of Saramito model (σ₀) to τ₀ as follows:

σ₀ = η[small gamma, Greek, dot above]_y + τ₀ (3)

where [small gamma, Greek, dot above]_y is the strain rate at which the slit is activated, which is dependent on frequency and strain amplitude. Therefore, depending on the [small gamma, Greek, dot above]_y, σ₀ changes. It should be noted that σ₀ is referring to the onset of plastic flow, not the flow from the rest state.

Donley et al.²¹ revisited Saramito model into a strain-separable version as in Fig. 2d, where the unrecoverable yielding-controlled dashpot is placed in series to a Kelvin section. Kamani et al.⁴⁹ replaced the linear unrecoverable dashpot in the same arrangement with a Herschel–Bulkley version and removed the slit. This model is known as KDR (Kamani–Donley–Rogers). Unlike the classical Saramito model, where yielding is governed by a fraction of total stress, the strain-separable versions define yielding based on the total stress, so the yield is constant in this model (in contrast to eqn (3)).

While the Saramito model is commonly used to represent EVP material, to the best of our knowledge, no study has successfully fitted this model to rheological data. For example, Cheddadi et al.⁵⁰ used the Saramito model to simulate foam flow by estimating the model's parameters according to prior experiments, rather than fitting data. This approach does not use the model to predict the flow behavior but rather assesses the model's capability to explain EVP behavior. In another study, Abdelgawad et al.⁵¹ used the Saramito model with a Herschel–Bulkley unrecoverable dashpot to simulate the flow of a Pluronic F127 solution through a single periodic wavy channel, and estimated constant parameters of the model from experiments. They observed only a qualitative match of flow fields and velocity profiles that does not fully align with the experimental data. Fraggedakis et al.²⁶ tried to fit a Carbopol data to the Saramito using a nonlinear regression approach by minimizing the error between experimental and simulated dynamic moduli, rather than directly fitting the model to stress data. However, this approach failed, and to address this issue, they replaced the linear slit-controlled dashpot with an isotropic-kinematic hardening model (IKH model).⁵² While they achieved a satisfactory fitting, this approach only considers the first harmonic of stress when fitting the model to G′ and G″ in the nonlinear region of the data. Even for the strain-separable Saramito and KDR, model parameters to predict LAOS curves were estimated from linear viscoelastic moduli and steady-shear flow curves,^21,53 which do not appear to accurately fit nonlinear stresses.

As we show in the results and discussion section, there is no unique single set of parameters for the model to represent the real type III materials at all strain amplitudes. This implies that the assumption of constant model parameters across different strain amplitudes is invalid, as the underlying microstructure formation and deformation may be strain dependent. There are some experimental findings discussing the structural changes under high strains/stresses. For example, the Payne effect^21,54 leads to a decrease in the elasticity and enhanced plastic flow with strain amplitude as the microstructure undergoes progressive breakdown. In the classic view, the Saramito model assumes a constant modulus for all elements to represent energy storage and loss, which is not consistent with microstructural changes. In other words, changes in the viscoelasticity of EVP materials cannot be solely attributed to variations in the contributions of linear springs and dashpots. The need for parameter adjustment reflects the model's inability to inherently capture the continuous evolution of viscoelastic and plastic contributions across different deformation amplitudes.

To extend the Saramito model to a more comprehensive version, we suggest a new approach in updating the model. First, we propose a physically meaningful correspondence between each component of the Saramito model and the underlying structural phenomena that govern rheology. Table 1 lists different phenomena for storing or dissipating energy in each element of the Saramito model. When the material structure is distorted reversibly in emulsions¹⁹ both energy and strain remain reversible upon removal of the stress load. This reflects elastic behavior, which can be represented by the spring element. In contrast, some processes involve energy dissipation without contributing to permanent deformation, like the rolling of particles over each other in colloidal materials,¹⁹ or the local vacancy diffusion of particles into neighboring voids.⁵⁵ These phenomena are captured with the recoverable-strain dashpot (i.e., the dashpot in parallel to the spring). The rest of the irreversible structural breakdowns lead to simultaneous energy dissipation and unrecoverable deformation that are represented by the unrecoverable dashpot, which is activated above the yielding point.

Table 1 Examples of deformation modes corresponding to the mechanical EVP model's elements

Element	Examples
Spring (G)	• Reversible colloidal deformation^56
	• Reversible structure distortion^19
	• Reversible stretching and recoiling of polymer chains^57
	• Swelling tendency of hydrogels due to osmotic pressure^58
Recoverable dashpot (η)	• Diffusion to void spaces^55
	• Rolling of particles/colloids^59
	• Polymer reptation^60
Unrecoverable dashpot (ηp)	• Irreversible structure breakdown^61,62
	• Mutual displacement of large blocks^63
	• Sliding^64
	• T1 neighboring switch in foams^65,66

---

Agrawal et al.⁶⁷ analyzed the effect of two different values of η = 500 Pa s extracted from the linear viscoelastic region and hypothetical η = 1 Pa s in the Saramito and KDR model, while keeping other parameters of the model fixed to study the accuracy of the model in predicting LAOS stress of Pluronic F127 hydrogel. The result showed that η = 500 Pa s is accurate in the linear region but overpredicts stress at large amplitudes, while with η = 1 Pa s it underpredicts the stress at low amplitudes. This can be interpreted in a way that dissipation mechanisms, such as droplet rolling, are more significant at low strain amplitudes or linear regimes. Such observations indicate that assuming a fixed value of η is a major disadvantage of these models, preventing them from accurately representing the reality of type III materials.

In this work, we suggest a nonlinear strain rate thinning recoverable dashpot (η([small gamma, Greek, dot above])). This modified model will remain thermodynamically consistent with respect to the second law of thermodynamics, regardless of how η([small gamma, Greek, dot above]) function is defined. Based on applied power (P) and stored energy (S) the model must satisfy:

D = P − Ṡ≥ 0 (4)

where γ_e is elastic strain, referring to the share of total strain that is associated with the spring. Eqn (2) was used to express the rate of stored energy as follows:

P = σ[small gamma, Greek, dot above] = τ[small gamma, Greek, dot above] + [small gamma, Greek, dot above]²η([small gamma, Greek, dot above]) (7)

By plugging eqn (6) and (7) into eqn (4), we have:Thus, for all positive values of η([small gamma, Greek, dot above]), the model remains thermodynamically consistent. Several candidate functions can offer shear rate thinning behavior, such as empirical power decay functions,⁶⁸ like Carreau fluid model, where the viscosity decreases smoothly from a low-rate plateau to a high-rate plateau.⁶⁹ In the power functions, the power index can adjust the curvature for accurate fitting. However, the trade-off is that this additional index makes the model more complicated. In contrast, an exponential decay form can be viewed as a simplified monotonic decay with fixed curvature, chosen as a lower-parameter alternative. Moreover, in natural processes, Euler's number e ≈ 2.71828 appears wherever growth or decay is continuous and proportional to the current state.⁷⁰ This proportionality in natural phenomena makes their mathematical descriptions commonly exponential, e.g., in biological growth,^71–73 physical decays,^74–76 and chemical reactions.^77,78 Accordingly, the exponential trend commonly exists in rheology, where the exponential stress relaxation is a well-established concept.^79,80

For modified Saramito model, we redefine η of the recoverable dashpot in Fig. 2c (black colored dashpot) as eqn (9) to capture the underlying continuous structure breakdown, which causes unrecoverable dashpot effectiveness to diminish.

In this equation, η₀, η_∞, and [small gamma, Greek, dot above]* are the zero-shear viscosity, infinite-shear viscosity, and critical shear rate respectively. The general form of the updated model is still similar to the classic Saramito model and can be fitted by PINN in a similar way. The only difference is that instead of fixed values of η, there are different corresponding values at each specific strain rate. For further simplicity and parameter reduction, we fix η_∞ = 0 since we expect the failed structure under strain to dominantly exhibit unrecoverable energy dissipation, and the recoverable dashpot contribution disappears. While the classic Saramito switches between Kelvin and Oldroyd below and above yield, imposing a decaying recoverable dashpot makes the model softly switch between Kelvin and Maxwell, with a gradual Oldroyd transition in between. However, there is still the sharp activation of constant η_p, decided by the constant τ₀. Even though in reality there is no sharp yield transition, the simplistic view is that at some point the structure fails irrecoverably, which implies yielding.^19,81 In the Results and discussion section, we show that this modification to the model enables it to represent oscillatory behavior from linear to nonlinear regions of strain amplitude.

Methodology

To check the developed model, we prepared oil-in-water nanoemulsions as an example EVP material. In the first step, we mixed 25% volume fraction of polydimethylsiloxane (PDMS from Clearco Products) in sodium dodecyl sulfate (SDS from Sigma-Aldrich) aqueous solution by magnetic stirring for 5 min at 700 rpm to form macroemulsions. The next step was homogenization to convert macroemulsions to nanoemulsions by ultrasonication (VCX500 ultrasonic processor, Sonics & Materials Inc.) at 40% amplitude until achieving an average droplet size of 125 nm. The process was carried out in an ice bath to maintain the temperature at approximately 15 °C and preserve the stability of the emulsion. Finally, we let the water in the aqueous phase evaporate at room temperature and constant humidity to concentrate the nanoemulsions to 40% volume fraction. The total SDS surfactant concentration in the system was adjusted to reach 180 mM after evaporation. After preparing the samples, oscillatory shears with amplitudes ranging from 1% to 100% applied using a parallel plate rheometer (Hybrid HR-3, TA Instruments) at 10 rad s⁻¹ frequency to capture both linear and nonlinear Lissajous curves. Since HR-3 rheometer is stress-controlled, it produces noisy data. We also tested other typical EVP samples such as kraft real mayonnaise sauce and pepsodent complete care toothpaste. These samples were measured with a strain-controlled rheometer (ARES-G2, TA Instruments), which resulted in cleaner, less noisy data.

We have implemented two customized PINN algorithms for fitting the Saramito model, which is detailed in the Appendix. The Python script is also available on GitHub, based on the PyTorch library. The first algorithm is the simple PINN fitting of the classical Saramito model to stress data at a single strain amplitude, shown in Fig. 3. The second one, presented in Fig. 4, is a two-step approach for the simultaneous fitting of the our modified Saramito model to stress data at several strain amplitudes. In the Results and discussion, we describe why this approach is essential to get meaningful fitting.

Fig. 3 Schematic workflow of PINN fitting for the Saramito model.

Fig. 4 Schematic workflow 2-step PINN fitting for the modified Saramito model.

Results and discussion

Modeling EVP materials with constant properties during LAOS

To validate the PINN fitting before fitting experimental data, we first attempted to fit synthetic data of the classical Saramito model as in eqn (1) and (2). Fig. 5a and b compares the fitted model with synthetic data generated from the parameter sets of , τ₀ = 75 Pa, η_p = 20 Pa s, and η = 10 Pa s for the strain amplitude of 0.5 and frequency of 10 rad s⁻¹. The fitting to stress vs. time is acceptable, and all fitting variables converge towards their corresponding expected original reference values. Additional plots are presented in the supporting information (SI). Fig. S1 shows the optimization of total loss with respect to epoch progress. Fig. S2 shows the optimization of fitting parameters of the Saramito model with respect to epoch progress.

Fig. 5 Difference between synthetic stress from the Saramito model and PINN fitting (after 90 000 epochs), shown as stress vs. time (a) and (c) and stress vs. strain (b) and (d). Parameters: (a) and (b) λ = 0.05 s, τ₀ = 75 Pa, η_p = 20 Pa s, η = 10 Pa s; (c) and (d) λ = 0.1 s, τ₀ = 25 Pa, η_p = 10 Pa s, η = 20 Pa s.

To investigate the influence of the weighting factor of α on the total loss function (data loss + α × physics loss) on PINN performance, we compared their fitting results after 10 000 epochs using four values of α in the range of 0.6 to 1.5 in Table 2. Despite variations in total loss, the recovered values of λ, η, η_p, and τ₀ remained relatively stable across different α values, showing only slight deviations. This suggests that the general framework of the PINN is sufficiently robust and stable for parameter estimation. Notably, α = 1.5 slightly compromises parameter accuracy (e.g., η_p = 16.10 vs. the original model value of 20), which can be attributed to the lower contribution of data loss in optimization, leading to a slightly off-fitted neural network. On the other hand, α = 0.6 produced the lowest loss with estimated parameters very close to the expected values, suggesting it may provide the most reliable trade-off between data fitting and physical model constraint enforcement. Fig. S3 compares the loss and elastic moduli of the original model and PINN fitted model in amplitude sweep results, confirming stable convergence of PINN regardless of various α values.

Table 2 The effect of the weighting factor α on PINN performance after 10 000 epochs

Predefined coefficients		α = 1.5	α = 1.0	α = 0.8	α = 0.6
		Loss = 4.83 × 10^−4	Loss = 1.486 × 10^−4	Loss = 6.036 × 10^−5	Loss = 5.177 × 10^−5
λ (s)	0.05	0.0405	0.0453	0.0470	0.0462
η (Pa s)	10	10.25	10.80	10.60	10.61
ηp (Pa s)	20	16.10	17.61	18.37	18.04
τ0 (Pa)	75	84.71	81.18	79.48	80.18

---

Aside from the advantages of PINN, some limitations of the conventional gradient descent approach persist. First, a small learning rate leads to slow convergence, and a large learning rate can skip the minimum value. The loss landscape may also have multiple local minima, making it challenging for gradient descent algorithms to locate the global minimum consistently.⁴³ Additionally, several possible fits of a complex model for a single set of experimental data may exist. Thus, the fitting algorithm may not accurately recover the original parameters.⁴⁵ As an example, for the synthetic data from the parameter sets of , τ₀ = 25 Pa, η_p = 10 Pa s, and η = 20 Pa s, the single Lissajous curve at the strain amplitude of 0.5 is not informative enough for recovering the fitting parameters because the fitting ends up with a simpler model expression where τ₀ converges to zero (see Table 3 and Fig. S4), even though the Fig. 5c and d suggests satisfactory stress prediction. In other words, the PINN fitting parameters with τ₀ = 0 Pa result in a stress function, which is numerically indistinguishable from the original synthetic data.

Table 3 The failure of parameter recovery of non-informative synthetic data for various epochs at α = 06

Predefined coefficients		Epoch = 10 000 (8 seconds)	Epoch = 50 000 (40 seconds)	Epoch = 90 000 (72 seconds)
		Loss = 2.531 × 10^−4	Loss = 1.441 × 10^−4	Loss = 1.428 × 10^−4
λ (s)	0.1	0.078	0.088	0.079
η (Pa s)	20	13.71	12.94	11.82
ηp (Pa s)	10	14.48	17.86	18.20
τ0 (Pa)	25	14.48	0.00	0.00

---

An important observation is when τ₀ is comparable to the induced stress in the spring, a more pronounced yield transition is observed during the oscillatory cycle at large strain amplitudes. The yielding produces stronger distortions in the stress–strain Lissajous curves. Physically, a higher τ₀ delays the onset of plastic flow, resulting in more abrupt yielding once the stress threshold is exceeded. In contrast, decreasing τ₀ reduces the nonlinearity and numerical resolution for fitting the parameters of the model.

This fitting failure is not an issue only associated with the PINN, but a general limitation of any fitting approach when the data is not informative for the corresponding model. As a solution, we suggest a two-step fitting to two separate linear and nonlinear oscillatory stress datasets of a sample. In this approach, G and η can be obtained by fitting the linear data with the Kelvin–Voigt model, which is equivalent to the unyielded version of the Saramito model. Then, by fixing these obtained parameters in the Saramito model, we proceed to PINN fitting of the nonlinear yielding oscillatory stress data for obtaining τ₀ and η_p. Table 4 compares the original and fitting parameters by this two-step approach for two examples with different values of λ. By capturing recoverable elements of the model separately from the linear rheological data, the approach is less prone to misleading fitting of nonlinear behavior and getting trapped in a local minimum. Furthermore, it exhibits faster performance in view of less calculation at each epoch and also the fewer required epochs for PINN to fit only τ₀ and η_p parameters.

Table 4 The improved parameter recovery of synthetic data after 30 000 epochs with a two-step approach for two sets of parameters

Predefined coefficients		α = 0.6	Predefined coefficients		α = 0.6
		Loss = 5.13 × 10^−6			Loss = 7.48 × 10^−5
λ (s)	0.1	0.114	λ (s)	0.05	0.051
η (Pa s)	20	20.00	η (Pa s)	20	20.00
ηp (Pa s)	10	11.42	ηp (Pa s)	10	10.10
τ0 (Pa)	25	24.46	τ0 (Pa)	25	25.11

---

Modeling EVP materials undergoing structural changes in LAOS

Thus far, we have dealt with clean synthetic data, where the precise values for the parameters were known to verify the credibility of the method. The next step is to assess Saramito model's credibility for real type III materials. In this study, we use the nanoemulsion sample. Fig. 6a displays the fitted Saramito model to the experimental stress at a frequency of 10 rad s⁻¹ and various strain amplitudes. For this real noisy data, the balance between physics loss and data loss prevented the neural network from overfitting the considerable noise, providing a smooth Lissajous curve suitable for parametric and comparative studies of different samples. Unlike synthetic data, experimental data lacks direct reference values for coefficients to verify the validity of the fitted model at a single strain amplitude. Even if we neglect this challenge and assume perfect fitting can be achieved at each strain amplitude, different sets of fitting coefficients are obtained. Therefore, single-step fitting (Fig. 6a) at best only represents the material's behavior within that particular deformation regime, and the algorithm can converge to various combinations of fitting parameters at each run.

Fig. 6 Comparison between model and experimental data of a typical nanoemulsion based on: (a) separate fitting at each strain amplitude, and (b) fitting the Saramito model at the strain amplitude of 0.5 and frequency of 10 rad s⁻¹.

To demonstrate the advantage of PINN fitting, it is compared with conventional fitting methods, such as gradient descent, in Table S2. The conventional fitting methods led to divergence due to severe noise in the data because it introduces drastically misleading values either by differentiation or integration during residual checking. It should be noted that by adjusting a good initial guess, Levenberg–Marquardt can converge for some noisy examples, but remains unstable, ending up with unreasonable infinite values. The comparison of fitting approaches for nanoemulsion sample at 10 rad s⁻¹ and the strain amplitude of 0.5 is provided in Table S1.

To assess the Saramito model's predictivity for real type III materials, we fit the data at the strain amplitude of 0.5 and then compared the model's prediction with experimental data at lower and higher strain amplitudes, as shown in Fig. 6b. The model significantly underpredicts the stress magnitude below the strain amplitude of 0.5 and overpredicts the stress magnitude above 0.5 strain. However, a meaningful fitted model should be able to predict unseen experimental data, with the same set of parameter constants. Otherwise, it is not useful for simulations and practical applications.

To come up with a meaningful fitted model for all strain amplitudes, the two-step fitting can be an appropriate approach, as worked for synthetic data. One other possible solution might be to fit the model to all strain amplitudes at the same time with a cumulative loss function. However, this approach does not lead to convergence during the fitting process, indicating that the Saramito model is inherently unable to accurately capture the flow behavior across the entire range of strain amplitudes.

In the next steps, we fit the modified Saramito model which assumes a decaying recoverable dashpot as in eqn (9). The increased number of parameters complicates the fitting due to the exponential decay function. Therefore, we use an enhanced version of the 2-step fitting framework as shown in Fig. 4. First, G and η₀ can be fitted to the model for very small strain amplitudes in the linear region, e.g. 0.01. As discussed in the previous section, this approach ensures fitting the unyielded regime with the unyielded model. The remaining parameters can be fitted using separate PINNs, but with an aggregated loss function for several strain amplitudes in the nonlinear regions. Here, we consider an equally weighted culminative loss function, so optimization sensitivity to each strain amplitude is equal for all of them. Obviously, based on this cumulative fitting approach and more complex model, a meaningful convergence requires more epochs and results in larger values of the loss function.

Fig. 7a shows the fitting of Saramito with the updated decaying recoverable dashpot to the samples of the nanoemulsion with 50 000 epochs (120 seconds). In addition to the enhanced fitting of the modified Saramito model, the 2-step approach offers a fair fit for the wide range of linear to non-linear regimes with a single set of parameters. However, there is a slight off prediction at the strain amplitude of 0.1, where the transition from elastic to plastic regime happens. This is due to the assumption of a constant unrecoverable dashpot and the sharp yielding transition, which is a trade-off for avoiding over-parametrization.

Fig. 7 Fitting of (a) modified Saramito model with decaying recoverable dashpot; (b) classical Saramito model, and (c) KDR model to nanoemulsion data.

In Fig. 7b and c, we also included fitting to the classical Saramito and KDR model as two important EVP models. We explained the equations of the KDR model in supporting information (SI). By adopting a similar 2-step fitting approach for them to fix the fitted parameter from linear data and fitting the rest of the parameters to nonlinear datasets, both classical Saramito and KDR models fail to converge meaningfully. The issue with these models is that enhancing the accuracy of fitting in the linear region decreases accuracy in the nonlinear region fitting, and vice versa.

To assess the uncertainty of the modified Saramito, a bootstrap analysis was employed. In this approach, we redistribute the error of each datapoint to generate new noisy data and check how the fitted parameter changes. Performing five bootstraps, the stage-1 parameters, G and η₀, remained unchanged across bootstrap replicates. Among the stage-2 parameters, the yield stress τ₀, plastic viscosity η_p, and rate-softening parameter [small gamma, Greek, dot above]* showed relatively small variability across the bootstrap runs. According to Table 5, τ₀ remained within a narrow range of approximately 204–218 Pa, while η_p varied between about 22.4 and 25.7 Pa s. These results suggest that the nonlinear constitutive parameters are reasonably robust to the noise level present in the experimental data, although this should be interpreted as a local uncertainty estimate rather than a full global identifiability analysis. It should be noted that this uncertainty is dependent on the state of the model. For example, if one dashpot dominates the stress signal, it will increase the uncertainty for recovering other parameters.

Table 5 The improved parameter recovery of synthetic data after 30 000 epochs with a two-step approach for two sets of parameters

Parameter	Mean	Standard deviation	95% Lower confidence interval	95% Upper confidence interval
G (Pa)	9664.82	0.00	9664.82	9664.82
η0 (Pa s)	160.38	0.00	160.38	160.38
η∞ (Pa s)	0.00	0.00	0.00	0.00
[small gamma, Greek, dot above]* (s^−1)	2.28	0.12	2.16	2.41
ηp (Pa s)	24.01	1.32	22.41	25.75
τ0 (Pa)	211.20	5.93	203.92	217.71

---

In Fig. 8, we showed two more examples of fitting the modified Saramito to the kraft real mayonnaise sauce and the pepsodent complete care toothpaste. The strain amplitudes of the input data for the two-step fitting must capture the linear to completely yielded regimes. Thus, the selection of the strain amplitudes depends on the material. For example, for the mayonnaise, we considered a strain amplitude of 2 to cover the fully yielded regime. The pronounced yield in data is the informative pattern for fitting the model. The model shows a fair fit for the mayonnaise as a type III soft matter. We also included the toothpaste example as a type I material, for which the model still provided an acceptable approximation, but with lower accuracy at higher amplitudes due to smaller yields.

Fig. 8 Fitting of the modified Saramito to (a) kraft real mayonnaise (50 000 epochs), and (b) pepsodent complete care toothpaste (80 000 epochs).

Having a single set of parameters for modeling EVP materials in all strains provides the opportunity for systematic studies and simplicity for interpretation. The current model is inherently frequency-dependent, due to the presence of combined spring and dashpots. However, it does not consider the direct functionality of dashpots to frequency. Accordingly, for further improvement, we need to embed the frequency dependence in the elements. In general, soft materials exhibit weak frequency dependence.^16,82 One possible approach for our future work is adopting the strain rate frequency superposition introduced by Wyss et al.⁸² By using this approach, we can approximate the response at other frequencies by altering strain amplitude.

Conclusion

In this work, we revisited the Saramito constitutive equation for modeling EVP materials to overcome key limitations of conventional EVP models and the identifiability challenges when it comes to fitting these models. Moreover, we discussed that a single set of fitting parameters cannot describe flow behavior across a wide range of strain amplitudes in both linear and nonlinear regions. Thus, we proposed a modified version of the Saramito model. To address fitting complexities of the updated model, including sensitivity to noise, reliance on numerically unstable gradient estimates, and difficulties imposed by discontinuous yield point, we have introduced a two-step PINN fitting approach. We showed that the modified Saramito model fitted by cumulative-loss PINN captures both linear and nonlinear regimes with a single parameter set. Although slight deviations persist near the elastic–plastic transition, the method substantially improves predictive capability while maintaining physical interpretability. This work highlights that incorporating strain-dependent recoverability is essential for modeling type III EVP materials and suggests future extensions incorporating frequency-dependence to achieve predictive rheology across both strain and frequency domains.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code for Physics-Informed Neural Network (PINN) for fitting modified Saramito model can be found at GitHub (https://github.com/babkvg/EVP-PINN). The version of the code employed for this study is version 1.0.

Appendix: workflow of the PINN fitting

One consideration with fitting models to LAOS experimental data is that t = 0 is assumed at the beginning of the recorded steady state cycle and does not take the transient start from zero stress into consideration, which means that some informative parts of the data are missing. However, rheological models will essentially lead to zero stress value when strain and strain rates are plugged in at t = 0 s as the real starting point. Additionally, it is common practice in PINN fitting to have the network learn the time derivative of the stress directly via the autograd graph. However, when the ODE constitutive model is conditional (for example, yielded and unyielded conditions), physics involves history-dependent accumulation, as seen in a yield-stress fluid, where elastoviscoplastic stress builds up while switching between unyielded to yielded regimes. Therefore, history dependency must be managed in the PINN approach. Otherwise, blindly embedding an ODE solver into the training graph will not converge to true cumulative stress values. To address these challenges, an updated step-by-step workflow of the PINN fitting is presented here:

Data acquisition and normalization

First, the algorithm will normalize all stress, strain, and strain rate data points between −1 to 1 to facilitate convergence of fitting.

Embedding neural networks into model calculations

Fitting is a reverse approach to check if the final response satisfies the equation with optimized sets of parameters. In this work, we implemented a fully connected feedforward architecture with an input layer of size 1, two hidden layers with 50 neurons each, and an output layer of size 1. Rectified linear unit (ReLU) activation functions are used in the hidden layers. The input is time (t), and the neural network predicts total stress (σ_NN) as the output. In addition to the network weights and biases, the model parameters λ, η_p, η, and τ₀ are treated as trainable variables within the PINN framework and are constrained to be non-negative through activation functions. First, there is no need to implicitly encode history dependence or conditional logic. It needs only to learn a scalar mapping from time to stress. Then, instead of directly plugging noisy σ into eqn (1), the algorithm plugs in σ_NN as total stress in the rewritten form of eqn (1), as in eqn (1A), to obtain τ, which we note as τ_M-NN since it is indirectly derived from the neural network.

τ_M-NN = σ_NN − η[small gamma, Greek, dot above] (1A)

At early epochs, the neural network output provides off predictions from experimental total stress, but it eventually converges to an accurate presentation of the total stress for efficient calculations. As mentioned before, in transient LAOS measurements performed by a rheometer, the reference time of recording stress is taken as a time where the system achieves steady-state Lissajous curves. To adjust the initial value at the starting point and embed the inherited stress when the model switches between Kelvin and Oldroyd, we need to calculate the rate of τ_M-NN by rewriting eqn (2) in the form of eqn (2A). To avoid dividing by zero, the ε = 10⁻⁸ or any other small value is considered in the denominator as follows: Here the term mimics yield-like behavior if |τ_M-NN| is bigger than τ_M0. Otherwise, this term becomes zero. Having At each datapoint, a cumulative summation results in τ_model values while complying with the initial value. Then, we can back-calculate total stress (σ_model), which is still entirely differentiable for straightforward gradient-based optimization after passing through elementary tensor operations. Thus, the model is not directly solved, but stress is being checked to see if it is satisfying the model regarding the initial value. This approach is faster than calculating from autograd command and checking the residual of the equation. Moreover, we do not need additional loss term to enforce the initial condition, but inherently it is embedded in the integration.

Defining loss function

For convergence with gradient descent fitting, a credible loss function must be defined as a metric for the accuracy of prediction. Thus, the difference between σ_NN and model-predicted total stress (σ_model) with experimental σ must be checked during optimization epochs. Furthermore, the algorithm should check if it can back-calculate τ_model from eqn (1A). This will introduce additional terms to the loss function. In this regard, we define data loss (Loss_data) and physics loss (Loss_phy) based on mean square error (MSE) to measure how well the PDE-type relationships are satisfied:

Loss_data = MSE(σ_NN, σ) (3A)

Loss_phy = MSE(τ_NN, τ_model) + MSE(σ_model, σ_NN) (4A)

In case of clean and noiseless experimental data, MSE(σ_model, σ_NN) leads to faster convergence as it provides an accurate reference of convergence from the first epoch. Finally, the total loss will be a weighted summation of data and physics losses. By considering the weight factor α for the physics loss, we can control how strongly the PINN enforces governing equations versus fitting observations during minimization of total loss, ultimately improving training stability.⁸³

Total loss = Loss_data + αLoss_phy = Loss_data + α(Loss_phy1 + Loss_phy2) (5A)

Considering , the PINN algorithm will update the τ_M0, λ, η, η_p, and σ_NN at each epoch with the Adam optimizer to minimize the Total Loss. In the case of simultaneous fitting, we add the total loss of all separate datasets together for the optimization.

Unscaling fitted parameters to the original scale

Finally, after reaching a satisfying threshold for the loss function (in this work, it was adjusted at 10⁻⁶) or completing the adjusted epochs, the fitted parameters must be denormalized based on their units and half-range values of stress and strain rate as follows:Early in training, the learning rate ramps up linearly to avoid large gradient steps too early. After warmup epochs, a cosine function is used to reduce the learning rate gradually. Thus, the learning rate (LR) equations are as follows before and after 500 epochs, considered as the threshold for warm-up with linear learning rate changes.This cyclical schedule can encourage better convergence and avoid getting caught up in a local minimum. In this approach, the network remains lightweight (predicting a single scalar per timestep), and the physics solver uses efficient tensor operations rather than iterative ODE solves in the autograd graph. Also, we avoid long chains of back-propagated operations that can vanish or explode over many timesteps. Moreover, the piecewise yielding law is implemented exactly outside the network, ensuring no approximation error in regime switching or integration.

The Python code is also available on GitHub. The code was run on the G4 GPU of the Google Colab service. In terms of training duration, we observed that every 10 000 epochs requires approximately 8 seconds for the single fitting and 24 seconds for the two-step fitting.

Acknowledgements

The authors would like to thank Dr. Zahra Abbasian Chaleshtari for providing the nanoemulsion rheological data. We also thank Dr. Iman Ghamarian (University of Oklahoma) for his insightful discussions on machine-learning approaches.

References

1. H. L. França, M. Jalaal and C. M. Oishi, Phys. Rev. Res., 2024, 6, 013226.
2. S. T. Chan, S. Varchanis, A. Q. Shen and S. J. Haward, PNAS Nexus, 2023, 2, pgad042.
3. M. Armstrong, A. Pincot, S. Rogers, T. Knight and D. Bailey, Front. Phys., 2022, 10, 889065.
4. M. Armstrong, K. Rook, W. Pulles, M. Deegan and T. Corrigan, Rheol. Acta, 2021, 60, 119–140.
5. R. R. Fernandes, D. E. V. Andrade, A. T. Franco and C. O. R. Negrão, J. Rheol., 2017, 61, 893–903.
6. I. Valizadeh, T. Tayyarian and O. Weeger, Addit. Manuf., 2023, 72, 103641.
7. F. Accetta and D. C. Venerus, Rheol. Acta, 2024, 63, 719–730.
8. N. Sanatkaran, M. Zhou and R. Foudazi, J. Rheol., 2021, 65, 453–461.
9. R. Agrawal and E. García-Tuñón, Soft Matter, 2024, 20, 7429–7447.
10. K. Suman, S. Shanbhag and Y. M. Joshi, J. Chem. Phys., 2023, 158, 054907.
11. V. V. Erramreddy, S. Tu and S. Ghosh, RSC Adv., 2017, 7, 47818–47832.
12. A. Gemant, Physics, 1936, 7, 311–317.
13. K. Hyun, M. Wilhelm, C. O. Klein, K. S. Cho, J. G. Nam, K. H. Ahn, S. J. Lee, R. H. Ewoldt and G. H. McKinley, Prog. Polym. Sci., 2011, 36, 1697–1753.
14. S. Rogers, J. Rheol., 2012, 56, 1129.
15. M. C. Rogers, K. Chen, M. J. Pagenkopp, T. G. Mason, S. Narayanan, J. L. Harden and R. L. Leheny, Phys. Rev. Mater., 2018, 2, 095601.
16. Z. A. Chaleshtari and R. Foudazi, Soft Matter, 2023, 19, 8337–8348.
17. M. E. Helgeson, S. E. Moran, H. Z. An and P. S. Doyle, Nat. Mater., 2012, 11, 344–352.
18. R. Foudazi, S. Qavi, I. Masalova and A. Ya. Malkin, Adv. Colloid Interface Sci., 2015, 220, 78–91.
19. B. Valipour Goodarzi and R. Foudazi, Langmuir, 2025, 41, 13688–13704.
20. P. Saramito, J. Non-Newton. Fluid Mech., 2007, 145, 1–14.
21. G. J. Donley, P. K. Singh, A. Shetty and S. A. Rogers, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 21945–21952.
22. G. Pagani, M. Hofmann, L. Govaert, T. Tervoort and J. Vermant, J. Rheol., 2024, 68, 155–170.
23. J. Oldroyd, Proc. R. Soc. London, Ser. A, 1950, 200, 523–541.
24. in Rheology Series, ed. H. A. Barnes, J. F. Hutton and K. Walters, Elsevier, 1989, vol. 3, pp. 37–54.
25. C. J. Dimitriou and G. H. McKinley, J. Non-Newton. Fluid Mech., 2019, 265, 116–132.
26. D. Fraggedakis, Y. Dimakopoulos and J. Tsamopoulos, Soft Matter, 2016, 12, 5378–5401.
27. M. Raissi, P. Perdikaris and G. E. Karniadakis, arXiv, 2017, preprint, arXiv:1711.10561 DOI:10.48550/arXiv.1711.10561.
28. M. Raissi, P. Perdikaris and G. E. Karniadakis, arXiv, 2017, preprint, arXiv:1711.10566 DOI:10.48550/arXiv.1711.10566.
29. G. P. P. Pun, R. Batra, R. Ramprasad and Y. Mishin, Nat. Commun., 2019, 10, 2339.
30. D. Chaffart, Y. Yuan and L. A. Ricardez-Sandoval, J. Phys. Chem. C, 2024, 128, 3733–3750.
31. M. R. Babaei, R. Stone, T. A. Knotts IV and J. Hedengren, J. Chem. Theory Comput., 2023, 19, 4163–4171.
32. J. Pateras, C. Zhang, S. Majumdar, A. Pal and P. Ghosh, Sci. Rep., 2025, 15, 7980.
33. W. Ji, W. Qiu, Z. Shi, S. Pan and S. Deng, J. Phys. Chem. A, 2021, 125, 8098–8106.
34. J. Lee, S. Shin, T. Kim, B. Park, H. Choi, A. Lee, M. Choi and S. Lee, Sci. Rep., 2025, 15, 16740.
35. A. Farkane, M. Ghogho, M. Oudani and M. Boutayeb, Int. J. Numer. Methods Fluids, 2024, 96, 381–396.
36. S. Thakur, M. Raissi and A. M. Ardekani, J. Non-Newton. Fluid Mech., 2024, 330, 105265.
37. H. Yan, J. Ding and C. Tao, Fluids, 2025, 10, 184.
38. A. Jain, R. Gurnani, A. Rajan, H. J. Qi and R. Ramprasad, npj Comput. Mater., 2025, 11, 42.
39. M. Mahmoudabadbozchelou and S. Jamali, Sci. Rep., 2021, 11, 12015.
40. M. Lardy, S. Tlili and S. Gsell, J. Non-Newton. Fluid Mech., 2025, 346, 105512.
41. Z. Zhang, G. Yuan, Z. Qin and Q. Luo, Expert Syst. Appl., 2026, 296, 129002.
42. M. Raissi, P. Perdikaris and G. E. Karniadakis, J. Comput. Phys., 2019, 378, 686–707.
43. P. Rathore, W. Lei, Z. Frangella, L. Lu and M. Udell, in Proceedings of the 41st International Conference on Machine Learning, JMLR.org, Vienna, Austria, 2024, vol. 235, pp. 42159–42191.
44. A. Satyadharma, M.-J. Chern, H.-C. Kan, Harinaldi and J. Julian, Phys. Fluids, 2024, 36, 103619.
45. M. Saadat, D. Mangal and S. Jamali, Digital Discovery, 2023, 2, 915–928.
46. M. Mahmoudabadbozchelou, K. M. Kamani, S. A. Rogers and S. Jamali, Proc. Natl. Acad. Sci., 2022, 119, e2202234119.
47. D. Izbassarov and O. Tammisola, Phys. Rev. Fluids, 2020, 5, 113301.
48. M. E. Rosti, D. Izbassarov, O. Tammisola, S. Hormozi and L. Brandt, J. Fluid Mech., 2018, 853, 488–514.
49. K. Kamani, G. J. Donley and S. A. Rogers, Phys. Rev. Lett., 2021, 126, 218002.
50. I. Cheddadi, P. Saramito, B. Dollet, C. Raufaste and F. Graner, Eur. Phys. J. E, 2011, 34, 1.
51. M. S. Abdelgawad, S. J. Haward, A. Q. Shen and M. E. Rosti, Phys. Fluids, 2024, 36, 113127.
52. C. J. Dimitriou and G. H. McKinley, Soft Matter, 2014, 10, 6619–6644.
53. K. M. Kamani, G. J. Donley, R. Rao, A. M. Grillet, C. Roberts, A. Shetty and S. A. Rogers, J. Rheol., 2023, 67, 331–352.
54. Y. Wang, G. Maurel, M. Couty, F. Detcheverry and S. Merabia, Macromolecules, 2024, 57, 3636–3646.
55. A. Ghorai, Diffus. Fundam., 2022, 36, 1226.
56. R. B. Reboucas, N. N. Nikolova and V. Sharma, Curr. Opin. Colloid Interface Sci., 2025, 77, 101904.
57. P.-C. Zhao, W. Li, W. Huang and C.-H. Li, Molecules, 2020, 25, 597.
58. R. Subramani, A. Izquierdo-Alvarez, P. Bhattacharya, M. Meerts, P. Moldenaers, H. Ramon and H. Van Oosterwyck, Front. Mater., 2020, 7, 212.
59. I. Masalova and A. Ya Malkin, Colloid J., 2007, 69, 185–197.
60. J. Klein, Nature, 1978, 271, 143–145.
61. X. Fan, H. Xu, Q. Zhang, D. Xiao, Y. Song and Q. Zheng, Polymer, 2019, 167, 109–117.
62. F. Yziquel, P. J. Carreau and P. A. Tanguy, Rheol. Acta, 1999, 38, 14–25.
63. A. Y. Malkin, V. G. Kulichikhin, S. Y. Khashirova, I. D. Simonov-Emelyanov and A. V. Mityukov, Polymers, 2024, 16, 442.
64. M. Cloitre and R. T. Bonnecaze, Rheol. Acta, 2017, 56, 283–305.
65. D. A. Reinelt, J. Rheol., 1993, 37, 1117–1139.
66. P. Rognon and C. Gay, Eur. Phys. J. E, 2009, 30, 291.
67. R. Agrawal, E. García-Tuñón, R. J. Poole and C. P. Fonte, J. Non-Newton. Fluid Mech., 2025, 338, 105407.
68. R. B. Kudenatti and N.-E. Misbah, Sci. Rep., 2020, 10, 9445.
69. U. Eberhard, H. J. Seybold, M. Floriancic, P. Bertsch, J. Jiménez-Martínez, J. S. Andrade and M. Holzner, Front. Phys., 2019, 7, 71.
70. E. Maor, e: The Story of a Number, Princeton University Press, 2011.
71. S. Cooper, Theor. Biol. Med. Model., 2006, 3, 10.
72. L. J. Curtis, Am. J. Phys., 1978, 46, 896–906.
73. W.-H. Lin, E. Kussell, L.-S. Young and C. Jacobs-Wagner, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 27795–27804.
74. G. Karczewski, A. K. Zakrzewski, L. Dobaczewski, W. Dobrowolski, E. Grodzicka, J. Jaroszyński, T. Wojtowicz and J. Kossut, Thin Solid Films, 1995, 267, 79–83.
75. A. A. Istratov and O. F. Vyvenko, Rev. Sci. Instrum., 1999, 70, 1233–1257.
76. F. Rocadenbosch, A. Comeron and D. Pineda, Appl. Opt., 1998, 37, 2199–2206.
77. S. Paul and G. Gangopadhyay, Chem. Phys. Lett., 2003, 369, 643–649.
78. A. I. Hanopolskyi, V. A. Smaliak, A. I. Novichkov and S. N. Semenov, ChemSystemsChem, 2021, 3, e2000026.
79. F. Meng, R. H. Pritchard and E. M. Terentjev, Macromolecules, 2016, 49, 2843–2852.
80. M. E. Cates and S. M. Fielding, Adv. Phys., 2006, 55, 799–879.
81. S. Jamali, G. H. McKinley and R. C. Armstrong, Phys. Rev. Lett., 2017, 118, 048003.
82. H. M. Wyss, K. Miyazaki, J. Mattsson, Z. Hu, D. R. Reichman and D. A. Weitz, Phys. Rev. Lett., 2007, 98, 238303.
83. Y. E. Chou, T. H. Liu and C.-A. Lin, arXiv, 2025, preprint, arXiv:2509.21393 DOI:10.48550/arXiv.2509.21393.

---