  We perform theoretical and numerical studies of the full relativistic
two-point galaxy correlation function, considering the linear-order scalar and
tensor perturbation contributions and the wide-angle effects. Using the
gauge-invariant relativistic description of galaxy clustering and accounting
for the contributions at the observer position, we demonstrate that the
complete theoretical expression is devoid of any long-mode contributions from
scalar or tensor perturbations and it lacks the infrared divergences in
agreement with the equivalence principle. By showing that the gravitational
potential contribution to the correlation function converges in the infrared,
our study justifies an IR cut-off $(k_{\text{IR}} \leq H_0)$ in computing the
gravitational potential contribution. Using the full gauge-invariant
expression, we numerically compute the galaxy two-point correlation function
and study the individual contributions in the conformal Newtonian gauge. We
find that the terms at the observer position such as the coordinate lapses and
the observer velocity (missing in the standard formalism) dominate over the
other relativistic contributions in the conformal Newtonian gauge such as the
source velocity, the gravitational potential, the integrated Sachs-Wolf effect,
the Shapiro time-delay and the lensing convergence. Compared to the standard
Newtonian theoretical predictions that consider only the density fluctuation
and redshift-space distortions, the relativistic effects in galaxy clustering
result in a few percent-level systematic errors beyond the scale of the
baryonic acoustic oscillation. Our theoretical and numerical study provides a
comprehensive understanding of the relativistic effects in the galaxy two-point
correlation function, as it proves the validity of the theoretical prediction
and accounts for effects that are often neglected in its numerical evaluation.
