Parina ; Patel

Terence Richard Hull
email: trhull@uclan.ac.uk

Anna A Stec

Richard E Lyon

Influence of Physical Properties on Polymer Flammability in the Cone Calorimeter

Keywords: cone calorimeter, retardant, flame, fire, polymer

  L'archive ouverte pluridisciplinaire

Introduction

As the utilisation of polymeric materials steadily embraces a wider variety of potentially hazardous applications, greater emphasis must be placed on mitigating the danger of fire. The physical characteristics of polymers and a better understanding of the behaviour of such materials when exposed to ignition sources is, therefore, a necessity. The ignitability and burning behaviour of polymers is a complex process involving interactions between a number of physical and chemical processes. Improved development of new fire safe materials would result from being able to understand the effect on burning behaviour of altering each variable independently. Unfortunately such studies are not practically feasible since any modification to the polymer, such as the incorporation of an additive, results in changes to a range of physical and chemical properties and processes. In many cases fire retardants (FRs) have chemical effects, such as intumescence, carbonisation, ceramicisation or stabilisation of the polymer. These effects are masked by changes in physical properties, resulting from the incorporation of FR additives, which are highlighted in this study. Where the behaviour of the processes can be reliably predicted, these can be incorporated into models of burning behaviour. Although state of the art models cannot yet make reliable predictions of time to ignition or heat release rate (HRR) history, such predictions are of great value in differentiating between expected, predictable behaviour and unexpected phenomena such as different chemical pathways leading to inhibition of decomposition and pyrolysis.

The development of calorimetric techniques based on the principle of oxygen depletion has greatly improved fire testing and research because it quantifies the heat release associated with real burning [START_REF] Hugget | Estimation of the Rate of Heat Release by Means of Oxygen Consumption Measurements[END_REF] [START_REF] Babrauskas | Development of the Cone Calorimeter[END_REF]. The cone calorimeter, [3] [START_REF] Astm E | Standard test method for heat and visible smoke release rates from materials and products using an oxygen consumption calorimeter[END_REF] developed at the National Bureau of Standards (NBS), now the National Institute of Standards and Technology (NIST), has been widely used for assessing the flammability of polymeric materials. This method was primarily developed for measuring the rate of heat release from a burning object as a function of incident flux [START_REF] Nelson | Polymer Ignition[END_REF]. The external radiation source is intended to simulate the effect of a burning object in close proximity. Interpretation of cone calorimeter data has not been adequately addressed within the fire retardant community [START_REF] Schartel | Application of Cone Calorimetry to the Development of Materials with Improved Fire Performance[END_REF], however, there is still a growing reliance on the instrument [START_REF] Babrauskas | Heat Release Rate: The Single Most Important Variable in Fire Hazard[END_REF] [START_REF] Sacristán | Cone Calorimetry Studies of Fire Retardant Soybean-Oil Based Copolymers Containing Silicon or Boron: Comparison of Additive and Reactive Approaches[END_REF] for the initial screening and subsequent assessment of new flammability modified materials. For example, the use of mineral fillers, such as aluminium hydroxide [START_REF] Zilbermann | Flame retardancy of some ethylenevinyl acetate copolymer-based formulations[END_REF][START_REF] Hull | An Investigation into the Decomposition and Burning Behaviour of EVA Nanocomposite Materials[END_REF] and nanofillers [START_REF] Hull | Fire Retardant Effects of Polymer Nanocomposites[END_REF][START_REF] Wang | Preparation and characterisation of a novel flame retarded PET/α zirconium phosphate nanocomposite (α-ZrP)[END_REF] will change both physical and decomposition behaviour.

A number of studies [START_REF] Staggs | A Theory for Quasi-Steady Single-Step Thermal Degradation of Polymers[END_REF] [14] [START_REF] Jia | The Numerical Simulation of the Noncharring Pyrolysis Process and Fire Development Within a Compartment[END_REF] have effectively demonstrated that a numerical pyrolysis model can be used to determine the relationships between the fundamental physical and chemical properties of polymeric materials and their gasification behaviour. Typically, the model is used to calculate the mass loss rate of a one-dimensional sample of solid fuel exposed to a uniform heat flux. ThermaKin is an example of such a model, which has been effectively utilised as a practical tool for the prediction and/or extrapolation of the results of fire calorimetry experiments [START_REF] Stoliarov | Thermo-Kinetic Model of Burning[END_REF] [17] [START_REF] Stoliarov | The Effect of Variation in Polymer Properties on the Rate of Burning[END_REF] [START_REF] Rahatekar | Effect of Carbon Nanotubes and Montmorillonite on the Flammability of Epoxy Nanocomposites[END_REF]. The model, which combines the absorption and transfer of thermal energy with Arrhenius kinetics for the decomposition of the polymer, predicts the overall behaviour of a pyrolysing object through mass and energy conservation equations. These equations are formulated in terms of rectangular finite elements, each element being characterised by component mass and temperature. Additionally, the model describes the transport of gaseous products through the condensed phase and follows changes in the volume of the bulk material.

For thermally thick solids (typically, thicknesses above 15 mm [START_REF] Quintiere | Principles of Fire Behaviour[END_REF]) the thermal inertia, kρc, the product of thermal conductivity (k), density (ρ) , and specific heat (c), of a material governs its ignition and flame spread properties. This determines the rate of rise in surface temperature and consequently, the time to ignition [START_REF] Drysdale | An Introduction to Fire Dynamics[END_REF]. The time to ignition (t ig ) of a thermally thick solid exposed to a constant net heat flux Q R = Q ext -CHF , where Q ext is the external heat flux from fire or radiant heater and CHF is the critical heat flux for ignition, has been expressed in Equation 1.

ݐ ݅݃ = ߨ 4 ݇ߩܿ ൫ܶ ݅݃ -ܶ 0 ൯ 2 ܳ ሶ ܴ ሖ 2 [Equation 1]
where T ig and T 0 are the ignition and ambient temperatures, respectively. The time to ignition of a thermally thin solid exposed to a constant net heat flux has also been expressed in Equation 2.

ݐ ݅݃ = ߩܿ߬ (ܶ ݅݃ -ܶ 0 ) ܳ ሶ ܴ ሖ [Equation 2]
Where τ refers to material thickness. Equations 1 and 2 follow from the concept of a constant ignition temperature T ign and temperature-independent thermal inertia. Once ignition has occurred and a flame is established on the surface, the net heat flux becomes

Q R = Q ext + Q flame -CHF b , where
Q flame is the additional heat flux supplied by the flame and CHF b ≈ σ T b 4 is the critical heat flux for burning in terms of the surface burning temperature T b and the Boltzmann radiation consant σ. It has been shown that T b ≈ T p where T p is the pyrolysis temperature measured in laboratory thermal analysis experiments using small samples and constant heating rates [START_REF] Lyon | The Role of Thermal Decomposition Kinetics in the Burning of Polymers[END_REF]. Thus, polymers with high pyrolysis temperatures reradiate more of the incident heat flux from the heater and flame back to the surroundings, and the net heat flux that drives the burning process is reduced accordingly.

The Arrhenius rate constant, k(T) = A exp[-E a /RT] is a reasonable descriptor of the temperature dependence of the rate of polymer thermal decomposition. The kinetic parameter A (s -1 ) represents the frequency of chemical bond breaking reactions in the polymer at temperature T while the activation energy E a represents the thermal energy barrier that must be overcome to break the chemical bonds and produce fuel gases. It has recently been demonstrated for a range of common polymers that the thermal decomposition temperature or peak pyrolysis rate temperature (T p ) measured in thermal analysis experiments such as pyrolysis combustion flow calorimetry (PCFC) [START_REF] Lyon | A Microscale Combustion Calorimeter[END_REF] or thermogravimetric analysis (TGA) has a large effect on ignition and burning [22][23]. Equation 3 is a derived result [START_REF] Lyon | The Role of Thermal Decomposition Kinetics in the Burning of Polymers[END_REF] that shows that the peak pyrolysis temperature in constant heating rate experiments such as TGA or PCFC is defined by an activation energy (E a ) and Arrhenius factor (A) that are not independent:

ܶ  = ܧ ܽ ܴ ln ܣ ݇  ൨ [Equation 3]
In Equation 3, R is the gas constant and

k p = k(T p ) = Aexp[-E a /RT p ] = βE a /R T p
2 is the value of the kinetic rate constant at T p measured for a milligram-size sample at a constant rate of temperature rise β = dT/dt derived from a semi-exact solution of the Arrhenius temperature integral [START_REF] Lyon | The Role of Thermal Decomposition Kinetics in the Burning of Polymers[END_REF]. Since surface heating rates of polymers burning in a cone calorimeter at 50 kW/m 2 external heat flux are comparable to those used to determine A and E a in PCFC or thermogravimetric analyses [START_REF] Lyon | The Role of Thermal Decomposition Kinetics in the Burning of Polymers[END_REF], β ≈ 1 K/s, and since R T p

2 /E a ≈ 20K (typically), a constant value, k p = (1 K s -1 )/(20K) = 0.05 s -1 was used in Equation 3 to calculate the pyrolysis temperatures in Table 1 for A and E a used in ThermaKin (also shown in Table 1). It has been proposed that uncertainty in the Arrhenius parameters manifests itself as uncertainty in modelling the fire response of polymers [START_REF] Lyon | The Role of Thermal Decomposition Kinetics in the Burning of Polymers[END_REF]. The processes modelled by ThermaKin have been summarised in Figure 1. For this study, radiant heat from above the sample is absorbed, emitted or reflected, and the condensed phase heat transfer process is modelled through the solid. The resulting temperature increases drives endothermic decomposition processes, leading to the gasification of volatile fuel components. When a critical mass flux for ignition is reached, ignition will occur, and the incident radiant flux is augmented by radiation from the flame. Thereafter, quasi-steady state conditions pertain, until the sample is so thin that it has no more capacity to absorb heat, and the rate of pyrolysis increases. This present study utilises ThermaKin as a means of relating the physical properties of a material to its HRR history in a cone calorimeter.

Modelling

The effect of physical properties on the fire behaviour of non-charring polymers was investigated using a one-dimensional pyrolysis model, ThermaKin; a complete description of the model's mathematical formulation and numerical algorithms has already been reported [START_REF] Stoliarov | Thermo-Kinetic Model of Burning[END_REF].

A sensitivity analysis has been carried out by the model's developers to determine the relative importance of individual properties [START_REF] Stoliarov | The Effect of Variation in Polymer Properties on the Rate of Burning[END_REF]. Within ThermaKin, changes to the fuel are accounted for as a change in the component. Each component is characterised by its physical state, density, heat capacity, thermal conductivity, gas transfer coefficient, emissivity and absorption coefficient. The model ignores changes in thermal conductivity resulting from changes in melt flow behaviour on heating. The chemical processes occurring are characterised through the reaction's activation energy and Arrhenius factor. The energy balance assumes that radiant heat may be absorbed or reflected by the sample, and then either conducted through it resulting in a localised temperature increase, or re-radiated from the surface. Higher temperatures will result in gasification forming vapour phase fuel which can b ignited and then burns, transferring some radiant heat back to the sample. From the sensitivity analysis [START_REF] Stoliarov | The Effect of Variation in Polymer Properties on the Rate of Burning[END_REF], it was established that for most synthetic polymers, the thermal, optical and chemical properties varied only within limited ranges. To demonstrate this, the lower, average and upper boundaries of each parameter were determined from experimental techniques and reported values [START_REF] Stoliarov | The Effect of Variation in Polymer Properties on the Rate of Burning[END_REF]. These are shown in Table 1. The study also determined that some parameters had a greater influence on time to mass loss, peak mass loss and average mass loss rate and subsequently time to ignition, peak heat release rate and average HRR when determined using the The ThermaKin model was set up as follows. The material (assumed to be an average polymer) was represented at all times by the average values -the property values and reaction parameters were varied within the boundaries of observed behaviour to give a lower, average and upper scenario for each. The material was specified to decompose via a first order reaction, defined by the Arrhenius parameters and heat of decomposition. It was assumed that the reaction had only gaseous products, which were instantly removed from the condensed phase, leaving no char residue.

