The \(M\)-dimensional unitary matrix \(S(E)\), which describes scattering of waves is in general a strongly fluctuating function of the energy, specially for complex systems such as ballistic cavities whose geometry induces chaotic ray dynamics. Its statistical behaviour can be expressed by means of the correlation functions of the kind \(\left \langle S_{ij}(E)S^\dag_{pq}(E)\right\rangle\). These correlations, involving an arbitrary number of matrix elements, have been much studied within the random matrix approach, in which \(S(E)\) is taken to be uniformly distributed either in the Circular Unitary or Circular Orthogonal ensembles, depending on whether or not time-reversal symmetry is present. In this work, we generalize these results to the situation when matrix elements taken at two different energies, e.g. \(\left \langle S_{ij}(E+\epsilon)S^\dag_{pq}(E-\epsilon)\right\rangle\). We express these functions as infinite series in \(1/M\), whose coefficients are rational functions of \(\epsilon\). From a mathematical point of view, this may be seen as a parametric generalization of the Weingarten functions of circular ensembles.