Performing Bayesian analyses with AZURE2 using BRICK: an application to the ${}^7$Be system

Phenomenological $R$-matrix has been a standard framework for the evaluation of resolved resonance cross section data in nuclear physics for many years. It is a powerful method for comparing different types of experimental nuclear data and combining the results of many different experimental measurements in order to gain a better estimation of the true underlying cross sections. Yet a practical challenge has always been the estimation of the uncertainty on both the cross sections at the energies of interest and the fit parameters, which can take the form of standard level parameters. Frequentist ($\chi^2$-based) estimation has been the norm. In this work, a Markov Chain Monte Carlo sampler, \texttt{emcee}, has been implemented for the $R$-matrix code \texttt{AZURE2}, creating the Bayesian $R$-matrix Inference Code Kit (\texttt{BRICK}). Bayesian uncertainty estimation has then been carried out for a simultaneous $R$-matrix fit of the $^3$He$(\alpha,\gamma)^7$Be and $^3$He$(\alpha,\alpha)^3$He reactions in order to gain further insight into the fitting of capture and scattering data. Both data sets constrain the values of the bound state $\alpha$-particle asymptotic normalization coefficients in $^7$Be. The analysis highlights the need for low-energy scattering data with well-documented uncertainty information and shows how misleading results can be obtained in its absence.


I. INTRODUCTION
Phenomenological R-matrix has been the standard analysis tool for cross section data that exhibit overlapping yet resolved resonances for many years [1]. It is used extensively to evaluate data for applications (e.g. the ENDF/B-VIII.0 evaluation [2]), to perform extrapolations to low, unobserved energies in nuclear astrophysics (e.g. Azuma et al. [3], Descouvemont et al. [4]), and to extract level parameters for nuclear structure. In all cases, it provides a reaction framework in which experimental information of various different types can be combined to improve estimates of the true cross sections. One challenging aspect of this type of analysis has been reliable uncertainty propagation.
Traditionally, data have been fitted using χ 2 minimization, with uncertainties being estimated using one of two methods. The first is using partial derivatives and the assumption that the quantity of interest is related linearly with the parameters of the model. The second is the assignment of confidence intervals based on some ∆χ 2 value. The assumption of linearity is often a poor one and the second method can become tedious or impossible to implement for a complicated model. Additional limitations are that one must assume Gaussian uncertainties on the input data and there is almost no ability to include prior information about the parameters. It is known that χ 2 methods may lead to biased results and/or underestimated uncertainties in data evaluations [5]. The reason for these issues is understood to be incomplete documentation or modeling of systematic uncertainties. While systematic uncertainties are a difficult subject in any approach, they are much easier to model and implement using the Bayesian methods described below. Finally, we would like to point out that a mixed approach is possible, where χ 2 minimization is combined with a Monte Carlo simulation of some uncertainties. This method was used by deBoer et al. [6] in a previous analysis of 3 He(α, γ) 7 Be and 3 He(α, α) 3 He.
Bayesian methods are increasingly becoming the standard for performing Uncertainty Quantification in physical sciences and engineering in general, and theoretical nuclear physics in particular . In contrast to a traditional χ 2 -minimization they offer the opportunity to examine the entire probability distribution for parameters of interest, rather than focusing on the values that maximize the likelihood. Perhaps equally important, in a Bayesian approach it is straightforwardmandatory even-to declare and include prior information on the parameters of interest. Bayesian methods, combined with the possibility to use Markov Chain Monte Carlo sampling to explore a high-dimensional parameter space, allow one to introduce additional parameters without fear of computational instabilities caused arXiv:2112.12838v2 [nucl-th] 31 Dec 2021 by shallow χ 2 minima. The use of MCMC sampling also makes uncertainty propagation straightforward, as we will demonstrate here. And a Bayesian framework is-to our knowledge-the only option if one wishes to incorporate a rigorous formulation of theory uncertainties into the statistical analysis. In this work, Bayesian uncertainty quantification is implemented by pairing the Rmatrix code AZURE2 [3,34] with the MCMC Python package emcee [35]. The pairing is facilitated by a Python interface BRICK (Bayesian R-matrix Inference Code Kit), enabling Bayesian inference in the context of R-matrix analyses.
To benchmark this code, it has been applied to the analysis of the 3 He(α, γ) 7 Be and 3 He(α, α) 3 He reactions. The 3 He(α, γ) 7 Be reaction is a key reaction in modeling the neutrino flux coming from our sun [36]. It also plays a role in Big Bang Nucleosynthesis (BBN). The reaction cross section is dominated by the direct capture process, but also has significant contributions from broad resonances (see Fig. 1). In recent years, high-precision measurements of this reaction have been performed, using direct γ-ray detection [37][38][39], the activation method [38][39][40][41][42], and a recoil separator [43]. Additional higher energy measurements have also been made recently by Szücs et al. [44], but are outside the energy range of the present analysis. Using these high precision measurements, several analyses have been made to combine these data sets and extrapolate the cross section to low energies using pure external capture [45], R-matrix [6], effective field theory [21,46], a modified potential model [47], and ab initio calculations [48][49][50][51]. These several recent analyses make this reaction an ideal case for benchmarking since they employ both more traditional and Bayesian uncertainty estimation methods.
As the energies pertinent to solar fusion and BBN the 3 He(α, γ) 7 Be cross section has a large contribution from external capture, 3 He(α, α) 3 He data, through its constraints on the scattering phase shifts, should also provide an additional source of constraint on the low-energy extrapolation. This type of combined analysis has been reported in deBoer et al. [6], but there it was found that the available scattering data of Barnard et al. [52] was inconsistent with the capture data, perhaps because of incomplete uncertainty documentation in the former. With this in mind, new measurements of the 3 He(α, α) 3 He cross section were recently reported by Paneru et al. [53].
In this work, a Bayesian uncertainty analysis is performed on an R-matrix fit to the low energy 3 He(α, γ) 7 Be [37][38][39][40][41][42][43] and 3 He(α, α) 3 He [52,53] data, including the new SONIK data. The sensitivity of the fit to the scattering data is the main focus, examining the differences resulting from the two different scattering data sets considered. The mapping of the posterior distributions of the fit parameters, cross sections, phase shifts, and scattering lengths gives new insights into the dependence of these quantities to the input scattering data.

