Integral formulas are your friends when you want to deal with estimates. This is because for many purposes, integrals behave better than derivatives. Let us recall how we can prove such a statement in single-variable complex analysis:

Suppose $f$ is entire and satisfies the bound $|f(z)|\leq A|z|^k$ (by modifying the proof below, you can consider bounds like $|f(z)|\leq \phi(|z|)$, and depending on the properties of the radial function $\phi$, deduce corresponding properties of $f$). We have a Taylor expansion $f(z)=\sum\limits_{n=0}^{\infty}c_nz^n$ on $\Bbb{C}$. Now, for each $n\geq 0$, and each $R>0$ we have that \begin{align} c_n&=\frac{1}{2\pi i}\int_{|z|=R}\frac{f(z)}{z^{n+1}}\,dz. \end{align} This immediately implies the bound \begin{align} |c_n|\leq \frac{1}{2\pi}\cdot\frac{AR^k}{R^{n+1}}\cdot 2\pi R=\frac{A}{R^{n-k}}. \end{align} As soon as $n-k>0$, we can take the limit $R\to\infty$ to find that $c_n=0$. Thus, we see that $c_{k+1}=c_{k+2}=\dots =0$, meaning that $f$ is a polynomial of degree at most $k$.

I hope you see how to adapt this to $n$-dimensions: write a Taylor expansion $f(z)=\sum_{\alpha}c_{\alpha}z^{\alpha}$, and write \begin{align} c_{\alpha}&=\frac{1}{(2\pi i)^n}\int_{|z_1|=\dots=|z_n|=R}\frac{f(z)}{z^{\alpha+1}}\,dz_1\dots\,dz_n, \end{align} where $z^{\alpha+1}$ is the shorthand multindex notation for $(z_1)^{\alpha_1+1}(z_2)^{\alpha_2+1}\cdots(z_n)^{\alpha_n+1}$. Now, bounding this will give you (letting $|\alpha|=\alpha_1+\dots+\alpha_n$) \begin{align} |c_{\alpha}|\leq \frac{1}{(2\pi)^n}\cdot\frac{AR^k}{R^{|\alpha|+n}}\cdot(2\pi R)^n=\frac{A}{R^{|\alpha|-k}}. \end{align} So, as soon as the indices are large enough: $|\alpha|-k>0$, we can take the limit $R\to\infty$ to conclude $c_{\alpha}=0$. Hence, $f$ is again a polynomial of degree at most $k$.

In fact, if you only assumed that $f$ was holomorphic on a region like $\{|z|>r_0\}$, then you could write a Laurent expansion centered at the origin, and the proof would show that if $f$ satisfied such a bound, then its Laurent expansion only contains terms up to power $k$ (by working in the limit $R\to 0^+$, you can prove that if you have bounds on $f$ near the origin, then you can place restrictions on the negative part of the Laurent coefficients).

So, purely from the integral representation formula of holomorphic functions, we are able to deduce how the Laurent-coefficients behave as soon as we know bounds on $f$ in the appropriate region (infinity or near the origin).

Going further, harmonic functions have the mean-value property, which is a form of an integral representation. Try to see if you can generalize these arguments for harmonic functions too (very similar statements hold).

Integral formulas are your friends when you want to deal with estimates. This is because for many purposes, integrals behave better than derivatives. Let us recall how we can prove such a statement in single-variable complex analysis:

Suppose $f$ is entire and satisfies the bound $|f(z)|\leq A|z|^k$ (by modifying the proof below, you can consider bounds like $|f(z)|\leq \phi(|z|)$, and depending on the properties of the radial function $\phi$, deduce corresponding properties of $f$). We have a Taylor expansion $f(z)=\sum\limits_{n=0}^{\infty}c_nz^n$ on $\Bbb{C}$. Now, for each $n\geq 0$, and each $R>0$ we have that \begin{align} c_n&=\frac{1}{2\pi i}\int_{|z|=R}\frac{f(z)}{z^{n+1}}\,dz. \end{align} This immediately implies the bound \begin{align} |c_n|\leq \frac{1}{2\pi}\cdot\frac{AR^k}{R^{n+1}}\cdot 2\pi R=\frac{A}{R^{n-k}}. \end{align} As soon as $n-k>0$, we can take the limit $R\to\infty$ to find that $c_n=0$. Thus, we see that $c_{k+1}=c_{k+2}=\dots =0$, meaning that $f$ is a polynomial of degree at most $k$.

I hope you see how to adapt this to $n$-dimensions: write a Taylor expansion $f(z)=\sum_{\alpha}c_{\alpha}z^{\alpha}$, and write \begin{align} c_{\alpha}&=\frac{1}{(2\pi i)^n}\int_{|z_1|=\dots=|z_n|=R}\frac{f(z)}{z^{\alpha+1}}\,dz_1\dots\,dz_n, \end{align} where $z^{\alpha+1}$ is the shorthand multindex notation for $(z_1)^{\alpha_1+1}(z_2)^{\alpha_2+1}\cdots(z_n)^{\alpha_n+1}$. Now, bounding this will give you (letting $|\alpha|=\alpha_1+\dots+\alpha_n$) \begin{align} |c_{\alpha}|\leq \frac{1}{(2\pi)^n}\cdot\frac{AR^k}{R^{|\alpha|+n}}\cdot(2\pi R)^n=\frac{A}{R^{|\alpha|-k}}. \end{align} So, as soon as the indices are large enough: $|\alpha|-k>0$, we can take the limit $R\to\infty$ to conclude $c_{\alpha}=0$. Hence, $f$ is again a polynomial of degree at most $k$.

In fact, if you only assumed that $f$ was holomorphic on a region like $\{|z|>r_0\}$, then you could write a Laurent expansion centered at the origin, and the proof would show that if $f$ satisfied such a bound, then its Laurent expansion only contains terms up to power $k$ (by working in the limit $R\to 0^+$, you can prove that if you have bounds on $f$ near the origin, then you can place restrictions on the negative part of the Laurent coefficients).

So, purely from the integral representation formula of holomorphic functions, we are able to deduce how the Laurent-coefficients behave as soon as we know bounds on $f$ in the appropriate region (infinity or near the origin).

Going further, harmonic functions have the mean-value property, which is a form of an integral representation. Try to see if you can generalize these arguments for harmonic functions too (very similar statements hold). Also, the fact that $f$ is $\Bbb{C}$-valued is irrelevant here; you can replace it with any complex Banach space $V$ (assume finite-dimensional for comfort if you like, and now the coefficients $c_{\alpha}$ will be elements of $V$).