domain: phys spi

text:
Janina Löffler
email: janina.loffler@hevs.ch

Mohamed Belhaj

Nenad Bundaleski

Juan J Diaz Leon

Jonathan Thomet

Samira Frey

Christophe Ballif

Nicolas Wyrsch

A comprehensive analysis of electron emission from a-Si:H/Al2O3 at low energies

Keywords: secondary electron emission, electron emission energy, a-Si:H, Al 2 O 3, high secondary electron emissive layers, microchannel

Recently developed microchannel plates based on amorphous silicon offer potential advantages with respect to glass based ones. In this context, secondary electron emission at very low energies below 100 eV has been studied for relevant materials for these novel devices. The aim of this work was to quantify the low energy electron emission -secondary emission and elastic scattering -from amorphous silicon and alumina and the dependence of the emission energy distribution on the primary electron energy, which was previously unknown. Secondary emission and energy distribution were both modelled and measured using equipment particularly designed for this energy range. The effects of roughness, angle of incidence and surface composition were analysed. We show crossover energies as well as the angular dependence of electron emission from amorphous silicon and alumina, with a maximum experimental emission yield value of 2 and 2.8, respectively, at an incident angle of 75°. A parameterization for the energy dependence of the emission energy spectrum at low energies was derived. This extensive analysis is fundamental for a comprehensive understanding of the performance of amorphous silicon-based microchannel plate detectors. It provides a complete model for secondary electron emission for a detailed description of the detector operation. The present results thus set the basis for a simulation framework, which is an essential element to increase the performance of these detectors and enable further developments.

Introduction

Secondary electron emission (SEE) is a fundamental physical phenomenon behind the working principle of electron multipliers. The latter represent a wide class of devices that provide sensing of low intensity radiation, down to the detection of single particles. An energetic electron impinging a dynode of a discrete electron multiplier produces several secondary electrons, which are then accelerated towards the next dynode. In the following collisions the electrons are further multiplied, resulting in an electron avalanche, eventually providing a current pulse per each particle impinging the detector.

We see from the previous example that the total electron yield (TEY), defined as the average number of electrons emitted per incident primary electron, is of extreme importance in electron multipliers. More precisely, dependence of the TEY on primary electron energy plays an essential role in their design and operation. The second important feature of SEE, from the perspective of electron multipliers, is the energy distribution of emitted electrons, since the initial electron energy determines both the electron trajectory and its impinging energy onto the next dynode.

Glass based microchannel plates (MCPs) are electron multipliers with high spatial resolution of the order of tens to hundreds of µm that were developed in the 1970s (Ladislas Wiza (1979)). However, both the collection efficiency, i.e. the probability that a particle will enter a channel, and temporal resolution of MCPs are limited as channel geometries and the vertical conductivity are largely fixed. State-of-the-art clean room processes have recently enabled fabrication of a novel electron multiplier with high resolution based on hydrogenated amorphous silicon (a-Si:H) (Franco et al. (2014)) with the potential to overcome the inherent limitations of MCPs. Amorphous silicon based microchannel plates (AMCPs) can be integrated monolithically and the vertical conductivity can be adjusted through doping. AMCP channel spacing and diameters can be customized down to narrow diameters below 3 µm (Löffler et al. (2018)). During the fabrication of AMCPs, channels with walls having different surface morphologies can be realized, which could additionally increase TEY. Another strategy in improving MCP performances is to deposit high TEY materials on the inner walls using atomic layer deposition (ALD), which is known to make conformal coatings even over nanostructured surfaces (Elam et al. (2006)). Although layers with higher secondary electron emissivity can be deposited by ALD, Al 2 O 3 layers are particularly suitable for this purpose due to the well-established and broadly available technology.

From the gain of AMCPs, the average number of collisions inside the channel has been calculated by Löffler et al. (2018). With a number of collisions of at least 6 and a bias voltage of 500 V, secondary electrons that are emitted with few eV gain an average of about 80 eV. An average energy of 40 eV has been found by [START_REF] Frey | IEEE Nuclear Science Symposium and Medical Imaging Conference[END_REF] for electrons that exit the channels. Knowing that electron paths between two collisions inside the channel become even shorter for high aspect ratio AMCPs with smaller channel diameters, electrons inside those channels hardly gain kinetic energies above 100 eV. Therefore, designing AMCPs with the most suitable architecture requires detailed understanding of SEE from this particular material in the energy range which is rarely investigated, i.e. below 100 eV. Thus, the present analysis is key to understanding the full potential and the limits of AMCPs, and to be able to assess their performance as a function of device geometries and surface morphology of the channels.

Low energy electron-matter interaction is a subject of great uncertainty with measurements of the TEY at 0 eV being controversial, as discussed by Andronov et al. (2013). Most of the commonly used physical models are inadequate for electron energies below 50 eV (Kieft and Bosch (2008)), implying that TEY modelling is also unreliable. While more accurate theories for low energy electron scattering are still to be explored, we have to rely on TEY measurements in this energy range. The latter, being a challenge on their own, have been performed mainly for materials that are used in space technologies by [START_REF] Balcon | [END_REF] including for metals such as Cu by Schulte et al. (2020), Ag, Au and crystalline silicon (c-Si) by Pierron et al. (2017). Further experimental difficulties originate when performing such measurements on nonconductive surfaces prone to surface charging.

In this paper, we extended the low energy SEE studies to a-Si:H and ALD-deposited Al 2 O 3 . A setup presented by Belhaj et al. (2010), highly suitable for TEY measurements from non-conductive surfaces, was used to determine the previously unknown TEY of a-Si:H and Al 2 O 3 surfaces for energies below 100 eV. In addition, the energy spectra of emitted electrons from the same surfaces for primary electron energies below 200 eV were measured using a setup particularly designed for such purposes by Villemant et al. (2017). The importance of this investigation is in the ability to follow the change of the overall spectra (e.g. the relative contribution of elastically backscattered electrons) with the primary electron energy. Furthermore, the study can be extended to other high TEY coatings, such as MgO and composite structures.

In parallel to the experimental investigations, we used Monte Carlo (MC) simulators and analytical models. Furthermore we developed an original approach to consider the influence of surface morphology on TEY and applied it to the materials of interest. This allowed us to evaluate their capabilities in describing SEE phenomena at low primary energies in the case of a-Si:H and Al 2 O 3 and thus their utility potential for future AMCP simulators.

Materials and methods

Sample preparation

The thin films of a-Si:H and Al 2 O 3 used in this study were deposited on HF cleaned polished c-Si wafer substrates by means of plasma enhanced chemical vapor deposition (PECVD) and ALD, respectively. Two different procedures were used for the deposition of a-Si:H, yielding distinct surface morphologies i.e. roughness. Some of the a-Si:H films were additionally coated by Al 2 O 3 .

The film thickness of both resistive materials studied here was kept as thin as possible in order to minimize charging problems, yet thick enough to measure electron emission from the film itself. Since the minimal film thickness is given by the secondary electron escape depth of few nm (Fitting et al. (2001)), all a-Si:H layers in this study were aimed to be about 15-20 nm.

The PECVD setup and deposition recipe of a-Si:H were described in detail by Franco et al. (2014). Rough surfaces were produced using the aforementioned recipe applied for the deposition of AMCP layers. Briefly said, the former deposition recipe is optimized for low intrinsic stress for the deposition of thick layers, resulting in a corrugated surface, when deposited as a thin film of few nm. Additionally, a procedure optimized for surface flatness was used to deposit a second set of samples. By comparing the SEE properties of the two types of films, we were able to estimate the influence of surface morphology on the overall secondary electron emission behaviour.

Al 2 O 3 was deposited in an Oxford FlexAl ALD setup using a mixture of Trimethylaluminum (Al 2 (CH 3 ) 6 ) and water. The thickness of the deposited films was about 5 nm on both the corrugated a-Si:H surface and the substrate surface (used as the reference sample).

Thickness of thin films was measured with an ellipsometer and interpreted using the dispersion model developed by Cauchy (1830).

Surface characterization

The two main characteristics of surfaces, i.e. their topography and their chemical composition, govern the process of secondary electron emission. We analysed the topography using an atomic force microscope (AFM) Bruker VT-102-2. The data analysis was performed using Gwyddion v.2.58. Surface chemical composition analysis was performed by means of X-ray photoelectron spectroscopy (XPS). The XPS setup consists of a monochromatic X-ray gun with Al anode, providing a photon energy of 1486.7 eV, and a hemispherical electron energy spectrometer. The spectra were acquired in the fixed analyzer transmission mode with the pass energy of 200 eV. The line intensities were determined after removal of the background of Shirley type. The composition was determined using appropriate sensitivity factors, under the assumption of in-depth uniform samples. Hydrogen concentration was not measured. In general, it is about 5-20% in a-Si:H, with a hydrogen rich surface (Blayo and Drévillon (1992)). Hydrogen on the surface can be present as adsorbed H 2 O and is usually removed by annealing in AMCPs, as described by Franco (2014).

SEE measurements

SEE measurements were performed on the ALCHIMIE and DESSE facilities located at Onera, both providing a base pressure of 10 -10 mbar. The TEY was measured at ALCHIMIE/DEESSE using a charge collection method with a pulsed electron beam, which has been proven to avoid charging effects and consequently enables an accurate measurement of the electron emission yield of insulators at low energies. The method and setup employed in this work were described in detail by Belhaj et al. (2010). The schematic of the setup is presented in figure 1a. The incident electron angle was varied from normal incidence (0°) up to 75°. The incident electron energies ranged from 10 to 1000 eV. The primary electron current was measured by a Faraday cup that was introduced between the electron gun and the sample holder. For each incident energy, 10 pulses of 6 ms and a current in the nA range were delivered. The TEY was calculated as the ratio between the sample charge and the incident charge.

The energy distribution of the emitted electrons was measured for electron irradiation with different incident energies using the hemispherical electron energy analyzer Omicron EA 125 of DEESSE, suitable for an energy range from few eV to 2 keV. The energy spectrometer was working in the fixed analyser transmission regime, with the electron pass energy of 200 eV. The schematic of the setup inside the vacuum chamber is shown in figure 1b. The samples were tilted by 45°with respect to the incident beam with monochromatic energy E 0 and the emission energy analyzer was facing the sample plane. The sample was kept at a bias of φ = -11.5 V. For an incident electron energy E 0 and flux Γ 0 , electrons with energies E e in the range from 0 eV to E 0 were emitted from the sample.

