text:
# Interplay between wetting and phase behavior in binary polymer films and wedges: Monte Carlo simulations and mean field calculations

## 1 Introduction

Confining a binary mixture, one can profoundly alter its miscibility behavior. The phase behavior of $`AB`$ mixtures in pores, slits and films has attracted abiding interest from both theorists and experimentalists alike. We study the interplay between (pre)wetting and phase behavior by self-consistent field (SCF) theory and Monte Carlo simulations. Particularly, we focus on situations where surfaces attract different components of the mixture.

The qualitative features of the miscibility in confined geometry are rather universal and shared by polymer mixtures as well as small molecules. Symmetric binary polymer blends are, however, particularly well suited to study the interplay between wetting and miscibility: (i) the wetting transition temperature typically is much lower than the critical temperature, where demixing occurs in the bulk, (ii) Fluctuations can be controlled by the degree of interdigitation: the more extended the molecule is, the larger is the number of neighbors it interacts with, and the smaller the effect of fluctuations. Therefore SCF calculations provide an accurate description for many properties except for the ultimate vicinity of critical points. The spatial extension of the molecules also sets the length scale of enrichment layers and facilitates experimental investigations. Indeed, wetting transitions have been studied in recent experiments. (iii) The vapor pressure of polymer films is vanishingly small, hence effects of evaporation can be neglected. (iv) Polymers tend not to crystallize easily. Therefore, wetting phenomena might not be preempted by crystal phases. Likewise there is no roughening transition of the interface as it occurs in Ising-like models.

Using a coarse-grained polymer model for a $`AB`$ binary melt we locate the first order wetting transition, the phase diagram in a symmetric slit pore (symmetric film), the phase diagram in a thin film where the substrate favor the $`A`$-component of the mixture with the same strength as the top surface attracts the $`B`$-component (antisymmetric film). Then we discuss the phase behavior in a quadratic pore where two neighboring surfaces favor the $`A`$-component and the other two neighboring surfaces favor the $`B`$-component (antisymmetric double wedge). We conclude with an outlook.

## 2 Model and techniques

We consider a binary polymer blend. Both species – $`A`$ and $`B`$ – contain the same number $`N`$ of monomers and have the same spatial extension $`R_e`$. They are confined into a thin film; the bottom substrate (W) might be a silicon wafer, while the other surface might be the interface to the vapor (vacuum, V). Depending on the ratio between the interface tension $`\gamma _{AB}`$ between the segregated bulk phases and the surface tension $`\gamma _{AV}`$, $`\gamma _{BV}`$ of the components and the vapor, the upper surface might be rough. The qualitative behavior is illustrated in Fig.1.

If the $`AB`$ interface tension is comparable to the liquid/vapor tension, it “drags” the film surface towards the substrate so as to reduce the length of the $`AB`$ interface. If the liquid/vapor tension exceeds the $`AB`$ interface tension by about two orders of magnitude, however, the surface is almost flat and the situation is equivalent to a binary mixture between two hard walls a distance $`\mathrm{\Delta }`$ apart. In the following we shall restrict ourselves to this limit $`\gamma _{AB}\gamma _{AV}`$ or $`\gamma _{BV}`$.

In the Monte Carlo simulations we use a computationally efficient, coarse-grained lattice model. The bond fluctuation model retains the universal features of polymers – connectivity, excluded volume of segments and a thermal interaction which leads to phase separation – but ignores details of chemical structure. Effective monomers prevent the corners of a unit cell of a 3D cubic lattice from double occupancy. We use chain length $`N=32`$ and $`R_e17u`$. Monomers along a chain are connected via bond vectors of length $`2,\sqrt{5},\sqrt{6},3`$ or $`\sqrt{10}`$ in units of the lattice spacing $`u`$. Different monomers repel each other by a square well potential of depth $`ϵ`$ which comprises the nearest 54 neighbors, like monomers attract each other. The strength of the repulsion is proportional to the Flory–Huggins parameter $`\chi =5.3ϵ/k_BT`$. Surfaces are structureless and impenetrable. They act on monomers in the two nearest layers ($`d_{\mathrm{wall}}=2`$) with strength $`ϵ_{\mathrm{wall}}`$.

In the SCF calculations we model the polymers as Gaussian chains. The repulsion between different species is quantified by the Flory-Huggins parameter $`\chi `$. Short-ranged interactions of strength $`\mathrm{\Lambda }_1`$ and $`\mathrm{\Lambda }_2`$ attract (repel) the $`A`$ ($`B`$) component in the vicinity of the surfaces. The total density profile of the film is imposed. It smoothly decays to zero at the surfaces in a boundary region of width $`0.15R_e`$. The blend is assumed to be incompressible. This standard Gaussian chain model is solved within mean field approximation.

## 3 Wetting transition

To accurately locate the wetting transition and calculate the contact angle of macroscopic $`A`$-drops we use Young’s equation $`\gamma _{AB}\mathrm{cos}\mathrm{\Theta }=\gamma _{WB}\gamma _{WA}`$. Computationally, this techniques has distinct advantages for locating first order wetting transitions: (i) The interface free energy $`\gamma _{AB}`$ and the difference $`\mathrm{\Delta }\gamma =\gamma _{WB}\gamma _{WA}`$ can be measured accurately in separate simulations thereby avoiding the need for huge simulation cells to simulate a thick $`A`$-layer at the surface in equilibrium with a $`B`$-rich bulk. (ii) Unlike observing the dependence of the thickness of the $`A`$-layer on temperature or monomer-surface attraction, one directly measures free energies. Therefore, we do accurately locate the transition, while the instability of the $`A`$-rich layer is located between the transition and the mean-field wetting spinodal. (iii) By virtue of the $`AB`$ symmetry, the difference $`\mathrm{\Delta }\gamma `$ can also be rewritten as the difference $`\mathrm{\Delta }\gamma =\gamma _{WB}\gamma _{WB}`$ of surface tensions of a wall that attract the $`A`$-component and a wall that attracts the $`B`$-component. This free energy difference can be measured by thermodynamic integration or expanded ensemble methods.

The results for our model are presented in Fig.2. From the crossing of $`\gamma _{AB}(ϵ)`$ and $`\mathrm{\Delta }\gamma (ϵ)`$ we locate the wetting transition. The fact that curves intersect under a finite angle indicates that the wetting transitions are of first order. As we reduce the monomer-surface attraction the wetting transition shifts to higher temperatures $`k_BT/ϵ`$ and become weaker. For all transitions studied, however, the wetting transition is of first order. This is also corroborated by SCF calculations, where we find first order wetting transitions for $`T/T_c<0.98`$.

If the wetting transition is of first order, then there will be only a small $`A`$-rich layer in the non-wet state. By virtue of the structural symmetry of the molecules, they loose the same amount of entropy as they pack against the surface. The surface free energy difference $`\mathrm{\Delta }\gamma `$ is mainly enthalpic. If we assume that the wetting transition is strongly first order, we can neglect the microscopic enrichment layer at the surface $`\mathrm{\Delta }\gamma =2ϵ_{\mathrm{wall}}d_{\mathrm{wall}}\mathrm{\Phi }`$, where $`d_{\mathrm{wall}}=2`$ denotes the range of the monomer-surface interaction and $`\mathrm{\Phi }=1/16`$ the monomer number density. Using the expression for the interface tension $`\gamma _{AB}=\mathrm{\Phi }b\sqrt{\chi /6}`$ ($`b=3.05`$: statistical segment length) in the strong segregation limit, we obtain: $`\chi _{\mathrm{wet}}=24\left(\frac{ϵ_{\mathrm{wet}}d_{\mathrm{wet}}}{bk_BT}\right)^2`$

This is in marked contrast to the value of the Flory–Huggins parameter at the unmixing transition in the bulk, $`\chi _c=2/N1/T_c`$. As both the interface tension $`\gamma _{AB}`$ and the difference in surface tension $`\mathrm{\Delta }\gamma `$ are chain length independent, so is the wetting transition temperature. The critical temperature $`T_c`$ of phase separation, however, increases linearly with chain length $`N`$. Therefore, critical phenomena associated with the bulk unmixing and wetting phenomena are well separated.

## 4 Thin films

### 4.1 Capillary condensation and interface localization/delocalization

If the mixture is confined into a film the surface interactions modify the phase behavior. As wetting is associated with the growth of an infinitely large enrichment layer, it is rounded-off in a thin film. If the wetting transition is of first order, there will be a pre-wetting transition: a coexistence between a thin and a thick (but microscopic) enrichment layer at a chemical potential which differs from the value at coexistence in the bulk. As pre-wetting transitions involve only enrichment layers of finite thickness, they might give rise to transitions in thin films.

First we regard a film with symmetric surfaces, i.e. both surfaces attract the $`A`$ component. The phase diagram as obtained from the simulations is presented in Fig. 3. Compared to the phase behavior in the bulk, the critical point is shifted to lower temperatures and larger composition of the species attracted by the surfaces. Moreover, the binodal in the vicinity of the critical points exhibit 2D Ising critical behavior in contrast to the 3D Ising behavior of the bulk unmixing transition.

Note the pronounced distortion of the $`B`$-rich binodal in the vicinity of the wetting transition. In the $`B`$-rich phase there are $`A`$-rich layers at the surfaces and the $`B`$ component prevails in the middle of the film. In the vicinity of the wetting transition the thickness of the $`A`$-enrichment layers grows as we increase the temperature. If we increased the film thickness this distortion would evolve into an additional two phase region, corresponding to a $`B`$-rich phase with thin and thick $`A`$-layers at the surface. This two phase region would correspond to the pre-wetting coexistence and it would join the $`B`$-rich binodal in a triple point.

The phase diagram of an antisymmetric film is also presented in Fig.3. In this case one surface attracts the $`A`$-component with exactly the same strength than the other surface the $`B`$-component. The phase diagram contains two critical points and a triple line. Around the critical temperature of the bulk, enrichment layers gradually form at the surfaces and stabilize an $`AB`$ interface that runs parallel to the surfaces. At the interface localization/delocalization transition this $`AB`$ interface becomes bound to one of the surfaces. In the case of a first order interface localization/delocalization transition this corresponds to a triple point of the phase diagram: an $`A`$-rich phase , a $`B`$-rich phase and a phase with symmetric composition coexist.