The burning was modelled to represent a scenario where the top surface of a one-dimensional material is exposed to a constant radiative heat flux. Natural convection was simulated by subjecting the top layer to a convection coefficient of 10 W m -2 K -1 with an ambient temperature of 300 K. The mass flux for ignition was set at 1 x 10 -3 kg m -2 s -1 corresponding to the onset of piloted ignition for polymers under convective conditions [START_REF] Lyon | Criteria for Piloted Ignition of Combustible Solids[END_REF]. The ignition process resulted in an additional 15 kW m -2 heat flux from the flame. The bottom surface remained insulated (to minimise heat loss effects) with a layer of glass wool set at a thickness of 1.5 x 10 -2 m, following the standard cone calorimeter set up [3]. The thickness of the material was set at 2.5 x 10 -3 m and the applied heat flux at 50 kW m -2 .

The mass loss rate for each scenario was obtained by numerical differentiation of the material mass versus time data -outputs of the ThermaKin model. To obtain HRR for a 'generic' polymer, the mass loss rate was multiplied by the heat of combustion for PMMA (a single decomposition product polymer which leaves no solid residue) obtained from literature to be 24.99 KJ g -1 [START_REF] Walters | Heats of Combustion of High Temperature Polymers[END_REF]. The choice of PMMA was somewhat arbitrary, but for the purpose of comparison, a single ... value must be used.

The values for the net heat of combustion were also varied, as shown in Table 1. Polyethylene (PE) was selected for the upper boundary and polyoxymethylene (POM) for the lower boundary. The physical properties input into the ThermaKin model are density, thermal conductivity and heat capacity. The effect of density on HRR histories is shown in Figure 2. The ThermaKin input is actually based on a defined volume and density, from which mass is implicit. For the purpose of this study, it is easier to compare samples of constant mass, so three samples of different density have been compared at constant mass. It is important to note that all the polymers whose parameters were identified were solids, not foams and other expanded materials containing voids of air or gas; which would have a more profound effect. The HRR corresponding to the lower limit for density gives a shorter time to ignition and the upper limit gives the longest time to ignition. The low density materials has a high surface temperature, reaching the critical pyrolysis flux first, but also a larger thermal gradient through its bulk, and hence requires more heat to get the remainder of the sample up to the pyrolysis temperature. This can be seen in the longer time to reach the peak HRR than the average or high-density polymer. The only effects of density (ρ) are on the time to ignition and the time to peak HRR. The effect of thermal conductivity on HRR histories is shown in Figure 3. Thermal conductivity relates to a material's ability to conduct heat and transfer energy through its bulk for scenarios where a temperature gradient exists. A low thermal conductivity indicates less ability to dissipate heat through the material. As a result, more heat remains concentrated at the surface and therefore the material ignites earlier. Again, the HRR history corresponding to the first sample to ignite, in this case also having the lower limit for thermal conductivity takes longer to reach its peak value in comparison to the HRR history for the upper limit of thermal conductivity. A material with a higher thermal conductivity allows heat to be dissipated effectively through the bulk and therefore, due to the build up of heat within the system, the material does ignite it will burn more quickly. A material with a low thermal conductivity is more likely to show thermally thick burning behaviour, indicating that the bulk of the fuel acts as an almost infinite heat sink during the early stages of burning. Similarly, a material with a greater thermal conductivity is more likely to show thermally thin burning. As a result, this gives a longer time to peak HRR as each layer ignites and burns as a separate entity. The effect of heat capacity on HRR histories is shown in Figure 4. The heat capacity of a material refers to the amount of energy required to raise the temperature of one kilogram of polymer by 1 K. A lower heat capacity indicates that less energy is required to raise the temperature of the surface to ignition temperature. A low heat capacity material will have a low thermal inertia (kρc) and hence the upper layers will reach the critical surface temperature for ignition more quickly than a material with higher heat capacity (or greater thermal inertia). At a constant external heat flux, this is governed by the time of exposure, which is balanced by the energy required for gasification (the heat of decomposition). The peak HRR and the time of test remain constant. The overall effect of the three parameters, taking the lowest, average and highest of each as the effect of thermal inertia on HRR histories is shown in Figure 5. This shows greater differences in both the shape and position of the respective HRR histories. The very short time to ignition for the lower thermal inertia sample, coupled with the poorer heat transfer results in ignition occurring in a material which is cooler underneath. As burning penetrates through the material, less heat is required to bring the remainder to pyrolysis temperature, so more is available for pyrolysis. This is seen as a steady increase in the HRR history. In contrast, the high thermal inertia sample a long ignition delay time, and a very steady, high HRR. The optical properties employed by the ThermaKin model are reflectivity and absorption coefficient. The process of absorption of radiant heat depends on the optical properties of the sample. A mirrorlike finish would not absorb any radiation, and would never ignite. Absorption is essentially a molecular process where a photon of radiation results in the excitation of a specific atom or molecule. The infrared radiation from the cone heater at 50 kW m -2 results in vibration and excitation which eventually relaxes down to increased thermal energy, or a higher temperature. The effect of reflectivity on HRR histories is shown in Figure 6. Reflectivity refers to the amount of infrared radiation that is reflected from the material's surface or essentially, how the surface interacts with external heat fluxes from the cone heater and then from both the cone heater and the flame, subsequent to ignition. The HRR history corresponding to the upper limit for reflectivity shows a longer time to ignition. This is due to the fact that more energy is reflected and as a result, less energy is absorbed by the surface. In a system with a constant external heat flux, this also results in a lower peak of HRR, as less heat is available for decomposition and volatilisation, so the overall process is slower. In contrast the low reflectivity surface shows a shorter time to ignition, higher peak HRR and an overall shorter burning time. Given that these predictions are from unadulterated polymers, the addition of an additive or a change in the sample preparation (surface roughening) could have a profound effect on the HRR history just because of changes in the reflectivity. The effect of absorption coefficient on HRR histories is shown in Figure 7. A material's absorption coefficient refers to the amount of energy absorbed by the material. Certain polar bonds are good absorbers of infrared radiation, and if these are present, most of the heat will be absorbed by the surface layers. The absence of such absorbers will allow more heat to penetrate within the bulk of the polymer, and possibly even pass straight through it. The surface layers of a material with a large absorption coefficient will heat up much quicker than those of a material with a lower absorption coefficient. However, the non-reflected radiation which penetrates the sample without being absorbed will cause heating in lower layers. Thus, highly absorbing materials have an early time to ignition however the lower absorbing materials will have the shortest time to peak HRR. The HRR corresponding to the lower limit for absorption shows a longer time to ignition and a more progressive decrease from its peak of HRR to 0 kW m -2 . This is because an increasing proportion of the radiation passes through samples as the thickness decreases. In practice, the absorption of polymers is likely to change during decomposition as the bonding changes, increasing the number of absorbing centres, so this type of behaviour is rarely observed. The chemical properties employed by the ThermaKin model are heat of decomposition, activation energy, Arrhenius factor and heat of combustion. The heat of decomposition, or heat of gasification describes the endothermic transition from condensed to gas phase fuel. It is analogous to the latent heat of vaporisation of a liquid, but prior to vaporisation the polymer must first breakdown into smaller, volatile fragments. The effect of heat of decomposition on HRR histories is shown in Figure 8. A simple energy balance exists between the heat of decomposition (or the heat required to produce the fuel) and the heat released by it. In thermally thick burning, at a constant heat flux, the rate of heat release is proportional to the rate of pyrolysis. If more heat is required to pyrolyse the fuel this will give a slower rate of pyrolysis and a lower, steadier HRR. The large differences in the shape of the heat release histories show the importance of the heat of decomposition parameter, which could also be described as the thermal stability of the polymer The values for T p in these simulations: 500, 700 and 900 K correspond to the variation of polymer properties used for the sensitivity analysis, although in this case, values are rounded off to the nearest 100 K. The effect of pyrolysis temperature expressed through the A/E a pair in Table 1 on the HRR histories calculated by ThermaKin is shown in Figure 9. The influence of high thermal stability (high T p ) on the burning behaviour is obvious in both the time to ignition and the burning rate after ignition. High thermal stability increases the time to ignition because the critical mass flux is not reached until the surface approaches T p . High thermal stability also reduces the burning (heat release) rate because a larger fraction of the incident energy is reradiated from the surface at high T p and the net heat flux is reduced accordingly. In the ThermaKin simulations, Q ext = 50 kW m -2 , Q flame = 15 kW m -2 and assuming CHF b = σ T p 4 the net heat flux would be

Q R = (1 -r)(Q ext + Q flame -CHF b ) = 54,
45 and 25 kW m -2 for T p = 500, 700 and 900 K, respectively. The maximum HRR in Figure 9 decreases in rough proportion to these CHF b because less energy is available to drive the burning process as the thermal stability of the polymer increases. The effect of heat of combustion on HRR histories is shown in Figure 10. The heat of combustion depends on the elemental composition of the material and its combustion efficiency. In many cases, the heat of combustion reported in literature is determined under rather artificial conditions of 25 atmospheres of pure oxygen in a bomb calorimeter. The sample with the highest heat of combustion, polyethylene (PE), contains just carbon and hydrogen, which release heat to form carbon dioxide and water. Lower heats of combustion are obtained from samples containing oxygen or other non-combustible elements. It is somewhat surprising to see that the heat of combustion has no influence on the modified HRR histories on the time of burning, and only affects the profile of HRR.

Conclusions

Cone calorimetry is probably the single most important tool for the assessment of flammability behaviours. A number of parameters govern the burning behaviour, which could be subdivided into: ignitability and flame spread (a sources of repeated ignitions); HRR (the larger and earlier the peak of HRR the faster the rate of fire growth); and smoke and toxic gas production. Adequate screening is essential for the development of fire safe materials -this requires proper understanding of the results -in order to correlate the burning behaviour to the underlying physical properties. By separating the effects of physical properties, the more specific chemical effects of fire retardants can then be identified.

The Thermakin modelling of cone calorimeter behaviour presented here predicts the individual effect of the optical, thermal and chemical properties of the fuel. Additive fire retardants, usually added at loadings of 10-70% have a profound effect on the physical properties, changing the absorption, transfer of heat and decomposition behaviour, as well as inhibiting the gas phase free radical processes. This study examines the effect of typical variations in polymer properties on the HRR history.

