  We develop a novel numerical scheme for the simulation of dissipative quantum
dynamics following from two-body Lindblad master equations. All defining
continuum properties of the Lindblad dynamics, hermiticity, positivity and in
particular trace conservation of the evolved density matrix are preserved. The
central ingredient is a new spatial difference operator, which not only fulfils
the summation by parts (SBP) property but also implements a continuum
reparametrization property. Using the time evolution of a heavy-quark
anti-quark bound state in a hot thermal medium as an explicit example, we show
how the reparametrization neutral summation-by-parts (RN-SBP) operator
preserves the continuum properties of the theory.
