Since each of the trajectories is independent in the PLDM approach, the simulation can be trivially parallelised 51 over a total of M trajectories. This drastically reduces the cost of the computation from an order of O(N 2 )O(T 3 )O(M ) to an order O(N 2 )O(T 3 ). Now, we can further make use of the simplified Hadamard products (Eq. 28 of PaperI, Ref. 1) to vectorize the code over each of the laser perturbations, which can, in principle, remove the scaling with respect to T , and our computation can thus, in principle, be further reduced. The vectorization is done by bunching different trajectories together by collecting the state vector from "m" different trajectories as different columns of a "batch" matrix, [|Ψ⟩], of size K × m,