Not yet reviewed by Pith; the record is open.

The third quantization (3rd Q) for bosons provides the exact steady-state solution of the Lindblad equation with quadratic Hamiltonians. By decomposing the interaction of the Bose Hubbard model (BHM) according to Hartree approximation, we present a self-consistent formalism for solving the open-system bosonic Lindblad equation with weak interactions in the steady state. The 3rd Q with Hartree approximation takes into account the infinite Fock space of bosons while its demand of resource scales polynomially with the system size. We examine the method by analyzing three examples of the BHM, including the uniform chain, interaction induced diode effect, and Su-Schrieffer-Heeger (SSH) Hubbard model. When compared with the simulations with capped boson numbers for small systems, the 3rd Q with Hartree approximation captures the qualitative behavior and suggests an upper bound of the steady-state value. Finite-size scaling confirms the results from the 3rd Q with Hartree approximation converge towards the thermodynamic limit. Thus, the manageable method allows us to characterize and predict large-system behavior of quantum transport in interacting bosonic systems relevant to cold-atom experiments.