These results hold for absolutely continuous random variables (see this question for the sample median behavior for discrete random variables). $F$ will be the cumulative distribution function, $f$ the probability density function, $\mu$ the mean. A hat will indicate sample quantities, and a bar sample means. Denote the quantile associated with probability $p$ by $\xi_p$, so that the median is $\xi_{1/2}\equiv \theta$.
The Bahadur representation of a sample quantile is
$$\hat \xi_p = \xi_p +\frac{1-\hat F_n(\xi_p)-(1-p)}{f(\xi_p)} +R_{n,p},\;\; R_{n,p} =o(1/\sqrt n)$$
$$\Rightarrow \sqrt n(\hat \xi_p - \xi_p) = \sqrt n\frac{1-\hat F_n(\xi_p)-(1-p)}{f(\xi_p)} +\sqrt nR_{n,p} \tag{1}$$
Define the indicator function $I_i\equiv I(X_i>\xi_p)$ (this is the complementary indicator function of the one the OP uses). We have that
$$E(I_i)= P(X_i>\xi_p) = 1-p,\;\; \operatorname{Var}(I_i)=p(1-p)\tag {2}$$
and also
$$\bar S_I \equiv \frac 1n\sum_{i=1}^nI_i = 1-\hat F_n(\xi_p), \;\; E(\bar S_I) = 1-p, \operatorname{Var}(\bar S_I) = p(1-p)/n \tag {3}$$
Using $(2)$ and $(3)$ we can write $(1)$ as
$$\sqrt n(\hat \xi_p - \xi_p) = \sqrt n\frac{1}{f(\xi_p)}\left(\bar S_I-E(\bar S_I)\right) +\sqrt nR_{n,p} \tag{4}$$
The quantity in the big parenthesis is a sample mean centered, the term in front is a constant, the last term asymptotically vanishes, so we are led to the CLT for sample medians. But this is not what we are after.
Consider now the bivariate random vector $(X_i, I_i)'$. The covariance between the two components is
$$\operatorname{Cov}(X_i,I_i) = E(X_iI_i) - \mu(1-p) \tag{5}$$
By an application of a multivariate CLT, $(5)$ will also be the covariance of the asymptotic distribution of the bivariate vector $\sqrt n\big(\bar X -\mu,\, \bar S_I-E(\bar S_I)\big)'$.
We are interested in the (asymptotic) quantity
$$\operatorname{Cov}\Big(\sqrt n(\bar X -\mu),\sqrt n(\hat \xi_p - \xi_p)\Big) = \operatorname{Cov}\Big(\sqrt n(\bar X -\mu),\sqrt n\frac{1}{f(\xi_p)}\left(\bar S_I-E(\bar S_I)\right)\Big) $$
$$=\frac{1}{f(\xi_p)}\operatorname{Cov}\Big(\sqrt n(\bar X -\mu),\sqrt n\left(\bar S_I-E(\bar S_I)\right)\Big)$$
$$= \frac{1}{f(\xi_p)}\left(E(X_iI_i) - \mu(1-p)\right) \tag {6}$$
For the sample median, $(6)$ becomes
$$\operatorname{Cov}\Big(\sqrt n(\bar X -\mu),\sqrt n(\hat \theta - \theta)\Big) = \frac{1}{2f(\theta)}\left(2\int_{\theta}^{\infty}xf(x)dx - \mu\right) \tag{7}$$
Eq. $(7)$ is the same as what the OP found (remember, we have defined the indicator function complementarily). Now, Ferguson's paper gives this quantity as $E(|X-\theta|)/2f(\theta)$. We have to verify that this is the same as $(7)$.
$$E(|X-\theta|) = \int_{-\infty}^{\infty}|x-\theta|f(x)dx = \int_{-\infty}^{\theta}(\theta-x)f(x)dx+\int_{\theta}^{\infty}(x-\theta)f(x)dx$$
$$=\theta\int_{-\infty}^{\theta}f(x)dx - \int_{-\infty}^{\theta}xf(x)dx + \int_{\theta}^{\infty}xf(x)dx - \theta\int_{\theta}^{\infty}f(x)dx$$
Bring the terms not involving $x$ together and add and subtract $\int_{\theta}^{\infty}xf(x)dx$ to obtain
$$E(|X-\theta|) =\theta F(\theta) -\theta[1-F(\theta)] + \\
+2\int_{\theta}^{\infty}xf(x)dx - \int_{-\infty}^{\theta}xf(x)dx - \int_{\theta}^{\infty}xf(x)dx$$
$$= \theta/2 - \theta/2 +2\int_{\theta}^{\infty}xf(x)dx -\int_{-\infty}^{\infty}xf(x)dx$$
$$\Rightarrow E(|X-\theta|) = 2\int_{\theta}^{\infty}xf(x)dx -\mu$$
So Ferguson's expression for the asymptotic covariance between sample mean and sample median is correct, without needing to impose $\mu = \theta$.