Franc ¸ois Varray Maxime Polichetti

Member, IEEE Bruno Gilles

Jean-Christophe Béra

Member, IEEE Barbara Nicolas

F Varray

Use of the cross-spectral density matrix for enhanced passive ultrasound imaging of cavitation

Keywords: Array processing, adaptive beamforming, passive ultrasound imaging, cavitation imaging I

Passive ultrasound imaging is of great interest for cavitation monitoring. Spatio-temporal monitoring of cavitation bubbles in therapeutic applications is possible using an ultrasound imaging probe to passively receive the acoustic signals from the bubbles. Fourier-domain beamformers have been proposed to process the signals received into maps of the spatial localization of cavitation activity, with reduced computing times with respect to the time-domain approach, and to take advantage of frequency selectivity for cavitation regime characterization. The approaches proposed have been mainly nonadaptive, and these have suffered from low resolution and contrast, due to the many reconstruction artifacts. Inspired by the array-processing literature, and in the context of passive ultrasound imaging of cavitation, we propose here a robust estimation of the second-order statistics of data through spatial covariance matrices in the Fourier domain, or cross-spectral density matrices. The benefits of such formalism are illustrated using advanced reconstruction algorithms, such as the robust Capon beamformer, the Pisarenko class beamformer, and the multiple signal classification approach. Through both simulations and experiments in a water tank, we demonstrate that enhanced localization of cavitation activity (i.e., improved resolution and contrast with respect to nonadaptive approaches) is compatible with the rapid and frequency-selective approaches of the Fourier domain. Robust estimation of the cross-spectral density matrix and the derived adaptive beamformers pave the way to the development of powerful passive ultrasound imaging tools.

last case). Cavitation corresponds to the formation, oscillation, and implosion of gas bubbles in a liquid that is subjected to pressure variations induced by HIFU. As functions of the excitation conditions, different effects are induced [START_REF] Wan | Cavitation in biomedicine[END_REF]: the oscillations of the bubbles can cause tissue heating and moderate mechanical effects (i.e., stable cavitation), or the implosions of the bubbles can induce destructive mechanical damage (i.e., inertial cavitation). The dynamics of these bubbles are complex, nonlinear, and hard to predict, especially due to the large number of tissue characteristics involved, such that cavitation effects might appear outside of the geometric focus of the instrument during the therapy [START_REF] Haworth | Using passive cavitation images to classify high-intensity focused ultrasound lesions[END_REF]. Thus, the spatio-temporal monitoring of cavitation bubbles is essential to control the actual treated area [START_REF] Wan | Cavitation in biomedicine[END_REF].

To localize cavitation, an ultrasound imaging probe composed of piezoelectric elements is used that passively receives acoustic signals from cavitation bubbles. These signals are processed via array-processing algorithms, which are also called beamformers, and which provide spatial mapping of the cavitation activity. The first method for passive cavitation imaging was proposed in [START_REF] Gyöngy | Passive spatial mapping of inertial cavitation during HIFU exposure[END_REF], which was inspired by the seismic imaging method known as 'time exposure acoustics' [START_REF] Norton | Time exposure acoustics[END_REF], [START_REF] Norton | Passive imaging of underground acoustic sources[END_REF]. At each point of the image, it focuses, combines, and integrates the signals received over time (thus here this is referred to as TD-DAS, for 'delay-and-sum in the timedomain'). In this way, the acoustic power emitted by the medium is mapped. However, TD-DAS has three major limitations that have been frequently reported in the literature for many applications: [START_REF] Escoffre | Therapeutic ultrasound[END_REF] it has very low resolution, and will not always localize a small spot of cavitation precisely; (2) it does not directly differentiate stable cavitation from inertial cavitation; and (3) the focusing of the signals in the time domain is very time consuming in terms of the calculations. As a result, several variants of TD-DAS have been proposed.

On the one hand, to improve spatial localization of cavitation, adaptive methods have been proposed [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF]- [START_REF] Polichetti | Advanced beamforming techniques for passive imaging of stable and inertial cavitation[END_REF], which were inspired by active ultrasound imaging [START_REF] Synnevag | Adaptive beamforming applied to medical ultrasound imaging[END_REF]- [START_REF] Polichetti | A nonlinear beamformer based on p-th root compression -Application to plane wave ultrasound imaging[END_REF]. In the time domain, robust Capon beamforming (TD-RCap) [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF], [START_REF] Stoica | Robust Capon beamforming[END_REF] has been adapted for passive cavitation imaging [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], [START_REF] Lu | Real-time monitoring of controllable cavitation erosion in a vessel phantom with passive acoustic mapping[END_REF], [START_REF] Lu | Passive acoustic mapping of cavitation using eigenspace-based robust Capon beamformer in ultrasound therapy[END_REF]. Based on pixel-wise optimization for combination of the delayed signals, this provides better resolution and better sidelobe suppression than the conventional approach, while being robust to model errors. However, this approach still operates in the time domain, where its optimizations makes it more time consuming than standard TD-DAS.

On the other hand, variants of TD-DAS have been proposed in the Fourier domain (FD) [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF]. The focusing of 0885-3010 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TUFFC.2020.3032345, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 2 signals is a simple phase shift, which significantly reduces the calculation time [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF]. In addition, this facilitates frequency analysis, to distinguish between different types of cavitation. Indeed, oscillation of the stable cavitation bubbles generates narrow-band harmonic acoustic signals, while implosion of inertial cavitation bubbles creates broadband noise [START_REF] Mclaughlan | A study of bubble activity generated in ex vivo tissue by high intensity focused ultrasound[END_REF]. These approaches do not require the challenging design of preprocessing frequency filter [START_REF] Lyka | Sum-ofharmonics method for improved narrowband and broadband signal quantification during passive monitoring of ultrasound therapies[END_REF], [START_REF] Jensen | Spatiotemporal monitoring of high-intensity focused ultrasound therapy with passive acoustic mapping[END_REF]. However, very few studies have been conducted to develop an adaptive beamformer in the FD for cavitation imaging.

