Found the problem. My noise budget code was wrong. So after I fixed it, the TE noise in substrate result from COMSOL agrees pretty well with the analytical result (within 20%).
The result from COMSOL is plotted in dashedblack line. The result from Cerdonio is plotted in dashed pink. Since my simulation uses the adiabatic assumption (used in BGV and Liu&Thorne paper), the results agree at high frequency. So I think the calculation is correct. I'll check some options (changing spot size, changing material) to see if TE noise can be made lower for AlAs/GaAs samples.
I attached my COMSOL file below. This is done in 3D model. It could have been done in 2D axis symmetric setting, but I used 3D for spacer sagging before, so I just used the same geometry I had.
Quote: 
I realized that the mesh size was too large, even with the finest mesh for default setting. So I reduced the mesh size around the beam area and the results got closer to the analytical prediction. It is still a factor of 2 below the prediction. I'll see if I can hunt down this problem . I think it will be a good idea to verify my model by using my model to calculate Brownian noise and comparing with the result reported by Braunschwig group.
When I defined mesh size in COMSOL, I used the predefined value provided by COMSOL. The finest mesh has maximum element size ~500 um, and minimum ~5um. Since the beam size is ~ 180 um, the maximum element size should be ~10 um. So I changed the values around, defined new area for smaller mesh until the results did not change much. I ran the simulation a few time to make sure that the solution converges. Right now my substrate has 3 regions
 a cylinder at the center where the beam hits the mirror, Radius = 2x180um, depth=2x180um , Min/Max mesh size = 12.6/27 um
 an outer region, radius = 5x180um, depth = 5x180 um, min/max mesh size = 10/50 um
 The rest of the substrate, fine mesh
I tried to change the mesh size/boundary size a bit to get the result accurate enough without taking too much time. The TE estimation still a factor of 2 below the analytical estimate.

