text:
# GUTPA/00/11/01WUB 01-01 Quark mass effects on the topological susceptibility in QCD

## I Introduction

One of the most intriguing features of QCD is its topological vacuum structure which results in phenomenological consequences with an important bearing on particle physics such as the breaking of the axial $`U_A(1)`$ and chiral symmetries.

In the context of pure gauge theories, lattice methods have by now matured to provide a valuable tool kit for accessing the continuum topological susceptibility, $`\chi `$. Various gluonic definitions of the topological charge lead to consistent estimates of $`\chi `$ , and remnants of the Atiyah-Singer index theorem are witnessed to hold on the lattice: the topological charge, $`Q`$, as obtained after cooling from the gluonic degrees of freedom is consistent with the fermionic index, as determined from counting the zero-level crossings of the eigenvalues of the Wilson-Dirac operator, $`\gamma _5(D/+m_0)`$, under variations of the bare quark mass, $`m_00`$ . Moreover, one might perceive the very value of $`\chi `$ as delivered by quenched lattice simulations , $`\chi _q^{1/4}=213(14)`$ MeV, to confirm the lattice approach to topological properties: it turns out to be in fairly good agreement with the large $`N_c`$ anticipation made by Witten and Veneziano , $`\chi _qf_\pi ^2(m_\eta ^{}^2+m_\eta ^22m_K^2)/6(180\text{MeV})^4`$. In conclusion, the lattice machinery appears to work for the study of topological aspects in gluodynamics.

Unfortunately, however, the situation regarding the QCD vacuum proper is by far less settled: for sheer cost reasons so far only the regime of intermediate sea quark masses, not much lighter than the strange quark, has been explored. In fact it is still under debate whether the data from full QCD simulations support the expected vanishing susceptibility at chiral sea quark masses, $`m0`$:

$$\chi =\frac{m}{n_f}\mathrm{\Sigma }+𝒪(m^2),\mathrm{\Sigma }=\underset{m0}{lim}\underset{V\mathrm{}}{lim}\overline{\psi }\psi .$$
(1)
Attempts to verify this prediction via lattice simulations have a long history. While the first lattice studies in this direction with Kogut-Susskind (KS) sea quarks  could not reach conclusive results, later investigations did yield some qualitative evidence in favour of the expected decrease of $`\chi (m)`$ .

Quite recently, this problem has been revisited from three different sides — but the debate is still open: the CP-PACS  and UKQCD  collaborations employed improved Wilson fermionic and different gluonic actions while the Pisa group  operated with two and four flavours of KS fermions. A further data point for $`n_f=2`$ KS flavours has been obtained by Hasenfratz  (for a recent review see e.g. Ref. ). While both the Pisa group (working at an inverse lattice spacing $`a^12`$ GeV and bare quark masses $`am_00.01`$) and the CP-PACS collaboration (at $`a^11.3`$ GeV, $`m_\pi /m_\rho >0.59`$) saw no evidence whatsoever in favour of the expected chiral behaviour, the UKQCD collaboration did verify (at $`a^12`$ GeV and $`m_\pi /m_\rho >0.57`$) a decrease of $`\chi (m)`$, consistent with theoretical expectations.

In this paper we shall present an analysis focused on this issue, based on the final statistics of our SESAM and T$`\chi `$L samples of QCD vacuum configurations. The article is organised as follows. In Sect. II we describe the details of our simulation and methodology, including a comparison between fermionic and gluonic definitions of the topological charge. Preliminary results on this comparison, based on smaller statistical samples, have been reported by us previously . In Sect. III we investigate correlations between the topological charge and the hadron spectrum. Finally, in Sect. IV, we present the lattice data on the topological susceptibility.

## II Determination of the topological charge

### A Measurements

We analyse ensembles of gauge configurations that have been generated by means of the hybrid Monte Carlo (HMC) algorithm using the Wilson fermionic and gluonic actions with $`n_f=2`$ mass degenerate quark flavours at the inverse lattice coupling, $`\beta =5.6`$, corresponding to an inverse lattice spacing $`a^1=(2.65_8^{+5}\pm 0.14)`$ GeV at physical sea quark masses. This was done on $`L_\sigma ^3\times L_\tau =16^3\times 32`$ as well as on $`24^3\times 40`$ lattices at five different values of the sea quark mass parameter. The corresponding chiralities can be quantified in terms of $`m_\pi /m_\rho `$-ratios, ranging between $`0.834(3)`$ and $`0.574(13)`$ . The relevant simulation settings are displayed in Table I . At each $`\kappa `$-value 4,000–5,000 thermalized HMC trajectories have been generated. In addition to the dynamical quark simulations, quenched reference measurements on $`16^4`$ lattices at $`\beta =6.0`$ were performed. The configurations that enter the analysis are separated by 25 (24 at $`\kappa =0.1565`$ and 10 in the quenched simulation) HMC trajectories.

Lower limits on the integrated autocorrelation times, $`\tau _{\text{int}}`$, are estimated by binning the data for the topological charge $`Q_i`$, $`i=1,\mathrm{},n_{\text{meas}}`$, into $`n_{\text{bin}}`$ blocks that contain $`m`$ successive measurements each. On each such block of length $`m`$ an average $`\overline{Q}_j^m,j=1,\mathrm{},n_{\text{bin}}`$, is calculated: $`\overline{Q}_j^m=(1/m)_{i=1}^mQ_{m(j1)+i}`$. We determine fluctuations between these bins,

$$\mathrm{\Delta }Q_m^2=\frac{1}{n_{\text{bin}}(n_{\text{bin}}1)}\underset{j=1}{\overset{n_{\text{bin}}}{}}\left(\overline{Q}_j^mQ\right)^2,$$
(2)
where the average $`Q0`$ is calculated on the first $`n_{\text{bin}}\times m`$ configurations. We then estimate the autocorrelation times,

$$\tau _{\text{int}}\frac{1}{2\nu _{\text{meas}}}\underset{m^2n_{\text{meas}}/2}{\mathrm{max}}\left(\frac{\mathrm{\Delta }Q_m^2}{\mathrm{\Delta }Q_1^2}\right),$$
(3)
that are included in the table. $`\nu _{\text{meas}}`$ denotes the measurement frequency. On the $`16^3\times 32`$ volume at $`\kappa =0.1575`$, in addition to the $`\nu _{\text{meas}}=1/25`$ time series, we also determined the topological charge with increased frequency $`\nu _{\text{meas}}=1/2`$, with a reduced number of 10 (as opposed to 60) cooling sweeps. The larger frequency enabled us to compute the autocorrelation time from the autocorrelation function itself, with full control over statistical errors . The resulting value, $`\tau _{\text{int}}=54(4)`$, is consistent with our estimate obtained in the way described above, $`\tau _{\text{int}}55`$. We take this as an indication that our estimates are reasonable. Interestingly, the autocorrelation times grow both with decreasing quark mass and increasing volume as already observed in Ref. .

We employ a gluonic as well as a fermionic definition of the topological charge. In the gluonic case we define a topological charge density,

$$q(x)=\frac{1}{16\pi ^2}\text{tr}F_{\mu \nu }(x)\stackrel{~}{F}_{\mu \nu }(x)=\frac{1}{16\pi ^2}F(x)F(x),$$
(4)
where we use the symmetric (clover leaf) definition ,

$$F_{\mu \nu }=\frac{i}{4a^2}\left(U_{\mu \nu }+U_{\nu ,\mu }+U_{\mu ,\nu }+U_{\nu ,\mu }4\right),$$
(5)
and $`\stackrel{~}{F}_{\mu \nu }=\frac{1}{2}ϵ_{\mu \nu \rho \sigma }F_{\rho \sigma }`$, which is correct up to order $`a^2`$ lattice artefacts. The topological charge,

$$Q=\underset{x}{}q(x)=\frac{1}{16\pi ^2}(F,F),$$
(6)
should then approach integer values as $`a0`$ in the infinite volume limit, on sufficiently smooth gauge configurations.

### B Renormalization

In a quantum field theory, both $`Q`$ and $`\chi `$ undergo multiplicative renormalization. Moreover, $`\chi `$ requires additive renormalization. Cooling  is meant to filter for (semi)-classical features; indeed, after cooling, both renormalization constants have been shown to be close to their trivial values . In the context of this article we will not investigate properties of the density distribution $`q(x)`$ itself but rather stick to the net topological charge $`Q`$ only. For this purpose the iterative application of cooling by simple minimisation of the Wilson plaquette action is appropriate. To keep the cooling update local we visit the lattice sites in an even-odd pattern, rather than in the sequential ordering that is usually employed in the Monte Carlo updating of pure gauge configurations. The inner-most loop within a cooling sweep runs across the directions $`\mu `$, and we carry out 60 such sweeps. In Fig. 1 we illustrate the resulting numerical pattern of topological charges which cluster nicely around integer values as anticipated.

The fermionic method consists of determining the topological charge à la Smit and Vink ,

$$Q=\nu _{}\nu _+=Z_Pm_0\text{tr}\gamma _5M^1,$$
(7)
that is inspired by the continuum axial divergence relation, $`_\mu [\overline{\psi }(x)\gamma _\mu \gamma _5\psi (x)]=2m\overline{\psi }(x)\gamma _5\psi (x)2n_fq(x)`$. The above trace is to be taken over Dirac, colour and space-time indices only (not over the flavours), i.e. $`\overline{\psi }\gamma _5\psi =n_f\text{tr}\gamma _5M^1`$ for $`n_f`$ mass degenerate quark flavours. $`Z_P=𝒪(1)`$ is a renormalization constant, $`m_0=(\kappa ^1\kappa _c^1)/(2a)`$ denotes the bare quark mass and $`M`$ is the lattice discretised version of $`D/+m_0`$,

$`2\kappa M_{xy}=\delta _{xy}`$ $``$ $`\kappa {\displaystyle \underset{\mu }{}}[(1\gamma _\mu )U_{x,\mu }\delta _{x+\widehat{\mu },y}`$ (8)
$`+`$ $`(1+\gamma _\mu )U_{x\widehat{\mu },\mu }^{}\delta _{x\widehat{\mu },y}].`$ (9)
We determine $`\text{tr}\gamma _5M^1`$ using $`Z_2`$ noisy sources with diagonal improvement as detailed in Ref. . On the $`\kappa =0.158`$ configurations 100 such estimates were performed while at all other $`\kappa `$ values we averaged over 400 estimates*** Approximating $`Q`$ by a finite number of noise vectors can result in an underestimated integrated autocorrelation time. Nonetheless, we find the $`\tau _{\text{int}}`$-estimates from the fermionic definition  to be consistent with the gluonic ones of Table I..

The renormalization constant $`Z_P`$ in Eq. (7) is unknown. We attempt to estimate the combination $`Z_Pm_0a`$ from the ratio of the gluonic \[Eq. (6)\] and of the fermionic \[Eq. (7))\] definitions. We determine $`m_0a`$ using the critical value, $`\kappa _c=0.15849(2)`$ . The results are displayed in Table II, assuming the multiplicative renormalization of the gluonic definition after cooling to be unity. The estimates of the fermionic traces are subject to statistical uncertainties, $`\mathrm{\Delta }Q/Z_P1`$ on the $`16^3\times 32`$ lattices and $`\mathrm{\Delta }Q/Z_P2`$ – 3 on the $`24^3\times 40`$ volumes. In addition, both definitions are expected to suffer from different $`𝒪(a)`$ lattice artefacts. We do not attempt to estimate the resulting statistical and systematic uncertainties on $`Z_P`$. We find $`Z_P`$ to be of order one and to depend monotonously on the quark mass, with the deviation from unity increasing when approaching the chiral limit.

The correlation between the two definitions is visualised in the scatter plot, Fig. 2, for the $`16^3\times 32`$ lattices at $`\kappa =0.1575`$. While the gluonic data points cluster around integer values, this is not the case for the (stochastically estimated) fermionic values. The data of the figure are normalised such that the points should collapse onto the line with slope one. Indeed, 137 out of 199 values (62 %) lie within the one $`\sigma `$ error band. The correlation between the two definitions is also visible from the Monte Carlo histories depicted in Fig. 3.

