Fix a Schwartz function $g_0 : {\bf R} \to {\bf C}$, say $g_0(x) = e^{-\pi x^2}$ (which looks like a natural choice in this context).

For $x_1,x_2 \in \bf R$ define $g(x_1,x_2)$ by the absolutely convergent sum $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n) e^{2\pi i n x_1}. $$ Clearly $g(x_1,x_2) = g(x_1+1,x_2)$ for all $x_1,x_2$. Also $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i {n+1} x_1} = e^{2\pi i x_1} \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i n x_1} = e^{2\pi i x_1} g(x_1, x_2 + 1). $$ Finally define $$ f(x_1,x_2) = e^{\pi i x_1 x_2} g(x_1,x_2). $$ Then $f$ is a smooth function satisfies the required identities $$ f(x_1,x_2) = e^{-i\pi x_2} f(x_1+1,x_2), \quad f(x_1,x_2) = e^{i\pi x_1} f(x_1,x_2+1) \tag{$\star$} $$ for all real $x_1,x_2$.

This construction was surmised by working backwards, observing that if $f$ satisfies the required quasiperiodicity then $g$ is periodic in $x_1$ and thus has a Fourier series $\sum_{n=-\infty}^\infty g_n(x_2) e^{2\pi i n x_1}$, and then the identity $g(x_1,x_2) = e^{2\pi i x_1} g(x_1, x_2 + 1)$ yields $g_n(x_2+1) = g_{n+1}(x_2)$, whence $g_n(x_2) = g_0(x_2+n)$ for each $n$.

The $L^2$ norm of $g$ on the unit square is the square root of $$ \sum_{n=-\infty}^\infty \int_0^1 |g_n(x_2)|^2 \, dx_2 = \sum_{n=-\infty}^\infty \int_n^{n+1} |g_0(x_2)|^2 \, dx_2 = \int_{-\infty}^\infty |g_0(x_2)|^2 \, dx_2, $$ which is the $L^2$ norm of $g$ on $\bf R$. So we have an injection (and probably an isomorphism) of $L^2({\bf R}, {\bf C})$ into the Hilbert space of functions satisfying ($\star$) which takes Schwartz functions to smooth functions.

The construction breaks the (anti)symmetry between $x_1$ and $x_2$. However, if we make $f$ periodic in $x_2$ by multiplying by $e^{\pi i x_1 x_2}$ instead of $e^{-\pi i x_1 x_2}$, we end up with much the same formula but with $x_1,x_2$ switched and $g_0$ replaced by its Fourier transform! (This calculation is similar to the one that gives the Poisson summation formula.) In particular, for our choice $g_0(x) = e^{\pi x^2}$ it is the same function. This suggests using the Hermite functions (orthogonal polynomials times $e^{\pi x^2}$) for $g_0$ to obtain an orthogonal basis for the space of functions satisfying ($\star$).

Fix a Schwartz function $g_0 : {\bf R} \to {\bf C}$, say $g_0(x) = e^{-\pi x^2}$ (which looks like a natural choice in this context).

For $x_1,x_2 \in \bf R$ define $g(x_1,x_2)$ by the absolutely convergent sum $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n) e^{2\pi i n x_1}. $$ Clearly $g(x_1,x_2) = g(x_1+1,x_2)$ for all $x_1,x_2$. Also $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i {n+1} x_1} = e^{2\pi i x_1} \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i n x_1} = e^{2\pi i x_1} g(x_1, x_2 + 1). $$ Finally define $$ f(x_1,x_2) = e^{\pi i x_1 x_2} g(x_1,x_2). $$ Then $f$ is a smooth function satisfying the required identities $$ f(x_1,x_2) = e^{-i\pi x_2} f(x_1+1,x_2), \quad f(x_1,x_2) = e^{i\pi x_1} f(x_1,x_2+1) \tag{$\star$} $$ for all real $x_1,x_2$.

This construction was surmised by working backwards, observing that if $f$ satisfies the required quasiperiodicity then $g$ is periodic in $x_1$ and thus has a Fourier series $\sum_{n=-\infty}^\infty g_n(x_2) e^{2\pi i n x_1}$, and then the identity $g(x_1,x_2) = e^{2\pi i x_1} g(x_1, x_2 + 1)$ yields $g_n(x_2+1) = g_{n+1}(x_2)$, whence $g_n(x_2) = g_0(x_2+n)$ for each $n$.

The $L^2$ norm of $g$ on the unit square is the square root of $$ \sum_{n=-\infty}^\infty \int_0^1 |g_n(x_2)|^2 \, dx_2 = \sum_{n=-\infty}^\infty \int_n^{n+1} |g_0(x_2)|^2 \, dx_2 = \int_{-\infty}^\infty |g_0(x_2)|^2 \, dx_2, $$ which is the $L^2$ norm of $g$ on $\bf R$. So we have an injection (and probably an isomorphism) of $L^2({\bf R}, {\bf C})$ into the Hilbert space of functions satisfying ($\star$) which takes Schwartz functions to smooth functions.

