was determined independently inside and outside of the phantom (the latter being referred to as world ). The RANMAR random number generator was chosen, as it provides 9 × 10 8 disjoint sequences with a length of 10 30 numbers [40]. The magnetic field was implemented using a corrected version [41] of the Geant4 class HadrontherapyMagneticField3D given in the HADRONTHERAPY application example of Geant4. As a verification for the thus defined beam model, depth-dose curves were scored in a water phantom (size x$_{wp}$ × y$_{wp}$ × z$_{wp}$ = 35 × 12 × 12 cm $^{3}$). For 24 energies, the range R$_{0}$ and energy spread σ$_{E}$ were obtained by fitting the analytical Bragg curve [31] to these depth-dose curves (as in Appendix B.1). These were then compared to those obtained from the experimental data. Furthermore, two-dimensional beam profiles were measured at the phantom entrance position behind the collimators using the pixelated scintillation detector. Profiles were scored in the simulation at the same position. Spatial spreads at the phantom entrance position σ$_{ph}$ yz were extracted by fitting a bivariate Gaussian function to the obtained experimental and simulated planar dose profiles. The physics model QGSP BERT HP EMZ showed the smallest deviation to the experimental data and was thus chosen for the simulation, although differences between the studied models in range and energy spread were smaller than 0 . 2 mm and 0 . 05 MeV, respectively. The minimum required production threshold for the world was found to be 5 mm and the minimum required production threshold for the phantom, step limiter and scorer resolution were found to be 0 . 1 mm. This corresponds well with values found in previous studies [36, 42]. Mean differences between the model and reference measurements of depth-dose curves and lateral beam profiles in water without magnetic field (see Appendix B.2), were 0 . 2 mm in range (max. 0 . 5 mm), 0 . 23 MeV in energy spread (max. 0 . 35 MeV) and 0 . 1 mm in spatial spread (max. 0 . 2 mm). The beam model was thus accepted for the simulations. Dose was scored in the PMMA phantom at the position of the film detector, on a 20 × 15 cm 2 plane with a thickness of 30 µ m corresponding to that of the active layer of a Gafchromic EBT3 film dosimeter (Ashland, Covington, USA). We note that the parameters defining the spatial and energy distributions were directly implemented as measured, and not adjusted iteratively to improve the match between simulated and measured data. This allowed an unbiased comparison of the different physics lists. For the deflection study (see main manuscript), the convergence study for the simulation resolution parameters was repeated with the range and the Bragg peak deflection as observables and a required accuracy 19 Appendix B.2. Beam model verification