diff --git a/hamilflow/models/harmonic_oscillator_chain.py b/hamilflow/models/harmonic_oscillator_chain.py
index 89e44eb..da3ece9 100644
--- a/hamilflow/models/harmonic_oscillator_chain.py
+++ b/hamilflow/models/harmonic_oscillator_chain.py
@@ -20,22 +20,6 @@ class HarmonicOscillatorsChain:
 
     This system can be solved in terms of _travelling waves_, obtained by discrete Fourier transform.
 
-    We can complexify the system
-    $$S_L[x_i] = S_L[x_i, \phi_j] \equiv S_L[X^\ast_i, X_j] = \int_{t_0}^{t_1}\mathbb{d} t \left\\{ \frac{1}{2}m \dot X^\ast_i \delta_{ij} \dot X_j - \frac{1}{2}m X^\ast_i A_{ij} X_j\right\\}\\,,$$
-    where $A_{ij} / \omega^2$ is equal to $(-2)$ if $i=j$, $1$ if $|i-j|=1$ or $|i-j|=N$, and $0$ otherwise;
-    $X_i \coloneqq x_i \mathbb{e}^{-\phi_i}$, $X^\ast_i \coloneqq x_i \mathbb{e}^{+\phi_i}$.
-
-    $A_{ij}$ can be diagonalised by the inverse discrete Fourier transform
-    $$X_i = (F^{-1})_{ik} Y_k = \frac{1}{\sqrt{N}}\sum_k \mathbb{e}^{i \frac{2\mathbb{\pi}}{N} k\mathbb{i}} Y_k\\,.$$
-
-    Calculating gives
-    $$S_L[X^\ast_i, X_j] = S_L[Y^\ast_i, Y_j] = \sum_{k=0}^{N-1} \int_{t_0}^{t_1}\mathbb{d} t \left\\{ \frac{1}{2}m \dot Y^\ast_k \dot Y_k - \frac{1}{2}m \omega^2\cdot4\sin^2\frac{2\mathbb{\pi}k}{N} Y^\ast_k Y_k\right\\}\\,.$$
-    Using the same transformation to separate the non-dynamic phases, we can arrive at a real action
-    $$S_L[y] = \sum_{k=0}^{N-1} \int_{t_0}^{t_1}\mathbb{d} t \left\\{ \frac{1}{2}m \dot y_k^2 - \frac{1}{2}m \omega^2\cdot4\sin^2\frac{2\mathbb{\pi}k}{N} y_k^2\right\\}\\,.$$
-
-    The origional system can then be solved by $N$ independent oscillators
-    $$\dot y_k^2 + 4\omega^2\sin^2\frac{2\mathbb{\pi}k}{N} y_k^2 \equiv 4\omega^2\sin^2\frac{2\mathbb{\pi}k}{N} y_{k0}^2\,.$$
-
     Since the original degrees of freedom are real, the initial conditions of the propagating waves need to satisfy
     $Y_k = Y^*_{-k \mod N}$, see [Wikipedia](https://en.wikipedia.org/wiki/Discrete_Fourier_transform#DFT_of_real_and_purely_imaginary_signals).
     """
@@ -46,6 +30,14 @@ def __init__(
         initial_conditions: Sequence[Mapping[str, float | tuple[float, float]]],
         odd_dof: bool,
     ) -> None:
+        """Instantiate an oscillator chain.
+
+        :param omega: frequence parameter
+        :param initial_conditions: a sequence of initial conditions on the Fourier modes.
+            The first element in the sequence is that of the zero mode, taking a position and a velocity.
+            Rest of the elements are that of the independent travelling waves, taking two amplitudes and two initial phases.
+        :param odd_dof: The system will have `2 * len(initial_conditions) - int(odd_dof)` degrees of freedom.
+        """
         self.omega = omega
         self.n_independant_csho_dof = len(initial_conditions) - 1
         self.odd_dof = odd_dof