[proofplan] Write $D=\Omega_1\times\cdots\times\Omega_n$, where each $\Omega_j\subset\mathbb C$ is simply connected, and first solve the scalar case of $(0,q)$-forms. The main analytic input is the one-variable Cauchy-Green solution operator for $\partial u/\partial\bar z=f$, with smooth dependence on all passive variables. We prove the scalar case by induction on the number of active complex variables: remove the part not containing $d\bar z_k$, analyze the remaining part containing $d\bar z_k$, and solve that part either by induction or by the one-variable solver. Finally, decompose a general $(p,q)$-form by its holomorphic multi-index and solve each scalar anti-holomorphic component, keeping track of the sign introduced by wedging past $dz_I$. [/proofplan]

[step:Record the one-variable $\bar\partial$ solver with smooth parameters] [claim:One-variable Cauchy-Green solver with smooth parameters] Let $\Omega\subset\mathbb C$ be a simply connected planar domain, and let $M$ be either a point or a product of planar domains. For every smooth function \begin{align*} F:M\times\Omega&\to\mathbb C, \end{align*} there exists a smooth function \begin{align*} S_\Omega F:M\times\Omega&\to\mathbb C \end{align*} such that \begin{align*} \frac{\partial}{\partial\bar z}(S_\Omega F)(w,z)=F(w,z) \end{align*} for every $(w,z)\in M\times\Omega$. Moreover, the solution may be chosen so that differentiation in passive variables commutes with the construction: if $\bar\partial_M F=0$, then $\bar\partial_M(S_\Omega F)=0$. [/claim] [proof] First suppose that $F$ has compact support in the $\Omega$-variable, locally uniformly with respect to the passive variable $w\in M$. Let $\widetilde F:M\times\mathbb C\to\mathbb C$ denote the extension of $F$ by zero outside $M\times\Omega$. Define the Cauchy transform \begin{align*} C F:M\times\Omega&\to\mathbb C\\ (w,z)&\mapsto \frac{1}{\pi}\int_{\mathbb C}\frac{\widetilde F(w,\zeta)}{z-\zeta}\,d\mathcal L^2(\zeta). \end{align*} The kernel $\zeta\mapsto (z-\zeta)^{-1}$ is locally integrable with respect to $\mathcal L^2$ because, in polar coordinates centered at $z$, the singular contribution is bounded by \begin{align*} \int_0^r\int_{\{|\theta|=1\}}\frac{1}{s}\,s\,d\mathcal H^1(\theta)\,d\mathcal L^1(s) =2\pi r. \end{align*} Thus the integral is well-defined on compact subsets of $M\times\Omega$. By the Cauchy-Green Formula and the distributional identity \begin{align*} \frac{\partial}{\partial\bar z}\left(\frac{1}{\pi(z-\zeta)}\right)=\delta_\zeta, \end{align*} where $\delta_\zeta$ is the Dirac measure at $\zeta$, we obtain \begin{align*} \frac{\partial}{\partial\bar z}(C F)(w,z)=F(w,z). \end{align*} Differentiation with respect to passive variables passes under the integral by dominated convergence on compact subsets, so \begin{align*} \bar\partial_M(CF)=C(\bar\partial_MF). \end{align*} In particular, $\bar\partial_MF=0$ implies $\bar\partial_M(CF)=0$. For a general smooth $F$, choose an exhaustion $(K_a)_{a\geq1}$ of $\Omega$ by compact Runge subsets and smooth cutoffs \begin{align*} \rho_a:\Omega&\to[0,1] \end{align*} such that $\rho_a=1$ on a neighbourhood of $K_a$ and $\operatorname{supp}\rho_a\subset K_{a+1}$. Choose also an exhaustion $(L_a)_{a\geq1}$ of $M$ by compact subsets when $M$ is not a point. Set \begin{align*} V_a:C^\infty(M\times\Omega)&\to C^\infty(M\times\Omega)\\ F&\mapsto C(\rho_aF). \end{align*} On a neighbourhood of $M\times K_a$, \begin{align*} \frac{\partial V_a(F)}{\partial\bar z}=F. \end{align*} Hence $V_{a+1}(F)-V_a(F)$ is holomorphic in the variable $z$ on a neighbourhood of $M\times K_a$. By the Runge Approximation Theorem with Parameters, for each $a$ there is a smooth function \begin{align*} G_{a+1}:M\times\Omega&\to\mathbb C \end{align*} which is holomorphic in $z$ and approximates $V_{a+1}(F)-V_a(F)$ in the $C^a$ norm on $L_a\times K_a$ to within $2^{-a}$. Define \begin{align*} U_1&:=V_1(F),\\ U_{a+1}&:=V_{a+1}(F)-\sum_{b=2}^{a+1}G_b. \end{align*} The estimates imply that $(U_a)_{a\geq1}$ is Cauchy in $C^\infty$ on every compact subset of $M\times\Omega$. Let \begin{align*} S_\Omega F:M\times\Omega&\to\mathbb C \end{align*} be its $C^\infty_{\mathrm{loc}}$ limit. Since the correcting functions $G_b$ are holomorphic in $z$, they do not change the $\bar z$-derivative, and on each compact subset all sufficiently large $a$ satisfy \begin{align*} \frac{\partial U_a}{\partial\bar z}=F. \end{align*} Passing to the limit gives \begin{align*} \frac{\partial}{\partial\bar z}(S_\Omega F)=F. \end{align*} The parameter version of Runge approximation preserves passive holomorphic dependence; therefore, if $\bar\partial_MF=0$, the functions $G_b$ may be chosen with $\bar\partial_MG_b=0$, and the limiting solution satisfies $\bar\partial_M(S_\Omega F)=0$. [/proof] [/step]

