How proton pulse characteristics influence protoacoustic determination of proton-beam range: simulation studies

In thermoacoustics, energy pulses cause pulsed heating, which creates expansion that emits longitudinal pressure waves. In the simplest conceptual case, an excited point source will emit a pressure wave proportional to the first time derivative of the excitation pulse. Thus, a Gaussian excitation pulse results in a bipolar acoustic emission consisting of a positive compression (increased pressure) peak—from the increase in the rate of heating as the excitation pulse rises—followed by a negative rarefaction (decreased pressure) peak—from the decrease in heating rate as the excitation pulse falls. The positive and negative pressure peaks are not due to heating and cooling. Rather, they are due to changes in the rate of heating: as the heating rate increases, there is a local expansion and positive emitted pressure, while during decreased heating contraction occurs and a negative pressure wave is observed. Conceptually, extrapolating beyond the point source, the pressure emitted from a complicated energy deposition is the sum of the individual responses that would be observed from decomposing the spatial deposition into point sources. Interestingly, because the acoustic emission depends on the excitation pulse time derivative, the acoustic wave amplitude will depend both on the deposited energy and the temporal shape of the excitation pulse. For a given excitation pulse shape, the thermoacoustic emission amplitude will be linear in the deposited energy (or number of protons for proton irradiation). But, if the excitation pulse has a small first derivative, a large amount of energy can be deposited without generating much thermoacoustic emission. For example, a long rectangular excitation pulse will deposit a large amount of cumulative energy, but acoustic emission will only be generated at the beginning rise (positive, compression acoustic peak) and the ending fall (negative, rarefaction peak). 1 Gy of deposited dose generates a ~240 μK temperature increase in water, but without knowing the shape of the exciting proton pulse, knowledge of the temperature change (integrated dose) is insufficient to predict the magnitude of the emitted protoacoustic pressure waves.

Ignoring heat diffusion and kinematic viscosity, the wave equation (Lai and Young 1982, Sigrist 1986, Tam 1986, Xu and Wang 2006) that describes the pressure, p, at a time, t, and position, $\vec{r}$, is:

$$\begin{eqnarray}&&{{\nabla}^{2}}p\left(\overrightarrow{\boldsymbol{r}},t\right)-\frac{1}{{{c}^{2}}}\frac{{{\partial}^{2}}}{\partial {{t}^{2}}}p\left(\overrightarrow{\boldsymbol{r}},t\right)=\frac{-\beta}{{{C}_{p}}}\frac{\partial H\left(\overrightarrow{\boldsymbol{r}},t\right)}{\partial t}\end{eqnarray} \tag{ 1 }$$

where c is the speed of sound (m s⁻¹), β is the thermal expansion coefficient (K⁻¹), C_p is the specific heat capacity (J kg⁻¹ K⁻¹) and H is the rate of input energy absorption that is converted to heat (J s⁻¹ m⁻³) (see appendix A for a list of symbols and nomenclature).

A Green's function solution (Morse and Feshbach 1953) to equation (1) is given by:

$$\begin{eqnarray}&&p\left(\overrightarrow{\boldsymbol{r}},t\right)=\frac{\beta}{4\pi {{C}_{p}}}\mathop{\int}^{}\mathop{\int}^{}\mathop{\int}^{}\frac{{{\text{d}}^{3}}{{r}^{\prime}}}{|\overrightarrow{\boldsymbol{r}}-{{\overrightarrow{\boldsymbol{r}}}^{\prime}}|}\frac{\partial}{\partial {{t}^{\prime}}}\left[H\left({{\overrightarrow{\boldsymbol{r}}}^{\prime}},{{t}^{\prime}}\right)\right]{{|}_{{{t}^{\prime}}=t-\frac{|\overrightarrow{\boldsymbol{r}}-{{\overrightarrow{\boldsymbol{r}}}^{\prime}}|}{c}}}\end{eqnarray} \tag{ 2 }$$

where the coordinate system is shown in figure 1. In equation (2), two shifted time variables are introduced: t is the time variable at the detector (at $\overrightarrow{\boldsymbol{r}}$) and t' is the time variable at irradiated point ${{\overrightarrow{\boldsymbol{r}}}^{\prime}}$. Given that the speed of protons and cascading electrons is much greater than the speed of sound, c_radiation » c, and the timescale of conversion from radiation to heat is on the order of

Zoom In Zoom Out Reset image size

10⁻¹² s, the local heating induced by the proton pulse is assumed to be instantaneous. t and t' share a common 0 defined by the arrival time of the proton pulse. The acoustic waves emitted from ${{\overrightarrow{\boldsymbol{r}}}^{\prime}}$ at t' will reach the detector at $t={{t}^{\prime}}+\frac{|\overrightarrow{\boldsymbol{r}}-{{\overrightarrow{\boldsymbol{r}}}^{\prime}}|}{c}$.

