Crystallographic Methods for Non-destructive Characterization of Mineral Inclusions in Diamonds

The mineralogy and chemical compositions of inclusions in diamonds are the primary source of information about the environment in which diamonds grow and help constrain the mechanisms of their growth. However, the vast majority of the information about inclusions has been gathered by extracting them from their diamonds, thus destroying all possibility of obtaining further information about the diamond–inclusion system as a whole with new experimental probes unavailable at the time of extraction. One such specific example is the recent discovery by X-ray tomography and in situ spectroscopy of the hydrous silicic fluid film that appears to be ubiquitous around silicate inclusions in lithospheric diamonds (Nimis et al. 2016); the films escaped detection in a multitude of analyses during more than 70 years of research involving the extraction of many thousands of such inclusions. Inclusions in diamond are under compressive stress as a result of their encapsulation at depth and ascent of the diamond to the Earth’s surface; extraction also destroys this stress and thus prevents the depth of entrapment from being determined from the stress state by elastic geobarometry. The stress release on extraction can also lead to the phase changes and/or conversion of the inclusion to a powder (e.g., Joswig 2011). For inclusions from diamonds suspected as being from super-deep sources, extraction therefore risks the loss of rare or possibly unique samples. Non-destructive characterization of inclusions in diamonds should therefore be made in situ whenever possible.

Extraction of inclusions has often been motivated by a lack of experimental probes to measure the inclusions in situ in their host diamonds. However, developments in technology mean that a wide range of spectroscopies is now available that can be used on inclusions in situ to offer rapid initial phase identification (e.g., Raman), and to determine certain aspects of chemical composition. Fourier transform infra-red spectrosocopy (FTIR) can be used to determine OH contents of minerals, synchrotron Mössbauer and µ-EXAFS to determine details about element oxidation states and structural environments in inclusions. None of these methods, however, provide a definitive identification of the mineral phase in the inclusion, nor the details of its crystallographic structure, nor its stress state. These require diffraction measurements. Single-crystal X-ray diffraction has been used since the work of Mitchell and Giardini (1953) to measure inclusions in diamonds. The presence of the diamond presents some challenges to obtaining the highest quality data and therefore the precision with which scientific conclusions can be derived from the data. But, by using methods developed over several decades to study crystals held in diamond-anvil pressure cells (DACs) and exploiting rapid developments in diffractometer technology and intensity of X-ray sources, diffraction measurements of inclusions can be of similar quality to those of crystals at ambient conditions, if the measurements are performed with care.

In this chapter we describe the methods that we have developed as best practice that, when used, will enable the reader to obtain the highest-quality data from most inclusions by using standard laboratory diffractometers. We focus on details of data collection and analysis methods and point out the pitfalls that can arise from poor set-up of the experiment. We show with examples how the crystallographic data measured from inclusions still buried deep within their diamonds can be used to provide simultaneously three key pieces of information to reconstruct the history of the diamond–inclusion pair. Intensity data collection and structure refinement confirms the identification of the mineral phase in the inclusion and allows the composition of the mineral to be determined, and the environment in which the host diamond grew to be characterized. Determination of the crystallographic orientations of inclusions with respect to their diamond hosts has recently provided new fundamental insights into how diamonds nucleate and grow in the mantle. The developments in methods to measure the stress state in inclusions and the theory of elastic geobarometry allows the depth of formation of diamonds to be inferred on the basis of a physical process that is independent from chemical equilibrium. The same principles and methods also apply to those inclusions that today are too small to be measured with laboratory X-ray sources and must still be measured at synchrotron beamlines. And they can be applied equally to all inclusions in diamond, whether the diamond hosts are believed to be lithospheric or sub-lithospheric (so-called super-deep diamonds). Indeed, the identification of the mineralogy of the inclusions and their stress state is the major determinant of the depth the formation of the diamond host.

The mineralogy of inclusions is the primary source of information about the environment in which diamonds grow within the Earth. Structure determination and refinement from single-crystal X-ray diffraction data is a non-destructive method of determining the mineral phase or phases present in the inclusion and, for most minerals that are solid solutions, the composition of the phase. It is especially valuable for colorless minerals which are difficult to identify by optical methods and phases that have not had their Raman spectra recorded previously (e.g., Joswig et al. 1999). Our own experience is that up to 20% of olivine and orthopyroxene inclusions previously examined by optical methods alone (i.e., without Raman) were misidentified. In addition, X-ray diffraction can also reveal the presence of hidden second phases in the inclusion, including inclusions within the main crystal (Joswig 2011; Howell et al. 2012b).

Further, in mineral structures such as olivines, garnets and pyroxenes, in which different cations are distributed over two or more cation sites, the structure refinement can provide information about the cation distribution and state of order which in turn may provide additional constraints on the conditions under which the host diamond grew. For example, cation order of Mg/Fe is believed to be relatively insensitive to pressure (e.g., Angel and Seifert 1999) and therefore could be used as a geothermometer, as has been done for crustal minerals. The specific cation content of a specific site in a mineral structure, such as the excess Si in the octahedral site in titanites (Knoche et al. 1998) or in garnet solid solutions (e.g., Wijbrans et al. 2016; Tao et al. 2018) can provide estimates of the pressure at the time of the growth of the inclusion. The key to the successful use of structural information from inclusions is a precise and accurate structure refinement, which in turn requires the measurement of high-quality X-ray intensity data from the inclusion embedded in the host diamond. This is not a trivial task.

When X-ray studies were first performed on inclusions within diamonds, data was collected on film (Mitchell and Giardini 1953; Futergendler and Frank-Kamenetsky 1961; Harris et al. 1970). The entire process of sample alignment and data collection was very time consuming, taking weeks to complete. Only a few inclusions could be studied which thereby limited the statistical significance of the results obtained. Further, the intensities of diffracted beams collected in this way were of low precision; the mineral species of an inclusion could be identified, but the composition could not be determined with any precision by structure refinement. The development of point detector diffractometers in the 1960s allowed more precise cell parameters to be measured through the use of peak profiling and diffracted-beam centering (Finger and King 1978), and quantitative measurements of diffracted beam intensities to be made. But point detectors are only able to collect one diffraction maxima at a time, meaning that data collections extended over days. And, because they only detect diffraction when the crystal is precisely oriented for a diffracted beam to enter the detector, finding the first few diffraction maxima from an inclusion crystal in a diamond was time consuming (e.g., Glinnemann et al. 2003; Joswig et al. 2003) because routine search procedures typically find the strong diffraction spots due to the diamond, and not the weaker ones of the inclusion (Fig. 1). For studies of inclusions on point detector instruments recourse was frequently made to using film methods to find the first diffraction spots, the same as was done for studies in diamond-anvil cells (Angel et al. 2000). This shows that an area detector such as film is far more efficient for finding diffraction maxima from a crystal of unknown structure and orientation. Therefore, the development of area detectors based on charge-coupled devices (CCDs) and image plates greatly simplified the process. These combined the intrinsic advantages of all area detectors, including film, with quantitative measurements of diffraction peak intensities. Some studies (e.g., Nestola et al. 2011; Howell et al. 2012a) have exploited the best aspects of both types of diffractometer in one study—using the area detector to find the peaks, the point detector to accurately center the inclusion and to determine its’ precise cell parameters, and the area detector diffractometer a second time to perform the intensity data collection.

Image plates and CCDs have, however, three disadvantages compared to point detectors. First, they cannot be collimated, meaning that the diffraction maxima from the inclusion crystal may be often overlapped by diffraction from the diamond which often extends over large areas of the detector (Fig. 1). Second, the read-out time with these detector technologies is slow, often significantly greater than the exposure times for a single frame measurement; a ‘frame’ being a single exposure while the crystal is rotated by a small angle in the beam, the angular change during the exposure being called the ‘frame width’. This prevents these detectors being used for data collections with small frame widths such as 0.2° to 0.5° that are necessary for peak profile measurements. Third, data collections with detectors of these two types are further slowed by the need to expose the detector and then to close the shutter while the data is read from the detector. The development of Hybrid Photon Counting (HPC) technology for area detectors (Forster et al. 2019) overcomes these last two problems of other area detectors by allowing the data signal to be read continuously from every pixel in the detector concurrently with the acquisition of diffraction data. This removes the need to stop the data collection to allow the data to be read from the detector, and so allows the data to be collected continuously during the rotation of the crystal, the same as was always done for diffractometers equipped with point detectors based on photo-multiplier tubes. HPC detectors also have effectively zero noise and very large dynamic range, making them ideal for collecting very weak diffraction intensities from small inclusions in the presence of strong diffraction spots from the diamond. These characteristics enables data to be collected with smaller frame widths (and thus shorter frame exposure times), the shortest time and smallest frame width only being limited by the frequency with which the detector is read, which is typically faster than 20 Hz.

These developments in technology, and the development of computing power and software algorithms to take advantage of the much-improved quality of data, have only become available within the last decade. In this section we describe the procedures for X-ray data collection and structure refinement that have been developed to exploit these advances and to obtain the highest-quality structure refinements of inclusion crystals trapped in diamonds. Many of these procedures have been modified from the techniques developed to measure high-quality structures from crystals under pressure in diamond-anvil pressure cells, and the reader is therefore encouraged to also study the methodological chapters in the Reviews in Mineralogy and Geochemistry volume 41 (Hazen and Downs 2000). The procedures described below for alignment of an inclusion on a diffractometer should also be followed even if only the orientation and cell parameters are required; failure to do so will often give results with significant systematic errors.

The most critical issue is that each inclusion must be accurately centered on the diffractometer prior to measurement. Optical centering alone is not sufficient because the high refractive index of diamond means that the image in the diffractometer telescope will be defocussed and also displaced if the face of the diamond is tilted even slightly with respect to the camera axis (Miletich et al. 2000). If an inclusion is not centered, the positions of the diffracted beams on the detector will include a contribution from the offset of the crystal. These will be interpreted as a wrong Bragg angle, giving rise to incorrect unit-cell parameters (Fig. 2). And the intensities of the diffracted beams will not be correct, at best leading to incorrect refined structures, at worst preventing structure refinement (e.g., Kunz et al. 2002). We therefore recommend the following procedures that have proved successful in a number of inclusion studies on laboratory sources (Joswig et al. 2003; Nestola et al. 2011; Kueter et al. 2016; Milani et al. 2016) and synchrotrons (Sobolev et al. 2000; Nestola et al. 2016a). A recent innovation at a synchrotron source is to use the much higher beam intensity to perform X-ray microtomography to map the inclusion locations in their host diamonds, and thereby center them in the beam for subsequent diffraction measurements (Wenz et al. 2019). In this section a working knowledge of single-crystal diffractometry and diffractometer terminology has been assumed.

Alignment. The exact orientation of the diamond on the goniometer head will be adjusted during the alignment procedure so the diamond must be mounted in a soft material like wax or plasticine, or on a goniometer head equipped with arcs. Mount the diamond so that a (111) face, or a polished window, through which you can best see the inclusion is vertical. In what follows, we call this the viewing face. If possible, mount the diamond in an orientation such that a ϕ rotation on the diffractometer will not bring a second inclusion into the X-ray beam. Then rotate the diamond to set the viewing face exactly perpendicular to the X-ray beam when the diffractometer is driven to zero.

Now drive the diffractometer to the position at which the viewing face is seen edge-on in the telescope. Adjust the physical orientation of the diamond until the viewing face is exactly vertical. Now drive the diffractometer to observe the inclusion through the viewing face, with this face exactly perpendicular to the telescope axis. Adjust the focus of the telescope to bring the inclusion into focus, and then center the inclusion in the telescope by using the translations on the goniometer head. If it is possible to see the inclusion through an opposite face on the diamond, rotate the diamond by 180°. The inclusion should still be centered, but it will probably be out of focus. If you believe that the inclusion is near the middle of the diamond, adjust both the video microscope focus and the goniometer head translation to bring it to focus. Repeat until the inclusion is in focus at both positions. If the inclusion is not in the middle of the diamond, you have to guess the correct position in depth for now.

This procedure should be sufficient to center the inclusion across the X-ray beam, but the centering along the X-ray beam at this point has relied on the optical focusing method, which in turn assumes that the inclusion lies in the middle of the diamond. That works for diamond-anvil cells in which the two anvils are of equal thickness. It does not work for inclusion crystals that lie away from the center of their host diamond. They will be off the center of rotation of the goniometer and thus as the ϕ or ω axes are rotated, the inclusion will be rotated off the centerline of the X-ray beam. On a conventional diffractometer with a sealed-tube source this is not a problem because the beam width is usually 300 µm or more, so the crystal stays in the beam during rotation and diffracted beam crystal centering (described below) can be used immediately to determine the crystal offsets which can then be corrected by moving the diamond with the translation slides on the goniometer head. Such a procedure works even if the inclusion cannot be seen optically (Joswig et al. 2003).

On a modern micro-focus laboratory diffractometer or a synchrotron beamline the full-width at half-maximum intensity (FWHM) of the beam at the crystal position may be 120 µm or significantly less. An inclusion that is not centered on the goniometer will be moved out of the beam during rotation, and it will either stop diffracting (if it is small compared to the beam diameter) or show reduced intensities if it is of similar size to the beam (e.g., see Angel et al. 2016). This phenomenon can be used to further improve the centering on diffractometers equipped with area detectors. Perform a rapid intensity data collection with a single ϕ-scan from 0 to 360°, with a step size of 3–6° per frame. This only takes a few minutes on a laboratory diffractometer equipped with a micro-source and a fast detector (e.g., HPC or CCD). The diffraction pattern should then be indexed and inspected for the disappearance of diffraction from the sample as it is rotated, which indicates that the inclusion has been moved out of the beam.

With the goniometer head translation slides, move the diamond across the beam by a reasonable amount (e.g., by the known FWHM of the beam) and repeat the fast data collection. Repeat this procedure until indexed diffraction spots from the inclusion are present throughout a 360° ϕ-scan. Final centering is then performed by integrating the intensities from the last of these ϕ-scans and observing the frame scale factors. The frame scale factors are determined by minimizing the internal agreement factor of the intensity data and they represent the relative deficit or surfeit in diffraction intensity on each frame (Angel et al. 2016). The measured intensities depend not only on how much of the inclusion crystal is in the beam, but also on the absorption path lengths in the diamond and the inclusion itself. For small silicate or graphite inclusions common in diamonds, the effects of self-absorption in the inclusion vary only slightly with orientation and can be ignored for this step. When an inclusion crystal is centered in the beam (Fig. 3a) the variation in scale factors as the crystal is rotated is only due to the changing absorption path lengths in the diamond and can be estimated quickly from the maximum and minimum dimensions of the diamond (Table 1). For extremely non-equant shapes like plates, the data collected in positions of high diamond absorption, and thus large-scale factors, should be ignored (Fig. 3b). When a mis-centered inclusion crystal moves out of the center of the beam during rotation on the goniometer, it diffracts less intensity so the frame scale factors are modified in combination with the effects of diamond absorption (compare Fig. 3c with 3a). The diamond should now be moved along the beam (when the diffractometer is at zero with the viewing face perpendicular to the beam) and the data collection repeated until the range in frame scale factors is reduced to that expected for the diamond shape (Table 1, and Fig. 3a,b). Normally, at most two iterations of this process of adjustment and 360° ϕ intensity data scan are required. We note at this point that it is critical to analyze the data with software that can immediately and automatically index the diffraction pattern of the inclusion even in the presence of the strong diffraction from the diamond. We use the program Crysalis Pro^TM (Rigaku-Oxford-Diffraction 2016) for this purpose, which is available free on request from Rigaku Oxford Diffraction at their discretion. It can be used to analyze diffraction images from a wide variety of area detectors on both laboratory diffractometers and synchrotron beamlines.

