To incorporate Lindblad decay dynamics into the PLDM method, we develop the following approach, which we refer to as the L-PLDM method, with a similar procedure documented in the work of L-MFE (Lindblad mean-field Ehrenfest) in Ref. 47. The current approach follows the original idea used in 47, of matching the exact time evolution of the density matrix governed by the Lindblad dissipator, with a deterministic change of the magnitude of the electronic expansion coefficients (in terms of mapping variables), but stochastic changes of the phases for these coefficients. The primary difference between MFE and PLDM, however, is that the PLDM method effectively uses the outer product of two separate coefficients (for forward and backward mapping variables, see Eq. 9) to calculate the estimator of the density matrix instead of a single set of coefficients in Ehrenfest dynamics. The use of two separate coefficients (associated with the forward and backward paths of mapping DOFs) allows for the possibility of populations being complex instead of purely non-negative and real. The L-MFE approach does not work correctly for complex populations, and thus a modified method for PLDM must be used.