@luc ,
Quickly looking at the weak form, it seems that you may be missing some terms. If your trial field p is given by
p = p_r + i p_i
and likewise the test function q with the inner product (see here) defined as
<p, q> = \int_\Omega p\overline{q}d\mathbf{x}
then
\int_\Omega \nabla p \cdot \nabla \overline{q}\ d\mathbf{x} = \int_\Omega(\nabla p_r \cdot \nabla q_r + \nabla p_i\cdot\nabla q_i) d\mathbf{x} + i\int_\Omega\left( \nabla p_i\cdot\nabla q_r - \nabla p_r\cdot\nabla q_i\right)d\mathbf{x}
similarly for the inner product <p, q> you will have two more terms. Assuming homogeneous Neumann data on the boundary, I would therefore change to
a = inner(grad(pr), grad(qr)) + inner(grad(pi), grad(qi)) - k**2 * inner(pr, qr) - k**2 * inner(pi, qi)
a += inner(grad(pi), grad(qr)) - inner(grad(pr), grad(qi)) - k**2 * inner(pi, qr) + k**2 * inner(pr, qi)
a *= dx
L = inner(Constant(0), qr) + inner(Constant(0), qi)
L *= dx