text:
# Dynamic Soft Elasticity in Monodomain Nematic Elastomers

## 1 Introduction

For the last two decades, thermotropic liquid crystalline polymers (LCP) have been a subject of increasing activity due to their fascinating potential in technical applications, as well as some fundamentally different properties in comparison to ordinary liquid crystals. Typical LCP consist of a polymer backbone and mesogenic groups. The backbone is usually made of siloxane or acrylic chains and the mesogenic groups are either grafted (end-on or side-on) using a flexible spacer, thus making a side-chain LCP, or incorporated directly into the polymer backbone (main-chain LCP). When the polymer is chemically crosslinked, it forms a three-dimensional permanent network connected with mesogenic moieties in what is called Liquid Crystalline Elastomers (LCE) fink81 , or gels if the network is swollen by a solvent. Many new physical results have been reported for LCE, but perhaps the most important novelty, unique for nematic LCE materials, is the effect of soft elasticity, where imposed strains could be completely accommodated without any (or with very low) elastic response due to the relaxation of internal microstructure. These unusual mechanical properties of LCE can be directly applied in technology, using LCE as anomalous damping materials, as well as artificial muscles or mechanical transducers.

The early research of LCP rheology, focused mostly on the viscosity of polymers, was summarized in wissbrun . Many interesting results have come to light once the mesogenic groups were incorporated into linear or branched polymer systems. The classical Leslie-Ericksen theory of low-molecular-weight liquid crystals is, in most of the cases, not applicable to LCP systems because of the underlying polymer dynamics making the response distinctly non-Newtonian. Hanus et al. hanus  measured dynamic-mechanical properties of a weakly crosslinked main-chain LCP having crystalline, smectic B, smectic A and isotropic phases. They only observed a clear step in the linear shear modulus $`G^{}`$ at the transition temperature from smectic A to the isotropic state. In the study of LCP having smectic C, smectic A, nematic and isotropic phases pakula , Pakula reported that decreases in the linear modulus were detected at the smectic/nematic transition as well as the glass/smectic transition. Also interestingly, although it was not discussed theoretically at that time, the increase of the modulus could be seen during the transition from the nematic to the isotropic phase. The issues of time-temperature superposition in LCP rheology have been actively discussed by several groups such as Colby et al. colby , Gallani et al. gallani2 ; gallani96 , Fourmaux-Demange et al. fourmaux-demange  and Weilepp et al. weilepp2 . All of these works also reported the characteristic power-law frequency dependence of the storage modulus $`G^{}(\omega )`$. Colby et al. used a sample possessing three phases (isotropic, nematic and smectic A) and concluded that the viscoelastic response is very sensitive to smectic-nematic and smectic-isotropic transitions, but insensitive to the nematic-isotropic transition. The time-temperature superposition worked only at the nematic-isotropic transition. Approaching the glass transition, $`G^{}(\omega )`$ obeys a power law of $`G^{}\omega ^{0.75}`$. The authors suggested that this power law might be connected to the large polydispersity of their sample colby , which was also argued by Fourmaux-Demange et al. fourmaux-demange  for their different LCP system. Gallani et al. measured dynamic-mechanical response of several LCP materials and also found a power law behaviour of the storage modulus $`G^{}`$, with the exponent of $`0.44`$, $`0.6`$, $`0.67`$ or $`0.75`$ depending on the sample. Here again, the time-temperature superposition worked only during the nematic-isotropic transitions. This was attributed fourmaux-demange ; gallani96  to the transient elastic clusters of macroscopic size (a few tens of micrometers), floating in a viscous medium. One may recall that a classical Rouse dynamics of polymer chains would lead to a square-root dependence $`G^{}\omega ^{0.5}`$ Doi .

In this paper we shall find that the power law dependence is also effective in crosslinked LCE materials. In their parallel study of LCE, Gallani et al. and Weilepp et al. gallani96 ; weilepp2  also clearly stated that the rheological response did not change significantly when crossing an isotropic-nematic phase boundary, whereas there has been a crucial change at the nematic-smectic A transition. The power-law frequency dependence of the elastic moduli has again been reported, with $`G^{}\omega ^{0.3}`$ in the smectic A phase. In the isotropic phase, the exponent was found to be 0.5, similar to the Rouse model prediction. A more detailed recent study of nematic LCE rheology by Stein et al. stein01  also reported $`G^{}\omega ^{0.5}`$ in a polydomain nematic elastomers, but could not obtain Master Curves for monodomain, aligned nematic materials. Working with such monodomain nematic rubbers, Clarke et al. clarkePRL01 ; clarke\_hotta2  reported a significant difference in response between polydomain samples (where the nematic director is disordered on the scale above few microns) and monodomain nematic rubbers (with the director is macroscopically aligned along the shear direction, or perpendicular to it). A substantial decrease in the low-frequency storage modulus $`G^{}`$ was registered as the temperature was decreasing through nematic-isotropic transition. This unusual effect has been theoretically described by a continuum theory of linear viscoelastic response in oriented monodomain nematic elastomers terentjev01 .

Continuous stress relaxation is a related subject, which offers a different angle of approach to the complex problem of polymer dynamics exploring the time, rather than the frequency domain. Chasset and Thirion, in their classical paper chasset , described a power law decay in ordinary crosslinked polymer networks such as the isoprene rubber (IR) and the butadiene-styrene copolymer rubber (SBR). They found that, after application of a constant strain, the stress response decreases with time, reaching at long times a characteristic power-law regime $`1/t^{0.10.15}`$ at very long times. In the linear regime (at small enough deformations), in an incompressible material (when the extension modulus is simply three times the shear modulus $`G`$), this translates directly into the relaxation of the effective modulus $`G(t)`$. Many subsequent experiments on stress relaxation, performed by different groups, have supported this long-time power law decay with a small exponent. Theoretical analysis of curro\_pincus , applying the reptation ideas for a network with long free dangling chain ends, has shown that the retraction of these chains could account for the observed stress relaxation law $`\sigma t^a`$ with a small exponent $`a0.10.15`$.

In liquid crystalline networks, the early results shown a different response. The problem in such experiments is always the need to access very long times, when the main transient processes in the polymer network have relaxed and only the slowest mode remains. A peculiar dual regime of long-time relaxation in a polydomain siloxane nematic LCE has been reported by Clarke et al. clarkePRL . A characteristic time $`t^{}3000`$ s separated a regime of fast power-law relaxation, $`Gt^{0.5}`$ or higher, and a very slow relaxation at the later stage, fitted in clarkePRL  with an inverse-logarithmic law $`\sigma [1+\mathrm{ln}(t/t^{})]^1`$. An even faster power-law decay was reported by Hotta et al. hotta\_terentjev2  in their study of acrylic nematic LCE, $`Gt^{0.67}`$ in the early stage of stress relaxation. At very long times, they find the different slow power law decay with the exponent $``$0.15, similar to the one observed in natural crosslinked rubber by Chasset and Thirion.

From this brief summary, it is clear that the power-law dependence of linear response function on time or frequency is a recurrent feature of LCE rheology. On the other hand, it is equally clear that the response is non-universal and the particular features of the response depend on the material composition and also on the nematic director texture and orientation with respect to the applied strain. In this paper we describe a simultaneous study of dynamic-mechanical and stress-relaxation responses for the two well-characterised nematic LCE materials. One sample is a “classical” siloxane side-chain nematic rubber first prepared and described in great detail by Finkelmann et al. and others kupfer ; kupfer2 ; stein01 ; clarkePRL01 . The other material had a high content of main-chain nematic polymer strands, which makes the effective chain anisotropy very high Loft , but also leads to a very slow dynamics. We shall examine in some detail the issues of time-temperature superposition, which is necessary to extrapolate the response function into the regions of time or frequency not practically accessible on experiment. We find that the results of two such different experimental techniques are, in fact, in a good agreement: the observed frequency power laws $`G^{}\omega ^a`$ are matched by the corresponding relaxation laws $`Gt^a`$, if the corresponding ranges of time and frequency are examined.

## 2 Material Characterisation

Two nematic elastomers were prepared in our group following the synthetic procedure pioneered by Finkelmann et al. kupfer . The polymer backbone was a standard poly-methylhydrosiloxane with approximately 60 $`Si`$-$`H`$ units per chain, obtained from ACROS Chemicals. The pendant mesogenic groups were 4’-methoxy phenyl-4-(1-buteneoxy) benzoate (MBB), as illustrated in Fig. 1, attached to the backbone end-on via the hydrosilation reaction in the presence of a commercial platinum catalyst COD, obtained from Wacker Chemie. The first monodomain nematic elastomer was prepared following the two-step crosslinking reaction kupfer , with 1,4 di(11-undecene) benzene (11UB), a small flexible difunctional molecule deemed to have relatively minor effect on the overall mesogenic properties of the liquid crystalline polymer (synthesized in the house).

The second network was prepared following the similar procedure, but using di-vinyl terminated polymer chains of $`\alpha `$-{4-\[1-(4’-{11-undecenyloxy} biphenyl)-2-phenyl\] butyl)-$`\omega `$-(11- undecenyloxy) poly-\[1-(4- oxydecamethyleneoxy)- biphenyl -2-phenyl\] butyl (MC) that themselves make a main-chain LCP percec . In our case, the crosslinking chains had $``$75 rod-like monomers between the terminal vinyl groups (determined by GPC, polydispersity $``$2). As a result, the properties of a nematic LCE are dominated by these long main-chain LCP strands connecting the siloxane side-chain molecules (as the Table 1 indicates, the overall MC content is $``$58%). In both cases, the crosslinking density was kept at 10 mol% of the reacting bonds in the siloxane backbone, so that on average, each siloxane chain has nine pendant mesogenic groups between crosslinking sites.

Both elastomers were prepared from the same batch of backbone and side-group mesogenic materials, with the same relative concentration of crosslinking groups.

Table 1 summarizes chemical and thermal characteristics of both materials. Equilibrium transition temperatures were determined on a Perkin Elmer Pyris 1 differential scanning calorimeter (DSC), extrapolating to low cooling rates, and the nematic phase identified by polarizing optical microscopy and X-ray scattering. The average degree of chain anisotropy was identified by observing the uniaxial thermal expansion of LCE samples as a function of temperature. Thermal expansion measurements were made by suspending the samples, without load, in a glass-fronted oven and measuring the variation in natural length of the samples with temperature, $`L(T)`$; see Loft  for further details.

## 3 Dynamic Mechanical Analysis

In this paper, we used the dynamic mechanical apparatus (DMA) Viscoanaliser VA4000, from Metravib RDS. In the previous work clarke\_hotta2 ; clarkePRL01  we have already examined the anisotropy of the shear response depending on the orientation of the nematic director, see the sketch in Fig. 2(a). Only the rotation-inducing geometry reflects the effects of soft elasticity and here we concentrate on this particular case. The way in which the samples were mounted in the DMA is shown in Fig. 2(b). It was important to choose the correct mounting for a given temperature range because of the great variation in the magnitude of response between the low-temperature (or high-frequency) glassy regime and the high-temperature rubber plateau. We used two different ways of sample mounting: The more stiff “sandwich” geometry, with each sample having the large area of contact with the clamps and the small thickness, was used to collect weaker signals from softened samples at higher temperature (in our DMA setup the thin-film sample had a round shape). The “tape” mounting geometry, in which the rectangularly cut thin-film sample is clamped sideways with only the small cross-section area, results in smaller forces and allows to measure the modulus of very rigid samples up to the glass plateau. The Table 2 gives the dimensions of typical sample mounting.

