Summation of $\sum_{n=0}^{\infty}a^nq^{n^2}$

I am trying to find the result for the sum of the form

$\sum_{n=0}^{\infty}a^nq^{n^2}$.

The special case for $a=1$ is easily given by $\vartheta(0,q)$, where $\vartheta(z,q)$ is the third Jacobi Theta function. So, whatever the answer is, it must collapse to $\vartheta(0,q)$ for $a=1$.

I tried the approach:

$\sum_{n=0}^{\infty}\exp(2niz)q^{n^2}$, where $z=-i\frac{\ln(a)}{2}$, to make the summation look like $\vartheta(z,q)$. However, Jacobi Theta functions include summation from $-\infty$ to $\infty$, while I need it to start at 0.

I tried to use the Parseval's relation as well as several other methods but I cannot proceed any further.

Edit: Following Somos' comment on partial theta functions, this relation indeed is one of the many partial theta functions

$\Theta_P(a;q)=\sum_{i=0}^n a^nq^{n^2}$ (Eq. I)

However, this is no good for computation as I cannot express this in an existing function in python. I know this is an addition to the question, yet I am still looking for a closed form expression that can be used for simulations.

Another partial theta function is given in the form

$\Theta_P(a;q)=\Pi_{k=1}^n (1+aq^{k})(1+a^{-1}q^{k})$ (Eq. II)

If there exists a relationship between Eq. I and Eq. II, it may be possible to find the summation I am looking for in terms of, including but not limiting to theta functions and q-pochhammer.

Another possibly promising solution I pursued was Abel-Poisson summation method, yet of knowledge my mathematics did not allow me to pursue this to the end.

Thank you.

A possible approach could be to search an approximated estimate of the summation. Applying the Euler-Maclaurin formula to $f(n)=a^nq^{n^2}$ in the interval between $0$ and $N$, we have

$$\displaystyle \sum_0^N a^nq^{n^2}=\int_0^N a^nq^{n^2}+ \frac{a^Nq^{N^2}+1}{2} \\ + \sum_{k=1}^{\infty }\,{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(N)-f^{(2k-1)}(0)\right)$$

where $B_{2k}$ denotes the $2k^{th}$ Bernoulli number and $f^{(2k-1)}$ is the $(2k-1)^{th}$ derivative.

To estimate the integral, we can make the substitution $$n=\frac{t}{\sqrt{\log q}}-\frac{\log a}{2 \log q}$$ We have

$$\displaystyle a^nq^{n^2}=a^{t/\sqrt{\log q} -\log a/(2\log q) } q^ {[t/\sqrt{\log q} -\log a/(2\log q) ]^2} \\ = e^{-\log^2 a/(4 \log q)}e^{t^2} $$

Also, we have

$$t=\frac{\log a +2n \log q}{2\sqrt{\log q}}$$ $$dt=\sqrt{\log q } \, dn$$ $$dn=\frac{1}{\sqrt{\log q}} dt$$

Then our integral can be rewritten as

$$\displaystyle \int_0^N a^nq^{n^2}dn\\= \frac{e^{-\log^2 a/(4 \log q)}}{\sqrt{\log q}} \int_{\frac{\log a}{2\sqrt{\log q}}}^{\frac{\log a +2N\log q}{2\sqrt{\log q}}} e^{t^2}$$

Setting by simplicity $$K_1=\log a/(2\sqrt{\log q})$$ and $$K_2=(\log a +2N\log q)/(2\sqrt{\log q})$$ we can further rewrite the integral as

$$\displaystyle \int_0^N a^nq^{n^2}dn =\frac{e^{-K_1^2 }}{\sqrt{\log q}} \int_{K_1}^{K_2} q^{t^2}\\ =\frac{e^{-K_1^2}}{\sqrt{\log q}} \left(\int_{0}^{K_2} q^{t^2} -\int_{0}^{K_1} q^{t^2} \right)\\ = \frac{\sqrt{\pi}\,e^{-K_1^2 }}{2\sqrt{\log q}} \bigg[\operatorname{erfi}(K_2) - \operatorname{erfi}(K_1)\bigg]\\$$

where $\operatorname{erfi}$ indicates the imaginary error function, defined as

$$\displaystyle \operatorname {erfi} (x)={\frac {2}{\sqrt {\pi }}}\int_{0}^{x}e^{t^{2}}\,dt$$

Alternatively, our integral could be written in terms of the Dawson function. In fact

