A Pan-European model of the Neolithic

K. Davison1, P. M. Dolukhanov2, G. R. Sarson1, A. Shukurov1 & G. I. Zaitseva3 1 School of Mathematics and Statistics, University of Newcastle upon Tyne, NE1 7RU, U.K. kate.davison@ncl.ac.uk< g.r.sarson@ncl.ac.uk< anvar.shukurov@ncl.ac.uk 2 School of Historical Studies, University of Newcastle upon Tyne, NE1 7RU, U.K. pavel.dolukhanov@ncl.ac.uk 3 Institute for the History of Material Culture, Russian Academy of Sciences, St. Petersburg, Russia


Introduction
The transition to the Neolithic was a crucial period in the development of Eurasian societies, defining to a large extent their subsequent evolution.The introduction of agro-pastoral farming, which originated in the Near East about 12 000 years ago and then spread throughout Europe, is usually considered to be a key feature of this transition (Zvelebil 1996).Yet the Neolithic was not a simple, single-faceted phenomenon.In his early definition of the Neolithic, Sir John Lubbock (1865) specified its main characteristics to be the growing of crops, the taming of animals, the use of polished stone and bone tools, and potterymaking.
Ceramic pottery is one of the defining characteristics of the Neolithic.It is true that there are examples of early farming communities apparently not involved in pottery-making.For example, aceramic Neolithic cultures have been identified in the Levant, Upper Mesopotamia, Anatolia (9800-7500 BC) and also in the Peloponnese (7000-6500 BC) and Thessaly Plain (7300-6300 BC).(All BC dates supplied are radiocarbon dates calibrated using OxCal v3.10 (Bronk Ramsey 2001) with calibration curve intcal04.14c.)Wheat, barley and legumes were cultivated at those sites; permanent houses with stone foundations were used.There is no widespread evidence of pottery (Perlès 2001) but recent excavations have revealed the occurrence of pottery in Thessaly, albeit in small quantities (J.K. Kozłowski, personal communication 27/03/2007).In contrast, the Neolithic in North-Eastern boreal Europe is identified with a sedentary ABSTRACT -We present a mathematical model, based on a compilation of radiocarbon dates, of the transition to the Neolithic, from about 7000 to 4000 BC in Europe.With the arrival of the Neolithic, hunting and food gathering gave way to agriculture and stock breeding in many parts of Europe; pottery-making spread into even broader areas.We use a population dynamics model to suggest the presence of two waves of advance, one from the Near East, and another through Eastern Europe.Thus, we provide a quantitative framework in which a unified interpretation of the Western and Eastern Neolithic can be developed.
In the present paper we attempt to develop a unified framework describing the spread of both the 'agro-pastoral' and 'boreal' Neolithic.Our quantitative model of the Neolithization is based on the large amount of relevant radiocarbon dates now available.