A geometry of uniaxial extension, more commonly used in studies of equilibrium stress-strain relaxation in rubbers, is less appropriate for an oscillating regime because of slow relaxation and incomplete sample recovery on each cycle. In the chosen simple shear setup, the middle part of the sample holder is forced to move in an oscillating fashion, maintaining the small shear strain – between $`0.003`$ and $`0.004`$ (that is, 0.3-0.4%). This value was chosen from a separate strain-ramping experiment, finding the range of strains that still return a linear response in the whole of the chosen frequency range.

### 3.1 Temperature Scans

The materials were studied over a wide temperature range encompassing the glass and the nematic/isotropic transitions. Typically the temperature scans ranged from $`50^\mathrm{o}C`$ to $`130^\mathrm{o}C`$, at a number of fixed frequency values. Figure 3 shows the results for the storage modulus $`G^{}`$ at just two frequencies, 0.01 and 100 Hz, for the SiH sample. The storage modulus $`G^{}`$, which is the real part of the linear complex modulus $`G^{}(\omega )`$ ferry  and the loss factor $`\mathrm{tan}\delta =G^{\prime \prime }/G^{}`$ (with $`G^{\prime \prime }`$ the imaginary part of $`G^{}(\omega )`$ – the loss modulus) were plotted against the temperature. As we mentioned above, it is technically difficult to obtain a clean mechanical signal from the sample that undergoes a change in stiffness of up to $`10^410^9`$Pa. Therefore, we used the two ways of mounting the samples, the “sandwich” at high and the “tape” at low temperatures. One can see in Fig. 3, as well as in Fig. 4 for the SiMC material, that the two portions of the full curve, for the tape and sandwich modes, are connecting quite successfully without any additional adjustment (the data were discarded if the vertical error was reaching $`15\%`$, on the assumption of a mounting fault). This indicates good reliability of DMA measurements for both the storage modulus $`G^{}`$ and loss factor $`\mathrm{tan}\delta `$.

The glass plateau modulus for both SiH and SiMC materials reaches up to nearly $`10^9`$ Pa, which is a reasonable value in polymeric materials. Between $`0^\mathrm{o}C`$ and $`10^\mathrm{o}C`$, both SiH and SiMC undergo the glass/nematic transition, showing a dramatic decrease in $`G^{}`$ on cooling, and the associated steep rise in $`\mathrm{tan}\delta `$. The storage modulus decreases down to $`10^5`$ Pa which is an expected order of magnitude in rubbery phases. However, both materials show a substantial drop in $`G^{}`$ in the nematic phase, below $`T_{ni}`$. This effect, most pronounced at low frequencies, is the dynamic soft elasticity: the reduction in the modulus and the rise in internal mechanical dissipation caused by the underlying director rotations induced by shear terentjev01 ; clarkePRL01 . The relative depth of the drop in $`G^{}`$ is higher in SiMC, the material with a higher chain anisotropy due to its main-chain LCP component. Above $`T_{ni}`$, the storage modulus returns to its (higher) value characteristic of ordinary isotropic rubbers. At very low frequencies $`G^{}`$ is slightly rising with temperature, revealing that the sample is in a crosslinked rubbery state with its equilibrium shear modulus $`Gn_xk_BT`$ (where $`n_x`$ is the crosslink density and $`k_B`$ is Boltzmann factor).

### 3.2 Frequency Scans

Figures 5 and 6 show the same data for the storage modulus $`G^{}`$, but this time plotted against frequency at a set of constant values of temperature. In our experiments, the frequency was changing between 0.01 to 900 Hz, measured 45 points in total. However, we were obliged to discard higher frequencies: the DMA data sets must pass a series of tests on internal consistency and failing these tests indicates one of several possible problems with acquisition (most common is the internal resonance of the mechanical frame). Only the reliable range of results is presented in Figs. 5 and 6.

The range of storage modulus variation with frequency matches well with the temperature scans, from $`10^9`$ Pa down to $`10^410^5`$ Pa. As in the temperature scans, the average value of $`G^{}`$ is rising on heating, when each sample approaches its isotropic rubbery phase and the mesogenic effects weaken. In each of the figures, the plot (b) shows the expanded version of the high-temperature results, clearly indicating the effect of dynamic soft elasticity. In addition, one can see for the SiMC material that the low-frequency rubber plateau is not achieved even at highest temperatures: the modulus continues to decrease with decreasing frequency indicating the on-going mechanical relaxation processes in the sheared rubbery network. It is in contrast with the SiH results, where one can clearly identify the frequency-independent rubber plateau in both the nematic and the isotropic phases.

Below we shall discuss this extremely slow relaxation of the MC-containing elastomers and match these results with the findings of static stress relaxation. We attribute this slowing down to the hairpin folds of the backbone of rods with flexible spacers, which has been demonstrated elias  to freeze the dynamics even of non-crosslinked polymer melts.

### 3.3 Time-Temperature Superposition

The frequency scans presented in Figures 5(a) and 6(a) invite the time-temperature superposition ferry . At the same time, one immediately discovers from the corresponding (b)-plots that such superposition will not be possible as the monodomain nematic rubber sheared in the “soft” $`𝑫`$-geometry undergoes the nematic transition. We assert that this conclusion can only be made for materials which show a clean and pronounced soft decrease in $`G^{}`$: in a frustrated polydomain elastomer, or in an overcrosslinked network, the required nematic director rotation may be impeded by internal constraints and, therefore, the characteristic nematic soft elasticity masked. In that case, one might expect the ordinary time-temperature superposition to be valid, same as in ordinary isotropic polymer systems.

Time-temperature (or, in this case, frequency-temperature) superposition is a very interesting, practically useful, but also theoretically ill-defined approach to data analysis in dynamic mechanical measurements. The empirical algorithm proposed by Williams, Landel and Ferry (WLF) ferry  allows extrapolating the results into the time or frequency ranges far beyond experimentally accessible. As traditional in this field, we attempted such a superposition of our data for nematic elastomers. The storage moduli $`G^{}`$ taken in the series of frequency scans at different temperatures are shifted along the frequency axis until each data sets superimposes with the previous (with no additional vertical shifting). The multiplicative shift factor, $`\omega ^{}=a_T\omega `$, is strongly temperature dependent. Another important factor is the choice of reference temperature $`T_{\mathrm{ref}}`$ at which the real frequency values are taken. WLF suggested using the glass transition temperature as $`T_{\mathrm{ref}}`$, however, we decided against this for the following reasons. First of all, the presented results show that the dynamic glass transition is, in fact, difficult to pinpoint as the transformation occurs over a range of temperatures. Secondly, no additional information is extracted from the choice of $`T_{\mathrm{ref}}`$ anyway; all it affects is the overall shift of the Master Curve along the frequency axis. Our choice of $`T_{\mathrm{ref}}`$ as the ambient room temperature was dictated by the possibility of determining the key parameter of polymer chains, the Rouse frequency $`\omega _\mathrm{R}=\tau _\mathrm{R}^1`$ at the relevant $`T_{\mathrm{ref}}`$, from the position of the dynamic glass transition on the frequency axis. In both materials, we find the Rouse time $`\tau _\mathrm{R}10^510^6`$s (which of course must be regarded as a very crude estimate).

Figures 7(a) and 8(a) show the results of superposition for the SiH and SiMC samples, respectively, in the form of the Master Curves. The procedure is only applicable at the temperatures below the nematic transition $`T_{ni}`$. Each graph also shows one data set $`G^{}(\omega )`$, appropriately shifted, at a high temperature, in the isotropic phase. The values are much higher than on the low-frequency nematic plateau and, of course, no shifting along the $`\omega `$-axis would superimpose them.

The very existence of Master Curves indicates that there is a unique physical process of mechanical relaxation, which spans the range of times, or frequencies that is far beyond experimentally accessible. Our measurements, at each temperature, are only able to detect a small portion of the whole process, but the time-temperature superposition suggests a way of building the whole picture from small pieces. An important part of this analysis is the temperature dependence of shift factors $`a_T`$, which is plotted in Figs. 7(b) and 8(b). The empirical WLF relationship is based on the ideas of glass transition and has been successful in describing the properties of simple homopolymer melts:

$$\mathrm{log}a_T=\frac{C_1(TT_{\mathrm{ref}})}{C_2+(TT_{\mathrm{ref}})},$$
(1)
where $`C_1`$ and $`C_2`$ are the constants, called the WLF coefficients. (In more complex polymers, and in crosslinked networks, deviations from this equation have been substantial). In our case, Figures 7(b) and 8(b), the shift factors follow the WLF relation very well at lower temperatures (below $`T_{\mathrm{ref}}`$), with the WLF coefficients $`C_117.4`$ and $`C_2102`$ for both materials,<sup>1</sup><sup>1</sup>1Recall that both SiH and SiMC have been analysed with nearly the same reference temperature, around $`20^\mathrm{o}`$C. very close to the values found in ordinary isotropic polymer systems. This is perhaps not surprising when one recalls that the empirical WLF relationship is “designed” to describe the dynamic glass transition region.

However, the shift factors deviate from the classical WLF behavior very strongly as soon as the higher temperatures (and, correspondingly, long relaxation times) are considered. Here, the two materials behave quite differently: SiH only allows the superposition of frequency scans with much smaller shifts $`a_T`$, while the Master Curve for SiMC demands much greater shifts. One should note that, in practice, the procedure of shifting of plots along the frequency axis to achieve superposition is very sensitive to the value of each shift factor $`a_T`$ and, accordingly, the errors in Figs. 7(b) and 8(b) are very small. It is more likely that the noise in the $`a_T`$ data is due to uncertainties in DMA temperature or modulus measurements.

Finally, we note that the region of the dynamic glass transition, in both Master Curves for SiH and SiMC, follows a clear power law $`G^{}\omega ^a`$. The plots indicate the non-universal values of the exponent, similar to many previous observations. The very low frequency (long-time) tail of continuing relaxation in SiMC may appear to follow another power law with a much smaller exponent of $`0.07`$ also labelled in Fig. 8(a). However, we are inclined to believe that this is an ambiguity of the complex transient regime: in the previous studies of static stress relaxation clarkePRL ; hotta\_terentjev2  we were distinctly unable to distinguish between a very low power-law exponent and an even slower logarithmic relaxation regime. Since the material never reached the equilibrium, whether in a real experiment or in extrapolation, the identification of the longest relaxation regime remained impossible.

## 4 Static Stress Relaxation

Stress relaxation curves were measured at several different temperatures at a constant extensional strain of $`0.09`$ for both SiH and SiMC, following the simple technique described in clarkePRL ; hotta\_terentjev2 . This strain is higher than the shear strain used in DMA experiments (for technical reasons we simply could not impose a smaller extension) and the question of whether the response is still linear may be raised. We have to stress two points. First of all, the strain of $``$9% is, in fact, on the border of the initial linear regime of stress-strain curves measured for both SiH and SiMC, cf. clarke\_hotta1 . This, however, has to be regarded with caution because it is impossible to obtain a truly equilibrium stress-strain relation (a certain rate of strain is always involved) and any extrapolation to zero strain rate would involve a theoretical model, which may or may not be ultimately successful. However, as a second point, we have observed that the key conclusion of this section, about the power laws of stress relaxation, remains valid even for higher strains – on and above the soft plateau when we are certain that the stress-strain is non-linear. Therefore, in practice, the choice of the strain step for the static relaxation experiments was not very important.

