Experimental variable effects on laser heating of inclusions during Raman spectroscopic analysis

Abstract Raman spectroscopy for fluid, melt, and mineral inclusions provides direct insight into the physicochemical conditions of the environment surrounding the host mineral at the time of trapping. However, the obtained Raman spectral characteristics such as peak position are modified because of local temperature enhancement of the inclusions by the excitation laser, which might engender systematic errors and incorrect conclusions if the effect is not corrected. Despite the potentially non-negligible effects of laser heating, the laser heating coefficient (B) (°C/mW) of inclusions has remained unsolved. For this study, we found B from experiments and heat transport simulation to evaluate how various parameters such as experimental conditions, mineral properties, and inclusion geometry affect B of inclusions. To assess the parameters influencing laser heating, we measured B of a total of 19 CO2-rich fluid inclusions hosted in olivine, orthopyroxene, clinopyroxene, spinel, and quartz. Our results revealed that the measured B of fluid inclusions in spinel is highest (approx. 6 °C/mW) and that of quartz is lowest (approx. 1 × 10−2 °C/mW), consistent with earlier inferences. Our simulation results show that the absorption coefficient of the host mineral is correlated linearly with B. It is the most influential parameter when the absorption coefficient of the host mineral (αh) is larger than that of an inclusion (αinc). Furthermore, although our results indicate that both the inclusion size and depth have little effect on B if αh > αinc, the thickness and radius of the host mineral slightly influence B. These results suggest that the choice of inclusion size and depth to be analyzed in a given sample do not cause any systematic error in the Raman data because of laser heating, but the host radius and thickness, which can be adjusted to some degree at the time of sample preparation, can cause systematic errors between samples. Our results demonstrate that, even with laser power of 10 mW, which is typical for inclusion analysis, the inclusion temperature rises to tens or hundreds of degrees during the analysis, depending especially on the host mineral geometry and optical properties. Therefore, correction of the heating effects will be necessary to obtain reliable data from Raman spectroscopic analysis of inclusions. This paper presents some correction methods for non-negligible effects of laser heating.


Introduction
Raman spectroscopy is now well established as a versatile tool for inclusion research because of its non-invasive analysis and high spatial resolution. This method has been applied to various fields of inclusion studies: ascertaining the depth provenance of rocks from mineral inclusions (Kohn, 2014), estimation of peak metamorphic temperatures from graphite inclusions (Beyssac et al., 2002), determination of the preeruptive volatile contents of magmas from melt inclusions (Moore et al., 2015), and characterization of paleo crust-mantle fluid activities from fluid inclusions (Frezzotti et al., 2012). To extract the physicochemical properties of unknown inclusions from the measured Raman spectra, experimental calibrations related to spectral features to the P-V-T-X properties of fluids or the strains and crystal orientations of minerals have been the subject of many studies (e.g., Frezzotti et al., 2012;Angel et al., 2019). Furthermore, the combined technique of microthermometry and Raman spectroscopy enables accurate estimation of the temperatures of phase transition in situ, even for phases that are difficult to assess under optical observation (e.g., Bakker, 2004). Typically, to obtain Raman spectra of an inclusion having a sufficient signalto-noise ratio for quantitative analysis, laser power of tens to hundreds of milliwatts is necessary, thereby necessitating focal intensities as high as several gigawatts per square centimeter. Even for analysis of transparent geomaterials, such high excitation energy can lead to nonnegligible temperature increments at the focal position because of light absorption by host minerals or inclusions themselves. Such undesired heating at the focus will change Raman spectral features and might engender misinterpretation of the data. However, this important experimental parameter, temperature, has been difficult to control and estimate during Raman spectroscopic analysis because local temperatures behave intricately depending on the sample geometry, experimental conditions, and sample optical and heat transport properties. To interpret the inclusion data obtained by Raman spectroscopy correctly, one must ascertain the temperature of inclusion under the analysis and clarify how the temperature is governed as functions of experimental, geometrical, and optical parameters.
Although the thermal damage caused by the excitation laser to samples has been a concern especially in the field of bioengineering (e. g., Català et al., 2017), no comprehensive work assessing laser heating effects of small objects hosted in a solid phase, i.e., inclusions, during Raman spectroscopic analysis has been reported to date. Fall et al. (2011) observed the phase transition of CO 2 fluid after laser irradiation. In their experiments, the CO 2 -bearing inclusions, for which homogenization temperatures are slightly higher than the ambient temperature, are heated by an excitation laser, being not homogenized during the experiment. Therefore, the effect of laser heating might be negligible (less than <1 • C), at least under their experimental conditions and when using minerals with very small visible light absorption coefficients, such as quartz. However, the absolute value of the temperature rises at a certain laser power (laser heating coefficient (B); • C/mW) has persisted as an open problem for inclusion in representative rock-forming minerals, including quartz.
In the field of Earth science, laser-heating effects on the sample are evaluated especially with respect to Raman spectroscopy of a carbonaceous material (RSCM) thermometer that is used widely for estimating peak metamorphic temperatures. Kagi et al. (1994) reported that an appreciable local temperature increase occurs because of absorption of the incident laser by graphite. According to data reported by Kagi et al. (1994), B is approx. 1 × 10 2 • C/mW. In addition, Umeda and Enami (2014) reported that laser power higher than 2 mW would change apparent spectral properties of CM and would provide apparently high peak metamorphic temperature by the RSCM thermometer (see review by Henry et al., 2019). Similar warming effects by transmitted light have been reported for near-infrared microthermometry for fluid inclusions in opaque to translucent minerals (e.g., Casanova et al., 2018). Therefore, depending on the experimental, geometrical, and optical parameters, the inclusions are heated by the laser during Raman spectroscopic analysis. Their local heating alters the Raman spectral properties. Nevertheless, little is known about which parameters strongly affect laser heating.
To ascertain a correct physical picture of laser heating in inclusion, we perform heat transport simulations of an inclusion hosted by a solid phase (mineral) using finite element method, taking special care of the effects of spherical aberration at the air-sample boundary and the entire three-dimensional spatial profile of the focused beam. Based on the model, the effects of various parameters were evaluated. We also report the experimentally determined value of B of CO 2 -rich fluid inclusions hosted by representative mantle minerals and quartz. Finally, by combining experimental and simulation results, we discuss the effects of the respective parameters on the laser heating of inclusions. Goals of this paper are to ascertain the most influential parameters on the laser heating of inclusions, and to propose favorable analytical procedures under which results obtained from Raman spectroscopy are reliable.

