What do physicists mean by an "integrable system"?

I don't know if I can state a clear single "definition", but hopefully I will be able to sort out some of the concepts and the confusion.

integrability is sometimes associated with having a closed form solution

This, I think, is categorically not true. At least in the usual sense of 'closed form'. If you take the Lieb-Liniger model, which is I believe one of the seminal examples of an integrable system, the solution obtained is in the form of a set of integral equations, that the authors proceed to solve numerically. This is not 'closed form'.

integrability is sometimes associated with having infinitely many conserved quantities

This is the definition I am familiar with, but it requires caution and there are some subtleties. Namely, every system at the thermodynamic limit has an infinite number of conserved quantities: the projectors onto the eigenstates of the Hamiltonian $|psi_nrangle langle psi_n |$. Therefore, this definition alone is not enough. One needs an infinite number of conserved quantities that are 'not trivial' in some sense. Sometimes they are defined by being with local support, but I am not sure that this is enough or unique. However, it usually guaranteed that if one has a solution of the system in terms of the 2-particle scattering matrix and the associated Yang-Baxter equation, one can construct this infinite number of conserved quantities.

integrability is sometimes kind of like the opposite of chaos

integrability is sometimes kind of like the opposite of thermalization

These two are related, as I understand them, and the notion is generally derived from the existence of the infinite number of conserved quantities. The idea is that if we have an infinite number of 'non trivial' conserved quantities, then we can describe the macroscopic observables using them, and then the observables keep their value throughout the time-evolution. This, of course, contradicts thermalization and chaos, in the sense that if a system is prepared in some state it will keep its initial observables, instead of thermalizing. However, this is a subject of a very lively debate, surrounding the questions of what is exactly the nature of the the conserved quantities, whether or not the 'eigenstate thermalization hypothesis' is true or not, and how can one generalize integrability to 'quasi-integrable' models.

I think that as in many other topics in contemporary physics, there is no clear definition of integrability. Once it was related to a system having an exact solution (usually via the Bethe-ansatz method or one of its relatives), and the infinite number of conserved quantities was a feature / definition depending on your point of view. Nowadays the term migrated and expanded, together with the interests of the community.

I am a bit late to the party, but I have had similar questions to yours in the past. I will summarise below what I know, which has been able to "quell" my dissatisfactions about integrable systems for the time being. Maybe it's just a placebo though...

Definition (at least one that I like)

A system with $n$ degrees of freedom, and $n$ constants of motion, such that the Poisson bracket of any pair of constants of motion vanishes, is known as a completely integrable system. Such a collection of constants of motion are said to be in involution with each other.

Completely integrable is opposed to just "partially" integrable, meaning that you cannot get a full analytical solution, and to superintegrable which are systems with $>n$ constants of motion. For example orbital motion is superintegrable (when neglecting inter-planetary interactions, or 'perturbations') because you not only have energy and angular momentum, but also the Runge-Lenz vector as conserved quantities..

For $F$ to be a constant of motion, you need ${F, H} = 0$, while for two constants of motion to be in involution you need ${F_1, F_2} = 0$. This essentially means that any two constants of motion are "compatible" with one another and can be employed at the same time. An example where this is not true is, jumping to quantum mechanics and hence replacing Poisson brackets with commutation relations, angular momentum. While each component of the angular momentum operator commutes with the Hamiltonian $[ hat L_i , hat H] = 0$, they are not in involution with each other as $[ hat L_i , hat L_j] = epsilon_{ijk} hat L_k$. Hence you cannot "use" all three of them at the same time, but usually opt for the conventional $hat{mathbf{L}}^2$ and $hat L_z$.

Practical definition

Integrable systems are nonlinear differential equations which ‘in principle’ can be solved analytically. This means that the solution can be reduced to a finite number of algebraic operations and integrations.

This definition is taken from here which also makes an interesting distinction between ODEs and PDEs, arguing that integrability is not fully defined for the latter.

