On the robustness of multilateration of ionoacoustic signals for localization of the Bragg peak at pre-clinical proton beam energies in water

Objectives. The energy deposited in a medium by a pulsed proton beam results in the emission of thermoacoustic waves, also called ionoacoustics (IA). The proton beam stopping position (Bragg peak) can be retrieved from a time-of-flight analysis (ToF) of IA signals acquired at different sensor locations (multilateration). This work aimed to assess the robustness of multilateration methods in proton beams at pre-clinical energies for the development of a small animal irradiator. Approach. The accuracy of multilateration performed using different algorithms; namely, time of arrival and time difference of arrival, was investigated in-silico for ideal point sources in the presence of realistic uncertainties on the ToF estimation and ionoacoustic signals generated by a 20 MeV pulsed proton beam stopped in a homogeneous water phantom. The localisation accuracy was further investigated experimentally based on two different measurements with pulsed monoenergetic proton beams at energies of 20 and 22 MeV. Main results. It was found that the localisation accuracy mainly depends on the position of the acoustic detectors relative to the proton beam due to spatial variation of the error on the ToF estimation. By optimally positioning the sensors to reduce the ToF error, the Bragg peak could be located in-silico with an accuracy better than 90 μm (2% error). Localisation errors going up to 1 mm were observed experimentally due to inaccurate knowledge of the sensor positions and noisy ionoacoustic signals. Significance. This study gives a first overview of the implementation of different multilateration methods for ionoacoustics-based Bragg peak localisation in two- and three-dimensions at pre-clinical energies. Different sources of uncertainty were investigated, and their impact on the localisation accuracy was quantified in-silico and experimentally.