### C Ergodicity

In Fig. 3 we survey the Monte Carlo histories of $`Q`$ for all our runs. These charts provide a gross evidence for the quality of our data in view of the decorrelation with respect to the topological sectors. The topological susceptibility,

$$\chi =\underset{x}{}q(x)q(0)=\frac{Q^2}{V},$$
(10)
should be independent of the volume $`V=L_\sigma ^3\times L_\tau a^4`$ to a first approximation. Therefore, the modulus of the topological charge $`|Q|`$ should scale in proportion to $`\sqrt{V}`$. Indeed, the topological charge distribution on the large lattice at $`\kappa =0.1575`$ is by a factor of about two wider than that on the small lattice. We also observe reduced fluctuations as we increase $`\kappa `$ at fixed $`\beta `$. While the $`16^3\times 32`$ time histories appear to tunnel ergodically through all topological sectors the total number of tunnellings observed for $`\kappa =0.1580`$ is not yet sufficient to achieve a symmetric distribution.

## III Topology and the hadron spectrum

We address the question whether the statistics presented in Fig. 3 suffices to expose a significant $`|Q|`$ dependency of hadronic states. To minimise statistical errors we subdivided each sample into two subsamples only, one containing configurations with $`|Q|1.5`$ and one with $`|Q|>1.5`$. For the smaller volume at $`\kappa =0.1575`$ this division results in two subsamples of approximately 100 configurations each.

The most interesting test case is given by the effective masses in the flavour singlet pseudoscalar channel, $`m_\eta ^{}^{\text{eff}}(t)`$, that should be particularly sensitive to the vacuum topology. This is confirmed by the data displayed in Fig. 4. We find the $`\eta ^{}`$ mass on vacuum configurations with $`|Q|1.5`$ to systematically lie below the corresponding mass obtained with the cut $`|Q|>1.5`$. Note that in the infinite volume limit one would expect such a sensitivity of correlation functions on $`|Q|`$ to disappear. For further details we refer to Ref. .

In Fig. 5 we show the corresponding flavour octet channel quantities, i.e. $`\pi `$ effective masses, where we anticipate no such correlation with topology. This is indeed borne out by the data: not only at $`t\mathrm{}`$ but time slice by time slice no sensitivity to the value of the modulus of the topological charge is detected, with rather high statistical accuracy. The horizontal line with error band indicates the asymptotic large $`t`$ results obtained from a fit to the entire data sample .

Apart from the $`\eta ^{}`$ none of the standard mesonic, baryonic and glueball-like states exhibit correlations between $`|Q|`$ and the respective effective masses. In Figs. 68 we illustrate this for the $`\rho `$, the nucleon and the $`\mathrm{\Delta }`$. In no case have we found any systematic effects on the static potential. Of course this does not exclude the possibility of correlations between hadronic properties and the distribution of instantons and anti-instantons and their relative orientations in colour and position space .

## IV The topological susceptibility

In Figs. 911 we display three histograms of topological charge distributions. Each bin with width $`\delta Q`$, centred around $`Q_j=j\delta Q`$, contains all measurements resulting in charges within the interval $`[Q_j\delta Q/2,Q_j+\delta Q/2]`$, where $`\delta Q=1`$ on the $`16^3\times 32`$ lattice depicted in Fig. 9 and $`\delta Q=3`$ for the $`24^3\times 40`$ lattices of Fig. 10 and Fig. 11. In addition to the data we display Gaussian distributions,

$$n(Q)=\frac{n_{\text{meas}}\delta Q}{\sqrt{2\pi Q^2}}\mathrm{exp}\left(\frac{Q^2}{2Q^2}\right).$$
(11)
In Fig. 9 we include the statistical uncertainties of the individual bins, $`\mathrm{\Delta }n=2\tau _{\text{int}}\nu _{\text{meas}}\sqrt{n}`$, while in Fig. 10 and Fig. 11 the error on the width of the distribution is reflected by the error band around the central curve.

The fitted topological susceptibilities $`\chi =Q^2/V`$ as well as $`Q`$ are displayed in Table III. For orientation, we convert $`\chi `$ into physical units in the last column of the table, using $`r_0^1=(394\pm 20)`$ MeV. The distribution at $`\kappa =0.158`$ is not symmetric around zero anymore, as reflected by the value $`Q=1.9\pm 1.0`$ and by Fig. 11. Therefore, in this case, the resulting value of $`\chi `$ should pass with a grain of doubt, notwithstanding the comfortable error bars. A comparison between the $`L_\sigma =16`$ and $`L_\sigma =24`$ results at $`\kappa =0.1575`$ reveals that the level of finite size effects on $`\chi `$ is below the statistical errors.

The Gell-Mann-Oakes-Renner (GMOR) relation connects the pion mass to the chiral quark mass, $`m`$, via the pion decay constantOur normalisation conventions are such that $`f_\pi 92`$ MeV. $`f_\pi `$ in the chiral limit and the chiral condensate, $`\mathrm{\Sigma }`$:

$$f_\pi ^2m_\pi ^2=2m\mathrm{\Sigma }+𝒪(m^2).$$
(12)
A direct lattice determination of $`f_\pi `$ from the appropriate current matrix element yields $`\sqrt{2}f_\pi =Z_A0|A_4|\pi /m_\pi 0.0496(34)a^1`$, in the limit $`m_\pi ^20`$ . Allowing for an additional systematic error of 10 % on $`f_\pi `$ to account for the fact that we have only determined the axial vector renormalization constants $`Z_A`$ perturbatively, we arrive at the lattice estimate

$$f_\pi =(0.238\pm 0.030)r_0^1=(94\pm 13)\text{MeV}$$
(13)
at $`\kappa =\kappa _c`$, using the value $`r_0^1a=0.1476_{30}^{+47}`$ . In combining Eq. (12) with Eq. (1) one expectsThis formula as well as the GMOR relation Eq. (12) only apply to $`n_f2`$ since for $`n_f=1`$ no light pion exists (for a finite number of colours, $`N_c<\mathrm{}`$).

$$\chi =\frac{f_\pi ^2m_\pi ^2}{2n_f}+𝒪(m_\pi ^4)$$
(14)
for small pseudoscalar masses and large Leutwyler-Smilga parameters, $`x=Vm\mathrm{\Sigma }1`$ .

In Table IV we display the spatial lattice extents, $`L_\sigma am_\pi `$ as well as $`x`$ and $`Q^2=\chi V`$. We estimate the (scheme independent) combination $`m\mathrm{\Sigma }`$ that appears within $`x`$ from Eqs. (12) – (13). The topological charge fluctuations $`Q^2`$ correspond to the average number of non-interacting instantons that can be accomodated within the simulated lattice volume . It has been argued  that $`x>10`$ is already a large volume in the Leutwyler-Smilga  sense, such that Eq. (14) should safely apply to all but the $`16^3\times 32`$ lattices at $`\kappa =0.157`$ and $`\kappa =0.1575`$. Even at the latter, more critical, $`\kappa `$-value, where we have two lattice volumes that correspond to $`x6`$ and $`x25`$ at our disposal, we do not find any volume dependence of $`\chi `$.

In the regime of large $`m_\pi `$ one would expect $`\chi `$ in units of some reference scale like $`r_0^1`$ to approach the quenched value as a smooth function of $`m_\pi `$. Prior to the comparison with Eq. (14) one should in principle extrapolate the lattice results on $`\chi `$ to the continuum limit, for instance along a trajectory of constant $`m_\pi r_0`$. Having only $`\beta =5.6`$ data at our disposal we cannot yet perform this extrapolation but conjecture that our results are already sufficiently close to the continuum limit for Eq. (14) to hold. This assumption is plausible in view of the fact that the topological charges obtained from the gluonic and fermionic definitions agree with each other reasonably well (with a renormalization constant $`Z_P`$ of order one) and of our previous studies of light hadronic quantities .

For the purpose of the present investigation we shall consider the dimensionless topological susceptibility, $`\chi r_0^4`$, as a function of the dimensionless combination $`m_\pi ^2r_0^2`$ in order to carry out the chiral extrapolation. The results are plotted in Fig. 12, together with the region that is allowed for by the leading order expectation, Eq. (14), for the $`f_\pi `$ value of Eq. (13). All unquenched data except for the point at the heaviest quark mass ($`m_\pi ^2r_0^25`$) are consistent with this expectation. Note that the most chiral data point stems from the run with $`\kappa =0.158`$ which, as can be seen from Figs. 3 and 11, is not of sufficient statistical quality for the purpose of topological studies. The right-most entry in Fig. 12 is the result of our quenched reference study.

A linear fit, excluding the largest mass point, to the parametrisation (14) renders $`f_\pi =0.255(11)r_0^1`$ with $`\chi ^2/N_{DF}=3.93/4`$. This value compares reasonably well with the expectation, Eq. (13). On the other hand our data are consistent with a mass independence of $`\chi `$ too: fitting them to a constant yields $`\chi r_0^4=0.0503(23)`$ with $`\chi ^2/N_{DF}=2.67/5`$, in agreement with the quenched reference point, $`\chi r_0^4=0.065(11)`$.

Our present study is based on data obtained at a fixed value of $`\beta `$ and the quark mass is varied by just tuning $`\kappa `$. It is worthwhile to relate our findings to those recently presented by the UKQCD collaboration§§§Unfortunately, none of the other collaborations have converted their results into units of $`r_0`$.  who simultaneously vary $`\beta `$ and $`\kappa `$ to change the quark mass while keeping the lattice spacing $`a`$ fixed in units of $`r_0`$. Apart from this they dispose of ensembles of gauge field configuration of comparable statistical sizes, using slightly coarser lattice spacings with an improved fermionic action.

In order to facilitate this comparison we have superimposed both data setsThe statistical errors quoted by the UKQCD collaboration are somewhat smaller than ours since no signs of autocorrelation effects have been found in their study . in Fig. 13. The data do not clearly contradict each other, however, we are unable to confirm the decrease of the topological susceptibility that has been reported in Ref. , at pion masses within the range of present-day lattice simulations.

## V Conclusions

We have demonstrated that Wilson fermions are suitable for lattice studies of topology related effects in the QCD vacuum. We have found agreement between the topological charge computed after cooling from the field strength tensor, $`Q=(F,F)/(16\pi ^2)`$, and the fermionic definition given by Smit and Vink with a renormalization constant $`Z_P`$ of order one. The trace has been estimated from diagonal improved $`Z_2`$ noisy sources. As the next step we plan to investigate correlations between the topological charge density distribution and eigenvectors of $`\gamma _5M`$.

The topological susceptibility at $`\beta =5.6`$ has been found to be consistent with that of quenched studies at large sea quark masses. Our data at smaller mass values is consistent with the asymptotic slope in $`m_\pi ^2`$ that is expected from our independently determined $`n_f=2`$ pion decay constant. However, we are unable to unambiguously confirm the decrease that has recently been reported by the UKQCD collaboration  with a different fermionic action. Results obtained by the CP-PACS collaboration  who cover a similar range of quark masses on somewhat coarser lattices do not clearly show this tendency either. We plan to clarify this issue in simulations at different $`\beta `$ values which will enable us to perform a continuum limit extrapolation, and at smaller sea quark masses.

###### Acknowledgements.

We thank S. Güsken and P. Überholz for their contributions at an earlier stage of this project. We gratefully acknowledge R. Burkhalter for detecting an important mistake in this article. Our European collaboration was funded by the EU network “Hadron Phenomenology from Lattice QCD” (HPRN-CT-2000-00145). G.B. received support from EU grant HPMF-CT-1999-00353 and PPARC grant PPA/G/O/1998/00559, N.E. from DFG grant Li 701/3-1. B.O. appreciates support from the DFG Graduiertenkolleg “Feldtheoretische und Numerische Methoden in der Statistischen und Elementarteilchenphysik”. The HMC productions were run on the APE100 systems at INFN Roma and NIC Zeuthen. We are grateful to our colleagues G. Martinelli and F. Rapuano for the fruitful T$`\chi `$L-collaboration. Part of the physics evaluation has been carried out on the ALiCE cluster computer at Wuppertal University.


