  A consistent microscopic theory for the response of strongly-coupled
superfluid fermionic systems is formulated. After defining the response as a
two-point two-fermion correlation function in the basis of the Bogolyubov's
quasiparticles, the equation of motion (EOM) method is applied using the most
general fermionic Hamiltonian with a bare two-body interaction, also
transformed to the quasiparticle space. As a superfluid extension of the case
of the normal phase, the resulting EOM is of the Bethe-Salpeter-Dyson form with
the static and dynamical interaction kernels, where the former determines the
short-range correlations and the latter is responsible for the long-range ones.
Both kernels as well as the entire EOM have the double dimension as compared to
that of the normal phase. Non-perturbative approximations via the cluster
decomposition of the dynamical kernel are discussed, with the major focus on a
continuous derivation of the quasiparticle-phonon coupling variant of the
latter kernel, where the phonons (vibrations) are composite correlated
two-quasiparticle states unifying both the normal and pairing modes. The
developed theory is adopted for nuclear structure applications, such as the
nuclear response in various channels. In particular, the finite-amplitude
method generalized beyond the quasiparticle random phase approximation, taking
into account the quasiparticle-vibration coupling, is formulated for
prospective calculations in non-spherical nuclei.