Literature review
Other research fields for which B has already been determined provide some examples. A field showing particular progress is optical tweezers for microscopic particle manipulation (e.g., Liu et al., 1995). Earlier studies measured temperatures surrounding optically trapped microspheres in water or glycerol during near infrared irradiation. Reportedly, temperature increases are 0.008-0.099 • C/mW. In most cases, determination of local temperature in an optical trap is based on measurement of the viscosity of the liquid (e.g., water) surrounding an optically trapped bead (e.g., silica) rather than a bead itself. By contrast, our target is a fluid inclusion (liquid CO 2 ) itself, located at the focal point of the laser, which is surrounded by a solid phase (host mineral).
According to earlier studies, one parameter that strongly influences laser heating is the absorption coefficient of the host phase, not of the trapped particle (Liu et al., 1995;Celliers and Conia, 2000;Peterman et al., 2003;Mao et al., 2005;Català et al., 2017). Because the absorption coefficient depends on the wavelength, it also depends on the laser wavelength (Haro-González et al., 2013). Regarding the bead size, two opposing results have been reported: Peterman et al. (2003) asserted that size has no effect on B. Also, Català et al. (2017) reported that larger particles suppress laser heating. The distance between the bead and the boundary of solvent and surrounding material, i.e. glass, also affects B. In fact, B becomes small only when the bead is close to the interface with the glass (Peterman et al., 2003;Català et al., 2017). Based on the model proposed by Mao et al. (2005), B is expected to show logarithmic growth with increasing sample thickness. In addition, B is correlated linearly with the reciprocal of the thermal conductivity of the liquid (Celliers and Conia, 2000;Peterman et al., 2003;Mao et al., 2005;Català et al., 2017). The specific B and measurement methods used for earlier studies are presented in Table S2 of the report by Català et al. (2017).

MODEL: Three-dimensional heat transport simulation using finite element method
To clarify the temperature distribution within and around an inclusion when illuminated by excitation laser, we performed heat-transfer simulation by finite element method using a software package: COM-SOL Multiphysics. Fig. 1 depicts a laser beam transmitted through the objective lens towards the inclusion. We assumed symmetric disc shape for the host mineral. The disc has H h thickness and R h radius. A spherical inclusion with radius of r inc exists at z inc below the top surface of the disc. The actual focal position (= inclusion depth) is given as z inc = (n h /n a )z 1 , where z 1 is the nominal focal position below the sample surface and where n a and n h respectively denote the refractive indices of the air and host mineral.
Temperature is in general terms governed by the heat equation where k h (m 2 /s) stands for the thermal diffusivity, K h (W/mK) expresses the thermal conductivity, and α h denotes the absorption coefficient of host mineral (1/m). Also, I ill (W/m 2 ) represents the spatial distribution of the focused intensity inside the sample as (e.g., Maruyama and Kanematsu, 2011) where q (m) represents the radial distance from the optical axis of an illuminated point, L (m) denotes the depth of an illuminated point, and w 0 (m) corresponds to the waist-radius of the laser-beam, where the laser-beam diameter converges to a minimum. The value of w 0 was estimated directly from the Raman instruments used for the measurements. The bitmap image of the laser spot reflected from the silicon wafer surface is used. The estimated w 0 was 1.46 μm for LU Plan ×50 (used at Hokkaido University) and 1.0 μm for TU Plan ×50 (used at Japan Agency for Marine-Earth Science and Technology (JAMSTEC)). Detailed procedures for estimating w 0 are presented in online Supplementary material S1. Also, I 0 (W/m 2 ) and w(L, z 1 ) respectively represent the laser intensity at the center of the beam waist and the beam-waist radius at the axial position, given as where P ill represents the laser power at the sample surface (W) and z 0 = πw 0 2 /λ (m) is the Rayleigh range. Also, λ denotes the excitation laser wavelength. In an earlier model, the boundary condition which determines heat flux at the interface was the Dirichlet condition (i.e., T = const. at boundaries), but in this study, in addition to the Dirichlet condition, we also performed simulations under convective heat flux boundary condition that is defined by h(T ext − T) at the glass-sample and air-sample interfaces. At the air-sample boundary, the convective heat transfer coefficient (h a/s ) was fixed to be 5 W/m 2 K at z = 0 plane and q = R h curved surface, assuming heat transfer by free convection (Whitelaw, 1997). In case of the boundary between the sample and the plate on which the sample is placed (z = H h ), it is difficult to settle the style of heat transport because it is affected by many parameters such as surface conditions, optical and heat transport properties of the plate, and the temperature distribution below the sample. Therefore, we considered the heat transport coefficient at the glass-sample interface (h g/s ) as a variable. Both the ambient (T ext ) and initial sample temperature (T sample ) are equal. Of the parameters throughout Eqs. (1) to (4), α h , k h , K h , and n h are optical and physical properties unique to the host mineral and are shown in Table 1. The influence of reflection and crystal orientation dependence of α h , k h , K h , and n h on the laser heating effect was not considered here.
No report of the relevant literature describes a study of the absorption coefficient of CO 2 at the specific conditions of our simulation, i.e. 532 nm and approx. 7 MPa. Thompson et al. (1963) estimated that absorption coefficient of CO 2 at 198 nm is 10 − 5 cm − 1 . According to Hansen (1997), absorption coefficient of solid CO 2 at 532 nm is about 10 − 4 cm − 1 . In both studies, the absorption coefficient decreases as the wavelength increases from approx. 180 nm. The absorption coefficient of CO 2 at 532 nm is so small compared to the host mineral (for the smallest case, e.g., 0.02 cm − 1 for quartz) that it was regarded as negligible in the present model. Therefore, we assumed it to be 0 (see section 7.1 for details). Thermal diffusivity of CO 2 is 0.02 mm 2 /s at the conditions of 0.7 g/cm 3 and 25 • C based on the NIST Chemistry WebBook (Linstrom and Mallard, 2020). The simulated inclusion temperature and B are an average over an inclusion domain in 3D models.
All the parameters used for this study are presented in Table 2. Hereinafter, we use an abbreviation for each parameter. As described in section 2, earlier studies have experimentally or theoretically examined the effects of λ, w 0 , α inc , z inc , r inc , K inc , K h , α h , H h , and R h on B. In addition to these parameters, we investigated the effects of n h , which can change the temperature distribution around the inclusion, k h , which affects the time to reach the steady state, and h a/s and h g/s , which determine the heat flow rate at the boundary, on B.

