Can I break the degeneracy of energy eigenstates if I know which irrep of the group they transform as? Say I have a quantum system with a symmetric potential, whose symmetry is described by a group $G$. I know the character table of $G$, its irreducible representations, can work out the projection operators $\Pi_j$ etc.
With imaginary time evolution, I can find the spatial part of the energy eigenstates $\phi_{E_i}$. If there are degeneracies, however, what I will get is the sum of the degenerate eigenstates at the same $E_i$: $\psi_{E_i} = \sum_j \phi^{(j)}_{E_i}$.
Question: if I know the group, the irrep, etc., can I decompose/break $\psi$ into the individual $\phi^{(j)}_{E_i}$?
Reason for the question: From this answer:

Suppose there is a  group of transformations $G$. Then it acts on the Hilbert space by some
set of unitary transformations $\mathcal{O}$. The Hilbert space is therefore a
representation of the group $G$, and it splits up into subspaces of
irreducible representations (irreps). The important thing is that if
$|\psi\rangle$ and $|\phi\rangle$ are in the same irrep iff you can get from one to the
other by applying operators $\mathcal{O}$.

So another way of phrasing my question would be: can I somehow get $\mathcal{O}$, such that $\phi^{(2)}_{E_i} = \mathcal{O}\phi^{(1)}_{E_i}$ and $\psi = \sum_j \phi^{(j)}_{E_i} = \sum_j \mathcal{O}^j\phi^{(1)}_{E_i}$ ?

Example
The 2D quantum harmonic oscillator, looking at the states with energy $E = 2\hbar\omega$ which are $\psi_1(x,y) = \phi_0(x)\phi_1(y)$ and $\psi_2(x,y) = \phi_0(y)\phi_1(x)$ where $0$ is the ground state and $1$ the first excited state.
I know that $|\psi_1|^2$ and $|\psi_2|^2$ should look like this:

but, from my code, I get the spatial distribution of the "overall" energy level $E=2$ so I get $|\psi_1+\psi_2|^2$:

The group and irrep information about the 2D harmonic oscillator is from here:

The set of states with total number $m$ of excitation span the irrep
$(m,0)$ of $SU(2)$. Thus the degeneracy is the dimension of this irrep
[...] this is just $m+1$.

With this info, can I get $\psi_1$ and $\psi_2$ from $\psi_1 + \psi_2$?
EDIT:
To clarify, I want to decompose $\psi$ into non-degenerate eigenstates, not $|\psi|^2$. I am just plotting $|\psi|^2$ instead of $\psi$ for simplicity.
 A: *

*It is not possible to “break the degeneracy” by combining degenerate states: any unitary transformation inside the degenerate subspace will produce a different set of eigenstates of $H$, but they will all have the same eigenvalues.


*This does not prevent you from organizing the states in your (degenerate) subspace into irreps of some group.  Presumably the elements of this group commute with the Hamiltonian so you will get legitimate eigenstates which also carry group labels.


*The reason you might want to use the group is that some perturbation might lift the degeneracy so that states in different irreps have different eigenvalues (or at least some irreps have different eigenvalues, as there is to guarantee that all degeneracies are lifted).


*In the example you give, you start with
$\phi_0(x)\phi_1(y)$ and $\phi_0(y)\phi_1(x)$.  Imagine you were to add an interaction between the two particles - say some sort of potential which would be of the type $\kappa (x-y)^2$.  A perturbative treatment of this interaction would depend, to first order, on the average of $(x-y)^2$.  This perturbation is invariant under $S_2$, the group of permutation of the two coordinates.  This group has 2 1-dimensional representation and (no surprise) the irreps are spanned (up to normalization) by
\begin{align}
\Psi_{\pm}(x,y)&=\phi_0(x)\phi_1(y)\pm \phi_1(x)\phi_0(y)\, ,\\
&\sim e^{-\lambda(x^2+y^2)}(y\pm x)
\end{align}
The states $\Psi_{\pm}(x,y)$ are still degenerate under the original 2D h.o. oscillator Hamiltonian but, in a perturbative treatment,  the effect of interaction in $(x-y)^2$ would depend, to first order, on the average of $(x-y)^2$ evaluated with either
$\Psi_{\pm}(x,y)$, and this average is different for the two combinations.
Moreover, it is easy to check using parity that
\begin{align}
\int dx\,dy\, \Psi_-(x,y)(x-y)^2\Psi_+(x,y)=0\, .
\end{align}
Thus, by organizing your state according to irreps of $S_2$, which leaves your perturbation invariant, the degeneracy is lifted by the perturbation and the basis states for the irreps of $S_2$ are eigenstates of the Hamiltonian plus perturbation, at least to first order.
To understand the splitting of terms there is a massive paper (translated in English) by Hans Bethe

Bethe, Hans A. "Splitting of terms in crystals." Selected Works Of Hans A Bethe: (With Commentary). 1997. 1-72.

(originally Ann.Physics 3 p.133 (1929))
and the other canonical source is

Tinkham, M., 2003. Group theory and quantum mechanics. Courier Corporation.

A: Assume you know the degeneracy of the level is $n$. Then you can run the imaginary time evolution $n$ times, orthogonalizing to all previously-found wavefunctions each time, to obtain a complete basis for that energy level. In this basis, the action of each group element corresponds to a finite-dimensional matrix. Computing the elements of these matrices is straightforward. The problem thus reduces to simultaneous block-diagonalization of the matrices.
Whether this is tractable depends on specific facts about the group. The easiest case is a finitely-generated Abelian group. Then by Schur's lemma, the problem reduces to simultaneous diagonalization of a finite set of finite-dimensional matrices, for which there are known algorithms.
For a finite group, see https://arxiv.org/pdf/1901.05274.pdf.
For a finitely-generated nonabelian group, simultaneous block-diagonalization of the matrices corresponding to the generators will suffice. I'm not sure if there are good known algorithms for this.
If your group is a Lie group, then probably it suffices to block-diagonalize every element of a basis for the corresponding Lie algebra (possibly with some restrictions that the group be simple or otherwise "nice", and note that this is only practical for finite-dimensional Lie groups).
