Deriving the addition formula for the lemniscate functions from a total differential equation

The lemniscate of Bernoulli $C$ is a plane curve defined as follows.

Let $a > 0$ be a real number.
Let $F_1 = (a, 0)$ and $F_2 = (-a, 0)$ be two points of $\mathbb{R}^2$.
Let $C = \{P \in \mathbb{R}^2; PF_1\cdot PF_2 = a^2\}$.
Then the equation of $C$ in the polar coordinates is:

$r^2 = 2a^2\cos 2\theta$

Let $P$ be a point of $C$ in the first quadrant.
Let $u$ be the arc length between $O = (0, 0)$ and $P$.

Then, by this question,
$u = \int_{0}^{r} \frac{2a^2dr}{\sqrt{4a^4 – r^4}}$

Let $2a^2 = 1$ and $x = r$.

Then
$u = \int_{0}^{x} \frac{dx}{\sqrt{1 – x^4}}$

$u = u(x)$ is defined in $0 \le x \le 1$.
However, the above integral can be defined on $[-1, 1]$.
So we extend the domain of $u(x)$ to $[-1, 1]$ by the above integral.
Since $\frac{1}{\sqrt{1 – x^4}}$ is invariant under the substitution $x \rightarrow -x$,
$u(-x) = -u(x)$ for every $x \in [-1, 1]$.

Since $u'(x) = \frac{1}{\sqrt{1 – x^4}} > 0$ on $(-1, 1)$, $u(x)$ is strctly increasing on $[-1, 1]$.
Hence there exists the inverse function of $u(x)$. We denote the inverse function of $u(x)$ by $s(u)$. We call $s(u)$ lemniscate sine.
Since arcsin $x = \int_{0}^{x} \frac{dx}{\sqrt{1 – x^2}}$, $s(u)$ is analogous to $\sin u$.
We denote $u(1) = \int_{0}^{1} \frac{dx}{\sqrt{1 – x^4}}$ by $\omega$.
$s(u)$ is defined on $[-\omega, \omega]$.
$\omega$ corresponds to $\frac{\pi}{2}$ in the analogy of $s(u)$ with $\sin u$.

Since $u(-x) = -u(x)$, $s(-u) = -s(u)$

We define a function $c(u)$ by $c(u) = s(\omega – u)$ and call it lemniscate cosine.
$c(u)$ is defined on $[0, 2\omega]$.

Pursuing the analogy with $\sin u$ and being motivated by this question, we consider the following total differential equation.

$$\frac{dx}{\sqrt{1 – x^4}} + \frac{dy}{\sqrt{1 – y^4}} = 0$$

Let $u = \int_{0}^{x}\frac{dx}{\sqrt{1 – x^4}}$
Then $x = s(u)$

Let $v = \int_{0}^{y}\frac{dy}{\sqrt{1 – y^4}}$
Then $y = s(v)$

Let $c$ be a constant.
Then $u + v = c$ is a solution of this equation.
Then we get

$$s(u + v) = \frac{x\sqrt{1 – y^4} + y\sqrt{1 – x^4}}{1 + x^2y^2}$$

Substituting $u = \omega$, $v = -u$, we get $x = s(\omega) = 1, y = s(-u) = -s(u)$.

Hence $s(\omega – u) = \frac{\sqrt{1 – y^4}}{1+y^2} = \sqrt{\frac{1 – y^2}{1 + y^2}}$

Hence $c(u) = \sqrt{\frac{1 – s^2(u)}{1 + s^2(u)}}$

Hence

$$s(u+v) = \frac{s(u)c(v) + s(v)c(u)}{1 – s(u)s(v)c(u)c(v)}$$

Since $c(u+v) = s(\omega – u – v) = s((\omega – u) + (-v))$,

$$c(u+v) = \frac{c(u)c(v) – s(u)s(v)}{1 + s(u)s(v)c(u)c(v)}$$

Remark
Since $c(u) = \sqrt{\frac{1 – s^2(u)}{1 + s^2(u)}}$ and $s(-u) = -s(u)$, $c(-u) = c(u)$.

My question
How do we prove the following equation?

$$s(u + v) = \frac{x\sqrt{1 – y^4} + y\sqrt{1 – x^4}}{1 + x^2y^2}$$

Solutions Collecting From Web of "Deriving the addition formula for the lemniscate functions from a total differential equation"