The density of most polymers lies within a fairly narrow range and has a small influence on the time to ignition and time to peak of HRR, but in these predictions has no effect on the actual height of the peak. The thermal conductivity varies by a factor of 4, and materials with the lowest thermal conductivity are predicted to have half the time to ignition and double the time to peak of heat release, as may be expected. The heat capacity of polymers varies over a narrow range, affecting the time to ignition but not the overall burning time or peak of heat release. The addition of most additives, particularly inorganic materials, will change the thermal inertia (the product of density, thermal conductivity and heat capacity) and this will influence the ignition and burning behaviour.

The optical characteristics of polymers, absorption and reflectivity vary by factors of 8 and 4 respectively, though surface coatings and additives could produce much larger differences. The absorption coefficient affects the shape of the heat release history, showing sharply decreasing absorbance by the lowest absorber as the sample becomes thinner. In practice, the absorbance of infrared radiation is likely to change during decomposition, with an increase in absorbing centres as the chemical structure becomes more heterogeneous, tending to mitigate this effect. The surface reflectivity affects both the time to ignition and the peak HRR, but not the total heat released (as may be expected), so highly reflective polymers take longer to ignite and burn more slowly.

The chemical composition of the polymer exerts the greatest influence on its burning behaviour. The heat of decomposition (or energy required for gasification) varies over more than two orders of magnitude, has no effect on the time to ignition, but a strong effect on the peak of heat release and the burning time. The pyrolysis temperature, which relates to the rate of gasification (or fuel production) has only been varied from 500-900K but has the strongest influence on both the time to ignition and the peak of HRR. This is the single most important parameter for predicting the burning behaviour of a polymer. Finally the heat of combustion of the gas phase fuel affects the peak HRR, but not the time to ignition (as may be expected) or the burning time (which is more surprising since presumably a lower heat of combustion is accompanied by a decrease in thermal feedback).

The effects of flame retardant chemicals on the burning rate of polymers are incorporated into ThermaKin by including reactions that generate products whose properties have the desired (observed) effect. For example, intumescence can be introduced into the model as a reaction product (component) having properties (density, thermal conductivity, heat capacity, mass transport coefficient, in-depth absorption of radiation, etc) that may differ from the polymer but which reproduce the desired (observed) effect on the burning rate. The trigger for reaction from one component to another is temperature. Thus, within the model of the burning material in the cone calorimeter, an intumescent foam could form first in the uppermost part of the material, and then the foam/non-foam interface will progressively penetrate downwards into the bulk of the material as the temperature rises. The only limit to the number of intermediate products that may be included is the availability of physical property data for them. Likewise, the optical properties of the flame retardant may be represented by an inert component with an emissivity that is higher (such as char or carbon nanotubes) or lower (such as titanium dioxide (TiO 2 ), mica or glass fibres) than the emissivity of the matrix polymer so as to absorb or reflect more of the incident radiant energy, respectively, at the surface or in depth. As the polymer matrix recedes, the solid component (char, TiO 2 etc.) will concentrate at the surface and change its emissivity. The effect of gas phase active flame retardants on burning rate can be captured empirically by adjusting (or measuring) the effective heat of combustion of gaseous products and/or the flame heat flux back to the material surface, which will be a function of the soot concentration and flame temperature.

Overall, the Thermakin model provides a useful tool for understanding how the physical properties of polymers influence their burning behaviour. The ability to isolate individual physical properties shows great potential to optimise formulations on a microscale, prior to screening on a cone calorimetry scale.

Figure 1 .

 1 Figure 1. Schematic of processes occurring in the cone calorimeter, as modelled by Thermakin

Figure 2 .

 2 Figure 2. Effect of Density (with Constant Mass) on Cone Calorimeter HRR Histories

Figure 3 .

 3 Figure 3. Effect of Thermal Conductivity on Cone Calorimeter HRR Histories

Figure 4 .

 4 Figure 4. Effect of Heat Capacity on Cone Calorimeter HRR Histories

Figure 5 .

 5 Figure 5. Effect of Thermal Inertia on Cone Calorimeter HRR Histories

Figure 6 .

 6 Figure 6. Effect of Reflectivity on Cone Calorimeter HRR Histories

Figure 7 .

 7 Figure 7. Effect of Absorption Coefficient on Cone Calorimeter HRR Histories

Figure 8 .

 8 Figure 8. Effect of Heat of Decomposition on Cone Calorimeter HRR Histories

Figure 9 .

 9 Figure 9. Effect of Pyrolysis Temperature on Cone Calorimeter HRR Histories

Figure 10 .

 10 Figure 10. Effect of Heat of Combustion on Cone Calorimeter HRR Histories

Table 1 .

 1 material's effective heat of combustion. The authors conclude that knowledge of the activation energy, Arrhenius factor, heat of decomposition and char yield was central to the accuracy of simulations, with knowledge of absorption coefficient and reflectivity also being fairly important. Sensitivity Analysis -Values for Lower, Average and Upper Parameter Boundaries for Common Polymers

	Parameter	Lower	Average	Upper
	Density (ρ)/ kg/m 3	830	1300	1830
	Thermal Conductivity (K)/ W m -1 K -1	0.11	0.24	0.42
	Heat Capacity (c)/ J kg -1 K -1	1700	2300	2900
	Reflectivity (r)	0.06	0.12	0.25
	Absorption Coefficient (α)/ m -1	1100	3800	9000
	Heat of Decomposition (h dec )/ kJ kg -1	0.1	1.3	2.5
	Heat of Combustion(h c )/ MJ kg -1	15.93	24.99	44.60
	Pre-Exponential Factor (A) /s -1 Activation Energy (E a )/J mol -1 F r o Pyrolysis Temperature (T p )/ K	1.0 x 10 8 8.90 x 10 4 500	1.0 x 10 14 2.05 x 10 5 700	1.0 x 10 20 3.67 x 10 5 900
	P
	e
	e r
		R
		e
		v i e
			w

John Wiley & SonsPolymers for Advanced Technologies

Acknowledgements

One of us (PP) would like to thank the Engineering and Physical Science Research Council (EPSRC) for the award of a CASE PhD Studentship (No. 07001444).


Enrique Regidor
email: enriqueregidor@hotmail.com

Paloma Astasio

Paloma Ortega

) David Martínez

Elisa Calle

Luis De La Fuente

Healthy and unhealthy migrant effect on the mortality of immigrants from wealthy countries residing in Spain

Keywords: immigrants, wealthy countries, unhealthy migrant effect, healthy migrant effect, early retirement, Spain

This study attempts to identify the possible existence of a healthy migrant effect and an unhealthy migrant effect on the mortality of immigrants from wealthy countries who move to Spain. Immigrants aged 35-64 years from France, Germany, Great Britain and 16 other wealthy OECD countries who resided in Spain were compared with respect to: 1) mortality from cancer, cardiovascular disease, and all other diseases and 2) employment status, duration of residence, and educational level, in two geographic areas: the -preferred destination area‖ --the Mediterranean coast, Balearic Islands and Canary Islands --and the rest of Spain. In general, cancer mortality was lower and mortality from cardiovascular disease was higher in immigrants who resided in the preferred destination area than in their countries of origin and than in immigrants who resided in the rest of Spain. Immigrants in the preferred destination area had a higher percentage of retired persons, longer time of residence and a lower percentage of persons with university education. The largest differences between the two areas in cardiovascular and all-disease mortality and in the frequency of the aforementioned sociodemographic characteristics were observed in British immigrants and those from the 16 OECD countries.

Possible explanations for these findings are suggested which are compatible with the presence of an unhealthy and/or healthy immigrant bias in the two areas.

Introduction

Studies have shown important differences in the mortality of immigrants with respect to mortality in their destination countries, differences that vary in magnitude and direction depending on the immigrants' place of origin and cause of death. [START_REF] Wild | Cross sectional analysis of mortality by country of birth in England and Wales[END_REF][START_REF] Singh | All-cause and cause-specific mortality of immigrants and native born in the United States[END_REF][START_REF] Bos | Ethnic inequalities in age-and cause-specific mortality in The Netherlands[END_REF][START_REF] Deboosere | Adult migrant mortality advantage in Belgium: evidence using census and register data[END_REF][START_REF] Gadd | Are there differences in all-cause and coronary heart disease mortality between immigrants in Sweden and in their country of birth? A follow-up study of total populations[END_REF] Various explanations for these results have been suggested, among them, characteristics of place or environment of origin, differences in genetic susceptibility, healthy migrant effect, time of residence in the host country, unhealthy re-migrant effect, differences in a wide variety of behavioural, cultural, socioeconomic and psychosocial characteristics, and adoption of the lifestyle in the host country. [START_REF] Abraido-Lanza | The Latino mortality paradox: a test of the -salmon bias‖ and healthy migrant hypothesis[END_REF][START_REF] Razum | The -healthy migrant effect‖not merely a fallacy of inaccurate denominator figures[END_REF][START_REF] Simth | Ethnic inequalities in health: a review of UK epidemiological evidence[END_REF][START_REF] Razum | Commentary: of salmon and time travellers -musing on the mystery of migrant mortality[END_REF] One aspect not considered to date in the study of immigrant mortality is the possibility of selfselection of sick individuals into migration. In wealthy countries many persons retire from the labour force because they develop diseases that make it impossible to work. [START_REF] O'reilly | Intra-European migration and the mobility-enclosure dialectic[END_REF][START_REF] Warnes | The diversity and welfare of older migrants in Europe[END_REF] Since this disability leads to reduced income, some people move to other countries where the cost of living is lower and the quality of life is similar to that of their country of origin. In these cases, any attempt to explain the mortality of immigrants from wealthy countries should consider the possibility that they are less healthy than the population in their country of origin.

Our hypothesis of an unhealthy migrant effect was motivated by our previous comparison of the mortality of immigrants residing in Spain with mortality in the Spanish population. We found that, in the early years of the present decade, mortality in immigrants aged 20-64 from western countries was 10% higher than mortality in the Spanish population (unpublished data). In contrast, a previous study conducted only in residents of the province of Madrid, located in the geographic centre of Spain, found that mortality in immigrants from western countries was 20% lower than mortality in the Madrid-born population. [START_REF] Regidor | Heterogeneity in cause-specific mortality according to birthplace in immigrant men residing in Madrid, Spain[END_REF] We then set ourselves the challenge of finding explanations for this paradox.

Spain has traditionally been one of the countries that offers the most attractions to persons in wealthy countries considering a change of residence. [START_REF] Rodríguez | European retirees on the Costa del Sol: a cross-national comparison[END_REF] According to a recent study, Spain is the European country with the highest European Quality of Life Index due to low taxes, cheaper essential goods and a larger number of hours of sunshine. [START_REF] Casado-Díaz | Northern European retired residents in nine southern European areas: characteristics motivations and adjustment[END_REF][START_REF]UK and Ireland rank bottom in Quality of Life Index www.uswitch[END_REF] However, the primary attraction is not just any location in Spain, but is concentrated in the provinces of the Mediterranean coast, the Balearic Islands and the Canary Islands. 16 Thus, the previously described finding of higher mortality in immigrants from wealthy countries than in the Spanish population may be due to higher mortality in immigrants who live in these areas than in immigrants who reside in the rest of Spain.

In this study, we attempt to obtain evidence to support the existence of a healthy migrant effect and an unhealthy migrant effect on mortality in immigrants from wealthy countries who move to Spain. For this purpose, we compare mortality in immigrants from wealthy countries who reside in two geographic areas --the preferred destination area and the rest of Spain --with mortality in the countries of origin, and we compare the frequency of various characteristics in immigrants who reside in each of these two areas of Spain.