Acoustic imaging in the FD has been widely investigated in the array-processing literature [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Van Veen | Beamforming: a versatile approach to spatial filtering[END_REF]. Theoretical developements are often illustrated with the well-known problem of the direction-of-arrival (DOA) estimation [START_REF] Stoica | Robust Capon beamforming[END_REF], [START_REF] Azimi-Sadjadi | Wideband DOA estimation algorithms for multiple moving sources using unattended acoustic sensors[END_REF], [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], and have been widely applied in aeroacoustics [START_REF] Merino-Martinez | A review of acoustic imaging methods using phased microphone arrays[END_REF]. Such studies usually represent data as a cross-spectral density matrix (CSM), which is the spatial covariance matrix of the received data at the probe expressed in the FD [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Van Veen | Beamforming: a versatile approach to spatial filtering[END_REF], [START_REF] Merino-Martinez | A review of acoustic imaging methods using phased microphone arrays[END_REF]. When appropriately estimated, the CSM allows access to the secondorder statistics of the data, to design optimal beamformers in terms of power (e.g., minimum variance robust adaptive beamforming for improved localization performance even in presence of error models [START_REF] Vorobyov | Principles of minimum variance robust adaptive beamforming design[END_REF]), and also to use further analysis to extract source information, such as in eigenvalue decomposition of the CSM (e.g., the so-called "high-resolution" beamformers like the multiple signal classification (Music) beamformer [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]).

The goal of this paper is not to establish a definitive ranking of beamforming methods, but rather to [START_REF] Escoffre | Therapeutic ultrasound[END_REF] demonstrate the compatibility of adaptive beamformers and the FD; and (2) provide useful concepts of array processing and their advantages within the passive cavitation imaging field. To do so, the paper introduces a robust estimation of the CSM in the context of passive ultrasound imaging of cavitation, and in this way, proposes adaptive beamformers in the FD. Four beamformers are illustrated through the CSM formalism. First, we illustrate how the beamformer design can be derived from the CSM by implementation of the nonadaptive FD-DAS approach as a reference, as this is close to the first FD nonadaptive approach used in [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Bader | Post hoc analysis of passive cavitation imaging for classification of histotripsyinduced liquefaction in vitro[END_REF], and then RCap in the FD (FD-RCap), which has already been implemented in the time domain [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF]. Then, the eigenvalue decomposition (EVD) of the CSM and its benefits are illustrated through the Pisarenko beamformer (FD-Pisa) [START_REF] Stoica | On spatial power spectrum and signal estimation using the pisarenko framework[END_REF]. Finally, taking advantage of the EVD, it is also possible to introduce another cavitation mapping approach that is not based on power estimation: the multiple signal classification beamformer (FD-Music) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF].

All of these methods have been commonly applied in many engineering fields, such as nondestructive testing, seismology, and underwater acoustics, and they are introduced in the context of passive ultrasound imaging. They are first illustrated in simulations. Then, under experimental conditions, we demonstrate their potential for spatio-temporal monitoring of cavitation for real cavitation bubbles.

II. MODEL AND ESTIMATION OF THE CROSS-SPECTRAL DENSITY MATRIX

This section first introduces the signal models in both the time and frequency domains. Then the CSM and its characteristics are defined. Finally, the robust estimation of the CSM is proposed.

A. Signal model

An acoustic source localized in r 0 emits a real signal s(t). This is passively received by a uniform linear array of N sensors, where n is the sensor index (n = 1, ..., N ). Assuming the source transmits spherical waves and the sound speed c is uniform in the propagation medium, the time of flight from the source to each element of the probe is

τ n ( r 0 ) = || r n -r 0 || c . ( 1
)
Assuming there is no attenuation of the wave, the data matrix that contains the received signals is given by [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF]:

y(t) = [s(t -τ 1 ( r 0 )) , ... , s(t -τ N ( r 0 ))] T + b(t), (2)
where . T is the transposition operator, and b(t) is the measurement noise. This b(t) should be white Gaussian noise, with zero-mean, and also spatially white. The data matrix y(t) is of size N × N t , where N t is the number of time samples (which corresponds to an acquisition duration of T rec ). In the FD, for a given frequency f , the received data are represented by an N -length column vector as [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Azimi-Sadjadi | Wideband DOA estimation algorithms for multiple moving sources using unattended acoustic sensors[END_REF]:

Y(f ) = S(f ) a( r 0 , f ) + B(f ), (3)
where S(f ) is the Fourier transform of the source signal s(t), B(f ) is the vector of the Fourier transforms of the noise contributions for each sensor, and a( r 0 , f ) is the propagation vector, such that:

a( r 0 , f ) = [e -i2πf τ1( r0) , ... , e -i2πf τ N ( r0) ] T . (4)
Acoustic imaging aims to build maps on the basis of received signals [START_REF] Khokhlova | Effects of nonlinear propagation, cavitation, and boiling in lesion formation by high intensity focused ultrasound in a gel phantom[END_REF], on which the peak values correspond to the positions of the sources. To do this, spatial filters, called beamformers, have been designed to detect the spatial signature a( r 0 , f ) of the array of an acoustic source placed at r 0 . More precisely, beamforming aims to detect the relative phase shifts between channels due to the delays τ n ( r 0 ).

B. Cross-spectral density matrix analysis

Beamforming design is often based on the spatial covariance matrix of the data received, for three main reasons [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Van Veen | Beamforming: a versatile approach to spatial filtering[END_REF], [START_REF] Mermoz | Imagerie, corrélation et modèles[END_REF]. First, this matrix directly expresses the phase shifts from one element to another of the signals received, and thus the spatial signature of the acoustic source. Secondly, if estimated in the correct way, this gives access to the secondorder statistics of the received data, which allows the optimal beamforming design (in terms of power or variance). Thirdly, the eigenstructure of the data can be used to extract the source information. In the first subsection II-B1, the CSM is defined. Then, subsections II-B2 and II-B3 introduce the notations and concepts derived from the CSM, which are required to understand the four beamformers considered.

0885-3010 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TUFFC.2020.3032345, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 1) Definition: Signals are considered broadband, zeromean, stationary, and ergodic, to deal with their second-order statistics. In the FD, the spatial covariance matrix is usually referred to as the CSM. The CSM is of size N × N , and it is defined as:

M(f ) = E lim Tint→+∞ 1 T int Y(f ) Y * (f ) , (5)
where E . is the expected value, T int is the integration time for the Fourier analysis, and . * is the transpose conjugate operator. The diagonal coefficients are the power auto-spectrum densities of the channels. The other coefficients are the power cross-spectrum densities, which reveal the relative phase-shifts between the channels.

2) Use of the cross-spectral density matrix for power estimation: Passive ultrasound imaging of cavitation was essentially developed through reconstruction of power maps [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF]. Generally speaking, the FD beamformers are spatial filters that consist of applying N complex coefficients, contained in the so-called steering vector h( r, f ), to the observations Y(f ), to compensate for the phase shifts induced by the propagation vector a( r 0 , f ) on the signals received, and then coherently sum them. The mean power spectral density for the random signal h * ( r, f )Y(f ) received at the probe and steered at r is expressed as:

P ( r, f ) = E lim Tint→+∞ 1 T int |h * ( r, f ) Y(f )| 2 . ( 6
)
This expression is looked upon as a power estimate, as commonly used for passive cavitation imaging, although with a statistical point of view, it can be rewritten using the CSM from Equation (5):

P ( r, f ) = h * ( r, f ) M(f ) h( r, f ). (7)
The advantage of this formalism compared to the methods used in [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF] is that power estimate at a given position r is expressed in terms of the second-order statistics of the data, which gives access to optimal beamforming in terms of power or variance.

To distinguish the source from the noise contributions, and thus to design the optimal steering vector, using Equations ( 3) and [START_REF] Pishchalnikov | Cavitation bubble cluster activity in the breakage of kidney stones by lithotripter shock waves[END_REF], Equation ( 7) can be developed as:

P ( r, f ) = |h * ( r, f ) a( r 0 , f )| 2 P s (f ) + h * ( r, f ) M b (f ) h( r, f ), (8)
where P s (f ) is the power spectral density of the source, and M b (f ) is the CSM of the noise. Here, the quantity G r ( r 0 , f ) = |h * ( r, f ) a( r 0 , f )| 2 is the receive spatial gain of the beamforming method. The benefits of this formalism are illustrated later for FD-DAS (III-A) and FD-RCap (III-B).

3) Eigenvalue decomposition of the cross-spectral density matrix: The CSM gives access to very interesting represen-tations of the received data through the EVD of the matrix M(f ), which is defined as [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Merino-Martinez | A review of acoustic imaging methods using phased microphone arrays[END_REF]:

M(f ) = U(f ) V(f ) U * (f ) =   u 1 . . . u N      v 1 . . . 0 . . . 0 . . . v N       oo u * 1 oo . . . oo u * N oo    = N n=1 v n .u n u * n , (9)
where V(f ) and U(f ) are the matrices that contain the eigenvalues (EVs) v n and the column eigenvectors u n , respectively. The EVs are sorted in decreasing order (v 1 > ... > v N ). In a noiseless case, J sources are said to be spatially decorrelated if they are orthogonal in the sensor space, i.e., the (normalized) propagation vector of the source j corresponds to a CSM eigenvector: u j = a( r j , f )/

√

N . In this case, each EV v j is N times the power spectral density P s,j (f ) of the source j, as a signal received on the array is recorded N times: v j = N P s,j (f ).

In the presence of spatially decorrelated background noise, e.g., M b (f ) = σ 2 b I, where σ 2 b (f ) is its power spectral density and I is the identity matrix, the decomposition is modified. The EVD cannot perfectly separate the sources j from the noise for finite-size observations (e.g., number of sensors N is finite). This is uniformly distributed over all EVs. The CSM is decomposed as the sum of a signal space E s that is composed of the first J eigenvectors, and a noise space E b that is composed of the remaining eigenvectors, such that:

M(f ) = E s + E b = J n=1 (N P s,n (f ) + σ 2 b ) u n u * n + N n=J+1 σ 2 b u n u * n .
(10) Each source j is decomposed on one of the subspaces u n u * n that forms the signal space E s , and is orthogonal to the other subspaces, and notably those of the noise space E b . The EVD gives access to new beamformers for passive ultrasound imaging, as illustrated in III-C and III-D.

C. Cross-spectral density matrix estimation

The quality of the cavitation maps greatly depends on how the matrix M(f ) is estimated, and the way that it is further processed. Assuming the received signals are ergodic and stationary, a natural way to estimate M(f ) is to consider K successive acquisitions y k (t) (usually called snapshots) of duration T snap :

M(f ) = 1 K K k=1 1 T snap Y k (f )Y * k (f ) (11)
where Y k (f ) is the Fourier transform of the k th snapshots y k (t). Beamformer performances increase with K [START_REF] Stoica | Robust Capon beamforming[END_REF]. However, T snap has to be high enough to ensure good spectral resolution. To artificially increase the number of realizations K without extensively reducing T snap , a single acquisition

K = T rec (1 -p ov )T snap . (12)

III. BEAMFORMING METHODS

A wide variety of beamformers has been proposed in the array-processing literature [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Van Veen | Beamforming: a versatile approach to spatial filtering[END_REF], [START_REF] Merino-Martinez | A review of acoustic imaging methods using phased microphone arrays[END_REF]. To illustrate the advantages of the proposed formalism, four beamformers that use the second-order statistics of the received data through the CSM are introduced within the framework of passive ultrasound imaging. First, we extend two already existing beamformers for passive cavitation imaging into the proposed formalism: FD-DAS is an extension of [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], but includes the robust estimate of CSM; and FD-RCap corresponds to the beamformer of [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], but transposed into the FD. In this way, we present nonadaptive and adaptive ways, respectively, to use the CSM for optimal steering-vector design. Then, the proposed formalism is used to introduce new beamformers for passive cavitation imaging. In particular, the benefits of exploitation of the eigenstructure of the received data are illustrated for two advanced beamformers. The FD-Pisa is an intermediate approach between FD-DAS and FD-RCap in terms of map quality and complexity. Unlike the other three beamformers, FD-Music is not a power-estimate-based method, although it is well known in the field of array processing for its high-resolution performance. A summary of the beamformer characteristics is given in Table I. Note that although the spectral signature of cavitation is broad, the following analyses consider single-frequency imaging, as is often the case in the literature [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF]. The perspectives for composite multifrequency imaging are further discussed in section VI-4.

A. The nonadaptive beamformer: FD-DAS

The formulation of FD-DAS can be described as a minimization of a cost function, as for DOA estimation in [START_REF] Stoica | Spectral analysis of signals[END_REF]. The design of the steering vector h( r, f ) is here independent of the signals received. First, its coefficients have to minimize the contribution of the noise, and in the basic case it is spatially uncorrelated, such that M b (f ) = σ 2 b I. Secondly, for each pixel in r, the coefficients in h( r, f ) must preserve the signal that would have been emitted with the spatial signature a( r, f ), with a unit gain G r ( r, f ) = 1. From Equation ( 8), the problem can be written as:

min h( r,f ) h * ( r, f )h( r, f ) subject to h * ( r, f )a( r, f ) = 1. (13)
The analytical solution is obtained by the Lagrange multipliers method:

h F D-DAS ( r, f ) = a( r, f ) N . ( 14
)
By introducing Equation ( 14) into the power spectral density estimation of Equation [START_REF] Arvanitis | Combined ultrasound and MR imaging to guide focused ultrasound therapies in the brain[END_REF], and replacing M(f ) by its estimate M(f ), this gives:

PF D-DAS ( r, f ) = a * ( r, f ) M(f ) a( r, f ) N 2 . ( 15
)
This expression provides access to a nonadaptive powerestimate-based beamformer with a robust estimation of the CSM using snapshot averaging. Note that when K = 1 (i.e., without block averaging for the CSM estimation; see Equation ( 11)), this corresponds to [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF]. The nonadaptive approach is very simple and quick, but it is well-known for its poor resolution and for numerous artifacts [START_REF] Haworth | Passive imaging with pulsed ultrasound insonations[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF].

B. The Capon adaptive beamformers: FD-Cap, FD-RCap

Capon adaptive beamforming [START_REF] Stoica | Spectral analysis of signals[END_REF], [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF] aims to optimize the design of the steering vector for each pixel using the information contained in the received data, via the matrix M(f ). In this way, acoustic sources are better resolved and artifact rejection is improved, with respect to the nonadaptive approach.

With a unit gain constraint for the signal that would have been emitted from r (i.e., the same constraint as for FD-DAS), FD-Cap proposes to find the optimal steering vector that minimizes the power of all types of noise (i.e., without assumptions as to the structure of its CSM M b (f ); unlike FD-DAS). Although it is not possible to estimate this directly, in practice M b (f ) is proportional to M(f ) [START_REF] Vorobyov | Principles of minimum variance robust adaptive beamforming design[END_REF]. The problem of minimization of FD-Cap is then written as:

min h( r,f ) h * ( r, f )M(f )h( r, f ) subject to h * ( r, f )a( r, f ) = 1. (16)
This problem has an analytical solution, which can be solved using the Lagrange multipliers method:

h F D-Cap ( r, f ) = M -1 (f ) a( r, f ) a * ( r, f ) M -1 (f ) a( r, f ) . ( 17
)
Introduction of Equation ( 17) into Equation ( 7), and replacing M(f ) by its estimate M(f ), gives:

PF D-Cap ( r, f ) = 1 a * ( r, f ) M-1 (f ) a( r, f ) . ( 18
)
The Capon approach is not used in practice, because it is very unstable in the presence of model errors [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF]. Many approaches for telecommunications or acoustic imaging have been developed to make this method more robust [START_REF] Synnevag | Adaptive beamforming applied to medical ultrasound imaging[END_REF], [START_REF] Vorobyov | Principles of minimum variance robust adaptive beamforming design[END_REF].

0885-3010 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TUFFC.2020.3032345, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control

The FD-RCap [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF], [START_REF] Stoica | Robust Capon beamforming[END_REF], [START_REF] Azimi-Sadjadi | Wideband DOA estimation algorithms for multiple moving sources using unattended acoustic sensors[END_REF] considers that if there are modeling errors (i.e., sound speed, probe calibration), then the true propagation vector a 0 ( r, f ) is included in a small region around the a( r, f ) vector, which is defined by a parameter . The FD-RCap looks for the optimal estimation â0 ( r, f ) of a 0 ( r, f ) from the received data. The optimization problem arises as

min â0( r,f ) â * 0 ( r, f ) M -1 (f ) â0 ( r, f ) subject to ||â 0 ( r, f ) -a( r, f )|| 2 ≤ . (19)
The top term here tries to maximize the power estimation in the sense of Capon (Equation ( 18)) and the bottom term defines the uncertainty region around a( r, f ). This problem has an analytical solution that can also be obtained by the Lagrange multipliers method:

â0 ( r, f ) = a( r, f ) -(I + λ( ) M(f )) -1 a( r, f ), ( 20
)
where I is the identity matrix, and λ( ) is the Lagrange multiplier that is optimally computed for each pixel using a Newton-Raphson numerical method that is controlled by the parameter. Introducing Equation ( 20) into the Capon solution of Equation ( 18), and simplifying [START_REF] Stoica | Robust Capon beamforming[END_REF], this gives

PF D-RCap ( r, f ) = ||â 0 ( r, f )|| 2 /N â * 0 ( r, f ) M-1 (f ) â0 ( r, f ) . (21)
The parameter depends on the received data and regulates the compromise between the robustness and the image quality. A too low value of defines a too small region of uncertainty, whereby FD-RCap tends towards FD-Cap, and thus loses robustness: at the slightest model error, the source signal is interpreted as noise and is rejected. A too high value of defines a too large region of uncertainty: the power of the source is overestimated, and the interference is amplified. The choice of this parameter is generally empirical, with the condition < ||a( r, f )|| 2 , and it is not a big issue as it can be found by a little trial and error [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF]. It should be noted that the inversion of M(f ) implies the matrix is well conditioned, with a sufficient condition of K > N [START_REF] Stoica | Robust Capon beamforming[END_REF], [START_REF] Stoica | Spectral analysis of signals[END_REF].

C. The Pisarenko class beamformer: FD-Pisa

An alternative to the pixel-wise optimal design of the steering vector is to artificially change the relative levels between the sources prior to steering. This is possible by taking advantage of the eigenstructure of M(f ) and applying nonlinear operations to its EVs. This approach is known as Pisarenko class beamforming [START_REF] Stoica | On spatial power spectrum and signal estimation using the pisarenko framework[END_REF], [START_REF] Pisarenko | On the estimation of spectra by means of non-linear functions of the covariance matrix[END_REF], [START_REF] Dougherty | Functional Beamforming[END_REF], or FD-Pisa. Namely, to sharpen the peaks that reflect the source positions (with respect to FD-DAS), the p th root is applied to the EVs of M(f ):

M 1/p (f ) = N n=1 v 1/p n .u n u * n . (22)
In this way, the gap between the source levels is compressed: the sidelobes generated by strong sources will interfere less with the lower source level when beamforming. The steering vector is the same as for FD-DAS, with a previous normalization such that no bias is introduced when using nonlinear operations on it,

h F D-P isa ( r, f ) = a( r, f ) √ N . (23)
After steering the compressed CSM M 1/p (f ), the p power is used on the resulting value to restablish the original dimension of an acoustic power:

PF D-P isa ( r, f ) = 1 N a * ( r, f ) √ N M1/p (f ) a( r, f ) √ N p ( 24
)
where the 1/N term scales the N factor within the EVD (see Equation ( 10) and [START_REF] Stoica | On spatial power spectrum and signal estimation using the pisarenko framework[END_REF]). The parameter p has to be empirically chosen: p = 1 corresponds to FD-DAS, and p > 1 is interesting to sharpen the peak values. The performances in terms of localization usually increase with p when uncorrelated sources are considered [START_REF] Stoica | On spatial power spectrum and signal estimation using the pisarenko framework[END_REF], [START_REF] Dougherty | Functional Beamforming[END_REF], [START_REF] Baron | Numerical and experimental assessment of functional beamforming for source quantification[END_REF].

D. The multiple signal classification beamformer: FD-Music

The FD-Music [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], [START_REF] Bienvenu | Optimality of high resolution array processing using the eigensystem approach[END_REF] is not power-estimation based, but it is well known as a high resolution beamformer. FD-Music localizes the sources using the contributions of the background noise, rather than the source signal directly. More precisely, FD-Music exploits the orthogonality between the signal subspace and the noise subspace obtained by EVD decomposition of M(f ), as described in Equation [START_REF] Khokhlova | Histotripsy methods in mechanical disintegration of tissue: towards clinical applications[END_REF].

For each pixel r of the map, FD-Music projects the steering vector h F D-M usic ( r, f ) = a( r j , f ) to the noise subspace: this is the sum of the scalar products of h F D-M usic ( r, f ) with each vector of E b . The steering vector is chosen proportional to that of the FD-DAS in Equation ( 14), without the 1/N factor, as the FD-Music does not seek to quantify the power. Two cases are then possible:

• The pixel in r coincides with the position r j of the source j. Then the pointing vector h F D-M usic ( r, f ) = a( r j , f ) is a eigenvector of the signal space E s . By construction of the EVD, it is therefore necessarily orthogonal to each eigenvector of the noise space E b , and therefore the sum of these scalar products is also zero:

N n=J+1 a * ( r j , f ).u n = 0. (25)
• The pixel in r = r j does not coincide with any of the J sources. So the steering vector h F D-M usic ( r, f ) = a( r, f ) is not an eigenvector of E s . There is, therefore, at least one eigenvector of the noise space for which its projection is not zero. So by necessity, one of the scalar products a * ( r, f ).u n (∀ n ∈ J + 1 ; N ) is not zero. Thus the sum of the scalar products is not zero:

N n=J+1 a * ( r, f ).u n = 0, if r = r j . ( 26
)
Taking advantage of the properties defined by Equations ( 25) and ( 26), the FD-Music map is defined by: where K s is an estimation of the signal space size, as the number of sources J is unknown in practice. The choice of size K s for E s is often guided by analysis of the decrease in EVs [START_REF] Marcos | Les méthodes à haute résolution[END_REF], [START_REF] Wong | On information theoretic criteria for determining the number of signals in high resolution array processing[END_REF]. The maxima of PF D-M usic ( r, f ) are reached when its denominator tends toward zero (equivalent to Equation ( 25)), and therefore indicates the source position. Indeed, in practice, on a source pixel, the denominator does not cancel itself out strictly, but tends toward 0, giving peak values on the map. For the sake of uniformity in localization map notation, it can be noted PF D-M usic ( r, f ), although it has to be kept in mind that it is not a power map.

PF D-M usic ( r, f ) = 1 N n=Ks+1 |a * ( r, f ).u n | 2 , (27)

IV. MATERIALS

To illustrate the relevance of this exploitation of the robust estimation of the CSM and the resulting FD adaptive beamformers within the passive cavitation imaging framework, single and multiple source configurations were carried out on simulations, and an experimental set-up was developed to create cavitation bubbles in a water tank. We first describe this set-up, and then the simulations using similar parameters. The performance metrics are also detailed.

A. Experimental set-up 1) High-intensity focused ultrasound settings: The experimental set-up is represented in Fig. 1. Cavitation bubble clouds were created using a 1-MHz single spherical transducer focused at 8 cm, with an active diameter of 10 cm, and a circular 4.5-cm diameter aperture in the center (Imasonic, France). This was immersed in a tank of filtered water (72 cm×25 cm×25 cm). The theoretical focal volume at -3 dB is ellipsoidal in shape, 8.3-mm long for the x-axis, and 1.1-mm in diameter. The HIFU transducer is connected to a power amplifier (GN1000; PR ÂNA, France). This can transmit bursts centered at 1 MHz, with each lasting 20 ms, and pulsed at a rate of 4 Hz. The amplitude delivered is adjustable, although it remains constant during each shot. The amplitudes can range from 10 mV to 150 mV, which generates peak negative pressures at the focus of the transducer of 0.72 MPa to 10.8 MPa. Here, we focus on the data at 7.2 MPa.

2) Passive ultrasound imaging: For passive ultrasound imaging, a linear ultrasound imaging probe (L7-4; Verasonics, USA) was used. It is composed of N = 128 elements with uniform spacing of 298 µm, for a total width of 38 mm. Its center frequency is 5 MHz and its bandwidth at -6 dB is 3.5 -7.5 MHz. It is connected to an ultrasound scanner (Vantage 256; Verasonics, USA). The probe makes continuous passive acquisitions with a duration of 24 ms at a sampling frequency of f s = 20.8 MHz. Although the therapy signal oscillates at a fundamental frequency of f HIF U = 1 MHz, outside the bandwidth of the imaging probe the images are generally reconstructed by exploiting higher frequencies than for cavitation, which are within the bandwidth of the imaging probe, to ensure better resolution [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Jensen | Spatiotemporal monitoring of high-intensity focused ultrasound therapy with passive acoustic mapping[END_REF], [START_REF] Bader | Post hoc analysis of passive cavitation imaging for classification of histotripsyinduced liquefaction in vitro[END_REF], [START_REF] Choi | Non-invasive and real-time passive acoustic mapping of ultrasoundmediated drug delivery[END_REF]. The imaging probe is placed 6 cm above the focus of the HIFU transducer. The lateral axis of the imaging probe is parallel to the x-axis, and points toward the HIFU transducer. The ultrasound images are reconstructed in the xz-plane.

3) Optical imaging: Ombroscopic imaging of the cavitation phenomenon was obtained using a fast optical camera (Phantom v12.1; Vision Research Inc., USA) that recorded at the rate of 10 000 frames per second, and an exposure time of 99 µs, over 24 ms. The images obtained are in the same xz-plane as the ultrasound images. The field of view is 1.5cm wide and 1.1-cm high (800×600 pixels), and contains the HIFU focus.

4) Cavitation generation: The objective was to compare optical images with the reconstructed passive ultrasound maps of cavitation events that consisted of a few bubbles nucleated with the HIFU bursts. Optical images are limited to a very small region compared to the broad area where the cavitation might be nucleated and within which this might move in the water tank (not necessarily near the focus). To ensure comparison with optical images, i.e., to trap the cavitation within the field of view of the camera, a polyester wire with asperities that provided nucleation sites along the wire was placed at the HIFU focus (Fig. 1). In this way, the cavitation was easier to image optically, as it occurs mainly in the neighborhood of the intersection of the wire and the HIFU focus. Also, the cavitation threshold was stabilized, and essentially depended on the gas concentration (estimated from the oxygen concentration). The cavitation threshold was measured using spectrogram analysis of the signal passively recorded with the central element of the ultrasound probe, as was done with a hydrophone in [START_REF] Mclaughlan | A study of bubble activity generated in ex vivo tissue by high intensity focused ultrasound[END_REF], [START_REF] Neppiras | Measurement of acoustic cavitation[END_REF]. Two oxygen concentrations were used: 2.5 mg/L and 6.0 mg/L, which corresponded to the measured cavitation thresholds of 5 MPa and below 0.7 MPa, respectively. This experimental set-up allowed comparison of the optical and ultrasound imaging modalities of a wide variety of cavitation events. In this paper, we focus on the cavitation events that consisted of the few nucleated bubbles that stayed at the wire surface. In a second configuration, a second wire (wire #2) was placed 3-mm downstream from the first one (wire #1); see Fig. 1. This wire #2 was still located within the -3 dB focal volume, although at a point where the HIFU field was less intense, and thus this produced a second weaker source in comparison to the first source created on wire #1, which constitutes a more challenging scenario to be imaged, with artifacts of the strong source spreading on the weaker one. In the results section, the presence of inertial cavitation through the whole HIFU burst was confirmed using spectrogram analysis, as in [START_REF] Mclaughlan | A study of bubble activity generated in ex vivo tissue by high intensity focused ultrasound[END_REF], [START_REF] Neppiras | Measurement of acoustic cavitation[END_REF].

5) Synchronization and registration of images: The control system for the HIFU transducer generates a synchronization signal 1 ms before the start of the HIFU shot. This signal is received by the ultrasound scanner and launches passive acquisition for a period of 24 ms. The same applies to the optical camera. The ultrasound propagation and device activation times are less than 100 µs. Thus, cavitation is imaged before, during, and after the entire shot.

