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}$$
 A: 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
A: 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.
