There is no solver for that exact equation, you will need to force it...
sesolve could be used with:
noise_array = np.random.randn(N_dts) # **0.5?
H = H0 - 0.5j * L.dag() @ L +1j * L * Coefficient(lambda t, noise: noise, args={"noise": noise_array[0])
solver = SESolver(H, options={ "method": "vern7", "interpolate": False}) # Default "adams" method remember data between steps.
solver.start(psi_0, t_0)
for i, t in enumerate(np.linspace(t_0, t_final, dt)):
state = solver.step(t + dt, args={"noise": noise_array[i]}) # evolve from t to t+dt and return the state at t+dt.
This will do the evolution with in chunks of `dt` with a different noise value for each steps.
`"interpolate": False` with verner method ensure each steps start from a clean state with the new noise value.
With just the limited info you shared, this seems a somewhat easy and promising method.