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:
% Analemma.m
clear all;
close all;
e = 0.01670332;
i = 23.43695091*pi/180;
omega = pi+102.99587119*pi/180;
nu = linspace(0,2*pi,300);
for k = 1:length(nu),
sinE = sqrt(1-e^2)*sin(nu(k))/(1+e*cos(nu(k)));
cosE = (cos(nu(k))+e)/(1+e*cos(nu(k)));
E = atan2(sinE,cosE);
m = E-e*sinE;
X = cos(m+omega)*cos(nu(k)+omega)+sin(m+omega)*sin(nu(k)+omega)*cos(i);
Y = -sin(m+omega)*cos(nu(k)+omega)+cos(m+omega)*sin(nu(k)+omega)*cos(i);
Z = sin(nu(k)+omega)*sin(i);
x(k) = asind(Z);
y(k) = -atand(Y/X);
end
plot(y,x);
set(gca,'Xdir','reverse');
title(['Analemma for e = ' num2str(e)]);
xlabel('Degrees West');
ylabel('Degrees North');
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.
