Alternate direction methods are based on the Lie formula (e(A/n)n e(B/n))n - e(A+B) = O(1/n) for square matrices A and B. A more precise formula is known: (e(A/2n)e(B/n)e(A/2n))n - e(A+B) = O(1/n2). The search for more accurate product formulae is negative, in die following sense: there is no integer k and no choice of nonnegative real numbers alpha(j) and beta(j), for 1 less-than-or-equal-to j less-than-or-equal-to k such that (e(alpha1A/n)e(beta1B/n) ... e(alphakA/n)e(betakB/n))n - e(A+B), is of order 3, unless the commutators [A, [A, B]] or [B, [B, A]] are linearly dependent. Replacing exponentials by approximation of exponentials does not improve the situation. However, linear combinations of product formulae can be constructed; they lead to stable numerical schemes of order 3 and more. Applications to domain decomposition methods are given, as well as some numerical experiments.