For which values of $a\in\mathbb{Q}$ does integer solutions to $x^2+x+1=a(y^2+1)$ exist? I am unable to determine for which values of $a\in\mathbb{Q}$ does integer solutions to
$$x^2+x+1=a(y^2+1)$$
in the form $(x,y)$ exist.
My initial idea was to set $a=\frac{c}{d}$ to get $dx^2+dx+d=c(y^2+1)$ for $c,d\in\mathbb{Z}$ but I am unsure how to convert this into Pell's equation nor do I know how to apply initial solutions. Any suggestions?
 A: I'll just speak to the question of how to get from $dx^2+dx+d=c(y^2+1)$ to (something like) a Pell equation: We can go
$d(x^2+x+1)=c(y^2+1)$
$d(4x^2+4x+4)=4c(y^2+1)$
$d(2x+1)^2+3d=4cy^2+4c$
$cd(2x+1)^2+3cd=(2cy)^2+4c^2$
$U^2-cdV^2=3cd-4c^2$, where $U=2cy$, $V=2x+1$. 
Now, that's not exactly a Pell equation, since the right side is $3cd-4c^2$ instead of $1$, but there are standard methods for solving this kind of equation (or showing it has no solution). 
A: Solutions of the equation: 
$$ax^2-by^2+cx-dy+q=0$$
you can record if the root of the whole: $k=\sqrt{(c-d)^2-4q(a-b)}$
Then using the solutions of the equation Pell: $p^2-abs^2=\pm1$
Then the formula of the solution, you can write:
$$x=\frac{\pm1}{2(a-b)}(((d-c)\pm{k})p^2+2(bk\mp(bc-ad))ps+b(a(d+c)-2bc\pm{ak})s^2)$$
$$y=\frac{\pm1}{2(a-b)}(((d-c)\pm{k})p^2+2(ak\mp(bc-ad))ps-a(b(d+c)-2ad\mp{bk})s^2)$$

Edited:
For the specific case 
$$x^2+x+1 = a(y^2+1)$$
let,
$$a = \frac{u}{v} = \frac{m^2+mn+n^2}{m^2+n^2}$$
Then,
$$x  = \frac{mp^2+2m\color{red}u\,pq+(m+n)uv\,q^2}{n} $$
$$y  = \frac{mp^2+(2m+n)\color{blue}v\,pq+muv\,q^2}{n} $$
and $p,q$ solve the Pell equation,
$$p^2-uv\,q^2=1$$
where one can easily set $n=1$ to get integer $x,y$.
A: Surely a zest of Galois theory could do no harm. 
1) Consider first the equation with rational parameter and variables : $x^2+x+1=a(y^2 +1), a, x, y \in \mathbf Q^*$  (1). With $i^2=-1$ and $j^3=1$, (1) is equivalent to $a= N_3(x-j).N_1(y-i)^{-1}$, where $N_n$ denotes the norm map of $\mathbf Q(\sqrt -n)/\mathbf Q$. The natural question is then how to characterize the elements of $\mathbf Q^*$ which are products of two norms as above. This is a purely Galois theoretic exercise (see e.g. Cassels-Fröhlich, ex. 4.4, p. 358)  : Let $K$ be any field of characteristic $\neq 2$. An element of $a \in K$ is a product of a norm from $K (\sqrt b)$ and a norm from $K(\sqrt c)$ iff $a$, as an element of $K(\sqrt bc)$, is a norm from $L=K(\sqrt b , \sqrt c)$  (2). In our case where $\mathbf Q(\sqrt -1)=\mathbf Q(i)$ and $\mathbf Q(\sqrt -3)=\mathbf Q(j)$, (2) is equivalent to "$a \in \mathbf Q (\sqrt 3)$ is a norm from $\mathbf Q (\sqrt 3)(i)$", or "$a$ is a sum of two squares in $\mathbf Q(\sqrt 3)$", say $a=(r+s\sqrt 3)^2 + (t+u\sqrt 3)^2$. Since $a \in \mathbf Q$, a straightforward calculation shows that the latter condition is equivalent to $a$ is of the form $r^2+t^2+3(s^2+u^2)$, with $r, s, t, u \in \mathbf Q$ such that $rs=-tu$  (3). Summarizing : (3) is equivalent to  (3b)  $a= (m^2+mn+n^2)/(\mu^2 +\nu^2)$, with rational parameter and variables. Cp. the edit of @Tito Pezias III.
2) Let us apply this approach to the special expression of $a$ given by (1). Elementary calculations show that the relation $rs=-tu$ in (3) boils down to $(2x+1)y=0$ and we have two different cases which correspond each to a single normic equation :  (i) $y=0$, $a= x^2+x+1= N_{3} (x-j)$ ; (ii) $x=-\frac 1 2$,  $3 a^{-1}=4(y^2 +1)= N_{1}(2(y-i))$. 
Let us finally come back to the original question, with $x, y \in \mathbf Z$. In case (i), $a$ is a norm from  $\mathbf Z[j]$, but this is only a necessary condition. In case (ii), $3 a^{-1}$ is a norm from $\mathbf Z[i]$, i.e. $3 a^{-1}$ is a sum of two squares of integers, but again this is only a necessary condition; actually the inequality $3/4a\ge1$ cannot hold if $a\in \mathbf Z$. In the end, the integrality condition imposed on the variables $x, y$ seems rather irrelevant. A more natural condition could be that $x, y\in \mathbf Z$ and $a\in\mathbf Q^*$ is defined up to a square, see (3b).
Edit : Because of the OP's uncertainty, I checked again the "elementary calculations" in part 2) of my answer... and I found a sign error which completely destroys the conclusion, so there is no reduction to two normic equations in one variable. I replace 2) by 3) and 4) below, following the same method but getting a different final result. I put in another answer because it seems that typo difficulties appear when the post is too long. 
A: Here are the announced corrections to my first answer.
3) Beforehand we must come back to the original question : "For which values of $a\in \mathbf Q^*$ do integer solutions to the equation $x^2+x+1=a(y^2+1)$ (1) exist ?" As pointed out by @Gerry Myerson, this formulation does not make much sense (just take $a$ to be a quotient of the desired form), so I guess that the OP actually meant : "What can we say about the integer solutions of this equation (1) when they exist ?"
 Let us first study the rational solutions of (1) as in part 1). The solutions to the more general equation (3b) are parametrized by equation (3) , which is just a way to express that $a=N(z)$, where $N$ is the norm map of $\mathbf Q (\sqrt 3)(i)/\mathbf Q (\sqrt 3)$ and $z=(r+s\sqrt 3) + i(t+u\sqrt 3)$. If we fix $a$ and a solution $z_0$, all the solutions to $a=N(z)$ will be of the form $z=Z z_0$, with $N(Z)=1$ (this amounts to determine the pre-image of $a$). The solutions of the latter normic equation are in turn parametrized by $R^2+T^2+3(S^2+U^2)=1$  (4), which expresses the representation of $1$ by a $4$-variable quadratic form. There is a well established theory for that (local- global principle, Hilbert symbols, etc., see e.g. CF, op. cit.), but let us use our previous parametrization. To simplify the notations, write $z=[r,s,t,u]$ when $z=(r+s\sqrt 3) +i(t+u\sqrt 3)$, and analogously $Z=[R,S,T,U]$. Then $z=z_0Z =[r_0R-t_0T+3(s_0S-u_0U), r_0s_0+s_0R-u_0T-t_0U, t_0R+r_0T+3(u_0S+s_0U),t_0S+u_0R+r_0U+r_0U+s_0T]  (5)$. A not very enlightening parametrization, but note that there is obviously an infinite number of solutions (in (4), just take $U=S=0$ and solve $R^2+T^2=1$) 