Simple geometrical considerations indicate that the precision with which the inclusion can be centered by observing the frame scales depends on the FWHM of the beam and the size, shape and orientation of the crystal on the goniometer head. Typically, these procedures are normally sufficient to reduce the inclusion offset from the middle of the diffractometer to less than half the full width at half maximum of the X-ray beam, thus maybe an offset of < 40 µm on a laboratory diffractometer with a micro-focus source. Further precision in unit-cell parameters and the integrated intensities can be attained by eliminating the effects of these residual offsets by using the method of ‘diffracted beam centering’ (Hamilton 1974; King and Finger 1979; Angel et al. 2000). On a conventional point detector diffractometer this consists of measuring the same Bragg reflection from a crystal in eight different orientations on the diffractometer.

From the differences in the peak positions (e.g., Fig. 2), the crystal offsets from the center of the diffractometer can be determined and these can be corrected by adjustment of the crystal prior to intensity data collection. These procedures have been standard on point detector diffractometers for many decades (e.g., Finger and King 1978; Angel and Finger 2011). On area detector diffractometers duplicate measurements of the same reflection appear on different image frames within the data collection. The corrections for small displacements of the crystal are therefore made during the analysis of the data after the full intensity data collection.

Intensity data collection. The general requirements for an intensity data collection on any diffractometer, whether equipped with a point detector or an area detector, are basically the same. Data should be collected, if possible, in orientations that minimize the absorption in the diamond host crystal. For approximately equant diamonds, such as octahedra, no useful optimization of diffractometer positions is possible. But data from inclusions within diamond plates (or polished mineral sections) should be collected in the same ‘fixed-ϕ mode’ that is used for DACs (Angel et al. 2000).

Each peak must be collected in a way that allows its profile shape to be determined, because this yields both more accurate unit-cell parameters, more accurate integrated intensities, and allows peaks affected by diffraction from the diamond to be identified. On a point detector instrument this is achieved by making the steps of the ω-scan through each peak significantly smaller than the intrinsic peak width (which depends on both the inclusion crystal and the divergence of the X-ray beam). Typically, a step size of 10–20% of the FWHM of the peaks is optimal, but this also depends on the methods used in the data reduction software to fit the data. With an area detector the minimum frame width is often restricted by motor speeds or detector read-out times, but it must be set so that each diffraction peak appears on several consecutive frames. For HPC detectors on laboratory micro-sources with significant beam divergence, we find that a frame width of 0.2–0.5° gives excellent results when processed in the Crysalis-Pro^TM (Rigaku-Oxford-Diffraction 2016) software.

The coverage of the data collection should be as high as possible. The resolution (the highest Bragg angle to which data is collected) must be sufficiently high that a change in the resolution of the data does not significantly affect the structural parameters of interest. For critical cases this can be determined from test refinements to test datasets from crystals of the same mineral in air, as we describe below. But the basic guideline is that data to a resolution of ~0.7 Å is sufficient for routine structure identification and determination of bond lengths and angles. For the determination of reliable site occupancies in olivines a resolution of ~0.46 Å is required, corresponding to 2q_max = 100° for Mo Ka radiation. In order to obtain reliable frame scales, we collect data with an average redundancy of 6, with the frame exposure times set to obtain an average I/σ(I) = 25 for all reflections. With a small detector the optimal strategy is usually to collect at both positive and negative Bragg angles, and to set the exposure time at high angles to four times that at low angles. On a commercial diffractometer the software can be used to determine an optimal data collection strategy that meets these requirements. As a guide, on an early-generation Oxford Diffraction SuperNova diffractometer equipped with a Mo micro-source X-ray tube and a Dectris Pilatus-200 K detector a data collection on an olivine inclusion to these requirements typically produces more than 2,500 frames of data in 20 hours that yield in excess of 10,000 diffraction peak intensities (Angel et al. 2016). Improvements in source and detector technologies mean that this data collection time is significantly reduced on more modern diffractometer systems, and the same is true when data are collected at a synchrotron beamline because of the higher X-ray intensities.

Intensity data reduction. In order to obtain the highest-quality intensity data it is important to eliminate from the data both the effects of Bragg diffraction from the diamond, and the diffuse background (Fig. 1). This is achieved for both point detector and area detector data by fitting the peak profiles of every single diffraction peak from the inclusion crystal. The low background of diffuse scattering and the tails of the diffraction peaks of the diamonds is eliminated with a combination of local and global background refinement when the peak profiles are fit. Then the profiles are examined and those with unusual shapes are deleted on the assumption that they are affected by overlap with diffraction from either the diamond or another inclusion or they are affected by some other problem. The remaining profile parameters are then used to determine the ‘normal’ peak profile parameters for the sample crystal, and these are then applied to the fitting of the weaker peaks to obtain more precise intensities for them. These procedures have been established for several decades as the ‘learned-profile’ method of intensity integration for point detector instruments (e.g., Diamond 1969; Angel 2003), and are incorporated into commercial software for area detectors. Experience shows that even with the advanced algorithms now available for eliminating peaks with poor profiles, it is necessary to check the results of the integration carefully. We find that the ‘profilecorrint’ plot provided in Crysalis-Pro^TM (Rigaku-Oxford-Diffraction 2016) is a sensitive indicator of residual problems in the peak profiles. The frame scale factors should also be inspected to see that they lie within the range expected given the shape of the diamond host (Table 1). For olivine inclusions, a good quality dataset has an overall R_int less than 0.04, with R_int < 0.02 for low-angle data and R_int < 0.10 for high-angle data, with the average F²/σ(F²) > 8 in all angle ranges.

On an area-detector diffractometer the orientation of an inclusion can be determined within 1–2° by using peak positions found by a simple peak hunt in the dataset but the unit-cell parameters from such a process can be in error by several %. On the other hand, the integration provides peak positions from the fitting of the peak profiles of only those reflections unaffected by overlaps with diamond reflections. These peak positions should then be used to determine not only the crystal offsets and other diffractometer aberrations, but also the unit-cell parameters and orientation unaffected by the offsets and aberrations.

The data should then be averaged in a way to eliminate reflections with intensities that are inconsistent with their symmetry equivalents. These may be due to poor frame scaling or overlaps in the data frames from diamond reflections or reflections from other inclusions; rejection on the basis of peak profiles does not always catch all weak overlaps. Further, the diamond can also reduce the intensity of the inclusion peaks by diffracting away part of the incident beam, or by diffracting the diffracted beam from the sample (Loveday et al. 1990). These types of events do not change the diffraction profile, but only lower the intensity, and they can only be eliminated by averaging. We recommend detecting and eliminating affected reflections by applying the Blessing (1987) criteria for outlier rejection first to the low-intensity outliers, prior to averaging the remaining data in the normal way. The Average program (Angel 2004) is freely available for this task.

With the procedures described above for the collection and reduction of X-ray intensity data from inclusions, the data quality can be similar to that from a crystal in air and significantly better than for data collected from crystals at high pressures in DACs for which the reciprocal space access restrictions of the DAC mean that the datasets from one crystal are typically incomplete (Miletich et al. 2000). Therefore, it is possible to proceed with structure solution (if required) and routine refinement in exactly the same way as for data collected from crystals in air. Such procedures have yielded the structures of inclusions in lithospheric diamonds of olivine (Nestola et al. 2011; Kueter et al. 2016), garnets (Joswig 2011; Kueter et al. 2016), pyroxenes (Joswig 2011; Nestola et al. 2016a) and a high-pressure phase of CaSiO₃ (Joswig et al. 2003; Brenker et al. 2021). The quality of refined structures is sometimes improved by the use of robust-resistant weighting (Prince and Collins 1992) to identify and eliminate the influence of any remaining outliers in the intensity dataset. Secondary extinction effects in inclusion crystals seem to be surprisingly strong and should be corrected with a physically-realistic model such as that of Becker and Coppens (1974).

The determination of the composition of an inclusion is often the motivation for performing a structure refinement. For simple solid solutions such as (Mg,Fe)O the increase in cell edge and volume with increasing Fe content gives an initial indication of composition. But this is a minimum estimate of Fe content because the unit-cell parameter of the inclusion is also decreased from its room-pressure value because the inclusion is under pressure. The true composition of the inclusion can be determined in situ by structure refinement, in which the Fe/Mg ratio of the only cation site in the structure is refined (Nimis et al. 2018). For more complex structures in which the cations can be distributed over two or more sites, extra care has to be taken to use the appropriate refinement models established for the mineral species, and to verify that the model returns valid results for intensity datasets collected from inclusions.

Olivine. We use data collected from a suite of olivine inclusions in diamonds from the Udachnaya kimberlite to illustrate the further precautions that may be required to determine by structure refinement the composition of a mineral that is a solid solution. These inclusions are a subset of those whose orientations were determined and analyzed by Nestola et al. (2014). The composition of these olivines is required in order to estimate their residual inclusion pressures from their unit-cell volumes measured while still trapped in their diamond hosts, because the room-pressure volume of mantle olivines varies significantly with composition.

The most important point to remember while performing structure refinement is that an artificial model of the scattering power of the atoms in the crystal structure is being refined to the intensities of the Bragg diffraction peaks. The model is artificial because it considers the crystal to be perfect, without defects, and that all unit-cells within the crystal are identical in structure and composition. This is not true in a solid solution, such as olivine or garnet. Further, the independent atom model (IAM) most used for conventional structure refinement considers the atoms as independent scatterers of radiation, so there is no representation of the bonding electrons which also contribute to the scattering of the X-rays. While the IAM will certainly provide the correct structural architecture for the inclusion and quite precise and accurate bond lengths and angles it often results in incorrect site occupancies for cations (e.g., Angel and Nestola 2016). Thus, any refinement model must be first tested by refinement of the structure of a crystal of the same mineral of similar and known composition. We have found that the most reliable structure model is different for different mineral groups, but extensive reading of the crystallographic literature will provide good indications of the most appropriate model to use on most mineral groups, and its limitations.

Chemical analysis of olivine inclusions in lithospheric diamonds has shown their chemistry to be dominated by Mg and Fe, with trace amounts of Ni, Ca, and Mn (e.g., Stachel and Harris 2008). However, the total scattering power of each site in the structure is proportional to the number of electrons associated with the cation or mixture of cations occupying the site. Thus, it is impossible to determine how the several cation species are distributed across the M1 and M2 sites in the olivine structure. Cations therefore have to be put together in groups of similar number of electrons and, hence, scattering power. Trace amounts of Ni and Mn can be modelled by Fe, while the Ca content is effectively modelled by Mg and Fe equally. Structure refinements were then performed with just Mg and Fe allocated to the M1 and M2 sites and the tetrahedral site was assumed to be fully occupied by Si. The effective Mg number from the structure refinement is then simply the average of the refined Mg occupancies of the M1 and M2 sites, and it should correspond to the value of

determined from independent chemical analysis of the olivines.

Refinements of the structures of a crystal of known composition that is similar to those of inclusions and a crystal of pure Mg₂SiO₄, both measured in air, and DFT simulations of pure forsterite were combined to develop a robust refinement model. Conventional structure refinements using scattering factors for spherical atoms (the IAM) yield cation site occupancies for all three cases that are systematically in error (Angel et al. 2016; Angel and Nestola 2016). This can be attributed to the failure of the IAM to account for scattering by the bonding electrons (Kroll et al. 1997) and it also leads to the refined cation site occupancies being strongly correlated with both data resolution and the values of the atomic displacement parameters. The ‘electron in bond’ model (Kirfel and Will 1980; Heinemann et al. 1999) which models the bonding electrons by an additional 1 electron scatterer placed in line between the Si and O positions removes much of these correlations. In combination with robust-resistant weights in the least-squares, the data collection procedures and structure refinements of olivine inclusions yield values of Mg* in agreement with microprobe analyses of the same crystals (Angel and Nestola 2016).

These sorts of refinement tests also allow the optimal data collection parameters including resolution, coverage and redundancy, for a given type of inclusion to be determined prior to measurement. Single-crystal X-ray intensity datasets were then collected from 13 olivine inclusions in 5 diamonds with an Oxford Diffraction SuperNova diffractometer equipped with a Mo micro-source and a Dectris Pilatus-200K area detector, with the procedures that we described above. The number of averaged reflections used in the structure refinement ranged from ~1,000 to ~1,600, the number being strongly correlated with the average ‘signal to noise’ of the data <F²>/<σ(F²)> which reflects differences in inclusion sizes and hence total scattering power. Longer data collections on the same instrument, or data collections with a more modern instrument or more intense synchrotron source would suppress these trends. Quality indices for the refinements ranged from R₁ = 0.017 and wR₂ = 0.026 for a dataset with average <F²>/<σ(F²)> = 35, to R₁ = 0.071 and wR₂ = 0.095 for a dataset with average <F²>/<σ(F²)> = 11. Uncertainties on bond lengths vary in the same way, for individual Si–O bonds from 0.0003 Å to 0.002 Å and those of M–O bonds are all better than 0.001 Å. The refined structures pass the most sensitive quality test for silicate structures; that the root-mean-square displacements of the Si and O atoms towards one another, that are calculated from refined anisotropic displacement parameters (Trueblood et al. 1996), are consistent with the strong bonding of the SiO₄ tetrahedra (Downs et al. 1990; Kunz and Armbruster 1990; Angel et al. 2016).

Figure 4 shows that the refined compositions of these olivine inclusions fall into the middle of the range of all chemically-analyzed olivine inclusions in cratonic diamonds (Stachel and Harris 2008). Further, while there are differences between the compositions of olivines in different diamonds, the compositions of all of the inclusions in any one diamond are the same within the estimated uncertainties (which come directly from the least-squares refinement of the structure). This result is consistent with the conclusion that some of these inclusions, such as all four in diamond 9 and inclusions 13B and 13C, are remnant pieces of the same original olivine crystals in the mantle at the time of diamond growth (Nestola et al. 2014). The geological reasonableness of the results summarized in Figure 4 demonstrates that the difficulties of data collection and structure refinement can be overcome to yield crystal structures of inclusions of a quality comparable to those obtained from samples measured in air.

