Monte Carlo method applied to the mechanical dating of the Turin Shroud

Alternative dating methods to radiocarbon have been recently developed to study the Turin Shroud (TS), the linen sheet that according to the tradition enveloped the dead body of Jesus Christ; among them, a mechanical one is based on the study of five mechanical parameters in reference to the breaking strength, the Young modulus and the loss factor. This method tests single TS linen fibers using a machine, built in for the purpose, capable to measure the stress-stain parameters during loading cycles. These parameters have been already shown to be dependent on the age of the linen under test, and a preliminary result has been obtained: the TS date is 400 A.D. ±400 years at 95% confidence level. A companion paper has combined this mechanical result with two chemical ones, coming from Raman and FT-IR spectra, and the combined result is 90 AD ±200 years at 95% confidence level, thus confirming the compatibility of the age of the TS with the period in which Jesus Crist lived in Palestine. As the evaluation of the uncertainty propagation is not simple in the case of the mechanical parameters, this paper tries to furnish a more reliable result applying the Monte Carlo method to the mechanical parameters directly measured. The obtained mechanical age of the TS of 279 A.D. ±216 years at 95% confidence level is compatible with the previous results, but with an uncertainty almost halved.


Introduction
The Turin Shroud [1,2,3] (TS) is a handmade 3-1 twill linen cloth, 4.4 m long and 1.1 m wide, on which the complete front and back images of a human body are indelibly impressed.According to the Catholic Christian tradition, the TS is the burial cloth in which Jesus Christ was wrapped before being placed in a tomb in Palestine about 2000 years ago; for this reason, it is the most important Relic of Catholic Christianity.Even if the Science has not demonstrated the contrary, the Catholic Christian Church does not impose any veneration of it.The TS enveloped the corpse of a scourged, thorn-crowned, and crucified man who was stabbed in the side by a lance.There are also many marks caused by blood, fire, water, and folding on the cloth that partially obscure the double, front and back, body image [3].The wounds are of great interest to forensic pathologists.In 1988 the TS has been radiocarbon dated in medieval age [4].The result has been widely criticized for many reasons [5,6,7,8], but no alternative dating has been proposed.For this, alternative dating methods have recently been suggested [9] both of chemical (based on Raman and FT-IT spectrometry) and of mechanical type [10].The scientific interest in the TS has been developed from its first photograph in 1898 where the image appeared as a negative.The luminance levels also revealed a three-dimensional (3D) image of a human body, interpreted as mapping cloth to body distance.The bloodstains originated from human blood left by direct body-cloth contact.
The most depth scientific analysis of the TS, performed in 1978 by the STURP (Shroud of TUrin Research Project) [1,2], yielded no explanations for the body image formation on the TS.As the breaking strength, the Young modulus and the loss factor are mechanical parameters that can be easily correlated with the age of a linen fabric, an ad hoc cyclic-loads machine able to measure these parameters of flax fibers has been built and used to date the TS [11,12].A preliminary result has been obtained: the TS is datable to 400 A.D. ±400 years at 95% confidence level, but the uncertainty has been perhaps overestimated because it is not easy to propagate the uncertainty of the measured mechanical parameters to the age of the TS.A companion paper [13] has combined this mechanical result with two chemical ones coming from Raman and FT-IR spectra [9], and the combined result is 90 A.D. ±200 years at 95% confidence level, thus confirming the compatibility of the age of the TS with the period in which Jesus Crist lived in Palestine.
In this paper, after a brief summary describing the cyclic loading machine and the measured mechanical parameters, the Monte Carlo method [14,15,16] has been applied to better evaluate the uncertainty relative to the TS age (the indirect parameter) starting from the mechanical measurements (the direct parameters) of the TS fibers with their uncertainties [11,12].

