On the origin of the negative energy-related contribution to the elastic modulus of rubber-like gels

TL;DR

This paper develops a coarse-grained polymer model to explain the negative energy contribution to the elastic modulus of rubber-like gels, revealing a temperature-dependent correction factor that improves understanding of gel elasticity.

Contribution

The study introduces an exact free energy expression for a polymer model that captures the negative energy contribution and proposes a modified elastic modulus formula incorporating a temperature-dependent factor.

Findings

The model accurately describes experimental data for tetra-PEG hydrogels.

The elastic modulus depends on temperature as G(T) ∝ k_B T w(T).

The correction factor w(T) relates to chain interactions and solvent effects.

Abstract

We consider a coarse-grained polymer model in order to investigate the origin of a recently discovered negative energy-related contribution to the elastic modulus $G(T)$ of rubber-like gels. From this model, we are able to compute an exact expression for the free energy of the system, which allows us to evaluate a stress-strain relationship that displays a non-trivial dependence on the temperature $T$. We validate our approach through comparisons between the theoretical results and the experimental data obtained for tetra-PEG hydrogels, which indicate that, although simple, the present model works well to describe the experiments. Importantly, our approach unveiled aspects of the experimental analysis which turned out to be different from the conventional entropic and energetic analysis broadly used in the literature. Also, in contrast to the linear dependence predicted by the traditional, {\it i.e.}, purely entropic, models, our results suggest that the general expression of the elastic modulus should be of the form $G(T) \propto k_BT w(T)$, with $w(T)$ being a temperature-dependent correction factor that could be related to the interaction between the chains in the network and the solvent. Accordingly, the correction factor allows the expression found for the elastic modulus to describe both rubber and rubber-like gels.

Equations

$$G(T)=n_e{k}_{B}T,$$

$$G(T)=a(T-T_0),$$

$$G(T)\propto k_{B}Tw(T),$$

$$\mathrm{\Omega }(N,n_b,n_s)=\frac{N!}{n_b!n_s!}.$$

$$\ell (n_s,n_b)=n_s{\ell }_s+n_b{\ell }_b,$$

$$U(n_s,n_b)=n_s{E}_s+n_b{E}_b.$$

$$p(N,n_b)=\frac{\mathrm{\Omega }(N,n_b)e^{-H(N,n_b)/k_{B}T}}{{\mathcal{Z}}_{\mathcal{G}}(N,T,f)},$$

$${\mathcal{Z}}_{\mathcal{G}}(N,T,f)={\left[e^{(f{\ell }_s-E_s)/k_{B}T}+e^{(f{\ell }_b-E_b)/k_{B}T}\right]}^N.$$

$$\mathcal{G}(N,T,f)=-N{k}_{B}T\ln \left[e^{(f{\ell }_s-E_s)/k_{B}T}+e^{(f{\ell }_b-E_b)/k_{B}T}\right],$$

$$\ell =-{\left[\frac{\partial }{\partial f}\mathcal{G}(N,T,f)\right]}_{N,T}=N\left[\frac{{\ell }_s+{\ell }_b{e}^{(\Delta E-f\Delta \ell )/k_{B}T}}{1+e^{(\Delta E-f\Delta \ell )/k_{B}T}}\right].$$

$${\ell }_0=N\left(\frac{{\ell }_s+{\ell }_b{e}^{\Delta E/k_{B}T}}{1+e^{\Delta E/k_{B}T}}\right).$$

$$\gamma =\frac{\ell -{\ell }_0}{{\ell }_0}.$$

$$f(T,\ell )=\frac{\Delta E}{\Delta \ell }+\frac{k_{B}T}{\Delta \ell }\ln \left(\frac{\ell -N{\ell }_b}{N{\ell }_s-\ell }\right).$$

$$\sigma (T,\gamma )=\frac{\Delta E}{A\Delta \ell }+\frac{k_{B}T}{A\Delta \ell }\ln \left(\frac{{\gamma }_b+\gamma }{{\gamma }_s-\gamma }\right),$$

$${\gamma }_s=\frac{N{\ell }_s-{\ell }_0}{{\ell }_0},$$

$${\gamma }_b=\frac{{\ell }_0-N{\ell }_b}{{\ell }_0},$$

$${\gamma }_b={\gamma }_s{e}^{-\Delta E/k_{B}T},$$

$$\sigma (T,\gamma )\simeq \frac{k_{B}T}{{\gamma }_{s}A\Delta \ell }\left(1+e^{\Delta E/k_{B}T}\right)\gamma .$$

$$G(T)=\underset{\gamma \to 0}{\lim }{\left[\frac{\partial }{\partial \gamma }\sigma (T,\gamma )\right]}_T.$$

$$G(T)=\frac{k_{B}T}{{\gamma }_{s}A\Delta \ell }\left(1+e^{\Delta E/k_{B}T}\right),$$

$$\sigma (T,\gamma )\simeq G(T)\gamma ,$$

$${\gamma }_s=\frac{({\ell }_s-{\ell }_b)e^{\Delta E/k_{B}T}}{{\ell }_s+{\ell }_b{e}^{\Delta E/k_{B}T}},$$

$$G(T)=\frac{k_{B}T}{A\Delta {\ell }^2}\left({\ell }_s+{\ell }_b{e}^{\Delta E/k_{B}T}\right)\left(1+e^{-\Delta E/k_{B}T}\right).$$

$$T_0^{\ast }=\frac{\Delta E}{k_B\ln \left(-{\ell }_s/{\ell }_b\right)},$$

$$\sigma (T,\gamma )={\sigma }_S(T,\gamma )+{\sigma }_E(T,\gamma ),$$

$$G(T)=G_S(T)+G_E(T).$$

$$G_S(T)=T\frac{dG(T)}{dT},$$

$$G_S(T)$$

$$G_E(T)=\frac{\Delta E}{A\Delta {\ell }^2}\left({\ell }_b{e}^{\Delta E/k_{B}T}-{\ell }_s{e}^{-\Delta E/k_{B}T}\right).$$

$${\sigma }_S(T,\gamma )\simeq G_S(T)\gamma ,$$

Full text

††preprint: arXiv:2302.12819

I Introduction

In contrast to many solid materials, rubber-like gels are semisolid viscoelastic materials formed by networks of disordered crosslinked polymeric chains immersed in liquid solvents which display isotropic and soft mechanical properties richtering2014softmatter . Examples annabi2014advmater ; thiele2014advmater ; kamata2015advhmater ; calo2015eurpolj ; basu2011macromol ; wen2012softmatter ; dai2021softmatter include not only biologically relevant but also synthetically designed gels, e.g., modified hyaluronic acids, polyethylene glycol (PEG), and polyacrylamide (PAM) gels. Due to the possibility of having a wide range of elasticity, mechanical strength, porosity, and swellability, hydrogels are materials that hold great promise in several applications gu2020angewchem ; danielsen2021chemrev such as, e.g., drug delivery systems and tissue engineering scaffolds.

Remarkably, because of the incompressibility that is generally observed in experiments yoshikawa2021prx , the mechanical response of rubber-like gels is usually characterized by a single elastic constant111For rubber-like gels, $G(T)$ corresponds to the shear modulus, which is usually much smaller than the bulk modulus $B(T)$, even so, the former governs their response to locally applied stresses and it also determines changes in their shape. $G(T)$ (see the discussions in Ref. treloarbook ). It is called the elastic modulus and it is usually assumed to be given by treloarbook ; toda2018aipadv

[TABLE]