The idea that the orientation of mineral inclusions within their host diamonds may imply something about whether the inclusion minerals grew prior to, or during, diamond growth dates back at least to the 1950s. The reasoning was that if an inclusion in diamond shows an epitaxial relationship with its diamond host, then it must have grown on a pre-existing free surface of the diamond. In this scenario, in order for the inclusion to be trapped, it is argued that diamond growth must have continued after growth of the inclusion crystal and therefore the diamond and inclusion growth occurred simultaneously. This would then constitute proof that the inclusions were syngenetic with the diamond and therefore geochemical analysis and dating of the inclusions would indicate the conditions at the time of the growth of the diamond and the age of the diamond. This logic is flawed; an epitaxial relationship could also result from a diamond nucleating and growing on a pre-existing mineral leading to a protogenetic inclusion.

A defining feature of ‘epitaxy’ is that the two minerals exhibit a very specific and crystallographically-simple orientation relationship that arises because of structural matching and bonding across the contact plane between them. The periodic nature of crystalline structures means that an orientation relationship that arises from epitaxy is exact. If epitaxy was an active process in the entrapment of a mineral phase in diamond, then it should result in a large number of inclusions of that mineral showing exactly the same crystallographic orientation with respect to their diamond hosts. Early studies reported the orientations of a small number of olivine, garnet and magnesiochromite inclusions (e.g., Mitchell and Giardini 1953; Futergendler and Frank-Kamenetsky 1961; Zyuzin 1967) relative to their host diamonds. Despite not meeting the strict requirements for epitaxy these orientations were then interpreted to reflect an epitaxial relationship (Hartman 1954; Futergendler and Frank-Kamenetsky 1961; Zyuzin 1967). However, further measurements showed that specific orientations between inclusions and diamond are not a common phenomenon (Harris et al. 1967; Harris 1968). And subsequent analysis showed that the claims made by Futergendler and Frank-Kamenetsky (1961) for structural matching between individual planes of the structures of olivine and diamond were crystallographically incorrect and physically impossible (Nestola et al. 2014). Nonetheless, despite the limited number of supporting data and the long-standing errors in crystallographic interpretation, the claims for epitaxy between diamond and its inclusions have long been used as evidence in favor of their syngenesis (e.g., Orlov 1977; Harris and Gurney 1979; Pearson and Shirey 1999). This was in large part due to the failure to clearly distinguish between the measurements of orientations and their interpretation in terms of the genetic mechanisms which may have created them, with the result that the interpretation in terms of epitaxy and syngenesis was promoted at the expense of clear and statistically-meaningful data evaluation.

In recent years, the advent of improved instrumentation has allowed many more inclusions to be measured in situ allowing statistically significant conclusions to be drawn about their crystallographic orientations relative to their host diamonds. We first describe the methods of measurement and data analysis, which yield the raw data on the orientation relationships exhibited by a given mineral when it occurs as an inclusion. Then we use several examples to show how different types of relative orientations can be used to suggest possible mechanisms for the entrapment of the inclusions by their host diamonds. These mechanisms often have implications for the environment in which the diamond grew. We also note that orientation information is also required if one wants to interpret the measured stress of inclusions of non-equant shapes in terms of entrapment conditions, as we discuss in the following section of this chapter.

Orientation measurements. The positions of the diffracted beams from a crystal (e.g., Fig. 1) on the detector of a diffractometer depend not only on the unit-cell parameters of the crystal and its symmetry, but also on the physical orientation of the crystal on the diffractometer; if the crystal is rotated the positions of the ‘spots’ on the detector will also be moved. Therefore, the orientation of an inclusion crystal relative to the diffractometer can be determined by measuring its diffraction pattern and indexing it. This yields the ‘orientation matrix’ (Busing and Levy 1967) which describes the position of the diffraction spots in reciprocal space relative to the diffractometer. The orientation matrix defines both the unit-cell parameters of the crystal and the orientation of the unit-cell axes with respect to the external reference system of the diffractometer. The orientation of the inclusion relative to its host diamond can therefore be determined by measuring the diffraction patterns of both crystals simultaneously. Their respective orientation matrices will each define their orientations relative to the diffractometer, and appropriate mathematical manipulation of these two matrices allows the relative orientation of diamond host and inclusion to be determined (Angel et al. 2015b).

In order to determine the orientations of inclusions by single-crystal diffraction it is important that the inclusion is well-centered on the diffractometer by the methods described in the previous section ‘Single-crystal diffraction experiments’. Failure to accurately center the inclusion will result in incorrect peak positions (Fig. 2) and thus small errors in the orientation measurements and, perhaps more importantly, it will lead to difficulties in indexing the diffraction patterns as the inclusion is rotated out of the beam. Once the inclusion is centered, the host diamond adjacent to the inclusion will also contribute peaks to the diffraction pattern (Fig 1). In order to obtain the orientation data we recommend performing a relatively fast data collection scan with a small frame width which allows the peak positions to be determined by peak-profile fitting, as is done for integrating the intensities for a structure refinement (see previous section). If enough data are collected, then the offset of the inclusion from the center of the diffractometer can be calculated from the data and its effects on the orientation matrix eliminated (Hamilton 1974; King and Finger 1979).

Diamond spots appear together with the diffraction from the inclusion, but the orientations of both crystals can normally be determined unambiguously by indexing. For a small crystal in air, well-centered on the diffractometer, the uncertainty in orientation determined in this way is significantly less than 0.1°. The relatively large size of the diamonds results in large diffraction spots, whose size may result in an uncertainty in orientation of the diamonds alone of up to 1.5°. A practical test of the precision of orientation measurements is simply to remount the diamond in several different orientations on the diffractometer and remeasure the same inclusion each time. Nestola et al. (2014) showed that the uncertainties in the relative orientations of olivine inclusions relative to their diamond hosts is no more than 2°. The orientation of magnesiochromite inclusions, which have smaller unit cells than olivine and therefore fewer diffraction spots, are reproducible to better than 4° (Nimis et al. 2019). If the inclusions are oriented such that their reflections mostly overlap with those of the diamonds (e.g., Nimis et al. 2018), then the determination of the exact peak positions is difficult, but the nature of the orientation relationship will be obvious even if it cannot be measured precisely.

EBSD measurements can be made on partially exposed inclusions to determine their orientation with respect to their diamond hosts, because the polishing of the diamond and exposure (and consequent partial stress release) of the inclusion is not expected to rotate the inclusion by more than 1–2°. A single EBSD measurement measures the orientation of a small volume of the inclusion corresponding to the size of the electron beam. This contrasts with an X-ray diffraction measurement that provides the average orientation of the entire single-crystal inclusion. EBSD measurements should therefore be performed on many points in the inclusion. If the EBSD measurements show that the entire inclusion has the same orientation, then it should be used as a single datum in subsequent analysis. Orientations of crystals measured by EBSD are normally expressed in terms of Euler angles which characterize the rotations required to bring the laboratory or instrument axial reference system into coincidence with the crystallographic axes. Or vice-versa. Care therefore has to be taken to use the appropriate conventions for the definition of the Euler angles when interpreting data from any instrument. Once these conventions are known, the Euler angles can be used to calculate an orientation matrix, equivalent to that determined by single-crystal X-ray diffraction measurements, if that is required (Angel et al. 2015b).

Orientation analysis. Both EBSD and single-crystal diffraction yield the orientation of the diamond and the inclusion relative to an instrument coordinate system. These descriptions of the orientations of crystals relative to the laboratory reference axes are not unique, because of the symmetry of the crystal. For example, olivine has orthorhombic symmetry, which makes the direction of the positive a-axis symmetrically-equivalent to the opposite direction. Because these directions are symmetrically equivalent it is impossible by diffraction or any other physical measurement method to distinguish them. Thus, for example, an olivine described as having its ‘a-axis’ vertically upwards on the instrument could be equally well described as having the same axis vertically downwards (Fig. 5). The orthorhombic symmetry of olivine therefore leads to 4 symmetrically equivalent and indistinguishable ways of indexing its diffraction pattern, and hence 4 symmetrically equivalent and indistinguishable descriptions of the orientation of one olivine crystal on the diffractometer or in an EBSD instrument. Figure 6a shows the directions of the a-, b- and c-axes of one olivine inclusion relative to its host diamond for all of these 4 symmetrically equivalent ways in which the olivine diffraction pattern can be indexed. It can be seen in Figure 6a that each crystallographic axis has two possible orientations, one in the upper hemisphere (closed symbol) closer to the +c direction than the −c direction of the host diamond, and one in the lower hemisphere (open symbol) closer to the −c direction of the diamond. Figure 6a therefore shows that even if all of the olivines in one diamond were measured with the diamond in the same orientation on the diffractometer, their orientations could appear different if this ambiguity caused by symmetry is not taken into account. A similar ambiguity occurs for indexing the diffraction patterns of the diamond host crystals; being cubic, the choice of how to label the three perpendicular crystallographic axes a, b, or c is arbitrary and this alone generates 24 equivalent relative orientations between the olivine inclusion and diamond (Fig. 6b). In combination the symmetry of the two crystals generates the 96 equivalent relative orientations shown in Figure 6c. Failure to account for this indexing ambiguity can result in olivine orientations relative to the diamond appearing more random than they really are. Data from multiple inclusions that resemble Figure 6c appear random, but in fact the inclusions have identical orientations when the symmetry is considered.

OrientXplot (Angel et al. 2015b) is a freely-available (www.rossangel.net) analysis program which addresses this problem by using the crystallographic symmetry of the inclusion and the diamond to generate all possible symmetry-equivalent relative orientations of the pair of crystals that are consistent with the diffraction patterns. All of the stereograms shown in this chapter were created directly with OrientXplot. The symmetry of both host and inclusion are entered into OrientXplot by the user, so that the effects of symmetry can be explored and accounted for. The program reads and correctly interprets the orientation matrices from all major commercial single-crystal diffractometers (Oxford Diffraction, Rigaku, Bruker, Stoe) so that data from different instruments can be analyzed together. User-defined conventions can also be used to analyze data collected on non-standard instruments such as often found on synchrotron beamlines. Software for EBSD measurements normally plots the data on a stereogram that is based on the instrument coordinate system (e.g., see Neuser et al. 2015) and we are not aware of any such software that accounts for the ambiguity due to symmetry. Therefore, OrientXplot accepts EBSD data as either Cartesian rotation matrices or as sets of Euler angles (with various conventions). Full details of input data formats, the mathematical methods used by the program, and worked examples are provided in Angel et al. (2015b) and the program documentation. Note that the common practice of plotting the poles of just the crystallographic axes of the diamond and the inclusion relative to the instrument coordinate system (e.g., Neuser et al. 2015; Davies et al. 2018) does not allow data from different diamond hosts to be quantitatively compared, and will often obscure the existence of anything but very specific orientational relationships. To clearly see the orientation relationships of a suite of inclusions to their hosts, OrientXplot can therefore plot the relative orientations of different inclusions in different diamond hosts simultaneously on the same stereogram.

For olivine inclusions, the ambiguity caused by symmetry means that the 96 symmetrically-equivalent relative orientations (Fig. 6c and Nestola et al. 2014) are all equally-valid descriptions of the orientation relationship between the diamond and the inclusion. There are several ways to choose which individual relative orientation to plot for analysis of the relative orientations of several inclusions of the same mineral. One is to choose the orientation which puts one particular axis of the inclusion within the asymmetric unit of the mutual orientation space of the two crystals. Another is to choose the relative orientation that puts one particular direction in the inclusion as close as possible to a direction in the diamond host. The criteria should be chosen carefully so as to test hypotheses about the possible relative orientations of inclusions and diamonds; the use of different choices is illustrated in the examples described below. Provided that the criteria always define one, no more and no less, of the possible symmetrically-equivalent orientations, the choice of which one to display does not affect the conclusions that are drawn from the resulting stereographic projections. If one choice indicates a strong orientation relationship, then all choices will do so. If there is no preferred orientation, then this will be evident in all plots of all valid choices of orientation pairs.

Crystallographic orientation relationships (CORs). The historical literature on the orientations of inclusions in diamonds has been confused for 70 years by the conflation of orientation measurements with their analysis and interpretation. For example, the word ‘epitaxy’ implies both a specific type of orientation relationship and the mechanism by which it was created. This has generated confusion over what the phrase ‘orientation relationship’ between inclusions and diamond means. Clarity has been brought to the subject by the concepts and terminology of Crystallographic Orientation Relationships (CORs) between a host and its inclusions (Griffiths et al. 2016, 2017; Habler and Griffiths 2017). A COR describes a type of orientation relationship between a suite of inclusions and their hosts. Different classes or types of COR are distinguished purely by geometrical criteria and, in particular, by their degrees of freedom. The COR concept is a data analysis tool and has no intrinsic genetic significance.

A specific COR describes pairs of crystals which have a specific and fixed orientational relationship to one another, for example as is found for true examples of epitaxial growth of one structure upon the substrate of another. A specific COR has no degrees of freedom, meaning that all inclusions belonging to the COR will have exactly the same orientation. Once allowance has been made for the ambiguity of indexing, poles of the crystallographic axes of the inclusions will all plot in exactly the same place in a stereogram based on the diamond host (Fig. 7a). Such CORs are well-known in minerals as the result of the solid-state exsolution of one mineral from another. One such example is the specific orientation of clinopyroxene grains within garnet inclusions in sub-lithospheric diamonds (Harte and Cayzer 2007) which are formed on decompression due to the decreasing capacity of the garnet to retain the Si associated with its majoritic component. Inclusions demonstrating a rotational COR have one specific direction of the inclusion parallel to one direction of the host diamond, but other directions will appear random as a result of rotation about the common axis. Crystallographic directions perpendicular to the common axis will plot on a single great circle in a stereogram (Fig. 7b). Rotational CORs therefore have one degree of freedom, the rotation around the common axis. Random CORs have two degrees of orientational variation, which is equivalent to saying that there is no specific orientational relationship between the inclusions and their host diamonds, which is immediately apparent in the scatter of the plots of crystallographic directions in a stereogram (Fig. 7c). Real measurements of inclusion orientations include experimental uncertainties, so the ideal cases of specific and rotational CORs are never obtained exactly, but only approximately. In addition, post-entrapment deformation of the diamond or modification of the inclusions may further scatter the orientation data.

The use of the COR concept does not imply anything about the mechanism by which the COR is developed for a particular type of mineral inclusion; it only describes the orientations of the inclusions. This approach, which we illustrate in the following examples, allows the analysis of inclusion orientation data to proceed in two well-defined and separated steps. First the type of COR or CORs is determined on the basis of statistical tests on the data, and these are generally not subject to significant uncertainty unless the number of data are not sufficient.

Then, once the statistical significance of the orientation measurements is clearly established, it can subsequently be interpreted in terms of the growth mechanisms of the diamonds and inclusions. If such interpretations subsequently prove to be incorrect, the underlying information about the COR remains valid and open to further alternative genetic interpretations. The COR concept also avoids the bias implicit in many earlier studies which often looked for evidence of specific CORs rather than describing the full statistical distribution of inclusion orientations present.

