Kohn-Sham DFT Calculations: Most of the bonding analyses were performed with scalar relativistic DFT calculations using functionals with varying fractions of exact exchange (eX). DFT calculations were performed in spin-unrestricted Kohn-Sham (UKS) fashion with the DKH2 Hamiltonian and the same basis sets as used in the MC-WFT calculations. Natural Bond Orbital (NBO) and Natural Localized Molecular Orbitals (NLMO) calculations and Natural Population analysis (NPA) were carried out with the NBO program version 6. The extent of LM donation reflecting the Pu 5f occupation in PuO 2 was examined by varying the percentage of eX between 0 and 100% in a hybrid based on the Perdew-Burke-Ernzerhof (PBE) 105 exchange-correlation functional. Note that a global hybrid functional combining PBE with 25% eX corresponds to the widely used PBE0 parametrization ('PBE1PBE' in the input syntax of some codes). The NPA identified considerable excess population of the Pu 5f shell as well as '6f' populations. We consider the '6f' occupations identified by the NPA as arising from the use of a large atom-centered basis set and interpret them as donation into diffuse regions of the 5f shell, labeled as 5f'. Considerable '6f' occupations from NPA were also noted in previous studies of actinide bonding and interpreted in a similar way as done here. Additional DFT calculations were performed for the PuO 2 cluster model and the aforementioned aqua ions to obtain the density of states (DOS). The latter was generated using Multiwfn. Additional DFT calculations with periodic boundary conditions (PBCs) were performed using a revised PBE functional for solids (PBEsol) 108 with the on-site Coulomb repulsion in the 5f shells treated via a Hubbard U correction 109 (DFT+U), using the Vienna ab-initio Simulation Package (VASP). Complete details about these periodic DFT calculations are given in the Supporting Information (SI), Section S9.