The ultrasound and optical images are spatially registered using two nodes on the wire as markers. Before the HIFU shot, these are localized by active ultrasound imaging. Their coordinates are then assigned to the nodes visible in the optical image. In this way, the optical images are described via the x and z axes, as with the ultrasound images. Comparison is then possible with maps obtained by passive ultrasound imaging.

6) Reconstruction parameters: For passive ultrasound imaging, spatio-temporal tracking of the cavitation is ensured by dividing the total acquisition of 24 ms into 120 subacquisitions of T rec = 200 µs, with each of these providing a localization map. For each map, snapshot averaging is achieved with K = 130 realizations and T snap 15 µs; i.e., with an overlap of p ov = 90 % (as the parameters that provided the best results in the simulation V-A3).

B. Simulation set-up

Simulations are conducted with probe parameters similar to those described in section IV-A2, except for the frequency response of the probe, which is not simulated. Four simulated configurations are considered, as summarized in Figure 2. All of the configurations consider sources that emit random Gaussian noise to mimic cavitation noise, as proposed in [START_REF] Lyka | Sum-ofharmonics method for improved narrowband and broadband signal quantification during passive monitoring of ultrasound therapies[END_REF], through the whole acquisition, which lasts for T rec = 200 µs. To mimic the measurement noise, white Gaussian noise is systematically added, independently to each sensor, with a varying signal-to-noise ratio (SNR).

