I am trying to derive the non-abelian Aharonov-Bohm effect by generalising Michael Berry's derivation to the case of non-abelian gauge field $A$.
My derivation so far
We require a degenerate eigenspace in order to achieve a non-abelian Berry phase therefore I take my Hilbert space to be $\mathcal{H} = \mathcal{H}_\text{spatial} \otimes \mathcal{H}_\text{internal}$, where $\mathrm{dim}(\mathcal{H}_\text{internal})=N$. The wavefunctions will take the form
$$\Psi(x,t) = \psi(x,t) \mathbf{v} ,$$
where $\psi(x,t) $ is the spatial wavefunction and $\mathbf{v} $ is the internal state vector of the system. I now take my Hamiltonian to be
$$ H(X) = - \frac{1}{2m } (\nabla \mathbb{I} - ie A)^2 + V(X-x)\mathbb{I}$$
where $V(X-x)$ is the confining potential which traps our particle inside a small box centred at position $X$, $A$ is our gauge field and $\mathbb{I}$ is the identiy on $\mathcal{H}_\text{internal}$. This Hamiltonian is almost identical to the Hamiltonian used in Berry's derivation, except now I have upgraded this to an operator on $\mathcal{H}$ by allowing $H$ to have internal indices too and allowing $A$ to be a non-abelian gauge field.
Generalising the result of Berry's paper, the $N$ eigenstates of the Hamiltonian with energy $E$ in a region where the curvature of $A$ vanishes is given by
$$ \Psi_j(X;x,t) =P \exp \left( - i \int_X^x A \cdot \mathrm{d} l \right) \psi_E(X;x,t) e_j $$ where $P$ stands for path ordering, $\psi_E$ is the spatial wavefunction with energy $E$ and $e_j$ are the basis vectors of $\mathcal{H}_\text{internal}$. This is easy to show as the differential operator $\nabla$ only acts on the spatial degrees of freedom, so we have one eigenstate for every basis vector $\mathbf{e}_j$ and hence our desired degeneracy required for a non-abelian Berry connection. The correseponding Berry connection is given by
$$ [\mathcal{A}_\mu]_{ij}(X) = i\langle \Psi_i(X) | \frac{\partial}{\partial X^\mu} | \Psi_j(X) \rangle \\ = i\int \mathrm{d}^n x e_i^\dagger \bar{P} \exp \left( i \int_X^x A \cdot \mathrm{d} l \right) (iA_\mu) P \exp \left( - i \int_X^x A \cdot \mathrm{d} l \right) e_j \psi_E^*(X;x,t) \psi_E(X;x,t)$$
where $\bar{P}$ is the anti-path ordering operator, which is due to taking the Hermitian conjugate. For the case of an abelian gauge field $A$, the exponentials would commute past everything and the Berry connection would reduce to $\mathcal{A} \propto A$, however I do not know how to evaluate this for the case of non-abelian connections.
My problem
Multiple sources suggest the non-abelian Aharonov-Bohm effect would yield a Wilson line of the gauge field,
$$ U = P \exp \left( -i \oint_C A \cdot \mathrm{d} l \right) $$ e.g. this and this, which suggests to me that the Berry connection is proportional to the gauge field, i.e. $\mathcal{A} \propto A$, however from my derivation I get stuck at the last line above where I am required to evaluate
$$ \bar{P} \exp \left( i \int_X^x A \cdot \mathrm{d} l \right) A_\mu P \exp \left( - i \int_X^x A \cdot \mathrm{d} l \right)=? $$
Is there some sort of generalised Baker-Campbell-Hausdorff formula for path ordered exponentials, i.e. something like $e^X Y e^{-X} = Y + [X,Y] + \frac{1}{2} [X,[X,Y]] + \ldots $?