Consider the basic one-dimensional wave equation
$$\phi_{tt}(x,t)=\phi_{xx}(x,t)$$
with boundary conditions taken to be the Sommerfeld radiation conditions:
\begin{align} \phi_t(0,t)&=\phi_x(0,t)\\ \phi_t(1,t)&=-\phi_x(1,t) \end{align}
These BCs demand that there be no waves coming from the left or right, respectively. The general solution of the wave equation is well known to be
$$\phi(x,t)=l(x+t)+r(x-t)$$
where $l(x+t)$ represents a "left-moving" wave and $r(x-t)$ a "right-moving" wave.
I am following some (very old) notes where this PDE is solved numerically using finite difference methods. They start by "recasting the wave equation in first-order form (first order in time, first order in space)".
The first step is to introduce auxiliary variables $\psi(x,t)$ and $\rho(x,t)$:
\begin{align} \psi(x,t)&=\phi_x(x,t)\\ \rho(x,t)&=\phi_t(x,t) \end{align}
The wave equation then becomes the following pair of first-order equations:
\begin{align} \psi_t(x,t)&=\rho_x(x,t)\\ \rho_t(x,t)&=\psi_x(x,t) \end{align}
and the boundary conditions become
\begin{align} \psi_t(0,t)&=\psi_x(0,t)\\ \psi_t(1,t)&=-\psi_x(1,t)\\ \rho_t(0,t)&=\rho_x(0,t)\\ \rho_t(1,t)&=-\rho_x(1,t) \end{align}
Then these equation can be solved using a second-order Crank-Nicolson approximation, and so on.
My problem: I do not understand this process of reducing the PDE to first-order form. For instance,
Why does the single equation $\phi_{tt}(x,t)=\phi_{xx}(x,t)$ turn into two equations $\psi_t(x,t)=\rho_x(x,t)$ and $\rho_t(x,t)=\psi_x(x,t)$ after introducing the auxiliary variables? To be more specific, where does the equation $\psi_t(x,t)=\rho_x(x,t)$ come from?
Why do the 2 boundary conditions become 4 boundary conditions? How are the new boundary conditions derived? I thought the transformations would just look like $$\phi_t(0,t)=\phi_x(0,t)\quad\rightarrow\quad \rho(0,t)=\psi(0,t)$$ $$\phi_t(1,t)=-\phi_x(1,t)\quad\rightarrow\quad \rho(1,t)=-\psi(1,t)$$ but this is apparently very wrong.
I am not able to find many examples online of this method. If you can recommend a useful resource on this topic, it would be appreciated. I know of the reduction-of-order method for ODEs but using it for PDEs has me confused.
To answer your first question: It appears you have already recognized that $$\psi:=\phi_x\quad\&\quad\rho:=\phi_t\\\implies\psi_x=\phi_{xx}\quad\&\quad\rho_t=\phi_{tt}\\\implies\psi_x=\rho_t\qquad(*)$$ To obtain the second equation, we need to factor the operator from the original PDE. $$\begin{align}\phi_{tt}-\phi_{xx}&=\left(\partial_t^2-\partial_x^2\right)\phi=\left(\partial_t-\partial_x\right)\left(\partial_t+\partial_x\right)\phi\\&=\left(\partial_t-\partial_x\right)\left(\psi+\rho\right)=\psi_t+\rho_t-\psi_x-\rho_x=0\end{align}$$ We know from $(*)$ that $\psi_x-\rho_t=0$, thus $$\psi_t-\rho_x=\psi_x-\rho_t=0\implies\psi_t=\rho_x$$
To answer your second question: Because we modified the PDE of a single function to be a system with 2 functions, we need twice the number of conditions to make sure we properly satisfy the original PDE. In order to obtain these boundary conditions, we need to assume smoothness in $\phi$, meaning $$\phi_{xt}=\phi_{tx}$$ on the closure. Now take derivatives of the radiation conditions: $$\partial_t\left(\phi_t-\phi_x\right)=\phi_{tt}-\phi_{tx}=\rho_t-\rho_x=0\\\implies\rho_t(0,t)=\rho_x(0,t)$$ $$\partial_x\left(\phi_t-\phi_x\right)=\phi_{tx}-\phi_{xx}=\psi_t-\psi_x=0\\\implies\psi_t(0,t)=\psi_x(0,t)$$ $$\partial_t\left(\phi_t+\phi_x\right)=\phi_{tt}+\phi_{tx}=\rho_t+\rho_x=0\\\implies\rho_t(1,t)=-\rho_x(1,t)$$ $$\partial_x\left(\phi_t+\phi_x\right)=\phi_{tx}+\phi_{xx}=\psi_t+\psi_x=0\\\implies\psi_t(1,t)=-\psi_x(1,t)$$ This is not my favorite, but one of many approaches to solve the wave equation.