The behavior can be analyzed qualitatively by looking at the interface potential $`g(l)`$ which describes the interaction between an $`AB`$ interface and a single surface. If the film is thick enough, the interface potential can be constructed as a superposition of the interface potentials emerging from each surface. The qualitative behavior in the vicinity of a first order wetting transition is depicted in Fig.4(a). Using a double-tangent construction we can construct the phase behavior in a thin film. At low temperatures there coexist an $`A`$-rich phase and a $`B`$-rich phase, in which the $`AB`$ interface is localized at the surface. Upon increasing the temperature, one encounters the triple point. This triple point is the thin film analogon of the first order wetting transition. As the film thickness increases the triple temperature converges towards the wetting transition temperature of the semi-infinite system. Above the triple temperature there are two phase coexistence regions, which correspond to thin and thick enrichment layers at the surfaces. This is the analogon of the pre-wetting transition in a thin film.

### 4.2 The tricritical interface localization/delocalization transition

If we reduce the film thickness, the interactions emerging from each surface interfere. The phenomenological considerations explain that this leads to a second order interface localization/delocalization transition at small film thicknesses. Both regimes are separated by a tricritical transition.

The scaled distribution functions prove convenient to locate the tricritical thickness accurately. To this end we have adjusted the temperature such that the central peak of the probability distribution of the order parameter $`m\varphi _A\varphi _B`$ is a factor 1.2 higher than the outer peaks. This corresponds to the behavior of the universal distribution of the 2D tricritical universality class. The results for various film thicknesses $`\mathrm{\Delta }`$ (in units of the lattice spacing) are presented in Fig.5. For $`\mathrm{\Delta }<\mathrm{\Delta }_{\mathrm{tri}}`$ the valleys between the three peaks are too shallow (cf. Fig.5), while they are too deep for $`\mathrm{\Delta }>\mathrm{\Delta }_{\mathrm{tri}}`$. In the latter case the transition is of first order and our estimate tends towards the triple temperature. At $`\mathrm{\Delta }_{\mathrm{tri}}14=0.89R_e`$ the distribution of our simulations is similar to the universal 2DT distribution, and this has been confirmed for larger lateral system sizes.

### 4.3 Crossover from capillary condensation to interface localization/delocalization

Realizing strictly (anti)symmetric surface interactions is often difficult in experiments. Varying the surface interaction $`\mathrm{\Lambda }_2N`$ of the top surface from attracting the $`A`$-component to attracting the $`B`$-component (while the bottom surface always attracts the $`A`$-component with fixed strength $`\mathrm{\Lambda }_1N`$) in the SCF calculations, we study the crossover from capillary condensation for symmetric surfaces to interface localization/delocalization. The dependence of the phase diagram on the surface interactions within the SCF calculations is presented in Fig.6. For symmetric surfaces (capillary condensation) the critical point is shifted towards lower temperatures similar to the simulation result. Of course, the binodals are parabolic in mean field theory independent from dimensionality. The coexisting phases have almost uniform composition across the film and differ in their composition. As we reduce the preference of the top surface for species $`B`$, the critical point and the critical composition tend towards their bulk values ($`\varphi =0.5`$, $`1/\chi N=0.5`$), i.e. the critical temperature increases and the critical composition becomes more symmetric. The coexistence curve in the $`1/\chi N`$-$`\mathrm{\Delta }\mu `$ plane approaches the symmetry axis. Upon making the top surface attracting the other component $`B`$, we gradually change the character of the phase transition towards an interface localization/delocalization transition. The critical temperature passes through a maximum and the critical composition through a minimum, respectively. For $`\mathrm{\Lambda }_2N<0`$ (surface attracting the $`B`$-component) there are enrichment layers of the $`A`$-component at the bottom and the $`B`$-component at the top, and the two coexisting phases differ in the location of the $`AB`$ interface which runs parallel to the surfaces. As the preferential interaction of the top surface increases, the critical temperature decreases and the critical composition becomes richer in $`A`$. When the coexistence curve intersects the pre-wetting line of the bottom surface at $`\mathrm{\Delta }\mu <0`$, a triple point forms at which an $`A`$-rich phase and two $`B`$-rich phases with a thin and a thick $`A`$-enrichment layer coexist. When the bottom surface attracts $`A`$ with exactly the same strength as the top surface $`B`$ (antisymmetric surfaces), the phase diagram becomes symmetric.

For symmetric surface fields the critical point occurs close to the bulk critical point (and converges towards it in the limit of infinite film thickness) while the critical points in antisymmetric films are associated with the wetting transition and converge towards the pre-wetting critical temperature of the semi-infinite system (if the wetting transition is of first order) for $`\delta \mathrm{}`$. In both cases, however, critical points belong to the 2D Ising universality class.

## 5 Interface localization/delocalization in an antisymmetric double wedge

### 5.1 Background

In the following we consider wetting (or rather filling) in a wedge geometry. Macroscopic considerations show that the wedge will be filled with liquid when the contact angle $`\mathrm{\Theta }`$ on a planar substrate equals the opening angle $`\alpha `$. Intriguingly, Parry and co-workers predict that the filling of a wedge is related to the strong fluctuation regime of critical wetting and that critical filling may even occur if the concomitant wetting transition of the planar surface is of first order. Specifically, they predicted the distance $`l_0`$ of the $`AB`$ interface from the bottom of a wedge to diverge as $`l_0(T_fT)^{\beta _s}`$ with $`\beta _s=1/4`$. Correlations along the wedge and in the other two directions are characterized by diverging correlation lengths $`\xi _y(T_fT)^{\nu _y}`$ and $`\xi _x\xi _{}(T_fT)^\nu _{}`$ with exponents $`\nu _y=3/4`$ and $`\nu _{}=1/4`$, respectively.

In the following we study a wedge with opening angle $`\alpha =\pi /4`$ of the wedge (c.f. Fig.7). Similar to the study of wetting an antisymmetric geometry is advantageous. Therefore we stack two wedges which attract different components on top of each other. This antisymmetric double wedge is a pore with quadratic cross-section of size $`L\times L`$. Let $`L_y`$ denote the length of the wedge (c.f. Fig.7). (i) If we used identical surface fields on all four free surfaces the analog of capillary condensation would occur in a wedge, i.e., phase coexistence would be shifted away from the bulk coexistence curve and the wetting layers would be only metastable (with respect to “wedge condensation”). (ii) As the wetting layer grew on all four surfaces in the case of symmetric boundaries, we would need larger system sizes to reduce interactions between the wetting layers across the wedge.

The phase behavior in such an antisymmetric double wedge geometry has been studied recently in the framework of an Ising model. When the wetting transition of the planar substrate was of first order, the wedge filling was also found to be of first order. When the wetting transition was of second order an unconventional scaling behavior was observed which is characterized by critical exponents $`\alpha =3/4`$, $`\beta =0`$, and $`\gamma =5/4`$. Those critical exponents can be related (see below) to the exponents of critical wedge filling, and the simulations of the Ising model confirm the predictions of Parry and co-workers.

In the following we corroborate these findings in the framework of the Ising model by our polymer simulations. Moreover, we present evidence for the unconventional second order transition in an antisymmetric double wedge even though the wetting transition on a planar substrate is of first order.

We present preliminary simulation data for two temperatures: $`ϵ/k_BT=0.025`$ ($`T/T_c=0.58`$) and $`ϵ/k_BT=0.05`$ ($`T/T_c=0.29`$). At both temperatures the wetting transitions, that occur at appropriate attractive strength $`ϵ_{\mathrm{wall}}`$ of planar surfaces, are of first order (c.f. Fig.2). In the former case it is weak first order wetting transition, in the latter case it is a strong first order transition.

### 5.2 First order transition in an antisymmetric double wedge

At the lower temperature $`ϵ/k_BT=0.05`$, the behavior is similar to a first order interface localization/delocalization transition. We consider here only the case $`\mathrm{\Delta }\mu =0`$ where phase coexistence in the bulk occurs. This excludes the rather interesting interplay between pre-wetting and pre-filling behavior studied in Ref.. At large surface interaction $`ϵ_{\mathrm{wall}}>ϵ_{\mathrm{wall}}^{\mathrm{trip},\mathrm{wedge}}`$ there runs an $`AB`$ interface along the diagonal which divides the two double wedges. This corresponds to the delocalized state. Upon decreasing $`ϵ_{\mathrm{wall}}`$ (or decreasing the temperature) the $`AB`$ interface becomes localizes in one of the wedges. In this case the composition of the double wedge is either $`A`$-rich or $`B`$-rich and we define as order parameter $`m\varphi _A\varphi _B`$. The two situations are separated by a triple point $`ϵ_{\mathrm{wall}}^{\mathrm{trip},\mathrm{wedge}}`$ at which the interface can be localized in either of the wedges or delocalized on the diagonal. The trimodal probability distribution in the vicinity of the tricritical point is presented in Fig.8. In analogy to the case of antisymmetric films we expect this triple point in a double wedge to correspond to a first order filling transition. In the inset we show the cumulant $`m^2/|m|^2`$. If the transition was of second order, these cumulants would depend monotonically on $`ϵ_{\mathrm{wall}}`$ and would exhibit a common intersection point. This is not at all what we observe, and we conclude that the interface localization/delocalization transition in the double wedge is of first order at the lower temperature $`ϵ/k_BT=0.05`$.

### 5.3 Critical behavior in an antisymmetric double wedge

Even though the wetting transition on a planar surface at $`ϵ/k_BT=0.025`$ is of first order, the behavior at the interface localization/delocalization transition in an antisymmetric double wedge at high temperature differs from the transition at low temperature. In the inset of Fig.9 (a) we present the dependence of the cumulant on the surface interaction strength for various system sizes. The cumulants depend monotonically on $`ϵ_{\mathrm{wall}}`$ and exhibits a common intersection point around $`ϵ_{\mathrm{wall}}^{\mathrm{crit}}0.035`$. In panel (b) we show the probability distribution of the composition $`\varphi `$ at this intersection point: the distribution is bimodal and the two largest system sizes collapse onto a master curve without any size dependent prefactor. Therefore we conclude that the interface localization/delocalization transition is of second order.<sup>1</sup><sup>1</sup>1The interface localization/delocalization transition might be of second order in a very thin antisymmetric film (c.f. Sec. 4.2) even if the wetting transition is of first order. Therefore still larger system sizes would be desirable to confirm this conclusion. We note however, that the thickness of the enrichment layer at the first order wetting transition of the planar substrate ($`ϵ/k_BT=0.0226`$, $`ϵ_{\mathrm{wall}}^{\mathrm{wet}}/k_BT=0.04`$) is only $`l_0445.2=L/\sqrt{2}`$. Therefore we believe that our conclusion is not affected by finite size effects.

