Here is the necessary and sufficient condition, in terms of $\det A$ and $\text{tr}(A)$, in order that a $2\times2$ matrix $A$ be the $k$-th power of some matrix with integer coefficients.

**Edit.Edit, 11.11.2020** The proof.The proof, is essentially routine; yet I post everything for convenience, before I forget all details.

It is convenient to introduce the polynomials

$$P_k(x,y):=\sum_{j\ge0}(-1)^j{ k-j\choose j }x^{k-2j}y^j;$$
their relevance in this context being that the polynomial $z^2-xz+y$ divides the polynomial $$z^{k+2}-xP_{k+1}(x,y)z+yP_{k}(x,y)$$ (even as elements of $\mathbb{Z}[x,y,z]$; see below for other properties we need).

**Characterization of the $k$-th powers in $ M_2(\mathbb Z)$.** Let $k\ge0$. *A matrix $A\in M_2 (\mathbb Z)$ is a $k$-th power of an element of $ M_2 (\mathbb Z)$ if and only if there are $t,d$ in $\mathbb Z$ such that*

**1.** $\det(A)=d^k$

**2.** $\text{tr}(A)=P_k(t,d)-P_{k-2}(t,d)d$

**3.** $P_{k-1}(t,d)$ divides $A+ P_{k-2}(t,d)d I $

*Precisely, if $B\in M_2 (\mathbb Z)$ verifies $B^k=A$ then $d:=\det(B)$ and $t:=\text{tr}(B)$ satisfy* **(1,2,3)**.

*Conversely, if $(d,t)$ satisfy ***(1,2,3)**, there exists a $B\in M_2 (\mathbb Z)$ with characteristic polynomial $z^2-tz+d$.
Precisely, if $P_{k-1}(t,d)\ne0$, it is unique, namely
$$B:=\frac1{P_{k-1}(t,d)}\Big(A+ P_{k-2}(t,d)d I \Big).$$
If $P_{k-1}(t,d)=0$, then $A=mI$, is an integer multiple of the identity, and all the infinitely many $B\in M_2 (\mathbb Z)$ with
with characteristic polynomial $z^2-tz+d$ satisfy $B^k=A$.

**Proof**. Assume $A=B^k$ and $B\in M_2 (\mathbb Z)$ and set $t:=\text{tr}(B)$ and $d:=\det(B)$. Then **(1)** is $\det(A)=\det(B^k)=\det(B)^k=p^k$. As seen above, the characteristic polynomial of $B$, $p_B(z):=z^2-tz+d $ divides the polynomial $z^k-P_{k-1}(t,d)z+P_{k-2}(t,d)d$, and since by Cayley-Hamilton $B^2-tB+d=0$, we also have
$$B^k-P_{k-1}(t,d)B+P_{k-2}(t,d)d I=0,$$
so, taking the trace, we have $\text{tr}(A)= P_{k-1}(t,d)t -2P_{k-2}(t,d)d$, which is **(2)**, while $P_{k-1}(t,d)B= A+P_{k-2}(t,d)d I$, is **(3)**.

Conversely, assume the above conditions **(1,2,3)** hold for integers $t,d$. Consider first the case $P_{k-1}(t,d)\ne0$. So one can define
$$B:=\frac1{P_{k-1}(t,d)}\Big(A+ P_{k-2}(t,d)d I \Big),$$
an element of $M_2 (\mathbb Z)$ thanks to **(3)**.
The trace and determinant of $B$ are then by **(1,2)**, hidding the variables $(t,d)$ in the $P_j$
$$\text{tr}(B)=\frac{\text{tr}(A)+ 2P_{k-2}d}{P_{k-1} } = \frac{P_k +P_{k-2} d }{P_{k-1} } =\frac{ t P_{k-1} }{P_{k-1} } =t $$
$$\det(B)= \frac{\det\Big(A+ P_{k-2} d I \Big)}{P_{k-1} ^2}
= \frac{P_{k-2}^2d^2 + \text{tr}(A)P_{k-2} d +\det(A) }{P_{k-1}^2}=$$ $$=\frac{P_{k-2}^2d^2 + \big(P_k-P_{k-2}d\big)P_{k-2}d +d^k }{P_{k-1}^2}= $$
$$=\frac{ P_kP_{k-2}d +d^k }{P_{k-1}^2}= d,$$
because $P_{k-1}^2 -P_kP_{k-2} =d^{k-1}$. Thus the characteristic polynomial of $B$ is $z^2-tz+d$, which implies $B^k=P_{k-1}B-P_{k-2}d\, I=A$.
Finally, consider the case $P_{k-1}(t,d)=0$. By **(3)** $A$ is then a multiple of the identity, $A=m I$, for $m:=-P_{k-2}(t,d)d$. If $m=0$, any nilpotent $B$ has the wanted properties. If $m\ne 0$, let $\lambda$ and $\mu$ be the roots of $z^2-tz+d$, so $t=\lambda+\mu$ and $d=\lambda\mu$. Then we have $\lambda\neq\mu$, otherwise $0=P_{k-1}(\lambda+\mu,\lambda\mu)=k\lambda^{k-1}$ and $\lambda=\mu=0=t=d$ and $A=0$. Also (see below)
$$0= P_{k-1}(\lambda+\mu,\lambda\mu)= \frac{\lambda^k-\mu^k}{\lambda-\mu} $$
whence $\lambda^k=\mu^k$, and
$$m =-\lambda\mu P_{k-2}(\lambda+\mu,\lambda\mu)=-\frac{\lambda^k\mu-\mu^k\lambda}{\lambda-\mu} = \lambda^k=\mu^k.$$
Let $B$ one of the infinitely many matrices in $M_2 (\mathbb Z)$ with characteristic polynomial $z^2-tz+d$. Since $\lambda\ne \mu$, $B$ is diagonalizable, $B=Q^{-1}\text{diag}(\lambda,\mu)Q$, so $$B^k=Q^{-1}\text{diag}(\lambda^k,\mu^k)Q=Q^{-1}mIQ=mI=A,$$
ending the proof.

**More details.** The sequence of polynomials $P_k(x,y)\in\mathbb{Z}[x,y]$ is defined by the two-term recurrence
$$\cases{P_{k+2}=xP_{k+1}-yP_k\\
P_0=1 \\ P_{-1}=0.}$$
One easily verifies by induction the expansion
$$P_k(x,y):=\sum_{j\ge0}(-1)^j{ k-j\choose j }x^{k-2j}y^j;$$
in fact $P_k$ may also be presented in terms of the Chebyshev polynomials of the first kind as $P_k(x,y^2)=y^kxT_k\big(\frac{x}{2y}\big)\in\mathbb{Z}[x,y^2]$.
They verify
$$P_k(u+v,uv)=\frac{u^{k+1}-v^{k+1}}{u-v}=\sum_{j=0}^{k} u^jv^{k-j},$$
and, related to that, for all $k\ge0$ one has:
$$z^{k+2}-P_{k+1}(x,y)z+yP_{k}(x,y)=\big(z^2-xz+y\big) \sum_{j=0}^kP_{k-j}(x,y)z^j, $$ both easily verified by induction. Finally, since they solve a two-term linear recursion, the Hankel determinant of order $2$ must be a $1$-term linear recurrence, and one finds
$$P_{k}(x,y)^2-P_{k+1}(x,y)P_{k-1}(x,y)=y^{k}.$$