Muography as a support technique for non-invasive research and three-dimensional localization of tombs in archaeological sites: a case study from Palazzone Necropolis (Perugia – Italy)

Transmission muography is a non-invasive imaging technique that exploits the penetrating power of atmospheric muons into matter to obtain two-dimensional and three-dimensional density images of the monitored structure. The detectors used are particle trackers. Muography enables the monitoring of large structures and it is also particularly useful in the archaeological field for a mapping of low-density underground anomalies potentially related to unknown or inaccessible tombs or tunnels. The Palazzone necropolis, located south of Perugia (Italy), dating back to Etruscan period, contains about 200 known tombs, some of which, such as the Volumni Hypogeum, can be visited thanks to a touristic route. The eastern area of the necropolis, on the other hand, does not have a touristic path and is partially unknown. The objective of the muographic measurement campaign is to support the re-evaluation of this archaeological area by searching for new anthropic cavities and identifying them three-dimensionally. One of the goals of this study is to obtain a three-dimensional localization of cavities starting from a single muographic measurement by exploiting an image focusing algorithm. For this purpose, an area that contains a known cavity was used as the reference cavity for the test of the three-dimensional reconstruction algorithm.


the muon transmission radiography technique
The muon transmission radiography technique (also called muography in the rest of the text) is an imaging technique that exploits the great penetrating power of secondary cosmic muons to create two-dimensional and three-dimensional density images of the internal parts of large volumes [1].
The technique finds many applications in various fields: geological, archaeological, industrial, civil engineering and nuclear safety and is currently known all over the world.The operating principle is similar to X-ray radiography, but uses muons as a source.Muons are secondary cosmic rays produced by the interaction between primary cosmic rays, made up mostly of protons (87%), with the upper layers of the atmosphere.When a primary cosmic ray penetrates the Earth's atmosphere, it interacts with its molecules composed mainly of nitrogen and oxygen, producing a shower of secondary particles.Protons have an 80% chance of interacting with the atmosphere before reaching 15 km altitude and therefore most of them generate showers [2].The component of the shower with higher probability to reach the ground level is made up of muons, thanks to their relatively long lifetime ( ≃ 2.2 μs).They have characteristics very similar to those of an electron  − [3] and have a mass of   = 105.6MeV/ 2 which corresponds to approximately 200 times the mass of the electrons.Muon flux at sea level in the vertical direction is of the order of 70 m −2 s −1 sr −1 for average energies of 4 GeV [3].This quantity varies with changes in latitude, longitude and altitude.The energy spectrum of muons on the ground can reach energies of the order of TeV.
The integral flux of muons [2] varies with the zenith angle  as ∝ cos  (), where  depends on the energy of the muons (i.e. = 2 per  = 3 GeV).The spectrum shifts to a higher momentum as the zenith angle increases and there are higher energy muons than in the vertical case.The dependence on the azimuth angle  is weaker and low energy muons are the ones more affected by this dependence.In -1 -particular, this dependence is related to the Earth's magnetic field and is called East-West asymmetry as more particles arrive from the west rather than the east [2].
Once the muons arrive at ground they interact with matter mainly by ionization and excitation.Given their large mass, the radiative losses (dominant at high energies) in matter are lower than in other secondary cosmic rays (such as electrons for example).Therefore muons are the most abundant charged component at ground level and the most penetrating component in matter.Highly energetic muons can pass through hundreds meters of rock.It is for this reason that they can be used for the imaging of large structures such as pyramids, volcanoes, geosites and structures of civil interest such as dams or buildings.
Historically, the first muographic application dates back to 1955 and was carried out by E.P. George in Australia to measure the rock overburden of a hydroelectric plant tunnel [4].The second application dates back to 1969 and was carried out on the pyramid of Chephren in Giza (Egypt) by Luis Alvarez (Nobel Prize) to search for a hidden chamber [5].With the muon transmission radiography technique, the transmission of muons through the target is studied by measuring the variation in the number of muons arriving on the detector in presence of the target and without the target.The detectors used are particle trackers developed in the field of high energy physics.Various detector technologies have been developed for muographic applications [1] which are based on the compactness of the detectors, low energy consumption, easy transport and installation and low maintenance requirements.
The main advantages of using muography compared to other prospecting techniques concern both the fact that the muon flux is an unlimited and free resource and that this technique is non-invasive.The main disadvantage concerns the long measurement time linked to the flux of cosmic rays on the ground, i.e. to the relatively low frequency of cosmic rays at ground level and on the dimensions of the structure to be studied.For the study of large structures such as glaciers, volcanoes and geosites, acquisition times can vary between a few months to a few years.

