Recently, we reported a series of global minima whose structures consist of carbon rings decorated with heavier group 14 elements. Interestingly, these structures feature planar tetracoordinate carbons (ptCs) and result from the replacement of five or six protons (H⁺) from the cyclopentadienyl anion (C₅H₅⁻) or the pentalene dianion (C₈H₆²⁻) by three or four E²⁺ dications (E = Si–Pb), respectively. The silicon derivatives of these series are the Si₃C₅ and Si₄C₈ clusters. Here we show that ptC persists in some clusters with an equivalent number of C and Si atoms, i.e., Si₅C₅, Si₈C₈, and Si₉C₉. In all these species, the ptC is embedded in a pentagonal C₅ ring and participates in a three-center, two-electron (3c-2e) Si-ptC-Si σ-bond. Furthermore, these clusters are π-aromatic species according to chemical bonding analysis and magnetic criteria.

A second-order many-body perturbation correction to the relativistic Dirac–Hartree–Fock energy is evaluated stochastically by integrating 13-dimensional products of four-component spinors and Coulomb potentials. The integration in the real space of electron coordinates is carried out by the Monte Carlo (MC) method with the Metropolis sampling, whereas the MC integration in the imaginary-time domain is performed by the inverse-cumulative distribution function method. The computational cost to reach a given relative statistical error for spatially compact but heavy molecules is observed to be no worse than cubic and possibly quadratic with the number of electrons or basis functions. This is a vast improvement over the quintic scaling of the conventional, deterministic second-order many-body perturbation method. The algorithm is also easily and efficiently parallelized with 92% strong scalability going from 64 to 4096 processors.

Density functional theory calculations support the possibility that Li ₁₀ Si ₄ and Li ₆ Si ₅ clusters can be combined to assemble structures, opening opportunities for silicon–lithium cluster assembled materials.