Designing optimal experiments: an application to proton Compton scattering

Abstract

Interpreting measurements requires a physical theory, but the theory’s accuracy may vary across the experimental domain. To optimize experimental design, and so to ensure that the substantial resources necessary for modern experiments are focused on acquiring the most valuable data, both the theory uncertainty and the expected pattern of experimental errors must be considered. We develop a Bayesian approach to this problem, and apply it to the example of proton Compton scattering. Chiral Effective Field Theory ($\chi$EFT) predicts the functional form of the scattering amplitude for this reaction, so that the electromagnetic polarizabilities of the nucleon can be inferred from data. With increasing photon energy, both experimental rates and sensitivities to polarizabilities increase, but the accuracy of $\chi$EFT decreases. Our physics-based model of $\chi$EFT truncation errors is combined with present knowledge of the polarizabilities and reasonable assumptions about experimental capabilities at HI$\gamma$S and MAMI to assess the information gain from measuring specific observables at specific kinematics, i.e.  to determine the relative amount by which new data are apt to shrink uncertainties. The strongest gains would likely come from new data on the spin observables $\Sigma _{2x}$ and $\Sigma _{2x^\prime }$ at $\omega \simeq 140$ to 200 MeV and $40^\circ$ to $120^\circ$. These would tightly constrain $\gamma _{E1E1}-\gamma _{E1M2}$. New data on the differential cross section between 100 and 200 MeV and over a wide angle range will substantially improve constraints on $\alpha _{E1}-\beta _{M1}$, $\gamma _\pi$ and $\gamma _{M1M1}-\gamma _{M1E2}$. Good signals also exist around 160 MeV for $\Sigma _3$ and $\Sigma _{2z^\prime }$. Such data will be pivotal in the continuing quest to pin down the scalar polarizabilities and refine understanding of the spin polarizabilities.

Appendices

Appendix A: Experimental design details

Suppose that our theoretical model $y_k({x};\vec {a})$ is related to measurements $y_{exp }({x})$ via additive theoretical and experimental noise, as in Eq. (5). We can linearize $y_k({x};\vec {a})$ about some point $\vec {a}_\star$ by keeping only the first order terms in its Taylor expansion, i.e.,

$$\begin{aligned} y_k({x};\vec {a})&\approx y_k({x}; \vec {a}_\star ) + \sum _i b_i({x}) [\vec {a}_i - \vec {a}_\star ] \nonumber \\&= c({x}; \vec {a}_\star ) + \vec {b}({x}) \cdot \vec {a}\,, \end{aligned}$$

(A1)

where $\vec {b}({x}) \equiv \partial y_k({x};\vec {a})/\partial \vec {a}$ evaluated at $\vec {a}_\star$ are our basis functions and $c({x};\vec {a}_\star ) \equiv y_k({x};\vec {a}_\star ) - \vec {b} \cdot \vec {a}_\star$ is constant with respect to the polarizabilities $\vec {a}$ but depends on the kinematic point ${x}$. Thus, the vector of N measurements $\mathbf {y}$ is related to the polarizabilities via the likelihood

$$\begin{aligned} \mathbf {y}\,|\,\vec {a}\sim {\mathcal {N}}[B\vec {a}+ {\mathbf {c}}, \Sigma ] \end{aligned}$$

(A2)

where $B \equiv \vec {b}({\mathbf {x}})$ is an $N \times 6$ matrix, ${\mathbf {c}} \equiv c({\mathbf {x}}; \vec {a}_\star )$ is a length N vector, and $\Sigma$ is the $N \times N$ covariance matrix due to theoretical and experimental error. That is, given some experimental covariance $\Sigma _{\mathrm{exp}}$ and a theoretical covariance ${\bar{c}}^2 R_{\delta k}$ from Eqs. (13) and (14), then

$$\begin{aligned} \Sigma = {\bar{c}}^2 R_{\delta k} + \Sigma _{\mathrm{exp}} \, . \end{aligned}$$

(A3)

Note that $R_{\delta k}$ depends on the values of the tuned $\ell _\omega$ and $\ell _\theta$, whose estimates from the order-by-order convergence pattern are given in Table 1.

The linear model of Eq. (A2) is well known in the statistics literature [68, 69], so here we will simply state the relevant results. If a Gaußian prior with mean $\vec {\mu }_0$ and covariance $V_0$ is placed on $\vec {a}$ as in Eq. (20), then the resulting posterior is also Gaußian, with mean and covariance given by

$$\begin{aligned} \vec {\mu }&= V\left[ V_0^{-1}\vec {\mu }_0 + B^\intercal \Sigma ^{-1} (\mathbf {y}- {\mathbf {c}})\right] \, , \end{aligned}$$

(A4)

$$\begin{aligned} V&= (V_0^{-1} + B^\intercal \Sigma ^{-1} B)^{-1} \, . \end{aligned}$$

(A5)

