How to find zero-point oscillations for this system? Consider the following Hamiltonian which is absolutely relativistic literally: only sensitive to absolute pairwise relative phase space variables of objects for a system of $N$ objects moving in one dimension:
$$H = \sum_{ij} M_{ij} \sqrt{ (p_i - p_j)^2 + (q_i - q_j)^2 }$$ where  $(q_i,p_i)$ are phase space variables, $i\in\{ 1,\dots, N\}$, and $M_{ij}$ is a constant symmetric matrix. The stationary equilibrium configuration of this system is $q_i=a$ and $p_i=b$ for all $i$ and arbitrary $a$ and $b$. (Static equilibrium requires $b=0$.) How does the system evolve after perturbed infinitesimally from the equilibrium state?
Coincidentally, this is the Hamiltonian of Newtonian vector resonant relaxation of a thin stellar disk. 
The equations of motions are
$$\dot{q}_i = \frac{\partial H}{\partial p_i} = \sum_j M_{ij} \frac{p_i - p_j}{\sqrt{ (p_i - p_j)^2 + (q_i - q_j)^2 }}\,,$$ 
$$\dot{p}_i = -\frac{\partial H}{\partial q_i} = \sum_j M_{ij} \frac{q_j - q_i}{\sqrt{ (p_i - p_j)^2 + (q_i - q_j)^2 }}\,.$$ 
These equations are remarkably simple in the complex plane. Define $z_i = q_i + i p_i$, then
$$H = \sum_{nm} M_{nm} |z_n - z_m| \\
\dot{z}_n = \sum_m i M_{nm} \frac{z_m - z_n}{|z_m - z_n|}\,.$$
Notice that the "effective force" is independent of the phase space distance, it depends on only the argument of the relative coordinates. This is similar to the interaction of charged infinite insulating plates. 
Is there a way to derive the normal mode oscillations of this nonlinear system for "small" perturbations around $z_i=0$? 
(Note: The Hamiltonian is non-relativistic in the usual sense of the word. However this is the general relativistic Hamiltonian of 1 object moving if the metric of spacetime is $g_{\mu \nu}=x^2(-1,1,1,1)$. The Lagrangian is then $L = -\sqrt{-g_{\mu \nu}\frac{d x^{\mu}}{d\tau}\frac{d x^{\nu}}{d\tau}} =-\sqrt{x^2( 1 - \dot{x}^2)} $ which yields $H=\sqrt{x^2+ p^2}$.)
 A: It is an interesting problem. Usually you would find the infinitesimal oscillations by setting $p_i = b + \delta p_i$, $q_i = a + \delta q_i$ and expanding the Hamiltonian to second order in the $\delta p_i$, $\delta q_i$. Here though, this doesn't work as you just get the same Hamiltonian,