Configuration 1 is designed to investigate the evolution of the common image quality metrics (i.e., lateral and axial main lobe widths, contrast, and error position, as detailed in section IV-C) with respect to the SNR varying between -15 dB and +15dB. To do so, a single source is placed at r 1 = [5, 60] mm.

Configurations 2 and 3 are carried out to illustrate the superiority of adaptive beamformers to resolve two closely spaced sources for the lateral and axial directions, respectively. To do so, a first source is placed at r 1 = [5, 60] mm, and a second one is added, for configuration 2 with 0.4 mm lateral spacing (about one wavelength λ = 0.31 mm), and for configuration 3 with 3.5 mm axial spacing (about 11 wavelengths). In both configurations, the two sources have the same transmitting amplitude, with a SNR of 10 dB.

Configuration 4 is simulated to highlight the influence of CSM averaging on cavitation maps; namely, in a complex scenario where three sources coexist. Moreover, configuration 4 investigates the interest in the use of FD adaptive beamformers to reduce the side lobes of a strong source that hides two weaker sources. To do so, three sources are simulated, with the first high-intensity source at r 1 = [5, 60] mm. The second source is placed in the side lobes generated by the FD-DAS at r 2 = [6, 63] mm, with a power of 10 dB less than the first one. The third source is placed close to the second source, at r 3 = [6.7, 63] mm (with an approximately two-wavelength spacing), and also with 10 dB less power than the first one. The SNR is 10 dB.

C. Beamformer performance metrics

The benefits of the CSM formalism are assessed through qualitative and quantitative comparisons of the imaging methods. Simple metrics commonly used in the context of ultrasound imaging are considered, to characterize the resolution and contrast of the images obtained for both the simulated and experimental data.

The full-width at half-maximum of the main lobe is measured in millimeters, as R ax in the axial dimension and R lat in the lateral dimension. This is a simple metric measured on single-source configuration (or isolated source), as commonly used to reflect the resolution of two sources by a beamformer.

The gap measurement between two close sources represents a more relevant metric to quantify the resolution performance of a beamformer, although this is more complex. It is computed as (before log-compression of the map values):

Gap = 10 log 10   max P ( r 1 ), P ( r 2 ) P ( r1+ r2 2 )   . (28)
The cavitation-to-artifact ratio (CAR) is defined as the ratio of the maximum value from the map over the averaged pixel values in an artifact region (before log-compression of the map values). The CAR measurement is a value with no physical meaning, but that aims to provide a relative comparison of the rejection of artifacts by each beamformer (e.g. side lobes, main lobe extension under -3 dB). The CAR reflects the contrast, and thus the unmasking of weak sources by the beamformer in the presence of stronger sources. The position and the extent of the artifact region do not favour a beamformer (except maybe in a very pathological case, as very far from the target, or including the main lobe above -3 dB).

The position error Err pos is the distance in millimeters between the actual position of a simulated point source and the maximum pixel value in the map (for a single source scenario).

The running time for the reconstruction is measured to reflect the complexities of the beamformers. All of the reconstructions are carried out in Matlab (R2016b, The MathWorks, Natick, MA, USA), on the same computer, which has 32 GB memory and an Intel(R) Core(TM) i7-6600 CPU at 2.60 GHz.

V. RESULTS

We first present the data for the simulations, to introduce the performances of the four beamformers and the influence of the robust CSM estimation for the different point source configurations (Fig. 2). Then, the experimental data for small clouds are described. According to the guidelines given in the literature, the adaptive parameters for FD-RCap [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF] and FD-Pisa [START_REF] Dougherty | Functional Beamforming[END_REF] are empirically chosen to satisfy the compromise between resolution and contrast (as indicated for each configuration), as they are not easily physically interpretable. For FD-Music, the choice of K s is more guided (as it corresponds to the number of uncorrelated sources, which is easy to know in simulation, and also due to the EV plots, as described in configuration 4). It should be noted that all of the maps do not represent the same quantity. On the one side, the power-estimate-based beamformers (i.e., FD-DAS, FD-RCap, FD-Pisa) are compared to each other. On the other side, the noise space projection-based beamformer FD-Music requires separate comment (namely, when comparing metrics on a dB scale, as for CAR or Gap).

A. Simulation results

The four configurations (Fig. 2) are analyzed with the same reconstruction frequency at 5 MHz (i.e., the center frequency of the probe used for the experiments), and with similar pixel spacing of 0.02 mm for the lateral dimension and 0.1 mm for the axial dimension.

1) Simulation, configuration 1: influence of the SNR: Configuration 1 investigates on a single source the evolution of R lat , R ax , CAR, and Err pos with respect to the SNR, for the four beamformers FD-DAS, FD-RCap with = 1 and = 10, FD-Pisa with p = 50 and p = 100, and FD-Music with K s = 1 and K s = 3. The parameters for the CSM estimation are the same through all of the results for this configuration: K = 130, with p ov = 90%, which corresponds to the optimal parameters described in V-A3. The image quality metrics are shown in Figure 3.

The error positions on Fig. 3d are negligible as they are equal to or smaller than 0.2 mm (i.e., mainly due to one or two pixel errors in the axial dimension, which is consistent with the observations in [START_REF] Abadi | Frequency-sum beamforming for passive cavitation imaging[END_REF]) for SNRs above -10 dB, except for FD-RCap with = 1. Indeed, = 1 corresponds to a relatively small tolerance value for error models. Thus, in the high-noise scenario (SNR below 0 dB), increasing up to 10 reduces the error position to 0.2 mm for SNRs above -10 dB.

