The correct bilinear form for the Helmholtz operator in a lossless medium is
\int\left (\nabla p_r \cdot\nabla q_r + \nabla p_i \cdot\nabla q_i - k^2(p_r q_r + p_i q_i)\right ) d\mathbf{x} .
The coupling between real and imaginary components comes at the exit boundaries (impedance boundary condition). You will have surface terms of the form
k \int \left (p_r q_i - p_i q_r\right ) d\mathbf{S}
on these wave exit boundaries.
Your error is that you appear not to have included any exit boundaries. All the boundary conditions are reflecting (Dirichlet).