A spatial and temporal characterisation of single proton acoustic waves in proton beam cancer therapy

: An in vivo range veriﬁcation technology for proton beam cancer therapy, preferably in real-time and with submilli-meter resolution, is desired to reduce the present uncertainty in dose localization. Acoustical imaging technologies exploiting possible local interactions between protons and microbubbles or nanodroplets might be an interesting option. Unfortunately, a theoretical model capable of characterising the acoustical ﬁeld generated by an individual proton on nanometer and micrometer scales is still missing. In this work, such a model is presented. The proton acoustic ﬁeld is generated by the adiabatic expansion of a region that is locally heated by a passing proton. To model the proton heat deposition, secondary electron production due to protons has been quantiﬁed using a semi-empirical model based on Rutherford’s scattering theory, which reproduces experimentally obtained electronic stopping power values for protons in water within 10% over the full energy range. The electrons transfer energy into heat via electron-phonon coupling to atoms along the proton track. The resulting temperature increase is calculated using an inelastic thermal spike model. Heat deposition can be regarded as instantaneous, thus, stress conﬁnement is ensured and acoustical initial conditions are set. The resulting thermoacoustic ﬁeld in the nanometer and micrometer range from the single proton track is computed by solving the thermoacoustic wave equation using k -space Green’s functions, yielding the characteristic amplitudes and frequencies present in the acoustic signal generated by a single proton in an aqueous medium. Waveﬁeld expansion and asymptotic approximations are used to extend the spatial and temporal ranges of the proton acoustic ﬁeld. V C 2022 Author(s). All article content,