The following proof is basically the same as the previous one, but some people may prefer this.

Let $u + v = c$, where $c$ is a constant.
It suffices to prove that

$$s(c) = \frac{s(u)c(v) + s(v)c(u)}{1 – s(u)s(v)c(u)c(v)}$$

Since $v = c – u$, the right hand side of this equation will be

$$f(u) = \frac{s(u)c(c – u) + s(c – u)c(u)}{1 – s(u)s(c – u)c(u)c(c – u)}$$

Suppose $f'(u) = 0$.
Then $f(u)$ is a constant.
Hence $f(u) = f(0)$.
On the other hand, since $s(0) = 0$ and $c(0) = s(\omega) = 1$, $f(0) = s(c)$.
Hence $f(u) = s(c)$ and we are done.

So it suffices to prove that $f'(u) = 0$.
This is only a tedious routine.

We follow the method of my answer to this question.

Let $u = \int_{0}^{x}\frac{dx}{\sqrt{1 – x^4}}$.
Then $x = s(u)$.

Let $v = \int_{0}^{y}\frac{dy}{\sqrt{1 – y^4}}$.
Then $y = s(v)$.

Let $c$ be a constant.
$u + v = c$ is a solution of the equation:

$$\frac{dx}{\sqrt{1 – x^4}} + \frac{dy}{\sqrt{1 – y^4}} = 0$$

It suffices to prove that $s(c) = \frac{x\sqrt{1 – y^4} + y\sqrt{1 – x^4}}{1 + x^2y^2}$.
Since $v = c – u$, the right hand side is a function of $u$.
We write this function by $\phi(u)$.
Namely, $$\phi(u) = \frac{x\sqrt{1 – y^4} + y\sqrt{1 – x^4}}{1 + x^2y^2}$$

Let us compute $\frac{d\phi}{du}$.

$\frac{dx}{du} = 1/\frac{du}{dx} = \sqrt{1 – x^4}$

$\frac{dy}{du} = -\frac{dy}{dv} = -1/\frac{dv}{dy} = -\sqrt{1 – y^4}$

$\frac{d^2x}{du^2}= \frac{d\sqrt{1 – x^4}}{du}\cdot\frac{dx}{du} = \frac{-2x^3}{\sqrt{1 – x^4}} \sqrt{1 – x^4} = -2x^3$

$\frac{d^2y}{du^2}= \frac{d^2y}{dv^2} = -2y^3$

Hence
$\frac{d\phi}{du} = \frac{d}{du}(\frac{-x\frac{dy}{du} + y\frac{dx}{du}}{1 + x^2y^2})$

Let $f = -x\frac{dy}{du} + y\frac{dx}{du}$

Let $g = 1 + x^2y^2$

Then $\frac{d}{du}(\frac{f}{g}) = \frac{\frac{df}{du}}{g} – \frac{f\frac{dg}{du}}{g^2}$

$\frac{df}{du} = -\frac{dx}{du}\frac{dy}{du} – x\frac{d^2y}{du^2} + \frac{dy}{du}\frac{dx}{du} + y\frac{d^2x}{du^2} = y\frac{d^2x}{du^2} – x\frac{d^2y}{du^2} = 2(xy^3 – yx^3)$

$\frac{dg}{du} = 2xy^2\frac{dx}{du} + 2x^2y\frac{dy}{du} = 2xy(y\frac{dx}{du} + x\frac{dy}{du})$

$f\frac{dg}{du} = 2xy(y^2(\frac{dx}{du})^2 – x^2(\frac{dy}{du})^2) = 2xy(y^2(1 – x^4) – x^2(1 – y^4))$

Hence
$\frac{d}{du}(\frac{f}{g}) = \frac{2(xy^3 – yx^3)}{1 + x^2y^2} – \frac{2xy(y^2(1 – x^4) – x^2(1 – y^4)}{(1 + x^2y^2)^2} = \frac{2(1 + x^2y^2)(xy^3 – yx^3) – 2xy(y^2(1 – x^4) – x^2(1 – y^4)}{(1 + x^2y^2)^2} = 0$

Hence
$\frac{d\phi}{du} = 0$

Hence $\phi(u)$ is constant.
Hence $\phi(u) = \phi(0) = y = s(v) = s(c)$ as desired