[proofplan]
We prove the $n$-variable [Cauchy integral formula](/theorems/345) by iterating the classical one-variable [Cauchy integral formula](/theorems/345). Starting from $f(z)$, we apply the one-variable formula in $z_1$ with parameters $z_2, \dots, z_n$ fixed, obtaining a single contour integral whose integrand is holomorphic in $z_2, \dots, z_n$. We then apply the one-variable formula again in $z_2$ to this integrand, and continue inductively until all $n$ variables have been integrated out. The validity of each step rests on the fact that the integrand at stage $k$ is holomorphic in $z_{k+1}, \dots, z_n$, which follows from differentiation under the integral sign justified by [uniform continuity](/page/Uniform%20Continuity) on the compact torus.
[/proofplan]
[step:Apply the one-variable Cauchy formula in $z_1$ with $(z_2, \dots, z_n)$ held fixed]
Fix $z = (z_1, \dots, z_n) \in \mathbb{D}^n(a, r)$. Since $|z_j - a_j| < r_j$ for each $j$, the point $z_1$ lies in the open disc $\mathbb{D}(a_1, r_1)$. For each fixed $(z_2, \dots, z_n)$ with $|z_j - a_j| < r_j$ for $j = 2, \dots, n$, the function
\begin{align*}
g_{z_2, \dots, z_n}: \mathbb{D}(a_1, r_1) &\to \mathbb{C} \\
\zeta_1 &\mapsto f(\zeta_1, z_2, \dots, z_n)
\end{align*}
is holomorphic in $\zeta_1$, because $f$ is holomorphic on $\Omega$ and the polydisc $\overline{\mathbb{D}^n(a,r)} \subset \Omega$. By the one-variable [Cauchy Integral Formula](/theorems/???), applied to $g_{z_2, \dots, z_n}$ on the disc $\mathbb{D}(a_1, r_1)$ with $|z_1 - a_1| < r_1$:
\begin{align*}
f(z_1, z_2, \dots, z_n) &= \frac{1}{2\pi i} \oint_{|\zeta_1 - a_1| = r_1} \frac{f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}\, d\zeta_1.
\end{align*}
[guided]
The idea is to treat $f$ as a function of $z_1$ alone, freezing the remaining variables. Since $f \in \mathcal{O}(\Omega)$ and $\overline{\mathbb{D}^n(a,r)} \subset \Omega$, for any fixed $(z_2, \dots, z_n)$ in the smaller polydisc $\mathbb{D}^{n-1}(a', r')$ (where $a' = (a_2, \dots, a_n)$ and $r' = (r_2, \dots, r_n)$), the map $\zeta_1 \mapsto f(\zeta_1, z_2, \dots, z_n)$ is holomorphic on an [open set](/page/Open%20Set) containing $\overline{\mathbb{D}(a_1, r_1)}$. In particular, it is holomorphic on a neighbourhood of the closed disc, so the classical one-variable [Cauchy integral formula](/theorems/345) applies directly.
Define
\begin{align*}
g_{z_2, \dots, z_n}: \mathbb{D}(a_1, r_1) &\to \mathbb{C} \\
\zeta_1 &\mapsto f(\zeta_1, z_2, \dots, z_n).
\end{align*}
The one-variable [Cauchy Integral Formula](/theorems/???) states: if $h$ is holomorphic on a neighbourhood of $\overline{\mathbb{D}(c, \rho)}$ and $w \in \mathbb{D}(c, \rho)$, then $h(w) = \frac{1}{2\pi i}\oint_{|\zeta - c| = \rho} \frac{h(\zeta)}{\zeta - w}\, d\zeta$. Applying this to $h = g_{z_2, \dots, z_n}$, $c = a_1$, $\rho = r_1$, $w = z_1$:
\begin{align*}
f(z_1, z_2, \dots, z_n) &= \frac{1}{2\pi i} \oint_{|\zeta_1 - a_1| = r_1} \frac{f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}\, d\zeta_1.
\end{align*}
This is the first step of the iteration. The integrand still depends on $z_2, \dots, z_n$, and we will apply the Cauchy formula to these variables next.
[/guided]
[/step]
[step:Verify that the integrand is holomorphic in $(z_2, \dots, z_n)$ for each fixed $\zeta_1$ on the circle]
For each fixed $\zeta_1$ with $|\zeta_1 - a_1| = r_1$, the function
\begin{align*}
h_{\zeta_1}: \mathbb{D}^{n-1}(a', r') &\to \mathbb{C} \\
(z_2, \dots, z_n) &\mapsto \frac{f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}
\end{align*}
is holomorphic in $(z_2, \dots, z_n)$, since $f(\zeta_1, \cdot)$ is holomorphic on $\mathbb{D}^{n-1}(a', r')$ (because $(\zeta_1, z_2, \dots, z_n) \in \overline{\mathbb{D}^n(a,r)} \subset \Omega$ and $f \in \mathcal{O}(\Omega)$), and the factor $(\zeta_1 - z_1)^{-1}$ does not depend on $z_2, \dots, z_n$.
The integrand $(\zeta_1, z_2, \dots, z_n) \mapsto \frac{f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}$ is continuous on the compact set $\{|\zeta_1 - a_1| = r_1\} \times \overline{\mathbb{D}^{n-1}(a', s')}$ for any $s' < r'$ componentwise, hence uniformly bounded there. Differentiation under the integral sign with respect to $z_j$ (for $j \geq 2$) is therefore justified by [uniform convergence](/page/Uniform%20Convergence) of the difference quotients, confirming that the integral representation preserves holomorphicity in the remaining variables.
[guided]
Why does this step matter? To iterate the Cauchy formula, we need the integral in $\zeta_1$ to produce a function that is holomorphic in the remaining variables. If the integral somehow destroyed holomorphicity, we could not proceed.
For each fixed $\zeta_1$ on the circle $|\zeta_1 - a_1| = r_1$, the point $(\zeta_1, z_2, \dots, z_n)$ lies in $\overline{\mathbb{D}^n(a,r)} \subset \Omega$ whenever $(z_2, \dots, z_n) \in \mathbb{D}^{n-1}(a', r')$. Since $f \in \mathcal{O}(\Omega)$, the map $(z_2, \dots, z_n) \mapsto f(\zeta_1, z_2, \dots, z_n)$ is holomorphic on $\mathbb{D}^{n-1}(a', r')$. The denominator $\zeta_1 - z_1$ is independent of $z_2, \dots, z_n$ and nonzero (since $|\zeta_1 - a_1| = r_1$ while $|z_1 - a_1| < r_1$), so
\begin{align*}
h_{\zeta_1}: \mathbb{D}^{n-1}(a', r') &\to \mathbb{C} \\
(z_2, \dots, z_n) &\mapsto \frac{f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}
\end{align*}
is holomorphic in $(z_2, \dots, z_n)$.
To pass holomorphicity through the integral, we use differentiation under the integral sign. Fix any compact polydisc $\overline{\mathbb{D}^{n-1}(a', s')} \subset \mathbb{D}^{n-1}(a', r')$ with $s_j < r_j$ for each $j \geq 2$. On the product $\{|\zeta_1 - a_1| = r_1\} \times \overline{\mathbb{D}^{n-1}(a', s')}$, the function $f$ is continuous (as a restriction of a [holomorphic function](/page/Holomorphic%20Function) to a compact subset of its domain), hence bounded, say by some $M > 0$. The denominator satisfies $|\zeta_1 - z_1| \geq r_1 - |z_1 - a_1| > 0$ uniformly for $z$ in the compact polydisc. Therefore the integrand is uniformly bounded, and differentiation under the integral sign with respect to the Wirtinger operator $\partial_{\bar{z}_j}$ for $j \geq 2$ yields
\begin{align*}
\partial_{\bar{z}_j} \left[\frac{1}{2\pi i} \oint_{|\zeta_1 - a_1| = r_1} \frac{f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}\, d\zeta_1\right] = \frac{1}{2\pi i} \oint_{|\zeta_1 - a_1| = r_1} \frac{\partial_{\bar{z}_j} f(\zeta_1, z_2, \dots, z_n)}{\zeta_1 - z_1}\, d\zeta_1 = 0,
\end{align*}
since $f \in \mathcal{O}(\Omega)$ implies $\partial_{\bar{z}_j} f = 0$. This confirms that the integral is holomorphic in $(z_2, \dots, z_n)$.
[/guided]
[/step]
[step:Iterate the one-variable Cauchy formula in $z_2, \dots, z_n$ to obtain the $n$-fold integral]
We apply the one-variable Cauchy formula to the integrand in $z_2$. For each fixed $\zeta_1$ on $|\zeta_1 - a_1| = r_1$, the function $z_2 \mapsto \frac{f(\zeta_1, z_2, z_3, \dots, z_n)}{\zeta_1 - z_1}$ is holomorphic in $z_2$ on $\mathbb{D}(a_2, r_2)$ (as established above). By the one-variable Cauchy formula in $z_2$:
\begin{align*}
\frac{f(\zeta_1, z_2, z_3, \dots, z_n)}{\zeta_1 - z_1} &= \frac{1}{2\pi i} \oint_{|\zeta_2 - a_2| = r_2} \frac{f(\zeta_1, \zeta_2, z_3, \dots, z_n)}{(\zeta_1 - z_1)(\zeta_2 - z_2)}\, d\zeta_2.
\end{align*}
Substituting into the integral from the first step:
\begin{align*}
f(z) &= \frac{1}{(2\pi i)^2} \oint_{|\zeta_1 - a_1| = r_1} \oint_{|\zeta_2 - a_2| = r_2} \frac{f(\zeta_1, \zeta_2, z_3, \dots, z_n)}{(\zeta_1 - z_1)(\zeta_2 - z_2)}\, d\zeta_2\, d\zeta_1.
\end{align*}
The interchange of the $\zeta_1$-integral with the $\zeta_2$-integral is justified by the continuity of the integrand on the compact torus $\{|\zeta_1 - a_1| = r_1\} \times \{|\zeta_2 - a_2| = r_2\}$ and [Fubini's theorem](/theorems/2961) for continuous functions on compact sets.
Repeating this procedure for $z_3, z_4, \dots, z_n$ in succession, at the $k$-th stage we apply the one-variable Cauchy formula in $z_k$ to the integrand (which is holomorphic in $z_k$ by the same argument as above: $f$ is holomorphic in all its arguments, and the Cauchy kernel factors $(\zeta_j - z_j)^{-1}$ for $j < k$ do not involve $z_k$). After $n$ iterations:
\begin{align*}
f(z) &= \frac{1}{(2\pi i)^n} \oint_{|\zeta_1 - a_1| = r_1} \cdots \oint_{|\zeta_n - a_n| = r_n} \frac{f(\zeta_1, \dots, \zeta_n)}{(\zeta_1 - z_1)\cdots(\zeta_n - z_n)}\, d\zeta_n \cdots d\zeta_1.
\end{align*}
Since the integrand is continuous on the compact torus $\mathbb{T}^n(a, r) = \{|\zeta_j - a_j| = r_j : j = 1, \dots, n\}$, the iterated integral equals the integral over $\mathbb{T}^n(a, r)$, and the order of integration is immaterial. This is the desired formula.
[/step]