Question
Consider two single qubit states $\left\{|\alpha_0\rangle,|\alpha_1\rangle\right\}$ which are not orthogonal or parallel, i.e. $\left|\langle\alpha_0|\alpha_1\rangle\right|\ne0,1$. Additionally, consider the unitary operation: $U|\alpha_i\rangle|0\rangle=|c_{ii}\rangle\forall i\in\left\{0,1\right\}$; where we want $|c_{ii}\rangle$ to approximate $|\alpha_i\rangle|\alpha_i\rangle$, i.e. approximate cloning - note that $|c_{ii}\rangle\ne|\alpha_i\rangle|\alpha_i\rangle\forall i$ due to the no-cloning theorem. We can define the fidelity $F_i\equiv\left|\langle\alpha_i|\langle\alpha_i|c_{ii}\rangle\right|^2$ as a measure of the quality of the approximation such that the best approximation is given by maximising $F_i$ with respect to $|c_{ii}\rangle$ subject to the following constraints:
- $|c_{ii}\rangle\in\operatorname{span}\left\{|\alpha_j\rangle|\alpha_j\rangle\right\}$
- $|\alpha_i\rangle|0\rangle\to|c_{ii}\rangle\forall i$ is unitary $\iff\langle\alpha_i|\alpha_j\rangle=\langle c_{ii}|c_{jj}\rangle\forall i,j$
- $F=F_0=F_1$
Contraint 1 allows us to express:
$$|c_{ii}\rangle\equiv\sum_ja_{ji}|\alpha_j\rangle|\alpha_j\rangle$$
Now if we let $A_{ij}=\langle\alpha_i|\alpha_j\rangle$ then constraint 2 can be expressed as:
$$\boldsymbol A=\boldsymbol a^\dagger\boldsymbol A^{\circ 2}\boldsymbol a$$
where $\boldsymbol A^{\circ 2}$ is the Hadamard power i.e. $\left[\boldsymbol A^{\circ 2}\right]_{ij}\equiv A_{ij}^2$. In this notation we can also express $F_i=\left|\left[\boldsymbol A\boldsymbol a\right]_{ii}\right|^2$.
We have already applied constraint 1 in the notation used leaving both constraints 2 and 3 to be applied.
As far as I see the problem of finding the values of $|c_{ii}\rangle$ that maximise $F$ from here can then be solved in one of three ways - each of which I have tried with no luck. I would be very grateful if someone was able to outline a solution for complex $\boldsymbol a\in\mathbb C^{2\times2}$. However, if no elegant solution for $\boldsymbol a\in\mathbb C^{2\times2}$ then a solution for $\boldsymbol a\in\mathbb R^{2\times2}$ where $|\alpha_i\rangle$ can be expressed in terms of $|0\rangle$ and $|1\rangle$ with only real coefficients would suffice to be accepted as the correct answer. I know there should be an elegant solution for $\boldsymbol a\in\mathbb R^{2\times2}$ as this was what was asked on a problem sheet I have been using for some self-learning, but I am interested in the more general solution $\boldsymbol a\in\mathbb C^{2\times2}$ if possible.
Finally, while I am an avid StackExchange user, this is my first post on Quantum Computing StackExchange and so I apologise if I have broken any site-specific rules and please do let me know. Additionally, if the moderators feel this is more appropriate for either Physics or Maths SE do migrate.
Attempts
Below I will outline three of my attempted methods in case these are of any help.
$\boldsymbol a\in\mathbb C^{2\times2}$ has 8 degrees of freedom, but as $F_i$ is independent of an overall phase factor the degrees of freedom reduce to 7 and we can express:
$$\boldsymbol a\equiv a\begin{bmatrix}1&be^{i\beta}\\ce^{i\gamma}&de^{i\delta}\end{bmatrix}$$
Using the method of Lagrange multipliers with a matrix constraint we must maximise:
$$L_i\equiv\left|\left[\boldsymbol A\boldsymbol a\right]_{ii}\right|^2+\operatorname{Tr}\left(\boldsymbol \mu\left[\boldsymbol A-\boldsymbol a^\dagger\boldsymbol A^{\circ 2}\boldsymbol a\right]\right)+\lambda\left(\left|\left[\boldsymbol A\boldsymbol a\right]_{ii}\right|^2-\left|\left[\boldsymbol A\boldsymbol a\right]_{\overline{ii}}\right|^2\right)\quad\forall i$$
where $\boldsymbol\mu\in\mathbb C^{2\times2}$ and $\lambda\in\mathbb C$ are the lagrange multipliers; and $\overline i\equiv\begin{cases}0,&i=1\\1,&i=0\end{cases}$. Thus,
$$\text{d}L_i=\left(1+\lambda\right)\left(\vec\pi_i^\dagger\boldsymbol a^\dagger\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol{A\text{d}a}\vec\pi_i+\vec\pi_i^\dagger\text{d}\boldsymbol a^\dagger\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol{Aa}\vec\pi_i\right)-\operatorname{Tr}\left(\boldsymbol{\mu a}^\dagger\boldsymbol A^{\circ2}\text{d}\boldsymbol a\right)-\operatorname{Tr}\left(\text{d}\boldsymbol a^\dagger\boldsymbol A^{\circ2}\boldsymbol{a\mu}\right)-\lambda\left(\vec\pi_\overline{i}^\dagger\boldsymbol a^\dagger\boldsymbol A\vec\pi_\overline{i}\vec\pi_\overline{i}^\dagger\boldsymbol{A\text{d}a}\vec\pi_\overline{i}+\vec\pi_\overline{i}^\dagger\text{d}\boldsymbol a^\dagger\boldsymbol A\vec\pi_\overline{i}\vec\pi_\overline{i}^\dagger\boldsymbol{Aa}\vec\pi_\overline{i}\right)$$
As $\left\{\text{d}\boldsymbol a\vec\pi_i,\vec\pi_i^\dagger\text{d}\boldsymbol a^\dagger\middle|\forall i\right\}$ are independent variables and $\operatorname{Tr}\left(\boldsymbol M\right)\equiv\sum_i\vec\pi_i^\dagger\boldsymbol M\vec\pi_i$ then for $\text{d}L_i=0$ the simultaneous equations are:
$$\begin{cases}\left(\delta_{ij}\left(1+2\lambda\right)-\lambda\right)\vec\pi_i^\dagger\boldsymbol a^\dagger\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol A-\pi_i^\dagger\boldsymbol{\mu a}^\dagger\boldsymbol A^{\circ2}=0&\forall i,j\\\left(\delta_{ij}\left(1+2\lambda\right)-\lambda\right)\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol{Aa}\vec\pi_i-\boldsymbol A^{\circ2}\boldsymbol{a\mu}\vec\pi_i=0&\forall i,j\\\vec\pi_i^\dagger\boldsymbol a^\dagger\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol{A\text{d}a}\vec\pi_i=F&\forall i\\\boldsymbol A=\boldsymbol a^\dagger\boldsymbol A^{\circ 2}\boldsymbol a\end{cases}$$
This gives $\lambda=-\frac{1}{2}$ and so:
$$\begin{cases}\frac{1}{2}\vec\pi_i^\dagger\boldsymbol a^\dagger\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol A-\pi_i^\dagger\boldsymbol{\mu a}^\dagger\boldsymbol A^{\circ2}=0&\forall i\\\frac{1}{2}\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol{Aa}\vec\pi_i-\boldsymbol A^{\circ2}\boldsymbol{a\mu}\vec\pi_i=0&\forall i\\\vec\pi_i^\dagger\boldsymbol a^\dagger\boldsymbol A\vec\pi_i\vec\pi_i^\dagger\boldsymbol{A\text{d}a}\vec\pi_i=F&\forall i\\\boldsymbol A=\boldsymbol a^\dagger\boldsymbol A^{\circ 2}\boldsymbol a\end{cases}$$
Then we can easily see that $F=\vec\pi_i^\dagger\boldsymbol{A\mu}\vec\pi_i=\vec\pi_i^\dagger\boldsymbol{\mu A}\vec\pi_i$ for all $i$. However, I do not know how I would proceed further than this.
Another approach would be to use constraint 3 to reduce the degrees of freedom to 6 by showing:
$$de^{i\delta}=-\alpha^*be^{i\beta}+\left|1+\alpha ce^{i\gamma}\right|e^{i\epsilon}$$
where $\alpha\equiv\langle\alpha_0|\alpha_1\rangle$. Now our remaining degrees of freedom are $\left\{a,b,c,\beta,\gamma,\epsilon\right\}$ and we need only maximise:
$$L_i'\equiv\left|\left[\boldsymbol A\boldsymbol a\right]_{ii}\right|^2+\operatorname{Tr}\left(\boldsymbol \mu\left[\boldsymbol A-\boldsymbol a^\dagger\boldsymbol A^{\circ 2}\boldsymbol a\right]\right)\quad\forall i$$
I won't write out the resulting simultaneous equations again as this time they are more involved as the reduction in degrees of freedom means $\left\{\text{d}\boldsymbol a\vec\pi_i,\vec\pi_i^\dagger\text{d}\boldsymbol a^\dagger\middle|\forall i\right\}$ are not all independent degrees of freedom.
Or finally, we could apply both constraints 2 and 3 before maximising; however, I have not managed to navigate my way through the swamp of simultaneous equations for this but in this case, we now only need to maximise $F$ with respect to the remaining degrees of freedom.
Now I have in all three cases managed to find the conditions under which these are maximised; however, I then have been struggling to solve the resulting simultaneous equations.