Yes, it is true (if we replace $2^d$ to $2d$). If the range of $P-\lambda$ is not dense, there exists a bounded linear functional $\eta\in \ell^\infty=(\ell^1)^*$ which vanishes on this range. In other words, there exists a not identically zero bounded function $\eta$ on $\mathbb{Z}^d$ which satisfies $P^*\eta=P\eta=\lambda\eta$. Consider two cases.

1) $|\lambda|>1$. Note that since $\eta(x)=\lambda^{-1}(2d)^{-1}\sum_{y\sim x} \eta(y)$, the sum of coeffcients on the right is $\lambda^{-1}$. Iterating this, we see that each $\eta(x)$ is a linear combination of other values of $\eta$ with sum of coefficients $\lambda^{-n}$ for any $n$. Since $\lambda^{-n}$ tends to 0 and $\eta$ is bounded, we conclude that $\eta\equiv 0$, a contradiction.

2) $\lambda$ is not real. Denote $\mu=2d \lambda$. Take large $N$ and denote $\Omega_N=[-N,N]^d\cap \mathbb{Z}^d$. Consider the sum $S_n=\sum_{x\in \Omega_N,y\sim x} \eta(y)\overline{\eta(x)}$. It equals, if we replace $\sum_{y\sim x} \eta(y)$ to $\mu \eta(x)$, to $S_n=\mu\sum_{x\in \Omega_n}|\eta(x)|^2$. On the other hand, all internal edges have real contribution to this sum. For boundary edges use an estimate $2|\eta(y)\overline{\eta(x)}|\leqslant |\eta(x)|^2+|\eta(y)|^2$. Then denoting $\mu=u+iv$ and taking the imaginary part, we get $|v|\sum_{x\in \Omega_n} |\eta(x)|^2\leqslant 2\sum_{x\in \Omega_{n+1}\setminus \Omega_{n-1}} |\eta(x)|^2$. Now denote $a_n=\sum_{x\in \Omega_n} |\eta(x)|^2$. This sequence is positive for large $n$, increases and satisfies an inequality $a_{n+1}-a_{n-1}\geqslant ca_n$ for certain positive $c=|v|/2$. It follows by induction that $a_n\geqslant c_0 \rho^n$, where $\rho >1$ is a root of $\rho ^2-1=c\rho$. On the other hand, since $\eta$ is bounded, $a_n=O(n^d)$. A contradiction.