Mechanical parameters of flax fibers
The age influences the mechanical properties of the vegetable fibers, because they undergo a spontaneous and irreversible degradation, in particular related to the polymerization degree of cellulose, with a corresponding increase of the amorphous cellulose [17].Therefore the mechanical analysis of the fibers can furnish parameters [18,19,20,21,22] that can be directly related to a specific historical age.
The mechanical parameters most suitable to study the aging of the plant fibers appear therefore to be the tensile strength, the longitudinal modulus of elasticity, or Young modulus, and the loss factor.The first parameter can be obtained from simple tensile tests, while the Young modulus and the loss factor must be extracted from cyclic loading tests.
The ancient flax fibers coming from the TS, for which the study of Refs.[10,11,12] were mainly aimed, have in fact a diameter ranging from 5 to 25 Pm and lengths between 1 and 5 mm.As the breaking strength of the modern fibers is much higher than that of ancient ones, it has been necessary to build a new cyclic-loads machine because of the impossibility to use traditional tensile tools, which are indeed unable to control loading cycles characterized by relatively low forces, of the order of few PN.An example of the stress-strain measurements furnished by the new machine is shown in Fig. 1.
The flax fiber has a quite complex geometry [23] that can be evaluated using a SEM (Scanning Electron Microscopy); nevertheless, this microscopic method is quite invasive and therefore not suitable for a subsequent mechanical analysis.On the other hand, a SEM performed on the prestressed fiber can furnish not reliable data because of the induced permanent deformations.
The presence of lumen (central void) in the flax fiber and the non-cylindrical geometry can both influence the cross section evaluation based on the measured "diameter" (intended as the maximum fiber width measured by means of an optical microscope connected to a vision system) also of 10%.The target uncertainty for the measurements of the mechanical parameters has been therefore chosen as ±10% considering a rectangular distribution, that, according to Ref. [14], must be divided by 3 1/2 to determine the corresponding standard uncertainty at 68% confidence level.
It is to remember here that Ref [14] explains how to evaluate the measurement uncertainty and tends to assign Gaussian distributions to the uncertainties of each parameter under analysis, especially if the number of parameters under evaluation is relatively high.This preference is based on the Central Limit Theorem that assures a near Gaussian distribution as the result of the combination of many distributions of different types.
The same Ref.[14] also reports a summary of the coefficients to be used to transform various kinds of distributions to the Gaussian one; in the present paper the coefficient 3 1/2 was used in the case of rectangular distributions.
Proceeding with the calculation of the measurement uncertainty Ref. [14] imposes to evaluate the standard uncertainty as the estimated standard deviation of the measured parameter and in the case of a Gaussian distribution this uncertainty corresponds to a confidence level of 68%.
Nevertheless the final result must be expressed using an expanded uncertainty, frequently evaluated at 95% confidence level, that in the case of Gaussian distributions corresponds to the standard uncertainty multiplied by a factor equal to 2. This paper proceeds just in this way.
A probabilistic distribution of the Gaussian type has been assigned for all the uncertainties that are expressed as standard uncertainty evaluated at 68% confidence level.
The following mechanical parameters have been measured for each fiber under test: -1.Breaking strength, σ R -2.Final increasing modulus, E f , regarding the last increasing part of the highest complete loading cycle among those performed.
-3.Initial decreasing modulus E i regarding the first decreasing part of the highest loading cycle among those performed.
-4. Loss factor η d regarding the last complete loading cycle; it represents the dissipated energy defined as: where D and U are respectively the dissipated and stored energies during a whole loading cycle.-5.Loss factor η i , regarding an inverse loading cycle; this was evaluated in reference to the last unloading phase of the cycle coupled with the last loading action that produces the fiber breaking.
Moduli E f and E i are determined using criteria similar to those adopted by the ASTM standards for the determination of the "initial modulus" [21,22].
As regards the determination of the mechanical parameters, standardization was made, because in the case of the flax the fiber's strength depends on diameter, length and environmental humidity.The "standardization" [10,11] concerning the measured stresses σ was done with reference to a standard fiber having length L=1.0 mm between the two clamps, diameter d=15 μm, relative humidity H=50%.The standard stress σ R applied to each fiber has been thus evaluated as: where K L , K d and K H are respectively the length, diameter and humidity coefficients defined [10,11] as follows: where L [mm] and d [Pm] are respectively the measured fiber length and diameter, and H [%] is the laboratory humidity.
The relationships between each mechanical parameter and the age of the fabric have been determined by measuring the corresponding parameter for nine reference samples having known age, obtained from historical information and/or radiocarbon analysis [10,11,12].Table 1 shows the mechanical parameters measured for each reference sample and Figure 2 displays the corresponding plots.The date of the fibers is shown in year unit, by considering it as a relative number (excluding zero that does not influence the assigned uncertainty) with a negative value, if the date is B.C., and otherwise with a positive one.
It is to note that three mechanical parameters (σ R , E f , E i ) increase exponentially with date, while the loss factors η d and η i decrease almost linearly with date.
In reference to the exponential increase of the first three parameters, this behaviour can be explained supposing that it is directly correlated to the natural cellulose degradation that is a firstorder reaction kinetics, such as the radiocarbon decay [24].The linear behaviour of the loss factor can be instead explained with a constant increment of the cellulosic chains rupture that produces a linear increase with time of the slippage among those chains that in turn increases the parameter in question, which is directly linked with energy dissipation.
Table 1.Mechanical parameters of the nine flax fibers used as reference samples measured by means of the cyclic loading machine.Sample's labels correspond to those of Refs.[10,11]; standard uncertainty corresponds to a confidence level of 68%.