$$\displaystyle \int_0^N a^nq^{n^2}dn\\ =\frac{e^{-\frac{\log^2 a}{4\log q} }}{\sqrt{\log q}} \left(\int_{0}^{K_2} q^{t^2} -\int_{0}^{K_1} q^{t^2} \right)\\ = \frac{e^{-K_2^2 +n \log a +2n^2\log q }}{\sqrt{\log q}} \int_{0}^{K_2} q^{t^2} \\ - \frac{e^{-K_1^2 }}{\sqrt{\log q}} \int_{0}^{K_1} q^{t^2} \\ =\frac{a^Nq^{N^2}}{\sqrt{\log q}} \, e^{-K_2^2} \int_{0}^{K_2} e^{t^2}\\ - \frac{e^{-K_1^1}}{ \sqrt{\log q}} \int_{0}^{K_1} e^{t^2} \\ =\frac{{a^Nq^{N^2} \,F_+(K_2) }- {F_+(K_1)}}{\sqrt{\log q}} $$

where $F_+(K_2)$ and $F_+(K_1)$ indicate the Dawson function, i.e. the one-sided Fourier–Laplace sine transform of the Gaussian function, usually expressed in the following two forms:

$${\displaystyle F_{+}(x)=e^{-x^{2}}\int _{0}^{x}e^{t^{2}}\,dt}$$

$${\displaystyle F_{-}(x)=e^{x^{2}}\int _{0}^{x}e^{-t^{2}}\,dt.\!}$$

Note that the Dawson function is linked to the error function and the imaginary error function by

$${\displaystyle \operatorname {erf} (x)={\frac {2}{\sqrt {\pi }}e^{-x^{2}}F_-(x)}}$$

$${\displaystyle \operatorname {erfi} (x)={\frac {2}{\sqrt {\pi }}e^{-x^{2}}F_+(x)}}$$

To proceed further, we have to note that the initial sum $\sum_0^\infty a^n q^{n^2}$ converges for $q<1$ and diverges for $q>1$ (the case $q=1$ reduces to the trivial sum $\sum_0^\infty a^n$, whose convergence or divergence depends on $a$). These two cases determine the presence of real or imaginary terms in the last equations. Since the OP asks for the result of infinite sums, I will focus only on the first case. However, the same following approach could be used to estimate partial sums in the second case.

If $q<1$, then $K_2$, $K_1$, and the term $\sqrt{ \log q}$ in the denominators of the RHS of the equations above are all imaginary. So if we write

$$K_2=\frac{\log a +2N\log q}{2\sqrt{\log q}}= \frac{\log a +2N\log q}{2\sqrt{-\log 1/q}}= -i\frac{\log a +2N\log q}{2\sqrt{\log 1/q}} =i\hat{K_2} $$

and similarly

$$K_1=\frac{\log a }{2\sqrt{\log q}}= \frac{\log a }{2\sqrt{-\log 1/q}}= -i\frac{\log a }{2\sqrt{\log 1/q}} =i\hat{K_1} $$

reminding that $\operatorname{erfi}(z i)= i\,\operatorname{erf}(z)$, we have

$$\displaystyle \int_0^N a^nq^{n^2}dn\\ = \frac{\sqrt{\pi} \, e^{\hat{K_1}^2} }{2\sqrt{\log q}}\left[ { \operatorname{erfi}(i\hat{K_2}) }- { \operatorname{erfi}(i\hat{K_1}) }\right] \\ = \frac{\sqrt{\pi} \, e^{\hat{K_1}^2} } {2\sqrt{\log (1/q)}}\left[ { \operatorname{erf}(\hat{K_2}) }- { \operatorname{erf}(\hat{K_1}) }\right] $$

Now we can consider that, as $N$ increases, the error function shows a fast convergence to $1$, satisfying $\operatorname{erf}(z)=1+O(e^{-z^2}z^{-1})$. Therefore, making the inverse substitution $\hat{K_1}= - (\log a )/(2\sqrt{\log 1/q}) $, we get

$$\displaystyle \int_0^\infty a^nq^{n^2}dn\\ = \frac{\sqrt{\pi} \, e^{ \frac{\log^2 a}{4\log 1/q}} } {2\sqrt{\log (1/q)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) }\right] $$

Turning back to the initial Euler-Maclaurin expansion, the degree of approximation for the sum $\sum_0^\infty a^nq^{n^2}$ depends on the number of terms considered after the integral. Among these, it is not difficult to note that the leading terms result from the values of the function and the derivatives in $n=0$. For most values of $a$ and $q$, a good level of approximation can be obtained considering these terms up to the third derivative. Because $$f^{(1)}(0)=\log a$$ $$f^{(3)}(0)=\log^2 a + 6\log a \log q$$

and $$B_2=\frac{1}{6}$$ $$B_4=-\frac{1}{30}$$

we have

