text:
# UV-Manipulation of Order and Macroscopic Shape in Nematic Elastomers

## I Introduction

The unique behaviour of liquid crystal elastomers derives from the intimate relationship between the elastic nature of the polymer network, and the ordering of their mesogenic, liquid crystalline moieties. Many interesting and unusual properties of these materials have been identified, and are summarised in recent review articles fink89 ; ober93 ; brand98 ; emt99 ; WT96 .

The property that is most relevant to this paper is the response of a liquid crystal elastomers macroscopic shape to temperature. A permanently, uniaxially aligned monodomain nematic liquid crystal elastomer network kupfer94  will exhibit a spontaneous contraction along its director axis when heated towards its nematic-isotropic phase transition temperature assfalg99 ; wermter01 . This phenomenon is due to the coupling between the average polymer chain anisotropy, and the nematic order parameter Q(T) WT96 . The large variation in magnitude of resulting spontaneous uniaxial deformation depends strongly on the extent of this coupling between the elastic and the liquid-crystal degrees of freedom in the material static01 ; thermal01 .

Accordingly, a mechanical action can be produced in response to a control signal which defines the underlying order parameter of the elastomer. Substantial macroscopic movement and shape changes are brought about by subtle variations in microscopic ordering. There are many methods of achieving mechanical motion in response to a control signal - or actuation - for many different applications fleck97 , ranging from precise piezoelectric positioning devices to artificial muscles degen97 ; ratna01 . Liquid crystal actuators predominantly rely on thermal wermter01 ; thermal01 ; ratna01  and, to a much lesser extent, electromagnetic meyer97  control signals. Recently, however, a less invasive method of manipulating the underlying nematic order parameter in an elastomer network has been reported for the first time pereira01 .

When exposed to ultraviolet (UV) radiation, liquid crystal molecules containing photo-isomerisable groups, such as azobenzene, experience a reduction in their nematic order and can therefore be induced to undergo a phase transformation from the nematic to the isotropic state ikeda87 ; ikeda99 ; legge92 . This is due to the UV-induced cis-trans isomerisation of the azo- (N=N) bond, whereby the photo-sensitive mesogenic molecules change from a rod-like shape, to a strongly kinked shape upon irradiation at a resonant wavelength of 365nm pereira01 ; ikeda99 ; janossy98 . The rod-like shape serves to stabilise the liquid crystal phase, whereas the kinked shape has the opposite effect, acting as an impurity and destabilising the nematic phase, reducing its order parameter ikeda99 , Fig. 1. The dynamics of this process, and the reverse trans$``$cis reaction occurring on heating or irradiation at 465nm, are well studied in ordinary low-molar weight liquid crystals ikeda99 ; kim95 .

With this method of non-invasive control over the degree of nematic ordering, there is now a new mechanism for actuation in liquid crystal elastomers. Such mechanical photo-induced actuation in nematic elastomers was first reported by Finkelmann et al. pereira01 . The degree of wide-spread, far-reaching interest in this new photo-mechanical effect is highlighted by recent popular articles in the New Scientist and The Financial Times ft , which appear even before any in-depth research had been conducted, and suggest exotic applications such as surgical micro-actuators.

In this article, a range of such photo-sensitive nematic liquid crystal elastomers is studied. Before investigating their varied response to ultraviolet light, their thermal responses are first characterised in section III.1. The differing UV responses for the different sample compositions are presented in section III.2. Section IV then focuses on the detailed quantitative analysis of the UV-mechanical response, comparing the results with theoretical expressions derived in section IV.2.

## II Experimental

All starting materials and resultant aligned, monodomain nematic liquid crystal elastomers were prepared in the Cavendish Laboratory. The procedure for forming the side-chain polysiloxanes by the hydrosilation of the terminal vinyl group in the mesogenic moiety with the Si-H bond of the polysiloxane chain, as well as the two-step crosslinking technique with a uniaxial stress applied after the first stage of crosslinking to produce and freeze the monodomain nematic alignment, has been developed over the years by Finkelmann et al. kupfer94 ; wermter01 ; greve01 . A number of minor modifications were made to the procedure, which are described in thermal01 . The azobenzene compounds were all synthesised according to standard literature techniques boden83 ; vogel ; moreuv . A summary of the different mesogenic moieties used is given in Fig. 2.

From the structure of the molecular moieties one can deduce that all the side-group compounds are mesogenic rod-like units. PBB and BMAzo have a four carbon atom spacer, which induces a parallel orientation of the rods to the siloxane backbone resulting in a prolate chain anisotropy in the nematic phase, with its principal radii of gyration $`R_{}>R_{}`$ . The third rod-like side group, UBAzo, has a flexible spacer of eleven carbon atoms and is thus only expected to be weakly coupled to the siloxane backbone. There are three distinct types of di-vinyl crosslinking agents: the flexible, non-mesogenic DiUB, the rod-like DiUAzo and the main-chain nematic polymer DiMC, which contains a total of approximately 75 rod-like monomer units between its reacting vinyl groups. The polymer backbone onto which all these compounds were grafted was polyhydromethylsiloxane (PHMS), containing approximately 60 SiH groups per chain, obtained from ACROS Chemicals.

The polymer networks were all crosslinked via the same hydrosilation reaction, in the presence of a commercial platinum catalyst COD, obtained from Wacker Chemie. The crosslinking density was 10 mol% (molar percent, relative to the total number of SiH units) in all resultant elastomers.

Phase sequences were established on a Perkin Elmer Pyris 1 differential scanning calorimeter, which correlated with the critical temperatures obtained by the thermal expansion measurements (section III.1). In all resulting elastomeric materials a broad nematic liquid crystalline phase below a certain critical temperature was observed, with optical microscopy and mechanical testing confirming that this was indeed the nematic phase. Some of the materials also possessed a second ordered phase at low temperatures, exhibiting high rigidity and little response to temperature or UV-irradiation, just above the glass transition. This additional phase was not a focus of the present work - all measurements were performed in the soft rubbery nematic phase of the elastomer samples.

Elastomer mechanical expansion measurements were obtained by freely suspending the samples in an active, glass-fronted Gallenkamp oven and measuring their variation in length with a travelling microscope. Cooling and heating rates were chosen such that thermal hysteresis effects were negligible (around 0.5oC/min, except for the sample containing DiMC, which required a much slower rate).

Samples (of approximate dimensions $`30\times 5\times 0.4`$mm) were uniformly irradiated with a Merck VL-4.L UV lamp, providing a narrow band at 365nm with an output of approximately 4W, at a distance of approximately 15cm. The temperature was determined with a thermocouple placed directly behind the samples.

## III Comparison of Material Properties

A wide range of elastomer compositions and topologies (summarised in the Table 1) was investigated, in order to probe different aspects of nematic elastomer response to UV-irradiation. It was universally found that the main consequence of UV-irradiation was the cis-trans isomerisation of the azobenzene groups, with a very low probability of side reactions. However, the differing compositions and topologies of the samples resulted in a range of substantially different magnitudes of the effect. One key factor to consider is the overall concentration of isomerisable azobenzene groups, which varied between 9 and 100mol% in our materials. Another crucial factor is the topological role of the crosslinker (cf. Fig. 2). The small, flexible crosslinking polymer, DiUB is deemed to have a minor effect on the overall mesogenic properties of the liquid crystal elastomer. This is quite unlike the other two crosslinkers \- DiMC is a highly anisotropic main-chain nematic polymer containing $``$75 rod-like monomer units, which can align itself and fold into hairpin defects in the nematic phase. DiUAzo contains a single, rod-like photo-isomerisable unit of its own; its effect on the mechanical response of the elastomer network to UV-irradiation was the subject of the pioneering study by Finkelmann et al. pereira01 .

### III.1 Uniaxial thermal expansion

Before their response to ultraviolet light was investigated, the spontaneous uniaxial thermal expansion of each of the samples was first characterised (Fig. 3). This determines the extent of the coupling between the local nematic order, described by the scalar order parameter $`Q(T)`$, and the average elastomer network anisotropy, measured as the spontaneous uniaxial extension $`\lambda =L/L_0`$, where $`L`$ is the current length of the sample, and $`L_0`$ its length in the isotropic phase (the latter is always shorter for an elastomer with a prolate backbone configuration, as is the case for all our materials).

Empirically, one can relate the uniaxial thermal expansion of a liquid crystal elastomer, $`\lambda `$, to the underlying nematic order parameter $`Q(T)`$, which can be independently determined by either wide-angle X-ray scattering, or from the optical birefringence of the elastomer. It has been shown static01 ; greve01  that the relationship between the macroscopic shape and the nematic order of liquid crystal elastomers can be very well described by a direct linear proportionality, $`\lambda =1+\alpha Q`$. The constant $`\alpha `$ is a sensitive function of the average network anisotropy, varying from $`\alpha 0.08`$ in weakly anisotropic elastomers to $`\alpha 3`$ or greater in materials containing main-chain crosslinkers static01 . Such a relationship between the spontaneous uniaxial extension and the underlying nematic order is well-understood theoretically within the molecular model of equilibrium nematic polymer networks WT96 .

However, a conceptual problem resides in the values of the nematic order parameter itself. The thermal expansion data sets for all our materials, as well as the analogous studies by other groups, e.g. greve01 ; ratna01 , were all found to follow $`\lambda =1+\alpha Q`$ very faithfully, but with $`Q`$ varying in the manner of a virtual second-order phase transition, $`Q|TT_{\mathrm{ni}}|^\xi `$ chaikin . In our present materials the critical exponent of the apparent critical behaviour, $`\xi `$, takes values between 0.14 and 0.31 (Table 2). There are several possible reasons for this observation of an apparent critical behaviour in a system that, by symmetry arguments, should undergo a discontinuous first order transition. One likely explanation is that the kinetic response of the rubber-elastic network is such that it is unable to perform the discontinuous jump required by the changing underlying nematic order parameter. Such an explanation, although appealingly simple, has several weaknesses, the main one being that practically no dependence on cooling/heating rate is observed, which should be expected if the kinetic retardation is at work.

Another possibility is that in forming the monodomain nematic materials, an effective aligning molecular field is frozen in by crosslinking. This makes the first-order nematic-isotropic transition supercritical and, essentially, eliminates the true isotropic phase altogether: even at high temperature there would still be a small residual order remaining. This could still be the case, even though in our materials we could not experimentally detect any residual order above $`T_{\mathrm{ni}}`$ and observe a clear and well-defined heat signature of the transition on DSC.

A perhaps more intriguing possibility is the presence of quenched sources of random disorder in the nematic field (the network crosslinking points), which genuinely transforms the weak first order transition into a continuous one, an effect already discussed in the context of spin glasses cardy99 . The present paper does not attempt to address this issue. We focus instead on a different physical property of a photo-sensitive nematic elastomer and simply use this experimentally observed, continuous “critical” variation of $`Q(T)`$ as the foundation for describing the UV-induced response of such an elastomer.

