A numerical inversion method for improving the spatial resolution of elemental imaging by laser ablation-inductively coupled plasma-mass spectrometry

Abstract

To improve the spatial resolution of the two-dimensional elemental images of solid organic and inorganic materials, a novel numerical correction method was developed for laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS). Diffusion and dilution of LA aerosol particles in the carrier gas during transportation from the LA cell to the ICP are the major cause of image diffusion along the axis of a one-dimensional line scan. The correction calculations use (1) numerical forward modelling of the diffused elemental signals and (2) a non-negative signal deconvolution inversion calculation technique to retrieve the original signal profiles. This method improves spatial resolution with a semi-quantitative determination of the ablated masses. We used this method to sharpen the spatial distribution images of rhodium particles contained within a meteorite sample.

---

1. Introduction

The distribution of nano- to micro-sized particles in organic and inorganic materials can provide information regarding their mechanisms of transport and concentration in space and time.^1,2 Metal nanoparticles may be toxic to living cells if they are overloaded, and hence, determination of their distributions and transport mechanisms is crucial for biological and medical applications.³ Nano- to micrometre particles are also present in natural inorganic materials and can provide clues to understanding the condensation mechanisms of particular elements such as precious metals in ore formations,^4,5 the effects of nanoscale mineral particles on the environment,⁶ and the origin of pre-solar grains in meteorites.^7,8 Several methods enable simultaneous investigation of the spatial distribution and elemental concentration of nano- to micro-sized particles. Transmission electron microscopy equipped with energy dispersive X-ray spectrometry (TEM-EDX) and secondary ion-microprobe mass spectrometry (e.g. Nano-SIMS) are commonly used to record the images of the two-dimensional distributions and elemental and isotopic compositions of nano- to micro-particles.⁹ However, these approaches are not usually applicable to macroscopic samples, such as human tissues or rock specimens with dimensions larger than a millimetre.

Laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) can generate both macro and microscopic images with extremely high sensitivity, and it has a high throughput capability (Fig. 1a and b).^3,10,11 Recently, LA-ICP-MS imaging systems with submicron- to micron-scale resolution have been developed.¹² Reportedly, the spatial resolution of these systems can be improved to about 1 μm when they are equipped with a time-of-flight (TOF)-ICP-MS^13,14 or even with conventional quadrupole ICP-MS instruments using a rapid sample aerosol transfer ablation cell and tubing.¹⁵ However, improving the spatial resolution of the LA-incorporated method requires more effort. The cause of deterioration in these images is twofold. One is the practical limit of the laser beam size, usually a few to a few dozen micrometres in diameter. This is partly because of the physical diffraction limit of the focusing objective lens of the laser and largely because of the trade-off between elemental sensitivity in the ICP-mass spectrometer and elemental abundances in samples. The second major cause of image deterioration is the aerodynamic diffusion–dilution of LA aerosols during transport from the sampling cell to the ICP ion source^15–17 particularly for a normal LA system using a large-volume sample cell and a transfer tubing that is ca. 4 mm in diameter (Fig. 1c).

Fig. 1 Schematic diagram of the laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) imaging technique. (a) Diagram of the LA-ICP-MS imaging setup. The samples are directly placed in a sampling cell and ablated with a nanosecond or femtosecond laser. The laser-induced sample particles are transported by a carrier gas and finally ionized in the plasma. (b) Diagram of the LA-ICP-MS imaging system. (c) Observed profile of a one-shot signal of ²³⁸U (the mean and standard deviation of 10 samples are shown as a dashed line and error bars, respectively) and fitted the profile using the response function g(t) (red solid line). The t = 0 point does not correspond to the time of the LA shot.

Examination of the fluid dynamics in the sample cell and transfer tubing has led to improvements in the geometric design of the LA-ICP-MS hardware (Fig. 1a and c). Recent hardware developments have improved sample transport with minimal stagnation and enabled fast washout of the sample cell to minimize the dilution–diffusion effects.^16–18 Nevertheless, the limited laser beam size and shape and the dilution–diffusion effect remain significant hurdles for nano- to micrometre-scale analyses in elemental imaging. This is particularly true when rapid imaging over a wide area is required, although this is one of the advantages of LA-ICP-MS over other micro-analytical techniques. To cope with this problem, numerical signal deconvolution approaches have been developed to improve the spatial resolution of elemental imaging. These approaches include the deconvolution of LA signals from overlapping ablation pits in grid ablation or line scanning using a circular laser beam^16,19–22 and the deconvolution of the diffusion–dilution effect in LA aerosol transport using continuous line scan ablation.^23,24