Samples
We selected a total of 19 CO 2 -rich fluid inclusions hosted in mineral separates of five species; olivine, orthopyroxene, clinopyroxene, spinel, and colorless quartz (Fig. 2). Fluid inclusions in olivine, orthopyroxene, clinopyroxene, and spinel are derived from a mantle xenolith designated as "En2A" from Ennokentiev, Sikhote-Alin, Far Eastern Russia. Petrological descriptions were presented by Yamamoto et al. (2012). Quartzhosted fluid inclusions are derived from pegmatites of Kadugannawa Complex from Owala-Kaikawala area, Matale, Sri Lanka. Seven mineral grains were selected and made into doubly polished chips. The fluid inclusion compositions are investigated using Raman spectroscopic analysis; volatile species other than CO 2 and H 2 O are not detected. Although no H 2 O liquid phase was visible under microscopic observation, trace amounts of H 2 O were confirmed in Raman spectra from some Fig. 1. Schematic of the simulation situation considered in the model. Refractive indices of the host mineral and air are presented respectively as n h and n a . The surface (i.e., air-host mineral boundary) of the host mineral is defined as the z = 0 plane. The focus position of a transmitted beam shifts from an intended position at z = z 1 to z = (n h /n a )z 1 . A blue sphere represents an inclusion with radius of r inc and depth of z inc . q is the radial distance from the optical axis of an illuminated point and L denotes the depth of an illuminated point. The host mineral shape is assumed to be a disc with radius R h and height H h . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)  Span and Wagner (1996), as shown in Table 3. The respective sizes and depths of the fluid inclusions were measured using a 2000× objective (VH-ZST; Keyence Co.) with a digital microscope (VHX-5000; Keyence Co.). The projected area of fluid inclusions was measured automatically by identifying the outline of a fluid inclusion using the brightness contrast. The inclusion radius (r inc ) was defined by the equivalent circle radius of the obtained projected area. We also measured the equivalent circle radius of projected area of host mineral (R h ). We calculated the depth of the inclusion from the distance between the mechanical position for focus at the top of the host mineral and at a position where we obtained a sharp image of a fluid inclusion. The actual depth (z inc ) is calculated approximately by the simple expression of z inc = refractive index of host mineral (n h ) × measured depth. As shown in Table 1, the values of n h of olivine, orthopyroxene, clinopyroxene, spinel, and quartz are, respectively, 1.65, 1.66, 1.65, 1.72, and 1.55 (Lucey, 1998;Palik, 1998;Ghosh, 1999). In addition to z inc , the host mineral thickness (H h ), which is defined by the distance between the mechanical position for focus at the surface of the host mineral and top of the slide glass, is shown in Table 3.

Overview of the hot band to Fermi diad intensity ratio thermometer
The CO 2 fluid temperature can be estimated using Raman spectra: a so-called hot band thermometer was developed using a high pressure optical cell (HPOC) (Rosso and Bodnar, 1995;Arakawa et al., 2008;Hagiwara et al., 2018). The experimental calibrations of earlier studies do not cover the P-T conditions of our interest. Therefore, we expanded the applicable P-T range using newly measured Raman spectra of CO 2 from 23 to 200 • C in approximately 20 • C increments under 10 temperature conditions (Fig. S1). Specific temperature, pressure, and density conditions are presented in Table S1. Details of the experimental apparatuses and procedures are provided in online Supplementary materials S1.   (5) is useful at 0.6 < ρ < 1.2 g/ cm 3 and 23 < T < 200 • C. The maximum absolute value of residual errors for all 131 measurements was 9.5 • C. The average of the absolute value of residual errors was 2.0 • C. Applicability of this thermometer for natural fluid inclusions was confirmed from natural fluid inclusions in quartz, which are shown as square (temperature increasing cycle) and triangle (decreasing cycle) symbols in Fig. 3. This test experiments were performed with a laser power of 7.7 mW because our experiments show that the hot band to Fermi diad intensity ratio of a CO 2 inclusion in quartz does not change for laser powers up to at least 15 mW ( Fig. S2 and Table S2).

Raman setup
To determine B, Raman spectra of fluid inclusions were measured using P ill ranging from 4.2 to 14.2 mW along both laser increasing and decreasing cycles. The values of P ill during the experiments were measured using a laser power meter (PD300-ROHS; Ophir Optronics Ltd.). The uncertainty of the measured P ill is ±3%. During the measurement, T sample was controlled by the room temperature, which was controlled by an air conditioner. The sample was placed on the Linkam stage. The actual T sample was monitored by the thermocouple of the Calibrated hot bands to the Fermi diad intensity ratio defined as ] were converted to temperature using Eq. (5). Table S3 presents the T sample and average fluid temperatures of five measurements and their standard deviation (1σ) at each P ill . Note that, T sample of quartz is not reported in Table S3 (see the paragraphs below for details). The duration of every Raman analysis is 300 or 600 s.  ] sens. calib + δ HB/FD ) and density obtained from a natural fluid inclusion at 26, 40, 60, 80, and 100 • C are plotted as square and triangle symbols, which are measured respectively along the temperature increasing and decreasing cycle (Table S4). δ HB/FD was calculated as − 0.0075 (Table S4).
Raman spectra of the CO 2 fluid were obtained using a micro-Raman spectrum analysis system set up at the Hokkaido University Museum, with 1200 grooves/mm grating. Raman spectra were acquired using diode-pumped solid-state laser (532 nm, Gem 532; Laser Quantum) excitation, a spectrometer with 75 cm focal length (Acton SP-2750; Princeton Instruments, Inc.), and a CCD camera (1650 × 200 pixels, 16-μm width, iVac; Andor Technology). A Nikon microscope (ECLIPSE E600-POL) with a × 50 (N.A. = 0.8, LU Plan; Nikon Corp.) objective was used to focus on the samples. Spectra were collected in a single window between 656 and 1481 cm − 1 .
In the case of fluid inclusions in quartz, which requires higher P ill for accurate determination of B, almost identical measurements were performed using a micro-Raman spectrum analysis system (RAMANtouch VIS-HP-MAST; Nanophoton) at JAMSTEC. The laser output was 1500-4000 mW at the source. Using neutral-density filters with variable transmission, we adjusted P ill at the sample surface from 5 to 351 mW. Raman spectra were acquired using semiconductor laser (532 nm) excitation, a 1200 grooves/mm grating, and a CCD camera (1340 × 400 pixels). In this experiment, the sample was placed on a slide glass, not on a Linkam stage. Therefore, T sample is not measured directly. The ambient temperature during the measurement was monitored using a thermometer (TR-73 U; T&D Corp.) placed next to a microscope and was almost constant at the range of 20.2-22.1 • C. A Nikon objective (N.A. = 0.8, TU Plan ×50) was used to focus on the samples. The duration of every Raman analysis is 30 to 300 s. Spectra were collected in a single window between 559 and 1780 cm − 1 . Sensitivity calibration, spectra fitting, and data correction procedures are described in Supplementary material S1 and Tables S4 and S5. The average of the sensitivity calibrated hot band to Fermi diad intensity ratio at each Table S3.