Methods

Sources of data

Mortality data by cause of death were obtained from the mortality register for 2002 to 2006, and population data were obtained from the population register for the same period. The data files from both registers were provided to the investigators by the National Institute of Statistics (NIS), which collects data on all deaths and demographic changes in the population. Information identifying individuals was eliminated to preserve confidentiality.

The mortality register includes province of residence, age, sex and birthplace of persons who die in Spain. These characteristics are collected by funeral directors based on information provided by relatives of the deceased. The records also include the underlying cause of death certified by a physician and coded according to the International Classification of Diseases, 10 th Revision. In the present study, external causes of death were excluded from the analysis because these are acute events. Besides deaths from all diseases, we also analysed deaths from cancer, from cardiovascular disease, and from all other diseases.

The population register includes information on all persons living in each municipality. The City Council is responsible for maintaining and updating the population register in each municipality.

Each City Council reports to the NIS any changes in the data contained in this register. The NIS produces an annual file from the population register with reference to 1 January of each year, based on which official estimates are made of the population residing in Spain. These files were used in the present investigation. The variables contained in this file, obtained from the information provided by each resident, are province of residence, age, sex and birthplace.

This study was restricted to immigrants from OECD countries that had a higher per capita income than that of Spain between 1996 and 2006. Only immigrants aged 35-64 years were included. The upper limit was established because, based on our previous findings, our hypothesis refers to subjects of working age. And the lower limit was chosen because the probability of developing chronic diseases is very low in persons under 35 years of age.

Age, sex and birthplace were recorded in the same way in both files. Age was categorised in 10year age groups, and birthplace was defined as follows: based on the population size, only persons from France, Germany and the United Kingdom could be differentiated in separate groups, since they represent 70% of the immigrants from wealthy countries residing in Spain.

The rest of the immigrants came from the following 16 countries: Switzerland, Italy, Belgium, The Netherlands, the United States of America, Sweden, Norway, Ireland, Denmark, Austria, Finland, Canada, Japan, Australia, Luxembourg and Iceland. Thus, the immigrants were grouped into four categories: those born in France, Germany, the United Kingdom and other OECD countries.

The National Survey of Immigrants, conducted between November 2006 and February 2007, was used to determine employment status, time of residence in Spain, and educational level of the immigrant population. The file with the survey microdata was downloaded from the NIS website. 17 The sample was designed so that the estimates would be representative of the immigrants from different countries, including the countries that were part of the European Union before its expansion in 2004, plus Switzerland and Norway, as well as the two countries of North America. For this purpose the sample was designed taking the 2006 population register as the reference. The sampling framework for this survey was the population aged 16 years and over born outside of Spain, who had resided in Spain for more than one year. The sample was selected in three stages: the first-stage units were the census sections stratified by municipal population size; the second-stage units were family households; and the third-stage units were household residents born outside of Spain. The non-response rate was 13%, and the number of immigrants interviewed was 15,645. Participants in this survey were asked about their employment status in the week before the interview, the year they arrived in Spain, and the highest educational level completed.

Statistical analysis

A separate analysis was made for each of the two study areas: the preferred destination area and the rest of Spain (figure 1). The preferred destination area consisted of 12 of the 50 Spanish provinces, where 65% of the immigrants included in this study resided. These provinces were: Girona, Tarragona, Castellón, Alicante, Murcia, Almería, Granada, Málaga and Cádiz (located on the Mediterranean coast), The Balearic Islands, and the two provinces of the Canary Islands (Santa Cruz de Tenerife and Las Palmas). Although Barcelona and Valencia are located on the Mediterranean coast, they were included in the -rest of Spain‖ area because they have a high level of commercial and industrial development, and it is likely many immigrants live there for work-related reasons.

Men and women were analysed separately. We first calculated the age-adjusted mortality rate per 100,000 person-years in each immigrant group for the selected diseases. The age distribution of the Spanish population was used as the standard population. The mortality rate in the two study areas was then compared for each immigrant group. Specifically, we estimated the mortality rate ratios in residents of the preferred destination area with respect to residents in the rest of Spain. The confidence intervals for the mortality rate ratios were calculated from the variance estimated by the Mantel-Haenszel method for person-time data. [START_REF] Greenland | Introduction to stratified analysis[END_REF] To determine whether the mortality of immigrants residing in each study area is different from mortality in the countries of origin, we compared mortality from different causes of death in French, German and British citizens residing in the preferred destination area and in the rest of Spain with mortality in France, Germany and the United Kingdom in the period 2002-2006.

Information on the deaths in those years and on the population aged 35-64 years from those countries was obtained from the mortality database of the World Health Organization. [START_REF]Health statistics and health information systems[END_REF] Mortality in the countries of origin of immigrants from "other OECD countries" was not estimated given that they comprise a heterogeneous group of subjects from many different countries.

In each immigrant group and destination area we also calculated the percentage of retired persons, the mean time of residence in Spain, and the percentage of persons with university education. Finally, we tested whether the magnitude of the mortality differences in the immigrant population with respect to mortality in the Spanish population varied depending on the area of residence. In each study area, we estimated the age-and sex-adjusted mortality rate ratios by birthplace, taking the mortality rate of persons born in Spain as the reference.

Results

Table 1 shows the population and number of deaths in immigrants residing in each study area, by place of birth. Most of the immigrants lived in the preferred destination area, except for the French immigrants who primarily resided in the rest of Spain. The mean age of immigrants in the preferred destination area was higher than in the rest of Spain: about 1.5 years higher in the French immigrants and 5-8 years in the other immigrant groups.

Table 2 presents the mortality rate ratios from the causes of death analysed in the preferred destination area with respect to the rest of Spain. Significant differences between the two areas for mortality from all diseases were seen only in British immigrants, whose mortality rate ratio was 1.44 (1.18 to 1.75) in men and 1.36 (1.04 to 1.77) in women, and in immigrants from other OECD countries, whose mortality rate ratio was 1.43 (1.24 to 1.64) in men and 1.52 (1.25 to 1.84) in women. In all immigrant groups the cancer mortality rate was lower in the preferred destination area than in the rest of Spain. In contrast, the cardiovascular disease mortality rate was higher in the preferred destination area than in the rest of Spain, except in French women. The largest differences were observed in British immigrants and those in other OECD countries, for whom the cardiovascular disease mortality rate ratio was 2.56 (95% confidence interval 1.71 to 3.83) and 2.56 (1.92 to 3.41) in men, and 4.32 (1.90 to 9.84) and 3.51 (2.08 to 5.91) in women, respectively. The highest mortality rate ratios for the rest of diseases were also observed in British immigrants and those from other OECD countries.

Table 3 presents the age-adjusted mortality rates and the mortality rate ratios for the causes of death analysed in each area and in the country of origin. In general, mortality rate ratios for cardiovascular diseases were significantly higher and mortality rate ratios for cancer were significantly lower in immigrants residing in the preferred destination area than in the country of origin. Likewise, the mortality rate ratios for cardiovascular diseases in British immigrants residing in the rest of Spain were lower than in the UK. With respect to mortality from the rest of diseases, French immigrants had lower mortality in the two study areas than in France, whereas British men residing in the preferred destination area had a higher mortality rate ratio than in the UK.

Immigrants in the preferred destination area had a higher percentage of retired persons and a lower percentage of persons with university educationexcept in German womenthan those in the rest of Spain (table 4). The largest differences between the two areas were seen in British immigrants and those from other OECD countries. Immigrants residing in the preferred destination area had lived in Spain a shorter time, with the largest difference between the two areas seen in British immigrants.

The mortality rate ratios in each immigrant group with respect to the Spanish population are presented in table 5. In all four immigrant groups the magnitude of the mortality rate ratios for cancer and cardiovascular disease was different in the two study areas; in contrast, the magnitude of the mortality rate ratios for the rest of diseases in the two areas was only clearly different in British immigrants and those from other OECD countries. The mortality rate ratios for all diseases were similar for French and German immigrants in the two study areas, and did not differ significantly with respect to the Spanish population. However, in British immigrants and those from other OECD countries these mortality rate ratios were, respectively, 1.14 (1.04 to 1.25) and 1.19 (1.07 to 1.32) in the preferred destination area, and 0.86 (0.64 to 1.16) and 0.86 (0.70 to 1.05) in the rest of Spain.

Discussion

Our comparison of mortality in immigrants from wealthy countries in the two study areas with mortality in their countries of origin yields the following three relevant findings: a) mortality from cancer is lower in the preferred destination area and in the countries of origin than in the rest of Spain; b) except in French women, mortality from cardiovascular disease is higher in the preferred destination area and in the countries of origin than in the rest of Spain; and c) in British immigrants and those from other OECD countries mortality from the rest of diseases is higher in the preferred destination area than in the rest of Spain. Because the largest differences between the two areas of Spain in mortality from cardiovascular disease and from the rest of diseases are seen in British immigrants and those from other OECD countries, mortality for all diseases in these immigrants is between 36% and 52% higher in the preferred destination area than in the rest of Spain. The comparison of immigrant mortality with mortality in the Spanish population also reflects these findings.

Several studies have shown that early retirement is associated with increased mortality. [START_REF] Morris | Loss of employment and mortality[END_REF][START_REF] Quaade | Mortality in relation to early retirement in Denmark: a population-based study[END_REF][START_REF] Tsai | Age at retirement and long term survival of an industrial population: prospective cohort study[END_REF] Some authors have noted that mortality from cardiovascular disease associated with early retirement may also be attributed to changes after retirementreduced income, abandoning healthy habits, or psychological consequencesthat increase the risk of developing these diseases. [START_REF] Bamia | Age at retirement and mortality in a general population sample: the Greek EPIC Study[END_REF] We found that immigrants residing in the preferred destination area had a higher proportion of retired people than immigrants residing in the rest of Spain. This could explain why mortality from cardiovascular disease was higher in the preferred destination area in most immigrants than in the rest of Spain and in their countries of origin.

Increased mortality associated with early retirement has also been attributed to illness and poor health before retirement. [START_REF] Morris | Loss of employment and mortality[END_REF][START_REF] Quaade | Mortality in relation to early retirement in Denmark: a population-based study[END_REF][START_REF] Tsai | Age at retirement and long term survival of an industrial population: prospective cohort study[END_REF][START_REF] Karpansalo | Perceived health as a predictor of early retirement[END_REF][START_REF] Mein | Predictor of early retirement in British civil servants[END_REF] Both cardiovascular disease and the rest of diseasesdigestive, respiratory, metabolic, renalthat constitute the leading causes of death are chronic diseases that generally have a prolonged natural history before retirement. Thus, the higher mortality from these causes of death in the preferred destination area suggests poor health before retirement and, consequently, a possible unhealthy immigrant bias in mortality in this area. The fact that the elevated mortality rate ratio for the rest of diseases by area of residence was only observed in British immigrants and those from other OECD countries may be to the large proportion of retired persons in these two immigrant groups in the preferred destination area.

In any case, it can not be ruled out that the elevated differences in mortality from cardiovascular disease and the rest of diseases in British immigrants and those from other OECD countries who reside in the preferred destination area and in the rest of Spain may be due to the presence of a healthy immigrant bias in the rest of Spain. Mortality from both causes of death in the United Kingdom was found to be much higher than that seen in British immigrants in the rest of Spain.