Selection of radiocarbon dates
The compilation of dates used in this study to model the spread of the Neolithic in Europe is available upon request from the authors; unlike all other similar studies known to us it includes dates from the East of Europe.We used data from Gkiasta et al. (2003), Shennan and Steele (2000), Thissen et al (2006) for Southern, Central and Western Europe (SCWE) and Dolukhanov et al. (2005), Timofeev et al. (2004) for Eastern Europe (EE).Our selection and treatment of the dates, described in this section, is motivated by our attempt to understand the spread of agriculture and pottery making throughout Europe.
Many archaeological sites considered have long series of radiocarbon dates: often with 3-10 dates, and occasionally with 30-50.Associated with each radiocarbon measurement is a laboratory error, which after calibration was converted into a calibration error σ i .The laboratory error characterises the accuracy of the measurement of the sample radioactivity rather than the true age of the archaeological site (Dolukhanov et al. 2005) and, thus, is often unrepresentatively small, suggesting an accuracy of 30 years on occasion.Therefore, we estimated an empirical minimum error of radiocarbon age determination of the archaeological age and then used it when treating sites with multiple dates.A global minimum error of σ min = 160 years is obtained from well explored, archaeologically homogeneous sites with a large number of tightly clustered dates.Such sites are: (1) Ilipinar, 65 dates, with the standard deviation σ = 168 years (and mean date 6870 BC); (2) Achilleion, 41 dates, σ = 169 years (mean 8682 BC); (3) Asikli Höyük, 47 dates, σ = 156 years (mean 7206 BC).Similar estimates are σ min = 100 years for LBK sites and σ min = 130 years for the Serteya site in North-Western Russia (Dolukhanov et al. 2005); the typical errors vary between different regions and periods but we apply σ min = 160 years to all the data here.
For sites with multiple radiocarbon date determinations, the dates are treated and reduced to two (and rarely more) dates that are representative of the arrival of multiple Neolithic episodes to that location.For the vast majority of such sites, the radiocarbon dates available can be combined, as discussed below, to just two possible arrival dates.Examples of sites with multiple radiocarbon measurements are Ilipinar and Ivanovskoye-2 where, respectively, 65 and 21 dates have been published.Figures 1a and b indicate that for these sites the series of dates form very different distributions; different strategies are used to process these different types of date series as described below (see Dolukhanov et al. 2005 for details).If a geographical location hosts only one radiocarbon measurement associated with the early Neolithic, then this is taken to be the most likely date for the arrival of the Neolithic.The uncertainty of this radiocarbon date is taken to be the maximum of the global minimum error discussed above and the calibrated date range obtained at the 99.7 % confidence level and then divided by six (to obtain an analogue of the 1σ error).There are numerous such sites in our collection, including Casabianca, Dachstein and Inchtuthil.
If only a few (less than 8) date measurements are available for a site and those dates all agree within the calibration error, we use their mean value and characterise its uncertainty with an error equal to the maximum of each of the calibrated measurement errors σ i , the standard deviation of the dates involved σ (t i ), 1 ≤ i ≤ n, and the global minimum error introduced above: n is the total number of dates in the cluster.An example of such a site is Bademagaci, where we have 4 dates, all within 60 years of one another; Figure 1c shows the histogram of radiocarbon dates of this site.The typical calibration error of these dates is approximately 30 years, thus Eq. (1) yields σ min as an uncertainty estimate.However, we apply a slightly different procedure for clusters of dates that do not agree within the calibration error.
For a series of dates that cluster in time but do not agree within the calibration error, we use different approaches depending on the number of dates available and their errors.Should the cluster contain less than 8 dates, we take the mean of the dates (as in the previous case), as any more sophisticated statistical technique would be inappropriate for such a small sample; the error is taken as in Eq. ( 1).An example of such a site is Okranza Bolnica -Stara Zagora with 7 measurements, and Figure 1f shows that the dates are tightly clustered around the mean value.
If however, the date cluster is large (i.e. more than 8 dates, such as Ilipinar, shown in Fig. 1a), the χ 2 statistical test can be used to calculate the most likely date T of a coeval subsample as described in detail by Dolukhanov et al. (2005): where σ i = max(σ i , σ min ).The coeval subsample is obtained by calculating the statistic: and comparing it with χ 2 .If X 2 ≤ χ 2 n-1 , the sample is coeval and the date T is the best representative of the sample.If X 2 > χ 2 n-1 , the sample is not necessarily coeval, and the dates that provide the largest contribution to X are discarded one by one until the criterion for a coeval sample is satisfied.This process is very similar to that implemented in the R_ Combine function of OxCal (Bronk Ramsey 2001).However, OxCal's procedure first combines the uncalibrated dates into one single radiocarbon measurement and then calibrates it.Our approach on the other hand first uses the calibration scheme of OxCal and then combines the resulting calibrated dates to give T. Furthermore, our procedure adds the flexibility of identifying and discarding dates with the largest relative deviation from T. Within R_Combine the minimum error is not used in the calculation of X 2 but is rather only incorporated into the final uncertainty estimate.We feel that it is more appropriate to include the minimum uncertainty into the calculation from the outset.As a check, we combined several set of dates using both OxCal and our procedure, and the results agree within a few years in most cases where such agreement could be expected.and 1e), here the oldest dates are 6950 and 8800 years BC and the associated errors are 217 and 167 years.
Apart from sites with either no significant peak or only one peak, there are sites whose radiocarbon dates have a multimodal structure which may indicate multiple waves of settlement passing through this location.Ivanovskoye-2 (with 21 dates) is a typical site in this category, and Figure 1b depicts two distinct peaks.In such cases multiple dates were attributed to the site, with the above methods applied to each peak independently.Admittedly our method of assigning an individual date to a specific peak could be inaccurate in some cases as appropriate stratigraphic and/or typological data are not invoked in our procedure.In future refinements to this technique we may consider fitting bimodal normal distributions to the data to avoid the rigid assignment of measurements to one peak or another.After selection and processing, the total number of dates in our compilation is 477.In our final selection, 30 sites have two arrival dates allocated and 4 sites have three arrival dates allocated, namely Berezovaya, Osipovka, Rakushechnyi Yar and Yerpin Pudas.