Importantly to our study of experimental design, the posterior covariance V depends on the kinematic points ${\mathbf {x}}$ where the experiment is performed, and on the specifics of the observable through $\Sigma$, but not on the exact results of the experiment $\mathbf {y}$.

Given that we choose to maximize the expected information gain in the polarizabilities, then the integrals of Eq. (22) must still be performed. The integral over $\vec {a}$ splits into the difference of two terms: the differential entropy of the prior ${{\,\mathrm{pr}\,}}(\vec {a})$ and of the posterior ${{\,\mathrm{pr}\,}}(\vec {a}\,|\,\mathbf {y}, {\mathbf {d}})$. The differential entropy of a Gaußian ${\mathcal {N}}(\mu , \Sigma )$ is well known to be $\frac{1}{2}\ln {|2\pi e\Sigma |}$. Therefore

$$\begin{aligned} U_{\text {KL}}({\mathbf {d}})&= - \int \ln [{{\,\mathrm{pr}\,}}(\vec {a})] {{\,\mathrm{pr}\,}}(\vec {a}) \,{\mathrm{d}}{\vec {a}} \nonumber \\&\quad + \int \ln [{{\,\mathrm{pr}\,}}(\vec {a}\,|\,\mathbf {y}, {\mathbf {d}})]{{\,\mathrm{pr}\,}}(\vec {a}\,|\,\mathbf {y}, {\mathbf {d}}) \,{\mathrm{d}}{\vec {a}} {{\,\mathrm{pr}\,}}(\mathbf {y}\,|\,{\mathbf {d}})\,{\mathrm{d}}{\mathbf {y}} \nonumber \\&= \frac{1}{2}\ln {|2\pi e V_0|} - \frac{1}{2}\ln {|2\pi e V|} \int {{\,\mathrm{pr}\,}}(\mathbf {y}\,|\,{\mathbf {d}})\,{\mathrm{d}}{\mathbf {y}} \nonumber \\&= \frac{1}{2} \ln \frac{|V_0|}{|V|} \, , \end{aligned}$$

(A6)

where we used the fact that V does not depend on $\mathbf {y}$ and then performed the trivial integration over all possible measurements $\mathbf {y}$.

Appendix B: Observable constraints and EFT truncation model details

Constraints on Compton observables are discussed in detail in Ref. [6]. Some of this is reproduced here, with particular attention paid to nth-order chiral corrections to observables $\Delta y_n$ rather than the value of the observable y itself. The $\Delta y_n$ impact the distribution for the $\chi$EFT uncertainty $\delta y_k$, but because we restrict the “experimentally accessible regime” in this study from small $\omega$, and forwards/backwards angles, these constraints are not as important as they otherwise would be. These constraints are summarized for particular $\omega _{\text {lab}}$ and $\theta _{\text {lab}}$ values in Table 5.

Table 5 The constraints on corrections to observables $\Delta y$ and their derivatives $\Delta y'$ at particular $\theta$ and $\omega$. The LO amplitude as well as all calculated higher orders fulfill them automatically, so these must only be enforced in the GP. The observables marked by a dagger $\dagger$ are zero below the pion-production threshold, but we impose no constraint on them at $\omega =\omega _\pi$, as discussed in the text

All observables that are nonzero below $\omega _\pi$ approach the Thomson limit as $\omega \rightarrow 0$ [1]. Thus, higher-order corrections must vanish there , and approach $\omega = 0$ as at least $\omega ^2$. Therefore, at least the first derivative of all corrections must vanish there as well.

The remaining observables must vanish for $\omega \le \omega _\pi$, but there is no constraint on the the derivative of corrections at $\omega = \omega _\pi$. We have found that the corrections approach 0 very quickly, so that imposing the constraint $\Delta y_n(\omega _\pi ,\theta ) = 0$ for all higher order terms is actually a worse approximation than not imposing the constraint at all; see, e.g., Fig. S.4 in the Supplemental Material. This comes back to the large cusps in the spin-observable $c_n$ found near $\omega _\pi$, discussed in Sect. 3.2.1, which remain an unresolved aspect of this model.

Due to the coordinate singularity at $\theta =0^\circ$ and $180^\circ$, observables or their derivative with respect to $\theta$ must vanish there [6]. But this does not preclude both the value and their derivatives from vanishing there. These constraints can be deduced by symmetry arguments, and are summarized in Table 5.

The hyperparameters ${\bar{c}}^2$ and $\ell _i$, shown in Table 1, are tuned to coefficients $c_n$ at the best known $\vec {a}$ (see Table 3) for $\Lambda _b = 650\, \text {MeV}$. The training data are on a grid with $\theta _{\text {lab}}= \{30^\circ , 50^\circ , 70^\circ , 90^\circ , 110^\circ , 130^\circ \}$ and $\omega _{\text {lab}}= \{200, 225, 250\}\, \text {MeV}$ for observables which are zero below $\omega _\pi$. For observables that are non-zero below $\omega _\pi$, the additional training points $\omega _{\text {lab}}= \{50, 75, 100, 125\}\, \text {MeV}$ are included, and common $\ell _\omega$ and $\ell _\theta$ are used between the two regions. The training region is well outside the kinematic endpoints where additional constraints arise on observables or their derivatives, and excludes the pion-production threshold region.