Specific procedures for laser heating coefficient determination and error analysis
Earlier reports have described linear correlation between sample temperature and P ill (Kagi et al., 1994;Liu et al., 1995;Haro-González et al., 2013;Català et al., 2017). Therefore, we can also expect the linear increase of inclusion temperature with P ill . Fig. 4 presents relations between P ill and the inclusion temperature obtained from simulations. Actual data obtained from spinel03 fi139 are also presented as an example. The open symbols show the simulated relation between P ill and inclusion temperature calculated using finite element method. Parameters used for the simulation are almost identical to the actual experimental conditions and listed in the caption of Fig. 4. The boundary conditions are calculated with varying h g/s from 5 to 50 W/m 2 K. As all measured results show the linear relation between inclusion temperature and P ill (Fig. 4 for an example), we define B as the slope of the linear fit of the measured values. The red and blue circles respectively present the measured data obtained from laser increasing and decreasing cycles. For estimating the best fit straight line to data, we applied bivariate least squares fitting method (York method) using OriginPro9.0 (solid line in Fig. 4). In this case, we obtained the slopes (B Raman ) and those of errors as 6.1 ± 0.6 • C/mW. To construct 68.3% confidence intervals around the calculated slopes, we used the two-tailed t-distribution at N − 2 degrees of freedom (t N− 2, 31.7% ) and the standard errors of the slopes (σ BRaman ). Estimated slopes and those of 68.3% confidence intervals (σ BRaman × t N− 2, 31.7% ) are presented in Table 4.

Simulation of temperature distribution inside the host mineral
In Fig. 5, we present a typical q-z section of the stationary temperature distribution around the inclusion (parameters are those of spinel03 fi139). We set h g/s as 10 W/m 2 K. Similarly to earlier studies, the radial temperature distribution from the optical axis at the laser focus showed logarithmic distance dependence (upper inset of Fig. 5) (Celliers and Fig. 4. Measured and simulated inclusion temperature as a function of laser power. Red and blue filled circles represent the data obtained from the fluid inclusion in spinel (spinel03 fi139) and measured respectively along the laser increasing and decreasing cycles. Error bars for temperature show the standard deviation (n = 5). Those for laser power are 3%. These data are presented in Table S3. A solid line shows the linear fit (York method) to eight data points. The laser heating coefficient determined from the linear fit is 6.1 ± 0.6 • C/mW. Open symbols show the simulated relations between inclusion temperature and laser power. Parameters used for simulation are the following: T sample = 25 • C, λ = 532 nm, P ill = 0-20 mW, w 0 = 1.46 μm, α inc = 0 m − 1 , z inc = 2.6 μm, r inc = 2.0 μm, K inc = 0.08 W/mK, k inc = 0.02 mm 2 /s, H h = 42 μm, R h = 0.534 mm, K h = 5.5 W/mK, k h = 1.79 mm 2 /s, α h = 2030 m − 1 , n h = 1.72, h a/s = 5 W/m 2 K, and h g/s = 5-50 W/m 2 K. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 4
Laser heating coefficients and errors obtained from natural fluid inclusions by using hot band thermometer.  Mao et al., 2005;Haro-González et al., 2013;Català et al., 2017). In addition, reports of earlier studies described that the temperature distribution along the optical axis also shows logarithmic decay (Català et al., 2017). However, our simulation results do not follow such a monotonic decrease (right inset of Fig. 5). Compared to results of earlier studies, our simulated temperature distributions along the optical axis show complex behavior: temperatures take a minimum value at z ~ z inc ; the temperature decreases gradually after taking a maximum value at z > z inc .

Laser heating coefficient as functions of experimental, geometrical, and optical variables
In Fig. 6, we present simulation of the effect of each variable on B by changing each parameter within the range that can occur in the study of inclusions. The parameters used in the simulation are presented in the caption in Fig. 6. Our simulation results showed a linear relation between B and α h (Fig. 6a). Although B and H h have a positive relation throughout whole simulated H h , it almost reaches a plateau when H h was greater than 200 μm (Fig. 6b). For H h < 100 μm, B decreases rapidly as H h decreases (Fig. 6b). Although R h and B show negative correlation, B remains nearly constant even as R h increases when R h exceeds 4 mm (Fig. 6c). The negative correlation between B and K h was confirmed by our simulations (Fig. 6d). k h and n h did not influence B (Fig. 6e and f). We investigated the r inc dependence of B over r inc = 0.1-10 μm at z inc = 10 μm, which revealed that B decreases slightly, concomitantly with increasing r inc (inset in Fig. 6g). Additionally, we simulated the size dependence of B at z inc = 100 μm from r inc = 0.1-50 μm. As a result, the relation between B and r inc was completely the same as that at z inc = 10 μm when r inc ≤ 10 μm, but the amount of change increased dramatically at r inc > 10 μm (inset in Fig. 6g). Our simulation results show an exponential decrease (i.e., B ∝ exp.(− r inc )), whereas those of Català et al. (2017) show logarithmic decay (i.e., B ∝ − log(r inc )). By contrast, in the case of α inc > α h , B clearly increases concomitantly with increasing r inc at constant w 0 (Figs. 6i and 7), probably because the inclusions generate heat more efficiently than the host minerals when α inc > α h . We investigated z inc dependence of B. It decreases near both the air-sample and glass--sample boundaries. However, the variation of B was very slight (Fig. 6h). We investigated variations of B by changing α inc from 10 − 4 to 10 2 cm − 1 and r inc from 0.1 to 2.5 μm while maintaining α h = 1.14 cm − 1 (Fig. 6i). As a result, the effect of α inc on B is negligible when α inc < approx. 0.1-1 cm − 1 , but it cannot be negligible when α inc > approx. 0.1-1 cm − 1 . The transition of the effect of α inc on B occurs around α inc = 1 cm − 1 , which is comparable to α h (1.14 cm − 1 ). Therefore, the impact of α inc on B cannot be negligible when α inc > α h (Fig. 6i). Fig. 6j presents evaluation of the effect of λ dependence of Rayleigh length on B, assuming constant α h (= 1.14 cm − 1 ). As a result, the difference in λ did not affect B between 300 and 1100 nm. Negative correlation between B and w 0 was observed when α inc is 100 and 1 cm − 1 (Fig. 6k). However, B was constant when α inc is 0 and 0.01 cm − 1 , irrespective of w 0 . If α inc < α h and if r inc is constant (each solid line in Fig. 7), then B is constant because heat generation is ineffective at the focus, which has the highest energy density and the highest potential for heat generation. The larger h a/s (or h g/s ) becomes, the smaller B is ( Fig. 6l and Fig. S3).
Among the parameters used for the simulations, α h , H h , R h , α inc , h a/s , and h g/s strongly affected B. In fact, B changed about one order of magnitude or more in the investigated range. However, k h , n h , z inc , and λ do not affect B. In addition, K h , r inc , and w 0 slightly affect B.