Likewise, we also can not rule out the possible influence of immigrants' educational level on the mortality from cardiovascular disease and the rest of diseases observed in British immigrants and those from other OECD countries who reside in the rest of Spain. Given the strong association between educational level and mortality from both causes of death, [START_REF] Kaplan | Socioeconomic factors and cardiovascular disease: a review of the literature[END_REF][START_REF] Smith | Education and occupational social class: which is the more important indicator of mortality risk?[END_REF][START_REF] Huisman | Educational inequalities in cause-specific mortality in middle-aged and older men and women in eight western European populations[END_REF] this low mortality may be due to the high educational level in British immigrants and those from other OECD countries, about half of whom had university education. These may be persons who have moved to Spain to continue their careers. This fact would help explain the large differences in cardiovascular mortality observed between the two areas in these immigrants.

Cancer also constitutes a collection of chronic processes, but the pattern of cancer mortality is different from the pattern seen in the other diseases. Various explanations can be suggested for this finding. One of them is related with the large proportion of retired people among immigrants residing in the preferred destination area. The probability of moving after early retirement may be lower if a person has been diagnosed with cancer than in the case of another disease. Since the prognosis for cancer is worse than for other diseases, people may wish to die in their country of origin. Another explanation may be that immigrants frequently return to their country of birth if they are diagnosed with cancer, but not if they are diagnosed with another disease. Nevertheless, it is unlikely that some immigrantssuch as the Britishwould have developed cancer after emigrating to the preferred destination area given that, on average, they had lived in Spain for quite a short time. Moreover, this explanation would entail the assumption that immigrants who reside in the rest of Spain do not return to their country of origin as a result of developing cancer.

And, finally, we can not rule out the possible existence of a healthy immigrant bias in relation to this disease in the preferred destination area. However, this would not explain why this bias is not observed in mortality from cancer in immigrants residing in the rest of Spain.

Mortality from the rest of diseases in immigrants from France and Germany was similar in both the preferred destination area and the rest of Spain. However, the reason for this finding was different in these two immigrant groups. Whereas a healthy immigrant bias was observed in the French immigrants, since mortality in both areas was lower than mortality in France, no significant differences in mortality were seen in German immigrants in the two areas in comparison to mortality in Germany.

In interpreting our study findings we must consider the possibility of a numerator/denominator information bias in relation with the country of birth, since the population data and mortality data were taken from different sources. However, the probability of this bias is small since the source of the population data is a continually updated register that receives data from the mortality register. Likewise, mortality rates could have been overestimated since not all immigrants in Spain are necessarily inscribed in the population register, and this overestimate could be greater in the preferred destination area. However, the possible influence of this bias on our findings is negligible, since it would not explain why the pattern of cancer mortality differs from that observed for cardiovascular mortality. On the other hand, neither the mortality register nor the population register includes information about employment status, educational level or time of residence in Spain. However, the effect of these factors on the results could only be evaluated indirectly, from the information provided by the National Survey of Immigrants.

In summary, immigrants from wealthy countries who reside in the preferred destination area, which consists of the Mediterranean coast, the Balearic Islands and the Canary Islands, have lower mortality from cancer and higher mortality from diseases other than cancer than those who reside in the rest of Spain and than in their countries of origin. These findings are compatible with the presence of an unhealthy and healthy immigrant bias in the two areas. We have suggested possible explanations here, but follow-up studies that permit estimation of mortality in both the migrant population and the population that remains in the places of origin will be needed to determine the reasons for these findings.

Figure 1 .

 1 Figure 1. Spanish provinces included in each study area: Preferred destination area and rest of Spain.

THE ATLANTIC OCEAN THE MEDITERRANEAN SEA Sta. Cruz de tenerife Las Palmas de Gran Canaria FRANCE PORTUGAL Madrid Valencia Barcelona

	Birthplace
		Person-years	Mean age		Person-years	Mean age
	Men France	74,721	43.9	201	119,372	42.4		256	Girona
	Germany	150,164	50.8	836	69,906	42.3		164
	United Kingdom	209,561	52.2	1,259	44,805	45.6		110
	Other OECD countries 2	172,701	49.9	901	122,430	44.3	Tarragona 286
	Women								The Balearic Islands
	France	81,364	44.7	106	134,603	Castellón 43.5	157
	Germany	163,477	50.6	403	71,199	43.4		120
	United Kingdom	235,983	51.8	694	41,821	46.6		62
	Other OECD countries 2	176,205	50.2	550	107,455	45.0		148
	1. The Mediterranean coast, The Balearic Islands and The Canary Islands		Alicante
	Murcia 2. Switzerland, Italy, Belgium, The Netherlands, the United States of America, Sweden, Norway, Ireland,
	Denmark, Austria, Finland, Canada, Japan, Australia, Luxembourg and Iceland
					Granada	Almería
			Cádiz	Málaga	ia
	The Canary Islands
			Preferred destination area	Rest of Spain.

Table 1 .

 1 Population and deaths in immigrants aged 35-64 residing in the two study areas, by birthplace.

	Spain

Table 2 .

 2 Number of deaths, and mortality rate ratios with 95% confidence intervals (95% CI) in immigrants aged 35-64 residing in the two study areas, by sex, cause of death, and birthplace.Spain 2002-2006

	Preferred destination area	Rest of Spain

1. Codes of the International Classification of Diseases, 10th Revision: A00-R99 (all diseases), C00 to C97 (cancer), I00 to I99 (cardiovascular diseases), and rest of codes from A00 to R99 (rest of diseases).

Table 3 .

 3 Age-adjusted mortality rate per 100,000 persons-years, and mortality rate ratios with 95% confidence intervals (95% CI) in immigrants resident in the two areas of Spain and in resident population in the county of birth, by age and cause of death,2002-2006

	Country	Preferred	Country
	of birth	destination area	of birth

Table 4 .

 4 Sample size and population characteristics of immigrants aged 35-64 residing in the two study areas, by sex and birthplace. Spain 2006

Mean time of residence in Spain (in years) Birthplace Men Women Sample size Percentage of retired population Percentage of population with university education

	Preferred destination area	Rest of Spain	Difference	Preferred destination area	Rest of Spain	Difference

All diseases


Pchelintsev Evgeny
email: evgen-pch@yandex.ru

Improved estimation in a non-Gaussian parametric regression

Keywords: Non-Gaussian parametric regression, Improved estimates, Pulse noise, Ornstein-Uhlenbeck process, Quadratic risk, Autoregressive noise Mathematics Subject Classification (2010) 62H12 -62M10

The paper considers the problem of estimating the parameters in a continuous time regression model with a non-Gaussian noise of pulse type. The noise is specified by the Ornstein-Uhlenbeck process driven by the mixture of a Brownian motion and a compound Poisson process. Improved estimates for the unknown regression parameters, based on a special modification of the James-Stein procedure with smaller quadratic risk than the usual least squares estimates, are proposed. The developed estimation scheme is applied for the improved parameter estimation in the discrete time regression with the autoregressive noise depending on unknown nuisance parameters.

Introduction

