How to find the critical point for this Coulomb field

Two equal positive charges are at distance $d$, $-d$ from the origin on the $y$ axis. What is the distance on the $x$ axis beyond which a small perturbation in $y$ will move a particle away from the $x$ axis?

My attempt is to consider the y components of electric field where $y<d$ $$\frac{y-d}{(x^2+(y-d)^2)^{3/2}}+\frac{d+y}{(x^2+(d+y)^2)^{3/2}}=0$$ and solve for $x$ as $y/d \to 0$. But this gives me nonsense unsolvable equations.

The answer should be $x = \sqrt{2}d$.

As illustrated above, We are finding the horizontal distance to where the lines start to curve away.

Due to symmetry, the vertical electric field $E_y$ is always zero at $y = 0$.

$$E_y \equiv \frac{y+d}{(x^2+(y+d)^2)^{3/2}} + \frac{y-d}{(x^2+(y-d)^2)^{3/2}}$$

We want $\partial E_y / \partial y = 0$, where a test particle right on the $x$-axis is no longer stable. Namely, the goal is to find the $x$ value along the $x$-axis ($y = 0$) where the field $E_y$ goes from "squeeze in" $\partial E_y / \partial y < 0$ to "repel" $\partial E_y / \partial y > 0$ (note that we are talking about the field not the potential).

\begin{align} \frac{\partial E_y}{\partial y} &= \frac1{(x^2+(y+d)^2)^{3/2}} - \frac{3(y+d)^2}{(x^2+(y+d)^2)^{5/2}} \\ &\hspace{36pt} + \frac1{(x^2+(y-d)^2)^{3/2}} - \frac{3(y-d)^2}{(x^2+(y-d)^2)^{5/2}} \\ \implies \left. \frac{\partial E_y}{\partial y} \right|_{y=0} &= \frac2{(x^2+d^2)^{3/2}} - \frac{6 d^2}{(x^2+d^2)^{5/2}} = 2\frac{x^2 - 2d^2}{(x^2+d^2)^{5/2}} \end{align} Setting $(\partial E_y / \partial y)|_{y=0} = 0$ we easily have $x^2 - 2d^2 = 0 \implies x = \pm \sqrt{2} d$ as needed.

Alternatively, perturbation (series expansion) as OP intended

The goal is to show the critical point where the vertical electric field $E_y$ at given finite tiny $y>0$ goes from negative (towards the $x$-axis) to positive (away from $x$-axis), where

\begin{align} E_y &= \frac{d+y}{(x^2+(d+y)^2)^{3/2}} - \frac{d-y}{(x^2+(d-y)^2)^{3/2}} \\ &= \frac{d}{ (x^2 + d^2)^{3/2} } \left( \frac{ 1 + y/d }{ \displaystyle \Bigl( 1 + \frac{ 2dy + y ^2}{ x^2 + d^2 } \Bigr)^{3/2} } -\frac{ 1 - y/d }{ \displaystyle \Bigl( 1 - \frac{ 2dy - y ^2}{ x^2 + d^2 } \Bigr)^{3/2} } \right)\end{align}

Since we are only interested solving $E_y = 0$, the overall factor in front can be dropped for convenience.

$$\begin{aligned} E_y &\propto \bigl(1 + \frac{y}d \bigr) \left( 1 + \frac{ 2dy + y ^2}{ x^2 + d^2 } \right)^{-3/2} - \bigl(1 - \frac{y}d \bigr) \left( 1 - \frac{ 2dy - y ^2}{ x^2 + d^2 } \right)^{-3/2} \\ &= \bigl(1 + \frac{y}d \bigr) \left( 1 + \frac{ 2d }{ x^2 + d^2 } \bigl( 1 + \frac{y}{2d}\bigr)y \right)^{-3/2} - \bigl(1 - \frac{y}d \bigr) \left( 1 - \frac{ 2d }{ x^2 + d^2 } \bigl( 1 - \frac{y}{2d}\bigr)y \right)^{-3/2} \end{aligned}$$

