Since the advent of laser-ion acceleration in 2000, the field has been steadily progressing^[1^–4^]. While previously ions were accelerated predominantly by target normal sheath acceleration (TNSA) to up to tens of MeV, state-of-the-art PW systems with improved temporal laser contrast reach new acceleration regimes. As such, radiation pressure acceleration and relativistic-induced transparency accelerate ions within sub-ps laser–plasma interactions, resulting in pronounced angular emission characteristics^[5^]. At these acceleration parameters, protons approaching the 100 MeV kinetic energy barrier as well as carbon ions exceeding 80 MeV/u have been observed^[6^–8^]. With improving parameters of this technique, the number of applications is increasing steadily. As such, the first radio-biological in vivo studies of tumor irradiation with laser-accelerated ions investigating the FLASH effect have been conducted^[9^], fundamental ion–matter interactions in high-energy deposition regions have been investigated^[10^–12^] and the application of laser-accelerated ion bunches for fuel ignition in inertial confinement fusion has been discussed^[13^–15^]. All of these applications require a precise ion spectrum determination for depth–dose verification, model input and cross-section maximization. Peak beam currents exceeding mega-amperes pose challenges to detectors applied for spectral measurements^[16^]. The harsh conditions encountered, especially in the aperture angle close to the laser–plasma interaction where strong electro-magnetic pulses (EMPs) disturb (or even damage) electronics, prevent the application of complementary metal oxide semiconductor (CMOS) detectors, such as RadEyes^[17^,18^]. On the other hand, when measuring further downstream of the laser–plasma interaction as, for example, done by Thomson parabolas, angular emission features such as energy-dependent cone narrowing are forfeited^[19^,20^]. Other non-invasive detection modalities, such as integrating current transformers, measure the total beam charge without a precise spectrum reconstruction^[21^]. Due to these issues, non-electronic methods such as radio chromic films (RCFs) or Columbia Resin #39 (CR39) are still utilized as widespread standard ion diagnostics in laser-ion-acceleration experiments^[22^,23^]. They suffer from none of the discussed drawbacks, but nevertheless they do not allow for online read-out. Particularly for highly energetic particles exceeding 100 MeV, the RCF stacks used in experiments can become quite thick and evaluation is tedious.

A rather new approach for energy density determination of laser-accelerated particles is ionoacoustics^[24^]. It capitalizes on the measurement of the ultrasound wave emitted by the heat-expansion of a medium when energy is deposited by ions in a short amount of time^[25^,26^]. The pressure trace $p\left({\mathit{\boldsymbol r}},t\right)$ at an ultrasonic transducer at position $\mathit{\boldsymbol r}$ , emerging from an arbitrary dose distribution $H\left(\mathit{\boldsymbol r}^{\prime },t\right)$ , can be described by the following:

$$\begin{array}{r}p(\text{boldsymbol}\mathit{r},t)=\frac{N_{\mathrm{p}}}{4\pi }\frac{\mathrm{\Gamma }}{c_{\mathrm{s}}^2}\frac{\partial }{\partial t}\int \mathrm{d}{\text{boldsymbol}\mathit{r}}^{\prime }\frac{1}{\mid \text{boldsymbol}\mathit{r}-{\text{boldsymbol}\mathit{r}}^{\prime }\mid }H({\text{boldsymbol}\mathit{r}}^{\prime },t-\frac{\mid \text{boldsymbol}\mathit{r}-{\text{boldsymbol}\mathit{r}}^{\prime }\mid }{c_{\mathrm{s}}}),\end{array}$$  (1)

The ion-bunch energy acoustic tracing (I-BEAT) detector measures this far-field pressure trace on a single-shot to reconstruct the initial energy deposition and the corresponding particle spectrum down to a minimum of ${10}^7$ protons/mm² per bunch^[27^]. The advantages of this method are, firstly, the EMP resistance as the pressure wave has a μs delay to the EMP when arriving at the ultrasound detector. The detector is a piezoelectric PZT-transducer with a saturation threshold orders of magnitude above the capabilities of modern laser accelerators. Secondly, as radiation damage or decay in water is negligible, it is radiation resistant in terms of re-usability. Further, water shows a low acoustic attenuation, allowing one to detect signals nearly unperturbed.