Intriguingly there are also marked differences between this second order transition in an antisymmetric double wedge and the second order transition in a thin film which belongs to the 2D Ising universality class. In the latter case, only the distribution of the scaled order parameter $`L^{\beta /\nu }m`$, where $`\beta =1/8`$ and $`\nu =1`$ are the critical exponents of the order parameter and the correlation length in the 2D Ising universality class, exhibits data collapse for different system size. Moreover, we present in Fig.9(a) the dependence of the absolute value of the magnetization in the vicinity of the transition. Curves for different system sizes exhibit a common intersection point which agrees well with the intersection point of the cumulants. The analogous curves at a Ising-like transition do not exhibit a common intersection point but monotonously converge towards $`|m||TT_c|^\beta `$ for $`T<T_c`$ and $`|m|0`$ for $`TT_c`$ upon increasing the system size.

To relate the critical behavior of the antisymmetric double wedge to the predictions of Parry et al., we regard the distance $`l_0`$ of the $`AB`$ interface from the corner of one wedge. Similar to an antisymmetric film (c.f. Sec. 4.1), we assume that we can approximate the distribution in a double wedge by the superposition of the distributions of single wedges $`P_{\mathrm{wedge}}(l_0)`$ via $`P(l_0)P_{\mathrm{wedge}}(l_0)+P_{\mathrm{wedge}}(\sqrt{2}Ll_0)`$ If the two distributions $`P_{\mathrm{wedge}}(l_0)`$ and $`P_{\mathrm{wedge}}(\sqrt{2}Ll_0)`$ do not overlap, the $`AB`$ interface will be located in either of the two wedges and the order parameter will not vanish. If the two distributions overlap, the interface fluctuates around the diagonal and the order parameter will be zero. Right at the transition the two distributions begin to overlap:

$$l_0+\xi _{}\stackrel{!}{=}\sqrt{2}Ll_0\xi _{}\text{(interface localization/delocalization in double wedge)}$$
(1)
where $`l_0`$ denotes the mean height in a single wedge and $`\xi _{}`$ its fluctuations. Importantly, Parry’s prediction of $`\beta _0=\nu _{}`$ in wedges (and also corners) means that the height and its fluctuations are of the same order. They diverge as we approach the critical filling transition.

The height of the interface $`l_0`$ is related to the order parameter $`m`$ of the localization/delocalization transition. Therefore we expect the distribution of the order parameter also to be bimodal. As $`l_0\xi L`$ at the transition and the order parameter is a function of $`l_0/L`$ the distribution of the order parameter will exhibit two peaks whose positions and widths will not depend on the system size. This is exactly what we observe in Fig.9(b). Using this observation and the standard finite size scaling assumption at a second order phase transitions

$$P(m)L^{\beta /\nu _{}}\stackrel{~}{𝒫}(L^{\beta /\nu _{}}m,L/\xi _{},L_y/\xi _y)L^{\beta /\nu _{}}𝒫(L^{\beta /\nu _{}}m,L^{1/\nu _{}}t,\eta )$$
(2)
where $`\stackrel{~}{𝒫}`$ and $`𝒫`$ are scaling functions, $`t=(TT_f)/T_f`$ denotes the relative distance to the filling transition, and $`\eta L_y/L^{\nu _y/\nu _{}}=L_y/L^3`$ denotes the generalized aspect ratio, we conclude $`\beta /\nu _{}=0`$. Due to the anisotropy of the fluctuations of the interface along the wedge with correlation length $`\xi _y`$ and perpendicular to the wedge with correlation length $`\xi _{}`$ the generalized aspect ratio appears as a scaling variable. In our simulations we have chosen the system geometry such that $`\eta `$ remains approximately constant to ensure that finite size finite-size rounding in the direction along the wedge and the rounding in the two other directions set in simultaneously.<sup>2</sup><sup>2</sup>2If we kept the ratio $`L_y/L`$ constant $`\eta 0`$ and the system would exhibit a behavior characteristic of a corner. In the limit $`L`$ fixed but $`L_y`$ the wedge becomes quasi-onedimensional and there is no transition. Hence, the scaling of the probability distribution not only confirms $`\beta =0`$ but also $`\nu _y=3\nu _{}`$.

Knowing the probability distribution of the order parameter we can calculate all its moments:

$$m^k=_k(L^{1/\nu _{}}t,\eta )$$
(3)
where $`_k`$ are scaling functions. As a special case, we calculate the susceptibility: $`\chi =L^2L_ym^2/k_BTL^2L_y\stackrel{~}{}_2(L/\xi _{},L_y/\xi _y)\xi _{}^2\xi _yt^{2\nu _{}\nu _y}t^\gamma `$ with $`\gamma =2\nu _{}+\nu _y=5/4`$. Gratifyingly these values for the exponents comply with the anisotropic hyperscaling relation $`\gamma +2\beta =(d1)\nu _{}+\nu _y`$. Using thermodynamic scaling $`2\alpha =\gamma +2\beta `$ we infer the critical exponent $`\alpha =3/4`$ for the specific heat. Another consequence of the absence of any $`L`$-dependent prefactor in Eq.(3) is the common intersection of moments of the order parameter at the transition. Again this is an agreement with our observation in Fig,8(a). As this intersection involves only the lowest moment of the order parameter it yields an accurate estimate of the location of the critical interface localization/delocalization transition in an antisymmetric double wedge.

It is interesting to relate the observation of first and second order interface localization/delocalization transition in a double wedge to the shape of the interface potential. Parry et al.  predict that the filling transition is second order if the interface potential between an $`AB`$ interface and a planar surface does not exhibit a free energy barrier between the minimum close to the surface and the behavior at large distances, i.e. if a macroscopically thick film is not even metastable.

In Fig.10 we present the interface potential obtained from the probability distribution of the composition in a simulation of an antisymmetric film at $`ϵ/k_BT=0.025`$. In the vicinity of the wetting transition the interface potential exhibits a maximum between the minimum close to the surface and the value at large distances. This fact confirms that the wetting transition is of first order. At the smaller value of $`ϵ_{\mathrm{wall}}`$, however, there is no such maximum within the statistical uncertainty of the Monte Carlo data and, in agreement with Parry’s predictions, we observe a second order transition in the double wedge.

## 6 Summary

We have investigated the interplay between wetting and phase separation of incompressible binary mixtures confined in thin films and wedges. In our polymer model the wetting transition if of first order and we can accurately locate it via Young’s equation. The concomitant pre-wetting behavior modifies the phase boundaries in thin films. If both surfaces attract the same component, capillary condensation occurs and the critical point is close to the critical unmixing transition in the bulk. If one surfaces attracts the $`A`$-component but the other attracts the $`B`$-component an interface localization/delocalization transition occurs. In this case there are two critical points which correspond to the pre-wetting critical points at each surface. If the film thickness is very small, however, the interface localization/delocalization transition might be of second order even if the wetting transition is of first order. The critical points in a thin film are characterized by Ising critical behavior.

In analogy to the interface localization/delocalization in an antisymmetric film, we have studied the transition in an antisymmetric double wedge and we relate the phase behavior to the filling transition in a single wedge. Importantly we present evidence that the analog of critical filling in an antisymmetric double wedge geometry gives rise to unconventional critical behavior characterized by an order parameter exponent $`\beta =0`$ and strong anisotropic fluctuations. We can relate the critical exponents to the predictions of Parry et al. on critical filling. In agreement with those predictions the filling transition can be critical even though the wetting transition on a planar substrate is of first order. This is practically important because there is no experimental realization of critical wetting on a solid substrate. Our findings suggests the polymer blends might be promising candidates to explore the filling behavior experimentally.

Acknowledgment: It is a great pleasure to thank E.V. Albano, D.P. Landau and A. Milchev for fruitful collaborations and A.O. Parry for stimulating discussions. Financial support by the DFG under grants Mu 1674/1-1, Bi314/17-4, and DAAD/PROALAR 2000 as well as computer time at NIC Jülich and HLR Stuttgart are acknowledged.


format: text

text:
# Forecasting characteristic earthquakes in a minimalist model

## I  Introduction

The earthquake process in seismic faults is a very complex natural phenomenon that present geophysics, in spite of its considerable efforts, has not yet been able to put into a sound and satisfactory status. However, in the crucial field of earthquake prediction, recent years have witnessed significant advances. For recent thorough reviews dealing with this issue, see Keilis-Borok (2002); Keilis-Borok and Soloviev (2002), and references therein, in particular chapter four by Kossobokov and Shebalin. See also Lomnitz (1994). The introduction of new concepts coming from modern statistical physics seems to add some light and put some order into the intrinsic complexity of the lithosphere and its dynamics. Thus, for example, references to critical phenomena, dynamical systems, hierarchical systems, fractals, self-organized criticality and self-organized complexity are now found very frequently in geophysical literature (Turcotte, 2000; Sornette, 2000; Gabrielov et al., 1999; Gabrielov et al., 2000). Hopefully, this conceptual framework will prove its usefulness sooner better than later.

We have recently presented a simple statistical model of the cellular-automaton type which produces an earthquake spectrum similar to the characteristic earthquake behaviour of some seismic faults (Vázquez-Prada et al., 2002). The largest earthquakes on a fault or fault segment (the events that break its complete length) are usually termed characteristic (Schwartz and Coppersmith, 1984; Wesnousky, 1994; Dahmen et al., 1998). For this reason, in the minimalist model the event of maximum size is called the characteristic one. Our model is inspired by the concept of asperity, i.e., by the presence of a particularly strong element in the system which actually controls its relaxation. This model presents some notable properties, some of which will be reviewed in Section II. In Section III, an algebraic procedure for the exact calculation of the probability distribution of the time of return of the characteristic earthquake is presented. The purpose of this paper is to quantify the forecasting of the characteristic earthquake occurrence in this model, using seismicity functions, which are observable, but not stress functions (Ben-Zion et al., 2003), which are not. In Section IV, we construct an error diagram (Molchan, 1997; Newman and Turcotte, 2002) based on the time elapsed since the occurrence of the last characteristic event. This permits a first assessment of the degree of predictability. In Section V, we propose a general strategy of classification of the seismic cycles which, adequately exploited in this model, allows a refinement of the forecasts. Finally, in Section VI we present the conclusions.

