Oh, sorry, I didn’t see your reply, yes it is almost the same solution as mine below.
Be sure to specify a range (i use it with slice) instead of just the two indices.
Sure, instead of one frequency index you could use a range, for example:
freq_limits = [8, 12]
freq_idx = [np.argmin(np.abs(freq - frequency)) for frequency in freq_limits]
freq_range = slice(freq_idx[0], freq_idx[1] + 1)
# average the frequency range:
psd_avg = psd[:, freq_range].mean(axis=1)
# plot
mne.viz.plot_topomap(psd_avg, epochs.info)
(BTW - for gamma you may prefer to use multitapering with some frequency smoothing instead of welch)