If you inspect the QED lagrangian, you'll find that the electromagnetic vector potential must have the units of [1] energy.
According to Peskin and Schroder's QFT text, the expansion of the electromagnetic field operator is
A_\mu(x)=\int\frac{d^3p}{(2\pi)^3}\frac{1}{\sqrt{2E_{\mathbf{p}}}}\sum_{r=0}^3\bigg(a_{\mathbf{p}}^r\epsilon_\mu^r(p)e^{-ip\cdot x}+a_{\mathbf{p}}^{r\dagger}\epsilon_\mu^{r*}(p)e^{ip\cdot x}\bigg).
3 factors of energy from d^3p, and -1/2 factor from 1/\sqrt{2E_{\mathbf{p}}}. If the polarization vectors, \epsilon_\mu^r, are unitless, then the ladder operators must carry units of 1½.
There is no problem at all. It is always possible to choose normalization in such away that the ladder operators become dimensionless;
Method(i): The simple box norlalization; expand the field as
A^{\mu}(x)=\sum_{r,\vec{k}}(2 L^{3}E)^{-1/2}\epsilon_{r}^{\mu}(\vec{k}) a_{r,\vec{k}} e^{-ikx} + ...
where the k summation is over all momenta allowed by the periodic boundary conditions (discretized momenta);
\vec{k}=\frac{2\pi}{L}\vec{n}
where n=0,\pm1,\pm2,....
Now, in the mass unit
[E]=[A]=-[L]=1, we have;
[A]=-1/2[E]-3/2[L]+[a], \Rightarrow[a]=0.
Now, do whatever you want with these dimensionless ladder operators.
Check Mandl & Shaw, they do QFT in a box.
Method(ii); If you want to work with the (continuum) integral expansion of the field, then you need to know the difference between
operator and
operator-valued distribution. In QFT, the solution of the field equation is an operator-valued distribution. There is a mathematically well-founded procedure in which the field distribution can be
smeared out with a "good" function in such away that it (the field) becomes a
bona fide operator linear with respect to this good function.
In QFT the relation
[a(p),a^{\dagger}(q)]=\delta^{3}(p-q)
seems to imply that the field consists of a noncountable set of oscillators labelled by the 3-momentum. However, it can be shown that we may use a countable set of oscillators instead of a continuum. This is achieved by using a set of good functions;
f_{k(i)}(\vec{p})\in L^2, i=1,2,..,
satisfying
1) orthonormality
\int (d^3p) f_{k(i)}(p) f_{k(j)}^{*}(p)=\delta_{ij}
2) completeness
\sum_{i} f_{k(i)}(p) f_{k(i)}^{*}(q)=\delta^{3}(\vec{p}-\vec{q})
Notice that [f]=-3/2 mass unit.
Now, we introduce (in your expansion) the smeared out operators (no longer distributions) by;
a_{k(i),r}=\int d^3pf_{k(i)}(p)a_{r}(p)
a_{k(i),r}^{\dagger}= \int d^3pf_{k(i)}^{*}(p)a_{r}^{\dagger}(p)
You can see that these operators are dimensionless, because
[a(p)]=-3/2 in your expansion.
Now, you can write your commutation relations in terms of these operators
[a_{k(i),r},a_{k(j),s}^{\dagger}]=-g_{rs}\delta_{ij}
and your(dimentionless) one photon state becomes
|1_{k,r}\rangle = \int d^3p f_{k}^{*}(p)a_{r}^{\dagger}(p)|0\rangle = a_{k,r}^{\dagger}|0\rangle,
and
\langle1_{k,r}|1_{k,r}\rangle = \langle0|a_{k,r}a_{k,r}^{\dagger}|0\rangle = -g_{rr}\langle0|0\rangle
(the norm of the state containing a scalar photon(r=0) is negative)
And finally, your number operator becomes
N=\sum_{k(i),r}a_{k(i),r}^{\dagger}a_{k(i),r} = \sum_{r} \int d^3p(-g_{rr})a_{r}^{\dagger}(\vec{p})a_{r}(\vec{p})
This is dimensionless as it should be.
good luck
regards
sam
