(another comment, not an answer.)

Experimentally, with the following code
```python
from sympy import *
var('x')
var('n', integer = True)

f = ( (2+(x+3) * sqrt(-x))/(2*(x+4)) *
      sqrt(-x)**n + ((2-(x+3) * sqrt(-x)) / (2*(x+4))) * (-sqrt(-x))**n +
      + (x+2+sqrt(x*(x+4)))/(2*(x+4)) *
        ( (x+sqrt(x*(x+4)))/2 )**n +
        (x+2-sqrt(x*(x+4)))/(2*(x+4)) *
        ( (x-sqrt(x*(x+4)))/2 )**n
        )
pprint(f)
for i in range(10):
   eq = simplify(f.subs(n,i))
   print ('========== n = ', i)
   pprint(eq)
   print ('have solutions')
   sols = solve(eq)
   pprint (sols)
   pprint ('approx. = ')
   pprint ([s.evalf() for s in sols])
```
we see the following:
$$
\begin{array}{l}
f_{0} = 1\\
f_{1} = 0\\
f_{2} = x^{2}\\
f_{3} = x^{2} \left(x + 2\right)\\
f_{4} = x^{2} \left(x^{2} + 2 x + 1\right)\\
f_{5} = x^{3} \left(x^{2} + 3 x + 1\right)\\
f_{6} = x^{4} \left(x^{2} + 4 x + 4\right)\\
f_{7} = x^{4} \left(x^{3} + 5 x^{2} + 7 x + 3\right)\\
f_{8} = x^{4} \left(x^{4} + 6 x^{3} + 11 x^{2} + 6 x + 1\right)\\
f_{9} = x^{5} \left(x^{4} + 7 x^{3} + 16 x^{2} + 13 x + 2\right)\\
\end{array}
$$

[![enter image description here][1]][1]

And it seems that already $f_9$ have some imaginary roots, albeit very small.
Do these imaginary roots really converges to 0 when $n\to \infty$ or they stay on the same magnitude ?


  [1]: https://i.sstatic.net/z91X1.png