## II  Some Properties of the Model

In the minimalist model (Vázquez-Prada et al., 2002), a one-dimensional vertical array of length $`N`$ is considered. The ordered levels of the array are labelled by an integer index $`i`$ that runs upwards from $`1`$ to $`N`$. This system performs two basic functions: it is loaded by receiving stress particles in its various levels and unloaded by emitting groups of particles through the first level $`i=1`$. These emissions that relax the system are called earthquakes.

These two functions (loading and unloading) proceed using the following four rules:

1. In each time unit, one particle arrives at the system.
2. All the positions in the array, from $`i=1`$ to $`i=N`$, have the same probability of receiving the new particle. When a position receives a particle we say that it is occupied.
3. If a new particle comes to a level which is already occupied, this particle has no effect on the system. Thus, a given position $`i`$ can only be either non-occupied when no particle has come to it, or occupied when one or more particles have come to it.
4. The level $`i=1`$ is special. When a particle goes to this first position a relaxation event occurs. Then, if all the successive levels from $`i=1`$ up to $`i=k`$ are occupied, and the position $`k+1`$ is empty, the effect of the relaxation (or earthquake) is to unload all the levels from $`i=1`$ up to $`i=k`$. Hence, the size of this relaxation is $`k`$, and the remaining levels $`i>k`$ maintain their occupancy intact.

Therefore, the size of the earthquakes in this model range from 1 up to $`N`$, being the event of $`k=N`$ the characteristic one. Note that the three first rules of this model are exactly those of the forest-fires model (Drossel and Schwabl, 1992). Our model has no parameter and, at a given time, the state of the system is specified by stating which of the $`(i>1)`$ $`N1`$ ordered levels are occupied. Each one of these possible occupation states corresponds to a stable configuration of the system, and therefore the total number of configurations is $`2^{(N1)}`$. These mentioned $`2^{(N1)}`$ stable configurations can be considered as the states of a finite, irreducible and aperiodic Markov chain with a unique stationary distribution (Durrett, 1999).

The evolution rules of the model produce an earthquake size-frequency relation, $`p_k`$, that is shown in Fig. 1a, where the results for $`N=10`$, $`N=100`$, and $`N=1000`$ are superimposed. Note that this spectrum has a distribution of the characteristic-earthquake type: it exhibits a power-law relationship for small events, an excess of maximal (characteristic) events, and very few of the intermediate size. Besides, the three superimposed curves of probability are coincident.

The result for the probability of return of the characteristic earthquake, $`P(n)`$, is shown in Fig. 1b for $`N=20`$. Here $`n`$ represents the time elapsed since the last characteristic event. During an initial time interval $`1n<N`$, $`P(n)`$ is null, then it grows to a maximum and then finally declines asymptotically to $`0`$. (In Section IV, $`P(n)`$ for $`N=20`$, will be usually denoted as curve $`a`$.) In Section III, we explain a general algebraic method for the exact computation of $`P(n)`$.

The configurations of the model are classified into groups according to the number of levels, $`j`$, that are occupied $`(0jN1)`$. Using the Markov-chain theory or producing simulations (Vázquez-Prada et al., 2002), one easily observes that in this model the system resides often in the configurations of maximum occupancy, i. e., in $`j=N2`$ and $`j=N1`$.

This last property can be observed in Fig. 1c, where we have represented, for $`N=100`$, the time evolution of the level of occupancy, $`j`$, in an interval long enough to observe the occurrence of several characteristic earthquakes. The typical pattern after a total depletion is a gradual recovery of $`j`$ up to a new high level of occupancy. Once there, the system typically presents a plateau before the next characteristic earthquake. Especially during the ascending recoveries, the level of occupancy $`j`$ suffers small falls corresponding to the occurrence of rather small earthquakes, that in this model are abundant. Of course, one also observes that occasionally $`j`$ falls in a significant way corresponding to the occurrence of a $`N>kN/2`$ intermediate earthquake.

Due to the fact that this model is not critical, it is reasonable to consider it as an example of self-organized complexity (Gabrielov et al., 1999).

## III  Algebraic Approach to P(n)

The function $`P(n)`$, for a minimalist system of size $`N`$, is obtained from the Markov matrix of the system, $`𝐌`$, following the following three steps: i) The element of the last row and first column of $`𝐌`$ is changed by a 0. After this pruning, the matrix will be called $`𝐌^{}`$. ii) The new matrix $`𝐌^{}`$ is multiplied by itself $`n1`$ times to obtain $`𝐌^{(n1)}`$ and the element of the first row, last column of this matrix is identified. iii) $`P(n)`$ is the product of this selected matrix element times $`1/N`$.

The whys of this recipe are explained in Vázquez-Prada et al. (2002), where $`P(n)`$ for $`N=2`$ is explicitly obtained by mere inspection. The result is

$$P(n)=\frac{n1}{2^n},N=2.$$
(1)
The explicit form of $`P(n)`$ for larger values of $`N`$, can be achieved by exploiting the Jordan decomposition of $`𝐌^{}`$,

$$𝐌^{}=𝐐𝐉𝐐^1,$$
(2)
and hence,

$$𝐌^{n1}=𝐐𝐉^{n1}𝐐^1.$$
(3)
The matrix $`𝐉`$ is formed by “Jordan blocks” in the diagonal positions, i.e., by square matrices whose elements are zero except for those on the principal diagonal, which are all equal, and those on the first superdiagonal, which are equal to unity. Thus, the task of obtaining an arbitrary power of $`𝐉`$ is simple because, as said, each Jordan block is the sum of two conmuting matrices: one is a constant times the unity matrix, and the other is nilpotent. Therefore, in the computation of any arbitrary power of $`𝐉`$, each block is independent and the corresponding Newton bynomial formula can be applied. As an example, we now present the calculation of the case $`N=3`$. In this case,

$$3𝐌^{}=\left(\begin{array}{cccc}1& 1& 1& 0\\ 0& 2& 0& 1\\ 1& 0& 1& 1\\ 0& 0& 0& 2\end{array}\right),$$
(4)
which is decomposed as

$$\left(\begin{array}{cccc}1& 1& 1& 0\\ 0& 2& 0& 1\\ 1& 0& 1& 1\\ 0& 0& 0& 2\end{array}\right)=\left(\begin{array}{cccc}\hfill 1& \hfill 1& \hfill 1& \hfill 2\\ \hfill 0& \hfill 0& \hfill 2& \hfill 1\\ \hfill 1& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 2\end{array}\right)\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 2& 1& 0\\ 0& 0& 2& 1\\ 0& 0& 0& 2\end{array}\right)\left(\begin{array}{cccc}\hfill 1/2& \hfill 1/4& \hfill 1/2& \hfill 3/8\\ \hfill 1/2& \hfill 1/4& \hfill 1/2& \hfill 3/8\\ \hfill 0& \hfill 1/2& \hfill 0& \hfill 1/4\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1/2\end{array}\right).$$
(5)
Therefore,

$$𝐉^{n1}=\left(\begin{array}{cccc}\hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 2^{n1}& \hfill (n1)2^{n2}& \hfill 1/2(n1)(n2)2^{n3}\\ \hfill 0& \hfill 0& \hfill 2^{n1}& \hfill (n1)2^{n2}\\ \hfill 0& \hfill 0& \hfill 0& \hfill 2^{n1}\end{array}\right)\left(\frac{1}{3}\right)^{n1},$$
(6)
and from eq. (3)

$$𝐌_{1,4}^{n1}=\frac{2^n}{32}(n2)(n+5)\left(\frac{1}{3}\right)^{n1}.$$
(7)
Thus, finally,

$$P(n)=\left(\frac{2}{3}\right)^n\frac{(n2)(n+5)}{32},N=3.$$
(8)
One could optimistically guess that $`P(n)`$, for an arbitrary $`N`$, can be deduced from the systematics observed in the previous low-$`N`$ cases. This is disproved by the following formula, which is the result of $`P(n)`$ for $`N=4`$.

$$P(n)=\left(\frac{1}{4}\right)^n\left[\frac{13}{16}+\frac{7n}{4}\frac{n^2}{2}+\frac{n^3}{32}+3^n\left[\frac{3}{16}+\frac{7n}{324}+\frac{n^2}{108}+\frac{n^3}{2599}\right]\right],N=4.$$
(9)
Although it is not apparent, this formula, as it should, vanishes for $`n=3`$. As in eqs. (1) and (8), $`P(n)`$ in eq. (9) is adequately normalized:

$$\underset{n=N}{\overset{\mathrm{}}{}}P(n)=1.$$
(10)
## IV  Error Diagram for the Forecasting of the Characteristic Earthquake

In the following paragraphs, we will stick to a model of size $`N=20`$ to make the pertinent comparisons. This size is big enough for our purposes here, and small enough to obtain good statistics in the simulations.

For $`n=20`$ the mean value of $`P(n)`$ is

$$<n>=\underset{i=20}{\overset{\mathrm{}}{}}P(i)i=121.05,$$
(11)
the standard deviation is

$$\sigma =\left[\underset{i=20}{\overset{\mathrm{}}{}}P(i)(i<n>)^2\right]^{1/2}=55.21,$$
(12)
and the skew of the distribution is

$$\gamma =\frac{1}{\sigma ^3}\underset{i=20}{\overset{\mathrm{}}{}}P(i)(i<n>)^3=0.10.$$
(13)
Now we enter into the matter of forecasting. As in any optimization strategy, we will try to achieve simultaneously the most in a property called A and the least in a property called B, these two purposes being contradictory in themselves. Here A is the (successful) forecast of the characteristic earthquakes produced in the system. Our desire is to forecast as many as possible, or ideally, all of them. B is the total amount of time that the earthquake alarm is switched on during the forecasting process. As is obvious, our desire would be that this time were a minimum. The maximization of A is equivalent to the minimization of an A that represents the fraction of unsuccessful forecasts.

