Pathway to Polyradicals: A Planar and Fully π-Conjugated Organic Tetraradical(oid)

In this work, we provide a general strategy to stabilize the ground state of polyradical(oid)s and make higher spin states thermally accessible. As a proof of concept, we propose to merge two planar fully π-conjugated diradical(oid)s to obtain a planar and cross-conjugated tetraradical(oid). Using multireference quantum chemistry methods, we show that the designed tetraradical(oid) is stabilized by aromaticity and delozalization in the π-system and has six thermally accessible spin states within 1.72 kcal/mol. Analysis of the electronic structure of these six states of the tetraradical(oid) shows that its frontier π-system consists of two weakly interacting subsystems: aromatic cycles and four unpaired electrons. Conjugation between unpaired electrons, which favors closed-shell structures, is mitigated by delocalization and the aromaticity of the bridging groups, leading to the synergistic cross-coupling between two diradical(oid) subunits to stabilize the tetraradical(oid) electronic structure.

Since tetraradicals can be seen at least conceptually as two diradicals, in the present case the diradical(oid) 2,2 ′ -(5,11-dihydroindolo [3,2-b]carbazole-3,9-diyl)dimalononitrile S is conceptually merged through the central benzene ring with the diradical(oid) 2,2 ′ -((1,4-phenylenebis(ethyne-2,1-diyl))bis(4,1-phenylene))dimalononitrile L, resulting in the potential tetraradical(oid) 2,2 ′ -(6,12-bis((4-(dicyanomethyl)phenyl)ethynyl)-5,11-dihydroindolo[3,2-b]carbazole-3,9-diyl)dimalononitrile LS .Thus, we refer the tetraradical(oid) described in this work as LS .The notations can be understood from the Figure S1 as follows: L represents longer chain of the molecule including the central benzenoid ring and two benzenoid ring connected to the central benzenoid ring, delimited by ethynyl residue, while S represents the shorter perpendicular chain in the molecule, composed of the indolocarbazole with dicyanomethylene substituents at the ends.Therefore, this notation will be used to help describe electronic structure and density distributions.Also, this notation can be adapted to describe 4 different singly occupied natural orbitals symbolically.There are five resonance structures given in Figure S1.In the top row, on the left, short (S ) chain is quinoidal and long chain (L) diradical, this resonance structure will be referred as LDSQ.In the , both systems short (S ) and long (L) chain are diradical in the combined system LS , thus, this resonance structure will be referred as LDSD.On the right, the short system (S ) is diradical and long system (L) is quinoidal, thus, this resonance structure will be referred as LQSD.In the bottom row resonance structures, both subunits L and S are radical, hence, we refer structure on the left as LRSR and structure on the right as LRSR ′ .In the main text of the manuscript, LRSR and LRSR ′ resonance structures are omitted for conciseness.Throughout this work, all orbitals are visualized as isosurfaces with value of 0.015 unless explicitly stated otherwise.In the top row, leftmost resonance structure is referred to as LDSQ, middle resonance structure as LDSD, and rightmost resonance structure as LQSD.In the bottom row, structure on the left will be referred to as LRSR and structure on the right as LRSR ′ .Rings with Clar's aromatic π-sextet are filled-in blue.

S2.1. Geometry optimization
For geometry optimizations of the molecules, Density Functional Theory (DFT) was used with Amsterdam Density Functional (ADF) module engine 1 within Amsterdam Modeling Suite (AMS) release of 2021 (ADF 2021.104version).Since the molecule presented in this manuscript had been determined to possess significant tetraradical character in the preliminary tests, the unrestricted Kohn-Sham (UKS) DFT formalism with quintet state (keyword line "SpinPolarization 4" in ADF engine) was used for the geometry optimization with Generalized Gradient Approximation (GGA) exchange-correlation functional BLYP 2,3 without any relativistic corrections due to absence of heavier atoms than the first row (H, C and N).The Slater-type all-electron triple zeta (TZP) 4 basis set available in ADF was used.The self-consistent field (SCF) convergence criteria were chosen as follows: the norm of the commutator between Fock and Density matrices ([F,P]) set to 1.0 • 10 −5 and maximum element of the [F,P] set to 1.0 • 10 −6 (defaults).Numerical quality for Becke's integration grid 5,6 was set to "VeryGood", which means Lebedev angular grid order varying in range [13, 41] and number of radial points varying in range [64,  72].Also, the numerical quality for density fitting approximation was set to "VeryGood", which means number of interpolation points varying in range [120,160] and maximum number of L-expansion varying in range [9,10].All the other technical parameters were left as default for ADF software package.In order to verify that the optimized geometry corresponded to the minimum on the Potential Energy Surface (PES), a Hessian with respect to nuclear coordinates was computed which turned out to have all eigenvalues positive, translating into all real calculated infrared absorption frequencies corresponding to the determined normal modes. 7It must be mentioned that the planar structure of the molecule LS results irrespective of the planarity of the starting geometry.For example, if in the starting geometry the angle between planes defining the indolo [3,2-b]carbazole backbone and the plane defined by the 2-(4-dicyanomethylene-phenyl)-ethynyl substituent moieties is about 45°and 30°(for the first and the second moiety, respectively), the optimized geometry is still planar.Furthermore, even if these angles are 65°and 50°, the resulting optimized geometry is still planar.Thus, the geometry used throughout the electronic-structure calculations is an optimized structure obtained from starting planar geometry belonging to the C 2h point group in order to take advantage of the molecular symmetry in the downstream more expensive ab-initio calculations (from Figure S1, the assignation of the point group is clear for the LS system).Finally, we remark that all computations were carried out for a single, isolated molecule in the gas phase.