Two-dimensional and three-dimensional muon imaging techniques
The muography technique, like X-ray radiography technique, generates a two-dimensional and three-dimensional image of muon transmission through the target.The objective of muography is to search not only for anomalies within the target, but also to give indications about the average density of the object studied.
To obtain a two-dimensional image of average density of the target, comparisons between measurements and simulations are needed and therefore the procedure follows these steps: the measurement in the target configuration, the measurement in the free-sky configuration and the simulation.The first measurement is carried out inside the measurement site by pointing the detector in the direction (, ) in order to observe the target.The second measurement is carried out in the open sky with the detector oriented with the same angles as the previous measurement.From the comparison between the two measurements we obtain the measured muon transmission  meas (, ) for each direction (, ) that falls within the detector's acceptance: where  meas (, ) represents the number of muons measured in the (, ) direction,  acq the acquisition time and  eff is  eff (, ) = (, ) •  with  the acceptance of the detector in the pointing direction -2 -and  its efficiency.This dimensionless quantity can vary between 0 and 1.The measured transmission is 0 if the target is completely opaque to the passage of muons, while it is 1 if it is completely transparent.The opposite quantity, i.e. muon measured absorption, can also be defined as  meas = 1 −  meas .From the simulation phase we obtain the expected muon transmission  simu (, ) in the case in which the target does not present internal anomalies.Therefore the simulation must contain the known geometric model of the target.The required precision of the target geometry depends on the size of the anomalies to be searched for.For small anomalies it is important to have a better precision of the known geometry than for the search for large anomalies.This geometry can be created through 3D modeling software or through the acquisition of point clouds with photogrammetry techniques for example.The simulated transmission is obtained following the relation (1.1) after a simulation with the target and a free-sky simulation.It is important that in the simulation the detector is placed at the same point as the measurement to avoid shifts between  meas and  simu which can lead to artifacts.The two-dimensional angular density distribution is obtained by comparing the simulated transmission by associating a homogeneous density  to the target with the measured one and searching for which density value we have  simu (, , ) =  meas (, ) for each observation direction.With this procedure it is possible to associate a density value (, ) with the measured muon transmission value  meas (, ).
From the individual angular images (transmission and density distributions) it is not possible to locate any density anomalies in space or have a 3D reconstruction of the internal structure of the target.In order to obtain three-dimensional information, specific algorithms have been developed, some of which have been taken from medical physics.These algorithms can be divided into two classes: • algorithms that are based on multiple muographic measurements performed in the same measurement site at different points.In this case the measurement site must allow installation in different points.It is possible to obtain a tomography of the target in terms of density; • algorithms that are based on a single muographic measurement, as the back-projection technique [6][7][8], applicable in all those cases in which, due to the conformation of the measurement site, it is not possible to install the detector in different points.
In the first case, the most used algorithms are the triangulation algorithm and inversion algorithms such as ART [9], SART [10] and DART [11].These algorithms give excellent results in the medical field where it is possible to surround the target with measurements.In the muographic case they are limited by the number of measurements that can be carried out within a measurement site and by the fact that the target cannot be completely surrounded as the flux of cosmic rays comes from above.Some applications can be found in [12,13].
The second type of algorithm is useful both for saving time, as only one measurement is necessary, and in all those applications where it is not possible to install the detector in multiple points.The backprojection technique [6][7][8], developed by the muography group of Florence, is based on the process of image focusing by which different points of the detector see the target from different positions of view and through a sort of internal triangulation allow to locate any anomalies in three-dimensions.The technique is applicable up to distances from the detector equal to ∼ /  where  is the distance between the tracking planes and   the detector angular resolution.This technique was applied in [14,15].

