Optimizing the cross-entropy is equivalent to optimizing the log-likelihood of the parameters given the data, $\ell(\theta)$, which is what we want, i.e. find the parameters that most likely generated the data.
So, the likelihood is defined as $$\mathcal{L}(\theta) = P(y \mid x; \theta),$$
i.e. a function of the parameters $\theta$.
The log-likelihood is just the logarithm of the likelihood
$$\ell(\theta) = \log \mathcal{L}(\theta)$$
To understand why we take the logarithm, read this answer.
Now, for a fixed $\hat{\theta}$, $\mathcal{L}(\hat{\theta}) = P(y \mid x; \hat{\theta})$ is a conditional probability distribution.
For simplicity, let's assume that $\mathcal{L}(\hat{\theta}) = P(y \mid x; \hat{\theta})$ is a (conditional) Bernoulli distribution, which is defined by a single parameter $\hat{p}$, whose probability mass function (pmf) is defined as
$$\hat{p}^y (1-\hat{p})^{1-y},$$
where $y = \{0, 1\}$.
Now, in the context of neural networks, we use a neural network to output $\hat{p}$, i.e. an estimate of the parameter of the Bernoulli, which we plug into the cross-entropy loss function. So, if $f_\theta: \mathcal{X} \rightarrow [0, 1]$ is your neural network, then $f_{\hat{\theta}}(x) = \hat{p}$, given the input-output pair $(x, y)$.
So, the Bernoulli pmf for a given pair $(x, y) \in \mathcal{D}$ can be written as
\begin{align}
P(y \mid x; \hat{\theta})
&= \hat{p}^y (1-\hat{p})^{1-y} \\
&= f_{\hat{\theta}}(x)^y (1-f_{\hat{\theta}}(x))^{1-y}
\end{align}
So, without fixing the parameters, we can write the likelihood (which is not a probability distribution wrt $\theta$) as follows
\begin{align}
\mathcal{L}(\theta)
&= f_{\theta}(x)^y (1-f_{\theta}(x))^{1-y} \\
\iff \ell(\theta) &= \log \left( f_{\theta}(x)^y (1-f_{\theta}(x))^{1-y}\right) \\
&= \log f_{\theta}(x)^y + \log \left(1-f_{\theta}(x)\right)^{1-y} \\
&= y \log f_{\theta}(x) + (1-y)\log \left(1-f_{\theta}(x)\right)
\end{align}
which is our well-known binary cross-entropy for a single pair $(x, y)$ (note that only one of the addends above is not zero). For multiple pairs, you can apply the same reasoning, and you will see why we also use the log, i.e. it will allow you to simplify a few expressions in the derivation (specifically, multiplications will turn into sums, which is nice numerically).
So, if we maximize $\mathcal{L}$, with respect to $\theta$ (e.g. the parameters of your neural network), then we will find an estimate $\hat{\theta}$, that produces the parameter $\hat{p}$, such that this parameter can be used to construct/define the Bernoulli distribution that most likely (hence the name likelihood) produced our given dataset $\mathcal{D}$.
This reasoning also applies to the case where $P$ is categorical or even Gaussian distribution. So, actually, when you're minimizing the MSE, you're actually also optimizing the log-likelihood. It just turns out that the expression is different because Gaussians are not Bernoullis.
I will not talk about the back-propagation, but it's just the chain rule. If the log-likelihood is differentiable, then you can use it in the context of back-propagation and GD. You can ask another question, if you're interested, but note that, in the general case, computing the partial derivatives can be a pain in the neck, so your question may be too broad. Find a very simple case.
To conclude, don't dwell on the word error. Think in terms of probabilities, probability distributions, likelihoods, and why we can optimize likelihoods. Read my answer multiple times, if you didn't understand something.
