Z Navrátil

P Dvořák

O Brzobohatý

D Trunec

Determination of electron density and temperature in capacitively coupled RF discharge in neon by OES complemented with CR model

Keywords: numbers: 52, 80, Hc, 52, 70, Kz

In this work a method of determination of electron temperature and electron density in plasma based on optical emission spectroscopy complemented with collisional-radiative modelling (OES/CRM) was studied. The radiofrequency (13.56 MHz) capacitively-coupled discharge in neon at 10 Pa was investigated by intensity calibrated optical emission spectroscopy. The absolute intensities of neon transitions between 3s -3p states were fitted with a collisional-radiative model in order to determine the electron temperature and electron density. Measuring techniques as imaging with ICCD camera were adopted for supplementary diagnostics. The obtained results were compared with the results of compensated Langmuir probe measurement and one dimensional Particle In Cell / Monte Carlo (PIC/MC) simulation.

The results of OES/CRM and PIC/MC method were in close agreement in the case of electron temperature in vicinity of a driven electrode. The determined value of electron temperature was about 8 eV. In bulk plasma, the measured spectra were not satisfactorily fitted. In the case of electron density only relative agreement was obtained between OES/CRM and Langmuir probe measurement, the absolute values differed by a factor of 5. The axial dependence of electron density calculated by PIC/MC was distinct from them, reaching the maximum values between the results of the other two methods. The investigation of power dependence of plasma parameters close to driven electrode showed a decrease in electron temperature and an increase in electron density together with increasing incoming RF power. The calculated spectra fitted very well the measured spectra in this discharge region.

Introduction

Radiofrequency (RF) capacitively coupled discharges are widely used in industrial applications for thin film deposition, etching, cleaning, surface modification, in high power lasers, as analytical sources, light sources etc [START_REF] Raizer | Radio-Frequency Capacitive Discharges[END_REF]. However, precise measurement of the basic plasma parameters such as electron density and electron temperature or electron distribution function is still demanding. Although the use of Langmuir probe in RF discharges is well established, it is still limited to low pressure plasmas, risking the plasma disturbance and suffering with low signal to noise ratio for highly energetic electrons, which have a low concentration but high impact on plasma processes.

On the other hand, optical emission spectroscopy (OES) seems to overcome all these drawbacks. However, in order to derive the plasma parameters from the measured optical emission spectra, a relatively complicated collisional-radiative modelling (CR modelling) is generally needed.

In case of rare gases, the capacitively coupled RF discharge has been extensively studied in helium and argon, but less often in neon. Electron distribution function (EDF) has previously been measured with Langmuir probe in argon and helium [START_REF] Godyak | Measurement of electron energy distribution in low-pressure RF discharges Plasma Sources[END_REF].

Anisotropy of electron energy distribution function was measured by electrostatic energy analyser in asymmetric helium RF discharge [START_REF] Okuno | Measurements of electron energy distribution function in an asymmetric radio-frequency discharge plasma[END_REF]. Thomson light scattering was adopted for low electron temperature and density measurement in helium RF discharge [START_REF] Wesseling | Thomson light scattering measurements of electron temperature and density in the α-γ transition of a capacitive rf discharge in helium[END_REF].

Transition between α and γ mode of RF discharge in argon was investigated in [START_REF] Deegan | Measurement of the electron energy distribution function in an argon radio-frequency discharge in the gamma mode[END_REF]. Timeaveraged axial potential profile was measured with the electrostatic probe in symmetric RF discharge in neon [START_REF] Kaneda | Time-averaged electric potential profiles in a capacitive-coupling parallel-plate electrode neon gas RF discharge plasma[END_REF]. Paper [START_REF] Kimura | Pressure dependences of electron energy distribution and power dissipation in symmetrical RF discharges of inert gases[END_REF] compared EDFs measured in a discharge centre at different pressures (0.1 -0.8 torr) in neon, argon and xenon. EDF in capacitive RF discharge in molecular gases was also investigated (see e.g. [START_REF] Turner | Anomalous sheath heating in a low-pressure RF discharge in nitrogen[END_REF][START_REF] Tatarova | Effects of nonlocal electron kinetics and transition from alpha to gamma regime in an RF capacitive discharge in nitrogen[END_REF] for nitrogen).

The cited papers show the necessity of knowledge of electron distribution function for understanding capacitive discharges, including understanding of heating mechanisms and prediction of reactions induced by electron collisions. Further important information concerning the physics of capacitively coupled discharges can be obtained by methods, that enable measurements with high temporal resolution. Uncompensated electric probes enable measurement of plasma potential changes occurring during one RF period [START_REF] Sobolewski | Electrical characteristics of argon radio frequency glow discharges in an asymmetric cell[END_REF][START_REF] Dvořák | Measurement of plasma potential waveforms by an uncompensated probe Plasma Sources[END_REF]. Optical emission spectroscopy with high temporal resolution brings an illustrative look into the processes occurring in each RF period, including discharge heating in α and γ modes [START_REF] Schulze | Different modes of electron heating in dual-frequency capacitively coupled radio frequency discharges Plasma Sources[END_REF], generation of electron beams by sheath expansion [START_REF] Schulze | Electron beams in asymmetric capacitively coupled radio frequency discharges at low pressures[END_REF] and field reversal during sheath collapse [START_REF] Schulze | Electric field reversals in the sheath region of capacitively coupled radio frequency discharges at different pressures[END_REF]. Despite the fact that these measurements brought numerous facts on time-development of capacitive discharges, they did not provide information concerning EDFs, electron concentration or temperature.

Therefore, measurement of EDFs based on OES, which is noninvasive and potentially enables high temporal resolution is still desired.