Muographic measurement: the case study of the Palazzone necropolis archaeological site
The Palazzone necropolis is located south of Perugia (in central Italy) and can currently be visited via a tourist route.It contains about 200 underground tombs (figure 1 left) from the Hellenistic period (3rd-1st century BC), with five of them from the Archaic period (6th century BC) and represents one of the most important Etruscan burial sites in the world.Most of the tombs are located at a depth of a few meters and consist of a single room.The access corridor has been excavated in the areas where the surface has changes in slope.The most significant tomb is the Volumni Hypogeum located about ten meters underground with a structure similar to that of a Roman-Italic house with 10 rooms branching off from the central room whose ceiling imitates a wooden ceiling.The numerous archaeological finds are found in the necropolis warehouse structure and in the Antiquarium used as a museum.The necropolis has also been the subject of geological studies and has been defined as an "Archeo-geosite " [16,17] as the tombs represent a three-dimensional environment in which it is possible to study the stratigraphy of the area.From geological studies, two dominant lithologies emerged: the Volumni Unit and the Palazzone Unit [17], both made up of conglomerates and sands with material density values ranging between 1.8 g/cm 3 and 2.3 g/cm 3 with higher values for the Volumni Unit but always in this range.Muography fits into this context, for the study of the eastern area of the necropolis not open to the public in which the presence of tombs not yet mapped is hypothesized.In this area, in fact, there are structures of archaeological importance such as the still intact Tomb n.11 and the columbarium (in which small urns were placed), but also the presence of archaic tombs not yet georeferenced and partially covered by vegetation and debris.Muography can therefore give an indication of the area in which these tombs may be located and possibly discover others.Furthermore, the information on the density of the terrain obtained with muography can be useful to improve the geological knowledge of the studied area.The first objective of the muographic measurement at the Palazzone necropolis was to test the back-projection algorithm on Tomb n.11 for three-dimensional reconstruction.hill not yet open to the public.It is oriented in NE-SW direction and is located inside a small raised area above the target hill.The entrance to the tomb is protected from atmospheric agents through a -4 -cover visible also from satellite images.The tomb has dimensions of about 11 m in length, 8-10 m in width and 4-5 m in height and has secondary corridors that branch off from the central room where niches are visible on the walls.The stratigraphy of the terrain is clearly visible on the walls.The geometry of the tomb was acquired by carrying out a terrestrial laser-scanner survey in order to compare it with the results obtained with muography.

The MIMA detector
The measurement at the Palazzone necropolis was carried out using the MIMA (Muon Imaging for Mining and Archaeology) muon tracker [18].MIMA appears as an aluminum cube of size (50 × 50 × 50) cm 3 housed in a mechanics that allows its altazimuth orientation.The design of the detector meets the lightness and transportability requirements necessary by the type of application.The total weight is 80 kg and it can be disassembled to allow passage even in the narrow places.
The tracker is made up of 3 module XY.Each plane is made up of plastic scintillator bars with a triangular section.The scintillator planes that form a module have the bars oriented orthogonal to each other.The signal developing in each bar is detected by a SiPM (Silicon PhotoMultipliers) which is placed on the triangular bases of the bars.The triangular shape of the bars allows the application of a centroid algorithm for the reconstruction of the impact point of the particle in the scintillator plane.The electronics allows the acquisition of data for each bar when the trigger condition occurs.The trigger condition is generated in the following cases: there has been a signal on all tracking planes, or on 5 tracking planes or on 4 adjacent tracking planes.Details on the electronics and the track reconstruction algorithm can be found in [12,18].The detector is also equipped with a temperature and pressure monitoring system.
MIMA has a spatial resolution of 1.5 mm and an angular resolution of 0.3 • and thanks to the choice to also select the tracks passing through 4 adjacent planes the angular acceptance of the detector with respect to the pointing direction reaches up to ±65 • .

Installation and measurement
The muographic measurement at the Palazzone necropolis involved the observation of the hill located to the east of the necropolis where the presence of tombs not yet discovered is not excluded.In fact, this area is not currently part of the tourist route and has only been partially explored.However, there is evidence of the Etruscan presence as there are structures of archaeological interest such as Tomb n.11 and the columbarium.
The optimal point for installing the MIMA detector was chosen after a preliminary survey and a software simulation.The selected place is located inside the warehouse corridor adjacent to the hill to be observed (figure 2 left).The warehouse is a partly underground structure built from a quarry for aggregates.From the chosen point of view, which is located about 12 m below the ground level, it is possible to observe the entire hill (figure 2 right).Tomb n.11 is located at a distance of approximately 56 m from the detector and it is visible at elevation angles of about 20 • .Therefore, to maximize the detector's acceptance on the known tomb area, it was decided to orient the detector at this elevation angle.The azimuth angle is 41.2 • Data collection began on May 13, 2022 and lasted approximately two months with a trigger rate of 1.9 Hz.During data collection, the environmental factors remained constant and the efficiency of the tracking planes were found to be constant.
-5 - The free-sky measurement was carried out at the INFN1 section of Florence as it was not possible to find the horizon free from distant reliefs at the site of the target measurement.Also in this case, the measurement, carried out in the following months, lasted approximately two months with a trigger rate of approximately 7 Hz.