The construction breaks the (anti)symmetry between $x_1$ and $x_2$. However, if we make $f$ periodic in $x_2$ by multiplying by $e^{\pi i x_1 x_2}$ instead of $e^{-\pi i x_1 x_2}$, we end up with much the same formula but with $x_1,x_2$ switched and $g_0$ replaced by its Fourier transform! (This calculation is similar to the one that gives the Poisson summation formula.) In particular, for our choice $g_0(x) = e^{-\pi x^2}$ it is the same function. This suggests using the Hermite functions (orthogonal polynomials times $e^{-\pi x^2}$) for $g_0$ to obtain an orthogonal basis for the space of functions satisfying ($\star$).

Fix a Schwartz function $g_0 : {\bf R} \to {\bf C}$, say $g_0(x) = e^{-\pi x^2}$ (which looks like a natural choice in this context).

For $x_1,x_2 \in \bf R$ define $g(x_1,x_2)$ by the absolutely convergent sum $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n) e^{2\pi i n x_1}. $$ Clearly $g(x_1,x_2) = g(x_1+1,x_2)$ for all $x_1,x_2$. Also $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i {n+1} x_1} = e^{2\pi i x_1} \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i n x_1} = e^{2\pi i x_1} g(x_1, x_2 + 1). $$ Finally define $$ f(x_1,x_2) = e^{\pi i x_1 x_2} g(x_1,x_2). $$ Then $f$ is a smooth function satisfies the required identities $$ f(x_1,x_2) = e^{-i\pi x_2} f(x_1+1,x_2), \quad f(x_1,x_2) = e^{i\pi x_1} f(x_1,x_2+1) $$ for all real $x_1,x_2$.

This construction was surmised by working backwards, observing that if $f$ satisfies the required quasiperiodicity then $g$ is periodic in $x_1$ and thus has a Fourier series $\sum_{n=-\infty}^\infty g_n(x_2) e^{2\pi i n x_1}$, and then the identity $g(x_1,x_2) = e^{2\pi i x_1} g(x_1, x_2 + 1)$ yields $g_n(x_2+1) = g_{n+1}(x_2)$, whence $g_n(x_2) = g_0(x_2+n)$ for each $n$.

Fix a Schwartz function $g_0 : {\bf R} \to {\bf C}$, say $g_0(x) = e^{-\pi x^2}$ (which looks like a natural choice in this context).

For $x_1,x_2 \in \bf R$ define $g(x_1,x_2)$ by the absolutely convergent sum $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n) e^{2\pi i n x_1}. $$ Clearly $g(x_1,x_2) = g(x_1+1,x_2)$ for all $x_1,x_2$. Also $$ g(x_1,x_2) = \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i {n+1} x_1} = e^{2\pi i x_1} \sum_{n=-\infty}^\infty g_0(x_2+n+1) e^{2\pi i n x_1} = e^{2\pi i x_1} g(x_1, x_2 + 1). $$ Finally define $$ f(x_1,x_2) = e^{\pi i x_1 x_2} g(x_1,x_2). $$ Then $f$ is a smooth function satisfies the required identities $$ f(x_1,x_2) = e^{-i\pi x_2} f(x_1+1,x_2), \quad f(x_1,x_2) = e^{i\pi x_1} f(x_1,x_2+1) \tag{$\star$} $$ for all real $x_1,x_2$.

This construction was surmised by working backwards, observing that if $f$ satisfies the required quasiperiodicity then $g$ is periodic in $x_1$ and thus has a Fourier series $\sum_{n=-\infty}^\infty g_n(x_2) e^{2\pi i n x_1}$, and then the identity $g(x_1,x_2) = e^{2\pi i x_1} g(x_1, x_2 + 1)$ yields $g_n(x_2+1) = g_{n+1}(x_2)$, whence $g_n(x_2) = g_0(x_2+n)$ for each $n$.

The $L^2$ norm of $g$ on the unit square is the square root of $$ \sum_{n=-\infty}^\infty \int_0^1 |g_n(x_2)|^2 \, dx_2 = \sum_{n=-\infty}^\infty \int_n^{n+1} |g_0(x_2)|^2 \, dx_2 = \int_{-\infty}^\infty |g_0(x_2)|^2 \, dx_2, $$ which is the $L^2$ norm of $g$ on $\bf R$. So we have an injection (and probably an isomorphism) of $L^2({\bf R}, {\bf C})$ into the Hilbert space of functions satisfying ($\star$) which takes Schwartz functions to smooth functions.

The construction breaks the (anti)symmetry between $x_1$ and $x_2$. However, if we make $f$ periodic in $x_2$ by multiplying by $e^{\pi i x_1 x_2}$ instead of $e^{-\pi i x_1 x_2}$, we end up with much the same formula but with $x_1,x_2$ switched and $g_0$ replaced by its Fourier transform! (This calculation is similar to the one that gives the Poisson summation formula.) In particular, for our choice $g_0(x) = e^{\pi x^2}$ it is the same function. This suggests using the Hermite functions (orthogonal polynomials times $e^{\pi x^2}$) for $g_0$ to obtain an orthogonal basis for the space of functions satisfying ($\star$).