$$\displaystyle \sum_0^N a^nq^{n^2}\approx \frac{\sqrt{\pi} \, e^{ \frac{\log^2 a}{4\log 1/q}} } {2\sqrt{\log (1/q)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) }\right] + \frac{1}{2} -\frac{\log a}{12}+\frac{\log a(\log^2 a +6 \log q)}{720} $$

For example, setting $a=1/2$ and $q=3/4$, we have

$$\sum_0^\infty a^nq^{n^2}=1.46413758...$$

as shown by WA here. Our approximated formula yields

$$\displaystyle \sum_0^N a^nq^{n^2}\approx \frac{\sqrt{\pi} \, e^{ \frac{\log^2 1/2}{4\log 4/3}} } {2\sqrt{\log (4/3)}}\left[ 1- { \operatorname{erf}\left(-\frac{ \log 1/2 }{2\sqrt{\log 4/3}}\right) }\right] + \frac{1}{2} -\frac{\log 1/2}{12}+\frac{\log 1/2 (\log^2 1/2 +6 \log 4/3)}{720} \\ =1.46407819...$$

as confirmed by WA here. The relative error, expressed in percentage, is $-0.004\%$.

Including further terms of the Euler-Maclaurin espansion, the approximation can be improved. For example, considering the terms up to the fifth derivative, taking into account that

$$f^{(5)}(0)\\=\log^5 a + 20 \log^3 a \log q \\ + 60 \log a \log^2 q $$

and $$B_6=-\frac{1}{42}$$

the approximated formula above would include the additional term

$$-\frac{1}{30240} \bigg[\log^5 (1/2) + 20 \log^3 (1/2) \log (3/4) + 60 \log(1/2) \log^2 (3/4) \bigg]=0.00005574... $$

as shown here. In this way, the estimated sum would become

$$=1.46407819 + 0.00005574 = 1.46413393 $$

with a relative error of $-0.00025\%$.

Lastly, to facilitate the computation of the approximated formula (e.g. in Python), the value of the $\operatorname{erf}$ function can be calculated using several numerical approximations. These allow to calculate the value of $$\operatorname{erf}\left(-\frac{ \log a }{2\sqrt{\log 1/q}}\right) $$ with good precision using elementary functions.

Sequence $\{a^n q^{n^2}\} = \{b^n\},\ $ where $\,\{b_n=aq^n\},\ $ converges to zero if and only if $|q|<1.$ If $aq>1,\ $ then it has increasing beginning part (see Wolfram Alpha plot for $q=3/4, a=1/2,1,2,4\ $ with the unit of the vertical axis in the decimal logarithms.).

The given series converges absolutely, if and only if $|q|<1.$

Known Laplace transform $$\mathcal L\left(e^{\Large -\,^{t^2}\!/\!_{4k^2}}\right)=k\sqrt\pi\, e^{k^2s^2}\operatorname{erfc}ks\,(k>0)\tag1$$ allows to apply the Euler-Maclaurin formlula with the same nodes, as the given sum.

On the other hand, the direct summation gives the better accuracy and calculation complexity.

At the same time, since $$\theta_3(z,q)=\sum\limits_{n=-\infty}^\infty q^{\large n^2}e^{\large 2inz},\tag2$$ then $$S=\sum\limits_{n=0}^{\infty} a^n q^{\large n^2} =\sum\limits_{n=-N}^{\infty} a^{n+N} q^{\large (n+N)^2} = \lim\limits_{N\to\infty}a^{\large N} q^{\large N^2}\sum\limits_{n=-N}^{\infty} q^{\large n^2}(aq^{2N})^{\large n}\\ = \lim\limits_{N\to\infty}a^{\large N} q^{\large N^2}\sum\limits_{n=-N}^{\infty} q^{\large n^2}e^{\large i n (-i(\ln a + 2N\ln q))},$$ $$S= \lim\limits_{N\to\infty}\left(aq^{\large N}\right)^{\large N} \theta_3\left(-i\dfrac{\ln a + 2N\ln q}2,q\right).\tag3$$

The worst case of summation is $aq^L = 1$ for the large values of $L.$ Then the greatest value of $M\approx a^nq^{n^2}$ achieves if $$\dfrac{a^{n+1}q^{(n+1)^2}}{a^nq^{n^2}} = 1, \quad n\approx \dfrac{-\log q}{2\log a}\approx\dfrac L2,\quad M\approx a^Lq^{L^2}\approx a^{\Large^L\!/\!_4},\quad a\approx \sqrt[\Large L]{M^4}.\tag4$$ If $a=1.335,\quad q=0.99887$, then $M\approx 10^8,\ L\approx 255.$

In this case, formula $(3)$ with $N=0$ gives $$S\approx 54764\,68327.05684,$$ and direct summation gives $$S = 54764\,683241.58768\dots$$