text:
# Spontaneous thermal expansion of nematic elastomers

## 1 Introduction

Liquid crystalline elastomers (LCE) have attracted a significant experimental and theoretical interest. Recent review articles summarise much of the recent-years research . The novel behavior of these materials arises from a coupling between the liquid crystalline ordering of mesogenic moieties and the elastic properties of the polymer network. Several unusual physical effects have been discovered in LCE: (a) Spontaneous, reversible shape changes on heating or illumination. (b) “Soft elasticity” – mechanical deformation without (or with very low) stress. (c) Mechanical instabilities and discontinuous stress-strain strain relations on switching of nematic director by mechanical fields . (d) Electrical switching of optical properties with accompanying mechanical strains in . (e) Solid phase nematohydrodynamics and unusual rheology and dynamics, leading to anomalous dissipation, mechanical damping and dynamic softness .

Theoretical work has been able to develop simple molecular models of ideal nematic network and qualitatively describe or predict nearly all observed phenomena with a description having no free parameters, based only on independently measurable quantities . This theory of an ideal nematic rubber is encompassed in the expression for the free energy density of a nematic rubber in response to a arbitrary general deformation expressed by the strain tensor $`\underset{¯}{\underset{¯}{𝝀}}`$:

$$F_{\mathrm{elast}}=\frac{1}{2}\mu \mathrm{Tr}\left[\underset{¯}{\underset{¯}{\mathbf{}}}_0\underset{¯}{\underset{¯}{𝝀}}^\mathrm{T}\underset{¯}{\underset{¯}{\mathbf{}}}^1\underset{¯}{\underset{¯}{𝝀}}\right],$$
(1)
where $`\mu n_\mathrm{x}k_BT`$ is the rubber modulus proportional to the crosslinking density of the network. Within the phantom Gaussian chain approximation, the matrix $`\underset{¯}{\underset{¯}{\mathbf{}}}`$ (as well as $`\underset{¯}{\underset{¯}{\mathbf{}}}_0`$) represents the anisotropic distribution of nematic chain segments after the affine deformation $`\underset{¯}{\underset{¯}{𝝀}}`$ is applied (or the initial anisotropy before deformation, in the case of $`\underset{¯}{\underset{¯}{\mathbf{}}}_0`$). The effective chain anisotropy is the crucial factor in the description. The matrix $`\underset{¯}{\underset{¯}{\mathbf{}}}`$ characterises the local nematic director $`𝒏`$ as its principal eigenvector, as well as the degree of local order expressed by the difference between its principal eigenvalues: in the component form $`\mathrm{}_{ij}=\mathrm{}_{}\delta _{ij}+[\mathrm{}_{}\mathrm{}_{}]n_in_j`$. Accordingly, the nematic polymer chain has the anisotropic radius of gyration, $`\overline{R}_{}=(\frac{1}{3}\mathrm{}_{}L)^{1/2}`$ and $`\overline{R}_{}=(\frac{1}{3}\mathrm{}_{}L)^{1/2}`$ for a Gaussian chain of contour length $`L`$. In the isotropic phase both principal values become equal to the single parameter, $`b`$, of chain persistence length: $`\mathrm{}_{ij}^{\mathrm{iso}}=b\delta _{ij}`$ and $`\overline{R}=(\frac{1}{3}bL)^{1/2}`$ along each of the three coordinate axes.

A more recent development allowed to combine the concepts of reptation theory of entangled networks with the anisotropic nature of nematic polymer strands , obtaining the “tube model” corrections to the ideal nematic rubber elastic free energy

$`F_{\mathrm{elast}}`$ $`=`$ $`{\displaystyle \frac{2}{3}}\mu {\displaystyle \frac{2M+1}{3M+1}}\mathrm{Tr}(\underset{¯}{\underset{¯}{\mathbf{}}}_0\underset{¯}{\underset{¯}{𝝀}}^\mathrm{T}\underset{¯}{\underset{¯}{\mathbf{}}}^1\underset{¯}{\underset{¯}{𝝀}})`$ (2)
$`+{\displaystyle \frac{3}{2}}\mu (M1){\displaystyle \frac{2M+1}{3M+1}}(|\underset{¯}{\underset{¯}{\mathbf{}}}^{1/2}\underset{¯}{\underset{¯}{𝝀}}\underset{¯}{\underset{¯}{\mathbf{}}}_0^{1/2}|)^2`$
$`+\mu (M1)\mathrm{ln}|\underset{¯}{\underset{¯}{\mathbf{}}}^{1/2}\underset{¯}{\underset{¯}{𝝀}}\underset{¯}{\underset{¯}{\mathbf{}}}_0^{1/2}|,`$
where $`M`$ is the average number of entanglements per network strand ($`M=1`$ for an ideal unentangled nematic rubber) and the notation $`\mathrm{}`$ stands for an angular averaging of a matrix converted with an arbitrary unit vector . In a real elastomer $`M`$ is expected to be rather high since the reptation diffusion and the resulting disentanglement are restricted in a crosslinked network, unlike in the case of uncrosslinked melt where one needs a certain chain length ($`N_e`$) to form a topological knot . Clearly, at high $`M`$ entanglements dominate the rubber modulus, as is indicated by the second and third lines in the Eq. (2). However, it turns out that, although the modified elastic free energy appears different, many key parameters and predictions of the ideal theory remain in force.

One such equilibrium effect stands out in the properties of liquid crystalline elastomers, being especially pronounced when a monodomain (single crystal, permanently aligned ) nematic network is prepared. It is the anomalous thermal expansion – the spontaneous elongation of the material along the director axis on cooling from the isotropic phase . It is rather straightforward to derive the equilibrium uniaxial extensional strain by minimisation of the Eq.(1) where the chain anisotropy changes from $`\mathrm{}_{}^{(0)}=\mathrm{}_{}^{(0)}=b`$, in the initial state, to a set of values $`\mathrm{}_{}\mathrm{}_{}`$ with $`\mathrm{}_{}`$ along the nematic director $`𝒏`$, which is also the principal axis of spontaneous extension (see  for detail). The result for the uniaxial expansion,

$$\lambda _{\mathrm{th}}=(\mathrm{}_{}/\mathrm{}_{})^{1/3},$$
(3)
is a sensitive function of temperature. Remarkably, the same expression for $`\lambda _{\mathrm{th}}`$ is obtained from the modified reptation theory of nematic elastomers (2). The reason for this spontaneous uniaxial thermal expansion is the direct coupling between the average chain anisotropy and the nematic order parameter. Depending on this coupling, the magnitude of $`\lambda _{\mathrm{th}}`$ could differ greatly (and even change sense from extension to contraction in a polymer with oblate backbone conformation), reaching a maximum in a main-chain nematic polymer systems which could expand, spontaneously and reversibly, by a factor of $`34`$ .

In this paper we focus on the role of the single model parameter in the above molecular theories – the effective chain anisotropy $`\underset{¯}{\underset{¯}{\mathbf{}}}`$. On the one hand, the strength of the theory is that the effective anisotropy parameter $`r=\mathrm{}_{}/\mathrm{}_{}`$ could be directly measured by examining the thermal expansion of a uniaxial nematic rubber, cf. Eq.(3), thus leaving no free parameters and allowing quantitative predictions. On the other hand, one appreciates that the value of $`r`$ measured in a macroscopic experiment, such as the mechanical deformation, is an effective average parameter of the whole network, not necessarily directly related to the persistent length anisotropy of a particular polymer forming the network.

We examine a sequence of materials, having essentially the same chemical structure and composition of side chain polysiloxane nematic polymer strands but different in the type of permanent chemical crosslinking. Maintaining the constant degree of crosslinking, at 10% (by reacting bonds), we continuously alter the proportion between the two crosslinking agents: (i) flexible di-alkeneoxybenzene units, which are thought to have only a minor effect on the anisotropy of nematic network strands, and (ii) long di-functional chains of main-chain nematic polymer, which create an additional (and very high) anisotropy in the composite material. In all cases we prepare uniformly aligned monodomain nematic networks – single-crystal liquid crystal elastomers (SCLCE) in terminology of . By measuring the spontaneous uniaxial thermal expansion $`\lambda _{\mathrm{th}}`$, we obtain the temperature-dependent effective network anisotropy $`r(T)`$. A parallel thermal study (DSC) identifies the points of of nematic phase transitions for all materials and correlates these with the thermo-mechanical response.

We also correlate the thermal expansion (and the effective mechanical anisotropy $`r`$) with the local nematic order parameter $`Q(T)`$, measured by X-ray diffraction. It is obvious that there has to be a direct relationship, with the limit $`r1`$ at $`Q0`$. However, the theoretical derivation of such a correspondence is difficult because it is based on the single-chain analysis and is strongly model-dependent, varying from a trivial $`r=(1+2Q)/(1Q)`$ for a freely-jointed chain model to a highly non-linear $`r\mathrm{exp}\left[3/(1Q)\right]`$ in the hairpin regime of main-chain homopolymers . The results show that a linear relationship between $`r(T)`$ and $`Q(T)`$ holds over a broad range of temperatures.