Olivine, clinopyroxene and garnet inclusions.Mitchell and Giardini (1953) in perhaps the first measurement of the orientations of inclusions in diamonds by diffraction showed that two elongated olivine inclusions in two diamonds had an orientation (010)olivine // (111)diamond and [101]olivine // [101]diamond. If generally true, this would constitute a specific COR and would probably imply some kind of epitaxial control on the orientations, as indeed claimed by Futergendler and Frank-Kamenetsky (1961) even though they also found three other olivines with significantly different orientations. Frank-Kamenetsky (1964) reported three more flattened olivines with the same orientation, but also found two other types of orientation. In recent years, individual measurements of at least 113 individual olivine inclusions in situ in diamond have been published from both EBSD (Neuser et al. 2015) and X-ray diffraction measurements (Nestola et al. 2014; Milani et al. 2016; Seryotkin et al. 2017; Sobolev et al. 2020). The as-published orientations appear random (Fig. 7c). When the randomization due to the indexing ambiguity is removed by plotting the b-axis of the olivine within the joint asymmetric unit of the diamond and olivine, they effectively fill the entire orientational space (Fig. 8a). This confirms that olivine inclusions in diamonds exhibit a random COR. It is also obvious from Figure 8a that a small proportion of olivines that are part of the random COR have orientations that are approximately in a specific alignment with the diamond, including the orientation reported by Mitchell and Giardini (1953). Such specific orientations will always be found, of course, as a subset of inclusions exhibiting a COR of more degrees of freedom.

The presence of a proposed specific COR can be easily tested with this data, by selecting the one symmetry-equivalent orientation of each measurement that is closest to the proposed COR. To test for the specific COR of Mitchell and Giardini (1953), the orientation with (010)olivine closest to (111)diamond is plotted on a stereogram. Figure 8b shows that the b-axes, equivalent to the (010) plane normals, do not plot systematically close to the [111] axis of the diamond; in fact, the angles between the two span the whole possible range of values from 0 to 54.7°. And the [101] directions of the olivine are also widely scattered and not coincident with [101] of the diamond. This excludes the proposed specific COR. Comparison of this plot with that expected for a rotational COR around [111]diamond (Fig. 7b) clearly also excludes the possibility that olivines in diamonds exhibit a rotational COR. This example demonstrates that distinguishing specific CORs from rotational CORs and rotational CORs from random CORs requires large numbers of inclusions to be measured, and conclusions must not be based on a handful of measurements.

The reason for olivines exhibiting a random COR, irrespective of the morphology of the inclusion, appears to be the recently-discovered ubiquitous presence of aqueous fluid films coating the inclusions, which must have been trapped along with the olivine at the time of diamond growth (Nimis et al. 2016). This fluid would prevent any direct interaction or bonding between the surfaces of the olivine and the diamond and thus prevent the development of any preferred orientation through a process such as epitaxy. The presence of the film also explains the absence of dislocations in the diamond adjacent to the inclusions (Agrosi et al. 2016). The incorporation of the fluid into the inclusion at the time of entrapment is due to the solid surfaces of olivine and diamond being mutually repulsive, as shown by DFT simulations (Bruno et al. 2016). In the presence of fluid, the olivine and diamond surfaces would therefore separate and fluid would be drawn in by capillary action to separate the olivine and diamond surfaces at the time of inclusion entrapment.

However, there are also many examples of two or more olivine inclusions in a single diamond that exhibit the same orientation as one another (Nestola et al. 2014; Milani et al. 2016; Seryotkin et al. 2017), and these inclusions also share the same composition as determined by in situ single-crystal diffraction (Fig. 4). Recent measurements of the orientations of garnet inclusions within lithospheric diamonds (Davies et al. 2018; Nestola et al. 2019b; Sobolev et al. 2020) show exactly the same characteristics as olivine inclusions; the overall COR between the garnets and diamond is also random (Fig. 8c), while iso-oriented groups of garnets occur in the same host diamond (Nestola et al. 2019b). Given that both olivine and garnet inclusions in diamond exhibit a random COR (Fig. 8), this grouping of orientations in one diamond requires another explanation that does not invoke interaction between the inclusion and the diamond or any influence of surface energies. Post-entrapment dissolution and re-precipitation could not lead to different isolated inclusions developing identical orientations. The only rational explanation proposed so far is that these iso-oriented inclusions were originally parts of one single crystal within the host rock prior to diamond growth (Nestola et al. 2014). Such single crystals would have been partially dissolved by the fluids from which the diamond grew in order to create space for the growth of the diamond, which then entrapped the remnants of the original crystals (Fig. 9). The mutual slight mis-orientations of iso-oriented inclusions, in the case of olivine, is consistent with the common observation of sub-grain boundaries in olivine within peridotite diamond host rocks (e.g., Baptiste et al. 2012).

Groups of two or three iso-oriented sets of inclusions suggest that diamonds nucleated on grain boundaries, or triple junctions, respectively, and then grew while dissolving the surrounding minerals (Fig. 9). Four differently-oriented groups of iso-oriented inclusions have not yet been found in a single diamond. The spacing and spatial arrangement of all such groups of iso-oriented inclusions is consistent also with the estimated grain size of minerals in mantle rocks, while the preferred orientation exhibited by olivines in a single diamond (but not with respect to the diamond) is consistent with the preferred orientation found in olivines in mantle peridotites (Nestola et al. 2014). There is insufficient published data to determine the nature of the COR of clinopyroxene inclusions with their host diamonds. But the observation that a single clinopyroxene inclusion with no specific orientation with respect to its diamond host did however have an identical orientation to an adjacent clinopyroxene in the garnet peridotite xenolith (Nestola et al. 2017), supports the model for diamond growth shown in Figure 9.

Such a scenario of dissolution of pre-existing mantle minerals to form the inclusions found in lithospheric diamonds means that, in strict application of the terminology, these inclusions are protogenetic with respect to their diamonds. Given that the datasets (Fig. 8) are from both inclusions with faces parallel to (111)-type planes of diamond and many other shapes, these data show that inclusion shape is not relevant for determining the genetic relationship of these inclusions to their host diamonds. The common occurrence of pseudo-octahedral shapes of inclusions is simply a reflection of the much lower surface energies of (111) faces of diamond compared to any other face of diamond. Thus, modern orientation data on olivine and garnet inclusions has shown that the morphology of the inclusions in diamond has no bearing on whether they are syngenetic or protogenetic. These conclusions open once again the question of whether ages and geochemistry of inclusions record the conditions and time of diamond growth. If the inclusion sizes are small enough then solid-state diffusion of elements could be sufficiently rapid to ensure re-equilibration of the inclusion with the diamond growth environment, including the resetting of isotopic systems used for dating (Davies et al. 2018; Nestola et al. 2019b). Whether this is also true for all inclusions is a matter for further research.

Ferropericlase inclusions. The limited published data on the orientations of these inclusions in sub-lithospheric diamonds suggest that they approximate a specific COR (Fig. 10a) in which the crystallographic axes of the cubic inclusions are very close to those of the diamonds (Nimis et al. 2018). If the COR was exact and the small scatter in orientations was just a consequence of measurement uncertainties, then all of the crystallographic axes would be scattered in a similar way when plotted on a stereogram. However, the [112] directions (and those close to it, such as [111]) of the ferropericlase inclusions and diamonds all coincide to within the 2° uncertainty in the measurement. This suggests that the ferropericlase inclusions exhibited a specific COR when they grew as a result of interfacial energy minimisation, in turn implying that no fluid film separated the inclusions and diamond during growth and that these inclusions are a product of the diamond-forming reaction (Nimis et al. 2018). This would mean that the compositions of such inclusions are not representative of bulk mantle ferropericlase, but only of the reaction itself. The interpretation is that a specific COR was subsequently transformed into a rotational COR as the result of post-entrapment plastic deformation of the diamonds (Nimis et al. 2018). These measurements are from just two diamonds, so many more measurements are required to determine whether all ferropericlase inclusions show this COR, or whether a range of CORs for ferropericlase are present in diamonds which would indicate a range of growth mechanisms and environments.

Graphite inclusions. Graphite is a common inclusion in diamond often associated with other minerals. This graphite is believed to be produced by internal graphitization (e.g., Harris and Vance 1972), a process perhaps enhanced by the presence of fluid films around silicate inclusions (Nimis et al. 2016). Graphite is more rarely found as isolated single-phase euhedral inclusions where evidence suggests that the diamonds may have nucleated on the graphite (Nasdala et al. 2005). Glinnemann et al. (2003) showed that out of 13 measured isolated graphite inclusions of the latter type, 12 have their c-axes within 4° of a <111> direction of the diamond. Insufficient data are given in the paper to completely define the relative orientations, so in drawing Figure 10b we have assumed that the reported misalignments are due to measurement uncertainties, and that the true relationship is exact. Graphite has hexagonal symmetry, so the a-axes of the graphite inclusions are perpendicular to the c-axes and therefore plot on a great circle around the [111]diamond, and appear to define a rotational COR. However, the 12 measured inclusions fall into two main groups, one with the <100> directions of the graphite within a few degrees of the <110> directions of the host diamond, and the other with the angle between these two directions being 34–38° (Fig. 10b). The fact that only one inclusion from 12 has an intermediate rotational angle suggests that these data might represent two distinct specific CORs, a conclusion supported by analysis of the structures of the interface planes of graphite and diamond which shows that these two orientations correspond to the best two structural matches (Glinnemann et al. 2003). Further careful measurements of graphite inclusions are needed to determine whether these conclusions are statistically significant, or instead whether graphite displays a full rotational COR. If only the specific CORs are found this would support the interpretation that the diamond nucleated on the graphite and subsequently entrapped it (Nasdala et al. 2005).

Magnesiochromite inclusions. Of the types of inclusions for which there are sufficient data to attempt a meaningful analysis, these show the most complex combination of CORs. The modern dataset consists of 58 inclusions in lithospheric diamonds from various Siberian kimberlites (Nimis et al. 2019; Sobolev et al. 2020) plus 7 from Damtshaa, Botswana and 6 from Panda, Canada (Nimis et al. 2019). A stereographic plot of the data (Fig. 11a) with the a-axis plotted in the common asymmetric unit appears, at first sight, very similar to the random CORs exhibited by garnets (Fig. 8c). But closer inspection shows that there is a clustering of inclusions with their a-axes close to the a-axis of the diamond (compare Fig. 11a carefully with Fig. 8a and 8c) and there appears also to be a clustering of the [111]mgchr directions towards [111]diamond. This is confirmed by replotting the data with the orientation that has [111]mgchr closest to [111]diamond; then 33 out of the 71 inclusions (46%) have the angle between these two directions less than the estimated measurement uncertainties of 4° (Fig. 11b). This means that their octahedral faces are parallel within error to an octahedral face of their host diamond. These inclusions form a rotational COR which is especially evident when the <110> directions which lie in the (111) plane of the inclusions are plotted on the stereogram (Fig. 11b). With the limited numbers of measurements, the group whose [111] mgchr are more than 4° from [111]diamond (Fig. 11c) appear to represent a random COR. The interpretation of these data is that magnesiochromite inclusions exhibit two CORs; one random and one rotational. The random COR is entirely consistent with the observation of fluid rims around these inclusions (Nimis et al. 2016; Nimis et al. 2019) that would, as for olivine and garnet inclusions, prevent formation of an exact COR through epitaxial growth. The only mechanism so far proposed for the creation of the rotational COR is that it is the result of mechanical impingement. Although spinels in mantle rocks do not normally display octahedral shapes, dissolution by the fluids present at the time of diamond growth will tend to dissolve the higher-energy faces of the magnesiochromites which will then tend to develop a habit that is octahedral. In the presence of fluid and under lithostatic load the tendency will be for the dominant octahedral faces of adjacent spinels and diamonds to become parallel to one another, simply by being pushed together mechanically. The apparent tendency for some of these inclusions to exhibit an approximate specific COR can then be explained as the mechanical interaction of a spinel remnant not with a single plane of the diamond, but with etch pits in, or growth steps on, the (111) surface of the diamond (Nimis et al. 2019).

Evidence for inclusion over-pressures in the form of induced birefringence haloes in diamond adjacent to inclusions has been known for many years (e.g., Sorby and Butler 1869; Meyer quoted in Cohen and Rosenfeld 1979). Inclusions are over-pressured because diamond is so much stiffer than the silicate minerals, or fluids, found as inclusions, so the diamond expands less on exhumation than would a free crystal of the inclusion mineral. Therefore, when the diamond is in the laboratory the inclusion crystal is constrained to have a smaller volume than a free crystal of the same mineral would have at room conditions. The inclusion will therefore exhibit an over-pressure. Elastic thermobarometry (or piezobarometry) is the process of determining the P and T of entrapment of inclusions from the over-pressure measured in them when their host mineral is at room conditions. It has recently provided new constraints on the P–T conditions attained during metamorphism (e.g., Zhong et al. 2018; Gonzalez et al. 2019; Korsakov et al. 2020; Kosminska et al. 2020). When applied to inclusions in diamond it can provide estimates of their depth of formation and growth beyond that available from identification of the mineral phases present as inclusions. Elastic thermobarometry has been applied to various types of inclusions in diamonds with variable degrees of apparent success (e.g., Harris et al. 1970; Cohen and Rosenfeld 1979; Izraeli et al. 1999; Sobolev et al. 2000; Howell and Nasdala 2008; Howell et al. 2010, 2012a; Nestola et al. 2011; Anzolini et al. 2016; Kueter et al. 2016; Anzolini et al. 2018) that reflect the concurrent development of the theory and techniques for characterizing the stress and strain states of inclusions in diamond.

The application of elastic thermobarometry relies on several fundamental assumptions. The first is that when the inclusions are trapped by the host diamond, both experience the same temperature and the same hydrostatic pressure. This is not unreasonable, since the evidence of the ‘fluid’ films around many inclusions (Nimis et al. 2016, 2019; Smith et al. 2016, 2018), the presence of fluid inclusions (e.g., Smith et al. 2015; Gubanov et al. 2019) and the resorption of inclusion minerals (Nimis et al. 2019) indicate that diamonds large enough to host macroscopic inclusions probably grow in fluid-saturated environments. The established theory also assumes that the inclusion space in the diamond is filled completely and exactly by the inclusion crystal without any void space, but the theory can be adapted to account for the presence of fluid as we describe below. The basic analysis further assumes that all volume changes in the system following entrapment are elastic. If there has been plastic deformation of the diamond after diamond growth and inclusion entrapment, some or all of the inclusion stress may be released. Similarly, the amount of inclusion stress released by cracking of the diamond around the inclusion depends on whether the inclusion also contains fluid in addition to the mineral. Cracks and plastic deformation can be detected by a variety of techniques, including X-ray tomography (e.g., Rustioni et al. 2015; Mazzucchelli et al. 2016; Nimis et al. 2016; Wenz et al. 2019) and topography (e.g., Agrosi et al. 2017). In all of these cases, the measured inclusion pressure will be reduced but this does not invalidate the use of elastic thermobarometry. Proper characterization of the inclusion often allows either the minimum entrapment pressure to be estimated, or the conditions of the resetting of the inclusion stress to be calculated. Phase changes in the inclusions after entrapment to lower-pressure polymorphs may result in the inclusion pressure being buffered (e.g., Gillet et al. 1984; Perrillat et al. 2003; Anzolini et al. 2016) to a higher value than expected for a single-phase inclusion of the higher-pressure polymorph present at the time of entrapment, and also requires more elaborate analysis.

