How does this FLL work?
The conjugate product of the two filters form a frequency discriminator with an error output that is proportional to a frequency offset at the input.

This is due to the relationship that the cross correlation in one domain is the conjugate product in the other domain (you can skip past the derivation unless interested):

$$c(\tau) = XCORR[x(t)y(t)]$$
$$c(\tau) = \int_{t=-\infty}^{\infty}x^*(t)y(t+\tau)dt$$
$$C(\omega) = \mathscr F\{y(t)\}= \int_{\tau=-\infty}^{\infty}\bigg[\int_{t=-\infty}^{\infty}x^*(t)y(t+\tau)dt\bigg]e^{-j\omega \tau}d\tau$$
$$=\int_{t=-\infty}^{\infty}x^*(t)\bigg[\int_{\tau=-\infty}^{\infty}y(t+\tau)e^{-j\omega \tau}d\tau\bigg]dt$$

From the shift property of the Fourier Transform, the terms inside the brackets are:

$$\int_{\tau=-\infty}^{\infty}y(t+\tau)e^{-j\omega \tau}d\tau = Y(\omega)e^{j\omega t}$$

Therefore:
$$C(\omega) =\int_{t=-\infty}^{\infty}x^*(t)\bigg[Y(\omega)e^{j\omega t}\bigg]dt$$

$$ = Y(\omega)\int_{t=-\infty}^{\infty}x^*(t)e^{j\omega t}dt  = Y(\omega)\bigg[\int_{t=-\infty}^{\infty}x(t)e^{-j\omega t}dt\bigg]^*$$

$$ =Y(\omega)X^*(\omega) $$

The above shows that the cross correlation of two functions x(t) and y(t) in the time domain is equivalent to the conjugate multiplication of their Fourier Transforms (in the frequency domain). Similarly, and this is the important point to explaining the operation of the frequency discriminator, is that the converse is also true: the cross correlation of two functions in the frequency domain is equivalent to the conjugate multiplication of the transform of the same two functions in the time domain.

Consider the following as example spectrums of the "Matched Filter" and "Frequency Derivative Matched Filter" in the block diagram given by the OP. Notably the Frequency Derivative Matched Filter is a Band-Edge Filter.



When there is no frequency error, the input spectrum would be aligned with both of these resulting in the following spectrums at the outputs of these filters. Notably, observe that convolving these two output spectrums will result in a zero result specifically due to the positive and negative balance in the Frequency Derivative Matched Filter Spectrum!



When there is a positive frequency error in the input signal as depicted in the graphic below, the input spectrum would be offset by the frequency error resulting in the following spectrums at the output of the two filters. In this case specifically, convolving the two output spectrums would result in a negative value corresponding to the frequency error.  Similarly a negative frequency error would result in a positive output value after the convolution in frequency of the two spectrums. 



The above shows how the convolution of a Matched Filter Spectrum and a Frequency Derivative Matched Filter Spectrum results in an output value that is the error term for a frequency discriminator. Additionally we showed that convolution in frequency is equivalent to the complex conjugate multiplication in time. Thus the complex conjugate multiplication of the two filters in the block diagram form a frequency discriminator as the error signal into the loop. 

The frequency discrimintor would operate as described without the scaling by the Signal to Noise Ratio (SNR) and Tanh() operation. 

The additional scaling by the SNR, and the Tanh() operation are done to perform maximum likelihood frequency estimation (MLE). This shows the ideal implementation (in a maximum likelihood sense) but to my understanding the SNR scaling is rarely done in practice (but could be if an SNR estimation is performed as part of the receiver operation), and the Tanh() operation is typically approximated with a piecewise linear operation to be the signal for the small signal case when the signal is less than ±1 or limited to the sign of the signal (+1/-1) otherwise.  So basically this path in the block diagram is hard limited for better performance with a limit that is dependent on the SNR of the received signal (a SNR dependent hard limiter is a typical operation within MLE implementations). 

I believe that fred harris (fh) points out an advantage to this approach is how the derivative matched filter response is concentrated at the band edge where as we see through the convolution operation most of the information regarding frequency offset resides (when using a complex conjugate multiplier as a discriminator). The filter naturally blanks out the center portion of the spectrum which is not contributing further the our frequency error answer but would be contributing noise otherwise. Other advantages to this approach are that the symbol timing clock (timing recovery) can also be derived. I have not studied it personally in detail nor implemented this approach, but I believe fh makes the case that this implementation results in much lower pattern noise versus other carrier recovery (and timing recovery) approaches.