The energy E a recorded at the analyzer was

E a = E e -eφ + (ϕ sample -ϕ spectrometer ) (1)
with the elemental charge e and the respective work functions of the sample and the spectrometer, ϕ sample and ϕ spectrometer . An extensive correction was elaborated for low energy emission spectra measurements by Villemant et al. (2017), as the electron energy analyzer is typically calibrated for energies higher than 300 eV. Additionally, geometrical constraints due to the limited field of view of the analyzer were taken into account. In summary, the emission spectra were obtained from the recorded spectra Γ a (E a ) at energies E a with an energy step ∆E = 0.2 eV as follows:

N (E e ) = 1 M i=1 Γa(E i ) H(E i ,φ) • Γ a (E a ) H(E e , φ)∆E , (2)
where N (E e ) is the normalized emission probability at energy E e . The transfer function

H(E e , φ) = H T (E e , φ) + H a (E a ) is composed of: H a (E a ) = A • exp (α • E a ) for 0 ≤ E a < 124 eV 1 for E a ≥ 124 eV, (3)
with A = 0.0192 and α = 0.0319 for the correction of the analyzer calibration and

H T (E e , φ) =    1 for 0 ≤ E e ≤ 3.25 eV arctan ( d 2•h ) arctan Ee |φ| for E e > 3.25 eV, (4)
with d = 0.03 m and h = 0.03 m for the geometrical correction.

SEE modelling

We used existing MC tools, Nebula by van Kessel and Hagen (2020) and Casino by Demers et al. (2011), to model SEE. Nebula was used with the modelling setup described by Löffler et al. (2019). The number of primary electrons in all Nebula results shown throughout this paper was 10 6 . The electron-material interaction simulator 3D Casino ver. 3.3 has been used with the standard physical models implemented, using the scattering cross section of the ELSEPA code by Salvat et al. (2005) and the energy loss model by [START_REF] Joy | [END_REF]. For each incident energy between 10 eV and 300 eV we simulated 3000 electrons for incident angles of 0°and 80°.

In addition, we used a novel simplified approach to model the energy dependence of total electron yield TEY(E) from corrugated surfaces. The latter represents an extension of the earlier work by Bundaleski et al. (2015), in which we predicted TEY(E) from a smooth surface for an arbitrary incident angle from the experimentally measured TEY at normal incidence. As an input to the model TEY(E) from the flat in depth uniform sample measured at normal incidence is used. Then, an AFM image of the actual emitting surface, representing the corrugation function z(x,y), is used to calculate the TEY(E) of that surface. The surface is considered as a mosaic of K flat areas, each having an arbitrary inclination with respect to the surface normal of a perfectly smooth sample. We used the AFM image to determine the local inclination at a point 'i' of the sample surface and the corresponding TEY i (E) according to the model by Bundaleski et al. (2015). The trajectory of each of the N electrons leaving the surface was monitored. These electrons were ejected from the point 'i', with random outgoing angles following a cosine distribution. If the number of electrons that managed to leave the surface without colliding with the surface protrusions is M i , the TEY(E) of the corrugated surface is estimated as

TEY(E) = K i=1 M i N • TEY i (E). ( 5
)
The major drawback of this approach in modelling TEY from realistic surfaces is the assumption that the sample is indeed in-depth uniform. Besides, we neglect emission of the second generation of electrons, emitted due to impingement of the emitted electrons hitting the surface protrusions.

Results and discussion

Thin film characterization

We considered in this work four different samples, summarized in Table 1 together with their surface morphology parameters:

• Sample 1 is the corrugated a-Si:H layer (rough);

• Sample 2 is the smooth a:Si-H layer (flat);

• Sample 3 is the corrugated a:Si-H layer coated by Al 2 O 3 (rough AlOx);

• Sample 4 is the smooth Al 2 O 3 layer (flat AlOx).

c-Si < 1 2 flat PECVD a-Si:H 12 < 2 polished c-Si < 1 3 rough AlOx ALD Al 2 O 3 5 10* PECVD a-Si:H 10 4 flat AlOx ALD Al 2 O 3 5 < 1* polished c-Si < 1
*Roughness given by the underlying substrate or layer

The typical surface topography of the samples fabricated with different roughness is shown in figure 2 and3. Both figures show a 2×2 µm 2 AFM scan of a sample with 5×5 cm 2 overall surface area. AFM measurements were done with a Bruker AFM Dimension Edge with the AFM probe TESPA-V2, which has a nominal tip radius of 7 nm. The topography was measured in tapping mode and image analysis was done with Gwyddion (Nečas and Klapetek (2012)). All surface parameters were calculated according to ISO standards.

For the corrugated a-Si:H layer in figure 2, we measured an average roughness of S a = 1.8 nm, with a mesokurtic shaped surface with normal peaks and a kurtosis of S ku = 3.01. The latter indicates height distribution close to a perfect Gaussian. Considering that ALD coatings are highly conformal, figure 2 and 3 illustrate the typical surface texture of the corrugated and smooth AlOx samples as well.

Since the surface composition controls the secondary emission, we used XPS to analyse the chemical composition of the samples as deposited. The samples were inevitably exposed to air to perform the post-fabrication analysis. Samples would need to be plasma cleaned inside the measurement setup in order to determine the secondary yield of the pure material a-Si:H. As AMCPs are not plasma cleaned before use, we measured the samples as deposited in order to determine a most comparable data set. An overview of their surface composition is given in Figure 4. Carbon is a standard surface contaminant present in all samples, typically in the form of adventitious carbon (saturated aliphatic hydrocarbons). We see that its relative amount is increased for both AlOx samples, being the result of the ALD coating process. Small quantities of fluorine in the same samples most probably have the same origin. No XPS analysis of the smooth sample was performed as it was expected to behave as the corrugated sample. Hydrogen content is not detected by XPS, both a-Si:H layers (corrugated and smooth) are typically composed of about 10% hydrogen.

Emission yield

The TEY was measured for all samples at incident electron angles ranging from 0°( normal incidence) to 75°and for incident electron energies starting at 10 eV up to 1000 eV. Furthermore, the number of secondary electrons per incident electron, the secondary emission yield (SEY) was modeled for the AMCP surface.

Hydrogenated amorphous silicon

The TEY of a-Si:H, measured for the corrugated and the smooth sample, is shown in figure 5. It was compared to experimental values of crystalline Si by Pierron et al. (2017); Bronshtein and Fraiman (1969) and from MC calculations by Pierron et al. (2017) and Kuhr and Fitting (1999). We also show the TEY of Si calculated with Nebula.

The measurements and simulations agree generally well for crystalline silicon, except for the low energy region where the more recent measurements by Pierron et al. (2017) show a higher yield than that predicted by the simulations. Of particular interest are the electron energies at which TEY = 1 (known as crossover energy), since the condition TEY > 1 is necessary to secure electron multiplication in MCPs. Compared to crystalline silicon, a-Si:H shows a higher emission yield and practically two times lower crossover energy of about 50 eV.

The incorporation of H into the atomic network of Si, yielding formation of a-Si:H, leads to H and OH termination. The overall effect is an increased electron emission yield of a-Si:H as compared to Si. This characteristic increase of the TEY due to H termination has also been observed for diamond by Shih et al. (1997). Indeed, electronic ab-initio calculations by Tao et al. (2015) have shown a decreased electron affinity of SiN through hydrogen termination, which effectively lowers the potential barrier for electron emission. The increase in emission yield observed here could thus be due to a lower potential barrier at the surface of a-Si:H. If this is not the case, the main reason for the increased TEY would be an increased energy gap, which forbids various energy loss processes and consequently extends the electron trajectories. The surface morphology of a-Si:H is also beneficial for secondary electron emission. With a maximum total emission yield over 1.6, the rough surface texture leads to a higher yield than the smooth surface. Additionally, we observe a further shift in the crossover energy of about 20 eV towards a lower energy for the corrugated a-Si:H sample compared to the smooth a-Si:H sample.

The TEY dependence on the incident angle is shown in figure 6. For higher incident angles we observe an increase towards higher yields. This is a well known phenomenon related with the deposition of the energy of primary electrons in a shallower region than at normal incidence. Therefore, the generated secondary electrons can be emitted with higher probability. This can be modeled for a number of pure conductive materials according to different approaches, such as the one presented by Bundaleski et al. (2015). It was demonstrated in the latter that the crossover energy shifts towards lower energies with increasing incident angle. In figure 6, we compare the crossover energy of the smooth layer at 75°with the measurements and literature values already shown in figure 5. The crossover energy for the smooth layer at 75°is about 20 eV lower than that at normal incidence. A decreased crossover energy is observed for the corrugated surface even at normal incidence. Indeed, this is due to local surface inclinations of corrugated surfaces, which cause an oblique incidence of primary electrons. This frequently results in amplification of the secondary electron emission. Due to the increase in elastic scattering for even higher incident angles, Bundaleski et al. (2015) predict a maximum yield at about 70°. Therefore, the TEY should decrease with further increasing angle.

In AMCPs an average emission yield of about 1.6 at an average incident energy of about 70 eV has been inferred from gain measurements of AMCP devices by Franco et al. (2014). Given the rough surface texture inside AMCP channels, this assumption is in good agreement with the observed trends. Electrons are accelerated with grazing trajectories in AMCPs and local electric field gradients could be the reason for an even further increase of their electron emission yield.

Aluminum Oxide

While the inherent TEY of a-Si:H layers can be enhanced by roughening its surface, highly emissive layers need to be used to further improve AMCP performances. Hence the use of materials that can be deposited via ALD, since it is the only method capable of uniformly coating the inside of a channel with high aspect ratio. We show the total emission yield of Al 2 O 3 in figure 7, where the measurements and MC calculations are compared to literature values.

Calculated TEY include results of MC simulations by Ivanov et al. (2018) at the respective incident angles of 0°and 80°. van 2017) for a 12.5 nm thick one. In figure 7, we show the experimental total emission yield of Al 2 O 3 at 0°incidence and at 75°incidence.

