There are excellent answers of kguler and belisarius. Not to compete with them, but to offer a different view, you might operate within a different paradigm. Here it is:
Step 1: Let us first find the solution of your equation for J
:
sl = Solve[16 x^3 - 4 x^2 + J^2/64 == 0, J]
(* {{J -> -16 Sqrt[x^2 - 4 x^3]}, {J -> 16 Sqrt[x^2 - 4 x^3]}} *)
There are two, and let us operate on the positive one. The negative may be treated analogously.
Step 2: Let us create a table entitled "lst" out of this second solution with the structure {J,x}
. It is easily seen that x=x(J)
is a two-valued function and I only select the smaller branch:
lst = Select[
Table[{sl[[2, 1, 2]], x}, {x, 0, 0.25, 0.01}], #[[2]] <= 0.17 &]
(* {{0., 0.}, {0.156767, 0.01}, {0.306933, 0.02}, {0.45028,
0.03}, {0.58657, 0.04}, {0.715542, 0.05}, {0.836909,
0.06}, {0.950352, 0.07}, {1.05552, 0.08}, {1.152, 0.09}, {1.23935,
0.1}, {1.31706, 0.11}, {1.38453, 0.12}, {1.44107, 0.13}, {1.48585,
0.14}, {1.51789, 0.15}, {1.536, 0.16}, {1.53866, 0.17}} *)
The upper one can be selected analogously.
Step 3: Let us now approximate (that is, fit) the obtained table with a polynomial:
model1 = a + b*J + c*J^2 + d*J^3;
ff1 = FindFit[lst, model1, {a, b, c, d}, J]
model2 = a + b*J + c*J^2 + d*J^3 + e*J^4;
ff2 = FindFit[lst, model2, {a, b, c, d, e}, J]
(* {a -> -0.00294399, b -> 0.106524, c -> -0.0890299, d -> 0.0569192}
{a -> 0.00137571, b -> 0.0172211, c -> 0.186767, d -> -0.219333,
e -> 0.0872711} *)
The first solution is for the first model (that is, the expansion up to the cube), while the second is for the second one (the expansion up to J^4).
Let us draw them to see the difference:
Show[{
ListPlot[lst, AxesLabel -> {"J", "x"}],
Plot[{model1 /. ff1, model2 /. ff2}, {J, 0, 1.5},
PlotStyle -> {Red, Green}]
}]
It returns the following plot:

The red here shows the fitting with the first model, while green - with the second one.
Since this fitting yields a polynomial this is though not equal to, but not really too far away from what kguler and belisarius proposed.
The real difference might come, if the expansion is not a final aim, but a step in some more general problem. Then one might think about a simple and exact enough approximation. Just to give a simple example, have a look at the following two models with a smaller amount of parameters:
model3 = (a*J)/(1 + d*J^2);
ff3 = FindFit[lst, model3, {a, d}, J]
model4 = (a*J + b*J^2)/(1 + d*J^2);
ff4 = FindFit[lst, model4, {a, b, d}, J]
(* {a -> 0.0596691, d -> -0.176778}
{a -> 0.0807892, b -> -0.0277843, d -> -0.269206} *)
The visualization is below:

Have fun!
Series[ToRadicals@Root[j^2 - 256 #1^2 + 1024 #1^3 &, 2], {j, 0, 4}]
? $\endgroup$Series[ToRadicals@ First[x /. First@Solve[16 x^3 - 4 x^2 + J^2/64 == 0, x, Reals][[2]]], {J, 0, 4}]
$\endgroup$