In the M-step of the EM algorithm, you have to maximize the expected log-likelihood of X with respect to z which is: $$\int d z P(Z \mid X, \theta^{old}) \ln P(X \mid Z, \theta)$$. Why do we maximize the expected log-likelihood and not the log of the expected likelihood which is $$\ln \left\{\int d z P(Z \mid X, \theta^{old}) P(X \mid Z, \theta)\right\}$$. Is it just an arbitrary choice for the algorithm or is there a deeper reason for that ?

Let us start with a more simple task: Given some data points $$x_1, x_2, ..., x_n$$ and given the assumption that they come from a normal distribution $$\mathcal{N}(\mu, \sigma)$$ with $$\sigma$$ known, let us try to figure out $$\mu$$. We assume that the $$x_i$$ were sampled independently (i.e. they are evaluations $$x_i = X_i(\omega)$$ for some 'true' $$\omega$$ and $$X_i \sim \mathcal{N}(\mu, \sigma)$$ and the $$X_i$$ are independent) write down the likelihood $$L(\mu) = p_\mu(x) = \prod_{i=1}^n p_\mu(x_i) = \text{const} \prod_{i=1}^n e^{(x_i - \mu)^2/2\sigma^2}$$
Why do people tend to put a log in front of that? The reason is simple: log has the magical property that $$\log(x*y) = \log(x) + \log(y)$$ i.e. it transforms products to sums. That is in fact the reason why we use the logarithm so often (and not sin, cos, any other crazy function: they all do not have this property). Why would one want that? The reason for that in turn is that almost the only optimization philosophy we have is: compute the gradient of the function at the current point and then step with the parameters into the direction of the gradient (because one can mathematically prove that at least locally the gradient points into the direction of the steepest ascend and we want to maximize). So let us compute the derivative of the likelihood... Uuh, there is an ugly product and the product rule $$(fg)' = f'g + fg'$$ makes it extremely difficult to take the derivative of the likelihood directly. How does it look like after taking the log? \begin{align} \log L(\mu) &= \text{const} + \sum_{i=1}^n \log e^{-(x_i - \mu)^2/2\sigma^2}\\ &= \text{const} + \text{const} \sum_{i=1}^n -(x_i - \mu)^2 \end{align}
and since we want to maximize we can ignore the constants, i.e. we have to maximize $$- \sum_{i=1}^n (x_i - \mu)^2$$ or minimize $$\sum_{i=1}^n (x_i - \mu)^2$$ and computing the gradient of that function is much more simple!
Now when dealing with the EM algorithm setup it is almost the same. We want to maximize (under the assumption that the random variables that sample $$(x_i, z_i)$$ are independent $$p_\Theta(x) = \int_{\mathcal{Z}_1} ... \int_{\mathcal{Z}_n} \prod_{i=1}^n p_\Theta(x_i, z_i) dz_1 ... dz_n = \prod_{i=1}^n \int_{\mathcal{Z}_i} p_\Theta(x_i, z_i) dz_i$$