format: text

text:
# Stress–softening and recovery of elastomers

## I Introduction

This study is concerned with modeling stress–softening and subsequent recovery (by annealing at an elevated temperature and by swelling) of elastomers. Stress–strain relations for unfilled and particle reinforced rubbery polymers have been the focus of attention in the past four decades. This may be explained by numerous applications of rubbery-like materials in industry (vehicle tires, seals, shock absorbers, flexible joints, etc.), as well as by some peculiarities in their mechanical behavior whose physical mechanisms are not quite clear. A renewal of the interest to the mechanical response of elastomers observed in the past decade, see Refs. , is associated with the development of constitutive equations that (i) account for the molecular structure of polymers at the micro-level and (ii) describe some features of their response induced by time-dependent changes in the internal structure (viscoplasticity, damage, micro-fracture, etc.).

Among these peculiarities, one of the most important for applications is the Mullins effect . This phenomenon is evidenced in tensile , compressive  and shear  tests as a noticeable difference between the stress–strain curve for a virgin specimen and that for the material reloaded after retraction. A conventional standpoint is that cyclic pre-stretching with 3 to 6 cycles and the maximum elongation ratio of about $`k=2`$ makes the response of an elastomer repeatable . However, the number of cycles, as well as their amplitude and frequency are chosen by the trial-and-error method and substantially depend on what the repeatability of a stress–strain diagram means for the experimentalist.

Although the physical mechanism of the Mullins effect in unfilled elastomers remains hitherto obscure, it is traditionally associated with material damage , non-affinity of deformation for a polymeric network , viscous drag excerted on the network chains by their environment  and mechanically-induced crystallization of rubbery polymers . These reasons may, however, be questioned, because observations reveal that with an increase in strains above the maximal level reached at pre-stretching, the stress–strain curve returns to that for a virgin sample.

The present study aims to explain stress-softening of elastomers under cyclic loading by an increase in the average size of globules formed by polymeric chains (uncoiling of macromolecules). For this purpose, we develop a phenomenological model for the stress–strain response of rubbery polymers and determine its parameters by fitting experimental data in uniaxial tensile tests for virgin, pre-stretched and recovered specimens.

An elastomer is taken as a network of linear macromolecules bridged by permanent (chemical crosslinks) and temporary (physical crosslinks, entanglements and van der Waals forces) junctions. For the sake of simplicity, it is presumed that junctions move affinely with the bulk material. The novelty of our approach is that the number of junctions is not fixed, but assumed to alter under loading. To explain this phenomenon, we ascribe two possible states to an entanglement: loose (passive) and tight (active). An entanglement in the loose state does not create a junction, but imposes only topological constrains on the number of available configurations for a chain. On the contrary, an entanglement in the tight state is thought of as the same junction as that formed by a chemical crosslink. The transition of an entanglement from its passive state to the active state (tightening of entanglements) is supposed to be driven by mechanical factors. Introducing a simple kinetic equation for the rate of this transition, we study the response of a temporary network with a time-varying number of junctions.

A chain that connects two neighboring junctions is modeled as a sequence of strands (statistically independent segments) bridged by bonds . With reference to the Kratky–Porod concept , strands are taken as inextensible (rigid) rods, whereas the mechanical energy of a chain is determined as the sum of strain energies for bonds. It is postulated that a bond may be in two stable conformations: flexed (trans) and extended (cis). A bond in the flexed conformation is thought of as a linear elastic solid, whereas the mechanical energy of a bond in the extended conformation vanishes. Constitutive equations for an elastomer and the kinetic equations for the transition of bonds from their flexed to extended conformations are derived by using the laws of thermodynamics. The development of stress–strain relations is based on the hypothesis that the characteristic time for changes in the concentration of active entanglements substantially exceeds that for mechanical deformation.

The paper is organized as follows. Deformation of a chain is discussed in Section 2. The mechanical energy of a temporary network is determined in Section 3. Stress–strain relations are developed in Section 4 using the laws of thermodynamics. Kinetic equations for transition of entanglements from their loose to tight state are introduced in Section 5. Uniaxial extension of a specimen is studied in Section 6. Results of numerical simulation are compared with experimental data for natural rubber vulcanizates in Section 7. Section 8 deals with the effects of cyclic stretching and recovery on experimental constants. Some concluding remarks are formulated in Section 9.

## II Deformation of a long chain

An elastomer is modeled as an ensemble of linear macromolecules bridged by crosslinks, entanglements and van der Waals forces. Active links between macromolecules are thought of as junctions (permanent in the case of chemical crosslinks and temporary for entanglements and van der Waals forces). A sequence of mers (belonging to a polymeric molecule) between two subsequent junctions is associated with a chain. Chains are divided into strands (statistically independent segments) bridged by bonds. In accord with the Kratky–Porod model, a chain is treated as an aggregate consisting of $`N+1`$ identical inextensible strands linked in sequel. The average number of bonds in a chain change under stretching of a specimen because of slippage of entanglements with respect to chains and transition of entanglements from their loose to the tight state. This implies that $`N`$ is taken as a function of time, $`N=N(t)`$.

We adopt a zigzag model  and assume that a bond bridging two neighboring strands is characterized by one of the two stable conformations: flexed and extended. In the stress-free state the angle between two strands linked by a bond in the flexed conformation equals $`\theta (0,\pi )`$ (the same for all bonds) and that for a bond in the extended conformation equals $`\pi `$ (see a sketch depicted in Figure 1 of Ref. ). In thermal equilibrium before deformation all bonds are in the flexed conformation. Mechanical loading activates chains, which implies that conformations of some bonds alter. The numbers of bonds in flexed and extended conformations, $`N_\mathrm{f}(t)`$ and $`N_\mathrm{e}(t)`$, obey the balance law

$$N_\mathrm{f}(t)+N_\mathrm{e}(t)=N(t).$$

Introducing the average concentration of bonds in the extended conformation, $`n=n(t)`$, we obtain

$$N_\mathrm{f}(t)=N(t)[1n(t)],N_\mathrm{e}(t)=N(t)n(t).$$
(1)
In a deformed state the angle between strands linked by a bond in the flexed conformation changes, whereas for a bond in the extended conformation (modeled as two rigid rods directed along a straight line) this angle remains unaltered. The strain, $`e`$, from the stress-free state of a chain to its deformed state equals the sum of strains for bonds in the flexed conformation,

$$e=N_\mathrm{f}e_\mathrm{f}.$$

This equality together with Eq. (1) implies that

$$e_\mathrm{f}=\frac{e}{N(1n)}.$$
(2)
Bonds in the flexed conformation are modeled as linear elastic solids with the mechanical energy

$$\frac{1}{2}\mu e_\mathrm{f}^2,$$

whereas the strain energy of bonds in the extended conformation vanishes. The mechanical energy of a chain, $`w`$, equals the sum of the mechanical energies for individual bonds. It follows from Eqs. (1) and (2) that

$$w(e)=\frac{\mu e^2}{2N(1n)}.$$
(3)
To express the strain in a chain, $`e`$, in terms of the macro-strain tensor for a network, $`\widehat{ϵ}`$, we consider a chain connecting neighboring junctions $`A_1(t)`$ and $`A_2(t)`$. Let $`\overline{r}_{01}(t)`$ and $`\overline{r}_{02}(t)`$ be radius vectors of these points in the stress–free state and $`\overline{r}_1(t)`$ and $`\overline{r}_2(t)`$ their radius vectors in the deformed state at time $`t0`$. The end-to-end length of the chain reads $`\delta _0(t)`$ in the stress-free state and $`\delta (t)`$ in the deformed state. Introducing the guiding vector in the stress-free state, $`\overline{l}`$, (the unit vector directed along the end-to-end vector for the chain), we obtain

$$\overline{r}_{02}(t)\overline{r}_{01}(t)=\delta _0(t)\overline{l}.$$

In the deformed state, the junctions occupy points with the radius vectors

$$\overline{r}_1(t)=\overline{r}_{01}(t)+\overline{u}(t,\overline{r}_{01}(t)),\overline{r}_2(t)=\overline{r}_{02}(t)+\overline{u}(t,\overline{r}_{02}(t)),$$

where $`\overline{u}(t,\overline{r})`$ is the displacement vector at point $`\overline{r}`$ for transition from the stress-free state to the deformed state at time $`t`$. The end-to-end vector for the chain in the deformed state is given by

$$\overline{R}(t)=\overline{r}_2(t)\overline{r}_1(t)=\delta _0(t)\overline{l}+[\overline{u}(t,\overline{r}_{01}(t)+\delta _0(t)\overline{l})\overline{u}(t,\overline{r}_{01}(t))].$$

Neglecting terms beyond the first order of smallness compared to $`\delta _0(t)`$, we find that

$$\overline{R}(t)=\delta _0(t)\overline{l}[\widehat{I}+\overline{}_0\overline{u}(t)],$$

where $`\overline{}_0`$ is the gradient operator in the stress-free state, $`\widehat{I}`$ is the unit tensor, the dot stands for inner product and the argument $`\overline{r}_{01}`$ is omitted. In terms of the radius vector in the deformed state, $`\overline{r}`$, this equality reads

$$\overline{R}(t)=\delta _0(t)\overline{l}\overline{}_0\overline{r}(t)=\delta _0(t)[\overline{}_0\overline{r}(t)]^{}\overline{l},$$
(4)
where $``$ stands for transpose. The end-to-end length of the chain in the deformed state, $`\delta (t)`$, is given by

$$\delta ^2=\overline{R}\overline{R}.$$

This equality together with Eq. (4) yields

$$\delta ^2(t,\overline{l})=\delta _0^2(t)\overline{l}\widehat{C}(t)\overline{l},$$
(5)
where

$$\widehat{C}(t)=\overline{}_0\overline{r}(t)[\overline{}_0\overline{r}(t)]^{}$$
(6)
is the Cauchy deformation tensor for transition from the stress-free state of the network to its deformed state at time $`t`$.

The extension ratio for the chain, $`\lambda (t,\overline{l})`$, is defined as the ratio of the current end-to-end length, $`\delta (t,\overline{l})`$, to that of the chain in its activated stress-free state, $`\delta ^{}(t,\overline{l})`$. The latter state is defined as a state in which the chain is unloaded, but the numbers of bonds in various conformations coincide with their current values, $`N_\mathrm{f}(t,\overline{l})`$ and $`N_\mathrm{e}(t,\overline{l})`$. It differs from the equilibrium stress-free state of the chain, where the number of bonds in the extended conformation vanishes. The difference between the end-to-end lengths of a chain in the activated state, $`\delta ^{}(t,\overline{l})`$, and in the equilibrium state, $`\delta _0(t)`$, determines the end-to-end elongation driven by transformation of bonds from their flexed conformation to the extended conformation. We assume that the transformation-induced end-to-end elongation is proportional to the number of bonds acquiring the extended conformation,

$$\delta ^{}(t,\overline{l})=\delta _0(t)[1+\eta n(t,\overline{l})],$$
(7)
where $`\eta >0`$ characterizes an increment of the end-to-end length driven by an individual transition. The quantity $`\eta `$ reflects the average size of a globule created by a polymeric chain: when $`\eta `$ is small, the chain is rolled into a tight coil, whereas an increase in $`\eta `$ is tantamount to unfolding the coil. Because the parameter $`\eta `$ describes the current state of a chain, it should be thought of as a function of time (strain). To simplify calculations, we treat it as a constant, but suppose that $`\eta `$ can change from one stress–strain curve to another in cyclic tests.

