I'll add background uniform charge to neutralise the point source so that the problem is well posed. I'll also rather consider the Green functions of $-\Delta$ which is positive. Actually, in dimensions $D\geq 2$, there is an ambiguity of the flat torus as you have many types of lattices up to a similarity (the Laplacian being invariant by similarities). Each equivalence class of lattices will correspondingly have its own fundamental solution which are not equivalent up to a simple change of variable.
In 1D, there is only one class, your domain is $\mathbb R/\mathbb Z$ with no loss of generality. You can calculate the fundamental solution directly:
$$
-\phi'' = \delta-1 \\
$$
or equivalently solve it in the unit interval $(0,1)$:
$$
\phi'' = 1
$$
with boundary condition:
$$
\phi'(0^+)-\phi'(0^-) = 1
$$
so:
$$
\phi = -\frac{x(1-x)}2
$$
if you set $\phi(0)=0$, which agrees with your formula on the line for $x\ll1$.
In 2D, you already have your answer.
For $N>3$, you'll start to get power law divergences at the origin, or in Fourier space, UV divergences. A general approach is to sum the Fourier representation of the Green's function.
I'll consider for simplicity the lattices of the form $\Lambda = \prod_{l=1}^DL_l\mathbb Z$, i.e. the domain is $\mathbb R^D/\Lambda$ and the dual lattice is $\Lambda^*=\prod_{l=1}^D\frac{2\pi}{L_l}\mathbb Z$. You problem is:
$$
-\Delta\phi = \delta-1
$$
In Fourier series:
$$
\phi(x) = \sum_{k\in\Lambda^*}\hat\phi(k)e^{ik\cdot x} \\
k^2\hat \phi(k) = 1 \\
\hat\phi = \sum_{k\in\Lambda^*-\{0\}}\frac{e^{ik\cdot x}}{k^2}
$$
The next steps are purely formal in order to reach a well defined expression which will define the appropriate regularised expression of the sum.
You can use the same trick as in the continuous case, you write:
$$
\frac1{k^2} = \int_0^{+\infty}e^{-k^2t}dt
$$
and permute sum and integral:
$$
\begin{align}
\hat \phi(x) &= \int_0^{+\infty}dt\sum_{k\in\Lambda^*-\{0\}}e^{ik\cdot x-k^2t} \\
&= \int_0^{+\infty}\left(-1+\prod_{l=1}^D\left(\sum_{n\in\mathbb Z}e^{i2\pi nx_l/L_l-(2\pi n/L_l)^2t}\right)\right)dt \\
&= \int_0^{+\infty}\left(-1+\prod_{l=1}^DK\left(\frac{x_l}{L_l},\frac t{L_l^2}\right)\right)dt \\
\end{align}
$$
I used the heat kernel of $\mathbb R/\mathbb Z$:
$$
K(x,t) = \sum_{n\in\mathbb Z}e^{i2\pi nx-(2\pi n)^2t}
$$
It is related to the Jacobi theta function by:
$$
K(x,t) = \theta(x,i4\pi t)
$$
This is an improvement as the sum is replaced by an integral. Actually, you can physically interpret this obscure mathematical trick. I am taking the time average of the heat kernel minus its average value (suggestively, $t$ represents time in this interpretation). Furthermore, using the asymptotic formula for $x\in(-1/2,1/2)$, $t\in\mathbb R_+^*$:
$$
\begin{align}
t&\to+\infty & K(x,t) &= 1+e^{-(2\pi)^2t}\cos(2\pi x)+... \\
t&\to0^+ & K(x,t) &\sim \frac{e^{-x^2/4t}}{\sqrt{4\pi t}}
\end{align}
$$
The integral converges absolutely for $x\in\mathbb R^D$ such that at most one $l\in\{1,...,D\}$ verifies $x_l\in L_l\mathbb Z$. For $x\in\Lambda$, you would expect a genuine divergence, but for the other cases, it means that the integral representation is insufficient. Thus the final expression can be taken as the regularised definition of the Green's function where it is well defined.
Note that the integral is defined everywhere in 1D, and you can check that it corresponds to the quadratic formula found at the beginning. Also note that for equal $x_l$ (along the diagonal), the expression even makes sense for non integer dimensions. The key to the simplification is the factorisation of the sum, which works due to the symmetries of the lattice and allows to revert the problem to 1D. You cannot expect this trick to work for a general lattice.
More rigorously, you can define the Fourier sum using zeta regularisation:
$$
\hat \phi(x,s) = \sum_{k\in\Lambda^*-\{0\}}\frac{e^{ik\cdot x}}{(k^2)^s}
$$
which converges absolutely for complex $s$ with $\Re s> \frac D2$. You can view it as the Green's function of $(-\Delta)^s$. By analytically continuation beyond the half complex plane and setting $s=1$, you recover the Green's function. To make the analytic continuation explicit, you can once again massage the sum into a single integral representation. The method is essentially the same as in Convergence of lattice sums and Madelung’s constant which can be viewed as a special case when $x_l=L_l/2$. You just need to generalize the previous trick to:
$$
\frac1{(k^2)^s} = \frac1{\Gamma(s)}\int_0^{+\infty}e^{-k^2t}t^{s-1}dt
$$
The steps are the same, and are now rigorous in the ROC of $s$:
$$
\phi(x,s) = \frac1{\Gamma(s)}\int_0^{+\infty}\left(-1+\prod_{l=1}^DK\left(\frac{x_l}{L_l},\frac t{L_l^2}\right)\right)t^{s-1}dt
$$
Unlike other functions like $\Gamma,\zeta$, you cannot extend the function by a Hankel contour since $K$ has a natural boundary at $\Re t=0$. However, the integral representation suffices to extend the function further depending on how many of the $x_l$ are in $L_l\mathbb Z$. If there are $m$, then you can extend the function to $\Re s>m/2$. In particular you can extend it to $s=1$ for $m=0,1$. In this case, it coincides with the previous formal derivation. For $m=D$, you expect a singularity, but for $2\leq m \leq D-1$, you would expect a finite value. For that, you can extend by continuity using the case $m=0$ for example.
The zeta regularisation would work also in the formula of the method of images. Using the free space Green's function:
$$
G(x) = \frac{\Gamma(D/2-1)}{4\pi^{D/2}r^{D-2}}
$$
Your Green's function is formally:
$$
\phi(x) = \sum_{\lambda\in\Lambda}\frac{\Gamma(D/2-1)}{4\pi^{D/2}|x-\lambda|^{D-2}}
$$
In a similar spirit, you can extend it to:
$$
\begin{align}
\phi(x,\sigma) &= \frac{\Gamma(D/2-1)}{4\pi^{D/2}}\sum_{\lambda\in\Lambda}\frac1{|x-\lambda|^{D-2+2\sigma}} \\
&= \frac{\Gamma(D/2-1)}{2\pi^{D/2}}\int_0^{+\infty} \left(\prod_{l=1}^D\sqrt{\frac{\pi}{u}}K(x_l/L_l,1/4u)\right)u^{D/2-1+\sigma}du \\
&= \frac{\Gamma(D/2-1)}4\int_0^{+\infty}\left( \prod_{l=1}^DK(x_l/L_l,1/4u)\right)u^{\sigma-1}du
\end{align}
$$
where you can recognise again $K$ from the method of images:
$$
\begin{align}
K(x,t) &= \sum_{n\in\mathbb Z}\frac{e^{-(x-n)^2/4t}}{\sqrt{4\pi t}} \end{align}
$$
This time the ROC is for $\Re\sigma > 1$ and the desired Green's function is obtained by setting $\sigma=0$.
You want to translate in real space the previous analytic continuation, then you can instead use the Green's function of $(-\Delta)^s$ in real space:
$$
G(x,s) =\frac{\Gamma(D/2-s)}{4^s\pi^{D/2}r^{D-2s}}
$$
so that the new regularised sum is rather:
$$
\phi(x,s) = \frac{\Gamma(D/2-s)}{4^s\pi^{D/2}}\sum_{\lambda\in\Lambda}\frac1{|x-\lambda|^{D-2s}} \\
$$
which is rigorously equal to the sum calculated using the Fourier series using the Poisson summation formula. While it is more physically transparent, it is a bit of a stretch at first to modify the prefactors as well. Notice that even for the apparently problematic case of $D=1$ in real space, the regularisation works. Actually, in $D=1$, if $x=0$, you recover the zeta function, so regularising your Green's function in this case is essentially equivalent to the infamous formula $\sum_{n=1}^\infty=-\frac1{12}$.
Hope this helps.