It can be observed that the different assumptions in the theoretical models lead to a great variation of the resulting TEY. The experimental emission yield data on the other hand agree well, especially in the low energy range. A higher maximum yield of the experimental result by van der Graaf et al. ( 2017) could be due to the layer thickness. Thicker layers of Al 2 O 3 are expected, theoretically, to exhibit a higher yield, when charging is avoided. As seen in the previous section, higher incident angles on flat layers also increase TEY. In reality the surface morphology, surface contamination and charging effects greatly affect experimental results. The TEY of the smooth of Al 2 O 3 layer at 75°is below the expected increase with incident angle. This and the flattened shape of the curve indicate that this particular measurement was subject to charging.

The increase of the TEY on a corrugated a-Si:H layer by using a highly emissive Al 2 O 3 coating is shown in figure 8. As the highly resistive a-Si:H layer leads to charging of the surface, the secondary yield of Al 2 O 3 cannot be fully determined when deposited on this rough surface. However, despite this charging effect and a total layer thickness of about 22 nm, we observe a higher secondary electron yield with the additional Al 2 O 3 coating compared to the yield of the 17 nm thick corrugated a-Si:H layer.

Typically, the ALD deposition temperature of Al 2 O 3 is between 200°C and 300° ( Guo et al. (2019)). We found that the deposition temperature has no significant effect on the TEY, as shown in figure 8. The smooth Al 2 O 3 layer, prepared at 200°C and a similar layer deposited at 300°C with otherwise identical parameters exhibit almost identical TEY values.

Total emission yield in AMCPs

How do the previous results allow one to deduce the TEY in AMCPs? As previously seen, the total emission yield largely depends on the surface morphology. Therefore, we considered the influence of the real surface structure in AMCP channels, shown in figure 9. In general, the surface inside the channel varies depending on the fabrication process and within one channel, along its depth. We also discuss the TEY for the limit of electrons with negligible energy approaching 0 eV.

Modelling TEY of the AMCP channel surface

Using the measured TEY(E) of smooth and corrugated a-Si:H layers shown in section 3.2.1, the expected emission yield of a corrugated AMCP channel surface has been calculated based on the model described in section 2.4. The aim of this calculation was to evaluate the possible effects of the AMCP channel surface morphology on the TEY. As a first step, the expected emission yields of the smooth and corrugated a-Si:H layers were calculated from their surface scans smooth sample at 0°was taken from the experiment to calibrate the TEY(E) model. Hence, the calculated TEY at 0°(flat calc) corresponds to the measured TEY (flat), as shown in figure 10 a. The calculated TEY for the corrugated surface (rough calc) shows a similar yield. In figure 10 b, we show the calculated TEY at 75°for the smooth (flat calc) and corrugated surface (rough calc) and the experimental TEY for the smooth sample (flat). The difference in surface roughness does not lead to a significant difference in TEY according to the calculations. The experimental yield however is significantly different from the calculated ones. Since the experimental TEY cannot be explained by the surface roughness, other effects, like hydro-carbonation of the surface or, charging effects might have caused a lower experimental yield than expected for the corrugated surface as well.

The surface corrugation function z(x,y) of the AMCP channel was reconstructed from an SEM image and is shown in figure 9. The expected TEY of this surface morphology was calculated (channel) and is shown in figure 11 alongside the calculated yield for the smooth surface (flat).

From the theoretical electron emission model, the surface morphology of the channel is expected to reduce the TEY, as electrons can be recaptured. In an operating AMCP on the other hand, the local electric field strength varies and could lead to re-emission of electrons and a lowered crossover energy. This means that, in general, surface morphology alone cannot be used to calculate secondary yield in AMCP channels, and while TEY measurements of a-Si:H layers can be used as a basis, the secondary yield of AMCPs must be deduced from AMCP gain measurements. Moreover, the maximum yield, and the yield curve at higher energies as just presented determine only a fraction of the electron dynamics, whereas the crossover energy plays a major role in determining the level of multiplication in the channel. The expected lowering of the crossover energy for corrugated surfaces could not be observed with the model.

Discussion

-Total emission yield at low energies Overall, we found in the previous section that the total emission yield decreases below 1 only for energies below 40 eV. Notably, the measured TEY does not approach 0 when the incident energy approaches 0 eV. In fact, elastic scattering is the predominant process at low energies, as will be confirmed in the following section. However, the TEY for incident electrons with energies close to 0 eV might be prone to the following measurement artefacts. Electrons may not reach the surface and therefore not interact with the surface at all. Also, electrons may be absorbed. The nature of the interaction and the resulting TEY depend as well on the local work function, which is not uniform for the samples presented here. Furthermore, the electron beams used to measure TEY consist of electrons with an energy spread of about 0.8 eV. In theory, it can be assumed that TEY(0 eV)=1 as discussed by Cazaux (2012), which has been shown experimentally for selected materials by Cimino et al. (2015). For AMCP modelling we assume hence a TEY(0 eV)=1, where all electrons are backscattered, when absorption and other effects are neglected.

With the previous results, we provided an overview of the magnitude and slope of the TEY at low energies as well as the crossover energies of a-Si:H and Al 2 O 3 and their dependencies on surface corrugation and incoming angle. Due to the influence of the surface chemistry and texture on the emission yield, the AMCP channel surface cannot be fully represented as a thin film and thus cannot be fully characterized by electron emission measurements. However, the present study is a valuable starting point for the parametrization of emission yield curves in order to describe the multiplication process, in the low energy regime, in electron multipliers.

Emission energy

The goal of this analysis is to describe the operation of an AMCP channel, where we are studying an electron 'cascade': the energy of each emitted electron has a direct consequence on its incident energy and angle, which then affects the number of electrons emitted and trajectories in the forthcoming steps. As presented above, the number of emitted electrons depends on the incident electron energy. At a fixed electric potential and AMCP channel geometry, and assuming a cosine function for the emission angle of electrons, the incident energy and angle of an electron i.e. the TEY and therefore the number of electrons in the cascade solely depends on the electron emission energy.

The emission energy distribution N(E e ) was measured for incident energies below 200 eV for the rough a-Si:H layer and the rough Al 2 O 3 layer. The emission energy distribution was determined with the electron energy analyzer shown in figure 1b for the effective incident electron energies of 12.5 eV-192.5 eV for a-Si:H and 3.5 eV-194.5 eV for Al 2 O 3 , corrected according to the equation 2 and according to their work function taken from Cazaux (2010). The spectra are shown in figure 12.

Emitted electrons can be categorized according to the physical process they undergo when interacting with the surface: secondary electrons (SE), inelastically backscattered electrons (IBS) and elastically backscattered electrons (EBS). Because SE and IBS can not be separated at low energies and for the purpose of comparing the measurements with literature values, we artificially divide the spectrum into three parts, although this does not have a physical meaning. EBS are characterized by an emission energy close to the incident energy (+/-2 eV). An emission energy both below 50 eV and 2 eV below the incident energy refers to a SE. The other electrons are considered as IBS. For primary electron energies below 50 eV we consider all electrons which are not EBS as secondary.

The TEY presented in section 3.2 is the sum of the secondary emission yield δ, the EBS yield η e and IBS yield η ie . Several empirical (Vaughan (1989)) and semi-empirical models (Bundaleski et al. (2015)) exist that describe δ. Therefore the resulting secondary yield curve can also be parameterized. The backscattering yield can also be described by existing models by Cazaux (2012); Cimino et al. (2004); Jablonski (2014). We used this categorization to obtain, with the help of the emission spectra measured at few energies shown above 12, probability functions for electron emission energies for all incident electron energies up to 500 eV. These functions could then be used for AMCP simulations. A comprehensive analysis of electron emission from a-Si:H/Al 2 O 3 at low energies 18 3.4.1. Relative backscattering probability At first, we used the energy spectra shown above to deduce the relative backscattering probability P . P is a useful measure to determine the relative proportion of high energy electrons versus low energy electrons that are emitted upon an incident electron with a particular energy. P is particularly useful in the event of elastic scattering, as all elastically backscattered electrons are emitted with the same as the incoming energy and their total emission yield is 1.

The elastic and inelastic backscattering yields are related to P in the following way: η e/ie = P EBS/IBS • T EY . We calculated P EBS as the integral of the elastic peak divided by the integral of the whole emission function. Figure 12 a) shows an inlet for 42.5 eV, where the area of the elastic peak is colored in blue, whereas the area of the secondary electron distribution function is colored in red. Similarly, P IBS is calculated from the integral of the remaining energies between 50 eV and the elastic peak.

As shown in the spectra in figure 12, we measured P EBS/IBS in four points for a-Si:H and in five points for Al 2 O 3 . In order to determine the backscattering probabilities for every possible incident energy in AMCPs, the experimental data points were fitted with an exponential decay function for the elastic contribution P EBS and with a linear function for the inelastic contribution P IBS . The probability of elastic backscattering at 0 eV was set to 1 P EBS (0) = 1 as in Cazaux (2012). The functions P EBS/IBS were arbitrary ones chosen to describe the experimental behavior:

P EBS (E in ) = (1 -a) exp (-b • E in ) + a for E in < 50 eV (6) P EBS+IBS (E in ) = P EBS (E in ) + c • E in + d for E in ≥ 50 eV. ( 7
)
The fit values are given in table 2.

Table 2: Fit values for the backscattering probability according to equations 6 and 7. We compared the incident electron energy dependence of the backscattering probability measured here for a-Si:H with existing literature values of Si, and with the results of MC calculations using Nebula and CASINO in figure 13. A similar comparison is shown for Al 2 O 3 in figure 14. An additional dependence of the backscattering probability on the incident angle can be observed for both materials, which was not taken into account in the parametrization. For both materials, variations in the relative backscattering are important in the low energy range below 50 eV. Consequently, the parametrization presented in this section can be used to separate the total emission yield into its respective contributions and calculate electron emission energies in AMCPs. 12) and fitted with equations 6 and 7, compared to the simulation results and values found in the literature.

avalanche process in an AMCP channel. In this section, we compare the low energy part of the a-Si:H emission energy spectrum to theoretical and experimentally measured secondary electron energy E s distributions ∂δ/∂E s (E in ) of a-Si:H and c-Si in figure 15. While secondary electrons are commonly defined to have an energy cutoff at 50 eV, most of the secondary electrons are emitted with energies below 10 eV.

The analytical functions include a physics based model by Henke et al. (1979) with

dδ/dE s (E s , χ, E g ) = C • E s (E s + χ -E g ) 2 • (E s + χ) 5/2 (8a)
for semiconductors and ) eq.( 8a), (Si), E g (Si)