Modelling
The mechanisms of the spread of the Neolithic in Europe remain controversial.Gordon Childe (1925) advocated direct migration of the farming population; this idea was developed in the form of the demic expansion (wave of advance) model (Ammerman and Cavalli-Sforza 1973).The Neolithization was viewed as the spread of colonist farmers who overwhelmed the indigenous hunter-gatherers or converted them to the cultivation of domesticated cereals and the rearing of animal stock (Price 2000).An alternative approach views the Neolithization as an adoption of agriculture (or other attributes) by indigenous hunter-gatherers through the diffusion of cultural novelties by means of intermarriages, assimilation and borrowing (Tilley 1994;Thomas 1996;Whittle 1996).Recent genetic evidence seems to favour cultural transmission (Haak et al. 2005).
Irrespective of the particular mechanism of the spread of the Neolithic (or of its various signatu- res), the underlying process can be considered as some sort of 'random walk', of either humans or ideas and technologies.Therefore, mathematical modelling of the spread (at suitably large scales in space and time) can arguably be based on a 'universal' equation (known as reaction-diffusion equation) with parameters chosen appropriately (Cavalli-Sforza and Feldman 1981).A salient feature of this equation is the development of a propagation front (where the population density, or any other relevant variable, is equal to a given constant value) which advances at a constant speed (Murray 1993) (in the approximation of a homogeneous, one-dimensional habitat).This mode of spread of incipient agriculture has been confirmed by radiocarbon dates (Ammerman and Biagi 2003;Ammerman and Cavalli-Sforza 1971;1973;1984;Gkiasta et al. 2003;Pinhasi et al. 2005).
In Figure 2a we plot the distance from a putative source in the Near East versus the 14 C dates for early Neolithic sites in SCWE; the linear interdependence is consistent with a constant propagation speed.Due to the inhomogeneous nature of the landscape we would not expect to see a very tight correlation between distance from source and time of first arrival, since there are many geographical features that naturally cause barriers to travel (e.g. the Mediterranean Sea).It is also suggested in a previous work (Davison et al. 2006) that there are local variations in the propagation speed near major waterways; this again detracts from the constant rate of spread.In spite of this, the correlation coefficient is found to be -0.80;reassuringly high given the above complications.There is also a tail of older dates that originate in early Neolithic sites in the Near East, where a Neolithic tradition began and remained until it saturated the area and subsequently expanded across the landscape.
In contrast to earlier models, we include the 'boreal', East-European (EE) Neolithic sites, which we present in the same format in Figure 2b.It is clear that the Eastern data are not all consistent with the idea of spread from a single source in the Near East.A correlation coefficient of -0.52 between the EE dates and distance to the Near East is sufficient evidence for that.Our modeling, discussed below, indicates that another wave of advance swept westward through Eastern Europe about 1500 years earlier than the conventional Near-Eastern one; we speculate that it may even have spread further to produce early ceramic sites in Western Europe (e.g. the La Hoguette and Roucadour groups).
Our population dynamics model, described in detail by (Davison et al. 2006), was refined for our present simulations.We thus solve the reaction-diffusion equation supplemented with an advection of speed V, arising from this anisotropic component of the random walk of individuals that underlies the large-scale diffusion (Davison et al. 2006;Murray 1993) where N is the population density, γ is the intrinsic growth rate of the population, K is the carrying capacity, and ν is the diffusivity (mobility) of the population.We solve Eq. ( 2) numerically in two dimensions on a spherical surface with grid spacing of 1/12 degree (2-10 km, depending on latitude).All the variables in Eq. ( 2) can be functions of position and time, as described by Davison et al. (2006).
We consider two non-interacting populations, each modelled with Eq. ( 2), but with different values of the parameters V, γ, K and ν; the difference is intended to represent differences between subsistence strategies (farmers versus hunter-gatherers) and/or between demic and cultural diffusion.We thus numerically solve two versions of Equation ( 2), one for each of two non-interacting populations with different origins of dispersal.The boundaries of the computational domain are at 75°N and 25°N, and 60°E and 15°W as shown in Figure 3, they are chosen to comfortably incorporate our pan-European area.The environmental factors included into the model are the altitude, latitude, coastlines and the Danube-Rhine river system.The equation describing the farming population also includes advection velocity V along the major waterways (the Danube, the Rhine and the sea coastlines; V ≠ 0 within corridors 10 km wide on each side of a river or 10 km inshore near the sea) which results from anisotropic diffusion in those areas.The prescription of the components of the advective velocity are given in Davison et al. (2006).
The focus of our model is the speed of the front propagation U, since this quantity can be most readily linked to the radiocarbon age used to date the 'first arrival' of the wave of advance.This feature of the solution depends only on the linear terms in Equation (2) and, in particular, is independent of the carrying capacity K.Moreover, to a first approximation U only depends on the product γν: Taking the intrinsic growth rate of a farming population as γ = 0.02 year -1 (Birdsell 1957), the mean speed of the front propagation of U ≈ 1 km/year for the population of farmers suggests the background (low-latitude) value of the diffusivity ν = 12.5 km 2 / year (Ammerman and Cavalli-Sforza 1971;Davison et al. 2006).For the wave spreading from Eastern Europe, U ≈ 1.6 km/year is acceptable as a rough estimate obtained from the EE radiocarbon dates (Dolukhanov et al. 2005); this estimate is confirmed by our model (see Fig. 2d).Analysis of the spread of Paleolithic hunter-gatherers yields U ≈ 0.8 km/year; the corresponding demographic parameters are suggested to be γ = 0.02-0.03year -1 and ν = 50-140 km 2 /year (Fort et al. 2004).These authors use an expression for U different from Eq. (3); it is plausible, therefore, that the intrinsic growth rate obtained by Fort et al. (2004) for hunter-gatherers is a significant overestimate; for ν = 100 km 2 /year and U ≈ 1.6 km/year, the nominal value of γ obtained from Eq. ( 3) is about 0.006 year -1 .A growth rate of γ = 0.01 year -1 has been suggested for indigenous North-American populations in historical times (Young and Bettinger 1992).The range γ = 0.003-0.03year -1 is considered in a model of Paleoindian dispersal (Ste-ele et al. 1998).Our simulations adopt γ = 0.007 year -1 and ν = 91.4km 2 /year for the hunter-gatherers.
For the wave that spreads from the Near East carrying farming, K and ν smoothly tend to zero within 100 m of the altitude 1 km, above which land farming becomes impractical.For the wave spreading from the East, K and ν are similarly truncated at altitudes around 1500 km as foraging is possible up to higher altitude than farming.The low-altitude (background) values of K adopted are 0.07 persons/km 2 for hunter-gatherers (Dolukhanov, 1979;Steele et al. 1998) and 3.5 persons/km 2 for farmers, a value 50 times larger than that for hunter-gatherers (Ammerman and Cavalli-Sforza 1984).The values of K do not affect any results reported in this paper.
In seas, for both farmers and hunter-gatherers, both the intrinsic growth rate and the carrying capacity vanish as seas are incapable of supporting a human population.The diffusivity for both farmers and hunter gatherers tails off exponentially as ν ∝ exp(-d/l), with d the shortest distance from the coast and l = 40 km, allowing the population to travel within a short distance offshore but not to have a sustained existence there.The value of l has been fine-tuned in this work in order to reproduce the delay, indicated by radiocarbon dates, in the spread of the Neolithic from the continent to Britain and Scandinavia.This provides an interesting inference regarding the sea-faring capabilities of the times, suggesting confident travel within about 40 km off the coast.
The inclusion of advection along the Danube-Rhine corridor and the sea coastlines is required to reproduce the spread of the Linear Pottery and Impressed Ware cultures obtained from the radiocarbon and archaeological evidence (see Davison et al. 2006 for details).The speed of spread of farming in the Danube-Rhine corridor was as high as 4 km/yr (Ammerman and Cavalli-Sforza 1971) and that in the Mediterranean coastal areas was perhaps as high as 20 km/yr (Zilhão 2001); we set our advective velocity in these regions accordingly.However, there are no indications that similar acceleration could occur for the hunter-gatherers spreading from the East.Thus, we adopt V = 0 for this population.
The starting positions and times for the two waves of advancei.e., the initial conditions -were selected as follows.For the population of farmers, we position the origin and adjust the starting time so as U = 2 γν.
(3) to minimize the root mean square difference between the SCWE 14 C dates and the arrival time of the modelled population at the corresponding locations; the procedure is repeated for all positions between 30°N, 30°E and 40°N, 40°E with a 1° step.This places the centre at 35°N, 39°E, with the propagation starting at 6700 BC.For the source in the East of Europe, we have tentatively selected a region centered at 53°N, 56°E in the Ural mountains (to the east of the Neolithic sites used here), so that the propagation front reaches the sites in a well developed form.We do not suggest that pottery-making independently originated in this region.More reasonably, this technology spread, through the bottleneck between the Ural Mountains and the Caspian Sea, from a location further to the east.The starting time for this wave of advance was fixed by trial and error at 8200 BC at the above location; this reasonably fits most of the dates in Eastern Europe attributable to this centre.For both populations, the initial distribution of N is a truncated Gaussian of a radius 300 km.