[step:Prove the scalar anti-holomorphic lemma by induction on active variables]For $1\leq k\leq n$, write \begin{align*} D_k&:=\Omega_1\times\cdots\times\Omega_k. \end{align*} Let $M$ denote either a point or a product of the remaining planar domains, used only as a passive parameter space. For $r\geq0$, let $\mathcal E^{0,r}_k(M\times D_k)$ denote the smooth forms on $M\times D_k$ which involve only the anti-holomorphic differentials $d\bar z_1,\ldots,d\bar z_k$ and have anti-holomorphic degree $r$. Define \begin{align*} \bar\partial_k:\mathcal E^{0,r}_k(M\times D_k)&\to\mathcal E^{0,r+1}_k(M\times D_k)\\ \sum_{|J|=r}a_J\,d\bar z_J&\mapsto \sum_{\ell=1}^k\sum_{|J|=r}\frac{\partial a_J}{\partial\bar z_\ell}\,d\bar z_\ell\wedge d\bar z_J, \end{align*} where each coefficient \begin{align*} a_J:M\times D_k&\to\mathbb C \end{align*} is smooth. We prove by induction on $k$ that, for every $r\geq1$, every $\alpha\in\mathcal E^{0,r}_k(M\times D_k)$ satisfying $\bar\partial_k\alpha=0$ has a primitive $\beta\in\mathcal E^{0,r-1}_k(M\times D_k)$ satisfying $\bar\partial_k\beta=\alpha$. For $k=1$, the only nonzero case is $r=1$. Write \begin{align*} \alpha=f\,d\bar z_1, \end{align*} where \begin{align*} f:M\times\Omega_1&\to\mathbb C \end{align*} is smooth. By the one-variable solver, define \begin{align*} \beta:M\times\Omega_1&\to\mathbb C \end{align*} by $\beta:=S_{\Omega_1}f$. Then \begin{align*} \bar\partial_1\beta=\frac{\partial\beta}{\partial\bar z_1}\,d\bar z_1=f\,d\bar z_1=\alpha. \end{align*} If $r>1$, then $\mathcal E^{0,r}_1(M\times D_1)=0$, so $\alpha=0$ and $\beta=0$ works. Assume the statement holds for $k-1$, with arbitrary passive parameter spaces. Let $\alpha\in\mathcal E^{0,r}_k(M\times D_k)$ satisfy $\bar\partial_k\alpha=0$. Decompose $\alpha$ uniquely as \begin{align*} \alpha=\alpha_0+\alpha_1\wedge d\bar z_k, \end{align*} where \begin{align*} \alpha_0&\in\mathcal E^{0,r}_{k-1}(M\times D_k),\\ \alpha_1&\in\mathcal E^{0,r-1}_{k-1}(M\times D_k). \end{align*} Here $z_k$ is treated as an additional passive variable for the operator $\bar\partial_{k-1}$. The component of $\bar\partial_k\alpha$ containing no factor $d\bar z_k$ is $\bar\partial_{k-1}\alpha_0$, so $\bar\partial_{k-1}\alpha_0=0$. By the induction hypothesis, choose \begin{align*} \gamma\in\mathcal E^{0,r-1}_{k-1}(M\times D_k) \end{align*} such that \begin{align*} \bar\partial_{k-1}\gamma=\alpha_0. \end{align*} Define \begin{align*} R:=\alpha-\bar\partial_k\gamma. \end{align*} The form $R$ has no component without $d\bar z_k$, so there is a unique \begin{align*} \eta\in\mathcal E^{0,r-1}_{k-1}(M\times D_k) \end{align*} such that \begin{align*} R=\eta\wedge d\bar z_k. \end{align*} Since $\bar\partial_k\alpha=0$ and $\bar\partial_k^2\gamma=0$, we have $\bar\partial_kR=0$. But \begin{align*} \bar\partial_k(\eta\wedge d\bar z_k) =\bar\partial_{k-1}\eta\wedge d\bar z_k, \end{align*} because the term containing $d\bar z_k\wedge d\bar z_k$ vanishes by alternatingness of the exterior product. The map $\theta\mapsto\theta\wedge d\bar z_k$ is injective on forms using only $d\bar z_1,\ldots,d\bar z_{k-1}$, hence \begin{align*} \bar\partial_{k-1}\eta=0. \end{align*} If $r\geq2$, the induction hypothesis applied to $\eta$ gives \begin{align*} \theta\in\mathcal E^{0,r-2}_{k-1}(M\times D_k) \end{align*} with \begin{align*} \bar\partial_{k-1}\theta=\eta. \end{align*} Set \begin{align*} \beta:=\gamma+\theta\wedge d\bar z_k. \end{align*} Then \begin{align*} \bar\partial_k\beta &=\bar\partial_k\gamma+\bar\partial_k(\theta\wedge d\bar z_k)\\ &=\bar\partial_k\gamma+\bar\partial_{k-1}\theta\wedge d\bar z_k\\ &=\bar\partial_k\gamma+\eta\wedge d\bar z_k\\ &=\bar\partial_k\gamma+R\\ &=\alpha. \end{align*} If $r=1$, then $\eta$ is a smooth function \begin{align*} \eta:M\times D_k&\to\mathbb C \end{align*} satisfying $\bar\partial_{k-1}\eta=0$. Apply the one-variable solver in the variable $z_k$, with $M\times D_{k-1}$ as the passive parameter space, to obtain \begin{align*} h:M\times D_k&\to\mathbb C \end{align*} such that \begin{align*} \frac{\partial h}{\partial\bar z_k}=\eta \end{align*} and, because $\bar\partial_{k-1}\eta=0$ and the solver preserves passive $\bar\partial$-closedness, \begin{align*} \bar\partial_{k-1}h=0. \end{align*} Set \begin{align*} \beta:=\gamma+h. \end{align*} Then \begin{align*} \bar\partial_k\beta &=\bar\partial_k\gamma+\bar\partial_{k-1}h+\frac{\partial h}{\partial\bar z_k}\,d\bar z_k\\ &=\bar\partial_k\gamma+\eta\,d\bar z_k\\ &=\bar\partial_k\gamma+R\\ &=\alpha. \end{align*} This completes the induction and proves the scalar $(0,r)$ exactness for all $r\geq1$.[/step]

