These are the likelihoods after each stage:
Phyml BioNJ -43819.695026
Raxml topology search -43686.158939
Phyml branch length optimization -43673.144059
The main thing is that the branch lengths are very small because the sequences are long (19kb). This means they are arbitrarily interacting with a lower bound in the optimization algorithm. The branches that are affected are those with the maximum likelihood value of zero.
Here is a plot of the equivalent terminal branch lengths:
The difference in a branch length of 1E-6 and 1E-8 does actually matter here because there are some large polytomies that result in a 'ladderised' affect.
Possibly there is an option to set the lower bound on branch lengths? Interestingly Phyml is clearly going hard against its lower bound of 1E-8 but Raxml has a lower bound of 1.47E-6 which seems unlikely to be a selected constant so I suspect this is the result of some other stopping criterion?
I suspect that a general solution would be to introduce a constant scaling factor for the branch lengths between the optimization search algorithm and the likelihood evaluation. A suggestion would be to use the ratio of the sequence length to the number of unique site patterns.