The direct quadrature method of moments (DQMOM) can be employed to close population balance equations (PBEs) governing a wide class of multivariate number density functions (NDFs) Such equations occur over a vast range of scientific applications, including aerosol science. kinetic theory, multiphase flows, turbulence modeling, and control theory, to name just a few. As the name implies, DQMOM uses quadrature weights and abscissas to approximate the moments of the NDF, and the number of quadrature nodes determines the accuracy of the closure. For nondegenerate univariate cases (i.e., a sufficiently smooth NDF), the N weights and N abscissas are uniquely determined by the first 2N non-negative integer moments of the NDF. Moreover, an efficient product-difference algorithm exists to compute the weights and abscissas from the moments. In contrast, for a d-dimensional NDF, a total of (l + d)N multivariate moments are required to determine the weights and abscissas, and poor choices for the moment set can lead to nonunique abscissas and even negative weights. In this work, it is demonstrated that optimal moment sets exist for multivartate DQMOM when N = n(d) quadrature nodes are employed to represent a d-dimensional NDF with n = 1-3 and d = 1-3. Moreover, this choice is independent of the source terms in the PBE governing the time evolution of the NDF. A multivariate Fokker-Planck equation is used to illustrate the numerical properties of the method for d = 3 with n = 2 and 3.