By making some account of astronomy I came across the parametric curve of the analemma:
$$ x(t) = \arcsin(s \sin(t)), y(t)= \arctan\left[\frac{(1-c)\tan(t)}{1+c \tan(t)^2}\right]$$ where $ s = \ sin (23,5 °) $, $ c = \ cos (23,5 °) $ (the sine and the cosine of the ecliptic angle) are constants, while the parameter $ t \in [0, 2 \pi] $. Now the symmetries with respect to the axes $ x, y $ are evident by virtue of the properties of the trigonometric functions. In addition, for the first derivative, we have:
$$ x'(t) = \frac{s \cos(t)}{\sqrt{1-s^2 \sin(t)^2}}, y'(t) = \frac{(1-c)(1-c\tan(t)^2)}{[(1+c\tan(t)^2)^2+(1-c)^2\tan(t)^2] \cos(t)^2}$$
The study of the y (x) function limited to the first quadrant $x,y>=0$ is immediate. In addition, the maximum $ t_0 = \arctan (1 / \sqrt {c}) $ is also well calculated. it is missing to study the curvature in the interval $ t \in [0, \pi / 2 [$, as well as the concavity of the function $ y (x) $ in the first quadrant. The negative curvature and therefore the concavity that gives us the classic eight form of the analemma should be proved. To derive the curvature with the formula of Serre is a difficult task. The same applies to the calculation of the second derivative of $ y (x) $. I tried to see the tangent expression of the tangent vector:
$$ y '(t) / x' (t) = \frac {(1-c) (1-c\ tan (t) ^ 2) \sqrt {1-s ^ 2 \sin (t) ^ 2} } {[(1 + c\tan (t) ^ 2) ^ 2 + (1-c) ^ 2 \tan (t)^ 2] s \cos (t)^3} $$
If I prove that the tangent, or $ \frac {y '} {x'} $, decreases , I show the negative curvature, the concavity. Numerator factors decrease, while denominator factors grow. There is only one problem: $ s \cos(t)^3$ decreases!!! Does anyone have any idea how to try the negative curvature, as well as the concavity?
The issue to me is that $x$ is the latitude north and $y$ is the longitude west of the sun from its mean equatorial position. That means that the coordinates describe a path on the unit sphere so I think a different metric may be in order. At least it's easier on the surface of the unit sphere!
Here's my analemma with the orbital elements taken from some JPL page: longitude of ascending node $\Omega=0$, longitude of perihelion $\varpi=102.99587119°$ thus for the sun $\omega=282.99587119°$, eccentricity $e=0.01670332$ and from Wikipedia $i=23.43695091°$ all for $2018$. If we point the $X$-axis at the first point of Ares, the $Z$-axis to celestial north and then get the $Y$-axis from $\hat Z\times\hat X=\hat Y$, then to locate a point on the analemma with true anomaly $\nu$, we start with $\hat X$, rotate it forward by $\nu+\omega$ about the $Z$-axis to get the angle from the ascending node, rotate forward by $i$ about the $X$-axis to get to where we should be on the ecliptic, the rotate backwards by $m+\omega$ about the $Z$-axis to get the position relative to the mean sun.
To get the eccentric anomaly $E$ from the true anomaly $\nu$ we use $$\begin{align}\sin E&=\frac{\sqrt{1-e^2}\sin\nu}{1+e\cos\nu}\\ \cos E&=\frac{\cos\nu+e}{1+e\cos\nu}\\ E&=\text{atan2}\left(\sin E,\cos E\right)\end{align}$$ And then we get the mean anomaly $m$ via $$m=E-e\sin E$$ Then we compute $$\begin{align}\begin{bmatrix}X\\Y\\Z\end{bmatrix}&= \begin{bmatrix}\cos(m+\omega)&\sin(m+\omega)&0\\-\sin(m+\omega)&\cos(m+\omega)&0\\0&0&1\end{bmatrix} \begin{bmatrix}1&0&0\\0&\cos i&-\sin i\\0&\sin i&\cos i\end{bmatrix} \begin{bmatrix}\cos(\nu+\omega)&-\sin(\nu+\omega)&0\\\sin(\nu+\omega)&\cos(\nu+\omega)&0\\0&0&1\end{bmatrix}\begin{bmatrix}1\\0\\0\end{bmatrix}\\ &=\begin{bmatrix}\cos(m+\omega)\cos(\nu+\omega)+\sin(m+\omega)\sin(\nu+\omega)\cos i\\ -\sin(m+\omega)\cos(\nu+\omega)+\cos(m+\omega)\sin(\nu+\omega)\cos i\\ \sin(\nu+\omega)\sin i\end{bmatrix}\end{align}$$ And the we can get $x=\pi/2-\theta=\pi/2-\cos^{-1}Z=\sin^{-1}Z$ and $y=-\phi=-\tan^{-1}(Y/X)$. Plot looks like this:
The original problem was easier in that $e=0$, so $\nu=E=m$ and $\nu+\omega=m+\omega=t$ so it works out to $$\begin{align}\begin{bmatrix}X\\Y\\Z\end{bmatrix}&=\begin{bmatrix}\frac12(1+\cos i)+\frac12(1-\cos i)\cos2t\\-\frac12(1-\cos i)\sin2t\\ \sin i\sin t\end{bmatrix}\\ \begin{bmatrix}\dot X\\ \dot Y\\ \dot Z\end{bmatrix}&=\begin{bmatrix}-(1-\cos i)\sin2t\\-(1-\cos i)\cos2t\\ \sin i\cos t\end{bmatrix}\\ \begin{bmatrix}\ddot X\\ \ddot Y\\ \ddot Z\end{bmatrix}&=\begin{bmatrix}-2(1-\cos i)\cos2t\\2(1-\cos i)\sin2t\\-\sin i\sin t\end{bmatrix}\end{align}$$ Here is a graph of the simpler analemma:
Now to get the sign of the eastward component of acceleration we just take $$\vec r\times\ddot{\vec r}\cdot\hat Z=\sin^2i\sin2t>0$$ For $0\lt t\lt\pi/2$. So maybe this wasn't what you wanted but it at least was more in line with what the curve physically represents as well as being more tractable.
EDIT: I'm not certain that I communicated my ideas sufficiently, so I'll make one last effort. The coordinate system of the original post is $$\begin{align}\vec r&=\langle x,y,z\rangle=\left\langle\sin^{-1}Z,-\tan^{-1}\left(\frac YX\right),0\right\rangle\\ &=\left\langle\sin^{-1}\left(\sin i\sin t\right),\tan^{-1}\left(\frac{(1-\cos i)\sin t\cos t}{\cos^2t+\cos i\sin^2t}\right)\right\rangle\\ &=\left\langle\sin^{-1}\left(\sin i\sin t\right),\tan^{-1}\left(\frac{\frac12(1-\cos i)\sin2t}{\frac12(1+\cos i)+\frac12(1-\cos i)cos2t}\right)\right\rangle\end{align}$$ Then $$\begin{align}\vec v&=\frac{d\vec r}{dt}=v\hat T=\left\langle\frac{\dot Z}{\sqrt{1-Z^2}},-\frac{\frac{\dot Y}X-\frac{Y\dot X}{X^2}}{1+\frac{Y^2}{X^2}},0\right\rangle\\ &=\left\langle\frac{\dot Z}{\sqrt{1-Z^2}},\frac{Y\dot X-X\dot Y}{X^2+Y^2},0\right\rangle\\ &=\frac1{1-Z^2}\left\langle\sqrt{1-Z^2}\dot Z,Y\dot X-X\dot Y,0\right\rangle\end{align}$$ We have use the fact that $\langle X,Y,Z\rangle$ lies on the unit sphere to pull out that factor of $(1-Z^2)$. Differentiating more, $$\begin{align}\vec a&=\frac{d\vec v}{dt}=a_T\hat T+a_N\hat N=a_T\hat T+\kappa v^2\hat N\\ &=\frac{2Z\dot Z}{(1-Z^2)^2}\left\langle\sqrt{1-Z^2}\dot Z,Y\dot X-X\dot Y,0\right\rangle\\ &\quad+\frac1{1-Z^2}\left\langle\sqrt{1-Z^2}\ddot Z-\frac{Z\dot Z^2}{\sqrt{1-Z^2}},Y\ddot X-X\ddot Y,0\right\rangle\end{align}$$ Now we can extract the curvature $\kappa$ via $$\begin{align}\vec v\times\vec a&=\kappa v^3\hat B\\ &=\frac1{(1-Z^2)^2}\left\langle0,0,\sqrt{1-Z^2}\dot Z(Y\ddot X-X\ddot Y)-\sqrt{1-Z^2}\ddot Z(Y\dot X-X\dot Y)+\frac{Z\dot Z^2(Y\dot X-X\dot Y)}{\sqrt{1-Z^2}}\right\rangle\\ &=\left\{\frac1{(1-Z^2)^{3/2}}\left[\dot Z(Y\ddot X-X\ddot Y)-\ddot Z(Y\dot X-X\dot Y)\right]+\frac1{(1-Z^2)^{5/2}}Z\dot Z^2(Y\dot X-X\dot Y)\right\}\hat k\end{align}$$ Considering that $$\begin{align}Y\dot X-X\dot Y&=\frac12(1-\cos i)^2+\frac12(1-\cos^2i)\cos2t\\ Y\ddot X-X\ddot Y&=-(1-\cos^2i)\sin2t\end{align}$$ The final expression for the curvature might not be too bad, but I won't pursue the issue further.