It seems that a consistent value of the “critical exponent” is obtained, $`\xi 0.2`$. The anomalously small value of $`\xi `$ for Azo18MC can be attributed to the significant difference in this sample’s response time from all the others: kinetic effects are markedly relevant and the long DiMC chains have been shown to slow down, and even completely freeze their dynamics, when in the hairpin regime of their nematic conformation elias99 . This would result in the non-equilibrated thermal expansion of the sample lagging the actual, equilibrated order parameter. This sample also shows a very significant pre-transitional effect, with the residual nematic order and apparent sample extension persisting far above the critical temperature. This phenomenon has been examined in a separate study thermal01  and is due to a composite nature of this material in which the side-chain and main-chain polymer chains are, crudely, in equal proportion. Another point of interest is the presence of highly ordered phases in DiAzo100 and Azo30 at lower temperatures, where the dependence of $`Q`$, and hence $`\lambda `$, on temperature deviates from the attributed nematic behaviour, and becomes essentially constant. However, this will not be dwelt upon, as this work focuses exclusively on the nematic phase and the nematic-isotropic phase transition.

The observed dependence of $`Q`$ and $`\lambda `$ on temperature is very important in analysing and modelling the response of the liquid crystal elastomers to UV-irradiation. In doing so, one sample (Azo18), which proved to have the most pronounced response to UV-irradiation, is studied in greater detail in Section IV.

### III.2 Response to UV-irradiation

In order to compare different compositions and crosslinking topologies of the nematic elastomers in their response to UV-irradiation, all samples were equilibrated before irradiation at approximately $`10^\mathrm{o}`$C below their critical temperature $`T_{\mathrm{ni}}`$, determined by the thermal analysis in section III.1. Such proximity to $`T_{\mathrm{ni}}`$ also aims to maximise the variation of $`\lambda `$, whilst avoiding the singularity at the phase transition.

During irradiation, the sample temperature was, without exception, found to increase in an approximately exponential manner, mirroring the behaviour of the spontaneous uniaxial extension, $`\lambda `$. Similarly, as the elastomers relaxed after the radiation was switched off, the measured temperature again mirrored $`\lambda `$, decreasing back to its pre-irradiated value. In both cases, this exponential variation in temperature was found to be modulated by a small-amplitude oscillatory behaviour. Possible explanations for this thermal phenomenon involve a trivial heating effect from the UV lamp interfering with the active temperature control of the oven. However, such an additional heating power should not result in the overall rise of saturated temperature value, as seen in Fig. 4: the oven would be expected to reduce its own heating and bring the temperature back to its set value. Also, it would be an inconceivable coincidence if the characteristic rate of the oven heating system are exactly the same as that of all our different materials.

By closely examining the Fig. 4 it becomes clear that we observe the effect of spontaneous sample heating, occurring due to the exothermic nature of isomerisation of azobenzene groups. In a following section we shall describe the population dynamics between the trans and cis N=N isomers. The relative proportion of these states is controlled by the rates of forward and reverse reactions, which reach the balance in the saturated state (after a long time of irradiation). If one assumes that the heat release on the return to the trans ground state is proportional to the rate of change of cis concentration, $`\dot{n}_{cis}`$, then the observation in Fig. 4 becomes sensible. Initially, all azobenzene groups are in their trans state and the sample is thermally equilibrated with the oven (at 59.5oC in this case). On irradiation, the increasing number of cis isomers is created, and also accompanied by an ever increasing relaxation flux back to the trans state. This heats the sample, at the rate exactly matching the rate of isomerisation itself. Finally, a saturation state is reached, when the UV-stimulated trans$``$cis flux is equal to the relaxation flux: at this stage the heat release in the sample is constant and maximal. On switching the UV-light off the exothermic cis$``$trans relaxation continues to heat the sample, but with a decreasing power as the overall concentration $`n_{cis}`$ decreases.

All the samples responded to UV-irradiation in a manner similar to Azo18, illustrated in Fig. 4, differing only in the magnitude and speed of their response. The resultant contraction upon irradiation for each sample, relative to the expected contraction due to the aforementioned increase in temperature, are summarised in Table 3. As can be seen, the proportion of UV-isomerisable azobenzene rods in the sample is by no means a direct indication of its response to UV-irradiation, interpreted as the maximal sample extension achieved on irradiation.

Azo18 gave the largest response to radiation, but the very similar Azo30, with its higher proportion of azo-rods, gave a significantly smaller response. This must be due to the position of the azo-group within the corresponding mesogenic moiety - with Azo30, it is only weakly coupled to the polymer backbone by a long spacer, in contrast to Azo18, where a strong parallel alignment between the rod and the backbone is induced. This could result in the kinked cis-form of the moiety having less of a destabilising effect on the backbone anisotropy, and so the elastomer response to UV-irradiation is curtailed.

DiAzo9 is very similar to the material studied by Finkelmann et al. pereira01 . Despite containing less than a third of the proportion of azobenzene rods found in Azo30, its response to UV-irradiation is comparable in magnitude, which should be attributed to the crucial position of the photo-isomerisable group, in the centre of the predominant, short crosslinker.

Azo18MC might have been expected to have a very large response compared to the other materials due to its large backbone anisotropy and large coupling constant to the underlying order parameter, illustrated by its large thermal expansion amplitude. However, this was definitely not the case for UV irradiation. One may speculate the cause may be not only that the actual proportion of azobenzene rods was quite low for this sample (much of the volume is occupied by the main-chain polymer crosslinker), resulting in a smaller absolute change in the order parameter, but more importantly, that large DiMC chains can essentially freeze the elastomer dynamics in response to changes in $`Q`$ elias99 , quenching its manifestation through $`\lambda `$.

All the mesogenic groups and crosslinkers in DiAzo100 were photo-sensitive, and so one might have expected an exceptional mechanical response to UV-irradiation. However, the critical temperature, $`T_{\mathrm{ni}}`$, of this sample was much higher than that of the other materials, and to maintain consistency between experiments on different samples this required an initial equilibrium temperature of approximately $`110^\mathrm{o}`$C. At such a temperature the thermally-activated relaxation taking the azobenzene groups from their metastable kinked cis-shape back to the rod-like trans-shape would be much accelerated. As a result, one is perhaps unable to isomerise a significant proportion of azobenzene rods at any given moment of time. A lower equilibrium population of cis-isomers would then result in a smaller overall change in $`Q`$, and hence in $`\lambda `$.

The solid curves in Fig. 4 are simple exponential fits to the data, for both the UV-on and the UV-off states. Apparently, it is a reasonably good fit. However, the parameters are rather arbitrary and appear to be quite unrelated to the corresponding spontaneous thermal expansion curves. It will shortly be seen that a different, universal theoretical dependence on temperature and irradiation time describes the data with a much greater degree of fidelity. Accordingly, the solid lines in Fig. 4 serve only as a rough guide to the eye.

## IV Manipulation of Nematic Order

### IV.1 Population dynamics

The internal energy of rod-like trans, and the kinked cis mesogenic moieties is adequately described as a two-level potential well, with the former configuration energetically more stable than the latter (Fig. 5). In the experiments conducted, there are three possible ways of moving between states, each characterised by a separate time constant:
$``$ Thermal excitations up, from trans to cis, with the characteristic time given by $`\tau _{tc}=\tau _0e^{U/kT}`$, with $`\tau _0`$ the bare attempt rate.
$``$ Thermally driven relaxation down, from cis to trans, with the time constant $`\tau _{ct}=\tau _0e^{\mathrm{\Delta }/kT}`$.
$``$ UV-isomerisation, occurring with the rate $`1/\tau _{\mathrm{uv}}\eta `$ radiation intensity.

One expects that the barrier height $`U`$ is much greater than the metastable state depth $`\mathrm{\Delta }`$, and hence one can neglect the rate of spontaneous, thermally driven trans to cis transitions, although certainly not its reverse.

Let us assign the total (fixed) number density of photo-isomerisable rod-like groups in a given elastomer as $`n_0=n_{trans}+n_{cis}`$, with the only independent variable defined as $`nn_{trans}`$, so that $`n_{cis}=n_0n`$. Thus, the rate of change in the population of trans isomers can simply be written as;

$$\frac{}{t}n(t,T)=\eta n\frac{1}{\tau _{tc}}n+\frac{1}{\tau _{ct}}(n_0n)$$
(1)
In the above approximation, the relaxation time $`\tau _{tc}`$ is assumed to be very large, and can therefore be neglected. The appropriate initial condition is then $`n_{t=0}=n_0`$ (that is, all azobenzene groups are in their rod-like trans conformation) and the dynamic solution to Eq. 1 is simply;

$`n(t,T)`$ $`=`$ $`{\displaystyle \frac{n_0}{1+\tau _{ct}\eta }}\left[1+\tau _{ct}\eta \mathrm{exp}\left(t/\tau _{\mathrm{eff}}\right)\right]`$
$`\mathrm{with}`$ $`\tau _{\mathrm{eff}}=\left(\eta +1/\tau _{ct}\right)^1`$
Hence, one can see that the increasing concentration of kinked cis-isomers, which act as impurities for the nematic ordering, will eventually reach a plateau due to the opposing thermal back-reaction, itself characterised by the rate $`1/\tau _{ct}`$. Both the rate and the final saturation value of the isomer density are crucially dependent on the irradiation rate $`\eta `$; this could be different for different materials, sample shapes and experimental conditions as it is affected, among other factors, by the radiation absorption in the medium.

The relaxation process, after switching the radiation off (that is, setting $`\eta =0`$), is the thermally driven flux of cis-isomers into the trans ground state over the barrier $`\mathrm{\Delta }`$. Assuming a fully saturated irradiated state, with $`n_{cis}=n_0n_{\mathrm{}}=\tau _{ct}\eta /(1+\tau _{ct}\eta )`$ at $`t=0`$ (and still neglecting the spontaneous thermal trans$``$cis process), results in the relaxation law

$`n(t,T)`$ $`=`$ $`n_0\left[1{\displaystyle \frac{\tau _{ct}\eta }{1+\tau _{ct}\eta }}\mathrm{exp}\left(t/\tau _{ct}\right)\right].`$ (3)
Naturally, the relaxation time scale is purely $`\tau _{ct}`$, a property of an individual azobenzene group, which should remain at least approximately universal between all different samples under study.

### IV.2 Macroscopic shape and nematic order

The influence of impurities on the nematic order is well studied degen . One of the simplest ways of representing their effect is via the Landau-de Gennes mean-field formulation, which describes the local free energy density as a power series in the magnitude of the scalar order parameter, $`Q`$. Such an approach is commonly used in the thermodynamics of second order phase transitions, as well as in application to liquid crystals. Although it is much less suitable for the nematic-isotropic transition in liquid crystals, which should be first-order by symmetry, many useful qualitative results are obtained within the Landau-de Gennes formalism.

However, in view of many limitations of such a phenomenological mean-field theory, we adopt a different approach. In section III.1, the actual form of the dependence of the order parameter on temperature was established experimentally as being a virtual, continuous transition with the apparent critical exponent $`\xi `$ quite far from the mean-field value predicted by the Landau-de Gennes approach (and certainly appears a continuous transition). Therefore, we simply take the empirical dependence, $`Q|TT_{\mathrm{ni}}|^\xi `$, and modify it by substituting in the impurity-shifted critical temperature, $`T_{\mathrm{ni}}=T_{\mathrm{ni}}[n(t,T)]`$.

