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

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). Di ﬀ usion 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 di ﬀ usion along the axis of a one-dimensional line scan. The correction calculations use (1) numerical forward modelling of the di ﬀ used elemental signals and (2) a non-negative signal deconvolution inversion calculation technique to retrieve the original signal pro ﬁ les. 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.


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. 3 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, 6 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 ionmicroprobe 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. 9 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. 12 Reportedly, the spatial resolution of these systems can be improved to about 1 mm when they are equipped with a time-of-ight (TOF)-ICP-MS 13,14 or even with conventional quadrupole ICP-MS instruments using a rapid sample aerosol transfer ablation cell and tubing. 15 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][16][17] particularly for a normal LA system using a largevolume sample cell and a transfer tubing that is ca. 4 mm in diameter (Fig. 1c).
Examination of the uid 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][17][18] Nevertheless, the limited laser beam size and shape and the dilution-diffusion effect remain signicant 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][20][21][22] and the deconvolution of the diffusion-dilution effect in LA aerosol transport using continuous line scan ablation. 23,24 Data science has contributed signicantly 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 articial intelligence. 25,26 Such approaches have been used to improve magnetic resonance imaging (MRI), 27 seismic tomography, 28 geodesic inversion, 29 etc. Following similar datascience 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 mm in a 10% valley denition 23 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 mm s À1 ). 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) 30 in a meteorite from northwest Africa.

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 mm-wide and 500 mm-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. 31 The meteorite was cut into pieces using a lowspeed 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 (hrt1@eqchem.s.u-tokyo.ac.jp).