Demonstrated by several experiments, this ionoacoustic approach is also capable of detecting mono-energetic bunches of 20 MeV protons down to minimum energy deposition of ${10}^{12}$ eV, equivalent to a particle number of $5\times$ ${10}^4$ with an accuracy in the range verification of $\pm 30\;\unicode{x3bc}$ m^[28^]. Similarly, ¹H, ²³⁸U, ¹²⁴Xe and ¹²C ions in excess of hundreds of MeV/u kinetic energies were characterized and applicability up to 1 GeV was proposed^[29^,30^]. Applying multiple transducers, it has been shown that ionoacoustic tomography of the ion beam’s Bragg-peak profile is possible in vivo and in real-time, with sub-mm accuracy^[31^].

Adversely, the inherent drawback of this method is the reliance on a spatial energy density gradient to produce the pressure amplitude. This gradient is induced by the Bragg-peak for mono-energetic ion beams. Laser-accelerated ion beams, however, typically exhibit an exponential ion energy spectrum with certain additional features, depending on the acceleration regime^[1^,5^,7^]. These spectra range over orders of magnitude and the corresponding energy densities deposited in matter are typically dominated by the low-energy particles. Therefore, the dynamic range of ionoacoustic measurements often is insufficient to measure and reconstruct the entire spectrum and is sensitive only to their low-energetic components^[27^,32^].

Here we present a simulation study of a new version of I-BEAT, solving this problem by artificially introducing gradients to the energy density distribution, that is, tracing ionoacoustic modulations of broad energy distributions (TIMBRE). The detector is designed to be placed within a few centimeter range behind the target to collect the majority of the accelerated particles. Similar to the spherical ionoacoustic waves with resonant frequency (SPIRE) technique^[33^,34^], we use special modulator foils to increase the electronic stopping power of the ions. This generates high-energy density gradients at the water–modulator interfaces and, hence, strong broadband frequency signals. In addition, placing the modulator foils equidistantly results in particularly strong resonant waves, allowing for reconstruction of the energy deposited in the modulator foils. This overcomes theoretical predictions of minimum detection thresholds of the Bragg-peak, given that media interfaces are spatially more pronounced in their electronic stopping gradients^[35^].

Figure 1(a) shows a cross-section scheme of the detector. It consists of a vacuum-compatible water tank that is placed within the centimeter range behind the laser target to collect most of the accelerated ions. A Kapton entrance window is used for the ions to enter the water volume. Inside the tank, multiple flat foils act as modulators to shape the dose distribution. By this method, spatially separated energy density distribution characteristics are encoded in the acoustic source signal. To find an optimal modulator material four quantities were optimized. First, the ion stopping power in the modulators should be high compared with water to generate large energy density gradients at the interfaces. Second, the modulator’s Grüneisen parameter $\Gamma$ (see Equation (1)), which defines the material’s property to translate a certain energy density into pressure, should be high to further increase the instantaneous pressure ${p}_0\equiv p\left(\mathit{\boldsymbol r},0\right)$ . Third, the acoustic impedance of the material ${z}_\mathrm{mod}$ needs to be matched to water, given that the reflectivity $R$ of a pressure wave is given by $R=\left({z}_\mathrm{w}-{z}_\mathrm{mod}\right)/\left({z}_\mathrm{w}+{z}_\mathrm{mod}\right)$ , where ${z}_\mathrm{w}$ is the water impedance ( $z=\rho {c}_\mathrm{s}$ ). Thus, pressure signals that have to percolate subsequent foils to propagate to the transducer are transmitted without extensive reflection of their energy content. Lastly, the speed of sound of the material should be as low as possible to keep generated frequencies $f$ low, because higher ones are attenuated more strongly in water, which shows a quadratic acoustic attenuation of $0.0022\;\mathrm{dB}\cdot\mathrm{cm}^{-1}\cdot{\mathrm{MHz}}^{-2}$ . Here, lead was found to show the best trade-off between all four conditions due to its balance among the speed of sound ${c}_\mathrm{s}$ , density $\rho$ and electronic stopping. The initial pressure in the lead modulators ( ${p}_{0,\mathit{\operatorname{mod}}}$ ) is more than two orders of magnitude higher than in water ( ${p}_{0,\mathrm{w}}$ ), that is, ${p}_{0,\mathit{\operatorname{mod}}}\approx 256\times {p}_{0,\mathrm{w}}$ , resulting from increased electronic stopping and a higher Grüneisen parameter, whilst the reflectivity $R$ at the interfaces is approximately 87%. Figure 1(b) shows a representative simulated depth–dose distribution for the reference case with 10 modulator foils and an input proton spectrum as accelerated by the petawatt high-energy laser for heavy ion experiments (PHELIX) laser at the GSI Helmholtzzentrum für Schwerionenforschung^[32^].

