We present in this paper a comprehensive account of an explicitly spin-free coupled cluster theory for treating energy differences of open-shell states relative to a closed-shell ground state, where the open-shell states of interest are dominated by a few simple configuration state functions. We develop a valence-universal coupled cluster formalism to achieve this via a novel cluster expansion ansatz for the valence part of the wave operator, where the orbital relaxation and the correlation relaxation accompanying ionization/excitation from the ground state are taken care of to all orders in compact, efficient, and explicitly spin-free manner. The essential difference of our proposed ansatz from the ordinary and the normal-ordered cluster ansatz in vogue is that (a) we allow the valence cluster operators to be connected among themselves with spectator valence lines only and (b) we use suitable combinatoric factors accompanying powers of cluster operators thus connected, which are equal to the number of ways the operators can be joined, leading to the same excitation (the automorphic factor). We emphasize that such an ansatz does not generate terms (diagrams) with chains of cluster operators joined among themselves via spectator lines only. Barring only a few, almost all the terms in the working equations determining the cluster amplitudes involve contraction of the Hamiltonian with the cluster operators via at least one nonspectator line, leading to what we call a "strongly connected" series. The structure of the working equation is remarkably similar to the single-reference closed-shell equation, with a few additional terms. The presence of contractions among cluster operators via spectator lines introduces the additional physical effects of orbital and correlation relaxation using low-body cluster operators. As an illustrative application of the new multireference coupled cluster (CC) theory, we consider in this paper computation of ionization potentials (IPs) of one-valence problem with only one active orbital. The numerical applications are made for both the core- and the inner- and outer-valence IPs for several molecular systems. The numerical values demonstrate the superiority of the relaxation-inducing CC theory, as compared to the normal-ordered ansatz.