II. WHAT IS BRICK?
BRICK is a python package that acts as an interface between the AZURE2 [3, 34] R-matrix code and an MCMC sampler. It is not a replacement for AZURE2 nor is it intended to be. The primary functionality that it provides is a user-friendly way to sample parameters that have already been set up with the AZURE2 GUI to be varied.
AZURE2 is a multilevel, multichannel, R-matrix code (open source) that was developed under the Joint Institute for Nuclear Astrophysics (JINA) [3,34]. While the code was created primarily to handle the added complexity of charged-particle induced capture reactions [54], also has capability for a wide range of other types of reaction calculations. The code is primarily designed to be used by way of a graphical user interface (GUI), but can also be executed in a command line mode for batch processes [34]. The code stores all of its setup information in a simple text input file. While this file is usually edited by way of the GUI, it can also be modified directly. This may be desirable for batch type calculations, as are being used here.
AZURE2 primarily uses the alternative R-matrix parameterization of Brune [55]. It has two main advantages. The first is that it eliminates the need for the boundary conditions present in the classical formalism of Lane and Thomas [1]. The second is that the remaining fit parameters become the observed level parameters. The remaining model parameters are the channel radii.
A key advantage in using the parameterization of Brune [55] for the fitting of low energy capture reactions is that level parameters for bound or near threshold resonances can be more directly included in the R-matrix analysis [56,57]. The use of the Bayesian uncertainty estimation further facilitates the inclusion of uncertainty information for these parameters. This provides an improved method for communicating the level structure information gained from transfer reaction studies into an R-matrix analysis in a statistically rigorous way.
B. Implementation a. Overview The role of BRICK in our R-matrix calculations is to act as a mediator. It maps proposed parameters -both R-matrix parameters and normalization factors -from an MCMC sampler to AZURE2 and R-matrix predictions from AZURE2 back to the sampler. First, it accepts proposed points in parameter space, θ, from the sampler -in this analysis we use emcee [35] and packages them into a format that AZURE2 can read. Then it reads the output from AZURE2 and presents it as a list. Each item of the list contains the predictions, 7 Be E x (MeV) 6 Li + p 3 He + α µ(θ), and data, y and σ, corresponding to a specific output channel configuration. The likelihood can then be calculated according to the user's choice, see (2) for one possibility. Accompanied by user-defined priors, one can readily construct a Bayesian posterior. This posterior value, ln P, is passed back to emcee. Finally, based on the ln P value, the MCMC algorithm proposes a new θ and the process repeats. A diagram is provided in Fig. 2 to illustrate the qualitative functionality of the different software packages. The process described above starts at the orange rectangle labeled "emcee".
b. Details BRICK is built such that different samplers can be used. The analysis presented in this paper uses emcee, so the details provided in this section will be somewhat specific to it.
When initializing an instance of an EnsembleSampler, the most relevant argument is log prob fn, the function that returns the logarithm of the probability. One of the advantages of emcee is that it allows the practitioner to perform arbitrary calculations inside that probability function. That function must meet only two requirements: (1) take an array of floating point numbers that represents the vector in parameter space and (2) return a floating point number that represents the logarithm of the probability associated with that array. In between those two steps, one is free to perform whatever calculations one needs. This can be seen on the left-hand side of Fig. 2. The parameter-space vector, θ, is output from emcee. The logarithm of the probability at that point, ln P, is subsequently input to emcee. In this sense, emcee is well-suited to the implementation of "black-box" physics models where one has limited access to the source code.
The primary tasks that BRICK accomplishes are (1) translating θ into a format that AZURE2 can read and (2) reading the output from AZURE2 such that a ln P value can be easily calculated. The means of accomplishing these tasks relies on the command-line interface (CLI) to AZURE2, which is accessible when installed on Linux machines. The CLI options available to AZURE2 are well documented in the manual [34]. The most critical argument is the input file, typically accompanied by the file extension .azr. This input file contains all of the necessary information to perform an R-matrix calculation with a given set of parameters. It is generated when the R-matrix and data models are built with the commonly used GUI, which AZURE2 provides. BRICK is not built to replace that GUI. It accompanies AZURE2 by allowing the user to bring their AZURE2-prepared R-matrix model over and sample what was previously optimized. Accordingly, the default behavior of BRICK is to respect the choices made by the user in the AZURE2 GUI. If a parameter is fixed in AZURE2, it is fixed in BRICK. If it is varied in AZURE2, it is sampled in BRICK.
BRICK accesses the AZURE2 CLI through the Python module subprocess. But prior to that, BRICK must map the values in θ to the proper locations in the input file. This is accomplished by reading the <levels> and <sectionsData> sections of the input file. BRICK reads the appropriate parameters and flags looking for varied parameters. As they are found, their locations are stored. When a new θ is proposed, BRICK creates a new input file and maps the values in θ to the varied parameter locations. Then AZURE2 is called with the newly generated input file. The output from AZURE2 is written to a sequence of files in the output directory by default. Those files are read and the predictions, µ, and experimental data, y, are extracted. A likelihood is then constructed. Under the assumption that the uncertainties associated with y are uncorrelated and normally distributed, this is a multivariate Gaussian distribution. Accompanied by a list of prior distributions corresponding to the preexisting knowledge of the sampled parameters, a posterior is finally constructed and passed back to emcee.
Initially, this process was built in a single-threaded manner. As emcee is a ensemble sampler, efficient exploration of the posterior relies heavily on many, simultaneous walkers. In order to scale this beyond the most basic calculations, we modified our implementation to allow each walker to write its own input file and read from its own output directory. Inside the log-probability function, there is no access to any kind of walker identifier, so each walker generates a file-space that is uniquely identified by an eight-character random string. This allows each walker to work independently, so on systems where many cores are available, each walker can have a dedicated core. Or at least the time spent waiting for CPU time is minimized. This also allows for an increased number of walkers, which is a common tactic used to decrease autocorrelation time.
The starting point for the R-matrix model used here is that of deBoer et al. [6]. In that work, ten levels were used with three particle pairs ( 3 He+α, 7 Be+γ 0 , and 7 Be+γ 1 ) for a total of 16 R-matrix fit parameters. Initial MCMC calculations showed that a 7/2 − background level was not statistically significant, and was thus dropped from the calculation. This already demonstrates one of the powerful feature of this type of MCMC analysis, it provides a clear identification of redundant fit parameters. The R-matrix model used here thus consisted of nine levels, three particle pairs, and 16 R-matrix fit parameters as summarized in Table I.

