Series expansion of the determinant for a matrix near the identity
Solution 1:
Note that the determinant of a matrix is just a polynomial in the components of the matrix. This means that the series in question is finite. This also means that the functions $f_1,f_2,\dots,f_N$ are polynomials.
We start by computing the term which is first order in $\epsilon$.
Let $A_1, A_2, \dots , A_N$ be the column vectors of the matrix $A$. Let $e_1,e_2,\dots, e_N$ be the standard basis; note that these basis vectors form the columns of the identity matrix $I$. Then we recall that the determinant is an alternating multi-linear map on the column space.
$$ \mathrm{det}( I + \epsilon A ) = \mathrm{det}( e_1 + \epsilon A_1, e_2 + \epsilon A_2, \dots , e_N + \epsilon A_N ) $$
$$ = \mathrm{det}( e_1, e_2, \dots , e_N ) + \epsilon \lbrace \mathrm{det}(A_1,e_2,\dots, e_N) + \mathrm{det}(e_1,A_2,\dots, e_N) + \cdots + \mathrm{det}(e_1,e_2,\dots, A_N)\rbrace + O(\epsilon^2)$$
The first term is just the determinant of the identity matrix which is $1$. The term proportional to $\epsilon$ is a sum of expressions like $\mathrm{det}(e_1,e_2,\dots, A_j, \dots, e_N)$ where the $j$'th column of the identity matrix is replaced with the $j$'th column of $A$. Expanding the determinant along the $j$'th row we see that $\mathrm{det}(e_1,e_2,\dots, A_j, \dots, e_N) = A_{jj}$.
$$ \mathrm{det}( I + \epsilon A ) = 1 + \epsilon \sum_{j=1}^N A_{jj}+ O(\epsilon^2) = 1 + \epsilon \mathrm{Tr}(A) + O(\epsilon^2)$$
$$\boxed{ f_1(A) = \mathrm{Tr}(A) } \qquad \textbf{(1)} $$
We have the first term in our series in a computationally simple form. Our goal is to obtain higher order terms in the series in a similar form. To do this we will have to abandon the current method of attack and consider the determinant of the exponential map applied to a matrix.
If $A$ is diagonalizable then we can define $\exp(A)$ in terms of its action on the eigenvectors of $A$. If $a_1,a_2,\dots , a_N$ are eigenvectors of $A$ with eigenvalues $\lambda_1, \lambda_2, \dots, \lambda_N$ then $\exp(A)$ is the matrix which satisfies,
$$ \exp(A) a_j = e^{\lambda_j} a_j.$$
It is not hard to show that $\det( \exp(A) ) = \exp(\mathrm{Tr}(A))$. Since $A$ is linear operator on a finite dimensional vector space it has a finite norm. This means that we can safely evaluate the exponential map as an infinite series. The infinite series is consistent with our definition in terms of the eigenbasis. Consider the following,
$$\det( \exp(\epsilon A) ) = \exp(\epsilon \mathrm{Tr}(A))$$
$$\det( I + \epsilon A + \frac{\epsilon^2}{2} A^2 +\cdots ) = 1 + \epsilon \mathrm{Tr}(A) + \frac{\epsilon^2}{2} (\mathrm{Tr}(A))^2 + \cdots \qquad (*)$$
On the left hand side of $(*)$ we can factor an $\epsilon$ an write,
$$ \det( I + \epsilon \lbrace A + \frac{\epsilon}{2} A^2 +\cdots \rbrace ) = 1 + \epsilon \mathrm{Tr}(A + \frac{\epsilon}{2} A^2 + \cdots)+ \epsilon^2 f_2( A + \frac{\epsilon}{2} A^2 + \cdots ) + O(\epsilon^3)$$
$$ = 1 + \epsilon \mathrm{Tr}(A) + \epsilon^2 \lbrace \frac{1}{2}\mathrm{Tr}( A^2)+f_2( A + \frac{\epsilon}{2} A^2 + \cdots ) \rbrace + O(\epsilon^3)$$
Now we compare the second order terms in $\epsilon$ and obtain,
$$ \frac{1}{2}\mathrm{Tr}( A^2)+f_2( A + \frac{\epsilon}{2} A^2 + \cdots ) = \frac{1}{2}(\mathrm{Tr}(A))^2, $$
Now allow $\epsilon \rightarrow 0 $,
$$ \boxed{f_2( A ) = \frac{\mathrm{Tr^2}(A)-\mathrm{Tr}(A^2)}{2}} \qquad \textbf{(2)}. $$
The higher order terms can be obtained systematically using the same trick as above though the computations become more involved.
Solution 2:
The full formula for the expansion of the determinant of an nxn matrix A in a polynomial of traces of powers of A is: $$\det(A)=\sum_{\pi\in\Pi(n)}(-1)^{|\pi|-n}\prod_{m=1}^n\frac{(Tr(A^m))^{\pi_m}}{m^{\pi_m}\pi_m!}$$ where symbols are as follows: $\Pi(n)$ is the set of partitions of $n$ in positive integers; its elements are sequences of non negative integer multipliers, but a finite number of them different from zero, such that $$1\cdot \pi_1+2\cdot\pi_2+...n\cdot \pi_n=n.$$ E.g. the set of partitions of 3 has 3 elements, namely $$\Pi(3)=\{(3,0,\ldots),(1,1,0,\ldots), (0,0,1,0,\ldots)\},$$ meaning that the number 3 can be written as a sum of distinct positive integers with multiplicities in these ways: $$3=1\cdot3, 3=1\cdot 1+2\cdot 1, 3=3\cdot 1.$$ $|\pi|=\sum_{m=1}^{\infty}\pi_m$ is the length of the partition $\pi$, (the total number of non negative integers to be summed) e.g.$|(3,0,\ldots)|=3$, $|(1,1,0,\ldots)|=2$, etc. The formula can be derived, without assuming that A is diagonalizable, but of course assuming characteristic 0, by expanding $\ln(\lambda+A)=\ln\lambda +\ln(1+\frac{A}{\lambda})$ in a Taylor series, which converges for $\lambda>||A||$, then applying $\det e^B=e^{Tr(B)}$, expanding again the double series which is actually a polynomial because $\det(\lambda+A)$ is, and collecting terms of order 0 in $\lambda$. I believe there are other ways (i.e. combinatorial, using Newton's identities or something of the kind) to derive it, but I have not seen an explicit proof. The formula given in previous answer for a 2x2 matrix follows using the partitions of 2 $$\Pi(2)=\{(2,0,\ldots),(0,1,0,\ldots)\}$$ to compute $$\det(A)=(-1)^{2-2}\frac{(Tr(A))^2}{1^2\cdot 2!}+(-1)^{1-2}\frac{Tr(A^2)}{2^1\cdot 1!}.$$
As another example, using the partitions of 3 given above, one would find for a 3x3 matrix: $$\det(A)=(-1)^{3-3}\frac{(Tr(A))^3}{1^3\cdot 3!}+(-1)^{2-3}\frac{Tr(A)\cdot Tr(A^2)}{1^1\cdot1!\cdot 2^1\cdot 1!}+(-1)^{1-3}\frac{Tr(A^3)}{3^1\cdot1!}$$ or $$\det(A)=\frac{1}{6}\left((Tr(A)^3-3Tr(A^2)\cdot Tr(A)+2Tr(A^3)\right).$$
For a 4x4 matrix one ends up with $$\det(A)=\frac{(Tr(A))^4}{24}-\frac{(Tr(A))^2\cdot Tr(A^2)}{4}+\frac{(Tr(A^2))^2}{8}+\frac{Tr(A)\cdot Tr(A^3)}{3}-\frac{Tr(A^4)}{4}$$ and so on.
EDIT I'm editing my post, because I had forgotten to mention, that by the same token one can obtain a general formula for all the coefficients of the characteristic polynomial of A, not just the determinant. For an nxn matrix A, by defining $$\Delta_A(\lambda)=\det(\lambda+A)=\sum_{r=0}^{n}P_r(A)\lambda^r$$ (notice the unusual sign though) and expanding as above, one finds $$P_r(A)=\sum_{\pi\in\Pi(n-r)}(-1)^{|\pi|-n+r}\prod_{m=1}^{n-r}\frac{(Tr(A^m))^{\pi_m}}{m^{\pi_m}\pi_m!}$$ of which $\det(A)=P_0(A)$ is a particular case.