S2.2. Remarks on the methodology
Investigation of the electronic structure has been done with wavefunction methods rather than DFT because, as it will be shown in Figure S7 in the Section S4, DFT shows heavy functional dependence on energy spectrum of the allowed spin states.Hartree-Fock (HF) with its restricted and unrestricted formalism (RHF and UHF, respectively) is a computationally cheap but definitive method to determine whether or not the electronic structure of the molecule in the ground state is open-shell.To achieve this, all that is required is to compute RHF and UHF solutions separately and compare results.If UHF solution has lower energy than RHF solution (it can have the same energy or lower only), it means that closed-shell electronic wavefunction is unstable and the ground-state wavefunction is better described by a UHF, an open-shell, solution.Nonetheless, monoconfigurational HF and DFT methods cannot give us an appropriate description of different spin states of such a multiconfigurational system as tetraradical(oid).
For an appropriate description of an electronic structure of tetraradical(oid), a qualitatively correct wave function can be obtained by multiconfigurational SCF (MCSCF) methods such as Complete Active Space SCF (CASSCF), also known as Full Optimized Reaction Space (FORS), [8][9][10][11][12][13][14][15] which describes the so-called static/nondynamic correlation.HF, CASSCF and CASCI (explained below) calculations were performed using The General Atomic and Molecular Electronic Structure System (GAMESS) version 2016 (R1). 16For UHF and UKS DFT benchmark given in Section S4 as well as for CASSCF and CASCI calculations, Dunning's correlation-consistent double zeta basis set cc-pVDZ was used. 17CAS Configuration Interaction (CASCI) is the single-variable mode of CASSCF, in which molecular orbitals coefficients are not being optimized but only CI coefficients.DFT and UHF benchmark given in Section S4 was done with GAMESS version 2023 (R2) 18 due to availability of building custom exchange-correlation functionals.Improved descriptions including dynamic electron correlation require expensive multi-reference perturbation (MRPT) or multi-reference configuration-interaction (MRCI) calculations, which considering the size of the system -CASSCF and CASCI with (16,16) active space -cannot be computed.Even for (4,4) minimal active space (considering that, by definition, only 4 unpaired electrons can be described) 681 external molecular orbitals that remain, considering the basis set and active space, are too many to include.Only including up to 30 or 40 external orbitals captures only negligible part of the correlation energy.This was checked by the corresponding Multireference Moller-Plesset second order Perturbation Theory (MRMP2) [19][20][21][22][23][24] calculations, which showed by more than 100 factors of correlation energy correction than what MRCI showed with only including 30-40 external orbitals.Therefore, the only computationally practical method for studying the system of the size that is presented in this work is MRMP2.Nonetheless, since the results from MRMP2 gave same qualitative picture as CASCI and CASSCF calculations, they are not included nor discussed in this work due to redundancy.S2.3.Check of the qualitative electronic structure of the ground state with HF At first, RHF and UHF calculations were done for a singlet state to determine whether or not the ground state electronic structure of this potential tetraradical(oid) is open-shell (justified by the resonance structures of the compound given in Figure S1).It turned out that the UHF broken-symmetry singlet has much lower energy (by 189.6 kcal/mol) than the RHF solution.Also, occupation numbers of Natural Orbitals (NOs) indicate that there are four Singly Occupied NOs (SONOs), which implies four unpaired electrons, meaning that molecule is a tetraradical(oid).NOs are eigenvectors of the first-order density matrix operator, and a corresponding eigenvalue to the particular NO is its occupation number 25 (denoted as n N O ).The Highest Occupied NO (HONO) is defined as the orbital which has the lowest occupation number n N O among NOs with n N O ≥ 1 .Lowest Unoccupied NO (LUNO) is defined as the orbital which has the highest n N O among NOs with n N O ≤ 1.In UHF, n HON O−i + n LU N O+i = 2.000, which is also manifested usually in CASSCF and CASCI natural orbitals.Since it is important to study the molecule as having potentially 4 persistent unpaired electrons, broken-symmetry solution was not simply obtained by mixing Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO), but with the following procedure.Firstly, UHF quintet calculation was run, then UHF NOs were taken as guess for Restricted Open-Shell Hartree-Fock (ROHF) quintet calculation.Then Pipek-Mezey localization 26 was done in order to localize converged ROHF open-shell/singly occupied orbitals.Having checked the localization, it showed that each Singly Occupied Molecular Orbital (SOMO) was localized on each particular spin center (see Figure S2).After localization, the orbitals were ordered in order to obtain a singlet of the type |α(1)α( 2)β(3)β(4)⟩, where orbitals are assigned as follows from the Figure S2: SOMO 1 ≡ 1, SOMO 2 ≡ 2, SOMO 3 ≡ 3, SOMO 4 ≡ 4. Thus, both subsystems S and L are initiated as having opposite spins for their respective two unpaired electrons.Resulting UHF singlet NOs with their occupation numbers are given in Figure S3.This showed that the molecular electronic structure of the LS is tetraradical(oid) with spin-opposite unpaired electrons for each substructure S and L. NOs occupation number complementarity and still not exactly being 1.000 (but very close to it) shows that the LS is indeed a tetraradical(oid) with very high tetraradical character.Although HF is sufficient to qualitatively gauge the ground state character of the electronic structure of LS , UHF states are not pure.To describe the electronic structure and distribution of states of the possible tetraradical(oid) qualitatively correctly, we need to employ multiconfigurational calculations such as CASSCF and CASCI.