In the very-high-noise scenario for SNR = dB, the maps are not reliable for FD-RCap and FD-Music any more, which leads to high position errors (i.e., greater than or equal to 0.7 mm).

For the main lobe size (Fig. 3a,b), the axial width is larger than the lateral width (i.e., about 10-fold), as is commonly seen in passive cavitation imaging [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF]. For all beamformers, both of the widths evolved with very similar trends with respect to the SNR. For SNRs above -5 dB, the three adaptive beamformers provide thinner main lobe widths (regardless of the adaptive parameter values), which are almost halved, at least. For FD-RCap, as reflects the tolerance to uncertainties, a larger value of is required for the low SNR, while a smaller value of in the high-noise scenario allows very thin main lobes to be assessed. However, FD-Pisa appears stable even for high p values, and also FD-Music, with slightly overestimated K s , whatever the SNRs. FD-Pisa (for both p values) provides intermediate performances between FD-DAS and FD-RCap, with the latter appropriately tuned with . Very similar observations are established for CAR in Fig. 3c (the artifact region is set as in Figure 2, to enclose the FD-DAS upstream side lobes).

This joint analysis on four image quality metrics in the single-source scenario leads to the conclusion that following the guidelines from the literature, the appropriate tuning of the adaptive parameters is not a difficult task for better image quality metrics than for FD-DAS, even for the high-noise scenario (SNR = -5 dB). The two following subsections illustrate the consequences of improved R ax , R lat and CAR, in terms of enhancement of source localization.

2) Simulation, configurations 2 and 3: closely spaced sources: Configurations 2 and 3 are designed as a proof of concept,with the aim to illustrate the superiority of FD adaptive beamformers for the resolution of two closely spaced sources. The four beamformers are considered: FD-DAS, FD-RCap with = 1 (i.e., low-noise scenario), FD-Pisa with p = 50, and FD-Music with K s = 2 (i.e., two sources). The parameters for the CSM estimation are the same for both configurations: K = 130 with p ov = 90%, which corresponds to the optimal parameters described in V-A3. Lateral slices (at z = 60 mm) for configuration 2 and axial slices (at x = 5 mm) for configuration 3, are shown in Fig. 4a,b, respectively.

In both configurations, FD-DAS does not separate the two closely spaced sources, but the adaptive beamformers do. Using FD-Pisa, a slight gap appears between the two sources: lateral Gap = 1.0 dB, and axial Gap = 0.4 dB. FD-RCap provides the best performance among these power-based beamformers, with lateral Gap = 6.1 dB and axial Gap = 3.0 dB. Also, FD-Music clearly detects both sources (as the map units are not the same as the power-based beamformers, we cannot indicate any resolution superiority for FD-Music). As commonly observed in passive cavitation imaging, the axial resolution is worse than the lateral resolution (with source spacing of 0.4 mm on the x axis, and 3.5 mm on the z axis). For FD-Pisa and FD-RCap, even if both of these sources nominally transmit the same power, a small difference in the power map estimates is observed for a single reconstruction frequency. Indeed, the assumption of white noise is true from a statistical point of view; thus at a single frequency, from one realization to the other, fluctuations of the power spectral density of each source up to 2 dB can be observed. Also, this is potentially due to interference between sources. It should be noted that these configurations are challenging cases that are chosen to highlight the potential of FD adaptive beamformers to resolve two closely spaced sources where the nonadaptive scenario fails.

3) Simulation, configuration 4: influence of estimation quality of the CSM: Configuration 4 is designed to highlight the necessity for CSM averaging when using FD adaptive beamformers; namely in a multiple-source scenario. The data for configuration 4 are presented in Figure 5, to provide a strategy to select the parameters for CSM averaging and also to illustrate their influence on the maps. The four beamformers are considered: FD-DAS, FD-RCap with = 5 (more sources, more noise, so increased compared to the three previous configurations), FD-Pisa with p = 50, and FD-Music with K s = 3 (i.e., three uncorrelated sources).

The strategy to choose the parameters for the CSM averaging is as follows. First, the recording duration T rec is usually imposed by the hardware capacity or the duration of the observed phenomenon. Moreover, K is constrained by the condition of inversion of FD-RCap (i.e., K > N ; see III-B). So, Fig. 5 first compares the use, or not, of snapshot averaging (i.e., K = 1 or K = 130), and then it illustrates the effects of p ov (directly linked to T snap via Equation ( 12)). Thus, four cases are considered: K = 1, i.e., without snapshot averaging (Fig. 5a-d); K = 130 without overlap, which leads to T snap = 1.5 µs (Fig. 5e-h); K = 130 with p ov = 50%, which leads to T snap = 3.1 µs (Fig. 5i-l); and K = 130 with p ov = 90%, which leads to T snap = 15.4 µs (Fig. 5m-p). a) Effects of cross-spectral density matrix averaging on cavitation maps: Without snapshot averaging (K = 1), as in Fig. 5a, FD-DAS detects the strong source, although the sidelobes it generates hide both of the weak sources. In Fig. 5b, FD-RCap does not detect any source, and the map is completely unstable, as M(f ) is ill-conditioned. In Fig. 5c, FD-Pisa only detects the strong source, and misses both of the weak sources, while in Fig. 5d, FD-Music barely detects the strong source, and it has many artifacts. The second row in Fig. 5 corresponds to a number of snapshots where K > N (without overlap), to ensure the inversion of M(f ) (see III-B). In Fig. 5f, FD-RCap is stabilized. Both of the weak sources emerge for FD-RCap and FD-Pisa (Fig. 5f,g,respectively). However, as snap is small, the spectral analysis is degraded, so the main lobe of the strong source is affected (Fig. 5e,g,h), in comparison to the first row (Fig. 5a,c,d). For the third row of Fig. 5, a snapshot overlap is introduced to maintain K > N , while increasing T snap with p ov = 50% (Fig. 5i-l): the main lobes are sharpened in the axial direction, with respect to the case where p ov = 0% (Fig. 5e-h). With p ov = 90% (Fig. 5mp), this trend is enhanced. In Fig. 5m, FD-DAS still fails to As FD-Music is not a power-based beamformer, it a specific y-axis. detect both of the weak sources, although all of the sources are detected with the three advanced adaptive beamformers: FD-RCap, FD-Pisa, and FD-Music (Fig. 5n-p). The following analyses compare the beamformers for enhanced weak source detection, with a focus on the set of parameters that gives the best performances (i.e., three detected sources for the three FD adaptive beamformers), such that K > N , p ov = 90 % and T snap = 15.4 µs, which correspond to Fig. 5m-p.

b) Power-estimate-based beamformer comparison:

The performance is in line with the previous configurations 1, 2, and 3. FD-RCap offers the best performance in terms of resolution (i.e., the thinnest main lobes) and contrast (i.e., the best artifact rejection) (Fig. 5n). However, FD-RCap suffers from the high computational cost (about 104.7 s versus 1.3 s for FD-DAS; to compute one map composed of 35,140 pixels). Again, FD-Pisa (Fig. 5o) offers intermediate performance between FD-DAS (Fig. 5m) and FD-RCap (Fig. 5n), although with an execution time of 1.4 s, so similar to FD-DAS. The two weak sources are not distinguished with FD-RCap either, although they are clearly revealed in comparison to FD-DAS. c) Influence of K s on FD-Music: The four maps obtained with FD-Music are shown in Fig. 6, as: K s = 1 (Fig. 6d) (with underestimated size of the signal subspace E s ); K s = 3 (Fig. 6e) (with correctly estimated size of E s ); and K s = 5 (Fig. 6f) or K s = 10 (Fig. 6g) (with overestimated size of E s ). The EVs are shown in Fig. 6h. When the real size of the signal subspace is underestimated, some of the sources are missing, as in Fig. 6d. When K s is the same as the effective size of the signal subspace, all of the sources are correctly detected (Fig. 6e) with the best performance metrics: the Gap between the two weak sources is 8 dB. However, when the effective size of the signal subspace is overestimated (K s = 5 Fig. 6f, or K s = 10 (Fig. 6g), the maps are degraded (i.e., with artifacts and main lobe extensions), with decreased localization performance (Gap between the two weak sources of 5.7 dB and 4.6 dB, respectively).

When the real size of the signal subspace is not known, this can be estimated considering the EV plot (Fig. 6h) [START_REF] Wong | On information theoretic criteria for determining the number of signals in high resolution array processing[END_REF]. A simple trial-and-error approach can be as follows. The sources are decomposed on the first EVs, then they are followed by a plateau for the EVs due to a spatially decorrelated noise; see Equation [START_REF] Khokhlova | Histotripsy methods in mechanical disintegration of tissue: towards clinical applications[END_REF] or [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF]. Thus, the plateaus (e.g., EV 2 and 3) or bends (e.g., EV 2 or 5) observed on the first eigenvalues before a broad plateau or a slow decrease in EV, are potentially the values for which the signal space and the noise space separate. As underestimation of the signal space can miss the source detection (Fig. 6d), it is better to slightly overestimate this space (Fig. 6f), although not by too much, to avoid loss of resolution or contrast (Fig. 6g). This approach is used to tune K s for the experimental results. More sophisticated methods can be investigated, although they are application dependent (e.g., see the exhaustive comparison of automated singular value thresholding for ultrasound Doppler imaging in [START_REF] Baranger | Adaptive spatiotemporal SVD clutter filtering for ultrafast doppler imaging using similarity of spatial singular vectors[END_REF]). Finally, whatever K s is, the FD-Music running time is 1.4 s, and so similar to that of FD-DAS. d) Joint qualitative comparison of the four beamformers: Figure 6 qualitatively compares the beamformers for enhanced weak source detection. The maps from Fig. 5m-p are shown in Fig. 6a-c, e with the specific dynamic ranges, from the lowest intensity in their own background, up to 0 dB. The displaying of these maps with different dynamic ranges is suitable to compare the localizations of the different sources. Indeed, it can be seen that FD-Music (Fig. 6e) provides a very smoothed background with respect to FD-RCap (Fig. 6b) and FD-Pisa (Fig. 6c). In this way, FD-Music is expected to be an excellent tool to distinguish sources from artifacts.

B. Experimental results