Sample's
Label Other possible relationships between the mechanical variables have also been tested using a linear model or a polynomial interpolation, but the exponential law for the first three parameters in all cases resulted the best fitting, confirming thus the agreement between experimental results and theoretical model.
The independency of the five mechanical parameters has been studied in a companion paper [13] using the Multiple Linear Regression Method based on the least squares estimate.A zero determinant value of a matrix product evidenced the existence of a multicollinearity condition for two couples of the five mechanical parameters.In particular, both the couples of the moduli E f and E i , and loss factors η d and η i , have been found respectively dependent, in agreement with the theoretical model, even if they were defined along different paths of the loading cycle.Therefore, the five measured mechanical parameters were reduced to three.
Table 2 therefore shows the three mechanical parameters defined for the present study where E m and η m are the arithmetical mean values respectively obtained by moduli E f and E i , and loss factors η d and η i .

The Monte Carlo method to determine the calibration curves
As shown in Section 2, the uncertainties of the sample's age and of the measured values using the cyclic loading machine, are characterized by a probabilistic distribution of Gaussian type, and are expressed, if not differently specified, as standard uncertainty that corresponds to a confidence level of 68%.
As it is not simple to propagate the uncertainties of the measured parameters to those of both the calibration curves and the age of the TS, it is useful to apply the Monte Carlo method to determine the probabilistic distribution of the coefficients of the interpolating lines, and in the following Section 4, the probabilistic distribution of the TS age.This method is in fact adopted in stochastic processes when the data vary randomly in a predetermined range.
According to Ref. [15] the three calibration curves of V R , E m , η m have been determined considering the statistically distribution on N curves generated by a proper code written in MatLab ® language.This code processes the nine mechanical parameters for each of the three calibration curves by adding an uncertainty to the corresponding measured mechanical values; the uncertainty is calculated by multiplying a pseudo-Gaussian random number with the corresponding assigned standard uncertainty.
As the first two parameters V R and E m , show an exponential behaviour, the following function has been considered: where m p is the mechanical parameter, y the age, while c 0 and c 1 are constants.Eq. ( 4) may be then represented in a linear form as: where m=c 1 is the slope, q=ln(c 0 ) is the intercept of an interpolating straight line.
In the case of the loss factors, considering their linear dependence on the age y, the calibration formula may be simply written as: where m is the slope, q the intercept.
It is known [14] that the probabilistic Gaussian distribution of the variables' values in the interpolation process, yields a probabilistic Gaussian distribution too, also for the coefficients m and q of the interpolating lines of Eqs.(5,6).
By applying the Monte Carlo method to the three mechanical parameters with N=100 000, the results plotted in Figs. 3, 4 and 5 have been obtained; the corresponding m and q coefficients with their standard deviations are shown in Table 3.As foreseen, the resulting distributions are pseudo-Gaussian.
In reference to Figs. 3 and 4, the sample K may be seen as an outlier, but the same can't be said in reference of Fig. 5. Instead of applying a mathematical method to decide the possible cancellation of sample K, the same Monte Carlo method has been applied to the series of 8 reference samples without the possible outlier in question.As the resulting curves show no appreciable variation and as the TS date varied of few tens of years, the decision to consider or not the sample K in the calibration appeared of second order with respect to the TS result, and therefore the sample K was not considered as an outlier in the present evaluation.

