Electromagnetic wave function of a photon

Question

I wanted to find the electromagnetic wave function of a photon, or at least what it might be. An interesting thing about photons is no matter how far they travel, their field stays concentrated instead of spreading out. Usually when one looks at electromagnetic waves being create by a source, they diffuse over distance. I started out my calculation with Maxwell’s equations and used them to find the electromagnetic wave equation in a homogeneous isotropic medium. $$\begin{matrix}\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial t}&\nabla\times\vec{H}=\frac{\partial\vec{D}}{\partial t}\\\end{matrix}$$ $$\nabla\times\nabla\times\vec{E}=\nabla\left(\nabla\cdot\vec{E}\right)-\nabla^2\vec{E}=\nabla\left(0\right)-\nabla^2\vec{E}=-\nabla^2\vec{E}=-\frac{\partial\left(\nabla\times\vec{B}\right)}{\partial t}=-\frac{\partial\left(\nabla\times\mu\vec{H}\right)}{\partial t}=-\mu\frac{\partial}{\partial t}\left(\nabla\times\vec{H}\right)=-\mu\frac{\partial}{\partial t}\left(\frac{\partial\vec{D}}{\partial t}\right)=-\mu\left(\frac{\partial^2\vec{D}}{\partial t^2}\right)=-\varepsilon\mu\left(\frac{\partial^2\vec{E}}{\partial t^2}\right)=-\frac{1}{c^2}\left(\frac{\partial^2\vec{E}}{\partial t^2}\right)$$ $$\nabla^2\vec{E}=\frac{1}{c^2}\left(\frac{\partial^2\vec{E}}{\partial t^2}\right)$$ I decided the wave would travel in the +z direction. As such, the field will be of form $$\vec{E}=\vec{A}\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\vec{B}\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}$$ $\vec A$ and $\vec B$ are functions of x and y only. Now to solve the equation $$\nabla^2\vec{E}=\left[\begin{matrix}\nabla^2E_x\\\nabla^2E_y\\\nabla^2E_z\\\end{matrix}\right]=-\frac{\omega^2}{c^2}\vec{E}+\left(\nabla^2\vec{A}\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\nabla^2\vec{B}\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}$$ $$\frac{\partial^2\vec{E}}{\partial t^2}=-\omega^2\vec{E}$$ $$\nabla^2\vec{E}=-\frac{\omega^2}{c^2}\vec{E}+\left(\nabla^2\vec{A}\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\nabla^2\vec{B}\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}=-\frac{\omega^2}{c^2}\vec{E}=\frac{1}{c^2}\left(\frac{\partial^2\vec{E}}{\partial t^2}\right)$$ $$0=\left(\nabla^2\vec{A}\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\nabla^2\vec{B}\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}$$ $$0=\nabla^2\vec{A}=\nabla^2\vec{B}$$ As such, the x, y, and z components of $\vec A$ and $\vec B$ need to be harmonic functions. We have another constraint though, the divergence must be 0. $$\nabla\cdot\vec{E}=\left(\frac{\partial A_x}{\partial x}\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\frac{\partial B_x}{\partial x}\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\frac{\partial A_y}{\partial y}\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\frac{\partial B_y}{\partial y}\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\frac{\omega}{c}\left(A_z\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}-B_z\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}\right)=\left(\left(\frac{\partial A_x}{\partial x}\right)+\left(\frac{\partial A_y}{\partial y}\right)-\frac{\omega}{c}B_z\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\left(\frac{\partial B_x}{\partial x}\right)+\left(\frac{\partial B_y}{\partial y}\right)+\frac{\omega}{c}A_z\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}=\left(\nabla\cdot\vec{A}-\frac{\omega}{c}B_z\right)\sin{\left(\frac{\omega}{c}\left(z-ct\right)\right)}+\left(\nabla\cdot\vec{B}+\frac{\omega}{c}A_z\right)\cos{\left(\frac{\omega}{c}\left(z-ct\right)\right)}=0$$ $$\begin{matrix}B_z=\frac{c}{\omega}\left(\nabla\cdot\vec{A}\right)&A_z=-\frac{c}{\omega}\left(\nabla\cdot\vec{B}\right)\\\end{matrix}$$ Now to find what limitations that imposes on $\vec A$ and $\vec B$. We’ll keep in mind that they are harmonic. $$\begin{matrix}\nabla^2B_z=\frac{c}{\omega}\nabla^2\left(\nabla\cdot\vec{A}\right)=\frac{c}{\omega}\nabla\cdot\left(\nabla^2\vec{A}\right)=\frac{c}{\omega}\nabla\cdot\left(0\right)=0\\\nabla^2A_z=-\frac{c}{\omega}\nabla^2\left(\nabla\cdot\vec{B}\right)=-\frac{c}{\omega}\nabla\cdot\left(\nabla^2\vec{B}\right)=-\frac{c}{\omega}\nabla\cdot\left(0\right)=0\\\end{matrix}$$ So no more limitations are imposed on them. As such the x and y components of $\vec A$ and $\vec B$ can be expressed as a linear combination of 2D harmonic functions. Their z components are calculated from there. In order for the wavefunction to give a result similar to photons, it must be finite everywhere, real, go to 0 away from the origin, and not be 0 everywhere. Unfortunately, harmonic functions reach their maximum and minimum values at the boundary. Since we defined the boundary to be 0 at infinity, it must b 0 everywhere. As such, we get a solution of 0. If we were to allow point discontinuities, such as $\ln{\left(r\right)}$ and $\frac{\sin{\left(n\varphi\right)}}{r^n}$ then we’d just get the function to diverge to $\pm\infty$, which would make the required energy per length spike to infinity. It seems an electromagnetic wave which meets all these criteria in a homogeneous isotropic medium doesn’t exist. I’m starting to wonder if photons don’t actually exist in homogeneous isotropic media. Maybe their energy warps the geometry of spacetime in and around them in a way which allows their electromagnetic energy to remain concentrated in a small volume. Any input/advice is appreciated.