I. INTRODUCTION
Proton therapy is an alternative for conventional photon therapy in radiation oncology.Proton beams have a theoretical advantage over photon beams by having a more localised dose distribution characterised by a peak at the end of their range, the Bragg peak, behind which no dose is present.Therefore, a more conformal therapeutic dose can be applied to the tumour while sparing the surrounding healthy tissue.Realising this advantage in the clinic relies on an accurate positioning of the Bragg peak within the tumour volume, which is a challenging demand.In particular, range uncertainties are characteristic of proton therapy.These are considered to be the greatest issue in current clinical practice and have induced an urgent need for accurate in vivo range verification techniques (Paganetti, 2011;Parodi and Polf, 2018).
An envisaged way of improving dose deposition accuracy in proton therapy is to measure the location of proton dose deposition within the anatomy of the patient, preferably in real-time and with a sub-millimeter resolution.Monitoring the actual delivered dose during or after the treatment provides the opportunity to either modify the treatment or compensate in later stages in case of deviations from the original treatment plan.Several methods for dosimetry in proton therapy are under investigation at the clinical or experimental level (Knopf and Lomax, 2013;Parodi, 2020), including prompt gamma imaging (Hueso-Gonz alez et al., 2015; Min et al., 2006), positron emission tomography (Parodi, 2015;Parodi et al., 2007), and ionoacoustic imaging (Hickling et al., 2018).In positron emission tomography and prompt gamma imaging, the correlation between the emerging nuclear induced secondary radiation and the Coulomb-induced dose deposition is not straightforward, so that a Bragg peak positioning accuracy better than a few millimeters is highly challenging in clinical situations.
Ionoacoustic imaging employs acoustic pressure signals generated by thermoelastic expansion of the confined area where protons deposit the majority of their energy as heat (Hickling et al., 2018).The frequency content of the emitted acoustic pressure wave is determined by the spatial dimensions of the heated area and the temporal characteristics of the heat deposition by the packet of protons.The spatial dimensions are mainly dictated by proton stopping physics while the temporal characteristics are mainly a result of the pulse characteristics of the proton beam.Stopping of protons is a stochastic process which makes that not all protons have the same range.While the axial length of the Bragg peak of an individual proton is several micrometers, the stochastic nature of proton stopping causes the axial dimension of the Bragg peak for a proton beam (consisting of roughly half a million protons per pulse for a 10 nA beam intensity and 10 ls pulse duration) to spread out to several millimeters.Consequently, ionoacoustic signals in clinical environments have frequencies on the order of kHz, hampering the desired submillimeter localization resolution.
To approach the problem from another perspective, one may consider the heat deposition of an individual proton.As this heat is deposited over a range of several micrometers, an impacting proton generates a broadband acoustic pressure wave.We will refer to this wave as the proton acoustic signal, in order to distinguish it from the previously mentioned ionoacoustic signal.The higher frequency of this proton acoustic signal will allow a higher localisation resolution than the ionoacoustic signal.However, the proton acoustic pulses cannot be detected at the surface of the body because their high frequency will give rise to very high acoustic attenuation.In the absence of these effects, the ionoacoustic signal would contain much higher frequencies than mentioned above.
Clinical ionoacoustic signals typically have a peakpressure on the order of tens of mPa in the kHz band (Jones et al., 2015).Averaging is necessary to detect these, as the state-of-the-art medical transducers have a noise-equivalent pressure on the order of hundreds of mPa (Matte et al., 2011;van Neer et al., 2010).Thus, it will now be clear that there is a challenge of measuring proton-induced pressure pulses at the surface of the human body in a higher frequency band, where the pressure content is even lower.To tackle this challenge, we assume that ultrasound contrast agents (UCAs) consisting of micron-sized bubbles may be useful (Lascaud et al., 2021).These microbubbles can be injected in the blood flow of patients and act as acoustic scatterers.If a single proton impacts in the vicinity of a microbubble, the proton acoustic wave might drive the microbubble into oscillation at its resonance frequency, which is typically on the order of a few megahertz.Moreover, monodisperse bubbles that are excited simultaneously might oscillate synchronously at the same resonance frequency.In this way, microbubbles might convert the broadband proton acoustic pulses into acoustic waves with amplitudes and frequencies detectable outside the body and with frequencies allowing submillimeter resolution.This idea was substantiated by experiments using broadband photoacoustic pulses with center frequency in the hundreds of MHz band, which drove microbubbles into oscillation at their natural frequency in the much lower MHz range (Daeichin et al., 2020;Lum et al., 2018).A similar concept could be based on the vaporization of nanodroplets by a proton beam, which has been recently observed in experiments (Carlier et al., 2020).
Since microbubble oscillations are driven by the pressure difference in the gas-liquid interface, the acoustic fields of single protons play a key role in the interaction with microbubbles and nanodroplets.The temporal evolution of the signal is relevant in the nanometer and micrometer distances, as this is the distance expected between the proton stopping position and the closest microbubble.To numerically simulate these interactions, models for the proton acoustic wave generation and the microbubble response to these waves are necessary but have not yet been developed.Until now, the pressure generation in proton beams has been described as a bulk phenomenon where in a certain region of interest the thermoacoustic effect of all protons in the beam was considered, resulting in the ionoacoustic signal.Baily (1992) has reviewed theoretical and experimental work on the generation of acoustic waves through the interaction between ionizing radiation, such as protons and materials.It was proven that the ionoacoustic signal generated by a proton beam fits with a theory based on the thermoacoustic expansion of the region where the proton beam deposits the majority of its energy.However, the properties of the acoustic waves generated by the individual protons in the beam are still underexposed.To bridge this theoretical gap, we have developed a framework that combines a single proton heat deposition model with numerical wave simulation algorithms, such that acoustic waves generated by single protons could be computed.Importantly, our model does not account for the possibility of a proton stopping within the gas core of a microbubble, as this is considered a highly unlikely event.An envisaged next step in this simulation study is to compute the microbubble response to the computed single proton acoustic waves, but the modelling of this physical process for this paper out of scope.
Localized heat deposition in a medium that is irradiated with swift ions has been studied for a variety of applications including proton beam cancer therapy (Toulemonde et al., 2009).In previous work, the inelastic thermal spike model has been used to study the temporal temperature distribution along the track of an individual proton stopping in water.The model that we will use inherently treats the proton heating process as a two step process, first the absorption of proton energy by electrons, creating secondary electrons or delta rays, and second the heating of the molecules through electron-phonon coupling.In an earlier study temperature spikes and pressure spikes in irradiated water have been discussed (Sun and Nath, 1993), but here all the energy lost by the ions was assumed to be immediately transferred to the molecules, thus neglecting the electronic component, resulting in the prediction of much higher temperatures and consequently pressure spikes.
We will use methods originating from the field of photoacoustics to predict the acoustic field generated due to thermoelastic expansion following localised radiation heating.The thermoacoustic wave equation will be solved in k-space using a wavenumber integration algorithm, yielding the proton acoustic wave (Cox and Beard, 2005).In addition, spatial and temporal wave extrapolation algorithms https://doi.org/10.1121/10.0009567will be used to stretch the computational limits on simulation ranges.In this way, the characteristic amplitudes and frequency content of the proton acoustic wave were obtained, providing a theoretical basis for continuing research on the behaviour of ultrasound contrast agents and nanodroplets in a proton beam.