The effect of impurities on the critical temperature is, in the first approximation, a linear function of the population of cis-isomers. The size of this population, $`n_{cis}`$, defined above as $`n_0n(t)`$, is itself a function of time;

$$T_{\mathrm{ni}}(n)=T_\mathrm{c}\beta [n_0n(t)],$$
(4)
where $`T_\mathrm{c}T_{\mathrm{ni}}^{(0)}`$ is the transition temperature in the “pure” material containing only the trans rod-like units (given in the Table 2 for our materials). Combining this with the known dependence of the spontaneous uniaxial sample extension, $`\lambda `$, on $`Q`$, one obtains the following relation to describe the observed mechanical response to UV-irradiation

$`\lambda _{\mathrm{uv}}(t,T)=`$ (5)
$`=1+\alpha \left[T_\mathrm{c}\beta n_0\left(1{\displaystyle \frac{\tau _{ct}\eta }{1+\tau _{ct}\eta }}(1+e^{t/\tau _{\mathrm{eff}}})\right)T\right]^\xi `$
Plotting $`\lambda `$ against temperature for different values of the intensity parameter $`\eta `$, Fig. 6, one sees that the expected consequence of UV-irradiation is simply to lower the critical temperature of the nematic-isotropic phase transition.

The relaxation of sample shape after switching the radiation off is determined by the rate of recovery of the initial population of trans isomers, according to Eq. (3);

$`\lambda _{\mathrm{off}}(t,T)=1+\alpha \left[T_\mathrm{c}\beta n_0{\displaystyle \frac{\tau _{ct}\eta }{1+\tau _{ct}\eta }}e^{t/\tau _{ct}}T\right]^\xi ,`$ (6)
where, as in the Eq. (3), we have assumed that a fully saturated state was achieved on preceding irradiation.

One important parameter has to be used in the further analysis: the relaxation time scale of the cis$``$trans process. From the data in Fig. 4 for the Azo18 at $`59^\mathrm{o}`$C it is estimated as $`\tau _{ct}1479`$ s, assuming that the dynamics is dominated by the rate of isomerisation (and not, for instance, the viscoelastic relaxation of a rubbery network). In the following analysis of experimental data, the effect of (Arrhenius) temperature variation of $`\tau _{ct}`$ within a given data set are, in the first approximation, neglected. This is justified by the relatively narrow range of absolute temperatures over which a given data set varies. However, across data sets, $`\tau _{ct}`$ is shown to vary significantly with temperature.

### IV.3 Nematic destabilisation by UV radiation

Applying the Eq. (5) to the data in Fig. 4, one obtains a remarkably improved fit (Fig. 7). In analysing the data we take the fixed values for the coefficient $`\alpha =0.22`$ and the exponent $`\xi =0.19`$, as determined by the thermal expansion measurements of section III.1 and the relaxation time $`\tau _{ct}=1479`$ s, as estimated from a crude exponential fit of Fig. 4. The values of current temperature $`T`$ are taken from the corresponding experimental data points. The free parameters are, thus, only the radiation rate $`\eta 2.210^4\text{s}^1`$ and the product $`\beta n_011.8`$, which is a measure of the overall concentration, and the efficiency, of azobenzene groups.

The result of such a fit for $`\lambda (t,T)`$ is very satisfying. The initial quasi-linear contraction and its subsequent rapid flattening and saturation are fully described, unlike the rough exponential fit in Fig. 4. More importantly, we now discover that what looked like noise in both the temperature and the strain $`\lambda (t)`$ in the initial graph, Fig. 4, is in fact a consistent small-amplitude temperature oscillation (probably due to the active oven control operation), which is also fully reproduced in the mechanical shape of the sample. This is a rather stringent test of the theoretical model and the parameters used to fit the experimental data.

All the other data sets for Azo18 at different temperatures are fitted equally well, as shown in Fig. 8. Once again, we fix the value of strain-order coupling coefficient $`\alpha =0.22`$ and the exponent $`\xi =0.19`$, as determined by the thermal expansion measurement for this sample, and the relaxation time $`\tau _{ct}=1479`$ s. On physical grounds, we also insist that the coefficient $`\beta n_0`$ remains constant ($`12`$) between these data sets, since it is determined by the intrinsic properties of the nematic order and the given concentration of azobenzene groups. With these stringent constraints, one obtains the values of the only remaining free parameter, the effective rate of induced isomerisation. The inverse of $`\eta `$, the characteristic time for the photo-induced trans$``$cis isomerisation, $`\tau _{\mathrm{uv}}=1/\eta `$, is plotted in Fig. 9(a) and has an average value $`\tau _{\mathrm{uv}}310^3`$s. Although there is a certain spread in these values, we find the outcome of this analysis quite satisfactory given the number of constraints imposed on the fitting, certain error in the data and, in spite of that, a clear evidence that not only the overall trend but also the small oscillations of temperature are well reproduced by the Eq. (5).

The values of the characteristic time for the reverse thermal cis$``$trans relaxation, $`\tau _{ct}`$, are consistent with those reported in the literature for azobenzene photochemistry moreuv , as well as within themselves, Fig. 9(b). The crude estimate of the cis$``$trans energy barrier from the Arrhenius fit of the relaxation time returns a reasonable value $`\mathrm{\Delta }410^{20}\text{J}10kT`$. This is in contrast with the isomerisation time $`\tau _{\mathrm{uv}}`$, which was found to be a constant, or very slowly decreasing function of temperature. However, the errors in determining $`\eta `$ were large, and there is much room for further studying its dependence on temperature.

### IV.4 The vicinity of critical point

Since the effect of UV radiation on a nematic rubber containing isomerisable azobenzene groups is clearly identified as an effective shift of critical temperature of the underlying nematic-isotropic phase transition, the question arises - what would be the response when the material is irradiated at a temperature close to the “bare” critical point $`T_\mathrm{c}`$? In this case, the increasing concentration of cis-impurities should take the system past the transition line into the isotropic phase. Continuing the study of Azo18, we now subject it to UV-radiation at a temperature $`T_{\mathrm{exp}}63^\mathrm{o}`$C, only 4o below $`T_\mathrm{c}`$.

Fig. 10 indicates that the mechanical response of the sample on irradiation, $`\lambda (t)`$, starts exactly in the same way as at all other temperatures. However, the initial steep drop in $`\lambda `$ (and in the underlying nematic order $`Q`$) rapidly reaches the point of $`\lambda =1`$ (or $`Q=0`$); the material becomes fully isotropic and no longer has any macroscopic response to the further increase of the population of cis-isomers. The irradiation continues for some time, as illustrated by the continuing rise in temperature, which nevertheless has never reached the “bare” $`T_\mathrm{c}`$. Then the UV-light is switched off and the temperature begins to fall, in the same way as all previous measurements have shown. However, nematic order is not established for some time: clearly the concentration of cis-impurities has to decrease sufficiently to allow the nematic phase to appear at a given temperature.

Accordingly, there are two simultaneous thermal effects in the material. One is the rise (in the UV-on state), or fall (in the off state), of the actual temperature of the sample. The other is the underlying shift of the critical temperature of the nematic phase transition, described by the Eq. (4) with the time-dependent concentration $`n(t)`$ given by Eqs. (IV.1) and (3), respectively (in the latter case of thermal cis$``$trans relaxation, Eqs. (3) and (6) assume that the state of full saturation has been reached on irradiation, otherwise a different initial condition has to be used, altering the prefactor to the exponential). It is straightforward to apply the corresponding analysis and determine the time when, in the UV-on state, the decreasing critical temperature $`T_{\mathrm{ni}}=T_\mathrm{c}\beta [n_0n(t)]`$ becomes equal to the increasing sample temperature. Equally, in the UV-off state, we find the time needed for the increasing $`T_{\mathrm{ni}}`$ to intersect the decreasing actual temperature $`T`$.

Applying this analysis to Azo18 in Fig. 10, one obtains that the nematic order disappears at $`t1100`$ s and subsequently re-enters the nematic phase, is $`t4850`$ s. Both times are marked with arrows in Fig. 10: clearly the agreement with the observed change in the sample length, $`\lambda =L/L_0`$, is excellent. This is another confirmation of the theoretical model, which assumes simple thermal activation laws for characteristic times, the linear dependence of critical temperature on the concentration of cis-isomers and, much less trivially, the empirical critical dependence of the nematic order parameter $`Q|TT^{}|^\xi `$.

## V Conclusions

The derivation of Eq. (5) in section IV.2, cumulating in the analysis of relaxation data sets in Figs. 7-8, makes it clear that, once the number of kinked moieties reaches a steady-state population, the resultant effect of UV-irradiation is equivalent to simply reducing the critical temperature of the nematic-isotropic phase transition by an amount

$$\mathrm{\Delta }T=\beta [n_0n_{\mathrm{}}(T)]=\beta n_0\frac{\tau _{ct}\eta }{1+\tau _{ct}\eta },$$

which is dependent on the intensity of the ultraviolet light used (through $`\eta `$).

Hence, a slightly modified form of $`\lambda =1+\alpha Q`$, making only the replacement $`T_{\mathrm{ni}}T^{}`$, should fully describe the final, steady-state length of the sample, $`\lambda _{\mathrm{uv}}(t\mathrm{})`$, as a function of temperature. All the other parameters in Eq. (5), in particular the exponent $`\xi `$, should remain unchanged. Indeed, one finds that this idea does describe $`\lambda _{\mathrm{uv}}(t\mathrm{})`$ very well, as shown in Fig. 11. The corresponding value of shift $`\mathrm{\Delta }T`$, which best fits the data in Fig. 11, is found to be $`(3.3\pm 0.7)^\mathrm{o}`$C which agrees well with the average asymptotic value of $`\beta [n_0n_{\mathrm{}}]`$ gleaned from extrapolating the data sets in Figs. 7-8, of $`(3.5\pm 1.9)^\mathrm{o}`$C.

This work, has studied a wide range of photo-sensitive materials, based on azobenzene molecular moieties, whose differing compositions and topologies are reflected in their varied mechanical response to irradiation with ultraviolet light. The fundamental details of this response were investigated, in particular its kinetics and the variation of the critical temperature near the nematic-isotropic phase transition. The obtained experimental results are very well described by a theoretical model assuming a simple linear shift in the critical temperature with the concentration of kinked cis groups, which results in a rather non-trivial model for the dependence of the underlying nematic order parameter $`Q`$ and macroscopic uniaxial mechanical extension $`\lambda `$ on irradiation time and temperature. Basic material parameters and constants were also determined, which will allow subsequent investigations to further probe the details of this fascinating phenomenon. With such a novel, non-invasive method of control over mechanical actuation, there is a great potential for application, but still very much is left unknown.

###### Acknowledgements.

We are grateful to Wacker Chemie for donating the platinum catalysts for the chemical preparations, and to H. Finkelmann and M. Warner for a number of useful discussions and advice. This work has been supported by EPSRC UK.


format: text

