I think the plot in question originally comes from Wallis (2021) so I recommend you have a look at that book for further explanation on the particulars of that graphical representation. In any case, the main reason why the Wilson score interval is superior to the classical Wald interval is that is is derived by solving a quadratic inequality for the proportion parameter that leads to an interval that respects the true support of the parameter. Contrarily, the Wald interval can go outside the true support, and it also has worse coverage properties (see Brown, Cai and DasGupta (2001) for further discussion).
It might help here to show you the derivation of the interval in algebraic terms. Suppose we have $n$ binary data values giving the sample proportion $p_n$ (which we will treat as a random variable) and let $\theta$ be the true proportion parameter. The classical Wald interval uses the asymptotic pivotal distribution:
$$\sqrt{n} \cdot \frac{p_n-\theta}{\sqrt{\theta(1-\theta)}} \overset{\text{Approx}}{\sim} \text{N}(0,1).$$
For the Wilson score interval we first square the pivotal quantity to get:
$$n \cdot \frac{(p_n-\theta)^2}{\theta(1-\theta)} \overset{\text{Approx}}{\sim} \text{ChiSq}(1).$$
Let $\chi_{1,\alpha}^2$ denote the critical point of the chi-squared distribution with one degree-of-freedom (with upper tail area $\alpha$). For any confidence level $1-\alpha$ we then have the probability interval:
$$\begin{align}
1-\alpha
&\approx \mathbb{P} \Big( n (p_n-\theta)^2 \leqslant \chi_{1,\alpha}^2 \theta(1-\theta) \Big) \\[6pt]
&= \mathbb{P} \Big( n (p_n^2 - 2 p_n \theta + \theta^2) \leqslant \chi_{1,\alpha}^2 (\theta-\theta^2) \Big) \\[6pt]
&= \mathbb{P} \Big( (n + \chi_{1,\alpha}^2) \theta^2 - (2 n p_n + \chi_{1,\alpha}^2) \theta + n p_n^2 \leqslant 0 \Big) \\[6pt]
&= \mathbb{P} \Bigg( \theta^2 - 2 \cdot\frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \cdot \theta + \frac{n p_n^2}{n + \chi_{1,\alpha}^2} \leqslant 0 \Bigg) \\[6pt]
&= \mathbb{P} \Bigg( \bigg( \theta - \frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \bigg)^2 \leqslant \frac{\chi_{1,\alpha}^2 (n p_n (1-p_n) + \tfrac{1}{4} \chi_{1,\alpha}^2)}{(n + \chi_{1,\alpha}^2)^2} \Bigg) \\[6pt]
&= \mathbb{P} \Bigg( \theta \in \Bigg[ \frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \pm \frac{\chi_{1,\alpha}}{n + \chi_{1,\alpha}^2} \cdot \sqrt{n p_n (1-p_n) + \tfrac{1}{4} \chi_{1,\alpha}^2} \Bigg] \Bigg), \\[6pt]
\end{align}$$
and substitution of the observed sample proportion (for simplicity I will use the same notation for this value) then leads to the Wilson score interval:
$$\text{CI}_\theta(1-\alpha) = \Bigg[ \frac{n p_n + \tfrac{1}{2} \chi_{1,\alpha}^2}{n + \chi_{1,\alpha}^2} \pm \frac{\chi_{1,\alpha}}{n + \chi_{1,\alpha}^2} \cdot \sqrt{n p_n (1-p_n) + \tfrac{1}{4} \chi_{1,\alpha}^2} \Bigg].$$
As you can see, solving the quadratic inequality in the probability interval leads to an interval that respects the true space of possible values of the proportion parameter (i.e., it is between zero and one). This is a major advantage of this method but it also has better coverage properties in general.