Because the first nonzero order often behaves differently than the corrections, we do not use it for induction on the $c_n$; that is we only train the hyperparameters on corrections. Hence, we train on $c_2$–$c_4$ for $\,{\mathrm{d}}{\sigma }$ and $\Sigma _3$, but otherwise we train on $c_3$ and $c_4$.

The coefficients for various observable slices are shown in Figs. 12 for the cross section, and Figs. S.1–S.6 of the Supplemental Material for the spin observables. These plots also include uncertainty bands for higher order coefficients, with the symmetry constraints given in Table 5 included. These constraints on both the coefficient functions and their derivatives propagate directly to the truncation error $\delta y_k$ by replacing $r(x,x';\ell _\omega , \ell _\theta )$ in Eq. (14) by its conditional form ${\tilde{r}}(x,x';\ell _\omega , \ell _\theta )$, see Refs. [27, 45]. For example, if the value of $c_n$ is known at the set of points ${\mathbf {x}}$, then one can compute its conditional GP, with covariance kernel given by

$$\begin{aligned} {\tilde{r}}(x, x') = r(x, x') - r(x, {\mathbf {x}}) r({\mathbf {x}}, {\mathbf {x}})^{-1} r({\mathbf {x}}, x') \, . \end{aligned}$$

(B1)

See Refs. [70,71,72] for details about adding derivative observations to GPs. Because the RBF kernel [Eq. (11)] is separable in $\omega$ and $\theta$, these constraints can simply be applied to each one-dimensional kernel separately, and multiplied to yield the total constrained kernel. We employ the gptools python package for easily implementing derivative constraints [73].

Fig. 11

A component of the standard deviation due to $\chi$EFT uncertainty at $\mathrm {N}^{4}\mathrm {LO} ^+$, see Eq. (14). The factor of ${\bar{c}}$ is unique to each observable, and is not included. See Table 1

Fig. 12

Coefficients for the differential cross section $\,{\mathrm{d}}{\sigma }$

For completeness, we also provide the profile for the truncation error standard deviation (up to factors of ${\bar{c}}$, which vary by observable); see Fig. 11. It assumes the form of Q provided in Eq. (8) along with the first omitted $\chi$EFT order given in Eq. (18).

This allows us to return to the discussion of the omitted constraints $\Sigma _i \in [-1, 1]$ on the spin observables in Sect. 4.1. Over the physically interesting kinematic range, the actual value of most spin observables lies in the much more narrow interval $[-0.7,0.7]$; see Fig. 5 in Ref. [6]. So, then the question becomes: are the mean prediction and its theory uncertainty contained in $[-1,1]$ with a high degree of probability? From Eq. (12), one can see that the $1\sigma$ interval for the truncation error $\delta y_k$ is $y_{\mathrm {ref}}{\bar{c}}$ times another factor $Q^{\nu _{\delta k}(\omega )}/\sqrt{1-Q^2(\omega )}$. Here $y_{\mathrm {ref}}=1$ and ${\bar{c}}\lesssim 0.7$ for most spin observables (see Table 1). The third factor is plotted in Fig. 11 and does not exceed $\approx 0.3$ at $\omega \lesssim 230\, \text {MeV}$, where our analysis shows the biggest sensitivities. Therefore, even for the spin observables with large magnitudes, the $1\sigma$ upper range of a GP will only give values about 0.2 larger than the established maximum of 0.7, namely about 0.9 in total. This is close but still below $|\Sigma _i|=1$. Therefore, a majority of our test functions in the GP will not probe, let alone exceed, the strict bounds on those spin observables. Furthermore, if observables and their truncation errors vanish at $\theta = 0^\circ$ or $180^\circ$, this will make the constraint even more trivially satisfied near these regions. We are therefore confident that implementing the constraint $\Sigma _i\in [-1, 1]$ would not impact our results for $\omega \lesssim 220\, \text {MeV}$, and cautiously optimistic that the impact would be small even at higher energies.

Though we likewise do not constrain the cross section to be non-negative, we are confident that within our constrained angle range, corrections are highly unlikely to be large enough for this to be a worry. According to Fig. 4 in Ref. [6], the $\mathrm {N}^{4}\mathrm {LO}$ cross section is small ($<10\,\mathrm {nb/sr}$) in a narrow region at forward angles around $\omega _\pi$. Figure 11 shows that the expansion parameter is small, and Fig. 12 shows that the coefficients $c_i$ are natural-sized. Therefore, the GP corrections are highly unlikely to exceed the size of the predicted cross section and create negative (unphysical) values.