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 that`c`

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+1right)^2 h^2+4}+c^2 h+h}{2 sqrt{left(c^2+1right)^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+1right) h}{2 sqrt{left(c^2+1right)^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-1right)^2+4}-h k^2+h}{2 sqrt{h^2 left(k^2-1right)^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-1right)^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[{[email protected],
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}+1right)^2}{2 left(c^4+c^2 left(sqrt{c^4+2 c^2+4 k^2+1}+2right)+sqrt{c^4+2 c^2+4 k^2+1}+4 k^2+1right)}$$

$$frac{2 k^2}{c^4-c^2 left(sqrt{c^4+2 c^2+4 k^2+1}-2right)-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}+1right)^2}{2

left(c^4+left(sqrt{c^4+2 c^2+1}+2right)

c^2+sqrt{c^4+2 c^2+1}+1right)}$$

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+1right)^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^2right)^2+4 h^2}+h^2right)^2}{4

left(frac{1}{4} left(-c^2

h^2+sqrt{left(h^2-c^2 h^2right)^2+4

h^2}+h^2right)^2+h^2right)}$$

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| cright| <1 \

{1}/{2} & left| cright| =1 \

0 & left| cright| >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 $htoinfty$ 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, [email protected][b, {h, Infinity, 0}]]
(* Out: 1 + O[1/h]^2 *)
Assuming[c <= 0, [email protected][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`

.