I'm going to assume you have some familiarity with linear algebra, as the math becomes much less tedious than trying to do three-dimensional trigonometry with the $x$, $y$, $z$-coordinates directly.
You are looking for a function describing the line of the horizon. Since it is a circle and thus a one-dimensional object, I'm going to call it $\vec{h}(t)$, where $t$ is a parameter variable such that $0 \leq t < 2\pi$. What I'm going to do is find the transformations that turn this horizon circle into a unit circle on the x-y plane centered at the origin (which I already know how to do), then apply those transformations in reverse.
First, define $\vec{s}$ as the vector indicating the position of the satellite with respect to the center of the Earth--the line of length $R+h$ in your diagram). It has components
$$\vec{s} = \begin{bmatrix}s_x\\s_y\\s_z\end{bmatrix}.$$The horizon circle $\vec{h}(t)$ is centered on the line from the center of the Earth ($\vec{0}$ in ECEF) to the satellite. This line is also perpendicular to the plane that contains the circle.
Imagine a line perpendicular to $\vec{s}$ that connects the horizon point. This line is a radius of the horizon circle, and where this line meets $\vec{s}$ is the center of the circle $\vec{h}(t)$. I'm going to call this point $\vec{h_c}$. To find out this quantity, I'm going to define the angle $\beta$ as the angle between the horizon point and the satellite, centered at Earth's center. It's equal to $90^o-\alpha$ ($\alpha$ should be given by $\sin^{-1}(d/(R+h))$) and calculated by
$$\beta = \cos^{-1}\left(\frac{R}{R+h}\right) = \cos^{-1}\left(\frac{R}{\left|\vec{s}\right|}\right).$$ See the picture below.
The distance from Earth's center to $\vec{h_c}$ is
$$\left| \vec{h_c} \right| = R \cos\beta.$$
So, the position of the center of $\vec{h}(t)$ is
$$\vec{h_c} = \left(R \cos\beta \right)\hat{s},$$
where $\hat{s}$ is the unit vector in the direction of $\vec{s}$. We can also immediately see that the radius of the horizon circle is $R \sin\beta$.
Now, we can turn the unknown $\vec{h}(t)$ into something simple.
$$\vec{h}(t) - \vec{h_c}$$ is the horizon circle moved so that it is centered at the center of the Earth.
$$\frac{\vec{h}(t) - \vec{h_c}}{R\sin\beta}$$
shrinks the horizon line so it has a radius of one.
Now, we need a transformation that rotates this unit circle onto the $x$-$y$ plane so we can write that circle as
$$\begin{bmatrix}\cos{}t\\\sin{}t\\0\end{bmatrix},$$ where this $t$ is the same as $\vec{h}(t)$. In other words, we need an orthogonal transformation $T$ such that
$$T \begin{bmatrix}0\\0\\1\end{bmatrix} = \hat{s}$$
I'm going to cheat a little bit here and refer to another of my answers where I derive one possible $T$.
$$T = \begin{bmatrix}
\frac{-\hat{s}_y}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \frac{-\hat{s}_x \hat{s}_z}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \hat{s}_x\\
\frac{\hat{s}_x}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \frac{-\hat{s}_y \hat{s}_z}{\sqrt{\hat{s}_x^2 + \hat{s}_y^2}} & \hat{s}_y\\
0 & \sqrt{\hat{s}_x^2 + \hat{s}_y^2} & \hat{s}_z
\end{bmatrix}$$
(Notice that this is actually $T^{-1}$ in the other answer since that question wanted the opposite transformation.) The variables $\hat{s}_{x,y,z}$ are the components of the unit vector $\hat{s}$.
Now, we have
$$\frac{\vec{h}(t) - \vec{h_c}}{R\sin\beta} = T\begin{bmatrix}\cos{}t\\\sin{}t\\0\end{bmatrix}.$$
Solving for $\vec{h}(t)$,
$$\vec{h}(t) = \vec{h_c} + \left( R \sin\beta \right) T\begin{bmatrix}\cos{}t\\\sin{}t\\0\end{bmatrix}.$$
With everything substituted in and multiplied together, I get
$$\vec{h}(t)
= R \left( \frac{R}{R+h} \begin{bmatrix}\hat{s}_x\\\hat{s}_y\\\hat{s}_z \end{bmatrix} + \sqrt{\frac{1-\left(\frac{R}{R+h}\right)^2}{\hat{s}_x^2 + \hat{s}_y^2}} \begin{bmatrix} -\hat{s}_y \cos{}t - \hat{s}_x \hat{s}_z \sin{}t \\ \hat{s}_x \cos{}t - \hat{s}_y \hat{s}_z \sin{}t \\ \left( \hat{s}_x^2 + \hat{s}_y^2 \right) \sin{}t \end{bmatrix} \right)$$
or
$$\vec{h}(t) = R \left( \frac{R}{\left|\vec{s}\right|} \begin{bmatrix}\hat{s}_x\\\hat{s}_y\\\hat{s}_z \end{bmatrix} + \sqrt{\frac{1-\left(\frac{R}{\left|\vec{s}\right|}\right)^2}{\hat{s}_x^2 + \hat{s}_y^2}} \begin{bmatrix} -\hat{s}_y \cos{}t - \hat{s}_x \hat{s}_z \sin{}t \\ \hat{s}_x \cos{}t - \hat{s}_y \hat{s}_z \sin{}t \\ \left( \hat{s}_x^2 + \hat{s}_y^2 \right) \sin{}t \end{bmatrix} \right)$$
$$ $$