text:
# Time dependent correlations in marine stratocumulus cloud base height records

## I Introduction

Cloud base height (CBH) profiles are difficult object to describe and analyze. The base of a cloud is determined by the local condensation level up to which unsaturated air parcels rise. CBH has a complex structure, that varies with the coordinates of the observation point and time. Such a complexity is due to a variety of processes that take place in the marine boundary layer (MBL) formed over the sea surface. Some of these processes are, e.g., turbulent motions, entrainment, radiative transfer and changes in the cloud microphysical structure. Thus the CBH looks like an erratic quantity in time and space. Therefore, advanced methods of computational statistical analysis are required in order to extract meaningful physical information from data series with such a high level of complexity that traditional analysis techniques fail to provide. The detrended fluctuation analysis (DFA) method is designed to search for long-range correlations in non-stationary, highly fluctuating signals. The method has been successfully applied in the investigation of turbulence, DNA sequences, ionic transport trough cell membrane, foreign currency exchange rates, meteorological phenomena, and stratus clouds dynamics. Using this technique we hereby sort out correlations and anti-correlations in stratocumulus cloud base height fluctuation data records. Our findings of the time evolution of the local $`\alpha `$-exponent (which is the Hurst exponent), i.e. describing a scaling law over a certain time range, are related to the physical processes in the cloud. These results reflect the ability of the method to indicate the type of physics that underlies the CBH profile phenomenon.

## II Marine Boundary Layer - structure and processes

A variety of physical processes takes place in the atmospheric boundary layer. At time scales of less than one day, significant fluxes of heat, water vapor and momentum are exchanged due to entrainment, radiative transfer, and/or turbulence. The turbulent character of the motion in the atmospheric boundary layer (ABL) is one of its most important features. The turbulence can be caused by a variety of processes, like thermal convections, or mechanically generated by wind shear, or following interactions influenced by the rotation of the Earth. This complexity of physical processes and interactions between them create a variety of atmospheric formations. In particular, in a cloudy ABL the radiative fluxes produce local sources of heating or cooling within the mixed-layer and therefore can greatly influence its turbulent structure and dynamics. Moreover the variations in the turbulent structure and dynamics of the clouds cause subsequent changes in the cloud boundaries, especially in the cloud base.

The atmospheric boundary layer is characterized by an inner (surface) layer at heights above the aerodynamic roughness length and below one tenth of the depth of the ABL, while the outer (Eckman) layer is at higher distances within the ABL. The outer region of an unstably stratified ABL is often referred to as the mixed layer because of the dominating convective motions that take place there, generated by strong surface heating from the Sun or by cloud-topped radiative cooling processes. In contrast, a stably stratified ABL occurs mostly at night in response to the surface cooling due to long-wave length radiation emitted into the space.

The top of the boundary layer in convective conditions is often well defined by the existence of a stable layer (capping or subsidence inversion), beneath which clouds form, the so-called cloud-topped boundary layer (CTBL). In presence of clouds (shallow cumulus, stratocumulus (Sc) or stratus (St)) the structure of the ABL is modified because of the radiative fluxes. Phase changes become more important. During cloudy conditions one can distinguish mainly : (i) the case in which the cloud and the subcloud layers are fully coupled; (ii) two or more cloud layers beneath the inversion, with the lower layer well-mixed with an upper elevated layer, decoupled from the surface mixed layer or (iii) a radiatively driven elevated mixed cloud layer, decoupled from the surface.

An interesting case is that of the marine boundary layer characterized by high concentration of moisture. It is wet, mobile and has a well expressed lower boundary. The competition between the processes of radiative cooling, entrainment of warm and dry air from above the cloud and turbulent buoyancy fluxes determine the state of equilibrium of the cloud-topped marine boundary layer. The Atmospheric Stratocumulus Transition Experiment (ASTEX) was designed to clarify the transition from stratocumulus to trade cumulus clouds in the MBL in the region of the Azores Islands. Several papers have been dedicated to the relevance of the main processes causing this transition. In Ref.  Betts et al. found a relationship between the thermodynamic structure of the marine boundary layer and the diurnal variation of the cloudiness.

## III DFA method and spectral analysis

A method that relaxes the requirement of stationarity of the investigated signal is the detrended fluctuation analysis (DFA) method. The DFA method is a tool used for sorting out long range correlations in a non-stationary self-affine time series with stationary increments. The method has been used previously in the meteorological field. It provides a simple quantitative parameter \- the scaling exponent $`\alpha `$, which is a signature of the correlation properties of the signal. The Detrended Fluctuation Analysis technique consists in dividing a random variable sequence $`y(n)`$ of length $`N`$ into $`N/t`$ (non-overlapping) boxes, each containing $`t`$ points ($`N/t`$ = 4,5, …). Then, the trend (assumed to be linear in this investigation, but it can be generalized without any difficulty) $`z(n)=an+b`$ in each box is computed using a linear least-square fit to the data points in that box. The detrended fluctuation function $`F^2(t)`$ is then calculated following :

$$F^2(t)=\frac{1}{t}\underset{n=kt+1}{\overset{(k+1)t}{}}\left[y(n)z(n)\right]^2\text{for}k=0,1,2,\mathrm{},\left(\frac{N}{t}1\right).$$
(1)
Averaging $`F^2(t)`$ over the $`N/t`$ intervals gives the mean-square fluctuations

$$<F^2(t)>^{1/2}t^\alpha .$$
(2)
If the $`y(n)`$ data are random uncorrelated variables or short range correlated variables, the behavior is expected to be a power law with an exponent $`\alpha `$=1/2 if the fluctuations are not correlated. An exponent $`\alpha `$ $``$ 1/2 in a certain range of $`t`$ values implies the existence of long-range correlations in that time interval as, for example, in fractional Brownian motion. A small value of $`\alpha `$ indicates antipersistence of correlations, i.e. a positive fluctuation is more likely to be followed by a negative one than a positive one. An $`\alpha `$ value between 0.5 and 1 indicates persistence of correlations, i.e. a positive fluctuation is more likely to be followed by a positive one than a negative one. The $`\alpha `$-exponent value that holds true for a certain time interval called scaling range, is a characteristic of the correlations in the fluctuations of a signal $`y(t)`$ defined between the beginning $`t_0`$ and end of the observations $`t_M`$, i.e. $`[t_0,t_M]`$. The value can vary in time if the process is not stationary and/or the scaling range is finite (see Sect.6). The $`\alpha `$ exponent is related to the usual fractal dimension $`D`$.

The concept of the invariance of the fractal shapes at different scales (self-similarity) can be transposed in the statistical analysis of a time series. One searches for a given scale range where the statistical properties of the signal are invariant. By definition the time series is self-affine if its spectrum has a power law dependence.

$$S(f)(1/f)^\beta .$$
(3)
It has been shown by Heneghan and McDarby that the relationship $`\beta =1+2\alpha `$ holds true for stochastic processes, i.e. for fractional Brownian walks. Depending on the value of the spectral exponent $`\beta `$, one has the following cases : if $`\beta <1`$, the process is stationary; if $`\beta >1`$, the process is non-stationary; if $`1<\beta <3`$ the process is non-stationary with stationary increments. In the case of white noise $`\beta =0`$ while $`\beta =2`$ for a Brownian walk signal.

## IV Data

Cloud base height data has been obtained during the Atlantic Stratocumulus Transition Experiment (ASTEX) in June 1992, at the Azores Islands using a laser ceilometer. The wavelength of the vertically transmitted into the atmosphere laser beam is 0.904 mm. The ceilometer works with a 15 m spatial resolution and reaches the maximum measurable height of 4 km. The time resolution of CBH records is 30 seconds. The cloud base height values were selected from the profiles of the backscattered signal intensity, using manufacturer’s algorithm based on a voltage threshold for the maximum value of the backscattered radiation at a height giving the local condensation level up to which unsaturated air parcels rise. We use data records measured on June 15, 1992 (Fig. 1(a)) and June 18, 1992 (Fig. 2(a)) to determine (i) the scaling properties of CBH fluctuations; (ii) to classify the type of time correlations, (iii) to investigate the diurnal changes in the CBH evolution, and (iv) to identify the type of the time dependent long-range correlation evolution.

According to, on June 15, 1992 an anticyclone was centered north of Santa Maria - one of the Azores Islands, where the measurements were taking place, producing a strong subsidence inversion beneath which were observed marine boundary layer clouds. The nocturnal MBL was found to be well-mixed, generally coupled to the surface, with a solid layer of stratocumulus, while the daytime MBL was observed to be decoupled from the surface, because of the coexistence of a stratocumulus cloud layer and a cumulus marine subcloud layer. After the sunrise the MBL became intermittently decoupled and CBH varied dramatically as a cumulus formed at the top of the marine subcloud layer. In this decoupled MBL, the cumulus convection, that arises from the latent instability generated in the marine subcloud layer, acts toward a recoupling of the formally decoupled cloud layers that persist beneath the trade inversion. Five sequences of strong decoupling were observed that day (Fig. 1(a)).

For June 18, 1992 we do not have any information about the meteorological situation, but in analogy with data in Fig. 1(a), the CBH profile on Fig. 2(a) leads one to distinguish periods of increased variability of CBH, followed by short CBH fluctuations.

## V Scaling properties of the CBH data

The CBH evolution, reported in Figs. 1(a) and 2(a) through $`S(f)`$, can be regarded as a non-stationary process with stationary increments, since their spectral exponents are $`\beta =1.28\pm 0.1`$, for June 15, 1992 (Fig. 3) and $`\beta =1.49\pm 0.08`$, for June 18, 1992 (Fig. 4). The data error bars are calculated from the r.m.s deviations in the fit according to standard techniques and have their $`95\%`$ confidence interval conventional meaning.

The type of correlations of the CBH fluctuations can be probed using the Detrended Fluctuation Analysis (DFA) method. Applying the DFA method, we find that long-range time correlations in CBH fluctuations exist and are of an anti-persistent type with $`\alpha =0.24\pm 0.002`$ for June 18, 1992 and $`\alpha =0.21\pm 0.005`$ for June 15, 1992 (Fig. 5). They hold for a time interval approximately equal to 140 min in the first case and 40 min for the second. In this scaling range, the relationship $`\beta =1+2\alpha `$ is reasonably satisfied. The flat like regions of both curves (Fig. 5) correspond to low values of the Hurst exponent $`\alpha `$ ($`\alpha =0.06\pm 0.003`$ for June 15, and $`\alpha =0.10\pm 0.003`$ for June 18) and indicate the presence of $`1/f`$-like noise.

For completeness, we show in Fig. 6, the distribution function of the CBH fluctuations for the 30 sec time lag interval. Notice the similarity between both cloud cases and the marked departure from the Gaussian behavior, as numerically exemplified in the $`S(f)`$ and DFA correlation methods in fact. A similar double triangular pyramid shape has been observed in other meteorological cases.

## VI Time correlation analysis

