Analytic continuation of the double sum $\sum_{n,m\ge0}x^ny^mt^{nm}$ Define function  $f(x,y,t)$ as the analytic continuation of the series
$$f(x,y,t)=\sum_{n,m\ge0}x^ny^mt^{nm}$$
This series definitely converges when all the arguments are small enough. I would like to understand the global properties of this function: presence of zeros, of poles or other singularities etc. In particular, I would like to locate the poles of this function w.r.t. to $y$ and find the corresponding residues.
Let us do a partial resummation
$$f(x,y,t)=\sum_{n\ge0}\frac{x^n}{1-yt^n}$$
which is possible when say $|y|,|t|<1$. 
Naively, from the above expression I expect that $f(x,y,t)$ has poles in $y$ for $y=t^{-n}, n\ge0$ with residues
$$\operatorname{Res}_{y=t^{-n}}f(x,y,t)=-x^n/t^n$$
However, there is a problem here for me. I will first formulate is loosely and then give more accurate description.
The problem
I. Loose formulation. Locations of poles of $f(x,y,t)$ are independent of $x$. Therefore, we don't expect them to change if $x$ changes. However, one can prove formula $f(x^{-1},y,t)=f(x,y^{-1},t)$ (modulo not-so-interesting terms). In the lhs we've only changed $x$, and hence we don't expect locations of poles to change. On the other hand in the rhs we've changed $y\to y^{-1}$ which causes poles to invert!
II. More accurate description.
 By formally operating with the series one can deduce three following properties


*

*$f(x,y,t)=f(y,x,t)$, obvious symmetry.

*$f(x,y^{-1},t)=-yf(xt^{-1},y,t^{-1})$ rule to invert an argument.

