How does this opamp phase shifter produce 2x phase shift at the output without changing the magnitude?

Well, we are trying to analyze the following circuit (assuming an ideal opamp):

^{simulate this circuit – Schematic created using CircuitLab}

When we use and apply KCL, we can write the following set of equations:

$$ \begin{cases} \text{I}_\text{x}=\text{I}_1+\text{I}_3\\ \\ \text{I}_1=\text{I}_2\\ \\ 0=\text{I}_3+\text{I}_4 \end{cases}\tag1 $$

When we use and apply Ohm's law, we can write the following set of equations:

$$ \begin{cases} \text{I}_1=\frac{\text{V}_\text{x}-\text{V}_1}{\text{R}_1}\\ \\ \text{I}_2=\frac{\text{V}_1}{\text{R}_2}\\ \\ \text{I}_3=\frac{\text{V}_\text{x}-\text{V}_2}{\text{R}_3}\\ \\ \text{I}_4=\frac{\text{V}_3-\text{V}_2}{\text{R}_4} \end{cases}\tag2 $$

Substitute \$(2)\$ into \$(1)\$, in order to get:

$$ \begin{cases} \text{I}_\text{x}=\frac{\text{V}_\text{x}-\text{V}_1}{\text{R}_1}+\frac{\text{V}_\text{x}-\text{V}_2}{\text{R}_3}\\ \\ \frac{\text{V}_\text{x}-\text{V}_1}{\text{R}_1}=\frac{\text{V}_1}{\text{R}_2}\\ \\ 0=\frac{\text{V}_\text{x}-\text{V}_2}{\text{R}_3}+\frac{\text{V}_3-\text{V}_2}{\text{R}_4} \end{cases}\tag3 $$

Now, when we have an ideal opamp we know that \$\text{V}_\alpha:=\text{V}_+=\text{V}_-=\text{V}_1=\text{V}_2\$. So we can rewrite equation \$(3)\$ as follows:

$$ \begin{cases} \text{I}_\text{x}=\frac{\text{V}_\text{x}-\text{V}_\alpha}{\text{R}_1}+\frac{\text{V}_\text{x}-\text{V}_\alpha}{\text{R}_3}\\ \\ \frac{\text{V}_\text{x}-\text{V}_\alpha}{\text{R}_1}=\frac{\text{V}_\alpha}{\text{R}_2}\\ \\ 0=\frac{\text{V}_\text{x}-\text{V}_\alpha}{\text{R}_3}+\frac{\text{V}_3-\text{V}_\alpha}{\text{R}_4} \end{cases}\tag4 $$

Now, for the output voltage we get:

$$\text{V}_3=\frac{\text{V}_\text{x}\left(\text{R}_2\text{R}_3-\text{R}_1\text{R}_4\right)}{\text{R}_3\left(\text{R}_1+\text{R}_2\right)}\tag{5}$$

---

Now, applying this to your circuit we need to use (from now on I use the lower case letters for the function in the 'complex' s-domain where I used Laplace transform):

$$\text{R}_1=\frac{1}{\text{sC}}\tag6$$

So, the output voltage will be:

$$\text{v}_3\left(\text{s}\right)=\frac{\text{v}_\text{x}\left(\text{s}\right)\left(\text{R}_2\text{R}_3-\frac{\text{R}_4}{\text{sC}}\right)}{\text{R}_3\left(\frac{1}{\text{sC}}+\text{R}_2\right)}\tag7$$

So, when we use the transformation \$\text{s}=\text{j}\omega\$ (where \$\text{j}^2=-1\$), we get:

$$\underline{\text{v}}_3\left(\text{j}\omega\right)=\frac{\underline{\text{v}}_\text{x}\left(\text{j}\omega\right)\left(\text{R}_2\text{R}_3-\frac{\text{R}_4}{\text{j}\omega\text{C}}\right)}{\text{R}_3\left(\frac{1}{\text{j}\omega\text{C}}+\text{R}_2\right)}\tag8$$

So, we can write a transfer function:

$$\underline{\mathcal{H}}\left(\text{j}\omega\right):=\frac{\underline{\text{v}}_3\left(\text{j}\omega\right)}{\underline{\text{v}}_\text{x}\left(\text{j}\omega\right)}=\frac{\text{R}_2\text{R}_3-\frac{\text{R}_4}{\text{j}\omega\text{C}}}{\text{R}_3\left(\frac{1}{\text{j}\omega\text{C}}+\text{R}_2\right)}=\frac{\text{R}_2\text{R}_3+\frac{\text{R}_4}{\omega\text{C}}\cdot\text{j}}{\text{R}_3\left(\text{R}_2-\frac{1}{\omega\text{C}}\cdot\text{j}\right)}\tag9$$

So, for the magnitude we will get:

$$\left|\underline{\mathcal{H}}\left(\text{j}\omega\right)\right|=\frac{\sqrt{\left(\text{R}_2\text{R}_3\right)^2+\left(\frac{\text{R}_4}{\omega\text{C}}\right)^2}}{\text{R}_3\sqrt{\text{R}_2^2+\left(\frac{1}{\omega\text{C}}\right)^2}}\tag{10}$$

And, for the phase we will get:

$$\arg\left(\underline{\mathcal{H}}\left(\text{j}\omega\right)\right)=\arctan\left(\frac{\text{R}_4}{\omega\text{C}}\cdot\frac{1}{\text{R}_2\text{R}_3}\right)-\left(\frac{3\pi}{2}+\arctan\left(\text{R}_2\omega\text{C}\right)\right)\tag{11}$$

---

Now, when \$\text{R}:=\text{R}_3=\text{R}_4\$ all this simplifies to:

$$\left|\underline{\mathcal{H}}\left(\text{j}\omega\right)\right|=\frac{\sqrt{\left(\text{R}_2\text{R}\right)^2+\left(\frac{\text{R}}{\omega\text{C}}\right)^2}}{\text{R}\sqrt{\text{R}_2^2+\left(\frac{1}{\omega\text{C}}\right)^2}}=1\tag{12}$$

And:

$$\arg\left(\underline{\mathcal{H}}\left(\text{j}\omega\right)\right)=-\pi-2\arctan\left(\text{R}_2\omega\text{C}\right)\tag{13}$$

How is the gain of this circuit ±1 for all frequencies?

Using superposition we can immediately say that the gain is: -

$$H(s) = (1+\dfrac{R'}{R'})\cdot (\dfrac{1}{1+sCR}) - \dfrac{R'}{R'}$$ $$ = \dfrac{2}{1+sCR} - 1$$ $$ = \dfrac{1- sCR}{1+sCR}$$

And, because the magnitude of \$1 - sCR\$ is the same as the magnitude of \$1+ sCR\$, the gain magnitude is 1.

The phase change in the numerator and the phase change in the denominator go in opposite directions so, given that they are disposed as "numerator" and denominator" any phase angle from the circuit is double what it would be for a single RC network.