II. NUMERICAL METHODS
As described above, in interactions between protons and the medium they are traversing, kinetic energy is converted to heat.As a result of this, a temperature distribution is induced along the proton track, which generates an acoustic wave through local thermoelastic expansion of the medium.Models for the heating process and acoustic wave generation are described in Secs.II A and II B, respectively.All presented calculations have been performed for a proton in water, but the framework can straightforwardly be extended to other media or ion species (Toulemonde et al., 2009;Walig orski et al., 1986).

A. Proton heat deposition model
A proton traversing a medium mainly slows down by transferring kinetic energy to electrons, creating a cascade of secondary electrons.These energetic secondary electrons heat the medium when becoming thermalized and bound.This radiation heating process can be modelled phenomenologically as a two-step process: (1) the absorption of proton energy by electrons and (2) the heating of the atoms through electron-phonon coupling (Qiu and Tien, 1992).Previously, for a proton in water this process was mathematically described by an inelastic thermal spike model consisting of two coupled equations of energy transfer (Toulemonde et al., 2009), given by where T e ðr; tÞ and Tðr; tÞ denote temperature increases of the electronic and molecular subsystems resulting from proton impact.The symbols C e , C(T), K e , and K(T) denote the specific heats and thermal conductivities of the electronic and molecular subsystems (Dufour et al., 1993;Toulemonde et al., 1996;Toulemonde et al., 2009).The electron-phonon coupling constant g is linked to the electron-phonon mean free path l by the relation g ¼ K e =l 2 (Meftah et al., 2005;Toulemonde et al., 1996;Toulemonde et al., 2000;Toulemonde et al., 2009).The symbol Dðr; tÞ denotes the energy density, or dose, supplied by the incident proton to the electronic subsystem and carried via the electron cascade over a distance r ¼ jrj from the proton track (Gervais and Bouffard, 1994;Toulemonde et al., 1996;Walig orski et al., 1986), where r is the position vector in a three-dimensional cylindrical domain with the longitudinal axis aligned with the proton track.Cylindrical symmetry is imposed such that the quantities are independent of the azimuthal coordinate.Equation ( 1) was solved numerically for the temperature distribution Tðr; tÞ resulting from proton impact.An overview of the parameters in the inelastic thermal spike model is displayed in Table I.