dδ/dE s (E s , χ) = C • E s (E s + χ) 3 (8b)
eq.( 8a), (a-Si:H), E g (a-Si:H) eq.( 9) E av (192.5 eV) eq.( 9) E av (12.5 eV)

eq.( 10) E av (192.5 eV) eq.( 10) E av (12.5 eV)

eq.( 11) E max = 0.9 eV exp Si E in = 500 eV (Bronstein)

exp eq.( 9) E av (194.5 eV) eq.( 9) E av (13.5 eV)

eq.( 10) E av (194.5 eV) eq.( 10) E av (13.5 eV)

eq.( 11 for insulators, where χ is the electron affinity, E g the band gap and C a normalization constant. We use the values listed in table 3 for χ and E g . The model is based upon the assumption that the electrons are ejected from the valence band. Empirical secondary

) = C • E -1.5 av • E s exp -1.5 E s E av . (9)
with the average emission energy E av . We use E av (E in = 12.5 eV) = 2.5 eV with a normalization factor of C = 5 and E av (E in = 192.5 eV) = 8.2 eV with C = 22.9. Furthermore, we tested the fit function of Price and Fraser ( 2001)

dδ/dE s (E s , E av ) = 3.3 • (E s /E av ) 0.47 • exp (-2.21 • E s /E av ) (10)
and the fit function of Kruschwitz et al. ( 2011)

dδ/dE s (E s , E mp ) = exp -(log (E s /E mp )) 2 2 • 1.33 2 (11)
with the most probable emission energy E mp = 0.9 eV. None of the theoretical curves can represent the measured curves perfectly, as shown in figure 15. We found the best fit with equations 11 and 8a for incident energies above 25 eV. Although the fit function of the equation 10 may be appropriate for the energy spectrum at E in = 12.5 eV, it does not hold for higher incident energies.

The low energy part of the Al 2 O 3 spectra were compared to the theoretical secondary electron distribution curves of Al 2 O 3 in figure 16. Here, the best fit was obtained with the theoretical curve of equation 8b for incident energies above 25 eV. None of the functions could be used to represent the energy spectrum at E in = 13.5 eV.

For incident energies below 25 eV, the emission energy distributions of both materials depend on the incident energy and cannot be represented by any of the functions above. In this energy range we use an energy dependent fit function for incident energies below 25 eV for a-Si:H and Al 2 O 3 using an energy dependence proposed for copper and stainless steel by [START_REF] Furman | [END_REF]:

∂δ/∂E s (E in ) = Γ((P 1 + P 2 )/2) • (P 1 /P 2 ) P 1 /P 2 • E P 1 /P 2 -1 s Γ(P 1 /2) • Γ(P 2 /2) • (1 + P 1 /P 2 • E s ) (P 1 +P 2 )/2 (12)
where Γ(x) is the Gamma function of x and x is an integer. In contrast to the original interpretation of P 1 and P 2 as probabilities of emitting 1 or 2 electrons, respectively, we consider these magnitudes as functions of the incident energies. The best fit was obtained for P 1 = 39 • exp(-E in /68) + 111 and P 2 = 2. The fit function is shown in figure 17, where the resulting probability density function (PDF) is compared with the measured secondary electron emission spectra of a-Si:H and Al 2 O 3 . In the same figure we also present corresponding cumulative distribution functions, obtained by integration of the PDFs, which are necessary for performing MC simulations of AMCPs.

From the energy distribution analysis in this chapter we derived a parameterization for the energy dependence of elastic backscattering. While secondary emission energy models from literature could represent the emission energy spectrum for high incident energies, we found that at low incident energies, the electron energy spectrum depends on the incident energy. Finally, we showed a suitable parameterization for the energy emission spectrum of both a-Si:H and Al 2 O 3 using equation 12 at low incident energies.

Conclusion

In this paper, we provide a complete analysis of electron emission characteristics, which is necessary for reliable AMCP modelling. The electron emission energy distribution and the emission yield of a-Si:H and Al 2 O 3 with two different surface textures have been characterized at low incident electron energies. A comparison to literature values confirmed the need for this detailed analysis, as there are significant variations between measurements, analytical expressions and Monte Carlo models. Furthermore, we defined parametrizations for the expected electron scattering process depending on the incident electron energy and determined analytical functions describing the energy dependence of the secondary electron energy distribution.

In general, amorphous silicon shows a higher emission yield than crystalline silicon which has a maximum yield below 1.2. The average emission yield of a-Si:H is about 1.6, deduced by Franco et al. (2014) from the operation of AMCPs. In the present study, a rough surface leads to a further increase of the electron emission yield of amorphous silicon up to a maximum of approximately 1.6 at normal and of 2 at 75 degree incidence. Curvature and trenches in a microchannel modify the actual surface texture in the channel and therefore influence the emission yield and the crossover energy such that it cannot be reconstructed based on thin film measurements. In addition, the local electric field strength in AMCP channels under operation can modify the secondary emission yield. Together with electron incidence at grazing angles, this could be the reason for the effective emission yield of 1.6 observed so far in AMCPs.

Since coatings with highly emissive ALD films are the key to improve the potentially high gain achievable with AMCPs, we explored the low energy electron emission properties of ALD-deposited Al 2 O 3 . We measured a maximum emission yield of up to 2.8. The TEY of Al 2 O 3 can be further increased by increasing the film thickness. Yield values up to 3.5 were measured and theoretic models by van der Graaf et al. ( 2017) predict even higher electron yields.

We determined the energy and angular dependencies of the emission energy distribution for a-Si:H and Al 2 O 3 . In accordance with literature, the probability for elastic backscattering is high at low incident energies. For incident energies below 25 eV, the secondary emission spectrum depends on the incident energy, which is not the case for higher incident energies. For electron multiplication processes at these low energies, such as in AMCPs, this emission energy variation needs to be considered.

With the present analysis, the determining factors of electron emission at low energies in a-Si:H were identified and extended for Al 2 O 3 . The angular dependence and energy dependence of the electron emission yield and energy distribution are therefore understood and can be used for more detailed modelling and thus further development of electron multipliers. This can be particularly relevant for developing electron multipliers that deviate from the cylindrical shape of microchannels.

This study could also be extended to other highly emissive layers that can be deposited via ALD, such as MgO, Ti doped MgO or multi layered structures as indicated by Tao et al. (2016). With a maximum yield of more than 10 for an 11 nm layer measured by [START_REF] Jokela | AIP Conference Proceedings[END_REF], MgO has a real potential to tremendously improve AMCP performance and could enable them to reach a multiplication factor well above 10 5 .

Figure 1 :

 1 Figure 1: (a) Experimental setup for total electron emission yield measurements. The sample was tilted at various angles for an electron beam incidence from 0°up to 75°, relative to the surface normal. (b) Electron energy analyzer for characterization of the emission energy distribution. The sample was mounted at 45°with respect to the incident beam. This analyzer was calibrated for low energies by Villemant et al. (2017).

Figure 3 :

 3 Figure 2: AFM scan of the rough a-Si:H sample and resulting surface parameters.

Figure 4 :

 4 Figure 4: XPS spectra and resulting surface composition (atomic percent) of the corrugated a-Si:H layer and the corrugated and flat AlOx layer. *in situ XPS measurement

  Figure 5: Total emission yield of a-Si:H measured for the smooth a:Si-H layer (flat) and the corrugated surface (rough) compared to literature values and simulations of crystalline Si.

