equation solving - Solve returns a wrong answer

I'm afraid the following is rather localized questions, but I don't know how to simplify and generalize it. I have the following two points (in the real plane):

p = {Sqrt[b y1^2 + 4 ϵ], y1}
q = {(3 b y1^2 + 12 ϵ + Sqrt[5] y1 Sqrt[b (b y1^2 + 4 ϵ)])/(2 Sqrt[b y1^2 + 4 ϵ]), 
     (3 b y1 + Sqrt[5] Sqrt[b (b y1^2 + 4 ϵ)])/(2 b)}

In particular, I have also the following assumptions:

$Assumptions = b > 0 && y1 ∈ Reals && ϵ > 0;

Then I try to solve the following:

Solve[p[[2]] == -q[[2]] && p[[1]] == q[[1]], y1]

and obtain the answer

{{y1 -> -(Sqrt[ϵ]/Sqrt[b])}, {y1 -> Sqrt[ϵ]/Sqrt[b]}}

However, one can easily test that the second answer is wrong:

In: {p, q} /. % // Simplify

returns: $\left\{\left\{\left\{\sqrt{5} \sqrt{\epsilon },-\sqrt{\frac{\epsilon }{b}}\right\},\left\{\sqrt{5} \sqrt{\epsilon },\sqrt{\frac{\epsilon }{b}}\right\}\right\},\left\{\left\{\sqrt{5} \sqrt{\epsilon },\sqrt{\frac{\epsilon }{b}}\right\},\left\{2 \sqrt{5} \sqrt{\epsilon },\frac{4 \epsilon }{\sqrt{b \epsilon }}\right\}\right\}\right\}$

What am I missing?

Answer

You can understand what is happening by comparing:

s = Solve[p[[2]] == -q[[2]] && p[[1]] == q[[1]], y1]
(*
{{y1 -> -(Sqrt[ϵ]/Sqrt[b])}, {y1 -> Sqrt[ϵ]/Sqrt[b]}}
*)

with:

s = Solve[p[[2]] == -q[[2]] && p[[1]] == q[[1]], y1, Reals]
(*    
{{y1 -> ConditionalExpression[-Sqrt[(ϵ/b)], ϵ > 0 && b > 0]}}
*)

or:

s = Solve[p[[2]] == -q[[2]] && p[[1]] == q[[1]], y1, MaxExtraConditions -> All]

Solve::useq: The answer found by Solve contains equational condition(s) {0==Sqrt[5] Sqrt[b] Sqrt[ϵ]-Sqrt[5] Sqrt[b ϵ],0==-Sqrt[5] Sqrt[b] Sqrt[ϵ]-Sqrt[5] Sqrt[b ϵ]}. A likely reason for this is that the solution set depends on branch cuts of Mathematica functions. >>

{{y1 -> ConditionalExpression[-(Sqrt[ϵ]/Sqrt[b]), b != 0 && -Sqrt[5] Sqrt[b] Sqrt[ϵ] + Sqrt[5] Sqrt[b ϵ] == 0 && ϵ != 0]}, 
 {y1 -> ConditionalExpression[Sqrt[ϵ]/Sqrt[b], b != 0 && Sqrt[5] Sqrt[b] Sqrt[ϵ] + Sqrt[5] Sqrt[b ϵ] == 0 && ϵ != 0]}}