In 1961, James and Stein proposed shrinkage estimates which outperform in mean square accuracy the maximum likelihood estimates in the problem of estimating the mean of a multidimensional Gaussian vector with unit covariance matrix (James and [START_REF] Stein | Estimation with quadratic loss[END_REF]. This result stimulated the development of the theory of improved estimation for different regression models with dependent errors. Fourdrinier and Strawderman and Fourdrinier and Wells solved the problem of improved parametric estimation in the regression with dependent non-Gaussian observations under spherically symmetric distributions of the noise [START_REF] Fourdrinier | A paradox concerning shrinkage estimators: should a known scale parameter be replaced by an estimated value in the shrinkage factor?[END_REF][START_REF] Fourdrinier | Comparaison de procédures de sélection d'un modèle de régression: une approche décisionnelle[END_REF]. Fourdrinier and Pergamenshchikov and Konev and Pergamenchtchikov investigated the problem of improved estimation in nonparametric setting [START_REF] Fourdrinier | Improved selection model method for the regression with dependent noise[END_REF][START_REF] Konev | General model selection of a periodic regression with a Gaussian noise[END_REF].

In this paper, we consider the problem of improved parametric estimation for a continuous time regression with dependent non-Gaussian noise of pulse type. The noise is specified by the Ornstein-Uhlenbeck process which is known to capture important distributional deviation from Gaussianity and to be appropriate for modelling various dependence structures (Barndorff-Nielsen and Shephard 2001).

Consider a regression model satisfying the equation

dy t = p j=1 θ j φ j (t)dt + dξ t , 0 ≤ t ≤ n, (1)
where θ = (θ 1 , ..., θ p ) ′ (the notation ′ holds for transposition) is the vector of unknown parameters from a compact set Θ ⊂ R p ; (φ j (t)) 1≤j≤p are oneperiodic [0, +∞) → R functions, orthonormal in the space L 2 [0, 1]. The noise (ξ t ) t≥0 in ( 1) is assumed to be a non-Gaussian Ornstein-Uhlenbeck process obeying the following stochastic differential equation

dξ t = aξ t dt + du t , (2)
where a ≤ 0, (u t ) t≥0 is a Levy process which is the mixture

u t = ̺ 1 w t + ̺ 2 z t (3)
of a standard Brownian motion (w t ) t≥0 and a compound Poisson process (z t ) t≥0 defined as

z t = Nt j=1 Y j , (4)
where (N t ) t≥0 is a Poisson process with the intensity λ > 0, and (Y j ) j≥1 is a sequence of i.i.d. Gaussian random variables with parameters (0,1). The noise parameters a, ̺ 1 , ̺ 2 and λ are unknown.

The problem is to construct improved estimates for the unknown vector parameter θ on the basis of observations (y t ) 0≤t≤n , which have higher precision as compared with the least squares estimates (LSE).

It will be observed that the regression model ( 1) is conditionally Gaussian given the σ-algebra G = σ{N t , t ≥ 0} generated by the Poisson process. In Section 3 it is shown that the problem of estimating parameter θ in (1) can be reduced to that of estimating the mean in a conditionally Gaussian distribution with a random covariance matrix depending on the unknown nuisance parameters. This enables one to construct a shrinkage estimates for unknown parameters (θ 1 , . . . , θ p ) in [START_REF] Anderson | The Statistical Analysis of Time Series[END_REF]. The main result is given in Theorem 3.1 which claims that the proposed estimate has less risk than the LSE.

In Section 2 we propose a special modification of James-Stein procedure for solving the problem of estimating the mean in a conditionally Gaussian distribution. This procedure allows one to control the mean square accuracy of estimates. It is shown (Theorem 2.1) that this estimate has less mean square risk than the usual LSE.

The rest of the paper is organized as follows. In Section 4 we apply Theorem 2.1 to the problem of parameter estimation in a discrete time regression under a Gaussian autoregressive noise with unknown parameters. Appendix contains some technical results.

On improved estimation in a conditionally Gaussian regression

In Section 3, we will shown that the initial problem of estimating the parameters (θ 1 , . . . , θ p ) in model [START_REF] Anderson | The Statistical Analysis of Time Series[END_REF] reduces to the following one. Suppose that the observation Y is a p-dimensional random vector which obeys the equation

Y = θ + ξ, (5)
where θ is an unknown constant vector parameter from some compact set Θ ⊂ R p , ξ is a conditionally Gaussian random vector with zero mean and the covariance matrix D(G), i.e. Law(ξ|G) = N p (0, D(G)), where G is some fixed σ-algebra.

The problem is to estimate θ. Consider a shrinkage estimate for θ of the form

θ * = 1 - c Y Y, (6)
where c is a positive constant which will be specified later.

The choice of this estimate ( 6) is motivated by the need to control the quadratic risk

R(θ, θ) = E θ θ -θ 2
in the case of conditionally Gaussian model ( 5). It will be observed that such control can not be provided by the ordinary James-Stein estimate (obtained from ( 6) by the change of Y to Y 2 ).

In order to get an explicit upper bound for the quadratic risk of estimate ( 6) we impose some conditions on the random covariance matrix D(G).

Assume that (C 1 ) There exists a positive constant λ * , such that the minimal eigenvalue of matrix D(G) satisfies the inequality

λ min (D(G)) ≥ λ * a.s. (C 2 )
The maximal eigenvalue of the matrix D(G) is bounded on some compact set Θ ⊂ R p from above in the sense that

sup θ∈Θ E θ λ max (D(G)) ≤ a * ,
where a * is a known positive constant.

Further we will introduce some notation. Let denote the difference of the risks of estimate [START_REF] Fourdrinier | Comparaison de procédures de sélection d'un modèle de régression: une approche décisionnelle[END_REF] and LSE θ = Y as

∆(θ) := R(θ * , θ) -R( θ, θ).
We will need also the following constant

γ p = p-2 j=0 2 j-1 2 (-1) p-j µ p-1-j Γ j+1 2 -(-µ) p I(µ) 2 p/2-1 Γ p 2 d
,

where µ = d/ √ a * , I(a) = ∞ 0 exp(-r 2 /2) a +
sup θ∈Θ ∆(θ) ≤ -[(p -1)λ * γ p ] 2 .
Proof. First we will find the lower bound for the random variable Y -1 .

Lemma 2.2. Under the conditions of Theorem 2.1

inf θ∈Θ E θ 1 Y ≥ γ p .
The proof of lemma is given in the Appendix.

In order to obtain the upper bound for ∆(θ) we will adjust the argument in the proof of Stein's lemma (James and Stein 1961) to the model ( 5) with a random covariance matrix.

We represent the risks of LSE and of (6) as

R( θ, θ) = E θ θ -θ 2 = E θ (E θ -θ 2 |G) = E θ trD(G); R(θ * , θ) = R( θ, θ) + E θ [E((g(Y ) -1) 2 Y 2 |G)] +2 p j=1 E θ [E((g(Y ) -1)Y j (Y j -θ j )|G)],
where

g(Y ) = 1 -c/ Y . Denoting f (Y ) = (g(Y ) -1)Y j and applying the conditional density of distribution of a vector Y with respect to σ-algebra G p Y (x|G) = 1 (2π) p/2 det D(G) exp - (x -θ) ′ D -1 (G)(x -θ) 2 ,
one gets

I j := E(f (Y )(Y j -θ j )|G) = R p f (x)(x -θ j )p Y (x|G)dx, j = 1, p.
Making the change of variable u = D -1/2 (G)(xθ) and assuming f (u) = f (D 1/2 (G)u + θ), one finds that

I j = 1 (2π) p/2 p l=1 D 1/2 (G) jl R p f (u)u l exp - u 2 2 du, j = 1, p,
where A ij denotes the (i, j)-th element of matrix A. These quantities can be written as

I j = p l=1 p k=1 E(< D 1/2 (G) > jl < D 1/2 (G) > kl ∂f ∂u k (u)| u=Y |G), j = 1, p.
Thus, the risk for an estimator (6) takes the form

R(θ * , θ) = R( θ, θ) + E θ ((g(Y ) -1) 2 Y 2 ) +2E θ p j=1 p l=1 p k=1 < D 1/2 (G) > jl < D 1/2 (G) > kl ∂ ∂u k [(g(u) -1)u j ]| u=Y .
Therefore, one has

R(θ * , θ) = R( θ, θ) + E θ W (Y ),
where

W (z) = c 2 + 2c z ′ D(G)z z 3 -2trD(G)c 1 z .
This implies that ∆(θ) = E θ W (Y ).

Since z ′ Az ≤ λ max (A) z 2 , one comes to the inequality

∆(θ) ≤ c 2 -2cE θ trD(G) -λ max (D(G)) Y .
From here, it follows that

∆(θ) ≤ c 2 -2c p i=2 E θ λ i (D(G)) Y .
Taking into account the condition (C 1 ) and the Lemma 2.2, one obtains

∆(θ) ≤ c 2 -2(p -1)λ * γ p c =: φ(c).
Minimizing the function φ(c) with respect to c, we come to the desired result, i.e.

∆(θ) ≤ -[(p -1)λ * γ p ] 2 .
Hence Theorem 2.1.

By applying the Stirling's formula for the Gamma function

Γ(x) = √ 2πx x-1/2 exp(-x) (1 + o(1))
one can check that r p → 0.5 as p → ∞. The behavior of the risk [START_REF] Stein | Estimation with quadratic loss[END_REF] for small values of p is shown in Fig. 1. It will be observed that in this case the risk of the James-Stein estimate θJS remains constant for all p ≥ 3, i.e.

R(0, θJS ) = 2

and the risk of the LSE θ is equal to p and tends to infinity as p → ∞.

Improved estimation in a non-Gaussian Ornstein-Uhlenbeck-Levy regression model

In this section we use the estimate (6) to a non-Gaussian continuous time regression model to construct an improved estimate of the unknown vector parameter θ. To this end we reduce first the initial continuous time regression model [START_REF] Anderson | The Statistical Analysis of Time Series[END_REF] to a discrete time model of the form (5) with a conditionally Gaussian noise. A commonly used estimator of an unknown vector θ in model ( 1) on the basis of observations (y t ) 0≤t≤n is the LSE θ = ( θ1 , . . . , θp ) ′ with the components θj = 1 n n 0 φ j (t)dy t , j = 1, p.

From here and (1), one has

θ = θ + n -1/2 ζ(n), (8)
where ζ(n) is the random vector with the coordinates

ζ j (n) = n -1/2 n 0 φ j (t)dξ t .
Note that the vector ζ(n) has a conditionally Gaussian distribution with a zero mean and conditional covariance matrix

V n (G) = cov(ζ(n), ζ(n) ′ |G) with the elements v ij (n) = E(ζ i (n)ζ j (n) | G).
Thus the initial problem of estimating parameter θ in (1) can be reduced to the that of estimating parameter θ in a conditionally Gaussian regression model [START_REF] Konev | General model selection of a periodic regression with a Gaussian noise[END_REF].

Theorem 3.1. Let the regression model be given by the equations (1)-( 4), ̺ 1 > 0. Then, for any n ≥ 1 and p ≥ 2, the estimator of θ

θ * = 1 - ̺ 2 1 (p -1)γ p n θ θ,
dominates the LSE θ:

sup θ∈Θ ∆(θ) ≤ - ̺ 2 1 (p -1)γ p n 2 .
To prove this theorem one can apply Theorem 2.1. To this end it suffices to check conditions (C 1 ), (C 2 ) on the matrix V n (G). The verification of conditions (C 1 ) and (C 2 ) is given in the Appendix.

Improved estimation in an autoregression

In this section we consider the problem of improved estimating the unknown mean of a multivariate normal distribution when the dispersion matrix is unknown and depends on some nuisance parameters. The models of autoregressive type are widely used in time series analysis [START_REF] Anderson | The Statistical Analysis of Time Series[END_REF][START_REF] Brockwell | Time Series: Theory and Methods[END_REF].

Let the noise ξ = (ξ 1 , . . . , ξ p ) ′ in (5) be described by a Gaussian autoregression process

ξ k = aξ k-1 + ε k , k = 1, p, (9)
where |a| < 1, Eξ 0 = 0 and ε 1 , . . . , ε p are independent Gaussian (0,1) random variables. Assume that the parameter a in ( 9) is unknown and belongs to interval [-α, α], where 0 < α < 1 is known number.

It is easy to check that the covariance of the noise ξ has the form

D(a) = 1 1 -a 2        1 a . . . a p-1 a 1 . . . a p-2
. . .

a p-1 a p-2 . . . 1        Proposition 4.1. Let the noise ξ in (5) be specified by equation (9) with a ∈ [-α, α]. Then, for any p > 1/(1 -α) 2 , the LSE is dominated by the estimate θ * = 1 -p - 1 (1 -α) 2 γ p Y Y in the sense that sup θ∈Θ ∆(θ) ≤ -p - 1 (1 -α) 2 2 γ 2 p .
Proof. We note that trD(a) = p/(1a 2 ). Now we will estimate of the maximal eigenvalue of matrix D(a). From the definition

λ max (D(a)) = sup z =1 z ′ D(a)z one has z ′ D(a)z = p i=1 p j=1 < D(a) > ij z i z j = 1 1 -a 2 1 + 2 p-1 i=1 p-i j=1 a j z i z j+i = 1 1 -a 2 1 + 2 p-1 j=1 a j p-j i=1 z j z i+j .
Applying the Cauchy-Bunyakovskii inequality yields

λ max (D(a)) ≤ 1 1 -α 2 1 + 2 ∞ j=1 α j = 1 (1 -α) 2 .
Thus,

trD(a) -λ max (D(a)) ≥ p - 1 (1 -α) 2 .
By applying Theorem 2.1 we come to the assertion of Proposition 4.1.

Appendix

5.1 Proof of the Lemma 2.2.

Proof. From ( 5), one has

J = E θ 1 Y = E θ 1 θ + ξ ≥ E θ 1 d + ξ .
Taking the repeated conditional expectation and noting that the random vector ξ is conditionally Gaussian with zero mean, one gets

J ≥ E θ 1 (2π) p/2 detD(G) R p exp(-x ′ D(G) -1 x/2) d + x dx.
Making the change of variable u = D(G) -1/2 x and applying the estimation u ′ D(G)u ≤ λ max (D(G)) u 2 , we find

J ≥ 1 (2π) p/2 R p exp(-u 2 /2) d + λ max (D(G)) u du.
Further making the spherical changes of the variables yields

J ≥ 1 2 p/2-1 Γ(p/2) E θ ∞ 0 r p-1 exp(-r 2 /2) d + λ max (D(G))r dr.
From here by applying the Jensen and Cauchy-Bunyakovskii inequalities and by the condition (C 2 ), we obtain

J ≥ µ 2 p/2-1 Γ(p/2)d ∞ 0 r p-1 exp(-r 2 /2) µ + r dr = γ p .
This leads to the assertion of Lemma 2.2.

5.2

The verification of the conditions (C 1 ) and (C 2 ) on the matrix V n (G).

Now we establish some properties of a stochastic integral

I t (f ) = t 0 f s dξ s , 0 ≤ t ≤ n
with respect to the process (2). We will need some notations. Let us denote

ε f (t) = a t 0 exp{a(t -v)} f (v) (1 + exp(2av)) dv ,
where f is [0, +∞) → function integrated on any finite interval. We introduce also the following transformation

τ f,g (t) = 1 2 t 0 2f (s)g(s) + ε * f,g (s) ds of square integrable [0, +∞) → R functions f and g. Here ε * f,g (t) = f (t)ε g (t) + ε f (t)g(t) .
E I t (f )I t (g) = ̺ * τ f,g (t) ( 10
)
where

̺ * = ̺ 2 1 + λ̺ 2 2 .
Proof. Noting that the process I t (f ) satisfies the stochastic equation

dI t (f ) = af (t)ξ t dt + f (t)du t , I 0 (f ) = 0 ,
and applying the Ito formula, one obtains (10). Hence Proposition 5.1.

Corollary 5.2. If f is function from L 2 [0, n] then E I 2 n (f ) ≤ 3̺ * n 0 f 2 (t)dt. (11)
Further, for an integrated [0, +∞) → R function f , we define the function

L f (x, z) = a exp(ax) f (z) + a x 0 exp(av) f (v + z) dv .
Let (T l ) l≥1 be the jump times of the Poisson process (N t ) t≥0 , i.e.

T l = inf{t ≥ 0 : N t = l}.

Proposition 5.3. Let f and g be bounded left-continuous [0, ∞) × Ω → R functions measurable with respect to B[0, +∞) G (the product σ algebra created by B[0, +∞) and G). Then This implies that where

E (I t (f )|G) = 0 and E (I t (f ) I t (g)|G) = ̺ 2 1 τ f,g (t) + ̺ 2 2 l≥1 f (T l )g(T l )1 (T l ≤t) + ̺ 2 2 l≥1 t 0 (f (s)L g (s -T l , T l ) + g(s)L f (s -T l , T l )) 1 (T l ≤s) ds.
z ′ V n (G)z = ̺ 2 1 z ′ z + z ′ F n z + z ′ B n (G)z ≥ ̺ 2 1 z ′ z, therefore inf z =1 z ′ V n (G)z ≥ ̺ 2
̺ * = ̺ 2 1 + λ̺ 2 2 .
Proof. We note that

E θ λ max (V n (G)) ≤ E θ tr(V n (G)) = p j=1 E θ ζ 2 j (n) = 1 n p j=1 E θ I 2 n (φ j ).
Applying (11) and n 0 φ 2 j (t)dt = n, we obtain the desired inequality. Hence Lemma 5.5.

Thus the matrix V n (G) is positive definite and satisfies for any compact set Θ ⊂ R p , the conditions (C 1 ) and (C 2 ) with λ * = ̺ 2 1 and a * = 3p̺ * .

  r dr and d = sup{ θ : θ ∈ Θ}. Theorem 2.1. Let the noise ξ in (5) have a conditionally Gaussian distribution N p (0, D(G)) and its covariance matrix D(G) satisfy conditions (C 1 ), (C 2 ) with some compact set Θ ⊂ R p . Then the estimator (6) with c = (p -1)λ * γ p dominates the LSE θML for any p ≥ 2, i.e.

Corollary 2 . 3 .Corollary 2 . 4 . 2 =

 23242 Let in (5) the noise ξ ∼ N p (0, D) with the positive definite non random covariance matrix D > 0 and λ min (D) ≥ λ * > 0. Then the estimator (6) with c = (p -1)λ * γ p dominates the LSE for any p ≥ 2 and compact set Θ ⊂ R p , i.e. sup θ∈Θ ∆(θ) ≤ -[(p -1)λ * γ p ] 2 . Remark 2.1. Note that if D = σ 2 I p then sup θ∈Θ ∆(θ) ≤ -[(p -1)σ 2 γ p ] 2 . If ξ ∼ N p (0, I p ) and θ = 0 in model (5) then the risk of estimate (6) is given by the formula R(0, θ * ) = p -(p -1)Γ((p -1)/2) √ 2Γ(p/2) : r p .

Figure 1 :

 1 Figure 1: Risk of θ * at θ = 0.

Proposition 5 . 1 .

 51 If f and g are functions from L 2 [0, n] then

Proof. 0 (̺ 2 1 f 0 ̺ 2 1 fa t 0 (

 01010 By the Ito formula one hasI t (f ) I t (g) = t (s)g(s) + a(f (s)I s (g) + g(s)I s (f ))ξ s )ds + ̺ 2 2 l≥1 f (T l ) g(T l ) Y 2 l 1 {T l ≤t} + t 0 (f (s)I s-(g) + g(s)I s-(f )))du s .Taking the conditional expectation E (•|G), on the set {T l > t}, yieldsE (I t (f ) I t (g)|G) = t (s)g(s)ds + ̺ 2 2 l≥1 f (T l ) g(T l ) 1 {T l ≤t} + f (s)E(I s (g)ξ s |G) + g(s)E(I s (f )ξ s |G)) ds .

