We know the formula of "family of circles intersecting at two points" as follows: $$x^2 + y^2 + D_1 x + E_1 y + F_1 + \lambda (x^2 + y^2 + D_2 x + E_2 y + F_2) = 0, \qquad \lambda \setminus \{-1\} \in \mathbb{R}$$ I could not find any rigorous derivation of the formula, which might involve theorems (linear independence, orthogonality of vectors by dot product) from linear algebra, vector algebra, etc., other than simply multiplying one of the circle equations known to intersect at two points by $\lambda$.
Question: How can we can derive the formula in a more rigorous way if possible?
You forgot to allow $\lambda=\infty$ to correspond to the circle $x^2+y^2+D_2x+E_2y+F_2=0$ as a candidate. A more convenient way is the affine linear combination $\lambda C_1+(1-\lambda)C_2=0$, $\lambda\in\mathbb{R}$ (the degenerate case $\lambda=\infty$ then corresponds to the radical axis).
An algebraic proof: A general circle is given by $x^2+y^2+Dx+Ey+F=0$, so in terms of coefficients, it makes sense to start by looking at the linear map $$ \begin{pmatrix} x(P_1)^2+y(P_1)^2 & x(P_1) & y(P_1) & 1\\ x(P_2)^2+y(P_2)^2 & x(P_2) & y(P_2) & 1\\ \end{pmatrix}\colon\mathbb{R}^4\to\mathbb{R}^2 $$ which obviously has rank 2. We are given the kernel containing both $$ \begin{pmatrix} 1\\D_1\\E_1\\F_1 \end{pmatrix} \text{ and } \begin{pmatrix} 1\\D_2\\E_2\\F_2 \end{pmatrix} $$ so every element of the kernel is a linear combination of them (since the kernel has dimension $2$). Now scalar multiplication doesn't change the circle, and obviously the zero vector doesn't give a circle. So projectivising gives $$ \lambda \begin{pmatrix} 1\\D_1\\E_1\\F_1 \end{pmatrix}+ (1-\lambda)\begin{pmatrix} 1\\D_2\\E_2\\F_2 \end{pmatrix} $$ i.e., $$ \lambda(x^2+y^2+D_1x+E_1y+F_1)+(1-\lambda)(x^2+y^2+D_2x+E_2y+F_2)=0 $$