For the constant-strain extension setup, each sample was cut with a razor blade into a rectangular shape of approximately $`6\times 2.5\times 0.2\text{mm}`$, measured by a micrometer (narrow elongated strips minimize the edge effects during the extensional deformation). The initial sample cross-section area of $`0.5\text{mm}^2`$ was used to calculate the nominal stress $`\sigma `$, in units of (mN/mm$`{}_{}{}^{2}`$ kPa), from the force measured by the temperature-compensated dynamometer (UF1/AD20 from Pioden Controls Ltd). The values of force, obtained in arbitrary units, were calibrated with weights at different experimental temperatures of 30, 50, 80, 105$`{}_{}{}^{\mathrm{o}}C`$. The samples were then extended with the engineering strain $`\epsilon =(LL_0)/L_0`$, increasing from 0 to 0.09 within a second. Both monodomain samples, SiH and SiMC, were stretched perpendicular to their uniform director orientation so that a soft-deformation geometry, comparable to the $`𝑫`$-shear geometry of DMA experiments, is achieved. In each case the stress relaxation was measured for two days $`(10^5\text{s})`$. Figure 9 shows typical stress relaxation curves taken in the series of these measurements. The data were then fitted with several possible stress relaxation equations, such as the power-law, the inverse logarithmic law and the exponentials.

In Fig. 9(a), for the “faster” SiH material, one can note that the stress relaxation was substantial at a lower temperature of 31$`{}_{}{}^{\mathrm{o}}C`$ and is almost unnoticeable at higher temperatures (in fact, the data set for 82$`{}_{}{}^{\mathrm{o}}C`$ is not shown because it is hard to distinguish from the 52$`{}_{}{}^{\mathrm{o}}C`$ case). The relaxation in SiMC continues even at high temperatures, but its extent is clearly smaller.

Instead of presenting and discussing the individual relaxation fits and their relations, we shall follow a different route of analysis. If we assume that the Master Curves, shown in Figs. 7(a) and 8(a), are genuine and represent the complete relaxation dynamics in our complex polymeric systems, then we should expect that the same laws must govern the static stress relaxation. Let us plot such a Master Curve against the inverse frequency, which would be a measure of real time. Such plots, Figs. 10(a) and 11(a) look reasonable, as one would expect the relaxing modulus to behave in real time, in double-logarithmic scale. Of course, one cannot assign anything more than a qualitative meaning to this representation because, let us not forget, the $`G^{}(\omega )`$ and $`G(tt^{})`$ are related by the Fourier transformation and merely mapping $`G(\omega )`$ onto the inverse frequency is not a rigorous operation. It does, however work for the power-law dependence: from the basic definitions

$`\sigma (t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^t}G(tt^{}){\displaystyle \frac{d}{dt^{}}}\epsilon (t^{})𝑑t^{};`$ (2)
$`\sigma (\omega )`$ $`=`$ $`\mathrm{i}\omega G^{}(\omega )\epsilon (\omega ),`$
and assuming the linear response function to decay as a power law, at long times $`G(t)G_0t^a`$, one can easily calculate the corresponding estimate of the complex shear modulus:

$`G^{}(\omega )G_0\mathrm{\Gamma }(1a)\left(e^{\frac{1}{2}\mathrm{i}a\pi }e^{\frac{3}{2}\mathrm{i}a\pi }\right)\omega ^a.`$ (3)
Apart from a constant coefficient, this is directly the inverse function of $`G(t)`$.

Let us now take the real-time stress relaxation data sets and attempt the time-temperature superposition on them, building a Master Curve by shifting along the time-axis. This operation works for all temperatures, but we only present three data sets in each of Figs. 10(b) and 11(b), to allow a clear view of the underlying line, which is the DMA Master Curve for inverse frequencies. The matching is quite remarkable, especially if we take a look at the shift factors (now for the time-shifts, and so with the opposite sign!) required to bring the relaxation data onto this Master Curve. Figures 7(b) and 8(b) show the frequency shift factors and their WLF analysis, but also the square symbols that represent the (inverse) time shift factors for the static relaxation data. (Recall that the Master Curves were built for the reference temperature of around 20$`{}_{}{}^{\mathrm{o}}C`$). The agreement and the trend in these respective shift factors is not bad.

## 5 Discussion

First of all, we must recognise the difference of our results from some literature data. Gallani et al. and Weilepp et al. gallani2 ; weilepp2  were both pointing out that the nematic-isotropic transition is trivial for the dynamic mechanical measurements. The reason could be that polydomain nematic elastomers were used, although we have earlier seen the effect of dynamic soft elasticity in a variety of polydomain elastomers clarke01 . The drop in $`G^{}`$ below $`T_{ni}`$ is smaller in polydomain samples, but it only disappears in the monodomain $`𝑽`$-geometry, as in fact the more recent study stein01  confirms. It is also interesting that many authors, e.g. gallani2 ; fourmaux-demange , reported their success in time-temperature superposition across the nematic-isotropic phase transition, which of course is related to the absence of pronounced soft drops in $`G^{}`$ values. We can only suppose that a different synthetic procedure of their materials led to the very non-soft (highly semi-soft) nematic networks.

Returning to our own results, the two main results of this work are the building of Master Curves and identifying the temperature regimes where the superposition follows the WLF and where it fails altogether, and establishing the direct analogy between the time and the frequency domains. In the first part, the most important observation is the substantial difference in rubber-plateau (low-frequency) moduli between the nematic and isotropic phase, in the shear geometry when the director rotation is induced. This effect of dynamic soft elasticity can lead to a number of technological advances, e.g., in the area of mechanical vibration damping clarke01  (exploiting the record-high values of the loss factor over a broad range of temperatures and frequencies), or even suggesting a novel field of polarised acoustics (using the fact that only certain, soft, shear-wave geometries are strongly attenuated).

Comparing the results of static stress relaxation with those of dynamic mechanical measurements is intriguing. Regarding the ability to quantitatively predict the static relaxation behaviour from the independently (and relatively quickly) acquired DMA data, or vice versa, by inverting the Master Curve at a specific reference temperature, we have to remain cautious. Although the possibility is highly attractive and the procedure has worked spectacularly well in our case, one needs to verify it on different materials, including a check of how the classical isotropic rubber behaves in this respect.

In our case, the Master Curve inversion has been successful and predicts that, for instance, a “traditional” kupfer ; kupfer2  siloxane side-chain nematic elastomer SiH, crosslinked in a way that allows the soft elasticity, would be able to relax under a constant strain and reach its mechanical equilibrium (at least approximately) after some $`10^4`$s at room temperature, with a very little residual relaxation after that. This is somewhat different from the results clarkePRL  on a very similar, but polydomain nematic elastomer, where the relaxation have not reached its end even after extrapolation to astronomical times. Analogous continuing relaxation was reported in other, chemically very different polydomain nematic rubbers hotta\_terentjev2 . The effect of internal mechanical constraints imposed by the domain interfaces is clearly dominating long-time regime, but is absent in our present monodomain network. In contrast, the composite nematic elastomer SiMC, albeit monodomain, also shows a non-ending relaxation. We assume, by analogy with the previous work on relaxation in non-crosslinked MC polymer elias , that internal constraints on “free” stress relaxation are now imposed by the hairpin folds of long nematic main-chain polymer strands.

In summary – dynamic-mechanical analysis have been performed using two different monodomain nematic elastomers having different degrees of anisotropy and internal microstructure. Master Curves have been built between the glassy state and the nematic-isotropic phase transition, where the modulus saturates on the low-level soft plateau. Above $`T_{ni}`$ the modulus rises substantially, since the internal relaxation is no longer able to reduce the elastic response – and the further time-temperature superposition fails. The dynamics of these elastomers is dominated by power laws (non-universal, depending on the material and temperature), which is confirmed by the successful operation of the Master Curve inversion to describe the static stress relaxation. Power laws (in time) of stress relaxation match very well with the corresponding frequency dependence of the dynamic modulus in the appropriate range of temperatures.

### Acknowledgements

We thanks EPSRC and Bridgestone Corporation for funding this work, and Wacker Chemie for the donation of the platinum catalyst. We are indebted to H. Finkelmann for the advice and assistance in chemical synthesis and to A.R. Tajbakhsh for the preparation of monodomain LCE materials.


format: text

text:
# Introduction

## Introduction

Ever since numerical lattice QCD computations have been done, the spectrum of light mesons has served as a benchmark problem. This is still true today, as actions respecting the global chiral symmetry among the light flavours are being tested and several groups have embarked on ambitious simulations of “full” QCD with two or possibly more dynamical flavours. These developments address one of the key issues in low-energy QCD. The fact that chiral symmetry is both spontaneously and explicitly broken generates pseudo-Goldstone bosons, i.e. particles that dominate (for small enough quark masses) the long-range behaviour of correlators between external currents and which are collectively called “pions”.

The lattice is not the only framework to address the low-energy structure of QCD. In the old days PCAC relations were exploited to predict the dependence of low-energy observables on shifts in the quark masses and external momenta. The one best known is

$$F_\pi ^2M_\pi ^2=(m_\mathrm{u}+m_\mathrm{d})|0|\overline{q}q|0|+O(m^2)$$
(1)
which connects the pion mass and decay constant to the product of the explicit and spontaneous symmetry breaking parameters. However, the Gell-Mann – Oakes – Renner relation (1) does not give a prediction how $`F_\pi `$ and $`M_\pi `$ separately depend on the quark mass. And in the real world the latter may be shifted only by a discrete amount (e.g. by replacing $`ds,\pi K`$ one gets a leading order prediction for the quark mass ratio $`m_s/m`$ with $`m(m_u+m_d)/2`$).

Today, the first limitation is overcome, since the successor of PCAC, Chiral Perturbation Theory , gives detailed predictions how either $`M_\pi `$ or $`F_\pi `$ individually depends on the quark mass $`m`$. And – in principle – the second problem is gone since one may vary $`m`$ continuously in a lattice computation. Therefore, it seems natural to combine the two approaches to benefit from their respective advantages. However, for that aim quarks need to be taken sufficiently light, and this is a numerical challenge on the lattice. Below, an elementary test is presented whether this is already achieved in present state-of-the-art results for $`M_\pi ^2`$ versus $`m`$, as published by the CP-PACS  and UKQCD  collaborations. The idea is to restrict the analysis to the “doubly degenerate” case with $`N_f=2`$ $`O(a)`$ improved Wilson quarks, i.e. to consider only the subset where both valence-quarks are exactly degenerate with the sea-quarks, as this avoids additional assumptions of the “partially quenched” framework.

## Prediction by Chiral Perturbation Theory

Here, I give a brief summary of the work by Gasser and Leutwyler (GL), with an application to lattice data in mind.

GL have calculated the substitute for the GOR relation (1) to NLO in Chiral Perturbation Theory (XPT) with $`N_f=2`$ quarks . To that order, two low-energy constants111Unlike $`F`$, the other LO “constant” $`B`$ is scheme- and scale-dependent, as is $`m`$. In the chiral representation, only the product $`mB`$ appears which is RG-invariant (cf. (5)). In order to determine $`B`$ and $`m`$ separately, a lattice computation is needed; the result of the CP-PACS study  along with (5) for the physical pion is
$$m_{\mathrm{phys}}(\overline{\mathrm{MS}},\mu 2\mathrm{GeV})3.5\pm 0.2\mathrm{MeV},B(\overline{\mathrm{MS}},\mu 2\mathrm{GeV})2.80.15\mathrm{GeV}.$$
(2)  from the LO Lagrangian appear ($`F=lim_{m0}F_\pi ,B=lim_{m0}0|\overline{q}q|0/F^2`$) together with the finite parts of the $`l_i`$ in the NLO Lagrangian. Restricted to the degenerate case their result  takes the form

