Here is an example that works with two bands. Up to you to generalize to four bands.
WITH raster AS (
SELECT ST_AddBand(ST_AddBand(ST_AddBand(ST_SetValue(ST_AddBand(ST_MakeEmptyRaster(2, 2, 0, 0 ,1), '8BUI'), 1, 1, 2), '8BUI'::text, 2), '8BUI'::text, 2), '8BUI'::text, 4) rast
), bands AS (
SELECT 1 band, ST_Band(rast, 1) rast FROM raster
SELECT 2 band, ST_Band(rast, 2) rast FROM raster
, max AS (
SELECT ST_Union(rast, 'MAX') rast FROM bands
), index1 AS (
SELECT ST_MapAlgebra(b1.rast, max.rast, 'CASE WHEN [rast1] = [rast2] THEN 1 ELSE NULL END') rast
FROM (SELECT rast FROM bands WHERE band = 1) b1,
), index2 AS (
SELECT ST_MapAlgebra(b2.rast, max.rast, 'CASE WHEN [rast1] = [rast2] THEN 2 ELSE NULL END') rast
FROM (SELECT rast FROM bands WHERE band = 2) b2,
), indexbands AS (
SELECT rast FROM index1
SELECT rast FROm index2
), indexmax AS (
SELECT ST_Union(rast, 'MIN') rast FROM indexbands
) SELECT (ST_PixelAsPolygons(rast)).* FROM indexmax;
The idea is to create a temporary raster being the union(max) and to compare each band to this raster to build indexes rasters and then to union again these indexes rasters to find the first matching index.
Might be slow on big rasters given the number of operations...
I agree this is overwhelmingly complex for a task that seems easy. We would need a ST_Union() taking a user custom aggregate function. That's how ST_Union() was designed at the very beginning but it's not the way it was implemented. The Pl/pgSQL version still exists if you want it. Might need a little refresh.