Finally, we conclude by discussing the possible application of spontaneous deformation of monodomain nematic rubber in thermo-mechanical actuators. Section 4 presents a study of a selected optimised sample, providing the data for deformation under load (the artificial muscle regime ) and for the stress-temperature relation at fixed geometry (the dynamic actuator regime). There are competing factors, for instance, the amplitude and the steepness of the response are in the opposite relation to the response rate. Comparing the materials that span the wide range of possibilities, we propose that highly anisotropic main-chain containing nematic elastomers are most suitable for large-amplitude actuation when the time is a less relevant factor, while the low-anisotropy nematic materials provide the low-amplitude but high-frequency response suitable, for instance, for acoustic applications.

## 2 Experimental

All starting materials and samples of side chain siloxane liquid crystalline elastomers were prepared in the Cavendish Laboratory by modification of the procedure of Finkelmann et al . The polymer backbone was a poly-methylhydrosiloxane with approximately 60 Si-H units per chain, obtained from ACROS Chemicals. The pendant mesogenic groups were 4’-methoxyphenyl-4-(1-buteneoxy) benzoate (MBB), as illustrated in Figure 1, attached to the backbone via the hydrosilation reaction. All networks were chemically crosslinked via the same reaction, in the presence of commercial platinum catalyst COD, obtained from Wacker Chemie, with di-functional crosslinking groups also shown in Fig. 1 (both synthesized in the house). In all cases the crosslinking density was calculated to be 10 mol% of the reacting bonds in the siloxane backbone, so that on average each chain has 9 mesogenic groups between crosslinking sites.

However, the samples were different in the proportion of two crosslinking components. The first crosslinking agent, 1,4 di(11-undecene) benzene (11UB), is a small flexible molecule deemed to have relatively minor effect on the overall mesogenic properties of the liquid crystalline polymer (apart from the trivial effect of added impurity, which is deliberately kept constant in our experiments). The second crosslinking agent was very different: di-vinyl terminated chains $`\alpha `$-{4-\[1-(4’-{11-undecenyloxy} biphenyl)-2-phenyl\] butyl)-$`\omega `$-(11-undecenyloxy) poly-\[1-(4-oxydeca methyleneoxy)- biphenyl-2-phenyl\] butyl chains that themselves make a main-chain liquid crystal . In our case the crosslinking chains were with $``$74 rod-like monomers between the terminal vinyl groups (determined by GPC, polydispersity $``$2). Each such long-chain crosslinker connects two siloxane polymer chains.

Calculating the crosslinking density by reacting bonds, as described above, can be quite different from the actual concentration of the species in resulting material. In particular, in the case of MC polymer crosslinker, the relative gram weight of it in the otherwise side-chain polymer matrix is very high: the molecular weight of an average MC crosslinking chain, cf. Fig. 1, is $`3453075\times 456`$ – which is much greater than the molecular weight of the corresponding 9-unit side-chain polymer strand between crosslinking points ($`26919\times 299`$). One may equally regard such a system as a main-chain nematic rubber network, end-linked with relatively small crosslinking groups made of side-chain nematic polysiloxane.

The series of materials, synthesized in this way, are referred to by the proportion of MC crosslinker in the total crosslinking composition. In this way, the sample labeled as MC0 is crosslinked purely by the 10 mol% of small-molecule group 11UB. At the opposite end, MC10 refers to the material crosslinked purely by the long main-chain nematic polymer. We find particularly interesting effects occurring at very small proportion of main-chain crosslinker. The borderline case, when the overall mass of MC material in the network composition is approximately equal to that of side-chain polysiloxane, is at approximately 1.75 mol% of MC. We find that compositions close to this are optimal in balancing the increasing amplitude of thermal expansion and the decreasing rate of response. The samples we have studied, MC0.1, MC0.25, MC0.5 and MC1, refer to the crosslinking composition listed in the Table 1, together with their characteristic transition temperatures. For comparison, the pure uncrosslinked melt of the same main-chain polymer, studied in some detail in , has $`T_g`$38 C and $`T_{\mathrm{ni}}`$ 112 C (for the longer-chain MC component).

Monodomain, aligned samples were made by preparing partially crosslinked films in a centrifuge, highly swollen in toluene (2-3 ml per 1 g of material), reacting for 25-35 minutes before evaporating the solvent and suspending the samples under load in an oven to complete the reaction for more than 5 hours at 120 C. A careful study of reaction kinetics ensured that approximately 50% of crosslinks were established in the first stage of this preparation. When a uniaxial stress is applied to such a partially crosslinked network, the aligned monodomain state in the resulting nematic elastomer is established with the director along the stress axis. This orientation is then fixed by the subsequent second-stage reaction, when the remaining crosslinks are fully established. Following the ideas and results of , in all cases we performed the second-stage crosslinking in the high-temperature isotropic phase: in this way a better alignment and mechanical softness are achieved (in contrast to the crosslinking in a stretched polydomain nematic phase , which results in a number of defects and domain walls frozen in the material).

