Consider the twisted Eisenstein series $$ \begin{align} E_{k \ge 1}\left[\begin{matrix} \phi \\ \theta \end{matrix}\right] := & \ - \frac{B_k(\lambda)}{k!} \\ & \ + \frac{1}{(k-1)!}\sum_{r \ge 0}' \frac{(r + \lambda)^{k - 1}\theta^{-1} q^{r + \lambda}}{1 - \theta^{-1}q^{r + \lambda}} + \frac{(-1)^k}{(k-1)!}\sum_{r \ge 1} \frac{(r - \lambda)^{k - 1}\theta q^{r - \lambda}}{1 - \theta q^{r - \lambda}} \ , \end{align} $$ where $\phi = e^{2\pi i \lambda}, q = e^{2\pi i \tau}$. We know that they all can be written in Fourier series, for example, $$ E_1 \begin{bmatrix} -1 \\ z \end{bmatrix} = \frac{1}{2i}\sum_{n \in \mathbb{Z}, n \ne 0} \frac{1}{\sin n \pi \tau} e^{2\pi i n \mathfrak{z}} \ , \qquad E_2 \begin{bmatrix} 1 \\ z \end{bmatrix} = - \frac{1}{12} -\frac{1}{4}\sum_{n \in \mathbb{Z}, n \ne 0} \frac{1}{\sin^2 n \pi \tau} e^{2\pi i n \mathfrak{z}} \ , $$ where $z = e^{2\pi i \mathfrak{z}}$.
Recently I run into the following series $$ \sum_{\substack{m,n \in \mathbb{Z}\\m,n,m+n \ne 0}}\frac{1}{\sin m \pi \tau\sin n \pi \tau \sin (m + n) \pi \tau}e^{2\pi i (m+n) \mathfrak{z}} \ . $$ I wonder if this can be rewritten in terms of Eisenstein series? Or is it known to be some other special functions?
============Update============
With some guess work and help from Mathematica, I think the answer is the following $$ \sum_{\substack{m,n \in \mathbb{Z}\\m,n,m+n \ne 0}}\frac{1}{\sin m \pi \tau\sin n \pi \tau \sin (m + n) \pi \tau}e^{2\pi i (m+n) \mathfrak{z}} \\ = - \frac{\vartheta'_1(\mathfrak{z})\vartheta''_1(\mathfrak{z})}{2\pi^3\vartheta_1(\mathfrak{z})^2} + \frac{\vartheta'''_1 (\mathfrak{z})}{6\pi^3\vartheta_1'(\mathfrak{z})} + 4 E_2(\tau) \frac{\vartheta_1'(\mathfrak{z})}{\pi\vartheta_1(\mathfrak{z})} \ , $$ where $E_2(\tau)$ is the standard quasi-modular Eisenstein series. In terms of twisted Eisenstein series, it reads $$ - 8 i \Bigg( E_3 \begin{bmatrix} 1 \\ z \end{bmatrix} + E_2 \begin{bmatrix} 1 \\ z \end{bmatrix} E_1 \begin{bmatrix} 1 \\ z \end{bmatrix} - E_2(\tau)E_1 \begin{bmatrix} 1 \\ z \end{bmatrix} \Bigg) $$
However, a proof is still lacking.
The guessed result is a quasi-Jacobi form of modular weight-three. If this were known beforehand, one can always make an anzatz in terms of the twisted Eisenstein series and work out the relative coefficients using Mathematica. So a relevant question would be: how to argue that the series gives rise to such a quasi-Jacobi form?