Computing the Lyapunov operator φ-functions, with an application to matrix-valued exponential integrators