It is of interest to test if the correlations remain the same for shorter intervals accommodated into $`[t_0,t_M]`$ or if they change with time, as should be anticipated for nonstationary time series data. In order to probe the existence of so called locally correlated and decorrelated sequences, one can construct an ”observation box” with a certain width, $`\omega `$, of size 7h, place the box at the beginning of the data, calculate $`\alpha `$ for the data in that box, move the box by $`\mathrm{\Delta }\omega `$= 30 min (or 60 points) toward the right along the signal sequence, calculate $`\alpha `$ in that box, a.s.o. up to the $`N`$-th point of the available data. A time dependent $`\alpha `$ exponent may be expected, and is given at the end of each box, thus for $`t`$ ranging from $`\mathrm{\Delta }\omega `$ to $`N`$.

Thus we obtain an ”instantaneous measure” for the degree of long-range correlations. Results for this instantaneous $`\alpha `$-exponents are shown in Figs. 1(b) and 2(b) for the cloud base height data on June 15, 1992 and June 18, 1992. Note that each $`\alpha `$ value represents the behavior resulting from all points in the box, which continuously overlap. This approach seems suitable to cloud base height data because it can be expected to reveal changes in the correlations dynamics of the clouds at various times for a given time lag $`\mathrm{\Delta }\omega `$.

The $`\alpha `$ evolution is within the limits of a quite anti-persistent behavior with a transition from weak to strong, at various instant, in other words, anti-persistence is found in the underlying dynamics. Five regions in the CBH evolution, corresponding to the description given in section 2, can be distinguished for the CBH data measured on June 15. The first region lasts from 0 to 6 h, the second from 6 to 13 h, the third holds from 13:00 to 16:00 h, the fourth - from 16:00 to 20:00 h and the fifth one from 20:00 to 24:00 h. The transition occurring from 08:00 to 10:30 in the evolution of the $`\alpha `$ exponent indicates the increasing role of a process that generates the second structure in the CBH that appears between 10-13 h in Fig. 1(b). This process is characterized by a mean value of the local $`\alpha =0.24\pm 0.01`$ and is understood as a process of $`destabilization`$, in the decoupled MBL. Between 14-18 h the downward trend of the $`\alpha `$ exponent indicates the increasing dominance of a process that generates the third part of CBH evolution that exists between 17-19 h in Fig. 1(b). The mean value of the local a exponent is $`\alpha `$ =0.20$`\pm `$0.02. From 19 h to 20:30 h a very slow upward trend leads to the process of the fourth region between 21:30 h - 24 h in Fig. 1(b), with a mean local $`\alpha `$ =0.26$`\pm `$0.02.

For June 18 we distinguish three main regions in the CBH evolution: 0-10 h, 10-19 h and 19-24 h. The trends in the $`\alpha `$ exponent evolution are as follows: 0-10 h - process, generating the first CBH structure, with mean local $`\alpha `$ =0.17$`\pm `$0.02. This corresponds to the dynamics of the MBL before its decoupling. For the period 10-16 h a large upward transition leads to a weak anti-persistent correlation in CBH fluctuations with a mean value of the local $`\alpha `$ =0.35$`\pm `$0.02. Hence the dominant process is destabilizing the MBL structure. For the time 16-20 h the dynamics of the process changes again, characterized by a mean local $`\alpha `$ =0.26$`\pm `$0.01 and a return to more stable MBL structure.

Thus the evolution of the Hurst exponent $`\alpha `$ in the cases that we have investigated, consists in the successive transitions from lower to higher values and vice versa. This corresponds to transitions from strong to weak anti-persistent correlations and reflects the changes in the underlying dynamics of the MBL.

Notice that the DFA function being a measure of the root mean square deviation from the linear best fit (a ”local trend”) is therefore a measure of a height-height correlation function. The so-called height is nothing else than a mass threshold. Therefore this DFA function informs on a corresponding physical property, i.e. the local compressibility (or bulk modulus) of the cloud fluid.

## VII Conclusions

We have presented an analysis of a marine stratocumulus cloud base height data series using the DFA method and the spectral density method in order to investigate the scaling range of time correlations in the CBH fluctuations searching for scaling laws and persistence effects. In all cases anti-persistent correlations hold with a well established day-night evolution measured by the instantaneous $`\alpha `$ index , - itself characterizing the dominant physical process. The transition, revealed by the trend in the $`\alpha `$ evolution, between processes inducing different levels of anti-persistent correlations of the fluctuations could be used as a trace parameter of the CBH dynamics. The accuracy of these statistics in CBH fluctuations suggests the usefulness of similar studies for the radiation transfer dynamics.

The DFA method allows one to obtain the scaling ranges and crossovers of the correlations between the fluctuations of a signal. The $`\alpha `$ exponent is in fact directly related to that for the auto correlation function. Incidently, if the relationship $`\gamma =2(1\alpha )`$ is supposed to hold for describing the power law dependence of the autocorrelation function, a value $`\gamma `$ = 1.56 is expected. This means that the CBH fluctuations can be empirically represented by a fractional Brownian motion process. The oscillations between different stability regimes could thus be mapped onto a Kramer problem, itself mimicking the basic Langevin equation for the signal increments.

Nevertheless this knowledge of the second moment of a signal is not usually enough in order to have a full understanding of the physical process nor of its underlying physical mechanisms governing the CBH. It would be of interest to couple these results to some analysis of the local (vertical and horizontal) velocity fluctuations, and of the temperature fluctuations in order to obtain additional information, whence the viscosity and thermal conductivity, for a better understanding of cloud physics and weather predictability.

## VIII Acknowledgments

This research was partially supported by Battelle grant number 327421-A-N4.


text:
# Human Sexual Contact Network as a Bipartite Graph

## I Introduction

Yet another human interaction network! *The web of human sexual contacts* is introduced in , where the authors make extensive analysis of 1996 Swedish survey of sexual behaviour and show that the cumulative distributions of number of sexual partners for males and females have power-law forms with exponents $`\alpha _f2.54`$ for females and $`\alpha _m2.31`$ for males, during a single year and $`\alpha _{f_{tot}}2.1`$, $`\alpha _{m_{tot}}1.6`$ for entire life time.

Human sexual behaviour is a subject area  in itself and beyond the scope of this work. Here we will concentrate on capturing the essential details needed to form such a social network for which we will make use of the developments from the previous network models  and general consensus about human sexual behaviour.

The model introduced in  has identified two basic principles that seem to govern the growth dynamics of human interaction networks, namely preferential attachment and continuous growth. This differs from the traditional complex network consideration  of connecting together a fixed number of network elements (nodes) which results in Poisson distributions of the connectivity.

The usual approach in most complex network models  that have human attributes is to employ the preferential attachment process; an incoming node is more likely to connect to a node with higher connections. This process yields power-law distributed connectivity in the network. A similar process takes place in a sociophysics model , where probability of winning a fight depends on the number of victories, consequence of which is a self-organised hierarchy in the system. Again a well studied biological evolution model , in which the power-law behaviour is due to increasing fitness of the species.

Contrary to previous human interaction network models, the model we consider here has two different types of nodes and the linking process is discriminative as well as preferential. In the proceeding section, the model is built by connecting new nodes to nodes with higher connections and analysed by a rate equation approach. The results and the discussions are presented in the final section.

## II The model

The growth process of the sexual contact network described by a random growing bipartite graph (illustrated in Fig. 1) considers the following:

* With probability $`q`$ a new (no previous connection) female is connected to an existing male, which contributes to the number of females with contact degree $`l=1`$ in the network and increases the degree $`k`$ of the connecting male.
* With probability $`p`$ a new male is connected to an existing female, which has the analogous effect as in (i).
* With probability $`r`$ a link is created between an existing male and an existing female, which increases the degrees of both sexes.

In addition to above rules, initially there must be at least one connected male-female pair and only heterosexual contact is allowed among the sexes. Also the process $`r`$ only counts once between the same pair; that is once a link between two individuals is created it remains and further contact between them does not increase their contact degrees.

We define male sexual degree distribution $`A_k`$ as the number of males with $`k`$ female connections. The rate at which a male attracts females is proportional to his number of connections and it grows linearly with increasing $`k`$; in other words the connection probability of a male increases as he becomes more experienced. And similarly does the female sexual degree distribution $`B_l`$ with attraction rate proportional to $`l`$. However, it seems necessary to consider some characteristics unique to both sexes so that we may underpin the social expectations about their sexual contacts. Here we introduce $`\lambda `$ and $`\mu `$ as the *charisma* parameters of males and females respectively. These parameters should not be confused with *charm* or *beauty*, which could be used as attractiveness of individuals. Here we are considering only the general attractiveness of the sexes, which are shaped by social expectations.

The charisma parameters are added into the rate of connections given $`k+\lambda `$ for males and $`l+\mu `$ for females. Since we need $`k+\lambda >0`$ and $`l+\mu >0`$ for the network to grow, $`\lambda `$ and $`\mu `$ are considered to be positive integers and small compared to the contact degrees of $`k`$ and $`l`$. Also, the general consensus is that males tend to exaggerate their number of sexual partners. This desire puts females in demand and increases female charisma. Conversely, the social expectation of females restricts them to report less and partly because of this they are more selective about their sexual partners. Using this information we argue that the rate at which females attract sexual partners should be greater than that for males and therefore $`\mu >\lambda `$.

In this model we are considering a large system where the random fluctuations in the linking process can be neglected when compared with the collective behaviour of the individuals. Therefore it is appropriate to use mean-field approximation. Then we can write the time evolution of the average number of males with $`k`$ contacts $`A_k(t)`$ as

$`{\displaystyle \frac{dA_k(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{q+r}{M_a}}\left[(k1+\lambda )A_{k1}(t)(k+\lambda )A_k(t)\right]`$ (2)
$`+p\delta _{k1}`$
and similarly for females with $`l`$ contacts,

$`{\displaystyle \frac{dB_l(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{p+r}{M_b}}\left[(l1+\mu )B_{l1}(t)(l+\mu )B_l(t)\right]`$ (4)
$`+q\delta _{l1}.`$
Since the structure of Eq. (2) and Eq. (4) are the same we will only analyse the former equation and apply results to the latter.

The first term in the square brackets indicates the contribution to $`A_k`$ when a male of $`k1`$ connections acquires a new contact, the corresponding losses are given by the second term. The growth of the network due to influx is expressed by the last term, with $`p`$ being a probability of addition of males with $`k=1`$ connection. The multiplicative factor $`M_a`$ is to ensure the appropriate normalization and is given by

$`M_a(t)`$ $`=`$ $`{\displaystyle \underset{k}{}}(k+\lambda )A_k(t).`$ (5)
The total number of males in the system is $`A(t)=_kA_k(t)`$ and the total number of females $`B(t)=_lB_l(t)`$. From initial conditions we have $`A(t)=1+pt`$ and similarly $`B(t)=1+qt`$ making an assumption that there is one of each sex representatives before the growth commences. In particular, the bipartite structure  of the graph ensures that the total number of connections in the system is the same for both males and females, which means that $`L(t)K(t)=_kkA_k(t)=t+1`$, with one link before the first time step. Also, the network does not contain monogamous individuals after the first time step.

We start our analysis by looking at the moments of $`A_k(t)`$, which are defined as

