In variational phase-field modeling of brittle fracture, the functional to be minimized is not convex, so that the necessary stationarity conditions of the functional may admit multiple solutions. The solution obtained in an actual computation is typically one out of several local minimizers. Evidence of multiple solutions induced by small perturbations of numerical or physical parameters was occasionally recorded but not explicitly investigated in the literature. In this work, we focus on this issue and advocate a paradigm shift, away from the search for one particular solution towards the simultaneous description of all possible solutions (local minimizers), along with the probabilities of their occurrence. Inspired by recent approaches advocating measure-valued solutions (Young measures as well as their generalization to statistical solutions) and their numerical approximations in fluid mechanics, we propose the stochastic relaxation of the variational brittle fracture problem through random perturbations of the functional. We introduce the concept of stochastic solution, with the main advantage that point-to-point correlations of the crack phase fields in the underlying domain can be captured. These stochastic solutions are represented by random fields or random variables with values in the classical deterministic solution spaces. In the numerical experiments, we use a simple Monte Carlo approach to compute approximations to such stochastic solutions. The final result of the computation is not a single crack pattern, but rather several possible crack patterns and their probabilities. The stochastic solution framework using evolving random fields allows additionally the interesting possibility of conditioning the probabilities of further crack paths on intermediate crack patterns.