How to correct phase of Diagonal Volterra kernels obtained with exponential sweep sine?
A potential pitfall in non-linear modeling is using different signal levels. The sweep test needs to be carried out at a sufficiently high signal level to cover the maximum possible real-world input levels, but not much higher as this will require more kernels for the same effective quality. Usually real-world signals to the system are much of the time at lower level than their peak level, with the 1st order linear kernel dominating.

To me your kernels look just fine. I don't know how common it is to see that a higher-order kernel (3rd order in your plots) is more prominent than the 1st order kernel. The saw-like phase comes from the large initial delay, which gives a steep linear dependency of phase on frequency. Wrapped phase of $\pi$ is the same as a wrapped phase of $-\pi$, and you are seeing wrapped phase, which gives the sawtooth pattern. If you remove the initial delay in the inspection of the results, you should be able to unwrap the phase more easily for inspection.

Most systems you would model by the technique will have right-sided kernels. The systems are causal and cannot store energy for very long after they have received an impulse as input. That would require a delay line. Instead, much of the energy passes through almost instantly as the initial spike, and the residual energy is released as the tail.

I suggest you test with a system you have complete control over, for example a memoryless non-linearity that is the sum of the $N$ first Chebyshev polynomials of the first kind, one for each of the $N$ harmonics you want to generate, including the 0th and 1st harmonics in the count. Then any phase shifts due to errors in the algorithms/software should be clearly visible. If the algorithms work correctly, each harmonic impulse response should be a delayed sinc, and the the coefficients of the memoryless nonlinearity polynomial should be resolved as the amplitudes of the scaled impulse-like kernels.

What I have done in my own implementation is to generate the sweep and its fade in and fade out in the frequency domain by programming the timing into the group delay, and amplitude into the magnitudes, rather than using a direct time domain equation for the sweep. Either way should work and will have its own problems. The frequency-domain approach gives amplitude oscillations near the fades. In the frequency domain it is easy to correct the non-whiteness of the sweep in the inverse filter. (I can't publish my implementation in detail at the moment.)

If there are phase problems, I'd try $\cos$ instead of $\sin$ in the time-domain sweep formula, and fade the sweep in and out. The fades will take care of eliminating transients at the ends. The initial and final phases of the sweeped sinusoid could be made such that the value of the sweep is zero both in the beginning and in the end, but with the fades in and out that is no longer needed for transient reduction, because the fades gently bring the value to zero in any case.

I don't see a problem with the sweep going all the way to the Nyquist frequency, if the non-linearity is analog and the analog-to-digital converter (ADC) has an anti-alias low-pass filter at the Nyquist frequency. This way the harmonic impulse responses and consequently the Hammerstein kernels have only non-aliased frequencies up to Nyquist frequency. It is sufficient information to model the non-linear system in the sub- Nyquist frequency band if the system is compatible with the Hammerstein model. Only the implementation has to be careful not to introduce aliasing. In any case, you'd typically want the linear kernel to have frequency response information up to or at least near the Nyquist frequency, which dictates the sweep end frequency.