Comparison of the model with radiocarbon dates
The quality of the model was assessed by considering the time lag ∆T = T-T m between the modelled arrival time(s) of the wave(s) of advance to a site, T m , and the actual 14 C date(s) of this site, T, obtained as described in Sect. 2. The sites were attributed to that centre (Near East or Urals) which provided the smallest magnitude of ∆T.This procedure admittedly favours the model, and the attributions have to be carefully compared with the archaeological and typological characteristics of each site.Such evidence is incomplete or insufficient in a great number of cases; we leave the laborious task of incorporating independent evidence in a systematic and de-

Fig. 3. Time lags, ∆T = T-Tm , between the actual and modelled arrival times for the early Neolithic sites shown against their geographical position: panels (a)-(c) refer to a model with a single source in the Near East, and panels (d)-(f) to our best model with two sources (with the second on the Eastern edge of Europe). The positions of the sources are shown in grey in panels (c) and (f). Sites with |∆T| < 500 are shown in (a) and (d), those with 500 yr < |∆T| < 1000 yr in panels (b) and (e), and those with |∆T| > 1000 yr in panels (c) and (f). There are 265, 132, 80 sites in panels (a)-(c) and 336, 113, 28 sites in (d)-(f), respectively. Many data points corresponding to nearby sites overlap, diminishing the apparent difference between the two models. The advantage of the two-source model is nevertheless clear and significant.
tailed manner for future work.Our formulaic method of attribution has inevitably failed in some cases, but our preliminary checks have confirmed that the results are still broadly consistent with the evidence available, (see below).
First, we considered a model with a single source in the Near East (see Fig. 4a for histogram of time lags).The resulting time lags are presented in Figure 3a-c.The best fit model with two sources is similarly illustrated in Figure 3d-f.The locations of the two sources are shown with grey ellipses in panels (c) and (f).
In Figure 3a the sites shown are those at which the model arrival date and the radiocarbon date agree within 500 years (55 % of the pan-European dates); Figure 3d gives a similar figure for the two source model (now 70 % of the pan-European dates fit within 500 years).The points in the EE area are significantly more abundant in Figure 3d than in Figure 3a, while the difference in the SCWE area is less striking.The SCWE sites are better fitted with the one source model, with |∆T| < 500 years for 68 % of data points, but the fit is unacceptably poor for EE, where only 38 % of the radiocarbon dates can be fitted within 500 years.A convenient measure of the quality of the fit is the standard deviation of the time lags with The standard deviation of the pan-European time lags here is s = 800 years.Outliers are numerous when all of the European sites are included (illustrated by the abundance of points in Figure 3c), and they make the distribution skewed, and offset from ∆T = 0 (see Fig. 4a).The outliers are mainly located in the east: for the SCWE sites, the distribution is more tightly clustered (s = 540 years), has negligible mean value, and is quite symmetric.In contrast, the time lags for sites in Eastern Europe (EE), with respect to the centre in the Near East, have a rather flat distribution (s = 1040 years), which is strongly skewed and has a significant mean value (310 years).

