A limit with determinants
Here is a conceptual proof that avoids expanding a complicated determinant:
The determinant of a linear operator (or of a square matrix) is the product of the eigenvalues, counting multiplicity. The trace of a linear operator (or of a square matrix) is the sum of the eigenvalues, counting multiplicity.
Now suppose that $\lambda_1, \dots, \lambda_n$ are the eigenvalues of $H$, counting multiplicity. Then the eigenvalues of $I + \epsilon H$ are $$ 1 + \epsilon \lambda_1, \dots, 1 + \epsilon \lambda_n, $$ counting multiplicity. Thus $$ \det(1 + \epsilon H) = (1 + \epsilon \lambda_1) \cdots (1 + \epsilon \lambda_n). $$ It is now clear that \begin{align*} \lim_{\epsilon\to 0} \frac{ \det(1 + \epsilon H) - 1}{\epsilon} &= \lambda_1 + \dots + \lambda_n\\ &= \text{trace } H, \end{align*} as desired.
As I suggested in my comment, we proceed by expanding $$|Id + \varepsilon H| = |A| = \left| \begin{array}{cccc} 1+\varepsilon h_{11} & \varepsilon h_{12} & \cdots & \varepsilon h_{1n}\\ \varepsilon h_{21} & 1+\varepsilon h_{22} & \cdots & \ \\ \vdots & \ & \ddots & \ &\\ \varepsilon h_{n1} & \ & \ & 1+\varepsilon h_{nn} \end{array} \right|$$ in powers of $\varepsilon$. First, we have: $$|A| = (1+\varepsilon h_{11}) \left| \begin{array}{cccc} 1+\varepsilon h_{22} & \varepsilon h_{23} & \cdots & \varepsilon h_{2n}\\ \varepsilon h_{32} & 1+\varepsilon h_{33} & \cdots & \ \\ \vdots & \ & \ddots & \ &\\ \varepsilon h_{n2} & \ & \ & 1+\varepsilon h_{nn} \end{array} \right| + \varepsilon\sum_{j=2}^n (-1)^{1+j} h_{1j} \det(A_{1j})$$ where $\det(A_{1j}) = O(\varepsilon)$ (here, I am using the usual notation $A_{1j}$ to be the matrix obtained by deleting the first row and $j$th column from $A$).
Justification of this part: The minimal power of $\varepsilon$ would occur when the maximal number of diagonal terms is included. In this case (dealing with $(n-1)\times (n-1)$ minors, this would mean $n-2$ diagonal terms, since all terms of the cofactor expansion along the first row (with the exception of the first one, which I am separating from the rest of the computation) would exclude the $1+\varepsilon h_{11}$, as well as the diagonal entry that sits in the $j$th column. Finally, if there are $n-2$ diagonal terms multiplied together (in the minimal case), then there must be one off-diagonal term, introducing the (claimed) factor of $\varepsilon$.
We hence conclude that $$|A| = (1+\varepsilon h_{11}) \left| \begin{array}{cccc} 1+\varepsilon h_{22} & \varepsilon h_{23} & \cdots & \varepsilon h_{2n}\\ \varepsilon h_{32} & 1+\varepsilon h_{33} & \cdots & \ \\ \vdots & \ & \ddots & \ &\\ \varepsilon h_{n2} & \ & \ & 1+\varepsilon h_{nn} \end{array} \right| + O(\varepsilon^2)$$
Continuing (inductively) to expand the determinant in this manner, we see that \begin{align} |A| &= \prod_{j=1}^n (1+\varepsilon h_{jj}) + O(\varepsilon^2)\\ &= 1+\varepsilon \sum_{j=1}^n h_{jj} + O(\varepsilon^2) \end{align} which exactly yields the desired equality, since this immediately says $$|A|-1 = \varepsilon \sum_{j=1}^n h_{jj} + O(\varepsilon^2),$$ hence \begin{align} \lim_{\varepsilon \to 0} \frac{1}{\varepsilon}(|A|-1) &= \lim_{\varepsilon\to 0} \frac{1}{\varepsilon} \left( \varepsilon \sum_{j=1}^n h_{jj} + O(\varepsilon^2) \right)\\ &= \lim_{\varepsilon\to 0} \left( \sum_{j=1}^n h_{jj} + O(\varepsilon) \right)\\ &= \sum_{j=1}^n h_{jj} = \text{Trace}(H) \end{align}
Hint:
$$\left|\begin{matrix} 1+\epsilon h_{11}&\epsilon h_{12}\\ \epsilon h_{21}&1+\epsilon h_{22}\\ \end{matrix}\right|=1+\epsilon h_{11}+\epsilon h_{22}+\epsilon^2\left(h_{11}h_{22}-h_{12}h_{21}\right).$$
Only the main diagonal generates terms in $\epsilon$. This generalizes to higher order, for instance using the expansion by minors.