Lyapunov Exponent of Nonlinear Schrödinger Operator

Question

I am interested in ideas for how to compute the maximal Lyapunov exponent of the nonlinear Schrödinger-Newton system given by $$ \partial_t s(t, \vec{x}) = L(s(t,\vec{x})) $$ where $$ L(s(t, \vec{x})) = i\left(\nabla^2 + \int_{\mathbb{R}^3} \frac{|s(t,\vec{y})|^2}{|\vec{x} - \vec{y}|} d^3\vec{y} \right) s(t, \vec{x}). $$ I have tried studying a trajectory $s'(t,\vec{x}) = s(t,\vec{x}) + \epsilon \eta(t,\vec{x})$ for small $\epsilon$ to see how $|\eta(t,\vec{x})|$ grows over time, but the calculations get rather intricate very fast. I am therefore looking for other ideas for how I might go about this.