hermitian matrix-vector product does not give real result

Background : If $H$ is a Hermitian matrix, then for any complex vector $\psi$ it holds that $\psi^* \cdot H\cdot \psi$ gives a real number.

The problem

I am creating a 16x16 Hermitian matrix in Mathematica (I am attaching code at the end of my question because it is quite large). Then, I construct some complex vectors and I compute the product stated above, e.i I do

FullSimplify[Conjugate[wp[x0, y0]].H0.wp[x0, y0]]


Where H0 is Hermitian matrix and wp[x0,y0] is a complex vector. My problem is that Mathematica cannot simplify the result to a real number, but complex exponentials are always present. Can somebody help me understand how I can use the simplification capabilities of Mathematica so that the code returns a real expression?

The code

Please note: I did not create a "minimal working example", because I don't care of simplifying an expression irrelevant to my problem and I can't be sure that whatever answer is given for a MWE would work for my case.

I "declare" some variables as real:

Im[\[Alpha]] ^:= 0; Conjugate[\[Alpha]] ^:= \[Alpha]
Im[x0] ^:= 0; Conjugate[x0] ^:= x0
Im[y0] ^:= 0; Conjugate[y0] ^:= y0
Im[B] ^:= 0; Conjugate[B] ^:= B
Im[t] ^:= 0; Conjugate[t] ^:= t
Im[kx] ^:= 0; Conjugate[kx] ^:= kx
Im[ky] ^:= 0; Conjugate[ky] ^:= ky
a1 = \[Alpha]*Sqrt[3]*{1, 0};
a2 = \[Alpha] * Sqrt[3]*{1/2, Sqrt[3]/2};
\[Delta] = {0, \[Alpha]}


Then I create the Hermitian matrix $H$ (Tight-Binding Hamiltonian for any experts)

pos = Range[16];
pos[[9]] = {x0, y0};
pos[[2]] = pos[[9]] - a1;
pos[[6]] = pos[[9]] - a2;
pos[[16]] = pos[[9]] + a1;
pos[[12]] = pos[[9]] + a2;
pos[[4]] = pos[[9]] - a1 + a2;
pos[[14]] = pos[[9]] + a1 - a2;
pos[[7]] = pos[[9]] + 2 a2 - a1;
pos[[10]] = pos[[9]] - 2 a2 + a1 + \[Delta]
Do[
pos[[i - 1]] = pos[[i]] + \[Delta],
{i, {2, 6, 4, 9, 14, 16, 12}}];

hmatrix[n_, m_] := Which[
n == m,
0,
Norm[pos[[n]] - pos[[m]]] === Abs[\[Alpha]],
Block[
{xn = pos[[n]][[1]], yn = pos[[n]][[2]],
xm = pos[[m]][[1]], ym = pos[[m]][[2]]},
-t*Exp[-I*B*((yn + ym)/2)*(xm - xn)]],
True,
0]

H = Table[hmatrix[i, j], {i, 16}, {j, 16}];
H0 = H /. B -> 0;


And indeed mathematica confirms that both H, H0 are Hermitian, as Transpose[Conjugate[H]] == H returns True.

I now create the complex vectors $\psi$:

pw = Table[Exp[I (kx*pos[[i]][[1]] + ky *pos[[i]][[2]])], {i, 16}];
gaussian[a_, b_] :=
Table[Exp[-(pos[[i]][[1]] - a)^2 - (pos[[i]][[2]] - b)^2], {i, 16}];
wp[a_, b_] := (g = gaussian[a, b]; Table[g[[i]]*pw[[i]], {i, 16}]);


Now, I compute the expression of interest, namely $\psi^* \cdot H\cdot \psi$. As expected, the first complex vector pw gives a real result,

FullSimplify[Conjugate[pw].H0.pw]


returns

$$-2 t \left(12 \cos \left(\frac{1}{2} \sqrt{3} \alpha \text{kx}\right) \cos \left(\frac{\alpha \text{ky}}{2}\right)+7 \cos (\alpha \text{ky})\right)$$

However, quite unexpectedly,

FullSimplify[Conjugate[wp[x0, y0]].H0.wp[x0, y0]]