Dating of the Turin Shroud by using the Monte Carlo method
The calibration curves (straight lines) that relate the three mechanical parameters V R , E m , and η m with the flax textile age, can be used to determine the age of the TS.Table 4 shows the mechanical parameters, with relative uncertainty, measured for the TS.
Table 4. Mechanical parameters of the TS; the standard uncertainty corresponds to a confidence level of 68%.Being the age y to be evaluated as a function of the mechanical parameter m p , eqs. ( 5) and ( 6) are respectively rearranged as: where the coefficients m and q have been previously determined by means of the interpolation processes described in Section 3.
The three ages (corresponding to each mechanical parameter) of the TS have been thus determined by applying the Monte Carlo method to Eqs. (5', 6') using the data of Table 4, with N = 100 000.The results are shown in Figs. 6, 7, 8, and in Table 5.
The final age of the TS has been evaluated, by using again the Monte Carlo method to average the three TS partial dates obtained by the corresponding mechanical parameters.
In particular for N = 100 000 times the average of the three mean TS values, to which the corresponding uncertainty defined by the product of a pseudo-Gaussian random number with the assigned standard uncertainty, has been evaluated; the result is plotted in the pseudo-Gaussian distribution of Fig. 9.The mean age of the TS results of 278 A.D. with an uncertainty of 216 years at 95% confidence level.The resulting age of 278 A.D. ±216 years at 95% confidence level obtained by applying the Monte Carlo method to the measured mechanical parameters, is compatible with a previous result reported in Refs.[10,11,12] as of 400 A.D. ±400 years at 95% confidence level, but the uncertainty has been here reduced.It is also compatible with a companion paper that has combined this mechanical result with two chemical ones coming from Raman and FT-IR spectra and the combined result is 90 AD ±200 years at 95% confidence level.
These results clearly show the incompatibility with the 1988 radiocarbon result of 1325 A.D. ±65 years at 95% confidence level.
This result here obtained is also in agreement with the tradition that considers the TS as the burial cloth of Jesus Christ in 30 A.D. or in 33 A.D..The reduced uncertainty obtained using the Monte Carlo method puts however the supposed age in proximity of the limit of the mean value minus two times the standard deviation, i.e. in proximity of the limit of 95% confidence level.
Thus, if the tradition of the TS is assumed as true, this fact leads to think that the mechanical results might be affected by a non-identified bias of about two centuries.This bias could be perhaps related to the fact that one or two more recent fibers not coming from the TS (not excluded coming from the sewed Holland Cloth) have been mixed with the TS ones thus slightly varying the averages of the corresponding mechanical parameters.

GAUSSIAN DISTRIBUTION VERIFICATION
The results obtained by the Monte Carlo simulation are here considered normally distributed both because the theoretical model states that the convolution of some Gaussian distributions gives again a Gaussian distribution, and because they qualitatively show pseudo-Gaussian histograms.
This qualitative evaluation can be confirmed by testing the goodness of fit of the empirical distribution under evaluation with the theoretical Gaussian one by using the classical test of Chisquare.
It is well known that the distribution of Figure 9 can be considered as Gaussian if, subdividing the resulting distribution in k classes and assuming a risk of error of 5%, its Chi-square χ 2 value is within the range whose limits are χ ଵ ଶ and χ ଶ ଶ of Table 6.As the pseudo-random generator of MatLab ® software used to build the pseudo-Gaussian distribution gives dissimilar results per each run, 7 different results have been compared in Table 6.It is evident that in the case of k=9 classes all the 7 runs give a χ 2 within the limits, but when k=103 classes (and therefore 100 degrees of freedom) the χ 2 value exceeds the upper limit in three cases.Therefore an 8 th case has been chosen for the analysis where the χ 2 results within the limits; this case is shown in Figure 9.

Robustness analysis
To perform a verification of the robustness of the Monte Carlo Method one of the nine test samples has been removed from the calibration procedure, and the data of the mechanical parameters of the sample has been used to evaluate its age.The data of the mechanical parameters relative to the eight remaining test samples were used for the calibration, while both the TS sample and the removed one have been again dated.
It was considered appropriate to remove from calibration a middle aged sample with date close to that supposed for the TS: the choice fell on the sample FII (A.D. aged 65 ±6 years at 68% confidence level).
The results are reported in Table 7, where also the data of the TS has been reported as a function of the eight samples.It is to observe that the new age of the TS differs only eleven years from the value calculated using the nine samples.
In correspondence of a known age of 65 ±6 years at 68% confidence level, the Monte Carlo Method furnishes a compatible date of 147 ±107 years at 68% confidence level, thus confirming the robustness of the method.