where $k_{B}$ is the Boltzmann’s constant, $T$ is the absolute temperature of the solvent, and $n_{e}$ is the effective number density of “elastic elements” in solution. Different models treloarbook ; sakumichi2021polymj have different interpretations for $n_{e}$, but virtually all theoretical approaches assume that it should be proportional to the number density of precursor molecules in solution, $n$, since both the number density of strands (i.e., bridged chains), $n_{f}$, and the number density of crosslinks, $n_{c}$, should be proportional to $n$.

Interestingly, most of the theoretical approaches agreed on the linear dependence of $G(T)$ on the temperature $T$ given by Eq. 1. Essentially, this occurs because energetic contributions were neglected in the development of the seemingly successful models treloarbook ; sakumichi2021polymj , e.g., those that are based on the classical affine flory1953book , phantom network james1953jcp , and junction affine flory1977jcp models. In other words, those models consider that only entropic effects of the chains in the system are responsible to the observed temperature-dependent elastic modulus. It seems that such approximation works fine for natural and synthetic rubber treloarbook ; anthony1942jpc , but when a purely entropic approach is used to describe rubber-like gels, inconsistencies have been found between the theoretical predictions and the experimental results.

Indeed, recent experiments have shown that the linear elasticity of some gels can be fully described only when a significant negative contribution to its elastic modulus is considered yoshikawa2021prx ; fujiyabu2021prl . This issue was analyzed experimentally for tetra-PEG hydrogels that form networks with well-controlled topological structures yoshikawa2021prx ; sakumichi2021polymj , and a phenomenological correction to the linear description was proposed in order to consider such a negative contribution. In particular, instead of assuming $G(T)$ to be proportional to the absolute temperature $T$, as in Eq. 1, the authors in Ref. yoshikawa2021prx added an ad hoc term to the elastic modulus, writing it as

[TABLE]

where $a>0$ and $T_{0}>0$ are parameters that depend on the experimental conditions that were used to create the gels such as the concentration of precursor molecules in solution ($c$), their molar mass ($M$), and the connectivity ($p$) of the network yoshikawa2021prx ; sakumichi2021polymj . Since the first term in this equation is proportional to $T$, i.e., $G_{S}(T)=aT$, the authors of Ref. yoshikawa2021prx termed it as an entropic contribution to the elastic modulus, whereas the second one, $G_{E}(T)=-aT_{0}$, is its negative energy-related contribution, so that the total elastic modulus is given by $G(T)=G_{S}(T)+G_{E}(T)$.

Even though the authors have used Eq. 2 to describe the experimental results presented in Ref. yoshikawa2021prx , its origin is still an open issue. In our view, the qualitative microscopic model discussed in Ref. yoshikawa2021prx correctly presumes that the interaction between the chains in the network and their neighbouring solvent molecules should play an important role in the viscoelastic behaviour of rubber-like gels as the tetra-PEG hydrogels. However, we believe that the energy $E_{s}$ of the chain in its elongated (e.g., swelled, s) conformation should not be smaller than the energy $E_{b}$ of its collapsed (e.g., blob, b) conformation, that is, $E_{b}>E_{s}$, as assumed in the model discussed in Ref. yoshikawa2021prx . It is worth noting that this feature was also included as an effective self-repulsive interaction in the statistical mechanical model presented in Ref. shirai2022arxiv that attempts to explain the origin of the negative contribution to $G(T)$. Nonetheless, such qualitative behaviour contrasts with results observed from numerical simulations on folding phenomena junghans2006prl ; chen2008pre ; liu2012jcp ; frigori2013jcp , where $E_{b}<E_{s}$, i.e., chains in their elongated states have higher energies than the collapsed ones. Indeed, as it can be argued from the results presented latter in our work, one should have $E_{s}>E_{b}$ in order to explain the origin of the negative contribution $G_{E}(T)$ in rubber-like gels.

Accordingly, here we consider one of the simplest mesoscopic models which could be explored to explain the negative contribution $G_{E}(T)$ to the elastic modulus $G(T)$ of rubber-like hydrogels. In addition, such a simple model might provides us an interpretation for the phenomenological expression given by Eq. 2. In fact, our results indicate that it is possible to obtain expressions for the elastic modulus for these kind of gels in the form

[TABLE]

where $w(T)$ is a correction factor that makes the classical models able to describe not only purely entropic rubbers but rubber-like gels as well. This result is validated for the whole set of experimental data presented in Refs. yoshikawa2021prx , where tetra-PEG hydrogels were studied.

The remainder of the manuscript is as follows. In Sec. II we present the model and discuss how a force-extension relationship can be obtained from its thermostatistics. Most of the results related to the stress-strain relationship and to the elastic moduli of the model, as well as their experimental validation, are presented in Sec. III. Also, in Sec. III we discuss how the theoretical results obtained from this simple model might be related to the phenomenological expression given by Eq. 2. Finally, our main conclusions are presented in Sec. IV.

II Model and its thermostatistics

For simplicity, we consider a coarse-grained model where a bridged chain is assumed to have $N$ segments that can be only in one of two possible conformational motifs: the elongated, or swelled (s), and the collapsed, or blob (b), states. As showed in Fig. 1, the swelled state corresponds to a longer conformation with length $\ell_{s}$ and energy $E_{s}$, while the blob state is characterized by a conformation with length $\ell_{b}$ and energy $E_{b}<E_{s}$, so that $\Delta\ell=\ell_{s}-\ell_{b}>0$ and $\Delta E=E_{s}-E_{b}>0$. As mentioned before, the assumption that $E_{s}>E_{b}$ is supported by numerical simulations performed with effective coarse-grained models that assume implicity interactions with the solvent to describe the collapse of the chain (see, e.g., Refs. junghans2006prl ; chen2008pre ; liu2012jcp ; frigori2013jcp ). The negative value222Actually, the negative value of $\ell_{b}$ is not an issue if one thinks about this model as representing a one-dimensional random walk with different energy costs for the backward and the forward steps of different sizes. In fact, the one-dimensional freely-jointed chain (FJC) model kubobook is a special case of the model presented here when one imposes $E_{s}=E_{b}$ and $\ell_{s}=-\ell_{b}$. of $\ell_{b}$ is explained in Fig. 1, and, as we will show later, it is a necessary feature for the model in order to have the possibility of having zero end-to-end distances, i.e., $\ell=0$. This is interesting because it seems to be the main qualitative behaviour observed for Gaussian chains that have been used to describe rubber-like gels nishi2015jcp . It is worth noting that, for simplicity, we assume that the given energies, $E_{s}$ and $E_{b}$, and lengths, $\ell_{s}$ and $\ell_{b}$, of the individual segments in the two conformational states do not depend neither on the force $f$ or on the temperature $T$ of the thermal bath (i.e., the portion of the solvent which does not interact directly with the chain). Obviously, the end-to-end distance of the chain will depend on both thermodynamic variables.

Now, by considering that the chain has $N=n_{s}+n_{b}$ segments, with $n_{s}$ and $n_{b}$ being the number of segments in the swelled and in the blob states, respectively, the number of possible configurations of the system is given by

[TABLE]

Thus, the end-to-end distance of the chain can be written as

[TABLE]

and its internal energy, which includes the interaction between the chain and their neighboring solvent molecules, is given by

[TABLE]

The thermostatistics of such simple model can be determined either on the constant-stress ensemble, where the force $f$ is fixed and the end-to-end distance $\ell$ can fluctuate, or on the constant-strain ensemble, where $\ell$ is fixed and $f$ can vary. Since both approaches lead to exactly the same force-extension relation, in the following we will present the results in terms of the constant-stress ensemble, while in the Appendix A we demonstrate how the same relation can be obtained from the constant-strain ensemble.

II.1 Relationship between force and extension

