This problem (or rather, statement that I cannot understand) has arisen in a paper I have been reading "Geometry of Integrable Billiards and Pencils of Quadrics" by Dragovic and Radnovic. I'd be most grateful for any explanations of it (it may be a simple fact, but I'm not sure).

Let $\Omega \subset \mathbb{R}^d$ be a bounded domain such that its boundary $\partial \Omega $ lies in the union of several quadrics from the (confocal) family $ \mathcal{Q}_{\lambda}: Q_{\lambda}(x)=1 $ where $$Q_{\lambda}(x)=\sum_{i=1}^{d}\frac{x_{i}^2}{a_{i}-\lambda}.$$ Then in elliptic coordinates, $\Omega$ is given by: $$\beta_{1}'\leq\lambda_{1}\leq\beta_{1}'', \ldots, \beta_{d}'\leq\lambda_{d}\leq\beta_{d}'' $$ where $a_{s+1}\leq \beta_{s}'\leq\beta_{s}''\leq a_{s}$ for $1\leq s \leq d-1$ and $- \infty < \beta_{d}'<\beta_{d}''\leq a_{d}.$

Define $P(x):= (a_1 -x)\ldots(a_d -x)(\alpha_{1} -x)\ldots(\alpha_{d} - x).$

Now, we consider a billiard system inside $\Omega$ with caustics $\mathcal{Q}_{\alpha_1}, \ldots, \mathcal{Q}_{\alpha_d-1}.$

Why does the system of equations: $$ \sum_{s=1}^{d}\frac{d\lambda_s}{\sqrt{P(\lambda_s)}}=0, \sum_{s=1}^{d}\frac{\lambda_{s}d\lambda_{s}}{\sqrt{P(\lambda_s)}}=0, \ldots, \sum_{s=1}^{d}\frac{\lambda_{s}^{d-2}d\lambda_{s}}{\sqrt{P(\lambda_s)}}=0,$$

(which are apparently due to Jacobi and Darboux - I'd appreciate a modern reference because the only version I can find is scanned page-by-page in German), where $\sqrt{P(\lambda_s)}$ is taken with the same sign in all expressions, represent a system of differential equations of a line tangent to all the caustics $\mathcal{Q}_{\alpha_1}, \ldots, \mathcal{Q}_{\alpha_d-1}$? Moreover, why does: $$\sum_{s=1}^{d}\frac{\lambda_{s}^{d-1}d\lambda_{s}}{\sqrt{P(\lambda_s)}}=2dl$$ where $dl$ is an element of ``the" line length?

I found similar looking equations on page 4 of another paper (by Buser and Silhol), but cannot understand them either.