logistically, I would suggest collecting the raster names from each of the two folders first, make the pairing and do the math
flder0 = ['a_x01', 'a_x02', 'a_x03', 'a_x04']
flder1 = ['b_c_y01', 'b_c_y02', 'b_yc_03', 'b_yc_04']
xs = [i.rsplit('_')[1] for i in flder0]
ys = [i.rsplit('_')[-1] for i in flder1]
xsys = list(zip(xs, ys))
names = ["{}{}".format(*i) for i in xsys]
for i, nme in enumerate(names):
eq = "{} = {} * {}".format(*[nme, flder0[i], flder1[i]])
print(eq)
x01y01 = a_x01 * b_c_y01
x02y02 = a_x02 * b_c_y02
x0303 = a_x03 * b_yc_03
x0404 = a_x04 * b_yc_04
Lines 1 and 2, would be a arcpy.ListRasters on each folder
Lines 3 and 4 would simply be splitting the names into their bits to
line 5 is just part of the process of formatting the output names
When assembled, lines 8 and on are a 'fake' run of the process.