Consider the series expansion for $\displaystyle \epsilon \equiv \frac{y}{2d} \ll 1$. Use the shorthand $\displaystyle K \equiv \frac{ 4d^2 }{ x^2 + d^2 }$ for this factor that doesn't involve $\epsilon$.

$$\begin{aligned} E_y &\propto \bigl(1 + \frac{y}d \bigr) \left( 1 + K \bigl( 1 + \frac{y}{2d}\bigr)\frac{y}{2d} \right)^{-3/2} - \bigl(1 - \frac{y}d \bigr) \left( 1 - K \bigl( 1 - \frac{y}{2d}\bigr) \frac{y}{2d} \right)^{-3/2}\\ &=(1 + 2\epsilon ) \bigl( 1 + K ( 1 + \epsilon ) \epsilon \bigr)^{-3/2} - (1 - 2\epsilon ) \bigl( 1 - K ( 1 - \epsilon ) \epsilon \bigr)^{-3/2}\\ &=(1 + 2\epsilon ) \left( 1 - \frac32 K (1+\epsilon)\epsilon + \frac1{2!}\frac32 \frac52 K^2 (1+\epsilon)^2\epsilon^2 - \frac1{3!}\frac32 \frac52 \frac72 K^3 (1+\epsilon)^3\epsilon^3 + \cdots \right) \\ &\hphantom{=} - (1 - 2\epsilon ) \left( 1 - \frac32 \bigl( -K (1-\epsilon)\epsilon \bigr) + \frac1{2!}\frac32 \frac52 \bigl( -K (1-\epsilon)\epsilon \bigr)^2 - \frac1{3!}\frac32 \frac52 \frac72 \bigl( -K (1-\epsilon)\epsilon \bigr)^3 + \cdots \right) \end{aligned}$$ Here we have $E_y \propto (1 + 2\epsilon) \cdot A_1 - (1 - 2\epsilon) \cdot A_2 = A_1 - A_2 + 2\epsilon \cdot ( A_1 + A_2)$ where $$\begin{aligned} A_1 &= 1 - \frac32 K (1 + \epsilon)\epsilon + \frac{15}8 K^2 (1 + \epsilon)^2\epsilon^2 - \frac{35}{16} K^3 (1 + \epsilon)^3\epsilon^3 + \cdots \\ A_2 &= 1 + \frac32 K (1-\epsilon)\epsilon + \frac{15}8 K^2 (1-\epsilon)^2\epsilon^2 + \frac{35}{16} K^3 (1-\epsilon)^3\epsilon^3 + \cdots \end{aligned}$$ so that $$\begin{aligned} A_1 - A_2 & = -3 K \epsilon + \frac{15}8 K^2 (4\epsilon)\epsilon^2 - \frac{35}{16} K^3 (2 + 6\epsilon^2)\epsilon^3 + \cdots \\ A_1 + A_2 & = 2 - 3 K \epsilon^2 + \frac{15}8 K^2 ( 2+2\epsilon^2)\epsilon^2 + \frac{35}{16} K^3 (6 \epsilon + 2\epsilon^3)\epsilon^3 + \cdots \end{aligned}$$

$$\bbox[5px,border:2px solid #4CAF50]{\implies\begin{aligned}[t] E_y &\propto -3 K \epsilon + O( \epsilon^3 ) + 2\epsilon \cdot \bigl( 2 + O(\epsilon^2) \bigr) \\ &= (4-3K)\epsilon + O( \epsilon^3 ) \end{aligned} }$$

Therefore the critical point $E_y = 0$ is given by \begin{align} 4 - 3K &= 0 &&\implies \frac{ 4d^2 }{ x^2 + d^2} = \frac43 && \implies 2d^2 = x^2 && \end{align} which is the desired $x = \sqrt{2} d$.