Impact of systematic and amplitude model correlations within and between systems of combined input: A case study with $\phi_2$ ($\alpha$)

The pursuit of experimental precision in the $CP$-violating weak phase $\phi_2$ ($\alpha$) is not without its challenges, in part due to the need to combine multiple physical observables from various related decay channels, and therein lies a fundamental issue. Similarities in analysis procedures give rise to systematic correlations between the measured inputs constraining $\phi_2$ that must be taken into account to avoid bias. Specifically, in the case of the irreducible model uncertainty accompanying analyses involving the $\rho$ meson, it is demonstrated that ignoring correlations derived from its pole parameters, or indeed even treating correlations individually contained within each decay channel, can ultimately lead to a bias in $\phi_2$ of $\mathcal{O}(1^\circ)$. Correct treatment on the other hand, markedly reduces wandering of its central value as a function of the model uncertainty strength with the added dividend of a further improved overall uncertainty. Bias in the combination of $B^0 \to (\rho \pi)^0$ and $B \to \rho \rho$ is also seen to depend on the statistical strength of the former in relation to that of the model uncertainty in the latter. This work can inspire other studies into the points at which systematic correlations beyond those determined in single measurements matter in combinations leading to other $CP$-violating weak phases such as $\phi_1$ ($\beta$), $\phi_3$ ($\gamma$) and $\phi_s$.


