We study chaoticity and thermalization in Bose-Einstein condensates in disordered lattices, described by the discrete nonlinear Schrödinger equation (DNLS). A symplectic integration method allows us to accurately obtain both the full phase space trajectories and their maximum Lyapunov exponents (mLEs), which characterize their chaoticity. We find that disorder destroys ergodicity by breaking up phase space into subsystems that are effectively disjoint on experimentally relevant timescales, even though energetically, classical localisation cannot occur. This leads us to conclude that the mLE is a very poor ergodicity indicator, since it is not sensitive to the trajectory being confined to a subregion of phase space. The eventual thermalization of a BEC in a disordered lattice cannot be predicted based only on the chaoticity of its phase space trajectory.Introduction In this Letter, we bring together the topics of disorder, nonlinearity, chaos, ergodicity, and Bose-Einstein condensation (BEC) in optical lattices. Inspired by the realisation of disordered optical potentials in experiments with ultracold atomic gases [1, 2], we explore the question of thermalisation in such systems. Previous theoretical works have, amongst others, addressed the topics of Bose and Anderson glasses [3], Anderson localisation [4,5], and Lifshits glasses [6]. In [6], various regimes of interaction strengths were investigated, and it was found that for sufficiently strong interactions, a disordered BEC is expected. We will focus on this regime, and study the Bose-Hubbard model [7,8], which describes a bosonic gas in a lattice, in the mean-field approximation. The resulting model can be obtained by discretising the Gross-Pitaevskii equation, and is also known as the discrete nonlinear Schrödinger equation (DNLS) [9]. In this Letter, we pose and answer the following question: what is the effect of disorder on thermalisation and ergodicity in the mean-field Bose-Hubbard model / DNLS?The connection between chaoticity and thermalisation was recently discussed in the disorder-free case [10]. Intuitively, one would expect chaotic trajectories to thermalise, since unlike regular ones, they are not confined to the neighborhood of stable periodic orbits. Not being confined, the expectation is that they are able to cover the available phase space, and that the system is therefore well-described by the microcanonical ensemble, since only the energy and particle number are conserved. We show that this expectation is not correct, by explicitly demonstrating the absence of equiprobability of states on the microcanonical shell. Since the energies involved are higher than the disorder potential, classical localisation cannot be responsible for this effect.The most commonly employed method of chaos detection, which quantifies the sensitive dependence on initial conditions, is the evaluation of the maximum Lyapunov exponent (mLE) [11][12][13][14]. In [10], a positive mLE, which