IRL ABooklet

Solution 1. When $n = 1$ this is obviously possible. We will show that this is possible for all integers $n \ge 2$ by proving a more general version using induction. The number 2021 will be replaced by an integer $k \ge 0$ on which we will carry out the induction. For this approach to work we will consider powers of numbers in a range that does not necessarily start with 1 and we will consider sums of powers lower than $k$.
To formulate this more general statement, we first introduce some notation. When $a < b$ are two integers, we denote by $[a, b]$ the set of integers $t$ that satisfy $a \le t \le b$. When $n \ge 2$ and $m \ge 0$ we let $d(m)$ be the digit sum of the base $n$ representation of $m$. Given two integers, $k \ge 0$ and $n \ge 2$, let $N = n^{k+1}$. For $i = 0, 1, \dots, n-1$ we define subsets of $[0, N-1]$ as follows
$$S_i = \{m \mid d(m) \equiv i \pmod n \text{ and } 0 \le m < N\}.$$
We also let $S_i(s) = \{m + s \mid m \in S_i\} \subset [s, N + s - 1]$ be the set obtained by shifting $S_i$ by the integer $s$. Finally, for $0 \le r \le k$ we define
$$\sigma_i^r(s) = \sum_{m \in S_i(s)} m^r.$$
Here, even when $m = 0$, we agree that $m^0 = 1$. We now fix an integer $n \ge 2$. Theorem. For all integers $n, k, r, s$ satisfying $n \ge 2$ and $0 \le r \le k$ we have
$$\sigma_0^r(s) = \sigma_1^r(s) = \dots = \sigma_{n-1}^r(s).$$
Remark. Setting $s = 1$ and $r = k = 2021$ solves the original problem.
Proof. We do induction on $k \ge 0$. We keep $n$ fixed throughout the proof. The base case, $k = 0$, is trivially true, because the only possible value for $r$ is zero, $m^0 = 1$ and so $\sigma_i^0(s) = n^k$ is the number of elements in $S_i(s)$. Let us now fix $k > 0$ and suppose that the statement of the theorem holds true when $k$ is replaced by $k-1$. To avoid confusion, we denote the sets $S_i(s)$ that correspond to $k-1$ by $\tilde{S}_i(s)$:
$$\tilde{S}_i = \{m \mid d(m) \equiv i \pmod n \text{ and } 0 \le m < \tilde{N}\}$$
where $\tilde{N} = n^k$. Accordingly, we define $\tilde{S}_i(s)$ and $\tilde{\sigma}_i^r(s)$. The base $n$ representation of the numbers in $[0, N-1]$ is obtained from those of the numbers in $[0, \tilde{N}-1]$ by prepending digits $0, 1, 2, \dots, n-1$. Because $d(m+cn^k) = d(m)+c$ for $m \in [0, \tilde{N}-1]$ and $c = 0, 1, \dots, n-1$, we see that
$$S_i = \tilde{S}_i(0) \cup \tilde{S}_{i-1}(n^k) \cup \tilde{S}_{i-2}(2n^k) \cup \dots \cup \tilde{S}_{i-(n-1)}((n-1)n^k). \quad (11)$$
Here we consider the subscripts modulo $n$. As we assume for $0 \le r < k$ that $\tilde{\sigma}_j^r(s) = \sum_{m \in \tilde{S}_j(s)} m^r$ does not depend on $j$, we immediately get that $\sigma_i^r(s)$ does not depend on $i$, as long as $0 \le r < k$.
We are left to prove that $\sigma_i^k(s)$ does not depend on $i$. Using (11) we obtain
$$\begin{aligned} \sum_{m \in S_i(s)} m^k &= \sum_{m \in \tilde{S}_i(s)} m^k + \sum_{m \in \tilde{S}_{i-1}(s+n^k)} m^k + \dots + \sum_{m \in \tilde{S}_{i-(n-1)}(s+(n-1)n^k)} m^k \\ &= \sum_{m \in \tilde{S}_i(s)} m^k + \sum_{m \in \tilde{S}_{i-1}(s)} (m+n^k)^k + \dots + \sum_{m \in \tilde{S}_{i-(n-1)}(s)} (m+(n-1)n^k)^k. \end{aligned}$$
We now expand the expressions $(m + jn^k)^k$ using the binomial theorem. What we obtain is of the form $m^k + c_1m^{k-1} + \dots + c_{k-1}m + c_k$ where each coefficient $c_g$ is a product of a binomial coefficient and a power of $jn^k$. The exact values of these coefficients are not relevant here, but it is important that they are the same for all $m$. Therefore,
$$\begin{aligned} \sum_{m \in \tilde{S}_{i-j}(s)} (m + jn^k)^k &= \sum_{m \in \tilde{S}_{i-j}(s)} m^k + c_1 \sum_{m \in \tilde{S}_{i-j}(s)} m^{k-1} + \dots + c_k \sum_{m \in \tilde{S}_{i-j}(s)} m^0 \\ &= \sum_{m \in \tilde{S}_{i-j}(s)} m^k + c_1 \tilde{\sigma}_{i-j}^{k-1}(s) + c_2 \tilde{\sigma}_{i-j}^{k-2}(s) + \dots + c_k \tilde{\sigma}_{i-j}^0(s). \end{aligned}$$
By inductive assumption, no term that comes with a coefficient $c_g$ depends on $i$. Therefore,
$$\begin{aligned} \sum_{m \in S_i(s)} m^k &= \sum_{j=0}^{n-1} \sum_{m \in \tilde{S}_{i-j}(s)} m^k + \text{terms not depending on } i \\ &= \sum_{m \in [s, s+\tilde{N}-1]} m^k + \text{terms not depending on } i \end{aligned}$$
and this does not depend on $i$ and the proof is complete. The last equality comes from
$$\tilde{S}_0(s) \cup \tilde{S}_1(s) \cup \dots \cup \tilde{S}_{n-1}(s) = [s, s + \tilde{N} - 1]$$
and $\tilde{S}_a(s) = \tilde{S}_{a+n}(s)$ because we considered subscripts modulo $n$. $\square$