B. Priors on R-matrix parameters
Because this is a Bayesian analysis, we must choose priors for all R-matrix parameters. We have chosen to employ uninformative uniform priors. However, the signs of the reduced width amplitudes (that is the interference solution), which are implemented in AZURE2 by the signs of the partial widths, were determined by the initial best χ 2 fit using AZURE2. In this case, a unique interference solution was found. This may not always be the case: sometimes other interference solutions may be possible. The emcee sampler may then not be able to easily find these other interference solutions in the parameter space. It seems is likely that in cases where different interference solutions are possible, each one will require a separate emcee analysis.
One common circumstance where a Bayesian analysis TABLE I. Sampled parameters in the R-matrix model. Numbers indicate that the level energies were fixed. A variable indicates that the parameter was sampled. The subscripts α and γ indicate the exit particle pair -scattering and capture, respectively. Capture particle pairs are distinguished by ground (0) and excited (1) 7 Be states.
J π E λ (Ex, MeV) Widths and ANCs Prior Distributions will improve on previous uncertainty estimates is in the ability to give priors for bound state level parameters determined from transfer studies. Unfortunately, in the case of the 7 Be system, there is limited information available for the bound state α-particle ANCs. A recent first measurement has been reported by Kiss et al. [58], but the ANCs are rather discrepant from those found from this and past R-matrix analyses of capture data. This inconsistency has not been investigated here, but needs to be addressed in future work. If reliable bound-state ANC determinations become available, that are independent of the capture and scattering data, it provides a path to further decrease the uncertainty in the low energy S-factor extrapolation. One could also adopt priors on the ANCs from ab initio calculations, although we have not done so. It is also tempting to implement more constraining priors into the R-matrix analysis from a compilation like the National Nuclear Data Center or the TUNL Nuclear Data Project [59]. However, great care must be taken to carefully understand the source of the values and uncertainties when weighted averages are used to determine adopted values for level parameters in these compilations. In particular, past analysis of the data being fit in the Rmatrix analysis may be a contributor to the evaluation values. Thus blindly using evaluation level parameters and uncertainties can lead to double counting and an erroneous decrease of uncertainties. It is for this reason that uniform priors on parameters are adopted in the present analysis. The posterior shapes then clearly stem solely from the data sets considered in the R-matrix analysis.
The priors for the R-matrix parameters used in this work are listed in Table I. In most cases, level energies are fixed. For details of the choices made in formation of the R-matrix model, see Paneru [60]. The distribution formed by the product of these R-matrix priors and priors on the parameters introduced in the next section is the overall prior π shown in Fig. 2.
C. Modeling systematic errors in the data a. Common-mode errors AZURE2 provides a method for the inclusion of a common-mode error for each data set using a modified χ 2 function where c α is the normalization fit parameter, n α is the starting normalization which is set to 1 in the present analysis, f (x α,j ) is the differential scattering cross section form the R-matrix, y α,j is the data point value, σ α,j is the combined statistical and point-to-point uncertainty of a data point, and δ cexp,α is the fractional commonmode uncertainty of the data set. The additional term in the χ 2 function is derived by making the approximation that the common-mode systematic uncertainty has a Gaussian probability distribution [61]. The accuracy of this approximation is often unclear.
Common-mode errors are implemented in the present analysis in BRICK, outside of AZURE2, i.e., the commonmode errors are applied to the AZURE2 output. In BRICK the R-matrix parameter set θ R is augmented by a set of normalization factors f α and energy shifts, ∆ E,α . (At present energy shifts are only implemented for scattering data.) The overall parameter set θ is then the union of the set θ R and {f α , ∆ E,α }. The likelihood L is formed as a product of standard Gaussian likelihoods for each data point, but with normalization factors applied to the AZURE2 predictions µ: (2) where we have omitted overall factors that do not affect the parameter estimation. Here x jα represents the kinematics of the jth data point in data set α, and σ jα is the combined statistical and point-to-point uncertainty of the corresponding datum, y jα . N α is the number of points in data set α, and the product over α runs over all the sets that have independent common-mode errors.
The priors on the f α 's are specified by the BRICK user. If a Gaussian prior centered at 1 with a width equal to the common-mode error reported in the original experimental publication is employed for the f α 's, then the product of that prior on the normalization factors and the likelihood (2) has the same maximum value as the "extended likelihood" corresponding to (1), that is used to estimate the f α 's in the frequentist framework implemented in AZURE2.
In our analysis of the 3 He(α,α) 3 He and 3 He(α,γ) 7 Be reactions, we adopted such a Gaussian prior, truncated to exclude negative values of the cross section. We used a different f α for each energy bin in the SONIK data, detailed in Section IV B, with the widths of the prior given by the common-mode errors stated in Table II. The common-mode error associated with the Barnard data, described in Section IV A, is taken to be 5%. The width of the priors for the f α 's to be applied to the capture data, discussed in Section IV C, are specified by the commonmode errors listed in Table III. All normalization-factor priors are of the form where and N (µ, σ 2 ) represents a Gaussian distribution centered at µ with a variance of σ 2 . b. Energy shifts BRICK also has the capability of estimating (overall) beam-energy shifts in a particular data set. This is implemented as another parameter to be estimated ∆ E,α . This parameter affects all the AZURE2 evaluations for data set α. BRICK implements the energy shift by generating a different input and data files for each value of ∆ E,α under consideration. The flowchart of Fig. 2 is thus not strictly accurate when this feature is included. Gaussian priors were defined, centered at zero, on possible energy shifts for the SONIK data and the Barnard data. The widths of the priors are based on information in the original papers, as summarized in Sections IV A and IV B. For the SONIK data, the energyshift parameter's prior has a standard deviation of 3 keV, based on the energy uncertainty quoted in Paneru et al. [53]. Barnard et al. [52] cites a much larger uncertainty of 20-40 keV, depending upon the energy. The standard deviation of the prior on the ∆E parameter is taken to be 40 keV for this data set, a much larger value than for the SONIK data. It should be noted that the energy uncertainty for the Barnard data set is not a constant, but it is not possible to improve our modeling of this uncertainty due to the lack of documentation of its origin. A new measurement of 3 He+α elastic scattering was performed at TRIUMF using the Scattering of Nuclei in Inverse Kinematics (SONIK) [62,63] target and detector system. SONIK was filled with 4 He gas maintained at a typical pressure of 5 Torr bombarded with 3 He with a beam intensity of about 10 12 pps. Elastic scattering cross sections were measured at nine different energies from E c.m. = 0.38-3.13 MeV. SONIK covers an angular range of 30 • < θ c.m. < 139 • -a markedly larger range than previous measurements. The detectors in SONIK were arranged such that they observed three different points, termed interaction regions, in the gas target along the beam direction. When the beam traversed the gas target it lost energy, so the bombarding energy, and therefore the scattering energy, was slightly different in each of the three interaction regions. As we will explore further below, the results for the differential scattering cross section from this measurement are consistent with previous determinations but have better precision. The data also extend to markedly lower energies. The uncertainties with this measurement are well quantified and are presented in Paneru et al. [53]. A separate normalization uncertainty is determined for each beam energy. These normalization uncertainties range from 4.1-9.8 %.
C. 3 He(α, γ) 7 Be data The data selection [37][38][39][40][41][42][43] for the 3 He(α, γ) 7 Be reaction for this work follows that of previous recent works [6,21,45,64]. Note that the LUNA measurements of Gyürky et al. [65] and Confortola et al. [66] are collected in Costantini et al. [39]. The combined data sets cover a wide energy range from E c.m. = 94 to 3130 keV, but still remain below the proton decay threshold. Older data are not included due to a long history of discrepancies, which manifested as differences between experiments that used either direct detection of γ-rays or the activation technique. More recent measurements have achieved consistency resulting from improved experimental techniques by performing consistency check measurements using both direct detection of γ-rays and the activation technique [45]. Details about the capture data sets, including common-mode errors for cross sec-tions, are listed in Table III.

