If I understand it correctly, you are asking how log likelihood in a multi-class classification problem relates to the cross entropy loss. So here is my try:
Assuming we have a multi class classification problem ($C$ different classes) where we estimated the conditional probabilities for each class given the data $\textbf{X} = \{ (\textbf{x}_i , y_i)_{i=1}^N \}$ (e.g., using a neural network) and where the classes are one-hot encoded:
$$
\hat{p}_i^{(k)} \in [0, 1] \quad \text{denotes the estimated probability of the } i\text{th sample for the } k\text{th class}\\
p_i^{(k)}\in \{0, 1\}\quad \text{denotes the true probability of the } i\text{th sample for the } k\text{th class}\\
y_i^{(k)} = p_i^{(k)} = 1_{y_{i}^{(k)}==1} \quad \text{hard labels, i.e.,} \begin{cases} 1 & \text{if }i \text{th sample belongs to } k\text{th class},\\ 0 &\text{else.}\end{cases}
$$
Note that $\hat{p}_i^{(k)} = \hat{p}_i^{(k)} (\textbf{x}_i, \boldsymbol{\theta})$ is actually a function that depends on parameters $\boldsymbol{\theta}$ (e.g., network weights) and the input data $\textbf{x}_i$. In the end, we want to optimize the parameters $\boldsymbol{\theta}$. For the sake of simplicity, we just write $\hat{p}_i^{(k)}$ (silently knowing that these estimated probabilities come from some kind of function approximator).
The idea of maximum likelihood (or rather the maximum likelihood estimator) is to find the optimal model parameters $\hat{\boldsymbol{\theta}}$ by maximizing a likelihood function which estimates the likelihood our model parameters $\boldsymbol{\theta}$ given the observed data $\textbf{X}$, i.e.,
$$
\hat{\boldsymbol{\theta}} = \arg\max_{\boldsymbol{\theta}} P \big(\boldsymbol{\theta} | \textbf{X} \big) \stackrel{Bayes}{=} \arg\max_{\boldsymbol{\theta}} \frac { P(\textbf{X} | \boldsymbol{\theta}) P (\boldsymbol{\theta}) } { P(\textbf{X}) } =
\arg\max_{\boldsymbol{\theta}} P\big( \textbf{X} | \boldsymbol{\theta}\big) P(\boldsymbol{\theta}),
$$
where $\hat{\boldsymbol{\theta}}$ is commonly referred to as maximum likelihood estimator, $P(\boldsymbol{\theta})$ can be used to encode prior knowledge (e.g., making some parameters more probable without any information about the data) and is therefore termed prior distribution. If we have no prior knowledge about the model parameters (which is usually the case when using neural networks), we can use a uniform distribution making each parameter vector equally likely, such that the maximum likelihood estimator reduces to
$$
\hat{\boldsymbol{\theta}} = \arg\max_{\boldsymbol{\theta}} P\big( \textbf{X} | \boldsymbol{\theta} \big) \stackrel{i.i.d}{=} \arg\max_{\boldsymbol{\theta}} \prod_{i=1}^N P\Big( (\textbf{x}_i, y_i) | \boldsymbol{\theta}\Big),
$$
where i.i.d. denotes the assumption that our observations are independent and identically distributed.
Lastly, we need to identify how we can express $P\Big( (\textbf{x}_i, y_i) | \boldsymbol{\theta}\Big)$ in terms of the true $\textbf{p}_i$ and estimated $\hat{\textbf{p}}_i$ probability vectors. Do we know something about the true probability distribution? .... Yes! It is a categorical distribution, i.e., each sample $\textbf{x}_i$ has exactly one class assigned (Thus, we have as as stated in the beginning $y_i^{(k)} = p_i^{(k)} = 1_{y_i^{(k)}==1}$). A very neat formulation for the likelihood of categorical distributions is as follows:
$$
P\Big( (\textbf{x}_i, y_i) | \boldsymbol{\theta} \Big) = \prod_{k=1}^C \left(\hat{p}_i^{(k)}\right)^{1_{y_i^{(k)}==1}}
$$
Let's take a moment to digest this formula. It's actually pretty easy what we are doing here, we want that the estimated conditional probabilities $\hat{p}_i^{(k)} \in [0, 1]$ correspond to the true labels $p_i^{(k)}\in \{0, 1\}$. Note that the true labels are hard labels, e.g., $\textbf{p}_0 = \begin{bmatrix}0 &1&0 \end{bmatrix}^{\text{T}}$, but the predictions are (probably) soft labels, e.g., $\hat{\textbf{p}}_0 = \begin{bmatrix}0.2 &0.7&0.1 \end{bmatrix}^{\text{T}}$. We optimize our likelihood by maximizing the predicted probability of the true label, i.e., where the label $y_i$ corresponds to the class $k$, this is denoted by the indicator function $1_{y_i^{(k)} == 1}$
The maximum likelihood estimator can be rewritten to
$$
\hat{\boldsymbol{\theta}} = \arg\max_{\boldsymbol{\theta}} P(\textbf{X}|\boldsymbol{\theta})= \arg\max_{\boldsymbol{\theta}} \prod_{i=1}^N \prod_{k=1}^C \left(\hat{p}_i^{(k)}\right)^{1_{y_i^{(k)}==1}} =
\arg\max_{\boldsymbol{\theta}} \prod_{i=1}^N \prod_{k=1}^C \left(\hat{p}_i^{(k)}\right)^{p_i^{(k)}}
$$
Now we take the negative log-likelihood (hence the maximization problem becomes a minimization problem):
$$
\hat{\boldsymbol{\theta}} = \arg\min_{\boldsymbol{\theta}} \Big(-\log P(\textbf{X}|\boldsymbol{\theta})\Big) = \arg\min_{\boldsymbol{\theta}} \left(- \sum_{i=1}^N \sum_{k=1}^C p_i^{(k)} \cdot \log \left(\hat{p}_i^{(k)}\right) \right)
$$
(This formula is the same as in your question).
Let's suppose, there is only one datapoint $N=1$, then we get
$$
- \sum_{k=1}^C p^{(k)} \cdot \log \left(\hat{p}^{(k)}\right) = H(p, \hat{p})
$$
This is known as the cross-entropy loss, i.e., minimization of the cross-entropy loss corresponds to maximum likelihood when hard labels are provided.