1 and we come to the assertion of Lemma 5. 4 . 5 . 5 .

 1455 Lemma Let (ξ t ) t≥0 be defined by[START_REF] Barndorff-Nielsen | Non-Gaussian Ornstein-Uhlenbeck-based models and some of their uses in financial mathematics[END_REF] with a ≤ 0. Then a maximal eigenvalue of the matrix V n (G) = (v ij (n)) 1≤i,j≤p with elements defined by (12), satisfy the following inequalitysup n≥1 sup θ∈Θ E θ λ max (V n (G)) ≤ 3p̺ * .

It is easy to check that

L f (t -T j , T j ) 1 {T j ≤t} .

From here one comes to the desired equality. Hence Proposition 5.3.

Thus, in view of ζ j (n) = n -1/2 I n (φ j ) and Proposition 5.3 the elements of covariance matrix V n (G) can be written as

Lemma 5.4. Let (ξ t ) t≥0 be defined by (2) with a ≤ 0. Then a matrix V n (G) = (v ij (φ)) 1≤i,j≤p with elements defined by (12), satisfy the following inequality a.s.

Proof. Notice that by (12) one can the matrix V n (G) present as

where F n is non random matrix with elements

+ n 0 φ i (t)L φ j (t -T l , T l ) + φ j (t)L φ i (t -T l , T l ) 1 (T l ≤t) dt].


Mehdi Lallouache
email: mehdi.lallouache@ens-cachan.fr

Anindya S Chakrabarti
email: aschakrabarti@gmail.com

Anirban Chakraborti
email: anirban.chakraborti@ecp.fr

Bikas K Chakrabarti
email: bikask.chakrabarti@saha.ac.in

Opinion formation in the kinetic exchange models: Spontaneous symmetry breaking transition

Keywords: numbers: 87.23.Ge 02.50.-r Econophysics, Sociophysics, kinetic theory

We propose a minimal multi-agent model for the collective dynamics of opinion formation in the society, by modifying kinetic exchange dynamics studied in the context of income, money or wealth distributions in a society. This model has an intriguing spontaneous symmetry breaking transition to polarized opinion state starting from non-polarized opinion state. In order to analyze the model, we introduce an iterative map version of the model, which has very similar statistical characteristics. An approximate theoretical analysis of the numerical results are also given, based on the iterative map version.

I. INTRODUCTION

Recently physicists have been studying social phenomena and dynamics leading to the growth of the interdisciplinary field of "Sociophysics" [START_REF] Eds | Chatterjee, Econophysics and Sociophysics: Trends and Perspectives[END_REF]. One of the problems is of "opinion formation", which is a collective dynamical phenomenon, and as such is closely related to the problems of competing cultures or languages [START_REF] Castellano | [END_REF]3]. It deals with a "measurable" response of the society to e.g., political issues, acceptances of innovations, etc. Numerous models of competing options have been introduced to study this phenomenon, e.g., the "voter" model (which has a binary opinion variable with the opinion alignment proceeding by a random choice of neighbors) [4], or the Sznajd-Weron discrete opinion formation model (where more than just a pair of spins is associated with the decision making procedure) [5]. There have been other studies of systems with more than just two possible opinions [6], or where the opinion of individuals is represented by a "continuous" variable [7][8][9] using real numbers. Also, since opinion formation in a human society is mediated by social interactions between individuals, such social dynamics has been considered to take place on a network of relationships (see [START_REF] Castellano | [END_REF] for recent review on such models).

A two body exchange dynamics has already been developed in the context of modelling income, money or wealth distributions in a society [10][11][12][13][14]. The general aim was to study a many-agent statistical model of closed economy (analogous to the kinetic theory model of ideal gases) [15], where N agents exchange a quantity x, that may be defined as wealth. The states of agents are characterized by the wealth {x i }, i = 1, 2, . . . , N , such that x i > 0, ∀i and the total wealth W = i x i is conserved. The question of interest is: "What is the equilibrium distribution of wealth f (x), such that f (x)dx is the probability that in the steady state of the system, a randomly chosen agent will be found to have wealth between x and x + dx?" The evolution of the system is carried out according to a prescription, which defines the trading rule between agents. The agents interact with each other through a pair-wise interaction characterized by a "saving" parameter λ, with 0 ≤ λ ≤ 1. The dynamics of the model (CC) is as follows [15]:

x ′ i = λx i + ǫ(1 -λ)(x i + x j ) , x ′ j = λx j + (1 -ǫ)(1 -λ)(x i + x j ) , (1)
where ǫ (0 ≤ ǫ ≤ 1) is a stochastic variable, changing with time. It can be noticed that in this way, the quantity x is conserved during the single transactions: x ′ i + x ′ j = x i + x j , where x ′ i and x ′ j are the agent wealths after the transaction has taken place. In general, the functional form for steady state distribution f (x) is seen to be close to the Γ-distribution [16,17]. As a further generalization, the agents could be assigned different saving propensities and the steady state distribution f (x) show Pareto-like power-law behavior asymptotically [18,19].

Earlier, Toscani [20] had introduced and discussed kinetic models of (continuous) opinion formation involving both exchange of opinion between individual agents and diffusion of information. Based on this model, During et al [21] proposed another mathematical model for opinion formation in a society that is built of two groups, one group of ordinary people and one group of strong opinion leaders. Starting from microscopic interactions among individuals, they arrived at a macroscopic description of the opinion formation process. They discussed the steady states of the system, and extended it to incorporate emergence and decline of opinion leaders.

Here, we report the studies of a minimal model for the collective dynamics of opinion formation in the society, based on kinetic exchanges.

II. MODEL FOR OPINION FORMATION AND RESULTS

A. Homogeneous multi-agent model

Following the CC model described in the earlier section, we present a minimal model [22] for the collective dynamics of opinion O i (t) of the i-th person in the society, consisting of N (N -→ ∞) persons. We assume that any particular person can discuss (interact) only with one other person each time (time increases discretely by unity after each such discussion). A two-person "discussion" is viewed here as a simple two-body scattering process in physics. Persons in the society may bump onto each other randomly and exchange opinions through such random two-person discussions. In general, a person i could have any opinion O i between two extreme polarities denoted by +1 and -1. In any discussion at time t + 1, a person retains a fraction of his/her older opinion O i (t), determined by his/her "conviction", parameterized by λ i . This parameter value is characteristic of a person and does not change with time t. Additionally, the person i is "influenced" stochastically by the other person j during the discussion having the "influence" parameter equal to his/her conviction parameter λ j . We further assume for simplicity that all agents are homogeneoushave the same conviction parameter λ. Mathematically the dynamics may be represented by