In Figure 2, the process of the pressure wave generation is shown. The initial pressure has a step gradient at the modulator interfaces, which comes with a wide frequency bandwidth in Fourier space. Off-resonant components of this frequency band interfere destructively in subsequent oscillations and frequencies close to the resonance ${f}_\mathrm{res}={nc}_{\mathrm{s},\mathit{\operatorname{mod}}}/\left(2{d}_\mathrm{f}\right)\;\left(n\in \mathrm{\mathbb{N}}\right)$ , defined by the modulator’s thickness ${d}_\mathrm{f}$ , are enhanced. This is represented in Figure 2(a) for $t\times {f}_\mathrm{res}\approx 1$ and $n=1$ . With each oscillation of these resonant standing waves within the modulator, a fraction of the initial amplitude is transmitted into the water volume. Figure 2(b) shows the corresponding spectra of the signals. Thus, a pulse train is generated with a carrier frequency at ${f}_\mathrm{res}$ and an envelope function according to the build-up and ring-down defined by the modulators. The entire energy contained in these oscillations is proportional to the spectral amplitude, which is increasing for ${f}_\mathrm{res}$ as the spectrum narrows down. Here, the first-order ${f}_\mathrm{res}={c}_{\mathrm{s},\mathit{\operatorname{mod}}}/\left(2{d}_\mathrm{f}\right)\;\left(n=1\right)$ contains most of the energy with higher orders rapidly decreasing in energy content, that is, spectral amplitude. Therefore, in the following only the first ( $n=1$ ) resonance will be considered.

Transmitting through a subsequent modulator foil will only weakly reduce the spectral amplitude of the pulse trains. This is due to the phase difference of the pulse train’s jth wave-cycle reflection ${R}_{j}$ from the back surface relative to the $\left(j+1\right)$ th wave-cycle reflection ${R}_{j+1}$ from the front interface. Due to the phase-shift at the backside reflection, the phase difference is exactly $\pi$ , meaning they interfere destructively. Since the initial amplitude of ${R}_j$ is larger compared to ${R}_{j+1}$ due to the ring-down shape of the pulse trains, their destructive interference is enhanced. Similar to the standing waves shown in Figure 2 (green), each oscillation of the incoming pulse train will excite a new pulse train. Superimposing all of them yields the transmitted pulse, which will be spread temporally compared to the incoming one, that is, its spectrum will be peaked more narrowly around ${f}_\mathrm{res}$ with a higher spectral amplitude. As a result, the transmission of a resonant pulse train through subsequent modulator foils will reduce its energy content by less than 20% in this setting. An ultrasonic piezoelectric transducer (PZT) that relies on the same mechanism is used to measure the signal. The thicknesses ${d}_\mathrm{f}$ of both the Kapton entrance window and the modulator foils were designed such that their fundamental ( $n=1$ ) resonance frequencies ${f}_\mathrm{res}$ coincide with the transducer’s response peak, here equal to 9.8 and 10 MHz, respectively.

For all simulations a reference differential TNSA angular and energy proton spectrum with a maximum energy of 60 MeV was assumed^[32^]. Both energy and angular distributions were included in a FLUKA source file to calculate the energy distribution in the detector^[36^,37^]. The detector was placed 6 cm behind the particle source. Considering the 1.7 cm thickness of the front flange, this defines the source to entrance window distance as 7.7 cm. The entire detector was modeled as shown in Figure 1(a), with an additional 2 cm thick plastic shielding (entrance window spared) and a 15 μm thick aluminum foil in front of the flange, to stop all ions of kinetic energies below approximately 2 MeV and reduce background noise by energy deposited in the detector’s casing.

Subsequently, the geometry and the deposited energy data were exported to a k-wave code to simulate the resulting pressure evolution within the detector using the deposition data according to Equation (1) as the source^[38^–40^]. For pressure measurement, an Olympus Videoscan V311-SU transducer with a resonance frequency of 10 MHz was modeled, while its electrical impulse response function was neglected.