Temporal evolution of temperature distribution
We simulated the temporal evolution of the inclusion temperature when the inclusions were irradiated with 0.5, 5, and 50 mW excitation laser, and found that the inclusion temperature approached the steadystate value asymptotically with time (Fig. 8). Simulations were performed by assuming geometrical and material parameters the same as those of spinel03 fi139. Details of parameters are described in the caption of Fig. 8. As a result, t 0.9 (= time required for temperature reaches 90% of its steady-state value) was reached in about 18.5 s. Furthermore, varying P ill from 0.5 to 50 mW does not affect t 0.9 (Fig. S4j). Therefore, even if we analyze the same inclusion with different P ill , the difference does not affect the time required for the system to reach its steady state. However, according to simulation of Català et al. (2017), the temperature of microbeads in water reaches 90% of its steady state value only after 1 ms, which is four orders of magnitude lower than our results. To clarify the causes of such differences in t 0.9 , we investigated the effects of experimental parameters on t 0.9 (Fig. S4). Additionally, we simulated t 0.9 as functions of these parameters under four boundary conditions: h g/s = 0, 10, 1000 W/m 2 K, and Dirichlet boundary condition (Fig. S5). As a result, the only parameters that can change t 0.9 as many as 4 orders of magnitude was the Fig. 5. Simulation of stationary temperature distribution inside the host mineral containing inclusion at z inc = 2.6 μm. Radial temperature distributions at the z inc = 2.6 μm plane following the ln(1/q) decay (upper inset graph). The temperature along the optical axis becomes minimum at z = 2.3 μm (right inset graph). The inset schematic shown at upper left is an enlarged view of the temperature distribution around and within a fluid inclusion. Parameters used for simulations are the following: T sample = 25 • C, λ = 532 nm, P ill = 5 mW, w 0 = 1.46 μm, α inc = 0 m − 1 , z inc = 2.6 μm, r inc = 2.0 μm, K inc = 0.08 W/mK, k inc = 0.02 mm 2 /s, H h = 42 μm, R h = 0.534 mm, K h = 5.5 W/ mK, k h = 1.79 mm 2 /s, α h = 2030 m − 1 , n h = 1.72, h a/s = 5 W/m 2 K, and h g/s = 10 W/m 2 K. It is noteworthy that, although simulations were performed assuming R h = 0.534 mm, this schematic depicts the radial distance up to around 0.42 mm because the isotherm farthest from the optical axis existed at around q = 0.42 mm. Fig. 6. Simulated laser heating coefficient as functions of a-f) host mineral, g-i) inclusion, and j-l) experimental parameters. Parameters used for simulations are fundamentally the following otherwise used as variables in each figure: T sample = 25 • C, λ = 532 nm, P ill = 5 mW, w 0 = 1.46 μm, α inc = 0 m − 1 , z inc = 10 μm, r inc = 2.5 μm, K inc = 0.08 W/mK, k inc = 0.02 mm 2 /s, H h = 200 μm, R h = 1 mm, K h = 5.5 W/mK, k h = 2.05 mm 2 /s, α h = 114 m − 1 , n h = 1.65, h a/s = 5 W/m 2 K, and h g/s = 10 W/m 2 K. boundary condition. Therefore, because Català et al. (2017) set the boundary condition to the Dirichlet condition, the extremely small t 0.9 in their simulations will be attributable to the assumption of high heat flux at the boundary.

Experiments
The total of 19 laser heating coefficients obtained from Raman spectroscopy are presented in Fig. 9 as functions of host mineral and inclusion parameters.
Our experimentally obtained results show linear correlation between α h and B (Fig. 9a). The mean B were quartz << olivine < opx ≈ cpx < < spinel ( Fig. 9a and Table 4). Additionally, we found negative correlations between B and each of H h , R h , and K h (Fig. 9b to d). However, these correlations are not necessarily causal because they might be apparent relations caused by parameters that are highly influential to B (see section 7.1 for details). In our experiments, we measured B Raman of fluid inclusions in olivine with a radius of 4.5-22.2 μm, but no remarkable r inc dependence of B was observed (Fig. 9e). From our simulations, B at r inc = 4.5 μm was only 1.2 times larger than B at r inc = 22.2 μm. Therefore, at least for r inc < 22 μm, we conclude from both the experimental and simulation results that r inc has little impact on B. Our experimentally obtained results show no systematic relation between z inc and B (Fig. 9f). For example, B Raman of the shallowest fluid inclusions for each host mineral species is identical to those of the deepest one, within the error (Fig. 9f). Mao et al. (2005) proposed a laser heating model that describes a temperature profile in the focal plane at the middle of the chamber and parallel to the coverslips. Based on the model, Haro-González et al. (2013) presented the following relation at radial distance w 0 (=0.6λ/ NA) from the optical axis of focus as