Figure 6 :

 6 Figure 6: (a) Total emission yield of a smooth a-Si:H sample, measured at incident angles of 0°, 30°, 60°and 75°, starting at 10 eV. (b) A zoom shows the crossover energy of a-Si:H with different surface textures (corrugated and smooth samples) and angles.

  der Graaf et al. (2017) used a low energy MC extension of a GEANT4 code to obtain the curves presented here. Additionally to Al 2 O 3 total emission yield measurements and simulations compared to experimental and calculated literature values. results from Nebula, we show results from Casino. The total emission yield of Al 2 O 3 has been measured experimentally by Jokela et al. (2012) for a 4 nm thick layer and by van der Graaf et al. (

  The total emission yield of Al 2 O 3 (flat AlOx) fabricated at temperatures of 200°C and 300°C respectively, was measured at incident angles of 0°and 75°. The total emission yield of the corrugated a-Si:H layer is improved by an Al 2 O 3 coating (rough AlOx), even though limited by charging when compared to the total emission yield of the smooth Al 2 O 3 layer.

Figure 9 :

 9 Figure 9: (a) SEM image showing the surface morphology inside an AMCP channel. (b) Height profile z(x,y) and height distribution of the AMCP channel for the area marked in (a).

Figure 10 :

 10 Figure 10: (a) TEY of the corrugated a-Si:H surface shown in figure 2 was calculated (rough calc) using TEY of the flat a-Si:H layer (flat exp) as a calibration and compared to the measurement (rough exp) at 0°. (b) TEY for the rough and the flat layer were calculated at 75°and compared to the experimental TEY of the flat layer.

Figure 11

 11 Figure 11: (a) TEY of the AMCP channel surface shown in figure 9 was calculated (channel calc) at 0°. The TEY of a perfectly smooth a-Si:H layer (flat exp) at 0°served as a calibration. (b) The calculated TEY for the flat sample and for the channel surface are shown at 75°.

3 Figure 12 :

 312 Figure 12: (a) Corrected a-Si:H emission spectra for incident energies of 12.5 eV, 42.5 eV, 92.5 eV and 192.5 eV. The inset shows the emission spectrum for 42.5 eV. We can identify the elastic peak at 42.5 eV, secondary electrons are all electron emitted with lower energies. (b) Corrected Al 2 O 3 emission spectra for incident energies of 3.5 eV, 13.5 eV, 43.5 eV, 94 eV, 194.5 eV. The inset shows the emission spectrum for 3.5 eV, where the elastic peak overlaps with emission of electrons with less than the incident energy.

  :H 0.01 0.05 4•10 -5 0.05 Al 2 O 3 0.005 0.06 2•10 -5 0.03

Figure 15 :

 15 Figure 15: Measured secondary emission energy distribution compared to theoretical curves from literature for a-Si:H and Si.

Figure 16 :

 16 Figure 16: Secondary emission energy distribution compared to theoretical curves from literature for Al 2 O 3 .

Figure 17 :

 17 Figure 17: Probability density function (PDF) and cumulative distribution function (CDF) of the secondary electron emission energy for low primary electron energies E in and their fit functions.

Table 1 :

 1 Thickness d and surface roughness S a of all samples from Ellipsometry, including the surface roughness S a of the underlying substrate or layer.

			top layer		substrate/lower layer
	sample	material	d [nm] S a [nm] material	S a nm
	1 rough	PECVD a-Si:H 17	10	polished

  3.4.2. Secondary emission energy distribution As secondary electrons drive AMCP performance, quantifying their emission energies is crucial for an understanding of the

	Backscattering probability P EBS	0.01 0.1 1					BS+IBS Backscattering probability P	0.01 0.1 1		Si MC 0° (Kuhr) Si MC 0° (Pierron) Si MC CASINO 0°S i MC CASINO 40°S i MC CASINO 80°S i MC Nebula 0°S i MC Nebula 40°S i MC Nebula 80°S i exp 0° (Bronstein) a-Si:H exp 45°a -Si:H exp 45° fit eq. (6,7)
		0	10	20	30	40	50	50	250	500
				Incident energy (eV)			Incident energy (eV)
	Figure 13: Experimental a-Si:H backscattering probability at 45°, measured at four
	points (see figure 12) and fitted with equations 6 and 7, compared to MC simulations,
	simulated and experimental literature values.
	Backscattering probability P EBS	0.01 0.1 1					EBS+IBS Backscattering probability P	0.01 0.1 1		MC CASINO 0°M C CASINO 40°M C CASINO 80°M C Nebula 0°M C Nebula 40°M C Nebula 80°e xp 45°e xp 45° fit eq. (6,7)
		0	10	20	30	40	50	50	250	500
				Incident energy (eV)			Incident energy (eV)
	Figure 14: Experimental backscattering probability of Al 2 O 3 measured at five points
	(see figure

Table 3 :

 3 Electron affinity and band gap of Si, a-Si:H and Al 2 O 3 . A comprehensive analysis of electron emission from a-Si:H/Al 2 O 3 at low energies 21 dδ/dE s (E s , E av

		Si	a-Si:H Al 2 O 3
	Electron affinity χ (eV) 3.2 3.92	1.35
	Band gap E g (eV)	1.12 1.8	7.0/7.6 Filatova and Konashuk (2015)

electron distribution functions are typically used for MCP modelling by

Yakobson (1966)

;

Price and Fraser (2001)

;

Kruschwitz et al. (2011)

. All of these functions are phenomenological. The fit function by

Yakobson (1966)

is

AUTHOR SUBMITTED MANUSCRIPT -JPhysD-132310.R1

Acknowledgements

We are grateful to Philippe Wyss for deposition of the flat a-Si:H layers, Mohammad Beygi for AFM imaging and Pierre Mettraux for the XPS analysis of the flat Al 2 O 3 layer. Results presented here were obtained with the financial support of the Swiss National Science Foundation project 200021 162525/1, the Swiss National Science Foundation through the Sinergia program (Grant No. 177165, 2018-2022) and the Portuguese National Funding Agency for Science, Research and Technology in the framework of the project CERN/FIS-TEC/0039/2019 and the grants UIDB/00068/2020 and UIDP/00068/2020. The authors have no conflicts of interest to declare.


Krishna Reddy

Rafael Estevez

Sumit Basu

Revisiting the mesoscopic Termonia and Smith model for deformation of polymers

Mesoscopic models for polymers have the potential to link macromolecular properties with the mechanical behaviour without being too expensive computationally. An interesting, popular and rather simple model to this end was proposed by Termonia and Smith (1987 Macromolecules 20 835-8). In this model the macromolecular ensemble is viewed as a collection of two-dimensional self-avoiding random walks on a regular lattice whose lattice points represent entanglements. The load is borne by members representing van der Waals bonds as well as macromolecular strands between two entanglement points. Model polymers simulated via this model exhibited remarkable qualitative similarity with real polymers with respect to their molecular weight, entanglement spacing, strain rate and temperature dependence. In this work, we revisit this model and present a detailed reformulation within the framework of a finite deformation finite element scheme. The physical origins of each of the parameters in the model are investigated and inherent assumptions in the model which contribute to its success are critically probed.

Introduction

Suitably designed mesoscopic models for the deformation and fracture of solid polymers form an important bridge between atomistic and continuum models. Atomistic methods such as molecular dynamics (MD) (e.g. [START_REF] Brown | Molecular dynamics simulation of an amorphous polymer under tension: 1[END_REF][START_REF] Capaldi | Molecular response of a glassy polymer to active deformation[END_REF][START_REF] Negi | A molecular dynamics study of the strength and ductility of high T g polymers Modelling Simul[END_REF] can simulate deformation behaviour at very high strain rates on nanoscopic specimens. On the other hand, continuum methods (e.g. Boyce et al 1988,WuandvanderGiessen1993) can model mechanical behaviour at macroscopic scales but are phenomenological in nature and as a result cannot directly establish connections between microscopic properties (such as molecular weight of the macromolecule) and macroscopic mechanical properties.

A simple lattice model for simulating the deformation and fracture of solid polymers has been proposed by Termonia andSmith (1987, 1988) and Termonia et al (1988). Henceforth we will refer to this as the TS model. This model treats entanglements between polymer chains as unbreakable (though in another version of the model slip of macromolecular chains through entanglement points and consequent occasional disentanglement is allowed) and a chain is thought of as being a self-avoiding random walk (SARW) through a regular grid of entanglement points. The portions of a chain between two entanglement points, called 'strands', behave as rubber elastic springs governed by the inverse Langevin function [START_REF] Treloar | Physics of Rubber Elasticity[END_REF]. Additionally, the non-bonded van der Waals interaction between the chains is represented by discrete linear elastic members. In the original model, both the van der Waals elements and the 'strand' elements can break leading to a loss in the stress carrying capacity of the material. Breakage rules for the van der Waals members are motivated by kinetic theory of fracture and the members are broken through a Monte Carlo lottery. Thus the breakage rule depends on a critical stress though no length scale characteristic of the mesh is included (Curtin andScher 1990, Jagota and[START_REF] Jagota | Element breaking rules in computational models for brittle fracture Modelling Simul[END_REF]. Models similar in spirit to the TS model have been used to deal with time dependent damage evolution (Curtin andPamel 1997, Vujosevic and[START_REF] Vujosevic | Creep rupture of polymers-a statistical model[END_REF] and fracture of fibres [START_REF] Gonzalez | A numerical model to simulate the deformation and fracture of polyethylene fibres Modelling Simul[END_REF]. Recently, variants of the Termonia and Smith model incorporating slippage of chains through entanglements have been used by [START_REF] Terzis | Entanglement network of the polypropylene/polyamide interface: 3. Deformation to fracture[END_REF] to study the deformation and fracture of bi-polymeric interfaces.

The original TS papers concentrated on the behaviour of polyethylene at high temperatures (>100 • C) but the modelling methodology is essentially applicable to even amorphous polymers in the glassy state. This is due to the fact that the stress-strain response predicted by the model for semi-crystalline polyethylene has qualitative features that resemble the behaviour of a range of polymers under a range of conditions. Comparison of the predictions from the TS model with experimental stress-strain behaviours of a number of commercial polymers has been performed by [START_REF] Bicerano | Microstructural model for prediction of stress-strain curves of amorphous and semicrystalline elastomers[END_REF]. The results for amorphous polymers showed good quantitative agreement with atactic polystyrene and ethylene-propylene-diene terpolymers at a range of temperatures.

Though the TS model is simple in principle, a systematic evaluation of the significance of the various parameters involved in the model has not been performed. In this work we present a clearly formulated finite element based numerical model for a two-dimensional polymer represented by SARWs on a lattice of entanglement points. The method adopted in this paper to deform the macromolecular ensemble is significantly different from thatusedbyT ermonia and Smith (1987) and [START_REF] Bicerano | Microstructural model for prediction of stress-strain curves of amorphous and semicrystalline elastomers[END_REF]. The method is based on the realization that the macromolecular ensemble is basically a two-dimensional truss-like structure amenable to a non-linear finite element analysis, where the individual members (the strands and the van der Waals members) can undergo large deformations and rotations. In this regard, the numerical algorithm used is similar in spirit to the rod finite elements adopted by [START_REF] Gonzalez | A numerical model to simulate the deformation and fracture of polyethylene fibres Modelling Simul[END_REF].

The TS model is able to reproduce the uniaxial stress-strain behaviour of polymers qualitatively. We show that the reason for this success lies in a changeover in the load bearing behaviour of the network with strain. In other words, the energetics of the network shift from being van der Waals dominated at low strains to rubber elasticity dominated at higher strains. Thus, parametrization of the discrete van der Waals members vis-a-vis the rubber elastic behaviour of the chains is the key to fitting the simulated results to the experimental ones. To this end, we take the parameters presented in the original papers asab a s i sa n d explore, through a systematic parametric study, the influence of the key parameters on the overall stress-strain response. This paper is organized in the following way. The next section describes the sample generation procedure, the finite element solution method and techniques used for breaking the van der Waals bonds and the macromolecular chains. In section 3 the influence of various model parameters is discussed. Conclusions and discussions on the efficacy of the present method are presented in section 4.

Computational procedure

The computational procedure consists of a few basic steps. First, a regular lattice of entanglement points (also called nodes) is filled up with SARWs of a given length N such that two and only two SARWs pass through each node. These SARWs represent macromolecular chains. Superposed on this ensemble is a regular lattice of van der Waals members obeying linear elasticity. The second step is to perform a finite element analysis of this two-dimensional structure under volume preserving uniaxial stretching at a particular rate. The finite element assemblage is composed of non-linear (strand elements) and linear (van der Waals members) truss elements. During deformation, both van der Waals bonds and the macromolecular chains can break leading to loss of stiffness. The analysis is stopped when the sample is separated into two parts.

Sample preparation

Figure 1(a) is a schematic showing the situation prevailing in a real polymer. The heavy black circles represent entanglements between the coiled polymer chains and the dotted lines denote the van der Waals bonds acting between different parts of a chain as well as between neighbouring chains. In figure 1(b) is shown the idealized TS model [START_REF] Termonia | Kinetic model for tensile deformation of polymers: effect of molecular weight[END_REF]. The heavy solid lines indicate chain vectors between entanglements and the individual van der Waals bonds have been replaced by discrete members (dotted lines) joining each entanglement to its nearest neighbours. In the TS model, the chains start or end at an entanglement point since loose chain ends do not bear any load and are thus superfluous in the present context.

The chain vectors in an actual, undeformed specimen are randomly oriented in a threedimensional space, i.e. cos 2 ψ = 1 3 , where ψ is the vector's angle with the draw axis y.I n our simplified, regular two-dimensional representation (figure 1(b)), we have taken all strands to be oriented at ψ = 54.7 • with the x-axis. We should note here that the TS model also used a similar orientation though in two-dimensions this does not lead to a strand that has equal orientation to both the axes. A value of ψ = 45 • is probably more suitable in two-dimensions. However, we found that the effect of the strand orientations on the overall results to be presented later is small and we have adopted ψ = 54.7 • to be consistent with the TS model.

The simplified model is generated with the following considerations.

(i) Each chain behaves like a SARW and thus a chain does not cross itself.

(ii) One chain does not overlap with another, i.e. no two chains have the same two consecutive entanglement points. (iii) Chains do not cross each other. (iv) Periodic boundary conditions apply in the x direction.

Ideally, the chain generation process stops after all entanglement points have been visited by two and only two chains. However, this may not always be possible satisfying all the above conditions (i)-(iv). Several trials may be required before a suitable sample is generated.

Selection of a suitable generation procedure for a single chain in a dense system is important for achieving complete occupancy of the lattice. We are interested in rather short SARWs (of 5-10 steps each) and so we adopt a simple algorithm suggested by [START_REF] Sokal | Monte Carlo methods for the self-avoiding random walk Monte Carlo[END_REF] for generating a single chain. It should be noted that more sophisticated schemes for single chain generation have been used for generating dense fully occupied systems by [START_REF] Meirovitch | A new simulation of branched polymers[END_REF] and [START_REF] Mansfield | Monte Carlo studies of polymer chain dimensions in the melt[END_REF] among others. However, for the short length SARWs that we use, the algorithm due to [START_REF] Sokal | Monte Carlo methods for the self-avoiding random walk Monte Carlo[END_REF] turns out to be adequate in the sense that it yields an almost monodisperse ensemble with all entanglement points occupied by two and only two chains.

Polymers live in three-dimensional spaces. The two-dimensional lattice used here is at best an approximation. The two-dimensional polymer sheet analysed can be thought of as being one out of an infinite number of sheets arranged like a stack of cards. This stack makes up the actual three-dimensional sample. Entanglements between chains are also basically three-dimensional entities and in two-dimensions are represented as nodes at which two selfavoiding walks are connected. The distance between the sheets at the onset in the z direction is given by T 0 . When the system is stretched in the y direction, it contacts in the other directions and consequently the distance between the sheets is also reduced. Assumptions used in determining T 0 are discussed in section 2.5.

Figure 2 shows the chain length distributions in the whole sample after filling it up with multiple strands. Some chains need to be truncated at the top and bottom surfaces of the domain while a few chains (especially when most of the system is filled up) cannot be grown to their desired length because of overcrowding. Both these contribute to t h ef a c tt h a tt h e average chain length p avg falls below the desired length p. In figure 2, the horizontal axis represents the size M of the M × M grid used. For the number of strands in a chain p = 1, 5, 10 and 15, we are able to generate p avg = 1, 4.3, 7.3 and 8.2, respectively. Also there is very little effect of M on the value of p avg . The polydispersity is somewhat high for SARW sizes greater than 5 steps produced with our simple algorithm for generating dense systems. For the model macromolecules considered in this study, a 5 step SARW represents a chain with about 700 monomers3 . The mesoscopic model of the polymer is obtained by taking into consideration both van der Waals bonds and the strands. This is shown in figures 3(a) and (b) for a small 25 × 25 grid containing 100 chains with p = 5. Figure 3(a) shows a sample with only the macromolecular strands. Each strand connecting two entanglement points is assumed to follow inverse Langevin based rubber elasticity [START_REF] Treloar | Physics of Rubber Elasticity[END_REF]. The van der Waals force acting between neighbouring atoms is lumped into a network of linear, breakable springs joining entanglement points and forming a regular two-dimensional triangular lattice as shown in figure 3

Finite element formulation for the assemblage

As mentioned earlier, the generated sample consists of non-linear members formed by the strands and linear van der Waals members. Under uniaxial stretching, the members can undergo large rotations and deformations.

In the finite element model we consider the general case of a truss member whose stiffness is a sum total of the stiffnesses of an inverse Langevin and a linear spring. Within the assemblage, a particular member may have either or both of these two contributions active at a point of time. Here we discuss the corotational formulation for such a general truss element.

We start with the equation of virtual power in terms of first Piola Kirchhoff stress P in the absence of body forces and surface tractions. Thus,

V 0 P : δ Ḟ dV 0 = 0, ( 1
)
where F is the deformation gradient tensor and for two second order tensors A and B,

A : B = A ij B ji .
In terms of the Cauchy stress σ,u s i n gP = J F -1 σ, L as the velocity gradient tensor with D as its symmetric part and J = det(F), we can write for volume preserving deformations (i.e. J = 0)

Ṗ = J F -1 ( σ -Lσ).
Thus, the incremental form of equation ( 1), namely, V 0 Ṗ : δ Ḟ dV 0 = 0, in the deformed configuration becomes

V δL :[ σ -Lσ]dV = 0, (2)
where δL = δ ḞF -1 . Now consider a one-dimensional 2-node truss element with a corotational axis x having a current cross sectional area A and a current length l (see figure 4). For a volume preserving deformation of this bar, along x, the virtual power equation simplifies to

V δ D x x ( σ -D x x σx x )A d x = 0,
where ()denotes quantities in the x coordinate system.

For the 2-node bar, the velocities û and v are interpolated in terms of their nodal values as

û = N 1 û1 + N 2 û2 and v = N 1 v1 + N 2 v2 ,
where N 1 = 1 -x/l,andN 2 =x/l. Thus,

D x x = [-1010 ] ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ û1 v1 û2 v2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ = B û. (3)
This immediately leads to

V δ ûT BT [ σ x x -B û σ x x ]A d x = 0. ( 4
)
In general an element supports two types of forces. The total force f on an element is given by

f = fs (λ) + fvd (λ), (5)
where fs is the 'strand force' arising out of the non-linear spring character of the strand and fvd is the van der Waals force. The quantity λ is the instantaneous stretch. Therefore

f . = f ′ s λ + f ′ vd λ. ( 6
)
The force on any member in the x direction is given by

f =σ x x A, ( 7
)
where A is the current area of the member, whereby,

f . = σ x x A + σx x Ȧ. ( 8
)
Using volume conservation the area and length in the reference configuration A 0 and l 0 can be related to A and l as A 0 l 0 = Al. The fact that the volume of each 'strand' element is conserved follows from the picture of the element that we have in mind. The radius of each element in the undeformed state is equal to r g ⊥ , the component of the radius of gyration perpendicular to the Gaussian end to end vector R ee . When the element is stretched uniaxially to λR ee , r g ⊥ scales as 1/ √ λ (see, [START_REF] Heymans | Radius of gyration, maximum extensibility and intrinsic crazing in thermoplastics[END_REF]). This implies volume conservation and we get

f . = σ x x A 0 λ - σx x A 0 λ 2 λ, (9)
where λ = l/l 0 . Using equation ( 6) in the LHS of equation ( 9)weget

f ′ s λ + f ′ vd λ = σ x x A 0 λ - σx x A 0 λ 2 λ, ( 10
) leading to σ x x = λ A 0 f ′ s + f ′ vd + σx x A 0 λ 2 λ = λ A 0 C(λ) λ, ( 11
)
where

C(λ) = f ′ s + f ′ vd + σx x A 0 /λ 2 . It can be easily shown that λ = λ B û. (12)
Using equation ( 12)in( 11), then inserting in equation ( 4) and using the usual finite element procedures yield

K û = 0, K = Kmat + Kgeo , (13)
where K mat is the material stiffness matrix and K geo is the geometric stiffness matrix. Further,

Kmat = V BT λ A 0 C(λ)λ B dV. ( 14
)
When written explicitly this becomes

Kmat = λ 2 A 0 C(λ) ⎛ ⎜ ⎜ ⎝
10 -10 00 00

-101 0 00 00 ⎞ ⎟ ⎟ ⎠ . (15)
Moreover, using the definition of B, it is easy to show that the expression for the geometric stiffness matrix is

Kgeo =- σ x x A l ⎛ ⎜ ⎜ ⎝
10 -10 00 00

-101 0 00 00 ⎞ ⎟ ⎟ ⎠ . (16)
The stiffness matrices in the ( x, ŷ) system K mat and K geo can be easily transformed to the global (X, Y ) system through

K mat = T T
Kmat T and

K geo = T T Kgeo T, ( 17
)
where the transformation matrix T is given by

T = ⎛ ⎜ ⎜ ⎝ cos θ sin θ 00 -sin θ cos θ 00 00 c o s θ sin θ 00 -sin θ cos θ ⎞ ⎟ ⎟ ⎠ . (18)
Therefore, in the global (X, Y ) system we have in terms of the angle θ the bar which it makes with the X axis (see figure 4),

K mat = λ 2 A 0 C(λ) ⎛ ⎜ ⎜ ⎝ cos 2 θ cos θ sin θ -cos 2 θ -cos θ sin θ cos θ sin θ sin 2 θ -cos θ sin θ -sin 2 θ -cos 2 θ -cos θ sin θ cos 2 θ cos θ sin θ -cos θ sin θ -sin 2 θ cos θ sin θ sin 2 θ ⎞ ⎟ ⎟ ⎠ . ( 19
)
The geometric stiffness is invariant to rotation and thus

K geo =- f s + f vd l ⎛ ⎜ ⎜ ⎝ 10-10 010-1 -1 010 0 -10 1 ⎞ ⎟ ⎟ ⎠ . (20)
The total stiffness in the global (X, Y ) system is given by the sum of the material and geometric stiffnesses as

K = K mat + K geo . (21)
According to classical theory of rubber elasticity for non-Gaussian chains, force on a chain strand having a vector length r is given by (see [START_REF] Treloar | Physics of Rubber Elasticity[END_REF])

fs = kT l s L -1 r Nl s , ( 22
)
where N is the number of statistical chain segments (Kuhn segments) and l s is their length.

The inverse Langevin function L -1 for an arbitrary variable s can be approximated as (see [START_REF] Perrin | Analytic stress-strain relationship for isotropic model of rubber elasticity[END_REF])

L -1 (s) = s(3 -s 2 ) (1 -s 2 ) . ( 23
)
Thus equation ( 22) becomes

fs = kT l s λ √ N (3N -λ2 ) N -λ 2 (24) and f ′ s = kT l s 1 √ N 3(N -λ 2 ) 2 +2λ 2 (3N -λ 2 ) (N -λ 2 ) 2 . ( 25
)
Further, the van der Waals force is assumed to lead to a linear elastic stress-stretch relationship and

fvd = K vd A vd 1 - 1 λ , (26)
where K vd is the elastic modulus. The area of the van der Waals member in the reference configuration is A vd , which is equal to A when the van der Waals member overlaps a strand. A different area is used for the horizontal van der Waals members which do not overlap with strand elements. Issues regarding estimation of K vd and A vd are discussed later. Thus,

f ′ vd = K vd A vd λ 2 . ( 27
)
Therefore the non-linear finite element problem can now be solved in principle with velocities applied at specific nodes.

Boundary conditions

The geometry chosen for analysis is as shown in figure 5. We start with a unit cell of dimension B 0 × H 0 × T 0 that deforms in a volume preserving way to a cell of dimensions B × H × T at time t. The boundary conditions in terms of the velocity components u and v applied on the edges of the unit cell are

u(0,y,t) = 0a n d v(x, 0,t) = 0. ( 28
)
In addition, at any point of time t v(x, H, t) = v 0 and In the above, the velocity of the top face v 0 is kept constant. Furthermore, in order to perform the analysis at constant volume, at any point of time we need to have the principal stretches in the y and x directions related as

u(y,B,t) = u 0 (t). ( 29
)
λ x = 1 √ λ y .
This can be achieved if

u 0 (t) =- v 0 2λ 3 2 y B 0 H 0 in equation (29).
The overall strain E is defined as

E y (t) = v 0 t H . ( 30
)
The overall stress y is obtained from

y = 1 V ∂U ∂E y , ( 31
)
where U is the total strain energy stored in the unit cell.

The unit cell also has a thickness T 0 such that only a single polymeric layer is enclosed in the cell. The initial thickness T 0 is assumed to be equal to the distance between the entanglement points √ Nl s , which equals the Gaussian end to end distance for a chain of length N. Thus, plane stress conditions are assumed to prevail where the thickness of the unit cell scales as

T = T 0 λ -1/2 y .

Sources of irreversibility

As the system of entangled polymeric chains analysed here is deformed, two debilitating processes occur. They are the breakage of van der Waals bonds and scission oft h e macromolecular chains. As in [START_REF] Termonia | Kinetic model for tensile deformation of polymers: effect of molecular weight[END_REF] van der Waals bonds break according to kinetic theory of fracture at the rate

ν = τ exp - (U 0 -β σ) k B T . ( 32
)
Here τ is the thermal vibration frequency and U 0 and β are the activation energy and volume, respectively. Also k B is Boltzmann's constant and T is the temperature. The simulation of van der Waals bond breaking is performed with the help of a Monte Carlo lottery that continues with a pre-set time step of t till a desired strain level E y is attained. At every time step, the calculations proceed in the following manner.

1. Apply a displacement increment v 0 t on the top face at y = H and the corresponding volume preserving displacements at the right edge x = B. 2. Solve the finite element system for the displacement increments applied in step 1 to obtain new displacements of all nodes and forces in all elements. 3. Start the bond breaking process. Initialize ζ = 0. 4. Calculate a probability

P i = ν i ν max , (33)
where ν i and ν max are the rates of breaking of the ith and the most strained bond in the ensemble. 5. Within the analysis time step t, visit the intact van der Waals bonds at random. 6. Generate a uniformly distributed random number R i on a visit to bond i. 7. Break bond i if R i <P i , else go to the next bond. 8. After each visit of a bond ζ is incremented by 1/[ν max n vd (t)], where n vd (t) is the total number of intact bonds at time t. 9. Terminate the process, increment time by t and go back to step 1 when ζ t. Additionally, in view of equation ( 23) it is evident that as the value of s → 1 the force fs →∞. We therefore break a strand when the draw ratio of a strand approaches √ N. In [START_REF] Washiyama | Fracture mechanisms of polymer interfaces reinforced with block copolymers: transition from chain pullout to crazing[END_REF], the minimum force required to cause scission of a polymeric chain is estimated to be around 2 nN. In our case, for a strand with N = 14, at T = 300 K, this would imply a stretch in the strand of the order λ ≃ 0.99 √ N. Thus, the probability of a strand breaking at lower stretches is negligibly small and unlike the van der Waals bonds, the strand breakage rule is assumed to be deterministic: a strand breaks when it attains a stretch of √ N.

The inherent assumption in the TS model therefore is that plasticity in polymers stems from the breaking of van der Waals bonds and scission of polymer chains. The TS model considers chain slippage and consequent disentanglement as possible additional sources of irreversibility. However, at low temperatures in the glassy regime disentanglement seems less likely compared with the mechanisms considered here.

Parameters

Several parameters are involved in the TS model. The most important parameters pertain to equation (32). In this work, following the original paper [START_REF] Termonia | Kinetic model for tensile deformation of polymers: effect of molecular weight[END_REF], we define a 'base' material for which the parameters are the same as in the TS model, i.e. U 0 = 1.9 × 10 -19 J, τ = 10 12 s -1 and β = 1.75 × 10 -26 m 3 . I ne q u a t i o n( 22), the number of statistical segments per strand is N = 14 with each segment being l s = 10 Å in length.

The simulations were conducted at T = 473 K unless otherwise mentioned. However, the exact value of the temperature is not very important (except to compare with the TS results) for our purpose since experimental uniaxial deformation behaviour at any other lower temperature can be fitted by a proper choice of parameters.

The van der Waals members have stiffness K vd = 4 MPa. There are two kinds of discrete van der Waals members in the assemblage. There are van der Waals members overlapping with the strand elements and horizontal ones. The horizontal van der Waals members have a cross sectional area A vd = 2Nl 2 s cos θ representing in a lumped manner the integrated non-bonded interaction between the neighbouring chains.

The overlapping van der Waals members and the strands have the same cross sectional area A 0 = πr g 2 ⊥ l s , assuming that these members are cylindrical rods having radii equal to the component of the radius of gyration r g ⊥ perpendicular to the strand axis. Assuming that the strand is randomly coiled at the onset, r g ⊥ = √ 1/3 √ N/6l s under the Gaussian approximation (see [START_REF] Treloar | Physics of Rubber Elasticity[END_REF]). The van der Waals forces along the strands represent the lumped effect of the intrachain non-bonded interactions. We use ∼10000 nodes (entanglement points) in our simulations representing an initial box size of about 0.8 × 0.8 μm 2 .

Results and discussions

In this section we present and discuss salient results from uniaxial deformation of the 'base' material. In particular we focus on the physical significance of the parametersinvolvedinthe TS model. The objective is to try and understand the inherent assumptions in the model so that at a later stage parameters for it can be derived from lower order simulations such as MD.

Significance of various parameters in the TS model

In figures 6(a)-(c), we illustrate the behaviour of the 'base' material as it is strained. In figure 6(a), the evolution of the number of intact van der Waals bonds n vd (normalized with n vd0 , the number of intact bonds at E y = 0) with the overall strain E y is demonstrated. Clearly, most van der Waals bonds break within the early part of the deformation. Beyond E y = 2no more van der Waals bonds break and the load is borne primarily by the strands. The transfer of load from the van der Waals bonds to the polymer strands manifest as a point of inflection on the total energy U versus E y curve shown in figure 6(b). The total energy is calculated by summing the strain energies of all elements in the system as

U(E y ) = i λ i 1 f i s l i 0 dλ i + λ i 1 f i vd l i 0 dλ i . ( 34
)
In the above equation, f i s and f i vd denote the strand and van der Waals forces on the element i, while l i 0 is its initial length. The quantity λ i is the stretch in element i when the overall strain is E y . At low strains, the total energy is dominated by the quadratic van der Waals energy. With further straining, the chains begin to stretch and van der Waals bonds also break at the same time. Thus the energy of the strands governed by the inverse Langevin function gains importance. This is also illustrated by the evolution of the individual energy components in figure 6(b). For the parameters chosen for the 'base' material, beyond E y ≈ 2.3, the energy of the strands dominates the total energy.

Since equation ( 31) is used to calculate the stress, the point of inflection in the total energy manifests as the yield drop in the overall stress-strain curve shown in figure 6(c). Thus, in this model the yield drop arises due to a changeover in the energetics of deformation from a van der Waals energy dominated situation to a strand energy dominated situation with a large scale breakage of the van der Waals bonds. The qualitatively realistic stress-strain curve obtained from the model exhibits linear elastic behaviour and low strains, softening at yield and finally strain hardening behaviour.

An important point to be noted about figures 6(a) through (c) is the gradual transition in the energetics which leads to a rounded yield drop in the stress-strain curve. The slope change at the point of inflection in the total energy depends on the parameters chosen and especially on the rate at which van der Waals bonds break. A sharp unrealistic drop in stress would result if these bonds broke too quickly. The yield drop in solid polymers and its influence on mechanical behaviour has been studied extensively (see, e.g., [START_REF] Lai | A numerical study of crack-tip plasticity in glassy polymers[END_REF], [START_REF] Smit | Predictive modelling of the properties and toughness of polymeric materials: I. Why is polystyrene brittle and polycarbonate tough?[END_REF]). In the TS model, the yield drop comes about as a consequence of a switch in energetics. The severity of the yield drop will consequently depend on the parameters governing van der Waals bond breakage (mainly U 0 and β) vis-a-vis the mechanical properties of the strand (mainly N). It is worth mentioning here that MD studies have shown (see, [START_REF] Negi | A molecular dynamics study of the strength and ductility of high T g polymers Modelling Simul[END_REF]) that the van der Waals energy does flatten with strain after the yield strain. However, further studies are needed to establish whether the physical reason for the yield drop is indeed linked to the early breakage of the van der Waals bonds.

As already discussed, the softening is a consequence of the early breakage of the van der Waals bonds. Thus the pre-softening elastic response is almost entirely dependent on the stiffness of the van der Waals bonds. This is shown in figure 7(c) where the parameter K vd in equation ( 27) is varied from 4 to 10 MPa. Increasing K vd has the effect of reducing the yield strain while keeping the yield stress almost unchanged. The yield stress is influenced more strongly by the activation energy U 0 and the activation volume β as shown in figures 7(a) and (b). The activation energy U 0 has been varied from 1.8178 × 10 -19 to 1.9178 × 10 -19 J (i.e. from the original TS values to somewhat higher ones: note that the simulations are very sensitive to the value of U 0 ). Higher U 0 makes the breaking of van der Waals bonds more difficult leading to a higher yield stress. On the other hand, increasing β makes the breaking of the bonds possible at lower levels of σ and consequently leads to a lowering of the yield

Effect of macromolecular and loading parameters

One of the major attractions of a mesoscopic model is that it allows us to connect molecular properties such as molecular weight M w of polymer chains and the entanglement spacing φ to the macroscopic mechanical response. In this section we demonstrate how M w and φ influence the stress-strain response of the polymer. The parameters chosen are identical to that of the 'base' material described in the previous section.

A single SARW step roughly corresponds to a chain of M w = 1900 g mol -1 .W eh a v e varied the molecular weight through 1900, 5700, 9500 and 19000. Thus a chainofM w = 19000 corresponds to 10 SARW steps. The overall stress-strain response of all these curves is shown in figure 8. The response is very brittle at M w = 1900 and gradually becomes ductile as M w increases. This behaviour is quantitatively similar to that observed by [START_REF] Termonia | Kinetic model for tensile deformation of polymers: effect of molecular weight[END_REF]. The stress-strain behaviour observed in these simulations is also qualitatively similar to the effect of M w on the ductility observed in experiments [START_REF] Nunes | Influence of molecular weight and molecular weight distribution on mechanical properties of polymers[END_REF]. In the experiments, the ultimate elongation as well as the impact strength of high density polyethylene (HDPE) has been shown to increase rapidly with molecular weight till a critical molecular weight of about 5 × 10 5 was attained. In real polymers, the effect of M w comes about due to increased intermolecular forces and entanglement effects at high M w . These effects are mimicked by the TS model successfully.

To systematically examine the effect on the deformation behavior of the entanglement spacing, i.e. of the number of statistical chain segments between entanglements, Termonia and Smith (1988) introduced the spacing factor φ defined as where M e is the molecular weight between entanglements. The numerical scaling factor 1900 is somewhat arbitrary but pertains to the approximate molecular weight between entanglements in melts of linear polyethylene macromolecules.

φ = 1900 M e , (35)
Figure 9 shows a series of nominal stress-strain curves calculated for monodisperse polyethylene of M w = 19 000 at different values of the entanglement spacing factor φ = 0.1, 0.2, 0.4, 0.8 and 1. A low value of φ corresponds to a lightly entangled polymer with strands longer than in linear polyethylene melts, while higher values simulate highly entangled polymers with shorter strands. This figure reveals that the post-yield modulus of highly entangled polymers (large φ) is higher than lightly entangled ones. Lightly entangled polymers have an almost perfect plastic response after yielding. Moreover, they bear much lower levels of stress though elongation at failure is somewhat larger for the lightly entangled samples. This again is in keeping with the results obtained by Termonia and Smith (1988) (see, figure 2 of their paper).

Polymers are known to be rate and thermally sensitive and Termonia et al (1988)h a d studied the effects of deformation rate and temperature on the TS model. The essential feature of the rate sensitivity is an increase in the yield stress at high rates of loading. On the other hand, temperature reduces the yield stress and weakens the hardening process. The weakening of the hardening response is due to thermally activated loss of entanglements(seeRahaand Bowden (1972), [START_REF] Arruda | Effects of strain rate, temperature and thermomechanical coupling on the finite strain deformation of glassy polymers[END_REF]). According to [START_REF] Raha | Birefringence of plastically deformed poly(methyl methacrylate)[END_REF], the density of strands obeys

n s (T ) = B -Dexp - E a RT , (36)
where E a and D are constants. Moreover since the total number of monomers in the ensemble is constant, n s N is fixed. This leads to the fact that the number of Kuhn segments between entanglements at temperatures T 1 and T 2 follows the relation Following [START_REF] Basu | A thermo-mechanical study of mode: I. Small scale yielding crack-tip fields in glassy polymers[END_REF] we have taken D/B = 1.648 × 10 3 and E a /R = 2.8 × 10 3 K in this work.

N 2 N 1 = 1 -D B exp -E a RT 1 1 -D B exp -E a RT 2 . ( 37
Figure 10 shows the effect of the strain rate on the stress-strain response of the 'base' material. The overall strain rates (i.e. rate of change in E y )usedare Ėy = 1, 5 and 10 min -1 . As expected, the material hardens with strain rate and, in particular, the yield stress increases with the increase in rate. In the TS model, the rate sensitivity arises from equation (32) which governs the breakage of the van der Waals bonds. This is illustrated in figure 11 where the fraction of the intact van der Waals bonds n s with respect to the initial number n s0 is plotted at various rates. At higher rates, owing to the short times involved in the deformation, larger number of van der Waals bonds survive large deformations allowing the material to carry higher stresses. Also note that the hardening response is not dependent on the deformation rate. In the TS model, the rate insensitivity of the large strain behaviour is due to thefactthatmostof the load at large strains is borne by strands which are governed by the rate-insensitive inverse Langevin function. This feature of the model is not entirely consistent with experimental results where the stress-strain curves are seen to shift upwards almost parallel to each other with the deformation rate (see, e.g., G 'Sell and Gopez (1981)).

Figure 12 shows the effect of temperature on the stress-strain response. At every temperature the number of Kuhn segments is adjusted according to equation (37). The large strain behaviour of the polymer at high temperatures is similar to that for low φ as in figure 8. This is not surprising since high temperatures lead to a larger number of Kuhn segments between entanglements.

Another interesting effect of high temperature is the change in the yield drop as the temperature changes. It should be noted that continuum models of polymers (e.g. Boyce et al 1988,WuandvanderGiessen1993) do not capture the change in yield drop with temperature. However, change in the yield drop with temperature is seen in experimental studies as in the work of G 'Sell and Gopez (1981). In this model, both the rate of breakage of the van der Waals bonds and the force stretch relationship of the strands depend on T . Thus the energetics of the van der Waals bond stretching as well as chain deformation are altered with temperature. For example, at 325 K a higher difference between the van der Waals energy and the energy of the strands at the yield point leads to a larger yield drop.

In summary, the above discussion reinforces the fact that the TS model is a viable mesoscale model for polymer deformation as it has intrinsic capability of capturing many of the essential features of the stress-strain response of polymers. These include the yield drop, the post-yield hardening response and the thermal sensitivity. However, in the present form rate sensitivity of these materials is not properly represented by the model.

The stress-strain responses shown in figures 8 and 9 match closely with the original results of Termonia andSmith (1987, 1988). This implies that the model formulated in this paper works as well as the original model proposed by [START_REF] Termonia | Kinetic model for tensile deformation of polymers: effect of molecular weight[END_REF]. However, the present approach systematically positions the method within the framework of the finite element method, the generality of which is firmly established. The results pertaining to the strain rate and thermal sensitivity of polymers are also realistic.

Conclusions

In this work, the Termonia andSmith (1987, 1988) model has been analysed critically. The primary objective is to ascertain the applicability of this mesoscopic model to predicting the mechanical response of polymers. The implementation of the model has been clearly explained, in our opinion, to an extent not done in the previous literature. The model has been implemented within the framework of the large deformation based finite element method wherein the two-dimensional lattice is viewed as a non-linear truss assemblage.

We have highlighted an important physical basis of the model, namely, that it implicitly assumes a switch in the energetics of deformation from being governed by the non-bonded van der Waals forces at low strains to the intrachain forces at higher strains. The parameters in the model basically tend to influence the competition between these two energies. The yield drop in the stress-strain response is a significant consequence of the switch in energetics. The post-yield hardening response of the polymer is on the other hand a consequence of the inverse Langevin nature of the strand response.

We have also looked at the significance of each parameter in the TS model and have shown that particular salient features in the stress-strain response of glassy polymers such as the elastic modulus, yield stress, yield drop and hardening are primarily governed by particular parameters in the model. A qualitative description of the stress-strain behaviour at low strain rates is thus given by the TS model.

Figure 1 .

 1 Figure 1. Schematic showing the two-dimensional TS idealization in (b)o fa na c t u a ld e n s e macromolecular system shown in (a).

Figure 2 .Figure 3 .

 23 Figure 2. Average chain lengths p avg of dense samples obtained using M ×M meshes. The desired chain lengths are given by p.

  (b). The final model consists of the triangular lattice (shown in figure 3(b)) superposed on the ensemble of SARWs representing the macromolecules in figure 3(a).

Figure 4 .

 4 Figure 4. A typical rod element used in this work. The symbols used are defined in the text.

Figure 5 .

 5 Figure 5. Schematic showing the problem analysed in this work and the applied boundary conditions. The quantity v0 is the constant velocity applied at the top face while u0 is the appropriate velocity applied on the right face in order to keep the volume of the unit cell constant.

Figure 6 .

 6 Figure 6. Variation of (a) the number of intact van der Waals bonds n vd as a fraction of the initial number of bonds n vd0 ,(b) van der Waals, strand and total energies and (c) overall stress y with imposed strain E y .

Figure 7 .

 7 Figure 7. Effects of the values of (a) U 0 ,(b) β and (c) K vd on the overall stress ( y ) and strain (E y ) response.

Figure 8 .

 8 Figure 8. Effect of molecular weight M w on the overall stress ( y ) versus the strain (E y ) response.

Figure 9 .

 9 Figure 9. Effect of entanglement spacing φ on the overall stress ( y ) versus the strain (E y ) response of the base material.

Figure 10 .

 10 Figure10. Overall stress ( y ) versus strain (E y ) responses at deformation rates of Ėy of 1, 5 and 10 min -1 for the base material.

Figure 11 .Figure 12 .

 1112 Figure11. Evolution of the fraction of the intact bonds (n vd /n vd0 ) with strain (E y ) at the same deformation rates as in figure10.

Assuming, as in the original TS paper, that each step comprises 14 Kuhn segments or strands per step and that each strand is made up of about 10 monomers.

Acknowledgments

One of the authors (SB) acknowledges the funding received for this work from the Department of Science and Technology, Government of India, through Project Number SR/S2/ CMP-35/2003.