The over-pressure in the inclusion leads to elastic deformation of the diamond host around the inclusion and this is what gives rise in otherwise optically-isotropic diamond to the strain-induced birefringence which has been used to characterize the inclusion stress state (Howell and Nasdala 2008; Campomenosi et al. 2020a). The deformation of the host diamond leads to elastic stresses within it which balance the over-pressure in the inclusion. The retention of stress or over-pressure inside inclusions is thus the result of a mechanical equilibrium between the host diamond around the inclusion and the inclusion itself. Removal of this confining stress by removal of the diamond by polishing it away, partially or completely, or by breaking out the inclusion, will lead to partial or complete loss of inclusion stress (e.g., Campomenosi et al. 2018; Mazzucchelli et al. 2018; Zhong et al. 2019). For example, Glinnemann et al. (2003) showed that graphite inclusions near the surface of their host diamonds exhibit less compression than those deep within their hosts, while garnet inclusions close to the outside surfaces of the diamond show no significant difference in Raman shifts compared to garnets measured at room pressure (Liu et al. 1990). Therefore, unlike the determinations of the orientation of inclusions or their crystal structure that we have discussed in the preceding sections, meaningful measurements of inclusion stress can only be made on fully-buried inclusions at least three inclusion radii from the diamond surface, other inclusions and any other defects (e.g., Zhang 1998; Tajčmanová et al. 2014). It is therefore not possible for EBSD measurements on partially exposed inclusions to return any meaningful data with respect to the true inclusion stress or strain state (Campomenosi et al. 2018; Zhong et al. 2019).

Forward problem.The development of stress and strain in an inclusion trapped within diamond is an elastic effect and therefore, by definition, reversible. We first describe the forward problem, i.e., what happens to the inclusion from encapsulation at depth until its diamond host is transported to the Earth’s surface. Diamond is cubic and extremely stiff compared to the minerals found within it as inclusions, so we can use the theory developed for elastically isotropic systems as an initial basis for most analyses of inclusions in diamonds. The various extensions and adaptations of the isotropic theory are described in the ‘examples’ below.

We start our analysis by considering an isotropic mineral inclusion in diamond, and we assume that the inclusion is a single mineral phase which does not undergo reactions or phase transformations. At the time of diamond growth and entrapment of the inclusion, both the diamond and the inclusion are at the same P and T. The inclusion fills the void in the diamond completely, so the diamond void and inclusion also share the same volume, V. If, after entrapment, the P and T changes in such a way that the natural expansion and contraction of a free crystal of the inclusion exactly matches that of the host diamond, the inclusion will continue to exactly fit the void space in the diamond without the application of any additional stress. This line in P–T space is called an isomeke (Rosenfeld and Chase 1961; Adams et al. 1975), and is defined by the difference in the volume thermal expansion coefficients (α_inc and α_dia) and compressibilities (β_inc and β_dia) of the inclusion phase and the diamond. The volume compressibility of a phase is the inverse of its bulk modulus, thus; $βinc=1Kinc=−1V(dVdP)T$⁠. The slope of an isomeke through any P–T point is:

Isomekes can therefore be calculated from the EoS of the diamond and of the inclusion phase. Isomekes have also been called ‘zero volume difference’ (ZVD) curves (Izraeli et al. 1999) and ‘isovolume loci’ (Barron 2003). An example of an isomeke of mantle-composition olivine with diamond is given in Figure 12. It has a positive slope because, as for most silicates, olivines have volume thermal expansion and compressibility coefficients larger than diamond (α_inc > α_dia and β_inc > β_dia), Table 2. Although diamond is very stiff, neither α_dia nor β_dia are zero, and the isomeke of the olivine with diamond is not the same as an isochor of either phase (Fig. 12).

The importance of the isomeke concept is that for P and T conditions on the entrapment isomeke, the inclusion will have the same pressure as the external pressure on the diamond host crystal. However, the typical eruption path of diamonds in kimberlites involves rapid pressure decrease and is obviously not along an isomeke (Fig. 12). When the external P and T drop below that of the diamond–inclusion isomeke through entrapment conditions, the inclusion crystal wants to expand more than it is allowed by the diamond because β_inc > β_dia, but it is restricted to the volume change of the stiffer diamond. Therefore, when the diamond is erupted the inclusion develops an over-pressure with respect to the external pressure applied to the diamond (Fig. 12) at all external pressures and temperatures that lie below the entrapment isomeke (e.g., Angel et al. 2014b; Ferrero and Angel 2018). Room conditions lie below the entrapment isomekes of silicate inclusions in diamonds (e.g., Fig. 12), which is why they exhibit positive remnant pressures.

The first step in the calculation of the over-pressure in the inclusion is to use the volume change of the host diamond from entrapment to room conditions (calculated from its P–V–T Equation of State or EoS) in the EoS of the inclusion mineral. This is just a thermodynamic calculation, so the resulting pressure has become designated as P_thermo (previously PI* in Angel et al. 2014b). It only depends on the EoS of the two phases, and does not depend on the size, shape or any other geometrical aspect of the system. As a simple example, consider entrapment of a mantle-composition olivine at 1230 °C and 6.2 GPa. At these conditions the molar volume of diamond calculated from its EoS (Table 2) is 3.420 cm³/mol, the same as its molar volume at room conditions; these entrapment conditions lie on the isochor of diamond that passes through room conditions (Fig. 12). Therefore, for this particular example, on recovery to room conditions the void in the diamond that contains the inclusion has the same volume as it had at entrapment. However, the olivine at entrapment is calculated to have V = 43.34 cm³/mol, 1.26% smaller than the volume of a free olivine at room conditions of 43.89 cm³/mol. Thus, because the olivine inclusion is restricted by the diamond to maintain the same volume as it had at entrapment, it is under 1.26% volume compression at room temperature, corresponding to a remnant pressure in the inclusion of ~1.64 GPa (Fig. 12).

This is a fictive state that is never attained in real systems because there is a pressure difference between the host, at ambient pressure, and the inclusion at P_thermo. The system is not in mechanical equilibrium. The pressure P_thermo is greater than ambient pressure for nearly all silicate inclusions in diamonds, so it drives expansion of the inclusion and reduces the inclusion pressure by an amount ΔP_relax (Fig. 13). By convention (Angel et al. 2014b) a relaxation that decreases the inclusion pressure has ΔP_relax < 0. The final inclusion pressure, P_inc, which is what is always measured, is therefore related to the thermodynamic pressure by:

The expansion of the inclusion crystal compresses the host in a radial direction, thereby increasing the radial stress in the host. Mechanical equilibrium is attained when the pressure in the inclusion becomes equal to the radial stress in the host at the wall of the inclusion (Fig. 13c). The radial expansion of the inclusion also causes expansion in the diamond in directions parallel to the inclusion surface. Thus, it is this elastic relaxation that creates the stress and strain fields around the inclusion which can be detected as the strain-induced birefringence (e.g., Sorby and Butler 1869; Howell and Nasdala 2008; Howell et al. 2010), or as changes in the Raman spectra of the host diamond (Nasdala et al. 2003; Campomenosi et al. 2020a). For most inclusions, the stresses drop off in proportion to the distance from the center of the inclusion cubed (Fig. 13c). At distances of three or more radii from the inclusion, the induced stresses in the diamond become insignificant. At shorter distances, these stresses can interact with any inhomogeneity in the host diamond, such as another inclusion, cracks, external surfaces etc. These interactions can lead to a modification of the elastic relaxation that therefore changes the final inclusion pressure away from that expected for an isolated inclusion.

Diamond is cubic and therefore expands and contracts isotropically in response to changes in P and T. But, it is elastically anisotropic because elasticity is a 4^th-rank tensor property (Nye 1957). This means that the elastic properties of diamond change with direction in the crystal. Therefore, the amount of elastic relaxation (Fig. 13) around an inclusion will be different in different directions, leading to small deviations in the inclusion stress and strain from isotropic. However, because diamond is so stiff the effect of its elastic anisotropy on P_inc is so small that it can be ignored and the diamond host can be treated as being elastically isotropic. The amount of relaxation of an isolated isotropic inclusion then depends on the contrast between the elastic properties of the inclusion and the diamond, and the shape of the inclusion. For a spherical inclusion in diamond, the fractional volume change of the inclusion due to relaxation (Zhang 1998; Angel et al. 2017) is:

when the diamond is at room pressure (considered to be effectively zero). The shear modulus of diamond at room conditions is G_dia = 535 GPa. The pressure change in the inclusion, due to the relaxation, ΔP_relax, can be calculated from this volume change and the EoS of the inclusion (Angel et al. 2017). For the example of olivine shown in Figure 12, ΔP_relax = −0.26 GPa for a spherical inclusion, and the final measurable inclusion pressure is P_inc = 1.38 GPa. For small inclusion pressures a rough estimate can be made by using the definition of the bulk modulus to obtain:

We have described the theory of piezobarometry in terms of recovery of the host diamond to room conditions. But the theory applies equally to the calculation of inclusion pressures at any other conditions, for example while the diamond is still deep in the Earth, or while the diamond and its inclusion are heated in the laboratory (Harris et al. 1970). For non-ambient conditions the volume relaxation is given by $3(Pinc−Pext)4Gdia$ and the value used for G_dia in Equation (4) must be changed to that for the external P and T. Equation (5) does not apply, so the ΔP_relax must be calculated from this volume change and the EoS of the inclusion at the temperature of measurement (Angel et al. 2017).

Inverse problem.For an isotropic spherical inclusion the determination of possible entrapment conditions consists of reversing the steps described in the preceding section. The inclusion pressure P_inc is first measured and the ΔP_relax is calculated through Equation (4) and the EoS of the inclusion, to obtain P_thermo. This is the one and only constraint on the entrapment conditions from the measurement of the inclusion. It is therefore not possible to determine a unique P and T of entrapment (two variables) from this single observation. But it is possible to determine the entrapment isomeke, the line in P–T space from which the inclusion would exhibit the same P_thermo and P_inc.

The point on the entrapment isomeke at room temperature can be calculated by considering the isothermal compression of the diamond host and the inclusion, starting from a common volume (Fig. 14) at which the diamond is at room pressure, and the inclusion is at P_thermo. Because the diamond is stiffer than the inclusion, the PV curve of the diamond is shallower than that of the inclusion and the two EoS curves converge with increasing pressure. At the point at which they cross both the diamond and the inclusion have the same P and V (and the same T). This is a point on the entrapment isomeke. This particular pressure on the entrapment isomeke at the final temperature T_end of the inclusion measurement has been called P_foot (Angel et al. 2014b). Once P_foot is determined, the full entrapment isomeke corresponding to the measured P_inc can be calculated from the EoS of the two phases, either step-wise by using Equation (2), or more directly by finding the P at each T at which the fractional volume change of both the inclusion and the diamond from P_foot, T_end are exactly the same (e.g., Angel et al. 2020).

Whether an inclusion phase makes a good geobarometer depends on the slopes of its isomekes with diamond. These are determined by the contrast between the elastic properties of the inclusion and those of diamond (Eqn. 2). Diamond is significantly stiffer than all reported inclusion phases, so the thermal expansion coefficient of the inclusion phase is the major determinant of the isomeke slope through its contrast to that of diamond, (α_inc – α_dia). Thus, phases with a small α_inc similar to that of diamond, such as coesite and graphite, have isomekes that are relatively flat in P–T space (Fig 15) that make them good geobarometers. Inclusions of these phases also therefore have P_foot similar to the pressures of entrapment (Fig. 15), although the measured P_inc is significantly lower because the elastic relaxation term is relatively large. Thus, a coesite inclusion trapped at the same time and conditions as an olivine or a garnet should exhibit a much higher residual inclusion pressure when measured in the laboratory. In principle, the best geothermometers would be stiff inclusion phases with a β_inc similar to that of diamond, with a large α_inc so as to make the isomekes with diamond as steep as possible (Eqn. 2). Unfortunately, stiff minerals tend to display small thermal expansion.

Calculation of entrapment conditions of an inclusion requires three steps. The first step is the determination of the stress or the strain state of the inclusion crystal when the host crystal is at known conditions of P and T, normally room conditions. This can be achieved by in situ diffraction or Raman spectroscopy. The second step is to calculate the amount of stress and strain relaxation that the inclusion has undergone as a consequence of the contrast between the stress in the inclusion and the external pressure on the host, and thereby calculate the strain or stress that the inclusion would have in the fictive unrelaxed state. The third step is to use this information about the unrelaxed state to calculate possible entrapment conditions of the inclusion by using the thermodynamic equations of state of the two minerals.

Measuring inclusion pressures.The fundamental requirement for success in thermobarometry is an accurate and precise measurement of the stress in the inclusion. It is important to understand that even in an isotropic inclusion its pressure is never measured directly as a pressure. Instead, the pressure is inferred from an in situ measurement of some physical property of the inclusion. Most commonly, the Raman spectrum of the inclusion is measured, and the shift in the frequencies of the Raman lines from their values at room pressure is used to calculate the inclusion pressure from their known rate of shift with pressure. This is strictly incorrect as we describe in detail below. More rarely, the unit-cell parameters and volume of an inclusion are determined with an X-ray diffraction measurement, and the inclusion pressure inferred from their known variation with pressure (e.g., Nestola et al. 2011). While the determination of orientation of the inclusion crystal by diffraction is not sensitive to small offsets of the crystal from the center of the diffractometer, the unit-cell parameters and volume values can be significantly changed by failure to center the inclusion exactly. It is therefore crucial to first follow the entire procedure for centering the inclusion, that is outlined in the sections above, before making the measurements to determine its unit-cell parameters and volume.

It is clear that diffraction measurements of the unit-cell parameters of an inclusion do not directly measure its stress state. But the difference from the values of the unit-cell parameters measured on a free crystal at room conditions can be turned into the average strain experienced by the inclusion trapped in the host relative to a free crystal of the same material at room conditions. For example, the strain along the crystallographic a-axis of an inclusion crystal measured as a_inc can be defined as:

where a₀ is the unit-cell parameter the inclusion crystal would have when it is free of the host and at room pressure and temperature. In general, strain is a second rank tensor with six components (e.g., Nye 1957). The three normal components (ε_1, ε₂ and ε₃ in Voigt notation) describe the change in dimensions in three mutually perpendicular directions in the inclusion, and for crystals of orthorhombic symmetry and higher can be defined purely in terms of the cell parameters a, b, and c, as in Equation (6). Three shear components (ε_4, ε₅ and ε₆) characterize the strains in three mutually-perpendicular planes. The definition of the individual strain components in monoclinic and triclinic crystals in terms of the unit-cell parameters depends on the choice of axial system for the tensor (e.g., Schlenker et al. 1978; Knight 2010). Equation (6) defines the Lagrangian infinitesimal strains, which are usually suitable for the analysis of the remnant strains in inclusions when the host is at room conditions because the strains are small.

The volume strain ε_V = ε₁ + ε₂ + ε₃ is defined as the fractional volume difference:

What is rarely made clear in the geological literature is that Raman shifts also do not measure stress or pressure. The change in the peak position or wavenumber of a Raman peak is determined not directly by the pressure or the temperature but is primarily a function of the strains applied to a crystal. A change in P or T changes the unit-cell parameters and this change (or strain, Eqn. 6) determines the wavenumber of the Raman mode. This concept was first established by Grüneisen (1926) for isotropic solids, who defined what is now known as the thermal or thermodynamic Grüneisen parameter which relates the change in internal energy of an isotropic solid to the change in pressure along an isochor. The quasi-harmonic approximation (QHA) for the thermodynamics of solids amounts to assuming that the thermal Grüneisen parameter is a constant which does not change with either P or T. For anisotropic crystals one can define a phonon or mode Grüneisen tensor γ^m (Ziman 1960; Key 1967; Cantrell 1980; Angel et al. 2019) for each Raman-active mode. This is a second rank tensor which defines the fractional change in the phonon wavenumber, Δω/ω, arising from the strain applied to the crystal. Written out in full in Voigt (1910) notation:

This equation means that the changes in the Raman peak positions depend on all of the strains in three dimensions experienced by the crystal, not just the volume. The values of $γim$ are different for different modes. If the crystal follows the QHA then the values of the six components, $γim$⁠, of a mode Grüneisen tensor will be constants and will not change with pressure and temperature, or with the magnitude of the strains. The use of the mode Grüneisen tensors has several advantages over using polynomials in P and T to describe the variation in peak shifts (e.g., Schmidt and Ziemann 2000; Ashley et al. 2017; Morana et al. 2020). First, the Grüneisen tensor has a sound physical basis, and reflects the underlying physics that the shifts due to P and T are not independent in materials that exhibit quasi-harmonic behavior. Second, they account for deviatoric strains, as required for the analysis of host–inclusion systems, with the same approach as used for hydrostatic conditions. And, third, less coefficients are needed to describe the peak shifts than are often required for polynomials in P and T.

For spherical inclusions the strain is homogeneous (Eshelby 1957) irrespective of the symmetry of the host or inclusion crystals. In a faceted inclusion the strains vary across the inclusion volume, reflecting stress concentration at edges and corners. As a consequence, the Raman shift of a single band varies across a shaped inclusion (e.g., Zhukov and Korsakov 2015; Campomenosi et al. 2018; Murri et al. 2018), and a single measurement at one random point in an inclusion cannot therefore be interpreted as being representative of the inclusion stress and strain state. But the shear strains at the center of a faceted inclusion are normally orders of magnitude less than the normal strains (e.g., Mazzucchelli et al. 2018). And, for spherical inclusion crystals in diamond with symmetries that are not triclinic or monoclinic, $γ4m=γ5m=γ6m=0$ by symmetry. In all of these cases the shifts of the Raman lines measured at the center of the inclusion, compared to a free crystal at room conditions, become governed by:

with the strains of the inclusion also being defined relative to a free crystal of the same mineral at room conditions. Equation (9) applies to triclinic, and monoclinic inclusions when they experience no shear strains, and to orthorhombic inclusion crystals. For higher symmetries, discussed below and in detail in Angel et al. (2019), some or all of the components of the mode Grüneisen tensor are equal.

Therefore, the strains in an inclusion can be determined if the values of the components of its mode Grüneisen tensors are known, and sufficient Raman lines are measured to enable the strains ε_i to be determined through Equation (9). The values of the tensor components are not easy to obtain experimentally, because this requires different strain states to be imposed on a crystal. In general, it is not sufficient to measure the Raman shifts under hydrostatic pressure, because this only generates a single strain state and not the independent ranges of strains required. A practical alternative is to use computer simulations of the mineral based on density-functional theory (DFT) to calculate the Raman spectra of a crystal under a range of deviatoric strains from which the tensor components are determined by a least-squares fit of Equation (9). Tensor components have been calculated in this way for quartz (Murri et al. 2018; Murri et al. 2019), zircon (Stangarone et al. 2019) and rutile (Musiyachenko et al. 2021). They can be used in the Strainman program (Angel et al. 2019) to determine the strains in the inclusions from the difference in the Raman line positions measured at the center of the inclusion from the positions of the same lines measured on a free reference crystal of the same composition in air (Murri et al. 2018; Alvaro et al. 2020; Campomenosi et al. 2020b; Musiyachenko et al. 2021). In particular, care must be taken to only measure isolated inclusions as the presence of surfaces, cracks or other inclusions within 3 inclusion radii will modify the elastic relaxation (Fig. 13) and thus the final inclusion pressure. Entrapment conditions calculated from non-isolated inclusions will therefore be incorrect.

Calculating relaxation and entrapment.The volume relaxation of a spherical isotropic inclusion in an isotropic host crystal at room pressure is given by Equation (4) and can be converted into a ΔP_relax via the EoS of the inclusion phase, from which P_thermo can be calculated. However, inclusions in diamonds are never spherical, and the shape of the inclusion changes the amount of relaxation. Finite element model (FEM) calculations have been used to model the shape effects. Mazzucchelli et al. (2018) defined a shape factor Γ that can be used to ‘correct’ the pressure measured by Raman spectroscopy at the center of the inclusion to the one expected for the same phase trapped at the same conditions but as a spherical inclusion:

If the $Pincmeasured$ is determined from X-ray diffraction measurements, this is an average over the entire inclusion volume, and a similar factor can be calculated from a FEM of the inclusion (Anzolini et al. 2019).

A shape factor of Γ = 0 represents an inclusion whose relaxation is the same as that of an isotropic sphere. The biggest influence on the shape factor Γ is the aspect ratio of the inclusion; platy inclusions have larger shape factors than more equant ones (Mazzucchelli et al. 2018). The shape factor for an octahedral inclusion of pyrope garnet in diamond (Mazzucchelli, pers. comm.) is −0.026, which gives a 2.6% correction to P_inc which is insignificant compared to the uncertainties associated with the measurement of inclusion pressures. However, for plate and needle-shaped inclusions the values of Γ become significant (Table 3). Values of the shape factor for other inclusions in diamond have not been calculated, but octahedrally-shaped inclusions of olivine, coesite and walstromite-structured CaSiO₃ (now known as breyite, Brenker et al. 2021) which are significantly softer than garnets will have smaller correction factors as indicated by the trends found by Mazzucchelli et al. (2018).

When the inclusion is a non-cubic mineral, or a cubic mineral with significant elastic anisotropy, then the relaxation of the stress itself is anisotropic and depends not only upon the contrast of the elastic properties but also on the orientation of the inclusion crystal with respect to both the host and the shape of the inclusion. Thus, the relaxation of a prismatic inclusion of olivine in diamond (e.g., Mitchell and Giardini 1953) will depend on which crystallographic axis olivine is parallel to the prism axis, and the amount of strain and stress relaxed in each crystallographic direction will be different. The characterization and quantification of this relaxation requires the numerical calculation of the relaxation tensor for each combination of orientation and shape for each host and inclusion pair of interest (Mazzucchelli et al. 2019). This tensor can be used to correct the measured inclusion stress or strain to that expected from a spherical inclusion and subsequent calculations can be performed based on the isotropic analysis for spherical inclusions (Bonazzi et al. 2019).

Therefore, after measurement of the inclusion pressure it should be corrected for shape effects, and the $Pincspherical$ can then be used in the EosFit-Pinc program (Angel et al. 2017) to calculate the relaxation, the P_thermo, the P_foot and the entrapment isomeke. The EosFit-Pinc program, available for free download at www.rossangel.net, has a simple graphical user interface and is not restricted to specific inclusion phases because the EoS can be chosen and loaded by the user. Most common forms of P–V–T EoS are supported by the EosFit program suite (Angel et al. 2014a), which gives the user the opportunity to explore the effect on the predicted entrapment conditions of issues such as uncertainties and differences between various published EoS.

Simple linear EoS models have been used in the past for inclusion calculations (e.g., Zhang 1998; Barron 2003) in which the properties such as thermal expansion and compressibility are assumed not to change significantly with temperature and pressure. This can be a good approximation for the conversion of the ΔV_relax of the inclusion into ΔP_relax, because the pressures involved are small compared to the bulk moduli. However, it results in calculated isomekes that are linear in P–T space, whereas they are actually curved (e.g., Fig. 12) as a result of the stiffening of solids and decrease in thermal expansion of solids with increasing pressure (and the counter effects with increasing temperature). This is especially important for inclusions in diamond, as its room-temperature thermal expansion coefficient is less than 20% of its value at lithospheric depths. This is one of the main factors leading to curved isomekes between diamonds and their inclusions. If the variation in properties was ignored and the isomeke shown in Figure 12 was extrapolated from P_foot as a straight line with a slope equal to its initial slope (i.e., with constant compressibilities and thermal expansivities for both phases), then it is clear that the entrapment pressure at 1230 °C would be overestimated by more than 1 GPa.

The calculation of an isomeke from a single inclusion with the isotropic model only provides a line of possible entrapment conditions in P–T space (e.g., Fig. 12). We emphasize again that the entrapment isomeke does not represent the path of exhumation of the diamond (Fig. 12), but it is only a set of P–T conditions at which the inclusion may have been entrapped; if it has not been modified by other processes following entrapment, the inclusion pressure P_inc holds no information about the exhumation path. Further information is required to obtain a pressure of entrapment and thus depth of diamond growth. One possibility is to use the intersection of the isomeke with a mantle geotherm (e.g., Cohen and Rosenfeld 1979; Anzolini et al. 2019) and a minimum estimate of entrapment pressure can be obtained from the point where the entrapment isomeke enters the stability field of diamond (Graham and Cybriwsky 1981). If two different phases are present as separate inclusions in a diamond and belong to the same growth stage, then the intersection of their isomekes provides an estimate of the P and T of entrapment, a calculation that can be readily performed in the EntraPT software (Mazzucchelli et al. 2021).

The EntraPT program (Mazzucchelli et al. 2021), freely available at www.mineralogylab.com, also provides further shape corrections based on finite element modeling (Mazzucchelli et al. 2018), and the calculation of the relaxation tensor for non-ellipsoidal elastically anisotropic inclusions in an isotropic host (Mazzucchelli et al. 2019; Morganti et al. 2020). EntraPT implements the methods applied by Bonazzi et al. (2019) to correctly calculate the mean stress in inclusions from the measured strains, and then to calculate entrapment isomekes with the isotropic model including the propagation of uncertainties, providing a clear workflow through the complex issues surrounding the handling of data from shaped anisotropic inclusions.

The key factor in measuring and interpreting inclusion stress states is the symmetry of the stress field in the inclusion when it is measured in the laboratory. It has important and different consequences for all three steps required to calculate the entrapment conditions of any inclusion crystal. Given that diamond is so much stiffer than the inclusions, the symmetry of the stress field depends only on the symmetry and the shape of the inclusion crystal. The following examples of common inclusions in diamonds have therefore been divided by the symmetry of the inclusion to illustrate how to use measurements of their remnant pressures to determine the depths of their entrapment by the growing diamond. In particular, we emphasize the analytical and computational procedures necessary to overcome, or work around, the fact that we are trying to use a model for inclusion pressures that is based on assumptions of isotropic properties and shape for inclusions that are neither spherical nor elastically isotropic.

Uncertainties in the EoS parameters of diamond and the inclusion phase contribute to uncertainties in isomeke, and thus entrapment, pressures. However, because diamond is so stiff, and its EoS parameters are relatively well determined (Table 2) the uncertainties in them contribute far less to uncertainties in entrapment pressures than do those arising from the uncertainties in measuring P_inc (Angel et al. 2015a). Uncertainties in published EoS of silicates vary considerably, as does their sensitivity to compositional variation in solid solutions, so the consequences for calculations of entrapment conditions must be evaluated carefully for each mineral system.

Cubic inclusions. When a cubic crystal is subject to hydrostatic stress or a change in temperature its size changes equally in all directions; the strain (Eqn. 6) is isotropic. Consequently, the application of an isotropic strain by a diamond on a spherical inclusion with cubic symmetry creates a hydrostatic stress in the inclusion, and the symmetry of the inclusion remains cubic. As explained above, the influence of the elastic anisotropy of diamond on the relaxation is extremely small and well below the resolution of measurements (Mazzucchelli et al. 2019), so the final stress state can be treated as hydrostatic. The pressure in the inclusion, P_inc, can therefore be determined by measuring its unit-cell volume by diffraction and using the known EoS for the mineral. Or the pressure can be determined by measuring the shift of Raman lines of an inclusion.

For cubic crystals the components of the mode Grüneisen tensors are $γ1m=γ2m=γ3m$⁠, and the other components are zero, so Equation (9) reduces to:

This says that the Raman shift depends only upon the total volume change. Therefore P_inc can be determined from the measured shifts of its Raman bands from those of a free crystal of the same composition at room pressure, by using calibrations from high-pressure hydrostatic experiments. Note that the positions of Raman bands of free crystals that are solid solutions, such as garnets, can depend significantly on composition, as do their dependence on pressure and this must be allowed for in the analysis as well as the effects of peak broadening on apparent peak positions (Gillet et al. 2002; Kalugina and Zedgenizov 2020). Hydrostatic conditions in a cubic inclusion can be verified by comparing the pressures indicated by different Raman bands; if they differ significantly this is an indication that non-hydrostatic stresses are significant, and the inclusion should be treated as one of lower symmetry.