An example would be a $1D$ system with Hamiltonian $H(p,q) = p^2/2 + V(q)$, obeying the usual Hamilton's equations $dot q = p, dot p = -mathrm{d}V/mathrm{d}q$. Using the conserved quantity (energy) $E = p^2/2 + V(q)$ coming out of ${H,H} =0$, you can write $p = pm sqrt{2(E-V(q))}$ and hence: $$ t = pm frac{mathrm{d}q}{sqrt{2(E-V(q))}},$$ which then you can invert to find $q(t)$. I only have one integral relating $q$ and $t$. In an $n$-dimensional system $q_1, q_2,... q_i$, if one can still break down the problem to $n$ integrals each involving a specific $q_i$ and $t$ only, then one would have fully integrated the system without needing to solve coupled differential equations simultaneously (and very probably numerically).

And I think this qualifies as an answer to why these systems are associated with being nice.

Closed form solution

Apart from having an 'analytical' solution, meaning you can write $N$ equations each with only one $q_i$ as outlined above (as opposed to $N$ coupled equations), the term 'closed form solution' may have a geometrical interpretation.

Given a system with variables $q_1, q_2, ... q_i, ... q_n$ and $p_1, p_2, ... p_i, ... p_n$, the parameter space is $mathbb{R}^{2n}$.

The actual trajectory that solutions to the equations of motion $(q_s,p_s)$ will follow live in the phase space $U subseteq mathbb{R}^{2n}$. For instance a 2D harmonic oscillator may give you closed Lissajous figures when $omega_1/omega_2 in mathbb{Q}$ (hence $U subsetmathbb{R}^{2n}$), but these curves will fill the whole space densely for $omega_1/omega_2 not in mathbb{Q}$ (hence $U = mathbb{R}^{2n}$).

A more topological explanation for this involves lines on invariant tori and Hopf fibration, which I don't know enough about.

Infinite conserved quantities

Take $N$ free particles in a closed container that are not interacting among each other, but only bounce back from the container walls. Each particle conserves momentum (and energy, though this is trivial in the absence of a potential). Make $Nrightarrow infty$, and you have infinitely many conserved quantities.

The system might look like a mess since you have a zillion particles going all over the place, but each particle is doing its own thing, following an equation of motion that is independent of what the other particles are doing.

Given the initial conditions, + the conservation of each particle's momentum, the system is completely integrable.

Thermalisation

Continuing the $N$ particle example from above. If the system is integrable, then it cannot thermalise.

Thermalisation means that, eventually, the velocity/momentum/energy distribution tends to a Maxwell-Boltzmann profile (or whatever for quantum degenerate gases). This can only be true if particles are allowed to interact (at least a bit) so as to "redistribute" the momentum. Total energy (and total momentum) will still be conserved, but by allowing inter-particle interactions you are not enforcing the $N$ ($rightarrow infty$) conservations of each momentum.

You might have heard of a phenomenon called Many-Body Localisation (MBL), where a (quantum) system in the presence of weak disorder seems to remain localised despite the presence of interactions, thereby not reaching thermal equilibrium. This is connected ot the Eigenstate Thermalisation Hypothesis (mentioned in the other answer) which seems to tbe one of the few criteria to classify this MBL phase, as you can't use any symmetries and stuff since it's not an equilibrium state.
Well, a mathematical model to explain how this might be possible (some maths for this in presented in this review) assumes a set of localised conserved charges that are constants of motion for the system, sometimes referred to as $ell$-bits, making the system (“locally”) integrable.

Relation to chaos

This question here is interesting about this, but to be honest I have yet to understand the precise definition of chaos. Sometimes it seems to be a deterministic system that heavily depends on initial conditions, while some other times it means that the 'approximate present cannot predict the future'.

I would like to add one "definition" of integrable quantum systems, which has not been mentioned above, and it is related to the statement

• integrability is sometimes kind of like the opposite of chaos.

Integrable quantum systems have Poissonian level statistics. Chaotic quantum systems usually obey random matrix theory which means that the many-body levels exhibit level repulsion -- it is very unlikely to find two eigenvalues of the Hamiltonian near each other. In the case of integrable systems degeneracies are allowed due to the large number of conserved quantities. Therefore levels do not repel but are distributed independently of each other in the spectrum. The level statistics is then Poissonian, similar to what you would find in a non-interacting system.