4) Switching from the general equation (3) to the original equation (1), we  adapt slightly the argument of 3) to exploit the particular expression of $a= N_3(x_0-j)/N_1(y_0-i)$ (notations of (1))= $N((x_0-j)/(y_0-i)$. Put $z_0=(x_0-j)/(y_0-i)$; note that by identification one can take $(y_o^{2}+1)z_0 = [y_0(x_0 -\frac 12), \frac 12, x_0+\frac 12,- \frac {y_0}2]$ (6). To get all the rational solutions $(x,y)$, the approach in 3) works without any change. To catch the integral solutions, use (5) and (6), although the computations are not very instructive.
Remark on the Pell-like equation.  Putting $z=x+\frac 12$, equation (1) becomes a "Pell-like" equation $z^2-ay^2=a- \frac 34$ (7). But the difficulty is that $a$ appears on both sides; this is related to the OP formulation as discussed at the beginning of part 3). To go on, take $a$ as a parameter given by  a solution  $(z_0,y_0)$ of (7). There are then two possible approaches : 
(i) Algebraically, as in part 3), we are led to solve a normic equation $N(z-y\sqrt d)=N(z_0-y_0\sqrt d)$, where $d=a-\frac 34$ and $N$ is the norm map of the quadratic extension $\mathbf Q(\sqrt d)/\mathbf Q$. By Hilbert's 90 (or direct computation) any element of norm $1$ in such an extension can be written under the form $(e-f\sqrt d)/(e+f\sqrt d)$, and we are brougt back to parametrize $(x,y)$ from $(z-y\sqrt d)/(z_0-y_0\sqrt d)=(e-f\sqrt d)/(e+f\sqrt d)$ by identification. This seems rather neat, but it's actually not, because $a$ is hidden under two disguises, $d$ and $(z_0, y_0)$. Besides, the symmetry in calculations is broken by the eviction of the fields $\mathbf Q (i)$ and $\mathbf Q(j)$, which played such a prominent role in parts 1) and 2).
(ii) Geometrically, the homogeneous form of (7), say $aY^2 + dT^2 = Z^2$, is the equation of a projective conic $C$ defined over $\mathbf Q$ and admitting a rational point $M_0$ (coming from $(x_0, y_0)$). Algebraic geometry then tells us that $C$ is isomorphic to the projective line $\mathbf P^1(\mathbf Q)$. The isomorphism is given by drawing the line joining $M_0$ to some fixed well chosen point of the projective plane $P^2(\mathbf Q)$, and then taking its intersection with $ P^1(\mathbf Q)$ (embedded in $P^2(\mathbf Q)$ as a coordinate axis). This is exactly the same kind of parametrization as for, say, the pythagorean triples. However the final parametrization does not look neater than (i). 
A: I just wanted to see just a few integral values of $a$ that work...( I know it asks for $a \in \mathbb{Q}$ )
$$x^2+x+1=a(y^2+1)$$
$$(2x+1)^2+3=4a(y^2+1)$$
$$\text{pell type equation} \qquad \to \qquad [2x+1]^2-a[2y]^2=4a-3$$
$$\alpha^2-a\cdot \beta^2=4a-3 \qquad \bigg{|} \qquad \alpha=2x+1 \qquad \beta = 2y$$
$$\implies$$
$$x=\frac{\alpha-1}{2}$$
$$y=\frac{\beta}{2}$$
$$\begin{align} 
& \underline{a \implies  (\alpha_n, \beta_n) \to (x_n,y_n)} \\ 
& 3 \implies (6,3) \to  \notin \mathbb{Z^2} \\ 
& \implies (21,12) \to (10,6) \\ 
& \implies (78,45) \to \notin \mathbb{Z^2} \\ 
& \implies (291,168) \to (145,84) \\ 
& \ \ \qquad \qquad \vdots \\ 
& 4 \implies (7,3) \to \notin \mathbb{Z^2} \\
& 7 \implies (635, 240) \to (317, 120) \\ 
& \ \ \qquad \qquad \vdots \\
& \vdots
\end{align}$$
With $a \in \mathbb{Q}$ as the OP suggested:
$$d(x^2+x+1)=c(y^2+1)$$
$$d(2x+1)^2+3d=4c(y^2+1)$$
$$\text{pell type equation} \qquad \to \qquad d \cdot [2x+1]^2-c \cdot [2y]^2=4c-3d$$
$$d \cdot \alpha^2-c\cdot \beta^2=4c-3d \qquad \bigg{|} \qquad \alpha=2x+1 \qquad \beta = 2y$$
$$\implies$$
$$x=\frac{\alpha-1}{2}$$
$$y=\frac{\beta}{2}$$
$$\begin{align} 
& \underline{(c,d) \implies  (\alpha_n, \beta_n) \to (x_n,y_n)} \\ 
& (7,1) \implies (635, 240) \to (317, 120) \\ 
& \vdots
\end{align}$$
A: You need the equation.  $$f(X^2+X+1)=e(Y^2+1)$$
To present in General form. To use this formula. For this we need to make the change.
$X=x+s$ ; $Y=y+t$
The equation then takes the form.
$$fx^2-ey^2+f(2s+1)x-2ety+f(s^2+s+1)-e(t^2+1)=0$$
You need to use the formula to solve this equation.
$$ax^2-by^2+cx-dy+q=0$$
Solvability is reduced to the solution of this equation.
$$k^2=(c-d)^2-4q(a-b)$$
$$k^2=((2s+1)f-2te)^2-4(f-e)((s^2+s+1)f-(t^2+1)e)$$
Then the solution can be written as.
$$e=p^2+(2s+1)pn+(s^2+s+1)n^2$$
$$f=p^2+2tpn+(t^2+1)n^2$$
$$k=(2t-2s-1)p^2+2(t^2-s(s+1))pn+((2s+1)(t^2+1)-2(s^2+s+1)t)n^2$$
The solution is given @Тито Piezas III a special case of the given formula.
A: This is just an oversized comment.
I certainly don’t disagree with any of the methods given, but I’ve just noticed some values of $a$ can be easily excluded by modular considerations. I expect this is just so obvious that it’s not worth mentioning, but here goes.
My assumption is $a=\frac{c}{d}$ where, critically, $gcd(c,d)=1$. If $x^2+x+1=u$ and $y^2+1=v$, then $a=\frac{u}{v}=\frac{c}{d}$.
Although there are a huge number of modular constrains, the most useful (IMHO) are when
$$c\not\equiv 0{\pmod j}$$
$$j=(2,5,11,17,23,29,41,47…)$$
I conjecture that the primes $(3n-1)$ follow http://oeis.org/A003627, but I’ve not had time to check.
$$d\not\equiv 0{\pmod k}$$
$$k=(4,3,7,11,19,23,31,43,47…)$$
I conjecture that the primes $(4n+3)$ follow http://oeis.org/A002145,  but again I’ve not had time to check.