$`M_i(t)`$ $``$ $`{\displaystyle \underset{k}{}}k^iA_k(t).`$ (6)
For the first few moments it is easy to show that

$`M_0(t)`$ $`=`$ $`M_0(0)+(q+r+p)t`$ (7)
and

$`M_1(t)=M_1(0)+pt.`$ (8)
For large times the initial values of the moments become irrelevant so that we get

$`M_a(t)`$ $`=`$ $`[q+r+p(1+\lambda )]t,`$ (9)
which is a linear function of time and indicates that Eq. (2) is also linear in time. Using this input we can write $`M_a(t)=m_at`$ implying $`m_a=q+r+p(1+\lambda )`$ and $`A_k(t)=a_kt`$. Substituting these relations in Eq. (2) leads to a time independent recurrence relation

$`\left({\displaystyle \frac{m_a}{q+r}}+k+\lambda \right)a_k=(k1+\lambda )a_{k1}+{\displaystyle \frac{m_a}{q+r}}p\delta _{k1}.`$ (10)
Iterating Eq. (10) we obtain

$`a_k={\displaystyle \frac{\mathrm{\Gamma }(k+\lambda )}{\mathrm{\Gamma }(\frac{m_a}{q+r}+k+\lambda +1)}}{\displaystyle \frac{\mathrm{\Gamma }(\frac{m_a}{q+r}+\lambda +1)}{\mathrm{\Gamma }(\lambda )}}{\displaystyle \frac{m_a}{q+r}}p.`$ (11)
Using the former relation for the case of females, by simply replacing the relevant parameters gives

$`b_l={\displaystyle \frac{\mathrm{\Gamma }(l+\mu )}{\mathrm{\Gamma }(\frac{m_b}{p+r}+l+\mu +1)}}{\displaystyle \frac{\mathrm{\Gamma }(\frac{m_b}{p+r}+\mu +1)}{\mathrm{\Gamma }(\mu )}}{\displaystyle \frac{m_b}{p+r}}q.`$ (12)
For large $`k`$ and $`l`$ Eq. (11) and Eq. (12) are asymptotically equivalent to

$`a_kk^{(1+\frac{m_a}{q+r})}`$ (13)
and

$`b_ll^{(1+\frac{m_b}{p+r})}.`$ (14)
The distributions $`a_kk^{\gamma _m}`$ and $`b_ll^{\gamma _f}`$ scale as power-laws, which is a vulnerable structure for social networks as far as the epidemics and the spreading of computer viruses are concerned . The scaling exponent

$`\gamma _m=2+{\displaystyle \frac{p(1+\lambda )}{1p}}`$ (15)
depends upon male charisma parameter $`\lambda `$ and the rate $`p`$ of addition of males into the network. Similarly, the exponent

$`\gamma _f=2+{\displaystyle \frac{q(1+\mu )}{1q}}`$ (16)
depends upon female charisma parameter $`\mu `$ as well as the rate of arrival of females. However, since we have the constraint $`p+q+r=1`$ implying one link at each time step, the growth of one network is effected by the growth of the other.

## III Discussion and conclusions

From the definition of the model, we have four free parameters, $`0p1`$, $`0q1`$, $`\mu >\lambda >0`$ and the fifth one is given by the relation $`(p+q+r)=1`$.

Starting with a simple example, let us ignore the process $`r`$ and take $`p+q=1`$. This yields

$`\gamma _f=2+{\displaystyle \frac{q(1+\mu )}{1q}}\text{and}\gamma _m=2+{\displaystyle \frac{(1q)(1+\lambda )}{q}}.`$ (17)
We get $`\gamma _f2`$, $`\gamma _m\mathrm{}`$ as $`q0`$ and $`\gamma _f\mathrm{}`$, $`\gamma _m2`$ as $`q1`$. Although this is an unrealistic assumption, it simply illustrates the strong coupling of the influxes of males and females into the system. A similar strong coupling is also observed in  where a more general coupled growing network is considered.

If we pick $`\lambda =2`$, $`\mu =3`$ and substitute into the scaling exponents in Eq. (15) and Eq. (16) respectively we get $`p0.30`$ and $`q0.28`$ for the corresponding power-law exponents in . Where $`\gamma _m=1+\alpha _m=3.31`$ and $`\gamma _f=1+\alpha _f=3.54`$ for the short time behaviour analysis. In this example choice of values for $`\lambda `$ and $`\mu `$ allowed us to obtain $`pq`$ meaning that the increase of links due to new males and females is approximately the same and about 40 % of the links are created among the pre-existing males and females. The average number of partners of the sexes $`_kkA_k/_kA_k=1/p3`$ is the same since we have a close system.

Let us now consider the distribution of contacts during the life time of the sexes. Here, we introduce another parameter $`\beta `$ to the attachment rates to describe bad experiences, aging or inactivity of the sexes and hence the new parameter is a negative input. It can be in the range $`[0,k+\lambda )`$ or $`[0,l+\mu )`$, which seems reasonable to assume for $`k,l>2`$. We can denote $`\lambda _1=\lambda \beta _m`$ and $`\mu _1=\mu \beta _f`$ where we assign a unique $`\beta `$ parameter to each kind to discriminate their sexual behaviour. The new scaling exponents

$`\gamma _f=2+{\displaystyle \frac{q(1+\mu _1)}{1q}}`$ (18)
and similarly

$`\gamma _m=2+{\displaystyle \frac{p(1+\lambda _1)}{1P}}.`$ (19)
Using the same values we used earlier for the free parameters gives $`\gamma _f=1+\alpha _{f_{tot}}=3.1`$ for $`\beta _f1.17`$ the same for all females and $`\gamma _m=1+\alpha _{m_{tot}}=2.6`$ for $`\beta _m1.6`$ for all males. Consequently we have $`\beta _f<\beta _m`$ implying that the sexual connections of males tend to suffer more from aging, bad experiences or inactivity than that of females. Perhaps this may be the reason why males resort to lies to maintain self image.

In conclusion, our simple model seems adequate to simulate the dynamics of the web of human sexual contacts. The gratifying feature of rate equation approach is that it gives flexibility to include diverse range of parameters to model any situation within the frame of this work. The charisma parameters used here are shaped by the western social expectations and they can be changed to reflect the views of other societies. Also, one can introduce *fitness*, multiplicative quenched disorder into the system as in  to account the activity level of different age groups. Finally, a challenging issue to address would be to study male-female degree correlations of this model.

###### Acknowledgements.

I would like to thank EPSRC for financial support, Dafang Zheng, Geoff Rodgers and Yan Fyodorov for useful discussions.


text:
# Improving the Prospects for Detecting Extrasolar Planets in Gravitational Microlensing Events in 2002

## 1 Introduction

The importance of well-aligned or high magnification microlensing events for detecting planetary companions to the lens star was first pointed out by Griest & Safizadeh (1998). They demonstrated that Jupiter-mass planets are readily detectable (if present) in events with maximum amplication, $`A_{\mathrm{max}}`$, as low as 10, and that Neptune-mass planets are detectable in events with $`A_{\mathrm{max}}=50`$, if they are monitored intensely around their times of peak magnification. Unfortunately, the vast majority of the more than 1000 microlensing events that have now been detected by survey groups \[Alcock et al. 2000, Derue et al. 2001, Udalski et al. 2000, Bond et al. 2001a\] were of low magnification. We are aware of only one event with high magnification that received intensive observational coverage at its peak. This event, MACHO 98–BLG–35, reached $`A_{\mathrm{max}}80`$. The observations yielded large exclusion regions for gas-giant planets surrounding the lens star, and also evidence for an Earth-mass planet near its Einstein ring \[Rhie et al. 2000, Bond et al. 2001b\]. Two other high magnification events, OGLE 00–BUL–12 and MACHO 99–LMC–2, received less intensive coverage yet still yielded large exclusion regions for gas-giant planets \[Bond et al. 2001b\]. An upper limit of the order of 30% on the abundance of Jupiter-like planets has also been obtained by the PLANET collaboration from a study of typical microlensing events \[Albrow et al. 2001, Gaudi et al. 2002\]. This study included the event MACHO 98-BLG-35 in which they excluded Jovian planets but were unable to draw any conclusions on the presence of terrestrial planets. This is not inconsistent with the conclusions of Bond et al (2001b) since the PLANET coverage of the peak of this event was less intensive and the data were analysed using a procedure that is generally less accurate.

Significant progress could clearly be made if one had a larger sample of high magnification events to work with. The purpose of this letter is to report the detection of 10 high magnification events thereby demonstrating that high magnification events can be detected with high efficiency, and to urge follow-up observations of future events.

## 2 Observations

During 2000–2001 a campaign of observations was undertaken by the MOA collaboration with the aim of improving the detection rate of high magnification events. A 0.6m telescope at the Mt John Observatory in New Zealand (170 E, 44 S) with a mosaic camera comprised of three 2k $`\times `$ 4k thinned CCDs was used. An area 17 deg2 towards the Galactic Bulge that is relatively unobscured by dust was monitored.

The MOA microlensing search procedure involves a combination of multiple observations (up to six times) per night and real-time difference imaging to pick up microlensing events. Images taken during the pilot year of 2000 were used to build a database of variable stars and to detect some microlensing events \[Bond et al. 2001a\]. The observations in 2001 were made primarily to search for microlensing events and to provide real-time alerts111http://www.phys.canterbury.ac.nz/physib/alert/alert.html to follow-up groups. A total of 53 possible microlensing events were detected in real time in 2001, of which 10 had $`A_{\mathrm{max}}>40`$. The details of these events are given in Table 1 and their light curves are shown in Fig. 1. Most of the events had $`A_{\mathrm{max}}100`$ and all stood out clearly as shown, for example, in Fig. 2.

We determined the parameters $`A_{\mathrm{max}}`$ and the the Einstein crossing time $`t_\mathrm{E}`$ by fitting the standard single lens microlensing profile given by Paczynski (1986). The constraints on these parameters were determined by a thorough examination over a range of values of $`A_{\mathrm{max}}`$ and $`t_\mathrm{E}`$. For some events only lower limits could be determined. The large uncertainty in $`A_{\mathrm{max}}`$ for some of the events in Table 1 was mainly due to the less than complete coverage of the peak of their profiles. The intrinsic faintnesses of the sources also contributed to the uncertainties. With the exception of the first event, all events shown in Table 1 had faint baseline intensities at or below the detection threshold. However, since difference imaging was used, the photometry and subsequent derivation of event parameters was unaffected by blending and thus free from systematic effects that may result using conventional profile fitting photometry techniques \[Alcock et al. 2000, Wozniak 2000, Bond et al. 2001a\]

The MOA strategy of nightly multiple observations of survey fields is different from that of earlier survey projects \[Alcock et al. 2000, Udalski et al. 2000, Derue et al. 2001\]. It is noteworthy that most of the events shown in Table 1 and Fig. 1 were observable only for a few days near their peaks. Events 32, 37, 41, 46, and 50 were especially rapid with most of the brightness changes ocurring in just one night. Such events had not been seen in previous surveys. If we had adopted the strategy of just a single observation of the survey fields per night, only events 2, 5, and possibly 7 would have been picked up. The multiple nightly sampling strategy would therefore appear to give 3–5 fold increase in the detection efficiency for rapid events.