Muon imaging at Palazzone necropolis
In figure 3 (left) and in the center the two angular distribution of muon counts obtained in the target and in the free-sky configuration are shown, respectively.By comparing the target measurement with the free-sky measurement, exploiting the eq.(1.1), the muon transmission measured at the Palazzone necropolis is obtained.Figure 3 (right) shows the angular distribution obtained in which the presence of the target hill is observed below elevation angles of 50 • .From this distribution it can be observed that the area in which the known tomb is located has higher measured transmission values and that even distant reliefs are visible.On the left and in the center, the angular distributions of muon counts in the target measurement and in the free-sky measurement, respectively.On the right, the measured muon transmission distribution.

Simulation geometry and simulated transmission
The hill under observation has a very thick layer of vegetation which must be eliminated from the geometry of the simulation as it represents a negligible contribution for the muons.The presence of vegetation in the simulation leads to an incorrect evaluation of the thicknesses of materials crossed by the muons and can lead to artifacts.
To acquire a geometry that could better approximate the actual shape of the hill surface, various point cloud acquisition campaigns were carried out: laser scanner, drone and GPS.These clouds were then compared with each other and with LiDAR [19] and DEM made available by the SILENE project [20] in order to find the optimal geometry.In [21] it was found that the geometry that best approximates the real development of the hill is obtained from a merge of several point clouds (figure 4

left).
The simulation carried out is based on a custom software that contains a muon flux model at the ground level obtained from the ADAMO [22] and DEIS [23] magnetic spectrometers data.From the simulation, associating to the target geometry an average density of 2.0 g/cm 3 , we obtain the simulated transmission angular distribution of figure 4 (center).The target hill is seen at elevation angles below about 50 • as in the measurement.In this distribution the acceptance is larger than the measurement, this is due to the fact that in the simulation the detector was taken as point-like.
Reference GPS points acquired on some areas of the hill were then compared with the geometry chosen for the simulation, finding an average deviation of  = |50| cm [21].Then quantitatively evaluating how much this difference could influence the muographic measurement, it was found that, for observation directions of the detector parallel to the hill surface, the difference in thickness between the simulation and a simulation with a geometry shifted upward by a length  are not negligible, reaching values of the order of 8-11 m in some areas.The result of this difference is show in figure 4 (right).These areas were classified as areas with a greater risk of finding artifacts due to the uncertainty in the description of the hill.

Two-dimensional angular distribution of average density and cavities identification
From the comparison between the measured and simulated transmission angular distribution, varying the target density in the simulation according to the procedure described in section 1.1, the average density angular distribution is found.Figure 5 (left) shows the result with the superposition of the areas at greatest risk of artifacts.It is observed that the density values obtained are compatible with -7 -the known densities of the materials that constitute the target hill (1.8 − 2.3 g/cm 3 ).At low elevations for azimuth between 60 • and 100 • an area with higher average density values than the rest of the distribution is visible, indicating the presence of the Volumni Unit which presents higher densities.
Being interested in the search for cavities, we search in this distribution for all the areas that have a density lower than the typical range of the materials that constitute the target.Therefore, excluding artifacts, it is possible to identify areas numbered from 1 to 6 in figure 5 (right), as possible signal areas associated with the presence of cavities.The signal identified with the number 1 is related to Tomb n.11 which is located in that angular region.
Before locating all the observed low-density angular regions in three dimensions, it was decided to apply the back-projection algorithm to localize Tomb n.11 in three-dimensions and estimate the goodness of the reconstruction starting from a single muographic measurement.

