  We present the implementation of quadratic response theory based upon the
relativistic equation-of-motion coupled cluster method. We showcase our
implementation, whose generality allows us to consider both time-dependent and
time-independent electric and magnetic perturbations, by considering the static
and frequency-dependent hyperpolarizability of hydrogen halides (HX, X = F-At),
providing a comprehensive insight into their electronic response
characteristics. Additionally, we evaluated the Verdet constant for noble gases
Xe and Rn, and discussed the relative importance of relativistic and electron
correlation effects for these magneto-optical properties. Finally, we calculate
the two-photon absorption cross-sections of transition ($ns^{1}S_{0}\to
(n+1)s^{1}S_{0}$) of Ga$^{+}$, and In$^{+}$, which are suggested as candidates
for new ion clocks. As our implementation allows for the use of
non-relativistic Hamiltonians as well, we have compared our EOM-QRCC results to
the QR-CC implementation in the DALTON code, and show that the differences
between CC and EOMCC response are in general smaller than 5\% for the
properties considered. Collectively, the results underscore the versatility of
our implementation and its potential as a benchmark tool for other approximated
models such as density functional theory for higher-order properties.