$`M_\pi ^2`$ $`=`$ $`M^2+{\displaystyle \frac{1}{F^2}}\left({\displaystyle \frac{M^4}{32\pi ^2}}\mathrm{log}({\displaystyle \frac{M^2}{\mu ^2}})2M^4l_3^\mathrm{r}(\mu )\right)`$ (3)
$`F_\pi `$ $`=`$ $`F+{\displaystyle \frac{1}{F}}\left({\displaystyle \frac{M^2}{16\pi ^2}}\mathrm{log}({\displaystyle \frac{M^2}{\mu ^2}})+M^2l_4^\mathrm{r}(\mu )\right)`$ (4)
$`M^2`$ $``$ $`2mB.`$ (5)
Note that (3, 4) has the typical structure of a NLO prediction: A chiral logarithm summarizing the contribution from the pion loops appears together with a counterterm. The $`l_i^\mathrm{r}`$ are the renormalized GL coefficients, i.e. they are the descendents of the $`l_i`$ which appear in the NLO Lagrangian and which are divergent quantities. As a result, the $`l_i^\mathrm{r}`$ depend on the (chiral) renormalization scale $`\mu `$. In the dimensionally regularized theory one has

$`l_i`$ $`=`$ $`l_i^\mathrm{r}(\mu )+\gamma _i\lambda (\mu )`$ (6)
$`\lambda (\mu )`$ $`=`$ $`{\displaystyle \frac{1}{16\pi ^2\mu ^{4d}}}\left({\displaystyle \frac{1}{4d}}+{\displaystyle \frac{\mathrm{log}(4\pi )+\mathrm{\Gamma }^{}(1)+1}{2}}\right)`$ (7)
$`l_i^\mathrm{r}(\mu )`$ $`=`$ $`l_i^\mathrm{r}(\mu ^{}){\displaystyle \frac{\gamma _i}{16\pi ^2}}\mathrm{log}({\displaystyle \frac{\mu }{\mu ^{}}})`$ (8)
and the $`\beta `$-function coefficients $`\gamma _i`$ are known. In the present context only $`\gamma _3=\frac{1}{2},\gamma _4=2`$  are relevant, and this suggests that one rewrites (3, 4) with the help of

$$l_i^\mathrm{r}=\frac{\gamma _i}{32\pi ^2}\left(\overline{l}_i+\mathrm{log}(\frac{M^2}{\mu ^2})\right),$$
(9)
where the $`\mu `$ dependence in $`l_i^\mathrm{r}=l_i^\mathrm{r}(\mu )`$ is traded for an $`M`$ dependence in $`\overline{l}_i=\overline{l}_i(M)`$, to get

$`M_\pi ^2`$ $`=`$ $`M^2\left(1{\displaystyle \frac{M^2}{32\pi ^2F^2}}\overline{l}_3+O(M^4)\right)`$ (10)
$`F_\pi `$ $`=`$ $`F\left(1+{\displaystyle \frac{M^2}{16\pi ^2F^2}}\overline{l}_4+O(M^4)\right).`$ (11)
The NLO part is given in terms of the LO parameters $`F,B`$ and the NLO coefficients $`\overline{l}_3,\overline{l}_4`$. While the former two are known quite accurately, for the $`\overline{l}_i`$ only their running in $`M^2`$ is known exactly and the phenomenological estimate of the integration constants has comparatively large error-bars. Gasser and Leutwyler give in their initial paper  the estimate

$$\overline{l}_3(m_{\mathrm{phys}})=2.9\pm 2.4,\overline{l}_4(m_{\mathrm{phys}})=4.3\pm 0.9$$
(12)
for real world quark masses. It turns out that even this limited information is useful, since it determines the curvature of $`M_\pi `$ as a function of the quark mass. Close to the chiral limit, both $`\overline{l}_i`$ are positive and as a consequence $`M_\pi ^2`$ does not rise strictly linear in $`m`$ but turns right, while $`F_\pi `$ has a positive first derivative in $`m`$. More specifically, (12) translates into

$$M_{\pi ,\mathrm{phys}}139\mathrm{MeV},F_{\pi ,\mathrm{phys}}92.4\mathrm{MeV}M_{\mathrm{phys}}141\mathrm{MeV},F86.1\mathrm{MeV}$$
(13)
which means that the physical pion is somewhat lighter than it would be, if the LO relation were exact, while is decay constant exceeds its value in the chiral limit by $`7\%`$.

A seemingly formal point which, in the end, proves convenient in analyzing the lattice data is the following. A naive look at (10) suggests that the typical structure of a NLO prediction is gone – rather than a $`M^4`$ and a $`M^4\mathrm{log}(M^2)`$ contribution, only the polynomial part seems left. The point is that this impression is entirely misleading; the IR divergencies (which are genuine to QCD in the chiral limit) are not gone, they are just hidden in the $`\overline{l}_i`$. The situation is, in fact, opposite – the $`M^4`$ part has been eliminated in favour of a pure $`M^4\mathrm{log}(M^2)`$ contribution, and the last step is to make this apparent. The quark mass dependence of $`\overline{l}_3`$ is given through

$$\mathrm{log}(\frac{\mathrm{\Lambda }_i^2}{M^2})=\mathrm{log}(\frac{\mathrm{\Lambda }_i^2}{M_{\mathrm{phys}}^2})+\mathrm{log}(\frac{M_{\mathrm{phys}}^2}{M^2})=\overline{l}_i(m_{\mathrm{phys}})+\mathrm{log}(\frac{M_{\mathrm{phys}}^2}{M^2})=\overline{l}_i(m)$$
(14)
and together with (5), the relation takes its final form (see e.g. )

$`M_\pi ^2`$ $`=`$ $`2mB{\displaystyle \frac{m^2B^2}{8\pi ^2F^2}}\mathrm{log}({\displaystyle \frac{\mathrm{\Lambda }_3^2}{2mB}})+O(m^3)`$ (15)
$`F_\pi `$ $`=`$ $`F+{\displaystyle \frac{mB}{8\pi ^2F}}\mathrm{log}({\displaystyle \frac{\mathrm{\Lambda }_4^2}{2mB}})+O(m^3)`$ (16)
where $`\mathrm{\Lambda }_4,\mathrm{\Lambda }_3`$ represent universal low energy scales. The estimates (12) translate into

$$\mathrm{\Lambda }_3=0.6\mathrm{GeV}\begin{array}{c}+1.4\mathrm{GeV}\hfill \\ 0.4\mathrm{GeV}\hfill \end{array},\mathrm{\Lambda }_4=1.2\mathrm{GeV}\begin{array}{c}+0.7\mathrm{GeV}\hfill \\ 0.4\mathrm{GeV}\hfill \end{array}.$$
(17)
It is worth emphasizing that all four parameters $`F,M^2=2Bm,\mathrm{\Lambda }_3,\mathrm{\Lambda }_4`$ do not depend on the QCD renormalization scheme , i.e. they are the proper physical low-energy parameters at order $`O(p^4)`$. Furthermore the $`\mathrm{\Lambda }_i`$ do not depend on the quark masses, and this means that the representation (15, 16) is perfectly suited to analyze the lattice data even if they are gotten at quark masses larger than $`m(m_\mathrm{u}+m_\mathrm{d})/2`$ in the real world – as long as they are not beyond the regime of validity of the chiral expansion itself.

The latter point is one of the key issues in the comparison we aim at. The chiral expansion is known to be asymptotic, and this means that increasing the order will enhance the accuracy near the chiral limit – at the price of worsening the prediction for heavier masses. What is the “critical scale” beyond which the chiral expansion “explodes” is, a priori, not known. From a formal point of view, one might think that (10, 11) indicate that the expansion is in $`M/(4\pi F)`$ and hence hope that it is good for pion masses up to $`1\mathrm{GeV}`$. In this paper I will argue that watching the convergence pattern at a fixed quark mass gives a more reliable estimate what is the permissible range. This is facilitated since the NNLO expression for $`M_\pi ^2`$ (with $`N_f=2`$) is known . The result reads (for the presentation I follow )

$$M_\pi ^2=2mB\left(1\frac{mB}{16\pi ^2F^2}\mathrm{log}(\frac{\mathrm{\Lambda }_3^2}{M^2})+\frac{m^2B^2}{64\pi ^4F^4}\left\{\frac{17}{8}\left(\mathrm{log}(\frac{\mathrm{\Lambda }_M^2}{M^2})\right)^2+k_M\right\}+O(m^3)\right)$$
(18)
where $`\mathrm{\Lambda }_M`$ is implicitly defined through $`51\mathrm{log}(\frac{\mathrm{\Lambda }_M^2}{\mu ^2})=28\mathrm{log}(\frac{\mathrm{\Lambda }_1^2}{\mu ^2})+32\mathrm{log}(\frac{\mathrm{\Lambda }_2^2}{\mu ^2})9\mathrm{log}(\frac{\mathrm{\Lambda }_3^2}{\mu ^2})+49`$ and the mass-independent $`k_M`$ accounts for the remainder at $`O(p^6)`$, in particular the new counterterms. Phenomenological values for $`\mathrm{\Lambda }_M`$ and $`k_M`$ will be mentioned below.

## Lattice Data

We are now in a position to esteem the results by the CP-PACS and UKQCD collaborations for the quark mass dependence of the pion mass. We shall consider, out of these $`N_f=2`$ data, only the two-fold degenerate case where both valence quarks have the same mass and where they are, at the same time, exactly degenerate with the sea quarks so that the theory is unitary.

The CP-PACS collaboration has simulated various ($`\beta ,\kappa `$) combinations with an RG improved gauge action and a mean-field improved clover quark action . They use a grid of size $`12^3\times 24,16^3\times 32,24^3\times 48,24^3\times 48`$ at $`\beta =1.8,1.95,2.1,2.2`$, respectively, which leads to a lattice spacing, if determined through the $`\rho `$ mass, between $`0.215\mathrm{fm}`$ and $`0.087\mathrm{fm}`$ and hence a spatial box size between $`2.58\mathrm{fm}`$ and $`2.08\mathrm{fm}`$. With so much information at hand, one could, in principle, attempt a continuum extrapolation for $`M_\pi ^2`$ versus the sum of the (degenerate) quark masses, $`2m`$. However, non-perturbative renormalization might be necessary, and/or the lattice spacing might be too large at the lower $`\beta `$ values. For this reason, I have decided to concentrate on the $`\beta =2.1`$ data, since here discretization effects are not supposed to be too large, and good statistics is available. At that $`\beta `$ value, even the lightest pion is unlikely to suffer from finite-size effects, since $`M_\pi L>7`$, and the lattice spacing determined via the $`\rho `$ mass is of order $`0.1\mathrm{fm}`$ and hence comparable with that in the UKQCD simulations.

The UKQCD collaboration works on a $`16^3\times 32`$ grid, using different actions: Wilson glue and non-perturbatively $`O(a)`$ improved clover quarks . Another point in which they differ from CP-PACS is a tactical one: they try to relax $`\beta `$ in pushing $`\kappa _{\mathrm{sea}}`$ up (i.e. the quark mass down) such that the lattice spacing, in units of $`r_0`$, stays constant. Numerically, it is $`a0.1\mathrm{fm}`$ for the data considered below (though the ensemble at $`(\beta ,\kappa _{\mathrm{sea}})=(5.2,0.1355)`$ is not matched any more), and the hope is that the size of discretization effects would be approximately constant. This choice implies that the physical box size stays constant, too, and the bound $`M_\pi L>4.5`$ maintained makes one feel comfortable that finite size effects are small.

The plan of this article is to ignore that the dynamical quark masses might be too heavy for the chiral prediction to be applicable and to ignore that in principle a continuum extrapolation is needed but instead to go ahead and simply compare the two datasets on a “as-is” basis with the LO/NLO/NNLO prediction from Chiral Perturbation Theory in the continuum. With the relevant low-energy constants on the chiral side given in physical units, $`r_0^1`$ must be so, too. In this article, this is done through the assumption that $`r_0^1`$ represents a universal low-energy scale, unaffected by unquenching effects; the numerical value used is $`r_0=0.5\mathrm{fm}`$ .