The P_inc should then be corrected by the appropriate shape factor (e.g., Table 3) to obtain $Pincspherical$ from which the entrapment isomekes can be calculated. Entrapment isomekes for several garnet endmembers are shown in Figure 16 which illustrates several issues critical for the correct interpretation of inclusion pressures as entrapment conditions. Most importantly, the P_inc expected for different garnet compositions trapped at the same pressure and temperature differs by up to 1 GPa. For example, a pure pyrope trapped at 6.1 GPa and 1200 °C will exhibit a P_inc = 0.1 GPa when measured at room conditions, whereas a pure almandine will have P_inc = 1.1 GPa. A measured P_inc = 0.5 GPa corresponds to entrapment pressures at 1200 °C of 6.9 GPa for pyrope, 5.6 GPa for grossular, and 4.8 GPa for almandine, a difference of more than 60 km in inferred depth of entrapment. Second, the thermodynamic calculations suggest that any inclusion trapped at temperatures higher (and pressures lower) than the zero-pressure isomeke should exhibit a negative P_inc. This would require bonding between the internal diamond surface and the inclusion to stretch the inclusion to a volume greater than that of a free crystal at room conditions. However, no negative P_inc values have ever been reported for inclusions in diamond. Therefore, such inclusions must expand within the diamond to their normal volume at room conditions, but no more, and they exhibit no pressure when measured at room conditions. Figure 16 therefore shows that pyrope, and pyrope-rich inclusions, cannot be used to determine entrapment conditions along geotherms corresponding to surface heat flows of more than ~37 mW/m². By contrast, grossular and almandine garnets will still exhibit P_inc > 0 for these entrapment conditions (Fig. 16).

The reason for these significant changes in P_inc with composition can be understood in terms of the slopes of the isomekes (Eqn. 2) shown in Figure 16. The total range of bulk moduli of garnets is < 20 GPa (Table 2), which is < 10% of the difference from the bulk modulus of diamond (444 GPa). Therefore, the difference in compressibilities β_inc – β_dia between different garnet inclusions and diamond varies by less than 10%. This contrast also controls the pressure spacing of isomeke lines, which are therefore very similar for all of the garnets with diamond (Fig. 16). In contrast, thermal expansion coefficients, a, of garnets vary between ~2.0 and 2.5 × 10⁻⁵ K⁻¹. Given that α_dia = 0.27 × 10⁻⁵ K⁻¹, this means that the slopes of the isomekes vary because of a 25% change in (α_inc – α_dia). Garnets with lower thermal expansion coefficients, such as almandine, therefore have shallower isomekes that result in higher P_foot values and thus P_inc pressures, than pyrope that has higher thermal expansion. In fact, the influence of bulk modulus on P_inc for pyrope–almandine–grossular garnets is insignificant compared to the influence of thermal expansion.

It is therefore critical to determine the composition of garnet inclusions in order to obtain correct estimates of their entrapment depths, for example by in situ X-ray diffraction as described above. The composition also allows the EoS to be determined from the known variation of EoS parameters of garnets with composition; for example, within the pyrope–almandine–grossular garnets, the bulk modulus can be considered constant for these purposes, and the thermal expansion coefficient appears to vary approximately linearly with composition (e.g., Bosenick and Geiger 1997). Once the appropriate EoS is established from the composition, the P_inc can be determined from the measured unit-cell volume, and the volume of a free crystal of the same composition determined from the systematics of how garnet cell volumes change with composition. In the calculation of P_inc the biggest effect of composition is on the determination of V₀. After correction for the shape effects the entrapment isomeke can be calculated with the EoS for that particular composition.

The effects of the uncertainties in the EoS parameters of end-member garnets and the differences between various EoS formulations for the thermal effects on P_trap and therefore estimated depths of entrapment are small for garnets, amounting to no more than 0.2–0.3 GPa in the worst cases, corresponding to less than 10 km uncertainty in depth. On the other hand, if the composition is unknown, then Figure 16 shows that depths cannot be determined from P_inc. Figure 16 also suggests that an alternative approach of calculating the P_trap for each of the endmembers and interpolating between these in proportion to the composition of the garnet might give reasonable estimates of entrapment pressures for eclogitic garnets in the absence of detailed EoS data for the composition of interest. Another alternative is to calculate the entrapment conditions for a mixture of the end-member components, as if they were acting as separate phases. The last caution to note is that at least some eclogitic garnets in diamond have a fluid rim (Nimis et al. 2016), whose influence on the inclusion pressure should be included in calculations, as we discuss below.

If the change in elasticity, the Raman bands, or even the unit-cell volume with composition is not known for the inclusion mineral, an alternative approach is to extract the inclusion after measurement of P_inc and use the free inclusion crystal itself to determine its EoS (Nestola et al. 2019a). This procedure is necessary for inclusions such as spinels whose elastic properties and unit-cell volumes are a complex function of not just composition but also state of cation order.

Non-cubic inclusions. These minerals have elastic properties that vary with direction. When subject to the isotropic strain imposed by the diamond as a result of recovery to room conditions, they should therefore develop anisotropic stresses. The mutual elastic relaxation of the host and the inclusion reduces the anisotropy of the stress, but not to the hydrostatic state (Mazzucchelli et al. 2019); if the diamond and inclusion have behaved perfectly elastically from entrapment, the inclusion should not be under hydrostatic stress when measured in the laboratory. Whether this is reduced to insignificant levels or not depends on the mineral properties and on their crystallographic orientation with respect to their morphology. This can be calculated by FEM calculations for each specific inclusion (Mazzucchelli et al. 2019) but is time consuming.

Therefore, a first step in analyzing non-cubic inclusions is to determine their stress state experimentally. This can be achieved by plotting the measured inclusion unit cell parameters, normalized to the room-pressure values for the mineral, against the fractional volume change of the inclusion normalized in the same way (e.g., Fig. 17). In the case of coesite, a monoclinic mineral, one of the principal axes of compression is the (100) plane normal, so this plane spacing is calculated from the unit-cell parameters and plotted along with b/b₀ and c/c₀. The measured inclusion data for both coesite (Sobolev et al. 2000) and graphite (Glinnemann et al. 2003) fall near to the lines calculated from hydrostatic compression measurements, indicating that these inclusions have a stress field that is close to hydrostatic. The one exception appears to be the lowest-pressure graphite inclusion measured by Glinnemann et al. (2003) which plots on the line of isotropic strain.

As noted above, the positions of Raman bands are dependent upon the strain applied to the crystal (Eqns. 8,9), and not to the mean stress or pressure. Therefore, another way to determine whether the stress state of an inclusion is hydrostatic is to convert the measured Raman shifts into pressures by using the appropriate hydrostatic calibration for the mineral. If different Raman bands give the same pressure, then the stress state in the inclusion is approximately hydrostatic and can be determined from the Raman shifts, and this pressure should agree with that determined by diffraction, as is the case for the coesite inclusions shown in Figure 17 (Sobolev et al. 2000).

When the stress state of a significantly anisotropic inclusion such as coesite or graphite is approximately hydrostatic, this indicates that some other physical process beyond mutual elastic relaxation has modified the inclusion stress state. Such processes may include the presence of fluid in the inclusion or plastic flow in the host diamond to reduce the inclusion pressure. Given that the graphite inclusions measured by Glinnemann et al. (2003) appear to retain a P_inc of 0.1 to 2.6 GPa, consistent with ‘entrapment’ in the thermodynamic stability field of graphite (Fig. 15), the observed hydrostatic stress may be the result of internal redistribution of stress due to internal plastic flow in the inclusion. One might then interpret the as indicating that the graphite formed internally by inversion from diamond, perhaps due to the presence of trapped fluids. However, Glinnemann et al. (2003) also report that the majority of these inclusions show a series of specific CORs with respect to the diamond (Fig. 10b) which they interpret in terms of a coherent interface; such coherency could result in strains that are significantly different from that expected for a spherical inclusion with no bonding to the diamond, in which case the estimates of entrapment conditions with the current simple model of inclusion pressures would be incorrect.

If the stress state of the inclusion is shown to be significantly non-hydrostatic, use of the Strainman program (Angel et al. 2019) allows the strain state of the inclusion to be calculated from the measured positions of the Raman lines. The strains measured in this way, or directly by diffraction, can be converted into stresses by use of the elastic modulus tensor of the inclusion crystal (Nye 1957); see Mazzucchelli et al. (2021) for applications of this approach to inclusions. Of course, minerals become stiffer with increasing pressure, so use of a room-pressure elastic tensor will lead to an under-estimate of the true stress values of the inclusion by a few %. This is insignificant for most inclusions, except those that retain pressures in excess of 1 GPa. For quartz inclusions in garnets Bonazzi et al. (2019) and Gonzalez et al. (2019) have demonstrated that the stress and strain state is non-hydrostatic, but that a good estimate of entrapment pressures can be obtained by using the mean normal stress as P_inc in subsequent calculations. A similar approach should work for non-cubic inclusions in diamonds.

For a number of low-symmetry inclusion phases neither the elastic tensors are known, nor have the mode Grüneisen tensors been determined, and therefore the detailed stress state in the inclusion cannot be determined. An alternative approach in these cases is to use DFT simulations to determine which Raman line or lines is least affected by non-hydrostatic stresses (Anzolini et al. 2018; Nestola et al. 2018a). The position of this line in the Raman spectrum of the inclusion is then measured, and the shift from the position of the same line in a free crystal is then converted to a pressure using the measured hydrostatic calibration. For a kyanite inclusion in a diamond from the Voorspoed kimberlite, without any detectable fluid present in the inclusion, this method gave a calculated entrapment pressure for the measured FTIR N-aggregation residence temperature that overlaps with P–T estimates obtained by traditional chemical geobarometry for diamonds from other kimberlites from the Kaapvaal craton (Nestola et al. 2018a). The calculated entrapment conditions for a CaSiO₃ inclusion (breyite) determined in this way are in the thermodynamic stability field of the phase (Anzolini et al. 2018).

The elastic anisotropy of olivine is significantly less than coesite or graphite and available measurements (e.g., Nestola et al. 2011) do not clearly distinguish whether the stress state of olivine inclusions is hydrostatic or not. Given the low elastic anisotropy and that the mean stress in olivine inclusions is typically less than 1 GPa, no significant error is introduced by determining the P_inc from the unit-cell volume of the inclusion as measured by diffraction. Figure 18 shows the measured volumes of 13 olivine inclusions from the Udachnaya kimberlite, plotted against their composition as determined by in situ X-ray data collection and structure refinement, as described in a previous section of this chapter. In order to determine the P_inc it is also necessary to know the unit-cell volume that the inclusion would have at room pressure, V₀. The latter varies with composition. If the inclusion cannot be freed from the diamond the V₀ of the inclusion can be calculated from the composition determined by in situ structure refinement, and the known variation of V₀ of olivines with Mg/Fe ratio (Schwab and Küstner 1977). Great care has to be taken in normalizing the unit-cell measurements because the exact values depend on diffractometer alignment, calibration and the peak-centering algorithms used (Angel et al. 2000); the cell parameters of the same crystal measured on different diffractometers will often differ by many times the estimated standard uncertainties. The simple solution is to measure in air a crystal or series of crystals of known composition with exactly the same protocols on the same diffractometer on which the inclusions are measured, at a similar time. Tests with data sets from other mantle olivine crystals collected in air show that this procedure yields the correct V₀ with uncertainties dominated by the uncertainty in composition from structure refinement (Fig. 18).

The equilibrium phase boundary between graphite and diamond (Day 2012) at mantle temperatures lies between the entrapment isomekes for olivine in diamond that would give rise to P_inc = 0.5 and 0.6 GPa (Fig. 19). Allowing for the small effects of shape on relaxation, a P_inc of approximately 0.5 GPa is the minimum pressure that should be found in any olivine inclusion. This corresponds to a volume compression of V/V₀ = 0.996, a boundary that is also shown in Figure 18. This type of graph then becomes a useful diagnostic tool because it immediately shows which inclusions have volumes, and hence P_inc values, consistent with entrapment during growth of the diamond within its thermodynamic stability field. At these minimum depths and pressures, uncertainties in the EoS of olivine (Angel et al. 2018) result in an uncertainty in P_inc of less than 0.03 GPa, corresponding to an uncertainty in unit-cell volume of 0.06 ų, about the size of the error bars in volume in Figure 18. It is therefore clear from this figure that these olivine inclusions have had their inclusion stress state modified by some other process that is not considered in the simple elastic model. The fact that the measured volumes are larger than the maximum expected from the elastic analysis is consistent with partial, but not complete, pressure release which could occur either by brittle fracture of the diamond at low temperatures or by plastic flow at higher temperatures. Certainly, a number of these inclusions have very fine fractures of a few µm width leading from them. If a crack of thickness t were planar and extended around the entire equatorial plane of a spherical inclusion of radius r it would create a fractional volume increase in the inclusion of 3t / 4r. A planar fracture 1 µm thick around an inclusion of radius 50 µm would increase the volume by 1.5%, which is more than enough to account for all of the inclusion volumes shown in Figure 18. However, microscopy and X-ray absorption tomography (e.g., Rustioni et al. 2015; Mazzucchelli et al. 2016; Nimis et al. 2016) show that cracks do not surround the inclusions in this way and therefore they can only contribute a much smaller increase in inclusion volume that is insufficient by orders of magnitude to account for the observed pressure release indicated by Figure 18.

The influence of fluids. It was recently discovered that olivine inclusions in gem-quality diamonds are surrounded by a thin film of hydrous silicic fluid that is no more than 1.5 µm thick (Nimis et al. 2016). Such films seem to be almost ubiquitous around silicate and oxide inclusions in lithospheric diamonds (e.g., Nestola et al. 2018a; Nimis et al. 2019) and can modify the stress state of the inclusion through several processes. First, because they are a fluid, they cannot support shear stresses which should ensure that the stress state of the solid part of the inclusion is hydrostatic. This may explain the hydrostatic stress state of some of the inclusions shown in Figure 17, and it allows the model of isotropic elasticity to be used to calculate entrapment conditions, provided that allowance is made for the presence and equation of state of the fluid. The latter is unknown because the exact composition of the fluid is not known. If the properties of the fluid can be modelled with those of water (Zhang and Duan 2005), then the final pressure of a coesite inclusion is decreased compared to a dry inclusion; a dry inclusion would exhibit P_inc ~3 GPa if it had been trapped at 4.2 GPa and 1000 °C, and 10 vol% water would decrease P_inc to ~2.1 GPa. This means that if the presence of water in a coesite inclusion is ignored the entrapment pressures will be under-estimated.

By contrast, in olivine the expected final inclusion pressures are increased by the presence of water for entrapment near to the lower stability limit of diamond, and the effect of water decreases with increasing depth of entrapment (Fig. 20). This sounds like a surprising and counter-intuitive result, because fluids are much more compressible than solid minerals, the bulk modulus of liquid water being less than 2 GPa at room conditions (Zhang and Duan 2005). However, their thermal expansion coefficients are also an order of magnitude higher than solids. As a consequence, the isomekes of water with diamond have similar slopes to those of olivine with diamond (but not coesite and diamond; Fig. 15) and hence the presence of pure water in an inclusion in olivine has a small effect on the inclusion pressures (Fig. 20).