It follows from Eqs. (5) and (7) that the Hencky strain for a chain, $`e=\mathrm{ln}\lambda `$, reads

$$e(t,\overline{l})=\mathrm{ln}\frac{\delta (t,\overline{l})}{\delta ^{}(t,\overline{l})}=\frac{1}{2}\mathrm{ln}\left[\overline{l}\widehat{C}(t)\overline{l}\right]\mathrm{ln}\left[1+\eta n(t,\overline{l})\right].$$
(8)
Equations (3) and (8) determine the mechanical energy of a chain, $`w`$, in terms of the Cauchy deformation tensor for the network, $`\widehat{C}`$.

## III Strain energy density of a network

We adopt the conventional hypothesis that the excluded-volume effect and other multi-chain effects are screened for an individual chain by surrounding macromolecules . This implies that the energy of interaction between chains is neglected (under the incompressibility condition for the network) and the mechanical energy of the network equals the sum of the strain energies for individual chains. Assuming the distribution of chains to be isotropic, we find the concentration of chains in a network (per unit mass) with guiding vector $`\overline{l}`$,

$$X(t,\overline{l})=\frac{\mathrm{\Xi }(t)}{4\pi }\mathrm{sin}\vartheta d\vartheta d\phi ,$$

where $`\vartheta `$ and $`\phi `$ are Euler’s angles which determine the position of the unit vector $`\overline{l}`$ with respect to some Cartesian coordinate frame and $`\mathrm{\Xi }(t)`$ is the number of chains (per unit mass) at time $`t0`$. To determine the strain energy of a network, we multiply the number of chains, $`X(t,\overline{l})`$, by their mechanical energy, $`w`$, and sum the results for various guiding vectors, $`\overline{l}`$,