Introduction
Violation of the combined charge-parity symmetry (CP violation) in the Standard Model (SM) arises from a single irreducible phase in the Cabibbo-Kobayashi-Maskawa (CKM) quarkmixing matrix [1,2]. Various processes offer different yet complementary insight into this phase, which manifests in a number of experimental observables over-constraining the Unitarity Triangle (UT). The measurement of such parameters and their subsequent combination is important as New Physics (NP) contributions can present themselves as an inconsistency within the triangle paradigm.
Decays that proceed predominantly through the b → uud tree transition (figure 1a) in the presence of B 0 -B 0 mixing are sensitive to the interior angle of the UT, φ 2 = α ≡  Figure 1. Leading-order Feynman diagrams shown producing B 0 → π + π − decays, though the same quark transition can also produce B 0 → ρ ± π ∓ , ρ + ρ − and a ± 1 π ∓ . (a) depicts the dominant (tree) diagram while (b) shows the competing loop (penguin) diagram. In the penguin diagram, the subscript x in V xb refers to the flavour of the intermediate-state quark (x = u, c, t).
More often than not, this involves the combination of several physics parameters extracted from related decay channels raising an important general question that hitherto has not yet been explored in any great detail. As experimental measurements become more and more precise, a crucial unknown is the point at which it will become necessary to consider systematic correlations, not arising simply within individual analyses, but rather in between the relevant analyses in order to avoid bias. As a specific case study of a broader issue, which includes but is not limited to combination-based approaches designed to measure other CP -violating weak phases such as φ 1 (β) [27][28][29][30][31][32][33][34], φ 3 (γ) [35][36][37] and φ s [28,31,34,[38][39][40], I open the discussion in this paper within the context of the φ 2 average. By and large, this problem is generally an internal matter for each collaboration, however there are irreducible systematic uncertainties that transcend experiment, warranting a more cooperative approach and thus is the primary focus here. Experience in amplitude analysis suggests that the model uncertainty of a dominant vector resonance tends to derive more significantly from its own pole parameters rather than the remainder of the model, which is converse to smaller contributions where the opposite trend appears to hold. This is because Breit-Wigner phases vary most rapidly at the poles, exacerbating the effect of small variations to manipulate interference patterns in the regions of greatest physical interest. In this specific consideration, these pole parameters are those of the ρ meson.
I open in section 2, with a description of the SU(2)-based approach for controlling distortions in experimental φ 2 measurements arising from the ever-present strong-loop gluonic penguin processes. Following this, I introduce the impact of systematic correlations in section 3 with a conceptually simpler example surrounding the branching fractions of the decay processes involved. In section 4, I move into the primary study on the bias in φ 2 caused by neglecting amplitude model correlations in the B → ρρ system arising from the ρ pole masses and widths. This bias, if left unchecked, can then go on to affect the otherwise immune B 0 → (ρπ) 0 analysis as explained in section 5. Finally, conclusions are drawn in section 6 along with some recommendations on how the community going forward can reduce φ 2 bias induced by systematic and amplitude model correlations.
2 Strong-penguin containment in φ 2 constraints In general, the extraction of φ 2 is complicated by the presence of interfering amplitudes that distort the experimentally determined value of φ 2 from its SM expectation and would mask any NP phase if not accounted for. These effects primarily include b → duu strongloop decays (figure 1b), although isospin-violating processes such as electroweak penguins, π 0 -η-η mixing, ρ 0 -ω-φ mixing [41] and the finite ρ width in B → ρρ [42] can also play a role.

Original approach
It is possible to remove the isospin-conserving component of this contamination by invoking SU(2) arguments. The original method considers the three possible charge configurations of B → ππ decays [43]. Bose-Einstein statistics rules out a total isospin I = 1 contribution, leaving just the I = 0, 2 amplitudes remaining. Strong penguins then only have the possibility to contribute an I = 0 amplitude, since the mediating gluon is an isospin singlet. However, in the specific case of B + → π + π 0 , the further limiting projection I 3 = 1 additionally rules out I = 0, thereby forbidding strong penguin contributions to this channel.  The complex B → ππ andB → ππ decay amplitudes obey the isospin relations respectively, where the superscripts refer to the combination of pion charges. The decay amplitudes can be represented as triangles in the complex plane as shown in figure 2. As B + → π + π 0 is a pure tree mode, its amplitude in isospin space is identical to that of its CP -conjugate and so these triangles lose their relative orientation to share the same base, A +0 =Ā +0 , allowing the shift in φ 2 caused by strong penguin contributions ∆φ 2 ≡ φ ± 2 −φ 2 , to be determined from the phase difference betweenĀ +− and A +− . These amplitudes can be constrained by 7 mostly independent physical observables for a two-fold discrete ambiguity in the range [0, 180] • , which are related to the decay amplitudes as where B, A CP and S CP are the branching fractions, CP violation in the decay and mixinginduced CP violation parameters, respectively. The superscript ij, represents the charge configuration of the final state pions and τ i+j B is the lifetime of the B + (i + j = 1) or B 0 (i + j = 0). Naturally for B + → π + π 0 , CP violation in the decay is forbidden by the isospin argument and mixing-induced CP violation is not defined. The ambiguity in φ 2 is also increased to 8-fold if S 00 CP of the colour-suppressed channel is not measured as is currently the case. This approach is also applied to the B → ρρ system analogously, substituting the ρ meson in place of each pion.

Next-generation approach
The B → ρρ system presents a greater theoretical and experimental challenge over B → ππ. It has already been pointed out that isospin-breaking (I = 1) ρ-width effects can be controlled by reducing the ρ analysis window of B 0 → ρ + ρ − and B + → ρ + ρ 0 according to the method outlined in ref. [44]. An open question to be studied is the extent to which this is systematically feasible in the presence of interfering and non-interfering backgrounds.
In this work, I espouse an alternate viewpoint in which the possibility to exploit the multi-body final state through directly modelling the structure of ρ 0 -ω mixing and I = 1 finite ρ-width effects is acquired in exchange for greater analysis complexity. To that end, I have already outlined the amplitude analysis framework by which this can be achieved, replacing the measured physical observables from eq. 2.2 by where λ ij CP is a CP -violation parameter and φ ij 2 is its effective weak phase. As these quantities are now related to the isospin triangles at amplitude level, the solution degeneracy in φ 2 for the range [0, 180] • is resolved [45] and as an added incentive, the 8-fold solution degeneracy in B 0 → a ± 1 π ∓ can also be lifted for the same range in the SU(3) approach [46]. Naturally, this method raises concerns regarding the potential impact on φ 2 coming from correlated amplitude model systematics which will be studied here.

Statistical method
In this paper, I employ the frequentist approach adopted by the CKMfitter Group [24] where a χ 2 is constructed comparing theoretical forms for physical observables expressed in terms of parameters of interest, µ, with their experimentally measured values, x. The most general form, is necessary here, where Σ is total covariance matrix composed of the statistical and systematic covariance matrices as Σ ≡ Σ Stat + Σ Syst . The statistical covariance matrix comes directly from the function minimisation procedure during the nominal fit to a sample, while the systematic covariance matrix is manually derived. Parameter variations are generated according to their uncertainties and the fit is repeated for each set of variations. Over N fits, the covariance between a pair of physical observables is given by where the barred quantities representing the means are obtained from the nominal fit. A scan is then performed, minimising the χ 2 to determine µ for each value of φ 2 fixed across a range. The value of ∆χ 2 from the global minimum is finally converted into a p-value scan, assuming it is distributed with one degree of freedom, from which confidence intervals can be derived.

Systematic correlations within systems
Before delving into the main point regarding amplitude model correlations, it may be advantageous to introduce this difficult topic by digressing to conceptually simpler systematic correlations that can be trivially accounted for here. One such example is the number of BB pairs produced at e + e − machines operating at the Υ(4S) resonance, N BB , that enters the absolute branching fraction calculations in B → ππ decays through where N ij is the extracted signal yield and ij is the reconstruction efficiency of that mode. Although equal production of B + B − and B 0 B 0 pairs is implicitly assumed here for simplicity, this will need to be evaluated at Belle II as the current uncertainties on their rates [56] would otherwise constitute the dominant systematic instead of those arising from N BB . It can immediately be seen that all three branching fractions are 100% systematically correlated in N BB , because as a quantity that is independent of the channel being studied, whatever direction it fluctuates in, all branching fractions must follow suit by the same factor. To illustrate, I repeat the φ 2 projection for Belle II with 50 ab −1 with and without accounting for systematic correlation arising from N BB . Input is borrowed from ref.
[47] where the first uncertainty is statistical and the second is systematic. matrix can be immediately written down as shown in table 2. The only known statistical correlation is between A CP (π + π − ) and S CP (π + π − ), set at +0.10 from the Belle result. The φ 2 scan is then conducted in the vicinity of the SM solution with and without systematic correlations, the results of which can be seen in figure 3. The leading edge of the solution consistent with the SM is seen to improve by 0.4 • when accounting for systematic correlations, a striking result within the context of the sub-degree precision anticipated at Belle II. At a first glance, this may seem counter-intuitive as some may recall the familiar summation of correlated uncertainties linearly over the more favourable summation in quadrature for uncorrelated cases. However, this is more applicable to instances of single physics parameters, whereas between physics parameters, correlations restrict statistical freedom. In this example, as the uncertainty in N BB is not allowed to nonsensically follow three independent statistical distributions, as would be encapsulated by the identity correlation matrix, the φ 2 constraint improves in consequence.
Although the knowledge of N BB is expected to dominate the systematic uncertainty on the branching fractions, ideally the full systematic covariance matrix will be constructed in future analyses considering all sources. For example, as common control samples provide the tracking, particle identification and π 0 reconstruction uncertainties, the branching fractions are again systematically correlated in these categories. Concurrently, the CP -violating parameters are also affected with the timing resolution and flavour-tagging performance being obtained from common studies. However, perhaps the most dangerous systematic here would be the shared method for evaluating tag-side interference from doubly-Cabibbosuppressed decays [48], which is also considered to be another irreducible systematic up until

Amplitude model correlations within systems
Although amplitude analysis has seen limited involvement [19] in the B → ρρ constraint of φ 2 , it stands to reason that this approach will become more attractive in controlling uncertainties as data samples increase. For the small cost of modelling one additional variable over current analyses, the necessary degrees of freedom are acquired to harness the full statistics particularly of the limiting colour-suppressed B 0 → ρ 0 ρ 0 decay, thereby improving φ 2 precision in this sector, and even opening the possibility to determine φ 2 separately for each of the three polarisation configurations of B → ρρ. Furthermore, amplitude analysis allows the direct modelling of interfering components such as the dipion I = 0, 1 resonant contributions and other non-resonant S-wave effects such as elastic and inelastic particle rescattering processes in the vicinity of the ρ, thus reducing model uncertainty estimations. Perhaps most importantly, isospin-breaking contributions known to bias φ 2 can also be accounted for in the amplitude model, such as with the ρ 0 -ω mixing lineshape of refs. [49,50] and the structure of I = 1 finite ρ-width effects suggested in ref. [42].
One aspect these three analyses have in common is fixed ρ pole parameters, so therein lies the potential for systematic model uncertainties to impact the φ 2 average. Unlike the N BB case mentioned in the previous section, these are not multiplicative factors to any physics parameter and as such the correlation matrix cannot immediately be written down. In order for the covariance matrix to be derived, repeated randomised systematic variations on a sample of each B → ρρ channel needs to be applied. Then to generate these samples, amplitude models are first required for which information is sparse, meaning that assumptions will have to be made on the magnitudes and relative phases between the B → ρρ polarisations. In previous works [45,46], unknown physics parameters were uniformly distributed in an ensemble test to give a sense of what to expect on average in their respective φ 2 studies. However in this case, applying systematic variations on top of all three amplitude analyses in an ensemble is not a practical endeavour and would be of unclear benefit, besides. Therefore, conclusions from this paper will be limited to identifying the scale of potential bias in φ 2 induced by neglecting amplitude model correlations, as opposed to providing a definitive range.

Amplitude model
Yields are set based on Belle results according to expectations for 50 ab −1 to be collected with Belle II. I consider rudimentary models with contributions to the 4-body phase space coming only from the channels known to exist in the analysis region. The amplitude for each intermediate state at position Φ 4 , is parameterised as where B L B represents the production Blatt-Weisskopf barrier factor [51] depending on the orbital angular momentum between the products of the B decays, L. Two resonances will appear in each isobar, denoted by R 1 and R 2 , for which respective decay barrier factors are also assigned. The Breit-Wigner propagators are represented by T , while the overall spin amplitude is given by S. Each isobar is Bose-symmetrised as necessary so that the total amplitude is always symmetric under the exchange of like-sign pions.
The Blatt-Weisskopf penetration factors account for the finite size of the decaying resonances by assuming a square-well interaction potential with radius r. They depend on the breakup momentum between the decay products q, and the orbital angular momentum between them L. Their explicit expressions used in this analysis are B 0 (q) = 1, Spin amplitudes are constructed with the covariant tensor formalism based on the Rarita-Schwinger conditions [52]. The spin S, of some state with 4-momentum p, and spin projection s z , is represented by a rank-S polarisation tensor that is symmetric, traceless and orthogonal to p. These conditions reduce the number of independent elements to 2S + 1 in accordance with the number of degrees of freedom available to a spin-S state. The sum over these polarisation indices of the inner product of polarisation tensors form the fundamental basis on which all spin amplitudes are built. Called the spin projection operator P , it projects an arbitrary tensor onto the subspace spanned by the spin projections of the spin-S state.
Another particularly useful object is the relative orbital angular momentum spin tensor L, which for some process R → P 1 P 2 , is the relative momenta of the decay products q R ≡ p 1 − p 2 projected to align with the spin of R, where the number of indices representing the tensor rank is equal to the value of L. Finally, to ensure that the spin amplitude behaves correctly under parity transformation, it is sometimes necessary to include the Levi-Civita totally antisymmetric tensor abcd p d R . Each stage of a decay is represented by a Lorentz scalar obtained by contracting an orbital tensor between the decay products with a spin wavefunction of equal rank representing the final state. Three spin topologies are necessary for B → ρρ as S-, P -and D-waves are permitted between the vector resonances, with total spin densities, In general, resonance lineshapes are described by Breit-Wigner propagators as a function of the energy-squared s, where M 2 (s) is the energy-dependent mass and Γ(s) is the total width which is normalised such that it represents the nominal width Γ 0 , at the pole mass, m 0 . For the ρ 0 resonance, the Gounaris-Sakurai parameterisation is used to provide an analytic expression for M 2 (s) and Γ(s) [53].

Pseudo-experiment generation method
The unknown strong complex couplings between contributions in the amplitude model are partly inspired by reverse-engineering the known branching fractions for each polarisation. The Monte Carlo (MC) is based on the decay rates in phase space, which for where q = +1(−1) for B + (B − ). On the other hand, the time-dependent decay rates of B 0 and B 0 decays to a self-conjugate final states are given by respectively, where A is the static decay amplitude, τ is the B 0 lifetime and ∆m d is the mass difference between the B H and B L mass eigenstates. This form assumes no CP violation in the mixing |q/p| = 1, and that the total decay rate difference between the two mass eigenstates is negligible. The total amplitude A, can be written in the typical isobar approach as the coherent sum over the number of intermediate states in the model with amplitude A i , as a function of 4-body phase space position Φ 4 , where a i is a strong complex coupling determined directly from the data. Incorporating a complex CP violation parameter λ i , for each weak contribution in the phase space, the totalĀ can be written asĀ for B + → ρ + ρ 0 , where the phase space of the CP -conjugated processΦ 4 , is set by convention to have the same sign as Φ 4 for all amplitude contributions, leaving A i to contain only strong dynamics blind to flavour. Conversely, for B 0 → ρ + ρ − and ρ 0 ρ 0 , the phase space of the CP -conjugated processΦ 4 , must be transformed relative to the elected particle ordering under C and P conjugation in order to achieve A i containing only strong dynamics. Complex couplings are then evaluated through a χ 2 fit relating the observed branching fractions for each isobar scaled to unity, to the fit fractions of each isobar calculated for the generated model in the 4-body phase space, The branching fractions for each polarisation are set mostly with HFLAV input [26] except where mentioned. As the longitudinal polarisation is known to dominate, the remainder is assigned exclusively to the P -or CP -odd P-wave for simplicity, while the longitudinal component is divided evenly between the P -or CP -even S-and D-waves for the flavourspecific and CP -conjugate final states, respectively. Naturally, there are 2 solutions for each free strong coupling, so whichever solution the fit converges to first is taken to generate the MC sample for each B → ρρ channel. Position in phase space is provided by the GENBOD algorithm [54] and qft++ gives the spin densities [55].
This amplitude analysis should be rather straight-forward as the phase space can be restricted to limit contributions from the a 1 (1260) resonances. Here, the yield is set to 100k events and the dipion range is restricted to be below the typical 1.1 GeV/c 2 . The input branching fractions with the determined couplings are shown in table 3. Note that here and throughout, the spin amplitudes given in eq. 4.4 are not normalised over the phase space, so there is no direct relation between the fitted couplings and their corresponding branching fractions. This also means that the relative strengths of each partial wave cannot be inferred from the couplings either as each spin factor contains different momentum scales by eq. 4.3, depending on the number of orbital angular momentum tensors involved.  Table 3. Branching fraction input with corresponding reverse-engineered couplings for B + → ρ + ρ 0 .
fractions are similar, dilution arising from imperfect flavour-tagging performance taken to be around the 30% mark for Belle II means that the yield is set to 30k events. The CPviolation parameter can be set from the known longitudinal quasi-two-body parameters assuming the solution closest to the SM expectation and uniformity between polarisations. The input branching fractions with the determined couplings and CP -violation parameters are shown in table 4. It should also be noted that despite an amplitude analysis being conducted here, it is very unlikely that a single solution for the effective φ +− 2 will emerge. As discussed in ref. [45], a lack of interfering contributions with a sizeable penguin contribution means that the CP -violation parameter will likely factorise in the isobar sum such that (ĀA * ) = (λAA * ) = (λ|A| 2 ) in eq. 4.7. As |A| 2 must be real-defined, the imaginary part of the aforementioned product evaluates to sin 2φ +− 2 , leaving two solutions remaining for φ +− 2 .

B
According to ref. [45], expanding the analysis space to include B 0 → a ± 1 π ∓ can ultimately resolve the φ 2 solution degeneracy in B → ρρ. Novelty aside, this strategy is prudent as the colour-favoured B 0 → a ± 1 π ∓ is otherwise either difficult to control systematically or statistically very expensive to remove. As such, the analysis range is defined as the dipion mass being less than 1.1 GeV/c 2 as before, or the 3-pion mass being below the production of open charm. The combined yield accounting for flavour-tagging performance is estimated in that previous work at 30k events.
An additional topology is necessary for B 0 → a ± 1 π ∓ , arising from the S-wave between the products of the a ± 1 → ρ 0 π ± decay, While a relative orbital angular momentum D-wave between the vector and pseudoscalar is possible, this has yet to be definitively seen, so is ignored at this time. Regarding the lineshape of the a ± 1 , potential dispersive effects are neglected, setting M 2 (s) to its polemass squared. The energy-dependent width of the a ± 1 is calculated from the integral over its phase space as a function of s, where A is the transition amplitude of the cascade, itself being comprised of barrier factors, a spin density and lineshape, with a coherent sum taken over the open polarisation indices of the initial state. Its exclusive decay to (ρπ) ± and isospin symmetry i.e., Γ a ± 1 →ρ 0 π ± (s) = Γ a ± 1 →ρ ± π 0 (s) is also assumed. The numerical form of the a ± 1 energy-dependent width can be seen in figure 4.  The input branching fractions with the determined couplings and CP -violation parameters are shown in table 5. In this case, only the LHCb results are used to set the B 0 → ρ 0 ρ 0 branching fractions [19], while the Belle result is used to determine the B 0 → a ± 1 π ∓ model [21].  Table 5. Branching fraction input with corresponding reverse-engineered couplings for B 0 → ρ 0 ρ 0 and B 0 → a ± 1 π ∓ .

Results
An amplitude fit is performed to the MC sample for each decay in order to determine nominal values for the physics observables of interest against which to compare in the ensuing systematic variations. Three scenarios are defined: "Current practice", in which systematic correlations are not considered at all, "Expected practice" for which each analysis independently handles their own systematic correlations and finally "Proposed practice", where systematic correlations are globally accounted for. Sets of ρ pole parameters are then generated by which to refit the MC samples and calculate covariance between physics observables. Although the ρ + and ρ 0 parameters are separately determined, they are essentially a manifestation of the same state with different charge, and so for the purposes of evaluating model uncertainties here, are also considered to be fully correlated at the theoretical level. Their parameters taken from ref. [56] are thus distributed with a multivariate normal distribution, recorded in table 6. Each MC sample is then refit applying each ρ variation upon which amplitude model covariance matrices are constructed for each scenario. Obviously for the Current practice, the amplitude model correlation matrix is set to the identity. In the Expected practice scenario, 100 ρ parameter variations are individually generated for each decay starting with a unique seed, while for the Proposed practice, a single final set of 100 variations is generated to be shared amongst the three analyses. The physics observables included are the branching fractions and where appropriate, the magnitudes of the CP -violating parameters and their effective weak phases for each polarisation, resulting in a 21 × 21 covariance matrix. For completeness, the model uncertainties and their corresponding correlation matrix are given in appendix A.
The φ 2 constraint summing over all polarisations is then conducted in each scenario. The model uncertainty is not taken to scale with statistics considering the large yields involved because the variations of the ρ parameters change the underlying interference pattern in the phase space in a predetermined albeit unknown way, which is why this systematic is irreducible without appreciable improvements to the pole properties themselves. As such, I repeat the φ 2 constraint in the context of a changing statistical uncertainty which is achieved by scaling the statistical covariance matrix obtained from the nominal fits, with the results shown in figure 5.
The raw drift of φ 2 can be seen in figure 5a. For reference, a statistical-only constraint ignoring the amplitude model covariance matrix is performed for which φ 2 is perfectly flat as a function of its own overall statistical uncertainty as expected. For these values of statistical uncertainty, the trends of each analysis scenario are then determined. When accounting for the model uncertainty, but not any corresponding correlations, the constraint has no problem at the statistical precision of 1 • , but rapidly deteriorates to become the worst performer once the model uncertainty begins to dominate. Surprisingly, the Expected practice scenario in which each analysis considers their own model correlations already has a visible bias at 1 • uncertainty, trending in the same direction as Current practice but without any intuitively discernible motion. Finally, the Proposed practice curve tightly follows the statistical-only curve, indicating that its covariance matrix by and large captures the correct structure of the model uncertainty. Perhaps the slight departure from the flat curve at low values indicates the point at which sensitivity to non-linear correlations begins to play a role. Figure 5b indicates the size of the model uncertainty as a percentage of the total uncertainty in φ 2 , which is obtained by quadratically subtracting the statistical uncertainty at that point, thereby assuming that all other systematics behave like the statistical error. Current practice performs well when the statistical uncertainty dominates, but becomes overtly large at the other end of the spectrum. Conversely, the model uncertainty fraction for Expected practice is noticeably larger even where the statistical uncertainty is supposed to be dominant, while its performance even improves to the point of the Proposed practice though with substantial bias as already seen.
The plots of figure 6 are variations of those already shown, though now as a function of the model uncertainty fraction in the Current practice scenario and relative to Proposed practice, which allows the difference in performance to be understood at the Belle II projection of 0.7 • total uncertainty with 50 ab −1 [47]. This is achieved by tuning the statistical covariance matrix scale factor, the result of which is shown in figure 7. Back in figure 6a, the bias in φ 2 with the full Belle II dataset would be only +0.03 • with Current practice, but −0.1 • with Expected practice for the studied model. Figure 6b indicates how much worse the model uncertainty scales as a function of its fraction, showing that at Belle II uncertainties, Current practice is worse-off than Proposed practice only by factor of 1.2, while Expected practice has a model uncertainty 3 times worse.
For the given amplitude models studied, it would appear that Proposed practice is objectively superior. Current practice may initially perform well, but is expected to become problematic, while Expected practice just looks dangerous from the onset, which is perhaps an unexpected outcome. A natural question arising at this point is whether this whole issue  can be sidestepped by releasing the ρ pole properties in the fit. However, the meaning of the φ 2 constraint would be unclear with multiple versions of these parameters present. In principle, it would be possible to release the ρ pole parameters in a simultaneous fit to the three B → ρρ modes, trading relative analysis simplicity for easier sytematics handling, though understandably this alternative is not very practical in terms of coordination and the short-term contractual nature of the field. In any case, there is more to each amplitude model than just the B → ρρ contributions as mentioned at the beginning of this section, so the full correlated model uncertainty accounting for other common lineshapes and additional contributions should also be studied. Minimal cooperation and overlap between analyses through the sharing of systematic variations and their signed fit residuals would seem to be the most sensible strategy moving forwards.

Amplitude model correlations between systems
Now that systematic bias in φ 2 and precision being left on the table is found to occur when not considering amplitude model correlations properly, or indeed even at all, attention turns to the other system also involving the ρ lineshape in an amplitude analysis, B 0 → (ρπ) 0 . While this involves a single analysis in which the complete set of systematic uncertainties are already considered as standard, the question remains whether the global treatment of amplitude model correlations in B → ρρ is sufficient, or whether its scope should also be expanded to include B 0 → (ρπ) 0 . Recalling the discussion of section 3, this is analogous to the realisation that all branching fractions of the B → ππ and ρρ systems are also systematically correlated through N BB .

Amplitude model
The overlap of the three B 0 → (ρπ) 0 charge combinations in the B 0 → π + π − π 0 phase space allows the tree (T ) and penguin (P ) processes involved to be distinguished, leading to the direct measurement of φ 2 in a single analysis [57]. Additional constraints from isospin pentagonal relations involving the charged B modes can also help to improve the constraint [58], though this extension will not be addressed further here, suffice it to say that a global approach to amplitude model correlations in the wider B → ρπ system will likely be needed as well in light of what has been seen so far. The single analysis applies the isospin symmetry argument to the penguin amplitudes rather than the tree as is the case in other systems, lessening the impact of theoretical uncertainties [41]. The decomposition of the complex amplitude couplings by charge and B flavour is given in eq. 5.1, A +− = T +− e −iφ 2 + P +− ,Ā −+ = T +− e +iφ 2 + P +− , A −+ = T −+ e −iφ 2 + P −+ ,Ā +− = T −+ e +iφ 2 + P −+ , A 00 = T 00 e −iφ 2 + 1 2 (P +− + P −+ ),Ā 00 = T 00 e +iφ 2 + 1 2 (P +− + P −+ ), (5.1) where the superscripts represent the ρ followed by the pion charge. It can be seen that sensitivity to the penguin amplitudes amongst the trees arises as a result of the isospin argument removing the independence of the colour-suppressed B 0 → ρ 0 π 0 penguin. However, by construction the tree and penguin couplings are highly correlated, so a parameterisation expanding the product of isobar sums for each of the three terms in eq. 4.7 was proposed [59]. For each resulting amplitude-squared-level form factor composed solely of strong dynamics, the expression for the corresponding coupling combination is then substituted with an independent free parameter. The parameter space thus increases dramatically in return for statistical stability of the fit, after which a χ 2 minimisation mapping the original tree and penguin amplitudes to the 27 bilinear coefficients can be executed to recover φ 2 . At this time, it is unclear if this method will persist going forwards due to the unphysical assumption that the higher ρ resonances are imposed to share the same P/T ratios and the related impracticalities in adding further structures to the amplitude model.
For the purposes of this study, since the tree and penguin amplitudes are known from MC generation, the fit will be performed with that paradigm for simplicity as φ 2 will be a free parameter. The input amplitudes themselves are still obtained from a χ 2 fit to the Belle bilinear coefficients, taking the solution consistent with SM expectations [11,12]. The parameters are given in table 7.

Results
The B 0 → (ρπ) 0 effective yield accounting for flavour-tagging dilution is set at 30k events for the 50 ab −1 expected at Belle II, with two scenarios considered for systematic variations. The first is that each system is treated independently, with B → ρρ already adopting the globally treated correlations proposed in the previous section as standard, while in the second, the systematic variations are common to both systems. The model uncertainty on φ 2 and how it correlates with the B → ρρ parameters are again given in appendix A. A φ 2 scan is performed over the lesser known quantities, which again includes the strength of the B → ρρ model uncertainty in that system and the total uncertainty in B 0 → (ρπ) 0 relative to B → ρρ. Figure 8a shows the scope of the bias in φ 2 , while Figure 8b demonstrates the loss in statistical power when treating each system independently. From these plots, it is clear that a coordinated approach to these analyses would be required if B 0 → (ρπ) 0 dominates the precision when at the same time a sizeable model uncertainty is present in the B → ρρ system.

A word on amplitude model correlations between experiments
So far, no mention has been made on the role of LHCb, which at this time is expected to provide only partial input to the φ 2 constraint, primarily from the B 0 → ρ 0 ρ 0 decay. Depending on the statistical power of this analysis, which will be considerable, it has already been shown that each analysis and therefore experiment by extension, when left to their The logical heresy is to combine Belle II and LHCb data in a single analysis and jointly handle the systematics. In lieu of this ideal scenario, which is understandably fraught with political difficulties, an unbiased outcome can still be achieved without the sharing of data sets through rigorous bookkeeping. These experiments can communicate with each other to define the ρ pole parameters amongst others and provide a standard set of variations for all analyses to use. In return, each analysis can report the signed fit residual obtained for each systematic variation on all physics observables measured so that systematic covariance matrices can be properly constructed.

Conclusion
One of the many challenges in φ 2 -sensitive studies is that similarity of the final states and overlapping analysis techniques inevitably leads to significant systematic correlations amongst the physics observables ultimately constraining φ 2 . A coordinated approach leads to an appreciably overall improved precision, while at the same time eliminating bias, which in the case of B → ρρ is shown to be at the level of 1 • in and of itself. Additional care may also need to be taken if B 0 → (ρπ) 0 turns out to dominate the average whilst at the same time, the amplitude model uncertainty in B → ρρ is significant. Physics parameter correlations in data models can cross experimental lines for which minimal cooperation through the sharing of signed fit residuals on commonly defined systematic variations offers a good compromise to combined data analyses. Hopefully, this work inspires other investigations into the role of systematic correlations beyond single measurements in the combinations of other CP -violating weak phases such as φ 1 (β), φ 3 (γ) and φ s .

A Model uncertainties and correlation matrix
For the considered B → ρρ and B 0 → (ρπ) 0 amplitude models, the systematic uncertainties on the extracted physics parameters constraining φ 2 that globally account for correlations arising from the ρ pole parameters are given in table 8. For comparison, how these uncertainties scale when ignoring such correlations within and between these systems is also shown. The associated correlation matrix in the global case is spread over tables 9-11, while the uncorrelated scenarios can be trivially inferred from these.

Parameter
Model  Table 11. Correlation matrix (3) under the Proposed practice scheme for each determined parameter entering the φ 2 constraint. Apart from the final row entry, square brackets indicate the partial wave in the B → ρρ sector.