What is the most natural way in which the Jarlskog invariant makes its appearance?

What is the most natural way in which the Jarlskog invariant can be introduced? I have seen ad hoc definitions given as $$J_{CP}=\text{Im}(V_{us}V_{cb}V^*_{ub}V^*_{cs})\tag{1}$$ which, to me, appears out of the blue (For example, Matt Schwartz's book of QFT, page 599). Does $J_{CP}$ appear in a natural way in calculating some mathematical expression for some physical observable? Ad hoc introduction of $J_{CP}$ via (1) doesn't seem to motivate why should such a non-trivial quantity be constructed.

I know that the Jarlskog invariant is a measure of the CP-violation but it's not quite clear why is this combination of elements gives the unique measure of CP-violation and there is no other measure.

Well, it is an ad hoc quantity, even though quite intuitive, and hence a popular one. It certainly featured first in $K^0-\overline {K^0}$ strangeness-oscillation amplitudes,

where the imaginary part of the 4 CKM matrix elements' product underlies CP violation. You need not rush to unobserved phenomena, when the archetype for CP violation in the quark sector all but shouts it at you: this CP-violating amplitude is computed to be linear in J. This is the very "hoc" of your "ad hoc".

As Matt takes care to detail, the encoder of this violation is the Hermitian and traceless commutator matrix of the up and down Yukawa matrices, and specifically its basis-invariant determinant--its imaginary part. One of the factors is the minuscule Jarlskog invariant, of the order of $3\cdot 10^{-5}$, twice the area of its CKM-matrix unitarity triangle.

It is so small by virtue of all 3 mixing sines of all 3 generations (so the t-quarks in the intermediate box propagators are crucial players, not an ignorable small correction, for this), and the smallness of the 2-3 and especially 1-3 angles: $J=s_{12}s_{23}s_{13}c_{12}c_{23}c_{13}\sin \delta$. (More dramatically, and intuitively, it is of $O(\lambda^6)$ of the Wolfenstein parameterization. Now, that's small.) So it follows that for 3 generations it is generic and systematic, and not ad hoc.

However, extension to 4 or more generations is meaningless or problematic; and it is not the only player in CP violation... That's why it is heuristic. It is only half of the full story.

The above determinant of the commutator is also proportional to the Vandermond matrix determinant of the actual masses of the 3 generations, as Matt soundly emphasizes, in his (29.90). This is the factor of the 6 pairwise mass differences among all 6 quarks, within the up and down sectors, normalized by the sixth power of the EW v.e.v. If the mass difference between any two flavors within an up or down flavor were small, then CP violation would be infinitesimal, regardless of how huge a hypothetical Jarlskog invariant were to be! Explicitly, this companion prefactor is $$\det \begin{pmatrix} 1 & m_u & m_u^2 \\ 1 & m_c & m_c^2 \\ 1 & m_t & m_t^2 \end{pmatrix} ~\cdot ~\det \begin{pmatrix} 1 & m_d & m_d^2 \\ 1 & m_s & m_s^2\\ 1& m_b& m_b^2 \end{pmatrix}=\\ (m_t-m_c)(m_c-m_u)(m_t-m_u)(m_b-m_s)(m_s-m_d)(m_b-m_d).$$

This is not an answer for the $J_{CP}$ in the quark sector but in the neutrino sector. This reference gives the oscillation probability $$P(\nu_e\to \nu_\mu)=4|U_{e3}|^2|U_{\mu 3}|^2\sin^2\Delta_{atm}-4\text{Re}\{U_{e1}U_{e2}^*U_{\mu 1}^*U_{\mu 2}\}\sin^2\Delta_{sol}-2J\sin(2\Delta_{sol})$$ ans similar expressions for $P(\nu_e\to \nu_\tau)$ and $P(\nu_\mu\to \nu_\tau)$. In these expressions, the Jarlskog invariant $J=J_{CP}$ appears quite naturally and is given by $$J=\text{Im}[U_{e2}U_{e3}^*U_{\mu 2}^*U_{\mu 3}].$$ Therefore, $J_{CP}$ makes its appearence while calculateing the observables such as neutrino and antineutrino oscillation probabilities. This to me, is a less ad hc way of introducing $J_{CP}$. I don't know whether there is such a natural way in which $J_{CP}$ appears in the quark sector.