Theta-function in the lower half-plane Standard theta function 
$$\vartheta(q)=\sum_{n=-\infty}^\infty q^{n^2} \qquad\qquad(1)$$
has a natural boundary of analyticity at $|q|=1$. This means that it can not be used to regularize expressions of the type $\sum_{n=-\infty}^\infty q^{n^2}$ with $|q|>1$. I however encounter this type of sums in some heuristic computations and would like to have some way to make sense of it (there is a finite answer available for comparison). In a similar context heuristic gives $\sum_{n} n^k$ and the correct answer agrees with the zeta-function regularization $\zeta(-k)$. Is there any similar go-to way to think about (1) but with $|q|>1$?
 A: If you are still active, I would be very interested to see the finite values you obtain in a physical context. For now, however, I will present two equal approaches we could take to regularize the sum.
Alternating Case
Ultimately, many series diverge because summation leaves some type of $x^N$ artifact remaining. The goal is to remove that artifact so that we end up with a sum that converges everywhere.
As an example, in the case of the geometric series the partial sums are $\sum_{n=0}^N x^n = \frac{1 - x^{N+1}}{1-x}$. If we let $S(N) = \sum_{n=0}^N x^n$, and define $S$ by the process $S(N+1) = S(N) + x^{N+1}$, then we realize we can take the sum both forwards or backwards. To take the sum backwards, write $S(N-1) = S(N) - x^N$, which gives rise to the sum $-\sum_{n=-1}^{-\infty} x^n = -\sum_{n=1}^\infty x^{-n} = \frac{1}{1-x}$ for $|x|>1$. This trick simply allows us to switch the artifact from $x^N$ to $x^{-N}$, which switches the convergent range from $|x|<1$ to $|x|>1$. Generally, we can use this trick to find representations for sums with only a single $x^N$ artifact. For instance $$\sum_{n=0}^{\infty}\sin\left(n+5\right)x^{n} \to -\sum_{n=1}^{\infty}\sin\left(-n+5\right)x^{-n}$$
or
$$\sum_{n=0}^{\infty}ne^{n}\cos\left(n-2\right)x^{n} \to -\sum_{n=1}^{\infty}\left(-n\right)e^{-n}\cos\left(-n-2\right)x^{-n}$$
The example is not too important, the main point is to say that finding a way to remove that $x^N$ artifact allows us to continue/regularize a function.
If we have an alternating function there's a nice approach to make this happen. We take $\sum_{n=0}^N (-1)^n f(n) x^n$ and transform it into
$$\sum_{n=0}^{N/2} f(2n) x^{2n} - \sum_{n=0}^{\frac{N-1}{2}} f(2n+1) x^{2n+1}$$
Next, we approximate $f(n)$ by a taylor series, so that we can define the sum for non-integer values. In the case of the sum $\sum_{n=0}^\infty (-1)^n q^{n^2}$ we get
$$\sum_{k=0}^{N/2} (-1)^n x^n q^{4n^2} =  \sum_{n=0}^{N/2} \sum_{k=0}^\infty (-1)^n x^n  \frac{(4\ln(q) n^2)^{k}}{k!} = \sum_{k=0}^\infty \frac{(4\ln(q))^{k}}{k!}\sum_{n=0}^{N/2} n^{2k} x^n (-1)^n =   \sum_{k=0}^\infty \frac{(4\ln(q))^{k}}{k!} \left[\left(x \frac{d}{dx}\right)^{2k} \frac{1-x^{\frac{N}{2}+1}}{1-x}\right]  $$
So, now we take $N=1$ and get the final sum as
$$\sum_{k=0}^\infty \frac{(4\ln(q))^{k}}{k!} \left[\left(x \frac{d}{dx}\right)^{2k} \frac{1-x^{\frac{1}{2}+1}}{1-x}\right] - f(1) x$$
This doesn't converge for $x=1$, and unfortunately becomes numerically unstable at around $|x-1| \approx .05$. However, I like this approach over some of the other ways of regularizing series because it converges even when $|q|<1$, which guarantees that this method will continuously extend all derivatives. When $x=1/2$, the graph looks like this:  where blue is $\sum q^{n^2} (-1/2)^n$, and orange is the continuation.
To extend this to the general complex case, we instead look at the roots of unity. If we are looking at a ray cast in the direction $\theta = \frac{p}{q} \tau$, then $e^{ i \theta}$ is periodic with period q. This leads us to the expansion:
$$\sum_{n=0}^N e^{\frac{p}{q} \tau i n } f(n) = \sum_{k=0}^{q-1} \sum_{n=0}^{(N-k)/q} e^{\frac{p}{q} \tau i k} f(qn+k)$$
Integral Approach
I originally saw this approach from Jorge Zuniga here, and I recommend you read his answer first. However, there are a few things I would add. The general way I view the approach is transforming
$\sum_{n=0}^\infty f(n)$ into $\sum_{n=0}^\infty \mathcal{M}\left\{\mathcal{M}^{-1}\left\{f(n)\right\}\right\} = \mathcal{M}\left\{\sum_{n=0}^\infty \mathcal{M}^{-1}\left\{f(n)\right\}\right\}$. Borel's method can be viewed as doing almost this same transformation, except with the Laplace instead of the mellin transformation. In particular, if we look at $\sum \mathcal{L}\left\{\mathcal{L}^{-1}\{f(n)(1/x)^n\}\right\}$ we get
$$\sum \mathcal{L}^{-1}\{f(n)(1/x)^n\} = \sum \frac{f(n)x^n}{n!}$$
and then
$$ \sum \mathcal{L}\left\{\frac{f(n)x^n}{n!} \right\} = \sum \int_0^\infty e^{-t} \frac{f(n)x^n}{n!} dt = \int_0^\infty e^{-t} \sum  \frac{f(n)x^n}{n!} dt$$
It makes sense to apply the inverse Mellin transform to the function $q^{n^2} = e^{\log(q) n^2}$ because the inverse Mellin transform takes a path that goes up/down the imaginary direction. While $e^{\log(q) n^2}$ grows quickly on the real line, traveling up the imaginary axis it converges to 0 very fast.
Anyway, once you have applied the method, you obtain
$$\sum_{n=0}^\infty q^{n^2} =\frac{1}{2\sqrt{\pi}}\int_{0}^{\infty}\frac{1}{t\left(1+xt\right)}\cdot\frac{1}{q^{\frac{1}{4}\log_{q}\left(t\right)^{2}}\ln\left(q\right)^{\frac{1}{2}}}dt$$
which converges only when $q>1$ (if $q<1$, you could obtain an integral representation by applying the mellin transform first, and then the inverse mellin transform, instead of the other way around).
The graph gives exactly the same result as the other approach. The integral approach is shown as black dots: .
This approach makes some interesting properties easier to see. For convenience, write
$$\theta(w,x) = \sum_{n=0}^\infty w^{n^2} x^n$$
Then the theta function satisfies the following equation inside the circle
$$\frac{\theta(w,x)-1}{x} = w\theta(w,w^2 x)$$
And it also formally holds this relationship outside the circle. Moreover, the partial sums asymptotically hold this relationship for $w>1$. We could show that the integral extension satisfies this functional equation.
Update
I was randomly thinking back on this question, and it occured to me to try to follow your approach using residues. The residue approach is well suited for this problem, because $e^{z^2}$ goes to zero very rapidly on the imaginary direction, so we can get a contour that converges even when the sum diverges. Thus, we can write
$$\vartheta(w,x) = \sum_{n=0}^\infty w^{n^2} (-1)^n x^n = \frac{1}{2 \pi i} \int_{-1/2 - i \infty}^{1/2 + i \infty} w^{n^2} x^t \Gamma(-t)\Gamma(t+1) dt$$
This has the benefit of now converging also when plugging in complex values for both $w$ and $x$. Moreover, previously, I had an issue with uniqueness because I only had the one functional equation $-\frac{\vartheta(w,x)-1}{x} = w\vartheta(w,w^2x)$. However, I have recently discovered the differential equation that is supposed to pair with that one, which is
$$\vartheta_w(w,x) = x^4 w^{11} \vartheta_{xx}(w, w^4 x) + 4 x^3 w^7 \vartheta_x(w, w^4 x) - x^2 w^2 \vartheta_x(w,w^2 x) + 2 x^2 w^3 \vartheta(w,w^4 x) - x \vartheta(w,w^2 x)$$
These two equations together with some initial conditions should uniquely define $\vartheta$, since the first equation allows you to extend the definition of $\vartheta(w,x)$ for any fixed $w$ into all values of $x$. The second equation allows you to use knowledge of $\vartheta(w,x)$ at all values of $x$ to get $\vartheta(w,x)$ at different values of $w$. These two functional equations actually naturally lead to the natural boundary that emerges at $|w|=1$, since values where $|w|<1$ will always depend on values of $\vartheta$ where $|w|<1$, and similarly points where $|w|>1$ will always depend only on values of $\vartheta$ where $|w|>1$. So I imagine the initial conditions would have to specify the behaviour on the boundary (probably in a distributional sense) to uniquely define the function.
A: The answer to the question being asked for, on Generalized Analytical Continuation (GAC) of $\theta_i, i=1,2,3,4$ functions can be found in Addenda 2 of this MO discussion that extends Borel´s summation (by means of Nachbin's Theorem) when nome $q$ holds $|q|>1$.
If you look such addenda with statistical eyes you realize that such formulation is not unique. In fact, the Integral for $s(\zeta,q)$ and their derivatives comes from an infinite sum of moments of a log-normal distribution that, as it is known, it cannot be uniquely determined by their moments.
So different measures can produce the same $s(\zeta,q)$ series. Thus a unique base $\theta_i$ function for $|q|<1$ can be analytically continued for $|q|>1$ to multiple (a truly infinite family of) integral functions. I rather prefer to talk about exo-analyticity since there is no "continuation" beyond the natural boundaries (the unit-$q$ disk in this case).
To illustrate this we use a classical example from Stieltjes
Stieltjes, T. -J., Investigations on continued fractions, Toulouse Ann. VIII, J1-J122 (1894); ibid. IX, A1-A47 (1895); C. R. CXVIII, 1401-1403 (1894). ZBL25.0326.01.
fitted to our notation with $\zeta\in\mathbb{C}\backslash(-\infty,0),\space q\in\mathbb{C}\wedge|q|>1$ and $|p|<1$,$$s(\zeta,q)=\sum_{n=0}^{\infty}(-\zeta)^n q^{n^2}=\frac{1}{2\sqrt{\pi\log q}}\int_{0}^{\infty}\exp(-\frac{1}{4}\frac{\log^2t}{\log q})\cdot\frac{1+p\sin(\pi\log_qt)}{1+t\zeta}\cdot\frac{dt}{t}$$
$p$ is an additional dummy parameter that does not contribute to the sum since all moments with the $\sin$ term vanish. Anyway the measure used in the mentioned MO discussion is a simple one (it corresponds to an unimodal distribution with $p=0$) that comes from basic Mellin Transforms. It allows to define a single equivalent exo-function for every $\theta_i$ Fourier $q$-series.
I leave a couple of references here
López-García, M., Characterization of solutions to the log-normal moment problem, Theory Probab. Appl. 55, No. 2, 303-307 (2011); translation from Teor. Veroyatn. Primen. 55, No. 2, 387-391 (2010). ZBL1225.60024.
Leipnik, R., The lognormal distribution and strong non-uniqueness of the moment problem, Theory Probab. Appl. 26, 850-851 (1982). ZBL0488.60024.