S2.4. CASSCF and CASCI calculation details
Having determined that LS is a potential tetraradical(oid), UHF quintet calculation was run to obtain UHF NOs, which includes 4 SONOs.Therefore, to study a potential tetraradical(oid), one of the best guess orbitals are UHF quintet NOs.The NOs were selected in the active space based on their NO occupation numbers with reasonable cutoff so as to include all the degenerate relevant orbitals and maintain the electronic problem tractable. 27ence, based on these limits, the largest active space that was used for MCSCF calculations was 16 orbitals with 16 electrons - (16,16).Note that multiconfigurational wavefunction can be obtained for symmetry-adapted (if guess orbitals are not localized) or C 1 point group states (if guess orbitals are localized).With S z = 0 determinant basis, (4,4) and (16,16) CASSCF calculations were run for 8 states from the UHF quintet NOs guess centered around HONO and LUNO, thus (4,4) representing minimal model space for 4 unpaired electrons (thus tetraradical(oid)).Also, to demonstrate the persistence of tetraradical spin states, the CASSCF(10,10) and CASSCF (14,14) calculations were run with UHF triplet NOs guess.We found with every size of active space that there are six states that are very close to each other in the ranges from which 603 cm −1 is a maximum.
After determining 6 spin states of tetraradical(oid) LS , and seeing that energies of these two singlets, three triplets and one quintet are very close, converged CASSCF (16,16) quintet optimized orbitals were taken, and (4,4) model frontier subspace (those that represent 4 SOMOs) was localized within (16,16) active space and used as a guess for CASCI (4,4) and CASCI (16,16) calculations.Having inspected that localization resulted each SOMOs being assigned to particular center as shown in Figure S4, we can also analyze the determinantal build of the corresponding wavefunctions of the spin states determined by CASCI.More than six states were determined with S z = 0 determinant basis.For (4,4) model space, these determinants are |α(1)α( 2 3)β(4)⟩.It turned out from CASCI that the ordering of states is the following: S 0 , T 0 , T 1 , S 1 , T 2 , Q 0 within range of 460.9 cm −1 .Orbital numbering is as follows from Figure S4: SOMO 1 ≡ 1, SOMO 2 ≡ 2, SOMO 3 ≡ 3, SOMO 4 ≡ 4. Thus, CI coefficients for each determinantal contribution can be translated into electron density distribution of spin centers.The existence of low-energy quasi-degenerate (magnetic scale difference) states with all the possible spin multiplicities for a tetraradical electronic structure shows that our molecule has 4 unpaired electrons.Thus, LS is a tetraradical(oid).