1) Single wire configuration: Fig. 7 shows the maps obtained with FD-DAS (Fig. 7a), FD-RCap (Fig. 7b) ( = 10), FD-Pisa (Fig. 7c) (p = 30), and FD-Music (Fig. 7d) (K s = 3), with the overlaying optical images. The configuration shown in Fig. 7 to illustrate the performances of the beamformers for the real cavitation events corresponds to an oxygen concentration of 2.5 mg/L, a measured cavitation threshold of 5 MPa (i.e., a little smaller than the peak-negative pressure at the (a,b,c,e) correspond to the maps of Fig. 5m,n,o,p, respectively, with different dynamic ranges for the displays. The maps are normalized and displayed with their own full dynamic ranges: from the lowest intensity in their own background, up to 0 dB. focus of 7.2 MPa) in the middle of the shot at t = 10 ms. Only a few bubbles are detected with the optical camera, and they are all located on the wire, in the small area that is 1-mm long and is circled by the white/ red dashed lines in Fig. 7. A reconstruction frequency of 4.7 MHz was chosen, close to the central frequency of the probe, to ensure that only inertial cavitation activity is mapped, and not the harmonic contribution that results from nonlinear distortion of the acoustic pulse due to the transducer itself, and potentially scattered by the wire [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF]. It should be noted that during the shot, the wire was pushed away from the transducer a little, with the maximum axial displacement of 0.3 mm. This wire displacement was not seen in the absence of cavitation, and results from the radiation pressure on the bubbles. For all of the maps, the beamforming methods succeed in localizing the cavitation activity with a single main lobe that coincides with the area where the bubbles are detected by the optical camera. This remains true over the entire shot, with the main lobe following the slight displacement of the wire. The characteristic hour-glass shape of FD-DAS (Fig. 7a), which has no physical meaning, is absent from the adaptive approach maps. FD-RCap (Fig. 7b) shows a smooth, oval lobe, as already seen in the literature for its time-domain version [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF]. This is also the case for FD-Pisa (Fig. 7c), while the main lobe of FD-Music (Fig. 7d) is less homogeneous and appears to agglomerate several small lobes. The localization performance statistics calculated on the 100 maps reconstructed during the 20 ms of shooting are shown in Fig. 7e for axial resolution R ax , Fig. 7f for lateral resolution R lat , and Fig. 7g for CAR. For the power-estimatebased beamformers, while the lateral resolution R lat is almost identical from one method to another, the axial resolution R ax is clearly better for the adaptive approaches, with an average of 7.5 mm for FD-DAS, versus only 3.7 mm for FD-RCap and 4.6 mm for FD-Pisa. The CAR (Fig. 7g) is also greatly improved with the adaptive methods. Once again, FD-Music cannot be directly quantitatively compared to the powerestimate-based beamformers, as it does not represent the same quantities. Nevertheless, the qualitative aspect for FD-Music described above is reflected in the quantitative metrics with relative considerations between the metrics. Indeed, the R lat metrics for FD-Music are very close to those of FD-DAS (less than 0.08 mm difference, for mean values), while R ax and CAR are clearly enhanced for FD-Music compared to FD-DAS, with an enhancement of more than 5 mm and 5 dB for the mean values, respectively.

2) Dual wire configuration: In the second configuration, which included the two wires, the oxygen concentration was increased to 6.2 mg/L, which reduces the cavitation threshold to below 0.7 MPa while keeping the excitation HIFU signal unchanged. This enabled the triggering of a second cavitation spot on the second wire that was placed downstream from the focus of the HIFU transducer. The maximum displacements of the wires due to the radiation pressure on the cavitation bubbles were 1.6 mm for the upstream wire and 0.6 mm for the downstream wire. These values are higher than for the first configuration, which is attributed to more sustained cavitation activity detected on the wires in optical images in the circled regions in Fig. 8a-d. As anticipated, the cavitation activity observed in the optical images appears to be more intense for the upstream wire, which was placed closer to the focus of the HIFU transducer than the downstream wire. Fig. 8 shows the representative maps of two cavitation sources on the two separate wires for FD-DAS (Fig. 8a), FD-RCap (Fig. 8b) ( = 10), FD-Pisa (Fig. 8c) (p = 30), and FD-Music (Fig. 8d) (K s = 4) during the first part of the shot, at t = 4 ms. The reconstruction frequency of 4.7 MHz was chosen. For all of the maps, both cavitating regions correspond to the two main lobes of the maps, which are centered on each wire at about 59 mm in depth, where the cavitation activity was detected optically. For the power-estimate-based beamformers, the source on the downstream wire is weaker. For FD-DAS, its intensity hardly reaches the level in the characteristic hourglass-shaped artifact from the strong source: this might have been hidden in another geometric configuration, such as in Fig. 5. For the adaptive beamformers FD-RCap and FD-Pisa, the better resolution allows the detection of the weak source more clearly. FD-Music gives a very different map, with two clear main lobes at the positions of the cavitating regions, whereby both emerge very clearly from a very smooth and homogeneous background, which is similar to what was noted in the simulation data. Again, the lobes appear to have a more complex internal structure than for the other beamformers.

3) Spatio-temporal monitoring: The projected maximum map values on the x-axis versus time are shown for each algorithm in Fig. 8e-h. This representation is interesting for measuring the ability of each method to bring out the main lobe of a source with respect to the background map, and to ensure the spatio-temporal monitoring of the cavitation with respect to the x-axis (similar to the spatio-temporal tracking in [START_REF] Choi | Non-invasive and real-time passive acoustic mapping of ultrasoundmediated drug delivery[END_REF]). The evolution of the local maxima of the maps was correlated with the displacement of the wires on which the cavitation bubbles were created (Fig. 8e-h, two dashed curved lines). In all of the cases, two trails are observed. The rightwing trace (at around x = 4 mm) was very energetic for the first half of the shot, and corresponds to the main lobe centered on the upstream wire (wire #1). The left-wing trace (at around x = 1 mm) was less energetic, and corresponds to the lobe centered on the downstream wire (wire #2). However, for FD-DAS (Fig. 8e) up to about t = 10 ms, this left-wing trace hardly emerges from artifacts from the more energetic right trace, while the adaptive methods FD-RCap (Fig. 8f) and FD-Pisa (Fig. 8g) efficiently reject these artifacts, and show clear identification of the left trace, thus enhancing the spatio-temporal monitoring of the cavitation activity. For FD-Music (Fig. 8h), the cavitation cloud on wire #2 is raised to almost the same level as the cavitation on wire #1. Again, this observation is consistent with the lack of preservation of the relative amplitudes of the sources using FD-Music, and makes it possible to highlight sources with little power. The absence of cavitation is detected on wire #1 beyond 10 ms for all of the systems here (Fig. 8e-h). The growth of the cavitation bubbles involves the recruitment of gas pockets trapped on the wire surface, as well as gas dissolved in the surrounding medium. As soon as larger-sized bubbles are obtained, they are pushed away by the radiation force of the acoustic wave, which degasses the wire and its surroundings, so that cavitation on the wire would require a higher magnitude of negative pressure. Finally, particularly with the adaptive beamformers, passive ultrasound imaging allows good spatiotemporal monitoring of complex scenarios of cavitation events. This successfully detects the extinction of the first cavitation spot around t = 10 ms, and the persistence of the second cavitation spot through the whole burst. The specificity of FD-Music can also be emphasized, as it unambiguously highlights the weak source from the begining to the end of the burst.

VI. DISCUSSION

In this paper, the concept of robust CSM estimation is introduced and its benefits for passive cavitation imaging are illustrated for four beamformers, through both simulation and experiments. In this way, we demonstrate that adaptive beamformers, and the enhanced image quality they provide, are compatible with low computational cost and frequency selectivity provided by FD approaches.

1) Robust cross-spectral density matrix estimation is required: Fig. 5 shows the need for robust CSM estimation to develop adaptive beamformers. As beamformer theory is derived from the theoretical CSM M(f ) expressed as an expected value E . , it is not surprising that averaging observations highly impacts upon the performance of the beamformers: this ensures the required matrix inversion for FD-RCap and a reliable EVD for the CSM for FD-Pisa and FD-Music. In this way, the maps are stabilized and the sources are correctly detected (Fig. 5). This result is well known in the array-processing field, and in particular for DOA estimation [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF], [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]. Averaging with snapshots is often referred to as the natural CSM estimation, although alternative or joint strategies can be considered, such as diagonal loading [START_REF] Abadi | Frequency-sum beamforming for passive cavitation imaging[END_REF] and spatial smoothing [START_REF] Synnevag | Adaptive beamforming applied to medical ultrasound imaging[END_REF]. Note that for the robust Capon beamformer in the time domain [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], the spatial covariance matrix is estimated from several observations, each of which corresponding to one acquisition time sample.

2) Power-estimate-based beamformers comparison: The CSM is the basis from which many beamformers can be derived [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF], [START_REF] Van Veen | Beamforming: a versatile approach to spatial filtering[END_REF], [START_REF] Merino-Martinez | A review of acoustic imaging methods using phased microphone arrays[END_REF]. Three beamformers are chosen here to illustrate the benefits of robust CSM estimation, and a preliminary comparison can be drawn between the powerestimate-based ones. The nonadaptive FD-DAS approach offers equivalent performance to [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF] (in the special case of FD-DAS for which the CSM estimation is based on a single observation). FD-DAS is used here as a nonadaptive reference method, as useful for comparing more advanced approaches. The FD-RCap method improves source localization in terms of resolution and contrast, as already demonstrated in its timedomain version (TD-RCap) [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF]. FD-RCap is therefore an excellent candidate for high-resolution localization of cavitation. Its main disadvantage is the computational overload that is required for the pixel-wise matrix inversion (see Equation ( 20)), with respect to FD-DAS. Note that for TD-RCap [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], the full set of the received signals y(t) has to be delayed for each pixel, and then used to construct the spatial covariance matrix, which is highly time consuming. In the FD, the CSM is obtained directly from raw signals and efficiently steered with fast matrix products. In this way, although the optimization operation with the inversion of the matrix remains with a huge amount of calculus, FD-RCap is a relative fast version for a given frequency of TD-RCap (similar comparison as in [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF]). As a point of comparison, and to give an order of magnitude to the time that can be saved, for the maps that are composed of 35,140 pixels, the low-complexity FD beamformers (i.e., FD-DAS with or without snapshot averaging, or FD-Pisa) run in the order of a few seconds, whereas the higher complexity FD beamformer (FD-RCap) requires one or a few minutes, and the time-domain beamformers (TD-DAS, TD-RCap) easily exceed an hour, mainly due to the very time consuming delaying step over 200 s of the integration time. Alternatively, the FD-Pisa, as the Pisarenko class beamforming, offers intermediate performance between FD-DAS and FD-RCap for a computing load equivalent to FD-DAS. Finally, from a statistical point of view, the CSM formalism developed for power-estimate-based beamformers addresses the issue of quantitative power estimation of cavitation, which is an active and challenging field of research in terms of both numerical [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF] and instrumentation [START_REF] Gray | Diffraction effects and compensation in passive acoustic mapping[END_REF] aspects.

