The question is a good one and has been asked also by [Ehrenfest (1932): "Einige die Quantenmechanik betreffende Erkundigungsfragen"][1]. The answer was given by [Pauli (1933): "Einige die Quantenmechanik betreffenden Erkundigungsfragen"][2]. Unfortunately I'm not aware of a english translation of these two publications. However one can find a slightly different form of the answer also in [Pauli's book "General Principles of Quantum Mechanics" p.16][3]. In that book Pauli writes

> a single real function is not sufficient in order to construct from
> wavefunctions of the form (3.1) a non-negative probability function
> that is constant in time when integrated over the whole space.

I will try to summarize his arguments here:

A wave packet to describe a single particle (basically deBroglie's idea) can be written generally like
$$
u(x,t) = \int U(k) e^{i(kx-\omega t)} dk = \int U(k) e^{ikx} dk \, e^{-i\omega t}
$$
where $U(k)$ is the Fourier transform of $u(x,0)$. The complex conjugate of this wave packet is
$$
u^*(x,t) = \int U^*(k) e^{-i(kx-\omega t)} dk = \int U^*(k) e^{-ikx} dk \, e^{i\omega t}
$$
One can also define such wave packets in Electrodynamics. But in Quantum mechanics we have an additional condition, namely that the probability $P(x,t)$ to find a particle must always be positive and the total probability to find a single particle somewhere must be one, so
$$
P(x,t) \geq 1 \\
\int P(x,t)\, dx = 1
$$
Pauli argues that the simplest ansatz to construct such a function $P(x,t)$ from $u(x,t)$ is a [definite quadratic form][4] from the functions $u$ and $u^*$, that means
$$
P(x,t) = a u^2 + b {u^*}^2 + c u u^*
$$
Now from the form of $u(x,t)$ and $u^*(x,t)$ we see that
$$
u^2 \sim e^{-2i\omega t}\  \text{and}\  {u^*}^2 \sim e^{2i\omega t}
$$
and an integral over space over these two functions can never be time independent. So the constants $a$ and $b$ must be zero in the ansatz for $P(x,t)$. Only the product of a wave packet and it's complex conjugate will yield a time independent total probability:
$$
1 = \int P(x,t)\, dx = \int uu^*\, dx = \iiint U(k)U^*(k') e^{i(kx-k'x)} \, e^{-i\omega t} e^{i\omega t} dk dx \\
= \iint U(k)U^*(k') \delta(k-k') dk dk' = \int \left|{U(k)}\right|^2 dk = \int P(k)\, dk
$$
Since the product $uu^* = Re[u]^2 + Im[u]^2$ it follows that - as Pauli said - in order to compute a meaningful probability from wave packets of the form $u(x,t)$ one needs the real and the imaginary part of $u(x,t)$ and the wave function in quantum mechanics must be complex.


  [1]: https://doi.org/10.1007/BF01331295
  [2]: https://link.springer.com/article/10.1007%2FBF01335695
  [3]: https://books.google.de/books?id=hVjsCAAAQBAJ&printsec=frontcover&hl=de#v=onepage&q&f=false
  [4]: https://en.wikipedia.org/wiki/Definite_quadratic_form