Given a real-valued matrix $A$, one can obtain its largest eigenvalue $\lambda_1$ plus the corresponding eigenvector $v_1$ by choosing a random vector $r$ and repeatedly multiplying it by $A$ (and rescaling) until convergence.
Once we have the first eigenpair, is there a similar way to estimate the second eigenpair?

Here is a theorem that you can use:
where $$\vec{x}=\frac{1}{\lambda_1 v_{1,k}}\begin{pmatrix} a_{k1}\\ a_{k2}\\ ...\\ a_{kn} \end{pmatrix}$$
$v_{1,k}$ is the $k$th component of $\vec{v}_1$, $a_{ki}$ is the $ki$th element of $A$. The row $k$ is smallest index such that $v_{1,k}$ is the infinity norm of $\vec{v}$, i.e., the largest component.
Then you can find the largest eigenvalue of $B$, which is the second largest of $A$.