*$f(xt,yt,t)=(xyt)^{-1}(f(x,y,t)-(1-x)^{-1}-(1-y)^{-1}+1)$ scaling rule.
Using their combination one shows that
$$f(x^{-1},y,t)=f(x,y^{-1},t)+(1-y)^{-1}-(1-x)^{-1}$$
Again, as stated in the loose formulation, the poles in $y$ do not seem to agree between the lhs and rhs: in the lhs they are at the points $y=t^{-n}, n\ge0$ while at the rhs at the points $y=t^n, n\ge0$.
Summary
I understand that most likely I messed up with divergent series somewhere, but I can't exactly find where. I would be grateful if somebody explained me where, but  my main question is the following: what are locations and residues of poles for function $f(x,y,t)$ wrt $y$-variable? Do they depend on $x$? If yes, how?
Appendix
Here I present formal derivation of properties 2, 3.
$$f(x, y^{-1},t)=\sum_{n\ge0}\frac{x^n}{1-y^{-1}t^n}=\sum_{n\ge0}\frac{x^nt^{-n}y}{yt^{-n}-1}=-yf(xt^{-1}, y, t^{-1})$$
And 
$$f(xt, yt,t)=\sum_{n,m\ge0}{x^ny^mt^{nm+n+m}}=(xyt)^{-1}\sum_{n,m\ge0}x^{n+1}y^{m+1}t^{(n+1)(m+1)}=(xyt)^{-1}\left(\sum_{n,m\ge0}-\sum_{n\ge0,m=0}-\sum_{n=0,m\ge0}+\sum_{n=0,m=0}\right)x^{n}y^{m}t^{nm}=(xyt)^{-1}(f(x,y,t)-(1-x)^{-1}-(1-y)^{-1}+1)$$
On Dmitry Vaintrob's answer
As far as I can the the factorization procedure suggested in this answer is not actually valid since the two factors has non-overlapping region of convergence. So, I'm still looking for a solution.
 A: I think you can do the following (I haven't checked convergence carefully, so don't trust me too much). Write $a = \frac{m+n}{2}, b = \frac{m-n}{2}.$ Then $mn = a^2-b^2.$ Now (let's ignore indeterminacy in the square root for now), $\ x^m y^n = (xy)^a(x/y)^b$. The advantage of this sketchy manipulation is that your sum can be expressed, formally, as $$\sum_{a, b\in \frac{\mathbb{N}}{2}a\equiv b\text{ mod }1} \left((xy)^a t^{a^2}\right)\left((x/y)^b t^{-b^2}\right),$$
which is almost a product of two theta functions. 
A little more rigorously, if you choose square roots $\alpha^2 = xy, \beta^2 = x/y$ and a fourth root $\tau^4 = t$ and write $A = 2a, B = 2b$ then under a suitable convergence hypothesis,
$$f(x,y,t) = \frac{G(\alpha,\beta,\tau)+G(-\alpha, -\beta, \tau)}{2},$$ where $$G(\alpha, \beta, \tau) = \sum_{A, B}  \alpha^{A}\beta^{B}\tau^{A^2-B^2} = \left(\sum_{A} \alpha^A\tau^{A^2}\right)\left(\sum_{B}\beta^B\tau^{-B^2}\right),$$
which is actually a product of two theta functions. Of course this product cannot converge either when $|\tau|>1$ or $|\tau|<1,$ but convergence on the circle $|\tau| = 1$ with small $\alpha,\beta$ (corresponding to small $x$ compared to $y, y^{-1}$) should be enough to get your analytic continuation. From this it should not be very hard to deduce the poles. I wonder if such a formula can also be obtained from an appropriate version of the triple product identity. A little bit more context for how this problem came up could be helpful to see if this is the case.
Edit: As Lev Borisov points out, this only works when the sum is taken over $n\in \mathbb{Z}$ (equivalently, when studying $f(x,y,t) + f(x,y^{-1},t^{-1})$). As such, it is not an answer to the OP. -DV
A: I will give an alternative expression for $f(x,y,t)$ that may be helpful for others trying to solve this problem.  Initially suppose $|x|, |y|< 1$, and $0 < t < 1$.  Then $t^{mn} = \exp(-mn \log(1/t))$, and so by the Mellin formula for $e^{-x}$, we have
$$t^{mn} = \frac{1}{2\pi i} \int_{(\sigma)} \Gamma(s) (mn \log(1/t))^{-s} ds,$$
where $\sigma > 0$.  
By a simple rearrangement, we have
$$f(x,y,t) = -1 + \frac{1}{1-x} + \frac{1}{1-y} + \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} x^m y^n t^{mn}.$$
The infinite double sum above can then be expressed as
$$\frac{1}{2 \pi i} \int_{(\sigma)} \Gamma(s) (\log(1/t))^{-s} Li_s(x) Li_s(y) ds,$$
where $Li_s(z) = \sum_{n=1}^{\infty} x^n n^{-s}$ is the polylogarithm.
I haven't thought much about if this integral representation can lead to a meromorphic continuation.  The polylogarithm does have analytic continuation but the above integral may not converge absolutely if $|x| > 1$ or $|y| > 1$.
A: Suppose that $|x|<1$ and $|t|<1$. Then the power series $\sum_{n\geq 1}\frac {x^n}{1-yt^n}$, considered as a function of $y$, converges uniformly on compact sets disjoint from $y=t^{-k}$. Indeed, the denominators approach $1$ uniformly as $n\to\infty$ so you can bound by the geometric series.
Moreover, the same argument shows that in a neighborhood of $y=t^{-k}$ the sum $\sum_{n\geq 1, n\neq k}\frac {x^n}{1-yt^n}$ converges uniformly. This shows that the poles are indeed simple with residues $-\frac {x^n}{t^n}$.
If you try to analytically continue beyond $|x|,|t|<1$ region, then you may run into trouble. There is still convergence if $1<|x|<|t|$, so I think that in this regime the pole locations and residues are the same. Note that you will now have an essential singularity at $y=0$. 
It is less clear what happens if you try to analytically extend away from these regions. Moreover, I am not at all sure whether such extension would be possible, and whether there would be some monodromy. My guess is that you would need to avoid the $x=t^{-n}$ locations, and there may be monodromy associated to loops around these.
A: I though it is worth pointing out a closely related case which appears to be solvable. The following formula is due to paper http://arxiv.org/abs/math/9904126 (formula (11) with slight relabeling and rearrangement)
$$\sum_{n\in\mathbb{Z}}\frac{x^n}{1-yt^n}=\prod_{k\ge1}\frac{(1-xyt^{k-1})(1-x^{-1}y^{-1}t^{k})(1-t^k)^2}{(1-xt^{k-1})(1-x^{-1}t^{k})(1-yt^{k-1})(1-y^{-1}t^{k})}$$
Series in the l.h.s. converges for $|t|<|x|<1$ (probably in some other region too). If $|t|<|x|,|x|^{-1},|y|,|y|^{-1}$ it can be presented as a more symmetric double expansion
$$\sum_{n\in\mathbb{Z}}\frac{x^n}{1-yt^n}=-1+\frac1{1-x}+\frac1{1-y}+\sum_{n,m>0}\left(x^ny^m-x^{-n}y^{-m}\right)t^{nm}$$
