Just because they don't have a covariance doesn't mean that the basic $x^t\Sigma^{-1} x$ structure usually associated with covariances can't be used. In fact, the multivariate ($k$-dimensional) Cauchy can be written as:
$$f({\mathbf x}; {\mathbf\mu},{\mathbf\Sigma}, k)= \frac{\Gamma\left(\frac{1+k}{2}\right)}{\Gamma(\frac{1}{2})\pi^{\frac{k}{2}}\left|{\mathbf\Sigma}\right|^{\frac{1}{2}}\left[1+({\mathbf x}-{\mathbf\mu})^T{\mathbf\Sigma}^{-1}({\mathbf x}-{\mathbf\mu})\right]^{\frac{1+k}{2}}} $$
which I have lifted from the Wikipedia page. This is just a multivariate Student-$t$ distribution with one degree of freedom.
For the purposes of developing intuition, I would just use the normalized off-diagonal elements of $\Sigma$ as if they were correlations, even though they are not. They reflect the strength of the linear relationship between the variables in a way very similar to that of a correlation; $\Sigma$ has to be positive definite symmetric; if $\Sigma$ is diagonal, the variates are independent, etc.
Maximum likelihood estimation of the parameters can be done using the E-M algorithm, which in this case is easily implemented. The log of the likelihood function is:
$$\mathcal{L}(\mu, \Sigma) = -{n\over 2}|\Sigma| - {k+1 \over 2}\sum_{i=1}^n\log(1+s_i)$$
where $s_i = (x_i-\mu)^T\Sigma^{-1}(x_i-\mu)$. Differentiating leads to the following simple expressions:
$$\mu = \sum w_ix_i/\sum w_i$$
$$\Sigma = {1 \over n}\sum w_i(x_i-\mu)(x_i-\mu)^T$$
$$w_i = (1+k)/(1+s_i)$$
The E-M algorithm just iterates over these three expressions, substituting the most recent estimates of all the parameters at each step.
For more on this, see Estimation Methods for the Multivariate t Distribution, Nadarajah and Kotz, 2008.