CONCLUSION
A cyclic-loads machine has been used for the measurement of micro-mechanical properties of single flax fibers of length in the range from 1 to 3 mm.
Five mechanical parameters relative to ancient fibers that have been successively reduced to three because of some dependence allowed defining a method for textiles dating alternative to the radiocarbon dating, which is costly and destructive.
Analyzing the stress-strain plots of loading cycles, the following mechanical parameters sensitive to ageing have been used: breaking strength, Young modulus and loss factors.A pre-selection of the flax fibers has allowed to considerably reduce the bias due to environmental factors like humidity and temperature.
The mechanical dating method has been applied to single flax fibers coming from the backside of the TS, picked up during the 1978 campaign in correspondence of the buttocks area [11,12] in order to obtain an age of the Relic, alternative to that furnished by the 1988 radiocarbon dating [4] that showed a questionable result [8].
To do this the calibration curves, obtained by measuring nine test samples have been employed in agreement with previous works [10,11,12] that dated the TS as 400 A.D. ±400 years at 95% confidence level.
To better evaluate the uncertainty of these results the Monte Carlo method has been applied in this work demonstrating the compatibility of this result with those obtained using other avenues, but reducing the uncertainty of the resulting age of the TS at 278 A.D. of an almost two factor; in particular of 216 years at 95% confidence level.
After controlling that the Chi-square value of χ 2 = 104 is within the limits of the Gaussian distribution, a verification of the robustness of the method has been done.The Sample FII has been removed from the calibration samples, and its known age of 65 A.D. ±12 years at 95% confidence level, has been compared with the calculated one that results of 147 A.D. ±214 years at the same confidence level, thus confirming the compatibility between the known and calculated ages.
The compatibility of the age of the TS with the period in which Jesus Christ lived in Palestine is therefore again confirmed, but the Monte Carlo method evidenced the possible presence of a systematic effect, perhaps due to the one or two more recent fibers measured within the TS ones, that increased the measured age of about two centuries.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 2015

Figure 1 .
Figure 1.Example of the result of a cyclic loading test on a fiber (on the left), and the same result (on the right) obtained after a proper processing addressed to highlight the sequence of cycles.

Figure 3 .
Figure 3. Calibration line relative to the breaking strength (on the top) with N=100 000 determined by means the Monte Carlo method.Lines relative to m±s m and q±s q have been added in the plot.The pseudo-Gaussian distribution of the m and q coefficients (on the bottom).

Figure 4 .
Figure 4. Calibration line relative to the Young modulus (on the top) with N=100 000 determined by means the Monte Carlo method.The lines relative to m±s m and q±s q have been added in the plot.The pseudo-Gaussian distribution of the m and q coefficients (on the bottom).

Figure 5 .
Figure 5. Calibration line relative to the loss factor (on the top) with N=100 000 determined by means the Monte Carlo method.The lines relative to m±s m and q±s q have been added in the plot.The pseudo-Gaussian distribution of the m and q coefficients (on the bottom).

Figure 6 .
Figure 6.On the left TS measured breaking strength (in red) on the calibration plot; the marker x indicates the mean value of TS.Resulting date of 603 A.D. with an uncertainty of 138 years at 95% confidence level in the pseudo-Gaussian plot resulting by Monte Carlo method on the right.

Figure 7 .
Figure 7. On the left TS measured Young modulus (in red) on the calibration plot; the marker x indicates the mean value of TS.Resulting date of 262 A.D. with an uncertainty of 204 years at 95% confidence level in the pseudo-Gaussian plot resulting by Monte Carlo method on the right.

Figure 8 .
Figure 8.On the left TS measured loss factor (in red) on the calibration plot; the marker x indicates the mean value of TS.Resulting date of 29 B.C. with an uncertainty of 596 years at 95% confidence level in the pseudo-Gaussian plot resulting by Monte Carlo method on the right.

Figure 9 .
Figure 9. Pseudo-Gaussian distribution of the age of the TS resulting from the application of Monte Carlo method to the average of the three ages evaluated in correspondence of the three mechanical parameters V R , E m , K m .The mean age of the TS results of 278 A.D. with an uncertainty of 216 years at 95% confidence level.

Table 2 .
[10,11]cal parameters of the nine flax fibers used as reference samples measured by means of the cyclic loading machine.Sample's labels correspond to those of Refs.[10,11]; standard uncertainty corresponds to a confidence level of 68%.

Table 3 .
Coefficients of the interpolating lines (mechanical parameter versus age) and corresponding standard deviations determined by means of the Monte Carlo method.

Table 5 .
Age of the TS resulting from the application of Monte Carlo method to the measured mechanical parameters.

Table 6 .
Results of the Chi-square test χ 2 applied to 7 different runs of the software using both k=9 and k=103 classes, and result of the distribution of Figure9(Run #8) using k=103 classes.The off limit values (external to χ ଵ ଶ and χ ଶ ଶ ) of χ 2 are underlined in italic.

Table 7 .
Age of the TS resulting from the application of Monte Carlo method to the measured mechanical parameters.