Abstract: Based on nonlocal nonlinear Schrodinger equations of double beams, the evolutional rules of parameters and critical powers of transmission of orthogonal polarization center-coincidence double Hermite-Gaussian beams are calculated with variational approximation method in strongly nonlocal medium by expanding the response function in Taylors series to the second order. The evolutional rules of the beam width and phase shift are numerically simulated with split-step Fourier algorithm. The orthogonal polarization center-coincidence double Hermite-Gaussian spatial optical solitons and their large phase shifts evolutional rules are obtained when the two beams are incident with critical powers. The beams can form breathers, but the breathers become unstable with the increasing of the order when the two beams are incident with the total critical power but two incident powers are different. For breathers of each order, the beam width with high power periodically compacted oscillates; the beam width with low power periodically extended oscillates. Moreover, the phase shift of higher power breather increases faster with the increasing of transmission distance. When the order of Hermite-Gaussian is under the fifth order, the variational approximate solution agrees well with the numerical solution.