The equations, material parameters, and boundary conditions for the chemomechanical phase field fracture model are outlined in supplementary information. The model geometry comprises an NMC811 agglomerate secondary particle with internal randomly generated primary particles. The primary particles were idealised as isotropic and linear elastic with E p = 150 GPa and Poisson's ratio of p = 0:3, while the lithiation induced strains were applied in a transversely isotropic manner -a direct consequence of the transversely isotropic lattice constants, " a and " c , as measured by Xu et al. The lattice strains were implemented such that their direction was transformed through random angles, representing that of a realistic polycrystalline particle. The carbon binder surrounding the particle is assumed to be a linear elastic, isotropic solid of Young's modulus E c = 0:3 GPa and Poisson's ratio of c = 0:3. It is typically a nano-porous solid that provides electronic conduction and structural rigidity, while allowing transport and conduction of lithium ions; thus, we assume that charge transfer reactions are permitted at the particle-binder interface despite the presence of a fully solid computational domain. We employ a single particle model that permits the solution of the interface over-potential and current density, as well as the particle voltage. The interface kinetics take the form of a typical Butler Volmer-type equation, while a simplification is made such that the electrolyte is not considered. This is suitable for calculations that involve low cycling rates, and therefore appropriate here due to the C/50 charging regime employed. Further information can be found on the single particle model in the review of Planella et al. A phase field formulation was used to describe fracture due to anisotropic lithiation-driven strains. Griffith's thermodynamics-based theory of fracture states that variation of total energy in a solid due to an incremental increase in the crack area is a function of the elastic strain energy density and the energy required to create to new surfaces. Here, Griffith's theory was recast in variational form, such that the unknown crack surface is calculated via a phase field variable, , which assumes a value of 1 for a fully cracked material point and 0 when the material is uncracked. Furthermore, a damage degradation function g can be defined that enables the diminishing stiffness (E H = gE) of the material to be quantified, where g = (1 ) 2 , and E and E' are the nominal and effective Young's modulus, respectively. A phase field length scale, l, which characterises the size of the fracture process zone was introduced and provided a relation with the material strength according to the AT2 model. In addition, to account for the different fracture energies of grain boundaries and bulk material, we introduce an auxiliary variable to characterise the transition between material properties (here fracture toughness) of the primary particles (when = 1) and that of the bond between the primary particles ( = 0). An additional auxiliary differential equation is solved, giving the distribution of , while the function h = (1 ) 2 defines the diffuse transition between the material properties; the equations, as well as the associated boundary conditions and material properties can be found in the supplementary information.