B. Proton acoustic wave model
When a region of a fluid is heated, a sound wave will be generated through thermoelastic expansion.The relation between the proton-induced temperature distribution Tðr; tÞ and the resulting proton acoustic wave pðr; tÞ is described by the thermoacoustic wave equation.
We consider a homogeneous, linear, and lossless acoustic medium, albeit that the signals have a high frequency.In general, the acoustic pressure pðr; tÞ generated through thermoelastic expansion obeys the following wave equation: Tðr; tÞ @t 2 ; (2) with pðr; t ¼ 0Þ ¼ 0 and @pðr; tÞ @t t¼0 ¼ 0: (3) The symbol t denotes the time coordinate, and r 2 indicates the Laplace operator.Tðr; tÞ denotes the temperature, v s is the ambient speed of sound, b V is the volume thermal expansivity, and j is the isothermal compressibility.Implicitly, it is assumed that thermal conductivity during heat deposition can be neglected, which condition holds for sufficiently short heating pulses and is known as thermal confinement (Wang, 2009).Proton heat deposition occurs on a timescale much shorter than the characteristic acoustic travel time of the resulting wave, thus may be regarded as instantaneous.Therefore, it may be assumed that all heat energy has been deposited before the medium has expanded and its equilibrium mass density q has changed.This condition is known as stress confinement and holds if the heating pulse duration t 0 satisfies t 0 < d c =v s with d c being the characteristic length of the temperature distribution.The concept of stress confinement originates from the field of photoacoustics, where short laser pulses are used as a heating source.In Sec.III, it is proved that the condition for stress confinement is satisfied for single proton acoustics.Under isochoric conditions, TABLE I. Parameters for the inelastic thermal spike model in water (Toulemonde et al., 2009).
the change in density Dqðr; t 0 Þ, just after the proton heating pulse, is related to the change in temperature Tðr; t 0 Þ and change in pressure pðr; t 0 Þ via the thermodynamic equation Dqðr; t 0 Þ ¼ qjpðr; t 0 Þ À qb V Tðr; t 0 Þ.Setting Dqðr; t 0 Þ to zero, which is dictated by stress confinement, yields a relation for the acoustic pressure at time t 0 , Equation ( 4) should be in agreement with the solution pðr; tÞ of Eq. ( 2) at the time instance t ¼ t 0 .Intuitively it can be deduced that, whereas Eq. ( 2) included a source term, in the case of instantaneous proton heating, the problem can be recast as an initial value problem with no explicit source term but with a given pressure distribution at the instant of the proton heating pulse pðr; t 0 Þ.This makes the two initial conditions required for a unique solution explicit.
Mathematically, the initial value problem for the specific case of proton acoustics is then defined as r 2 pðr; tÞ À 1 v 2 s @ 2 pðr; tÞ @t 2 ¼ 0; (5) with pðr; t ¼ t 0 Þ ¼ p 0 ðrÞ and @pðr; tÞ @t The first condition in Eq. ( 6) defines the acoustic pressure distribution at a time t 0 after proton impact, when the temperature increase has reached its maximum value.The second initial condition is equivalent to zero particle velocity ṽp ðrÞ everywhere at t 0 , which can be safely assumed.Solving Eq.
(5), using these initial conditions and the acoustic parameters from Table II, yields the proton acoustic wave pðr; tÞ in water.Three algorithms were used for full spatiotemporal characterisation of the proton acoustic field, which will be introduced in the following.The algorithm in Sec.II B 1 calculates the proton acoustic field in a three-dimensional Cartesian domain for a single time instant.As such, spatial characteristics of the field can be visualised efficiently but studying dynamic behavior requires successive evaluation of the algorithm for a series of time instants and therefore is, particularly for large domains, computationally involving.
With the applied algorithm the maximum distance to the proton track jrj and time t at which pðr; tÞ can be evaluated is limited by the amount of available computational memory.
To tackle this, wave extrapolation algorithms were used as a complement to the propagation algorithm.The algorithm in Sec.II B 2 was used to efficiently calculate the proton acoustic wave for large radial distances.The long time behaviour was approximated by using the temporal wave extrapolation algorithm from Sec. II B 3.
1. Wave propagation: Field for single time instant This scheme relies on the Green's function in k-space for computing a pressure field pðr; tÞ from an arbitrary initial pressure distribution (Cox and Beard, 2005).For the acoustic wave following the proton-induced pressure distribution p 0 ðr 0 Þ, it provides an exact solution that is given by with r 0 a position vector given in Cartesian coordinates.
Using the definition of the Fourier transformation, Eq. ( 7) may be written as  and Millero, 1973;Huang et al., 2009).pðr; tÞ ¼ F À1 x;y;z cos ðv s kðt À t 0 ÞÞF x;y;z p 0 ðrÞ where F x;y;z and F À1 x;y;z denote the three-dimensional spatial Fourier and inverse Fourier transformation, respectively.Numerical computation of the acoustic field at any time t requires only two three-dimensional fast Fourier transformations and a multiplication with the exact time propagator cos ðv s kðt À t 0 ÞÞ.Because the changes of pðr; tÞ over time are calculated using an exact propagator, it is not necessary to calculate the field at intermediate times, as would be required with, for example, finite difference methods.Since the pressure is calculated on a grid of points, the grid spacing must meet the usual Nyquist criterion to avoid aliasing in the spatial domain.

Spatial wavefield extrapolation
A scheme is introduced that may be used to extrapolate the pressure field pðr p ; z; tÞ recorded along a line parallel to the proton track at an arbitrary radial distance r p , to larger radial distances r > r p .It exploits the cylindrical symmetry of the proton acoustic problem, thereby reducing the number of spatial dimensions from 3 to 2. This makes the extrapolation algorithm suitable for swift computation of the pressure field for micrometer radial distances.The algorithm is based on the equation (Candel and Chassaignon, 1984) pðr; z; tÞ ¼ F À1 0 denoting the zeroth order Hankel function of the first kind and l defined as For lr greater than about 2 or in other terms if the distance to wavelength ratio r=k is greater than 1=p, using an asymptotic expansion of the Hankel function is justified to accelerate numerical computations.The wavelength is inherent to the proton acoustic signal but by choosing r sufficiently large the condition can be satisfied For the scheme to yield accurate results, the grid spacing in the axial direction should be sufficiently small compared to the typical wavelength k of the acoustic signal recorded at r p .Typically, Dz should be less than k=2, and much smaller than k=2 for increased accuracy.A smart interpolation method can be used to adjust the sampling frequency of pðr p ; z; tÞ (Verweij and Huijssen, 2009).