[guided]The induction separates the last anti-holomorphic differential $d\bar z_k$ from all earlier ones. Define $D_k:=\Omega_1\times\cdots\times\Omega_k$, and let $M$ be a passive parameter space. A form in $\mathcal E^{0,r}_k(M\times D_k)$ is allowed to have coefficients depending on all variables in $M\times D_k$, but its exterior factors are only among $d\bar z_1,\ldots,d\bar z_k$. The active Dolbeault operator is \begin{align*} \bar\partial_k:\mathcal E^{0,r}_k(M\times D_k)&\to\mathcal E^{0,r+1}_k(M\times D_k)\\ \sum_{|J|=r}a_J\,d\bar z_J&\mapsto \sum_{\ell=1}^k\sum_{|J|=r}\frac{\partial a_J}{\partial\bar z_\ell}\,d\bar z_\ell\wedge d\bar z_J. \end{align*} The base case $k=1$ is exactly the one-variable $\bar\partial$ problem. If \begin{align*} \alpha=f\,d\bar z_1, \end{align*} with $f:M\times\Omega_1\to\mathbb C$ smooth, the one-variable Cauchy-Green solver gives a smooth function $\beta:=S_{\Omega_1}f$ satisfying \begin{align*} \frac{\partial\beta}{\partial\bar z_1}=f. \end{align*} Therefore \begin{align*} \bar\partial_1\beta=f\,d\bar z_1=\alpha. \end{align*} There are no nonzero $(0,r)$-forms in one active variable for $r>1$, so the base case is complete. Now assume exactness has been proved for $k-1$ active variables, with arbitrary passive parameters. Let $\alpha\in\mathcal E^{0,r}_k(M\times D_k)$ satisfy $\bar\partial_k\alpha=0$. We split $\alpha$ according to whether a term contains $d\bar z_k$: \begin{align*} \alpha=\alpha_0+\alpha_1\wedge d\bar z_k, \end{align*} where $\alpha_0$ uses only $d\bar z_1,\ldots,d\bar z_{k-1}$ and has degree $r$, while $\alpha_1$ uses only $d\bar z_1,\ldots,d\bar z_{k-1}$ and has degree $r-1$. The condition $\bar\partial_k\alpha=0$ has a component with no $d\bar z_k$ factor. That component is precisely $\bar\partial_{k-1}\alpha_0$, so \begin{align*} \bar\partial_{k-1}\alpha_0=0. \end{align*} The induction hypothesis applies because $z_k$ is merely a passive parameter for $\bar\partial_{k-1}$. Hence there is \begin{align*} \gamma\in\mathcal E^{0,r-1}_{k-1}(M\times D_k) \end{align*} such that \begin{align*} \bar\partial_{k-1}\gamma=\alpha_0. \end{align*} Subtract the exact form $\bar\partial_k\gamma$ from $\alpha$: \begin{align*} R:=\alpha-\bar\partial_k\gamma. \end{align*} The purpose of this subtraction is to remove every term not involving $d\bar z_k$. Indeed, the no-$d\bar z_k$ part of $\bar\partial_k\gamma$ is $\bar\partial_{k-1}\gamma=\alpha_0$. Thus $R$ has the form \begin{align*} R=\eta\wedge d\bar z_k \end{align*} for a uniquely determined \begin{align*} \eta\in\mathcal E^{0,r-1}_{k-1}(M\times D_k). \end{align*} The residual $R$ is still closed: \begin{align*} \bar\partial_kR =\bar\partial_k\alpha-\bar\partial_k^2\gamma =0. \end{align*} Computing $\bar\partial_kR$ from $R=\eta\wedge d\bar z_k$ gives \begin{align*} \bar\partial_k(\eta\wedge d\bar z_k) =\bar\partial_{k-1}\eta\wedge d\bar z_k, \end{align*} because the part involving $\partial\eta/\partial\bar z_k$ contains $d\bar z_k\wedge d\bar z_k$ and vanishes by alternatingness. Since wedging by $d\bar z_k$ is injective on forms using only the earlier differentials, we get \begin{align*} \bar\partial_{k-1}\eta=0. \end{align*} If $r\geq2$, then $\eta$ has positive degree $r-1$ in the first $k-1$ variables. The induction hypothesis applies again, giving \begin{align*} \theta\in\mathcal E^{0,r-2}_{k-1}(M\times D_k) \end{align*} with \begin{align*} \bar\partial_{k-1}\theta=\eta. \end{align*} Then \begin{align*} \beta:=\gamma+\theta\wedge d\bar z_k \end{align*} satisfies \begin{align*} \bar\partial_k\beta &=\bar\partial_k\gamma+\bar\partial_k(\theta\wedge d\bar z_k)\\ &=\bar\partial_k\gamma+\bar\partial_{k-1}\theta\wedge d\bar z_k\\ &=\bar\partial_k\gamma+\eta\wedge d\bar z_k\\ &=\bar\partial_k\gamma+R\\ &=\alpha. \end{align*} If $r=1$, then $\eta$ is a function rather than a positive-degree form, so the induction hypothesis in the first $k-1$ variables cannot solve $\bar\partial_{k-1}\theta=\eta$. This is exactly where the one-variable solver in the last variable is needed. Since $\bar\partial_{k-1}\eta=0$, the function $\eta$ is $\bar\partial$-closed in the passive variables $z_1,\ldots,z_{k-1}$. Applying the one-variable solver in $z_k$ produces a smooth function \begin{align*} h:M\times D_k&\to\mathbb C \end{align*} such that \begin{align*} \frac{\partial h}{\partial\bar z_k}=\eta \end{align*} and the parameter-preserving property gives \begin{align*} \bar\partial_{k-1}h=0. \end{align*} With $\beta:=\gamma+h$, we compute \begin{align*} \bar\partial_k\beta &=\bar\partial_k\gamma+\bar\partial_{k-1}h+\frac{\partial h}{\partial\bar z_k}\,d\bar z_k\\ &=\bar\partial_k\gamma+\eta\,d\bar z_k\\ &=\bar\partial_k\gamma+R\\ &=\alpha. \end{align*} Thus every closed scalar $(0,r)$-form in $k$ active variables has a primitive, and the induction is complete.[/guided]

