The problem is to find the second order term in the series expansion of the expression $\mathrm{det}( I + \epsilon A)$ as a power series in $\epsilon$ for a diagonalizable matrix $A$. Formally, we will write the series as follows
$$ \det( I + \epsilon A ) = f_0(A) + \epsilon f_1(A) + \epsilon^2 f_2(A) + \cdots +\epsilon^N f_N(A), $$
In particular, we are looking for an expression in terms of the trace of the matrix $A$.
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.