It is essentially a discrete convolution between $\frac{1}{\cosh}$ and $\frac{1}{\cosh}$. It is interesting to point out that $\frac{1}{\cosh}$ is more or less a fixed point of the Fourier transform:
$$ \mathscr{F}\left(\frac{1}{\cosh x}\right)=\frac{\pi}{\cosh\frac{\pi s}{2}},$$
and its inverse Laplace transform is given by
$$ \mathscr{L}^{-1}\left(\frac{1}{\cosh x}\right) = 2\sum_{n\geq 0}(-1)^n \delta(s-2n-1)=f(s),$$
so:
$$ \frac{1}{\cosh\left(x\left(\frac{1}{2}-n\right)\right)}=\int_{0}^{+\infty}f(s) e^{-sx(n-1/2)}\,ds$$
and:
$$ \frac{1}{\cosh\left(x\left(n-\frac{1}{2}\right)\right)\cosh\left(x\left(n+p-\frac{1}{2}\right)\right)}=\int_{0}^{+\infty}\int_{0}^{+\infty} f(s)\,f(t)\,e^{-ptx}\,e^{-(s+t)x(n+1/2)}\,ds\,dt$$
so:
$$ \sum_{n\geq 0}\frac{1}{\cosh\left(x\left(n-\frac{1}{2}\right)\right)\cosh\left(x\left(n+p-\frac{1}{2}\right)\right)}=\int_{0}^{+\infty}\int_{0}^{+\infty}\frac{ 2\, f(s)\,f(t)\,e^{-ptx}}{\sinh\left((s+t)\frac{x}{2}\right)}\,ds\,dt$$
equals:
$$ \sum_{a\geq 0}\sum_{b\geq 0}\frac{8(-1)^{a+b} e^{-px(2b+1)}}{\sinh((a+b+1)x)}=e^{-px}\sum_{h=0}^{+\infty}\frac{8(-1)^h}{\sinh((h+1)x)}\cdot\frac{1-e^{-2(h+1)px}}{1-e^{-2px}} \tag{1}$$
that is a pretty fast convergent series.
Addendum. We may notice that the degenerate case $p=0$ is way easier to deal with.
Given the Weierstrass product for the $\cosh$ function, we have:
$$ \frac{1}{\cosh^2(x)}=\frac{d^2}{dx^2}\log\cosh(x) = -\sum_{m\geq 0}\left(\frac{1}{\left(x+\pi i\frac{2m+1}{2}\right)^2}+\frac{1}{\left(x-\pi i\frac{2m+1}{2}\right)^2}\right)$$
so we may recognize in the RHS a meromorphic function with equally-spaced and equally-behaving double poles along the imaginary axis. It follows that the meromorphic function defined through
$$ f(x)=\sum_{n\in\mathbb{Z}}\frac{1}{\cosh^2(x+n)} $$
is a constant multiple of a Weierstrass $\wp$ function, whose connection with elliptic integrals is well-known. It is interesting to point out that the RHS of $(1)$, as $p\to 0$, gives a fast-converging expansion for a Weierstrass $\wp$ function. On the other hand, the meromorphic function
$$ g(x)=\sum_{n\in\mathbb{Z}}\frac{1}{\cosh(x+n)\cosh(x+n+t)} $$
has only simple poles, due to:
$$\frac{1}{\cosh x}=\sum_{n\geq 0}(-1)^n\left(\frac{i}{x+(2m+1)\frac{\pi i}{2}}-\frac{i}{x-(2m+1)\frac{\pi i}{2}}\right)$$
so $g(x)$ is related with a Jacobi elliptic function; a useful identity (already exploited in a slightly disguised form through the inverse Laplace transform) is:
$$ \frac{1}{\cosh x} = \frac{2}{\pi}\int_{0}^{+\infty}\frac{\cos\left(\frac{2x}{\pi}u\right)}{\cosh u}\,du\tag{2} $$
that makes us aware of the following fact: a discrete convolution of $\frac{1}{\cosh}$ is the integral over $\mathbb{R}^+$ of $\frac{1}{\cosh u}$ times a discrete convolution of $\cos$, way easier to deal with. It turns out that the RHS of $(1)$ is related with $\frac{1}{\cosh(px)}$, too.
The first one is a pretty stupid one, but it's bugging me. In eq (1), how did you do the last step? I'm assuming its just by defining h=a+b, but then I've got issues with the limits of h...
Secondly, you don't happen to know if there's a closed expression for the series in (1)? I want to continue manipulating for quite a while and that would be very useful! As you pointed out, we know these series are related to the jacobi theta-functions, so my guess would be it can be expressed in terms of elliptic integrals or something..
– Lou May 09 '16 at 10:49