Exact global phase of the Gross-Pitaevskii wavefunction

I will follow closely the book by Kaspar Sakmann chapter 3.2 about derivation of the time-dependent Gross-Pitaevskii equation. Starting from the time-dependent variational principle we arrive at the equation: $$N\left[h(x) + \lambda_0(N-1)|\phi(x,t)|^2 \right]\phi(x,t) = \mu(t)\phi(x,t) + Ni\hbar\dot{\phi}(x,t)$$ where $$h(x) = -\hbar^2/(2m)d^2/dx^2$$, $$N$$ - number of particles and $$\lambda_0$$ is a constant parameter. Lagrange multiplier $$\mu(t)$$ is introduced here to ensure unit normalization (it was always there). Now, one can redefine the orbital $$\phi(x,t) = \tilde{\phi}(x,t)e^{\frac{i}{\hbar}\int\limits_{0}^{t}dt' \mu(t')/N}$$ and get the Gross-Pitaevskii equation that always appears in the textbooks: $$\left[h(x) + \lambda_0(N-1)|\tilde{\phi}(x,t)|^2 \right]\tilde{\phi}(x,t)= i\hbar\dot{\tilde{\phi}}(x,t)$$ Of course global change of the phase does not change physical quantities. However, I am dealing with more sophisticated problem and the true evolution of the function $$\phi(x,t)$$ is needed and change of the phase will lead to different result.

As K. Sakmann wrote in his book one can determine $$\mu(t)$$ bu multiplying the equation by $$\int dx' \phi(x',t)^*$$ from the left and putting expression for $$\mu(t)$$ back to the equation. At the end we get integro-differential equation:

$$i\hbar \left[\dot{\phi}(x,t) - \phi(x,t)\int dx'\phi(x',t)^* \dot{\phi}(x',t) \right] = \left[h(x) + \lambda_0(N-1)|\phi(x,t)|^2 \right]\phi(x,t) - \phi(x,t)\left[ \int dx' \phi(x',t)^* h(x') \phi(x',t) + \lambda_0(N-1)\int dx' |\phi(x',t)|^4\right]$$

which is very hard to solve even numerically. Equation for $$\tilde{\phi}(x,t)$$ is much simpler and can be solved e.g. numerically. My question is, is it possible to determine $$\mu(t)$$ by having the solution for $$\tilde{\phi}(x,t)$$? In other words, by having the solution for $$\tilde{\phi}(x,t)$$ I want to recover the solution for $$\phi(x,t)$$ through time-dependent phase factor multiplication. What should be the phase?

Is it possible to give expression for $$\mu(t)$$ just in terms of $$\tilde{\phi}(x,t)$$ and $$\tilde{\phi}^*(x,t)$$, but without time derivatives?

UPDATE

To give you a flavor what the real problem is I will refer to the following paper by A. Sinatra and I. Castin. They consider two-component condensate and want to find time-evolution of the initial particle number state $$|N_a : \phi_0, N_b : \phi_0\rangle$$. The ansatz for the time-evolution is the following (like in the MCTDHB) $$|\psi(t)\rangle = e^{-i/\hbar A(t)}|N_a : \phi_a(t), N_b : \phi_b(t)\rangle$$ Orbitals $$\phi_a,b$$ are solutions of the coupled GPEs: $$i\hbar \dot{\phi_a}(x,t) = \left[h(x) + N_ag_{aa}|\phi_a(x,t)|^2 + N_b g_{ab}|\phi_b(x,t)|^2 \right]\phi_a(x,t)$$ $$i\hbar \dot{\phi_b}(x,t) = \left[h(x) + N_b g_{bb}|\phi_b(x,t)|^2 + N_a g_{ab}|\phi_a(x,t)|^2 \right]\phi_b(x,t)$$ Note that there are no Lagrange multipliers there so these equations correspond to the second equation of the single component GP. Lagrange multipliers were removed by the phase transformation of orbitals which results in global phase $$A(t)$$. Now, the next step is enigmatic for me and this is the point I don't understand. In order to find the phase $$A(t)$$ the say that the state $$|\psi(t)\rangle$$ solves Schrodinger equation $$i\hbar \langle \psi(t)|\partial_t |\psi(t)\rangle = \langle \psi(t) |\hat{H}|\psi(t)\rangle$$ $$\dot{A} = \langle N_a : \psi_a(t), N_b : \psi_b(t) |\hat{H} - i\hbar \partial_t|N_a : \psi_a(t), N_b : \psi_b(t)\rangle$$

Is it really correct that the true state, the one you put into Schrodinger action with Lagrange multipliers solves also the Schrodinger equation?

After some thinking and looking at the derivation of the MCTDHB equations I think I get the answer. Let assume we have the following ansatz for the wavefunction (single-component): $$|\psi(t)\rangle = C(t)|N: \psi(t)\rangle$$
Schrodinger action is equal to: $$S[C, \psi] = \int dt \langle \psi(t) | \hat{H} - i\hbar \partial_t | \psi(t)\rangle - \mu(t)\left[\langle\psi(t)|\psi(t)\rangle - 1 \right] - \epsilon(t)\left[ |C(t)|^2 - 1\right]$$
As the following paper says we get at the end the equations: $$N\left[h(x) + \lambda_0(N-1)|\phi(x,t)|^2 \right]\phi(x,t) = \mu(t)\phi(x,t) + Ni\hbar\dot{\phi}(x,t)$$ $$\langle N:\phi(t)| \hat{H} - i\hbar\partial_t | N:\phi(t)\rangle C(t) = i\hbar \dot{C}(t)$$ Now, we can change orbital $\psi$ and reversely coefficient $C$ (we don't want to change initial function we have) by a phase factor to eliminate the Lagrange multiplier: $$\phi(x,t) = \tilde{\phi}(x,t)e^{i\hbar \int dt'\mu(t')}$$ $$C(t) = \tilde{C}(t)e^{-Ni\hbar \int dt'\mu(t')}$$
One can easily check that $$\langle N:\tilde{\phi}(t)| \hat{H} - i\hbar\partial_t | N:\tilde{\phi}(t)\rangle \tilde{C}(t) = i\hbar \dot{\tilde{C}}(t)$$ and we get $$\tilde{C}(t) = e^{-i/\hbar\int dt' \langle N:\tilde{\phi}(t')| \hat{H} - i\hbar\partial_{t'} | N:\tilde{\phi}(t')\rangle}C_0$$ which can be now calculated because we don't have Lagrange multiplier in the equation for $\tilde{\phi}$. This is what I was looking for!