Goal: Ultimately, I would like to find a trustworthy approximation for the ratio of the functional determinant of two differential operators using the formula
$$ frac{text{Det}(hat{D}_0)}{text{Det}(hat{D})} =frac{prod^infty_i lambda_{0,i}}{prod^infty_j lambda_j} approx frac{prod^N_i lambda_{0,i}}{prod^N_j lambda_j} $$
where $lambda$ are the eigenvalues of the respective operators and $N$ is the number of eigenvalues used in the approximation.
Step 1: defining the operators.
To define these two operators I first must solve the following coupled ODE boundary value problem in $phi(t)$ and $p_phi(t)$. The system depends on the parameter $tau > 0$.
bvp = {D(p(Phi)(t), t) == ((1/(Tau)))*p(Phi)(t)^2*(Sin(2*(Phi)(t))/2) +
(1/(Tau))*p(Phi(t)*Cos(2*(Phi)(t)) + (1/(Tau))*(Sin(2*(Phi)(t))/2),
D((Phi)(t), t) == 2*Pi  ((1/(Tau))*Sin(2*(Phi)(t)))/2 +
(1/(Tau))*p(Phi)(t)*Sin((Phi(t))^2,
(Phi)(0) == 2*ArcTan((1 + Sqrt(1 + 4*Pi^2*(Tau)^2))/(2*Pi*(Tau))),
(Phi)(1) == 2*ArcTan((1 + Sqrt(1 + 4*Pi^2*(Tau)^2))/(2*Pi*(Tau))) + 2*Pi}
Which is readily solved with
sol = ParametricNDSolveValue(bvp, {(Phi)(t), p(Phi)(t)}, {t, 0, 1}, {(Tau)})
Using these solutions I can define the two differential operators (which depend on a choice of parameter $tau$ and the numerical solution of $phi(t)$) with Dirichlet boundary conditions
$hat{D}_0$
{D((Tau)*Csc(sol((Tau))((1)))^2*D((Psi)(t), t), t), (Psi)(0) ==
0, (Psi)(1) == 0}
and $hat{D}$
{D((Tau)*Csc(sol((Tau))((1)))^2*D((Psi)(t), t), t) + (2*(Tau)*Cot(sol((Tau))
((1)))*Csc(sol((Tau))((1)))^2*D(sol((Tau))((1)), {t, 2}) 
(Tau)*Csc(sol((Tau))((1)))^4*(Cos(2*sol((Tau))((1))) + 2)*D(sol((Tau))((1)), t)^2 +
2*Pi*Csc(sol((Tau))((1)))^4*(2*Pi*(Tau)*(Cos(2*sol((Tau))((1))) + 2)  Sin(2*sol((Tau))
((1)))))*(Psi)(t), (Psi)(0) == 0,
(Psi)(1) == 0}
Step 2: Numerically calculating the eigenvalues.
First selecting a particular value of the parameter $tau$ to be $0.1$, I use the inbuilt function NDEigenvalue to numerically calculate N=39 Eigenvalues for $hat{D}_0$ and multiply them together
Times @@ NDEigenvalues({D(0.1*Csc(sol(0.1)((1)))^2*D((Psi)(t), t), t), (Psi)(0) == 0,
(Psi)(1) == 0}, (Psi)(t), {t, 0, 1}, 39)
Similarly for $hat{D}$
Times @@ NDEigenvalues({D(0.1*Csc(sol(0.1)((1)))^2*D((Psi)(t), t), t) +
(2*0.1*Cot(sol(0.1)((1)))*Csc(sol(0.1)((1)))^2*D(sol(0.1)((1)), {t, 2}) 
0.1*Csc(sol(0.1)((1)))^4*(Cos(2*sol(0.1)((1))) + 2)*D(sol(0.1)((1)), t)^2 +
2*Pi*Csc(sol(0.1)((1)))^4*(2*Pi*0.1*(Cos(2*sol(0.1)((1))) + 2)  Sin(2*sol(0.1)((1)))))*(Psi)(t), (Psi)(0) == 0,
(Psi)(1) == 0}, (Psi)(t), {t, 0, 1}, 39)
Then I simply take the ratio. I do this also for additional choices of the parameter $tau$.
The issue: The accuracy of this approximation is not good enough for small values of the parameter $tau$.
From the underlying physics, I strongly expect that
 The ratio will be a monotonically increasing function of $tau$
Unfortunately, after calculating a few of these ratios in the parameter range $ 0<tau<1 $, I find that the ratios do not appear to be monotonically increasing. Additionally, the order of magnitude seems to change significantly as I change $N$, the number of eigenvalues used in the approximation. These two things make me conclude the numerical solutions I have so far are untrustworthy.
If I try to increase the number of Eigenvalues used in the approximation beyond 39 I get an error message saying “A maximum number of 39 eigenvalues and functions can be computed for this discretized system.”
Solution Attempt:
My main guess at how to solve this issue is to find a way to increase the number of eigenvalues that NDEigenvalue returns.

I have tried to increase the accuracy of the parametric solutions for $phi(t)$, but this did not seem to help very much

I have also tried to change methods used in NDEigenvalues (see the code below)
NDEigenvalues({D(0.1*Csc(sol(0.1)((1)))^2*D((Psi)(t), t),t), (Psi)(0) == 0, (Psi)(1) == 0}, (Psi)(t), {t, 0, 1}, 70, Method > {"SpatialDiscretization" > {"FiniteElement","MeshOptions" > {"MaxCellMeasure" > 0.0001}}})
This will return about 60 eigenvalues before returning another error message (see image below) but has the same issue as the original method and will also sometimes find eigenvalues that when compared with the set of the first 39 eigenvalues found by the default method are sometimes an order of magnitude different!
So far all of this code runs fairly quickly, so I feel that there should be a way to get Mathematica to think for longer and find more eigenvalues but so far I have not found how to do this in a satisfactory way. I do not know if there are other ways I should be trying to solve my problem and would very much appreciate any advice I can get on how to make progress.