OK, here goes the "annulus".
Let's say that a domain $\Omega$ has a good approximation property at a point $w\in\Omega$ if for every $\rho>1$ there is a compact set $K=K(\rho)\subset \Omega$ and a constant $C=C(\rho)>0$ such that for every function $f$ analytic in $\Omega$ and every $m\ge 0$ there is a function $g$ analytic in $\Omega$ such that $g(z)-f(z)=O(|z-w|^{m+1})$ as $z\to w$ and $\sup_\Omega|g|\le C(\rho)\rho^m\sup_K|f|$.
Claim 1: The unit disk has a good approximation property at $0$.
Proof: Just let $K$ be the circle of radius $r\in(\rho^{-1},1)$ and take for $g$ the Taylor polynomial of $f$ of degree $m$.
Claim 2: Good approximation property is a conformal invariant on the Riemann sphere (with the understanding that if $w=\infty$, then the condition $g(z)-f(z)=O(|z-w|^{m+1})$ is replaced with $g(z)-f(z)=O(|z|^{-m-1})$, of course).
Proof: Obvious (map everything).
Claim 3: The Riemann sphere $\widehat{\mathbb C}$ with finitely many (reasonable) cuts $\Gamma_j\subset \mathbb C$ has good approximation property at $\infty$ (I assume $\Gamma_j$ are reasonable so that they are compact, pairwise disjoint, and $\widehat{\mathbb C}\setminus\Gamma_j$ is conformally equivalent to the unit disk for every $j$)
Proof: This is just a version of the simplest lemma on the separation of singularities. Fix $\rho>1$. First, each $\widehat{\mathbb C}\setminus\Gamma_j$ has good approximation property at $\infty$, so we can surround $\Gamma_j$ by some very close to it smooth contour $C_j$ that has the property of the compact $K$ in the definition of the good approximation property for $\widehat{\mathbb C}\setminus\Gamma_j$. Next, the function $f$ can be written as the sum $f(\infty)+\sum_j f_j$ where $f_j$ is the (clockwise) Cauchy integral of $f$ on any contour surrounding $\Gamma_j$ and having the evaluation point and other $\Gamma_k$ outside it. We can now readily estimate $f_j$ on $C_j$ by $\max_{C_j}|f|+|f(\infty)|+\sum_{k\ne j}\frac{\ell(C_k)}{2\pi\operatorname{dist}(C_j,C_k)}\max_{C_k}|f|$ by the triangle inequality, which yields $\max_{C_j}|f_j|\le C\max_K|f|$ where $K=\cup_j C_j$. Now, given $m$, just find an appropriate $g_j$ for each $f_j$ separately and put $g=\sum_j g_j$.
Assume now that $\Omega$ is a bounded domain with good approximation property at $0$ (we can always map the Riemann sphere with cuts and a fixed point conformally to this configuration). Let $H$ be the Hilbert space of functions analytic in $\Omega$ and square integrable with respect to the area measure (you can also introduce some reasonable weight, if you feel like it). Consider the subspaces $H_m=z^mH$ (those are just subspaces of functions in $H$ vanishing to order $m$ at the origin).
Clearly, we have $H=H_0\supset H_1\supset H_2\supset\dots$ and each $H_{m+1}$ is a closed subspace of $H_m$ of codimension $1$. Let $h_m$ be the function in $H_m$ orthogonal to $H_{m+1}$ of unit norm. Then, since $\cap_m H_m$ consists of analytic functions vanishing at $0$ with all derivatives, i.e., of just $0$, the functions $h_m$ form an orthonormal basis in $H$, so every function in $H$ can be written uniquely as its Fourier series in $h_m$, which converges in $H$ and, thereby, uniformly on compact subsets of $\Omega$ with all derivatives.
Notice now that the coefficient at $h_0$ is uniquely determined by $f(0)$ (all other terms vanish at $0$), the coefficient at $h_1$ is thus determined by $f(0)$ and $f'(0)$, and so on. This allows to write a formal Fourier decomposition into $h_m$ for any formal Taylor series at $0$. The only task is to show that the resulting series converges uniformly on compact sets in $\Omega$ if that Taylor series represents a function analytic in $\Omega$.
Note that due to the high order of $0$ of $h_m$ at the origin and the uniform bound for the $L^2$ norm with respect to the area measure, for every compact set $Q\subset\Omega$, there exists $B=B(Q)>0$ and $r=r(Q)<1$ such that $\max_Q|h_m|\le Br^m$. Now choose $\rho>1$ such that $\rho r<1$ and use the good approximation property. It yields a compact $K\subset\Omega$ and, for every $m$, a function $g$ that is bounded by $C\rho^m\max_K|f|$ and has the same first $m$ derivatives as $f$ at $0$. Then the formal Fourier coefficient of $f$ at $h_m$ is the same as of $g$, but $g\in H$ with norm bounded by the square root of the area of $\Omega$ times the uniform norm, so its Fourier coefficient is bounded by $C\rho^m$. This immediately implies the geometric convergence of the formal Fourier series of $f$ on $Q$.
Now it remains to note that $Q$ was arbitrary and that the sum of the formal Fourier series of $f$ (which we now know to converge) has the same derivatives at $0$ as $f$, so it must converge to $f$ itself.
Formally this doesn't work for general domains, but I still wonder if one can find a counterexample to this approach. Of course, one needs to replace the area measure by some weight very fast decreasing near the boundary, but it is not at all clear to me what the convergence properties of the corresponding formal Fourier series will be. The simplest case to consider is the punctured unit disk. There is no hope to get the good approximation property there (no matter what reasonable norm of $g$ you want to bound geometrically), so the last trick that yields a good bound for Fourier coefficients out of nothing won't work anymore, but can we still do something? Or can somebody come up with a completely different idea?