Temporal wavefield extrapolation
A scheme is presented that may be used to extrapolate the solution pðr; z; t p Þ at a sufficiently large time instant t p to later times t > t p .In doing so, the cylindrical symmetry of the proton heat deposition and the temporal characteristics of the heat deposition will be exploited.The algorithm enables obtaining the energy content of the proton acoustic wave in the MHz bands, which is desirable because the resonance frequency of ultrasound contrast agents is typically in the MHz range.Given the symmetry of the proton acoustic problem, the thermoacoustic equation from Eq. ( 2) may be written in cylindrical coordinates as 1 r @ @r r @pðr; z; tÞ @r !þ @ 2 pðr; z; tÞ Tðr; z; tÞ @t 2 : (13) Figure 1 shows the local temperature increase following proton impact as functions of time and radial distance to the proton track, for a location in the axial plane of the Bragg peak (z ¼ -1.6 lm).It is found that proton heat deposition occurs within a time t 0 ¼ 70 fs, which is several orders of magnitudes faster than typical acoustic time scales.Therefore, it is justifiable to approximate the temporal behaviour of the temperature distribution by a step function, such that Tðr; z; tÞ % Tðr; z; t 0 ÞHðt À t 0 Þ; (14) with Hðt À t 0 Þ representing the Heaviside step function.
Next, combining Eqs. ( 13) and ( 14) yields 1 r @ @r r @pðr; z; tÞ @r !þ @ 2 pðr; z; tÞ Now, we assume that the pressure at an arbitrary location behaves as if it were generated by a uniform line source.In that case, the solution for the acoustic pressure field in Eq. (15) will behave as the time derivative of the cylindrical Green's function (Verweij et al., 2014).The cylindrical Green's function is given by Taking the time derivative yields @Gðr; tÞ @t ¼ dðt À r=cÞ The asymptotic behaviour of the acoustic pressure field for long times is thus pðr; z; tÞ $ À 1 t 2 ; for t !1: (18) The corresponding asymptotic behaviour for frequencies approaching zero can be found by obtaining the frequency domain counterpart of the time derivative of the Green's function, and consider its behavior for x !0, which is similar to making a low frequency approximation.The Fourier transform of Eq. ( 17) is where J 0 and Y 0 represent the zeroth order Bessel functions of the first and second kind, respectively.As the frequency approaches zero, the Bessel functions can be approximated as (Abramowitz, 1974) By plugging Eqs. ( 20) and ( 21) into Eq.( 19) and taking the modulus, it is found that the low frequency behaviour of the pressure signal is characterised by In deriving Eqs. and ( 22), the dependence of the initial temperature distribution Tðr; z; t 0 Þ on the axial position z has been implicitly neglected.In Fig. 2, this initial temperature distribution Tðr; z; t 0 Þ for t 0 ¼ 70 fs is displayed, which shows that in fact there is a dependence of initial temperature on position.It can be seen that the temperature distribution has a dependency on the axial position, with a characteristic length scale on the order of a micrometer.For positions close to the proton track, say for radial distances on the order of a nanometer, the distribution acts as a uniform line source.The validity of the assumption weakens for larger radial distances, however, it may be expected that the derived asymptotic behaviour at least yields information about the order of magnitude of the low frequency energy content of the proton acoustic wave.

III. RESULTS AND DISCUSSION
A. Proton heat deposition The inelastic thermal spike model from Eq. ( 1) was solved using the Partial Differential Equation Toolbox TM from MATLAB V R , yielding the cylindrically symmetric temperature increase distribution Tðr; tÞ.First, the transient characteristics of Tðr; tÞ are of interest for determining the duration t 0 of the proton heating pulse.Figure 1 shows the temperature increase as functions of time and radial distance to the proton track for a location in the axial plane of the Bragg peak (z ¼ -1.6 lm).It is found that proton heat deposition occurs within a time t 0 ¼ 70 fs.From an acoustical perspective, the spatial temperature distribution Tðr; t 0 Þ is of interest as it sets the acoustical initial conditions.A snapshot of the spatial temperature distribution at t 0 ¼ 70 fs after proton impact is displayed in Fig. 2.

