The ellipse case is exactly the Klein model of hyperbolic geometry

Along a diameter. Take $x=0$ and $y=(r,0,\dots,0)$ with $0<r<1$. The line through them meets $\partial\mathbb{B}$ at $a=(-1,0,\dots,0)$ and $b=(1,0,\dots,0)$, so $|ay|=1+r$, $|ax|=1$, $|bx|=1$, $|by|=1-r$, giving $d_{\mathbb{B}}(0,y)=\frac12\ln\frac{1+r}{1-r}=\operatorname{artanh}(r)$. This matches exactly the standard formula for hyperbolic distance from the center to a point at Euclidean radius $r$ in the Beltrami–Klein model of curvature $-1$.

Off a diameter. For general $x,y\in\mathbb{B}$, choose a rotation $R\in O(n)$ (an isometry of the round ball fixing $0$) carrying the line through $x,y$ to a coordinate axis; since $R$ preserves $\mathbb{B}$ and maps the four points $a,x,y,b$ to another quadruple in the same cross-ratio, $d_{\mathbb{B}}(x,y)=d_{\mathbb{B}}(Rx,Ry)$, reducing to the diametral case worked out above after also translating one point to $0$ via a hyperbolic isometry of $\mathbb{B}$ (a Möbius transformation preserving $\mathbb{B}$, which likewise preserves cross-ratios). Since the formula obtained this way is exactly the Beltrami–Klein distance formula, $(\mathbb{B},d_{\mathbb{B}})$ is isometric to hyperbolic $n$-space.

General ellipsoids. Let $T:\mathbb{R}^n\to\mathbb{R}^n$ be an invertible linear (or affine) map with $T(\mathbb{B})=E$. An affine map is a projective transformation of $\mathbb{RP}^n$ that fixes the hyperplane at infinity, and cross-ratios of collinear points are invariant under every projective transformation; since $T$ sends the line through $x,y\in\mathbb{B}$ to the line through $Tx,Ty\in E$, and sends the boundary pair $a,b$ to the boundary pair $Ta,Tb$ of $E$, it sends the defining quadruple to the defining quadruple, so $[a,x,y,b]=[Ta,Tx,Ty,Tb]$ and hence $d_E(Tx,Ty)=d_{\mathbb{B}}(x,y)$ for all $x,y$. Thus $T:(\mathbb{B},d_{\mathbb{B}})\to(E,d_E)$ is an isometry, so $(E,d_E)$ is again hyperbolic $n$-space.