For the reconstruction of the energy density from a pressure trace, a semi-3D algorithm was specially designed. It requires only the geometry and a measured pressure trace $p(t)$ as input. In the first calculation step, the geometry is reduced to a 1D-array in the direction $x$ normal to the transducer surface, that is, axial to the tank. On this array, for every interface point, all possible propagation paths for a normalized source signal to the transducer are calculated, including reflection and transmission coefficients up to a certain threshold. This is done by multiplying the signal alternating with a propagation matrix and a transducer array, similar to the ABCD-matrices in ray optics. Next, all the detector signals with the same origin are stitched to one array $S(x)$ , which is then scaled accordingly. Scaling happens by the acoustic attenuation $\mathrm{att}(x)$ and the fraction of the signal that actually reaches the transducer $\mathrm{dif}(x)$ when the normalized source signal is diffracted either in two or three dimensions. For the diffraction, the Fresnel-integral is solved. To cover the transverse dimension of the signal shapes, all direct modulator signals $C$ in $p(t)$ are isolated, normalized and convoluted into $S(x)$ to create an array ${S}_\mathrm{c}(x)=S(x)\ast C$ that contains all contributions originating from a single modulator foil. Thus the entire signal at the transducer can be expressed as follows:

$$\begin{array}{r}p(t)=\sum _i{S}_{\mathrm{c},i}(x)\times \mathrm{dif}(x)\times \mathrm{att}(x)\times p_{0,i}.\end{array}$$  (2)

$$\begin{array}{r}{\epsilon }_{\mathrm{rec},i}=max\left\{\frac{\mathcal{F}[S_{\mathrm{c},i}(t)\times p_{0,i}]}{{\mathrm{\Gamma }}_{\mathrm{mod}}}\right\}.\end{array}$$  (3)

The average reconstruction time of a depth–dose profile scales linearly with the temporal window that is considered for the isolated signals $C$ . Considering the ring-down function for lead modulator foils, a time frame of 2 μs was found to work best for the 10 modulator foils case. The temporal windows embedded the signals 0.3 μs prior and 1.7 μs post to their arrival timings at the transducer. Here the mean reconstruction time was 0.42 s, which can be considered online given that most state-of-the-art high-power laser systems operate in the Hz repetition rate regime.

Figure 3(a) presents the pressure trace $p(t)$ that results from the reference source distribution in Figure 1(b) (blue ). Here, 10 modulator lead foils separated by 3 mm were used. The resonance frequencies of the entrance window, modulators and transducer response were tuned to approximately 10 MHz. The last modulator foil with the signal occurring at 46 μs was chosen right behind the Bragg-peak volume of the highest energetic ions ( $\approx$ 60 MeV). For comparison, the pressure signal ${p}_\mathrm{ref}(t)$ for the same ion spectrum without any modulator foils in the water reservoir is shown as well (orange). The Fourier-transformed signals $p(f)=\mathrm{\mathcal{F}}\left[p(t)\right]$ in Figure 3(b) reflect the corresponding spectra. In the modulated case, three features can be distinguished. Firstly, the overall modulation spanning the entire spectrum is produced by the overlap of individual, temporally separated pulse trains and consequently shows a comb with $\Delta f={c}_\mathrm{s}/\left({md}_\mathrm{f}\right)$ , where $m\in \left[1,10\right]$ describes the space numbers in between foils. The second feature is the low frequency peak at $f<$ 2 MHz. This signal is generated by the overall broad deposition region excluding the modulations, that is, a DC background. It contains wavelengths ${\lambda}_\mathrm{DC}$ that are much larger than the thickness of the foils and the transducer’s coating and therefore are not suppressed. Lastly, the actual first- ( $n=1$ ) and second-order modulator ( $n=2$ ) signals are at 9.8 and 19.6 MHz, respectively. The first-order signal contains the majority of the radiated energy and is the basis for the reconstruction.

Comparing both traces in Figure 3, modulated and unmodulated, shows the inherent asset of this detector approach. In the temporal traces it can be seen that the dynamic range in the unmodulated case (orange) is more than ${10}^6$ , that is, it is not possible to reconstruct the signal generated by the initially high energetic ions since the entrance window signal dominates the trace. In contrast, the modulated curve (blue) shows a total dynamic range of ${10}^4$ and for the isolated modulator signals only ${10}^3$ . Thus, common ultrasound transducers can record the entire signal with a high-pressure resolution. Further, from the spectra it becomes clear that the modulation causes a sharp narrowing of the frequency band around the resonance, while the total energy content is slightly increased (spectral amplitude). This may be due to a diffraction effect, as the foils predominantly emit in the axial dimension, whereas the unmodulated instantaneous pressure shows a larger transverse component.