On the lattice, there are two definitions of the “quark mass”, one through the vector Ward-Takahashi identity (VWI mass), one through the axial identity (AWI mass). The bare masses (open symbols in Fig. 1) need not agree, while after renormalization (filled symbols in Fig. 1) they would (up to $`O(a^2)`$ effects) if the renormalization factors were computed non-perturbatively. For unquenched ($`N_f=2`$) data the relevant non-perturbative factors are not yet available. For this reason, I have decided to renormalize both the CP-PACS and the UKQCD data at one-loop order (the details being given in the appendices). In the case of the CP-PACS data this basically repeats their calculation , albeit with two notable differences: First, the scale is set through the measured $`r_0`$, since this is the only possibility with the UKQCD data. Second, the “boost factor” $`u_0`$ is derived from the measured plaquette rather than the plaquette in the chiral limit, since in the UKQCD data there is no uniquely defined chirally extrapolated version. This perturbative calculation then lays the ground on which the CP-PACS and UKQCD data may be compared to each other and to the chiral prediction.

Fig. 1 displays $`(M_\pi r_0)^2`$ versus the renormalized quark masses with filled symbols, open symbols indicate the bare data to visualize the shift. In the upper part the renormalization was performed with the $`c_{\mathrm{SW}}`$ values as they were used in the simulations (suggested by a mean-field analysis  or the Alpha study ). In the lower part the renormalization factors were computed with $`c_{\mathrm{SW}}=1`$ throughout, which is a consistent choice at one-loop order. Obviously, the overall consistency of the data is much better with this latter choice, as is particularly obvious from the quadratic fits (constrained to go through zero) to the AWI data. Henceforth, we shall stick to the latter choice ($`c_{\mathrm{SW}}=1`$), but it is useful to keep in mind that the difference to the upper part is a clear indication of the size of inherent perturbative uncertainties. Finally, it is worth mentioning that all error-bars in this article represent only statistical errors.

We now turn to the physics content. To convey a feeling for the scales, I mention that at $`2mr_0=1`$ the sum of the valence quark masses is of order $`400\mathrm{MeV}`$, i.e. about four times as much as in a physical $`K`$, and the corresponding “pion” weighs about $`1\mathrm{GeV}`$. The physical kaon weighs $`496\mathrm{MeV}`$, i.e. $`(M_Kr_0)^21.58`$. If the pion would satisfy $`(M_\pi r_0)^21.58`$ this would mean that its $`u`$\- and $`d`$-quarks are about half as heavy as the $`s`$-quark in the real world. The lightest pion in the CP-PACS and UKQCD simulations have $`(M_\pi r_0)^2`$ values 1.75 and 2.20, respectively, and from this we conclude that their lightest $`u`$\- and $`d`$\- quarks (in the unitary theory) have about 55% and 70% of the mass of the physical strange quark.

With such heavy “light” quarks it is a priori not clear whether XPT is of any use to extrapolate them to the physical $`u`$\- and $`d`$-quark masses. In an attempt to shed a light on this issue, Fig. 2 shows the renormalized data (as in the bottom part of Fig. 1) together with the predictions from XPT at tree/one-loop/two-loop level (LO/NLO/NNLO). The low-energy constants are taken from phenomenology, i.e. these curves represent parameter-free predictions. At LO the chiral prediction is a straight line with the slope parameter $`B`$ taken from (2); the $`\pm 1\sigma `$ bounds are indicated by dotted lines. At NLO the prediction is curved, and the numerical values of the additional parameters are taken from (13) and (17). $`\mathrm{\Lambda }_3`$ is varied within its phenomenological $`\pm 1\sigma `$ bound (full versus dotted lines – lowering $`\mathrm{\Lambda }_3`$ to $`200\mathrm{MeV}`$ yields the upper, increasing it to $`2\mathrm{GeV}`$ yields the lower dotted line), with $`B,F`$ fixed at their central values. At NNLO the parameters $`\mathrm{\Lambda }_M,k_M`$ in (18) are determined as follows. With $`\mathrm{\Lambda }_1=0.11\pm \genfrac{}{}{0pt}{}{0.04}{0.03}\mathrm{GeV},\mathrm{\Lambda }_2=1.2\pm 0.06\mathrm{GeV}`$  and (17) the relation beneath eqn. (18) gives

$$\mathrm{\Lambda }_M=0.60\pm 0.03\mathrm{GeV}$$
(19)
where errors have been added in quadrature. This means that over the range considered the uncertainty in $`\mathrm{\Lambda }_M`$ is neglibigle compared to that coming from $`k_M`$. Phenomenological arguments indicate $`|k_M|2`$ , and the sum-rule estimates for the NNLO counterterms given in  may be converted into a more accurate estimate. For the purpose of the present article it is sufficient to use $`k_M=0\pm 2`$, and the associate curves (holding $`B,F,\mathrm{\Lambda }_3,\mathrm{\Lambda }_M`$ fixed at their central values) are included in Fig. 2.

This information is sufficient to asses the chiral convergence behaviour in $`M_\pi ^2`$ versus $`2m`$. To that aim we compare the LO/NLO/NNLO predictions for a given quark mass. Fig. 3 shows the relative shift in $`M_\pi ^2`$ when one more order is included or the new low-energy constants are varied within reasonable bounds. At first sight, the uncertainties at higher orders due to the error bars of the counterterms seem large, but one should keep in mind that the associate shifts are 100% correlated over the whole range. For instance, if the CP-PACS VWI point at $`2mr_00.2`$ sits on the $`1\sigma `$ curve (upper dotted NLO line in Fig. 2), then the one at $`2mr_00.4`$ should too, if only $`\mathrm{\Lambda }_3`$ needs to be adapted. This means that precise lattice data are ideally suited to reduce the error on $`l_3^\mathrm{r}(\mu )`$ (or $`\mathrm{\Lambda }_3`$). Fig. 3 shows that the crossover where the uncertainty due to the NLO exceeds that due to the LO contribution is at $`2mr_0=0.45/0.02/0.10`$ when $`\mathrm{\Lambda }_3`$ is taken at its central/+1$`\sigma `$/-1$`\sigma `$ value. By averaging with weights 2/1/1 one arrives at the estimate that in the case of $`M_\pi `$ the chiral expansion is sufficiently well-behaved that the NLO functional form is useful for quark masses up to

$$2mr_00.25(M_\pi r_0)^22M_\pi 560\mathrm{MeV},$$
(20)
and maybe more. This, if correct, means that current state-of-the-art simulations make contact with the regime where NLO-XPT holds, but so far a non-trivial “lever arm” which is needed to make model independent predictions in the deeply chiral regime, is likely missing.

While it is clear that future simulation data will allow to test the prediction (20) and provide, if it is correct, the lever arm needed, one might, already at this time, go ahead and try what comes out if one assumes that the estimate (20) is too pessimistic and hence uses the NLO chiral ansatz to fit the data over an extended range. This is what we shall do below, but it is clear that this attempt is rather speculative and results should be taken with care.

Since everything below is about the deviation from a linear relationship, it is useful to make the curvature optically visible. This is conveniently done by dividing out a factor $`2mr_0`$. The result is displayed in Fig. 4 where the predictions from XPT at LO/NLO/NNLO are included for completeness. In this representation, the genuine feature of the NLO curve is that it lies below the LO constant for light pions, but above if the (LO-) pion mass is larger than $`\mathrm{\Lambda }_3`$.

To fit the data, one replaces $`2mB`$ by $`M_\pi ^2`$, since in this form the data are plotted against something which is directly measured and the change is of yet higher order in the chiral counting. In other words, the statement is that one should use the representation

$$\frac{\stackrel{~}{M}_\pi ^2}{2\stackrel{~}{m}}=\stackrel{~}{B}\left\{1\frac{\stackrel{~}{M}_\pi ^2}{32\pi ^2\stackrel{~}{F}_\pi ^2}\mathrm{log}\left(\frac{\stackrel{~}{\mathrm{\Lambda }}_3^2}{\stackrel{~}{M}_\pi ^2}\right)\right\}$$
(21)
and adjust the 3 dimensionless parameters $`\stackrel{~}{B}Br_0,\stackrel{~}{F}F_\pi r_0,\stackrel{~}{\mathrm{\Lambda }}_3\mathrm{\Lambda }_3r_0`$.

Attempts to fit the data with either the VWI or the AWI definition of the quark mass to formula (21) are summarized in Tab. 1 and – for the AWI mass – illustrated in Fig. 5. Since the parameter $`\stackrel{~}{F}_\pi =F_\pi r_0=92.4\mathrm{MeV}\mathrm{\hspace{0.17em}0.5}\mathrm{fm}0.234`$ is known, it is held fixed at its physical value in (a, b) and fitted only in case (c). Fit (b) differs from (a) in that the data point with the heaviest quark mass has been omitted. From the overall spread it is clear that there are considerable systematic uncertainties and this is why I refrain from quoting statistical errors.

It is clear that with the quality of the present data one cannot claim direct evidence for chiral logs. Nonetheless, it is worth noticing that fits (a) and (b) indicate that the data are consistent with the logarithmic form (21) suggested by NLO chiral perturbation theory. Moreover, that the fit with $`F_\pi `$ held fixed at its physical value does not dramatically change if the heaviest data point is omitted (a$``$b), suggests that the estimate (20) for the permissible range at NLO level is indeed reasonable. What one can learn from Table 1 is that in QCD with $`N_f=2`$ the low-energy parameter $`\mathrm{\Lambda }_3`$ does not differ by orders of magnitude from the original estimate by Gasser and Leutwyler ; all $`\mathrm{\Lambda }_3`$ values in Tab. 1 lie between the central and the $`+1\sigma `$ level in (12). This is important, since it excludes – at least for $`N_f=2`$ – an alternative scenario of the chiral symmetry breaking in which $`B`$ would not be the adequate order parameter  (for an introduction to this topic see ).

## Discussion

The aim of the present note has been to compare the CP-PACS and the UKQCD data for $`M_\pi ^2`$ as a function of the quark mass to each other and to the prediction from Chiral Perturbation Theory (XPT). We have seen that after renormalizing with 1-loop lattice perturbation theory the two sets are reasonably consistent (maybe with the exception of the CP-PACS data with the VWI definition of the quark mass). After dividing out a factor $`2mr_0`$ to “zoom in” on the deviation from a linear behaviour, both sets show very little curvature. Postponing the issue whether or not it is permissible to compare non-continuum-extrapolated lattice data to the prediction from XPT in the continuum we found that the current size of the error-bars assures consistency with the logarithmic form (21), albeit one cannot claim evidence for chiral logs.

The situation might change as more precise data are being released. In fact, the JLQCD collaboration claims that their data with $`N_f=2`$ non-perturbatively $`O(a)`$ improved clover fermions and standard (Wilson) glue are inconsistent with NLO chiral perturbation theory, if they keep the parameter $`F_\pi `$ fixed at its physical value .

At this point it is mandatory to think about possible reasons why the data might not coincide with the chiral prediction. An incomplete list is the following one:

1. They need not, if the pions are too heavy. Most of the data have been obtained in a regime of quark masses where there is no guarantee that XPT works. Phenomenological experience tells us that this expansion works for mesons as heavy as the physical $`K`$, i.e. for $`(M_\pi r_0)^21.58`$. I have argued that this condition might be relaxed to (20). Even if this is true, at best 1 point from either set survives, and hence there is, strictly speaking, no room for trying the NLO-XPT ansatz which has 3 parameters.
2. Scaling violations, in particular implications of the broken chiral symmetry might be so severe that the chiral logs – even if they exist in the mass range considered – might be buried. To check one would have to perform a continuum extrapolation or attempt a dynamical simulation with fermions which obey the Ginsparg-Wilson relation.
3. For unknown reasons (e.g. an algorithmic flaw), the data might represent a partially quenched rather than the fully unquenched situation. Obviously, this is a rather remote possibility. The point is that the NLO prediction $`M_\pi ^2/M^2=1+\mathrm{const}M^2\mathrm{log}(M^2/\mathrm{\Lambda }_3^2)`$ has no genuine $`M^2`$ term. On the other hand, in the partially quenched case a true $`M^2`$ contribution exists ; hence even perfectly linear behaviour in a plot analogous to Fig. 4 or 5 does not necessarily imply unreasonable values of the low-energy constants.