Thus, in practice, our goal in this paper is to obtain simultaneously a minimum value for the two following functions, $`f_e`$ and $`f_a`$. The first represents the fraction of unsuccessful forecasts, or fraction of failures; the second represents the fraction of alarm time. These two functions, in this first one-dimensional strategy of forecasting, are dependent only on the value of $`n`$, that is, the time elapsed since the last main event, and to which the alarm is connected. Using the function $`P(n)`$ previously defined, they read as follows:

$$f_e(n)=\underset{n^{}=1}{\overset{n}{}}P(n^{}),$$
(14)
$$f_a(n)=\frac{_{n^{}=n}^{\mathrm{}}P(n^{})(n^{}n)}{_{n^{}=0}^{\mathrm{}}P(n^{})n^{}}.$$
(15)
These two functions are plotted in Fig. 2a. By eliminating $`n`$ between $`f_e(n)`$ and $`f_a(n)`$, we obtain Fig. 2b, which is the standard form of representing the so-called error diagram. The diagonal straight line would represent the result of a random forecasting strategy. The curved line is the result of this model.

Error diagrams were introduced in earthquake forecasting by Molchan who contributed with rigorous mathematical analysis to the optimization of the earthquake prediction strategies (Molchan, 1997). In his papers Molchan used $`\tau `$ and $`n`$ to represent the alarm fraction and the error fraction respectively; and put $`\tau `$ in the horizontal axis.

To fix ideas, it is convenient to define a so-called loss function, $`L`$, which expresses the trade-off between costs and benefits in the forecasting (Keilis-Borok, 2002). Among all the possible loss functions, we will choose the simple linear function

$$L=f_a+f_e.$$
(16)
$`L(n)`$ is also drawn if Fig.2a. The position $`n_a=66`$ provides the minimum value of $`L(n)`$. $`L(n_a)=0.578`$. Note that $`n_a`$ does not coincide either with the $`n`$ that maximizes $`P(n)`$, or with $`<n>`$.

## V  Improving the Forecasts

In Section IV we adopted the strategy of connecting the alarm at a fixed time, $`n`$, after the occurrence of a characteristic event. The evaluation of this strategy leads to the conclusion that for $`n=n_a=66`$, the loss function has a minimum value $`L(n_a)=0.578`$. The question now is: Can we think up other strategies that render better results? To answer this question, we now return to our previous comments on Fig. 1.

If we define a medium-size earthquake as an event with a size between $`N/2`$ and $`N1`$, i.e. $`N>kN/2`$, by observing the graphs in Fig. 1, one is led to the conclusion that in this model the occurrence of a medium-size earthquake is not frequent but when it actually takes place, the time of return of the characteristic quake in that cycle is increased.

This qualitative perception can be substantiated by numerically obtaining the probability of having cycles where no medium-size earthquake occurs, i.e., $`k<N/2`$. This information is completed by the distribution of cycles where the condition $`N>kN/2`$ does occur. These two distributions are shown in Fig. 3a as lines $`b`$ and $`c`$ respectively. Here, line $`a`$ represents the total distribution of the times of return of the characteristic earthquake in this model (the same as plotted in Fig. 1b). Note that, as it should, the distribution $`a`$ covers both distributions $`b`$ and $`c`$. The mean time $`<n>`$ for the three distributions is $`<n>_a=121.05`$, $`<n>_b=107.57`$ and $`<n>_c=166.84`$. The fraction of cycles under $`b`$ is $`0.77`$ and the fraction under $`c`$ is $`0.23`$. A splitting of this type, in which the $`a`$ distribution separates into $`b`$ and $`c`$, will be denoted henceforth as $`a=bc`$.

To check if $`a=bc`$ is potentially useful for our purposes, we will now analyze independently these two sets of cycles, $`b`$ and $`c`$, with the method used in Section IV for curve $`a`$. The result is the following: the best working $`n`$ for dealing with the cycles under distribution $`b`$ is $`n_b=60`$. And with respect to the cycles belonging to the distribution under $`c`$, the best $`n`$ is $`n_c=124`$.

Therefore, we will now study again the whole set of cycles, i.e. those under $`a`$, by means of a retarding strategy, which is based on the splitting $`a=bc`$. We will adopt the following steps: in any cycle, we will wait until an $`n`$, named $`n_{ret1}`$ (which is near to $`n_b`$), before taking any decision. If no medium earthquake has occurred so far, then the alarm is connected at $`n_{ret1}`$. If, on the contrary, a medium quake has occurred before $`n_{ret1}`$, then we move the alarm to $`n_{ret2}`$ (which is close to $`n_c`$). This notation $`n_{ret1}`$ and $`n_{ret2}`$ comes from the retarding strategy that we are exploring now. This two-dimensional strategy is implemented by varying $`n_{ret1}`$ and $`n_{ret2}`$ looking for the best value of $`L`$. This is illustrated in Fig. 4a. The best option is $`n_{ret1}=61`$ and $`n_{ret2}=101`$, with $`L(n_{ret1},n_{ret2})=0.549`$.

Now we look for a similar property that can classify the cycles from another point of view. This new property consists in identifying the cycles where the sum of the sizes of all the earthquakes before the characteristic one is less than $`N/2`$. This condition will be represented by $`SUM<N/2`$. The reason for this choice is that if $`SUM<N/2`$, the system, statistically speaking, tends to reach more rapidly the configurations of maximum occupancy, $`j=N2`$ and $`j=N1`$, and the time of return of the characteristic quake in that cycle tends to be smaller (see Fig. 1c). In Fig. 3b, line $`a`$ represents, as in Fig. 3a, the distribution of return intervals of the characteristic earthquake for all the cycles of the model. And lines $`f`$ and $`g`$ represent, respectively, the separation of line $`a`$ according to the fulfilment, or not, of the $`SUM<N/2`$ condition, $`a=fg`$. The mean value of the $`f`$ and $`g`$ distributions is $`<n>_f=88.78`$ and $`<n>_g=151.69`$ respectively. The fraction of events under the $`f`$ and $`g`$ lines is $`36.96`$ and $`63.04`$ respectively.

This second splitting of the whole set of cycles in the model, $`a=fg`$, can be used as an advancing strategy in parallel to what we did with the retarding strategy. Thus the independent analysis of curve $`f`$ leads to $`n_f=60`$, and the similar analysis of curve $`g`$ leads to $`n_g=90`$.

We will now study again the whole set of cycles (under $`a`$) by means of the advancing strategy, which is based on the splitting $`a=fg`$. Therefore, we proceed as follows: In any cycle, we wait until $`n_{adv1}`$ (which is close to $`n_f`$) before taking any decision. If the condition $`SUM<N/2`$ has been fulfilled, then the alarm is connected at $`n=n_{adv1}`$. If, on the contrary, this condition has not been fulfilled, then we move the alarm to $`n_{adv2}`$ (which is close to $`n_g`$). This two-dimensional strategy is implemented by varying $`n_{adv1}`$ and $`n_{adv2}`$ looking for the lowest $`L`$. This is illustrated in Fig. 4b. The search for the best option leads to $`n_{adv1}=61`$ and $`n_{adv2}=90`$, with $`L(n_{adv1},n_{adv2})=0.537`$. This value of $`L`$ is slightly better than that obtained using the retarding strategy.

Inspired by these results, we will now analize the possibilities of a mixed strategy which contains conceptual elements of the two partial strategies discussed so far. Here we will explore a 3-dimensional grid of points $`(n_1,n_2,n_3)`$ looking for the minimization of $`L`$. The first coordinate, $`n_1`$,will be explored in the neighbourhood of $`n_{adv1}`$, the second coordinate $`n_2`$ in the neighbourhood of $`n_{adv2}`$, and finally $`n_3`$ near $`n_{ret2}`$. The two succesive key decisions to be taken are:

1. In any cycle, we wait until $`n=n_1`$. If $`SUM<N/2`$ IS fulfilled, we connect the alarm at $`n_1`$ and leave it there. If at $`n=n_1`$, $`SUM<N/2`$ is NOT fulfilled, we move the alarm to $`n_2`$. And,
2. (We are now at $`n_2`$). If no medium-size event has occurred between $`n_1`$ and $`n_2`$, we leave the alarm connected at $`n_2`$. If, on the contrary, one or more medium-size events have occurred in this interval, then we move the alarm to $`n_3`$.

The search for the triplet $`(n_1,n_2,n_3)`$ that makes $`L`$ minimum is illustrated in Fig. 5. The result corresponds to $`(n_1=61,n_2=84,n_3=104)`$, and there, $`L=0.528`$. This is the best result obtained in this work.

Thus, the improvement obtained in $`L`$, when passing from $`L(n_a)`$ to $`L(n1,n2,n3)`$ is around $`8.5\%`$.

## VI  Conclusions

In this paper, we have analyzed the behaviour of the minimalist model in relation to a quantitative assesment of the forecasting of its successive characteristic earthquakes. We have chosen a simple loss function, $`L=f_a+f_e`$. Our first try, based on a one-dimensional search in $`n`$, produces a minimum result of $`L`$ around $`0.578`$. This was illustrated in Fig. 2a. With the aim of improving the forecasts, we then explored two modes of a common strategy that divides the probability distribution of the time of return of the characteristic earthquake into two distinct distributions. The first mode consists in using the occurrence of intermediate-magnitude earthquakes as a sign that the characteristic earthquake would likely return at a time later than usual in that cycle. This is based on the fact that medium-size events significantly deplete the load in the system and its recovery induces a retardation. This effect takes place in any system of the sand-pile type. The exploitation of this idea leads to a two-dimensional search that finally renders an $`L`$ value around $`0.549`$. (Fig. 4a). The second idea consists in using the fact that a significant absence of small earthquakes during a sizeable lapse of time in the cycle is a sign of imminence of the next characteristic event, or at least of a shortening of its period of return. This strategy is similar to the old wisdom in seismology that links a steady absence of earthquakes in a fault with the increase in the risk of occurrence of a big event. The exploration of this idea proceeds similarly to what we did with the retarding strategy: this also leads to a two-dimensional search. It renders a minimum $`L`$ around $`0.537`$. (Fig.4b).

Finally, a mixed strategy that tries to incorporate the information acquired is implemented by means of a three-dimensional search, and provides a value of $`L=0.528`$. The identification of the three optimum parameters is illustrated in Fig. 5.

