To conclude, the arithmetic averaging between S||H L -||S and S||H L + ||S * restores the correct symmetry of these elements. Usage of the triplet transition density does not lead to contamination of S||H Lz ||S matrix element with the unphysical singlet contribution. We presented a theory for calculating matrix elements of the BP SOC operator based on the application of Wigner-Eckart's theorem to reduced one-particle density matrices. The key equations are given in the second-quantized form. The algorithm is ansatz-agnostic and can be used with any electronic stricture method for which state and transition density matrices are available. The current implementation is based on the EOM-CCSD suite of methods. The approach allows one to compute the SOC matrix for the entire multiplet from just one transition density matrix, which solves the problem of accessibility of states of different spin projections. It also addresses, by constriction, the phase problem due to Born-Oppenheimer's separation of the nuclear and electronic degrees of freedom. We also highlighted an important aspect of calculation of SOCCs within a non-Hermitian theory and proposed using arithmetic averaging of the EOM-CC matrix elements, which leads rotationally invariant SOCCs.