S2.5. Broken-symmetry calculations with HF and DFT
For broken-symmetry calculations with HF and DFT, converged CASSCF (16,16) quintet orbitals were taken and (4,4) frontier subspace localized as shown in Figure S4.These orbitals are assigned as for previous CASCI calculations in Section S2.4.Then α and β sets of localized guess orbitals were reordered to model the determinants for each spin states such that S = S z .Since, in principle, |αααα⟩ is equivalent to |ββββ⟩ because both describe all-spin-parallel determinants, this equivalence and molecular symmetry have been used to define unique determinants for each S = S z subspace.Unique singlet determinants are |α(1)α( 2  For all the calculation results, the 8 types of orbitals are obtained for frontier singly occupied molecular (or natural) orbitals.If we assign symbols, we can categorize each spin state with these orbitals, their symmetry and

S3.1. CASSCF and CASCI results
Table S4.CASSCF (16,16) states determined by state-specific calculations.For all the calculations, UHF quintet NOs were used as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to Figure S5 and NO occupation number.

State Symmetry
Energy (a.u.) CASCI NOs occupancy, identity and symmetry Overall electron density corresponds to the same as if they were from Table S6.CASCI (16,16) from CASSCF (16,16) quintet converged orbitals as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to Figure S5 and NO occupation number.

State Symmetry
Energy (a.u.) CASCI NOs occupancy, identity and symmetry.∆E from G. S. (cm −1 ) HONO -1 HONO LUNO LUNO + 1 Table S7.CASCI(4,4) determinant contributions for 6 near-degenerate tetraradical(oid) spin states.These states are obtained for Sz = 0 basis from CASSCF (16,16) converged quintet localized (in (4,4) SOMOs subspace) initial guess orbitals.their NO occupation numbers.These assignations are given in Figure S5.CASSCF(4,4) calculations results are given in Table S1, CASSCF(10,10) calculations results are given in Table S2, CASSCF (14,14) calculations results in Table S3, and CASSCF (16,16) calculations results in Table S4.As we see from these results, we have 6 states which are very close to each other in energy.Also, the properties of corresponding spin states are very similar in calculations across different sizes of active space.The maximum range is, as mentioned before, 603 cm −1 .This means that LS is a magnetic molecule as the transition energies between its different spin states are very small.To better understand the properties of the molecule, we ran CASCI(4,4) and CASCI (16,16) from CASSCF(16,16) quintet converged orbitals with 4 SOMOs localized.CI coefficients of each determinantal contribution can be translated into electron density distribution of spin centers, because of localized guess.Results for CASCI (4,4) are given in Table S5 and results for CASCI (16,16) are given in Table S6.To gain further insight into the electronic structure of these spin states of the tetraradical(oid), we can analyze determinantal composition of wavefunction for each spin state.This kind of analysis for CASCI (4,4) calculation is given in Table S7 and for CASCI (16,16)  calculation is given in Table S8.

S3.2. Results from Heisenberg-Dirac-van Vleck model and effective Hamiltonian theory
Heisenberg model describes open-shell systems as particle-per-site model of spin centers. 28The general structure of this model is well known in the chemistry and (especially) molecular or solid-state physics community.The model spin Hamiltonian is called Heisenberg-Dirac-van Vleck (HDVV) Hamiltonian [28][29][30] and corresponding energy expressions of its eigenvalues are expressed in terms of two-body magnetic coupling constants J 12 , J 13 , J 14 , J 23 , J 24 , J 34 .We can make a matrix representation of the HDVV hamiltonian as follows from Eq. S1: 31 ĤHDV For a 4-site model, this Hamiltonian can be expanded as follows from Eq. S2: Where J ij is the exchange-coupling constant between Ŝ i and Ŝ j localized spin moments.According to the adopted definition of this model Hamiltonian in Eq.S1, the positive value of the exchange-coupling constant corresponds to an antiferromagnetic interaction between Ŝ i and Ŝ j magnetic moments, while negative value corresponds to the ferromagnetic interaction between these magnetic moments i and j (anti-parallel and parallel spin alignments, respectively).
In order to build up this Hamiltonian matrix in the basis that is appropriate, we can choose S z = 0 subspace for which all spin states that a tetraradical(oid) could have (singlet, triplet and quintet) can be represented with expansion of 6 determinants, which are all neutral determinants in (4,4) active space with S z = 0. Matrix elements can be found by applying the following relation expressed in Eq.S3 31 as follows: Where Ŝ i+ and Ŝ j− are ladder operators applied to particular spin center i and j, respectively.Therefore, by application of these operators, we can find matrix elements of the HDVV Hamiltonian in the basis of the S z = 0 determinants.The basis is ordered as follows: |ααββ⟩, |ββαα⟩, |αβαβ⟩, |βαβα⟩, |αββα⟩, |βααβ⟩.For example, to find the part of element ĤHDV V (3,3) for J 12 term, we can apply this expansion as shown in Eq.S4: Because 6 determinants of the S z = 0 form an orthonormal set.By application of these operators, we obtain the following HDVV Hamiltonian matrix as given in Eq.S5.
If we note that the molecule has C 2h point group symmetry and refer to the notation of SOMOs in Figure S4, we can note that interactions between magnetic centers 1 and 2 is equivalent to those between 3 and 4, and interactions between magnetic centers 1 and 3 is equivalent to the interaction between centers 2 and 4 (due to C 2 rotational symmetry).Therefore, we create the following notation for further simplification: J 12 = J 34 = J a , J 13 = J 24 = J b .In addition, considering the previous notation of S and L, we create the following notations: J 14 = J S and J 23 = J L .Hence, the HDVV Hamiltonian is simplified as given in Eq.S6: We can diagonalize this Hamiltonian to find its eigenvalues and eigenvectors.This results in the following eigenvalues as given in Eq.S7 and Eq.S8, with unnormalized eigenvectors in Eq.S9. Where In order to demonstrate how HDVV Hamiltonian could be used as a correspondent to an effective Hamiltonian ( Ĥeff ), if this system is to be appropriately modeled by Heisenberg Hamiltonian, we need to assemble Ĥeff for the problem and make one-to-one correspondence between ĤHDV V and Ĥeff .In this case, an appropriate matrix representation of the effective spin Hamiltonian is constructed from CASCI (4,4) wave functions because, as we see in Table S7, more than 99% of the wavefunction on this level of theory is composed of complete set of (4,4) neutral determinants that belong to S z = 0 subspace.Hence, effective Hamiltonian theory can be used as a rigorous mathematical scheme to reduce the electronic Hamiltonian to a spin Hamiltonian of the HDVV type by means of projection techniques. 32][35] The procedure of this approach starts by selecting the spin space (which is the model space) determining the low-lying energy spectrum.To construct the effective Hamiltonian, we need to project the approximate solutions of the exact nonrelativstic electronic Hamiltonian (in our case, (4,4) and (16,16) CASCI solutions) onto the N -dimensional model space S that contains the information of interest. 36This formalism ensures one-toone correspondence between the eigenvalues of the exact Hamiltonian and those of the effective Hamiltonian, and coincidence of the eigenfunctions of Ĥeff with the projection of the approximate solutions of the exact Hamiltonian onto the model space as closely as possible.Let us denote E m and ϕ m as energy and the corresponding eigenfunction of state m.It follows that: We can further define a projector that targets the model space S as Where {|I⟩} is an orthonormal basis of this model space, which in our case is 6 S z = 0 determinants |ααββ⟩, |ββαα⟩, |αβαβ⟩, |βαβα⟩, |αββα⟩, |βααβ⟩ given in Table S7 which span ĤHDV V shown in Eq.S5.In application, this requires previous localization step (note that our (4,4) and (16,16) CASCI were run with localized 4 frontier SOMOs).Hence, an effective Hamiltonian can be defined as given in Eq.S12: where |ϕ S,m ⟩ = PS |ϕ m ⟩.Since the basis in which ϕ S,m is expressed is not necessarily orthonormal, we use orthonormalized projections as proposed by des Cloizeaux given in Eq.S13 or orthonormalize this set using Gram-Schmidt procedure.
where S m,q are elements of the overlap matrix between states m and q as shown in Eq.S14: Hence, finally, an effective Hamiltonian is obtained from its spectral decomposition, and its matrix representation in the model space spanned by, in our case, 6 S z = 0 neutral determinants as given in Eq.S15.
Hence, the elements of Ĥeff can be directly compared to the elements of ĤHDV V , which in our case is given Eq.S5.By application of this procedure, we have computed the effective Hamiltonian for (4,4) CASCI results given in Eq.S16 and (16,16) CASCI results projected onto the (4,4) model space given in Eq.S17.
From comparing Eq.S6, Eq.S16 and Eq.S17, we see that it is trivial to determine value of each exchangecoupling constant from ĤHDV V .We see that largest exchange-coupling constant value is J S = 312.06cm −1 and the second largest is J a = 59.92 cm −1 , followed by J L = 43.34cm −1 and J b = 12.92 cm −1 .This would mean that both subsystems L and S have their electrons anti-ferromagnetically (spin anti-parallel) coupled and in turn this would mean that we would have a singlet ground state, which is indeed a case.Ovchinnikov's rule 37 predicts a singlet ground state for L, S , and LS .The qualitative picture is similar for the effective Hamiltonian obtained from (4,4) and (16,16) CASCI calculations.As we see, the Heisenberg model gives us correct qualitative picture of what all CASSCF and CASCI calculations showed.Nonetheless, the quantitative picture of predicted by the Heisenberg model should be analyzed thoroughly.It is very important to point out that collective weight of the neutral determinants that span (4,4) S z = 0 subspace for (16,16) CASCI is about 72% for all 6 magnetic spin states given in Table S8.This means that other determinants that lie out of (4,4) S z = 0 and neutral determinants subspace have significant contribution (about 28%) into the wavefunction obtained from (16,16) CASCI calculation.Hence, even though qualitatively correct ground state is predicted, the wavefunction of this molecule cannot be fully spanned by (4,4) active space and consequent 6 × 6 model Hamiltonian ĤHDV V built in the basis of six neutral determinants belonging to (4,4) active space.It is important to remark that our investigation through effective Hamiltonian reveals exchange-coupling interaction J a between spin centers 1 and 2 (refer to Figure S4 for assignations) to be greater than between 2 and 3, which is J L .This can be partially rationalized from the resonance structures given in Figure S1.Electrons in centers 1 and 2 could pair-up antiferromagnetically and close the shell.Nonetheless, it would not be obvious to expect that J a is greater than J L .Hence, ĤHDV V gives us additional insight into the dominant interactions in this magnetic system.This significant interaction J a provides synergistic effect that keeps every open shell faithfully from diradical(oid)s L and S to create tetraradical(oid) LS .These interactions can schematically be represented as shown in Figure S6a and low-energy spin spectrum is depicted in Figure S6b.