It is important to remark that the information we have used in our forecasts is based only in the observed systematics of earthquake occurrence in the model, i.e., only seismicity functions have been used. Thus, for example, in Section V we have not used the state of occupancy of the system $`j`$, which would have given much more accurate predictions. In real life, the use of this information would be equivalent to knowing, in real time, the value of the stress level and the failure threshold at any point in a fault.

### Acknowledgements

We are grateful to the two Referees of the first version of this paper for their thorough and helpful reviews. This work was supported by the project BFM2002-01798 of the Spanish Ministry of Science. Miguel Vázquez-Prada and Álvaro González are respectively supported by the PhD research grants B037/2001 (funded by the Autonomous Government of Aragón and the European Social Fund) and AP2002-1347 (funded by the Spanish Ministry of Education).

## References

Ben-Zion, Y., Eneva, M., and Liu, Y.: Large Earthquake Cycles and Intermittent Criticality On Heterogeneous Faults Due To Evolving Stress And Seismicity, J. Geophys. Res., 108 (B6), 2307, doi: 10.1029/2002JB002121, 2003.
Dahmen, K., Ertas, D., and Ben-Zion, Y.: Gutenberg-Richter and characteristic behaviour in simple mean-field models of heterogeneous faults, Phys. Rev., E58, 1444-1501, 1998.
Drossel, B., and Schwabl, F.: Self-organized critical forest-fire model, Phys. Rev. Lett., 69, 1629-1632, 1992.
Durrett, R.: Essentials of Stochastic Processes, Chapter 1, Springer, 1999.
Gabrielov, A., Newman, W. I., and Turcotte, D. L.: Exactly soluble hierarchical clustering model: inverse cascades, self-similarity and scaling, Phys. Rev., E60, 5293-5300, 1999.
Gabrielov, A., Zaliapin, I., Newman, W., and Keilis-Borok, V. I.: Colliding cascades model for earthquake prediction, Geophys. J. Int., 143, 427-437, 2000.
Keilis-Borok, V.: Earthquake Prediction: State-of-the-Art and Emerging Possibilities, Annu. Rev. Earth Planet. Sci., 30, 1-33, 2002.
Keilis-Borok, V. I., and Soloviev, A. A. (Eds.): Nonlinear Dynamics of the Lithosphere and Earthquake Prediction, Springer, 2002.
Lomnitz, C.: Fundamentals of Earthquake Prediction, J. Wiley & Sons, 1994.
Molchan, G. M.: Earthquake Prediction as a Decision-making Problem, Pure. Appl. Geophys., 149, 233-247, 1997.
Newman W. I., and Turcotte, D. L.: A simple model for the earthquake cycle combining self-organized complexity with critical point behavior, Nonlin. Proces. Geophys., 9, 453-461, 2002.
Sornette, D.: Critical Phenomena in Natural Sciences, Springer Verlag, Berlin, Germany, 2000.
Schwartz, D.P. and Coppersmith, K. J.: Fault behaviour and characteristic earthquakes: Examples from the Wasatch and San Andreas Fault Zones, J. Geophys. Res., 89, 5681-5698, 1984.
Turcotte, D. L.: Fractals and Chaos in Geophysics, 2nd Edit., Cambridge Univ. Press, 2000.
Vázquez-Prada, M., González, Á., Gómez, J. B., and Pacheco, A. F.: A Minimalist model of characteristic earthquakes, Nonlin. Proces. Geophys., 9, 513-519, 2002.
Wesnousky, S. G.: The Gutenberg-Richter or Characteristic-Earthquake Distribution, Wich It Is?, Bull. Seism. Soc. Am., 84, 1940-1959, 1994.


text:
# Elliptic flow measurements from STAR

## 1 Introduction

Elliptic flow characterizes the anisotropy in particle emission “in” and “out” of the reaction plane. The word flow is used to describe collective behavior but does not necessarily imply a hydrodynamic interpretation. Elliptic flow is commonly characterized by the second harmonic coefficient $`v_2`$ of an azimuthal Fourier decomposition of the the particle momentum distribution versus the reaction plane.

Based on general arguments it is thought that elliptic flow develops mostly in the first few fm/$`c`$ ($`<`$ radius of the nucleus) and thus provides information about the amount of thermalization achieved early in the collision. In fact, the observed elliptic flow for charged and identified particles at RHIC is interpreted as:

* one needs very strong interactions between the quarks and gluons at very early times in the collision .
* a well developed quark-gluon plasma

## 2 Integrated elliptic flow and non-flow contributions

Experimentally the reaction plane is not known, and elliptic flow is often reconstructed from two-particle azimuthal correlations. Two-particle azimuthal correlations can be affected by many other sources besides elliptic flow. These so called non-flow effects could be large  and this would change the interpretation of strong re-interactions of the constituents early in the collision. STAR has estimated the contribution of non-flow in the first elliptic flow paper from RHIC .

Figure 2 taken from  shows the correlation between two event plane angles for two independent subevents for the first and second harmonic as a function of centrality. In the case of flow the correlation of the second harmonic event plane angles will be proportional to M$`v_2^2`$, where M is the multiplicity of particles used in the determination of the event plane. The observed peaked shaped as a function of centrality is characteristic of elliptic flow. This peaked shape originates from the fact that M increases as a function of centrality while elliptic flow decreases. Every model description of a heavy-ion collision which includes final state interactions will yield such a shape, the magnitude and peak position as a function of centrality will depend on the amount of re-interactions. Non-flow contributions will be monotonic or almost constant for this quantity, which is also true in the specific case of non-flow due to mini-jets as calculated in . Based on this correlation the estimated maximum non-flow contribution, taken constant as a function of centrality, was 0.05. The propagation of this estimate to the measured $`v_2`$ values is shown in Fig. 2. In this figure the $`v_2`$ values versus centrality are shown as open circles and the asymmetric uncertainties are the non-flow estimates .

More recently a new analysis method based on cumulants was proposed  which utilizes the fact that true flow is a multi-particle correlation. The obtained elliptic flow values using this method for four particle correlations, $`v_2`${4}, is also shown in Fig. 2 as solid stars. A detailed description of the cumulant analysis in STAR for the $`\sqrt{s_{_{NN}}}=130`$ GeV data is given in Ref. . The reduction in the integrated elliptic flow between the two methods shows the (possible) contribution of non-flow effects at this energy. There are however two important caveats associated with calculating $`v_2`$ using two or multi-particle correlations. First the four particle cumulant is reliable when the magnitude of $`v_2`${4} $`>1/N^{3/4}`$ , where N is the number of “clusters” contributing to the measurement. However, this is an important improvement over the two particle analysis which is reliable when $`v_2`${2} $`>1/N^{1/2}`$. Secondly, event by event fluctuations could affect $`v_2`${2} differently than $`v_2`${4}. The two particle analysis is sensitive to $`v_2^2`$ while the four particle cumulant analysis is sensitive to the difference between $`v_2^4`$ and $`v_2^2`$. The $`v_2`$ value is obtained by averaging over events and due to event by event fluctuations in general $`v_2^2v_2^2`$ and $`v_2^4v_2^4`$. This can also lead to a reduction in the $`v_2`$ obtained from the four particle compared to the two particle correlation methods .