Analytical methods
The analytical conditions of the ICP-MS and LA are detailed in Table 1. Briey, for the resolution test of a line prole 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 proling of Cr. The ICP-MS instrument was a magnetic sector-based ICP-MS (Nu AttoM, Nu Instruments, Wrexham, UK) that was tuned for 53 Cr. A 5 mm 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 mm s À1 . 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 53 Cr was continuously monitored with a dwell time of 100 ms.
An ArF excimer-based laser ablation system (NWR193, ESI New Wave Research, Montana, USA) was used for twodimensional (imaging) Rh mapping of the CAI meteorite sample. The ICP-MS instrument was quadrupole based (iCAP Qc, Thermo Fisher Scientic, Bremen, Germany). The LA-ICP-MS experimental parameters were optimized to obtain the maximum signal intensities of 63 Cu and 64 Zn by ablating an SRM 610 synthetic glass standard from the National Institute of Standard and Technology (NIST). Single line proles were obtained by moving a circular laser spot 10 mm in diameter with a scanning velocity of 25 mm s À1 . The repetition rate of the laser pulse was set at 20 Hz. A two-dimensional map was obtained by repeated line-proling analyses across the sample. The distance between two adjacent analysis lines was 10 mm. The dwell time of the data acquisition was 10 ms for all analytes (in total 26 elements/isotopes, covering 19 F to 238 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 103 Rh.

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 prole. The resulting line-prole data were visualized in two-dimensional space using MATLAB (The MathWorks, Inc.). The ND algorithm implemented in MATLAB is available from the rst author (TA).

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. 14 The generated aerosols were 10-500 nm in diameter (100-150 nm median), 32 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 ow regime. 33 The dilution completely depends on the conguration of the LA cell and tubing. 16,18 Our two LA-ICP-MS systems exhibited a typical transient signal increase and decay prole 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 ow regime. 33 The element signal y(t) evoked by a single pulse is modeled as follows: 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 prole of the aerosol dilution and washout during the transport process. The temporal prole is shown in Fig. 1c and can be described by second-order and dead-time lags as follows: where s 1 and s 2 are the decay constants, and D is the dead time.
It should be noted that we set an amplitude of g(t) in eqn (2) to satisfy ð N 0 gðtÞdt ¼ 1. This restriction is necessary because the amount of an element in the aerosol sample must be conserved. A least squares tting method using the single-shot ablation data indicated an excellent t (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 x 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: where y consists of M sampling points of the observed signal, because y ¼ [y(x 1 ), y(x 2 ),., y(x M )] T , x consists of N sampling points of the original aerosol mass, i.e., , 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 x 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 s 1 ¼ 0.3401 s, s 2 ¼ 0.0730 s, and D ¼ 0.1107 s for Cr patterns sampled by fs LA and s 1 ¼ 0.3886 s, s 2 ¼ 0.3665 s and D ¼ 0.2072 s for CAI sampled by ns LA. The numbers below decimals indicate signicant digit gures not for measurement accuracy but for the numerical calculation of the inversion analysis. Note that the dead time D 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.

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][35][36] Recently the non-negative estimation methods are attracting researchers' attention, and especially the efficacy of non-negativity is highly regarded in the eld of neuroscience. [34][35][36] We can estimate the original aerosol mass x by solving the simultaneous linear eqn (4), referred to as the generative model. However, nding 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: 34 subject to x $ 0.
Here, the rst 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 l 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 nonnegative, 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 ofx obtained with this method is constrained to have a non-negative value, which reects 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 lter (WF) with the smooth constraint:x ¼ (G T G + lL T L) À1 G T y 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 l in the objective function J. 25 The kth observation was retained as validation data, and the remaining observations were used as training data for estimating x denoted asx \k . Then, we quantied the prediction error by calculating the mean-squared error between the kth validation observation and the corresponding prediction based on the estimate ofx \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 l by searching for the minimum GE for each LA line prole.

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 articial signal with Gaussian noise mimicking the quantized signal in the micrometre-scale particle mass measurements (Fig. 2).
To quantitatively evaluate the proposed method, we used the following synthetic data. First, we synthesized the original aerosol mass x that accurately reects 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 mm in a 10% valley denition. 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 x 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 s 1 ¼ 0.3886 s, s 2 ¼ 0.3665 s, and D ¼ 0.2072 s, as noted above. Furthermore, 1 s corresponded to 25 mm because the stage velocity was 25 mm s À1 in the actual equipment for CAI. The articial 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 conrmed that the cross-validation method could select a nearly optimal l value ( Fig. 2a and c). The ND algorithm with the chosen l could retrieve the original aerosol mass x accurately reecting the particles with a resolution of 30 mm in a 10% valley denition (Fig. 2b), and the ND algorithm outperformed the WF method in terms of accuracy (Fig. 2d). The experiments with synthetic data conrmed that the ND algorithm worked well (Fig. 2a and c). Furthermore, the experiments using synthetic data with Poisson noise conrmed 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).

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 x 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 s 1 ¼ 0.3401 s, s 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 l. 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 l that minimize the TE and GE, respectively. (b) and (d) Estimates of the ND and WF methods, respectively, with an optimal l. Original aerosol mass x (blue lines) and estimate x^(red lines). 0.0730 s, and D ¼ 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.
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, 37 the particle size distribution of the LA aerosols that arises from heterogeneous diffusion during aerosol transfer (see Fig. 1), icker noise from the ICP, 38 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, 38 but the rst 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 proles.
When the Gaussian noise was applied, the prole 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 prole in response to the s-factors and showed rounded rectangular shapes aer 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 prole 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 rst 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 mm length scale, which includes two non-overlapping 5 mm LA craters (Fig. 3c). The time intervals between the base level to the rst peak were ca. 2.5 s (25 mm) in contrast to 1.5 s (15 mm) in the ND model results, indicating an improvement in spatial resolution. Although the ND model offset did not signicantly reduce tailing during the fall of the signal aer the ablation, it eliminated spiky noise signals on the falling slopes that were above the resolution in a 10% valley denition threshold (Fig. 3d). Overall, the ND model improved the resolution of the edges of Cr patterns from 25 to 15 mm measured based on the average length between the rst contact of LA to the Cr patterns and the rst 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.

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-prole analysis for Rh was performed on a straight line approximately 1200 mm in length (Fig. 4a). As shown in Fig. 4d, the ND algorithm with a l that minimizes the GE (inset in Fig. 4c) succeeded in extracting sharp peaks with a resolution of less than 30 mm compared to the original diffuse peaks with a resolution of ca. 50 mm, 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 micrograins 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).
The spatial resolution of the individual peaks is better than 30 mm in a 10% valley denition, so that each identiable peak may originate from a single Rh nugget with a diameter <30 mm or singular to aggregated Rh particles of submicron size dispersed within a region <30 mm in diameter. Measurements using a scanning electron microscope did not nd 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 prole 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 uctuations in ionization efficiency are minor based on the instantaneous mass load if the ablated mass is small assumed. A semiquantitative 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 mm (Fig. 4d). Note that peak deconvolution using a Wiener lter is not appropriate for such a quantitative approach (Fig. 4e).
The line prole 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 signicantly, ranging from <1 to 5 mm in CAI. 7,30,39 This suggests that the size of the metal grains examined in this study is signicantly smaller than 30 mm 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 mm was not achieved. This is simply because of the large LA beam with a diameter of 10 mm scanned at a velocity of 25 mm s À1 with a 20 Hz laser pulse repetition rate. The use of a smaller crater (e.g. 2 mm) with a lower repetition rate (e.g. 2 Hz) and a slower scanning velocity (e.g. 2 mm s À1 ) 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 19 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. 40 Addressing these issues will be the goal for future development of the method described herein.

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 mm resolution (in a 10% valley denition) during rapid macroscopic elemental imaging in a 2.5 Â 2.5 mm area. The ND method improves the spatial resolution by overcoming the diffusiondilution effect of LA aerosol transport with the help of data science methods. Unlike hardware-based approaches, this numerical correction algorithm represents a more rened 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 conicts to declare.