Why does the Cauchy distribution model resonance? I can see that being the quotient of two independent normal distributions it may model things that spin (although I don't understand how it can be well defined at zero).
But what is the connection to the solution of differential equations modelling resonance?
 A: Given the absence of more knowledgeable answers, I will share what I found on the topic - hopefully better or more more comprehensive answers will come along.
In short, free-induction decay (FID), which results from interactions among dipoles in MRI (magnetic resonance imaging) is influenced by inhomogeneities in the magnetic field which have long tails (not exponential or Gaussian), and can be modelled using a Cauchy distribution:

... The distribution of microscopic inhomogeneous fields is not necessarily Gaussian. In some systems of importance in areas as diverse as oilwell
logging and medical magnetic resonance imaging (MRI) the distribution of inhomogeneous fields is approximated by the Cauchy-Lorentz (C-L, Cauchy in statistics, Lorentz in physics) distribution. This is the
field due to randomly located point dipoles within a sample. Again, the phase distribution width at an effective exchange time $t_{ex}$ is proportional to $\chi \nu t_{ex}.$ However, the distribution width for $n = t/t_{ex}$ random samples is $t/t_{ex}$ times the original width, so the width of
the phase distribution at time $t$ does not depend on $t_{ex}.$ Here, the Fourier transform of the $n$-fold convolution of a C-L distribution is the $n$-fold product of its simple exponential FT, not an exponential with the width squared as in the Gaussian. That is, the FID decay does not depend on $a$ or $D.$ Furthermore, the additional FID decay rate is proportional to the width of the CL distribution of phases, not to the square of the width, as for the Gaussian.

where $\chi\nu$ (susceptibility difference times frequency), $D$ (diffusion coefficient), $a$ (a pore dimension), echo-time $t$ for single echoes and half-echo-spacing $\tau$ for CPMG (Carr-Purcell-Meiboom-Gill (CPMG) spin-echo measurements), $t_{ex}$ for exchange time.

More relevant, and in reference to a more robust correspondence between time and frequency domains in MRI than the FFT (e.g. non-parametric fast Padé transform)

Considered strictly, a general line spectrum is associated only with for purely stationary, stable bound states (e.g. ground states of atoms/molecules) of practically infinite lifetimes. On the other hand, peaks in spectra are resonances that represent non-stationary, transient states of molecules. Such metastable states last for a given lifetime, after which they decay. This is reflected in distributions of frequencies in the spectral lineshapes, instead of the stick spectra. These distributions are mainly the Lorentzian functions that are in mathematics known as the Cauchy distributions. In physics, they are also called Breit-Wigner lineshape profiles from the theory of resonant scattering of particles and photons on general targets. The Cauchy, Lorentz and Breit-Wigner distributions coincide with one of the simplest forms of the Padé approximant (PA), symbolized as [0/2], and given by:
$$\begin{aligned} Y_{\mathrm{A}} \left( \nu \right) =\frac{d_0 }{2\pi }\frac{\varGamma _0 }{\left( {\nu -\nu _0 } \right) ^{2}+\left( {\varGamma _0 /2} \right) ^{2}}\quad \hbox {(Absorptive Lorentzian)} \end{aligned}$$
It is seen that the Padé spectrum $Y_A(\nu)$ depends on $\nu$ in a continuous way. If $Y_A(\nu)$ is written in the general form $P_L(\nu)/Q_K(\nu),$ it would follow that $P_0(\nu)=d_0Γ_0//(2π)$ and $Q_K(\nu)=Q_2(\nu)=(\nu−\nu_0)^2+(Γ_0/2)^2.$ Here, $d_0$ is the amplitude, $\nu_0$ and $Γ_0$ are the resonant frequency and the full width at half-maximum (FWHM), respectively. All three parameters $\nu_0,$ $Γ_0$ and $d_0$ are real-valued quantities, with the understanding that the phase of the generally complex amplitude is equal to zero. The notation [L/K] is the usual symbol for the general PA as a ratio of two polynomials $P_L$ and $Q_K$ of degrees $L$ and $K,$ respectively. The physical meaning of a stick or bar spectrum is that it consists of lines associated with mono-energetic (mono-chromatic) absorption or emission of radiations. These latter processes result from transitions between discrete rotational, vibrational or rovibrational molecular levels for the ground state as well as for excited states. As noted, such line spectra belong to pure bound, stable states of molecules. Metastable states are prone to decay and their spectra are not lines or sticks. Rather, due to the underlying dissipative dynamics (there is a loss of energy), the absorptive radiations cannot be mono-energetic. This implies that the resulting spectra do not have sharp lines, in lieu of which a distribution of energies (or frequencies) appears. The most typical example of such frequency distributions is the absorptive Lorentzian function. In other words, the Lorentzian lineshape centered at $\nu_0$ is a broadened line of breadth or full width $Γ_0.$ These notions are also familiar by reference to a vibrating dipole which radiates energy according to classical physics. As a result, the vibrational amplitudes decrease. In other words, the dipole vibrations are attenuated, and the lineshape of the emitted radiation energy or frequency is the same Lorentzian distribution instead of a simple mono-chromatic phenomenon.

In comparing to the pdf of the Cauchy distribution, the equation for the Absorption Lorentzian above would be encapsulate the density of the random variable resonant frequency, $\nu,$ with scale parameter $\gamma = 2,$ and location $x_0=\nu_):$
$$f_X(x) = \frac 1 {\pi \gamma\left ( 1 + (\frac {x-x_0}{\gamma}) \right)^2}$$