By considering that the chain is subjected to a constant force $f$ at its ends, and by suppressing the dependence on $n_{s}$ in Eqs. 4, 5, and 6, since $n_{s}=N-n_{b}$, one can write the probability density function related to finding a specific state of the system as

[TABLE]

where the function $H(N,n_{b})=U(N,n_{b})-f\ell(N,n_{b})$ can be thought as an enthalpy. The normalization constant corresponds to the partition function of the system and is given by $\mathcal{Z}_{\mathcal{G}}(N,T,f)=\sum_{n_{b}=0}^{N}\mathrm{\Omega}(N,n_{b})\,e^{-H(N,n_{b})/k_{B}T}$, which yields kubobook

[TABLE]

Hence, the above result can be used to compute the following free energy (i.e., Gibbs thermodynamic potential) $\mathcal{G}(N,T,f)\equiv-k_{B}T\,\ln[\mathcal{Z}_{\mathcal{G}}(N,T,f)]$, which leads to

[TABLE]

from where one can obtain all the thermostatistical properties of the system. In particular, one can estimate the average value of the end-to-end distance of the chain as333For the free energy defined by the potential $\mathcal{G}$ one has the differential form $d\mathcal{G}=-SdT-\ell df$ where $S$ and $\ell$ are the entropy and the end-to-end distance of the system, respectively.

[TABLE]

It is worth noting that, at finite temperatures, this result indicates that one may have a non-zero end-to-end distance even in the absence of forces. So, we use Eq. 10 with $f=0$ to define the natural (i.e., relaxed) end-to-end distance as $\ell_{0}\equiv\ell(N,T,0)$, that is,

[TABLE]

Thus, Eqs. 10 and 11 can be used to define the strain $\gamma\equiv\gamma(T,\ell)$ through the following relation

[TABLE]

We note that, since $\ell\equiv\ell(N,T,f)$ and $\ell_{0}\equiv\ell(N,T,0)$, the strain $\gamma$ is a function of the temperature and the force applied to the chain. As expected, this definition is consistent with having $\gamma=0$ if $f=0$ for any value of the temperature $T$. Also, expression 12 indicates that the strain $\gamma$ can be either positive or negative, as it only depends whether the chain is being stretched ($\ell>\ell_{0}$) or compressed ($\ell<\ell_{0}$).

In the constant-stress ensemble, one can evaluate the force in terms of the end-to-end distance by inverting Eq. 10, which yields the seemingly known result given by kubobook

[TABLE]

It is worth noting that, as we mentioned before, here we assume that $\ell_{b}<0$, so that $\ell$ corresponds to the end-to-end distance and not the total (i.e., contour) length of the chain. Besides its different interpretation, we indicate that, in that case, Eq. 13 is physically well-defined only for $0\,\leq\,\ell<N\ell_{s}$. In addition, since we assume that $\ell_{b}<0$ and require that the natural end-to-end distance of the chain $\ell_{0}$ to be positive by definition, one has to consider that Eqs. 11, 12, and 13 will have a precise physical meaning only for temperatures $T$ higher than a characteristic temperature $T_{0}^{*}$, which is found from the condition that $\ell(N,T_{0}^{*},0)=0$. The consequences of this condition are better discussed in the Sec. III.4.

III Results

In this section we discuss how the results obtained for the present model can be used to provide a suitable framework to describe the experimental data extracted from Ref. yoshikawa2021prx for tetra-PEG hydrogels. In particular, we will focus on the mechanical response of the gel networks described by stress-strain curves and also on the temperature-dependent behaviour of the elastic modulus.

III.1 $(T,\gamma)$-dependent stress

Accordingly, in order to link the force-extension behaviour of a single chain to the response of the whole system, we assume that the force should be related to the stress as $f=A\sigma$, where $A$ is an effective cross-sectional area444Here, the effective cross-sectional area $A$ related to the volume occupied by a bridged chain is assumed to be independent of the temperature $T$, however, as it will be argued later, it might depend on other quantities such as the concentration of precursor molecules ($c$) in solution, the molar mass ($M$) of such molecules, and the connectivity ($p$) of the network. Also, since we are dealing with an implicit averaging procedure that should be done over different orientations of the chains of a disordered isotropic system, possible differences related to the way one applies the stress and, consequently, to the way one defines the elastic modulus, are immaterial and should be incorporated into $A$. perpendicular to the direction of the force. Hence, by considering Eqs. 12 and 13 one can write the temperature-dependent stress in terms of the strain as555As shown in Appendix B, this result can be also obtained directly from the derivative of the Helmholtz free energy density and it corresponds to an engineering stress treloarbook .

[TABLE]

where

[TABLE]

and

[TABLE]

with $\ell_{0}$ given by Eq. 11.

We note that the above definitions of $\gamma_{s}$ and $\gamma_{b}$ are consistent with $\sigma(T,0)=0$, thus, from Eq. 14, one finds that

[TABLE]

which can be also verified through Eqs. 15 and 16 with the definition of $\ell_{0}$ given by Eq. 11. Also, it is worth noting that $\gamma_{s}\equiv\gamma_{s}(T)$ and $\gamma_{b}\equiv\gamma_{b}(T)$ will both present a dependence on the temperature because $\ell_{0}\equiv\ell(N,T,0)$, however, since they are defined as ratios between lengths, the stress $\sigma(T,\gamma)$ will not depend on the value of the number of segments $N$. As we will see later, the stress-strain relation given by Eq. 14 leads $\sigma(T,\gamma)$ to display a non-linear behaviour as a function of the temperature.

Interestingly, for small strains, one can assume that $\gamma/\gamma_{s}\ll 1$ and $\gamma/\gamma_{b}\ll 1$ (as a consequence of Eq. 17). Hence, by considering the linearization $\ln(1+x)\simeq x$, the stress given by Eq. 14 can be approximated by

[TABLE]

III.2 $T$-dependent elastic modulus

By considering the stress $\sigma(T,\gamma)$, one can evaluate the temperature-dependent elastic modulus as yoshikawa2021prx

[TABLE]

Hence, using the result 14 one finds that

[TABLE]

which is exactly the same result that would be obtained if one had replaced the linearized expression for stress, Eq. 18, in the definition of $G(T)$ given by Eq. 19. Accordingly, the above expression can be identified as the temperature-dependent pre-factor in Eq. 18, which can be rewritten as

[TABLE]

demonstrating that, at least in the linear regime, the limit of taking a small value for $\gamma$ is interchangeable with the partial derivative in Eq. 19.

Now, by considering Eqs. 11 and 15, one may write the full temperature dependence of $\gamma_{s}\equiv\gamma_{s}(T)$ as

[TABLE]

so that the explicit temperature-dependent expression for the elastic modulus given by Eq. 20 can be written as

[TABLE]

III.3 Stress $\sigma(T,\gamma)$ vs. temperature

In order to validate our approach, we include in Fig. 2 a comparison between the stress-strain relation given by Eq. 14 and the experimental data extracted from Ref. yoshikawa2021prx for tetra-PEG hydrogels. By trying several sets of parameters, we found that the parameters $\Delta E=3.4\,$pN.nm, $\ell_{s}=29\,$nm, $\ell_{b}=-4\,$nm, and $A=11.575\,\textrm{nm}^{2}$ are the best ones in order to describe, simultaneously, the whole set of experimental data. Indeed, as it can be seen from the results presented in Fig. 2 for $\sigma(T,\gamma)$ as a function of the temperature $T$ and several different values for the strain $\gamma$, the stress-strain relation given by Eq. 14 obtained here displays a good agreement when compared to the data extracted from the experiments.