Our use of difference imaging analysis has permitted the detection of events that would have otherwise been missed due to some combination of intrinsic faintness of the source and the degree of image crowding. Previous offline analyses using difference imaging have been carried out on the image databases produced by MACHO \[Alcock et al. 2000\] and OGLE \[Wozniak 2000\] and have also shown that microlensing events with faint source stars can be detected. However, the sampling rate over the target fields in these surveys was only of the order of once per night and events like the rapid ones presented appear to have been missed.

Our survey procedure also differs from previous ones in that the difference imaging analysis is carried out on-line. This crucial feature enables events to be alerted to other observatories in real-time for follow-up. Without this, the potential for planet detection described in the following section could not be realized. The real-time feature of our procedure was also a factor in the scheduling of our survey fields. For some of the events shown in Fig. 1, the observation frequency of those fields containing the events was increased following the event alert. Unfortunately, bad weather and interruptions due to daylight prevented us from obtaining complete coverage of their peaks.

In general, the MOA survey procedure, which combines multiple observations of each of the survey fields nightly with on-line difference imaging analysis, evidently allowed us to tap into a huge reservoir of faint microlensing sources that had previously not been fully utilized. Our procedure enables detection of these events with small telescopes during the brief time when they are highly magnified, i.e. when the lens and the source are well-aligned. The fitted Einstein crossing times, $`t_\mathrm{E}`$, for these events are not, however, unusually short. The crossing time for an event depends upon the mass of the lens and the transverse velocity of the lens with respect to source, and these parameters are not significantly biased by our detection procedure. With the exception of event 5, which has an unusually long duration of $`200`$ days222Six similarly long timescale events with $`t_\mathrm{E}>200`$ days have been previously detected in the MACHO survey. These may be due to microlensing by stellar mass black holes \[Bennett et al. 2001\], all events had $`t_\mathrm{E}`$ in the range 10–40 days, consistent with the distribution found using the conventional survey technique \[Udalski et al. 2000\]. Thus we see no evidence that the MOA survey probes a different population of lens stars from that in previous surveys. However, it does probe a different population of source stars. The data in Table 1 imply the source stars for all the events except the first are sun-like when allowance is made for reddening caused by dust \[Schlegel, Finkbeiner, & Davis 1998\]. These dwarf stars are ideal for planet hunting because source-size effects that tend to wash out planetary signals are minimised \[Griest & Safizadeh 1998\].

## 3 Discussion

It is evident from the above that microlensing events of high magnification with solar-type sources can be readily detected. Griest & Safizadeh (1998) considered the detection of planets with masses ranging from that of Neptune up to that of Jupiter in microlensing events with $`A_{\mathrm{max}}`$ up to 50. We have considered the detectability of lower mass planets in events with $`A_{\mathrm{max}}50`$ such as those reported here. It is beyond the scope of this letter to consider in depth, the range of sampling rates, detection criteria, telescopes, etc, that could be used when monitoring high magnification events. We present here some calculations as a benchmark as to what could be achievable in future events similar to those presented here.

We simulated events with $`A_{\mathrm{max}}`$ of 50, 100, and 200 that reached a peak I magnitude of 15. We assumed they were monitored intensively during the time interval $`[0.5t_{\mathrm{FWHM}},0.5t_{\mathrm{FWHM}}]`$ where $`t_{\mathrm{FWHM}}`$ is the full-width at half-maximum of the microlensing light curve. This is given by $`t_{\mathrm{FWHM}}=3.5t_\mathrm{E}/A_{\mathrm{max}}`$ and is typically in the range 10–30 hours. We assumed 300 measurements were made in this time interval with accuracy two times worse than the photon statistical limit achievable with a 1-m class telescope, consistent with accuracies attained using the difference imaging technique. Full details of these simulations are reported elsewhere (Rattenbury et al 2002, submitted to MNRAS). Simulated light curves corresponding to a range of planetary positions were generated and the quantity $`\mathrm{\Delta }\chi ^2=\chi _{\mathrm{single}}^2300`$ was calculated for each light curve. Here, $`\chi _{\mathrm{single}}^2`$ denotes the value derived by fitting the light curve for a single, planet-less lens to the simulated data for a lens with a planet. We adopted a detection criterion of $`\mathrm{\Delta }\chi ^2>60`$ corresponding to a $`<`$1% probability for statistical fluctuations causing a planet-like signal.

Fig. 3 shows the zones of detectability for an Earth-mass planet under the above conditions. It is seen that Earth-mass planets can be detected with high efficiency if their orbital radii projected onto the lens plane lies within the range 0.7–1.5$`r_\mathrm{E}`$ or $``$1.3–2.9 AU in high magnification events that are monitored intensely during $`t_{\mathrm{FWHM}}`$. For larger planets heavier than Neptune, the zones of detectability extend far beyond the region around the Einstein ring especially when $`A_{\mathrm{max}}100`$. Jupiter-mass planets are detectable almost anywhere in these events.

During the 2002 southern winter viewing season, we expect to detect around 10 high magnification microlensing events in real-time. We propose that these events be monitored with a sampling as dense as possible during the the time interval given by $`t_{\mathrm{FWHM}}`$ for the event. Given the rapid nature of many of these events, follow-up observations will need to commence within 0–24 hours following the event alert. For events with $`t_{\mathrm{FWHM}}24`$ hours, a world-wide network of 1-m class telescopes is required to carry out the peak measurements. For more rapid events, fewer telescopes are required, but a commensurate increase in telescope aperture up to 2-m is required to achieve the sensitivity shown in Fig. 3. Difference imaging should be applied to all observations so as to achieve photometry with uncertainties approaching that given by the photon noise and free from systematic effects caused by blending of the source star with neighbouring stars. Looking further ahead, continuous monitoring of the peaks of high magnification events could be carried by follow-up observations from facilities in the Antarctic or in space.

Follow-up measurements using ground or space based facilities with deep imaging capability would be required for all events in order to determine to determine accurately the baseline intensity of the source star. This information is required in order to determine the value of $`A_{\mathrm{max}}`$ accurately and hence the planet:star mass ratio for any detected planet. In the case of non-detections, an accurate value of $`A_{\mathrm{max}}`$ is required for the calculation of planetary exclusion regions. We note that these deep follow-up observations are not time-critical and can be performed at any time after the events. Further follow-up measurements a few years later with the Next Generation Space Telescope should enable the lens star in any event to be observed directly as it begins to diverge from the source star. This would enable the absolute value of the mass of a planet to be determined, and also its absolute instantaneous projected radius at the time of the microlensing event.

## 4 Summary

We have demonstrated that high magnification microlensing events can be readily detected in microlensing surveys with a strategy of high frequency sampling of survey fields and real-time difference imaging analysis. We have presented 10 events with $`A_{\mathrm{max}}40`$ that were detected by the MOA microlensing survey during 2001. The purpose of this letter is to bring this capability to the attention of the microlensing community. We encourage intensive follow-up observations of future high magnification events detected by MOA and similar experiments.

Microlensing probes distant planets and is thus complementary to the radial velocity and transit techniques. Ultimately, a future space-based mission such as GEST \[Bennett et al. 2002\], dedicated to observing thousands of microlensing events will be required in order to detect substantial numbers of low mass planets (including those with masses as low as that of Mars). In the meantime however, it would appear that intensive observations of high magnification events would provide the best chances of detecting extra-solar planets in microlensing events observed from the ground.

## Acknowledgments

We thank P. Dobcsanyi, A. Gilmore, M. Honda, J.Jugaku, T. Nakamura, G. Nankivell, N. Rumsey, H. Sato, M. Sekiguchi, and Y. Watase for assistance, and the Marsden Fund of New Zealand and the Ministry of Education, Science, Sports and Culture of Japan for financial support.


text:
# Numerical Convergence of Hydrodynamical SPH Simulations of Cooling Clusters

## 1. Introduction

Galaxy clusters are the largest virialized structures known in the universe and also bright X-ray sources. Useful cosmological information can be obtained from the statistical properties of the ensemble of X-ray clusters. X-ray observations of cluster number counts, the X-ray temperature function (Henry & Arnoud 1991; Edge et al. 1990; Henry 1997) and the X-ray luminosity function (Rosati et al. 1998; Ebeling et al. 1998) are powerful probes to constrain the values of the cosmological parameters $`\mathrm{\Omega }_0`$ and $`\sigma _8`$ (Henry & Arnoud 1991; White, Efstathiou, & Frenk 1993; Bahcall & Fan 1998; Eke, Cole, & Frenk 1996; Kitayama, Sasaki, & Suto 1998)

Hydrodynamical simulations have been widely used to predict for different theoretical frameworks the time evolution of the gas and temperature distributions. The numerical methods used are either Eulerian (Cen & Ostriker 1994; Anninos & Norman 1996; Bryan & Norman 1998; Kang et al. 1994; Bryan et al. 1994, Cen 1997; Cen et al. 1995) with a fixed or adaptative grid, or Lagrangian (Evrard 1988; 1990; Thomas & Couchman 1992; Navarro, Frenk, & White 1995; Eke, Navarro, & Frenk 1998; Katz & White 1993; Yoshikawa, Itoh, & Suto 1998; Valdarnini, Ghizzardi, & Bonometto 1999) In these simulations the gas component is treated as a single adiabatic fluid, without taking into account the effects of radiative cooling, and the physical processes which can be modelled are merging, substructure formation, shocks and compressional heating of the gas. A comparison between different numerical simulations shows that they are successful in reproducing the gross features of the cluster properties (Frenk et al. 1999). With increasing availability of computational power, numerical hydrodynamical simulations have been attempted to model the effects of radiative cooling on the gas in the formation and evolution of cluster galaxies (Katz & White 1993; Suginohara & Ostriker 1998; Anninos & Norman 1996; Yoshikawa, Jing, & Suto 2000; Pearce et al. 2000; Lewis et al. 2000) The numerical problems posed by the inclusion of gas cooling are challenging, mainly because the required increased spatial resolution also requires that one keeps under control two-body heating mechanisms. Previous simulations have produced some conflicting results (Yoshikawa, Jing, & Suto 2000; Pearce et al. 2000), thus the question of the minimum resolution in this kind of simulations is still to be fully settled.

In this contribution I will show the preliminary results that have been obtained from a series of hydrodynamical SPH simulations of galaxy clusters. The simulations have different numerical resolution and include the effects on the gas component of radiative losses, star formation and energy feedback from SN. Final profiles are compared in order to assess the effects of numerical resolution, or different star formation prescriptions, on the cluster X-ray variables.

## 2. Simulations