Influential variables on the laser heating of inclusions
where K h denotes the thermal conductivity, H h represents the sample thickness, NA stands for the numerical aperture of an objective, and λ shows the laser wavelength.
Although earlier studies have argued that B is proportional to α h /K h (e.g., Català et al., 2017), we found that such conditions are limited to those in which the boundary condition is a Dirichlet condition. Also, the effect of K h on B decreases as the heat transport at the interface decreases (Fig. S3d). The host minerals used for this study have K h of 5.5-9.0 W/mK. Therefore, the expected difference in B is only about twice as large at most. Therefore, K h is not the dominant factor in determining B, at least for the host minerals used for this study. However, if a large difference exists in K h between materials (e.g., diamond and barite), then such K h differences cause differences in B.
As one might infer from Eq. (6), the logarithmic relation between H h and B is also confirmed by our simulation (Fig. 6b). To verify the effect of H h on B from experimental data, we should compare the B of fluid inclusions in the same host mineral species with different thicknesses. Fig. 7. Simulated laser heating coefficient as a function of waist-radius (w 0 ), inclusion radius (r inc ), and absorption coefficient of the inclusion (α inc ). When α inc is 0 cm − 1 , the laser heating coefficient is constant independent of the w 0 , but the laser heating coefficient becomes smaller with increasing r inc . However, when α inc is 10 cm − 1 , the laser heating coefficient increases with decreasing w 0 . When α inc = 10 cm − 1 , the r inc is larger with a larger laser heating coefficient.
From the discussion presented above, the effects of H h and K h on B are small. From our simulations, the λ and w 0 effects on B are negligible, at least when α inc < α h (Fig. 6j and k). Therefore, as expected from Eq.
(6) and earlier studies, the dominant factor which determines B will be α h (Català et al., 2017). Indeed, our simulation and experimentally obtained results showed a linear relation between them (Figs. 6a, 9a, and   Fig. S3a). The spinel has about 10 3 times larger α h than that of quartz: the measured B Raman of spinel is approx. 6 • C/mW, which is 600 times Fig. 9. Relations between material parameters and the laser heating coefficient of fluid inclusions measured by hot band thermometer: a) absorption coefficient of the host mineral, b) host mineral thickness, c) circle equivalent radius of the host mineral, d) thermal conductivity of the host mineral, e) radius of the fluid inclusion, and f) the distance of the fluid inclusion from the sample surface. Each graph shows a total of 19 data obtained using a hot band thermometer (B Raman ). Sample geometry, absorption coefficient, thermal conductivity, and laser heating coefficient are presented in Tables 1, 3, and 4. larger than that of quartz (approx. 0.01 • C/mW). In addition, B of orthopyroxene, clinopyroxene, and olivine, which have an intermediate α h , fall between approx. 0.1 and 1 • C/mW, which lies between those of spinel and quartz. This order is equal to the order of α h except for the order of clinopyroxene and orthopyroxene. Therefore, although it remains possible that parameters other than α h affect B to some extent, our experiments and simulations support that the experimental parameter with the most influence on the B would be α h , within the parameter range that could occur in the inclusion study, as suggested by results of earlier studies (Haro-González et al., 2013;Català et al., 2017).
According to Català et al. (2017), the boundary distance from the optical axis has little effect on B at least boundary distance >80 μm (i.e., R h > 80 μm). However, unlike their results, our simulation results showed negative correlation between R h and B (Fig. 6c and Fig. S3c). This discrepancy is expected to result from the difference in boundary conditions because Català et al. (2017) assumed a Dirichlet boundary condition, whereas we assumed a convective heat flux. In fact, when the glass-sample boundary condition is changed to the Dirichlet boundary condition in our simulation, R h at which the relation between R h and B reaches a plateau decreases drastically from approx. R h = 4 mm to 0.2 mm (red circles in Fig. S3c). Therefore, the discrepancy of the simulation result is explainable by the difference of boundary conditions. To evaluate the effect of R h on B from the experimental data, we need host minerals with different R h . However, we have not prepared such samples in this study. Therefore, the effect of R h on B cannot be discussed further from our experimental data. Nonetheless, earlier reports describe that even quartz, which has very small α h , is heated to a considerable degree by a laser when R h is small (Chio et al., 2003). Therefore, it would be effective NOT to make the host mineral too small during sample preparation to reduce the damage and spectral changes caused by the excitation laser.
Regarding the effects of r inc on B, Peterman et al. (2003) concluded that only a small effect derives from the size of the trapped beads. By contrast, Català et al. (2017) concluded that B logarithmically decreases concomitantly with increasing bead size. It is related directly to the radial temperature profile. However, their data cannot be compared directly to those of the present study because both their experimental and simulation results have been based on the temperature of the buffer (either water or glycerol) rather than in the beads. Both our experimental and simulation results, which evaluated B of the inclusion itself, reveal that r inc has little impact on B, when the inclusion radius is smaller than 22 μm. By contrast, our simulations indicate that the extremely large inclusion (r inc >~30 μm) shows B decrease concomitantly with increasing size (r inc ) because the temperature rise from the focused laser is restricted near the focal point (inset in Figs. 5 and 6g). Such a restriction of the temperature rise region has also been reported by Català et al. (2017). When we assume the case with α inc > α h , possible case for solid inclusions, the larger inclusions are more likely to be heated. The effect of laser heating must be considered when comparing results obtained from inclusions with highly different size (Fig. 7).
The experiments reported by Català et al. (2017) demonstrated that heating was decreased significantly when the laser was focused on a particle in water and placed 10 μm from the bottom of glass-water boundary; whereas B remained almost constant irrespective of depth in the case where the laser focus was distant from the bottom by >10 μm.
They argued that the temperature rise from the focused laser around the focal point was restricted, with volume as small as 10 μm. Therefore, when the focus is close to the glass with a lower α/K, the heat produced by the laser is dissipated efficiently to the glass. In contrast, our simulation results show the B dependence on z inc only when the glass-sample boundary condition was in a Dirichlet condition (Fig. S3h). Even in such a condition, B is almost constant near the z = 0 plane. Furthermore, the experimentally obtained results showed no z inc dependence of B. Raman spectroscopy is usually applied to shallow inclusions rather than the bottom of the sample because the increase in inclusion depth generally makes it difficult to obtain the Raman scattered light of the inclusions efficiently (e.g., Everall, 2010). Therefore, z inc dependence of B is negligible in a normal analysis environment. For that reason, no bias is found in measurements of inclusions with various depths. However, it is noteworthy that the proximity of inclusion to the thin-section surface might modify the residual pressure of inclusions (Mazzucchelli et al., 2018;Zhong et al., 2020). Català et al. (2017) and Peterman et al. (2003) asserted that heating is governed mainly by laser absorption in the surrounding buffer, but not in the particles. In earlier studies, absorption coefficients of beads at their excitation wavelength were 0.06 cm − 1 (polystyrene) and 5 × 10 − 5 cm − 1 (silica), which were sufficiently lower than those of the surrounding buffer: water (0.14 cm − 1 ) and glycerol (0.21 cm − 1 ). Their results are consistent with our simulation results, which suggest that α inc did not affect B when α inc was lower than α h . Here we consider two typical fluid species, CO 2 and H 2 O, as examples. In the case of CO 2 , α inc at 532 nm is expected to be much smaller than that of quartz, which has the smallest α h used for this study. Therefore, it is reasonable to assume α inc = 0 for the evaluation of laser heating in inclusions from Fig. 6i. Similarly, α inc of pure water at 532 nm is approx. 2 × 10 − 4 cm − 1 . Its temperature and salinity dependence are negligible; the assumption of α inc = 0 is also valid for brine fluid inclusions (Pegau et al., 1997;Röttgers et al., 2014). Conversely, the impact of α inc cannot be negligible when α inc > α h (Fig. 6i). For that case, B is controlled strongly by r inc , α inc , and w 0 (Figs. 6i and 7). Therefore, for α inc > α h , bias attributable to laser heating can occur when measuring inclusions of different r inc and α inc or with different objective lenses.
In summary, α h , R h , and h g/s have significant effects on B, but r inc and z inc do not affect B. Therefore, to reduce the temperature rise during the analysis of inclusions, it is effective to make the host mineral larger and thinner and to increase the heat flow rate at the boundary. In addition, α h depends on the wavelength. Therefore, if the optimal excitation laser setting and sample geometry are selected, then the inclusion temperature during the analysis will be suppressed; the obtained data will approach the true value. To help researchers to reach such optimal analytical conditions, α h of representative rock-forming minerals and their wavelength dependence are presented in Table S6.