Figure 4(a) shows the energy density reconstructed from the modulated pressure trace in Figure 3 with an SNR of ${10}^8$ . The agreement between ground truth input and reconstruction is excellent. The ratio of ${\varepsilon}_\mathrm{rec}$ to the input values of ${\varepsilon}_\mathrm{in}$ equals $1\pm 0.003$ for all 11 data points. Figure 4(b) shows this ratio $q={\varepsilon}_\mathrm{rec}/{\varepsilon}_\mathrm{in}$ for each modulator position when artificially reducing the SNR. Because the added noise applies to the complete monotonically decaying curve, the SNR is smaller for a larger depth. Therefore, the SNR was normalized to the signal level at the respective modulator. The reconstruction yields very good results down to a normalized SNR $\gtrsim$ 10. For SNR $<$ 10, the noise amplitude impacts the optimization of ${p}_0$ . As explained above, the signal duration of individual foils increases upon transmission through subsequent foils because their spectrum is narrowed, while their spectral amplitude is reduced only minimally. Regarding the temporal trace this means that the peak-to-mean ratio decreases and hence is lower for modulator signals originating from shallower foils. Whilst in the reconstruction the mean is important, also the peaks contribute notably, whereas the normalized SNR only considers the mean. Therefore, a trend is observable of shallower foils with a worse peak-to-mean ratio starting to increase in the performance ratio at higher SNR values compared to the signals originating from deeper modulators. For noise levels equal to the signals, the optimization predominantly starts to reconstruct the noise instead of the pressure trace. Here $q$ starts to diverge rapidly. For the entrance window this effect sets in at SNR $\approx$ 0.35, whereas the last modulator shows this behavior for SNR $\approx$ 0.08. Hence, the sensitivity of the method to the high-energy part of the spectrum increases by almost an order of magnitude.

The assignment of a reconstructed energy density to a specific modulator foil limits the spatial resolution to the thickness of this foil, that is, 100 μm, which is similar to what an integrated RCF stack provides. However, the intermediate signal originating from in between foils cannot be reconstructed. Stacking modulators closer and increasing their number will shift the $q$ –SNR curve of shallower sources to higher SNR values (as explained above). Hence, there is a trade-off between the reconstruction of arbitrary distribution sample points and reconstruction performance based on the SNR level. A random noise signal in time as introduced in this study also shows a random noise distribution in Fourier space. This allows for frequency filtering during post-processing to reduce the noise contribution in off-resonant parts of the spectrum. Figure 4(c) exemplifies the result of this method for a Gaussian filter around ${f}_\mathrm{res}$ with a full width at half maximum (FWHM) of $\Delta f=$ 2 MHz. The overall shift to lower SNR values by approximately ${10}^2$ suggests a clear improvement of the reconstruction algorithm as the noise signal is reduced by a fraction of approximately $\Delta f/{f}_{\mathrm{max}}=\Delta f\times \mathrm{d}t/2$ , where $\mathrm{d}t$ is the temporal resolution of the measurement. In addition, the performance ratio $q$ of deeper modulators starts to increase for higher SNR values as compared to signals originating from more shallow foils. It is assumed that this effect and the higher variability are caused by distortions due to the interference between the remaining resonant noise and the actual signal, since the envelope is used for optimization. This overall worsens reconstruction performance for temporally shorter modulator signals or smaller numbers of foils.

The simulation study suggests that TIMBRE is capable of measuring broadband, laser-accelerated ion bunches. The increased energy density gradient and the enhancement at the resonance frequency due to the modulation of the water reservoir increase the sensitivity by at least $2.5\times {10}^2$ and approximately ${10}^3$ with additional frequency filtering. This suggests that less than ${10}^5$ protons/mm² even in low-ionizing regions can be detected. An earlier version without modulators has been successfully tested in close vicinity to the laser target and showed little disturbance to EMPs or other noise sources^[32^].

In particular, due to its capability of immediate feedback and simplicity, we expect the presented detector approach to become a valuable addition to other diagnostics, if not a primary ion diagnostic for laser–plasma-based ion sources, once tested extensively in experiments.