The nature of magnetic orbitals is analyzed on a series of five Cu(d(9))...Cu(d(9)) complexes, comparing various theoretical approaches. The magnetic orbitals are usually defined from ab initio mean-field calculations. It is shown that they are practically identical for the lower and upper multiplets. Diagonalization of the density matrices obtained from accurate configuration interaction wave functions provides natural magnetic orbitals, which should be considered as the reference information. No appreciable differences between the natural orbitals of both states are observed. The natural magnetic orbitals are significantly more delocalized on the ligands than the mean-field ones. It is shown that the definition of magnetic orbitals from spin-unrestricted density functional theory (DFT) calculations is not straightforward. When carefully determined, the DFT magnetic orbitals appear to strongly overestimate the metal-ligand delocalization. The consequences on the spin density are discussed.