Correction methods of laser heating effects on Raman spectra
As the inclusion temperature during Raman spectroscopy increases concomitantly with increasing P ill , the effect of laser heating on the Raman spectrum can be minimized by lowering P ill to the greatest extent possible. However, the lowest P ill does not mean that the laser heating effect on the Raman spectrum is completely absent. Therefore, depending on the situation, it is necessary to make some corrections to the measured data. In the following, we discuss how corrections should be made to the Raman spectral features when inclusions are suspected to be heated by a laser.
Therefore, ρ can be found from the measured Δ. However, because the relation between Δ and ρ depends on temperature, if the inclusions are heated to a considerable degree by the laser, then the obtained Δ must be corrected to ascertain a more reliable density specifically (e.g., Wang et al., 2019).
To derive a correction method, the relation between P ill and Δ was investigated using the fluid inclusion in quartz (quartz01 fi56) and that in olivine (En2A oli02 fi13). The analysis was conducted using a grating of 1800 grooves/mm; in all other respects, the analytical procedures are the same as those described in section 5.2. Analyses were conducted along both P ill increasing (circle symbol) and decreasing (triangle symbol) cycles with four P ill of 4.2, 7.7, 10.9, and 14.3 mW (Fig. 10). The measured Δ (Δ raw ) are first calibrated using measured and known distances between the Ne lines according to the explanation presented by Lamadrid et al. (2017) (Δ Ne calib in Table 5). Subsequently, discrepancies in the calibration curves between laboratories were corrected by adding the correction term δ Δ to Δ Ne calib . Finally, corrected Δ are defined as Δ Ne & Qtz calib = Δ Ne calib + δ Δ . δ Δ was calculated using Eq. 10 of Hagiwara et al. (2020), Eq. 8 of Sublett et al. (2020), and the 8 Δ Ne calib obtained from fluid inclusions in quartz. Detailed calculation processes of δ Δ are shown in Table S7. The values of ρ Raman calculated using substituting T sample and Δ Ne & Qtz calib into Eq. 8 of Sublett et al. (2020) and P ill at that time are presented in Table 5. The spectra were obtained five times at each P ill . The mean values of Δ Ne & Qtz calib were presented in Fig. 10b. Their standard deviations (1σ) were shown as error bars.
The relations between P ill and ρ Th − ρ Raman are presented in Fig. 10a.
The solid lines in Fig. 10a show a linear fit (York method) to eight data points. In the case of the inclusion in quartz, the slope was − 0.00019 ± 0.00036 g/cm 3 /mW, showing no significant positive or negative correlation with P ill , at least in the studied P ill range. As determined in this study, the B of fluid inclusion in quartz is about 0.01 • C/mW. Therefore, the expected laser-induced temperature rise is at most 0.2 • C. The Δ drop caused by laser heating is only 0.0006 cm − 1 , which is well below the measurement precision of Δ. If no correlation exists between P ill and the spectral features, then no correction is needed because the laser heating effect is negligible, at least within the studied P ill range.
To ascertain whether spectral features of inclusions in a given mineral are unaffected by P ill , it is desirable to select the inclusions that are most likely to be heated by the laser from all the inclusions to be analyzed. In other words, the fluid inclusions in the host mineral having the thickest and smallest diameter are optimal for evaluation. Some minerals have α h varying with composition. The measurement of inclusions in such minerals must consider B differing in different compositional zones even in a single grain. Indeed, infrared studies of inclusions have revealed that even synthetic fluid inclusions in a single mineral grain, which are expected to have similar fluid composition and density, show different phase-transition temperatures depending on the compositions of the surrounding host phase (Casanova et al., 2018). Therefore, it is important to repeat this work (i.e., obtain the spectra of each inclusion while varying P ill ), even in a single crystal, if the host phase has significant compositional variation.
In the case of inclusion in olivine, the relation between ρ Th − ρ Raman and P ill shows a significant positive correlation; the slope was 0.0012 ± 0.0006 g/cm 3 /mW (Fig. 10a). Earlier reports have described that even for fluids with constant density, ρ Raman is underestimated as the temperature increases (Wang et al., 2011;Yuan et al., 2017;Wang et al., 2019;Sublett et al., 2020), which is consistent with our experimentally obtained results (Fig. 10a). Therefore, to obtain a more accurate ρ Raman , the data must be measured using lower P ill or corrected to account for the laser heating effect.
If the target spectral properties (here Δ) and temperature are almost linearly correlated at constant ρ, then the intercept of the linear fit between Δ and P ill can be regarded as Δ at 0 mW. To verify the linearity of the relation between Δ and temperature at constant ρ, we calculated the relation between them at ρ = 0.698 g/cm 3 using Eq. 8 of Sublett et al. (2020). For the calculation, we assumed B = 1.1 • C/mW so that the obtained slope of the line in the calculations corresponds with that obtained from measured relation between Δ and P ill . As a result, although their equation is not a linear function, it shows an almost linear relation at least within P-T of our interest (black line in Fig. 10a). Therefore, the intercept at the linear fit of the measured Δ (Δ = 104.283 ± 0.016 cm − 1 ) is regarded as being the true value corrected for the laser heating effect Fig. 10. a) Laser power dependence of the difference between the density of CO 2 calculated from the homogenization temperature (ρ Th ) and the density calculated from the Fermi diad splits (ρ Raman ). ρ Raman was calculated by substituting Δ Ne & Qtz calib and the temperature monitored in Linkam stage into Eq. 8 of Sublett et al. (2020). The stage temperature and ρ Raman are presented in Table 5. We measured fluid inclusions in olivine (olivine02 fi13) and quartz (quartz01 fi56), with data plotted respectively as green and red symbols. The measurements were made with the laser power increasing cycle (circle) and decreasing cycle (triangle). In the case of inclusions in olivine, ρ Raman is underestimated with increasing laser power. However, in the case of inclusions in quartz, ρ Raman was kept constant with increasing laser power. b) Relation between Δ of fluid inclusions in olivine (olivine02 fi13) and laser power. The green solid line shows the linear fit (York method) to eight data points. The line slope is − 0.0031 ± 0.0014 cm − 1 /mW. The black line represents calculated relations between laser power and Fermi diad split using Eq. 8 of Sublett et al. (2020). The sample temperature was assumed to be 27.2 • C with laser power of 0 mW. The calculated slopes of the relations between Δ and laser power at laser heating coefficients of 1.1 • C/mW were − 0.0031 cm − 1 /mW. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) because Δ and P ill is linearly correlated (Fig. 10a). Assuming that the error of the correction term δ Δ determined in Table S7 is also 0.016 cm − 1 , the final error of the obtained intercept is ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 0.016 2 + 0.016 2 √ = 0.023 cm − 1 . In addition, because the standard error of Eq. 8 of Sublett et al. (2020) is reported as ±0.006 g/cm 3 , the error of Δ at 0.698 g/cm 3 is expected to be ±0.015 cm − 1 , which is the difference of Δ obtained by Eq. 8 of Sublett et al. (2020) applied for the conditions of ρ = 0.698 g/ cm 3 and ρ = 0.704 (or 0.692) g/cm 3 at the temperature of 27.2 • C.
Therefore, deviation of measured intercept (Δ = 104.283 cm − 1 ) from calculated one (Δ = 104.321 cm − 1 ) is almost equal to the sum of the errors of the measured and calculated Δ. By this series of analyses, the laser heating effect on Δ of a single fluid inclusion can be corrected. However, if the relation between the Raman spectral properties and temperature is not linear in the P-V-T-X range of interest, then the intercept should be obtained by fitting with a function that best describes the relation.
However, it is very time-consuming to perform a series of operations (i.e., to obtain the spectra of each inclusion while varying P ill ) to estimate the intercept values for dozens of inclusions. Fortunately, however, if the temperature dependence of the calibration curve is known, then there is no need to repeat this series of operations if only B can be determined. This is true because, once B is known, the temperature of the inclusions in the analysis is calculable. We can correct the laser heating effect using a calibration curve with known temperature dependence. In the following, we explain how to estimate B from the relation between the Raman spectral features and P ill of an inclusion when the temperature dependence of the calibration curve is known.
For a Raman CO 2 densimeter, the slope of the relation between Δ and P ill (∂Δ/∂P ill ) depends on ρ and B. Consequently, from the relation between Δ, ρ, and T obtained in an earlier study, B satisfying the measured (∂Δ/∂P ill ) ρ is obtainable. According to Eq. 8 of Sublett et al. (2020), B satisfying (∂Δ/∂P ill ) ρ = − 0.0031 ± 0.0014 at ρ = 0.698 g/cm 3 is B = 1.1 ± 0.5 • C/mW (Fig. 10a). This result is consistent with the actual B of oli02 fi13, B Raman = 0.6 ± 0.3 • C/mW (Table 4). Therefore, if the temperature dependence of the calibration curve is known, but B is unknown, then one can find B using the procedure described above. The precision of the measurement of B depends on the degree of temperature dependence of the calibration curve, but it is possible to determine B more accurately by performing the above series of analyses over a wider P ill range. However, if the relation between temperature and target spectral feature at desired P-V-T-X conditions is unknown, then B cannot be estimated by the procedure described above. Therefore, under such conditions, the raw data should be corrected by estimating the value of the intercept from the relation between the Raman spectral features and P ill .
To evaluate the effect of laser heating on the geological application of inclusion analysis, as an example, we consider an example of pressure estimation from a CO 2 fluid inclusion in a mantle xenolith. Oglialoro et al. (2017) showed the existence of a magma storage region located in the lower oceanic crust at a depth of 0.26-0.34 GPa beneath the El Hierro Volcano, Canary Islands, based on the density of CO 2 -rich inclusions in mantle xenoliths. For example, when measuring Δ of a CO 2 fluid inclusion of 0.7 g/cm 3 in spinel assuming P ill = 10 mW and B = 6.0 • C/mW, the inclusion temperature during the analysis reaches 85 • C; the Δ would drop by 0.16 cm − 1 lower than the true value according to Eq. 8 of Sublett et al. (2020). Results show that the density obtained from the measured Δ is 0.636 g/cm 3 , which is 9% lower than the true value. If we calculate the CO 2 trapping depth based on this density and assuming 1000 • C, then the result would be 0.28 GPa, which is 0.06 GPa lower than the true depth (0.34 GPa). Pressures are calculated following the method described by Pitzer and Sterner (1994). Because this deviation is much larger than the uncertainty of depth derived from the error in the temperature estimation (±30 • C corresponding to ±0.01 GPa), laser heating can be a major source of error. However, for inclusions with B = 1.0 • C/mW, such as those hosted in olivine and pyroxene, the underestimate of depth provenance is only 0.01 GPa, which is in similar range to the error derived from the temperature estimation. Therefore, the effect of laser heating is not so great. Therefore, although correcting the effects that laser heating will give data closer to the true values, such corrections are not necessarily required for the geological application when inclusions have small B.
Results of this study show that the use of the lowest P ill can reduce the deviation of the measured Raman spectral feature from the true value. Therefore, the Raman data obtained under such conditions will be approximately equal to the true value. Corrections for the effects of laser heating are not always necessary because those effects depend on the quality of the data which the researcher desires. However, at least when comparing data obtained from different mineral species, the presence or absence of laser heating effect should be described. In addition, if the laser heating effect on the measured values is not negligible, we recommend reporting of how the ideal measurement conditions were reached and what corrections were made to the raw data to interpret them. is defined as "Δ Ne calib + δ Δ " cm − 1 . δ Δ can be found by optimizing Eq. 10 of Hagiwara et al. (2020) to the minimum and is determined as 0.087 cm − 1 .
b ρ Raman is density calculated by substituting T sample and Δ Ne & Qtz calib into Eq. 8 of Sublett et al. (2020). c σρ Raman is density error calculated by substituting T sample and Δ Ne & Qtz calib + σ Δ into Eq. 8 of Sublett et al. (2020).

