In a many-body quantum system, local operators in Heisenberg picture O(t) = e iHt Oe −iHt spread as time increases. Recent studies have attempted to find features of that spreading which could distinguish between chaotic and integrable dynamics. The operator entanglement -the entanglement entropy in operator space -is a natural candidate to provide such a distinction. Indeed, while it is believed that the operator entanglement grows linearly with time t in chaotic systems, we present evidence that it grows only logarithmically in generic interacting integrable systems. Although this logarithmic growth has been previously established for non-interacting fermions, there has been no progress on interacting integrable systems to date. In this Letter we provide an analytical upper bound on operator entanglement for all local operators in the "Rule 54" qubit chain, a cellular automaton model introduced in the 1990s [Bobenko et al., CMP 158, 127 (1993)], and recently advertised as the simplest representative of interacting integrable systems. Physically, the logarithmic bound originates from the fact that the dynamics of the models is mapped onto the one of stable quasiparticles that scatter elastically. The possibility of generalizing this scenario to other interacting integrable systems is briefly discussed.Understanding the out-of-equilibrium dynamics of isolated quantum many-body systems has been a prominent challenge since the early days of quantum mechanics [1].A key recurring idea is that, at long times, local properties are captured by statistical ensembles [1-4], despite the global dynamics being unitary. This suggests the possibility of a huge compression of information. In one dimension (1d) it implies that the reduced density matrix of a subsystem goes to a steady state well approximated by a Matrix Product Operator (MPO) [5-10]. This contrasts with the intermediate time behavior, where one faces an "entanglement barrier" [11,12] reminiscent of the generic linear growth of the entanglement entropy of a pure state after a quantum quench [13].In the late 2000s, the physical intuition that it could sometimes be more efficient to simulate the dynamics of operators -e.g. density matrices-rather than the one of pure states spurred another idea [8,[14][15][16][17]: that local observables in Heisenberg picture, O(t) = e iHt Oe −iHt , could also be approximated that way. In an insightful paper, Prosen andŽnidarič [8] observed numerically that there was a crucial distinction to be made between chaotic [18,19] and non-interacting dynamics: the bond dimension necessary for an MPO representation of O(t) was apparently blowing up exponentially with t in the former case and polynomially in the latter.