Systems undergoing simultaneous conservative and dissipative dynamics are prevalent across science and engineering. For example, $\partial_t\phi + u(y)\partial_x\phi = D(\partial_x^2 + \partial_y^2)\phi$ describes passive scalar transport in a laminar shear flow, and consists of the unitary advection dynamics combined with dissipative diffusion dynamics. Any linear system of differential equations $d\vec{\phi}/dt = M\phi$ can be separated into a Hermitian (dissipative) component $H_1$ and anti-Hermitian (conservative) component $iH_2$ such that $M=H_1 + iH_2$, by defining $H_1 = (M+M^\dagger)/2$ and $iH_2 = (M-M^\dagger)/2$. Quantum computers are naturally suited to simulating conservative, unitary dynamics given by $e^{iH_2 t}$, and the simulation of purely dissipative dynamics given by $e^{H_1t}$ has been well-studied in the context of imaginary time evolution algorithms~[1]. However, $e^{H_1t + iH_2t} \ne e^{H_1 t}e^{iH_2t}$, but it can be approximated to order 2 using operator splitting by $e^{(H_1 + iH_2)t} = \left(e^{iH_2\Delta t/2}e^{H_1 \Delta t}e^{iH_2\Delta t/2}\right)^{t/\Delta t} + O(\Delta t^2)$, where the overall evolution time $t$ is separated into $t/\Delta t$ steps. Improving accuracy beyond order 2 is non-trivial, since negative coefficients necessarily appear in the evolution. This is acceptable for the implementation of $e^{iH_2\Delta t}$, as passing a negative $\Delta t$ retains a unitary evolution. However, for the implementation of $e^{H_1 \Delta t}$, passing a negative $\Delta t$ results in an amplification, which can no longer be encoded into a unitary operation and destabilises the overall evolution. We overcome this by considering splitting methods with complex coefficients. A scheme that satisfies our constraints was derived by Castella et al.~[2], consisting of real coefficients for the unitary evolution ($a=1/4$), and complex coefficients for the dissipative evolution with a positive real part ($b_1 = 1/10-i/30$, $b_2 = 4/15+2i/15$ and $b_3 = 4/15-i/5$), which retains non-amplifying steps. The operator is then split as \begin{equation} e^{(H_1 + iH_2)t} \approx \left(e^{H_1b_1\Delta t}e^{iH_2a\Delta t}e^{H_1b_2\Delta t}e^{iH_2a\Delta t}e^{H_1b_3\Delta t}e^{iH_2a\Delta t}e^{H_1b_2\Delta t}e^{iH_2a\Delta t}e^{H_1b_1\Delta t}\right)^{t/\Delta t} + O(\Delta t^4) \nonumber \end{equation} The method applies to arbitrarily higher orders of accuracy that meet the required conditions, such as