Also, the results presented in Fig. 2 indicates that $\sigma(T,\gamma)$ may not display a linear behaviour as a function of the temperature $T$, as proposed in Ref. yoshikawa2021prx and expected from the phenomenological model defined by Eq. 2. In fact, Eq. 14 indicates that this dependence is slightly non-linear since the parameters $\gamma_{s}$ and $\gamma_{b}$ inside the logarithm term are defined by Eqs. 17 and 22, with both displaying a non-linear dependence on the temperature.

III.4 Temperature $T_{0}^{*}$

Interestingly, Fig. 2 indicates that the stress-strain relation given by Eq. 14 displays a characteristic temperature $T_{0}^{*}$ where $\sigma(T_{0}^{*},\gamma)=0$. Although it might be unphysical to reach that temperature in practice because many solvents like water will be frozen at such a low temperature, we discuss the condition where the temperature $T$ is equal to $T_{0}^{*}$ since it might shed light on the phenomenological temperature $T_{0}$ defined in Eq. 2 (see also Ref. yoshikawa2021prx ).

We first note that having $\sigma(T_{0}^{*},\gamma)=0$ is the same as having that the force $f$ equal to zero at this specific temperature, since the force and stress are related as $f=A\sigma$. Then, as suggested in the inset of Fig. 2, where $\ell_{0}/N$ is plotted as a function of the temperature $T$, we argue that the condition $f\approx 0$ (i.e., $\sigma\approx 0$) will be satisfied when the natural end-to-end distance tends to zero at the temperature $T_{0}^{*}$, that is, $\ell_{0}^{*}\equiv\ell(N,T^{*}_{0},0)\approx 0$. In particular, by imposing that $\ell_{0}=0$ and considering Eq. 11, one finds the temperature

[TABLE]

which will be positive if $\Delta E>0$, $\ell_{b}<0$, and $\ell_{s}>|\ell_{b}|$, as in the case we are considering here. At this temperature, one has that the ratio determined by Eq. 17 is given by $\gamma_{s}(T_{0}^{*})/\gamma_{b}(T_{0}^{*})=-\ell_{s}/\ell_{b}$. Also, from Eq. 12, one has that $\ell^{*}\equiv\ell(N,T^{*}_{0},f)=(1+\gamma)\ell^{*}_{0}\approx 0$. Hence, the force computed by Eq. 13 should be given by $f(T_{0}^{*},\ell^{*})=(\Delta E/\Delta\ell)-(k_{B}T_{0}^{*}/\Delta\ell)\ln\left(-\ell_{s}/\ell_{b}\right)$, which is zero according to the $T_{0}^{*}$ defined by expression 24. This means that, at the temperature $T^{*}_{0}$, the stress $\sigma(T_{0}^{*},\gamma)$ must be zero as well. In fact, the condition where $\sigma(T_{0}^{*},\gamma)=0$ can be possible if the stress in Eq. 14 is independent of strain. It can happen, in particular, if the argument inside the logarithm term in expression 14 tends to a finite constant value while $\gamma_{s}(T_{0}^{*})\gg\gamma$ and $\gamma_{b}(T_{0}^{*})\gg\gamma$. In that case, one will have $(\gamma_{b}(T_{0}^{*})+\gamma)/(\gamma_{s}(T_{0}^{*})-\gamma)\approx\gamma_{b}(T_{0}^{*})/\gamma_{s}(T_{0}^{*})\approx-\ell_{b}/\ell_{s}$, where in the last step we make use of Eq. 17. As one may check, by replacing this result and $T_{0}^{*}$ given by Eq. 24 in Eq. 14, it yields the behaviour observed in Fig. 2, that is, $\sigma(T_{0}^{*},\gamma)=0$ at $T_{0}^{*}=124.3\,$K.

III.5 $S$ and $E$ contributions

As in Ref. yoshikawa2021prx , we assume that both the stress and the elastic modulus can be defined as a sum of two contributions, that is,

[TABLE]

and

[TABLE]

In principle, the authors of Ref. yoshikawa2021prx used the subscripts $S$ and $E$ to denote entropic and energetic contributions, respectively. Although we keep that notation to facilitate the comparison to experiments, their approach indicate that, in contrast to the conventional definitions used in the literature flory1953book , both $\sigma_{E}(T,\gamma)$ and $G_{E}(T)$ should include not only terms which are related to the derivative of the internal energy but also terms related to the derivative of the entropy, that is, $f_{E}=A\sigma_{E}\neq(\partial U/\partial\ell)_{N,T}$ (see Appendix C for more details). Below we explicitly derive the expressions for both contributions.

III.5.1 Elastic moduli

The $S$ contribution to the elastic modulus can be evaluated from Eq. 23 using the relation yoshikawa2021prx

[TABLE]

which yields

[TABLE]

The first term in Eq. 28 is just the elastic modulus $G(T)$ given by Eq. 23, thus, by considering Eq. 26, the complementary contribution is easily obtained as

[TABLE]

Clearly, if $\Delta E=0$, there will be no energy-related contribution to the elastic modulus from $G_{E}(T)$, however, if one assumes that $\Delta E>0$, $\ell_{s}>0$, and $\ell_{b}<0$, Eq. 29 indicates that $G_{E}(T)$ is a negative quantity for all values of $T>T^{*}_{0}$. In addition, by considering relatively small strains, linearized response functions can be obtained for the $S$ and $E$ contributions to the stress defined by Eq. 25, that is,

[TABLE]

and

[TABLE]

Hence, if $G_{E}(T)<0$, one will observe a negative contribution to the stress, i.e., $\sigma_{E}(T,\gamma)<0$, while the $S$ contribution to the stress, $\sigma_{S}(T,\gamma)$, will be positive for all $T>T^{*}_{0}$, since $\sigma(T,\gamma)\simeq G(T)\gamma$ must be also positive. In order to verify this behaviour for the general case, that is, not only at the linear response regime but also at higher values of $\gamma$, next we calculate the full expressions of $\sigma_{S}(T,\gamma)$ and $\sigma_{E}(T,\gamma)$.

III.5.2 Positive ($S$) contribution to the stress $\sigma(T,\gamma)$

Following the approach presented in Ref. yoshikawa2021prx , the $S$ contribution to the stress can be evaluated as

[TABLE]

so that Eq. 14 yields

[TABLE]

with $\gamma_{b}$ and $\gamma_{s}$ given by Eqs. 17 and 22, respectively, whereas

[TABLE]

is the coefficient of thermal expansion of the chains in the gel, with $G(T)$ being the elastic modulus given by Eq. 23.

It is worth mentioning that, under tension, i.e., for $f\neq 0$, the thermal expansion behaviour of the chains is governed by $\alpha(T,f)=\ell^{-1}[\partial\ell/\partial T]_{N,f}$, with the end-to-end distance $\ell$ given by Eq. 10, which can be explicitly evaluated as

[TABLE]

Thus, if $f>\Delta E/\Delta\ell$, it is possible that the chains react as predicted by the Gough-Joule effect with $\alpha(T,f)$ being a negative function, i.e., it is possible that tensioned chains contract when the temperature is increased, as it seems to occur in tetra-PEG hydrogels explored in Ref. yoshikawa2021prx .

III.5.3 Negative ($E$) contribution to $\sigma(T,\gamma)$

Now, by considering Eq. 25, that is, $\sigma_{E}(T,\gamma)=\sigma(T,\gamma)-\sigma_{S}(T,\gamma)$, one can consider Eqs. 14 and 33 to identify the negative contribution to the stress as

[TABLE]