S4. Broken-symmetry UKS DFT and UHF benchmark results
This DFT and UHF benchmarks was performed in GAMESS version 2023 (R2). 18Benchmark results for DFT PBE 38,39 functional hybrid with variable percentage of HF exchange (PBEX) are given in Table S9.In order to properly visualize the variance of the results for particular spin states with different approximate Hamiltonian, we can look at the Figure S7.It is apparent that there is quite significant variability between energies of states and their relative distribution of minimum and maximum energies depending upon the % of HF exchange in their build.Also, there are significant differences between UHF and UKS results.Thus, this is an empirical demonstration that single-reference methods like HF and DFT are not sufficient to characterize all relevant spin states of LS .It is noteworthy that some of the states calculated end up in a closed-shell solutions that does not describe the qualitative electronic structure of this molecule, which is indeed, according to all CASSCF and CASCI calculations, a tetraradical(oid).Hence, in description of LS , we can rely only on CASCI and CASSCF results as qualitatively trustworthy.Nonetheless, it must be mentioned that quintet state is always described qualitatively correctly with respect to 4 frontier SOMOs that show its tetraradical(oid) identity.Hence, quintet DFT for a tetraradical(oid) can be partially relied upon for some basic qualitative properties, but only describing 4 SOMO part of the larger relevant active space as shown in Table S7 and Table S8.Furthermore, there are some spin states for which some functional properly captured the main aspects of the electronic structure.These cases can be used for the calculation of several properties of this tetraradical(oid).