3) Focus on FD-Music, the nonpower-estimate-based approach: Apart from the power-estimate-based beamformers, FD-Music is well known as a high-resolution localization method [START_REF] Krim | Two decades of array signal processing research: the parametric approach[END_REF]. It has been demonstrated that FD-Music highlights the presence of sources, but does not provide a pixelwise power estimate: absolute quantification or relative quantification of the power of a source is not possible. When monitoring cavitation, weak power sources have to be detected very early in their development if they are located in a nontargeted area, to avoid undesirable effects like healthy tissue ablation. The simulation results in Fig. 6 and the experimental results in Fig. 8 show that FD-Music efficiently highlights weak sources. A specific aspect of this beamformer is that the background is smoother than that of the power-estimatebased beamformers (i.e., FD-DAS, FD-RCap, FD-Pisa), and therefore weak power sources are clearly identifiable. FD-Music is thus an interesting tool for localization of weak power sources, a situation where power-estimate-based beamformers are not appropriate. Furthermore, FD-Music has a running time similar to FD-DAS. Finally, unlike for FD-RCap or p for FD-Pisa, which are chosen empirically [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], [START_REF] Dougherty | Functional Beamforming[END_REF], the adaptive parameter K s of FD-Music can be adaptively chosen in line with the number of sources. For this experimental study, a very simple trial-and-error approach is used, as described in simulation, as a trade-off between the correct number of expected sources and the image quality metrics. However, the deduced K s values are slightly above our expectations (i.e., K s = 3 for a single spot; K s = 4 for two spots). This would suggest that a single bubble spot should instead be looked upon as several sources. Note also that for the sake of simplicity, we consider a constant value for K s that is relevant through the whole HIFU shot duration. Automated and more sophisticated EV thresholding can be investigated [START_REF] Wong | On information theoretic criteria for determining the number of signals in high resolution array processing[END_REF], [START_REF] Baranger | Adaptive spatiotemporal SVD clutter filtering for ultrafast doppler imaging using similarity of spatial singular vectors[END_REF].

4) Methodological perspectives :

The proposed formalism is based on second-order statistics of CSM and its subsequent assumptions of ergodic and stationary signals. This assumption is often assessed with respect to acquisition duration. This preliminary study illustrates that in some simple cavitation scenarios, these assumptions are verified, although further investigations are needed. For example, a higher-order statistics beamformer has been proposed for passive cavitation imaging in the time domain [START_REF] Lyka | Passive acoustic mapping using data-adaptive beamforming based on higher order statistics[END_REF] that overcomes the TD-RCap resolution. Also, the assumption for uncorrelated sources should be investigated for extended bubble clouds.

These results demonstrate that single-frequency imaging provides accurate localization of cavitation, even if its spectral signature is broad. In the passive cavitation imaging literature [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF], [START_REF] Patel | Closed loop spatial and temporal control of cavitation activity with passive acoustic mapping[END_REF], due to the frequency selectivity of FD approaches, single-frequency imaging is often considered to enable discrimination of stable cavitation from inertial cavitation. However, composite multifrequency images can also be considered by superimposing maps at different frequencies (i.e., summing) to further reject incoherent noise on maps from one frequency bin to another [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF], [START_REF] Patel | Closed loop spatial and temporal control of cavitation activity with passive acoustic mapping[END_REF]. Note that the computing load increases almost proportionally with the number of reconstructed bins [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF]. Moreover, as the beamwidth is frequency dependent [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], constant-beamwidth approaches can be considered for better consistency of acoustic maps, even for adaptive beamformers [START_REF] Wang | Constantbeamwidth and constant-powerwidth wideband robust Capon beamformers for acoustic imaging[END_REF]. A specific aspect to adaptive beamformers would be that the adaptive parameters can be frequency dependent, and as the error uncertainties might not be constant over the full bandwidth of the probe, this would potentially influence and p [START_REF] Dougherty | Functional Beamforming[END_REF], [START_REF] Wang | Constantbeamwidth and constant-powerwidth wideband robust Capon beamformers for acoustic imaging[END_REF].