with the usual definitions of $\gamma_{b}$ (Eq. 17) and $\gamma_{s}$ (Eq. 22), and $\alpha(T)$ given by Eq. 34.

III.6 Stresses $\sigma$, $\sigma_{S}$, and $\sigma_{E}$ vs. strain

Figure 3 includes the comparison between our theoretical results and the experimental data extracted from Ref. yoshikawa2021prx determined for tetra-PEG hydrogels at a constant temperature. In Fig. 3 one sees how well the non-linear relations found for the stresses, given by Eqs. 14, 33 and 36, describe the experimental results obtained as a function of the strain $\gamma$ at $T=288\,$K. This corroborates expression 14 for the total stress and validates its consequences, i.e., Eqs. 33 and 36, which were determined for its $S$ and $E$ contributions, respectively. Importantly, the results in Fig. 3 show that, even when the full dependence of the contribution $\sigma_{E}(T,\gamma)$ to the stress as a function of $\gamma$ is considered, a negative elastic contribution is observed for strains $\gamma$ as higher as 1.6, where the system is already at the non-linear response regime.

Furthermore, the results presented in Fig. 3 indicate that, when the positive ($S$) and the negative ($E$) contributions to the elastic modulus are given by Eqs. 28 and 29, respectively, the linearized relations (dash-dotted lines) evaluated through Eqs. 30 and 31 also describe reasonably well the experimental data. Even so, it is worth noting that Eqs. 21, 30, and 31 are valid only if $\gamma<\gamma_{s}$, since the non-linear expressions for $\sigma(T,\gamma)$, $\sigma_{S}(T,\gamma)$, and $\sigma_{E}(T,\gamma)$, all diverge for $\gamma=\gamma_{s}$. Thus, the value of $\gamma_{s}$ establishes the maximum strain that the chains in the network can be submitted at a temperature $T$.

III.7 Stresses $\sigma$, $\sigma_{S}$, and $\sigma_{E}$ vs. temperature

Now, in order to illustrate how the stress and its positive ($S$) and the negative ($E$) contributions behave as functions of the temperature but at a fixed strain, we include in Fig. 4 a comparison between our theoretical results and the experimental data obtained for tetra-PEG hydrogels extracted from Ref. yoshikawa2021prx . Accordingly, one can observe a good agreement between the experimental data (filled symbols) and our non-linear expressions (continuous lines), i.e., Eqs. 14, 33 and 36. Although qualitatively, a similar behaviour is observed for the linearized relations (dash-dotted lines), that is, Eqs. 21, 30, and 31. Besides the results obtained from our model (continuous and dash-dotted lines), Fig. 4 also includes the results obtained with the phenomenological model (dashed straight lines) proposed in Ref. yoshikawa2021prx , where $\sigma(T,\gamma)=G(T)\gamma$, with $G(T)$ defined by Eq. 2.

We note that, in contrast to what is expected from the phenomenological model, the positive contribution to the stress, i.e., $\sigma_{S}(T,\gamma)$ given by Eq. 30, is not a linear function of the temperature. In that model, the $S$ contribution to the stress is a bi-linear function of the temperature and the strain, that is, $\sigma_{S}(T,\gamma)=G_{S}(T)\gamma=aT\gamma$. As a consequence, the $S$ contribution to the elastic modulus is simply $G_{S}(T)=aT$, which is clearly different from what is found from our model, where both $\sigma_{S}(T,\gamma)$ and $G_{S}(T)$ are more complicated functions of the temperature. Furthermore, as shown in Fig. 4, the negative contribution to the stress, $\sigma_{E}(T,\gamma)$, and to the elastic modulus, $G_{E}(T)$, in our model are also temperature-dependent, while in the phenomenological model those contributions are constant through temperature variations, i.e., $\sigma_{E}(T,\gamma)=G_{E}(T)\gamma=-aT_{0}\gamma$.

Yet, one might expect that, even with such differences, for small strains and relatively high temperatures, both models might describe the experimental data well. As we discuss next, such equivalence can be used to establish the dependence of the parameters defined in our model and the experimental conditions.

III.8 Relationship between our model and the

phenomenological model of Ref. yoshikawa2021prx

To understand how the parameters of our model might be related to the model and results presented in Refs. yoshikawa2021prx , we first consider Eq. 14 for small strains and relatively high temperatures, i.e., the case where $\gamma/\gamma_{s}\ll 1$ and $\Delta E/k_{B}T$ $\ll 1$. In this case, one can assume Eq. 21 to be valid, so that the elastic modulus, given by Eq. 23, can be linearized by assuming that $e^{x}\simeq 1+x$. Hence, by neglecting second-order terms in $x=\Delta E/k_{B}T$, one finds that

[TABLE]

Remarkably, with such linearization, the approximated elastic modulus $G(T)$ obtained from our model seems to exhibit the same linear behaviour of the phenomenological model that was proposed in Ref. yoshikawa2021prx . Thus, by comparing Eq. 2 and Eq. 37, one finds the following relations

[TABLE]

and

[TABLE]

Now, it is instructive to inspect what happens if the parameters obtained from the heuristic non-linear fit of Eqs. 14, 33, and 36 to the data presented in Figs. 2, 3, and 4 are plugged into the above expressions. In particular, by assuming $\Delta E=3.4\,$pN.nm, $\ell_{s}=29\,$nm, $\ell_{b}=-4\,$nm, and $A=11.575\,$nm2, one finds from the above relations that $a=5.47\times 10^{-2}\,$kPa/K and $T_{0}=162.6\,$K, which overestimates the values $a=(5.0\pm 0.2)\times 10^{-2}$ kPa/K and $T_{0}=(131\pm 8)\,$K obtained from the fit of the data to the linear expressions of the phenomenological model yoshikawa2021prx , i.e., those determined from Eq. 2. Nevertheless, because this correspondence is not valid in general, but only in a limiting case, one may not expect that these relations will give consistent values for $a$ and $T_{0}$. Indeed, the temperature $T_{0}=162.6\,$K (Eq. 39) estimated from the parameters obtained with the non-linear expressions is different not only from the characteristic temperature $T_{0}^{*}=124.3\,$K (Eq. 24) of our non-linear model, but also from the value determined from the phenomenological model, i.e., $T_{0}=(131\pm 8)\,$K. Interestingly, the stresses obtained through the linear approximation of $G(T)$ given by Eq. 37 seems to fit the experimental data much better if small changes are made to the parameters. In particular, by considering $\Delta E=2.7\,$pN.nm, $\ell_{s}=30\,$nm, $\ell_{b}=-4.35$ nm (i.e., $\Delta\ell=34.35\,$nm), and $A=11.575\,\textrm{nm}^{2}$, the relations 38 and 39 yield $a=5.18\times 10^{-2}\,$kPa/K and $T_{0}=131\,$K, which are closer to the parameters obtained in Ref. yoshikawa2021prx and that were used to obtain the curves for different strains presented in Fig. 2 as the dashed lines. However, it is important noting that, if this last set of parameters is plugged into the non-linear expression, Eq. 14, one gets $T_{0}^{*}=101.32\,$K (Eq. 24) and the agreement with experimental data is poor, which indicates that the equivalence between the $G(T)$ given by Eq. 37 and the elastic modulus of the phenomenological model (Eq. 2) is somewhat fortuitous due to the linearization procedure.

As it can be inferred from Refs. yoshikawa2021prx ; sakumichi2021polymj , the experimental results presented in Figs. 2, 3, and 4 were all measured at the same experimental condition, which is characterized by gels with network connectivity $p=0.915$ assembled from tetra-functional precursor molecules with molar mass $M=20\,$kg/mol at a concentration $c=60\,$kg/m3. In the following we will explore an analysis based on the experimental data extracted from Refs. yoshikawa2021prx ; sakumichi2021polymj in order to understand how the parameters $\ell_{s}$, $\ell_{b}$, $\Delta E$, and $A$ might depend on $p$, $M$, and $c$.