S5. Study of Aromaticity of Tetraradical(oid) LS
Anisotropy of Induced Chemical Current (ACID) 40 calculations were done in order to explore the aromaticity of this tetraradical(oid).Since the qualitative electronic structure of several spin states modeled with PBE-X is aberrant, we need to select the combination of functional build and modeled spin states that reproduces correct qualitative structure of LS .We could select PBE-X with 25% of HF exchange, which corresponds to PBE0 hybrid exchange-correlation functional. 41,42For PBE0, the singlet (|ααββ⟩) and quintet (|αααα⟩) show qualitatively correct tetraradical(oid) structure, as they manifest nearly 4 SONOs.Thus, exploring the aromaticity for these states with PBE0 is justified.
We can note that aromaticity can help explain the stability of the tetraradical(oid) electronic structure.Indeed, in the resonance structure in Figure S1, in which the molecule has 4 unpaired electrons (LDSD), all of its benzenoid rings possess 6 π electrons, making them Clar's aromatic sextets. 43Other resonance structures (LQSD, LDSQ and LRSR/LRSR ′ in Figure S1) have only 2 Clar's π-sextets, making them less aromatic cumulatively.Moreover, these resonance structures preclude "global aromaticity" in the indolocarbazole aromatic backbone.The singlet and quintet states were calculated with PBE0/cc-pVDZ level of theory in Gaussian 09 program.From those results, ACID plots were computed separately and are given in Figure S8. Figure S8.ACID plots for singlet (Figure S8a) and quintet (Figure S8b) states of LS calculated with PBE0/cc-pVDZ level of theory.Isosurface value is 0.030.
Since we see on the ACID plot that there is a continuous ring current across the indolocarbazole backbone, the compound possesses "global aromaticity" that is a further stabilizing factor than merely local aromaticity of Clar's π-sextets.Thus, in the resonance structure LDSD we have 5 Clar's π-sextets (three more than other resonance structures) and "global aromaticity" that stabilizes tetraradical(oid) electronic structure even more.From these results, we can conclude that the active space of this molecule contains two subsystems, which interact weakly but appreciably enough to affect the spin spectrum, because rest of the correlation is important for defining the spin states of a tetraradical(oid).The first of these subsystems is aromatic core(s) in the indolocarbazole and two other phenyl rings which act as an almost inner part of the active space, which appears in Table S8 as determinants of type S-L*, L-S*, S-S*, L-L*, S,S-S*,S*, L,L-L*,L*, L,S-L*,S*.These are excited determinants which contain excitations in the aromatic moieties of this molecule.This proves that the electron correlation in the aromatic part of the molecule is indeed important.The second of these subsystems is 4 SOMOs that define tetraradical(oid) with all possible spin states in low-energy spectrum within 603 cm −1 .It must be emphasized that conjugation between radical centers affects the spin spectrum of this tetraradical(oid), because the interaction/coupling between radical centers is evidently non-zero.Nonetheless, this conjugation effect between radical centers is somewhat diminished by the presence of aromaticity of phenyl rings and "global aromaticity" of indolo [3,2-b]carbazole backbone.Hence, this shows that the nature of the tetraradical(oid), the role of the aromaticity in its stabilization and the degree of interaction between these aromatic and tetraradical subsystems define spin states of this fully planar π-conjugated tetraradical(oid) molecule with multiple thermally accessible singlet, triplet and quintet spin states.
In addition to the argument from the "global aromaticity" standpoint, it is also interesting to measure local aromaticity of the rings within LS .For those, we use Multicenter Index (MCI) of aromaticity, 44,45 which was calculated for each ring given in Figure S9.For a given system, the MCI shows the simultaneous electron sharing between all centers, irrespective of the position of the atoms in the molecule. 45MCI results for each ring of singlet, triplet and quintet states are given in Table S10.Notably, MCI values are hardly affected by changing the exchange-correlation functional, even by using CAM-B3LYP.To obtain these states, we start with a quintet state calculation.For both states, we use the quintet converged guess orbitals, and flip the spin on two of the radical-bearing atoms (using "SpinFlip" keyword in ADF engine) to obtain a singlet state, while for a triplet state, we flip the spin on one of the radical-bearing atoms.All of these converged states were verified to qualitatively approximate the tetraradical(oid) electronic structure.These results show that tetraradical(oid) spin states of LS involve aromatic rings as depicted in the LDSD resonance structure in Figure S1.The values of MCI given in Table S10 are close to MCI values of the benzenoid rings in other aromatic compounds. 44,45