The energy deposition rate is separable into a time-dependent function, E(t') [J s⁻¹], and a spatial energy deposition, E(s,z) (m⁻³), where, due to the cylindrical symmetry of the proton pencil beams considered here, s  =  (x²  +  y²)^1/2 is used:

$$\begin{eqnarray}&&H\left(s,z,{{t}^{\prime}}\right)=E\left({{t}^{\prime}}\right)E(s,z)\end{eqnarray} \tag{ 3 }$$

E(s,z) is a relative quantity whose volume integral is normalized to 1; for a total input energy ($\int_{0}^{\infty}E(t)\text{d}t$), the fraction that is deposited within a voxel is given by the integration of E(s,z) over that volume. Because the simulation is computed assuming a homogeneous volume, the cumulative dose deposition, D(s,z) [Gy], is given by integrating over the full proton pulse and dividing by the density: $D(s,z)=\frac{E(s,z)}{\rho}\int_{0}^{\infty}E(t)\text{d}t$. Therefore the spatial dose deposition is proportional to the spatial energy deposition. Because the heating is assumed to be instantaneous, then E(t')  =  E(t) is equal to the proton pulse (in units of J s⁻¹). Proton sources operate by modulating a stream of nanosecond proton 'micro'-pulses, whose periods (often  <10 ns) are much faster than the transit time of sound across dose deposition length scales of interest (1 mm/c  =  700 ns » 10 ns). Therefore, the temporal 'micro'-structure is not relevant to this discussion, and, instead, all discussion of proton pulse widths or characteristics refers to the 'macro'-structure modulating envelope.

Equation (2) can be simplified into two equations (see appendix B) that separate the contributions from the spatial distribution, E(s,z), and the time-dependence of the proton pulse, E(t):

$$\begin{eqnarray}&&{{P}_{\delta}}\left(\overrightarrow{\boldsymbol{r}},t\right)=\frac{c\beta}{4\pi {{C}_{p}}}\int_{0}^{2\pi}\int_{0}^{\pi}{{(tc)}^{2}}\sin {{\theta}^{\prime\prime}}\text{d}{{\theta}^{\prime\prime}}\text{d}{{\phi}^{\prime\prime}}\frac{E(s,z)}{tc}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&p\left(\overrightarrow{\boldsymbol{r}},t\right)=\frac{\partial}{\partial t}\left[E(t)*{{P}_{\delta}}\left(\overrightarrow{\boldsymbol{r}},t\right)\right]=\left[\frac{\partial E(t)}{\partial t}*{{P}_{\delta}}\left(\overrightarrow{\boldsymbol{r}},t\right)\right]=\left[E(t)*\frac{\partial {{P}_{\delta}}\left(\overrightarrow{\boldsymbol{r}},t\right)}{\partial t}\right]\end{eqnarray} \tag{ 5 }$$

P_δ is introduced because it is at the heart of calculating the pressure. Equation (4) is written without simplifying and pulling tc out of the integral to highlight the physical meaning of P_δ—it is proportional to a spherical surface integral of the spatial energy deposition divided by the radius, E(s,z)/R'', where R''  =  tc is the radius of a sphere centered at the detector. The 1/R'' in equation (4) is due to the radial dependence of the pressure amplitude, which drops off with the inverse of the radius (the energy/s is proportional to pressure squared, which, in accordance with the inverse-square law, scales with 1/R''²). The symbol P_δ is chosen to highlight its relationship to the pressure. Working back from equation (5), when a Dirac-delta-function (δ-function) proton pulse is substituted into equation (5), it collapses to ${{p}_{\delta}}\left(\overrightarrow{\boldsymbol{r}},t\right)=\frac{\partial}{\partial t}{{P}_{\delta}}\left(\overrightarrow{\boldsymbol{r}},t\right)$. Thus, P_δ is the indefinite integral (with respect to time) of the detected pressure generated by an impulsive, proton pulse (with infinitely small pulse-width). Based on this description we can build an intuitive understanding of P_δ. When a δ-function proton pulse irradiates a homogeneous volume, it creates a source pressure proportional to the energy deposition, E(s,z). Each irradiated source pressure volume element emits a micro bipolar pressure wave. The pressure measured at the detector (p_δ) is the sum of the micro pressure waves emitted by each volume element. Because of the constant sound speed throughout the medium, from the detector's perspective, the pressure that arrives at a certain time t is related to the pressure waves emitted at a common radius (R''  =  tc). Thus, P_δ, a spherical surface integral over the deposited energy (or source pressure), correlates with the expected p_δ that arrives at the detector; due to the derivative relationship, peaks in P_δ translate into bipolar pressure waves in p_δ. This discussion reduces to a simple observation: if the spherical surface integral (P_δ) peaks at a radius that does not pass through the maximum Bragg peak volume, then the p_δ pressure wave arrival time will be offset from the expected arrival time (distance to BP divided by c), which translates into TOF range error.

The above discussion was focused on the pressure waves that result from δ-function proton pulses. Equation (5) is general to any proton pulse. For non-δ-function proton pulses, the pressure measured at the detector is equal to the first time derivative of the convolution of the proton pulse (E(t)) with P_δ.

As the first step in simulating the protoacoustic pressure waves, to determine the manifestations of irradiating with different proton pulses, equation (4) was used to discretely, numerically calculate the P_δ in water (the double integral was summed over θ" and φ'' with steps of dθ'' and dφ''  ⩽  0.01 rad) for a grid of detector positions: s  =  0–100 mm in 5 mm steps and z  =  4–309, 334–709 mm in 5 mm and 25 mm steps, respectively. Given the cylindrical symmetry of the dose deposition, this 2D planar grid of detectors is equivalent to a 3D array of detectors covering a cylinder of radius s  =  100 mm and length z  =  709 mm. An ideal, point detector was assumed. P_δ was calculated with Δt  =  100 ns time steps (Δt*c  =  148.1 μm) with c  =  1481 m s⁻¹, which was chosen to match the sound speed in water at 20 °C. For E(t), different Gaussian proton pulses with full-width half-max (FWHM) of $\sigma _{t}^{\text{FWHM}}$ and a maximum instantaneous amplitude of 300 nA proton current (1.87  ×  10¹² protons s⁻¹, 44.9 J s⁻¹) were then substituted into equation (5) to calculate the simulated pressure, p. Single Gaussian pulses of $\sigma _{t}^{\text{FWHM}}$  =  1.4 and 56 μs deposit 1.1 and 43.3 cGy (2.8  ×  10⁶ and 1.1  ×  10⁸ protons) per pulse. Calculation of each P_δ took ~100 s on a personal computer. Thus, the pressure simulations were performed by first calculating P_δ (equation (4)) and then substituting it into equation (5) with the specified proton pulse, E(t).

The 300 nA instantaneous proton current limit was chosen to assess the pulse length effect on protoacoustic amplitude because we anticipate that proton sources will be limited by maximum current, either due to safety restrictions (which are really intended to limit average proton current rather than instantaneous current (Jones et al 2015)) or instrument capacity (there is an upper limit to the proton current output). The simulations are therefore meant to reflect the situation in which the operator has control over the pulse shape but is limited to an instantaneous proton current ceiling. Another logical approach is to simulate various pulses with equivalent protons per pulse. Conversion to this framework is straightforward: for a particular proton pulse, the pressure amplitude is proportional to the number of protons per pulse (which is also proportional to $\sigma _{t}^{\text{FWHM}}$: the dose deposited at the BP by a single 300 nA pulse  =  0.77 cGy/μs *$\sigma _{t}^{\text{FWHM}}$). This analysis is also included in the plot of the amplitude as a function of dose (see figure 4 below).

For clinical relevance, the spatial energy deposition for a 150 MeV proton pencil beam was used. The energy deposition, E(s,z), was calculated analytically based on empirical measurements using

$$\begin{eqnarray}&&E(s,z)=\frac{1}{\int_{-\infty}^{\infty}\text{IDD}(z)\text{d}z}\frac{\text{IDD}(z)}{\left(\frac{\pi \sigma _{s(z)}^{\text{FWHM}}}{4\ln 2}\right)}\exp \left(\frac{-{{s}^{2}}4\ln 2}{\sigma _{s(z)}^{\text{FWHM}}}\right)\end{eqnarray} \tag{ 6 }$$

This equation separates the dose deposition into the product of the lateral profile and the depth profile. Along the depth, the characteristic proton Bragg peak is described by IDD(z), the integrated dose deposition of a pristine Bragg peak with maximum dose deposited at 14.9 cm (Bortfeld and Schlegel 1996, Bortfeld 1997). The lateral profile is given by a cylindrically-symmetric Gaussian profile along s. The investigated proton beam's diameter increases with penetration (until the BP depth) due to scattering, which is captured by $\sigma _{s(z)}^{\text{FWHM}}$, the depth dependent lateral Gaussian FWHM. IDD(z) and $\sigma _{s(z)}^{\text{FWHM}}$ are both drawn from experimental measurements(Lin et al 2014), and they are shown in figure 2. The $\sigma _{s(z)}^{\text{FWHM}}$ is beam-line specific. We use the parameters measured at the Roberts Proton Therapy Center at the University of Pennsylvania(Lin et al 2014). The denominator terms in equation (6) are included for normalization so that $\iiint E(s,z)\text{d}x\text{d}y\text{d}z=1$. The proton energy distribution is shown in the xz plane in figure 3(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To summarize the above, using the energy distribution from figure 3(a), P_δ(t) was calculated for various detectors (at positions $\overrightarrow{\boldsymbol{r}}$). Then, from P_δ(t), p(t) was calculated for each detector as a function of variable proton pulses ($\sigma _{t}^{\text{FWHM}}$). From these data traces, different measures were collected to assess the amplitude, frequency, and range verification accuracy. Table 1 lists the data and features measured for analysis. To understand the amplitude dependence, the maximum pressure, p_max, was measured from each p(t). To assess the frequency dependence, the time-domain pressure traces were Fourier Transformed, and the spectral peak frequency, ω_max, and width, $\sigma _{\omega}^{\text{FWHM}}$, were measured. To understand the range verification accuracy, the arrival time of different features was measured as described below.

Table 1. The simulation data and extracted information for each study. Simulation data was calculated at various detector positions ($\overrightarrow{\boldsymbol{r}}$) and for Gaussian excitation proton pulses of various widths ($\sigma _{t}^{\text{FWHM}}$)

Study	Simulation data	Extracted information
Amplitude	$p(t)$	pmax
Frequency	$|FT\left[\,p(t)\right]|$	ωmax, $\sigma _{\omega}^{\text{FWHM}}$
Range verification accuracy	Pδ(t)	$\tau \left(P_{\delta}^{\max}\right)$
	$p(t)$	τcenter, τmax, τmin, τzero-cross

The TOF calculation of distance is based on the arrival time, τ, of the pressure wave of interest. In TOF calculations, the distance is calculated as the product of the speed of sound and τ. Therefore, the error in determining the distance to the Bragg peak from a TOF calculation, ε^TOF, is given by:

$$\begin{eqnarray}&&{{\varepsilon}^{\text{TOF}}}\left(\overrightarrow{\boldsymbol{r}},\sigma _{t}^{\text{FWHM}}\right)=\tau _{{}}^{\gamma}c-l_{\text{BP}}^{\det}\end{eqnarray} \tag{ 7 }$$

where ε is expected to depend on the detector position, $\overrightarrow{\boldsymbol{r}}$, and the proton pulse width, $\sigma _{t}^{\text{FWHM}}$. To calculate the distance to the beam propagation axis, $s_{\text{beam}}^{\det}$, then τ^α is used in place of τ^γ. Because simulation data is noiseless, the error (ε) is a standard error. Initial analysis of the results plots ε as a function of detector position, $\overrightarrow{\boldsymbol{r}}$, for a particular $\sigma _{t}^{\text{FWHM}}$. Later analysis attempts to map a more complicated exponential relationship between τ^γ and $l_{\text{BP}}^{\det}$ for particular $\sigma _{t}^{\text{FWHM}}$ without considering $\overrightarrow{\boldsymbol{r}}$.

Because each protoacoustic pressure wave consists of a bipolar combination of a positive compression peak followed by a negative rarefaction peak, there are multiple metrics to determine τ. Considered arrival times metrics are indicated in figures 3(d) and (e). τ_max and τ_min are the arrival times of the compression and rarefaction peaks, respectively. τ_zero-cross is the time when the pressure wave crosses zero between τ_max and τ_min. τ_center is the average of τ_max and τ_min.

The simulation calculates a noiseless dataset. In experiments, however, averaging is essential to remove random electrical, acoustic, and thermal noise because of the low levels of protoacoustic signals. To understand the manifestations of noise, two analyses were performed, one to investigate the signal-to-noise ratio (SNR) and another to assess how noise affects the calculation of arrival times. The effects of random noise on SNR are well established, and the noiseless simulation data was analyzed using the analytic equations below. The effect of random noise on determining arrival times is more complicated, however, and measurements of τ were repeated with noise added to the simulation data to empirically assess the manifestations.

2.4.1. SNR.

Using p_max from the simulations, the SNR was assessed for a hypothetical experimental random noise level using the analytical equation (8). If the random noise has a single measurement standard deviation of N, then the SNR resulting from averaging the protoacoustic signal resulting from n averages is given by:

$$\begin{eqnarray}&&\text{SNR}=\sqrt{n}\frac{{{p}_{\max}}}{N}\end{eqnarray} \tag{ 8 }$$

where p_max is the maximum protoacoustic signal amplitude resulting from a single proton pulse measured at the considered detector. For n proton pulses, the total dose at the Bragg peak, $D_{\text{total}}^{\text{BP}}$, is the product of n with the dose at the Bragg peak from a single proton pulse, $D_{\text{total}}^{\text{BP}}=nD_{\text{single}}^{\text{BP}}$. Clinically, the number of protoacoustic traces that can be averaged is limited by the total dose deposited at the Bragg peak. The ratio of SNR to $D_{\text{total}}^{\text{BP}}$ is given by:

$$\begin{eqnarray}&&\frac{\text{SNR}}{D_{\text{total}}^{\text{BP}}}=\left[\left(\frac{{{p}_{\max}}}{D_{\text{single}}^{\text{BP}}}\right)\sqrt{D_{\text{single}}^{\text{BP}}}\right]\frac{1}{N\sqrt{D_{\text{total}}^{\text{BP}}}}\end{eqnarray} \tag{ 9 }$$

From a detection stand point, the ideal proton pulse sequence is the one that maximizes the SNR while delivering the prescribed dose $\left(D_{\text{total}}^{\text{BP}}\right)$. p_max and $D_{\text{single}}^{\text{BP}}$ depend on the Gaussian proton pulse width ($\sigma _{t}^{\text{FWHM}}$) while $D_{\text{total}}^{\text{BP}}$ is a constant specific to the patient's prescribed dose (typically between 1.8 and 2.5 Gy per fraction).

2.4.2. Arrival times in the presence of noise.

The p_max simulated here (t_10–90%  =  18 μs, p_max  =  70 mPa,) is on the same order of magnitude as a recent experimental report with similar beam characteristics (15–20 mPa) (Jones et al 2015). The difference in simulation and experiment may be due in part to the finite bandwidth of the detector (Ahmad et al 2015), although further experiments are necessary for further comparison.

The proton pulse with a rise time of 2.9 μs (4.0 μs $\sigma _{t}^{\text{FWHM}}$) generates the maximum amplitude per pulse for the beam characteristics and distal detector position considered here (figure 4(a)). At shorter Gaussian $\sigma _{t}^{\text{FWHM}}$, the compression and rarefaction features destructively interfere, which causes the drop in maximum pressure at $\sigma _{t}^{\text{FWHM}}$  <  4.0 μs. For the ~600 μs FWHM rectangular proton pulses, there is no destructive interference because of the time separation between the compression and rarefaction waves (generated by the rise and fall of the proton pulse), and the drop in maximum pressure with decreasing rise time is not observed. Agreement between the maximum amplitude for the rectangular and Gaussian pulses at long rise times indicates that the amplitude is directly related to the rise time; despite delivering different numbers of protons, rectangular and Gaussian pulses with the same rise times generate protoacoustic emissions with the same amplitudes. For example, for a proton pulse with t_10–90%  =  6.5 μs, a ~600 μs rectangular pulse (1.1  ×  10⁹ protons, 4.29 Gy at BP) and a Gaussian proton pulse (1.8  ×  10⁷ protons, 7.0 cGy at BP) generate a protoacoustic pulse with the same p_max (112 and 117 mPa, respectively). Moreover, ~600 and ~800 μs rectangular proton pulses with the same rise time generate the same p_max. Mathematically, based on the convolution shown in equation (5), the maximum protoacoustic amplitude per pulse is achieved when the Gaussian proton pulse best overlaps the compression peak of p_δ. Intuitively, a δ-function-like proton pulse initiates the emission of the p_δ pressure waves from the irradiated deposition volume. If the proton pulse is longer than the p_δ pressure wave peak-widths, then stress confinement is not satisfied, and the heating (induced by the proton pulse) extends the pressure emission time period. If the proton pulse is shorter than the peak-widths of the p_δ pressure wave, then less energy is deposited and the amplitude is lower than ideal. If the proton pulse matches the peaks of the p_δ pressure wave, then the protoacoustic emission amplitude is maximized without broadening the pressure signal. Plotting p_max as a function of $\sigma _{t}^{\text{FWHM}}$ for different detector positions (not shown) reveals the same trend as in figure 4(a), but p_max peaks at different $\sigma _{t}^{\text{FWHM}}$ (for figure 3(a) detectors #1 and #2, p_max peaks at $\sigma _{t}^{\text{FWHM}}$  =  8 and 7 μs, respectively). This is because the p_δ pressure waves have different widths depending on the detector's orientation relative to the proton beam, which will match different width proton pulses.

Although the amplitude per pulse (p_max) has a non-monotonic dependence on rise time, the ratio of protoacoustic amplitude per deposited dose (${{p}_{\max}}/D_{\text{single}}^{\text{BP}}$) monotonically increases as the Gaussian proton pulse shortens (equivalently, as the rise time becomes faster) (figure 4(b)). Because both p_max and $D_{\text{single}}^{\text{BP}}$ are proportional to the maximum instantaneous proton current (300 nA here), their ratio, ${{p}_{\max}}/D_{\text{single}}^{\text{BP}}$, is independent of the proton current magnitude (and the number of protons per pulse). Fast rise time proton pulses, however, may not lead to the highest signal-to-noise ratios per total deposited dose. In previous studies, averaging of the protoacoustic signal has been necessary to remove random noise (Sulak et al 1979, Hayakawa et al 1988, 1989, Tada et al 1991, Hayakawa et al 1995, Albul et al 2001, Albul et al 2005, 2004, Capone and De Bonis 2006, Graf et al 2006, Terunuma et al 2007, De Bonis 2009, Bychkov et al 2010, Assmann et al 2015). Based on the analytical equation (8), the SNR depends on the square root of the number of averages while the signal per Gaussian proton pulse peaks at $\sigma _{t}^{\text{FWHM}}$  =  4.0 μs. Once both of these effects are considered (equation (9)), the SNR per total dose peaks for a Gaussian proton pulse with $\sigma _{t}^{\text{FWHM}}$ of 2 μs (figure 4(b)). Clinically, if the signal from many proton pulses is averaged, this study indicates that a 2 μs Gaussian proton pulse will provide the highest SNR per total dose relative to other Gaussian pulses (for this detector position). Because the $\sqrt{D_{\text{single}}^{\text{BP}}}$ term depends on the proton current magnitude, the $\text{SNR}/D_{\text{total}}^{\text{BP}}$ will also exhibit a square-root proton current dependence (higher proton currents translate into higher $\text{SNR}/D_{\text{total}}^{\text{BP}}$), as is shown in figure 4(c). For other beam characteristics, detector positions, and non-Gaussian proton pulse shapes the $\text{SNR}/D_{\text{total}}^{\text{BP}}$ may peak at different values. For detectors #1 and #2 from figure 3(a), $\text{SNR}/D_{\text{total}}^{\text{BP}}$ peaks at $\sigma _{t}^{\text{FWHM}}$  =  4 μs.

Based on the pressure trace simulated at (s,z)  =  (0,709 mm), approximately ~10⁻¹¹ of the total deposited proton beam energy is converted into protoacoustic pressure waves. This protoacoustic energy is calculated by assuming the pressure is a spherical wave originating at the BP through $4\pi {{r}^{2}}\mathop{\int}^{}\text{d}t\left({}^{p{{(t)}^{2}}}/{}_{2\rho c}\right)$. The ratio of protoacoustic pressure energy to deposited proton energy is largest at $\sigma _{t}^{\text{FWHM}}$ ~ 2 μs. This is consistent with the above discussion: when the proton pulse matches the p_δ(t) acoustic peaks, energy will be most efficiently transferred

For a mono-energetic proton pencil beam, the dose deposition can be approximated by a pre-Bragg peak cylinder and a Bragg peak disc(Jones et al 2014). The cylinder segment emits pressure waves that exhibit a s^−0.5 amplitude fall off in the radial direction, which is expected for a cylindrical antenna(Sulak et al 1979)(figure 5(c)). Along the beam axis, distal to the Bragg peak, the disc segment emits pressure waves that exhibit an $l_{\text{BP}}^{\det}$^−0.96 amplitude fall off, which is close to the r⁻¹ behavior expected for a point source(figure 5(c)). Because the pressure waves are a linear combination of the contributions from each separate segment, the amplitude falloff for the full dose deposition will exhibit a more complicated behavior.

Characterizing the frequency spectrum through simulation allows for proper selection of detectors. In particular, cutting out high frequency components either through inappropriate filtering or the use of bandwidth-limited detectors will broaden time-domain protoacoustic peaks, which may skew the measurement of arrival times (and BP positions). Exclusion of low frequency components is expected to be less drastic. Understanding the spectrum also allows for future assessment of protoacoustic wave propagation through tissue, which is frequency-dependent.

As the proton pulse rise time lengthens, the average protoacoustic spectrum decreases in frequency. The protoacoustic frequency spectra for proton pulses with $\sigma _{t}^{\text{FWHM}}$  ⩽  1.4 μs (⩽1.0 μs rise times) are indistinguishable, which indicates that the these shorter pulses are δ-function-like. For proton pulses with $\sigma _{t}^{\text{FWHM}}$  ⩾  5.6 μs (⩾4.0 μs rise times), the longer rise times are blurring out sharper features in P_δ. Even for the δ-function-like proton pulses, there is negligible amplitude above 500 kHz, and the frequency spectrum peaks at 45.5 kHz.

The frequency, ω, of the protoacoustic spectrum is expected to be defined by the lesser of either the proton pulse FWHM (ω ~ 1/$\sigma _{t}^{\text{FWHM}}$) or the transit time of sound across the dose deposition: ω  ≈  1/ ($\sigma _{s(z)}^{\rm FWHM}/c$) for detectors placed in the pre-Bragg peak region or ω  ≈  1/ ($\sigma _{\text{BP}}^{\text{FWHM}}/c$). From figure 6(b), the frequency spectrum from the protoacoustic signal simulated at a detector placed at (s  =  0, z  =  309 mm) peaks at 100–200 kHz for a δ-function-like proton pulse, which is consistent with the prediction (ω($\sigma _{\text{BP}}^{\text{FWHM}}$  =  7.8 mm)  ≈  190 kHz). As the proton pulse $\sigma _{t}^{\text{FWHM}}$ increases at and above 5.6 μs, the ω_max from the single detector decreases, which is also consistent with prediction (ω_max($\sigma _{t}^{\text{FWHM}}$  =  5.6 μs)  ≈  179 kHz); once the proton pulse $\sigma _{t}^{\text{FWHM}}$ increases beyond 1.4 μs, the maximum of the frequency spectrum becomes limited by the proton pulse time characteristics rather than the transit time of sound across the dose deposition structures. In photoacoustic measurements, short photon pulses allow for high spatial resolution of heterogeneous structures within the illuminated volume. In protoacoustic measurements, the current goal is to determine the characteristics of the dose deposition, which has millimeter-length-scale features. δ-function-like proton pulses will not limit the frequency spectrum, but using short proton pulses with $\sigma _{t}^{\text{FWHM}}$  <  1.4 μs (rise time  <1.0 μs) will provide no advantage in dose deposition spatial resolution. If, on the other hand, protoacoustics was applied similarly to photoacoustics to image heterogeneous structures in the illuminated volume, then faster rising proton pulses would be required to achieve sub-millimeter resolution.

Given that there is negligible spectral amplitude at  >500 kHz in both the summed spectra and the individual spectra (figures 6(a) and (b)), low frequency (<500 kHz) transducers are ideal for protoacoustic measurements of the simulated proton beam. If slower rising proton pulses are used, then lower frequency transducers are needed (or are sufficient). If lower energy protons are used, the spatial dimensions of the dose deposition $\sigma _{\text{BP}}^{\text{FWHM}}$ will be reduced, and proton pulses with faster rise times will be required along with transducers having higher matching bandwidths (Assmann et al 2015).

4.3.1. δ-function-like proton pulses.

4.3.2. Non δ-function-like proton pulses.

4.3.3. Simulation-based mapping.

4.4.1. Comparison of different metrics.

4.4.2. Noise.

4.4.3. Deconvolutions.

Although the analyzed simulations were restricted to Gaussian proton pulses and a 150 MeV beam with a beam-line specific diameter, the results can be generalized more broadly. Although the proton pulse was primarily identified by its FWHM, $\sigma _{t}^{\text{FWHM}}$, the rise time (τ_90−10%  =  0.716$\sigma _{t}^{\text{FWHM}}$) is the important quantity that affects the protoacoustic signal properties because the acoustic waves are related to the time derivative of the proton pulse, and the rise time is a measure of the first derivative. The proton pulses should be sorted into two categories: δ-function-like and longer. The δ-function-like proton pulses have rise times less than the transit time of sound across the dose deposition features. This restriction is known in thermoacoustics as stress confinement (Xu and Wang 2006). Within the δ-function-like range, the protoacoustic frequency spectrum and TOF-calculated range accuracy are independent of proton pulse rise time. In this case, the accuracy of calculating distance to the BP from the TOF is determined by the inherent systematic error. For the considered 150 MeV beam, the cross-over from δ-function-like to longer proton pulses occurs between $\sigma _{t}^{\text{FWHM}}$ of 1.4 and 5.6 μs (τ_90-10% of 1.0 and 4.0 μs) and the systematic error limits the accuracy to  ±2.6 mm (figure 7). Because of the strong angular dependence of the error, however, knowledge of the detector relative to the beam axis will provide corrections that improve that error (along the beam axis and distal to the BP, the inherent systematic error is ~0 mm).

For non δ-function-like proton pulses, increasing the rise time results in a lower frequency protoacoustic spectrum and larger errors in TOF calculations of the BP distance. Interestingly, tabulating the data from simulations allows for more accurate range calculations by fitting an exponential function that maps the BP distance to the arrival time of the protoacoustic wave. For example, for the proton pulse with $\sigma _{t}^{\text{FWHM}}$ of 56 μs, mapping allows for more accurate determination of the BP distance (<6 mm systematic error) compared to simple TOF-calculation (up to  >20 mm). Although mapping requires the complication of interpreting clinical results with matching simulations, other proton range verification methods such as PET and prompt gamma also rely on supporting simulations (from Monte Carlo) to interpret data (Moteabbed et al 2011, Sterpin et al 2015).

Examining the range verification accuracy in light of the low frequencies reveals an interesting observation highlighting the difference between ultrasound and protoacoustics. In ultrasound, due to the reflections inherent to measurement, the best axial resolution is defined by ~c/2ν, where ν is the frequency. In protoacoustic range verification, although lower frequency will generally translate into decreased accuracy, the same relationship between signal frequency and range verification accuracy does not apply: the arrival time of a 50 kHz protoacoustic wave can be determined to within a standard deviation of 0.5 μs (0.7 mm) for SNR  =  10. Stated differently, even though the peaks in the protoacoustic signal are broad and the frequencies are low, the arrival time can be measured accurately.

Aside from TOF accuracy, the other important factor for clinical application of protoacoustic range verification is the SNR per total deposited dose. $\text{SNR}/D_{\text{total}}^{\text{BP}}$ peaks at Gaussian proton pulse widths close to, but shorter than, the turnover point between δ-function-like and longer proton pulses. Although the shape of $\text{SNR}/D_{\text{total}}^{\text{BP}}$ plotted in figure 4(b) will not change with proton current (protons/pulse), the $\sqrt{D_{\text{single}}^{\text{BP}}}$ term of equation (9) imparts a square root proton current dependence on the magnitude of $\text{SNR}/D_{\text{total}}^{\text{BP}}$. Therefore, a maximum protons/pulse should be used to maximize $\text{SNR}/D_{\text{total}}^{\text{BP}}$.

Therefore, to retain the TOF accuracy and maximize the $\text{SNR}/D_{\text{total}}^{\text{BP}}$, the ideal Gaussian proton pulse for generating protoacoustic emissions has a width slightly less than the width of pressure peaks in p_δ. For the dose deposition considered here, for a detector placed on the beam axis and distal to the Bragg peak, the ideal $\sigma _{t}^{\text{FWHM}}$  =  2 μs.

For clinical sources that have longer-than-ideal proton pulses ($\sigma _{t}^{\text{FWHM}}$  >  2 μs), there will be a reduced $\text{SNR}/D_{\text{total}}^{\text{BP}}$ and larger range calculation errors, but simulation based mapping has the potential to reduce the error down to 2 mm in the region beyond the BP. All current clinical proton sources operate with proton pulse rise times that are classified as larger than δ-function like. Synchrocyclotron sources ($\sigma _{t}^{\text{FWHM}}$  =  7 μs(Kleeven et al 2013)) are the most promising for protoacoustic measurements given they can generate proton pulses closest to the ideal pulses identified here. It is expected that their TOF accuracy will display a ~1–2 mm systematic offset, as in the $\sigma _{t}^{\text{FWHM}}$  =  9.8 μs case, and the $\text{SNR}/D_{\text{total}}^{\text{BP}}$ is ~0.5 of the ideal case. Specially modulated ($\sigma _{t}^{\text{FWHM}}$  =  18 μs (Jones et al 2015)) and standard cyclotrons ($\sigma _{t}^{\text{FWHM}}$ ~ 50 μs (Measurements 2007)) will display TOF errors that require mapping (see table 2), their $\text{SNR}/D_{\text{total}}^{\text{BP}}$ will drop to 0.2 and 0.07 of the ideal proton pulse, respectively, and the high dose per pulse (19 and 54 cGy, respectively) limits the number of averages before the prescribed treatment dose is delivered. The low amplitude, high dose per proton pulse (2.2 Gy for maximum proton current of 300 nA), and low expected TOF accuracy for standard synchrotron pulses ($\sigma _{t}^{\text{FWHM}}$ ~ 200 μs (Measurements 2007)) most likely makes them untenable for clinical protoacoustic measurements. Potentially, these identified proton pulse rise times may be reduced through modifications such as the addition of deflectors or choppers (Yokota et al 1997) (isochronous cyclotrons) or pulse-picking mechanisms (Chuvilo et al 1984) (synchrotrons).

Analysis of the frequency spectra reveals that there are negligible frequency components above 500 kHz. Thus, a  <500 kHz low frequency hydrophone is ideal for experimental measurements of clinical proton beams. Any proton pulse temporal structure less than the δ-function-like timescale is blurred out by the spatial dimensions of the dose deposition. Hence, the ns 'micro' structure of clinical proton beams is not apparent in the protoacoustic signal (see section 2.1).

Although clinical application of protoacoustic range-verificaiton will have to consider the manifestations of tissue heterogeneity, the presented simulations lay the ground work for understanding the effects of using different proton pulses. Tissue heterogeneity introduces attenuation, reflections, and variable sound speeds. Attenuation will reduce the amplitude, although the low (<100 kHz) frequency ensures better transmission than (>1 MHz) ultrasound waves, and tissue's higher (37° C) temperature and higher Grüneisen parameter are expected to increase initial protoacoustic amplitudes. The reflections will create complicated echoed pressure waves, but, because it is measured from the earliest arriving peak, τ_max is expected to be unaffected. As in ultrasound, an average sound speed can be used to minimize the effects of variable sound speeds, although this will be a source of additional error.

The presented work considers the protoacoustic pressure wave generated by a single PBS beam with a pristine Bragg peak. During treatment, the irradiated volume is raster-scanned in layers by pristine BP beams with varying energies (depths) and intensities. In theory, given high enough SNR, the protoacoustic signal generated by each individual pristine BP can be measured by synchronizing acoustic detection with delivery. Without synchronization, the detected protoacoustic signal is the sum of the signals generated by each individual beam. Extracting information from the resulting sum is complicated (arguably impossible if the exciting proton pulse shape is variable), and would most likely require 3D arrays to map dose(Alsanea et al 2015). Given the large beam sizes and long beam-on times (10s of ms), protoacoustic detection is less promising for monitoring double-scattered proton therapy beams compared to PBS.