Short answer As stated, the identity need not hold in general: the operator S need not be of trace class. A counterexample and a correction are given by Lyons in the erratum to his paper. The identity does hold (and S is trace class) under additional structure on the marginal metric spaces, namely when both metrics admit an isometric embedding into a Hilbert space (i.e., are of negative type so that the metric equals a squared Hilbert distance). Under those assumptions one can factor S and compute its trace, which equals the integral of the diagonal kernel; consequently the sum of eigenvalues equals that integral.

Details and a proof under a natural sufficient condition

Assume:

• There are separable Hilbert spaces H_X, H_Y and measurable maps Φ_X: X → H_X, Φ_Y: Y → H_Y such that d_X(x,x') = ||Φ_X(x) − Φ_X(x')||^2 and d_Y(y,y') = ||Φ_Y(y) − Φ_Y(y')||^2. (This is the “isometric embeddability into a Hilbert space” condition; equivalently, d_X and d_Y are of negative type.)
• μ, ν have finite first moments with respect to d_X, d_Y. For squared Hilbert distances this means ∫ ||Φ_X(x)||^2 dμ(x) < ∞ and ∫ ||Φ_Y(y)||^2 dν(y) < ∞.

Define m_X = ∫ Φ_X dμ, m_Y = ∫ Φ_Y dν and the centered maps g_X(x) = Φ_X(x) − m_X, g_Y(y) = Φ_Y(y) − m_Y.

A direct expansion shows the “double centering” identity d_μ(x,x') = d_X(x,x') − ∫ d_X(x,z) dμ(z) − ∫ d_X(x',z) dμ(z) + ∬ d_X(z,z') dμ(z) dμ(z') = −2⟨g_X(x), g_X(x')⟩, and similarly d_ν(y,y') = −2⟨g_Y(y), g_Y(y')⟩.

Let A_X: L^2(μ) → L^2(μ) be the integral operator with kernel d_μ, and A_Y analogously. Then S is the tensor product S = A_X ⊗ A_Y on L^2(μ × ν) ≅ L^2(μ) ⊗ L^2(ν), since its kernel is d_μ(x,x') d_ν(y,y').

Factor A_X as follows. Define T_X: H_X → L^2(μ) by (T_X v)(x) = ⟨g_X(x), v⟩. Since ∫ ||g_X(x)||^2 dμ(x) < ∞, T_X is Hilbert–Schmidt with ||T_X||_HS^2 = ∫ ||g_X(x)||^2 dμ(x). The adjoint T_X*: L^2(μ) → H_X is given by T_X* h = ∫ h(x) g_X(x) dμ(x). Then (T_X T_X* h)(x) = ⟨g_X(x), ∫ h(x') g_X(x') dμ(x')⟩ = ∫ ⟨g_X(x), g_X(x')⟩ h(x') dμ(x'). Hence A_X = −2 T_X T_X*. Therefore A_X is trace class and tr(A_X) = −2 tr(T_X T_X*) = −2 tr(T_X* T_X) = −2 ||T_X||_HS^2 = ∫ d_μ(x,x) dμ(x) (the last equality uses d_μ(x,x) = −2 ||g_X(x)||^2).

The same construction gives that A_Y is trace class with tr(A_Y) = ∫ d_ν(y,y) dν(y).

Standard facts about tensor products of trace-class operators (see, e.g., Simon, Trace Ideals and Their Applications) say that A_X ⊗ A_Y is trace class, that its nonzero eigenvalues are precisely the pairwise products {α_i β_j} of the eigenvalues of A_X and A_Y (counted with multiplicity), and that tr(S) = tr(A_X ⊗ A_Y) = tr(A_X) tr(A_Y).

Putting these together, sum_i λ_i = tr(S) = tr(A_X) tr(A_Y) = (∫ d_μ(x,x) dμ(x)) (∫ d_ν(y,y) dν(y)) = ∫ d_μ(x,x) d_ν(y,y) d(μ × ν)(x,y), which is exactly the desired identity.

Remarks and references

• Without the Hilbert-embedding assumption, S need not be trace class and the identity can fail. See R. Lyons, Errata to “Distance covariance in metric spaces” (2013), which provides a counterexample to S being trace class in general.
• Under the Hilbert-embedding assumption, the above argument is essentially what is given in that erratum: A_X and A_Y are trace class via the factorization A = −2 T T*, and S = A_X ⊗ A_Y is trace class with trace equal to the stated integral.
• Two useful references for the operator-theoretic facts used above are:
  • B. Simon, Trace Ideals and Their Applications, 2nd ed., especially the sections on Hilbert–Schmidt and trace-class operators and tensor products.
  • M. Reed and B. Simon, Methods of Modern Mathematical Physics I: Functional Analysis, sections on compact, Hilbert–Schmidt, and trace-class operators.

A small identity you might find handy: for any μ with finite first moment, ∫ d_μ(x,x) dμ(x) = − ∬ d_X(x,x') dμ(x) dμ(x'), so the right-hand side of your desired identity equals the product of those two double integrals. This is consistent with the trace computation above.