Introduction
Owing to their characteristic inverse depth dose profile, proton beams deliver their maximal dose to tissue inside a confined volume (so-called Bragg peak, BP), offering better sparing of healthy tissues compared to other external beam therapies such as, electron or photon beam therapies (Newhauser andZhang 2015, Baumann et al 2016). However, there are still physical and biological challenges in ion beam therapy limiting its full potential exploitation. One challenge is represented by range uncertainties which are, among others, present because of anatomical changes between the treatment fractions (Li et al 2015), inaccurate knowledge of the tissue stopping power relative to water, and patient misalignment (Andreo 2009, Paganetti 2012, Grün et al 2013, Knopf and Lomax 2013. Therefore, the possibility of in vivo range verification during proton therapy is expected to help reducing the safety margins applied during the treatment (irradiated volume typically increased by 3% of the range along the beam axis) and hence decreasing the dose to the healthy tissues.
When energy is deposited in a medium by a pulsed proton or ion beam, it generates thermoacoustic waves, also called ionoacoustics (IA) (Jones et al 2014, Assmann et al 2015. Indeed, assuming sufficiently short proton pulses (thermal and stress confinement 4 ), the irradiated volume briefly heats up due to the energy loss along the proton beam path, mostly caused by Coulomb inelastic interactions with the electrons in matter. The resulting initial pressure distribution is proportional to the deposited energy multiplied by a material-dependent conversion factor (Grüneisen parameter) (Hickling et al 2018). In homogeneous media, the acoustic emissions predominantly emerge from the BP volume allowing to retrieve the maximum of the deposited energy (or dose) from ToF estimation using acoustic localisation methods such as triangulation or multilateration. Alternatively, setting up a larger number of sensing elements, the full two-or three-dimensional dose can in principle be accurately reconstructed (Kellnberger et al 2016, Lascaud et al 2021b, Yao et al 2021 as the wave characteristics provide three-dimensional information on the underlying energy distribution. In recent years, acoustic BP localisation based on linear least square approach (LS) using time of arrival (TOA) algorithm has been proposed as a method for range verification , Patch et al 2018. These preliminary simulation studies have shown the first evidence of acoustic source localisation feasibility for in vivo range verification during prostate and liver irradiation. Acoustic BP localisation in a homogeneous medium based on nonlinear least square optimisation and the time difference of arrival (TDOA) algorithm has been proposed recently (Otero et al 2019). Contrary to TOA, which takes as inputs the ToFs estimated at each sensor position, TDOA relies on estimating the ToF difference between a reference sensor and all remaining sensors included in the network. Based on simulated and experimental data, the authors (Otero et al 2019) highlighted that the localisation accuracy depends on the number of sensors and their arrangement relative to the source position. In a complementary in-silico study from the same authors (Otero et al 2020), it was shown that BP could be located in the brain with an accuracy of 1mm. For all studies conducted by Otero and co-workers, the ToF was evaluated using generalized cross-correlation (Knapp and Carter 1976), but the effect of the ToF estimation method on the localisation accuracy was not discussed.
Since the ionoacoustic signal encodes information on the dose distribution, its temporal evolution (the signal shape) depends on the portion of the proton beam seen by a given acoustic sensor. Two main wavefronts are typically observed in homogeneous media . A quasi-spherical wave propagating from the BP region (γ-wave) and a cylindrical wavefront (α-wave) from the plateau region. The energy variation at the position where the proton beam enters the phantom leads to the emission of an additional plane wave (entrance signal). The superposition of these different wavefronts results in a variation of the signal shape depending on the sensor position, which hampers accurate ToF estimation . Other sources of uncertainties on the ToF determination in homogeneous media result from the detector response, the noise and an inaccurate knowledge on the measurement starting time (Lehrack et al 2017). In heterogeneous media, the variation of the speed of sound in the different tissues may cause additional errors in the ToF estimation if not properly taken into account .
The present study investigates the robustness of BP localisation using TOA and TDOA multilateration algorithms at pre-clinical proton beam energies (20 and 22 MeV), in the context of the SIRMIO project . A compact portable system is developed along the SIRMIO project, allowing for precision imageguided small animal proton irradiation at existing clinical proton centres. The platform will include a dedicated beamline for degradation and refocusing of the clinical beam to relevant pre-clinical energies with submillimeter beam dimensions (Gerlach et al 2020). The anatomical images will rely on proton computed tomography of the animal in the treatment position before irradiation (Schnürle et al 2023), combined with ultrasound imaging for enhanced tumour visualisation. IA will be used for proton range verification with pulsed beams delivered by clinical synchro-cyclotron accelerators, to complement an in-beam positron emission tomography system foreseen for range monitoring at isochronous cyclotrons and synchrotrons, as well as for biological guidance. The compactness of the SIRMIO platform and co-integration of the different modalities limit the space available around the animal. Consequently, the number of IA sensors and their position should be optimised to not interfere with the treatment delivery and anatomical imaging, while allowing for accurate BP localisation. To this aim, the robustness of different localisation methods and sensor arrangements is assessed. The choice of the reference sensor for TDOA and the impact of the ToF uncertainties on the source localisation precision and accuracy are first estimated in ideal two-and three-dimensional scenarios. To this aim, the ToFs are geometrically derived from the Euclidean distances, and ToF errors are introduced. Thereafter, the accuracy of the ToF extracted from IA signals generated by a pulsed 20 MeV proton beam in water is assessed in-silico for 4 The heating time is short enough to neglect the thermal diffusion in the irradiated volume and shorter than the time required for an acoustic wave to propagate across the heated volume. several ToF estimation methods and different sensor positions. Lastly, BP localisation is investigated experimentally for two different sensor arrangements and proton beam energies.

Material and methods
2.1. Simulation study 2.1.1. Robustness analysis in ideal scenarios The geometries used to assess the influence of the ToF error on the multilateration outputs (robustness analysis) are illustrated in figure 1. Multilateration algorithms were first investigated in two-dimensions. An ideal source was virtually moved over 126 positions onto an orthogonal grid of 56 × 36 mm 2 by steps of 4 mm (see figure 1(a)). The ToFs were determined from the Euclidean distance between the source and the sensors, assuming a speed of sound of 1500 m s −1 . Here only 3 sensors (namely Bs 1 , Bs 2 and Bs 3 ) were used, as illustrated in figure 1(a). The sensor positions were chosen based on the setup available for the ionoacoustic experiments described in section 2.2.1. The simulation geometry was later extended to three-dimensions, as shown in figure 1(b) (grid of 56 × 36 × 48 mm 3 with a grid spacing of 4 mm), and a fourth sensor (Bs 4 ) was added on the top. To facilitate the result analysis, the field-of-view (FOV) of the sensor network is introduced and defined as the surface (or volume in 3D) enclosed by the sensors, as annotated in figure 1.
Random uncertainties in the ToF estimation were considered aiming to model speed of sound variations, inaccurate knowledge of the sensor's spatial location and error on the ToF. For that reason, the ToF error for each sensor was sampled from a multivariate normal distribution ( , ( ) m S  ) with mean zero and a standard deviation (σ random ) equal to 5% of the absolute ToF to reproduce in vivo variation of the speed of sound. There was no correlation between random uncertainties of individual sensors hence, off-diagonal elements in , ( ) m S  were zero. Systematic uncertainties (i.e. coming from an inaccurate knowledge of the measurement starting time) were also considered and modelled by a multivariate Normal distribution. The standard deviation (σ syst ) was set to 1 μs according to previous studies from our group (Lehrack et al 2017, Wieser et al 2021. The systematic uncertainties on the ToF were assumed to be perfectly correlated between all the sensors (off-diagonal elements in , ( ) m S  were >0). In the results presented hereafter, random and systematic uncertainties were modeled at the same time. The number of samples for each type of uncertainty was fixed to 150 (N rand = 150 and N syst = 150), which gives a total of 22 500 samples for each source position, chosen to have a representative sampling of the multilateration error distribution.

Ionoacoustic simulations
The performances of the multilateration algorithms (TOA and TDOA) for realistic ionoacoustic signals were evaluated in-silico in a homogeneous water phantom. The study's objective was to quantify the ToF error dependency on the sensor position. The three-dimensional dose distribution (D(r)) of a 20 MeV monoenergetic proton beam was obtained from the FLUKA Monte Carlo code (version 2020.0.4, using PRECISIOn defaults) (Ferrari et al 2005, Böhlen et al 2014. Additional information on the simulation setup can be found in the supporting material of Lascaud et al (2021a). The initial pressure (p 0 (r)) was deduced from the deposited energy (D(r) × ρ, with ρ the water density equals to 998 kg m −3 ) multiplied by the Grüneisen parameter (Γ=0.11). The initial pressure was propagated in three dimensions onto an anisotropic grid (spacing of 75 μm along the proton beam axis and lateral spacing of 150 μm) using the MATLAB k-Wave toolbox (MATLAB R2019a, MathWorks, Natick, MA) (Treeby and Cox 2010). Note that higher frequencies are propagated along the proton beam axis due to the sharpness of the BP (300 μm FWHM) compared to the lateral beam dimensions (2.5 mm FWHM at the phantom entrance was experimentally estimated from a Gafchromic film measurement), which motivated the choice of the anisotropic grid spacing. A 1 mm thick air gap positioned downstream to the proton beam was included in the k-Wave simulation setup, as illustrated in figure 2(b), to properly model the acoustic emission at the interface between the air and the water. The pressure waves were recorded using ideal point sensors (843 individual sensors) arranged in a semi-circular configuration with a diameter of 60 mm. The sensor network was positioned such that the center of the arc corresponds to the maximum of the proton dose. The arc dimension was chosen to be comparable with the setup used in the experiment shown in figure 2(a). After propagation, each signal was convolved with a 200 ns square pulse to account for the proton pulse shape observed experimentally.

Ionoacoustic experiments 2.2.1. 3-sensor configuration
The experiment was conducted at the Maier-Leibnitz-Laboratory in Munich (Germany) with a 20 MeV monoenergetic pulsed proton beam in a water-filled aluminium box. For this first experiment, the beam was delivered by pulses of 200 ns at a repetition rate of 4.9 kHz. The beam current was set to 3 nA, corresponding to a dose per pulse of 1.69 Gy deposited at the BP. As illustrated in figure 2(a), the proton beam entered the setup by an air-channel terminated by a 50 μm thick polyimide foil (entrance window). The ionoacoustic wavefront was recorded using 3 detectors positioned on the proton beam axis (axial sensor) and on an axis perpendicularly to it (lateral sensors). The sensors were directly mounted into the aluminium box designed such that the BP was at the focus of all the detectors through water-sealed apertures in the wall, ensuring a fixed position of the detectors relative to the entrance window during the experiments. All the single element detectors were focused piezoelectric transducers manufactured by Olympus-Parametric (3.5 MHz centre frequency with a 73% fractional bandwidth and a diameter of 12.7 mm). The axial sensor Bs 1 had a focal length ( fl) of 50.8 mm, whereas the focal length was equal to 25.4 mm for the two lateral detectors (Bs 2 , Bs 3 ).
The setup was mounted onto a motorised three-axis stage used for alignment. To this purpose, ionoacoustic measurements were performed with the setup moved at different locations along the x-and y-axis (as defined on figure 2(a)) by steps of 1 mm. The setup was assumed to be aligned with the proton beam (hereafter referred to as the on-axis position) at the location for which the highest signal amplitude was detected on Bs 1 . Additional measurements were performed after shifting by 5 mm the phantom along the x-axis (off-axis position).
The ionoacoustic signals were amplified by 60 dB using a low-noise amplifier (HVA-10M-60-B, FEMTO Messtechnik GmbH, Germany) before being acquired with a digital oscilloscope (6404D PicoScope, Pico Technology Ltd., GB) at a sampling frequency of 156.25 MHz. The signal acquisition was triggered using a synchronisation signal delivered by the chopping system of the particle accelerator. Such a trigger results in a time offset of 1.43 μs between the starting time of the measurement and the time at which the first protons hit the target, as previously characterised (Wieser et al 2021), which was later subtracted to the estimated ToF. The measurements were systematically composed of 1000 consecutive acquisitions (1000 proton pulses). The temperature was monitored during the experiment using a PT1000 probe immersed in the water. The mean temperature was found to be 21.92 • C, which corresponds to a speed of sound in water of 1488.1 ms −1 , later used to estimate the BP position from the signal ToF. For each acquisition, the signal was denoised in postprocessing using the Daubechies wavelet family with 5 decomposition levels (Ayat et al 2009, Dautov and Özerdem 2018, Baldazzi et al 2020, Vallicelli et al 2021 and the BP was localised with a single shot (1.69 Gy).

5-sensor setup
Additional experiments were performed in a water tank using a 5-sensor arrangement mechanically assembled together using a dedicated holder (horseshoes), as represented in figure 3. For this setup, the beam energy was increased to 22 MeV. The beam current was set to 3.5 nA at a repetition rate of 10 kHz, corresponding to a dose per pulse of 0.58 Gy delivered at the BP. Similarly to the previous experiments, 12.7 mm diameter piezoelectric single elements with a 3.5 MHz center frequency were used and arranged as follows: one axial transducer Bs 2 ( fl = 50.8 mm), two transducers inclined by an angle of 28°w.r.t. the beam axis, i.e. Bs 1 and Bs 3 ( fl = 25.4 mm), and two lateral transducers Bs 0 (unfocused) and Bs 4 ( fl = 25.4 mm), as shown in figure 3. Due to the limited number of channels available on the picoScope, the ionoacoustic signals were acquired sequentially (first with Bs 0 to Bs 2 , and second with Bs 2 to Bs 4 ). As earlier, the measurements were repeated over 1000 consecutive proton pulses. Due to the orientation of the transducers Bs 1 and Bs 3 with respect to the beam axis and lower dose per pulse, the signal-to-noise ratio (SNR) for the IA signals recorded by the two tilted transducers was lower compared to the previous experiments. Therefore, the signals were averaged over 50 acquisitions (total peak dose of 29 Gy deposited at the BP) before performing the wavelet filtering.
The spatial location of the transducers (c 0 to c 4 ) were estimated with a 3D optical tracking system (Hybrid Polaris Spectra System, Northern Digital Inc., Waterloo, ON, Canada) using a passive marker tool. For focused detectors, the sensor location was defined as the projection of the outer rim's centre onto the curved surface, as depicted in figure S9 of the supplementary material. To this aim, the position of the rim's center was estimated by locating 3 points on the circle, and its projection on the curved surface was retrieved from a fourth measurement on the curved surface. The transducer localisation's accuracy was improved by performing a 3D rigid registration (coherent point drift algorithm) between a template derived from the holder drawing (moving image) and the position obtained from the optical measurements. The position where the proton beam enters the water tank (location of the entrance window, c ew ) was determined by irradiating a Gafchromic film. The marked location of the irradiated area centre was later acquired with the optical tracking system. Figure 4 illustrates the multilateration workflow used for the BP localisation. The ToF was estimated for each sensor Bs 1 to Bs n , where n is the number of sensors, located at a spatial location c 1 to c j , respectively, and converted to a distance knowing the medium speed-of-sound (v medium ). Multilateration was performed afterwards either using TOA or TDOA algorithms. To this aim, the cost functions described in equations (1) and (2) for TOA and TDOA, respectively, were minimised using the Levenberg-Marquardt algorithm available from the Matlab software (Matlab R2021a) with 500 iterations and a termination tolerance set to 10 −9 .

Data analysis 2.3.1. Multilateration workflow
Two types of reference sensors were considered for the study based on ideal scenarios (point sources at multiple locations) using TDOA, hereafter referred to as static or dynamic reference sensors. For the former, the same reference sensor is used independently of the source location, whereas for the latter, the reference sensor was changed depending on the ToF reads-out on each sensor of the network. Three dynamic allocation strategies were employed to define the reference sensor set, namely, the minimum ToF (l Ref,min ), the maximum ToF ). For the latter, the reference sensor was selected to correspond with the detector for which the ToF is the closest (minimal difference) to the mean ToF across all the sensors.

Ionoacoustic time-of-flight extraction
For the studies with IA signals, the ToF was extracted from the direct signal (γ-wave), and different methods were compared. Figures 5(a) and (b) depict simulated IA signals recorded for different sensor positions (axial and lateral sensors, respectively) on which the several signal features used to extract the ToF are annotated. Those correspond to: • The zero-crossing between the compression and rarefaction pulses for the axial sensor, hereafter referred to as zero-crossing. As seen in figure 5(b), for lateral sensors, there are multiple zero-crossing. In that case, the first zero-crossing point was chosen.
• The maximum of the signal amplitude (max amplitude).
• The minimum of the signal amplitude (min amplitude).
• The maximum of the signal envelope is obtained from the absolute value of the Hilbert transform (max envelope).

Multilateration robustness assessment
Without uncertainties, the accuracy was defined by the error in position, which is expressed as: . Where x, y and z are the known source positions and x s , y s , z s are the reconstructed source position from the optimisation. The 3D metric is reducible to 1D (ò z ) or 2D When uncertainties were modelled, the first metric was changed to be the root mean square error RMSE, i.e. the RMSE between the actual and estimated source position. From the RMSE other metrics were derived as, for instance, the mean root mean square error, i.e. RMSE m , and the standard deviation of the root mean square error, The accuracy of various ToF extraction methods was assessed by evaluating the error in ToF estimation, i.e. ò ToF , which was defined as the product of the speed of sound in the medium, i.e. v medium , and the absolute difference between the known ToF (ground truth) and the ToF estimated for a particular extraction method. For all the results reported afterwards, the ground truth for the 20 and 22 MeV proton beam were determined from the corresponding FLUKA Monte Carlo simulations. Although the simulations were not compared to another independent measurement, the general good agreement obtained between the experimental and simulated IA signals gives confidence in the reliability of the model used, which was also previously validated in a time-offlight spectrometry study (Würl et al 2018). Moreover, for the experimental studies the SNR for each single sensor position, was defined as the ratio between the peak-to-peak amplitude and the noise standard deviation, defined as: SNR Bs

Comparison between TOA and TDOA in 2D and 3D ideal scenarios
The performance of 2D multilateration in ideal scenarios was first assessed without uncertainties. In this case, the error obtained for TOA is almost independent of the source localisation and its position relative to the sensor network FOV (error of ò xy = 0.81 ± 0.46 μm and 0.74 0.59 m xy FOV m =   for all the sources position and restricting to the FOV, respectively (see figure S1)). Conversely, the error rapidly increases for sources outside the FOV when using TDOA and depends on the reference sensor. For the geometry investigated, the lowest errors are obtained using a static sensor (Bs 1 ) as a reference (error of ò xy = 2.0 ± 1.8 μm and 0.78 0.05 xy FOV =   μm, (see figure S1)). The corresponding mean errors and standard deviations are summarized in table 1.
Similarly to the previous results, with uncertainties the source localisation is more accurate at the center of the array independently of the method employed. For TOA, the error progressively increases inside the FOV when going from the center (RMSE of 1.63 mm) toward the edges (RMSE of 4.03 mm), leading to a FOV RMSE m of 2.68 mm and a FOV RMSE s equals to 0.63 mm for the range of uncertainty modelled (see figure S4(a)). Using TDOA, the error is more homogeneous inside the FOV and slightly lower than for TOA at the center (RMSE = 1.61 mm). However, the error rapidly increases in the close vicinity of Bs 2 and Bs 3 , going up to 4.03 mm for the sources closest to Bs 2 . The mean root mean square error and its standard deviation were found to be 2.32 mm For the multilateration performed in 3D, only the results in presence of random and systematic uncertainties are reported. For TOA, the mean root mean square error inside the FOV is equal to 3.74 mm with a standard deviation of 0.54 mm (see figure S8(a)). The mean error reduces when using TDOA ( 3.35 mm FOV RMSE m = , whereas the standard deviation increases to 0.74 mm (see figure S8(b)). For both algorithms, as it was the case in 2D, the error increases when the source is located at the edges of the FOV with RMSE FOV going up to 5.35 mm and 8.37 mm for TOA and TDOA, respectively. The lowest errors are observed at the centre where RMSE FOV is as low as 2.58 mm and 2.47 mm for TOA and TDOA, respectively. All the results are reported in table S6. Figure 6(a) shows the error of the ToF determined from simulated ionoacoustic signals depending on the extraction method and for all the sensors (Bs 1 to Bs 843 ) in the arc array. For all the ToF extraction methods, the error in the ToF is minimum for axial sensors (Bs 422 ) and increases when going to lateral positions, i.e. Bs 1 , Bs 843 . The lowest errors are obtained by extracting the ToF from the signal max envelope, for which the error is lower than 0.1 μs at 4 positions near the beam axis (Bs 322 , Bs 403 , Bs 441 ,B 523 ). On average, for all the sensors positions investigated, the lowest error is obtained using zero-crossing method. Figure 6(b) shows the performances of TDOA depending on the reference sensor, i.e. from Bs 1 to Bs 843 , when the multilateration is carried out in 2D using all the 843 sensors of the array and extracting the ToF with the zero-crossing method. Similarly to the ToF error pattern, the localisation error is minimal for the reference sensor positioned on the beam axis (ò xz = 0.18 mm with Bs 422 as a reference) and increases when selecting a lateral sensor as a reference (up to ò xz = 1.60 mm when Bs 1 is the reference). It is worth noticing that the minimum error on the ToF estimation (for Bs 103 and Bs 743 ) does not correspond with the lowest localisation error.

Accuracy of ionoacoustic localisation depending on the sensor position
Shown in table 2 are the errors in BP localisation using TOA and TDOA (with Bs 422 as a reference) multilateration for various ToF extraction methods. Using 843 sensor positions, the ToF method resulting in the lowest error in the reconstructed BP position (error of 0.55 mm and 0.18 mm for TOA and TDOA, respectively) corresponds with those offering the lowest average ToF error (zero-crossing). Optimally selecting the sensors used for the multilateration and ToF extraction method to minimise the error on the ToF, the accuracy of the multilateration can be improved while decreasing the number of sensors required. For instance, using only 3 sensors (Bs 1 , Bs 422 and Bs 843 ) and evaluating the ToF from the maximum of amplitude of the signal envelope, 2.0 × 10 −3 6.9 × 10 −3 1.9 × 10 −3 3.8 × 10 −3 TOA none 8.1 × 10 −4 7.4 × 10 −4 4.6 × 10 −4 5.9 × 10 −4 which allows for the lowest error for axial sensors, the localisation error reduces to 0.09 mm and 0.35 mm for TOA and TDOA with Bs 422 as a reference, respectively. Figures 7(a) to (c) shows the IA signals recorded by the three detectors (axial sensor Bs 1 and the two lateral sensors Bs 2 and Bs 3 ) when the setup was positioned on the beam axis (on-axis) or after shifting it laterally by 5 mm (off-axis). As depicted on figure 7(a), the signal measured axially is made of three distinct pulses, namely the direct signal (from which the ToF was extracted), an entrance signal produced at the position where the proton beam enters the water phantom, and the reflection of the direct signal at the interface between air and water. As the energy gradients are sharper along the proton beam axis, the axial signal is of higher frequency and amplitude as the pressures measured laterally (figures 7(b) and (c) using Bs 2 and Bs 3 , respectively). Considering the 5 mm shift, with regards to the 12.7 mm diameter of the transducers used, moving the setup laterally mostly affects the ToF derived from the lateral sensors (increased ToF from Bs 2 and reduced for Bs 3 ).

Ionoacoustic experiments with 3 sensors
The distribution of the reconstructed lateral and axial locations of the BP for the 1000 consecutive measurements with the phantom positioned on the beam axis is presented in figures 7(d) and (e), respectively, in which both results obtained from TOA and TDOA are compared. The corresponding error on the BP position in  Table 2. Multilateration of the BP position in 2D for different ToF extraction methods. Two cases were considered, one for which all the 843 sensors were used, and a second, where only 3 sensors optimally selected were considered. The ground truth corresponds to the error on the multilateration when performed using a ToF derived from the Euclidean distance between the BP and the given sensors. For the multilateration performed with 843 sensors, the average error in ToF is also reported for all the sensor positions. ) arises from the reconstruction of the lateral position ( x ŝ =  0.37 ± 0.17 mm for TOA and TDOA, respectively), which is attributed to an inaccurate knowledge of the transducer position relative to the beam entrance (i.e. transducer location derived from the phantom center which was determined from a 1 mm step scan). With the phantom positioned off-axis, the total error in position increases to 1.26 ± 0.13 mm and 1.25 ± 0.14 mm for TOA and TDOA, respectively, resulting from inaccurate ToF estimation (i.e. the error on Bs 1 increases compared to the previous setup due to its almost lateral position) and lower SNR, as reported in table 4. Figures 8(a) and (b) show the results of BP multilateration in 3D performed with TOA and TDOA of the setup shown in figure 3. For this setup, for which particular attention was paid to the determination of the sensor locations, the total error on the determination of the 3D BP position reduces to 1.00 ± 0.72 mm and 0.82 ± 0.23 mm for TOA and TDOA, respectively. It is worth noticing that, without performing the proposed co-registration step to improve the accuracy of the sensor localisation, the BP multilateration error increases to 2.48 ± 0.30 mm for TOA and 8.82 ± 1.52 mm for TDOA.

TOA and TDOA algorithms in ideal scenarios
In an ideal scenario, i.e. without uncertainties, both TOA and TDOA perform similarly, although the accuracy of TDOA is more sensitive to the source location. Generally, the source should be inside the sensor array FOV, preferably in the centre. For the selected sensor geometry (36 mm height triangle) and range of uncertainties investigated, if the reference sensor is appropriately chosen, TDOA results in a lower localisation error (

Starting time of the signal acquisition
For the experiments reported in this study, the measurements were triggered using a stable synchronisation signal provided by the chopping system of the tandem accelerator. The time difference between the rising edge of the synchronisation signal and the actual moment when the protons enter the irradiated medium was previously determined with a precision better than 0.1 ns (corresponding to a distance smaller than 150 nm in water). The known offset on the acquisition starting time was thereafter systematically corrected during the data postprocessing. By doing so, the systematic error is negligible compared to the error on the ToF estimation, and therefore, similarly to the simulations, the best results are obtained with TOA. Depending on the accelerator and facility, the access to a stable synchronisation signal directly from the machine may not be possible, as pointed out in a previous study for which a jittering of up to 850 ns was observed (1.3 mm in 22 • C water) at a clinical facility (Lehrack et al 2017). Alternatively, the acquisition can be triggered using an external detector to sense secondary emissions such as prompt gammas resulting from the interactions between the proton beam and the irradiated medium. In that case, the correlation between the triggering events and the proton pulse time profile should be known as it may introduce an additional time offset on the estimation of the ToF. All in all, if the relation between the synchronisation signal and the proton pulse cannot be determined or if the signal is not stable enough, multilateration using TDOA algorithm should be favoured.

Multilateration from ionoacoustic signals in homogeneous media
For ionoacoustic-based BP localisation in homogeneous medium of known speed of sound and assuming accurate knowledge of the sensor locations, the most relevant sources of uncertainties are the starting point of the measurement (systematic error on the determination of the exact time at which the first protons enter the irradiated phantom) and the error on the ToF extraction. As illustrated by the in-silico study with an arc sensor arrangement, the error on the ToF depends on the detector location relative to the BP. It is hence not correlated between all the sensors composing the array. Consequently, TOA provides more accurate BP localisation, which can be optimised by adequately selecting the sensor position to minimise the ToF error. Moreover, as shown in the simulation studies, the performance of TDOA is intrinsically related to the choice of the reference sensor. In particular, for ionoacoustic applications, the error on the localisation can be up to 8 times higher depending on the reference sensor position and thus the error on the ToF estimation. However, the accuracy of the multilateration is not directly related to the ToF error observed on a given reference sensor (see figure 6(b)). The optimal reference sensor rather seems to correspond with positions for which the error on the ToF is close to the mean error, on average compensating for the ToF errors on the other sensors used for the multilateration.

Multilateration in heterogeneous media
The present study was limited to the irradiation of a homogeneous water phantom. In this case, the speed of sound in the medium is constant. However, the speed of sound varies in vivo depending on the tissue type, which is expected to reduce the BP localisation precision. Assuming a conservative 5 % error on the average speed of sound on the acoustic path, we showed that the error after multilateration increases by about 2 mm for the investigated geometry. Therefore, further implementation of multilateration for in vivo range verification should account for the speed of sound aberration, similar to the strategy employed to track wireless capsule endoscopy using radio-frequency signals (Mohan et al 2016, Khan et al 2018. To this end, the speed of sound could be empirically derived from the x-ray pre-treatment imaging (Mast 2000) or imaged with ultrasound (Rau et al 2021) shortly before the treatment delivery. It is worth mentioning that a previous in-silico study has shown that the over-and under-estimation of the speed of sound along the propagation path tends to cancel out for the investigated clinical applications (liver and prostate). As a result, similar localisation errors were obtained for speed of sound assumed to be either the one in water or the correct average speed of sound estimated from x-ray images .

Optimal sensor arrangement
Our results highlight the requirements for a preliminary study from the treatment plan (pre-treatment ionoacoustic simulation) to assess the ideal sensor position and reference sensor to mitigate the error on the ToF estimation. The adjustment of the sensor positioning should however be further studied, notably the robustness of the optimised arrangement when the BP is moved through the treatment field. In fact, as shown for the 3-sensor configuration experimentally investigated, similar localisation errors were obtained for TOA and TDOA when the proton beam was shifted by 5 mm laterally due to a larger error on the ToF. Therefore, the determination of the ideal sensor positions, their number and reference sensor for TDOA should be assessed for all the single beams composing the treatment plan (or at least for the ones with the highest intensity that are more reliable for the monitoring). Consequently, it also points out that a fixed sensor array geometry may not be optimal for multilateration, and sparse detector arrangement for which the sensor positions would be optimised for a given treatment plan will likely give the highest accuracy in localising the BP.
Furthermore, the sensor positioning should also consider the specific constraints of a given treatment plan. In particular, the locations where the beam enters the patient should be detector free to prevent the influence of the usually high-density, high-Z (atomic number) piezoelectric sensors on the beam quality and degradation of the IA detector performance. In the SIRMIO project, the small animal will be maintained in a sterile environment during the irradiation using a specially designed holder . The IA sensors should ideally be in direct contact with the animal skin to maximise the SNR of the detected signals. Therefore, the detectors should be positioned before enclosing the sterile space and setting it onto the SIRMIO platform. Bulky commercial ultrasonic detectors, as used in this study, cannot be inside the proton imaging FOV. Indeed, their large water equivalent thickness (WET >> 1 cm) and atomic number Z would considerably increase proton scattering and consequently reduce the image resolution. This additional constraint further limits the locations where the IA detector can be positioned. As an alternative, sensors with a low material budget could be employed. Emerging ultrasonic transducer technologies such as micromachined transducers (Lascaud et al 2019) and optical hydrophones (Sueyasu et al 2022), or more conventional piezoelectric polymers with a dedicated acoustic design (Lascaud et al 2021c) that have been proposed in the recent years for ionoacoustics are possible candidates. However, further investigation is needed to determine whether such a sensor could be sensitive enough to detect the very low pressure of a few mPa (i.e. at least one order of magnitude lower than in clinical scenarios) expected in SIRMIO after degradation of the clinical proton beam (Lascaud et al 2022).

Transition to clinical applications
Assuming optimal synchronisation of the signal acquisition, proper speed of sound correction, and accurate positioning of the sensors, multilateration is fundamentally limited by the capability to extract the ToF. As illustrated by our results, the shape of the signal detected varies depending on the sensor position. This phenomenon is particularly enhanced at pre-clinical energies for which the proton range is relatively small compared to the propagation length and sensing aperture. In these conditions, the signals detected laterally are the superposition of the acoustic waves emerging from the BP, the plateau region, and potentially the energy discontinuity at the phantom entrance, which leads to a large error in the estimated ToF. Illustratively, for clinical applications, the proton range increases from about 4 mm at 20 MeV to 122 mm at 130 MeV in water, allowing for a better separation between the different wavefronts. The error on the ToF extraction for a setup similar to the one used in this study in the case of 130 MeV proton beam delivered at clinical synchrocyclotron accelerator is presented in figure S11b of the supporting material. The maximum error in the ToF is 0.75 μs, corresponding to a distance of 1.11 mm which is in good agreement with previous studies , Lehrack et al 2017. This gives confidence that the knowledge gained at low proton beam energies can be applied to clinical scenarios.
It should be noted that the irradiation of a patient will give rise to additional acoustic emissions at the interfaces between the different tissues due to energy discontinuities and variation of the Grüneisen parameter. However, these signals of rather high spatial frequencies are expected to be filtered out when irradiating soft tissues by the microsecond proton pulse clinically available. Therefore, for favourable locations without bones or air cavities in the acoustic and beam path, the millimetre accuracy seems achievable, as suggested by a recent experimental study in an anthropomorphic phantom (Patch et al 2021).

Conclusion
Acoustic localisation of the Bragg peak position in water was investigated in-silico and experimentally for preclinical scenarios. The robustness of two multilateration algorithms (TOA and TDOA) was compared, as well as the influence of the sensor position with respect to the proton beam on the localisation accuracy. The simulation study showed that the position of the Bragg peak in water could be located in two dimensions with an accuracy of 90 μm (2% error) when the sensors are optimally positioned to minimise the error on the time-of-flight estimation. Experimentally, the localisation error was found to vary between 0.4 mm to 1 mm depending on the sensor arrangement and beam energy, which was attributed to inaccurate knowledge of the sensor positions and/or intrinsic low signal-to-noise ratio resulting from the weak ionoacoustic emissions and their directivity. Further studies are required to assess the robustness of the multilateration method during the delivery of a (pre-) clinical treatment plan for which the moving position of the Bragg peak relative to the sensor array may result in a large variation of the acoustic localisation accuracy. Moreover, high doses were used during our experiments to improve the signal-to-noise ratio of the acquired ionoacoustic signals, this was done to investigate the method's performance under ideal statistical conditions, which could in future be achieved with improvements in detector technologies.