Finally, I would like to highlight 4 points:

($`i`$) So far, we have ignored systematic uncertainties. Going back to Fig. 1 one sees that they are far from being negligible. What we see exemplifies the standard wisdom that perturbative renormalization factors which turn out to deviate from 1 by, say, 10% call for either repeating the exercise with non-perturbatively determined $`Z`$-factors or at much smaller lattice spacing. In the present context this means that the figures in Tab. 1 should not be taken as lattice determinations of the low-energy constants $`F,B,\mathrm{\Lambda }_3`$. Nonetheless, even this preliminary attempt shows that $`\mathrm{\Lambda }_3`$ cannot differ by orders of magnitude from the original estimate by Gasser and Leutwyler . Moreover, as the data get more precise – and under the proviso that consistency is maintained – lattice determinations of the scales $`\mathrm{\Lambda }_3`$ and $`\mathrm{\Lambda }_4`$ in conjunction with chiral perturbation theory will yield precise predictions of intrinsically Minkowskian quantities, e.g. the $`I=0,2`$ pi-pi scattering lengths (see  for the connection).

$`(ii)`$ For a thorough comparison with chiral perturbation theory, one needs to get control over the lattice artefacts. A safe way to do this is to first extrapolate all data to the continuum. Then one can determine the permissible mass range and extrapolate in a second step to the chiral limit. If it turns out that this order cannot be sustained in practice, the chiral framework may be extended to account for the main discretization effects. For the case of unimproved fermions this has been done , and it is clear that this approach could be generalized to improved actions. The only disadvantage is the added number of counterterms that need to be fixed from the data.

($`iii`$) In the present note the scale is set through $`\widehat{r}_0^1`$ throughout. In the CP-PACS studies, $`\widehat{r}_0`$ has been found to strongly depend on $`\widehat{m}^{\mathrm{VWI}}`$ or $`\widehat{M}_\pi ^2`$ . If $`r_0`$ itself depends on the sea-quark mass, normalizing all masses through the measured $`\widehat{r}_0`$ may not be the best way to compare to XPT; maybe one should use $`\widehat{r}_{0,0}`$ (the chirally extrapolated version) instead.

($`iv`$) Strictly speaking, the “phenomenological” value of $`l_3^\mathrm{r}`$ (or $`\mathrm{\Lambda }_3`$) that we used in comparing the lattice data to the prediction from XPT is not quite adequate, since the phenomenological determination is with $`m_\mathrm{s}`$ fixed at its physical value, while the lattice studies are $`N_f=2`$ simulations, in other words here $`m_\mathrm{s}`$ is sent to infinity. The connection to the $`N_f=3`$ GL coefficients $`L_i^\mathrm{r}(\mu )`$ is given by

$$l_3^\mathrm{r}=8\left(2L_6^\mathrm{r}L_4^\mathrm{r}\right)+4\left(2L_8^\mathrm{r}L_5^\mathrm{r}\right)\frac{1}{576\pi ^2}\left(\mathrm{log}\left(\frac{\underset{m_{\mathrm{u},\mathrm{d}}0}{lim}M_\eta ^2}{\mu ^2}\right)+1\right)$$
(22)
\[the limit on the r.h.s. refers to a situation with $`m_\mathrm{s}`$ held fixed at its physical value\], and it is tempting to use it as the starting point of a little gedanken experiment: Assume the $`s`$-quark mass would be such that even when it is doubled the chiral expansion would not break down. Eqn. (22) tells us that under these circumstances enhancing $`m_\mathrm{s}`$ by a factor $`2`$ (and hence roughly doubling $`M_\eta ^2`$, too) would lower $`l_3^\mathrm{r}`$ by $`0.0001`$. Since (12) translates into $`l_3^\mathrm{r}(\mu M_\rho )=0.0008\pm 0.0038`$ such a shift would be negligible compared to the error-bar, and it seems therefore reasonable to assume that the effect of an added third flavour in the simulations, fixed at the physical $`s`$-quark mass (which would fully justify the comparison with $`N_f=2`$ XPT), would be small compared to the inherent theoretical uncertainties.

Formula (22) is interesting in yet another respect, since it tells us that lattice studies which pin down $`M_\pi ^2`$ versus $`2m`$ (and hence $`\mathrm{\Lambda }_3`$ or $`l_3^\mathrm{r}`$) have a say in another issue. In the 3-flavour theory there is the famous “Kaplan Manohar ambiguity”, i.e. the chiral Lagrangian stays invariant under a simultaneous transformation of the quark masses

$$m_\mathrm{u}m_\mathrm{u}+\lambda m_\mathrm{d}m_\mathrm{s},m_\mathrm{d}m_\mathrm{d}+\lambda m_\mathrm{s}m_\mathrm{u},m_\mathrm{s}m_\mathrm{s}+\lambda m_\mathrm{u}m_\mathrm{d},$$
(23)
and an appropriate modification of $`L_6^\mathrm{r},L_7^\mathrm{r},L_8^\mathrm{r}`$. This would, in principle, allow to tune $`m_\mathrm{u}=0`$ and hence provide a simple and elegant solution to the strong CP problem . In XPT terminology it is the low-energy combination $`2L_8^\mathrm{r}(\mu )L_5^\mathrm{r}(\mu )`$ that decides whether this is a viable option . In the past, XPT has been augmented by theoretical assumptions or model calculations to exclude $`m_\mathrm{u}=0`$. More recently it has been proposed to determine the $`L_i^\mathrm{r}`$ from lattice simulations  and important steps in this program have been taken . Formula (22) tells us that a lattice determination of $`l_3^\mathrm{r}`$, augmented by knowledge about $`L_4^\mathrm{r}`$ and $`L_6^\mathrm{r}`$, helps to constrain $`m_\mathrm{u}=0`$.

The present analysis certainly emphasizes the need to compute renormalization factors non-perturbatively and to perform a continuum extrapolation with dynamical data.

### Acknowledgements

This work has taken its origin in Seattle, during the INT program “Lattice QCD and Hadron Phenomenology”. Discussions with Maarten Golterman and Rainer Sommer are gratefully acknowledged, as well as useful correspondence with Derek Hepburn and Akira Ukawa. I am indebted to Heiri Leutwyler for providing me with phenomenological estimates of the NNLO constant $`k_M`$ and to Stefano Capitani for a check of the generic formulas in appendices A,B.

## Appendix A: Renormalization of the VWI quark mass

Setup with $`m^{\mathrm{VWI}}a=\mathrm{log}(1+\frac{1}{2}(\frac{1}{\kappa }\frac{1}{\kappa _\mathrm{c}}))\frac{1}{2}(\frac{1}{\kappa }\frac{1}{\kappa _\mathrm{c}})`$ and $`u_0=P^{1/4}=(\frac{1}{3}\mathrm{Tr}U_{\mathrm{}})^{1/4}`$ :

$`m_{\overline{\mathrm{MS}}}^{\mathrm{VWI}}(\mu )`$ $`=`$ $`Z_m(\mu a)\left(1+b_mam\right)m`$ (24)
$`m_{\overline{\mathrm{MS}}}^{\mathrm{VWI}}(\mu )`$ $`=`$ $`\stackrel{~}{Z}_m(\mu a)\left(1+\stackrel{~}{b}_m{\displaystyle \frac{am}{u_0}}\right){\displaystyle \frac{m}{u_0}}`$ (25)
From , naive versus tadpole-improved, with $`g^2=6/\beta `$ and $`\stackrel{~}{g}^2=\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2(2\mathrm{GeV})`$:

$`Z_m(\mu a)`$ $`=`$ $`1+{\displaystyle \frac{g^2}{4\pi }}\left({\displaystyle \frac{z_m}{3\pi }}{\displaystyle \frac{1}{\pi }}\mathrm{log}(a^2\mu ^2)\right)`$ (26)
$`\stackrel{~}{Z}_m(\mu a)`$ $`=`$ $`1+{\displaystyle \frac{\stackrel{~}{g}^2}{4\pi }}\left({\displaystyle \frac{\stackrel{~}{z}_m}{3\pi }}{\displaystyle \frac{1}{\pi }}\mathrm{log}(a^2\mu ^2)\right)`$ (27)
From  (key to orig. literature: their \[21-28\]), for generic actions:

$$\frac{1}{g_{\overline{\mathrm{MS}}}^2(\mu )}=\frac{1}{g^2}+d_g+d_fN_f+\frac{112N_f/3}{16\pi ^2}\mathrm{log}(a^2\mu ^2)$$
(28)
Wilson/Clover: $`d_g=0.4682,d_f=0.0314917`$ (for $`c_{\mathrm{SW}}=1`$) and $`P=11/3g^2`$

$`{\displaystyle \frac{1}{g_{\overline{\mathrm{MS}}}^2(\mu )}}`$ $`=`$ $`{\displaystyle \frac{1}{g^2}}0.4682+0.0314917N_f+{\displaystyle \frac{112N_f/3}{16\pi ^2}}\mathrm{log}(a^2\mu ^2)`$ (29)
$`{\displaystyle \frac{1}{\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2(\mu )}}`$ $`=`$ $`{\displaystyle \frac{P}{g^2}}0.1349+0.0314917N_f+{\displaystyle \frac{112N_f/3}{16\pi ^2}}\mathrm{log}(a^2\mu ^2)`$ (30)
$`z_m`$ $`=`$ $`12.953+7.738c_{\mathrm{SW}}1.380c_{\mathrm{SW}}^2`$ (31)
$`\stackrel{~}{z}_m`$ $`=`$ $`z_m\pi ^2=\{\begin{array}{cc}\hfill \text{13.0}& \text{(}c_{\mathrm{SW}}2\text{)}\\ \hfill \text{9.44}& \text{(}c_{\mathrm{SW}}=1\text{)}\end{array}`$ (34)
Iwasaki/Clover: $`d_g=0.1000,d_f=0.0314917`$ (for $`c_{\mathrm{SW}}=1`$) and $`P=10.1402g^2`$,

$`R=10.2689g^2`$, thus $`1=P+0.1402g^2=R+0.2689g^2=3.648P80.331R0.2006g^2`$

$`{\displaystyle \frac{1}{g_{\overline{\mathrm{MS}}}^2(\mu )}}`$ $`=`$ $`{\displaystyle \frac{1}{g^2}}+0.1000+0.0314917N_f+{\displaystyle \frac{112N_f/3}{16\pi ^2}}\mathrm{log}(a^2\mu ^2)`$ (35)
$`{\displaystyle \frac{1}{\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2(\mu )}}`$ $`=`$ $`{\displaystyle \frac{P}{g^2}}+0.2402+0.0314917N_f+{\displaystyle \frac{112N_f/3}{16\pi ^2}}\mathrm{log}(a^2\mu ^2)`$ (36)
$`{\displaystyle \frac{1}{\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2(\mu )}}`$ $`=`$ $`{\displaystyle \frac{R}{g^2}}+0.3689+0.0314917N_f+{\displaystyle \frac{112N_f/3}{16\pi ^2}}\mathrm{log}(a^2\mu ^2)`$ (37)
$`{\displaystyle \frac{1}{\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2(\mu )}}`$ $`=`$ $`{\displaystyle \frac{3.648P2.648R}{g^2}}0.1006+0.0314917N_f+{\displaystyle \frac{112N_f/3}{16\pi ^2}}\mathrm{log}(a^2\mu ^2)`$ (38)
$`z_m`$ $`=`$ $`4.858+5.301c_{\mathrm{SW}}1.267c_{\mathrm{SW}}^2`$ (39)
$`\stackrel{~}{z}_m`$ $`=`$ $`z_m0.4206\pi ^2=\{\begin{array}{cc}\hfill \text{5.76}& \text{(}c_{\mathrm{SW}}=1.47\text{)}\\ \hfill \text{4.74}& \text{(}c_{\mathrm{SW}}=1\text{)}\end{array}`$ (42)
In  one finds:

$`\mathrm{Wilson}/\mathrm{Clover}:b_m`$ $`=`$ $`1/20.09623g^2+O(g^4)`$ (43)
$`\stackrel{~}{b}_m`$ $`=`$ $`b_m(11/12\stackrel{~}{g}^2)1/20.05456\stackrel{~}{g}^2`$ (44)
$`\mathrm{Iwasaki}/\mathrm{Clover}:b_m`$ $`=`$ $`1/20.0509g^2+O(g^4)`$ (45)
$`\stackrel{~}{b}_m`$ $`=`$ $`b_m(10.03505\stackrel{~}{g}^2)1/20.0334\stackrel{~}{g}^2`$ (46)
## Appendix B: Renormalization of the AWI quark mass

Setup with $`m^{\mathrm{AWI},\mathrm{imp}}a=\frac{1}{2}_\mu A_\mu ^{\mathrm{a},\mathrm{imp}}(x)O^\mathrm{a}(0)/P^\mathrm{a}(x)O(0)`$ and $`u_0=P^{1/4}=(\frac{1}{3}\mathrm{Tr}U_{\mathrm{}})^{1/4}`$ :

$`m_{\overline{\mathrm{MS}}}^{\mathrm{AWI}}(\mu )`$ $`=`$ $`{\displaystyle \frac{Z_A}{Z_P(\mu a)}}{\displaystyle \frac{1+b_Aam}{1+b_Pam}}m^{\mathrm{AWI},\mathrm{imp}}`$ (47)
$`m_{\overline{\mathrm{MS}}}^{\mathrm{AWI}}(\mu )`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{Z}_A}{\stackrel{~}{Z}_P(\mu a)}}{\displaystyle \frac{1+\stackrel{~}{b}_Aam/u_0}{1+\stackrel{~}{b}_Pam/u_0}}m^{\mathrm{AWI},\mathrm{imp}}`$ (48)
From , naive versus tadpole-improved, with $`g^2=6/\beta `$ and $`\stackrel{~}{g}^2=\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2(2\mathrm{GeV})`$:

$`Z_A`$ $`=`$ $`1+{\displaystyle \frac{g^2}{4\pi }}{\displaystyle \frac{z_A}{3\pi }},Z_P(\mu a)=\mathrm{\hspace{0.33em}1}+{\displaystyle \frac{g^2}{4\pi }}\left({\displaystyle \frac{z_P}{3\pi }}+{\displaystyle \frac{1}{\pi }}\mathrm{log}(a^2\mu ^2)\right)`$ (49)
$`\stackrel{~}{Z}_A`$ $`=`$ $`1+{\displaystyle \frac{\stackrel{~}{g}^2}{4\pi }}{\displaystyle \frac{\stackrel{~}{z}_A}{3\pi }},\stackrel{~}{Z}_P(\mu a)=\mathrm{\hspace{0.33em}1}+{\displaystyle \frac{\stackrel{~}{g}^2}{4\pi }}\left({\displaystyle \frac{\stackrel{~}{z}_P}{3\pi }}+{\displaystyle \frac{1}{\pi }}\mathrm{log}(a^2\mu ^2)\right)`$ (50)
From , for generic actions:

Wilson/Clover:

$`z_A`$ $`=`$ $`15.7970.248c_{\mathrm{SW}}+2.251c_{\mathrm{SW}}^2`$ (51)
$`\stackrel{~}{z}_A`$ $`=`$ $`z_A+\pi ^2=\{\begin{array}{cc}\hfill \text{2.581}& \text{(}c_{\mathrm{SW}}2\text{)}\\ \hfill \text{-3.924}& \text{(}c_{\mathrm{SW}}=1\text{)}\end{array}`$ (54)
$`z_P`$ $`=`$ $`22.596+2.249c_{\mathrm{SW}}2.036c_{\mathrm{SW}}^2`$ (55)
$`\stackrel{~}{z}_P`$ $`=`$ $`z_P+\pi ^2=\{\begin{array}{cc}\hfill \text{-16.372}& \text{(}c_{\mathrm{SW}}2\text{)}\\ \hfill \text{-12.513}& \text{(}c_{\mathrm{SW}}=1\text{)}\end{array}`$ (58)
Iwasaki/Clover:

$`z_A`$ $`=`$ $`8.1920.125c_{\mathrm{SW}}+1.610c_{\mathrm{SW}}^2`$ (59)
$`\stackrel{~}{z}_A`$ $`=`$ $`z_A+0.4206\pi ^2=\{\begin{array}{cc}\hfill \text{-0.746}& \text{(}c_{\mathrm{SW}}=1.47\text{)}\\ \hfill \text{-2.556}& \text{(}c_{\mathrm{SW}}=1\text{)}\end{array}`$ (62)
$`z_P`$ $`=`$ $`10.673+1.601c_{\mathrm{SW}}1.281c_{\mathrm{SW}}^2`$ (63)
$`\stackrel{~}{z}_P`$ $`=`$ $`z_P+0.4206\pi ^2=\{\begin{array}{cc}\hfill \text{-6.936}& \text{(}c_{\mathrm{SW}}=1.47\text{)}\\ \hfill \text{-6.202}& \text{(}c_{\mathrm{SW}}=1\text{)}\end{array}`$ (66)
In  one finds:

$`\mathrm{Wilson}/\mathrm{Clover}:b_A`$ $`=`$ $`1+0.15219(5)g^2+O(g^4)`$ (67)
$`\stackrel{~}{b}_A`$ $`=`$ $`b_A(11/12\stackrel{~}{g}^2)\mathrm{\hspace{0.33em}1}+0.06886(5)\stackrel{~}{g}^2`$ (68)
$`b_P`$ $`=`$ $`1+0.15312(3)g^2+O(g^4)`$ (69)
$`\stackrel{~}{b}_P`$ $`=`$ $`b_P(11/12\stackrel{~}{g}^2)\mathrm{\hspace{0.33em}1}+0.06979(3)\stackrel{~}{g}^2`$ (70)
$`\mathrm{Iwasaki}/\mathrm{Clover}:b_A`$ $`=`$ $`1+0.0733(5)g^2+O(g^4)`$ (71)
$`\stackrel{~}{b}_A`$ $`=`$ $`b_A(10.03505\stackrel{~}{g}^2)\mathrm{\hspace{0.33em}1}+0.0383(5)\stackrel{~}{g}^2`$ (72)
$`b_P`$ $`=`$ $`1+0.0744(12)g^2+O(g^4)`$ (73)
$`\stackrel{~}{b}_P`$ $`=`$ $`b_P(10.03505\stackrel{~}{g}^2)\mathrm{\hspace{0.33em}1}+0.0394(12)\stackrel{~}{g}^2`$ (74)
In all cases the relationship between $`z_X`$ and $`\stackrel{~}{z}_X`$ as well as $`b_X`$ and $`\stackrel{~}{b}_X`$ ($`X\{m,A,P\}`$) is mine.

## Appendix C: Implementation with tadpole resummation

The renormalization with $`u_0`$ defined via the plaquette is traced in Tab. 2 for the UKQCD data. Alternatively, one could define $`u_0`$ as $`1/(8\kappa _\mathrm{c})`$, but this would bring different perturbative coefficients than those listed in App. A/B. Notice the effect of the tadpole resummation, i.e. the difference of the coupling $`g_{\overline{\mathrm{MS}}}^2`$ computed via (29) and $`\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2`$ via (30). The statistical uncertainty for the VWI quark mass is larger than that for the AWI definition due to a limited accuracy of $`\kappa _\mathrm{c}`$, which is defined in a partially quenched sense .

For the CP-PACS data, both the plaquette $`P`$ and the rectangle $`R`$ are published , and this gives, in principle, the option to define the tadpole resummation (besides the usual option $`1/(8\kappa _\mathrm{c})`$) via (37) or (38), where the latter choice reflects the specific combination used in the action. All these options would, however, imply different perturbative coefficients than those listed in App. A/B, i.e. we restrict ourselves to the choice (36) along with $`u_0P^{1/4}`$. It is interesting to see that the CP-PACS RG improved action achieves agreement of $`g_{\overline{\mathrm{MS}}}^2`$ with the “standard” $`\stackrel{~}{g}_{\overline{\mathrm{MS}}}^2`$ (via (36)).

In the argument of the logarithm that converts to the $`\overline{\mathrm{MS}}`$ scheme, the lattice spacing is multiplied with a physical scale, here $`2\mathrm{GeV}`$. This means that $`a`$ must be assigned a physical value, as well. To compare like with like this is done via the measured $`r_0`$ in both sets, assuming $`r_0=0.5\mathrm{fm}`$ in physical units. As a consequence, our $`Z`$ factors for the CP-PACS data depend slightly on the quark mass, even though we work in a mass-independent scheme.


text:
# AUTOMATED STELLAR SPECTRAL CLASSIFICATION AND PARAMETERIZATION FOR THE MASSES

## 1 A Brief History of Automated Classification

Current or planned large-scale surveys, such as the Sloan Digital Sky Survey (York et al. 2000) or the GAIA mission (scheduled for launch around 2011), have led to increased interest in automated spectral classifiers (e.g., Bailer-Jones 2000). There are many other reasons to develop automated classifiers, not the least of which are the homogeneity of the results and the repeatability of the process. Automated spectral classification of stars goes back decades. For instance, in an early attempt Jones (1966) fit a few major stellar lines and correlated these indexes with MK type (Morgan, Keenan, & Kellman 1943). Malyuto & Shvelidze (1994) later developed this technique further. Unfortunately, the line fitting technique suffers from the disadvantage that one has first to know the approximate stellar type before determining which lines to fit, otherwise very different features will be found at the same wavelengths. Kurtz (1983; see also LaSala 1994) developed a minimum vector distance technique that matched spectra to a library of standards weighting the comparison to different spectral regions for different types of stars. The minimum vector distance technique has had some success – classifying stellar spectral types to within $`\sigma `$=2.2 spectral subtypes – but refining this technique is cumbersome since the weighting vectors need to be carefully established, yet they vary as a function of spectral type and luminosity class.

In the middle of the last decade four independent groups (Gulati et al. 1994; von Hippel et al. 1994; Vieira & Ponz 1995; Weaver & Torres-Dodgen 1995) began to successfully apply artificial neural networks (ANNs) to the spectral classification problem. Neural networks are trained to yield classifications that are identical to those previously assigned, and thus have the advantage that the ANN builder need not become a classification expert, but rather can rely on the true experts in the field, via their many previous classifications. For example, von Hippel et al. built an ANN using the objective prism spectra and classifications of Houk (1982, and references therein), who had at that time already classified more than 105 stars. This first generation of neural networks were applied to low resolution ultraviolet (Vieira & Ponz 1995) or optical (the other three studies) spectra and achieved $`\sigma `$=0.6 spectral subtypes and $`\sigma `$=0.35 luminosity classes for A stars and returned E($`BV`$) (Weaver & Torres-Dodgen 1995), or $`\sigma <2`$ subtypes over the broad range of O3 to M4 stars. These studies validated the ANN approach and indicated its tremendous potential.

The second generation of ANN spectral classification studies (Weaver & Torres-Dodgen 1997; Bailer-Jones, Irwin, & von Hippel 1998; Singh, Gulati, & Gupta 1998; Weaver 2000) focused on increasing sample size and moving on to two-dimensional classification. At this point spectral type classification became very good, with $`\sigma `$=0.5–0.7 subtypes. Luminosity classification quality varied from $`\sigma 0.3`$ luminosity classes (Weaver & Torres-Dodgen 1997) to a statistically reliable luminosity classification for dwarfs and giants, though not sub-giants (Bailer-Jones et al. 1998). Interestingly, Weaver (2000) also showed that he could provide two dimensional classification for both components of artificial binaries!

A few ANN studies (Bailer-Jones et al. 1997; Snider et al. 2001) moved from MK stellar classification to parameterization of astrophysical parameters ($`T_{\mathrm{eff}}`$, log($`g`$), \[Fe/H\]). Although the philosophy of classification and parameterization are different, they strive to serve the same community. Classification seeks to place any program star within the framework defined by a series of standards. Since the MK classification system is so widely used and the connection to stellar parameters such as absolute magnitude, surface temperature, and mass, are in general well known, this has been a productive route for studying individual stars, local stellar populations, and Galactic structure. Stellar parameterization seeks to skip the initial step of MK classification and directly determine atmospheric parameters. The cost is that the results are model dependent, but in many parts of the HR diagram model spectra look very much, though not exactly, like real stellar spectra. In addition, model atmospheres can be easily constructed for subsolar metallicity, and therefore these can be applied to spectra which would not be possible to classify on the MK system. Bailer-Jones et al. (1997) passed their objective prism spectra through ANNs trained on stellar atmosphere models and derived a detailed mapping between MK classifications and the Kurucz (1979, 1992) model atmosphere set they used (as implemented with the program SPECTRUM by Gray & Corbally 1994). They also reported that the mean metallicity in the solar neighborhood, as represented by Houk’s (1982) objective prism spectra, is slightly subsolar, at \[Fe/H\] = $`0.2`$. Snider et al. (2001) turned the problem around and trained ANNs on real stellar spectra using atmospheric parameters previously derived from fine abundance analysis work in the literature, now in three-dimensional parameter space. They achieved $`\sigma `$($`T_{\mathrm{eff}})150`$ K, $`\sigma `$(log($`g`$)) $``$ 0.33 dex, and $`\sigma `$(Fe/H) $``$ 0.2 dex.

## 2 A Few Lessons Learned

Here we offer a few comments on lessons we have learned in applying ANNs to the problems of spectral classification and parameterization:

* Classification in spectral type and parameterization in $`T_{\mathrm{eff}}`$ are easy, the ANNs have little trouble finding a good global solution, and the results are generally precise to less than a spectral subtype or 200 K.
* Luminosity classification and log($`g`$) parameterization are possible, but more difficult. Both require spectra with an adequate combination of resolution, wavelength coverage, and S/N. This spectral quality is just achieved with classical MK objective prism resolution and wavelength coverage.
* ANNs are best suited to stellar classification/parameterization when a single wavelength range and resolution are used for a particular ANN.
* If low S/N spectra are used, besides the higher random errors, ANNs may make systematic errors unless they have been trained on low S/N spectra (Snider et al. 2001).
* Supervised ANNs treat spectra as patterns with a known correlation between those patterns and answers, and attempt to learn that relationship. The better one can homogenize the training data so that the ANN does not find spurious correlations between the input catalog and the answers the better results one achieves. Larger catalogs always help in this regard, as do uniform data sets. Spurious correlations caused by real astrophysics can also be a problem. As an example of this phenomenon, the Houk catalog is magnitude-limited with $`V_{\mathrm{limit}}`$ = 10–11. This magnitude limit creates a correlation between spectral type and luminosity class, i.e., the catalog contains mostly early type dwarfs and late type giants. Without proper scrutiny one might believe a trained ANN has achieved true luminosity classification using such a catalog when the ANN could have learned to classify luminosity statistically based on spectral type.
* Principle Component Analysis (PCA) can be used to compress stellar spectra or to remove spurious signals (e.g., Storrie-Lombardi et al. 1995). PCA also bears a strong resemblance to ANNs (Lahav 1995), and is a good pedagogical tool to gaining a heuristic understanding of ANNs. By creating a series of vectors, a linear combination of which will recreate any star in the library, PCA recasts a stellar spectral library in much the same way as the hidden nodes in a single hidden layer ANN.
* Spectral classification errors are a function of spectral type, based largely on the number of examples and variance within a given type or range of types. For example, Bailer-Jones et al. (1998) found $`\sigma `$(SpT)=0.5 for B3 to A0 stars and $`\sigma `$(SpT)=0.8 for F3 to G1 stars. The lower errors for the B and A stars are probably the result of reduced sensitivity of their spectra to abundance differences.
* It is easiest to build multiple ANNs, each specializing in a specific dimension, when solving multi-dimensional spectral classification or parameterization problems. Not only are the ANNs more likely to converge on a good global solution, but less data are required for this approach. For example, Snider et al. (2001) built ANNs specializing in each of $`T_{\mathrm{eff}}`$, log($`g`$), and \[Fe/H\]. Weaver (2000) has also found it helpful to use one ANN for initial rough classifications, followed by specialist ANNs, trained on a limited spectral type ranges, to refine the classifications.

## 3 Proposed Classification and Parameterization for the Masses

Can we build automated spectral classifiers for general use? For some time Bob Garrison has pointed out that we could build stellar classifiers for particular spectrographs. Users of such spectrographs might have a near real-time reduction pipeline, immediately following which they would receive a stellar classification. This is certainly possible. In fact, if the data were obtained in a standard manner, the only reduction steps required prior to ANN classification would be spectral extraction and wavelength calibration.

In practice, neither we nor, to the best of our knowledge, any other group has taken this approach. The difficulties are not technical, but rather the time-consuming nature of building multiple such classification machines for the many possible spectrographs in use by stellar spectroscopists. Certainly, if a stellar spectroscopic survey of sufficient size were to be undertaken which would create a uniform data set of sufficient quality, we and others would be motivated to build a tailor-made stellar classification or parameterization machine for that instrument/survey.

We propose instead a thematically related approach. Instead of building ANN classification/parameterization machines for particular spectrographs, we propose to build such machines for particular combinations of wavelength coverage, resolution, and S/N, and make these available for use via the web. It would be up to the user to process their data onto a linear flux and wavelength scale at one of the resolutions and wavelength ranges supported by our web site algorithms. Users would upload their spectra, run the classification or parameterization ANNs, and receive a spectral type and luminosity class and/or the stellar astrophysical parameters, along with the associated uncertainties.

We hope to develop such a tool first by beginning with stellar parameterization based on model atmospheres. Our entire approach would be modular and would initially support a single resolution and wavelength range, while covering the parameter range 4500 $`T_{\mathrm{eff}}8000`$, 2 $``$ log($`g`$) $``$ 5, $`4.5`$ \[Fe/H\] $`+0.5`$, and 0 $``$ microturbulence $``$ 1/2 the stellar rotation break-up velocity. Our initial resolution and wavelength range have not been finalized, but would probably be $`R`$ 2000 and 150 Å around H$`\beta `$, respectively. This would allow anyone with higher resolution spectra or spectra with a broader range of wavelengths to take advantage of our first automated parameterization algorithms. As we move to support a wider range of spectral resolutions and wavelength coverages we also intend to add modules to increase our effective temperature range to $`T_{\mathrm{eff}}`$ = 50,000, decrease our surface gravity range to log($`g`$) = 1, include different relative O, Mg, and Ca abundances, and begin spectral classification on the MK system for stars of near solar metallicity. Eventually we hope to push the stellar parameters into the M, L, and T dwarf regimes. We recognize that model atmospheres are never perfectly accurate and that they improve with time. From time to time, where meaningful advances have been made for a particular range of stellar atmospheres, we will upgrade our parameterization modules.

Our anticipated users are spectroscopists doing fine-abundance analysis who want a starting point or a sanity check on their result, those conducting surveys who need classifications or parameterizations for statistical or pre-selection purposes, and those wanting independent determinations of the classifications/parameters for their program stars for a wide variety of studies.

ACKNOWLEDGMENTS

We wishes to thank Bob Garrison, Richard Gray, and Coryn Bailer-Jones for thought-provoking discussions. We gratefully acknowledge major financial support for this work from the State of Texas Advanced Research Program.

REFERENCES

* Bailer-Jones, C. A. L. 2000 A&A 357, 197

Bailer-Jones, C. A. L., Irwin, M., Gilmore, G. & von Hippel, T. 1997 MNRAS 292, 157

Bailer-Jones, C. A. L., Irwin, M. & von Hippel, T. 1998 MNRAS 298, 361

Gray, R. O. & Corbally, C. J. 1994 AJ, 107, 742

Gulati, R. K., Gupta, R., Gothoskar, P. & Khobragade, S. 1994 ApJ, 426, 340

Houk, N. 1982 University of Michigan Catalogue of Two-Dimensional Spectral Types for the HD Stars, Vol. 3

Jones, D. H. P. 1966 R. Obs. Bull. 126, 219

Kibblewhite, E. J., Bridgeland, M. T., Bunclark, P. S. & Irwin, M. J., 1984 (ed.) Klinglesmith D.A. NASA-2317, Astronomical Microdensitometry Conference. NASA, Washington D.C., p. 277

Kurtz, M. J. 1983 (ed.) Garrison, B. F., The MK Process and Stellar Classification, David Dunlop Observatory, Toronto, p. 136

Kurucz, R. L. 1979 ApJS, 40, 1

Kurucz, R. L. 1992 (ed.) B. Barbuy & A. Renzini, Stellar Populations of Galaxies, IAU Symp. 149, Dordrecht, Kluwer, 225

Lahav, O. 1995 Vistas in Astronomy 38, 251

LaSala, J. 1994 (ed.) C. J. Corbally, R. O. Gray, & R. F. Garrison, The MK Process at 50 Years: A Powerful Tool for Astrophysical Insight, ASP Conf. Ser. 60, San Francisco, ASP, 312

Malyuto, V. & Shvelidze, T. 1994 (ed.) C. J. Corbally, R. O. Gray, & R. F. Garrison, The MK Process at 50 Years: A Powerful Tool for Astrophysical Insight, ASP Conf. Ser. 60, San Francisco, ASP, 344

Morgan, W. W., Keenan, P. C. & Kellman, E., 1943, An Atlas of Stellar Spectra, University of Chicago Press, Chicago

Singh, H. P., Gulati, R. K. & Gupta, R. 1998 MNRAS 295, 312

Snider, S., Allende Prieto, C., von Hippel, T., Beers, T. C., Sneden, C., Qu, Y. & Rossi, S. 2001 ApJ 562, 528

Storrie-Lombardi, M. C., Irwin, M. J., von Hippel, T. & Storrie-Lombardi, L. J. 1995 Vistas in Astronomy 38, 331

Vieira, E. F. & Ponz, J. D. 1995 A&AS 111, 393

von Hippel, T., Storrie-Lombardi, L. J., Storrie-Lombardi, M. & Irwin, M. J. 1994 MNRAS 269, 97

Weaver, W. B. 2000 ApJ 541, 298

Weaver, W. B. & Torres-Dodgen, A. V. 1995 ApJ 446, 300

Weaver, W. B. & Torres-Dodgen, A. V. 1997 ApJ 487, 847

York, D. G. et al. 2000 AJ 120, 1579