S6. Study of Tetraradical(oid) LSH
Since cyano groups are somewhat stabilizing factor for radical centers via π-conjugation, it is pertinent to study the molecule which has hydrogens instead of cyano groups as terminal substituents in LS .We will refer to this molecule as LSH , for which resonance structures are given in Figure S10.Geometry optimization, CASSCF(4,4), CASSCF (16,16), CASCI(4,4) and CASCI(16,16) calculations were run for LSH with exactly the same procedure and level of theory as for LS .Results of CASSCF(4,4) calculations are given in Table S11 and results of CASSCF (16,16) calculations are given in Table S12.Since LS and LSH overlap in most of their structures, it turns out it is possible to make same assignations to frontier NOs of LSH as it was for LS .Hence, we can have a common assignation of frontier NOs as given in Figure S5 and we can analyze spin states of LSH accordingly.From CASSCF(4,4) and CASSCF (16,16) results for LS and LSH , we can see that LSH has its diradical and tetraradical character close to and slightly lower than those of LS .Thus, since we have previously shown that LS is a persistent tetraradical(oid), we can state the same about LSH .This comparison shows that even though cyano groups stabilize radical centers, contribution of this effect in overall stabilization of the tetraradical electronic structure is not determinative.The determining factors for the stability of tetraradical electronic structure of LS and LSH are cross-conjugated π-topology and aromaticity.To explore the wavefunction of states of LSH , CASCI (4,4) and CASCI (16,16) were run from CASSCF (16,16) Q 0 converged state with (4,4) frontier localized subspace (see the orbital numbering assignations in Figure S4).CASCI (4,4) results are given in Table S13 and CASCI (16,16) results are given in Table S14.
Table S11.CASSCF(4,4) states of LSH determined by state-specific calculations.For all the calculations, UHF quintet NOs were used as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to Figure S5 and NO occupation number.

State Symmetry
Energy (a.u.) CASSCF NOs symmetry, identity and occupation number ∆E from G. S. (cm −1 ) HONO -1 HONO LUNO LUNO + 1 1.000 1.000 312.73 Overall electron density corresponds to the same as if they were Table S12.CASSCF (16,16) states of LSH determined by state-specific calculations.For all the calculations, UHF quintet NOs were used as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to Figure S5 and NO occupation number.

State Symmetry
Energy (a.u.) CASSCF NOs symmetry, identity and occupation number ∆E from G. S. (cm −1 ) HONO -1 HONO LUNO LUNO + 1 Overall electron density corresponds to the same as if they were Table S13.CASCI(4,4) of LSH from CASSCF (16,16) quintet converged orbitals as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to Figure S5 and NO occupation number.

State
1.000 A D We can also determine exchange-coupling constants between the radical centers with exactly the same procedure as we did for LS , as explained in Section S3.2.Since the procedure of CASSCF and CASCI calculations have been the same including the ordering of localized frontier orbitals from which CASCI calculations were run, we Table S15.CASCI(4,4) determinant contributions for 6 near-degenerate tetraradical(oid) spin states of LSH.These states are obtained for Sz = 0 basis from CASSCF (16,16) converged quintet localized (in (4,4) SOMOs subspace) initial guess orbitals.can use the same model Hamiltonian for LSH as we did for LS .Determinant contributions from CASCI(4,4) and CASCI (16,16) results for LSH are given in Table S15 and Table S16, respectively.If we compare composition of the wavefunction for each state for LSH , they parallel those of LS for the corresponding states very well.From CASCI(4,4) results, we find that J 12 = J 34 = J a = 141.26cm −1 , J 13 = J 24 = J b = 17.86 cm −1 , J 23 = J L = 90.74cm −1 , J 14 = J S = 627.02cm −1 .From CASCI (16,16) results, we find that J 12 = J 34 = J a = 116.92cm −1 , J 13 = J 24 = J b = 19.02cm −1 , J 23 = J L = 75.56cm −1 , J 14 = J S = 521.08cm −1 .The ordering of strengths of these exchange-coupling constants is same for LSH as for LS .Hence, the scheme for LSH would resemble the scheme of LS as shown in Figure S11a with moderated color shades to emphasize difference from LS .The low-energy spin spectrum is for LSH is given in Figure S11b.We find that the range of energy of states in lowenergy spectrum increases by about 390 cm −1 and exchange-couplings are also stronger in LSH compared to LS .Nonetheless, even though cyano substituents increase the tetraradical character slightly and decrease the range of energy of states in low-energy spectrum, their effect on the electronic structure of the tetraradical(oid) is rather small.Therefore, the tetraradical(oid) electronic structure of both LS and LSH is stabilized by cross-conjugation and aromaticity.
S7. Study of Diradical(oid)s L and S In order to gain further insight in the properties of this tetraradical(oid), it could be useful to analyze its conceptual constituent diradical(oid)s.The structure for diradical(oid) L is given in Figure S12 and the structure of S is given in Figure S13.Their geometries have been optimized with BLYP/TZP level of theory in ADF in the triplet state because of their diradical(oid) nature.Note that L has D 2h point group symmetry and was optimized as such, and S has C 2h symmetry and was optimized accordingly as well.CASSCF (14,14) and CASCI (14,14) (or CASCI (16,16)) calculations were run for each diradical(oid) molecule.All calculations were initiated from UHF triplet NOs guess as it is one of the best initial guess for diradical systems. 27The results for S are given in Table S17 for CASSCF (14,14) and in S18 for CASCI (14,14).As we see, we have two spin states, S 0 with A g symmetry as a ground state and T 0 with B u symmetry as a first excited state with singlet-triplet gap (∆E S−T = E T − E S ) amounting to only 1.63 kcal/mol according to CASSCF (14,14) and 0.77 kcal/mol according to CASCI (14,14).Higher excited states are too distant from these two states by more than 75 kcal/mol for each case.Hence, this molecule S is a persistent diradical(oid).The results for molecule L are given in Table S19 for CASCI (16,16) calculations and in Table S20 for CASSCF (14,14).We find that L has two very close-lying spin states T 0 with B 3u symmetry and S 0 with A g symmetry with singlet-we described two subsystems S and L, which are independently diradical(oid)s and through their conceptual union we obtained tetraradical(oid) LS , which has 6 spin states in the range of 603 cm −1 .The merge of these diradical(oid)s is synergistic towards maintaining unpaired electrons to form a tetraradical(oid) as reflected by significant cross-couplings between the diradical(oid) subunits.