[step:Apply the scalar result to each holomorphic multi-index component] Let $\mathcal I_p$ denote the set of increasing multi-indices \begin{align*} I=(i_1,\ldots,i_p),\qquad 1\leq i_1<\cdots<i_p\leq n, \end{align*} and write \begin{align*} dz_I:=dz_{i_1}\wedge\cdots\wedge dz_{i_p}. \end{align*} Every form $\omega\in\mathcal E^{p,q}(D)$ has a unique decomposition \begin{align*} \omega=\sum_{I\in\mathcal I_p}dz_I\wedge\alpha_I, \end{align*} where each \begin{align*} \alpha_I\in\mathcal E^{0,q}(D) \end{align*} is a smooth anti-holomorphic $q$-form. Since $\bar\partial dz_I=0$ and $dz_I$ has degree $p$, \begin{align*} \bar\partial(dz_I\wedge\alpha_I)=(-1)^p dz_I\wedge\bar\partial\alpha_I. \end{align*} The forms $dz_I\wedge d\bar z_J$ are linearly independent over $C^\infty(D)$, and $\bar\partial\omega=0$, so \begin{align*} \bar\partial\alpha_I=0 \end{align*} for every $I\in\mathcal I_p$. By the scalar result with $k=n$ and no passive parameter space, for each $I\in\mathcal I_p$ there exists \begin{align*} \beta_I\in\mathcal E^{0,q-1}(D) \end{align*} such that \begin{align*} \bar\partial\beta_I=\alpha_I. \end{align*} Define \begin{align*} u:=(-1)^p\sum_{I\in\mathcal I_p}dz_I\wedge\beta_I. \end{align*} Then $u\in\mathcal E^{p,q-1}(D)$, and \begin{align*} \bar\partial u &=(-1)^p\sum_{I\in\mathcal I_p}\bar\partial(dz_I\wedge\beta_I)\\ &=(-1)^p\sum_{I\in\mathcal I_p}(-1)^p dz_I\wedge\bar\partial\beta_I\\ &=\sum_{I\in\mathcal I_p}dz_I\wedge\alpha_I\\ &=\omega. \end{align*} Thus every $\bar\partial$-closed $(p,q)$-form with $q\geq1$ is $\bar\partial$-exact. [/step]

[step:Identify the Dolbeault cohomology groups as zero] For $q\geq1$, the $q$-th cohomology of the Dolbeault complex $\mathcal E^{p,\bullet}(D)$ is \begin{align*} H^{p,q}_{\bar\partial}(D) =\frac{\ker\left(\bar\partial:\mathcal E^{p,q}(D)\to\mathcal E^{p,q+1}(D)\right)} {\operatorname{im}\left(\bar\partial:\mathcal E^{p,q-1}(D)\to\mathcal E^{p,q}(D)\right)}. \end{align*} The previous step proves that every element of the numerator belongs to the denominator. The reverse inclusion follows from $\bar\partial^2=0$. Hence the numerator and denominator are equal, and \begin{align*} H^{p,q}_{\bar\partial}(D)=0 \end{align*} for every $q\geq1$. This is exactly the asserted exactness of the Dolbeault complex $\mathcal E_D^{p,\bullet}$ in positive anti-holomorphic degree. [/step]