This is a partial answer. I'm going to derive the equations satisfied by the center of desired ellipses. As of this moment,, solving the equations by anything non-numerical is out of reach. The main result is
The center lies on the intersection of two algebraic curves, one cubic and another quartic.
To avoid conflict with other uses of the variable $a,b$. I will use $\alpha,\beta$ to denote the semi-major/semi-minor axis for the desired ellipse.
main result - center lies on intersection of a cubic and a quartic curve.
First, let us determine the condition for the origin to be such a center. If $\theta$ is the angle between the semi-major axis with the $x$-axis.The equation of the ellipse will be
$$\begin{align}
& {\small \frac{(x\cos\theta + y\sin\theta)^2}{\alpha^2} + \frac{(-x\sin\theta+y\cos\theta)^2}{\beta^2} = 1}\\
\iff & {\small \left(\frac{\cos^2\theta}{\alpha^2} + \frac{\sin^2\theta}{\beta^2}\right)x^2
+ \left(\frac{\sin^2\theta}{\alpha^2} + \frac{\cos^2\theta}{\beta^2}\right)y^2 +
\left(\frac{1}{\alpha^2}-\frac{1}{\beta^2}\right)\sin(2\theta)xy = 1}\\
\iff & {\small \frac12\left(\frac{1}{\alpha^2} + \frac{1}{\beta^2}\right)(x^2 + y^2)
+ \frac12\left(\frac{1}{\alpha^2}-\frac{1}{\beta^2}\right)\left((x^2-y^2)\cos(2\theta) + 2xy\sin(2\theta)\right) = 1 }
\end{align}
$$
Let $\epsilon = \frac{\sqrt{\alpha^2-\beta^2}}{\alpha}$ be the eccentricity of the ellipse and define $\sigma$ and $\lambda$ by
$$\frac{1}{\sigma^2} = \frac12\left(\frac1{\alpha^2} + \frac{1}{\beta^2}\right)
\quad\text{ and }\quad
\lambda = \frac{\alpha^2-\beta^2}{\alpha^2+\beta^2} = \frac{\epsilon^2}{2-\epsilon^2}$$
$\sigma$ and $\lambda$ will be an alternative measure of the size and eccentricity of the ellipse. In particular, for the conic above to be an ellipse, we need $0 \le \epsilon < 1 \iff 0 \le \lambda < 1$. In terms of them, above equation becomes
$$\frac{1}{\sigma^2}(x^2+y^2) - \frac{\lambda}{\sigma^2}\left[(x^2-y^2)\cos(2\theta)+2xy\sin(2\theta)\right] = 1$$
For any $p = (x,y), q = (x', y') \in \mathbb{R}^2$, define
- $\tilde{p} = (-y,x)$, the image of $p$ after a $90^\circ$ counterclockwise rotation.
- $p \cdot q = xx' + y y'\;$ be ordinary 2-d dot product.
- $p \times q \stackrel{def}{=} \tilde{p}\cdot q = x y' - y x'\;$ be the "cross product" for 2-d numbers.
- $U(p) = \begin{bmatrix} x^2 + y^2\\ x^2 - y^2\\ 2xy\end{bmatrix}$
and $V(p) = U(\tilde{p}) = \begin{bmatrix} x^2 + y^2\\ -x^2 + y^2\\ -2xy\end{bmatrix}$
In terms of them, the equation for 3 points $A,B,C$ to lie on the ellipse centered at origin is
$$\Lambda \cdot U(A) = \Lambda \cdot U(B) = \Lambda U(C) = 1\quad\text{ where }\quad
\Lambda = \frac{1}{\sigma^2}\begin{bmatrix}
1\\
-\lambda\cos(2\theta)\\
-\lambda\sin(2\theta)
\end{bmatrix}$$
With a little bit of vector algebra, we can solve above three equations to get
$$\Lambda = \frac{U(A)\times U(B) + U(B)\times U(C) + U(C)\times U(A)}{U(A)\cdot (U(B) \times U(C))}$$
It is not hard to verify the $U(\cdot), V(\cdot)$ satisfy following identities
- $U(P)\cdot U(Q) = V(P)\cdot V(Q) = 2(P\cdot Q)^2$.
- $U(P)\cdot V(Q) = V(P)\cdot U(Q) = 2(P\times Q)^2$.
- $U(P)\times U(Q) = 2(P\times Q)\left[V\left(\frac{P+Q}{2}\right)-V\left(\frac{P-Q}{2}\right)\right]$
Using these identities, the denominator simplifies to
$$\small \begin{align} & U(A)\cdot(U(B)\times U(C))\\
= & 2(B\times C)\,U(A)\times\left[V\left(\frac{B+C}{2}\right) - V\left(\frac{B-C}{2}\right)\right]\\
= & 2(B\times C)\left[ 2\left(A \times \frac{B+C}{2}\right)^2 - 2\left(A \times \frac{B-C}{2}\right)^2\right]\\
= & -4(B \times C)(C\times A)(A\times B)
\end{align}$$
Let $\Delta$ be the area of $\triangle ABC$ and $u,v,w$ be the barycentric coordinates of the origin with respect to $\triangle ABC$. We have
$$2\Delta u = B \times C,\quad 2\Delta v = C \times A\quad\text{ and }\quad 2\Delta w = A \times B$$
We find the denominator equals to $-32\Delta^3 uvw$.
Doing similar thing to the numerator and let $D, E, F$ be the midpoint of $BC, CA$ and $AB$, the numerator becomes
$$\small 4\Delta \left[ u (V(D) - V(E-F)) + v(V(E) - V(F-D)) + w(V(F)-V(D-E)) \right]$$
This is a little bit clumsy to write down. Let us use the notation $\sum_{cyc}$ to indicate a cyclic sum of over parameter set
$$u \to v \to w,\quad A \to B \to C\quad\text{ and }\quad D \to E \to F$$
In this new notation, the equation for $\Lambda$ becomes
$$\Lambda = - \frac{1}{8\Delta^2 uvw} \sum_{cyc} u(V(D) - V(E-F))$$
Switch to other coordinate system where the center, let's call it $Z$, is no longer the origin, the equation for the center becomes
$$\Lambda = -\frac{1}{8\Delta^2 uvw}\sum_{cyc}u(V(Z-D) - V(E-F))\tag{*1}$$
To proceed further, we will do two things.
We will switch to a coordinate system where $N$, the nine-point center of $\triangle ABC$ is the origin. In this coordinate system, $D,E,F$ will be lying on a circle centered at origin $N$ with radius $\frac{R}{2}$ ($R$ is the circumradius of $\triangle ABC$).
We will identify the Euclidean plane with the complex plane, we will use the lower case letter to denote the complex number corresponds to a point. e.g. $Z$ becomes $z$ and $D, E, F$ become $d, e, f$.
After this, we can repress $(*1)$ as two equations
$$\begin{align}
\Lambda_1 = \frac{1}{\sigma^2}
&= -\frac{1}{8\Delta^2 uvw} \sum_{cyc} u(|z-d|^2 - |e - f|^2)\tag{*2a}\\
\Lambda_2 + i\Lambda_3 = -\frac{\lambda}{\sigma^2} e^{2i\theta}
&= +\frac{1}{8\Delta^2 uvw} \sum_{cyc} u((z-d)^2 - (e - f)^2)\tag{*2b}
\end{align}
$$
Let $o = d + e + f$ (the point corresponding to it is the circumcenter $O$ of $\triangle ABC$) and $u',v',w'$ the barycentric coordinate of $Z$ with respect to $\triangle DEF$. We have
$$u = \frac{1-u'}{2},\cdots
\implies \sum_{cyc} ud = \frac12\sum_{cyc}(1-u')d = \frac12(o-z)$$
Notice
$$\begin{align}
|z-d|^2 - |e-f|^2 &= |z-d|^2 + |e+f|^2 - 2|e|^2 - 2|f|^2\\
&= |z-d|^2+|o-d|^2 - 4|d|^2\\
&= |z|^2 -d(\overline{z+o}) - \bar{d}(z+o) + |o|^2 - 2|d|^2
\end{align}$$
Multiply by $u$ and take cyclic sum, we get
$$\begin{align} \sum_{cyc}u(|z-d|^2 - |e-f|^2)
= & |z|^2 - \frac{(o - z)(\overline{z+o})}{2} - \frac{(\overline{o-z})(z+o)}{2} + |o|^2 - 2|d|^2\\
= & 2|z|^2 - \frac{R^2}{2}\end{align}$$
Equation $(*2a)$ becomes
$$\bbox[border:1px solid blue;padding: 1em;]{z\bar{z} - \frac{R^2}{4} + \frac{4\Delta^2 uvw}{\sigma^2} = 0} \tag{*3a}$$
Since $u,v,w$ are linear in $(x,y) = (\Re z,\Im z)$. This describe a cubic curve
in the Euclidean plane and it is the locus of the center when the size parameter $\sigma$ is held fixed.
Let $\Omega = \frac{def}{|d|^2} = \frac{4def}{R^2}$. By a similar procedure, we have
$$\begin{align}(z-d)^2 - (e-f)^2
&= (z-d+e+f)(z-d-e-f) = (z-o+2e)(z-o+2f)\\
&= (z-o)^2+2(e+f)(z-o) + 4ef\\
&= z^2 - o^2 - 2d(z-o) + 4\Omega \bar{d}
\end{align}
$$
Multiply by $u$ and take cyclic sum, we get
$$\begin{align}\sum_{cyc} u((z-d)^2 - (e-f)^2)
&= z^2 - o^2 - (o-z)(z-o) + 2\Omega(\bar{o} - \bar{z})\\
&= 2\left(z(z-o) - \Omega(\bar{z} - \bar{o})\right)\end{align}$$
Equations $(*2b)$ becomes
$$z(z-o) - \Omega(\bar{z} - \bar{o})
= -\lambda e^{2i\theta}\frac{4\Delta^2 uvw}{\sigma^2}$$
Compare this with equation $(*3a)$, we obtain
$$z(z-o) - \Omega(\bar{z} - \bar{o}) = \lambda e^{2i\theta} \left(z\bar{z} - \frac{R^2}{4}\right)$$
Taking absolute value and square, we obtain a quartic curve
$$\bbox[border:1px solid blue;padding: 1em]{|z(z-o) - \Omega(\bar{z}-\bar{o})|^2 = \lambda^2 \left|z\bar{z} - \frac{R^2}{4}\right|^2}\tag{*3b}
$$
This is the locus of the center when the eccentricity $\epsilon$ is help fixed.
observation I - locus for $\sigma$.
All locus of $\sigma$ passes through $D,E,F$ and the three foots of $\triangle DEF$.
When orthocenter $H$ of $\triangle ABC$ belongs to interior of $\triangle DEF$. $\sigma$ takes a local minimum at $H$. This should corresponds to a special circumellipse associated with $\triangle ABC$. I'm unable to figure out what that is. Anyone has any idea?
observation II - locus for $\lambda$.
When $\epsilon = 0 \implies \lambda = 0$, the locus reduces to a single point $O$. The circumcenter of $\triangle ABC$.
When $\epsilon = 1 \implies \lambda = 1$, the quartic terms in Eq. ($*3b$) cancel out. The corresponding curve reduces a union of $3$ lines through the midpoints. i.e. the lines $DE$, $EF$ and $FD$.
These $3$ lines split the plane into $7$ regions. When one varies $\epsilon$ (and hence $\lambda$) from $0$ to $1$. The locus sweep across the interior of
$4$ out of the $7$ regions ( $\triangle DEF$ or the $3$ cone with apex at $D, E, F$). This means in order for a point to be center for an ellipse, it need to be
either $D, E, F$ or belongs to the interior of above $4$ regions.
Infinitely many loci of $\lambda$ passing through $D, E, F$. Combine with the observation $I_1$, $D, E, F$ are centers for infinitely many circumellipses of $A,B,C$.
To be continued???
Best Answer
In a purely geometric fashion, given two tangents and two tangency points on them,
the wanted ellipse can be constructed as the Steiner inellipse of $ABC$, where $A$ is given by the intersection of the given tangents and $B,C$ are given by the symmetric of $A$ with respect to the tangency points. In particular the center of the wanted ellipse lies at the centroid of $ABC$.