Even with these caveats, the most reliable way of calculating $`v_2`$ when the reaction plane is not known is via a higher order cumulant method. Figure 4 shows the measured $`v_2`${4} values as a function of centrality for $`\sqrt{s_{_{NN}}}=`$ 130 and 200 GeV. The elliptic flow as a function of centrality for both energies is very similar. The maximum difference is about 20%, however note that only the statistical uncertainties are shown. For the preliminary 200 GeV results the systematic uncertainties are also about 20%. Also shown are the measurements from NA49 at the CERN SPS. The increase in the integrated elliptic flow between SPS and RHIC energies can be caused by an increase in the mean transverse momentum or due to a higher slope of the differential $`v_2`$($`p_t`$. Figure 4 shows the differential $`v_2`$($`p_t`$) for SPS and RHIC energies. While the slope increases, which indicates the increase in elliptic flow at RHIC, the dominant contribution to the integrated elliptic flow comes from the value of $`v_2`$ around mean $`p_t`$ of the particles. In the $`p_t`$ range of $`350500`$ MeV/$`c`$, the $`v_2`$($`p_t`$) values at SPS and RHIC energies are very close. However the mean $`p_t`$ difference between pions at the SPS and charged particles at RHIC is about 150 MeV/$`c`$ which already accounts for most of the difference in integrated $`v_2`$.

## 3 Elliptic flow for identified particles

The elliptic flow as a function of transverse momentum, $`v_2`$($`p_t`$), depends on the temperature, radial flow velocity, azimuthal variation of the transverse flow velocity and the spatial anisotropy of the system at freeze-out . The measurement of $`v_2`$($`p_t`$) for different particle masses constrains these parameters with reasonable precision . These freeze-out parameters in combination with the initial conditions strongly constrain the Equation Of State (EOS). In Fig. 6 the measured $`v_2\{2\}`$ versus $`p_t`$ at $`\sqrt{s_{_{NN}}}=`$ 130 GeV for pions and protons + antiprotons is shown  together with hydrodynamical model predictions  for two different EOS. This clearly illustrates that the heavier particles are more sensitive to the underlying EOS. This can be understood from the fact that the lighter particles are very sensitive to the freeze-out temperature while the heavier particles more directly reflect the flow. From Fig. 6 it is clear that for these model calculations the data prefer the quark-gluon plasma EOS. Figure 6 shows in solid symbols the $`v_2\{2\}`$ versus $`p_t`$ for pions, $`K_S^0`$’s, protons + antiprotons and $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ at $`\sqrt{s_{_{NN}}}=`$ 200 GeV. For comparison the $`v_2`$($`p_t`$) for pions and protons + antiprotons at $`\sqrt{s_{_{NN}}}=`$ 130 GeV in open symbols are included. The lines shown are the results from the blastwave fits to the $`\sqrt{s_{_{NN}}}=`$ 130 GeV data . From this comparison it is clear that the $`v_2`$($`p_t`$) for different particles at these different energies is very similar. The pions do show however a slightly higher slope at $`\sqrt{s_{_{NN}}}=`$ 200 GeV.

## 4 Elliptic flow at intermediate transverse momentum

Elliptic flow measurements can quantify the possible modifications of the created medium on the particle yields as a function of $`p_t`$. A medium modification like the predicted mechanism of parton energy loss, jet quenching, will inevitably lead to a finite $`v_2`$ at high-$`p_t`$.

Figure 8 shows the measured $`v_2`$($`p_t`$) up to 7 GeV/$`c`$ obtained by two, $`v_2`${2} and $`v_2`${RP}, and four, $`v_2`${4}, particle correlation methods. The measurement of $`v_2`${4} shows that there is a significant amount of elliptic flow out to $`p_t=56`$ GeV/$`c`$. The medium modification inferred from this observable in addition to the suppression of the single inclusive particle yield , characterized by RAA, and the disappearance of high-$`p_t`$ angular back to back correlations, as seen in the centrality dependence of IAA , is in qualitative agreement with the predictions of jet quenching. At high-$`p_t`$ the measurement of RAA, IAA and $`v_2`$ for the different particle species can provide a better constraint on the underlying mechanism responsible for the modification of the particle yield. That RAA depends on particle species was observed by PHENIX  for $`\pi ^0`$’s and charged particles. The comparison of the single inclusive particle yield for $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ and $`K_S^0`$ in central and peripheral collisions as measured by STAR also shows this species dependence . The suppression in the single inclusive particle yield at intermediate $`p_t`$ for $`K_S^0`$ is stronger than for $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$. The measured $`v_2\{2\}`$, shown in Fig. 8, for $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ and $`K_S^0`$ at the same intermediate $`p_t`$ is also particle dependent. The $`v_2`$($`p_t`$) above 2 GeV/$`c`$ for the $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ is larger than the $`v_2`$($`p_t`$) of $`K_S^0`$. It is assumed that the non-flow contributions in this measurement are approximately equal for both particle species and therefore the difference is a real difference in flow.

Explaining the origin of this suppression of single inclusive particle yield and elliptic flow for identified particles at intermediate-$`p_t`$ due to jet quenching alone seems not to work. To understand this observed behavior better a measurement of the contribution of the initial state Cronin effect in dA as well as a measurement of IAA for $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ and for $`K_S^0`$ is needed. Another possible explanation is that the hydro-like behavior extends further in $`p_t`$ for $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ then for $`K_S^0`$.

Another puzzle is the rather large values of elliptic flow at intermediate $`p_t`$. The interpretation that this is caused by radiative energy loss or due to inelastic interactions in a parton cascade would lead to rather large initial parton densities. Corrections for non-flow would already reduce the elliptic flow values and therefore the parton densities needed. In addition a mechanism like parton coalescence would require much smaller parton elliptic flow to account for the measured particle elliptic flow and would also explain the different $`v_2`$ values observed at intermediate-$`p_t`$ for $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ and $`K_S^0`$ .

## 5 Conclusions

Large elliptic flow values have been measured by STAR both at $`\sqrt{s_{_{NN}}}=`$ 130 and 200 GeV Au+Au collisions. The large magnitude of the charged and identified particle elliptic flow at low-$`p_t`$ is interpreted as due to strong interactions between the partons at early times in collision and even approaches the ideal hydrodynamical limit. While non-flow is not negligible at higher-$`p_t`$, elliptic flow extends at least up to 6 GeV/$`c`$. This unambiguously shows the effect of medium modification on the particle yield in this transverse momentum range. The dependence of the elliptic flow and the single particle inclusive yield on particle species at intermediate-$`p_t`$ shows that radiative energy loss can not be the only medium induced modification of the particle yield in this $`p_t`$-range.


text:
# Chapter 1 Identification of Showers with Cores Outside the ARGO-YBJ Detector

## Chapter 1 Identification of Showers with Cores Outside the ARGO-YBJ Detector

### Abstract

In any EAS array, the rejection of events with shower cores outside the detector boundaries is of great importance. A large difference between the true and the reconstructed shower core positions may lead to a systematic miscalculation of some shower characteristics. Moreover, an accurate determination of the shower core position for selected internal events is important to reconstruct the primary direction using conical fits to the shower front, improving the detector angular resolution, or to performe an efficient gamma/hadron discrimination.

In this paper we present a procedure able to identify and reject showers with cores outside the ARGO-YBJ carpet boundaries. A comparison of the results for gamma and proton induced showers is reported.

### 1. Introduction

Showers of sufficiently large size can trigger a detector even if their core is located outside its boundaries. The corresponding core positions are generally reconstructed not only near the carpet edges but also well inside the boundaries. As a consequence, sofisticated algorithms able to reduce the contamination of external events are needed. The goal is to identify and reject a large fraction of external events before exploiting any reconstruction algorithm, only by using some suitable parameters.

In this paper we present a reconstruction procedure able to identify and reject a large fraction of showers with cores outside the ARGO-YBJ detector.

### 2. Identification of external events

The ARGO-YBJ detector consists of a single layer of RPCs with dimensions of $`74\times 78m^2`$. The area surrounding this central detector (carpet), up to $``$ 100 $`\times `$ 110 $`m^2`$, is partially ($`50\%`$) instrumented with RPCs (guard-ring). The basic element is the logical pad ($`56\times 62cm^2`$) which defines the time and space granularity of the detector. The detector is divided in 6 $`\times `$ 2-RPC units (clusters): the central carpet contains 10 $`\times `$ 13 clusters. For a detailed description of the ARGO-YBJ detector see .

Various parameters based on particle density or time information are under investigation to identify showers with core position outside a given fiducial area. The most interesting ones are the following: (1) position of the cluster with the highest particle density, (2) position of the cluster row/column with the highest total particle density, (3) mean distance $`R_p`$ of all fired pads to the reconstructed shower core.

To perform these calculations we have simulated, via the Corsika code , gamma and proton induced showers with energy spectra ($`E^{2.5}`$ and $`E^{2.75}`$, respectively) ranging from 100 GeV to 50 TeV. The detector response has been simulated via a GEANT-3 based code.

As an example, in Fig. 1 we show the distributions of the position of the cluster with the highest particle density for $`\gamma `$-induced showers. In the plots we compare the events with the core really external to a 80 $`\times `$ 80 m2 fiducial area (solid histograms) and the truly internal ones (dashed histograms). To investigate the discrimination power of this particular parameter we have simulated a detector completely instrumented up to $``$ 100 $`\times `$ 110 $`m^2`$, i.e., containing 14 $`\times `$ 17 clusters. Therefore, the cluster coordinates run from 1 to 14 (X view) and from 1 to 17 (Y view) starting from the lower left corner of the carpet.

The $`R_p`$ distribution for showers reconstructed inside a 80 $`\times `$ 80 m2 fiducial area is shown in Fig. 2 (solid histogram). The dashed line refers to truly IN events while the dotted histogram refers to OUT showers erroneously reconstructed as internal. The shower cores have been calculated by means of the simple center of gravity method. As can be seen, the parameter $`R_p`$ identifies quite well the events with core outside the carpet. Large distances between the truly and the reconstructed shower axis lead to larger $`R_p`$ values. This fact offers the possibility to define a cut in $`R_p`$ to identify these events. A conservative choice is to reject showers with $`R_p>25`$ m.

From these studies it follows that the identification of a large fraction of external events can be achieved by defining a suitable fiducial area togheter with a combination of cuts in the parameters discussed above.

### 3. Maximum Likelihood Method (LLF)

Different algorithms have been investigated to reconstruct the shower core position in the ARGO-YBJ experiment . The most performant is the Maximum Likelihood Method. We point out that expression for -LLF of  refers to the case of a Poisson distribution in which the pads are not fired with probability $`P_i(0)`$ or fired with probability $`P_i(>0)=1P_i(0)`$ (hereafter ’LLF1 method’). In our study almost always the fired pads have particle multiplicity 1, and therefore such a simple discrimination can be made. However, if we consider a larger area as the whole RPC, the multiplicity can be $`>`$ 1, and the proper Poisson distribution on the fired RPCs appears more adequate. In this case the sum on fired elements is:

$$\mathrm{\Sigma }_jlnP_j(>0)=\mathrm{\Sigma }_jN_jln(\rho _j)ln(S_{RPC})\mathrm{\Sigma }_jN_j+\mathrm{\Sigma }_jln(N_j!)+S_{RPC}\mathrm{\Sigma }_j\rho _j$$
(1)
where $`N_e\rho _j`$ is the particle density expected on the j-th RPC at a distance $`R_j`$ from the core, $`N_j`$ is the recorded particle number and $`S_{RPC}`$ is the RPC area. The shower size can be calculated via the equation

$$N_e=\frac{\mathrm{\Sigma }_jN_j}{S_{RPC}\mathrm{\Sigma }_j\rho _j}.$$
(2)
We define this calculation the ’LLF2 method’. As a consequence, we expect that the differences between LLF1 and LLF2 increase with the particle density, for a fixed area. In Fig. 3 we compare the shower core position resolution calculated by applying the LLF1 and LLF2 methods on the RPCs for $`\gamma `$-induced showers with the core randomly sampled inside a 80 $`\times `$ 80 m2 area. As expected, the resolution worsens with multiplicity if the LLF1 approach is applied when the number of particles hitting the RPC is $`>`$ 1. We note that for very low multiplicities ($`N_{hit}<80`$) the method LLF1 is more performant than LLF2. In fact, the algorithm based on RPC occupancy (LLF1) provides a better representation of the hit distribution in very poor showers.

For very high multiplicities ($`N_{hit}>10^3`$) the shower core position is determined by LLF2 with an uncertainty $`<`$ 1 m.

### 4. Results

A possible procedure to reject external events in the ARGO-YBJ detector is the following one: (1) Rejection of the events whose highest density clusters are on the guard ring (X = 1, 14; Y = 1, 17) or on the boundaries of the central carpet (X = 3, 12; Y = 3, 15). (2) Rejection of the events whose highest total density rows or columns are respectively in positions $`\{1,3,15,17\}`$ or $`\{1,3,12,14\}`$. (3) Reconstruction of core coordinates $`\{X_c,Y_c\}`$ using the Maximum Likelihood Method. (4) Further rejection of events with $`R_p>25`$ m.

In Fig. 4 the fraction of events (internal and external to an area of 80 $`\times `$ 80 m2, respectively) rejected after the steps (1) - (4) is reported. As can be seen, this procedure is able to identify and reject a large fraction of external events. For low multiplicities ($`N_{hit}<100`$) a significative fraction of internal events is erroneously rejected, especially in proton-induced showers.

### 5. References

1. Heck D. et al. 1998, Report FZKA 6019 Forschungszentrum Karlsruhe.

2. Martello D., Bleve C., Di Sciascio G. 2001, 27th ICRC, Hamburg, 7, 2927.

3. Surdo A. et al. 2003, in this proceedings.


text:
# Two-way shape memory effect and mechanical properties of Pulse Discharge Sintered Ni2.18Mn0.82Ga

## I Introduction

Beginning from 1996 when Ullakko and co-worker 1-u  reported on a magnetic-field-induced strain of about 0.2% in a non-stoichiometric Ni-Mn-Ga single crystal, research in this field has attracted a considerable attention due to the great technological potential of this effect. Later this phenomenon was found in other compounds exhibiting shape memory effect when in the ferromagnetic state, such as Fe-Pd and Fe-Pt. 2-j ; 3-k  Whereas Fe-Pt and Fe-Pd show magnetic-field-induced strains not exceeding 1%, in Ni-Mn-Ga single crystals the values of magnetic-field-induced strains can be as large as 6% (Ref. 4). The mechanism of this phenomenon is believed to be the redistribution of different twin variants under action of a magnetic field. 5-o  Since these giant strains are attained in comparatively low magnetic fields (less then 1.5 T), they are easy suppressed by an external stress. For instance, the 6% field-induced strain was completely blocked by a compressive stress of order 2 MPa. 4-m  This suggests that the giant magnetostrains arising from the process of twin-boundary motion might be useful for large stroke and small force applications.

Another way to attain a large magnetic-field-induced strain in the ferromagnetic shape memory alloys is a shift of the martensitic transition temperature caused by a magnetic field. 6-v ; 7-d  However, in this case the magnitude of the applied field must be high in order to overcome the temperature hysteresis of martensitic transformation. From general consideration, it can be expected that the maximum achievable strain in this case is equal to the striction of the transition or even more if material is trained for the two-way shape memory effect. The work output has to be large (as in conventional shape memory alloys), which can be useful for mediate stroke and large force applications.

From the point of view of widespread use of magnetic-field-induced strains, observed in the ferromagnetic shape memory alloys, there is a need to investigate polycrystalline materials. Previous studies of a series of Ni2+xMn1-xGa ($`x=0.160.20`$) showed that, among the composition studied, Ni2.18Mn0.82Ga is characterized by a considerable striction of transition. 8-t  Because of that we have studied Ni2.18Mn0.82Ga prepared by a Pulse Discharge Sintering (PDS) process.

## II Samples preparation and measurements

Ingots of the Ni2.18Mn0.82Ga composition were prepared by arc-melting of high-purity initial elements. The ingots were annealed at 1100 K for 9 days and quenched in ice water. A part of the arc-melted ingots was used to fabricate PDS samples. For this aim the arc-melted ingots initially were crushed into particles and ground into fine powder with a particle size less than 53 $`\mu `$m. Meshed powder was filled in a graphite die with two graphite punches. The die was set in a pulse discharge system (Sodick Co., Ltd). The pulse discharge system was evacuated to a vacuum of 3 Pa prior to the sintering process. Maximum pressure and temperature during the PDS process were equal to 80 MPa and 1173 K, respectively. Disc-shaped billets with thickness of 6 mm were sintered. X-ray diffraction measurements of the samples, performed in a wide temperature range, showed that the high-temperature austenitic structure has a cubic modification whereas the low-temperature martensitic phase has a complex tetragonally based crystal structure. Samples with dimensions of $`3\times 3\times 6`$ mm3 were spark-cut from the billets. Temperature and magnetic field dependencies of strain were measured by a strain gage technique. For this aim a non-magnetic strain gage with a compensated temperature range from 273 to 423 K was attached along the longest dimension of the samples. The configuration of the experimental setup allowed detecting the relative change in the length of a specimen with an accuracy of 0.005%. The samples were inserted into a variable temperature chamber of a superconducting magnet. Temperature was monitored by a Lake Shore calibrated platinum resistance thermometer with an accuracy of 0.1 K. Stress-strain measurements and compression of samples were done at room temperature by an Instron machine.

## III Results and discussion

Shown in Fig. 1 are the temperature dependencies of Young’s modulus of a Ni2.18Mn0.82Ga PDS sample during heating and cooling. These dependencies were obtained from ultrasonic measurements performed in a temperature interval from 173 to 373 K. The value of Young’s modulus was evaluated through the formula

$$E=\rho \frac{v_s^2(3v_l^24v_s^2)}{(v_l^2v_s^2)},$$

where $`\rho `$, $`v_l`$ and $`v_s`$ are the density of the material and the velocities of longitudinal and shear waves, respectively. Marked dips on the temperature dependencies of Young’s modulus correspond to the direct and reverse martensitic transformations. As evident from Fig. 1, the Young’s modulus of Ni2.18Mn0.82Ga at room temperature is equal to 95 GPa.

The results of compression tests for Ni2.18Mn0.82Ga prepared by the PDS method and by a conventional arc-melting technique are shown in Fig. 2. The compression tests were done at room temperature with the same speed of compression for both the samples. The comparison of these curves clearly indicates that the PDS sample shows higher yield strength than the arc-melted one. This characteristic is of importance for some practical applications, and it makes the PDS materials more attractive in this sense.

Figure 3 shows the temperature dependencies of strain measured in a Ni2.18Mn0.82Ga PDS sample upon cooling and heating in zero and 5 T magnetic fields. In zero magnetic field the sample length monotonously increases upon heating up to the onset of the reverse martensitic transformation, $`A_s=343`$ K. The martensite \- austenite transformation is accompanied by a rapid increase in the sample length, which flattens out at austenite finish temperature $`A_f=357`$ K. Subsequent cooling down results in the direct martensitic transformation at $`M_s=342`$ K, which is accompanied by a shortening of the sample. As evident from Fig. 3, the striction of the transition is approximately 0.18%. The temperature dependencies of strain measured in a magnetic field of 5 T revealed that the striction of the transition remains essentially the same as in the case of the measurements without magnetic field. The results of the measurements performed in zero and 5 T magnetic fields leads to the conclusion that the application of the magnetic field results in an upward shift of the characteristic temperatures of the martensitic transition with a rate of about 1 K/T. This value agrees very well with the results reported for polycrystalline Ni2+xMn1-xGa prepared by arc-melting method. 7-d

An interesting feature of the PDS samples is that the two-way shape memory effect can be induced in these materials by a simple loading - unloading cycle. This feature is presented in Fig. 4. Indeed, the increase in the sample length caused by the martensitic transition is approximately 0.18% in the case of the sample which was not subjected to compression. Another sample cut from the same ingot was compressed for 2% at room temperature in the martensitic state. After unloading the residual deformation was approximately 1.2% (Fig. 2). The sample recovered approximately 75% of its initial length upon the first heating, showing shape memory effect. The subsequent cooling - heating process revealed that for this sample the change in the length associated with the martensitic transformation increased twofold as compared with the uncompressed sample and reached 0.4%. It is also seen from Fig. 4 that the change in slope of the curves at the characteristic temperatures of martensitic transformation in the compressed sample becomes less pronounced than in the compression-free sample. Further themocyclings demonstrated that this compression-induced two-way shape memory effect does not degrade and the 0.4% change in the length of the sample is very well reproducible at least up to the tenth heating - cooling cycle. It is interesting to note that a well-defined two-way shape memory effect has also been found recently in Ni-Mn-Ga thin films. 9-o

To study this compression-induced two-way shape memory effect in more detail, several Ni2.18Mn0.82Ga PDS samples were compressed for values of strain, ranging from 1 to 6%. After unloading the residual deformation in these samples was from 0.4 to 4%, respectively. The behavior of the samples upon the first heating process was found to be dependent on the value of residual deformation. Thus, the sample with 0.4% residual deformation showed a perfect shape memory effect whereas the sample with the largest residual deformation did not exhibit shape memory effect, which means that the residual deformation in this sample was essentially plastic one. The two-way shape memory effect was found only in the samples compressed for 2 and 3%. Together with the observation that the sample compressed for 2 and 3% did not revert the original shape completely after the first heating, these facts evidence that the two-way shape memory effect appears when the applied stress exceeds some critical limit which is enough for the occurrence of an irreversible slip. We suggest that it arises from the strain field of dislocations induced upon compression.

The measurements of magnetic-field-induced strain in the temperature interval of the martensitic transformation in the PDS-sintered Ni2.18Mn0.82Ga samples showed that these samples exhibit rather small values of magnetic-field-induced strain even in high (up to 5 T) magnetic fields. For instance, in the compression-free Ni2.18Mn0.82Ga sample this value is equal to 0.02% in a 5 T magnetic field. However, in the case of the samples which demonstrate the compression-induced two-way shape memory effect the magnetic-field-induced strain is six times greater as compared to the samples without two-way shape memory effect. In fact, such a tendency could be expected, since in this case the magnetic-field-induced strain is due to the shift of the martensitic transition temperature caused by the applied magnetic field. The magnetic-field-induced strain is proportional to the relative change in the dimension of the sample per 1 K and this characteristic is much better in the samples with the two-way shape memory effect. It can be suggested that an appropriate training procedure of the PDS material will results in enhancement of the magnitude of the two-way shape memory effect, leading to an increase in the value of magnetic-field-induced strain. It should be noted, however, that the strain of about 0.12% is not perfectly recovered. Figure 5 shows that the reversible magnetic-field-induced strain is equal to 0.06%. This value of magnetic-field-induced strain has been reversible for many cycles of application and removal of the magnetic field.

## IV Conclusion

In conclusion, the most interesting findings of this study are that the two-way shape memory effect can be induced in the Ni2.18Mn0.82Ga PDS materials by an ordinary compression of the samples in the martensitic state. The samples with the two-way shape memory effect show a significant enhancement in the magnitude of magnetic-field-induced strain observed in the temperature interval of martensitic transformation.

## Acknowledgements

This work was partially supported by the Grant-in-Aid for Scientific Research (C) No. 11695038 from the Japan Society of the Promotion of Science.