I am a little late but I post my answer anyway so that someone having the same question will find it interesting and discuss.
The discrete baseband multipath channel can be modeled as a FIR, i.e.
$$y[n] = \sum_{l=0}^{L-1} x[n-l] h_l + w[n]$$
where $L$ is the number of channel taps. $L$ depends on the relation between the bandwidth of basis waveform and the delay spread of channel.
The term "Rayleigh fading" channel implies that the channel taps $h_l$ can be modeled as i.i.d circular symmetric Gaussian complex random variables because :
- $h_l$ is the sum of a large number of small independent circular symmetric random variables, each random variable is channel of a physical path. This the rich scattering assumption, which is typically vaid for urban environment.
- there is not a particular path having gain much significant than others. Otherwise, we have Rician fading.
Let me call this random variables "Rayleigh".
With a sufficient cyclic prefix ("sufficient" means greater than delay spread, thus OFDM receiver captures all delayed versions of OFDM symbol, proof can be found at single-tap OFDM regardless subcarrier spacing), the demodulated data at subcarrier $k$ is
$$z[k] = x[k] \times \sum_{l=0}^{L-1}h_l e^{-j2\pi \frac{l}{N} k} = x[k] \times H[k]$$
where $N$ is DFT size.
The channel taps $h_l$ are i.i.d circular symmetric Gaussian random variables, $H[k]$ are circular symmetric Gaussian random variables, but they are in general not iid.
As pointed out by Maximilian Matthé in the comment, the covariance matrix is $F \mathrm{diag}(\vec{p}_0) F^H$ where $\vec{p}_0$ is Power Delay Profile zero-padded to size $N$. The frequency bins at $k = u \times N/L, u \in \mathbb{N}$ are independent, if $N/L$ is integer. Other ones are sinc-interpolated thus they are correlated. Note that $N/L \times \Delta f = 1 / L T_s \approx 1/\tau_m$ can be seen as coherence bandwidth.