Computer simulation, especially collisional-radiative modelling, has become widely used for diagnostics of various types of discharges [START_REF] Vlček | A collisional-radiative model applicable to argon discharges over a wide range of conditions. I: Formulation and basic data[END_REF][START_REF] Golubovskii | Nonlocal electron kinetics and radiation of a stratified positive column of discharge in neon[END_REF][START_REF] Zhu | Measurement of the electron density in atmospheric-pressure low-temperature argon discharges by line-ratio method of optical emission spectroscopy[END_REF]. The method of electron temperature measurement based on a collisional-radiative model for argon was reported e.g. in [START_REF] Kano | Spectroscopic measurements of electron temperature and density in argon plasmas based on collisional-radiative model Plasma Sources[END_REF]. Trace rare gases optical emission spectroscopy (TRG-OES), based on the addition of a small admixture of rare gas into the studied plasma and evaluation of plasma parameters from the best fit between the measured and with CR model calculated relative emission intensities, was reviewed in [START_REF] Donnelly | Plasma electron temperatures and electron energy distributions measured by trace rare gases optical emission spectroscopy[END_REF]. Comparative diagnostics of helium and argon microwave plasmas at a pressure of 1 -5 Torr utilizing CR model, OES and probe measurement was published in [START_REF] Mizuochi | Evaluation of electron energy distribution function in microwave discharge plasmas by spectroscopic diagnostics with collisional radiative model[END_REF]. Statistical analysis of reconstruction of electron distribution function in the positive column of neon dc discharge from the emission spectra was presented in [START_REF] Dodt | Reconstruction of an electron energy distribution function using integrated data analysis[END_REF]. Application of Particle In Cell/Monte Carlo method (PIC/MC) for investigation of plasma was published in [START_REF] Birdsall | Plasma Physics Via Computer Simulation[END_REF][START_REF] Hockney | Computer simulation using particles[END_REF][START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, molecule-molecule, and Coulomb collisions for particle modeling of materials processing plasmas and gases[END_REF].

In this work a method of the determination of electron temperature and density from temporally averaged optical emission spectra is developed for neon plasma and compared with the results from Langmuir probe measurement and PIC/MC simulation.

The paper is organized as follows: details of experimental apparatus, optical emission spectroscopy and Langmuir probe diagnostics are given in section 2. Collisional-radiative model and Particle in Cell / Monte-Carlo approach are described in sections 3 and 4, respectively. The results obtained with the above mentioned methods are presented, compared and discussed in section 5. The conclusions are given in section 6.

Experimental set-up

The schematic of the experimental set-up is displayed in figure 1. The RF discharge was generated in a grounded stainless-steel vacuum chamber with an inner diameter of 33 cm. Two parallel plate circular electrodes with a diameter of 80 mm were placed 40 mm apart. The upper electrode (located at position 0 mm), embedded in a grounded ring, was capacitively coupled to a RF sinusoidal voltage generator (13.56 MHz, max. power 50 W, max. applied voltage 470 V pp ). The lower electrode (located at position 40 mm) was grounded. The DC self-bias on the powered electrode varied within -130 V to -300 V showing an asymmetry of the discharge.

The vacuum chamber was pumped down using a turbomolecular pump to an ultimate pressure 5 × 10 -5 Pa and then filled with neon gas of research purity 5.0 up to pressure of 5 -15 Pa. The discharge was sustained in flowing regime with gas flow rate approx. 5 sccm. The butterfly valve between the chamber and the pump controlled the pumping speed, allowing the pressure to be set independently of the flow rate.

The pressure in the chamber was measured by MKS Baratron. A large window, with diameter 7 cm, enabled the observation of the discharge.

Langmuir probe measurement

Electron energy distribution function, electron concentration, mean electron energy and mean (DC) plasma potential, were measured using a compensated Langmuir probe (ESPion, Hiden Analytical). The probe consisted of a 1 cm long platinum wire, radius harmonic frequencies was passively compensated [START_REF] Hopkins | Langmuir probe measurements in the Gaseous Electronic Conference RF reference cell[END_REF]. Eventual disturbance of plasma potential due to high sheath resistivity was compensated by means of a reference probe [START_REF] Kleber | Sheath resistance measurements in the GEC reference reactor Plasma Sources[END_REF]. Electron energy distribution functions were calculated from second derivative of the probe characteristics by the well-known Druywestein formula. If EDFs had Maxwellian shape, the electron concentration and mean energy were calculated from a linear fit to a logarithm of measured probe characteristics in order to exclude the impact of noise produced by numerical calculation of the second derivative. Analogically, when it was possible to describe EDF by so-called general distribution, appropriate function was fitted to the logarithm of the measured data. The general distribution f g is described by relation

f g = C g n e E 1/2 exp - E κ 2κE κ p
where n e , E and E p are electron concentration, energy and the most probable energy, respectively. C g is a constant depending only on E p and κ. κ describes the shape of the EDF. For κ = 1 or κ = 2 the general distribution is identical to Maxwellian or Druywestein distribution, respectively. When EDFs had another shape, electron concentration and mean energy were calculated by integration of EDFs obtained by the Druywestein formula. At the pressure 10 Pa the measured EDFs were similar to Maxwellian EDF. The statistical error of electron temperature and density determined from repeated probe measurement was below 5 % and 4 %, respectively.

Optical emission spectroscopy

The emission of radiation from the discharge was studied by optical spectroscopy in UV/VIS spectral range (300 -750 nm). The spatial resolution of the measurement was achieved by two iris diaphragms, mounted with optical fibre on a movable table in front of the chamber window. The light coming from a narrow cone going through the whole discharge diameter, but with a distinct axial position between the electrodes, was taken out with the diaphragms and optical fibre and analysed with Jobin Yvon FHR 1000 spectrometer with CCD detector (focal length 1 m, grating 2400 gr mm -1 ). The spectra were integrated over many periods of RF signal. An Oriel Tungsten halogen lamp was used to calibrate the fibre irradiance. Original uncertainty of the lamp intensity was 2 %. A new recalibration with 5 % uncertainty revealed as much as 10 % difference in the calibrations.

The emission coefficient ǫ ij of specific transition i → j, integrated over the lineshape, was derived from the fibre irradiance. If the cylindrical symmetry of plasma with uniform radial profile for region between the electrodes can be assumed

ǫ ij (r) = ǫ ij , r ≤ R 0, r > R,
where ǫ ij is the emission coefficient of unit discharge volume defined as

ǫ ij = 1 4π n i A ij Λ ij hω ij ,
where n i is the population of excited state i, A ij and Λ ij Einstein coefficient and socalled escape factor of spontaneous emission, respectively, and hω ij is the corresponding photon energy, then the irradiance of the optical fibre ("detector") is

I ij = 1 S det V pl S det ǫ ij (r) ρ 2 Acc(θ)dV pl dS det ,
where ρ is the distance of the plasma element dV pl from the detector part dS det , θ is the angle of incidence and Acc(θ) represents the fibre acceptance function, determined experimentally. The escape factor Λ ij aproximates the reduction of emitted photons due to reabsorption [START_REF] Holstein | Imprisonment of resonance radiation in gases[END_REF][START_REF] Molisch | Radiation-trapping in cylindrical and spherical geometries[END_REF]. Under the above mentioned considerations the emission coefficient can be taken out from the integral. In our case, plasma region contributing through the iris diaphragms to the irradiation of fibre aperture was assumed to have cylindrical shape with length 110 mm and diameter 5 mm. Constant value of emission coefficient along this length and diameter was assumed and taken into calculations.

Supplementary temporally and spatially resolved, but spectrally unresolved measurement of discharge emission of radiation was carried out with an intensified CCD camera PI-MAX 1024RB-25-FG-43, controlled by a ST-133 controller and working in image mode. The temporal resolution was 5 ns and the spatial resolution 0.3 mm.

Collisional-radiative model

The collisional-radiative model (CR model) used in this work was described in detail in [START_REF] Navrátil | Collisional-radiative model of neon discharge: determination of E/N in the positive column of low pressure discharge[END_REF], where a method of determination of reduced electric field strength in the positive column of neon DC glow discharge was developed. Basically, comparing the measured spectra of the discharge with the spectra calculated with the CR model by least-squares method, the plasma parameters as reduced electric field strength, electron temperature or electron density may be determined from the best fit.

In this work, Maxwellian electron distribution function (Maxwellian EDF) was assumed for calculation of rates of collisional processes. Other distribution functions, such as, so-called, bi-Maxwellian EDF, general EDF f g , or product of PIC/MC simulation could be included in the calculations. The comparison between the measured and calculated spectra was performed on the basis of 1/I-weighted sum of squares, where I stands for intensity. Since the measured spectra were determined in absolute values of emission coefficient (in W m -3 sr -1 ), no scaling factor between the measured and calculated intensities was needed. On the contrary, a strong dependence of absolute intensity on the electron density increased the sensitivity of the fit to the electron density.

26 transitions of neon were taken for fitting.

In the CR model thirty excited states of neon, arising from 2p 5 3s, 2p 5 3p, 2p 5 4s and 2p 5 3d configurations, were considered. Various collisional, radiative and also diffusion processes were taken into account: electron impact excitation out of the ground state and 2p 5 3s states and deexcitation to these states, electron impact excitation transfer between 3s levels, spontaneous emission of radiation, radiation trapping, electron impact ionization of the ground-state and metastable atoms, chemoionization, associative ionization, 2-body and 3-body collision with ground-state neon atoms or diffusion of metastable-state atoms to the wall. The radiation trapping was incorporated by means of so-called trapping factors (or escape factors, respectively). In their calculations, Voigt profile with Doppler, Stark, natural and resonance broadening components was assumed for lineshapes.

The Einstein coefficients, cross sections, rate coefficients and other data were taken primarily from [START_REF] Navrátil | Collisional-radiative model of neon discharge: determination of E/N in the positive column of low pressure discharge[END_REF]. The rate constants for electron impact induced excitation transfer between 3s states were adopted from [START_REF] Ivanov | Electron-impact-induced excitation transfer between 3s levels of the neon atom[END_REF]. The original sources of Einstein coefficients and oscillator strengths were [START_REF] Sugar | NIST Atomic Spectra Database[END_REF] and [START_REF] Seaton | Oscillator strength in Ne I[END_REF], respectively. For 3s-3p transitions, the accuracy of Einstein coefficients was below 10 % (class B), the calculated oscillator strengths differed from experimental data no more than 10 %, mostly they were 3 % higher [START_REF] Seaton | Oscillator strength in Ne I[END_REF]. This is in agreement with our comparison of both data sources. However, for 3p-3d transitions, with only indirect influence on fitted 3s-3p lines, the differences between the data sources reach 30 -40 %. Accuracy up to 50 % (class D) is reported in [START_REF] Sugar | NIST Atomic Spectra Database[END_REF].

Particle in Cell / Monte Carlo model

One-dimensional (1D) electrostatic PDP1 code [START_REF] Vahedi | Capacitive RF discharges modelled by particle-in-cell Monte Carlo simulation. I: analysis of numerical techniques Plasma Sources[END_REF][START_REF] Verboncoeur | Simultaneous potential and circuit solution for 1D bounded plasma particle simulation codes[END_REF] was modified and employed for the study of various plasma parameters of the RF discharge in neon, e.g. of spatial distribution of plasma density and electron distribution function. The PDP1 code applied the approach of Particle in Cell and Monte Carlo (PIC/MC) methods. A motion of the charged particles in the electric field was solved using the PIC method [START_REF] Birdsall | Plasma Physics Via Computer Simulation[END_REF][START_REF] Hockney | Computer simulation using particles[END_REF] incorporating a self-consistent solution of the Newton motion equations

m i a i = q i E, i = 1, . . . , N,
where m i , a i , q i , E and N denote the particle mass, the acceleration, the charge, the electric field strength and the number of charged particles in the system, respectively, together with the Poisson equation

∇ 2 ϕ = -ρ/ε 0 ,
in which ϕ and ρ stand for the electric field potential and the charge density, respectively.

The equations of motion were solved numerically employing the well-known leap-frog (explicit) method. To reduce a huge number of charged particles in the studied plasma the so-called computer particles were brought into play [START_REF] Birdsall | Plasma Physics Via Computer Simulation[END_REF]. The number of such computer super-particles was then ≈ 10 4 .

Collisions of charged particles with the neutral atoms were considered in our model using a MC method [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, molecule-molecule, and Coulomb collisions for particle modeling of materials processing plasmas and gases[END_REF][START_REF] Vahedi | Capacitive RF discharges modelled by particle-in-cell Monte Carlo simulation. I: analysis of numerical techniques Plasma Sources[END_REF]. Free flight time of charged particles was calculated using the null-collision method. The initial velocities of particles, e.g. after ionization collisions, were also generated using the MC method. Elastic, excitation and ionization collisions were considered for the electrons and the elastic and charge transfer collisions for the ions. The methods briefly described above are nowadays well utilized in computer simulations of plasma [START_REF] Vahedi | Capacitive RF discharges modelled by particle-in-cell Monte Carlo simulation. II: comparisons with laboratory measurements of electron energy distribution functions Plasma Sources[END_REF][START_REF] Longo | Monte Carlo models of electron and ion transport in non-equilibrium plasmas Plasma Sources[END_REF][START_REF] Hagelaar | A Monte Carlo modelling study of the electrons in the microdischarges in plasma addressed liquid crystal displays Plasma Sources[END_REF]].

The model system was set in accordance with the experimental set-up. It consisted of two parallel electrodes having the same diameter of 8 cm. The distance between the electrodes was set at 40 mm. One electrode was powered by RF voltage, the second electrode was grounded. The frequency f of the power source was 13.56 MHz, the voltage amplitude V 0 was 300 V. The pressure p of the neutral gas was 10 Pa or 15 Pa, respectively. Since with the 1D model a negative self-bias could not be modelled, it was introduced into the model as a boundary condition [START_REF] Nagayama | Particle simulation of radio-frequency plasma discharges of methane for carbon film deposition[END_REF].

The boundaries of our computer model were given with the electrodes, where secondary electron emission and recombination occurred. Generally, secondary electron yields for incident electrons and ions, γ e and γ + , respectively, differ. A secondary electron yield γ e depends on electron energy and incident angle. On the other hand a secondary electron yield γ + for low energetic ions, i.e. in case of Auger emission, is dependent only on the type of ion [START_REF] Lieberman | Principles of Plasma Discharges and Material Processing[END_REF][START_REF] Lide | CRC Handbook of Chemistry and Physics[END_REF]. Assuming the a constant values γ for ions and electrons (independent on projectile energy and incident angle) the problem was simplified in our model. The distribution of initial velocities of the secondary electrons was chosen arbitrarily as half-Maxwellian with the temperature kT = 1 eV. The performed tests proved that this choice of initial secondary electron temperature (kT ∈ (0.5 -3) eV) did not influence the results significantly (not presented in this paper).

Results and discussion

Axial dependence

The axial electron temperature dependencies determined in discharge at pressure 10 Pa simulation predicted generally non-Maxwellian EDF, both effective temperature T eff calculated as 2/3 of mean electron kinetic energy, and T high from a slope of EDF at electron energies above 16 eV are displayed. There is obviously a good agreement between OES/CRM method and PIC/MC simulation close to the driven electrode (axial position 0 -7 mm in the plots). Both dependencies of OES/CRM method gave similar values of electron temperature within the error, since calculated and fixed metastable densities were of the same order. There was also a perfect agreement between the spectra measured with OES and spectra calculated with the CR model (see subsection 5.2 for more information). Langmuir probe could not be used in RF sheath. However, going into the bulk plasma, considerable discrepancy between the measurement and simulations was observed. According to PIC/MC simulations, T eff decreases to about 2 eV at positions 20 -25 mm, while T high stays constant at approx. 7 eV. Neither result of OES/CRM nor of Langmuir probe were in agreement with this prediction. Although Langmuir probe registers particularly the development of T eff , the temperature determined from the probe was changing only slowly around 5 eV. The optical spectra depend mostly on T high , but the electron temperature values determined from OES/CRM method were considerably higher. Also the calculated spectra did not satisfactorily fit the measured spectra in the bulk plasma. Although PIC/MC simulation predicts electrons with a temperature about 10 -20 eV at similar discharge conditions due to stochastic heating (see figure 3 and also next), they are expected just conversely close to the driven electrode.

The differences between Langmuir probe measurement and the PIC/MC model can be understood comparing directly the EDFs (see figure 4). PIC/MC simulation generally predicts bi-Maxwellian, or a more complicated EDF. Closely to the electrode, PIC/MC expects only electrons with higher temperature, T high ≈ 7 eV, and the EDF is Maxwellian. Going into the bulk plasma a group of low energetic electrons with temperature ≈ 0.7 eV appears, but T high remains the same, though the density of electrons with this temperature decreases. This behaviour of EDF is responsible for the development of T eff observed in figure 2: T eff approaches T high at the electrodes, but decreases in the middle. The EDF measured by Langmuir probe was almost Maxwellian in all measured positions (see figure 4b). An only weakly steeper low-energy part of EDF was observed in the discharge centre (positions 20 -28 mm). Consequently, the probe measurements did not confirm the bi-Maxwellian character of the EDF. However, the determination of the low-energy part of EDF by Langmuir probe is known to suffer from several problems [START_REF] Godyak | Measurement of electron energy distribution in low-pressure RF discharges Plasma Sources[END_REF]. A high amount of low-energy electrons, present in bulk plasma, may be unregistered by the probe. For these reasons, the temperature determined from the probe measurement (see figure 2) does not follow the axial dependence of T eff from PIC/MC simulation with pronounced temperature minimum in the bulk plasma.

The bi-Maxwellian EDF was observed firstly by Godyak et al [START_REF] Godyak | Abnormally low electron-energy and heating-mode transition in a low-pressure argon RF discharge at 13.56[END_REF] in argon, employing the Langmuir probe. The origin of bi-Maxwellian EDF is relatively simple.

The high-temperature group of electrons gain their energy in so-called stochastic heating process, in which interactions of electrons with the oscillating plasma sheaths are crucial. This group of electrons has enough energy to enter into inelastic collisions with neutral atoms, furthermore these electrons can easily overcome the ambipolar potential barrier interacting therefore more frequently with the oscillating sheath edges [START_REF] Godyak | Abnormally low electron-energy and heating-mode transition in a low-pressure argon RF discharge at 13.56[END_REF].

On the other hand, electrons created in the ionization collisions inside the bulk of the plasma belong to the group of low-temperature electrons. These electrons oscillate in very weak electric field presented in the bulk of the plasma and hence cannot gain energy to overcome the ambipolar potential barrier on the boundary of plasma sheath and plasma bulk, where there exists a maximal electric field and where stochastic heating takes place [START_REF] Lieberman | Principles of Plasma Discharges and Material Processing[END_REF][START_REF] Godyak | Abnormally low electron-energy and heating-mode transition in a low-pressure argon RF discharge at 13.56[END_REF][START_REF] Brzobohatý | Influence of the Ramsauer minimum on the plasma characteristics studied via computer simulation Czech[END_REF][START_REF] Brzobohatý | Influence of substrate material on plasma in deposition/sputtering reactor: experiment and computer simulation[END_REF]. Thermalization of these electrons occurs only due to the elastic collisions with neutral atoms, which become important at higher pressures. However, no substantial improvement was obtained. Bi-Maxwellian EDF gave a similar value of electron temperature in high energy region to Maxwellian EDF. Moreover, it was not possible to let the fitting algorithm change low energy temperature independently, since optical spectra are less sensitive to this parameter through stepwise excitation and deexcitation processes. General distribution function was found to describe the electrons better (evaluated by sum of squared differences in measured and fitted spectra), but only in electrode vicinity, where the deviation from Maxwellian distribution was small. As expected, the use of EDF from PIC/MC simulation did not produce optical spectra similar to the spectra observed by OES.

Influence of metastables

In agreement with our previous work [START_REF] Navrátil | Collisional-radiative model of neon discharge: determination of E/N in the positive column of low pressure discharge[END_REF], the calculated intensities of 3p -3s transitions depended considerably on the populations of metastables. The influence of metastables (and of all states of 3s configuration) on the populations of 3p states is through a stepwise electron impact excitation and it is crucial at low electron energies. At higher electron energies the direct excitation from the ground state becomes the most important. This provides a sensitivity of the spectra to the electron temperature, although excitation energies of 3p states differ by 0. at most. On the other hand, this requires a reliable cross-sectional data for electron impact excitation out of the ground state to 3s states and out of 3s states. In this work a cross section for 1s 5 -2p 9 transition from [START_REF] Boffard | Measurement of electronimpact excitation cross sections out of the neon 3 P 2 metastable level[END_REF] was scaled with oscillator strength to obtain a consistent cross-section set for all optically allowed excitations out of 3s states [START_REF] Navrátil | Collisional-radiative model of neon discharge: determination of E/N in the positive column of low pressure discharge[END_REF]. The uncertainty of the collision cross sections is about ± 30 % as determined in [START_REF] Boffard | Measurement of electronimpact excitation cross sections out of the neon 3 P 2 metastable level[END_REF], but the results of various authors somewhen differ by as much as a factor of 3 (compare [START_REF] Boffard | Measurement of electronimpact excitation cross sections out of the neon 3 P 2 metastable level[END_REF][START_REF] Behnke | Investigation about stepwise excitation cross sections in rare gases Contrib[END_REF][START_REF] Leveau | Destruction par chocs électroniques des atomes métastables et pseudo-métastables de néon dans une colonne positive[END_REF]). Moreover, reliable calculation of self-absorption of resonant lines and also of lines originating from transitions between 3p -3s states, is needed.

The dependence of optical spectra on axial position is shown in figure 5. Going into the bulk plasma the intensities of individual lines increase at different growth rates.

Whilst at axial position 4 mm line at 640 nm is the most intensive, the line at 585 nm dominates in spectrum at larger axial distances. Since upper state of line at 585 nm (2p 1 in Paschen notation) is populated mainly by direct electron impact excitation from the ground state and upper state of line at 640 mm (2p 9 ) is populated substantially also by stepwise excitation, the increase of intensity ratio 585 nm/640 nm is interpreted by the CR model as enhancement of electron energies in the discharge middle. However, this is in strong contrast with the decrease of electron temperature predicted by PIC/MC simulation.

The high values of electron temperature determined by OES/CRM method in the discharge middle suggest, that the relative intensities in spectra are influenced by other processes. E.g., the spectra may be influenced by the absorption of light in the outer regions of the discharge and in the vacuum chamber. Metastable atoms diffuse to the chamber wall and can selectively absorb intensities of transitions ending in metastable states. Since these transitions in opposite direction are also stimulated by electron impact, their absorption may be misinterpreted as enhancement of direct excitation against the stepwise excitation. A simulation of spatial distribution of metastables in the discharge chamber is needed to analyse the effect of absorption. However, this distribution cannot be determined using the local CR model. For measured metastable densities the escape factor of line 640 nm is ≈ 0.7. The local CR model predicts even higher metastable densities and the escape factor can reach values ≈ 0.1. Such small values of escape factors show a large influence of the metastables on the line intensities.

Limitations of PIC/MC code The 1D PIC/MC model is generally valid in cases

where the separation between the electrodes is much smaller than electrode diameters.

In our case the separation was 40 mm and the diameters of electrodes were only 80 mm. Furthermore, the 1D model is symmetric (both electrodes have the same area) and therefore no self-bias can be observed. This drawback of our PIC/MC simulation code was overcome by introducing the bias voltage artificially. The value of the bias was set to -100 V for pressure 10 Pa in simulations, which was below the experimental value. At -220 V the simulation was unstable. However, as can be seen from figure 3, the bias mostly affects the maximal electron temperature ‡ of beam electrons, but the temperature in the bulk remains approximately the same, the minimum only shifts in axial direction.

The discrepancy between measured and simulated EDFs in bulk plasma (the measured EDF was Maxwellian, but the simulated EDF was bi-Maxwellian) may be due to simple incorporation of inelastic collisions. Only total excitation into 3s state was considered in the model from excitations. E.g., superelastic collisions can reasonably modify the shape of EDF. However, the influence of higher excited states on EDF is usually negligible. For better agreement with the experimental data the more complex 2D or 3D model should be employed.

The comparison of determined electron density in discharge at pressure 10 Pa is shown as a function of axial position in figure 6. The density dependence is different from the intensity development due to spatially varying electron temperature. The relative axial dependencies of electron density determined by OES/CRM and Langmuir probe are rather similar. However, they differ approx. 5× in their absolute values. This can be due to complicated estimation of spatial dependence of emission coefficient. According to optical and probe measurement of radial dependencies of involved parameters (electron temperature and density, discharge intensity) a constant radial profile was assumed for the emission coefficient of plasma region with 110 mm in diameter (electrode diameter was 80 mm). However, this assumption is only approximate and at 10 Pa there is still some radiation produced outside this region, which increases the measured intensity ‡ We keep using effective temperature (k b T = 2/3E mean ), although E mean would be more appropriate in this discussion. The electron density simulated by PIC/MC method reaches values between OES/CRM and Langmuir probe results in the bulk plasma. The spatial profile of electron density simulated by PIC/MC method is sensitive to the self-bias, which is given due to the asymmetry of the electrodes. Better agreement would be achieved by employing a more complex 2D code.

Influence of temporal development

The axially and temporally resolved, but spectrally unresolved, light emission of the discharge measured with the ICCD camera is displayed in figure 7. The light emission is axially and temporally inhomogeneous, although the variation of light intensity through the period is substantially smaller in neon at 10 Pa than at lower pressure [START_REF] Schulze | Electron beams in asymmetric capacitively coupled radio frequency discharges at low pressures[END_REF][START_REF] Schulze | Electric field reversals in the sheath region of capacitively coupled radio frequency discharges at different pressures[END_REF] or in other gases [START_REF] Schulze | Electric field reversals in the sheath region of capacitively coupled radio frequency discharges at different pressures[END_REF]. In our case the light emission lasted nearly the whole period. Since the light emitted from the discharge was accumulated on the CCD detector over many periods of RF signal, the obtained results of OES/CRM were averaged over the whole period. E.g. the values of electron temperature and density measured at positions close to electrode do not reflect the conditions in sheath only, but rather the plasma flooding the region when the sheath collapses. However, the electron temperature and density calculated from temporally integrated intensities need not be simply equal to their temporal mean values, determined e.g. in PIC/MC simulations.

Power dependence

The agreement of measured and calculated spectra and OES/CRM and PIC/MC method at positions close to the driven electrode was further analysed by decreasing the incoming RF power. The optical emission spectra were measured at position 4 mm from the driven electrode, the RF power was changed in the range 5 -50 W. The metastable density was calculated from the CR model.

The intensities of spectral lines increased monotonously with increasing power. However, they increased in a different way. The dependencies of two lines at 585 nm and 640 nm on RF power are displayed for pressure 10 and 15 Pa in figure 8. It can be seen, that whilst the 640 nm line increases nearly linearly with power, the increase of 585 nm line slightly decelerates. This behaviour can be understood considering figure 9 and figure 10, in which dependencies of electron density and electron temperature determined from the spectra are shown (taken at the same axial position 4 mm). It can be seen that with increasing power the electron density increases. This is why intensities of spectral lines increase. Due to the increasing conductivity of the plasma, electrons gain less energy, which probably causes the decrease in electron temperature. Then, decreasing electron temperature more affects the excitation of those higher excited states, which are not populated by stepwise excitation through the metastable levels but mostly by direct excitation. This is the case of e.g. 2p 1 state, from which 585 nm line originates. The examples of fitted spectra for pressure 10 Pa and input RF power 10 W and 50 W, are shown in figure 11 and figure 12, respectively. Obviously, the agreement between the calculated and measured spectra is appreciably good. This is evidence of a valid assumption of Maxwellian electron distribution function and correct description of kinetics in this region.

Conclusion

In

Figure 1 .

 1 Figure 1. The scheme of experimental set-up. DE & GE -driven and grounded movable electrode, G -13.56 MHz generator, O -oscilloscope, OF -optical fibre, ID -iris diaphragms, OES -Jobin Yvon FHR 1000 spectrometer, LP -Langmuir probe.

Figure 2 .

 2 Figure 2. The electron temperature vs. axial position for pressure 10 Pa as determined by various methods: PIC/MC T eff , T high -effective temperature determined from mean electron kinetic energy and slope of EDF at electron energies above 16 eV, respectively, Langmuir probe measurement, CR model with metastable density determined from the model, CR model with metastable density set according to absorption measurement. Axial profile of intensity of 585 nm line (dotted line) is shown to clarify the axial position: 0 mm -driven electrode, 40 mm -grounded electrode.

Figure 3 .

 3 Figure 3. The electron temperature determined from PIC/MC model for pressure 10 and 15 Pa and bias -100 and -150 V, respectively. The temperature is only effective and it is calculated from the mean electron kinetic energy.

Figure 4 .

 4 Figure 4. Electron distribution function at different axial positions, as determined from (a) 1D PIC/MC simulation for pressure 10 Pa and bias -100 V and (b) Langmuir probe measurement. The EDF determined in PIC/MC simulation was Maxwellian at the electrodes and departured from this distribution going into the bulk plasma. The measured EDF was Maxwellian in all measured positions.

5. 1 . 1 .

 11 Influence of EDF type It is interesting that OES/CRM and PIC/MC methods gave similar temperature values only in electrode region, where EDF was Maxwellian according to PIC/MC simulations. Considerable attention was therefore paid to describe the electrons of the bulk plasma in CR model by more sophisticated EDFs. Bi-Maxwellian, general distribution function and also directly the EDF produced by PIC/MC were tested in the CR model to explain the discrepancy between the measured optical spectra and the spectra calculated for temperature expected by other methods.

Figure 5 .

 5 Figure 5. The dependence of measured optical spectra of neon RF discharge at pressure 10 Pa and incoming power 50 W on axial position. The line intensities were integrated over their lineshapes. Note the different scales of vertical axes.

Figure 6 .

 6 Figure 6. The electron density vs. axial position as determined by various methods: PIC/MC simulation, OES/CRM method with metastable density set according to the absorption measurement and Langmuir probe measurement. Axial profile of intensity of 585 nm line (dotted line) is shown to clarify the axial position: 0 mm -driven electrode, 40 mm -grounded electrode.

Figure 7 .

 7 Figure 7. The axially and temporally resolved light emission of RF discharge in neon at 10 Pa, measured with ICCD camera during RF voltage period. The intensity was measured in the discharge centre. The driven electrode is at position 0 mm, grounded electrode at position 40 mm.

Figure 8 .

 8 Figure 8. The intensity dependences of two lines at 585 nm and 640 nm on RF power at position 4 mm from the driven electrode.

Figure 9 .

 9 Figure 9. The electron density determined from the spectra taken at position 4 mm from the driven electrode as a function of power.

Figure 10 .

 10 Figure 10. The electron temperature determined from the spectra taken at position 4 mm from the driven electrode as a function of power.

Figure 11 .

 11 Figure 11. The example of spectra fit at position 4 mm and input RF power 10 W. The electron temperature and density determined from the fit were 12.2 eV and 1.2 × 10 8 cm -3 , respectively.

Figure 12 .

 12 Figure 12. The example of spectra fit at position 4 mm and input RF power 50 W. The electron temperature and density determined from the fit were 8.4 eV and 6.2 × 10 8 cm -3 , respectively.

Acknowledgments

The work was supported by research project MSM0021622411 of Ministry of Education of the Czech Republic, grants No. 202/09/0800 and No. 202/07/1669 of Czech Science Foundation and by project No. KAN 101630651 of Academy of Sciences of the Czech Republic.

with increasing incoming RF power. The measured and calculated spectra agreed very well in this region.

The obtained results suggest that the determination of electron density based on OES with absolute intensity measurement and CR modelling suffers from serious problems. Complicated geometry, unknown radial intensity profile, reflection of light in the chamber, etc. -all of these attributes make the estimation of absolute emission coefficient and thus also of electron density difficult. On the other hand, the electron temperature is determined mainly from the relative shape of optical spectra, which is sensitive mostly to high energetic region of the EDF. Although this sensitivity in low energetic region of EDF is increased by the stepwise excitation, fitting problems appear when more free parameters are used to describe the shape of EDF.

A very important topic is the correct determination of absolute value and spatial profile of neon metastables (or all 3s states). They are not only the initial states, out of which the low energetic electrons populate the studied 3p levels, but they can also alter the spectrum measured by spectrometer due to the radiation trapping. This is in contrast e.g. with the TRG-OES method, in which the studied emissions may be assumed to be under optically thin conditions [START_REF] Donnelly | Plasma electron temperatures and electron energy distributions measured by trace rare gases optical emission spectroscopy[END_REF].


Iva Dostálková

Marek Špinka
email: spinka@vuzv.cz

synchronization of all-or-nothing activity transitions in grouped animals

Keywords: group behaviour, modelling of behaviour, synchronization, group cohesion, group size

For groups of animals to keep together, the group members have to perform switches between staying in one place and moving to another place in synchrony. However, synchronization imposes a cost on individual animals, because they have to switch from one to the other behaviour at a communal time rather than at their ideal times. Here we model this situation analytically for groups in which the ideal times vary quasi-normally and grouping benefit increases linearly with group size. Across the parameter space consisting of variation in the grouping benefit/cost ratio and variation in how costly it is to act too early and too late, the most common optimal solutions are full synchronization with the group staying together and zero synchronization with immediate dissolution of the group, if the group is too small for the given benefit/cost ratio. Partial synchronization, with animals at the tails of the distribution switching individually and the central core of the group in synchrony, occurs only at a narrow stripe of the space. Synchronization cost never causes splitting of the group into two as either zero, partial or full synchronization is always more advantageous. Stable solutions dictate lower degree of synchrony and lower net benefits than optimal solutions for a large range of the parameter values. If groups undergo repeated synchronization challenges, they stay together or quickly dissolve, unless the animals assort themselves into a smaller group with less variation in the ideal times. We conclude with arguing that synchronization cost is different from other types of grouping costs since it does not increase much with increasing A c c e p t e d m a n u s c r i p t 2 group size. As a result, larger groups may be more stable than smaller groups. This results in the paradoxical prediction that when the grouping benefit/grouping cost ratio increases, the average group sizes might decrease, since smaller groups will be able to withstand synchronization challenges.

INTRODUCTION

Many animals live in stable or semi-permanent groups or at least spend shorter or longer periods of time in loose aggregations. In order to stay in proximity, group members need to coordinate when to stay in one place and when to travel. In other words, each animal has to strike a compromise between switching behavioural activities in close synchrony with others and switching the activities as close as possible to its individual ideal time (e.g, [START_REF] Ruckstuhl | To synchronise or not to synchronise: a dilemma for young bighorn males?[END_REF]Meldrum and Ruckstuhl, 2009).

One example of such situation is fleeing from a predator. In many instances, the predator is noted well ahead of possible attack (either directly by the animal itself or through the alert behaviour of other flock/herd members, [START_REF] Sirot | Coordination and synchronization of vigilance in groups of prey: the role of collective detection and predators' preference for stragglers[END_REF] and then an "economics of fleeing" comes into effect [START_REF] Ydenberg | The economics of fleeing from predators[END_REF]. Initially, the costs of staying are minute (because the risk of a successful attack from large distance is negligible) and therefore lower than the costs of fleeing (including energy for the flight and lost opportunity to continue with foraging or resting). As the predator approaches, costs of staying increase at an accelerating rate (because of the climbing probability of actual attack and risk of being caught) while the costs of fleeing remain stable or increase only slightly (e.g. the animal has to run faster or for longer). At the moment when the costs of staying overcome the costs of fleeing, the animal should depart. Because many factors affecting both types of costs (such as individual ability to evade attack, current position in the group, satiation/energy reserves) vary between individuals, the ideal departure time is different for each animal. If everybody sticks to its ideal timing, animals will depart individually, the group will disintegrate and the benefits of staying and behaving together will be lost. On the other hand, should the group stay fully together and enjoy the associated benefits, everybody has to make a smaller or larger concession from its ideal timing. This dilemma is present in every situation where grouped animals can draw a benefit from switching synchronously from one behaviour

A c c e p t e d m a n u s c r i p t

3 to another such as departure from a resting place, stopping migration to rest, switching from monitoring a predator to fleeing. However, departure from a place is perhaps the most typical example and therefore the model will be described in terms of who leaves when.

Synchronization can be viewed as a consensus decision making [START_REF] Conradt | Consensus decision making in animals[END_REF][START_REF] Conradt | Group decisions in humans and animals: a survey[END_REF]. Several theoretical and empirical studies have documented that animals can "agree" on the common solution using simple behavioural rules [START_REF] Couzin | Effective leadership and decisionmaking in animal groups on the move[END_REF][START_REF] Amé | Collegial decision making based on social amplification leads to optimal group formation[END_REF][START_REF] Martín | Collective detection in escape responses of temporary groups of Iberian green frogs[END_REF]. If there are just two or a few disjunct possible solutions (e.g., a choice between distinct spatial goals or shelters) then simple local rules can be found that guide the group either to one of the goals, or lead to splitting of the group [START_REF] Amé | Collegial decision making based on social amplification leads to optimal group formation[END_REF][START_REF] Couzin | Effective leadership and decisionmaking in animal groups on the move[END_REF]. Also, for groups of two animals and a continuous range of possible preferences (e.g., times of a behavioural switch), models have been developed that combine the mechanisms of coordination (e.g., communicating about inner state, simple reaction to the behaviour of the other animal) with their fitness outcomes [START_REF] Rands | Spontaneous emergence of leaders and followers in foraging pairs[END_REF][START_REF] Dostálková | Synchronization of behaviour in pairs: the role of communication and consequences in timing[END_REF]. Situations where there is a continuum of preferences in larger groups have been much less investigated. The only two previous studies to address these questions in terms of costs and benefits are those by Conradt and Roper. In the first paper, [START_REF] Conradt | Group decision-making in animals[END_REF] show that agreeing democratically on the time of the switch is usually more beneficial for the group than yielding to the time enforced by one leader animal, even if that animal has more experience or better judgement than the others. In their next study, [START_REF] Conradt | Democracy in animals: the evolution of shared group decisions[END_REF] document that for groups of three animals, agreeing on any of the three times preferred by the three group members is evolutionary stable. They then expand this analysis to larger groups and argue that also for those, any shared decision (e.g., even departing when the first animal wants to depart) is evolutionary stable. This is due to the fact that if a small subset of the group chooses to depart at another time, it will acquire some benefit from more appropriate timing (i.e. it will decrease its consensus costs), but it loses the whole benefit of acting together with the original large group. This argument holds only if the costs are small relative to grouping benefits. [START_REF] Conradt | Democracy in animals: the evolution of shared group decisions[END_REF] state this verbally but stop short of analysing the balance between consensus costs and grouping benefits quantitatively and that is where our study carries on.

The trade-off can result in four types of situations: either all animals "agree" on a common departure time (full synchronization); or each animal leaves individually (no synchronization); or something in between occurs with some individuals departing on their own while the rest leaves at a common time (partial synchronization); and lastly, animals may

A c c e p t e d m a n u s c r i p t

4 be also leaving in smaller subgroups (group splitting). Our model aims at answering three questions. First, how probable are the cases of full, partial and zero synchronization, and of group splitting? How do these probabilities depend on the benefit of being synchronized, and the costs departing too early or too late? Second, are the animals in a group able to cooperate on an optimal synchronization or are the evolutionary stable strategies for this synchronization problem suboptimal in terms of resulting net benefits? Third, the sychronisation of group departures may have an important influence on group stability and therefore stable group size. Therefore we also asked: Will repeated synchronization challenges lead to group stabilization or disintegration?

MODEL I DESCRIPTION: ONE COMMON DEPARTURE VERSUS INDIVIDUAL

DEPARTURES

The detailed elaboration of the model and derivation of its results is given in Appendix A. We suppose that a group of N animals is about to switch from one behaviour to another (e.g., to leave a resting place or to switch from monitoring a predator to fleeing). The animals in the group have different ideal times for the switch. These individual ideal times * i t are symmetrically distributed around zero and constrained to the interval (-t max , t max ). For simplicity we assume (-t max , t max ) = (-1, 1) (Fig. 1). For most of the animals, the individual ideal times are not far from t * = 0 while only a negligible proportion of animals prefer times immediately after t * = -1 or just before t * = 1. In other words, the probability of the position of ideal time t * for a randomly chosen animal is the highest for t * 0 and it is close to zero for t * near -1, 1 . . We assume that concrete ideal times t * are realizations of continous random variable "ideal time" defined by its density function f.

The model is setup in terms of grouping benefits and synchronization costs.

An animal that departs at a moment t i when a fraction p(t i ) (0 ≤ p(t i ) ≤ 1) of the N-sized group depart gets a grouping benefit proportional to the number of animals in the departing group, p(t i ) NQ, where Q is a constant. If an animal departs too early, i.e. at a time t i before its ideal

time * i t , it is penalised by -(1-ω)L ( * i t -t i ); if it leaves too late at t i > * i t , it receives the penalty of -ωL (t i -* i t
). The parameter ω (0 ≤ ω ≤1) describes whether it is more costly to leave prematurely (ω < 0.5) or rather with a delay (ω > 0.5). The product NR (where R = Q/L and is to be in synchrony with others, relative to the importance of departing close to the individual ideal time. The ratio R is multiplied by the group size N because the model assumes that the larger the group, the higher the grouping benefit. Thus the product NR quantifies how much an animal can get from synchronization: the higher the parameter NR, the higher net potential benefit can the animals cash from being staying with the group.

The model examines how an animal should behave (i.e., when to leave) in order to accrue the highest possible net benefit, given its specific ideal time * i t and the parameters NR and ω. Parameters NR and ω are equal for all animals in the group. Group size N is considered fixed for this model. We consider two strategies: either to depart individually at its own ideal time, or to leave together with other animals in one mass departure (Fig. 1). The solitary departure avoids any penalization for improper timing but gets no benefit from synchronization, hence its net benefit is zero. The animals participating in the mass departure benefit from synchronization (in direct proportion to the percentage of animals that leave with them), but pay for deviating from their ideal times. The cost of participating in the mass departure are higher for animals at the tails of the distribution and therefore, animals at the tails should be more inclined to use the solitary strategy while animals with their t* i closer to the centre of the distribution should rather join the mass departure (Fig. 1). We suspect that when ω < > 0.5, the time of mass departure t D will be different from zero and also the times t 1 and t 2 (breaking points where the solitary changes into the mass departure strategy and back) will be asymmetrically positioned. In the solution of the model, we are searching for such combination of times t 1 , t D and t 2 that brings the highest net benefit, depending on the parameters NR and ω. If t 1 = t 2 , then each animal departs on its own and there is no synchronization. If -1 < t 1 < t 2 < 1, then some animals depart solitary and some in the mass departure, and we label this partial synchronization. And finally full synchronization occurs when t 1 = -1 and t 2 = 1.

The net benefit function B(t 1 , t D , t 2 ) is the sum of the grouping benefit function B G (t 1 , t D , t 2 )

and timing cost function B C (t 1 , t D , t 2 ).

B(t 1 , t D , t 2 ) = B G (t 1 , t D , t 2 ) + B C (t 1 , t D , t 2 )
, where

B G (t 1 , t D , t 2 ) = NR ( 2 1 ) ( t t dx x f ) 2
A c c e p t e d m a n u s c r i p t

6 B C (t 1 , t D , t 2 ) = - D t t D dx x f x t 1 ) ( ) ( -(1-) 2 ) ( ) ( t t D D dx x f t x
Thus, the individually departing animals do not get any grouping benefits and do not pay any timing (synchronization) costs, whereas for animals leaving in the mass departure, the net benefit equals the grouping benefit minus the synchronization costs due to un-ideal timing of the departure.

We look for two types of solutions: an optimal and a stable one. The optimal solution is such that an animal acquires, on average, the highest summary net benefit, whereas the stable solution is a solution whose net benefit cannot be improved by unilateral adoption of another solution by a small subset of the animals.

Therefore optimization procedure means to find triplets [t 1 , t D , t 2 ] depending on NR and so that the value of net benefit function B is maximal. Using the classical metod for maxima finding we obtain the set of maxima ( Similarly to the optimal solutions, there is a croissant-shaped area (delimited by dashed lines in Fig 3) within which partial synchronization is a stable strategy. Above the croissant, full synchronization is stable. The lower boundary for partial synchronization and the boundary between partial and full synchronization are both shifted upwards (specifically, to twice as high NR) in comparison with the optimal solution. That is, the stable solution dictates a lower degree of synchronization than the optimal solution in a large area of the parameter space, including a place where the optimal solution prescribes full synchrony while the stable solution is no synchrony at all (Fig. 4). In this area, the stable strategy brings substantially lower net benefit than the optimal one (Fig. 4). Why this is so? An animal should participate in the mass departure only if the benefit outweighs the cost. The cost of improper timing is always paid solely by the specific individual, but there are two kinds of benefits: one goes to the animal itself (as it is part of the synchronized subgroup) and another goes to every synchronized animal because through the contribution of the focal animal, the synchronized

A c c e p t e d m a n u s c r i p t

8 subgroup increases in size. In the optimal solution, where the costs and benefits are averaged over the group, both types of benefits enter the equation and they are able to pay for higher timing cost, thus drawing more of the animals from the tails into the synchronization. In the stable solution, only the first type of benefit counts and thus more animals will opt out the synchrony.

The second difference between the stable and the optimal solutions is that zero synchronization is a stable strategy across all the parameter space, and not only below the croissant. In the areas of higher and high NR, i.e. within and above the croissant, two stable strategies coexist: the partial/full synchronization and the zero synchronization, of which the former brings higher net benefit. Zero synchronization is stable because once it becomes established, each animal departs at a different time and therefore it is difficult for a small fraction of animals to set foot with an alternative strategy, as their attempts to synchronize get no support from the majority.

Splitting into subgroups

Splitting into two subgroups is never an optimal strategy because for any given NR, the optimal strategy with one central departure always brings higher net benefit than the best twodeparture option (Fig. 5).

Repeated synchronizations

What happens when the synchronization challenge is repeated? If the grouping benefit/consensus cost ratio (NR) is high, then no animals are shed and the group remains of the same size across a single or repeated challenges. Low NR leads to zero synchronization and therefore to immediate dispersal of the group on the first challenge. The interesting case occurs with partial synchronization during which the group size is trimmed (Fig. 6). What happens on the next encounter depends on what is the distribution of ideal times in the downsized group. To simplify the solution, we restrict the analysis of this problem to symmetric costs of time (corresponding to ω = 0.5 in the previous model) and to the optimal solution.

If each animal keeps its ideal time from one challenge to the next then during partial synchronization both the group size and the variability of the density distribution of ideal

A c c e p t e d m a n u s c r i p t

9 times is also trimmed (Fig. 6b), but the latter more so because the costly tails of the distribution are shed. Therefore, the grouping benefit/consensus cost ratio increases and the group size stabilizes at a reduced size. The derivation of this result is given in Appendix C.

We call this option the Assortment scenario as it results in a subgroup where members are more alike in terms of their time preferences than in the original group.

If partial synchronization leads to a smaller group in which the variation of ideal times remains as large as before the reduction (i.e. individuals do not keep their individual times from one challenge to another, Non-Assortment scenario), then the consensus costs remain unaltered while grouping benefits fall (Fig. 6c). Consequently, the group is reduced further during next challenge or challenges until it disintegrates (see Appendix C).

Minimal group size

It follows from the analysis of repeated synchronizations that in Fig. 2, that the upper border of the partial synchronization area delineates the minimal R value for which the group of size N remains cohesive in the Non-Assortment scenario, and the lower border demarcates group cohesion boundary for the Assortment scenario. In the analysis of Models I and II, we have assumed that the group size was fixed, and the relative "per animal" grouping benefit R was varied. Conversely, the chart can be understood as illustrating how the stability of a group depends on its size N, for a given "per animal" grouping net benefit R (Fig. 7).

DISCUSSION

Our models yield three main insights into the links between the grouping-benefitversus-consensus-costs trade-off, behavioural synchronization and group size and stability.

First, the results indicate that synchronization of departures will often be an all-ornothing phenomenon. Model predicts that if the benefit of staying in a group of animals depends linearly on the size of the group (which is the assumption of the model) and the cost of staying in the group varies quasinormally among the group members, then the group will most of the time either switch behaviour synchronously and thus stay together, or all individuals will behave independently and the group will disintegrate. The cases when animals at the tails of the distribution will act solitarily and the rest of the group in synchrony are predicted to be rare. In other words, the overall grouping-benefit-to-consensus-costs ratio

A c c e p t e d m a n u s c r i p t

10 is predicted to decide, in most cases, for all animals to either depart synchronously or individually, with individual differences such as boldness/fearfulness [START_REF] Reale | Integrating animal temperament within ecology and evolution[END_REF] or satiation [START_REF] Rands | State-dependent foraging rules for social animals in selfish herds[END_REF] playing only a marginal role in this situation. Moreover, Model II shows that splitting into two subgroups brings lower net benefit than the central synchronization in one group, indicating that synchronization challenges will not lead to separation into smaller units. Finally, the analysis of repeated synchronization events with no assortment shows that even when partial synchronization occurs, the group will quickly diminish and disintegrate.

Second, Model I shows that when individuals in groups behave selfishly, a much weaker synchronization, higher probability of group disintegration and lower net benefits are predicted than in groups that behave optimally from the overall fitness point of view. This is due to the fact that the cost of staying with the group is carried individually by each animal, whereas benefits from the presence of the individual in the group are cashed both by the animal in question (it enjoys the company of others) and by all other group members (because the groups stays larger). Therefore, if the animal counts in only its own benefits (the stable solution), it is prepared to expend less costs to stay in the group than if everybody's benefits are taken into account. Moreover, Model I documents that zero synchronization is a stable solution all across the parameter space. Why zero synchronization is stable everywhere and full synchronization is not? Full synchronization is unstable with low NR since any individual can, independently of what others do, default from the synchrony, leave at its individual ideal time and cash in the benefit of not paying any penalty for improper timing. In the contrary, skipping individually from zero synchronization to a partial or full synchronization just by individual action is impossible, since in order to achieve synchrony, you need pals to go with you. The existence of zero synchronization as a solution for the entire parameter space in the stable, but not in the optimal solution further underscores the general implication of our models that evolutionary stable regime of synchronized departures will often bring substantially lower net benefit than the optimal one.

The optimal-stable solution dichotomy shows that there is a Prisoners Dilemma type of game underlying our models [START_REF] Fletcher | The Evolution of Altruism: Game theory in multilevel selection and inclusive fitness[END_REF]. If an animal plays Stable against a group that plays Optimal, it gets a better pay-off than if it would play Optimal. However, as animals shift from the Optimal to the Stable strategy, the whole group gradually starts playing Stable, and everybody ends with a worse pay-off than it would get if the whole group kept playing Optimal. The combination of our two models thus predicts that sub-optimally low

A c c e p t e d m a n u s c r i p t

11 degree of synchronization will be frequently observed in departures and other major cohesion influencing behavioural switches unless the optimal strategy could be maintained by altruismsupporting mechanisms such as relatedness, reciprocity or assortment (Fletcher andZwick, 2006, 2007;[START_REF] Frommen | shoaling decisions in threespined sticklebacks (gasterosteus aculeatus) -familiarity, kinship and inbreeding[END_REF][START_REF] Fischhoff | Reproductive status influences group size and persistence of bonds in male plains zebra A c c e p t e d m a n u s c r i p t 14 (Equus burchelli)[END_REF].

Third, and perhaps most interestingly, the current study brings to light the phenomenon of minimal stable group size. To our knowledge, this is the first study pointing to this phenomenon and identifying synchronization of behaviour as its driving force. Our results indicate that the need to synchronize the starts and stops of movements will in many cases set a minimal threshold group size below which group cohesion will no longer be beneficial. The situation is especially clear in the non-assortment scenario, where groups above a size that makes full synchronization stable are predicted to stay together while smaller groups should disintegrate. Thus the model predicts that when per capita grouping benefits are high, the average group size should decrease because small groups will be able to withstand synchronization challenges.

How does this seemingly paradoxical prediction arise? The synchronization cost differs substantially from other types of grouping costs considered so far. In most of the models and reviews published to date, the costs to an individual of being in a group were thought in terms of foraging interference (Stillmann et al., 2000;Sansom et al., 2007), reproductive interference [START_REF] Westneat | Density and extra-pair fertilizations in birds: a comparative analysis[END_REF], proximity-transmitted pathogens (Côté and Pohlin, 1995) or similar factors. For instance, in their influential book, Krause and Ruxton (2002, page 41) list twelve different types of grouping costs, but do not mention synchronization cost. The costs listed by Krause and Ruxton were assumed to increase with increasing group size, often at an accelerating rate (Fig. 8a). In contrast, the synchronization cost in our model remains constant, irrespective of group size (Fig. 8b). This is because the synchronization cost to an individual is given by the time distance between its ideal time and the time of common departure and this distance neither shrinks nor expands if more or fewer other animals also pay their costs. Previous models assumed that even very small grouping is better than solitary life and sought to find the optimal and the maximal stable group size (Fig. 8a). Our model shows that when the need to synchronize is the major cost of grouping, then net benefit may be lower in small groups than in solitary animals and a certain group size is needed for the grouping benefits to offset the synchronization costs (Fig. 8b). It remains to be investigated how important role the movement synchronization costs play in determining animal group sizes in nature. In some cases, it is possible that the synchronization and other

A c c e p t e d m a n u s c r i p t

12 types of grouping costs will combine so that the resulting cost function, together with depreciating grouping benefits, will set both minimum and maximum profitable group size (Fig. 8c). Our results also show that smaller groups are stable if animals assort according to a feature that determines their preferred departure times such as boldness/shyness [START_REF] Croft | Behavioural trait assortment in a social network: Patterns and implications[END_REF] or body size [START_REF] Ruckstuhl | To synchronise or not to synchronise: a dilemma for young bighorn males?[END_REF][START_REF] Hoare | Body size and shoaling in fish[END_REF].

The predictions of our models are linked to two crucial assumptions, namely that the grouping benefit increases linearly with the group size, rather than according to a depreciating function, and that the heterogeneity in ideal departure times is unimodally distributed. If the first assumption is relaxed, regular group splitting might occur in the form of a fission-fusion dynamics [START_REF] Conradt | Activity synchrony and social cohesion: a fission-fusion model[END_REF][START_REF] Pays | The dynamics of group formation in large mammalian herbivores: an analysis in the european roe deer[END_REF]. If ideal times are not unimodally distributed, more permanent splitting will probably happen, such as in the case of sexual segregation [START_REF] Ruckstuhl | Modelling sexual segregation in ungulates: effects of group size, activity budgets and synchrony[END_REF]. Also, the total constancy of synchronization cost over all group sizes is an idealization since in very small groups with quorum decisionmaking, individual animals will have some influence on the time of common departure and cost will be consequently somewhat lower. Further, our model only addresses possible decrease in size, splitting or dissolution of already existing groups, but not their formation or increase through individuals joining. Finally, the models only address major behavioural transitions in a world where being member of a group is an all-or-nothing state that has to be defended or abandoned during a challenge, rather than a state defined by a continuous proximity in space [START_REF] Ruckstuhl | Modelling sexual segregation in ungulates: effects of group size, activity budgets and synchrony[END_REF][START_REF] Rands | State-dependent foraging rules for social animals in selfish herds[END_REF] or continuous degree of synchronization in time [START_REF] Šárová | Synchronization and leadership in switches between resting and activity in a beef cattle herd a case study[END_REF][START_REF] Michelena | Social cohesion in groups of sheep: effect of activity level, sex composition and group size[END_REF]. Nevertheless, Models I and II apply not only to timing, but also to any situation where a group, in order to stay together, must reach a consensus along a continuous variable such as gradient in space, escalation of conflict with other group etc.

In conclusion, this study documents that synchronization of activity switches such as departures brings a specific type of grouping costs that has specific consequences. First, it may most often cause groups to either stay together or completely dissolve during synchronization challenges. Second, it implies that groups will often react to synchronization challenges with lower-than-optimal synchronization because of the Prisoners-Dilemma type of payoff inherent to the situation. And third, synchronization cost will often set a limit of a smallest group size below which grouping is no longer advantageous over the solitary condition. Let it be continuous function to the exclusion of points -1+b, -1+b+c, -1+b+c+d. Let it be defined by the following Table A1.

Table A1. where strategy J is present can be expressed by

* i t (-1, -1 + b) non-synchronized departures, each at individual * i t 0 * i t (-1 + b, -1 + b + c) synchronized departure at -1 + b + c NR ( d c b b dx x f 1 1 ) ( ) -((-1 + b + c) -x) * i t (-1 + b + c, -1 + b + c + d) synchronized departure at -1 + b + c NR ( d c b b dx x f 1 1 ) ( ) -(1-)(x -(-1+b + c)) * i t (-1 + b + c + d,
dx x f 1 1 ) ( = NR 2 ) 1 ( ) 1 ( c. (A.4) d c b c b dx x f 1 1 ) ( = NR 2 ) 1 ( ) 1 ( d, (A.5) d c b b dx x f 1 1 ) ( = NR 2 ) 1 ( ) 1 ( (c + d), (A.6) c = d ) 1 ( ( A .
B(I, J) = -(-1 + p + q) q p p dx x f 1 1 ) ( + q p p dx x xf 1 1 ) ( + (1-)(-1 + p + q) r q p q p dx x f 1 1 ) ( -(1- ) r q p q p dx x xf 1 1 ) ( + NR ( d c b b dx x f 1 1 ) ( )( r q p p dx x f 1 1 ) ( )
If I is replaced by J, then

B(I, I) = -(-1 + b + c) c b b dx x f 1 1 ) ( + c b b dx x xf 1 1 ) ( + (1-)(-1 + b + c) d c b c b dx x f 1 1 ) ( -(1- ) d c b c b dx x xf 1 1 ) ( + NR ( d c b b dx x f 1 1 ) ( ) 2
We define a function

p + q) q p p dx x f 1 1 ) ( + q p p dx x xf 1 1 ) ( + (1-)(-1 + p + q) r q p q p dx x f 1 1 ) ( -(1-) r q p q p dx x xf 1 1 ) ( + (-1 + b + c) c b b dx x f 1 1 ) ( - c b b dx x xf 1 1 ) ( -(1-)(-1 + b + c) d c b c b dx x f 1 1 ) ( + (1-) d c b c b dx x xf 1 1 ) ( + NR ( d c b b dx x f 1 1 ) ( )( r q p p dx x f 1 1 ) ( ) -NR ( d c b b dx x f 1 1 ) ( ) 2 G(b, c, d) = 0.
Strategy I is stable if and only if G (p, q, r) < 0 for p b, q c, r d and J = I is local maximum of G, G(b, c, d) = 0.

Conditions for the local stability

By the same procedure as in point 2 we will try to find maxima of G inside of simplex p + q + r 2, p 0, q 0, r 0, where in Here strategy I is stable, if, in agreement with procedure in point 2, This equation is identical to (D.2) and therefore e assort * = e*, which means that all animals will depart together during the second challenge and the group size will be stabilized.

2 (1 -) NR ) 0 ( ) 1 ( f ) 1 ( , c = d ) 1 ( ,

Non-Assortment scenario

The density distribution of ideal times after the first challenge and partial synchronization under the non-assortment scenario f nonassort remains the same as before the challenge as follows f nonassort (x) = f(x) for x (-1, 1).

From Appendix B we have for the partial synchronization during the first challenge For any given value of parameter R, the optimal strategy of the reduced one-departure model brings higher net benefit than the optimal strategy for the two-departure model II.

  is the second parameter of the model. The ratio R = Q/L quantifies how important it

  zero, partial and full behavioural synchronization.. After comparison of all possible maximal points for given parameters NR and we obtain the point of global maximum -the optimal strategy I=[t 1 , t D , t 2 ]. The strategy I = [t 1 , t D , t 2 ] is then stable, if the mixed net benefit is locally maximal. We consider an alternative strategy J "sufficiently" close to the strategy I, defined as J = [u 1 , u D , u 2 ]. The mixed net benefit can be defined as B(I, J) = B G (I,J) + B C (J) = NR

  Fig.2shows which of the three basic scenarios (full, partial and no synchronization) brings the highest net benefit for specific combinations of NR and ω. Most of the parameter space is occupied by either no synchronization (the area below the lower full line) or full synchronization (the area above the upper full line) while partial synchronization is superior only within the narrow croissant-shaped area. As ω deviates more and more from 0.5, full synchronization becomes the most prominent strategy. This is because for ω near zero or one, departing after or before the individual ideal time, respectively, gets almost no penalty and thus animals can cheaply achieve a common departure by waiting very late or by acting very prematurely.Fig.2also maps (in grades of gray shading) the net benefits achieved by the optimal strategies. Higher net benefits can be harvested (for the same value of parameter NR) for extreme values of ω, where either premature action or delay are cheap and full synchronization is therefore easily achieved by agreeing to a very early or very late common departure

Fig. 1

 1 FIGURE CAPTIONS

Fig. 2

 2 Fig. 2 Optimal strategies and their net benefits as dependent on the parameters ω and NR. The thick lines demarcate boundaries between the areas of zero synchronization, partial synchronization and full synchronization The net benefit increases with the darkness of the shading. The curved lines are isolines of equal levels of net benefit. The net benefit level is described by the number on the isolines.

Fig. 4

 4 Fig. 4 Differences between benefits of stable and optimal strategies in the parameter space ω, NR. The thick and dotted lines are those given in Figs.2 and 3. The darker the shading, the

Fig. 5 Fig. 6 .

 56 Fig.5Net benefits of the one-departure (solid curve) and two-departure (dashed curve) strategies in Model II, as dependent on NR.

Fig. 7 .

 7 Fig. 7. Group stability as dependent on the group size. For this example, the R parameter was fixed at R = 0.02.

  have the set of N animals. 2. There exist the ideal switch points denominated as * i t for the animal i, 1 i N that are the realizations of the random variable T*, which is defined by its continuous density function f. 3. The function f has the following properties: a. f is symmetric in 0.b. f(t) > 0 in the interval ( , ), f(t) = 0 for t (-1, ).

  c. f'(t) ≥ 0 in the interval ( , ), f'(t) < 0 in the interval ( I = [b, c, d] is defined by costants b, c, d b = t 1 -(-1) = t 1 + 1 c = t D -t 1 d = t 2 -t D (see Fig. 1), Where b + c + d 2, b 0, c 0, d 0. The behaviour of an individual i with the ideal x, b, c, d) be the utility function for ideal switch point placed in position x on time axis.

  Definitions of behaviours and of the utility function. For synchronized departures, the first part of the utility function quantifies the grouping benefit, the second quantifies the synchronization costs. R, are parameters, R 0, 1 0 Intervals defining position of ideal switch time point Behaviour Utility function U(x, b, c, d)

  B for given b, c, d and for given position of ideal switch point (e.g. for given x on time axis) is defined B(x, b, c, d) = f(x)U(x, b, c, d) where f(x) is the density function defined above. The net benefit B for the whole group is defined as B(b, c, d) Conditions for the local maxima of net benefit. To obtain internal maxima we have the following conditions inside of the simplex b + c + d 2, b 0, c 0, d 0 : dB/db = dB/dc = dB/dd = 0. Therefore the following relations must be valid: dB/dd = f(-1 + b + c + d)( As f(-1 + b + c + d) > 0 inside the simplex, we can derive the following relations c b b

  conditions for the function f and according to the mean value theorem we can see that 1As slopes of the lines in the right sides of relations (A.4) -(A.6) are the same, the solution of this system is unique. That means there exists unique triplet [b, c, d] inside of the simplex b + c + d 2, b 0, c 0, d 0 that can be maximum of B. This fact is valid for each [ Other maxima. Other maxima are placed on the border of the simplex. The possible points of maxima are b 0, c = d = 0, and b = 0, c + d = 2. The first case describes the set of non-synchronized switches, the second case describes the set of totally synchronized switches. 4.4. Conclusion -the points of maxima of optimal strategies ZERO SYNCHRONIZATION. If NR ≤ < (1 -) , then partial synchronization strategy [b, c, d] , b + c + d < 2, b > 0, c > 0, d > 0, c = d If NR ≥ (1 -), then the full synchronization strategy [0, c, d] , c + d = 2, such that c = d ) 1 ( , is the optimal strategy and brings the net benefit B(0, c, 2 -c) = NR -Derivation of the benefit functions. Let I = [b, c, d] and J = [p, q, r] are two strategies. Benefit of the strategy I, in the group,

G

  (p, q, r) = B(I, J) -B(I, I) = -(-1 +

  addition b + c + d 2, b 0, c 0, d 0. Conditions similar to (A.1) -(A.3) are dG/dr = f(-1 + p + q + r)of maxima the function G must be concave, that means in our case 2 2 dp G 2 G (b, c, d) = (1-) f(-1 + b + c + d) -f(-1 + b + c ) + f(-1 + b) < 0. That means 0 < b < 1, 2 > b + c + d > 1.

  )d, b < 1 and b + c + d > 1 5.3. Other stable strategies On the border of the simplex p + q + r 2, p 0, q 0, r 0 other stable strategies are placed. Here in addition b + c + d 2, b 0, c 0, d 0. If the point [b, c, d] is point of local maximum of the function G and if we have normal vector n = (n 1 , n 2 , n 3 ) to the border of the simplex p + q + r 2, p 0, q 0, r 0 at that point [b, c, d], then there exist positive numbers 1 1 , 2 2 , 3 3 so that dG/dr -1 1 n 1 = 0 dG/dq -2 2 n 2 = 0 dG/dp -3 3 n 3 = 0 Therefore we can conclude results about all maxima of G. 5.4. Conclusion -stable strategies The following three strategies (and no other strategy) are locally stable: ZERO SYNCHRONIZATION. The zero synchronization strategy [b, 0, 0] , is stable for any combination of values R, and brings the net benefit B(b, 0, 0) = 0. b, c, d] , b + c + d < 2, b > 0, c > 0, d > 0, c = d If NR 2 (1-), then the full synchronization strategy[0, c, d] , c + d = 2, such that c = d ) 1 ( , is stable and brings the net benefit B(0, c, 2 -c) = NR -DETAILS AND SOLUTION OF THE MODEL II 1. ASSUMPTIONS ABOUT IDEAL SWITCH TIMES The assumptions are the same as in Model I in Appendix A. 2. ASSUMPTIONS ABOUT STRATEGIES A strategy I = [r,s] is defined by costants r,s where r 0, s 0, r+s 1. The behaviour of an individual i with the ideal switch time * i t is determined by the position of * i t on the interval (-1, 1), according to Table B.1 3. ASSUMPTIONS ABOUT THE NET BENEFIT FUNCTION Let U(x, r, s) be the utility function for ideal switch point placed in position x on time axis. Let it be continuous function to the exclusion of points -r-s, -s, 0, s, r+s. Let it be defined by the following Table B.1. Table B.1. Definitions of behaviours and of the utility function. For synchronized departures, the first part of the utility function quantifies the grouping benefit, the second quantifies the synchronization costs. R is a parameter, R 0. intervals defining position of ideal switch time point Behaviour benefit -cost function U(x, b, c, d) Conditions for the local maxima of benefit. To obtain internal maxima we have the following conditions inside of the simplex r + s 1, r 0, s 0: dB/dr = dB/ds = 0. Therefore the following relations must be valid: r B r B2 = f(r + s) 4 NR ( placed on the border of the simplex. The possible points of maxima are r = s = 0, and r + s = 1. The first case describes the set of non-synchronized switches, the second case describes the set of totally synchronized switches. 2.4. Conclusion -the points of maxima of optimal strategies so that the following points are valid: ZERO SYNCHRONIZATION. For each NR , NR ≤ B (0, 0) is maximum of B. B( 0, 0) = 0. PARTIAL SYNCHRONIZATION. For each NR , < NR < , there exists [r, s], r + s < 1, r > 0, s > 0 so that relations (1) -(2) are valid and B(r, s) is maximum of B. B (r, s) SYNCHRONIZATION. For each NR, NR ≥ , there exists [r, s,] , r + s = 1B (1 -s, s) is maximum of B. B (1 -s, s) Reduced model I. The model described in Appendix A can be reduced to symmetric case, if = 1/2. We defined the following transformations of variables in the model I: -1 + b = -e, -1 + b + c = 0, -1 + b + c + d = e. Then benefit function is function of one variable maxima of B are similar to general model I: ZERO SYNCHRONIZATION. For each NR , NR ≤ For each NR , NR ≥ 1/2, B(1) is maximal. But that means B1(e) B2(r, s), where e is maximal point for B1 and [r, s] is maximal point for B2. 5.3. Conclusion -models I and II comparison.

B

  challenge, the benefit function is, according to Appendix B, part 5.1 e* is the value of e bringing the optimal benefit during the first challenge.After the first challenge and partial synchronization group size is thus reduced from N to pN where p is the proportion of animals that have departed together on the first challenge and therefore remain in the group p= their optimal times t* i within the interval (-1, -e*) or the interval (e*, 1) disappear from the group. The the new density function f assort (x) is defined as f assort (x) = 0 for x (-1, -e*

  lower proportion of p nonassort < p will depart together on the second challenge, until eventually e nonassort * = p nonassort = 0 and the group will dissolve.

  Figure 1

  Figure 2

Figure 8

 8 Figure 8

ACKNOWLEDGEMENTS

A c c e p t e d m a n u s c r i p t

13 Supported by grant IAA601410612 from GA AV ČR and grant MZE0002701404.