Let $(x,y)^R=(-y,x)$ represent rotation by $\pi/2$ counterclockwise and
$$
\gamma(t)=(1-t)^3p_0+3t(1-t)^2p_1+3t^2(1-t)p_2+t^3p_3
$$
define a cubic bezier with control points $\{p_0,p_1,p_2,p_3\}$.
Suppose $p_0=(x_0,y_0)$, $p_3=(x_1,y_1)$, and $c=(c_x,c_y)$ are given so that $|p_0-c|=|p_3-c|=r$ ($p_3$ is counterclockwise from $p_0$). Then $p_1=p_0+\alpha(p_0-c)^R$ and $p_2=p_3-\alpha(p_3-c)^R$ where
$$
\alpha=\frac43\tan\left(\frac14\cos^{-1}\left(\frac{(p_0-c)\cdot(p_3-c)}{r^2}\right)\right)
$$
For a quarter of a circle, $\alpha=\frac43(\sqrt2-1)$, and $\gamma$ is no more than $0.00027$ of the radius of the circle off.
Here is a plot of $\gamma$ in red over the quarter circle in black. We really don't see the circle since it is no more than $0.1$ pixels off from $\gamma$ when the radius is $400$ pixels.
$\hspace{3.5cm}$
Computation of $\boldsymbol{\alpha}$
Looking at an arc with an angle of $\theta=\cos^{-1}\left(\frac{(p_0-c)\cdot(p_3-c)}{r^2}\right)$
$\hspace{1.5cm}$
we see that the distance from $c$ to the middle of the arc is
$$
r\cos(\theta/2)+\frac34\alpha r\sin(\theta/2)
$$
we wish to choose $\alpha$ so that this is equal to $r$. Solving for $\alpha$ gives
$$
\begin{align}
\alpha
&=\frac43\frac{1-\cos(\theta/2)}{\sin(\theta/2)}\\
&=\frac43\tan(\theta/4)
\end{align}
$$
A Slight Improvement
Using a circle of radius $1$, the maximum error in radius produced using $\alpha=\frac43\tan(\theta/4)$ is approximately
$$
0.0741\cos^4(\theta/4)\tan^6(\theta/4)
$$
and the error is always positive; that is, the cubic spline never passes inside the circle. Reducing $\alpha$ reduces the midpoint distance by $\frac34\sin(\theta/2)=\frac32\tan(\theta/4)\cos^2(\theta/4)$ times as much, so to distribute the error evenly between the positive and negative, a first guess, assuming that the amplitude of the radius is unchanged, would be to reduce $\alpha$ by $0.0247\cos^2(\theta/4)\tan^5(\theta/4)$.
A bit of investigation shows that, when equalizing the positive and negative swings of the radius, the amplitude increases and that
$$
\alpha=\frac43\tan(\theta/4)-0.03552442\cos^2(\theta/4)\tan^5(\theta/4)
$$
gives pretty even distribution of the error between positive and negative for $\theta\le\pi/2$. The maximum error, both positive and negative, is approximately
$$
0.0533\cos^4(\theta/4)\tan^6(\theta/4)
$$
When $\theta=\pi/2$, this agrees with the article mentioned by bubba in comments.
Note however, that in minimizing the radial error from the circle, the actual variation in radius is increased. Using the simple formula for $\alpha$, which puts the cubic bezier outside the circle, the radius varies by $0.0741\cos^4(\theta/4)\tan^6(\theta/4)$. However, when we minimize the error, the radial variation increases to $0.1066\cos^4(\theta/4)\tan^6(\theta/4)$.
For a circular arc with angular span $\theta$ and unit radius, described as $C(t) = (cos(t), sin(t))$, where $t=[0, \theta]$, a good cubic Bezier curve approximation can be obtained with the following control points:
$P_0 = (1, 0)$,
$P_1 = (1, 0) + L(0,1)$,
$P_2 = (cos\theta, sin\theta) - L (-sin\theta, cos\theta)$,
$P_3 = (cos\theta, sin\theta) $
where $L = \frac43tan(\frac\theta4)$
The cubic Bezier curve approximation obtained this way always interpolates the two end points and the mid-point of the circular arc and the approximation error is always positive, which means the cubic Bezier curve is always "outside" the circular arc. The maximum radial error $(x(t)^2 + y(t)^2 - 1)$ from this approximation can also be computed analytically as
$E_{max} = \frac4{27} ( sin^6(\frac\theta4)/cos^2(\frac\theta4) )$
You can use this maximum error formula to decide how many Bezier segments you want to break up your arc into and then use the control points formula to obtain control points for each Bezier segment.
Best Answer
Wikipedia has it all worked out.
The control points between $(1,0)$ and $(0,1)$ are at $(1,k)$ and $(k,1)$ with $k=\frac43(\sqrt2-1)$.