D. Data Models
Two distinct data models are analyzed here, D CS and D CSB , where C indicates the inclusion of the capture data described in Sec. IV C, S indicates the inclusion of the SONIK data described in Sec. IV B, and B indicates the inclusion of the Barnard data described in Sec. IV A. D CSB is a more complete data model in the sense that it includes more data and would naively be considered the "best" data model. But, there are notable effects when the data of Barnard et al. [52] are included that are highlighted and discussed in Section V.

V. RESULTS
The results of our analysis are presented here in two subsections. The first discusses results in the energy regime of the data that was analyzed. The second computes extrapolated quantities -observables that lie in energy regimes outside those covered by the analyzed data.

A. Fits to data
First we examine the extent to which our results match experimental data. We do this by comparing predicted and measured observables. Figure 3 shows the total capture S-factor data alongside bands representing 68% intervals from the analyses of both data models, D CSB and D CS . For energies above 400 keV both analyses give very similar results. However, below that energy, the D CS analysis provides much better agreement with data. The LUNA data in particular discriminate between the two data models. The fit to the CSB data includes a normalization factor for the LUNA data that differs from 1 by about three times the stated common-mode error, cf. below. The normalization factors are not applied to the data in Fig. 3, which is why the CSB band sits well below the LUNA data.

Capture Data
Branching ratio results for both data models-D CS and D CSB -are shown in Fig. 4. The most prominent differences between the D CSB and D CS results occur near the upper and lower ends of the energy range. However, in the context of the experimental uncertainties, these differences are not significant. Over the entire energy range, the predictions from D CS and D CSB overlap at the 1-σ level. , LUNA [39], ERNA [43], and Notre Dame [37]. Colors, symbols, and line styles are the same as Fig 3. Bands indicate 68% intervals.