Solution 2. We will again replace 2021 by an integer $k \ge 0$. This solution takes a more structured approach using polynomials. Like in Solution 1, we denote by $[a, b]$ the set of integers $t$ that satisfy $a \le t \le b$, where $a < b$ are two given integers. When $n \ge 2$ and $m \ge 0$ we let $d(m)$ be the digit sum of the base $n$ representation of $m$. Given any two integers, $k \ge 0$ and $n \ge 2$, we let $N = n^{k+1}$ and for $i = 0, 1, \dots, n-1$ we define subsets of $[1, N]$ as follows
$$S_i = \{m \mid d(m-1) \equiv i \pmod{n} \text{ and } 1 \le m \le N\}.$$
We are going to prove that for any two integers, $k \ge 0$ and $n \ge 2$,
$$\sum_{a \in S_i} a^k = \sum_{b \in S_j} b^k \quad \text{for all } i, j.$$
Define $f(x) = 1+x+x^2+\dots+x^{n-1}$ and let $f_j(x,y) = f(x^{nj}y)$ for $0 \le j \le k$. We consider the polynomial $F(x,y) = x \cdot \prod_{j=0}^k f_j(x,y)$ and we think of $F$ as a polynomial in $x$ that has coefficients which are polynomials in $y$. Note that setting $x=1$ gives $f_j(1,y) = f(y)$ for all $j$. This will be of importance towards the end of the argument.
When we multiply out completely the product that defines $F$, we obtain an expression of the form
$$F(x,y) = \sum_{m=1}^{N} c_m(y)x^m$$
with coefficients $c_m(y)$ being polynomials in $y$. The exponents of $x$ that appear in $f_j(x,y)$ are $0, n^j, 2n^j, 3n^j, \dots, (n-1)n^j$ and the coefficient in front of $x^{cn^j}$ is $y^c$. Therefore, $c_m(y) = y^{d(m-1)}$ where $d(m-1)$ is the digit sum of the base $n$ representation of $m-1$.
Differentiating $p(x) = x^m$ with respect to $x$ gives $p'(x) = m x^{m-1}$. When we differentiate $i$ times, we get $p^{(i)}(x) = (m-i+1)(m-i+2)\cdots(m-1)m x^{m-i}$. This is conveniently written using the falling factorial power notation (see for example [1, p. 47])
$$m^i = m(m-1)(m-2)\cdots(m-i+1)$$
which allows us to write $p^{(i)}(x) = m^i x^{m-i}$ for the $i$-th derivative of $p(x) = x^m$. This holds true even if $i > m$, because then both sides are equal to zero. Substituting $x=1$ we obtain
$$p^{(i)}(1) = m^i. \qquad (12)$$
Usual powers $m^k$ and the falling powers $m^0, m^1, \dots, m^{k-1}, m^k$ are related in a beautiful way involving the Stirling numbers of the first kind, [1, p. 248]:
$$m^k = \sum_{i=0}^{k} \binom{k}{i} m^i. \qquad (13)$$
The Stirling numbers $\binom{k}{i}$ are defined combinatorially as the number of ways to partition a set that contains $k$ elements into $i$ non-empty subsets. For $k \ge 0$ and $i \ge 0$, these numbers are determined by the recurrence relation
$$\binom{k}{i} = i \binom{k-1}{i} + \binom{k-1}{i-1} \quad \text{for } k, i > 0$$
together with the initial values $\binom{0}{0} = 1$ and $\binom{r}{0} = \binom{0}{r} = 0$ if $r > 0$. In our situation, it is sufficient to know that such a formula exists, the exact values of the coefficients do not matter.
Combining equations (12) and (13) we get $m^k = \sum_{i=0}^{k} \binom{k}{i} p^{(i)}(1)$ if $p(x) = x^m$. Recalling that we wrote $F(x, y) = \sum_{m=1}^{N} c_m(y)x^m$, we now see that
$$\sum_{m=1}^{N} c_m(y)m^k = \sum_{i=0}^{k} \binom{k}{i} F^{(i)}(y, 1) \quad (14)$$
where $F^{(i)}(y, 1)$ is the polynomial in $y$ obtained by setting $x = 1$ in the $i$-th derivative of $F(x, y)$ with respect to $x$.
Since we defined $F(x, y) = x \cdot \prod_{j=0}^{k} f_j(x, y)$, the product rule for derivatives shows that, for $i \le k$, $F^{(i)}(y, 1)$ can be written as a sum of products, each having at least one factor of the form $f_j(1, y)$. Since $f_j(1, y) = f(y)$, as we noted earlier, the polynomial (14) is a multiple of $f(y)$. In other words
$$\sum_{m=1}^{N} c_m(y)m^k \equiv 0 \pmod{f(y)}.$$
Because $(y-1)f(y) = y^n-1$, we can replace $y^n$ by 1, or more generally, $y^d$ by $y^r$ whenever $d \equiv r \pmod{n}$, when we work (mod $f(y)$). When $c_m(y) = y^d$ and $0 \le r < n$ is the remainder of $d$ on division by $n$, we let $\bar{c}_m(y) = y^r$. We then have $\sum_{m=1}^{N} \bar{c}_m(y)m^k \equiv 0 \pmod{f(y)}$. This means there exists $C \in \mathbb{Z}$ such that
$$\sum_{m=1}^{N} \bar{c}_m(y)m^k = C f(y) = C(1 + y + y^2 + \dots + y^{n-1}).$$
As we have seen above, $c_m(y) = y^{d(m-1)}$. Thus, $\bar{c}_m(y) = y^r$ exactly when $d(m-1)$ leaves remainder $r \in [0, n-1]$ on division by $n$ and so $C = \sum_{m \in S_i} m^k$ for all $i = 0, 1, \dots, n-1$.