O i (t + 1) = λ(O i (t) + ǫ t O j (t)) , O j (t + 1) = λ(O j (t) + ǫ ′ t O i (t)) , (2)
where the opinion -1 ≤ O i (t) ≤ 1 for all agents i and time t, the conviction parameter 0 ≤ λ ≤ 1 is quenched 2). (Inset) Numerical results for the variation of the variance (∆O) 2 ≡ (O -Ō) 2 against λ, following dynamics of Eq. ( 2).

(does not change with time), and the stochastic parameters ǫ t and ǫ ′ t are annealed variables (change with time)uncorrelated random numbers uniformly distributed between zero and unity. Note that the equations are linear, but non-linearity is introduced in this model by imposing that -1 ≤ O i (t) ≤ 1 for all agents i and times t.

The question we are interested is that if such social dynamics continually take place, can any consensus be reached or polarity evolve after a long time? Mathematically, we are interested in the steady state distribution of O and other statistical properties. It is noteworthy that unlike in the market models, here we have no conservation of opinion. Rather, the steady state of value of Ō(t) = (1/N )|Σ i O i (t)| represents the order of the average opinion in the society after a long time t. We study the relaxation dynamics in the society: the relaxation and fluctuation of Ō, the steady state value of Ō(t) for t > τ , the relaxation time.

Remarkably, we find there is appearance of "polarity" or consensus, starting from initial random disorder (where O i 's are uniformly distributed with positive and negative values). In the language of physics, there is a "spontaneous symmetry breaking" transition in the system: starting from Ō(0) = 0 the system evolves either to the "para" state with Ō ≡ Ō(t > τ ) = 0 (where all individual agents have the opinion 0) for λ ≤ 2/3, or (continuously) to the "symmetry broken" state Ō ≡ Ō(t > τ ) = 0 (where all individuals have either positive or negative opinions) for λ ≥ 2/3 (see Fig. 1) for times t > τ . We note, however, that the fluctuation in Ō does not diverge, and shows a cusp near λ c (see inset of Fig. 1). We also study the relaxation behavior of Ō(t) and the critical divergence of the relaxation time τ near λ = λ c = 2/3 (see subsection II C, Fig. 4).

B. Random multiplier map

The basic nature of transition produced by Eq. ( 2), can perhaps be reproduced by the following simple iterative map

O(t + 1) = λ(1 + ǫ t )O(t) (3)
with the restriction that O(t) ≤ 1, which is ensured by assuming that if O(t) ≥ 1, O(t) is set equal to 1. As usual, ǫ t is a stochastic variable ranging between 0 and 1 (assumed to be uniformly distributed in our case). In a mean-field approach, the above equation reduces effectively to a multiplier map like O(t + 1) = λ(1 + ǫ )O(t), where ǫ = 1/2. Clearly for λ ≤ 2/3, O(t) converges to zero. The initial value O(0), is assigned either a positive or negative value. If it starts from a positive (negative) value, O(t) remains positive (negative). We note that there are subtle differences in the dynamics of Eq. ( 2) and Eq. (3). Apart from the absence of "spontaneous symmetry breaking" of the multi-agent model (from ±O i (0) values to all positive or all negative transition beyond λ c ), the nature of the phase transition (singularity) in the iterative map is also slightly different. The critical value λ c = exp {-(2 ln 2 -1)} ≈ 0.6796 has an analytical derivation [23], but for most numerical studies done here, we take λ c = 0.68. The time variation of the average opinion Ō(t) = (1/N )Σ i |O i (t)|, where i refers to different initial realizations and N refers to the total of all such realizations, and its fluctuations are studied numerically (see Fig. 2). We study the relaxation behavior of Ō(t) and the critical divergence of the relaxation time τ near λ = λ c = 0.68 (see subsection II C, Fig. 4). We note again that the fluctuation in Ō does not diverge, and shows a cusp near λ c = 0.68 (see inset of Fig. 2). We also note that the steady state fluctuation ∆O near λ c , is generally much higher in magnitude for the map case.

C. Results and analyses

For both the multi-agent model and the iterative map, we study the variation (with λ) of the fraction p of the agents having O i = ±1 at any time t in the steady state (t > τ ). This parameter p gives the average "condensation" fraction (of people in the society having extreme opinions |O i | = 1) in the steady state. The growth of p, as shown in Fig. 3, is seen to be similar to that of Ō. The inset shows that the growth behavior for p above (respective) λ c , for both the multi-agent model and map, are identical.

We studied the relaxation behavior of Ō and p, for both the multi-agent model and map. In each case, the relaxation time is estimated numerically from the time value at which Ō or p first touches the steady state value within a pre-assigned error limit. We find diverging growth of relaxation time τ (for both Ō and p) near λ = λ c (see Fig. 4). The values of exponent z for the divergence in τ ∼ |λ -λ c | -z have been estimated numerically for both the multi-agent model and the map (for both λ > λ c and λ < λ c , wherever accurate data were obtained). For the multi agent model, the fitting values for exponent z corresponding to Ō and p, respectively, are z ≈ 1.0 ± 0.1 and z ≈ 0.7 ± 0.1. For the map case, the fitting values for exponent z corresponding to both Ō and p, turn out to be the same: z ≈ 1.5 ± 0.1.

For the iterative map Eq. ( 3), we study carefully the time evolution of the condensation fraction p of |O| = 1 in different realizations at different values of λ. The variation of the steady state value p against λ is shown in Fig. 2. It may be noted that while the steady state value of Ō starts to grow from λ ≈ 2/3 (see Fig. 2) , the steady state value of p starts growing at λ ≈ 0.68 (see Fig. 3). Numerical results for the growth of the relaxation time τ for both Ō and p, against λ are shown in Fig. 4. Both diverge at λ ≈ 0.68. This clearly indicates that p, rather than Ō, is the order parameter for the transition.

An approximate analysis of the above transition for λ closer to unity can be done for the iterative map Eq. (3) as follows. In Fig. 5, we give the numerical results for the steady state distribution opinion, P (|O|) for three different values of λ; we observe roughly a bi-modal nature of the distribution as λ → 1: one mode is the uniform distribution within the range |O min | < |O| < 1 (and |O min | ≈ λ) and another a delta function at |O| = 1. We therefore approximate the steady state distribution of opinion by assuming that opinion O(t) is distributed uniformly starting from a minimum O min upto unity with (integrated) probability (1 -p), and a δ-function at exactly unity with probability p.

Then Ō = (1 -p)O av + p.1, (4)
where O av = (O min + 1)/2. We have assumed that the value O(t) stays in those two regions (from λ to 1 and exactly at 1) with probability (1 -p) and p. Hence, the corresponding equations are

O(t + 1) = λ(1 + ǫ)O(t)
with probability 1 -p, Note that the first equation is realized

only if λ(1 + ǫ)O(t) < 1 or ǫ < ǫ max = 1 λOav -1. This cut-off implies that p = ǫmax 0 dǫ = 1 λOav -1, since ǫ ∼uni[0,1]
. By substituting O av and p in Eq. ( 4), we derive the result that

Ō = 5λ + 2λ 2 -λ 3 -2 2λ(1 + λ) (5)
which is compared with the numerical simulations for λ → 1 in Fig. 6. It is evident that the approximation holds well, only for λ → 1.

III. SUMMARY AND DISCUSSION

In summary, we proposed a minimal model for the collective dynamics of opinion formation in the society, by modifying kinetic exchange dynamics studied in the context of markets. The multi-agent model (dynamics given by Eq. ( 2)) and its map version (dynamics given by Eq. (3)) have kinetic exchange like linear contributions from random two-person discussions or scattering processes, though the saturation of |O i | ≤ 1 induces nonlinearity in the dynamics. This model has an intriguing spontaneous symmetry breaking transition to polarized opinion state starting from non-polarized opinion state. Specifically, in the multi-agent model, we see that for λ > λ c = 2/3, starting from random positive and negative O i values (or for that matter any arbitrary state), at t = 0, the system eventually evolves to a state at t > τ where all O i are either positive or negative, with | Ō| determined by the λ value! This is similar to the growth of spontaneous magnetization in Ising magnets (where starting from arbitrary up and down spin states, a preferred direction is chosen by fluctuation), with magnetization determined by the temperature below its transition value. The appearance of spontaneous symmetry breaking in this simple kinetic opinion exchange model is truly remarkable. It appears to be one of the sim- plest collective dynamical model of many-body dynamics showing non-trivial phase transition behaviour. Indeed, it may be noted that for λ ≤ λ c , at t > τ , all O i 's become identically zero (without any fluctuation), while for λ > λ c , O i 's have fluctuations but the average has a steady state value depending on the value of λ. As such the nature of the phase transition in this model is quite different and does not fit to the commonly studied three absorbing state models studied (see e.g., [24]). In order to understand the nature of the transition, we also studied a simple iterative map and derived approximate result for the order parameter variation under certain limits, which compares quite well with the numerical simulations. Specifically, we find that the fraction p of people with extreme opinion |O i | = 1, and its fluctuations determine the nature of the phase transition in our model and locate the critical point accurately (from numerical studies). With the two mode distribution (uniform and delta) of O, valid close to λ → 1 (see Fig. 5), we could develop an approximate analysis of the variation of the steady state mean opinion | Ō| against λ as in Eq. ( 5). In any case, further investigations are necessary for understanding this phase transition. Additional studies for the heterogeneous conviction factors λ i 's, in influence of "field terms" that represent the external influence of media, etc. will be reported elsewhere [25]. Also, the study of this phase transition behavior for an extended model with separate"conviction" and "influence" parameters in Eq. ( 2), has recently been reported [26].

2 λFIG. 1 .

 21 FIG.1. Numerical results for the variation of the average opinion Ō(t) for large t (steady state value of Ō) against λ, following dynamics of Eq. (2). (Inset) Numerical results for the variation of the variance (∆O) 2 ≡ (O -Ō) 2 against λ, following dynamics of Eq. (2).

2 λFIG. 2 .

 22 FIG.2. Numerical results for the variation of the average opinion Ō(t) for large t (steady state value of Ō) against λ, following dynamics of Eq. (3). (Inset) Numerical results for the variation of the variance (∆O) 2 ≡ (O -Ō) 2 against λ, following dynamics of Eq. (3).

FIG. 3 .

 3 FIG.3. Numerical results for the variation of the average condensate fraction p(t) for large t (steady state value of p) against λ, following dynamics of Eq. (2) in black diamonds, and dynamics of Eq. (3) in red circles. (Inset) Numerical results for the growth of p, following dynamics of Eq. (2) in black diamonds, and dynamics of Eq. (3) in red circles, close to λc.

FIG. 4 .

 4 FIG. 4. Numerical results for relaxation time behaviors τ versus λλc, for (a) Multi-agent model with Ō (b) Multiagent model with p (c) Map with Ō (d) Map with p . (Insets) Determination of exponent z from numerical fits of τ ∼ |λ -λc| -z .

FIG. 5 .FIG. 6 .

 56 FIG. 5. Numerical results for for the steady state distribution opinion, P (|O|), for three values λ = 0.8, 0.85, 0.9 showing bimodal distributions in each case. (Inset). The same steady state distribution P (|O|), for three values λ = 0.8, 0.85, 0.9, but close to |O| = 1.

ACKNOWLEDGMENTS

The authors acknowledge F. Abergel, A. Chatterjee, A. Jedidi, K. Kaski, S.N. Majumdar, M. Marsili, N. Millot, M.A. Muñoz, M. Patriarca and S. Pradhan for useful discussions, comments and correspondences.