Data science has contributed significantly to these numerical approaches. In particular, recent interdisciplinary studies have focused on inverse problems such as deconvolution for extracting pertinent information from large sets of experimental data by using statistics, machine learning, and other forms of artificial intelligence.^25,26 Such approaches have been used to improve magnetic resonance imaging (MRI),²⁷ seismic tomography,²⁸ geodesic inversion,²⁹etc. Following similar data-science methods, we developed a novel non-negative deconvolution (ND) calculation algorithm to increase the spatial resolution of LA-ICP-MS. We focused on two-dimensional elemental imaging of micrometre-sized metal particles using line scans where the acquired data were heavily distorted because of dilution–diffusion effects during aerosol transport. The metal particles of interest are sufficiently smaller than the laser beam size and are occasionally sampled by the LA. The ND method proposed here improves the spatial resolution from ∼50 to ∼30 μm in a 10% valley definition²³ used for roughly estimate of resolution with a semi-quantitative estimation of the individual particle signals in rapid line-scan mode (laser scan velocity of 25 μm s⁻¹). Herein, we describe the ND algorithm and demonstrate improved imaging of rhodium metal nuggets in an approximately 2.5 × 2.5 mm area of a calcium–aluminium-rich inclusion (CAI)³⁰ in a meteorite from northwest Africa.

2. Sample preparation and analytical methods

2.1. Patterned and meteorite samples

We used a printed metal sample in the evaluation of the spatial resolution. The sample was a high-resolution 1951 USAF Test Target (Edmund Optics, Karlsruhe, Germany) for testing the resolution of optical microscopes. Positive patterns of thin Cr metal were coated on a glass base. We used 100 μm-wide and 500 μm-long rectangular patterns and analysed them in the short-side direction for resolution tests at their edges.

