Calculate the determinant of the matrix under the condition $x \neq y$

Calculate the determinant of the matrix:

$$\left|\begin{array}{ccccc} a_{1} & x & x & \ldots & x \\ y & a_{2} & x & \ldots & x \\ y & y & a_{3} & \ldots & x \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ y & y & y & \ldots & a_{n} \end{array}\right|.$$

Attempts and thoughts:

One of my good "colleagues" noticed that if $D_n$ is the determinant of an $n×n$ matrix, then

$$D_{n}=P_{n}-x y \sum_{i=0}^{n-2}(-1)^{i} P_{n-2-i} X_{i}.$$

Here $X_k$ is the sum of all products, $x^iy^j$, where $i+j=k$.

$P_k$ is the sum of all possible products involving exactly $k$ distinct $a_i$ values.

Manually counting the determinants, you can see a certain pattern:

$$ \begin{aligned} &D_{0}=P_{0} \\ &D_{1}=P_{1} \\ &D_{2}=P_{2}-x y[P_{0} X_{0}] \\ &D_{3}=P_{3}-x y[P_{1} X_{0}-P_{0} X_{1}] \\ &D_{4}=P_{4}-x y[P_{2} X_{0}-P_{1} X_{1}+P_{0} X_{2}]. \end{aligned} $$

This is true up to $n=10$. I did not check further. However, a better option appeared on the horizon, namely $\displaystyle\frac{x f(y)-y f(x)}{x-y}$, where $f(z)=(a_{1}-z)(a_{2}-z) \ldots(a_{n}-z)$. And verily, these two expressions agree with each other. Option $1$ is more complicated and requires additional calculations, while this formula is much simpler and fully satisfies the condition of the problem.

Question: I can prove the second formula by the method of mathematical induction. It is not very difficult, much easier than deriving this formula. And here's the question. How do you derive this formula?

Here is a proof without induction that $$\tag1 D_n=\frac1{x-y}\left(x\prod_{j=1}^n(a_j-y)-y\prod_{j=1}^n(a_j-x)\right)\mbox{if }x\neq y.$$ It could have been found, I think, without knowing the formula in advance, but of course, I had seen it in the question. This "derivation" is based on the fact that $D_n$ can be evaluated "easily" if $x=a_i$ for some $i$.

I give the derivation for distinct $a_i$, $i=1,...,n$ and $y\not\in\{a_1,\ldots,a_n\}$; for general $a_i$ and $y$ it follows by continuity.

Claim: If $1\leq i\leq n$ and $a_i=x$ then $D_n=a_i\displaystyle\prod_{1\leq j\leq n,j\neq i}(a_j-y)$.

Proof: I give it for $i=3$, for other $i$ it is similar. First, subtract row $3$ of $D_n$ from rows $1,2$. The resulting determinant has a $2\times(n-2)$-block of zeros in the upper right. It is therefore the product of the two diagonal blocks of size $2\times 2$ and $(n-2)\times(n-2)$, respectively. The determinant of the $2\times2$-block is $(a_1-y)(a_2-y)$. The lower right block is $$\left|\begin{array}{ccccc} x & x & x & \ldots & x \\ y & a_{4} & x & \ldots & x \\ y & y & a_{5} & \ldots & x \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ y & y & y & \ldots & a_{n} \end{array}\right|.$$ Factoring $x$ out of the first row and subtracting $y$ times the resulting first row from the other rows, we obtain an upper triangular matrix with diagonal elements $1,a_4-y,\ldots,a_n-y$. The determinant of the lower right block is therefore $x(a_4-y)\cdot\ldots\cdot(a_n-y)$. The product of both subdeterminants gives the result because $x=a_3$.

So we can regard $D_n=D_n(x)$ as a polynomial of $x$ of degree $\leq n-1$ (which follows from the expansion formula of determinants) with the values $$D(a_i)=a_i\prod_{j\neq i}(a_j-y), i=1,\dots,n.$$ This is enough to determine $D_n(x)$. The fact that there is a single missing factor, namely $a_i-y$, in the above product suggests to introduce the polynomial $P_n(x)=(x-y)D_n(x)$ of degree $\leq n$. It has the values $$P_n(a_i)=a_i\prod_{j=1}^n(a_j-y), i=1,\ldots n, P_n(y)=0.$$ The first $n$ values imply that $P_n(x)-x\prod_{j=1}^n(a_j-y)$ is a polynomal of $x$ of degree $\leq n$ with the $n$ zeros $a_i$. Here it is used that $a_1,\ldots,a_n$ are distinct. Hence it is a "constant" multiple of $\prod_{j=1}^n(a_j-x)$, i.e. $$P_n(x)=x\prod_{j=1}^n(a_j-y)+F(y)\prod_{j=1}^n(a_j-x)$$ with some $F(y)$ independent of $x$. The facts that $P_n(y)=0$ and, by assumption, $\prod_{j=1}^n (a_j-y)\neq0$ now give that $F(y)=-y$ and we have proved (1).

Remarks: 1. $D_n$ is a polynomial of $x,y$ and there is no problem with $x=y$. Formula (1), however, is only valid if $x\neq y$. If $x=y$ then $D_n$ equals the limit of formula (1) for $x\to y$. This limit is the derivative of the big parenthesis with respect to $x$ at the point $x=y$. The result is $D_n=f_n(y)-yf_n'(y)$, where $f_n(z)=\prod_{i=1}^n(a_i-z)$.
2. Using a single summation, formula (1) can be rewritten as an expression valid for all $x,y$. First $$D_n=\prod_{j=1}^n(a_j-y)+\frac y{x-y}\left(\prod_{j=1}^n(a_j-y)-\prod_{j=1}^n(a_j-x)\right)\mbox{if }x\neq y.$$ Now $$\prod_{j=1}^n(a_j-y)-\prod_{j=1}^n(a_j-x)=\sum_{i=0}^{n-1}\left(\prod_{j=1}^{i+1}(a_j-y)\prod_{j=i+2}^n(a_j-x)-\prod_{j=1}^i(a_j-y)\prod_{j=i+1}^{n}(a_j-x)\right)$$ where empty products are considered to be equal to 1. We find $$\prod_{j=1}^{i+1}(a_j-y)\prod_{j=i+2}^n(a_j-x)-\prod_{j=1}^i(a_j-y)\prod_{j=i+1}^{n}(a_j-x)=(x-y)\prod_{j=1}^{i}(a_j-y)\prod_{j=i+2}^n(a_j-x)$$ and finally $$D_n=\prod_{j=1}^n(a_j-y)+y\sum_{i=0}^{n-1}\prod_{j=1}^{i}(a_j-y)\prod_{j=i+2}^n(a_j-x)$$ where again empty products are considered to be equal to 1.