The function you are looking at is the LogSumExp function:
$$f(\mathbf{a}) = \ln \Big( \sum_{i=1}^n \exp(a_i) \Big)
\quad \quad \quad
\text{for all } \mathbf{a} \in \mathbb{R}^n.$$
Its gradient vector and Hessian matrix are given respectively by:
$$\begin{equation} \begin{aligned}
\nabla f(\mathbf{a})
&= \frac{1}{\sum_{i=1}^n \exp(a_i)} (\exp(a_1),...,\exp(a_n)), \\[12pt]
\nabla^2 f(\mathbf{a})
&= \text{diag}(\nabla f(\mathbf{a})) - \nabla f(\mathbf{a}) \nabla f(\mathbf{a})^\text{T}. \\[6pt]
\end{aligned} \end{equation}$$
(Here we have written the Hessian matrix in terms of the gradient vector. This is useful for the next step.) For any $\mathbf{z} \in \mathbb{R}^n$ we have the quadratic form:
$$\begin{equation} \begin{aligned}
\mathbf{z}^\text{T} (\nabla^2 f(\mathbf{a})) \mathbf{z}
&= \mathbf{z}^\text{T} \Big[ \text{diag}(\nabla f(\mathbf{a})) - \nabla f(\mathbf{a}) \nabla f(\mathbf{a})^\text{T} \Big] \mathbf{z} \\[6pt]
&= \mathbf{z}^\text{T} \text{diag}(\nabla f(\mathbf{a})) \mathbf{z}
- \mathbf{z}^\text{T} \nabla f(\mathbf{a}) \nabla f(\mathbf{a})^\text{T} \mathbf{z} \\[6pt]
&= \mathbf{z}^\text{T} \text{diag}(\nabla f(\mathbf{a})) \mathbf{z}
- (\nabla f(\mathbf{a}) \cdot \mathbf{z})^\text{T} (\nabla f(\mathbf{a}) \cdot \mathbf{z}) \\[6pt]
&= \mathbf{z}^\text{T} \text{diag}(\nabla f(\mathbf{a})) \mathbf{z}
- || \nabla f(\mathbf{a}) \cdot \mathbf{z} ||^2 \\[6pt]
&= \sum_{i=1}^n \bigg( \frac{\exp(a_i)}{\sum_{i=1}^n \exp(a_i)} \bigg) z_i^2
- \sum_{i=1}^n \bigg( \frac{\exp(a_i)}{\sum_{i=1}^n \exp(a_i)} \bigg)^2 z_i^2 \\[6pt]
&= \frac{1}{\sum_{i=1}^n \exp(a_i)} \sum_{i=1}^n \exp(a_i) z_i^2 \Bigg[ 1 - \frac{\exp(a_i)}{\sum_{i=1}^n \exp(a_i)} \Bigg] \\[6pt]
&= \frac{\sum_{i=1}^n \sum_{j \neq i} \exp(a_i) \exp(a_j) z_i^2}{(\sum_{i=1}^n \exp(a_i))^2} \geqslant 0. \\[6pt]
\end{aligned} \end{equation}$$
This establishes that the Hessian matrix is non-negative definite, which means that the LogSumExp function is (weakly) convex.