How to calculate a delay between two signals in frequency domain?
Calculating a delay between two signals is not as simple as you described and it cannot be generalized as easily either.


  
  The case for two sinusoidal signals seems trivial, as the delay then agrees with the phase (divided by the angular frequency) of the Fourier component (of this pointwise multiplication) for the respective frequency.
  


This is not true in the general case. Calculating the delay from a single sinusoidal component is inherently ambiguous because of the signal's periodic nature. For a single phase $\phi$, there are infinitely many valid delays $\tau$.

$$
\tau = \frac{2\pi \cdot n + \phi}{2\pi f}, n \in \mathbb{Z}
$$


  
  This can be generalized to two identic signals with a fixed delay, as the phase response equals then the angular frequency multiplied by the delay, for all Fourier components.
  


This isn't true either. As noted above, the phase response plus a multiple of 2π equals the angular frequency multiplied by the delay, for all Fourier components. Even worse, this is a different multiple for each Fourier component in the general case.

You want to find the delay $\tau$ that fulfills this equation for all Fourier components, or at least for most of them, if the signals are not exacly identical.

The inverse fourier transform solves this problem in a brute-force manner. For each candidate delay, it calculates how strong the measured phase response agrees with it. Just pick the maximum and you're done. There is no better algorithm for the general case because the problem is inherently ambiguous and many (or even all) solutions might agree equally strong. Obviously, there cannot be a method to calculate the solution of such a problem directly.

Things change if you want to calculate the delay with higher precision than one sample. You could just zero-pad the phase response before calculating the inverse fourier transform, but this is computationally inefficient. A better approach is to calculate a coarse solution first, using an unpadded inverse fourier transform and then refine it using Newton's algorithm. Newton's algorithm is perfecly fine here because


it never requires the calculation of a delay from a phase response (which would be ambiguous). You can get away with calculating agreement scores for a given phase response and a candidate delay.
convergence to a local maximum is perfectly fine if all you wanted to do is to refine the global maximum.


There is no non-brute-force algorithm to find the global maximum of an arbitrary function. You are trying to find the global maximum of the cross correlation function, so unless you can guarantee certain properties of the signals (such as convexity or the like) or you're happy with a local maximum, there will be no way around an inverse fourier transform.