An apparently "simple" limit?
There's a bit more to the story. Mathematica treats variables as complex by default, and I for one have had trouble figuring out how Limit figures out how to treat variables such as c in this case.
Some analysis
First, let's examine a0 (= a in OP) with the assumption thatc is real:
a0 = (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2 /
(4 (h^2 + 1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2));
aR = FullSimplify[a0, h > 0 && c \[Element] Reals]
(* -> (h + c^2 h + Sqrt[4 + (1 + c^2)^2 h^2])/(2 Sqrt[ 4 + (1 + c^2)^2 h^2]) *)
$$\frac{\sqrt{\left(c^2+1\right)^2 h^2+4}+c^2 h+h}{2 \sqrt{\left(c^2+1\right)^2 h^2+4}}$$
It simplifies a little if we split the fraction up:
1/2 + Factor[aR - 1/2]
(* -> 1/2 + ((1 + c^2) h)/(2 Sqrt[4 + (1 + c^2)^2 h^2]) *)
$$\frac{\left(c^2+1\right) h}{2 \sqrt{\left(c^2+1\right)^2 h^2+4}}+\frac{1}{2}$$
And let's look at it if c is an imaginary number k I:
aI = FullSimplify[a0 /. c -> k I, h > 0 && k \[Element] Reals]
(* -> (h - h k^2 + Sqrt[4 + h^2 (-1 + k^2)^2])/(2 Sqrt[4 + h^2 (-1 + k^2)^2]) *)
$$\frac{\sqrt{h^2 \left(k^2-1\right)^2+4}-h k^2+h}{2 \sqrt{h^2 \left(k^2-1\right)^2+4}}$$
Again it simplifies if we split the fraction up:
1/2 + Factor[aI - 1/2]
(* -> 1/2 - (h (-1 + k) (1 + k))/(2 Sqrt[4 + h^2 (-1 + k^2)^2]) *)
$$\frac{1}{2}-\frac{h (k-1) (k+1)}{2 \sqrt{h^2 \left(k^2-1\right)^2+4}}$$
We see if c is +/-I (if k is +/-1), the fraction is 0, so the limit is 1/2. If k > 1 or k < -1, then the limit of the fraction is -1/2, so the limit of a0 is 0. And if -1 < k < 1, then the limit of the fraction is +1/2, so the limit of a0 is 1. If c is real, then up above one can see the limit of the fraction is 1/2, so the limit of a0 is 1.
Here are some examples:
Table[Print["c = ", c, " --> ", HoldForm[Limit[a0, h -> Infinity]],
" = ", Limit[a0, h -> Infinity]], {c, {1, I/2, I, 2 I}}];
c = 1 --> Limit[a0, h -> Infinity] = 1 c = I/2 --> Limit[a0, h -> Infinity] = 1 c = I --> Limit[a0, h -> Infinity] = 1/2 c = 2 I --> Limit[a0, h -> Infinity] = 0
How to find the limit sought
One should probably use assumptions:
a1 = Simplify[a];
a2 = FullSimplify[a];
Block[{$Assumptions = c \[Element] Reals}, Limit[a0, h -> Infinity]]
Block[{$Assumptions = c \[Element] Reals}, Limit[a1, h -> Infinity]]
Block[{$Assumptions = c \[Element] Reals}, Limit[a2, h -> Infinity]]
(* -> 1, 1, 1 *)
Then the answers are correct.
EDIT - Addendum
Prompted to investigate further by various comments, I have something to add, including what the actual limit for a generic complex $c$ is (which one might expect Mathematica could return as the correct answer). This goes beyond the current version of the OP's question, which stipulates that $c$ be real, but I hope it sheds some light on the issues Limit has dealing with this function.
What should the limit be?
Working it out by hand, I get the limit for a complex $c$ to be
$$\begin{cases} 1 & {\rm Im}(c)^2-{\rm Re}(c)^2<1 \\ {1}/{2} & {\rm Im}(c)^2-{\rm Re}(c)^2=1 \\ 0 & {\rm Im}(c)^2-{\rm Re}(c)^2>1 \end{cases}$$
Here is a plot of the real part of a0 vs. c when h = 100. The imaginary parts are nearly zero, and the real part is nearly equal to the limit. You can see the discontinuity forming along the red cliff. The limits at six test values are shown by the green spheres. We will use these test values for c again later.
testvalues = {1, I/2, I, 2 I, 1 + 2 I, 2 + 2 I};
Show[ParametricPlot3D[
Evaluate[{Re[c], Im[c], Re[a0]} /. {c -> r E^(I q), h -> 100.}],
{r, 0, 3}, {q, 0, 2 π}, PlotPoints -> {50, 150}, MaxRecursion -> 5,
Exclusions -> None, MeshFunctions -> {#2^2 - #1^2 &, #4 &, #5 &},
MeshStyle -> {Red, GrayLevel[0.5], GrayLevel[0.5]},
Mesh -> {{-1, 0, 1, 2}; {0.97, 1.03}, 11, {0, π/2, π, (3 π)/2}},
MeshShading -> {{{Yellow, Red, Lighter[Purple, 0.5]}}},
BoundaryStyle -> GrayLevel[0.5], Lighting -> "Neutral", PlotRange -> Full,
AxesLabel -> {Re[c], Im[c], ""}],
Graphics3D[{Darker@Green,
Sphere[{Re[#], Im[#], Limit[a0 /. {c -> #}, h -> Infinity]}, 0.07] & /@
testvalues}]
]
There are spikes (poles) where the denominator of a0 is zero:
Solve[((4 (h^2 + 1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2)) /. h -> 100) == 0, c] // N
(* -> {{c -> -0.0099995 + 1.00005 I}, {c -> 0.0099995 - 1.00005 I}, {c -> -0.0099995 - 1.00005 I}, {c -> 0.0099995 + 1.00005 I}} *)
How could Mathematica get different answers?
First, it is not a problem with Simplify (a1 above) or FullSimplify (a2 above). All give the same correct limits with definite values for c:
Limit[a0 /. {c -> #}, h -> Infinity] & /@ testvalues
Limit[a1 /. {c -> #}, h -> Infinity] & /@ testvalues
Limit[a2 /. {c -> #}, h -> Infinity] & /@ testvalues
(* -> {1, 1, 1/2, 0, 0, 1}
{1, 1, 1/2, 0, 0, 1}
{1, 1, 1/2, 0, 0, 1} *)
I came up with ways to arrive algebraicallly at the different answers for the limits of a0 and a2 for a generic c. I doubt this is how Mathematica arrives at its answers, but it shows how the mistake might arise.
To find a limit at infinity, we can replace h by 1/k and let k approach 0. In calculus we learn to simplify the expression until we can plug in k = 0.
b0 = Simplify[a0 /. h -> 1/k, k > 0]
b2 = Simplify[a2 /. h -> 1/k, k > 0]
(* -> (1 + c^2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2])^2/(2 (1 + c^4 + 4 k^2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2] + c^2 (2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2]))) *)
(* -> (2 k^2)/(1 + c^4 + 4 k^2 - Sqrt[1 + 2 c^2 + c^4 + 4 k^2] - c^2 (-2 + Sqrt[1 + 2 c^2 + c^4 + 4 k^2])) *)
b0 // TeXForm // Print
b2 // TeXForm // Print
$$\frac{\left(c^2+\sqrt{c^4+2 c^2+4 k^2+1}+1\right)^2}{2 \left(c^4+c^2 \left(\sqrt{c^4+2 c^2+4 k^2+1}+2\right)+\sqrt{c^4+2 c^2+4 k^2+1}+4 k^2+1\right)}$$
$$\frac{2 k^2}{c^4-c^2 \left(\sqrt{c^4+2 c^2+4 k^2+1}-2\right)-\sqrt{c^4+2 c^2+4 k^2+1}+4 k^2+1}$$
So far, so good:
b0 == b2 // Simplify
(* -> True *)
For b0, setting k to zero yields 1 the same as Limit[a0, h -> Infinity]
b0 /. k -> 0 // TeXForm
$$\frac{\left(c^2+\sqrt{c^4+2 c^2+1}+1\right)^2}{2 \left(c^4+\left(\sqrt{c^4+2 c^2+1}+2\right) c^2+\sqrt{c^4+2 c^2+1}+1\right)}$$
which we can see by multiplying out the numerator and simplifying:
b0 /. k -> 0 // Expand // Together
(* -> 1 *)
Surprisingly, Simplify does something else and returns
b0 /. k -> 0 // Simplify // TeXForm
$$\frac{c^2+\sqrt{\left(c^2+1\right)^2}+1}{2 c^2+2}$$
whose value depends on c:
% /. {c -> #} & /@ testvalues
(* -> {1, 1, Indeterminate, 0, 0, 1} *)
Turning to b2, setting k to zero yields 0, the same as Limit[a2, h -> Infinity]
b2 /. k -> 0
(* -> 0 *)
As with a0 and a2, the limits of b0 and b2 at definite values of c are correct:
Limit[b0 /. {c -> #}, k -> 0] & /@ testvalues
Limit[b2 /. {c -> #}, k -> 0] & /@ testvalues
(* -> {1, 1, 1/2, 0, 0, 1}
{1, 1, 1/2, 0, 0, 1} *)
A related example
Change 1+c^2 to 1-c^2:
aIm = a0 /. c -> I c
(* -> (h^2 - c^2 h^2 + Sqrt[4 h^2 + (h^2 - c^2 h^2)^2])^2/(4 (h^2 + 1/4 (h^2 - c^2 h^2 + Sqrt[4 h^2 + (h^2 - c^2 h^2)^2])^2)) *)
$$\frac{\left(-c^2 h^2+\sqrt{\left(h^2-c^2 h^2\right)^2+4 h^2}+h^2\right)^2}{4 \left(\frac{1}{4} \left(-c^2 h^2+\sqrt{\left(h^2-c^2 h^2\right)^2+4 h^2}+h^2\right)^2+h^2\right)}$$
Now I'm primarily interest in real c and the limit at h -> Infinity.
It is similar to previous one above:
$$\begin{cases} 1 & \left| c\right| <1 \\ {1}/{2} & \left| c\right| =1 \\ 0 & \left| c\right| >1 \end{cases}$$
Mathematica gives one of the answers by default
Limit[aIm, h -> Infinity]
(* -> 1 *)
At definite values of c, we get the right answers
Limit[aIm /. {c -> #}, h -> Infinity] & /@ {0, 1, 2}
(* -> {1, 1/2, 0} *)
With various assumptions, one cannot get an answer that depends on c.
Limit[aIm, h -> Infinity, Assumptions -> c \[Element] Reals]
(* -> 1 *)
The correct answer appears in some cases:
Limit[aIm, h -> Infinity, Assumptions -> c > 1 || c < -1]
Limit[aIm, h -> Infinity, Assumptions -> c == 1]
Limit[aIm, h -> Infinity, Assumptions -> c == -1]
Limit[aIm, h -> Infinity, Assumptions -> -1 < c < 1]
(* -> 0, 1/2, 1/2, 1 *)
These are wrong; the limit is 1/2 in both cases:
Limit[aIm, h -> Infinity, Assumptions -> c \[Element] Reals && Abs[c] == 1]
Limit[aIm, h -> Infinity, Assumptions -> c == -1 || c == 1]
(* -> 0, 0 *)
Conclusion
Whatever definite value I set c to, I always got the correct limit, but Limit does not handle c as a symbol correctly. With some careless algebraic operations, I was able to derive the answers Limit gave for original function and the FullSimplified version. Even with Assumptions that yield a unique limit, Limit does not always give the right answer. I would be surprised if it were mathematically impossible to develop an algorithm for finding a large class of limits including the OP's, given what can be done for integrals and by standard calculus techniques. I've never found a lot of use for Limit, but perhaps it is because of its limitations.
It is a bug in Series. Note that
a := (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/( 4 (h^2 + 1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2))
b = FullSimplify[a]
Series[a,{h,Infinity,0}]
(* Out: 1 + O[1/h]^2 *)
Series[b,{h,Infinity,0}]
(* Out: O[1/h]^2 *)
The fact is that for $h\to\infty$ there are two terms cancelling each other in the denominator of $b$, one of which is a square root.
In fact further assumptions make it work
Assuming[c >= 0, Simplify@Series[b, {h, Infinity, 0}]]
(* Out: 1 + O[1/h]^2 *)
Assuming[c <= 0, Simplify@Series[b, {h, Infinity, 0}]]
(* Out: 1 + O[1/h]^2 *)
This is a known limitation in Series and Limit. Series does not handle roots in a flawless manner. For example, here is an expansion at the branch point (zero) that is only "half" correct.
In[4]:= Series[Sqrt[x^2], {x,0,2}]
3
Out[4]= x + O[x]
This is fine for re(x)>0, but not so good for re(x)<0 in a neighborhood of the origin.
Limit has the ability to handle this situation if provided sufficiently strong assumptions. For this example one could do as follows.
a = (h^2 + c^2 h^2 +
Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/(4 *(h^2 +
1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2));
a2 = FullSimplify[a];
Limit[a2, h -> Infinity, Assumptions -> Element[c, Reals]]
(* Out[80]= 1 *)
This next one is actually not so good.
Limit[a, h -> Infinity]
(* Out[63]= 1 *)
That is to say, while the original result was what the poster wanted, it is only by accident of haw some symbolic radicals were handled. Here is why it is not fully correct: for "many" values of c the limit is not 1.
Limit[{a, a2} /. c -> 2 I, h -> Infinity]
(* Out[81]= {0, 0} *)
Indeed, there is an interplay between real and imaginary parts that determines the limiting behavior.
Limit[{a, a2} /. c -> 1 - 2 I, h -> Infinity]
(* Out[85]= {0, 0} *)
Limit[{a, a2} /. c -> 3 - 2 I, h -> Infinity]
(* Out[86]= {1, 1} *)
As for bug reports, this one is in the data base and has been for quite some time. It is unlikely to change unless and until there is some redesign both of Series and Limit.