III.8.1 Energy difference $\Delta E$ and temperature $T_{0}^{*}$

We first consider how the energy difference $\Delta E$ changes for different values of $p$, $M$, and $c$. In fact, in the following we replace the concentration of precursor molecules $c$ by its number density using the relation

[TABLE]

where $N_{A}$ is the Avogadro’s constant. It was established from experimental data analyses in Refs. yoshikawa2021prx ; sakumichi2021polymj that $T_{0}\equiv T_{0}(c,M)$, i.e., the phenomenological temperature seems to be independent of the connectivity $p$. Also, it was suggested in Ref. yoshikawa2021prx that $T_{0}\propto M^{-1/2}\,n^{-1/3}$ for lower values of $M$ and $n$. Now, since Eq. 39 establishes that $\Delta E\propto T_{0}$, one needs to discuss how $\delta\ell/\Delta\ell$ depends on the experimental conditions, with $\delta\ell=\ell_{s}+\ell_{b}$ and $\Delta\ell=\ell_{s}-\ell_{b}$. In particular, by assuming that both $\ell_{s}$ and $\ell_{b}$ display the same dependency666The specific dependence of $\ell_{s}$ and $\ell_{b}$ on $M$ might be inferred from known polymer models rubinsteinbook such as, e.g., the Gaussian chain, the worm-like chain (WLC), the and freely-jointed chain (FJC). For instance, one can assume $M_{f}=2(M/z)$ as the molar mass of bridged chains between crosslinks of functionality $z$, so that the lengths $\ell_{s}$ and $\ell_{b}$ would follow the end-to-end distance dependence on $M_{f}$, which is a power-law of the type $\ell\propto(M_{f})^{\nu}$, with $\nu$ being the Flory exponent flory1953book . on the molar mass $M$ and no dependency on the number density of precursor molecules $n$, one finds that the ratio $\delta\ell/\Delta\ell$ will be independent of $M$ and $n$. Hence, one can infer that the energy difference behaves in the same manner as $T_{0}$, that is, $\Delta E\propto M^{-1/2}n^{-1/3}$. In Fig. 5(a) we include the right side of Eq. 39 as a function of the concentration $c$ for different values of $M$, while in Fig. 5(b) we indicate that $\Delta E$, just as $T_{0}$, should be indeed proportional to $M^{-1/2}n^{-1/3}$, as suggested in Ref. yoshikawa2021prx .

Finally, it is worth noting that, because of the definition of $T_{0}^{*}$, Eq. 24, which involves only the ratio $\ell_{s}/\ell_{b}$, one will also have $T_{0}^{*}\propto\Delta E$, so that this temperature is expected to display the same dependence on the molar mass and number density of the phenomenological temperature $T_{0}$, that is, $T_{0}^{*}\propto M^{-1/2}n^{-1/3}$. Although systematic procedures based on numerical simulations of polymer systems suggest that, indeed, coarse-grained potentials should depend on the number density reith2003jcomputchem , we are not aware of any quantitative relationship between $\Delta E$ and either $M$ or $n$ that is obtained from first principles.

III.8.2 Pre-factor of $G(T)$ for different connectivities $p$

Accordingly, Ref. yoshikawa2021prx presents an extensive study where experiments with gels formed by networks with tetra-PEG chains were performed in a broad range of experimental conditions. In particular, the study considered different connectivities $p$ while keeping the same molar mass $M$ and number density $n$, which allowed us to explore their results in order to infer how the pre-factor related to $A$ in Eq. 23 depends on $n$, $M$, and $p$.

From the theoretical point-of-view, diluted networks with different connectivities $p$ can be described by the effective medium theory based on regular networks presented in Ref. nishi2015jcp . In this case, the elastic modulus can be written as

[TABLE]

with $G_{c}(T)$ being the elastic modulus of the fully connected network, i.e., $p=1$, and $g_{z}(p)$ being a function that depends on the connectivity $p$ of the network. For sufficiently highly connected networks, the latter can be written as a linear function as nishi2015jcp

[TABLE]

where $p^{*}_{z}=2/z$ is a parameter related to the coordination number (or functionality) $z$ of the crosslinks in the idealized complete network (e.g., $z=4$ and $z=6$ for square and cubic lattices, respectively). Since the experiments were performed with the tetra-PEG gels sakai2008macromo ; yasuda2009macromo where the precursor molecules present a coordination equal to $z=4$, one expects that $p_{z}^{*}=1/2$ and $g_{z}(p)=2p-1$. Furthermore, as suggested in Refs. yoshikawa2021prx ; sakumichi2021polymj , the elastic modulus of a fully connected gel network can be approximated by

[TABLE]

where $n$ is the number density given by Eq. 40, $T_{0}$ is the phenomenological temperature, and $w^{*}$ is a positive phenomenological parameter that may be found from experimental analysis (in Refs. yoshikawa2021prx ; sakumichi2021polymj , the authors found its value to be $2.4$). By considering the elastic modulus $G(T)$ given by the phenomenological model, Eq. 2, and the definitions 41 and 43, one finds that

[TABLE]

which is corroborated by the linear behaviour observed for $a=2k_{B}\delta\ell/A\Delta\ell^{2}$ presented in Fig. 5(c). Hence, by taking the data from Ref. yoshikawa2021prx and reanalyzing it consistently with $g_{z}(1)=1$, we found from Fig. 5(c) that the phenomenological parameter $w^{*}$ is closer to $2.17$. In Fig. 5(d) we plot the function $g_{z}(p)$ found from the fit as a function of $p$, from where one can verify that the linear expression given by Eq. 42 describes the experimental data of Ref. yoshikawa2021prx well. Also, the fit of Eq. 42 to the data presented in Fig. 5(d) imposing that $g_{z}(1)=1$ yields $p^{*}_{z}\simeq 0.47$, which is remarkably close to the value that is expected for tetra-PEG networks with $z=4$.

It is worth noting that, through the identity 44 and the relation based on the phenomenological model given by Eq. 38, one finds that

[TABLE]

As discussed before, we assume that $\delta\ell/\Delta\ell=(\ell_{s}+\ell_{b})/(\ell_{s}-\ell_{b})$ is independent of $p$ and $M$, thus Eq. 45 implies that $A\Delta\ell\propto[g_{z}(p)n]^{-1}$. Even though one cannot tell exactly how the effective area $A$ changes with the different values of $p$, $n$, and $M$, since it is not possible to know exactly how $\Delta\ell$ changes with those quantities, one can infer, at least, that the product of those parameters is inversely proportional to $g_{z}(p)$ and the number density $n$.

In addition, one should note that, although we have considered specific values for $\ell_{s}$, $\ell_{b}$, and $A$, there are, in fact, only two parameters (besides $\Delta E$) that are relevant to the behaviour of the stresses and the elastic moduli, which are the ratio $q_{\ell}^{\,}=-\ell_{s}/\ell_{b}$ and the product $r_{\ell}^{\,}=A\Delta\ell$. For instance, by setting $q_{\ell}^{\,}=7.25$ and $r_{\ell}^{\,}=381.975\,$nm3, one could have used half of the lengths, that is, $\ell_{s}=14.5\,$nm and $\ell_{b}=-2\,$nm, then twice the value of the effective area, i.e., $A=23.150\,$nm2, in order to obtain exactly the same curves showed in Figs. 2, 3, and 4, since the values of $q_{\ell}^{\,}$ and $r_{\ell}^{\,}$ remain the same. As it happens, Eq. 45 establishes an interesting relation between $q_{\ell}^{\,}$ and $r_{\ell}^{\,}$ that can be conveniently explored in further studies to describe other experimental data sets. It is worth mentioning that the values $\ell_{s}=29\,$nm and $\ell_{b}=-4\,$nm are consistent with the estimates obtained from PEG chains used in single-molecule force-extension experiments kolberg2019JACS .