B. Proton acoustic wave simulation
It can be derived from the temperature distribution in Fig. 2 that the characteristic radial dimension of the heat heterogeneity d c is less than 10 nm.Considering the derived heating pulse duration of t 0 ¼ 70 fs and a speed of sound in water v s ¼ 1481 m=s, the condition t 0 < d c =v s for acoustic confinement is well-satisfied for the case of single proton heating.Thus, Eq. ( 4) may be used to obtain the proton-induced acoustic pressure distribution.Since the pressure field in this equation sets the acoustical initial conditions, it will be referred to as the initial pressure distribution, see Fig. 3.

Wave propagation: Field for single time instant
In Fig. 4, snapshots of the acoustic field are shown for three time instants.are computed with a scheme based on Eq. ( 8).Note that the axial dimension is three orders of magnitude larger than the lateral dimensions, such that relative axial wave propagation is marginal over the simulated timescale.
The temporal behavior of the proton acoustic signal, recorded on a line at a radial distance of 30 nm from the proton track, is shown in Fig. 5.The proton acoustic signal in the axial plane of the Bragg peak consists of a bipolar spike with a center frequency of 86.7 GHz.At 30 nm radial distance, the pressure in the Bragg peak has dropped from 45.2 to 3.1 MPa.The dropping of the peak pressure is attributed to the geometrical spreading of the acoustic field, considering that attenuation was not incorporated.The shape of the proton acoustic signal is similar to the signature of the source term in Eq. ( 2).This signature is the double time derivative of the almost stepwise increase in the temperature in Fig. 1, yielding the bipolar behavior of the pressure pulse.Considering the absence of dispersion in the numerical scheme, the shape of the simulated signal necessarily remains unchanged during propagation.The red crosses in Fig. 5(c at the corresponding axial position.Lower frequencies are observed for positions before the Bragg peak.This can be explained by the more widespread initial pressure distribution at these locations, see Fig. 3(a).
To study possible interactions with ultrasound contrast agents, in particular the acoustic energy in the megahertz frequency bands is of interest.As such, it is desirable to simulate the signal for time scales on the order of nanoseconds.Due to memory limitations this was not possible with the current scheme using Eq. ( 8) and therefore a wave extrapolation method was developed.

Spatial wavefield extrapolation
Figure 6 shows the local positive peak pressure for the last five micrometers of the proton track.Peak values up to approximately 5 nm radial distance are set by the initial pressure distribution.Hereafter, the values show a pressure falloff according to an inverse square root relationship with radial distance.For locations close to the proton track, the energy spreading is essentially cylindrical.For cylindrical waves, acoustic energy spreads inversely proportional with radial distance.Since acoustic energy scales with the square of pressure, an inverse square root proportionality for pressure could intuitively be expected.
The spatial wave extrapolation scheme using Eq. ( 12) was used to simulate the proton acoustic wave for radial distances from 150 nm to 5 lm. Figure 7 shows the local positive peak pressures for these radial distances, for the last five micrometers of the proton track.The positive peak pressure decreases further, according to the inverse square root relationship in Eq. ( 12), from 1.41 to 0.24 MPa between 150 nm and 5 lm radial distance.Note that in Fig. 7 at a radial distance of 5 lm a discrepancy, though small, arises with the inverse square root relationship.This is presumably because of the nonuniform initial pressure distribution in the axial direction, which makes that for very large radial distances the proton may be seen as a point source rather than a cylindrical source.

Temporal wave extrapolation
The temporal wave extrapolation scheme was used to approximate the proton acoustic wave for long times and to approximate its frequency spectrum in the low frequency regime.In Fig. 8 the asymptotic behaviour from Eqs. ( 18) and ( 22) has been plotted and compared to the simulated data from Fig. 5.These analytical expressions were plotted by fitting the relation to a single simulation data point.The agreement is good for short times and high frequencies, and therefore the analytical expressions may be used to estimate the pressure for long times and for low frequencies.In Fig. 8(a) it is shown that the time signal runs until 160 ps and is cut off afterwards.The simulated frequency spectrum in Fig. 8(b) was calculated from the numerical time signal using a finite number of samples, which unavoidably affected the low frequency band.The analytical approach taken to derive the asymptotic behaviour in Fig. 8(b) does not suffer from such a procedure and is therefore considered to be a better reflection of reality.The discrepancy between the simulation results and the asymptotic behaviour in Fig. 8(b) arose from the finite length of the numerical time signal.According to the inverse quadratic behaviour with time, it is found that the simulated negative pressure of À2.6 kPa after 160 ps at a radial distance of 30 nm from the Bragg peak will decay to À68 Pa after 1 ns.

IV. CONCLUSIONS
A theoretical model for describing the acoustic field caused by a single proton has been developed, laying a foundation for characterisation of possible acoustic interactions between protons and ultrasound contrast agents or nanodroplets.The model consists of a set of equations describing the heat deposition of a proton and schemes to solve the thermoacoustic wave equation.The numerical simulations yielded the shape, frequency content, and amplitude of the proton acoustic wave on nanometer and micrometer scales.It was found that the proton acoustic wave consists of a bipolar spike with a center frequency of around 86.7 GHz.Positive peak pressures of 45.2 MPa were found in the single proton Bragg peak, falling off according to an inverse square root relationship with radial distance to 0.24 MPa at 5 lm radial distance from the Bragg peak.The temporal behaviour was simulated up to 160 ps and showed that at 30 nm from the Bragg peak a negative pressure of À2.6 kPa was present at that time.An analytical approximation illustrated that the pressure falls off according to an inverse square relationship with time.
All in all, this work has provided a full spatiotemporal characterisation of the acoustic field of a single proton.As next steps, it would be interesting to use this model for predicting how ultrasound contrast agents and nanodroplets would respond to the combined acoustical field generated in a proton beam, where one has an abundant amount of protons, each generating an acoustical field.Comparing these results to recent observations of ultrasound contrast agents in a proton beam could yield valuable insight in the physical processes happening at nanometer and micrometer scales (Lascaud et al., 2021).However, a theoretical gap currently still exists in describing the microbubble response to submicrometer length acoustic pulses, which should first be filled to enable further characterisation of the interaction between protons and ultrasound contrast agents or nanodroplets.

FIG. 1
FIG. 1. (Color online)Temperature increase as a function of time and radial distance, in the axial plane of the Bragg peak (z ¼ -1.6 lm).The temperature peaks in around 70 fs after proton impact, due to the swift energy deposition into the electronic subsystem and subsequent electron-phonon interactions.Since heat diffusion in the atomic subsystem occurs on considerably longer time scales, the temperature remains approximately constant over a picosecond until finally relaxing to the ambient value again.
FIG. 2. (Color online) Temperature distribution at t 0 ¼ 70 fs after proton impact.The Bragg peak spans a region of approximately 3 lm axially and 5 nm radially.Highly localised temperature increases approaching 120 K are observed.
FIG. 3. (Color online) Initial pressure distribution induced by proton impact.Due to the applicability of acoustic confinement, this pressure distribution is essentially a scaled version of the temperature distribution at 70 fs after proton impact.A peak pressure of 45.2 MPa was found in the Bragg peak, of which the axial location is denoted by the red dashed line.
FIG. 4. (Color online) Snapshots of the propagating proton acoustic wave.The images are made with 3.4 ps time intervals, equivalent to a snapshot for every 5 nm of propagation distance.Note the cylindrical symmetry of the propagating wave and the dropping of the peak pressure due to geometric spreading.
FIG. 5. (Color online) Temporal and spectral behaviour of the proton acoustic signal, at a radial distance of 30 nm from the proton track.(a) Time signal.(b) Time signal in the axial plane of the Bragg peak.Although hardly visible, at the position of the red cross a negative pressure of approximately 25 kPa is present.(c) Frequency spectrum of the time recording, normalised against the maximum of the spectrum.The red crosses indicate the center frequency of the signal recorded at the corresponding axial detector position.(d) Frequency spectrum in the axial plane of the Bragg peak.Note that the low frequency part of the spectrum is highly influenced by the truncation of the time domain signal, and therefore is not reliable.
FIG. 6. (Color online) Local positive peak pressures for radial distances up to 150 nm.Peak values up to approximately 5 nm radial distance are set by the initial pressure distribution.Hereafter, the values show a pressure falloff according to an inverse square root relationship with radial distance.FIG. 7. (Color online) Local positive peak pressures for radial distances up to 5 lm.This figure is complementary to Fig. 6.