I found in the following paper: Comments on ‘‘A note on a three-term recurrence for a tridiagonal matrix’’ that we can compute the determinant of a block tridiagonal matrix A via a recursion.
In my particular case A is $4n\times4n$,
$$\textbf{A}=\begin{pmatrix} \textbf{B}_L-h\textbf{R} & J\space\textbf{R} & \textbf{0} & \cdots & \textbf{0} \\ J\space \textbf{R} & -h\textbf{R} & J\space\textbf{R} & & \textbf{0} \\ \textbf{0} & J\space \textbf{R} & -h\textbf{R} &\ddots &\vdots \\ \vdots & &\ddots &\ddots & J\space\textbf{R}\\ \textbf{0} & \textbf{0} & \cdots & J \space\textbf{R} & \textbf{B}_R-h\textbf{R} \end{pmatrix}$$
where $$\textbf{R}=\begin{pmatrix} 0&0&1&0\\ 0&0&0&1\\ -1&0&0&0\\ 0&-1&0&0\\ \end{pmatrix},\space \textbf{B}_{L,R}=\begin{pmatrix} 0&\frac{i}{2}\Gamma_{+}^{\text{L,R}}&-\frac{i}{2}\Gamma_{-}^{\text{L,R}}&\frac{1}{2}\Gamma_{-}^{\text{L,R}}\\ -\frac{i}{2}\Gamma_{+}^{\text{L,R}}&0&\frac{1}{2}\Gamma_{-}^{\text{L,R}}&\frac{i}{2}\Gamma_{-}^{\text{L,R}}\\ \frac{i}{2}\Gamma_{-}^{\text{L,R}}&-\frac{1}{2}\Gamma_{-}^{\text{L,R}}&0&\frac{i}{2}\Gamma_{+}^{\text{L,R}}\\ -\frac{1}{2}\Gamma_{-}^{\text{L,R}}&-\frac{i}{2}\Gamma_{-}^{\text{L,R}}&-\frac{i}{2}\Gamma_{+}^{\text{L,R}}&0\\ \end{pmatrix} $$ and $$J,h,\Gamma_{+}^{\text{L,R}},\Gamma_{-}^{\text{L,R}} \in \mathbb{R}$$
Now let me state the recursion mentioned in the above paper,
$$\text{det}(\textbf{A})=\prod_{k=1}^{n}\text{det}(\Lambda_{k})\space\space\space\space(1)$$
where (in my case),
$$ \Lambda_{1} = \textbf{B}_L-h\textbf{R}\\ \Lambda_{k} = -h\textbf{R}-J^{2}\textbf{R}\Lambda_{k-1}^{-1}\textbf{R}\\ \Lambda_{n}=\textbf{B}_R-h\textbf{R}-J^{2}\textbf{R}\Lambda_{n-1}^{-1}\textbf{R} $$ Now according to (1), the set of eigenvalues of $\textbf{A}$ should contain the eigenvalues of $\Lambda_{1} = \textbf{B}_L-h\textbf{R}$ since applying (1) to $\textbf{A}-\lambda I_{4n}$ gives $\Lambda_{1}^{'} = \textbf{B}_L-h\textbf{R}-\lambda I_{4}$.
Now here comes my issue.
I computed the spectrum of A in Mathematica for $n=50$ for the values $h=1,J=1.5,\Gamma_{+}^{\text{L}}=1.6,\Gamma_{+}^{\text{R}}=1.3,\Gamma_{-}^{\text{L}}=-0.4,\Gamma_{-}^{\text{R}}=-0.7$
I found that the spectrum did not contain the eigenvalues of $\textbf{B}_{L}-h\textbf{R}$. My question is: Is this recursion not applicable to my case, or am I incorrect in stating that the set of eigenvalues of A should contain the eigenvalues of $\Lambda_{1}$?
The recursion in the paper is a generalization of the classical determinant relation for block matrices (see Wiki): if $A$ is invertible then $$ \det\begin{bmatrix}A & B\\C & D\end{bmatrix}=\det A\cdot\det(D-CA^{-1}B). $$ The block matrix does not share eigenvalues with the $A$ block in general, since $\lambda I-A$ becomes not invertible for an eigenvalue $\lambda$ of $A$, hence, the formula cannot be applied.
Edit: there is an alternative formula for evaluation of the determinant of a tridiagonal matrix here (look Theorem 2). It uses only inversions of off-diagonal blocks that do not contain $\lambda$ in your case.