$$
H = \sum_{i,j} M_{ij} \sqrt{(\delta p_i - \delta p_j)^2 + (\delta q_i - \delta q_j)^2}
$$
I guess you already knew this. However, the result does tell you something important: the problem has a scale invariance. Mapping $p_i \to \alpha p_i$, $q_i \to \alpha q_i$ and $t \to \alpha t$ for scale factor $\alpha$ (where $t$ is time) reproduces exactly the same equations of motion.
If a problem has some natural length scale $l$ then you can look for infinitesimal solutions around the equilibrium: i.e. solutions where $\delta q_i \ll l$ for all $i$ and which are correct to leading order in $\delta q_i /l$. The normal modes of the system are infintesimal solutions which have a simple harmonic oscillator behaviour. The lack of a natural scale in this problem means it is not possible to talk of 'infinitesimal solutions', because there is no way to define what 'infinitesimal' means. For the same reason you cannot expect the system to have normal modes of oscillation.
It may still be possible to find approximate solutions, but it's hard to know what kind of approximation to make without knowing what kind of behaviour/questions you are interested in (for example, some forms of $M_{ij}$ may be easier than others).
[Edit:]
Note that the system has some conserved quantities:
$$
\sum_{n} z_n = \mathrm{constant}
$$
$$
\sum_{m , n} M_{mn}|z_n-z_m| = \mathrm{constant}
$$
and
$$
\sum_n |z_n|^2 = \mathrm{constant}
$$
The first two are related to the symmetries of translation in $z$ and in time, respectively. The third follows from the symmetry under $z_n \to z_n e^{i\theta}$. The scale invariance doesn't seem to have any connected conservation law.
Proof of third conservation law:
$\frac{d}{dt}\sum_n z_n {z_n}^\ast = \sum_n (\dot{z}_n {z_n}^\ast + z_n {\dot{z}_n}^\ast) = \sum_{n,m} i M_{mn} \frac{(z_m-z_n){z_n}^\ast - z_n({z_m}^\ast - {z_n}^\ast)}{|z_m-z_n|} = \sum_{n,m} i M_{mn} \frac{z_m{z_n}^\ast - z_n{z_m}^\ast}{|z_m-z_n|}$
The summand is anti-symmetric under permutation of $m$ and $n$, and therefore the total sum is zero.
A: This is as far as I got: The following are exact solutions for the phase space evolution of a system of $N$ objects. These solutions may be considered stationary equilibria. Is there a way to do perturbation theory relative to these solutions?
MarkA has shown that there are three conserved quantities
$$\sum_n z_n = \rm const\,,\quad \sum_n |z_n|^2 = {\rm const}, \quad \sum_{nm} M_{nm} |z_n - z_m| = {\rm const} $$
MarkA has also noticed a spatial scale invariance. The first two implies the motion is bounded, the scale invariance implies infinite long-range interactions, which suggest that symmetric lattice or symmetic polygon mass distributions may represent stable normal modes as in solid state physics.
First we consider the case of $M_{nm}=M$ for all $n\neq m$. In all cases the solutions exhibit rigid body rotation in phase space.
Equilateral polygon (vertices)
A particular solution is
$$z_n(t) = a + b \,e^{2\pi i n/N}e^{i \omega t} \quad n\in\{1,2,\dots, N\}$$
where $a$ and $b$ are arbitrary constants which respectively set the position of the center of mass and the size and initial orientation of the polygon and
$$\omega= \sum_{n=1}^{N-1}M \frac{1 - e^{2\pi i n/N}}{ | 1 - e^{2\pi i n/N} |} = \sum_{n=1}^{N-1} M \frac{1 - e^{2\pi i n/N}}{ 2 \sin (\pi n/N)  } = -i M\sum_{n=1}^{N-1} e^{\pi i n/N} = M\cot \frac{\pi}{2N}$$
We may add an object to the center of mass of the polygon, which would maintain this configuration. 
Line segment Initially equidistant objects along a line segment also exhibit stable rotations. The configuration spins around its origin with an angular velocity independent of spatial scale
$$z_n(t) = a + b \frac{n-1}{N-1} e^{i \omega t (n-N)/2} \quad n\in\{1,2,\dots, N\}$$
$$\omega =\sum_{n=1}^{N-1} M(n-1) = M \frac{N(N-1)}{2} $$
The angular velocity is always independent of the size of the perturbation $b$. Each objects behaves as a harmonic oscillator.
It is conceivable that the the uniform rotation of rigid objects form a wider class of solutions for which $M_{ij} = m_im_j$. Is this true?
Infinite systems
$N$ need not be finite in the previous examples. The $N\rightarrow \infty$ limit is a circular wire of circulating objects or a straight rod spinning around its center of mass.
Symmetric infinite symmetric lattices 
$$z_n = a n + b m \quad n {\rm ~ and ~}m\in Z$$
are static if $M_{nm}=M$ or more generally if $M_{nm}$ measures distance along the mesh with an arbitrary translation- and rotation-invariant metric. 
Inhomogeneous couplings
If $M_{nm}$ are finite block-diagonal matrices with all elements equal within a given block but different across different blocks, then we can superpose the previous solutions. The rigid body polygons associated with each block rotate independently around their center of mass. 
If $M_{nm}$ is "almost" block diagonal, with small but nonzero offdiagonal blocks (relative to the diagonal blocks), then the different polygons sufficiently far from each other in phase space will spin around their centres, and the centres will orbit one another slowly corresponding to the coupling coefficients in the offdiagonal blocks. Note that "sufficiently far" means a large angular diamater distance, so that the different objects within the two distinct structures are separated by approximately the same complex angle in phase space.
The weak off-diagonal coupling will distort the polygons on very long timescales.