$$W(t)=\frac{\mu \mathrm{\Xi }(t)}{8\pi N(t)}_0^{2\pi }𝑑\phi _0^\pi \frac{e^2(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}\mathrm{sin}\vartheta d\vartheta .$$
(9)
The quantities $`N(t)`$ and $`\mathrm{\Xi }(t)`$ are connected by the mass conservation law

$$mN(t)\mathrm{\Xi }(t)=1,$$
(10)
where $`m`$ stands for the average mass of a strand. It follows from Eqs. (9) and (10) that

$$W(t)=\frac{\mu m\mathrm{\Xi }^2(t)}{8\pi }_0^{2\pi }𝑑\phi _0^\pi \frac{e^2(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}\mathrm{sin}\vartheta d\vartheta .$$
(11)
Let $`\mathrm{\Xi }_{}`$ be the average concentration of chains (per unit mass) in a (hypothetical) totally disentangled network and $`M(t)`$ the current number of active entanglements per unit mass. Because any tight entanglement doubles the number of chains involved in the knot, the concentration of chains in the entangled network reads

$$\mathrm{\Xi }(t)=\mathrm{\Xi }_{}+2M(t)=\mathrm{\Xi }_0[1+\nu (t)],$$
(12)
where $`\mathrm{\Xi }_0=\mathrm{\Xi }_{}+2M(0)`$ is the average number of chains (per unit mass) in an equilibrium network (before a test) and the dimensionless function

$$\nu (t)=2\frac{M(t)M(0)}{\mathrm{\Xi }_0}$$
(13)
characterizes the ratio of the number of loose entanglements that have transformed into the active state within the interval $`[0,t]`$ to the initial number of chains. Combining Eqs. (11) and (12), we obtain

$$W(t)=\frac{\mu m\mathrm{\Xi }_0^2}{8\pi }[1+\nu (t)]^2_0^{2\pi }𝑑\phi _0^\pi \frac{e^2(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}\mathrm{sin}\vartheta d\vartheta .$$
(14)
The average mechanical energy per chain is given by

$$w_0(t)=\frac{W(t)}{\mathrm{\Xi }(t)}.$$

Substitution of Eqs. (12) and (14) into this equality yields

$$w_0(t)=\frac{\mu m\mathrm{\Xi }_0}{8\pi }[1+\nu (t)]_0^{2\pi }𝑑\phi _0^\pi \frac{e^2(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}\mathrm{sin}\vartheta d\vartheta .$$
(15)
Our objective now is to calculate the derivative of the function $`W`$ with respect to time. For this purpose, we suppose that the characteristic rate for transition of bonds from their flexed to extended conformation substantially exceeds that for changes in the concentration of active (tight) entanglements. This assumption may be explained by the fact that the characteristic length of a chain (the length between neighboring entanglements) dramatically exceeds that for a strand (the length-scale associated with a bond). This implies that the function $`\nu (t)`$ may be treated as a constant when expression (14) is differentiated with respect to $`t`$. Bearing this hypothesis in mind, we find from Eq. (14) that

$`{\displaystyle \frac{dW}{dt}}(t)`$ $`=`$ $`{\displaystyle \frac{\mu m\mathrm{\Xi }_0^2}{4\pi }}[1+\nu (t)]^2{\displaystyle _0^{2\pi }}𝑑\phi {\displaystyle _0^\pi }{\displaystyle \frac{e(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}}{\displaystyle \frac{e}{t}}(t,\vartheta ,\phi )\mathrm{sin}\vartheta d\vartheta +J_1(t),`$
$`J_1(t)`$ $`=`$ $`{\displaystyle \frac{\mu m\mathrm{\Xi }_0^2}{8\pi }}[1+\nu (t)]^2{\displaystyle _0^{2\pi }}𝑑\phi {\displaystyle _0^\pi }{\displaystyle \frac{e^2(t,\vartheta ,\phi )}{[1n(t,\vartheta ,\phi )]^2}}{\displaystyle \frac{n}{t}}(t,\vartheta ,\phi )\mathrm{sin}\vartheta d\vartheta .`$ (16)
It follows from Eq. (8) that

$$\frac{e}{t}(t,\overline{l})=\frac{1}{2[\overline{l}\widehat{C}(t)\overline{l}]}\left[\overline{l}\frac{d\widehat{C}}{dt}(t)\overline{l}\right]\frac{\eta }{1+\eta n(t,\overline{l})}\frac{n}{t}(t,\overline{l}).$$
(17)
Differentiation of Eq. (6) with respect to time yields

$$\frac{d\widehat{C}}{dt}(t)=\overline{}_0\overline{v}(t)[\overline{}_0\overline{r}(t)]^{}+\overline{}_0\overline{r}(t)[\overline{}_0\overline{v}(t)]^{},$$

where $`\overline{v}(t)=d\overline{r}(t)/dt`$ is the velocity vector for the network. Taking into account that

$$\overline{}_0\overline{v}(t)=\overline{}_0\overline{r}(t)\overline{}(t)\overline{v}(t),$$

where $`\overline{}(t)`$ is the gradient operator in the deformed state at time $`t`$, we obtain

$$\frac{d\widehat{C}}{dt}(t)=2\overline{}_0\overline{r}(t)\widehat{D}(t)[\overline{}_0\overline{r}(t)]^{},$$
(18)
where

$$\widehat{D}(t)=\frac{1}{2}\left[\overline{}(t)\overline{v}(t)+(\overline{}(t)\overline{v}(t))^{}\right]$$

is the rate-of-strain tensor for the network. It follows from Eq. (18) that

$$\overline{l}\frac{d\widehat{C}}{dt}(t)\overline{l}=2\overline{l}\overline{}_0\overline{r}(t)\widehat{D}(t)[\overline{}_0\overline{r}(t)]^{}\overline{l}=2\widehat{F}(t,\overline{l}):\widehat{D}(t),$$
(19)
where

$$\widehat{F}(t,\overline{l})=[\overline{}_0\overline{r}(t)]^{}(\overline{l}\overline{l})\overline{}_0\overline{r}(t)$$
(20)
is the generalized Finger tensor , the colon stands for convolution and $``$ denotes tensor product. Substitution of Eqs. (17) and (19) into Eq. (16) implies that

$`{\displaystyle \frac{dW}{dt}}(t)=\widehat{\mathrm{{\rm Y}}}(t):\widehat{D}(t)J(t),`$
$`\widehat{\mathrm{{\rm Y}}}(t)={\displaystyle \frac{\mu m\mathrm{\Xi }_0^2}{4\pi }}[1+\nu (t)]^2{\displaystyle _0^{2\pi }}𝑑\phi {\displaystyle _0^\pi }{\displaystyle \frac{e(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}}{\displaystyle \frac{\widehat{F}(t,\vartheta ,\phi )}{\overline{l}\widehat{C}(t)\overline{l}}}\mathrm{sin}\vartheta d\vartheta ,`$
$`J(t)={\displaystyle \frac{\mu m\mathrm{\Xi }_0^2}{8\pi }}[1+\nu (t)]^2{\displaystyle _0^{2\pi }}𝑑\phi {\displaystyle _0^\pi }{\displaystyle \frac{H(t,\vartheta ,\phi )}{[1n(t,\vartheta ,\phi )]^2}}{\displaystyle \frac{n}{t}}(t,\vartheta ,\phi )\mathrm{sin}\vartheta d\vartheta ,`$
$`H(t,\vartheta ,\phi )=e(t,\vartheta ,\phi )\left\{{\displaystyle \frac{2\eta [1n(t,\vartheta ,\phi )]}{1+\eta n(t,\vartheta ,\phi )}}e(t,\vartheta ,\phi )\right\}.`$ (21)
## IV Constitutive equations

Observations evidence that under cyclic loading with small frequency (less than 10 Hz), the temperature increment is negligible and temperature $`T`$ remains close to its reference value $`T_0`$ . This means that the effect of temperature on material parameters, as well as thermal expansion of the network may be disregarded. It is assumed that the deformation process is rather slow, which implies that at any instant $`t0`$ thermodynamic potentials are correctly defined. For affine deformation of an incompressible network, the Clausius–Duhem inequality reads

$$T\frac{dQ}{dt}=S\frac{dT}{dt}\frac{d\mathrm{\Psi }}{dt}+\frac{1}{\rho }(\widehat{\sigma }_\mathrm{d}:\widehat{D}\frac{1}{T}\overline{q}\overline{}T)0.$$
(22)
where $`\rho `$ is mass density, $`\overline{q}`$ is the heat flux vector, $`\widehat{\sigma }_\mathrm{d}`$ is the deviatoric component of the Cauchy stress tensor $`\widehat{\sigma }`$, $`\mathrm{\Psi }`$ is the free (Helmholtz) energy, $`S`$ is the entropy and $`Q`$ is the entropy production per unit mass. We accept the following expression for the free energy:

$$\mathrm{\Psi }=\mathrm{\Psi }_0+(cS_0)(TT_0)cT\mathrm{ln}\frac{T}{T_0}+W,$$
(23)
where $`S_0`$ and $`\mathrm{\Psi }_0`$ are the entropy and the free energy in the equilibrium stress-free state at the reference temperature $`T_0`$ and $`c`$ is the specific heat. The second and third terms on the right-hand side of Eq. (23) characterize the energy of thermal motion. Unlike conventional theories of rubber elasticity , the terms associated with configurational entropy of chains are disregarded in Eq. (23). According to Ref. , this approximation is acceptable, provided that

$$\mu k_BTN,$$

where $`k_B`$ is Boltzmann’s constant. In the sequel, we suppose that this inequality is satisfied.

Substituting Eqs. (14) and (23) into Eq. (22) and assuming changes in the concentration of junctions to be rather slow (which means that Eq. (21) for the derivative of the mechanical energy $`W`$ with respect to time may be used), we arrive at the formula

$$T\frac{dQ}{dt}=\left(\frac{\widehat{\sigma }_\mathrm{d}}{\rho }\widehat{\mathrm{{\rm Y}}}\right):\widehat{D}\left(SS_0c\mathrm{ln}\frac{T}{T_0}\right)\frac{dT}{dt}+J\frac{1}{\rho T}\overline{q}\overline{}T0.$$
(24)
Applying the conventional reasoning  to Eq. (24), we find that the expressions in braces vanish. This assertion results in the standard formula for the entropy

$$S=S_0+c\mathrm{ln}\frac{T}{T_0},$$
(25)
where the configurational entropy of the network is neglected, and, together with Eqs. (20) and (21), the constitutive equation

$$\widehat{\sigma }(t)=P(t)\widehat{I}+G[1+\nu (t)]^2[\overline{}_0\overline{r}(t)]^{}_0^{2\pi }𝑑\phi _0^\pi \frac{e(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}\frac{\overline{l}\overline{l}}{\overline{l}\widehat{C}(t)\overline{l}}\mathrm{sin}\vartheta d\vartheta \overline{}_0\overline{r}(t),$$
(26)
where $`P(t)`$ is pressure and $`G=\rho \mu m\mathrm{\Xi }_0^2/(4\pi )`$. We substitute Eqs. (25) and (26) into Eq. (24) and find that the rate of entropy production is nonnegative for an arbitrary loading program, provided that

1. the heat flux vector $`\overline{q}`$ obeys the Fourier law $`\overline{q}=\kappa \overline{}T`$ with a positive thermal diffusivity $`\kappa `$,
2. the function $`n(t,\vartheta ,\phi )`$ satisfies the kinetic equation

$$\frac{n}{t}(t,\vartheta ,\phi )=\alpha (t)e(t,\vartheta ,\phi )\left\{\frac{2\eta [1n(t,\vartheta ,\phi )]}{1+\eta n(t,\vartheta ,\phi )}e(t,\vartheta ,\phi )\right\},n(0,\vartheta ,\phi )=0,$$
(27)
where $`\alpha `$ is a nonnegative function of time.

We postulate that the rate of trans–cis transformation is proportional to the average number of strands in a chain,

$$\alpha (t)=\alpha _1N(t),$$
(28)
where $`\alpha _1>0`$ is a material constant. Equation (28) is explained by the fact that the force opposing transition of a bond from its flexed to extended conformation is driven by the action of neighboring bonds, whereas stresses in these bonds are inversely proportional to the number of strands in a chain, see Eq. (3). Combining Eqs. (27) and (28) and using Eqs. (10) and (12), we arrive at the nonlinear differential equation

$$\frac{n}{t}(t,\vartheta ,\phi )=\alpha _0\frac{e(t,\vartheta ,\phi )}{1+\nu (t)}\left\{\frac{2\eta [1n(t,\vartheta ,\phi )]}{1+\eta n(t,\vartheta ,\phi )}e(t,\vartheta ,\phi )\right\},n(0,\vartheta ,\phi )=0$$
(29)
with $`\alpha _0=\alpha _1/(m\mathrm{\Xi }_0)`$.

## V Evolution of the concentration of junctions

To close the model, the kinetics of stress–induced changes in the concentration of temporary junctions should be described. Denote by $`M_0`$ the total concentration (per unit mass) of entanglements (both loose and tight). Adopting the first-order kinetics for the mechanically-induced evolution of the number of tight entanglements, we assume that the rate of increase in the concentration of active entanglements is proportional to the current number of loose ones,

$$\frac{dM}{dt}(t)=\beta _{}(t)[M_0M(t)],$$
(30)
where $`\beta _{}(t)`$ is a positive function. It follows from Eqs. (12) and (30) that

$$\frac{d\nu }{dt}(t)=\beta _{}(t)[\nu _0\nu (t)],$$
(31)
where the quantity

$$\nu _0=2\frac{M_0M(0)}{\mathrm{\Xi }_0}$$

equals the ratio of the initial number of loose entanglements to the initial number of chains. We postulate that the rate of growth for the number of active entanglements, $`\beta _{}`$, is proportional to the mechanical energy per chain and inversely proportional to the number of strands in a chain,

$$\beta _{}(t)=\beta _0\frac{w_0(t)}{N(t)},$$
(32)
where $`\beta _0>0`$ is a material constant. The proportionality of $`\beta _{}`$ to the average mechanical energy, $`w_0`$, reflects the fact that tightening of loose entanglements is driven by mechanical factors. In general, any measure of straining may be chosen instead of $`w_0`$ on the right-hand side of Eq. (31), see a discussion of this issue in . We introduce the strain energy density, $`w_0`$, by analogy with conventional models for damage of elastomers . The inverse proportionality of $`\beta _{}`$ to the average number of strands in a chain, $`N`$, reflects slowing down of the process of tightening entanglements with an increase in the average length of chains. The latter is driven by the growth of the chains’ mobility (estimated in terms of the number of available configurations).

It follows from Eqs. (10), (12), (31) and (32) that the function $`\nu (t)`$ obeys the kinetic equation

$$\frac{d\nu }{dt}(t)=\frac{\beta }{4\pi }[\nu _0\nu (t)]_0^{2\pi }𝑑\phi _0^\pi \frac{e^2(t,\vartheta ,\phi )}{1n(t,\vartheta ,\phi )}\mathrm{sin}\vartheta d\vartheta ,\nu (0)=0,$$
(33)
where $`\beta =\beta _0\mu /2`$. Governing equations (26), (29) and (33) are determined by 5 adjustable parameters: an analog of the shear modulus $`G`$, the rate of trans–cis transition $`\alpha _0`$, the rate of stress–induced increase in the concentration of active entanglements $`\beta `$, the constant $`\eta `$ which characterizes the end-to-end elongation of a chain driven by transformation of a bond from its flexed to extended conformation, and the ratio, $`\nu _0`$, of the initial number of loose entanglements to the initial number of chains. To determine these quantities, we analyze uniaxial extension of a rod.

## VI Uniaxial tension of a specimen

At uniaxial tension of an incompressible medium, Cartesian coordinates in the deformed state, $`x_i`$, are expressed in terms of the Cartesian coordinates in the stress-free state, $`X_i`$, by the formulas

$$x_1=k(t)X_1,x_2=k^{\frac{1}{2}}(t)X_2,x_3=k^{\frac{1}{2}}(t)X_3,$$

where $`k=k(t)`$ is the extension ratio. It follows from these equalities and Eq. (6) that

$`\overline{}_0\overline{r}(t)`$ $`=`$ $`k(t)\overline{e}_1\overline{e}_1+k^{\frac{1}{2}}(t)(\overline{e}_2\overline{e}_2+\overline{e}_3\overline{e}_3),`$
$`\widehat{C}(t)`$ $`=`$ $`k^2(t)\overline{e}_1\overline{e}_1+k^1(t)(\overline{e}_2\overline{e}_2+\overline{e}_3\overline{e}_3).`$ (34)
The unit vector $`\overline{l}`$ is given by

$$\overline{l}=\mathrm{cos}\vartheta \overline{e}_1+\mathrm{sin}\vartheta (\mathrm{cos}\phi \overline{e}_2+\mathrm{sin}\phi \overline{e}_3),$$
(35)
which implies that the tensor $`\overline{l}\overline{l}`$ is determined by the matrix

$$\overline{l}\overline{l}=\left[\begin{array}{ccc}\mathrm{cos}^2\vartheta & \mathrm{sin}\vartheta \mathrm{cos}\vartheta \mathrm{cos}\phi & \mathrm{sin}\vartheta \mathrm{cos}\vartheta \mathrm{sin}\phi \\ \mathrm{sin}\vartheta \mathrm{cos}\vartheta \mathrm{cos}\phi & \mathrm{sin}^2\vartheta \mathrm{cos}^2\phi & \mathrm{sin}^2\vartheta \mathrm{sin}\phi \mathrm{cos}\phi \\ \mathrm{sin}\vartheta \mathrm{cos}\vartheta \mathrm{sin}\phi & \mathrm{sin}^2\vartheta \mathrm{sin}\phi \mathrm{cos}\phi & \mathrm{sin}^2\vartheta \mathrm{sin}^2\phi \end{array}\right].$$
(36)
Equations (34) and (35) yield

$$\overline{l}\widehat{C}(t)\overline{l}=k^2(t)\mathrm{cos}^2\vartheta +k^1(t)\mathrm{sin}^2\vartheta .$$
(37)
Because of the axial symmetry of deformation, the functions $`n=n(t,\vartheta )`$ and $`e=e(t,\vartheta )`$ are independent of $`\phi `$, which implies that all terms but $`\overline{l}\overline{l}`$ on the right-hand side of Eq. (26) are independent of $`\phi `$. It follows from Eq. (36) that

$$\frac{1}{2\pi }_0^{2\pi }\overline{l}\overline{l}𝑑\phi =\mathrm{cos}^2\vartheta \overline{e}_1\overline{e}_1+\frac{1}{2}\mathrm{sin}^2\vartheta (\overline{e}_2\overline{e}_2+\overline{e}_3\overline{e}_3).$$

Substitution of this expression and Eqs. (34) and (35) into Eq. (26) implies that

$$\widehat{\sigma }(t)=\sigma (t)\overline{e}_1\overline{e}_1+\sigma _0(t)(\overline{e}_2\overline{e}_2+\overline{e}_3\overline{e}_3),$$

where

$`\sigma (t)`$ $`=`$ $`P(t)+2\pi G[1+\nu (t)]^2{\displaystyle _0^\pi }{\displaystyle \frac{k^2(t)\mathrm{cos}^2\vartheta }{k^2(t)\mathrm{cos}^2\vartheta +k^1(t)\mathrm{sin}^2\vartheta }}{\displaystyle \frac{e(t,\vartheta )}{1n(t,\vartheta )}}\mathrm{sin}\vartheta d\vartheta ,`$
$`\sigma _0(t)`$ $`=`$ $`P(t)+\pi G[1+\nu (t)]^2{\displaystyle _0^\pi }{\displaystyle \frac{k^1(t)\mathrm{sin}^2\vartheta }{k^2(t)\mathrm{cos}^2\vartheta +k^1(t)\mathrm{sin}^2\vartheta }}{\displaystyle \frac{e(t,\vartheta )}{1n(t,\vartheta )}}\mathrm{sin}\vartheta d\vartheta .`$ (38)
The boundary condition on the lateral surface of the specimen implies that $`\sigma _0(t)=0`$. Combining this equality with Eq. (38), we find the longitudinal stress

$$\sigma (t)=\pi G[1+\nu (t)]^2_0^\pi \frac{2k^2(t)\mathrm{cos}^2\vartheta k^1(t)\mathrm{sin}^2\vartheta }{k^2(t)\mathrm{cos}^2\vartheta +k^1(t)\mathrm{sin}^2\vartheta }\frac{e(t,\vartheta )}{1n(t,\vartheta )}\mathrm{sin}\vartheta d\vartheta .$$

Introducing the notation

$$z=\mathrm{cos}\vartheta ,\stackrel{~}{e}(t,z)=e(t,\vartheta ),\stackrel{~}{n}(t,z)=n(t,\vartheta )$$

and bearing in mind that $`\stackrel{~}{e}`$ and $`\stackrel{~}{n}`$ are even functions of $`z`$, we obtain

$$\sigma (t)=E[1+\nu (t)]^2_0^1\frac{2k^2(t)z^2k^1(t)(1z^2)}{k^2(t)z^2+k^1(t)(1z^2)}\frac{\stackrel{~}{e}(t,z)}{1\stackrel{~}{n}(t,z)}𝑑z,$$
(39)
where

$$E=\frac{1}{2}\rho \mu m\mathrm{\Xi }_0^2.$$
(40)
It follows from Eqs. (8), (34) and (35) that

$$\stackrel{~}{e}(t,z)=\mathrm{ln}\frac{[k^2(t)z^2+k^1(t)(1z^2)]^{\frac{1}{2}}}{1+\eta \stackrel{~}{n}(t,z)}.$$
(41)
In the sequel, we focus on stretching with a constant rate of engineering strain, $`\dot{ϵ}_0>0`$,

$$k(t)=1+\dot{ϵ}_0t.$$
(42)
Combining Eqs. (29) and (42), we arrive at the kinetic equation

$$\frac{\stackrel{~}{n}}{k}=a\left[\frac{2\eta (1\stackrel{~}{n})}{1+\eta \stackrel{~}{n}}\stackrel{~}{e}\right]\stackrel{~}{e},\stackrel{~}{n}(1,z)=0$$
(43)
where $`a=\alpha /\dot{ϵ}_0`$. Substituting expression (42) into Eq. (33) and calculating the integral over $`\phi `$, we find that

$$\frac{d\nu }{dk}=b(\nu _0\nu )_0^1\frac{\stackrel{~}{e}^2dz}{1\stackrel{~}{n}},\nu (1)=0,$$
(44)
where $`b=\beta /\dot{ϵ}_0`$. Equations (39), (41), (43) and (44) are determined by 5 adjustable parameters: $`E`$, $`a`$, $`b`$, $`\eta `$ and $`\nu _0`$. This number is quite comparable with the number of experimental constants employed in conventional stress–strain relations in finite elasticity of elastomers .

## VII Comparison with experimental data

To determine adjustable parameters, we fit observations for natural rubber vulcanizates in uniaxial tensile tests with the strain rate $`\dot{ϵ}_0=2.0`$ min<sup>-1</sup> at room temperature. For a description of specimens and the experimental procedure, see Refs. .

We begin by matching experimental data in cyclic tests with the maximal elongation $`k_{\mathrm{max}}=6.0`$ for a vulcanizate with an unspecified composition. First, we approximate the stress–strain curve for a virgin sample. Given a Young’s modulus, $`E`$, the constants $`a`$, $`b`$, $`\eta `$ and $`\nu _0`$ are found by the steepest-descent algorithm. The parameter $`E`$ is determined by the least-squares technique. Afterwards, we fix the values $`b`$ and $`\nu _0`$ (which are responsible for the evolution of the concentration of junctions) and match stress–strain curves measured after retraction by using only three constants: $`E`$, $`a`$ and $`\eta `$. Experimental data in tests with various numbers of cycles, $`i`$, are depicted in Figure 1 together with results of numerical simulation. Material constants are listed in Table 1. The parameters $`E`$, $`a`$ and $`\eta `$ are plotted versus the number of cycles $`i`$ in Figures 2 and 3, which show that observations are fairly well approximated by the stretched exponential functions

$$E=E_010^{i^{\kappa _E}},a=a_010^{i^{\kappa _a}},\eta =\eta _010^{i^{\kappa _\eta }},$$
(45)
where the constants $`E_0`$, $`a_0`$, $`\eta _0`$ and $`\kappa _E`$, $`\kappa _a`$, $`\kappa _\eta `$ are found by the least-squares algorithm.

We proceed by fitting observations for a natural rubber vulcanizate with an unspecified composition in tensile loading–unloading tests with an increasing maximal amplitude of stretching, $`k_{\mathrm{max}}`$. First, we approximate experimental data for a virgin specimen and determine adjustable parameters using the steepest-descent procedure. Afterwards, we fix the constants $`b`$ and $`\nu _0`$ and match stress–strain curves measured after retraction with the help of three adjustable parameters: $`E`$, $`a`$ and $`\eta `$. Figure 4 demonstrates good agreement between experimental data and results of numerical simulation with adjustable parameters collected in Table 2. The quantities $`E`$, $`a`$ and $`\eta `$ are plotted in Figures 5 to 7 versus the number of cycles $`i`$. These figures reveal an acceptable quality of fitting experimental data by Eq. (45).

To evaluate the influence of annealing (24 h at the temperature $`T=100`$ C) on adjustable parameters in the constitutive equations, we approximate observations in tensile tests for three vulcanizates: A (polysulfide crosslinks), B (monosulfide crosslinks) and C (carbon–carbon crosslinks) with various amounts of crosslinkers (sulfur and dicumyl peroxide, respectively). A detailed description of the composition of specimens in presented in . For any type of vulcanizates, we begin by fitting the stress–strain curve for a virgin sample with the minimum content of crosslinker, $`\varphi `$, and find experimental constants by using the steepest-descent algorithm. Afterwards, we fix the quantities $`b`$ and $`\nu _0`$ and match other observations by using only 3 adjustable parameters: $`E`$, $`a`$ and $`\eta `$. Experimental data together with results of numerical simulation are depicted in Figures 8 to 10 for vulcanizate A, in Figures 11 to 13 for vulcanizate B and in Figures 14 and 15 for vulcanizate C. Experimental constants are collected in Tables 3 to 5.

To compare the effects of thermal recovery and recovery by swelling, we match experimental data in tensile tests for a virgin specimen and for the same specimen after retraction (vulcanizates A, B and C with unspecified compositions), and repeat the same procedure for a recovered sample (swolen in benzene for 24 h and dried in vacuum). Adjustable parameters are found by using the same numerical procedure as for the specimens annealed at an elevated temperature. Observations and results of numerical simulation are depicted in Figures 16 to 18, whereas material constants are listed in Table 6.

## VIII Discussion

Figures 1, 4, 8 to 18 demonstrate fair agreement between experimental data and results of numerical analysis, which implies that the model may be applied to fit observations in uniaxial tensile tests with the axial elongation up to $`k_{\mathrm{max}}=8.0`$.

Figures 3 and 7 show that the parameter $`\eta `$ increases with the number of cycles, $`i`$. This conclusion is also confirmed by the data listed in Tables 3 to 6 (results for the first and second stretching). On the contrary, recovery of specimens (by annealing at an elevated temperature and by swelling) reduces $`\eta `$. This means that the quantity $`\eta `$ may be treated as a parameter responsible for material damage at the micro-level. Because $`\eta `$ reflects the average end-to-end elongation of a chain driven by an individual trans–cis transition (which implies that it may be thought of as an average measure of coiling of a chain), an increase in $`\eta `$ under cyclic stretching is tantamount to an increase in the average size of a globule formed by a chain. This picture results in the conclusion that recovery of natural rubber vulcanizates may be treated as coiling of long chains. Tables 3 to 5 evidence that this process occurs more intensively for rubbery polymers with high concentrations of crosslinkers, and vulcanizates with dicumyl peroxide are annealed more effectively than those with sulfide crosslinks. These tables also show that the parameter $`\eta `$ monotonically decreases with the growth of the content of crosslinker, $`\varphi `$. Because the growth of the concentration of chemical crosslinks results in an increase in the volume fraction of junctions (which is confirmed by the growth of Young’s modulus with $`\varphi `$), and, as a consequence, a decrease in the average length of a chain, $`N`$, we may conclude that short chains are coiled up relatively more strongly \[the qualification “relatively” is important here, because $`\eta `$ is multiplied by the number of strands in a chain, $`N`$, in Eq. (7)\] than longer ones.

Tables 3 to 6 show that at the second stretching the rate of transition from the flexed to extended conformation of bonds, $`a`$, exceeds that for a virgin specimen (the data for vulcanizate B in Table 6 are the only exception from this rule). This observation may also be treated as some kind of material damage: if we take a bond as a hinge, then straightening two rods linked by the hinge occurs more easily for a preloaded sample than for a virgin one. Tables 3 to 5 reveal that $`a`$ substantially grows with the concentration of crosslinker, $`\varphi `$. This increase may be associated with appropriate changes in $`\eta `$: an increase in the degree of coiling makes trans–cis transitions easier. In most cases, annealing of specimens (thermal and by swelling) results in a decrease in $`a`$ which is accompanied by a decrease in $`\eta `$. However, for vulcanizate A with high concentrations of crosslinker, an inverse tendency is observed: despite a decrease in $`\eta `$, the rate of transition from the flexed to extended conformation increases, see Tables 3 and 6.

Young’s modulus $`E`$ strongly increases with the content of crosslinker. This conclusion is in accord with Eq. (40) which implies that $`E`$ is proportional to the square of the initial number of chains. Tables 3 to 6 demonstrate that the elastic modulus for a recovered specimen exceeds that for a prestrained material. The growth of $`E`$ during recovery seems quite natural, provided that the recovery process is associated with an increase in the concentration of entanglements.

Unlike experimental data collected in Table 2, Table 1 demonstrates that cyclic loading results in a pronounced decrease in $`a`$ and $`E`$ with the number of cycles, $`i`$. Figure 2 reveals that this decrease may be fairly well approximated by the stretched exponential law. To explain the difference between the data listed in Tables 1 and 2, we refer to Figures 5 to 7, which evidence that $`E`$ and $`a`$ decrease with $`i`$ at small amplitudes of stretching and strongly increase at large amplitudes. Changes in adjustable parameters observed in the cyclic test with an increasing amplitude are similar to those presented in Table 1 during the first 4 cycles and are analogous to those presented in Tables 3 to 5 when the maximum amplitude of stretching, $`k_{\mathrm{max}}`$, exceeds 4. One can only speculate about a physical mechanism for this phenomenon. As a possible explanation, we suppose that cyclic loading affects the internal structure of a rubber vulcanizate in different ways at small and large amplitudes of stretching. Periodic loading with amplitudes which do not exceed some threshold results in material training similar to that occurring at isothermal annealing (Table 1). On the contrary, cyclic stretching with amplitudes exceeding the threshold strain results in damage of the internal structure which is reflected by the data collected in Tables 2 to 5.

For each series of tests, the parameters $`b`$ and $`\nu _0`$ are found from the condition of the best fit for one stress–strain curve and remain fixed in matching observations for other curves. Excellent agreement between results of numerical simulation and experimental data leads to the conclusion that these parameters are not affected by mechanical factors.

Table 6 reveals that the differences between the parameters $`a`$ and $`E`$ found by fitting observations on recovered specimens at the first and second stretching are rather small compared with those for virgin samples. This means that recovery by swelling is equivalent to equilibration of rubbery polymers. On the contrary, analogous differences in $`\eta `$ remain practically constant, which confirms the main hypothesis of this study that stress–softening of elastomers may be associated with uncoiling of long chains.

## IX Conclusions

Constitutive equations have been derived for the isothermal mechanical response of elastomers at finite strains. Stress–strain relations are developed using the laws of thermodynamics and are applied to fit experimental data in uniaxial tensile tests for several natural rubber vulcanizates. Fair agreement is demonstrated between experimental data and results of numerical simulation.

We study the effects of cyclic loading, annealing at elevated temperature and recovery by swelling on material constants. The following conclusions are drawn:

1. Stress-softening of elastomers is reflected by changes in the parameter $`\eta `$ which characterizes a mechanically-induced increase in the average size of globules formed by long chains. The quantity $`\eta `$ is increased under periodic stretching with large amplitudes (uncoiling of chains) and is decreased at recovery (by annealing or swelling).
2. In contrast to conventional models for stress-softening of rubbery polymers, fitting of observations demonstrates that Young’s modulus is altered relatively weakly. This implies that damage should be associated rather with changes in the topology of a temporary network than with rupture of long chains.
3. Cyclic straining results in a noticeable increase in the rate of trans–cis transition of bonds, $`a`$, with the number of cycles, $`i`$ (which is fairly well described by the stretched exponential law). On the contrary, recovery of specimens (by thermal annealing and swelling) induces a decrease in $`a`$.
4. Young’s modulus, $`E`$, and the rate of transformation of bonds from their flexed to extended conformation, $`a`$, substantially increase with the content of crosslinker, $`\varphi `$, whereas the parameter $`\eta `$ is a decreasing function of $`\varphi `$.
5. Adjustable parameters change in a similar way during thermal annealing of specimens at an elevated temperature and during recovery by swelling.
6. Cyclic loadings with relatively small and large amplitudes affect the internal structure of elastomers in different ways. Periodic straining with large amplitudes causes damage to the internal structure, whereas the influence of that with small amplitudes is tantamount to recovery of vulcanizate and may be thought of as a kind of training of rubbery polymers.
7. The material parameters $`b`$ and $`\nu _0`$ which reflect the kinetics of mechanically-induced transition of loose entanglements into the active state are rather robust. They weakly depend on the intensity of loading, as well as on the composition of vulcanizates and on the content of crosslinkers.
8. Recovery by swelling reduces the increments of Young’s modulus and the rate of trans–cis transition driven by a cycle of loading. This implies that swelling of a rubber vulcanizate induces equilibration of its mechanical properties.

The above assertions are obtained in the framework of the phenomenological model, which implies that they should be treated with caution. These conclusions may, however, be taken as rather plausible assumptions (mechanically-induced uncoiling of chains, training of an elastomer by sub-threshold periodic stretching, acceleration of trans–cis transformations under loading, disentanglement at retraction, etc.) which have to be checked experimentally. Validation of these hypotheses will be the subject of a subsequent study.

### Acknowledgement

ADD acknowledges financial support by the Israeli Ministry of Science through grant 1202–1–00.

## Figure legends

Figure 1: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber with an unspecified composition. Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: $`i=1`$; curve 2: $`i=2`$; curve 3: $`i=3`$; curve 4: $`i=8`$

Figure 2: Young’s modulus $`E`$ MPa (unfilled circles) and the dimensionless parameter $`a`$ (filled circles) versus the number of cycles $`i`$. Symbols: treatment of observations . Solid lines: approximation of the experimental data by Eq. (45). Curve 1: $`E_0=12.3705`$, $`\kappa _E=0.0716`$; curve 2: $`a_0=9.2561`$, $`\kappa _a=0.1219`$

Figure 3: The dimensionless parameter $`\eta `$ versus the number of cycles $`i`$. Cirles: treatment of observations . Solid lines: approximation of the experimental data by Eq. (45) with $`\eta _0=4.1975`$ and $`\kappa _\eta =0.2055`$

Figure 4: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate with an unknown composition. Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: a virgin specimen; curve 2: the 4th stretching; curve 3: the 5th stretching; curve 4: the 6th stretching; curve 5: the 7th stretching; curve 6: the 8th stretching; curve 7: the 9th stretching

Figure 5: Young’s modulus $`E`$ MPa versus the number of cycles $`i`$. Circles: treatment of observations . Solid lines: approximation of the experimental data by Eq. (45). Curve 1: $`E_0=21.5638`$, $`\kappa _E=0.0125`$; curve 2: $`E_0=3.3807`$, $`\kappa _E=0.5003`$

Figure 6: The dimensionless parameter $`a`$ versus the number of cycles $`i`$. Circles: treatment of observations . Solid lines: approximation of the experimental data by Eq. (45). Curve 1: $`a_0=9.9095`$, $`\kappa _a=0.0863`$; curve 2: $`a_0=3.7108`$, $`\kappa _a=0.3952`$

Figure 7: The dimensionless parameter $`\eta `$ versus the number of cycles $`i`$. Circles: treatment of observations . Solid lines: approximation of the experimental data by Eq. (45). Curve 1: $`\eta _0=4.6700`$, $`\kappa _\eta =0.0035`$; curve 2: $`\eta _0=2.7620`$, $`\kappa _\eta =0.2524`$

Figure 8: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate A ($`\varphi =1.25`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 9: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate A ($`\varphi =2.50`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 10: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate A ($`\varphi =4.17`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 11: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate B ($`\varphi =0.2`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 12: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate B ($`\varphi =0.4`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 13: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate B ($`\varphi =0.6`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 14: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate C ($`\varphi =0.5`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 15: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate C ($`\varphi =3.5`$ phr). Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: the first stretching; curve 2: the second stretching; curve 3: stretching after annealing

Figure 16: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate A with an unspecified composition. Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: a virgin sample, the first stretching; curve 2: a virgin sample, the second stretching; curve 3: a recovered sample, the first stretching; curve 4: a recovered sample, the second stretching

Figure 17: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate B with an unspecified composition. Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: a virgin sample, the first stretching; curve 2: a virgin sample, the second stretching; curve 3: a recovered sample, the first stretching; curve 4: a recovered sample, the second stretching

Figure 18: The true stress $`\sigma `$ MPa versus the engineering strain $`ϵ`$ for natural rubber vulcanizate C with an unspecified composition. Circles: experimental data . Solid lines: results of numerical simulation. Curve 1: a virgin sample, the first stretching; curve 2: a virgin sample, the second stretching; curve 3: a recovered sample, the first stretching; curve 4: a recovered sample, the second stretching


text:
# Caustic Crossing in the Gravitational Lens Q2237+0305

##

The monitoring of the gravitational lens Q2237+0305 carried out by the OGLE group during 1997–1999 is analyzed. The significant light amplifications in the C and A quasar components with maxima in mid- and late 1999, respectively, are interpreted as the crossing of microlens caustics by the source. A constraint on the emitting-region size $`R10^{15}`$ cm has been obtained from the light-curve shape by assuming a power-law quasar brightness distribution $`(r^2+R^2)^p`$. To estimate the exponent $`p1.2`$ requires refining the standard model for the quasar continuum formation in an optically thick accretion disk with $`p=1.5`$.

##

Key words: quasars, gravitational lenses, accretion disks

## 1 INTRODUCTION

Quasar variability under the effect of microlenses depends both on parameters of the mass distribution for compact bodies and on the appearance of the emitting region. For instance, the larger is the quasar, the smaller are the amplitudes of light variations in its images. Thus, it becomes possible to formulate the inverse problem of determining the sizes and structure of quasars from their observed light curves. The main difficulty in solving this problem is that the distribution of the amplification produced by microlenses is not known in advance and exhibits a fairly complex pattern with many randomly located caustic lines (see, e.g., Zakharov 1997 and Zakharov and Sazhin 1998). In general, the specific form of this distribution is not known and can be analyzed only statistically. Exceptions are only those portions of the light curve that correspond to caustic crossing by the quasar. The amplification of a point source during caustic crossing obeys a simple law: it remains approximately constant as the caustic is approached, then abruptly increases to infinity at the caustic itself, and subsequently falls off as $`x^{1/2}`$ with increasing distance $`x`$ from the caustic (Chang and Refsdal 1984; Blandford and Narayan 1986).

Microlenses must show up most clearly in the multiple quasar images produced by the gravitational effect of galaxy macrolenses. First, the microlensing probability is rather high in such situations. Second, intrinsic quasar variability can, in principle, be separated from microlensing variability. Of all the gravitationally lensed quasars, Q2237+0305 ($`z_s=1.675`$) is undoubtedly the most promising object for microlensing analysis. Because the unique proximity of a lensing galaxy ($`z_l=0.039`$), microlensing variability in this object must take place faster than in other objects by an order of magnitude and with a large amplitude.

The quasar Q2237+0305 has been monitored virtually since its discovery, and it actually proved to be the first object in which microlensing variability was detected (Irwin et al. 1989; Corrigan et al. 1991). The observations by Østensen et al. (1996) showed that virtually all four quasar images were more or less variable.

The regularity and quality of the Q2237+0305 monitoring have improved markedly when the OGLE group joined it in the last four years (Wozniak et al. 2000) (see Fig. 1). Measurements are made in the V band approximately once a week during the observing season from May through December. The latest observations (http://www.astro.princeton.edu/$``$ogle/ogle2/huchra.html) show that image C passed its intensity peak in mid-1999, while image A peaked in late 1999. Interpreting the light-curve maxima as resulting from caustic crossing allows the size and structure of the quasar emitting region in the object under study to be determined.

## 2 AMPLIFICATION OF AN EXTENDED
SOURCE IN THE CAUSTIC REGION

When a source crosses a caustic line, an additional pair of images appears (or disappears). The total intensity of this pair depends on the distance to the caustic as $`x^{1/2}`$. Therefore, the intensity of a point source in the caustic region can be represented as (Schneider and Weiss 1987)

$$I_p(x)=I_0+\theta (x)a_0x^{1/2}.$$
(1)
Here, $`I_0`$ is the intensity of all the remaining microimages except the additional pair, $`\theta (x)`$ is the Heaviside unit function, and $`a_0`$ is the caustic strength. The amplification of an extended source with a brightness distribution $`P(𝒓)`$ can be calculated by ordinary summation over the set of infinitely small sources with individual amplification factors:

$$I=\left[𝑑𝒓P(𝒓)I_p(𝒓)\right]/\left[𝑑𝒓P(𝒓)\right].$$
(2)
To describe the brightness distribution in the source, we use a power-law model

$$P(r)=(1+\frac{r^2}{R^2})^p(p>0),$$
(3)
which is determined by the source radius $`R`$ and by the rate of brightness decline $`p`$.

Let the source center be at a distance $`D`$ from the caustic line. In the normalized coordinates $`\xi =x/R`$ and $`\eta =y/R`$ and using normalized distance $`d=D/R`$, we obtain

$$I=I_0+a_0R^{1/2}J(d),$$
(4)
where the function $`J(d)`$ is

$$J(d)=\left[_{\mathrm{}}^+\mathrm{}𝑑\eta _0^+\mathrm{}𝑑\xi P(\xi ,\eta )\xi ^{1/2}\right]/\left[_{\mathrm{}}^+\mathrm{}𝑑\eta _{\mathrm{}}^+\mathrm{}𝑑\xi P(\xi ,\eta )\right]$$
(5)
or, taking into account the symmetry in $`\eta `$ and the total intensity

$$P(r)𝑑𝒓=\frac{\pi R^2}{p1},$$
(6)
we obtain:

$$J(d)=2\frac{p1}{\pi }_0^{\mathrm{}}𝑑\xi _0^{\mathrm{}}𝑑\eta \xi ^{\frac{1}{2}}\left[1+(\xi d)^2+\eta ^2\right]^p.$$
(7)
Taking the internal integral over $`\eta `$ yields

$$J(d)=\frac{p1}{\pi }_0^{\mathrm{}}𝑑\xi B(\frac{1}{2},p\frac{1}{2})\xi ^{1/2}\left[\xi ^22\xi d+d^2+1\right]^{\frac{1}{2}p},$$
(8)
where $`B`$ is the beta function. The subsequent integration over $`\xi `$ yields

$$J(d)=\frac{\mathrm{\Gamma }(p\frac{1}{2})\mathrm{\Gamma }(2p\frac{3}{2})}{\mathrm{\Gamma }(p1)\mathrm{\Gamma }(2p1)}(1+d^2)^{\frac{3}{4}p}{}_{2}{}^{}F_{1}^{}(\frac{1}{2},2p\frac{3}{2};p;\frac{1}{2}\left(1+\frac{d}{\sqrt{1+d^2}}\right)).$$
(9)
Here, $`\mathrm{\Gamma }`$ and $`{}_{2}{}^{}F_{1}^{}(a,b;c;z)`$ are the gamma function and the Gauss hypergeometric function, respectively. Figure 2 shows the function $`J(d)`$ for several values of $`p`$. We clearly see from the figure that the sharpness of the jump in amplitude increases with increasing source brightness concentration toward the center during caustic crossing, tending to an infinite point-source limit for very large $`p`$.

For some particular $`p`$ values, for example, for $`p=3/2,2`$ and $`3`$, the hypergeometric function can be expressed as a combination of elementary and other special functions. Thus, we have for $`p=3/2`$

$$J(d)_{|p=3/2}=\left[2(1+d^2)((1+d^2)^{1/2}d)\right]^{1/2}.$$
(10)
## 3 VARIABILITY ANALYSIS

We see from Fig. 1 that the intensities of all four quasar images have varied during the last four years. The largest variations were observed in image C, which passed its maximum in mid-1999, and in image A, which reached its maximum in late 1999. The approach to interpreting the variability of these two images is the same. Let us assume that, in both cases, the source crosses the caustic; the brightness distribution must follow the law (4). The general form of the curve depends on five parameters:

1. Contribution $`I_0`$ from the remaining microimages, which is assumed to be approximately constant;
2. Caustic strength $`a_0`$;
3. The time it takes for the source to cross its radius $`\mathrm{\Delta }t`$ of the source, which is proportional to the source size $`R`$;
4. The time $`t_0`$ of caustic crossing by the source center;
5. Exponent $`p`$ in the brightness distribution (3).

Estimating the five parameters reduces to minimizing the sum of the squared differences between model and observed light curves

$$\chi ^2=\underset{i=1}{\overset{N}{}}\frac{1}{\sigma _i^2}\left[I_{model}(I_0,a_0,\mathrm{\Delta }t,t_0,p)I_{obs}\right]^2.$$
(11)
The summation is performed over all $`N`$ points of observations with weights that are inversely proportional to the squares of the observational errors $`\sigma _i`$. The best set of parameters is sought by the Marquardt method (see, e.g., Chapter 15.5 in the book by Press et al. 1992).

### 3.1 Image A

The model parameters estimated from the light curve of image A with 182 data points during the entire observing period 1997-1999 are given in the Table 1. The combination $`a_0/\mathrm{\Delta }t`$ is more convenient to calculate than the caustic strength $`a_0`$. The listed formal accuracies of the parameter estimates should be considered only as their lower limits.

Figure 3 shows the model light curve together with measured values. The 1997 observations are poorly fitted by a single curve, implying that approaching the caustic did not show up in the first year. Such a behavior is characteristic of numerical microlensing models for Q2237+0305, in which one amplification event is often superimposed on another to form complex light curves. Excluding the 1997 data from our analysis, while significantly improving the total $`\chi ^2`$ residual, affects the parameter estimates only slightly.

When the caustic is crossed in image A, an additional pair of images emerges. A characteristic feature of this direction of motion is a steep rise in light amplification followed by a gentler decline. The dashed line in Fig. 3 represents the expected behavior of image A in the immediate future. We assume that the brightness will continue to decline and (if no additional causes of amplification arise) will reach the original 1997 level in two to three years. The latest observations for 2000 are added for comparison.

### 3.2 Image C

Attempts to fit the light curve of image C over the entire observed period have failed. At the same time, excluding 80 data points for 1997 from our analysis results in quite reasonable estimates, which are given in the table and shown in Fig. 4. Interestingly, the best solution corresponds to caustic crossing in the negative direction, with the pair of images disappearing. The local minimum of $`\chi ^2`$ corresponding to the motion in the positive direction is several times greater than the absolute minimum reached during the motion in the negative direction.

The inability to fit the entire observed period by a single curve becomes more understandable in light of the recent results by Wyithe et al. 2000. These authors argue that the intensity variations in image C have been composite in pattern during the last four years, and, what is quite possible, another caustic crossing was overlooked during the observing seasons in 1997 and 1998.

The dotted and dashed lines indicate the best extrapolation computed in the model with single caustic crossing in 1997 and our forecast until the end of 2000, respectively. No appreciable intensity variations in image C are expected in the immediate future.

## 4 THE SIZE AND STRUCTURE OF
THE QUASAR EMITTING REGION

For the quasar radius estimates $`\mathrm{\Delta }t`$ in time units to be transformed to linear sizes $`R`$, we must know the apparent quasar velocity; only the velocity component $`v_{}`$ is perpendicular to the caustic line is important. Of course, the exact velocity is unknown. Nevertheless, a statistical analysis of the time derivatives of brightness variations by Wyithe et al. 1999 yielded an upper limit of $`v<500`$ km/s. The perpendicular velocity component can only be lower than this value. The most probable value of $`v_{}`$ was found to be $`300`$ km/s. Given the angular distances from the observer to the quasar $`D_{os}`$ and to the lens $`D_{ol}`$ in the model of a flat Universe with $`H_0=60`$ km/s $`Mpc^1`$, the most probable source radius is estimated to be

$$R=v_{}\mathrm{\Delta }t\frac{D_{os}}{D_{ol}}2.010^{13}\left(\frac{v_{}}{300\text{km/s}}\right)\left(\frac{\mathrm{\Delta }t}{1\text{day}}\right)\text{cm}.$$
(12)
Using the crossing time $`\mathrm{\Delta }t_A=90`$ days for image A gives the most probable quasar radius $`R=1.810^{15}`$ cm, while substituting the crossing time for image C, $`\mathrm{\Delta }t_C=30`$ days, reduces this estimate to $`R=610^{14}`$ cm. At the same time, the velocity constraint $`v<500`$ km/s together with the crossing time $`\mathrm{\Delta }t_C`$ give an upper limit on the source radius, $`R10^{15}`$ cm.

Another parameter that characterizes the mass distribution is the exponent $`p`$. It follows from the table that its computed value lies between 1.1 and 1.25. It is interesting to note that different models for the quasar structure can lead to different dependences of the emissivity on distance from the center. For instance, a power-law dependence with $`p0.5`$ follows from the model of an optically thin accretion disk (Manmoto et al. 1997). At the same time, the standard model of an optically thick accretion disk yields an $`r^3`$ dependence (Shakura and Sunyaev 1973), which changes to $`(r^2+R^2)^{3/2}`$ for a finite radius. Our estimate $`p1.2`$ favors the standard model, but more complex accretion-disk models should be analyzed to achieve better agreement.

## 5 DISCUSSION

Let us consider the legitimacy of some assumptions made here. The main assumption is associated with the hypothesis of caustic crossing. Numerical calculations show that there are two effects capable of causing a significant increase in amplification. Apart from fold-caustic crossing, the source can also pass near the caustic cusp. However, events of the second type for the images of Q2237+0305 are several times less probable than those of the first type (Wambsganss et al. 1992; Lewis and Irwin 1996) and, in general, are more symmetric. These two properties can serve as a statistical justification for using the hypothesis of caustic crossing.

A power-law model with a finite core radius is used to calculate the quasar size and structure. Three models of a caustic-crossing source are encountered in the literature: a homogeneous disk, a Gaussian source (Schneider and Weiss 1987), and a $`(1r^2/R^2)^{1/2}`$ distribution (Schneider and Wagoner 1987). All of them are completely determined by their radii. The $`(1+r^2/R^2)^p`$ model differs radically in that it is a two-parameter model. Estimating the rate of brightness decline $`p`$ allows us to choose between different models for the quasar structure.

Allowance for the spatial orientation of the emitting region appreciably complicates the analysis. Two additional parameters associated with the orientation-angle components of the initial circular source appear. The first and second parameters determine, respectively, the apparent-ellipse eccentricity and the angle of motion with respect to the perpendicular to the caustic (in this case, allowance for the direction of motion does not reduce to a simple substitution of the transverse velocity). Besides, accretion-disk rotation can cause an additional asymmetry of the emitting region in the spectral range considered due to the Doppler effect.

In view of many influencing factors, the possibility of reconstructing, at least in principle, the one-dimensional quasar profile (along the $`x`$-axis) as the solution of an integral equation (Grieger et al. 1991; Mineshige and Yonehara 1999; Agol and Krolik 1999) seems of great interest. However, this is a separate, independent problem, which is yet to be applied to an actual monitoring.

The simple amplification behavior $`x^{1/2}`$ for a point source located in the caustic region is possible only when several conditions are satisfied. The radius of caustic curvature must be considerably larger than the source size (see Fluke and Webster 1999 for curvature allowance). The caustic must be isolated lest the source cover several caustics simultaneously. In addition, it is implied that there is no large-scale time gradient in amplification variations. Thus, for example, introducing a constant slope as an additional free parameter allows a model curve of image C to be easily fitted to the entire observing period 1997–1999. However, since there were no such large gradients throughout the entire 15-year-long monitoring history of Q2237+0305 and since their physical origin is not completely understood, we have to abandon the additional parameter.

Finally, the problem can be further complicated by intrinsic quasar variability. In general, the latter is rather difficult to separate from microlensing variability. However, given that the delays between the images in our cases do not exceed several days, intrinsic quasar variability must be repeated in all quasar images (with individual amplification factors) virtually simultaneously. The fact that the brightness variations in all four components are not synchronous provides circumstantial evidence that intrinsic quasar variability is negligible in this source.

Despite the possible complicating factors listed above, we have every reason to believe that the simple model used here is capable of faithfully reproducing the observational data, and that its implications deserve a careful analysis.

After this paper was mainly complete, Wyithe et al. 2000 independently published a preprint where they also interpreted the OGLE-group observations of Q2237+0305. These authors analyzed the light curves by using statistical methods based on conditional probability distributions.

Wyithe et al. 2000 focused mainly on computing the probability of occurrence of brightness bursts with observed parameters and on estimating the possibility of subsequent bursts. The conclusion that there is an additional overlooked brightness burst in image C associated with caustic crossing between the observing seasons of 1997 and 1998 seems of considerable interest. Such an event can account for the difficulties of modeling the light curve for image C in terms of the model with a single caustic crossing. Wiythe et al. 2000 also expect another caustic crossing in image C 500 days after the 1999 summer maximum (with a large uncertainty of $`1002000`$ days, though).

The amplification in image A in late 1999 is interpreted as caustic crossing in the negative direction, just as we did here. However, the peak of image C in mid-1999 is assumed to be caused by the passage of the source near the caustic cusp. The choice between the two interpretations of image C variability could be made by analyzing color variations of the source, which are much larger during caustic crossing than during cusp passage. A color analysis could be performed by invoking additional data of the simultaneous monitorings at the Apache Point and Maidanak Observatories (Bliokh et al. 1999) through various filters during the same observing period.
Acknowledgments. I wish to thank the OGLE group for organizing and carrying out the monitoring of Q2237+0305 and for providing access to the data before their publication.

## REFERENCES

Agol E. and Krolik J., Astrophys. J. 524, 49 (1999)
Blandford R. D. and Narayan R., Astrophys. J. 310, 568 (1986)
Bliokh P. V., Dudinov V. N., Vakulik V. G. et al., Kinematika Fiz. Nebesnykh Tel 15, 338 (1999)
Chang K. and Refsdal S., Astron. Astrophys. 132, 168 (1984)
Corrigan R. T., Irwin M. J., Arnaud J. et al., Astron. J. 102, 34 (1991)
Fluke C. J. and Webster R. L., Mon. Not. Roy. Astron. Soc. 302, 68 (1999)
Grieger B., Kayser R., and Schramm T., Astron. Astrophys. 252, 508 (1991)
Irwin M. J., Webster R. L., Hewett P. C. et al., Astron. J. 98, 1989 (1989)
Lewis G. F. and Irwin M. J., Mon. Not. Roy. Astron. Soc. 283, 225 (1996)
Manmoto T., Mineshige S., and Kusunose M., Astrophys. J. 489, 791 (1997)
Mineshige S. and Yonehara A., Publ. Astron. Soc. Jpn. 51, 497 (1999)
Østensen R., Refsdal S., Stabell R. et al., Astron. Astrophys. 309, 59 (1996)
Press W, Teukolsky S., Vetterling W. et al. Numerical Recipes in FORTRAN (Cambridge Univ. Press, Cambridge, 1992, 2nd ed.), p. 675
Schneider P. and Wagoner R. V., Astrophys. J. 314, 154 (1987)
Schneider P. and Weiss A., Astron. Astrophys. 171, 49 (1987)
Shakura N. I. and Sunyaev R. A., Astron. Astrophys. 24, 337 (1973)
Wambsganss J., Witt H. J., and Schneider P., Astron. Astrophys. 258, 591 (1992)
Wozniak P. R., Alard C., Udalski A. et al., Astrophys. J. 529, 88 (2000)
Wyithe J. S. B., Webster R. L., and Turner E. L. , Mon. Not. Roy. Astron. Soc. 309, 261 (1999)
Wyithe J. S. B., Turner E. L., and Webster R. L. , astro-ph/0001307 (2000)
Zakharov A. F., Gravitational Lenses and Microlenses (Yanus, Moscow, 1997)
Zakharov A. F. and Sazhin M. V., Usp. Fiz. Nauk 168, 1041 (1998) (Phys. Usp. 41, 945 (1998))