The meteorite sample used for the Rh imaging analysis was a NWA (Northwest Africa) 7678 carbonaceous chondrite (CV3). The dark weathered exterior of the saw-cut sample revealed many chondrules smaller than 1 mm, with a few as large as 5 mm, and scattered CAIs larger than 1 cm that were set in a dark grey groundmass. Nearly all of the Fe–Ni was oxidized by weathering.³¹ The meteorite was cut into pieces using a low-speed diamond cutter (Isomet 1000, Buehler, Lake Bluff, Illinois). Prior to the image analysis using LA-ICP-MS, the surface of the sample was polished using diamond paste (#8000). Neither coatings of conductive materials nor additional matrix components were prepared in this study. Information and access to the NWA7678 meteorite sample are available from one of the authors TH (E-mail: hrt1@eqchem.s.u-tokyo.ac.jp).

2.2. Analytical methods

The analytical conditions of the ICP-MS and LA are detailed in Table 1. Briefly, for the resolution test of a line profile using the 1951 USAF Test Target pattern, a 266 nm femtosecond (fs) laser ablation system (CyberProbe UV, Cyber Laser Inc., Tokyo, Japan) was used for line profiling of Cr. The ICP-MS instrument was a magnetic sector-based ICP-MS (Nu AttoM, Nu Instruments, Wrexham, UK) that was tuned for ⁵³Cr. A 5 μm diameter nanosecond (ns) pulse LA beam was used at a repetition rate of 2 Hz with a line scanning velocity (of the translation stage) of 10 μm s⁻¹. This resulted in 2-shot LA pulses delivered in 1 s to form non-overlapping consecutive craters on the sample surface. The mass peak of ⁵³Cr was continuously monitored with a dwell time of 100 ms.

Table 1 Instrumentation and operational settings

ICP-mass spectrometer
Instrument	Nu AttoM SF-ICP-MS
ICP incident power	1350 W
He carrier gas	740 mL min^−1
Ar makeup gas	850 mL min^−1
Scan mode	Deflector scan
Monitored isotopes	^53Cr
Data acquisition	Time resolved analysis (TRA)
Dwell time	100 ms per isotope
Instrument	iCAP Qc (Thermo Fisher Scientific)
ICP incident power	1550 W
He carrier gas	736 mL min^−1
Ar makeup gas	850 mL min^−1
Scan mode	Standard mode
Monitored isotopes	^19F, ^23Na, ^24Mg, ^27Al, ^29Si, ^31P, ^32S, ^35Cl, ^44Ca, ^57Fe, ^59Co, ^60Ni, ^103Rh, ^139La
	^140Ce, ^146Nd, ^147Sm, ^175Lu, ^178Hf, ^182W, ^193Ir, ^204Pb, ^206Pb, ^208Pb, ^232Th, ^238U
Data acquisition	Time resolved analysis (TRA)
Dwell time	10 ms per isotope
Laser ablation system
Instrument	CyberProbe UV Ti:S femtosecond laser THG
Wavelength	266 nm
Pulse duration	230 fs
Fluence	2.4 J cm^−2
Repetition rate	2 Hz
Ablation pit size	5 μm
Scan speed	10 μm s^−1
Instrument	NWR193 (ESI), ArF excimer laser
Wavelength	193 nm
Pulse duration	5–10 ns
Fluence	2.4 J cm^−2
Repetition rate	20 Hz
Ablation pit size	10 μm
Scan speed	25 μm s^−1

---

An ArF excimer-based laser ablation system (NWR193, ESI New Wave Research, Montana, USA) was used for two-dimensional (imaging) Rh mapping of the CAI meteorite sample. The ICP-MS instrument was quadrupole based (iCAP Qc, Thermo Fisher Scientific, Bremen, Germany). The LA-ICP-MS experimental parameters were optimized to obtain the maximum signal intensities of ⁶³Cu and ⁶⁴Zn by ablating an SRM 610 synthetic glass standard from the National Institute of Standard and Technology (NIST). Single line profiles were obtained by moving a circular laser spot 10 μm in diameter with a scanning velocity of 25 μm s⁻¹. The repetition rate of the laser pulse was set at 20 Hz. A two-dimensional map was obtained by repeated line-profiling analyses across the sample. The distance between two adjacent analysis lines was 10 μm. The dwell time of the data acquisition was 10 ms for all analytes (in total 26 elements/isotopes, covering ¹⁹F to ²³⁸U: Table 1), and the total time for a single mass scan was 0.287 s (∼287 ms), corresponding to the time intervals of the analysis cycles on ¹⁰³Rh.

3. Numerical modelling

The numerical modelling consisted of (1) a forward model of dilution–diffusion processes of LA aerosols during transportation in the tubing between the LA cell and the ICP and (2) ND of elemental signals from a single line profile. The resulting line-profile data were visualized in two-dimensional space using MATLAB (The MathWorks, Inc.). The ND algorithm implemented in MATLAB is available from the first author (TA).

3.1. Generative model

The formulation of the generative model requires a single laser pulse. This pulse was applied for ∼230 fs or 5–10 ns and ablated the sample at the LA site during this time interval.¹⁴ The generated aerosols were 10–500 nm in diameter (100–150 nm median),³² which became suspended in He carrier gas. The aerosol transfer tubing (4 m long with a 4 mm inner diameter) is where most aerosol dilution and washout occur in the laminar flow regime.³³ The dilution completely depends on the configuration of the LA cell and tubing.^16,18 Our two LA-ICP-MS systems exhibited a typical transient signal increase and decay profile from a single laser pulse over 3000 ms (Fig. 1c), although they differed slightly in their details (see below).

A simple linear dynamic model mathematically consists of a dead-time lag and a higher-order lag representing the aerosol transport through the tubing in the laminar flow regime.³³ The element signal y(t) evoked by a single pulse is modeled as follows:

y(t) = xg(t), (1)

where x is the total aerosol mass (counts) ablated by a single laser pulse at t = 0, and g(t) is a single-pulse response function that describes the temporal profile of the aerosol dilution and washout during the transport process. The temporal profile is shown in Fig. 1c and can be described by second-order and dead-time lags as follows:where τ₁ and τ₂ are the decay constants, and Δ is the dead time. It should be noted that we set an amplitude of g(t) in eqn (2) to satisfy . This restriction is necessary because the amount of an element in the aerosol sample must be conserved. A least squares fitting method using the single-shot ablation data indicated an excellent fit (Fig. 1c).

The LA-ICP-MS technique uses multiple line scans to generate images like those shown in Fig. 1b. Under linear conditions, the element signal evoked by N laser pulses can be given by a linear summation of single-pulse response functions (eqn (1)), as follows:

where x_i is the total aerosol mass ablated by the i-th laser pulse at t = t_i (i = 1,…, N). Denoting the discrete time of the data acquisition as ξ_k (k = 1,…, M), the generative model for the element signal evoked by scanning with N laser pulses can be formulated using the following matrix representation:

y = Gx + s, (4)

where y consists of M sampling points of the observed signal, because y = [y(ξ₁), y(ξ₂),…, y(ξ_M)]^T, x consists of N sampling points of the original aerosol mass, i.e., x = [x₁, x₂,…, x_N]^T, G is an M × N Green's function matrix consisting of elements G_ki = g(ξ_k − t_i), and s is the observation noise.

In the measurement of the Cr patterns printed on a glass substrate, the fs laser pulse irradiation interval was set to 0.5 s and the signal acquisition period was 0.10 s, while in the measurement of Rh in the CAI inclusions, the ns laser pulse irradiation interval was set to 0.05 s and the signal acquisition period was 0.287 s (Table 1). Thus, the discrete intervals of t_i and ξ_k in the M × N Green's function matrix G were set to the individual interval values in each of the two cases. The values of the single-pulse-evoked g(t) parameters were determined to be τ₁ = 0.3401 s, τ₂ = 0.0730 s, and Δ = 0.1107 s for Cr patterns sampled by fs LA and τ₁ = 0.3886 s, τ₂ = 0.3665 s and Δ = 0.2072 s for CAI sampled by ns LA. The numbers below decimals indicate significant digit figures not for measurement accuracy but for the numerical calculation of the inversion analysis. Note that the dead time Δ calculated in this model does not correspond to the actual delay time of the equipment, as t = 0 in the abscissa does not correlate to the time of the LA shot. This discrepancy does not affect following discussions.

3.2. Non-negative deconvolution

We employed the ND method, which was previously shown to be effective in the detection of neuronal spikes from noisy calcium imaging data of multiple neuronal activities in the nervous system.^34–36 Recently the non-negative estimation methods are attracting researchers' attention, and especially the efficacy of non-negativity is highly regarded in the field of neuroscience.^34–36 We can estimate the original aerosol mass x by solving the simultaneous linear eqn (4), referred to as the generative model. However, finding a unique solution to eqn (4) is impossible because matrix G is not full rank (i.e. the number of time points of x, N, is different from that of y, M). This problem is thus categorized as ill-posed. To address this issue, we formulated the ND method as the following optimization problem:³⁴

J = ‖y − Gx‖² + λ‖Lx‖² (5)

subject to x ≥ 0.

Here, the first term of J is the square error that measures the residual between the generative model and observed signal. The second term in eqn (5) is an L2-norm regularisation term to overcome the ill-posed problem, where λ is a regularisation parameter. L in the regularisation term is an N × N matrix given by

L gives smoothness to x and thus makes the estimates robust against noise. Furthermore, to improve the quality of the solution, we introduced a strict constraint that x is always non-negative, which is consistent with the physical constraint that the ion signal intensity must be non-negative.

We were able to determine a suitable value for the regularisation parameter by using the cross-validation method. This optimization problem can be solved using the quadratic programming method. The estimate of x̂ obtained with this method is constrained to have a non-negative value, which reflects the physical constraint that the ion signal intensity must be non-negative. It should be noted that if the objective function J is minimized without non-negative constraints, the conventional Wiener filter (WF) with the smooth constraint: x̂ = (G^TG + λL^TL)⁻¹G^Ty can be obtained. To determine the efficacy of the non-negative constraint, we compared the performance of the ND algorithm with that of the WF with the smooth constraint.

We performed a leave-one-out cross-validation to determine the value of the regularisation parameter λ in the objective function J.²⁵ The kth observation was retained as validation data, and the remaining observations were used as training data for estimating x denoted as x̂_\k. Then, we quantified the prediction error by calculating the mean-squared error between the kth validation observation and the corresponding prediction based on the estimate of x̂_\k. By repeating the calculation of the prediction error with each of the observations and averaging them, we obtained the generalization error (GE). We then determined the value of the regularisation parameter λ by searching for the minimum GE for each LA line profile.

3.3. Numerical experiments with synthetic data

Numerical experiments with synthetic data were performed, as is common practice in data science to evaluate the validity and effectiveness of a method. Here, we synthesized an artificial signal with Gaussian noise mimicking the quantized signal in the micrometre-scale particle mass measurements (Fig. 2).

Fig. 2 Evaluation of the performance of the ND and WF methods using synthetic data. (a) and (c) True errors (TEs) and generalization errors (GEs) of ND and WF methods as a function of λ. The TEs were quantified by calculating the mean square error between the original aerosol mass x and the estimate x̂. Arrows on blue and red lines indicate the values of λ that minimize the TE and GE, respectively. (b) and (d) Estimates of the ND and WF methods, respectively, with an optimal λ. Original aerosol mass x (blue lines) and estimate x̂ (red lines).

To quantitatively evaluate the proposed method, we used the following synthetic data. First, we synthesized the original aerosol mass x that accurately reflects the typically measured distributions. The positions of the particles were randomly generated using a Bernoulli process, and the spatial elemental distribution of each of the particles was represented by a Gaussian distribution with a resolution of 30 μm in a 10% valley definition. The peak value of each of the elemental distributions was randomly generated from a Gamma distribution. Next, using the generative model described by eqn (4), we synthesized the observed signal y with Gaussian noise.

The time scale of the synthetic data was adapted to mirror that of the real data. In the actual equipment, the interval of the laser pulse was 0.05 s, and the signal acquisition period was 0.287 s for the CAI. The discrete intervals of t_i and ξ_k in the generative model (eqn (4)) were set to 0.05 s and 0.287 s for ns LA. The parameter values of g(t) were set as τ₁ = 0.3886 s, τ₂ = 0.3665 s, and Δ = 0.2072 s, as noted above. Furthermore, 1 s corresponded to 25 μm because the stage velocity was 25 μm s⁻¹ in the actual equipment for CAI. The artificial data used in the numerical experiments were composed of an 80 s time sequence (Fig. 2a and b).

To determine the effect of the non-negative constraint on the model performance of particle detection in CAI, we compared the performance of the ND with that of the WF. The experiments confirmed that the cross-validation method could select a nearly optimal λ value (Fig. 2a and c). The ND algorithm with the chosen λ could retrieve the original aerosol mass x accurately reflecting the particles with a resolution of 30 μm in a 10% valley definition (Fig. 2b), and the ND algorithm outperformed the WF method in terms of accuracy (Fig. 2d). The experiments with synthetic data confirmed that the ND algorithm worked well (Fig. 2a and c). Furthermore, the experiments using synthetic data with Poisson noise confirmed that the ND algorithm could work better than the WF when counting rates are low at small signal sizes or small laser spot sizes used (data not shown).

3.4. Resolution test using a 1951 USAF test target

As a further resolution test, we performed numerical experiments and actual measurements on Cr patterns printed on a glass base of a 1951 USAF test target. The interval of the fs laser pulse irradiation was 0.5 s, and the signal acquisition period was 0.10 s. The discrete intervals of t_i and ξ_k in the generative model (eqn (4)) were set to 0.5 s and 0.01 s, respectively. The g(t) parameter values were set at τ₁ = 0.3401 s, τ₂ = 0.0730 s, and Δ = 0.1107 s for the fs LA. Fig. 3a and b show the results of numerical experiments using a LA repetition rate of 2 Hz without overlap of LA craters. The dwell time of the data acquisition was set at 100 ms, and ablation was applied for 12 s between 10 and 22 s, etc. Note that the abscissas in Fig. 3a and b are not a length scale because of the dependence on the crater diameter used for the non-overlap measurement.

Fig. 3 Evaluation of the performance of the ND method using synthetic and measured line profile data. (a) Synthetic data with Gaussian noise. (b) Synthetic data with Poisson noise. (c) Photomicrographs of Cr positive print patterns on the 1951 USAF test target. The central part of the crater is effectively ablated more than the perimeters because of the Gaussian laser beam profile of the fs-LA. The circles indicate the crater size after >10 laser pulses. Single pulse craters formed during the experiment look smaller. (d) Measured and modelled data from Cr patterns. The red lines are synthetic and measured data. The blue lines are results modelled by the ND method. The shaded areas show signals acquired from samples. Thin horizontal lines show average intensities and 10% valleys. (e) Calculation results of the regularisation parameter λ for (a), (b), and (d).

In practice, measurements on LA signals tend to include high-frequency noise; such noise originates from irregular ablation of the sample surface due to LA–sample surface coupling,³⁷ the particle size distribution of the LA aerosols that arises from heterogeneous diffusion during aerosol transfer (see Fig. 1), flicker noise from the ICP,³⁸ and quantisation noise at the analogue–digital converter. The last two sources are known to affect analytical precision when the LA signals are close to the limit of detection,³⁸ but the first two affect even strong LA signals. To simulate their effects, we added both Gaussian (Fig. 3a) and Gamma noise (Fig. 3b) to the synthetic LA profiles.

When the Gaussian noise was applied, the profile obtained by the synthetic experiment reproduced well the onset and offset of the ablation of a hypothetical Cr pattern having very sharp edges (Fig. 3a). The synthetic LA signal with Gaussian noise (Fig. 3a, red) strictly followed the profile in response to the τ-factors and showed rounded rectangular shapes after the onset/offset of the sample ablation. The intensity curve calculated with the ND model (Fig. 3a, blue) reproduced quite well the original rectangular signal profile with almost vertical edges at the onset and offset of Cr ablation. When gamma noise was applied (Fig. 3b, red), the calculated ND slopes were slanted during rise/fall at the onset/offset of ablation (Fig. 3b, blue). Time intervals during the rise and fall of the synthetic signal were ca. 1.5 s between the base level to the first peak and between the last peak to the base level (red), whereas they were reduced to ca. 1 s by the ND model (blue).

The measured LA signals from the 1951 USAF test target (Fig. 3c) were more irregular than the synthetic signals with Poisson noise. In Fig. 3d, the horizontal scale of 1 s correlates with a 10 μm length scale, which includes two non-overlapping 5 μm LA craters (Fig. 3c). The time intervals between the base level to the first peak were ca. 2.5 s (25 μm) in contrast to 1.5 s (15 μm) in the ND model results, indicating an improvement in spatial resolution. Although the ND model offset did not significantly reduce tailing during the fall of the signal after the ablation, it eliminated spiky noise signals on the falling slopes that were above the resolution in a 10% valley definition threshold (Fig. 3d). Overall, the ND model improved the resolution of the edges of Cr patterns from 25 to 15 μm measured based on the average length between the first contact of LA to the Cr patterns and the first peaks of the Cr signal. It should be noted that no resolution improvement can be expected from the model when the signal-to-noise ratio is close to the lower limit of detection.

3.5. Extraction of rhodium signals from CAI

As a practical application, we measured metal particles containing Rh in CAI inclusions in a Northwest African meteorite (Fig. 4). First, a line-profile analysis for Rh was performed on a straight line approximately 1200 μm in length (Fig. 4a). As shown in Fig. 4d, the ND algorithm with a λ that minimizes the GE (inset in Fig. 4c) succeeded in extracting sharp peaks with a resolution of less than 30 μm compared to the original diffuse peaks with a resolution of ca. 50 μm, based on the average widths of peaks measured at 10% peak heights. Some of the originally broad single peaks were clearly separated into two or three sharp discrete peaks, thereby showing an improvement in spatial resolution (Fig. 4d). These results indicate that the resolution is improved for spiky signals from Rh-rich micro-grains in the CAI, although high-frequency noise signals were smoothed out by the model and signals from minor Rh particles may not have been detected because they were close to the lower limit of detection. We also tried a ND model excluding the smoothing effect of the high-frequency noise; the solutions of GE did not differ between the models (not shown).

Fig. 4 (a) Optical image of the analyzed CAI. (b) A corresponding uncorrected LA-ICP-MS image of the Rh distribution. (c) Image obtained from the ND correction. Comparison of the (d) ND-corrected and (e) WF-corrected line profiles of the Rh distribution using LA-ICP-MS imaging of calcium–aluminum-rich inclusion (CAIs) in NWA at a scan speed of 25 μm s⁻¹. The insets (d) and (e) show the GE of the ND and WF methods as a function of λ.

The spatial resolution of the individual peaks is better than 30 μm in a 10% valley definition, so that each identifiable peak may originate from a single Rh nugget with a diameter <30 μm or singular to aggregated Rh particles of submicron size dispersed within a region <30 μm in diameter. Measurements using a scanning electron microscope did not find any large Rh particles of micrometre size (not shown); the appropriate interpretation is thus sub-micrometre-scale Rh-rich particles. Given that the normalized single-pulse response function conserves the total ion count and the non-negative constraint maintains the non-negativity of the estimated mass (Fig. 4d), the aerial integral of the signal intensity profile closely reproduces the sub- to tens of micrometre-scale particle mass volume when the LA-ICP-MS sensitivity factor is provided for various matrices and 100% transfer of the ablated sample to the ICP and the fluctuations in ionization efficiency are minor based on the instantaneous mass load if the ablated mass is small assumed. A semi-quantitative treatment can be applied using the areas of the deconvolution peaks, even without the sensitivity factor. Therefore, the area of each peak corresponds to the integrated signal intensity from the aggregated particles. This indicates that smaller amounts of Rh metal particles result in smaller peak areas and they are spatially resolvable when they are separated by >30 μm (Fig. 4d). Note that peak deconvolution using a Wiener filter is not appropriate for such a quantitative approach (Fig. 4e).

The line profile analysis was followed by Rh imaging analysis of a 2.5 × 2.5 mm square area (Fig. 4a–c). Typical sizes of metal grains in refractory siderophile elements, including Rh, Ir, and Os, can vary significantly, ranging from <1 to 5 μm in CAI.^7,30,39 This suggests that the size of the metal grains examined in this study is significantly smaller than 30 μm and the grains are dispersed as single or aggregated metal particles within the regions of bright spots of Rh concentrations scattered over the 2.5 × 2.5 mm measured area in the CAI (Fig. 4a). The resulting Rh distribution images with and without ND corrections are compared in Fig. 4b and c. These images demonstrate the considerable sharpening of peaks corresponding to the Rh nuggets.

Under the ablation conditions used in this study, a spatial resolution better than 10 μm was not achieved. This is simply because of the large LA beam with a diameter of 10 μm scanned at a velocity of 25 μm s⁻¹ with a 20 Hz laser pulse repetition rate. The use of a smaller crater (e.g. 2 μm) with a lower repetition rate (e.g. 2 Hz) and a slower scanning velocity (e.g. 2 μm s⁻¹) in combination with a rapid aerosol transfer cell and tubing can further improve the spatial resolution. The improvement can only be realized when the transient sampling of the circular laser ablation pits is properly decomposed¹⁹ in tandem with the ND method for diffusion–dilution in aerosol transport, as proposed in this study. This double deconvolution can deal with the two major causes of deterioration of spatial resolution in LA imaging (see the Introduction). As such, the model calculations will be the ultimate challenge in the numerical modelling approach for LA-ICP-MS imaging. For this purpose, further hardware development for high-speed signal integration and data processing in the ICP-MS will be required to minimize the quantization problem of the acquired LA signals.⁴⁰ Addressing these issues will be the goal for future development of the method described herein.

4. Conclusions

A novel numerical modelling method, called non-negative deconvolution (ND), that uses LA-ICP-MS was shown to improve the spatial resolution of micrometre-sized metal particle distribution imaging and semi-quantitative measurements of individual particle/aggregate mass volumes at ∼30 μm resolution (in a 10% valley definition) during rapid macroscopic elemental imaging in a 2.5 × 2.5 mm area. The ND method improves the spatial resolution by overcoming the diffusion–dilution effect of LA aerosol transport with the help of data science methods. Unlike hardware-based approaches, this numerical correction algorithm represents a more refined way in the deconvolution of the LA-ICP-MS signal. The method is also applicable to organic materials of various sizes ranging from cells to tissues.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors are grateful for support from the Japan Society for Promotion of Science, KAKENHI JP15H02148, JP16H01123, and JP18H04372 to J.-I. K., JP25120005, JP25280090, and JP15K20864 to T. K., JP16H01539 and JP16H01269 to T. A., and JP26247094 to T. H. This research was also supported by a PRESTO grant (Grant Number JPMJPR1676) from the Japan Science and Technology Agency. We would like to thank comments from the journal editor, two anonymous and an adjudicative reviewers which greatly improved the manuscript. We also thank to comments from an anonymous reader's given to the early view version of this paper.

References

1. K. Galler, K. Brautigam, C. Grosse, J. Popp and U. Neugebauer, Analyst, 2014, 139, 1237–1273,  10.1039/c3an01939j.
2. R. Trouillon, M. K. Passarelli, J. Wang, M. E. Kurczy and A. G. Ewing, Anal. Chem., 2013, 85, 522–542,  DOI:10.1021/ac303290s.
3. H. A. O. Wang, D. Grolimund, C. Giesen, C. N. Borca, J. R. H. Shaw-Stewart, B. Bodenmiller and D. Günther, Anal. Chem., 2013, 85, 10107–10116,  DOI:10.1021/ac400996x.
4. H. M. Helmy, C. Ballhaus, R. O. C. Fonseca, R. Wirth, T. Nagel and M. Tredoux, Nat. Commun., 2013, 4, 2405,  DOI:10.1038/ncomms3405.
5. J. M. González-Jiménez and M. Reich, Ore Geol. Rev., 2017, 81, 1236–1248,  DOI:10.1016/j.oregeorev.2016.06.022.
6. M. F. Hochella, S. K. Lower, P. A. Maurice, R. L. Penn, N. Sahai, D. L. Sparks and B. S. Twining, Science, 2008, 319, 1631–1635,  DOI:10.1126/science.1141134.
7. A. Bischoff and H. Palme, Geochim. Cosmochim. Acta, 1987, 51, 2733–2748,  DOI:10.1016/0016-7037(87)90153-0.
8. C. Floss, F. J. Stadermann, J. P. Bradley, Z. R. Dai, S. Bajt, G. Graham and A. S. Lea, Geochim. Cosmochim. Acta, 2006, 70, 2371–2399,  DOI:10.1016/j.gca.2006.01.023.
9. Y. Lin, R. Trouillon, G. Safina and A. G. Ewing, Anal. Chem., 2011, 83, 4369–4392,  DOI:10.1021/ac2009838.
10. L. Mueller, H. Traub, N. Jakubowski, D. Drescher, V. I. Baranov and J. Kneipp, Anal. Bioanal. Chem., 2014, 406, 6963–6977,  DOI:10.1007/s00216-014-8143-7.
11. I. Konz, B. Fernández, M. L. Fernández, R. Pereiro and A. Sanz-Medel, Anal. Bioanal. Chem., 2012, 403, 2113–2125,  DOI:10.1007/s00216-012-6023-6.
12. C. Giesen, L. Waentig, T. Mairinger, D. Drescher, J. Kneipp, P. H. Roos, U. Panne and N. Jakubowski, J. Anal. At. Spectrom., 2011, 26, 2160–2165,  10.1039/c1ja10227c.
13. C. Giesen, H. A. O. Wang, D. Schapiro, N. Zivanovic, A. Jacobs, B. Hattendorf, P. J. Schuffler, D. Grolimund, J. M. Buhmann, S. Brandt, Z. Varga, P. J. Wild, D. Gunther and B. Bodenmiller, Nat. Methods, 2014, 11, 417–422,  DOI:10.1038/nmeth.2869.
14. R. N. Straus, A. Carew, D. Sandkuijl, T. Closson, V. I. Baranov and A. Loboda, J. Anal. At. Spectrom., 2017, 32, 1044–1051,  10.1039/c6ja00420b.
15. S. J. M. Van Malderen, J. T. van Elteren and F. Vanhaecke, Anal. Chem., 2015, 87, 6125–6132,  DOI:10.1021/acs.analchem.5b00700.
16. S. J. M. Van Malderen, A. J. Managh, B. L. Sharp and F. Vanhaecke, J. Anal. At. Spectrom., 2016, 31, 423–439,  10.1039/c5ja00430f.
17. S. J. M. Van Malderen, J. T. van Elteren and F. Vanhaecke, J. Anal. At. Spectrom., 2015, 30, 119–125,  10.1039/c4ja00137k.
18. D. N. Douglas, A. J. Managh, H. J. Reid and B. L. Sharp, Anal. Chem., 2015, 87, 11285–11294,  DOI:10.1021/acs.analchem.5b02466.
19. M. Fox, A. Tripathy-Lang and D. L. Shuster, Chem. Geol., 2017, 467, 30–41,  DOI:10.1016/j.chemgeo.2017.07.001.
20. A. Plotnikov, C. Vogt and K. Wetzig, J. Anal. At. Spectrom., 2002, 17, 1114–1120,  10.1039/b202258c.
21. A. Plotnikov, C. Vogt, K. Wetzig and A. Kyriakopoulos, Spectrochim. Acta, Part B, 2008, 63, 474–483,  DOI:10.1016/j.sab.2007.12.012.
22. J. T. van Elteren, A. Izmer, V. S. Šelih and F. Vanhaecke, Anal. Chem., 2016, 88, 7413–7420,  DOI:10.1021/acs.analchem.6b02052.
23. J. T. van Elteren and F. Vanhaecke, J. Anal. At. Spectrom., 2016, 31, 1998–2004,  10.1039/c6ja00234j.
24. J. Triglav, J. T. van Elteren and V. S. Šelih, Anal. Chem., 2010, 82, 8153–8160,  DOI:10.1021/ac1014832.
25. C. M. Bishop, Pattern Recognition and Machine Learning, Springer, New York, 2006.
26. Y. Igarashi, K. Nagata, T. Kuwatani, T. Omori, Y. Nakanishi-Ohno and M. Okada, J. Phys.: Conf. Ser., 2016, 699, 012001,  DOI:10.1088/1742-6596/699/1/012001.
27. M. Lustig, D. Donoho and J. M. Pauly, Magn. Reson. Med., 2007, 58, 1182–1195,  DOI:10.1002/mrm.21391.
28. T. Kuwatani, K. Nagata, M. Okada and T. Toriumi, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2014, 90, 042137,  DOI:10.1103/physreve.90.042137.
29. R. Nakata, H. Hino, T. Kuwatani, S. Yoshioka, M. Okada and T. Hori, Sci. Rep., 2017, 7, 6129,  DOI:10.1038/s41598-017-06185-0.
30. G. J. MacPherson, in Treatise on Geochemistry, ed. K. K. Turekian, Pergamon, Oxford, 2007, pp. 1–47.
31. The Meteoritical Society, https://www.lpi.usra.edu/meteor/about.php, 2017.
32. J.-I. Kimura and Q. Chang, J. Anal. At. Spectrom., 2012, 27, 1549–1559,  10.1039/c2ja10344c.
33. J. G. Balchen and K. I. Mummé, Process Control : Structures and Applications, Van Nostrand Reinhold, New York, 1988.
34. J. T. Vogelstein, A. M. Packer, T. A. Machado, T. Sippy, B. Babadi, R. Yuste and L. Paninski, J. Neurophysiol., 2010, 104, 3691–3704,  DOI:10.1152/jn.01073.2009.
35. E. A. Pnevmatikakis and L. Paninski, Adv. Neural Inf. Process. Syst, 2013, 26, 1250–1258.
36. E. A. Pnevmatikakis, D. Soudry, Y. Gao, T. A. Machado, J. Merel, D. Pfau, T. Reardon, Y. Mu, C. Lacefield, W. Yang, M. Ahrens, R. Bruno, T. M. Jessell, D. S. Peterka, R. Yuste and L. Paninski, Neuron, 2016, 89, 285–299,  DOI:10.1016/j.neuron.2015.11.037.
37. M. Guillong and D. Gunther, J. Anal. At. Spectrom., 2002, 17, 831–837,  10.1039/b202988j.
38. R. Thomas, Practical Guide to ICP-MS, Marcel Dekker, Inc., New York, 2004.
39. J. T. Armstrong, A. E. Goresy and G. J. Wasserburg, Geochim. Cosmochim. Acta, 1985, 49, 1001–1022,  DOI:10.1016/0016-7037(85)90315-1.
40. M. D. Montaño, J. W. Olesik, A. G. Barber, K. Challis and J. F. Ranville, Anal. Bioanal. Chem., 2016, 408, 5053–5074,  DOI:10.1007/s00216-016-9676-8.

---