does not return a real result, in the sense that Mathematica cannot simplify it:

$$t \left(-e^{-\frac{1}{2} \alpha \left(32 \alpha +i \sqrt{3} \text{kx}+2 i \text{ky}\right)}\right) \left(2 e^{6 \alpha ^2+\frac{1}{2} i \sqrt{3} \alpha \text{kx}}+2 e^{9 \alpha ^2+\frac{1}{2} i \sqrt{3} \alpha \text{kx}}+2 e^{12 \alpha ^2+\frac{1}{2} i \sqrt{3} \alpha \text{kx}}+e^{15 \alpha ^2+\frac{1}{2} i \sqrt{3} \alpha \text{kx}}+2 e^{6 \alpha ^2+\frac{1}{2} i \alpha \left(\sqrt{3} \text{kx}+4 \text{ky}\right)}+2 e^{9 \alpha ^2+i \sqrt{3} \alpha \text{kx}+\frac{3 i \alpha \text{ky}}{2}}+2 e^{9 \alpha ^2+\frac{1}{2} i \alpha \left(2 \sqrt{3} \text{kx}+\text{ky}\right)}+2 e^{9 \alpha ^2+\frac{1}{2} i \alpha \left(\sqrt{3} \text{kx}+4 \text{ky}\right)}+2 e^{12 \alpha ^2+\frac{1}{2} i \alpha \left(2 \sqrt{3} \text{kx}+\text{ky}\right)}+e^{15 \alpha ^2+\frac{1}{2} i \alpha \left(2 \sqrt{3} \text{kx}+\text{ky}\right)}+e^{\frac{1}{2} i \alpha \left(2 \sqrt{3} \text{kx}+\text{ky}\right)}+e^{\frac{1}{2} i \alpha \left(2 \sqrt{3} \text{kx}+3 \text{ky}\right)}+2 e^{\frac{1}{2} \alpha \left(24 \alpha +2 i \sqrt{3} \text{kx}+3 i \text{ky}\right)}+2 e^{\frac{1}{2} \alpha \left(24 \alpha +i \sqrt{3} \text{kx}+4 i \text{ky}\right)}+e^{\frac{1}{2} \alpha \left(30 \alpha +2 i \sqrt{3} \text{kx}+3 i \text{ky}\right)}+e^{\frac{1}{2} \alpha \left(30 \alpha +i \sqrt{3} \text{kx}+4 i \text{ky}\right)}+2 e^{9 \alpha ^2+\frac{i \alpha \text{ky}}{2}}+2 e^{9 \alpha ^2+\frac{3 i \alpha \text{ky}}{2}}+2 e^{12 \alpha ^2+\frac{i \alpha \text{ky}}{2}}+e^{15 \alpha ^2+\frac{i \alpha \text{ky}}{2}}+e^{15 \alpha ^2+\frac{3 i \alpha \text{ky}}{2}}+e^{\frac{i \alpha \text{ky}}{2}}+e^{\frac{3 i \alpha \text{ky}}{2}}+2 e^{\frac{3}{2} \alpha (8 \alpha +i \text{ky})}\right)$$

I have noticed this happening in many different cases. For some vectors the simplification works, for some others not. How should I proceed?

I feel kind of dumb, but I found the answer to my own question in Mathematica's documentation. The function ComplexExpand is what I should be using.
Simplify[ComplexExpand[Conjugate[wp[x0, y0]].H0.wp[x0, y0]]]

$$-2 e^{-16 \alpha ^2} t \left(2 \left(e^{9 \alpha ^2} \left(2 e^{3 \alpha ^2}+e^{6 \alpha ^2}+2\right)+1\right) \cos \left(\frac{1}{2} \sqrt{3} \alpha \text{kx}\right) \cos \left(\frac{\alpha \text{ky}}{2}\right)+e^{6 \alpha ^2} \left(2 e^{3 \alpha ^2}+2 e^{6 \alpha ^2}+e^{9 \alpha ^2}+2\right) \cos (\alpha \text{ky})\right)$$ which is real as expected!
Most importantly, this also correctly simplifies to real much more complex expressions than the ones mentioned here and it can get away with using Simplify instead of FullSimplify!