Scattering Data
The differential cross sections from the SONIK [53] and Barnard et al. [52] measurements are shown in Figs. 5 and 6, respectively, with the predictions from our analyses. In all cases, both analyses reproduce the data to high accuracy. However, the D CS analysis results in a much lower χ 2 /datum at max ln P : 0.72 for the SONIK [53] data vs. 0.95 for the D CSB analysis of the SONIK + Barnard [52] data sets.

B. Parameter Distributions
Separate corner plots for each data model are provided in the Supplemental Material [67]. There are notable differences in several R-matrix parameters. In particular, the D CS ANCs are significantly larger and their posterior distributions are noticeably wider. The D CS analysis also produces a significantly smaller ratio of ANCs, C 1 /C 0 . This is consistent with the smaller branching ratios at low energies shown in Fig. 4.
The D CS partial α widths in the 1/2 + , 3/2 + , and 5/2 + channels are smaller and separated by more than two standard deviations from the D CSB widths. The distributions for Γ (5/2 + ) γ,0 seem to indicate opposite signs. The posterior is markedly smaller and narrower, and the constraints on Γ (7/2 − ) α from D CSB are dramatically tighter. This is presumably due to the much larger amount of data in the vicinity of the 7/2 − resonance that is present in the Barnard et al. [52] data set. It is also worth noting the "non-Gaussian" behavior of several of these distributions-a characteristic that would be difficult to identify in a typical analysis that assumed linear propagation of uncertainties around a minimum of the posterior pdf. Using Gaussian approximations and linearizing would likely underestimate uncertainties in the 70 [53] are shown relative to the Rutherford prediction with grey x's and error bars. Each panel includes the measurements from three interaction regions [60]. Bands indicate 68% intervals. Green bands are generated for the analysis of DCS . Blue bands correspond to DCSB . case of Γ (3/2 + ) γ,0 , for example. All parameters shown in Fig. 7 are well-constrained. By comparing to the prior distributions listed in Table I, one can see the dominance of the data's influence over the information in the prior: all posterior distributions are markedly narrower than the priors chosen. As discussed in Sec. III B, several R-matrix-model iterations were taken to remove redundant parameters.
The correlation matrix of the R-matrix parameters is shown in Fig. 8. The figure represents an approximation of the full information contained in the corner plot given in the Supplemental Material [67]. There, significant, often-nonlinear, correlations are observable between several pairs of R-matrix parameters. In particular, the influence of the ANCs over the entire R-matrix parameter space, either directly or indirectly, means that it is very important for scattering data to have well-defined uncertainties over its full energy range.
The normalization factors applied to the theory predictions for each of the total capture data sets are shown for both data models in Fig. 9. The comparison reveals good agreement between D CS and D CSB for all but the LUNA data set [39]-the lowest-energy capture data set in our analysis. The D CS analysis yields a normalization factor for these data that is very close to 1. In contrast, the D CSB analysis requires that the LUNA data be shifted by nearly 10%. (Recall from Eq. (2), that f is applied to the theory prediction, and so an f > 1 corresponds to a systematic error that reduces the experimental cross section and uncertainties.) To put this in perspective, the LUNA collaboration estimates their common-mode error at 3.2%. Because the LUNA data set is the lowest capture data set, this disagreement between the D CS and D CSB analyses corresponds to a significant difference in the extrapolated S(0) of these two analyses.
The normalization factors applied to the theory predictions for each of the SONIK energies are shown in Fig. 10. When the data of Barnard et al. [52] are included in the analysis, the SONIK normalization factors are significantly larger. This effect is systematically apparent at lower energies. In more than half the cases, the D CSB and D CS results are inconsistent with each other. For eight out of ten SONIK energies, the normalization factor obtained from the fit is within the common-mode error estimated by the SONIK collaboration. Note that the common-mode error in this experiment was estimated to be different at different beam energies [60] 1 . This is IG. 6. Differential cross section as a function of energy as reported in Barnard et al. [52], shown as grey x's with error bars. Blue bands represent the 68% intervals generated from the DCSB analysis.  represented in Fig. 10 by the varying heights of the grey bands, which are priors in accord with these experimentally assigned common-mode errors, see Table II.
The posteriors for f Barnard and the energy shifts for both the Barnard et al. [52] and SONIK [53] data sets (see Sec. III C) are shown in Fig. 11. The result for f Barnard is 1.002 +0.003 −0.002 : well within the estimated systematic uncertainty of 5% given in Barnard et al. [52]. A shift of 19.26 +2.90 −2.51 keV in the energies reported in Barnard et al. [52] is found, but this result is consistent with the energy uncertainty estimates ranging from 20-40 keV given in that paper. However, even such a clearly nonzero shift does not seem to significantly impact extrapolated quantities. Finally, the SONIK energy shift indicated by our analyses is 1.59 +2.43 −1.81 keV. This result matches very well with the reported energy uncertainty estimate of 3 keV. The prior for this parameter was a normal distribution centered at 0 keV with a 1-σ width of 3 keV. The primary difference between the posterior and the prior for this parameter is the loss of probability in the negative energy region. If any energy shift in the SONIK data [53] is necessary, it is positive, but since 0 keV is well within one standard deviation, there is strong evidence for no shift.
The ANCs corresponding to the two bound 7 Be states are of particular interest for extrapolating threshold quantities. First, we point out that the inclusion of scattering data significantly reduces the uncertainty of the ANCs. Our posterior is much narrower than that obtained using capture-only data in Ref. Zhang 9. The normalization factors applied to the total cross section predicted by our R-matrix model are compared for each of the total capture data sets (Seattle [2] Weizmann [40], LUNA [39], ERNA [43], Notre Dame [37], and ATOMKI [42]). DCSB (blue) and DCS (green) results are shown together for each data set. Color online. This highlights the importance of scattering data in constraining bound-state properties and the amplitudes associated with transitions to them.
Second, the choice of scattering data set matters. The C 1 results from analyzing D CS and D CSB are discrepant at the 1-σ level. The C 0 results disagree by approximately 2-σ. The contrast is highlighted in Fig. 12 where the squares C 2 1 and C 2 0 are compared. The differing values directly impact the S-factor extrapolations discussed below.