Thus, the presence of water alone therefore cannot explain the low inclusion pressures that are uniformly observed in olivine inclusions, which would be expected to exceed 1 GPa if they were entrapped significantly deeper than the depths of the graphite–diamond phase boundary (Figs. 15, 19, 20). But, the presence of apparently amorphous solid on the diamond surface and in contact with the fluid rim in these inclusions (Nimis et al. 2016) suggests that the fluid may have precipitated solids that were originally dissolved in the fluid at the time of entrapment, just as brucite is found as a precipitate in some super-deep diamonds (Palot et al. 2016). If the precipitation is associated with an overall reduction in the volume of the fluid plus precipitate it will contribute to a further decrease in inclusion pressure. The same reasoning applies to the precipitation of additional olivine onto the olivine, or carbon onto the diamond during ascent of the diamond. The evaluation of the magnitude of these effects awaits the full characterization of the chemical composition of the various components of the inclusion rims, and thermodynamic modelling.

The presence of fluid also greatly magnifies the effects of cracks around the inclusion compared to cracks around purely solid inclusions. The formation of the crack allows the fluid to flow out from the rim and into the crack, effectively increasing the volume of the fluid rim by the volume of the crack. The bulk modulus of fluid water at shallow lithospheric depths is of the order of 10 GPa (Zhang and Duan 2005). Therefore, the drop in pressure of an inclusion fluid as the result of the opening of a crack of just 1/3 of the volume of the fluid rim would be enough to drop the inclusion pressure by ~3 GPa. This is enough to reduce to zero the P_inc of any olivine trapped at any depth in the lithosphere. Smaller crack volumes are therefore sufficient to explain the range of pressures and unit-cell volumes measured for olivine inclusions (Fig. 18 and Nestola et al. 2011) and why olivine P_inc values are often much less than expected. The conclusion is that without the characterization of the fluid in detail and the volumes of the cracks and fluid rim, olivine inclusion pressures should not be used to infer growth conditions of the diamonds. Both the sulfide-carbide and silicate inclusions in CLIPPIR and blue diamonds (Smith et al. 2016; Smith et al. 2018) also contain both CH₄ and H₂ as fluid phases, and calculations similar to those that we have outlined here will have to be made in order to calculate entrapment conditions from the observed inclusion stress states.

Phase transformations. Many mineral phases that are stable in the deep Earth are thermodynamically unstable at lower pressures and, in the absence of constraints, transform to lower density polymorphs or unmix into a less dense assemblage of phases. For example, the density changes associated with the transformations of ringwoodite to wadsleyite and to olivine are responsible for the well-known seismic discontinuities that define the Earth’s transition zone. In the mantle, such phase changes can occur because the materials are not constrained. By contrast, an inclusion phase is constrained by the surrounding diamond and is not free to expand. As a consequence, the state of thermodynamic equilibrium is changed to one in which the boundary conditions are no longer P and T, but the T and the volume imposed by the diamond (Ferrero and Angel 2018). This concept is also the origin of inclusion pressures that we measure as P_inc. If the inclusion pressure drops into the stability field of a lower pressure polymorph and the temperature is sufficiently high, whether or not the transformation takes place depends on whether the volume change due to the transformation can be accommodated by the diamond expansion on exhumation, unless the necessary volume is provided by processes such as plastic flow or brittle fracturing of the diamond host.

If a mineral transformation is accompanied by a large density change, then the volume is not available to accommodate the volume change due to the transition, and the high-pressure phase can be preserved. This is why it is possible to recover single-phase inclusions of phases such as jeffbenite, a high-pressure polymorph of garnet (Armstrong and Walter 2012; Nestola et al. 2016b) and ringwoodite (Pearson et al. 2014), as well as calcium silicate perovskite with minor exsolution (Nestola et al. 2018b). The reverse argument has been made to conclude that CaSiO₃ inclusions with the walstromite structure (the mineral breyite) cannot be inverted from CaSiO₃ perovskite, because the 30% volume change could not be accommodated by the diamond host (Anzolini et al. 2016). Indeed, the measured P_inc of 4.3 GPa is consistent with entrapment of the CaSiO₃ inclusion within the stability field of breyite (Anzolini et al. 2018; Woodland et al. 2020). The P_inc is also higher than the lower-pressure limit to breyite stability, which explains why it did not invert to its lower-pressure wollastonite polymorph. Fractured inclusions of breyite exhibit much lower pressures (Anzolini et al. 2016) which indicates that the fracturing occurred at temperatures that were too low to allow breyite to invert to wollastonite.

Transformations that involve smaller volume changes can be accommodated by the expansion of the diamond as the pressure of the inclusion drops below the equilibrium phase boundary. Then, the high-pressure phase will start to transform to the low-pressure polymorph. However, because the low-pressure polymorph is less dense, the pressure in the inclusion will increase as the transformation proceeds, until the pressure returns to that of the phase boundary. At this point the transformation will stop. Further exhumation may lead to further expansion of the inclusion volume, allowing further transformation to occur, so the pressure will be maintained at approximately the pressure of the phase boundary until the transformation of the inclusion is complete. This is the concept of pressure buffering. The various possibilities for the subsequent evolution of the inclusion pressure, which depend on whether the transformation proceeds to completion, are discussed for the inversion of coesite to quartz in garnet host crystals by Ferrero and Angel (2018). These principles will also apply to inclusions in diamonds. For coesite inclusions in diamonds, the evidence of the measurements of Sobolev et al. (2000) and Smith et al. (2015), which are in agreement with EoS calculations (Fig. 15) is that sealed and isolated coesite inclusions can remain at pressures in the stability field of coesite and no transformation to quartz occurs.

Mineral inclusions that are single-phase solid solutions at the time of entrapment often exsolve other phases due to the reduced mutual solubility of chemical components at lower temperatures and pressures. Examples include clinopyroxene exsolution from majoritic garnets (Harte and Cayzer 2007), spinel exsolution from magnesiowustite (Wirth et al. 2014; Palot et al. 2016) and the exsolution of titanite perovskite from silicate perovskite (Nestola et al. 2018b). As for simple isochemical phase transformations, the density changes associated with these reactions mean that it is possible for the pressure in the inclusion to be buffered to a higher value than expected for a single phase. The interpretation of the residual pressures of such multiphase inclusions is further complicated by the need to know the volume fractions of all of the phases present, their individual equations of state and, in the case of coherent exsolution, the requirement to consider the effects of the coherency between the phases on the stress and strain evolution of the composite inclusion. While calculations of coherency strains have been made for crustal minerals at ambient pressure (e.g., Robin 1974; Williame and Brown 1974), we are not aware of any being applied to mantle phases trapped in diamonds.

Plastic deformation. Diamond can accommodate stresses by plastic deformation generated by the motion of dislocations. This is a thermally activated process whose rate also depends on the applied stress. Plastic deformation is not significant in diamond below 1000 °C (Weidner et al. 1994) but at 1700 °C it is quite rapid even when the applied deviatoric stresses are low (Howell et al. 2012b). Therefore, inclusion pressures can be reduced by plastic flow driven by the pressure difference between P_inc and the external pressure during diamond ascent at high temperatures; the driving force is zero at the time of entrapment because the inclusion is trapped at the external pressure but increases while the diamond travels towards the crust. In general, lower temperatures and rapid exhumation rates will allow less plastic deformation to occur. As a consequence, there is little evidence for significant plastic flow around inclusions in lithospheric diamonds (Agrosi et al. 2016), and consideration of plastic flow is not needed to interpret the inclusion pressures measured in these diamonds. However, pervasive dislocation networks are present in sub-lithospheric diamonds (Cayzer et al. 2008; Agrosi et al. 2017; Smith et al. 2018). It is clear that this plastic deformation is the mechanism by which these diamonds can accommodate the volume changes associated with the inversion of high-pressure phases to the lower-pressure polymorphs sometimes found in these diamonds (Cayzer et al. 2008). The effect on inclusion pressures can, in principle, be calculated either by explicit finite element modelling of the inclusion (Anzolini et al. 2019) or by a 1-D visco-elasto-plastic mechanical model which enables the effects of different exhumation rates on the inclusion pressure to be explored (Zhong et al. 2018, 2020) and, in principle, constraints to be placed on the speed of ascent of the kimberlite hosting the diamond.

The history of the study of inclusions in diamonds shows that it is extremely important to characterize mineral inclusions in situ so as to preserve them for measurement by techniques that are either not yet capable of the sensitivity or resolution required, or perhaps not yet even developed. An obvious example is the discovery of new mineral phases as inclusions in diamond which has been made possible by the far greater X-ray intensities available at the synchrotron compared to the laboratory (e.g., Meyer et al. 2019). The scarcity of examples discussed in this chapter shows how few data have been collected in situ from inclusions in diamonds, compared to the thousands of ex situ measurements. Yet these examples have shown the new discoveries (e.g., the fluids found by Nimis et al. 2016) that can be made, and the far-reaching questions that can be answered, by in situ measurements with equipment available in any crystallography laboratory. There are thus many opportunities to use modern analytical methods and computational power to answer very specific questions about inclusions in diamonds and thus how they were formed and encapsulated, and at what conditions in the Earth.

The techniques that we have described here in some detail for certain mineral inclusions in diamond can be readily applied to other mineral inclusions in diamonds. In particular, a large number of orientation measurements will enable the CORs of inclusions with respect to diamonds to be unambiguously identified and tests performed to determine whether there is a correlation between COR types and the type of host rock or other aspects of the diamond growth environment. Such measurements should continue to be applied to inclusions such as olivine, garnet or graphite, for which there is insufficient data to determine with statistical significance whether multiple CORs exist as sub-sets of the data; for example, it is not clear whether the graphite inclusions shown in Figure 10b represent two specific CORs, or a single rotational COR. Many more determinations of graphite orientations in diamond would answer this simple question and help determine how such inclusions, and their host diamond, were formed. Similarly, there is only published data for ferropericlase inclusions in two diamonds; many more measurements are required to determine whether the specific COR observed is a general phenomenon and whether it has genetic significance. Orientation measurements of inclusion clusters can be combined with imaging techniques such as X-ray tomography to provide clear insights into the mechanism by which the inclusions were trapped and, therefore, how the diamond grew (Nimis et al. 2019). Similar insights may well be obtained by further measurements of inclusions not only in situ in the diamond, but while the diamond is still in contact with its current host rock (Nestola et al. 2017).

The major element chemistry of an inclusion can be determined by in situ X-ray diffraction and spectroscopy. Doing so preserves the inclusion for future research while providing useful data. For inclusions in sub-lithospheric diamonds that may represent the only known terrestrial example of a mineral, in situ measurements avoid the risk of the destruction of unique samples either through phase transformation or amorphization, or by simple loss of a sample much smaller than its host diamond. Thus the discovery of the ringwoodite as a single inclusion in diamond (Pearson et al. 2014) was proven by a combination of in situ single-crystal diffraction and Raman spectroscopy. The in situ measurement of its infra-red absorption spectrum was the first ever proof that the mantle transition zone is significantly hydrated, all obtained without risk of sample loss.

Despite studies published over the last half-century (e.g., Harris et al. 1970) the use of inclusion stress states to determine depths of diamond growth is only now becoming feasible through concurrent development of theory with the improvement in the precision of both the in situ measurements that we have described and the elastic properties of minerals. The presence of fluid around many lithospheric inclusions (Nimis et al. 2016) both complicates and simplifies the analysis of their stress state. It complicates analysis because it is now necessary to determine the phase equilibria of the trapped fluid and the associated amorphous rims in order to know how the total volume of the rim varies with P and T, not only from their respective EoS but also as a result of reactions within the fluid and the precipitation of the amorphous material (Palot et al. 2016). If that is done, then the problem of inclusion stress is simplified in two ways. The mineral component of the inclusion will be under hydrostatic stress and thus the shifts of its Raman lines will directly yield the inclusion pressure via calibrations of Raman shifts with pressure. The elasticity of the composite fluid + mineral inclusion can be treated as isotropic and the entrapment conditions calculated from P_inc using the Reuss bound on the elastic properties of the mineral plus those of the rim (Musiyachenko et al. 2021).

The current models for determining inclusion entrapment conditions from the inclusion stress assume free-slip boundary conditions between the host and the inclusion. Therefore, if fluid is absent then orientation measurements must first be performed to determine whether or not the inclusions have an exact COR to the diamond. If they do, this may indicate epitaxial bonding between the diamond host and the inclusion which may modify the stress state significantly through coherency, as is probably the case for graphite (Glinnemann et al. 2003) and ferropericlase (Nimis et al. 2018) inclusions. Interpretation of the inclusion stress state then requires explicit calculations by FEM of the contribution from the coherency strains and how they vary as the external P and T evolve during diamond ascent. In the absence of coherency the mean stress in the inclusion can be treated as a pressure in the isotropic model to obtain a good estimate of the entrapment isomeke (Bonazzi et al. 2019; Gonzalez et al. 2019). Beyond that approach, the anisotropy of the elastic properties of an inclusion should allow the isomekes of different directions in a single inclusion crystal to be used to determine the entrapment conditions uniquely, as recently demonstrated for quartz inclusions in garnets (Alvaro et al. 2020).

We have not discussed the characterization of the strain in the diamond created by the pressure difference between the inclusion and the external ambient pressure on the diamond because its quantitative analysis is far more complicated and at the limits of current capabilities in mineral physics. The difficulty arises from the fact that the symmetry of the diamond in the vicinity of the inclusion is broken by the combination of radial compressional strains and tensional tangential strains arising from elastic relaxation (Fig. 13) and it is therefore no longer cubic. The reduction in the symmetry of the diamond means that its physical tensor properties change, including the optical properties which lead to the optical birefringence haloes around inclusions. Similarly, the values of the mode Grüneisen tensor components may be expected to change from their values in cubic diamond. Once these have been determined by DFT simulation and confirmed to agree with experimental measurements (Grimsditch et al. 1978), the shifts of the diamond Raman modes in the host diamond have the potential to be used to infer the inclusion pressure (Nasdala et al. 2005). Similarly, it is possible that the intensities of the Raman line of diamond could be used to map the strain field around inclusions, similar to the method applied to garnets that host inclusions (Campomenosi et al. 2020a).

The general conclusion to be drawn from this review is that rapidly-developing instrumental technologies along with new concepts from mineral and crystal physics and the application of advanced computational methods such as DFT and FEM, can be used in combination to provide new ways to characterize inclusions in diamonds and yield new insights into how and where they nucleate and grow within the Earth.

Our own research in this field has been supported by funding from the European Research Council under the European Union’s Horizon 2020 research and innovation program through grant agreements 307322 to Fabrizio Nestola and 714936 to Matteo Alvaro, and by the SIRMIUR Grant MILE DEEp (no. RBSI140351) to Matteo Alvaro. We would like to thank Mattia Mazzucchelli for providing calculations on shape correction factors, Paolo Nimis for recalculating geotherms, and both for discussions and comments. For more than a decade Mathias Meyer (Rigaku Oxford Diffraction) has given us much valuable advice about collecting and processing X-ray data and has willingly developed new features in the Crysalis^TM software for our research programs. We thank Michelle Wenz, Martin Kunz, Daria Pasqual, and the editors Karen Smit and Graham Pearson for their very helpful and thorough reviews and comments.