The NEA phosphorescence spectrum [Eq. ( )] was simulated with 1000 geometries sampled from a quantum harmonic oscillator Wigner distribution in the T1 state taking  = 0.05 eV. For each geometry, the oscillator strength between 1 T and 0 S was computed with Eq. ( ). The perturbative expansion [Eqs. ( ) and ( )] included 20 singlet and 20 triplet states. S0 is determined with density functional theory (DFT) with the B3LYP functional and the 6-31G(d,p) basis set . All other states were computed with linear-response time-dependent DFT (TDDFT) , with triplet states treated as excitations of the singlet reference. Supplementary calculations were done with TDDFT with the CAM-B3LYP and B97X functionals. The results with TD-B3LYP were systematically superior, particularly predicting the 3 n*-3 * energy splitting in MeBA and MoBA in excellent agreement with the experiments . Moreover, despite multiple attempts, we have not found the 3 n* minimum for BA and MeBA with TD-B97X, for which T1 was always 3 *. With TD-CAM-B3LYP, the T1 optimization of BA was unsuccessful, systematically overstretching the carboxyl group and yielding negative excitations.