Hi Dimitris,
You can use the script below.
Consider that the script runs on windows only and is for only one zonal statistic type (no eg Min_Max), if you would like for more, a small adjustment is needed.
I am new in python, so the code is unconventional... but the result is correct.
The script loop through a list of band defined by a multiband raster and executed zonal statistics table, join only the field with the statistics from the table to the specified polygon layer, add a new field with the band name and values equal to the joined field and then delete the joined field.
The polygon is saved to the desired location when the loop run as many times as the band number is.
Let me know if this works.
P.S. Raki is expected when I will be at Heraklion