III.9 Temperature-dependent correction factor $w(T)$

It is important noting that the authors in Refs. yoshikawa2021prx ; sakumichi2021polymj did not explained the origin of the phenomenological correction $w^{*}$ in Eq. 43, and so in the expressions 44 and 45. Here we argue that it could be effectively described as a temperature-dependent correction factor to the elastic modulus $G(T)$. In particular, we define

[TABLE]

which should be constant if no temperature-dependent correction is needed to describe the experimental data, as in the case of the elastic modulus $G(T)$ obtained considering only the entropic contributions to the system, i.e., $G(T)\propto k_{B}T$. As one can see in Fig. 5(e), we present the plot of the ratio defined by Eq. 46 as a function of the temperature extracted from the data of Fig. 2 (here, only the data that correspond to strains $\gamma=0.2$ and $\gamma=0.4$ were considered to calculate the experimental elastic modulus $G(T)\simeq\sigma(T,\gamma)/\gamma$). In contrast to what is expected from the purely entropic models, Fig. 5(e) shows that the ratio $w(T)$ is not constant, indicating that a temperature-dependent factor must be present in any modelling approach that aims to describe rubber-like gels.

We have shown that both the stress-strain relation obtained from our model, Eq. 14, and its linearized approximation, Eq. 21, with the elastic modulus given by Eq. 37, describe the experimental data of Figs. 2, 3, and 4, reasonably well. Hence, by assuming that the relationship given by Eq. 45 is also valid for our model, one has that

[TABLE]

which can be replaced in the elastic modulus $G(T)$ given by Eq. 23, so that, in principle, the temperature-dependent correction factor can be explicitly written as

[TABLE]

As one may observe from Fig. 5(e), the above expression describes the experimental data of Ref. yoshikawa2021prx well. For relatively high temperatures, i.e., where $x=\Delta E/k_{B}T\ll 1$ and $e^{x}\simeq 1+x$, Eq. 48 leads to a temperature-dependent correction factor which is given by

[TABLE]

where we used the definition of $T_{0}$ given by Eq. 39. Accordingly, by considering the definition of $w(T)$, Eq. 46, the above approximation leads to the phenomenological relation for $G(T)$ given by Eq. 2, or equivalently, to the elastic modulus $G_{c}(T)$ given by Eq. 43, as expected. Interestingly, if $\Delta E=0$, one finds from both above expressions that $w(T)=w^{*}$ for any temperature $T$, so that Eq. 46 gives $G(T)=w^{*}g_{z}(p)nk_{B}T$, which is precisely the elastic modulus of a purely entropic gel.

As mentioned before, the emergence of the phenomenological factor $w^{*}$ was not explained in Refs. yoshikawa2021prx ; sakumichi2021polymj , so here we suggest that it might be related to a correction on the counting of the number of effective elastic elements in solution. In particular, if one considers that the number density of chains between crosslinks can be defined as $n_{f}=cN_{A}/M_{f}$, and that the molar mass $M_{f}$ of bridged chains between crosslinks is related to the molar mass $M$ of $z$-coordinated precursor molecules as $M_{f}=2(M/z)$, one finds that $n_{f}=(z/2)n$, where $n$ is given by Eq. 40. Hence, if one consistently assumes that the number density of effective elastic elements in solution is, in fact, given by $n_{e}=g_{z}(p)n_{f}\approx w^{*}g_{z}(p)n$, so that $n_{f}\approx w^{*}n$, these relations lead to

[TABLE]

which, in the case of tetra-PEG gels where crosslinks have functionality equal to $z=4$, is indeed close to $2.17$. In practice, by considering that the above equality is correct, it can be used to recast the function $w(T)$ that is given by Eq. 48 together with its general definition (Eq. 46) in order to yield a final expression for the elastic modulus of our model, that is,

[TABLE]

where we recall that $n$ corresponds to the number density of $z$-coordinated precursor molecules (which is also proportional to the number density of crosslinks $n_{c}$), $p$ is the connectivity of the chains in the network (which should be within the range $2/z\ll p\leq 1$ for the function $g_{z}(p)$ given by Eq. 42 to be valid), and $w(T)$ is the temperature-dependent correction factor, which is found to also depend on $z$ and is given by

[TABLE]

Accordingly, for gels with tetra-functional crosslinks of coordination equal to $z=4$, the term inside the left parenthesis in the above expression will ensures that $w(T)=2$ for any value of $T$ if $\Delta E=0$. Even so, it is worth noting that, in general, when $\Delta E=0$, Eq 52 implies that $w(T)=z/2$ for any temperature.

IV Conclusions

In this study we explored a simple coarse-grained model in order to investigate the origin of the negative, $G_{E}(T)$, and positive, $G_{S}(T)$, contributions to the elastic modulus $G(T)$ of rubber-like gels. By computing an exact expression for the free energy, we were able to find the non-linear temperature-dependent stress-strain relation given by Eq. 14, which provided expressions for the positive, $\sigma_{S}(T,\gamma)$, and the negative, $\sigma_{E}(T,\gamma)$, contributions to the stress $\sigma(T,\gamma)$ given by Eqs. 33 and 36, respectively.

The theoretical expressions obtained for the elastic moduli $G(T)$, $G_{S}(T)$, and $G_{E}(T)$, given by Eqs. 23, 28, and 29, respectively, strongly support the idea that the negative contribution to the elastic modulus of rubber-like gels arises from the effective interaction energy between the chains in the gel network and their neighboring solvent molecules, which is characterized in the present model by a positive energy difference $\Delta E$ and a negative value for $\ell_{b}$. Indeed, our theoretical results indicate that $G_{E}(T)<0$ if $\Delta E>0$, $\ell_{s}>0$, and $\ell_{b}<0$, while the condition $\Delta E=0$ implies that there will be no negative (i.e., energy-related) contribution to the elastic modulus, that is, $G_{S}(T)=G(T)$, and thus $G_{E}(T)=G(T)-G_{S}(T)=0$. This last conclusion can be made directly from the equations deduced from our model, i.e., Eqs. 29 and 36, which yields $G_{E}(T)=0$ and $\sigma_{E}(\gamma,T)=0$ when $\Delta E=0$ for any values of temperature $T$ and strain $\gamma$. Also, from Eqs. 24 and 39, one finds that having $T^{*}_{0}=0$ (or, similarly, $T_{0}=0$) is an equivalent way of concluding that no negative energy-related contribution is present in the system, i.e., $\Delta E=0$.

