  We present the development and implementation of the relativistic coupled
cluster linear response theory (CC-LR) which allows the determination of
molecular properties arising from time-dependent or time-independent electric,
magnetic, or mixed electric-magnetic perturbations (within a common gauge
origin), and take into account the finite lifetime of excited states via damped
response theory. We showcase our implementation, which is capable to offload
intensive tensor contractions onto graphical processing units (GPUs), in the
calculation of: \textit{(a)} frequency-(in)dependent dipole-dipole
polarizabilities of IIB atoms and selected diatomic molecules, with a emphasis
on the calculation of valence absorption cross-sections for the I$_2$
molecule;\textit{(b)} indirect spin-spin coupling constants for benchmark
systems such as the hydrogen halides (HX, X = F-I) as well the H$_2$Se-H$_2$O
dimer as a prototypical system containing hydrogen bonds; and \textit{(c)}
optical rotations at the sodium D line for hydrogen peroxide analogues
(H$_{2}$Y$_{2}$, Y=O, S, Se, Te). Thanks to this implementation, we are able
show the similarities in performance--but often the significant
discrepancies--between CC-LR and approximate methods such as density functional
theory (DFT). Comparing standard CC response theory with the equation of motion
formalism, we find that, for valence properties such as polarizabilities, the
two frameworks yield very similar results across the periodic table as found
elsewhere in the literature; for properties that probe the core region such as
spin-spin couplings, we show a progressive differentiation between the two as
relativistic effects become more important. Our results also suggest that as
one goes down the periodic table it may become increasingly difficult to
measure pure optical rotation at the sodium D line, due to the appearance of
absorbing states.