Localization and three-dimensional reconstruction from a single muographic measurement of a known tomb as a test of the technique
In order to reconstruct the known cavity (Tomb n.11) in three-dimensions, the back-projection technique was applied [6][7][8]: ninety back-projection maps were generated starting from the center of the detector up to beyond the observation hill.The figure 6 (left) shows a diagram of the procedure.The content  (, ) of each back-projection map, as explained in more detail in [6], consists of the difference between the number of muon counts observed in the measurement  target meas in the direction (, ) and the expected one  simu (, , ) freesky meas (, ) by choosing a target density  to insert into the simulation.In this case  = 2.0 g/cm 3 was chosen obtained by measuring the density of some soil samples.The values  target meas and  freesky meas are normalized to the acquisition times.In figure 6 (center), a Tomb n.11 back-projection map is shown.The focusing algorithm was applied to vertical section of the signal, as indicated in figure.The distance from the detector at which the angular width of the signal is minimum was sought, finding a point cloud that identifies the face of the cavity closest to the detector (figure 7 left, in black).
The air thickness of the cavity was estimated, as explained in [14], by combining the knowledge of the depth of the target hill and the density values found in the angular directions that identify Tomb n.11.The result is shown in figure 6 (right), observing that the thicknesses found (7-12 m) are compatible with the real air thickness of the cavity.Thanks to this estimate, it was also possible to create the -8 - point cloud representing the other faces of the cavity (figure 7 left, in red).By creating a mesh with the CloudCompare [24] software, we observe, in figure 7 (right), how the volume reconstructed with muography well localized the tomb even if a little lower than the real volume.The result is very good if we consider that the tomb is located approximately 56 m from the installation point of the detector.The reconstructed volume presents an error of 1 m on the points focused with the back-projection technique and of 1.3 m on the points obtained with the estimate of the air thickness and is a little larger than the actual volume of the cavity.This may be due to the phenomenon of multiple scattering which produces a smearing of the image which at the distances where the cavity is located results in an error of approximately 1 m.

Conclusions and future developments
Muography is a non-invasive imaging technique which, exploiting the penetrating power of cosmic muons, generates density images of the structure under study.This technique fits into the archaeological context of the Palazzone necropolis for the localization of cavities attributable to tombs in an area not yet completely inspected.This study aimed to identify the presence of cavities and to reconstruct a known tomb three-dimensionally using an innovative algorithm that uses a single muographic measurement.The volume of the known tomb has been reconstructed with good accuracy.In the future the technique will also be applied to the other identified signals.

Contents 1 Introduction: the muon transmission radiography technique 1 1 . 1 2 2 2 4and future developments 9 1
Two-dimensional and three-dimensional muon imaging techniques Muographic measurement: the case study of the Palazzone necropolis archaeological site 4 -dimensional angular distribution of average density and cavities identification 7 3.3 Localization and three-dimensional reconstruction from a single muographic measurement of a known tomb as a test of the technique 8 Conclusions Introduction:

Figure 1 .
Figure 1.On the left: the Palazzone necropolis with the plans of the known tombs.On the right: the laser-scanner survey of Tomb n.11 with the stratigraphy of the materials visible on the walls.

Figure 2 .
Figure 2. On the left: the MIMA tracker installed inside the necropolis warehouse.On the right: the area of the necropolis observable from the detector point of view.The y axis coincides with the north direction.

Figure 3 .
Figure 3. On the left and in the center, the angular distributions of muon counts in the target measurement and in the free-sky measurement, respectively.On the right, the measured muon transmission distribution.

Figure 4 .
Figure 4. On the left: the geometry for the simulation.Center: the expected simulated transmission from the detector's point of view with a target density of 2.0 g/cm 3 .On the right: the distribution indicating the areas with higher probability to find artifacts due to the imperfect geometry of the simulation.

Figure 5 .
Figure 5. On the left and right, respectively, the two-dimensional angular map of average density with the overlap of the areas at greatest risk of artifacts and with the identified low-density signal areas.

Figure 6 .
Figure 6.On the left: the scheme of the generation of the back-projection maps with steps .In the center is the back-projection map at a distance of 50 m from the detector.On the right the air thickness of Tomb n.11.

Figure 7 .
Figure 7. Three-dimensional reconstruction of Tomb n.11 with muography technique (black and red points) compared with the real development obtained with the laser scanner (green point cloud).On the right, in blue, the mesh reconstructed with muography.