Finally, it is worth mentioning that, by comparing our theoretical results and the results obtained from experiments, we show that, although simple, the present model yields expressions that display a good agreement with the experimental data of Ref. yoshikawa2021prx obtained for tetra-PEG gels. In particular, we validate our approach by considering the experimental data obtained for $\sigma(T,\gamma)$ at different strains (Fig. 2), and its positive ($\sigma_{S}$) and negative ($\sigma_{E}$) contributions as functions of both the strain (Fig. 3) and the temperature (Fig. 4). In addition, we inferred from the experimental results presented in Fig. 5 that the general expression of the elastic modulus should be in the form of Eq. 3, i.e., $G(T)=g_{z}(p)nk_{B}Tw(T)$, with $w(T)$ being a temperature-dependent correction factor, $n$ the number density of $z$-coordinated precursor molecules in solution, and $g_{z}(p)$ a connectivity-dependent function given by Eq. 42. Interestingly, the final expression for $G(T)$ obtained from our model, i.e., Eq. 51, should allow one to describe the mechanical response of both rubber materials ($\Delta E=0$) and rubber-like gels ($\Delta E>0$). Therefore, we believe that the general discussions presented here may provide a useful theoretical framework which could be used in future studies that want to probe and interpret the elastic modulus of rubber-like gels taking into account its negative energy-related contribution, including those based on numerical simulations toda2018aipadv ; sugimura2013polymerj ; rizzi2016jcp .

Acknowledgements. The authors acknowledge the Brazilian agencies FAPEMIG (Process APQ-02783-18), CNPq (N∘ 312999/2021-6 and 426570/2018-9), and CAPES (code 001).**

Appendix A Force-extension relationship in the

constant-strain ensemble

In this Appendix we indicate how the force-extension relation given by Eq. 13 can be also obtained from the free energy $\mathcal{F}(N,T,\ell)$ evaluated in the constant-strain ensemble. First, by considering that $n_{b}=N-n_{s}$, one finds that Eq. 5 yields $\ell=(N-n_{s})\ell_{b}+n_{s}\ell_{s}$, which should be within the range $0\leq\ell<N\ell_{s}$, hence

[TABLE]

The above expressions can be used to obtain both the internal energy $U(N,\ell)$ and the entropy $S(N,\ell)$ that are used to evaluate the free energy $\mathcal{F}(N,T,\ell)=U(N,\ell)-TS(N,\ell)$. Thus, from such free energy one can obtain the mean force $f$ (i.e., the mean stress $\sigma=f/A$) that should be applied to the chain in order to keep its length equal to $\ell$, that is,

[TABLE]

Since the internal energy is given by Eq. 6, one have that the energy-related contribution to the force is simply

[TABLE]

Additionally, the entropy can be estimated as $S(N,\ell)=k_{B}\ln\mathrm{\Omega}(N,\ell)$, with the density of states, Eq. 4, given by

[TABLE]

Here one can consider the Stirling’s approximation, i.e., $\ln(n!)\approx n\ln(n)-n$, which yields the derivative of $S(N,\ell)$ with respect to $\ell$ equals to

[TABLE]

Hence, by replacing the expressions 55 and 57 in Eq. 54, one can verify that the result for the force $f$ is equivalent to the expression 13 derived from the constant-stress ensemble.

Appendix B Stress-strain relationship from the Helmholtz free energy density

We include this Appendix in order to indicate that the stress-strain relationship given by expression 14 can be also evaluated as

[TABLE]

where $\mathcal{F}/V_{0}$ is the Helmholtz free energy density of the system, with $\mathcal{F}=\mathcal{G}+f\ell$ and $V_{0}=A\ell_{0}$. In particular, one can consider that $\mathcal{G}$ and $f$ are given, respectively, by Eqs. 9 and 13, replace $\ell$ by $\gamma$ through Eq. 12, and then use Eq. 11 for $\ell_{0}$, which yields

[TABLE]

where $\Delta\tilde{E}=(\ell_{s}E_{b}-\ell_{b}E_{s})/\Delta\ell$. It is worth emphasizing that, as we have mentioned in Sec. III.1, the response of the whole rubber-like material is obtained by analyzing the properties of its fundamental building block, i.e., a single bridged chain immersed in a incompressible fluid solvent. Arguably, since there is not a preferable orientation for the chains and the material is expected to be disordered and isotropic at mesoscopic length scales, we assume that this approach leads to the same result as if we had analyzed the whole system and the stress depends only on a scalar Helmholtz free energy density which is given by Eq. 59.

Appendix C Definitions of $f_{S}$ and $f_{E}$

In this Appendix we explicitly discuss the alternative definitions of the $S$ and $E$ contributions to the stress that we and the authors of Ref. yoshikawa2021prx explored, since it seems to differ from the usual definitions adopted in the literature flory1953book . In general, for isotropic rubber-like materials, one usually considers the following Helmholtz free energy

[TABLE]

which differential form is given by

[TABLE]

where $V$, $T$, and $\ell$ are the volume, the temperature, and the axial length of the material, while $P$, $S$, and $f$ are the hydrostatic pressure, the entropy, and the axial load, respectively. Hence, from Eq. 61, one can write the equation of state of such material as

[TABLE]

so that Eq. 60 leads to

[TABLE]

Usually, in order to determine the change in the entropy in the second term at the right side of the above equation, one could think of using the Maxwell relation yoshikawa2021prx $-(\partial S/\partial T)_{T,V}=(\partial f/\partial T)_{T,\ell}$. As it happens, the experimental setup to measure the derivative $(\partial f/\partial T)_{V,\ell}$ is not very practical, so to overcome this issue one is led to the following approximation flory1953book

[TABLE]

which turns Eq. 63 into

[TABLE]

where $\lambda=\ell/\ell_{0}=1+\gamma$ is the elongation of the material, and $\gamma$ is the strain defined as in Eq. 12. It is worth metioning that, usually, the definitions of the energetic and the entropic contributions to the force are given, respectively, by the first and the second terms in the above expression flory1953book .

Interestingly, although the experiments in Ref. yoshikawa2021prx were carried out in a fixed strain setup and at a constant pressure, the subscripts $\gamma$ and $P$ were suppressed in the notation of their thermodynamic derivatives. Even so, since we have decided to put forward a modelling approach aiming to describe their experimental data, we avoid approximations by considering that, since $\gamma$ is defined through Eq. 12, so $\lambda$ must also depend on the temperature $T$. In that scenario the differential form of $\ell$ should be written as

[TABLE]

Thus, through Eq. 61, and recalling the fact that the PEG hydrogel can be considered as an incompressible fluid yoshikawa2021prx , i.e., $dV=0$, one may rewrite the differential form of the free energy as

[TABLE]

so that

[TABLE]

Now, from Eq. 62, one has that

[TABLE]

thus, through the identity between second derivatives and by considering the result given by Eq. 68, one finds that

[TABLE]

Accordingly, in this work we assume that the positive ($S$) contribution to the force is defined as

[TABLE]

with the derivative of $f$ given by Eq. 70, that is,

[TABLE]

Furthermore, with that definition, the energy-related ($E$) contribution to the force can be evaluated as $f_{E}(T,\ell)=f(T,\ell)-f_{S}(T,\ell)$, which yields

[TABLE]

Hence, in contrast to the usual definitions that have been adopted in the literature flory1953book ; yasuda2009macromo , one have that $f_{S}\neq-T(\partial S/\partial\ell)_{T,V}$ and $f_{E}\neq(\partial U/\partial\ell)_{T,V}$. Even though neither Eq. 72 nor Eq. 73 correspond to the usual expressions for the energetic and entropic contributions, respectively, the sum of the two expressions lead to the very same equation of state, i.e., Eq. 63. Also, from the above results for $f_{S}(T,\ell)$ and $f_{E}(T,\ell)$, and by considering $\ell=\lambda\ell_{0}=(1+\gamma)\ell_{0}$ (see Eq. 12), one have that the expressions for $\sigma_{S}(T,\gamma)$ and $\sigma_{E}(T,\gamma)$ can be evaluated as

[TABLE]

and

[TABLE]

It can be easily verified that these definitions are not only fully consistent with the expressions derived in Sec. III.5, but also that they are also consistent with the way that the experimental results were obtained in Ref. yoshikawa2021prx .