C. Extrapolated quantities
The Coulomb-modified effective range function is given in Hamilton et al. [68] and van Haeringen [69] as where k is the relative momentum, is the angular momentum, η is the Sommerfeld parameter, Γ is the gamma function, u (η) is given by  with and Ψ representing the digamma function [70]. This effective range function is an analytic function of E (or k 2 ) near E = 0. From the phase shifts, obtained with BRICK, calculated over a range of low momenta, one can fit the scattering length, a 0 , and effective range, r 0 , according to the low-energy expansion Our calculation involves 70 equally spaced phase shifts over a range of energies from 0.57 keV to 3.93 MeV. The results are used to evaluate the effective range function defined by (5). The energy dependence is then fit to (9) using a non-linear least squares fit. The results from D CSB and D CS are shown in Fig. 13. As in the ANC comparison, they are strikingly discrepant.
The naive expectation would be that D CSB distributions would be smaller subsets of the D CS distributions. For many relevant quantities, this is not the case. Figure 15 shows a comparison of the scattering lengths obtained from the D CS and D CSB analyses. A comparison to Zhang et al. [21], also included in Figure 15, reveals the impact of including scattering data: the inclusion of scattering data drives the median downward and constrains the uncertainties significantly. A summary of these posteriors is given in Table IV.
The D CSB scattering length and effective range are both smaller and more tightly constrained. One might have expected that with more data-and more data at lower energies-this extrapolated quantity would become more tightly constrained. The two-dimensional posteriors shown in Fig. 13 seem to lie on the same line or band that defines the correlation between a 0 and r 0 , though two extended posteriors is not sufficient to define such a line.
The total capture S factor at zero energy was extrapolated by evaluating the S factor at 100 evenly 32.5 35.0 37.5 a 0 (fm) spaced points between 1 to 100 keV, constructing a cubicpolynomial interpolation function to represent the calculations, and evaluating that function at zero energy. Errors from the interpolation/extrapolation process are negligible when compared to contributions from parameter uncertainties. The results are shown alongside previous results in Fig. 14. As expected from the different low-energy behaviors shown in Fig. 3, the D CS and D CSB results are discrepant, only overlapping at the 2-σ level. The inclusion of the Barnard et al. [52] data reduces the uncertainty in S(0) and pulls the entire distribution downward, outside the uncertainties of the D CS analysis. This effect is not seen in [6] because the ANCs in that analysis were not varied freely. The D CSB result is discrepant with the D CS results and those reported in [6] and [21]. A summary of these posteriors is given in Table IV.
Insights into the relevance of parameters can be obtained by examining the correlations between them. In Fig. 16, the correlations between S(0) and a 0 , C 2 1 and C 2 0 are shown. While the D CS and D CSB results are discrepant in several astrophysically relevant cases, the discrepancy is consistent, and this figure exposes, to a large extent, why: the ANCs, particularly the groundstate ANC, strongly correlates with S(0). The Barnard et al. [52] data more tightly constrain these parameters at smaller values, and this directly lowers the predicted S(0) extrapolation.

