  We study the numerical approximation of stochastic evolution equations with a
monotone drift driven by an infinite-dimensional Wiener process. To discretize
the equation, we combine a drift-implicit two-step BDF method for the temporal
discretization with an abstract Galerkin method for the spatial discretization.
After proving well-posedness of the BDF2-Maruyama scheme, we establish a
convergence rate of the strong error for equations under suitable Lipschitz
conditions. We illustrate our theoretical results through various numerical
experiments and compare the performance of the BDF2-Maruyama scheme to the
backward Euler--Maruyama scheme.