Conclusions
Experiments and heat transport simulations were combined to determine the laser heating coefficient ( • C/mW) of the inclusions hosted in olivine, orthopyroxene, clinopyroxene, spinel, and quartz. We identify the influential parameters on the laser heating of the inclusions during Raman spectroscopic analysis. The experimental methods estimates B from the slope of the linear fitting of the inclusion temperature and P ill obtained with a hot band thermometer.
The average B of inclusions in each host mineral was in the order of quartz << olivine < opx ≈ cpx < < spinel. Because this order is equal to the order of α h , except that clinopyroxene and orthopyroxene are in reverse order, the most influential experimental parameter on B would be α h , as suggested by earlier studies and our simulations. Aside from the α h , the host mineral thickness and size and the heat transport at the boundaries can be the main factors that control B. Therefore, it is possible to reduce B to some extent by adjusting the shape of the host mineral during the sample preparation process. However, the radius and depth of inclusions have little effect on B if α h > α inc . Therefore, if researchers find that the physicochemical properties of inclusions measured by Raman spectroscopic analysis depend on the inclusion size or the distance from the sample surface, then laser heating can be excluded from the cause of that dependence.
For small host minerals and translucent minerals such as Cr-spinel, the temperature of the inclusions can rise tens or hundreds of degrees, even with the laser power typically used for inclusion analysis. For such situations in which the effects of laser heating are unavoidable, we propose some recommendations for the correction of laser heating effects on measured Raman spectra features. With such treatment, even if the Raman spectral features are altered by laser heating, results obtained from Raman spectroscopy become accurate.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.