VI. CONCLUSIONS
We have described and applied the Bayesian R-matrix Inference Code Kit (BRICK), which facilitates communication between the phenomenological R-matrix code AZURE2 [3] and a Markov Chain Monte Carlo (MCMC) sampler such as emcee [35]. It thereby enables MCMC sampling of the joint posterior probability density function (pdf) for the R-matrix parameters and normalization factors. With samples that represent such a posterior in hand, the computation of the pdf for any quantity that can be calculated in the R-matrix formalism is straightforward.
While BRICK is a general tool, we have also provided an example of its application to an R-matrix fit of 3 He-α scattering and the 3 He(α, γ) 7 Be capture reaction data, in order to make inferences about the 7 Be system. This application was partly motivated by the availability of a new 3 He-α scattering data set obtained using the SONIK detector at TRIUMF [60] following the suggestion of de-Boer et al. [6]. These data have more carefully quantified uncertainties than a previous measurement by Barnard et al. [52]. Our study shows this motivation was well justified, finding discrepant values for extrapolated quantifies when the data of Barnard et al. [52] were included. Our analysis of the SONIK data shows consistency between them and capture data, producing an S factor in accord with analyses of capture data alone: our final D CS (capture + SONIK data) result for the S-factor at zero energy is S(0) = 0.539 +0.011 −0.012 keV b. When the   Barnard et al. [52] data were included in the analysis, the D CSB results produced significantly lower ANCs and S(0) extrapolation. Indeed, the D CSB analysis produces values for S(E) at c.m. energies of 10-20 keV that can only be reconciled with the LUNA data [39] if the normalization of these data is adjusted by 2-3 times the quoted common-mode error.
This emphasizes the importance of detailed uncertainty quantification when data sets are to be used for accurate inference of extrapolated quantities, where Barnard et al. [52] does not include these kinds of details regarding the experiment. This makes the tension between the Barnard et al. [52] and SONIK data regarding S(0) difficult to resolve, thus the Barnard et al. [52] data may need to be omitted from future evaluations. We emphasize, though, that these previous data were invaluable in advancing our understanding of the 7 Be system to its current state, but data with more well defined uncertainties are needed for current applications.
Zhang, Nollett, and Phillips pointed out that the s-wave 3 He-α scattering length is correlated with this result [21]. The D CS analysis produces a 0 = 36.59 +0.55 −0.53 fm. Premarathna and Rupak simultaneously analysed capture data and 3 He-α phase shifts in EFT and found a 0 = 40 +5 −6 fm (Model A II of Ref. [20])-in good agreement with this number. However, it disagrees by 2σ with the a 0 extracted using EFT methods from capture data alone by Zhang et al. [21]: a 0 = 50 +6 −7 . Recently Poudel and Phillips [71] performed an EFT analysis of the SONIK data, using priors on the 7 Be ANCs from the capture analyis of Ref. [21], and extracted a 0 = 60 ± 6 fm-even further away from the results of this R-matrix analysis.
Improvements in the analyses presented here could occur if there were: • Better documentation of the energy dependence of systematic uncertainties in published data sets. The Bayesian formalism that underlies BRICK allows systematic uncertainties with any correlation structure to be incorporated into the analysis.
• Improved understanding of the way theory uncertainties in the phenomenological R-matrix formalism affect the extrapolation of data.
• Detailed modern data with full uncertainty quantification in the vicinity of the 7/2 − resonance. This may help resolve some of the ambiguities in results between the D CS and D CSB analyses.
• Ab initio constraints, e.g., on ANCs could be incorporated in the analysis.
• Data from transfer reactions that provided complementary information on the 7 Be ANCs.
Future applications of BRICK could include posteriors for astrophysical reaction rates. This would enhance BRICK's utility as a tool for performing detailed uncertainty quantification on nuclear reactions, especially those of astrophysical interest. AZURE2 already includes the necessary functionality. Implementing this feature ought to be a straightforward process.   2.