This is an interesting problem if one insists on a solution constructable using purely geometric operations.
For simplicity in describing the solution, we will choose a coordinate system such that
the two foci $A, B$ are $(1,0)$ and $(0,1)$ respectively. We will also assume the circle
$\mathscr{C}$ we are dealing with is centered at $C = (\alpha,\beta)$ in the first
quadrant and the radius of $\mathscr{C}$ is $r$.
It is known that the family of ellipses and hyperbolas having $A$ and $B$ as foci are given by:
$$\frac{x^2}{1+\lambda} + \frac{y^2}{\lambda} = 1\quad\text{ and is }
\begin{cases}
\text{ an ellipse, }&\text{ if } \lambda > 0\\
\text{ a hyperbola, } &\text{ if } 0 > \lambda > -1
\end{cases}\tag{*0}$$
For a point $p = (x,y)$ on this ellipse/hyperbola, the normal is in the direction $(\frac{x}{1+\lambda}, \frac{y}{\lambda})$. i.e. the normal line passing through $(x,y)$ is given by:
$$\left\{\;\left( x ( 1 + \frac{t}{1+\lambda} ),\,y ( 1 + \frac{t}{\lambda} ) \right) : t \in \mathbb{R}\;\right\}$$
If this ellipse/hyperbola is tangent to any circle centered at $C$ at point $p$, the above normal line with pass through $C = (\alpha,\beta)$ and one can show that $(x,y)$ lies on
a cubic curve $\mathscr{P}$:
$$(\alpha y - \beta x )\left(x(x-\alpha) + y(y-\beta)\right) - (x-\alpha)(y-\beta) = 0\tag{*1}$$
This is the red curve depicted in above picture.
For $r$ not too large, $\mathscr{C}$ and $\mathscr{P}$ in general will intersect at 4 points. Two of them are "tangent" points to two ellipses and the other two are tangent points to two hyperbolas. All these four ellipses/hyperbolas are having $A, B$ as their foci. For these four points, one can simplify $(*1)$ and shows that they lie on a conic
$\mathscr{Q}$:
$$(\alpha y - \beta x)(\alpha x + \beta y - \gamma) - (x-\alpha)(y-\beta) = 0\tag{*2}$$
where $\gamma = \alpha^2 + \beta^2 - r^2$. This is the cyan curve containing
$C, F, G, J, K, L$ in above picture.
Some of these points are relatively easy to construct.
- $C$ - of course we know where $C$ is. In fact, one can show that $\mathscr{Q}$ is tangent
to the line $OC$ at point $C$. Counting multiplicity, this gives "two" points to determine $\mathscr{Q}$.
- If one draws a circle using $OC$ as diameter, it will intersect $\mathscr{C}$ at two points $D$ and $E$. The line $DE$ in turn intersects the horizontal line containing $C$ at $F$ and the vertical line containing $C$ at $G$.
A conic is determined by any five points on it. If for whatever means we get an extra
point on $\mathscr{Q}$, then we can construct $\mathscr{Q}$ and its intersections with $\mathscr{C}$, e.g. the point $L$ in the above picture. These are the tangent points we want.
One can show that the conic $\mathscr{Q}$ intersects the $x$-axis at $J = (\mu,0)$
and $K = (1/\mu,0)$ where $\mu$ is a root of the quadratic equation: $$\alpha x^2 - (\gamma +1) x + \alpha = 0\tag{*3}$$
If one doesn't mind a little bit of algebra, we are done. Otherwise, we can convert this
formula to a geometric construction of $J$ and $K$ in 3 steps:
- Construct an auxiliary hyperbola $xy + 1 = 0$ (the olive green hyperbola in above picture). The foci of this hyperbola are $(\sqrt{2},-\sqrt{2})$ and $(-\sqrt{2},\sqrt{2} )$ and it contains the point $(1,-1)$. All of them can be constructed using compass and ruler start from $O$ and $A$.
- Construct an auxiliary horizontal line $\beta y + 1 = 0$ (the blue horizontal line containing $I$ in above picture). One first constructs a line segment joining $Y = (0,\beta)$ to $B = (-1,0)$. One then constructs another line segment parallel to it
start at $Y' = (0,1)$. This new segment will intersect the $x$-axis at $B' = (-\frac{1}{\beta},0)$. Another quarter circle will bring us to the point $I = (0,-\frac{1}{\beta})$ and the horizontal line at $I$ is the one we need.
- Let $H$ be the intersection of the line $DE$ with the blue line in previous step.
Project it onto $x$-axis as $H_1$. Now construct a ray start from $H_1$ at $-135^{\circ}$ with respect to the $x$-axis. This ray will intersect the "olive green" hyperbola at $H_2$ and $H_3$. If one projects $H_2$ and $H_3$ back to the $x$-axis, we get $K$ and $J$ respectively.
Update
Since Stefan asks, here are some details how to derive the various equations.
If the normal passes through $C = (\alpha,\beta)$, we have a $t$ such that:
$$x (1 + \frac{t}{1+\lambda}) = \alpha \quad\text{ and }\quad y (1+ \frac{t}{\lambda}) = \beta$$
Eliminating $t$ give us:
$$\begin{align}(1+\lambda)(\frac{\alpha}{x}-1) = \lambda (\frac{\beta}{y} - 1 )
\iff & (1+\lambda)(x-\alpha)y - \lambda(y-\beta)\alpha = 0\\
\iff & (x-\alpha)y - \lambda( \alpha y - \beta x ) = 0
\end{align}$$
This gives us
$$
\lambda = \frac{(x - \alpha)y}{\alpha y - \beta x}
\quad\text{ and }\quad
1 + \lambda = \frac{(y - \beta)x}{\alpha y - \beta x}
$$
Substitute this back into $(*0)$, we get $(*1)$:
$$\begin{align}
& (\alpha y - \beta x)\left(\frac{x}{y-\beta} + \frac{y}{x-\alpha}\right) = 1\\
\iff & (\alpha y - \beta x)\left(x(x-\alpha)+y(y-\beta)\right) - (x-\alpha)(y-\beta) = 0
\end{align}$$
If $(x,y)$ also lies on $\mathscr{C}$, the $2^{nd}$ factor in above expression can be simplified as:
$$\begin{align}
x(x-\alpha)+y(y-\beta) = & (x-\alpha)^2 + (y-\beta)^2 + \alpha (x-\alpha)+\beta(y-\beta)\\
= & \alpha x + \beta y - (\alpha^2 + \beta^2 - r^2)\\
= & \alpha x + \beta y - \gamma
\end{align}$$
Substituting this into $(*1)$ gives us the equation of $\mathscr{Q}$ in $(*2)$:
$$(\alpha y - \beta x)(\alpha x + \beta y - \gamma) - (x-\alpha)(y-\beta) = 0$$
Now $(*3)$ is really the equation of intersecting $\mathscr{Q}$ with the $x$-axis.
We just set $y = 0$ in $(*2)$, eliminate the common factor $\beta$ and we are done.
Best Answer
The answer depends on the number of prime factors of $k$ when we consider a factorization over the ring of Eisenstein integers $\mathbb{Z}[\omega]$. Have a look at this answer, too. For instance, there are just trivial solutions if $k$ is a square-free number and a product of primes of the form $3m-1$, that do not split over $\mathbb{Z}[\omega]$. So the number of integer points on such ellipses has a very irregular arithmetic behaviour, but a quite regular behaviour on average, just like in the Gauss circle problem.