Fig. 4. Time lags, ∆T = T-Tm, between the actual and modelled arrival times for the early Neolithic sites. (a)-(c): Histograms of the time lags, with a normal distribution fit (solid), for a model with a single source in the Near East (a), for a single source in the Urals (b) and for a two-source model (c). (d)-(f): The cumulative probability distribution of the time lags from panels (a)-(c), respectively, rescaled such that a normal probability distribution corresponds to a straight line (known as a normal probability plot). The straight lines show the best-fitting normal distribution, and the 95 % confidence interval. A significant reduction in the number of outliers can be seen in (f) or (c) as compared to (d) or (a) and (e) or (b). The distributions of panels (d) and (e) fail the Anderson-Darling normality test, while (f) passes the test confirming that ∆T is normally distributed (p-value = 0.149).
The failure of the single-source model to accommodate the 14 C dates from Eastern Europe justifies our use of a more complicated model that has two sources of propagation.Attempts were made at locating the single source in various other locations, such as the Urals, but this did not improve the agreement (see Fig. 4b for the histogram of time lags for the model with single source in the Urals).
Adding another source in the East makes the model much more successful: the values of the time lag, shown in Fig. 3d-f, are systematically smaller; i.e. there are significantly fewer points in Fig. 3f (5 %) compared to Fig. 3c (17 %).The resulting ∆T distribution for all the sites is quite narrow (s = 520 years) and almost perfectly symmetric, with a negligible mean value (40 years), see Fig. 4c.The distributions remain similarly acceptable when calculated separately for each source (with s = 490 and 570 years for the sites attributable to the Near East and Urals, respectively).The improvement is especially striking in EE, where the sites are split almost equally between the two sources.
We tentatively consider a model acceptable if the standard deviation, s, of the time lag ∆T is not larger than 3 standard dating errors σ, i.e., about 500 years, given our estimate of σ close to 160 years over the pan-European domain.This criterion cannot be satisfied with any single-source model, but is satisfied with two sources.While we would never expect a large-scale model of the sort proposed here to accurately describe the complex process of the Neolithization in fine detail (and so the resulting values of ∆T cannot be uniformly small), the degree of improvement in terms of the standard deviation of ∆T clearly favours the two-source model.The reduction in s is statistically significant, and cannot be explained by the increase in the complexity of the model alone.The confidence intervals of the sample standard deviations s for one-source and two-source models do not overlap (740 < σ < 840 and 480 < σ < 550, respectively); the F-test confirms the statistical significance of the reduction at a 99 % level.
It is also instructive to perform some further basic statistical analysis of the time lags ∆T.We use the Anderson-Darling test to assess if the sample of time lags can be approximated by the Gaussian probability distribution (i.e., in particular, have a symmetric distribution with an acceptably small number of outliers).The null hypothesis of the test is that the time lags have a Gaussian distribution with the sample mean and standard deviation, while the alternative hypothesis is that they do not.This test leads us to accept the null hypothesis in the case of the twosource model (p-value = 0.149) while rejecting the null hypothesis for both one source models.Figure 4d-f show the cumulative probability distributions of the time lags for each model studied, rescaled such that a normal probability distribution corresponds to a straight line (known as a normal probability plot).The straight lines show the best-fitting normal distribution together with its 95 % confidence interval.As quantified by the test, the time lags more closely follow the straight line in (f) than in (d) or (e); the number of outliers is reduced very significantly in (f).Table 1 shows those sites that have |∆T| > 1000 years, i.e., where the disagreement between the data and the best-fit, two-source model is the strongest.There are 28 such sites: 14 of these have not undergone any statistical treatment, while the remaining 14 are a result of date combination or selection as described in Section 2. Five of the dates in Table 1 arise from the four sites (Berezovaya, Osipovka, Rakushechnyi Yar and Yerpin Pudas) where we have been unable to isolate less than three representative dates (see Section 2).This may suggest that a reinvestigation of these sites in particular is required and improved stratigraphic and typological data are required for these sites.
As further quantification of the quality of the model, the χ 2 statistic has been calculated for each model: The results are shown in Table 2.
The values of X 2 , given in Table 2, may then be compared to the χ 2 value at the 5 % level with N-1 degrees of freedom (χ 2 N-1 = 527.86).On all occasions the value of X 2 significantly exceeds χ 2 N-1 (at 5 % level) this is not surprising given the simplicity of our model.The χ 2 statistical test would be satisfied if we discard about one third of the sites.It should be highlighted however that there is an approximate threefold increase in the accuracy of the model with two sources with respect to a single-source model.Some increase in the fit would be expected since we have increased the complexity of the model, but an increase of this magnitude surpasses what we believe could be attributed simply to the increase in model complexity.Development and application of further statistically robust techniques for comparison of our model with archaeological evidence is subject to our ongoing study.
(4) It is instructive to represent the data in the same format as in Figure 2a, b, but now with each date attributed to one of the sources, as suggested by our model.This has been done in Figure 2c, d, where the close correlation of Figure 2a is restored for the pan-European data.Now, the dates are consistent with constant rates of spread from one of the two sources.Using straight-line fitting, we obtain the average speed of the front propagation of 1.1 ± 0.1 km/year for the wave originating in the Near East (Fig. 2c), and 1.7 ± 0.3 km/year for the source in the East (Fig. 2d); 2σ values are given as uncertainties here and below.The spread from the Near East slowed down in Eastern Europe to 0.7 ± 0.1 km/year; the dates from the west alone (as in Fig. 2a) gives a higher speed of 1.2 ± 0.1 km/year.The estimates for the data in both western and eastern Europe are compatible with earlier results (Dolukhanov et al. 2005;Gkiasta et al. 2003;Pinhasi et al. 2005).Care must be taken when using such estimates, however, since the spread occurs in a strongly heterogeneous space, and so cannot be fully characterised by a single constant speed.The rate of spread varies on both pan-European scale and on smaller scales, e.g., near major waterways (Davison et al. 2006).
Our allocation of sites to sources suggested and used above requires careful verification using independent evidence.Here we briefly discuss a few sites.Taking Ivanovskoye-2 (56.85°N, 39.03°E) as an example, the data form two peaks (Fig. 1b 3 and 4 should be reassessed both in terms of the statistical processing of multiple measurements and in terms of the agreement with independent archaeological data.

Conclusions
Our model has significant implications for the understanding of the Neolithization of Europe.It substantiates our suggestion that the spread of the Neolithic involved at least two waves propagating from distinct centres, starting at about 8200 BC in Eastern Europe and 6700 BC in the Near East.The earlier wave, spreading from the east via the 'steppe corridor', resulted in the establishment of the 'eastern version' of the Neolithic in Europe.A later wave, originating in the Fertile Crescent of the Near East, is the better-studied process that brought farming to Europe.
It is conceivable that the westernmost extension of the earlier (eastern) wave of advance produced the pre-agricultural ceramic sites of La Hoguette type in north-eastern France and western Germany, and Roucadour-type (also known as Epicardial) sites in western Mediterranean and Atlantic France (Berg and Hauzer 2001;Jeunesse 1987).The available dates for the earlier Roucadour sites (7500-6500 BC) (Roussault-Laroque 1990) are not inconsistent with The nature of the eastern source needs to be further explored.The early-pottery sites of the Yelshanian Culture (Mamonov 2000) have been identified in a vast steppe area stretching between the Lower Volga and the Ural Rivers.The oldest dates from that area are about 8000 BC (although the peak of the culture occurred 1000 years later) (Dolukhanov et al. 2005).Even earlier dates have been obtained for pottery bearing sites in Southern Siberia and the Russian Far East (Kuzmin and Orlova 2000;Timofeev et al. 2004).This empirical relation between our virtual eastern source and the earlier pottery-bearing sites further east may indicate some causal relationship.
According to our model, the early Neolithic sites in Eastern Europe belong to both waves in roughly equal numbers (56 % to near-eastern wave and 44 % to eastern wave).Unlike elsewhere in Europe, the wave attributable to the Near East does not seem to have introduced farming in the East.The reason for this is not clear and may involve the local environment where low fertility of soils and prolonged winters are combined with the richness of aquatic and terrestrial wildlife resources (Dolukhanov 1996).
Regardless of the precise nature of the eastern source, the current work suggests the existence of a wave which spread into Europe from the east carrying the tradition of early Neolithic pottery-making.If confirmed by further evidence (in particular, archaeological, typological, and genetic), this suggestion will require serious re-evaluation of the origins of the Neolithic in Europe.
Financial support from the European Community's Sixth Framework Programme under the grant NEST-028192-FEPRE is acknowledged.

Fig. 1 .
Fig. 1.Histograms of calibrated radiocarbon ages from archaeological sites in kyr BC, binned into 200 year intervals representing various temporal distributions.(a) The 65 dates from Ilipinar (40.47°N, 29.30°E) are approximately normally distributed, so the χ 2 criterion can be employed to calculate the age of this site as described by Dolukhanov et al (2005).The resulting Gaussian envelope is shown solid.(b) Ivanovskoye-2 (56.85°N, 39.03°E) has 21 dates showing a multimodal structure where each peak can be treated as above.(c) The 4 dates from Bademagaci (37.40°N, 30.48°E) combine into a single date when their errors are taken into account.(d) The 6 dates from Mersin (36.78°N, 34.60°E) are almost uniformly distributed in time, so the oldest date can be used as representative of the arrival of the Neolithic.(e) The 9 dates from Halula (36.40°N, 38.17°E) are treated as in (d).(f) The 7 dates from Okrazna Bolnica -Stara Zagora (42.43°N, 25.63°E) are not numerous enough to justify the application of the χ 2 test, but they form a tight cluster, so the mean date can be used for this site.

Fig. 2 .
Fig. 2. Radiocarbon dates of early Neolithic sites versus the great-circle distance from the assumed source.Inset maps show the location of the sites plotted, and the straight lines correspond to spread at a constant speed given below.(a) Sites from Southern, Central and Western Europe (SCWE) with respect to a Near Eastern source (Jericho).The linear correlation (cross-correlation coefficient C = -0,80) suggests a mean speed of advance of U = 1.2 ± 0.1 km/year (2σ error).(b) Sites from Eastern Europe (EE) show very poor correlation with respect to the same Near-Eastern source (C = -0,52), so that straight-line fitting is not useful.(c) Sites attributed, using our two-source model, to the Near-Eastern source (note a significant number of EE sites clearly visible in the inset map) show a reasonable correlation (C= -0,77) and a mean speed U = 1.1 ± 0.1 km/year.(d) Sites attributed to the Eastern source (from both EE and SCWE) show a correlation similar to that of Panel (c) (C = -0,76), and a mean speed U = 1.7 ± 0.3 km/year.
We further consider those sites which are geographically in the west (i.e., to the west of a boundary set to join the Baltic Sea to the Black Sea) but are allocated to the source of pottery making in the Ural mountain area.These sites are shown in Table3.There are 40 such sites (i.e., 14 % of sites in the west); they deserve further analysis in order to verify the attribution suggested by the model and, if necessary, to further refine the model to improve the agreement with the archaeological data.There are also 104 sites in the east of the above boundary that are allocated to the source of farming in the Near East (i.e.56 % of data points in the east).These sites are listed in Table4.Where a site is characterised by a combined date obtained as described above, only the final age estimate is given (see entry in the column labelled 'Note' for the selection technique applied).All sites in Tables