In a previous paper (Valdarnini et al. 1999, hereafter VGB) a large set of hydrodynamical simulations was used to study the global X-ray cluster morphology and its evolution. The simulations were run using a TREESPH code with no gas cooling or heating. In order to assess the numerical reliability of the numerical integrations, four different clusters were selected as a representative sample of all the simulation clusters. For this cluster sample a large set of different integrations was performed by varying two numerical input parameters: the number of particles and the softening parameters. I refer to VGB for a detailed description of the simulations. In Valdarnini (2001) I have used the same cluster sample to study the effects of including in the simulations additional physics such as gas cooling and star formation. I will report here the results for the cluster with label $`\mathrm{\Lambda }`$CDM$`00`$ in VGB. This is the most massive cluster ($`M_v1.510^{15}M_{}`$) extracted from a cosmological $`\mathrm{\Lambda }`$CDM N-body simulation with size $`L=200h^1Mpc`$, matter density $`\mathrm{\Omega }_m=0.3`$ and Hubble constant $`H_0=70Kmsec^1Mpc^1`$. In order to check the effects of radiative cooling for this cluster a set of TREESPH simulations was performed, with initial conditions provided by the cosmological simulation, and different values of the gas softening parameter $`\epsilon _g`$ and number of gas particles $`N_g`$. For this cluster Table 1 reports the values of $`\epsilon _g`$ and $`N_g`$ for the simulation runs. The generic simulation has cluster index cl$`00j`$, with $`j=00,01,\mathrm{},05`$. The cluster cl$`0000`$ is the reference case without cooling. The effects of radiative cooling are modelled in these simulations by adding to the SPH thermal energy equation an energy-sink term. The total cooling function includes contribution from recombination and collisional excitation, bremsstrahlung and inverse Compton cooling.

Allowing the gas to cool radiatively will produce dense clumps of gas at low temperatures ($`10^4{}_{}{}^{}K`$). In these regions the gas will be thermally unstable and will likely meet the physical conditions to form stars. In TREESPH simulations star formation (SF) processes have been implemented using different approaches. According to Katz, Weinberg, & Hernquist (1996, hereafter KWH) a gas particle is in a star forming region if the flow is convergent and the local sound crossing time is larger than the dynamical time (i.e. the fluid is Jeans unstable). In a simplified version Navarro & White (1993 , hereafter NW) assume as a sufficient condition that a gas particle must be in a convergent flow and its density exceeds a threshold ($`\rho _g>\rho _{g,c}=710^{26}grcm^3`$).

If a gas particle meets these criteria then it is selected as an eligible particle to form stars. The local star formation rate (SFR) obeys the equation

$$d\rho _g/dt=c_{}\rho _g/\tau _g=d\rho _{}/dt,$$
(1)
where $`\rho _g`$ is the gas density, $`\rho _{}`$ is the star density , $`c_{}`$ is a characteristic dimensionless efficiency parameter, $`\tau _g`$ is the local collapse time and it is the maximum of the local cooling time $`\tau _c`$ and the dynamical time $`\tau _d`$. Gas particles with $`T\stackrel{<}{}10^4{}_{}{}^{}K`$ have long cooling times and $`\tau _g=\tau _d`$. At each step the probability that a gas particle will form stars in a time step $`\mathrm{\Delta }t`$ is compared with a uniform random number. If the test is successfull then a mass fraction $`\epsilon _{}`$ of the gas mass is converted into a new collisionless particle. This star particle has the position, velocity and gravitational softening of the original gas particle. Typical assumed values are $`\epsilon _{}=1/3`$ and $`c_{}=0.1`$ (KWH). For the local SFR, NW adopt Eq. 1 with $`c_{}=1`$, $`\tau _g=\tau _d`$ and $`\epsilon _{}=1/2`$ when a gas particle can convert part of its mass into a star particle.

The numerical tests have been performed following the NW prescriptions for selecting gas particles which can form stars. Once a star particle is created it can release energy into the surrounding gas through supernova (SN) explosions. All the stars with mass above $`8M_{}`$ end as a SN, leaving a $`1.4M_{}`$ remnant. The SN energy ($`10^{51}erg`$) is released gradually into the gas according to the lifetime of stars of different masses.

## 3. Results

The final radial gas density and temperature profiles are shown in Fig. 1. As can be inferred from Table 1, the numerical strategy has been to perform runs of increasing spatial resolution in order to resolve the core radius of the gas density profile. The most important result is that the inclusion of a star formation prescription has been effective in removing the cold gas particles ($`T\stackrel{<}{}10^4{}_{}{}^{}K`$) from the cluster center. In all the simulations the gas density profiles have a well-defined core radius, with size $`r_c50100Kpc`$, approximately $`0.05`$ of the virial radius $`r_v`$ ($`2.9Mpc`$). The shapes of the temperature profiles show that convergence is achieved for $`N_g\stackrel{>}{}20,000`$. All the central values for the gas temperatures at $`r=10Kpc`$ are within a factor $`1.5`$ for the highest resolution simulation runs. The profiles increase inwards from the virial radius up to $`100200Kpc`$. Thereafter the profiles stay almost flat, or with a modest decline in $`T(r)`$ towards $`r=0`$. There is not a strong drop of the temperature in the very central region. Final X-ray luminosities are quite stable ($`L_X210^{45}ergsec^1`$) versus the simulation numerical resolution. Compared to the non-radiative run the luminosities increase by a factor $`2`$.

A different question is the sensitivity of the estimated cluster properties to the numerical method used to describe star formation in the hydrodynamical simulations. To this end three simulations were performed with the same numerical parameters as for the cl$`0005`$ run, but using different SF methods or parameters. The NW simulation (I) is the standard case showed in Fig. 1. In the other two runs (II and III) the KWH prescription is adopted for converting gas particles into stars, but with different values of the star formation efficiency parameter $`c_{}`$$`0.1`$ and $`1.`$ (see Fig. 2). The most important differences are found between the KWH simulations with different $`c_{}`$. The differences are dramatic in the final X-ray luminosities, which differ by a factor $`40`$. The source of this discrepancy relies in the different gas density profiles, which have substantial differences in the cluster core regions for $`r\stackrel{<}{}100Kpc`$. These differences are localized at the cluster center, beyond $`r100Kpc`$ all the profiles converge, as it is shown in the plots of Fig. 2. The temperature profiles have a peak value of $`10^8{}_{}{}^{}K`$ at $`r100Kpc`$ and thereafter decline outward by a factor $`2`$ out to $`r_v`$. Below $`100Kpc`$ the profiles instead show large differences. Compared to the NW run the simulation with $`c_{}=0.1`$ has gas temperatures which decline inwards by a factor $`10`$ from $`100Kpc`$ down to $`r10Kpc`$. These radial decays follow because of the less efficient conversion of the cooled gas into stars compared to the NW run. There is a remarkable agreement for the ratio of cluster mass locked into stars to the gas mass, which is $`10\%`$ at $`r_v`$ for all the runs considered. In order to assess the effects of numerical resolution upon final results simulations I, II and III have been run again but with a number of particles increased by a factor $`3`$ (cl$`0005`$H). These simulations will be referred as IH, IIH and IIIH, respectively. The simulation results are shown in Fig. 2. For simulations IH there are not appreciable differences in the radial profiles. This confirms the previous results, that is the NW runs have reached numerical convergence in the physical variables for the numerical parameters of the cl$`0005`$ simulation of Table 1. The profiles of simulation IIH are instead different from those of run II at $`r\stackrel{<}{}50Kpc`$. The strong drop in $`T(r)`$ has been removed and the gas density profile is much closer to the NW one. Simulation IIIH yields final profiles very similar to the ones of the parent simulation. The bottom left panel of Fig. 2 shows that high resolution runs have final X-ray luminosities which can differ within a factor $`2`$ from the parent simulations.

To summarize, the above results demonstrate that simulations I and III have an adequate numerical resolution to reliably predict X-ray cluster properties, such as the X-ray luminosity. For simulation II ( KWH with $`c_{}=0.1`$ ) there are large differences at the cluster core between the final profiles when the numerical resolution is increased. For the simulation runs I and III $`L_X10^{45}ergsec^1`$, while $`L_X410^{46}ergsec^1`$ for the run with $`c_{}=0.1`$. Thus simulations I and III are consistent at the $`1\sigma `$ level with the $`L_{bX}T_X`$ relation derived from a sample of cooling flow clusters (Allen & Fabian 1998). For simulation II the predicted cluster temperature is outside the $`2\sigma `$ limits.

## References

Allen, S.W., & Fabian, A.C. 1998, MNRAS, 297, L57

Anninos, P., & Norman, M.L. 1996, ApJ, 459, 12

Bahcall, N.A., & Fan, X. 1998, ApJ, 504, 1

Bryan, G.L., Cen, R., Norman, M.L., Ostriker, J.P., & Stone, J.M. 1994, ApJ, 428, 405

Bryan, G.R., & Norman, M.L. 1998, ApJ, 495, 80

Cen, R. 1997, ApJ, 485, 39

Cen R., & Ostriker J., 1994, ApJ, 429, 4

Cen, R., Kang, H., Ostriker, J.P., & Ryu, D. 1995, ApJ 451, 436

Ebeling, H. et al. 1998, MNRAS, 301, 881

Edge, A.C., Stewart, G.C., Fabian, A.C., & Arnaud, K.A. 1990, MNRAS, 245, 559

Eke, V.R., Cole, S., & Frenk, C.S. 1996, MNRAS, 282, 263

Eke, V.R., Navarro, J.F., & Frenk, C.S. 1998, ApJ, 503, 569

Evrard, A.E. 1988, MNRAS, 235, 911

Evrard, A.E. 1990, ApJ, 363, 349

Frenk, C.S. et al. 1999, ApJ, 525, 554

Henry, J.P., & Arnoud, K.A. 1991, ApJ, 372, 410

Henry, J.P. 1997, ApJ, 489, L1

Kang, H. , Cen, R. , Ostriker, J.P., & Ryu, D. 1994, ApJ, 428, 1

Katz, N., & White, S.D.M. 1993, ApJ, 412, 455

Katz, N., Weinberg, D.H., & Hernquist, L. 1996, ApJS, 105, 19

Kitayama, T., Sasaki, S., & Suto, Y. 1998, PASJ, 50, 1

Lewis, G.F., Babul, A., Katz, N., Quinn, T., Hernquist, L., & Weinberg, D.H. 2000, ApJ 536, 623

Navarro, J., & White, S.D.M. 1993, MNRAS, 265, 271 (NW)

Navarro, J., Frenk C.S., & White, S.D.M. 1995, MNRAS, 275, 720

Pearce, F.R., Thomas, P.A., Couchman, H.M.P., & Edge, A.C. 2000, MNRAS 317, 1029

Rosati, P., Della Ceca, R., Norman, C., & Giacconi, R. 1998, ApJ, 492, L21

Suginohara, T., & Ostriker, J.P. 1998, ApJ, 507, 16

Thomas, P.A., & Couchman, H.M.P. 1992, MNRAS, 257, 11

Valdarnini, R., Ghizzardi, S., & Bonometto, S. 1999, New Astr, 4, 71 (VGB)

Valdarnini, R. 2001, in preparation

White, S.D.M., Efstathiou, G., & Frenk, C.S. 1993, MNRAS, 262, 1023

Yoshikawa, K., Itoh, M., & Suto, Y. 1998, PASJ, 50, 203

Yoshikawa, K., Jing, Y.P., & Suto, Y. 2000, ApJ, 535, 593 (YJS)