S8. Optimized Geometries
The optimized geometry of tetraradical(oid) LS with BLYP/TZP quintet state is given below.This molecule has C 2h point group symmetry and has been optimized as such.As mentioned in the methodology, planar geometry results irrespective of whether the starting geometry is planar or not.Total bonding energy: -18.797140365861537 a.u.

Figure S2 .
Figure S2.Localized singly occupied molecular orbitals of converged solution of ROHF quintet obtained from starting with converged UHF quintet NOs guess.

Figure S3 .
Figure S3.UHF SONOs of converged solution of UHF singlet obtained from starting with converged ROHF quintet localized frontier orbitals guess.

Figure S4 .
Figure S4.Localized singly occupied molecular orbitals of converged solution of CASSCF(16,16) quintet obtained from starting with converged UHF quintet NOs guess.

Figure S5 .
Figure S5.Symbolic assignations to singly cccupied natural orbitals that appear for frontier orbitals in the CASSCF and CASCI solutions.

Figure
Figure S6.(a) Scheme of strength of exchange-coupling interactions in the tetraradical(oid) LS.Red is the strongest, blue is the second strongest closely followed by green and the weakest is gray.Rods sticking out from the horizontal line (1 -4) correspond to the indole nitrogens.(b) low-energy spectrum of the tetraradical(oid) LS.

Figure
Figure S11.(a) Scheme of strength of exchange-coupling interactions in the tetraradical(oid) LSH.Dark red is the strongest by more than factor of 4 relative to dark blue that is the second strongest, followed by moderate green and the weakest is gray.Rods sticking out from the horizontal line (1 -4) correspond to the indole nitrogens.(b) low-energy spectrum of the tetraradical(oid) LSH.

Table S1 .
CASSCF(4,4)states determined by state-specific calculations.For all the calculations, UHF quintet NOs were used as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to FigureS5and NO occupation number.S + 1.155A u L + 1.038 B g L − 0.961 B g S − 0.846 0.00 T 0 B u -2295.6575974661A u S + 1.155 B g L − 1.001A u L + 0.998 B g S − 0.846 24.34 S 1 A g -2295.6565732806A u C 0 1.046 B g C 1 1.011 B g C

Table S3 .
CASSCF(14,14)states determined by state-specific calculations.For all the calculations, UHF triplet NOs were used as initial guess.Under HONO, LUNO, etc. three columns represent symmetry of the particular orbital, its symbolic representation according to FigureS5and NO occupation number.

Table S9 .
DFT with variable HF % in PBE hybrid and UHF benchmark for all possible qualitatively distinct determinants for Sz = 0, 1, 2. Energies under each custom hybrid functional are given in atomic units.All quintets (|αααα⟩) have proper qualitative structure.Results with closed-shell qualitative structure are in bold, results with diradical(oid) structure are in bold italic, and results with proper tetraradical(oid) qualitative structure are in italic fonts.
The triplet-optimized D 2h geometry of molecule L is given below (BLYP/TZP level of theory).Total bonding energy: -11.131136314102742 a.u.The triplet-optimized C 2h geometry of molecule S is given below (BLYP/TZP level of theory).Total bonding energy: -10.315426748220721 a.u.