optimization problem, any solution?
there is a solution to your problem. First of all, let us consider the problem $$\mathcal{P}:\\ z : = \max_{\mathbf{F}} Re\{\mathbf{b}\mathbf{F}^H \mathbf{C} \mathbf{F} \mathbf{d}\} {\rm ~~s.t.~~~} Tr(\mathbf{F} \mathbf{F}^H) \leq p. $$
Let $\mathbf{F}_\star$ the solution to the above problem. It is clear that if $ z < a Tr(\mathbf{F}_\star \mathbf{F}_\star^H) $, then $\mathbf{F} = \mathbf{0}$ is the solution to your problem. Otherwise, $\mathbf{F}_\star$ is the also the solution to your problem and $Tr(\mathbf{F}_\star {\mathbf{F}_\star}^H) = p$.
Now have a deeper look at $\mathcal{P}$. As $\mathbf{C}$ is a Hermitian matrix, there exists the eigendecomposition $\mathbf{C} = \mathbf{U} \Sigma \mathbf{U}^H $, where $\mathbf{U}$ is an orthonormal matrix and $\Sigma $ is a diagonal matrix with real-valued entries $\sigma_1,\dots,\sigma_N$. For simplification, and without loss of generality we assume that $\mathbf{F} = \mathbf{U} \mathbf{H}$. Then, we can rewrite $\mathcal{P}$ as
$$ \mathcal{Q}: \\
\max_{\mathbf{H}} Re\{\mathbf{b}\mathbf{H}^H \Sigma \mathbf{H} \mathbf{d}\} {\rm ~~s.t.~~~} Tr(\mathbf{H} \mathbf{H}^H) = p.
$$
here we have exploited the fact that $ \mathbf{U}^H \mathbf{U} = \mathbf{I}$ and $ tr( \mathbf{FF}^H)= tr( \mathbf{U}\mathbf{H} ( \mathbf{U} \mathbf{H})^H) = tr(\mathbf{HH}^H$). As $\Sigma$ is diagonal we can reformulate the cost function as
$ Re\{\mathbf{b}\mathbf{H}^H \Sigma \mathbf{H} \mathbf{d}\} = Re\{\sum_{n=1}^N \sigma_i \mathbf{b} \mathbf{h}_n^H \mathbf{h}_n \mathbf{d} \} = \sum_{n=1}^N Re\{\sigma_n \mathbf{h}_n \mathbf{d} \mathbf{b} \mathbf{h}_n^H \}$, where $\mathbf{h}_n$ is the $n$th column of $\mathbf{H}$. Further, we see that $Re\{\sigma_n \mathbf{h}_n \mathbf{d} \mathbf{b} \mathbf{h}_n^H \} = \mathbf{h}_n(\sigma_n\frac{1}{2}(\mathbf{d} \mathbf{b} +\mathbf{b} \mathbf{d})) \mathbf{h}_n^H$. Defining the matrix $\mathbf{Q}_n:= \sigma_n\frac{1}{2}(\mathbf{d} \mathbf{b} +\mathbf{b}^H \mathbf{d}^H)$, we can finally rewrite $\mathcal{Q}$ as
$$ \mathcal{R}: \\
\max_{\{\mathbf{h}_n\}_{n=1}^N} \sum_{n=1}^N \mathbf{h}_n\mathbf{Q}_n \mathbf{h}_n^H {\rm ~~s.t.~~~} \sum_{n=1}^N \| \mathbf{h}_n\|^2 = p.
$$
From the definition it is clear that each $\mathbf{Q}_n$ is a Hermitian matrix of rank two. All of these matrices have the (same) two eigenvectors $\mathbf{q}_1$ and $\mathbf{q}_2$ with (different ) eigenvalues $\lambda_{1,n}$ and $\lambda_{2,n}$. It is clear that vector components in $\mathbf{h}_n$ that are orthogonal to $\mathbf{q}_1$ and $\mathbf{q}_2$ have no influence on the cost function of $\mathcal{R}$ but increase the constraint value. Therefore, we assume without loss of optimality that $\mathbf{h}_n = c_{1,n}\mathbf{q}_1 +c_{2,n} \mathbf{q}_2 $. Then $\mathcal{R} $ boils down to
$$
\max_{\{{c}_{1,n},{c}_{2,n}\}_{n=1}^N} \sum_{n=1}^N (\lambda_{1,n}|c_{1,n}|^2+\lambda_{2,n}|c_{2,n}|^2) {\rm ~~s.t.~~~} \sum_{n=1}^N (|c_{1,n}|^2+|c_{2,n}|^2) = p.
$$
Solution
The solution to the above problem is simple: Let $i_\star,n_\star = \arg \max \lambda_{i,n}$, then $|c_{i_\star,n_\star}|^2= p $ and $|c_{i,n}|^2= 0 $. Consequently, the solution matrix $\mathbf{H}_\star$ has only one non-zero row (the $n_\star$th row) which is given by $(p)^{0.5} \mathbf{q}_{i_\star}$.
Consequently, $\mathbf{F}_\star =\mathbf{U} \mathbf{H}_\star $ and $z = p \lambda_{i_\star,n_\star} $.
All the best, Ert