This paper focuses on the enhancement of detection and localization of possibly simultaneous cavitation events, but to go further, the integration of maps over frequencies from power-estimate-based beamformers might provide access to power quantification within a given band, as in [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF], [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Baron | Numerical and experimental assessment of functional beamforming for source quantification[END_REF]. Nevertheless, this requires further challenging considerations, which include fine numerical scaling [START_REF] Haworth | Quantitative frequency-domain passive cavitation imaging[END_REF], [START_REF] Harris | On the use of windows for harmonic analysis with the discrete Fourier transform[END_REF], and also the appropriate calibration and acoustic modeling [START_REF] Coviello | Passive acoustic mapping utilizing optimal beamforming in ultrasound therapy monitoring[END_REF], which are beyond the scope of this paper.

Cavitation is a three-dimensional physical phenomenon. Imaging this with linear arrays limits its localization into a two-dimensional plane. The subsequent issues of this aspect are that off-plane cavitation events might not only be missed, but might also generate artifacts on the two-dimensional maps. On the one hand, the potential for adaptive beamformers to remove such artifacts has to be investigated. On the other hand, the CSM-derived beamformers are good candidates to deal with three-dimensional imaging. As these are intrinsically formalised in the FD, they are very computationally efficient [START_REF] Arvanitis | Passive acoustic mapping with the angular spectrum method[END_REF]. This property is essential to efficiently process the large volumes of data involved in three-dimensional imaging, for which more than 1000 sensors can be used [START_REF] Petrusca | Fast volumetric ultrasound b-mode and doppler imaging with a new high-channels density platform for advanced 4d cardiac imaging/therapy[END_REF].

VII. CONCLUSIONS

This study has shown the potential of CSM-derived beamformers for spatio-temporal monitoring of cavitation. These provide: (1) high quality maps in (2) shorter computational times than needed for time-domain approaches, and with (3) direct frequency selectivity. These three attractive properties are gathered within the proposed CSM-derived formalism, which is therefore expected to be a powerful tool to develop passive ultrasound imaging. In the long term, these CSMderived beamformers can be expected to contribute to the enhancement of spatio-temporal regulation of cavitation for therapeutic applications [START_REF] Patel | Closed loop spatial and temporal control of cavitation activity with passive acoustic mapping[END_REF].

Fig. 1 :

 1 Fig. 1: Experimental set-up showing the xz-plane view (a) and the yz-plane view (b). Wire #1 is placed in the xz-plane, to cross the acoustic axis at the theoretical focus of the HIFU transducer. HIFUs propagate in the opposite direction to the x-axis. When the HIFU transducer fires, a cavitation cloud is created: the acoustic signal emitted by the bubbles is passively received by the imaging probe. In the second configuration, a second wire #2 is added 3 mm downstream of wire 1, and parallel to it. In (a), white dots indicate the portion of the xz-plane that is imaged by the high-speed camera.

Fig. 2 :

 2 Fig. 2: Descriptions of the four simulated configurations and their respective goals. Note that the wavelength is λ = 0.31 mm.

Fig. 3 :

 3 Fig. 3: Performances of the four beamformers, considering the four image quality metrics of (a) R ax , (b) R lat , (c) CAR, and (d) Err pos versus the SNR, for the single-source simulation (configuration 1). Two adaptive parameters are used for each adaptive beamformer. As FD-Music is not a power-based beamformer, it is set apart for the CAR plots.

Fig. 4 :

 4 Fig. 4: Lateral slices at z = 60 mm (a) and axial slices at x = 5 mm (b) for two closely spaced sources, as in configurations 2 and 3, respectively. Four beamformers are represented: FD-DAS, FD-RCap with = 1, FD-Pisa with p = 50, and FD-Music with K s = 2. The actual source positions are plotted as the dashed lines. The plots are normalized by their maxima.As FD-Music is not a power-based beamformer, it a specific y-axis.

Fig. 5 :

 5 Fig. 5: Influence of the cross-spectral density matrix estimation on the cavitation maps obtained with the four different beamformers, in simulation. The maps are reconstructed at 5 MHz. Each column corresponds to a beamformer, as indicated. For FD-RCap, = 5, for FD-Pisa, p = 50, and for FD-Music, K s = 3. Different parameters for CSM estimation are tested for each row: (a-d) K = 1 (i.e., without snapshot averaging); (e-h) K = 130 without overlap, which leads to T snap = 1.5 µs; (i-l) K = 130 with p ov = 50 %, which leads to T snap = 3.1 µs; and (m-p) K = 130 p ov = 90 %, which leads to T snap = 15.4 µs. The actual bubble positions are the centers of the white circles. All of the maps are normalized and displayed with the same dynamic range, of 20 dB.

Fig. 6 :

 6 Fig. 6: Beamformer comparisons in simulation. The maps are reconstructed at 5 MHz. All of the maps are reconstructed with K = 130 and p ov = 90 % for CSM averaging. The first row corresponds to the power-estimate-based beamformers with optimal settings for FD-DAS (a), FD-RCap with = 5 (b), and FD-Pisa with p = 50 (c). The second row corresponds to the maps obtained for FD-Music, with K s = 1 (d), K s = 3 (e), K s = 5 (f), and K s = 10 (g). (h) The normalized eigenvalues are plotted logarithmically. The actual source positions are the centers of the white circles.(a, b, c, e) correspond to the maps of Fig.5m, n, o, p, respectively, with different dynamic ranges for the displays. The maps are normalized and displayed with their own full dynamic ranges: from the lowest intensity in their own background, up to 0 dB.

Fig. 7 :

 7 Fig. 7: Representative experimental localization maps for a single cavitation spot on a wire immersed in degassed and filtered water, for FD-DAS (a), FD-RCap (b) ( = 10), FD-Pisa (c) (p = 30), and FD-Music (d) (K s = 3), in the middle of the HIFU shot (10 ms). The HIFU transducer is on the right in the images. The maps are reconstructed at 4.7 MHz. The bubbles detected with the optical camera are located within the white/ red dashed lines. The dotted rectangle indicates the artifact region for the CAR measurement, with respect to the maximum amplitude of the map. The performance statistics are calculated on the 100 maps reconstructed during the 20 ms of shooting, and are shown for axial resolution (e), lateral resolution (f), and CAR (g), with the median in red, the mean in green, the first and third quartiles in blue, and the total extent in black.

Fig. 8 :

 8 Fig. 8: Representative experimental localization maps for the two cavitation spots on the two wires immersed in the filtered water, for FD-DAS (a), FD-RCap (b) ( = 10), FD-Pisa (c) (p = 30), and FD-Music (d) (K s = 4), at t=4 ms. The HIFU transducer is on the right of the images. The maps are reconstructed at 4.7 MHz. The bubbles detected with the optical camera are located within the white/red dashed lines. Projected maximum map values on the x-axis versus time is shown for each map, respectively (e-h). The projection considers the depth range 54 mm to 68 mm. The vertical dashed lines are the positions of the wire centers throughout the shot.

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TUFFC.2020.3032345, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control

0885-3010 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

TABLE I :

 I Summary of the beamformer characteristics. of length T rec is divided into partially overlapping snapshots y k (t). The overlap p ov is expressed as a percentage of T snap . The parameters to estimate the CSM are then linked by the following relationship:

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TUFFC.2020.3032345, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control

Authorized licensed use limited to: University of Grenoble Alpes. Downloaded on December 10,2020 at 15:00:38 UTC from IEEE Xplore. Restrictions apply.

0885-3010 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

Authorized licensed use limited to: University of Grenoble Alpes. Downloaded on December 10,2020 at 15:00:38 UTC from IEEE Xplore. Restrictions apply.

This study was supported by the LABEX CELYA (ANR-10-LABX-0060) and was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


E Mathurin

C Dartois

J Engrand

A Duprat

A Deniset-Besseau

Y Dazzi

T Kebukawa

David Noguchi

Troadec

Samples J Mathurin

E Dartois

C Engrand

J Duprat

A Deniset- Besseau

A Dazzi

Y Kebukawa

T Noguchi

D Troadec

NANOMETRE-SCALE INFRARED CHEMICAL IMAGING (AFM-IR) OF ORGANIC MATTER IN ULTRA-CARBONACEOUS ANTARCTIC MICROMETEORITES (UCAMMS) AND FUTURE ANALYSES OF HAYABUSA

Nanometre-scale infrared chemical imaging (AFM-IR) of organic matter in ultra-carbonaceous Antarctic micrometeorites (UCAMMS) and future analyses of

Hayabusa 2 samples.

Introduction:

The chemical composition of organic matter (OM) in interplanetary samples (meteorites and micrometeorites) is suitably characterized by the distribution of the different chemical bonds using infrared (IR) vibrational spectroscopy (see e.g. [1]). Classical IR microscopy provides a global view of the dust grain chemical structure content but remains limited by the diffraction, with typical spot sizes sampling a few micrometers in the mid-IR range. This spatial resolution limitation is well above that of complementary techniques such as isotopic imaging with NanoSIMS or transmission electron or X-ray microscopy techniques. These techniques reveal mineralogical, chemical and isotopic heterogeneities at the sub-micron scale but do not give full access to the distribution of the various chemical bonds. The IR diffraction limitation can be circumvented by using AFM-IR microscopy. This technique opens a new window for studies of OM at ten to tens of nanometer scales and will be of importance for studies of the samples from carbonaceous asteroid Ryugu, returned by the Hayabusa 2 space probe in December 2020.

AFM-IR is now a well-established microscopy technique in the vibrational field. It combines an atomic force microscope (AFM) and a tunable IR source to detect photo-thermal effect and access chemical information down to a nanoscale resolution [2]. This technique is now applied in a wide diversity of scientific fields [3], and was recently used to analyze extraterrestrial OM [4,5]. We report here on recent results obtained on imaging two UltraCarbonaceous Antarctic MicroMeteorites (UCAMMs) using AFM-IR [5]. A small fraction of the Antarctic micrometeorites from the Concordia collection consists in UCAMMs, particles with extreme concentrations in OM, most of them exhibiting large deuterium excesses [6]. UCAMMs are also found in Japanese interplanetary dust collections [7][8][9]. These UCAMMs most likely originate from the surface of small icy bodies in the outer regions of the solar system [1,6,7,10]. The large OM fraction of UCAMMs (considerably higher than in the most carbon-rich meteorites) enables direct analyses without the pre-treatment generally applied to extract the OM from other meteoritic samples, and give access to unaltered chemical maps of the intimate association of minerals and organics.

Samples and methods: Two UCAMM fragments from the Concordia collection were analyzed by AFM-IR: DC06-07-18 (hereafter DC06-18) and DC16-14-309 (hereafter DC16-309). DC06-18 was prepared by crushing a small fragment on a diamond window. DC06-309 was ultramicrotomed after embedding in sulfur; sulfur was then sublimated from the sections overnight in a desiccator at about 50°C. After AFM-IR analysis, complementary high resolution scanning electron microscopy images and EDX spectral maps of elements from C to Fe where obtained on the samples. In contact mode, the tip is in close contact with the surface and explores mainly the repulsive potential. In tapping mode, the cantilever oscillates and approaches the surface only for a short time, oscillating back and forth on the potential curve.

We used two different AFM-IR setups: i) contact mode setup where the AFM tip stays in contact with the sample (Fig. 1 left panel), and ii) tapping mode setup where the AFM tip is oscillating at the surface of the sample (Fig. 1 right panel). The latter mode is more adapted to the AFM-IR imaging on loosely bound samples. AFM-IR maps were acquired in contact mode for DC06-18 and tapping mode for DC16-309. In both setups, two acquisition modes are possible. In the first mode, the emission wavenumber is fixed and the tip moves along the surface allowing to acquire simultaneously AFM topography and IR absorption map of the sample at a selected wavenumber. In the second mode, the tip position's is fixed and the IR laser is tuned to explore the whole spectral range (here 1900 cm -1 -1350 cm -1 ) giving access to local absorption spectra at the point of interest.

Results: Figure 2a shows the local absorption at different wavenumbers on a microtomed section of UCAMM DC06-309, each wavenumber corresponding

1983.pdf 51st Lunar and Planetary Science Conference (2020)

to the specific absorption of a given chemical bond. The AFM-IR maps reveal the heterogeneity of the absorption of the C=O (1710 cm -1 , in green) and C=C (1600 cm -1 , in red) bonds. Based on these maps, ratio of different absorptions can be obtained to map the functional group diversity. Figure 2b shows the ratio of the C=O bond (green) versus the C=C bond (red). This ratio shows heterogeneities appearing as patches of tens to hundreds of nanometers in size. The size of these heterogeneities is consistent with complementary results obtained by STXM-XANES for UCAMMs [9,11,[START_REF] Guérin | pdf 51st Lunar and Planetary Science Conference[END_REF]. Local spectra were also acquired (Figure 2c) to allow a spectral identification of vibrational signatures at very fine scale. The heterogeneity of such signatures escape detection in a classical diffraction-limited IR analysis with typical sampling of several microns in size. AFM-IR also allowed the identification of a submicrometric contaminant brought during the sample preparation of DC16-309, as observed at 1460 cm -1 (light blue in Fig. 2a and spectra lines S2 & S3 in Fig. 2c) (see [5]).

Conclusions and outlook for the analysis of Hayabusa 2 samples:

The chemical functional group diversity of OM in UCAMMs without pre-treatment can be quantitatively characterized at a sub-micron scale by AFM-IR [5]. Using this technique with an extended spectral range should increase the outcome of the chemical characterization of the OM and will better define the association (and potential interaction) between the OM and the minerals in the sample.

In addition to the preliminary work by Kebukawa et al. [4], this first study paves the way for the development of AFM-IR studies of different extraterrestrial objects such as carbonaceous chondrites. More recently, we successfully applied AFM-IR in tapping mode to carbonaceous chondrites for the first time. To prepare the analysis of Hayabusa 2 samples, a preliminary study of the Jbilet Winselwan (heated CM2) carbonaceous chondrite was performed, and heterogeneous distribution of C=O and C=C was observed (Figure 3). One of the main objectives of these investigations is to develop an upfront and dedicated expertise on complex materials such as extraterrestrial samples, to fully apply the AFMIR technique on the samples soon returned by the Hayabusa 2 mission.

Figure 1 :

 1 Figure 1: Schematic view of the AFM-IR technique in contact (left panel) and tapping (right panel) modes. The interaction potential normal to the surface is illustrated by the colored line. The approximate excursion ranges on this interaction potential for the two AFM modes are shown in yellow.In contact mode, the tip is in close contact with the surface and explores mainly the repulsive potential. In tapping mode, the cantilever oscillates and approaches the surface only for a short time, oscillating back and forth on the potential curve.

Figure 2 :

 2 Figure2: AFM-IR study of UCAMM DC16-309. a) AFMIR maps recorded at 1710 cm -1 (green), 1600 cm -1 (red), and 1460 cm -1 (blue), composite (red, green, blue) image of the three maps, and AFM height and back-scattered electron (BSE) images, b) ratio of the AFM-IR maps at 1710 and 1600 cm -1 . The white isocontour delineates the height corresponding to 5% of the maximum height measured by the AFM tip, c) AFM-IR spectral maps labelled S1 to S4, performed in selected locations shown on the b) map in the 1900-1350 cm -1 range. Vertical green and red lines on c) mark the 1710 cm -1 and 1600 cm -1 wavenumber positions.

Figure 3 :

 3 Figure 3: AFM-IR study of Jbilet Winselwan (heated CM2) carbonaceous chondrite -AFMIR maps recorded at 1710 cm - 1 (green), 1600 cm -1 (red) and AFM height

Acknowledgments: This work was supported by the French Agence Nationale de la Recherche (COMETOR, ANR-18-CE31-0011) as well as INSU (PNP, PCMI), IN2P3, CNES, DIM-ACAV (Région Île de France), CNRS, the "Défi Origines" from the MITI (CNRS) and Université Paris-Sud. We are grateful to the French and Italian polar institutes IPEV and PNRA, for their financial and logistic support to the micrometeorites collection at the vicinity of the CONCORDIA station (Dome C). A part of this study was supported by JSPS Core-to-Core program "International Network of Planetary Sciences."