Equilibrium transition temperatures given in the Table 1 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. 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)`$. Precise length measurements were made using a traveling microscope, as the samples were slowly cooled (at a rate 0.5o/min) after the complete annealing in the isotropic state.

The force measurements were performed on a custom built device consisting of a temperature compensated stress gauge and controller(UF1 and AD20 from Pioden Controls Ltd) in a thermostatically controlled chamber (Cal 3200 from Cal Controls Ltd). The rectangular shaped samples (approximately $`30\times 5\times 0.3`$ mm) were mounted with clamps, at room temperature, in such a way that the sample length $`L`$ remains fixed throughout the experiment. The samples were then heated (at a rate $`0.5^\mathrm{o}`$/min) to the isotropic phase. The stress and temperature values were acquired by connection to a Keithley multimeter (2000 series) and stored on a PC over an IEEE interface. The increasing force was measured by the stress gauge and later correlated with the effective strain in the sample, due to the decreasing underlying natural length $`L(T)`$. Since the effective strain reached high values, $`\lambda _{\mathrm{th}}33.5`$ when $`T_{\mathrm{ni}}`$ was reached, many samples broke or tore. However, a good reproducibility of force-temperature, as well as stress-strain relations allowed one to extrapolate the data over the whole temperature range.

## 3 Thermal and mechanical equilibrium

### Phase transitions

Figure 2 presents the data of DSC measurements of phase transformations in our materials. All samples exhibit a single nematic phase above their glass transition, however, the position and the morphology of phase transformations strongly depends on the composition. It is expected, and well-known, that liquid crystalline polymers show a much broader and diffuse signature of the weakly first-order nematic-isotropic transition. This effect is especially pronounced in crosslinked elastomers, where the quenched sites of random disorder introduced by crosslinks lead to a more diffuse transition and stronger hysteresis. Fig. 2(a) shows an example of heating and cooling scans for a selected sample, MC0 and MC10, at the opposite ends of the series. One can clearly recognize the effect of main-chain polymer component, which increases the glass and the nematic-isotropic transition temperatures.

Figure 2(b) focuses only on the region of nematic transition and brings together the DSC scans for all materials, obtained at the same cooling rate of 10 C/min after a long annealing at 120 C. Comparison between MC10 and MC1 tells that even at the relative proportion of MC crosslinking chains of 1% the properties of the resulting network are dominated by their effect. This observation will be supported below by the comparison of spontaneous thermal expansion of MC10 and MC1, which has the same overall character and the transition point, but rather different from the other members of the series. One should expect this, since the calculation of the relative mass of the side- and the main-chain components in our materials becomes approximately equal at the level of 1.75 mol% of MC (see Table 1).

The second key result of this DSC study, evident in all materials between MC0.1 and MC1, is the presence of two sequential nematic transitions. We have verified by X-ray and optical observations that no additional phase exists in the materials, which are optically transparent and homogeneous at all temperatures. Clearly, the higher-temperature transition is related to the main-chain nematic component of the network: its temperature hardly shifts on composition (remember that the data in Fig. 2(b) is obtained on cooling after the complete annealing), but its overall latent heat decreases in proportion with the decreasing percentage of MC. The second, lower-temperature transition is, of course, that of the side-chain polysiloxane polymer. Its temperature continuously shifts up on increasing proportion of the MC crosslinker and, accordingly, the stronger nematic field provided by the more strongly anisotropic MC polymer (which has already undergone the phase transformation, as indicated by thermal scans). We shall now investigate the related mechanical response and return to the discussion of the possible reasons for such a behavior.

### Uniaxial thermal expansion

Changing the temperature affects the nematic order parameter $`Q(T)`$ of a nematic elastomer and, in a permanently aligned monodomain state, alters the effective chain anisotropy of the rubber, $`r=\mathrm{}_{}/\mathrm{}_{}`$. Figure 3 shows the variation of natural length of free, unloaded samples, measured along the nematic director. Obviously, in an incompressible rubber the volume is conserved and the two perpendicular directions experience the symmetric contraction by $`1/\sqrt{\lambda }`$ when the principle direction is elongated by $`\lambda =L/L_0`$. Here the length $`L(T)`$ is measured with respect to the length $`L_0`$ the corresponding sample has in the isotropic phase. Elongation on cooling into the nematic phase means that the network is characterised by an effective prolate anisotropy, that is $`r>1`$, below $`T_{\mathrm{ni}}`$. The spontaneous uniaxial thermal expansion described by the strain $`\lambda _{\mathrm{th}}`$, Eq.(3), is an equilibrium reversible process. Fig. 3(a) illustrates this by showing the evolution of sample length for one selected material, MC0.25, on heating and on cooling. There is an evident hysteresis which, however, is rather small at the cooling rate used in this experiment. Without going into greater detail, we must mention that the characteristic response time, or equivalently – the width of hysteresis, increase dramatically with the increasing MC concentration. In the material MC0 we find no trace of a hysteresis, while in MC10 the response to changing temperature is rather slow and the hysteresis in $`L(T)`$ is very noticeable at an already slow cooling/heating rate of $`0.5^\mathrm{o}`$/min. This slow response and relaxation of main-chain nematic polymers , the likely reason being the freezing of chain motion in narrowly confining hairpins.

Analyzing the sequence of $`\lambda _{\mathrm{th}}(T)`$ data for the full series of samples, one is impressed by the high value of strain achieved, especially in the MC-containing materials. The arguments already presented above, about the relative mass of the long-chain MC material in the overall network composition becoming approximately equal to that of the siloxane side-chain polymer at the level of $``$MC-1.75%, are confirmed by the observation that no significant further increase in the effective network anisotropy $`r=\lambda _{\mathrm{th}}^3`$ occurs in high-MC containing samples; the increase in $`\lambda _{\mathrm{th}}`$ at a given $`T`$ in going from MC1 to MC10 is much less than from MC0 to MC1. One could argue that the ambient-temperature value of $`r42`$ found in MC10 is close to the natural anisotropy of our MC polymer. This would make a prediction for the radius of gyration anisotropy $`R_{}/R_{}6.5`$ (since $`R^2\mathrm{}L`$), which is a high but not unexpected value for this type of semiflexible polymer, as seen by neutron scattering from the melt (see, e.g.,  and references therein).

Another remarkable aspect of the results shown in the Fig. 3(b) is the dual-transition feature of the transition, similar to that seen in the DSC study. Networks with the low MC-concentration, when the main-chain component is the minority in the composition, show their first sign of nematic transition (indicated here by the initial small uniaxial expansion) at the higher $`T_{\mathrm{ni}}`$, followed by a subsequent rapid increase at the lower $`T_{\mathrm{ni}}`$. The amount of average network anisotropy generated at the first transition increases with the concentration of MC, and is consistent with the overall latent heat of this higher-temperature transition in the Fig. 2(b). The second transition, seen as a near-critical jump in $`\lambda _{\mathrm{th}}`$ for MC0.1, moves towards higher temperatures as the MC-concentration increases, again consistently with the DSC data.

One fundamental question that arises from the results shown above is the nature of the first, higher $`T_{\mathrm{ni}}`$, nematic transition. A natural explanation would be the phase separation of very long MC chains during the synthesis, when they are in the minority composition. It would then be natural that small volumes with the high concentration of MC undergo the nematic transition at $`108`$ C. The high uniaxial elongation of these small MC-regions would, however, be absorbed by the still isotropic network of polysiloxane chains and the overall mechanical strain would be small. This would also explain the shift of the $`T_{\mathrm{ni}}`$, because the uniaxially aligned nematic MC-regions would provide a driving field for the rest of the mesogenic material. However, it is not straightforward to reconcile this hypothesis about the phase separation with the fact that our materials show no optical inhomogeneity (which could simply mean that the MC-regions are too small), as well as with the fact that each MC-chain has to be crosslinked at both ends to a siloxane strand. Since the first-step crosslinking takes place in the presence of solvent in a well-mixed centrifuge environment, one expects a relatively homogeneous distribution of junction points, which should then be preserved after the second crosslinking stage, in spite of de-swelling. An intriguing alternative possibility could be that we observe the precursor to the nematic transition in essentially separate chains of MC. Theoretical studies  have long ago predicted the existence of nematic-nematic phase transition lines when two competitive ordering trends are present in the system; in that terminology, we could be observing the discontinuous $`N_{\mathrm{III}}`$-$`N_{\mathrm{III}}^{}`$ transition at a temperature below $`T_{\mathrm{ni}}`$ for the change between the isotropic and the $`N_{\mathrm{III}}`$ phases. We have not pursued this question in any more depth here, instead concentrating on the mechanical aspects of nematic ordering in rubber.

## 4 Thermo-mechanical actuation

An important practical question of spontaneous thermal expansion of monodomain nematic rubbers, relevant to practical use of this effect in thermo-mechanical actuators, is the amount of stress the material can sustain. The work cycle of a corresponding artificial muscle depends on the amplitude of motion, but also on the force the nematic rubber exerts – or the load it can work under. In this section we focus on the sample MC1, which shows an almost maximal strain $`\lambda _{\mathrm{th}}(T)`$, but still has a reasonably fast response and small hysteresis.

Figure 4(a) shows the difference in thermal expansion curves when the sample was loaded with increasing weights. For the sample of cross-section area $`1.5\text{mm}^2`$ each 5 g of weight corresponds to the stress $`\sigma 30`$ kPa; the curves in the plot are for the stress increasing from $`0`$ to $`90`$ kPa. Clearly, the qualitative nature of the effect is preserved when changing the temperature under load. An increase in the strain values is the natural consequence of rubber elasticity: at any given temperature and natural length, the applied load would cause an extra stretching of the rubber band. Several interesting effects are seen in this plot, in particular, the small but unambiguous shift of the nematic transition temperature, $`\mathrm{\Delta }T_{\mathrm{ni}}5^\mathrm{o}`$, with increasing load. Such an effect is expected on theoretical grounds in a weak first-order transition system in an external field – however, in ordinary liquid crystals, electric and magnetic fields are usually too small to create a noticeable shift. On the other hand, there is no trace of supercritical behaviour or paranematics phase above $`T_{\mathrm{ni}}`$ : the transition remains abrupt and well-defined; the small increase of $`\lambda `$ above $`T_{\mathrm{ni}}`$ is purely due to the rubber-elastic response to the applied load.

Another practically important experiment is reported in Fig. 4(b). Here we measure the force that a nematic rubber band can exert when thermally contracting. A strip of MC1 elastomer with a cross-section of $`5\times 0.3`$ mm (standard in all experiments described in this paper) is rigidly clamped in a dynamometer device at ambient temperature, so that its length $`L`$ remains constant. On heating, the underlying natural length decreases, $`L_0(T)`$, according to the data of the previous section, cf. Fig. 3(b) for MC1. Accordingly, the effective extensional strain increases in the material and the sample exerts an increasing force on its clamped ends. This force is measured and transformed into the nominal stress by dividing by the sample cross-section area.

The plot in Fig. 4(b) shows two separate experiments on heating the constrained rubber band MC1, circles and diamonds data points (the reproducibility of results is comforting). In both cases the sample broke before reaching the isotropic phase – clearly, the stress of the order of 100 kPa has been the limit this particular material could bear (one could conceivably increase the strength by tailoring the sample shape and in particular, its thickness). We performed a number of such measurements and the results are very reproducible for a given material and sample shape. The second $`y`$-axis on this plot gives the variation of the underlying effective strain. We calculate it determining the actual gradient of underlying displacement with respect to the “virtual” equilibrium length, $`\epsilon =[L_{\mathrm{fix}}L(T)]/L(T)`$, where $`L_{\mathrm{fix}}`$ is device-fixed, while $`L(T)`$ follows the data measured and reported in the previous section. One can clearly see that the internal stress and the effective strain follow exactly the same graph and, therefore, are linear functions of each other over a large region of deformations.

## 5 Conclusions

In this article we report the systematic study of uniaxial thermal expansion of monodomain nematic elastomers. Composite materials combined main-chain and side-chain nematic polymers in their networks and the results are sensitive functions of concentration of the MC component. Increasing the proportion of MC in the otherwise side-chain nematic network has a great effect on the average effective anisotropy of polymer chains and, as a result, on the magnitude of spontaneous strain $`\lambda _{\mathrm{th}}`$. At the same time, there is a pronounced slowing down of all response processes: the stress relaxation, the equilibration of natural length and the shape recovery. One does expect such slow dynamics in MC-dominated materials, which is consistent with other studies of main-chain nematic polymers. The isotropic-nematic transition at $`T_{\mathrm{ni}}`$ is abrupt in all our materials, with no supercritical effects or paranematic phase. However, both the mechanical and the thermal data indicate that at a small concentration of MC, the network experiences a sequence of two nematic transitions. A possible explanation for this effect could be a microphase separation of MC-rich regions, however, there factors against it and the question is far from being fully understood.

The study of mechanical actuation, that is the contraction under load, or exerting a force on changing the temperature shows that uniformly aligned nematic elastomers have a great potential in practical applications, from stress and deformation gauges, to artificial muscles and micromanipulators. On the map of material properties of different materials used in actuator design , the nematic elastomers would occupy the region of low stress (unlikely to exceed a few MPa) and extremely high deformations. In comparison, the present “large-amplitude” actuators are based on shape-memory alloys, which reach strains of up to 10% – something that pales in comparison with over-300% deformations described in this work.

### Acknowledgements

This work has been supported by EPSRC UK. We are grateful to S.M. Clarke, M. Warner and A. Hotta for a number of useful discussions and help with the measurements. Special thanks to H. Finkelmann and A. Greve for advice and help with material synthesis procedures.


format: text

text:
# Calculations of neutrino spectra with a program based on the GEANT library.

## 1 Introduction

In 1988 - 1989 the new generation of neutrino experiments started at IHEP (Protvino). It was a result of the physical start of the new IHEP-JINR neutrino detector, increased intensity of the IHEP 70 - GeV proton accelerator and somewhat changed situation in the neutrino physics, which required better accuracy of measurements and inspired interest to experiments with electron neutrinos.
Most of these experiments were carried out with a nonstandard mode of work of the neutrino facility. It concerns the beam-dump experiment, the experiment with focusing system switched off and the new experiment with short decay base. These experiments required a new program for neutrino spectra calculations which could take into account additional effects that become important in these modes of the neutrino beam formation, for example hadronic interactions in the elements of the neutrino channel other that the neutrino target.
For this purpose a program based on the GEANT library  was chosen. It has the following advantages:

* developped, multiply checked and supported package for a simulation of passage of particles through matter, including hadron showers;
* the way of describing a geometry is sufficiently developped and familiar to most of physicists.

The disadvantage of the full GEANT - like simulations of processes in the neutrino channel is that it is rather time consuming. This required usage of tricks that would increase the number of useful decays, but in spite of them, only after appearance of sufficiently powerful computers an accurate calculation of neutrino spectra by such a program became possible.

## 2 The program for the neutrino spectra calculation

In the GEANT library the procedures for random decays of particles in flight are provided. If GEANT accomplishes a decay, the particle is stopped (no more traced). Thus, a particle can decay only once, which required the simulation of many particles and a lot of CPU time. For this reason routines different from the ones of the GEANT library were used for the simulation of decays. These routines are called several times on every step of a particle. The GEANT table of branching ratios was changed (some branching ratios are zeroed), so that GEANT does not accomplish decays to neutrinos of the type for which the spectra are being calculated. To take into account correctly the zeroed branching ratios the GEANT routine GDECAY is slightly changed. Namely, if no decay mode is chosen (it is now possible since the sum of branching ratios is not equal to unity), a particle is not stopped. Thus, a particle makes more steps, which also improves the statistical accuracy of the calculations. The effect of decreasing number of particles along the beam due to a finite lifetime is taken into account through weights ascribed to particles.
Unfortunately, although all these tricks certainly decrease the overall statistical errors of calculations, they do not decrease by the same factor all components of these errors. For example, the fluctuations due to hadron showers are not decreased. This makes very difficult a direct estimation of statistical uncertainties of spectra. The only obvious way to estimate them is to compare results of several runs of the program. This consideration also suggests that there is no sense in very small steps or many decays in one step.
Some branching ratios, such as the branching ratios to hadrons, are not changed in order to simulate more exactly the beam composition (for example if a kaon decays to pions, these pions are also traced later). For the same purpose and in order to take into account decays of muons during the calculation of electron neutrino spectra the program has 2 modes of work:

1. The calculation of muon neutrino and antineutrino spectra. The branching ratios of decays to electron neutrinos are not changed since in some of them charged pions are produced.
2. The calculation of electron neutrino and antineutrino spectra. The branching ratios of decays to muons are not zeroed.

In each of these modes it is possible to simulate directly hadron showers in the neutrino target induced by initial protons or (not for the beam dump) to use parameterized yields of particles from the standard neutrino targets.
As most of programs using GEANT library the program consists of the following parts:

1. The initialization. Here beside usual GEANT calls the GEANT table of branching ratios is changed.
2. The geometry definition. It is a usual GEANT geometry definition  except that some volumes should have standard numbers (the target, lenses and some volumes for the calculation of control distributions). Besides, the constants for the spectra calculations, such as a distance to a possible detector and its maximal sizes, are set.
3. The routine for the action before each event. Here either a proton is put into the GEANT initial bank (direct calculation) or several hadrons with random momentum and angle (parameterization).
4. The routine for the action after each step. It calls the two main routines of the program:
a) The routine that ascribes a weight to a new particle or calculates a weight after each step. A weight 1 is ascribed to primary protons. In the mode of work with a parameterization a weight to primary hadrons at their first appearance is ascribed according to this parameterization. For other particles, at their first appearance, their creation vertex is searched for in the array of vertices (this is not in the standard GEANT banks and common blocks, see below). From this array the initial weight, the type and momentum of a parent particle, the spin direction (for a muon) are taken. If there is a correction factor for the process of creation of this particle by its parent particle the weight is corrected. After each step the weight is decreased according to the life time of particle divided by the sum of all zeroed branching ratios.
b) The calculation of the neutrino spectra. For two body decays the limits on the angles in the center of mass system are calculated, so that neutrinos from a decay can be inside the maximal detector volume. For three body decays, if a neutrino energy is above threshold and its direction does not hit the detector maximal volume, its azimuthal angle is resimulated withing the limits that contain the maximal detector volume. The probability of a decay in case of additional limits is multiplied by the corresponding phase space factor. The matrix element of decays is taken into account. For decays of muons it depends on the spin direction.
The probabilities of decays are put into double precision arrays that correspond and will be rewritten in the end of run into one- and two-dimensional histograms (energy - radius, energy - angle).
5. The routine for the action after each act of creation of new particles. As usual, particles below cut and those that do not contribute to the neutrino spectra, are ignored. Here the array of vertices mentioned above is filled

## 3 The spectra of hadrons.

The first calculations of WBB spectra were made in the mode with the full simulation of the hadron shower in the neutrino target, which is a 60 cm long aluminium cylinder. Although the final calculations were made using the parameterization, it is interesting to compare the GEANT yields with the ones measured in . These measurements are rather detailed, they include secondary particles from the primary beam of 67 - GeV protons in the momentum range from 5 - 7 GeV to 60 GeV and the angle range from 0 to 50 mrad, subdivided into energy bins of 1 GeV and angle bins of 2 mrad.
The comparison for GEANT 321 with GHEISHA is shown in Fig. 1.
The comparison for GEANT 321 with FLUKA is shown in Fig. 2.

From these fugures one can see that GEANT does not simulate correctly the yields of hadrons, the difference reaching a factor of 2 for the case of GHEISHA. For this reason in the mode of work with a full simulation the correction by weights (previous section, 4a) should be used.
Special attention should be paid to the spectrum of $`K^0`$, for which there are no reliable measurements at 70 GeV. Often it was assumed that there is the following relation between the charged and neutral K mesons:

$$K_L^0=K_S^0=\frac{K^++K^{}}{2}$$
(1)
It is clear from the following considerations that the relation should be different. We are usually interested in spectra of rather fast particles, which means that they should be ”leading” particles, i. e. should contain a quark from the incident proton. Simple counting of quarks shows that the spectrum of $`K^0`$ should be closer to the spectrum of $`K^{}`$ than to the spectrum of $`K^+`$.
In order to find the best description the following measurements were analysed: the $`K^0`$ measurement at FNAL  at 300 GeV, the charged kaon measurements at CERN at 400 GeV  and the charged kaon measurements  at IHEP (Protvino) at 67 GeV. All measurements were made on thin targets. The method  of comparison of spectra at different energies was used. It uses the variable $`x_R=E^{}/E_{max}^{}`$, the ratio of the energy in the center of mass system to the maximal energy at a given $`p_t`$. The best description was obtained with the following formula:

$$K_L^0=K_S^0=\frac{K^++(2n1)K^{}}{2n}$$
(2)
where $`n`$ is the ratio of the distributions u/d in a proton. This value is a function of a scaling variable $`x`$, which we substitute here by $`x_R`$, which does not introduce a large error since the function does not change rapidly with $`x`$. It can be found for example in . $`n(x_R)`$ changes from 1 to 5 as $`x_R`$ changes from 0 to 1 (Fig. 3). The following approximation was used afterwards in calculations of spectra:

$$n=1.435x_R^3+5.3883x_R^2+0.22269x_R+1.0189$$
(3)
The uncertainty of formula (2) is estimated to be 15%, which is a quadratic sum of the uncertainty of the FNAL measurements (10%), the uncertainty of the scaling from 70 to 400 GeV (10%) and a smaller uncertainty of the CERN measurements.
In Fig. 4 one can see that the formula (2) is in a much better agreement with the copper data than $`(K^++K^{})/2`$. The agreement of the formula (2) on berillium is also quite reasonable. Apparently the measurements  on berillium slightly underestimate the yields.

## 4 The beam dump neutrino spectra.

Practically the only possible way of spectra calculation for beam dump experiments is the mode with a full simulation of hadron showers in the target. It is difficult here to use parameterizations. However, at least the GEANT (GHEISHA) yields of particles from the first interaction of protons were compared with the results of measurements on the thin copper target . The difference was similar to that seen in Fig. 1. A corresponding correction was used in the calculations of spectra (see section 2, 4a).
The systematic uncertainties of the calculated spectra were determined by the systematic uncertainties of the measurements . Possible uncertaintied in the simulation of the secondary interactions by GEANT were suppressed since corresponding secondary particles contribute 20 - 30% to the neutrino spectra interesting for the beam dump experiment. These contributions were so small due to the fact that the experiment  was devoted to the search of hard neutrinos from the charm production which required rather high energy cuts.
It is possible to reduce these uncertainties for the electron neutrino and antineutrino spectra. For this purpose the reconstructed muon neutrino events were used. The high energy region of the muon neutrino and antineutrino spectra is dominated by decays of charged kaons (the contribution from charm estimated from the beam-dump extrapolation method was small). For this reason, comparing the number of reconstructed muon neutrino events with the calculated one it is possible to check the charged kaon spectra and to reduce systematic uncertainties.
The final accuracy of the spectra calculations were 10% for the electron neutrinos and 16% for the electron antineutrinos.
Besides, for the beam-dump experiment such important value as a ratio of effective densities of the two targets can be calculated with a very good accuracy. Since the dimensions of the half-density target were finite and there were constructive gaps, it is not exactly 2 as it was designed. Furthermore, this ratio for muon neutrinos is different from that for electron neutrinos since for the latter there is a contribution from $`K_S^0`$. For such calculations mainly a correct simulation of the shower shape is needed, which is reproduced by GEANT more reliably (this was checked, for example, during the calibration of the neutrino detector). The calculated ratio of the effective densities is 2.0 for muon neutrinos and 1.9 for electron neutrinos (accuracy 1%, mainly statistical).
For additional control of the beam and spectra during the beam-dump experiment the muon fluxes in 9 gaps of the muon filter were measured . In Fig. 5 the calculated fluxes are compared with the measured ones (the accuracy of measurements 3%, to which the uncertainty of the prompt muon pair production subtraction 7% should be added). Unfortunately, such control does not allow to improve the accuracy of the electron neutrino spectra calculations since the main source of muons (more than 80%) is the decay of pions.

## 5 The neutrino beam with the focusing system turned off.

Two rather long runs of the IHEP-JINR neutrino detector were carried out with the neutrino channel focusing system turned off. Such mode of work has the advantage of a higher fraction of electron neutrinos, sufficient to perform some studies. These runs were also used for the measurement of the total cross section of neutrino and antineutrino interactions (simultaneously). The corresponding total neutrino fluxes are significantly lower than with the focusing system turned on, but the loss of statistics for the neutrino detector was not dramatic because due to the lack of a resonant beam extraction (only a fast extraction is available in IHEP Protvino) the event aquisition rate is limited by the maximal allowed occupancy of the neutrino detector. In these measurements it is also important that neutrino and antineutrino spectra are less steep, which decreases the effects of smearing.
Since in this neutrino beam (as well as in all the beams described below) the standard neutrino target was used, the calculations were accomplished in the mode of work of the program that uses the parameterization of secondary particle spectra from the target. But the data of  could not be used directly since the measurements  were carried out with a 6 cm diameter aluminuim target, whereas in the neutrino runs a 1 cm diameter target was used. The neutrino spectra in this case are better (harder), but the parameterization should be corrected. The correction was obtained in the following way: the program was run in the mode of full shower simulation and the correction was obtained by dividing the corresponding yields of hadrons from the targets.
In the mode of work with the parameterization the GEANT correction procedure mentioned in section 2, 4a, is still needed since part of protons traverse the target without interaction and interact in the lenses, which is simulated by GEANT. What is significantly better in this mode of work is that statistical errors at high energy are much lower since the unweighted spectrum of hadrons from a target to be traced by GEANT is flat (is not a function of momentum).
It turned out that without focusing the secondary effects, such as interactions of hadrons in the channel elements are much more important than with normal focusing. This happens because a bigger fraction of particles hit the walls of the decay tunnel. In addition, when the focusing system is turned on, the secondary particles from the interactions in the neutrino channel elements are not focused or are not in optimal conditions for a focusing. For this reason they contribute less to the neutrino spectra.
Only after calculations with a program based on GEANT the agreement between the calculated and measured  muon fluxes was obtained. This is shown in Fig. 6.

The systematic errors of the muon neutrino spectra in the first approximation are those of , but they are additionally controlled, mainly below 10 GeV (pionic part), by the muon flux measurements (the accuracy of the measurements about 3%).
In the electron neutrino spectrum for this neutrino channel a contribution from muon decays is significant. For the calculation of this contribution it is important to take into account the polarization of muons. These decays increase the electron neutrino spectrum by up to 20% at 3 GeV (thay give softer neutrinos so this difference is the smaller the higher is the energy). This effect is even more important for the neutrino beam with focusing.

## 6 The neutrino beam with the focusing system on.

Here some difficulty was to introduce into the program the parabolic lenses used in IHEP because there is no paraboloid shape in GEANT. If define these lenses as cylinders, the magnetic field in them is not only non-uniform but non-contiguous and GEANT does not trace particles correctly. For this reason a routine for the step limitation in the lenses was introduced. The step is limited so that it is either entirely in a region with field or entirely in a region without field or is very small.
The electron neutrino spectrum with a contribution from decays of muons is shown in Fig. 7.
The systematic uncertainty of the neutrino flux calculation consists of the same uncertainties as for the beam with the focusing system off and additional uncertaintied of the focusing system such as a possible inaccurate target and lenses positioning and an uncertainty in the simulation of starting points of particles inside the target. No studies of the focusing system uncertainties were made with this program, but the previous studies  show that, taking into account muon flux measurements, they can be kept as low as 7%

## 7 The neutrino beam with a short decay base.

For this experiment , with rather complicated and asymmetrical geometry of the decay region, it was very convenient to use GEANT tools for the description of the geometry.
The decay cavity can be described as 4 adjacent (with steps in the vertical projection) trapezoids, asymmetric in the vertical projection. It is narrow near the target to provide a better radiation shielding and wider near the end closest to the detector. It was made asymmetric because of the danger of ”underground” muons that can go around the iron muon filter, in the earth under it due to a multiple scattering.
In this experiment a contribution from muon decays is also rather significant, including decays of muons in the muon filter (the path length of muons in it is comparable or even bigger than the length of the decay region).

## 8 Conclusion

A program for the neutrino spectra calculations based on the GEANT library is created. The range of possible applications of the program is very wide.

## 9 Acknowledgements

The author expresses the gratitude to A.S.Vovenko, A.I.Mukhin, Yu.M.Sviridov, Yu.M.Sapunov for helpful discussions, V.I.Kravtsov, Yu.M.Sapunov, K.E.Shestermanov for provided routines.
The developments of the program for the neutrino spectra calculations have been possible thanks to the computer power and other support offered by the Physics institute of the University of Urbino (Italy).


text:
# Nucleon momentum distribution in deuteron and other nuclei within the light-front dynamics method

## I Introduction

The most recent experimental data on the deuteron structure functions and tensor polarizations obtained at Thomas Jefferson Laboratory have been reported in  (and references therein). The new data on tensor polarization observables measured in elastic electron-deuteron scattering for the squared momentum transfer between 0.66 and 1.7 (GeV/c)2  make it possible to determine quite precisely the deuteron charge form factors for momenta comparable with the deuteron mass. Thus these data are important to probe the relativistic dynamics inside the deuteron. It was shown in  that the empirical results for the tensor polarization $`t_{20}`$ and both monopole ($`G_C`$) and quadrupole ($`G_Q`$) charge form factors are in very good agreement with the calculations within two relativistic and covariant models: i) a model developed in the framework of the explicitly covariant version of light-front dynamics (LFD) (e.g. ) and, ii) a model using a three-dimensional reduction of the Bethe-Salpeter equation . At the same time the non-relativistic impulse approximation approach  (even with inclusion of meson exchange currents and relativistic corrections) cannot explain, for instance, the behaviour of the monopole charge form factor in the region of the first node and the secondary maximum shown for the first time in .

It is well known  that one of the main feature of the realistic nuclear models (beyond the mean-field approximation) is that they have to describe simultaneously both important ground state characteristics, the nucleon momentum and density distributions of the system. This is not the case in the Hartree-Fock approximation where the dynamical short-range and tensor nucleon-nucleon (NN) correlations are only partly incorporated. The main characteristic feature of the data on the nucleon momentum distribution $`n(q)`$ extracted by the nuclear $`y`$-scaling analysis  (we shall call them all along this paper $`y`$-scaling data (YSD)) is the existence of high-momentum components for momenta $`q>2`$ fm-1 due to the short-range and tensor NN correlations. This has been shown also for particular nuclei such as 2H, 3,4He, 12C, 16O, 40Ca and others and in nuclear matter within various theoretical correlation methods (see e.g.  and reviews in ). In addition, both experimental and theoretical analyses confirm the conclusion that the high-momentum behaviour of the nucleon momentum distribution ($`n(q)`$/A at $`q>2`$ fm-1) is similar for nuclei with mass numbers $`A`$=2, 3, 4, 12 and 56 and for nuclear matter.

The first aim of the present work is to apply the LFD method to calculations of the momentum distribution in the deuteron. Thus this study is an important complementary test of the possibility of the LFD to describe simultaneously both deuteron charge form factors (that has been shown in ) and momentum distribution. In this line the results of calculations using wave functions of existing NN potentials are also given and compared with the LFD results and the YSD.

The second aim of the work is to develop a method which would allow us to calculate the nucleon momentum distribution for a wide range of $`(A,Z)`$-nuclei on the basis of the knowledge of $`n(q)`$ in the deuteron and especially of its high-momentum components. The latter are important because in various processes with high transferred momenta only relativistic nucleons can take part . Such analyses require the knowledge of the wave function (WF) at momenta with magnitude of the nucleon mass. We make an attempt to relate the relativistic effects in $`n(q)`$ in the deuteron with the explanation of the high-momentum components of $`n(q)`$ in heavier nuclei which are important for calculations of the cross sections of various nuclear reactions. This question is discussed together with the effects of short-range NN correlations on the nucleon momentum distributions.

A brief review of the LFD method is given in Sect. II. The results of the calculations of $`n(q)`$ in deuteron within the LFD and the comparison with calculations using wave functions corresponding to existing NN potentials and the YSD are given in the same Section. In Sect. III a method to calculate the nucleon momentum distribution in $`(A,Z)`$-nuclei is developed using the natural orbital representation and the results on the high-momentum components of the nucleon momentum distribution in the deuteron. Sect. IV contains the conclusions of the present work.

## II Nucleon Momentum Distribution in the Deuteron Within the LFD Method

The results of the experimental and theoretical studies of the structure and form factors of the deuteron have been presented in the most recent review . In it special attention has been paid to the different theoretical models for investigating the deuteron electromagnetic form factors and nucleon momentum distribution. Among them, the LFD method has been discussed as one of the promising approaches. We start this Section with a brief description of the LFD (e.g. ) which is applied to the analyses of relativistic bound systems at high relative momenta. The LFD is a self-consistent method in which the relativistic wave functions are the Fock column components of the state vector defined on an arbitrary light-front (LF) hypersurface given by invariant equation $`\omega x=0`$. The four-vector $`\omega =(\omega _0,\stackrel{}{\omega })`$, $`\omega ^2=0`$, determines the position of the LF surface. The state vector satisfies a dynamical equation following from rotations of the light front. In LFD the transformations of the reference system and of the light front are independent of each other. The dynamical dependence of the state vector on the LF position is parametrized by the explicit dependence of the wave functions on the four-vector $`\omega `$. The dependence on $`\omega `$, however, is not a property of the observable amplitude. The form factors depend on four-momentum transfer squared only and do not depend on $`\omega `$, though they are related to the $`\omega `$-dependent wave functions . The relativistic deuteron WF on the LF $`\mathrm{\Psi }(\stackrel{}{q},\stackrel{}{n})`$ depends on two vector variables: i) the relative momentum $`\stackrel{}{q}`$ and ii) on a unit vector $`\stackrel{}{n}`$ along $`\stackrel{}{\omega }`$. Due to this, the WF is determined by six invariant functions instead of two ($`S`$\- and $`D`$-waves) in the non-relativistic case. Each one of these functions depends on two scalar variables $`q`$ and $`z=cos(\widehat{\stackrel{}{q},\stackrel{}{n}})`$. In LFD these six functions are calculated within the relativistic one-boson-exchange model. The kernel in the dynamical equation for $`\mathrm{\Psi }(\stackrel{}{q},\stackrel{}{n})`$ has been calculated  using Lagrangians of interaction of nucleon with pseudoscalar $`(\pi ,\eta )`$-, scalar $`(\sigma ,\delta )`$\- and vector $`(\rho ,\omega )`$ mesons using the parameters of the Bonn model . Thus in the LFD calculations the deuteron structure is determined by the relativistic nucleon-meson dynamics, supposing that the nucleons in the deuteron interact by exchanging relativistic mesons , without using potential approximations.

The relations of the covariant LFD method to other relativistic approaches, such as the standard LFD , the Bethe-Salpeter formalism  and its three-dimensional reductions, the quasi-potential equations  and the Gross equation method  are given in . For instance, it is shown in  that after the projection of the Bethe-Salpeter amplitude on the light front, the six components of the LFD deuteron WF are expressed through integrals over the eight components of the deuteron Bethe-Salpeter amplitude. Provided the NN interaction is the same, these approaches incorporate by different methods the same relativistic dynamics. The WF’s in LFD are the direct relativistic generalization of the non-relativistic ones in the sense that they are still the probability amplitudes. Therefore they can be used in the relativistic nuclear physics (e.g. ). Another advantage of LFD is that due to the explicit relativistic covariance of the method, the form factors do not depend on the system of reference and on the direction of the $`z`$-axis. One should emphasize that the LFD WF’s have been successfully used also in the hadron physics .

As mentioned in the Introduction, the LFD results give a good description of the data  for the tensor polarization $`t_{20}`$ and for the structure function $`A(Q^2)`$ (up to $`Q^2`$$`3`$ (GeV/c)2). Concerning the structure function $`B(Q^2)`$ its minimum occurs below the position indicated by the data. The latter is related  to the insufficient accuracy of the perturbative approach in calculating the components of the LFD wave functions. We would like to emphasize, however, that the LFD method describes very well the charge and quadrupole form factors which are expressed by the structure functions A(Q2), B(Q2) and the tensor polarization $`t_{20}`$.

In this Section we firstly calculate the nucleon momentum distribution in deuteron $`n(q)`$ within the LFD using the WF $`(\mathrm{\Psi }(\stackrel{}{q},\stackrel{}{n}))`$ which is normalized generalizing the non-relativistic normalization condition :

$$\frac{1}{3}\frac{m}{(2\pi )^3}\mathrm{\Psi }^2(\stackrel{}{q},\stackrel{}{n})\frac{d^3\stackrel{}{q}}{\epsilon (\stackrel{}{q})}=\frac{m}{(2\pi )^3}F(\stackrel{}{q}^2,\stackrel{}{n}\stackrel{}{q})\frac{d^3\stackrel{}{q}}{\epsilon (\stackrel{}{q})}=1,$$
(1)
where $`F`$ is expressed by the six scalar functions $`f_{16}`$ depending on the scalar variables $`\stackrel{}{q}^2`$ and $`\stackrel{}{n}\stackrel{}{q}`$:

$$F=f_1^2+f_2^2+f_3^2+(3+z^2)f_4^2+(1z^2)(f_5^2+f_6^2)+(3z^21)f_2f_3+4zf_4(f_2+f_3).$$
(2)
In (1) $`\epsilon (\stackrel{}{q})=\sqrt{\stackrel{}{q}^2+m^2}`$, where $`m`$ is the nucleon mass. As shown in , in the non-relativistic limit ($`q<<m`$) the functions $`f_{36}`$ become negligible, $`f_{1,2}`$ do not depend on $`z`$ and turn into $`S`$\- and $`D`$-waves ($`f_1u_S`$, $`f_2u_D`$) and the WF $`\mathrm{\Psi }(\stackrel{}{q},\stackrel{}{n})`$ becomes the usual non-relativistic wave function. One of the most important properties of the functions $`f_{16}`$ found in  is that for $`q`$ 2$`÷`$2.5 fm-1 the component $`f_5`$ (being related mainly to $`\pi `$-exchange) exceeds sufficiently the $`S`$\- and $`D`$-waves. This fact is very important in the calculations of $`n(q)`$ in deuteron as it will be shown below.

In  a Legendre expansion of $`f_{16}(q,z)`$ relative to $`z`$ has been used:

$$f_i(q,z)=\underset{l}{}(2l+1)f_i^l(q)P_l(z)$$
(3)
($`l`$=0,2 for $`i`$=1,2,3,5 and $`l`$=1,3 for $`i`$=4,6) with given values of the coefficients $`f_i^l(q)`$.

We calculate in our work the angle-averaged nucleon momentum distribution in the deuteron defined as:

$$n(q)=C(q)\overline{F(q)},$$
(4)
where

$$C(q)=\frac{m}{(2\pi )^3\epsilon (q)}$$
(5)
and $`\overline{F(q)}`$ is the angle-averaged function $`F(q,z=\mathrm{cos}\theta )`$:

$$\overline{F(q)}=\frac{1}{4\pi }F(q,\mathrm{cos}\theta )𝑑\mathrm{\Omega }.$$
(6)
In accordance with Eq. (1) the normalization of $`n(q)`$ is given by:

$$n(\stackrel{}{q})d^3\stackrel{}{q}=1.$$
(7)
The LFD calculations have shown that, as expected, the most important contributions to the total $`n(q)`$ give terms related to the $`f_1`$, $`f_2`$ and $`f_5`$ functions

$$n(q)n_1(q)+n_2(q)+n_5(q),$$
(8)
where

$$n_1(q)=C(q)\overline{f_1^2(q)},n_2(q)=C(q)\overline{f_2^2(q)},n_5(q)=C(q)\overline{(1z^2)f_5^2(q)}.$$
(9)
Here we would like to note that the use of the functions $`f_1`$, $`f_2`$ and $`f_5`$ averaged over $`z=cos(\widehat{\stackrel{}{q},\stackrel{}{n}})`$ can be justified by their smooth (almost constant) behaviour as functions of $`z`$ at different values of $`q`$ (see Figs. 10 and 11 of ). The contributions of $`n_1`$, $`n_2`$, $`n_{12}=n_1+n_2`$ and $`n_5`$ are compared in Fig. 1. It can be seen that, while the functions $`f_1`$ and $`f_2`$ give a good description of the YSD of $`n(q)`$ for $`q<2`$ fm-1 (like the $`S`$\- and $`D`$-wave functions in the non-relativistic case), it is impossible to explain the high-momentum components of $`n(q)`$ at $`q>2`$ fm-1 without the contribution of the function $`f_5`$. We note that the deviation of the total $`n(q)`$ from the sum $`n_{12}=n_1+n_2`$ starts at $`q`$ around 1.8 fm-1. All this shows the important role of NN interactions which incorporate exchange of relativistic mesons in the case of the deuteron. We emphasize the consistency of the LFD method obtained in the simultaneous description of both deuteron momentum distribution (made in this work) and of charge form factors . Here we would like to discuss this point. Considering the LFD results on the nucleon momentum distribution in the deuteron we do not imply that the LFD method and the contribution of the $`f_5`$ function are the only way to understand the high-momentum components of $`n(q)`$. In the following of this Section we will pay attention to the results of other approaches considering, however, this question in close connection to the possibility of simultaneous description of both momentum distribution and form factors. It is known, for instance, that the deuteron wave function related to the Argonne v18 potential for the NN interaction  gives a realistic description of $`n(q)`$ in the deuteron . However, as can be seen from Figs. 2 and 4 of Ref. , the recent non-relativistic impulse approximation calculations using the Argonne v18 potential and even those including meson exchange currents and relativistic corrections do not give a good agreement with the new data of  on the tensor polarization $`t_{20}`$ and the charge form factors $`G_C`$ and $`G_Q`$.

We calculate in the present work also the angle-averaged nucleon momentum distribution in the deuteron using $`S`$\- and $`D`$-wave functions $`\mathrm{\Psi }_S(q)`$ and $`\mathrm{\Psi }_D(q)`$ corresponding to various NN potentials, such as the charge-dependent Bonn potential , the Argonne v18 , the Nijmegen - I ,- II and - Reid 93  and Paris 1980  in the expression:

$$n(q)=\frac{1}{4\pi }[\mathrm{\Psi }_S^2(q)+\mathrm{\Psi }_D^2(q)]n_S(q)+n_D(q)$$
(10)
with

$$n(q)d^3\stackrel{}{q}=1.$$
(11)
In Fig. 2 the result for $`n(q)`$ using the charge-dependent Bonn potential  is given and compared with the YSD. As can be seen, the $`D`$-component of $`n(q)`$ is important but even its inclusion does not give a very good agreement with the data for $`q2`$ fm-1. In the next Fig. 3 we present the LFD result for $`n(q)`$ compared with the calculations using the WF’s corresponding to Nijmegen-I,-II, -Reid 93, Argonne v18 and Paris 1980 NN potentials and with the $`y`$-scaling data. We would like to note that: i) the results of the calculations using the NN potentials, such as Nijmegen-II, -Reid 93, Argonne v18 and Paris 1980 (shown in Fig. 3) are in better agreement with the YSD than those using the charge-dependent Bonn potential (Fig. 2). This might be related to the fact that these potentials describe NN phase shifts up to larger energies (e.g. the Nijmegen-II potential gives reasonable $`pp`$ phase shifts up to 1.2 GeV, while the charge-dependent one-boson exchange Bonn potential fits the phase-shift data below 350 MeV); ii) It can be seen from Fig. 3 that there are small differences between the curves corresponding to different NN potentials for $`q`$ 3 fm-1 (which give a good description of the YSD and almost coincide with the LFD result) and larger ones for $`q>3`$ fm-1. Large differences take place, however, between all of them and the LFD result for $`q>3`$ fm-1. An important conclusion of this fact is that data at higher momentum transfers are needed (and probably a more refined analysis of the currents used for the coupling of the electron with the nucleons) in order to distinguish the properties of the covariant LFD method from the potential approaches. The momentum region of $`q`$ just larger than 3 fm-1 will be decisive in this sense. The main advantage of the LFD method is, however, as mentioned already, that it describes simultaneously the deuteron charge form factors as well, which is not the case with e.g. the Argonne v18 NN potential method (see Fig. 4 in ).

## III Nucleon Momentum Distributions in 4He, 12C and 56Fe

As mentioned in the Introduction, realistic many-body calculations (e.g. ) have shown the existence of high-momentum components in the nucleon momentum distributions for studied nuclei at $`q2`$ fm-1, due to presence of short-range and tensor correlations. This conclusion agrees well with the $`y`$-scaling data . Secondly, it has been shown that all nuclear momentum distributions $`n_A(q)`$ for $`(A,Z)`$-nuclei ($`A>2`$) at these high momenta are simply rescaled versions of the nucleon momentum distribution in the deuteron $`n(q)`$ :

$$n_A(q)\alpha _An(q),$$
(12)
where $`\alpha _A`$ is a constant.

Since the magnitude of the high-momentum tail of $`n_A(q)`$ is proportional to the mass number, this effect is associated with the nuclear interior rather than with the nuclear surface. In the present paper we develop a model to calculate $`n_A(q)`$ for nuclei with $`A>2`$ from that of the deuteron. The method has similarities to that suggested in  on the basis of 4He. However, in the present work we incorporate the analysis performed in Sect. II for the deuteron as the simplest two-nucleon bound system to the consideration of $`n_A(q)`$ for heavier nuclei, namely because the high-momentum components are related to the specific features of the two-nucleon interaction at short distance (around the repulsive core). The latter determine the corresponding behaviour of the nucleons in the central part of the nucleus, where the density is higher. There the two nucleons could be closer to one another and the short-range correlations could be operative. In this sense, we consider the present method as a next step and improvement on that of .

In this Section we suggest a practical method to calculate the nucleon momentum distribution for the general case of $`(A,Z)`$-nuclei ($`A>2`$), while as we mentioned already, the correlation methods consider only particular nuclei (e.g. 3,4He, 12C, 16O, 40Ca and nuclear matter) and due to their complexity the applications to other nuclei is often a very difficult task.

We use the natural orbital representation (NOR)  of the one-body density matrix. In our case we consider its diagonal elements in the momentum space, which represent the nucleon momentum distribution $`n_A(q)`$. The NOR allows us to use the transparency of the single-particle picture even being beyond the mean-field approximation and to account for the correlation effects in a proper way. Using the common theoretical ground of NOR, the method makes it possible to combine the mean-field predictions for $`n_A(q)`$ (which are expected to be realistic at $`q<2`$ fm-1) with that part of the momentum distribution which includes correlation and relativistic effects.

In the natural orbital representation the momentum distribution in a nucleus with $`A`$ nucleons can be written in the form:

$$n_A(q)=N_A[n_A^h(q)+n_A^p(q)],$$
(13)
where $`N_A`$ is the normalization constant and

$$n_A^h(q)=\frac{1}{A}\underset{nlj}{\overset{F.L.}{}}2(2j+1)\lambda _{nlj}C(q)R_{nlj}(q)^2$$
(14)
$$n_A^p(q)=\frac{1}{A}\underset{F.L.}{\overset{\mathrm{}}{}}2(2j+1)\lambda _{nlj}C(q)R_{nlj}(q)^2$$
(15)
are the hole- and particle-state contributions, respectively. In (14) and (15) F.L. denotes the Fermi level, $`C(q)`$ is given by (5), $`\lambda _{nlj}`$ are the natural occupation numbers ($`0\lambda _{nlj}1`$) for a state with quantum numbers $`nlj`$ and $`R_{nlj}(q)`$ are the natural orbitals. We call hole-state natural orbitals those natural orbitals for which the numbers $`\lambda _{nlj}`$ are significantly larger than the remaining ones, called particle-state natural orbitals . As shown in , the high-momentum components of the total $`n_A(q)`$ caused by short-range correlations are almost completely determined by the contributions of the particle-state natural orbitals. This fact, together with the approximate equality of the high-momentum tails of $`n(q)/A`$ for all nuclei allows us to make the main assumption of the method that the particle-state contributions to $`n_A(q)`$ are almost equal for all nuclei. The deuteron $`n(q)`$ from (8) can be written in a form similar to (13):

$$n(q)=N_d[n_d^h(q)+n_d^p(q)]$$
(16)
with

$$n_d^h(q)=\frac{1}{2}[n_1(q)+n_2(q)],$$
(17)
$$n_d^p(q)=\frac{1}{2}n_5(q),$$
(18)
and $`N_d=2`$. The assumed equality of the particle-state contributions for all nuclei implies the equality of the particle-state contribution $`n_A^p(q)`$ for $`(A,Z)`$-nuclei from (15) with that to the deuteron momentum distribution $`n_d^p(q)`$ from (18)

$$n_A^p(q)n_d^p(q).$$
(19)
The right-hand side of (19) can be taken from different sources. This could be the YSD or results of calculations within given theoretical correlation models. Having the results from Sect. II we take $`n_d^p(q)`$ to be related to $`n_5(q)`$ (18) because within the LFD namely this function is responsible for the high-momentum tail of $`n(q)`$ for deuteron:

$$n_A^p(q)\frac{1}{2}n_5(q).$$
(20)
Then the resulting momentum distribution for $`(A,Z)`$-nucleus normalized to unity has the form:

$$n_A(q)=N_A[\underset{nlj}{\overset{F.L.}{}}2(2j+1)\lambda _{nlj}C(q)R_{nlj}(q)^2+\frac{A}{2}n_5(q)],$$
(21)
where

$$N_A=\{4\pi 𝑑qq^2[\underset{nlj}{\overset{F.L.}{}}2(2j+1)\lambda _{nlj}C(q)R_{nlj}(q)^2+\frac{A}{2}n_5(q)]\}^1.$$
(22)
It is known from  that the hole-state natural orbitals are almost unaffected by the short-range correlations and, therefore, the functions $`R_{nlj}(q)`$ in (21) and (22) can be replaced by single-particle wave functions from the mean-field approximation (shell-model or Hartree-Fock WF). In our calculations we use Woods-Saxon s.p. wave functions for protons and neutrons. The hole-state occupation numbers $`\lambda _{nlj}`$ are close to unity and we set them equal to unity with good approximation. These properties lead to a similarity of our model to that suggested for calculations of the spectral function .

The calculations of $`n_A(q)`$ for 4He,12C and 56Fe are given in Figs. 4-6 and are compared with the $`y`$-scaling data from . We note that, as expected from (12) the value of $`\alpha _A`$ must be almost constant. In our case $`\alpha _A=N_AA/2`$ and its numerical values are 8.335, 8.352 and 8.372 for 4He, 12C and 56Fe, respectively.

The dashed line in Figs. 4-6 gives the calculated normalized hole-state contribution (the mean-field component) to $`n_A(q)`$ only, i.e. without the contribution of $`n_5(q)`$ in Eq. (21). The comparison between the total $`n_A(q)`$, its hole-state contribution and the $`y`$-scaling data clearly shows the role of $`n_5(q)`$ contribution. It can be seen that the difference between the results with $`n_5(q)`$ and without it starts for all three nuclei at $`q1.5`$ fm-1 and increase rapidly at $`q>1.5`$ fm-1. This shows the role of the high-momentum components of the deuteron momentum distribution $`n_5`$ which can be considered as an ”ingredient” of the high-momentum components of the momentum distribution of $`(A,Z)`$-nucleus. In our opinion, the related function $`f_5`$ of the LFD method incorporates the main part of the short-range features of the NN interactions which determine the correlation effects seen in $`n_A(q)`$ when calculated within the non-relativistic correlation methods.

The good agreement of the results on $`n_A(q)`$ with the data for the three nuclei achieved, including also the region of the high momenta, proves the assumption that the nucleon momentum distribution in nuclei can be extracted on the basis of the realistically described properties of the deuteron as the smallest bound system of two nucleons.

## IV Conclusions

The results of the present work can be summarized as follows:

(i) The light-front dynamics method is applied to calculate the nucleon momentum distribution in the deuteron. A very good description of the $`y`$-scaling data is obtained. It is shown that, while the functions $`f_1`$ and $`f_2`$ describe well the YSD of $`n(q)`$ for $`q<2`$ fm-1, it is impossible to explain the high-momentum components of $`n(q)`$ at $`q>2`$ fm-1 without the contribution of the function $`f_5`$. Thus, this study gives a positive answer to the question about the possibilities of the LFD method to describe simultaneously both deuteron momentum distribution and charge form factors (the good agreement with the latter as well as with the data for the tensor polarization $`t_{20}`$ and for the structure function $`A(Q^2)`$ up to $`Q^23`$ (GeV/c)2 was shown in ).

(ii) The nucleon momentum distribution in deuteron was calculated also in this work using the $`S`$\- and $`D`$-wave functions corresponding to various NN potentials (charge-dependent Bonn, Paris 1980, Nijmegen-I and -II, -Reid 93, Argonne v18). We show that the results using these potentials explain very well (with the exception of the charge-dependent Bonn potential) all the available data for $`n(q)`$ up to $`q3`$ fm-1 exactly like the LFD method. However, as shown in , the recent non-relativistic impulse approximation calculations (e.g. with Argonne v18 potential) and even those including meson exchange currents and relativistic corrections do not describe well the new data  on the tensor polarization $`t_{20}`$ and the charge form factors $`G_C`$ and $`G_Q`$. This is the main difference between the results of the LFD method and of the calculations using WF’s corresponding to NN potentials for the important deutron characteristics considered. Secondly, we should mention also that for $`q>3`$ fm-1 the LFD results for $`n(q)`$ deviate strongly from those of the calculations using NN potentials. Thus data at higher $`q`$ are needed in order to distinguish the abilities of the LFD from those of the potential methods. Here we would like to mention also other relativistic methods such as the method in  (using a three-dimensional reduction of the Bethe- Salpeter equation) and the method in  (based on the Gross equation ) which could give a simultaneous description of both deuteron form factors and momentum distribution.

(iii) A correlation method for calculating the nucleon momentum distribution in nuclei with $`A>2`$ $`n_A(q)`$ is proposed. The method uses the natural orbital representation of the one-body density matrix and combines the mean-field component of $`n_A(q)`$ with the correlated high-momentum components taken from the deuteron momentum distribution. The method is based on the well-known results of realistic many-body calculations which have shown that, due to the short-range and tensor NN correlations, all nuclear momentum distributions $`n_A(q)`$ at high momenta ($`q>2`$ fm-1) are rescaled versions of the nucleon momentum distribution in the deuteron $`n(q)`$. In the present work we incorporate in the method the LFD results for $`n(q)`$. The role of the wave function $`f_5`$ is clearly shown by the comparison of the calculations of $`n_A(q)`$ using the hole-state (mean-field) contribution without and with an inclusion of the $`n_5`$($`f_5`$) contribution. In our opinion, the LFD function $`f_5`$ incorporates the main part of the short-range features of the NN interactions which can be seen in calculations of $`n_A(q)`$ within non-relativistic correlation methods. The suggested method gives an easy practical way for calculations of $`n_A(q)`$ for any nucleus because only hole-state wave functions and occupation probabilities enter its main relationships (Eqs.(21) and (22)). The method gives a good description of the $`y`$-scaling data for $`n_A(q)`$ in 4He, 12C and 56Fe including the high-momentum region, showing a clear difference from the mean-field predictions for $`q>1.5`$ fm-1.

## V Acknowledgments

One of the authors (A.N.A.) thanks Prof. V.A. Karmanov and Prof. J. Carbonell for the helpful discussions. He is also grateful to the Royal Society in London and the Bulgarian Academy of Sciences for support during his visit to the University of Oxford. The work was partly supported by the Bulgarian National Science Foundation (Contracts No.$`\mathrm{\Phi }809`$ and $`\mathrm{\Phi }905`$).