GDAL supports a mask layer, but it also supports a "no data" colour which I tend to use because it cuts the size of the final image down by 25%. So the steps I would use with your data, after cutting off the collar is:
gdal_translate -expand rgb -gcp <pixel/coordinate pairs> -co TILED=YES -a_srs <EPSG code for the map's projection> -a_nodata 255 255 255 original.tif translated.tif
Note, I've put in Ground Control Points (GCPs) and a projection which you'll have to determine empirically. Also note I've assigned pure white to the nodata value.
These images are still not properly georeferenced however, as they still only have GCPs rather than proper coordinates. So you will need to run gdalwarp:
gdalwarp -t_srs EPSG:4326 -co TILED=YES -srcnodata 255 255 255 -dstnodata 255 255 255 translated.tif warped.tif
Note again the explicit setting of the nodata values. This is probably overkill, but it doesn't hurt, and I've had gdalwarp ignore what I thought was implicit nodata vales.
Next, you need to run gdalbuildvrt like you have it, but I usually explicitly specify the nodata values again, to be on the safe side:
gdalbuildvrt -srcnodata 255 255 255 -vrtnodata 255 255 255 temp.vrt <warped files>.tif
Edit
This may be closer to what you want:
gdalbuildvrt -srcnodata 255 255 255 -hidenodata temp.vrt <warped files>.tif
Finally, you can call gdal2tiles.py as you have it.
Best Answer
Sounds like you'd want to use the image mosaic plugin: http://docs.geoserver.org/stable/en/user/data/raster/imagemosaic/index.html
The configuration of the mosaic can be scary, but if all you want is to merge all the images together visually, then put them in a directory and just point the GeoServer image mosaic at it with no other setups, it should self-configure.
Then you might find it's slow.... the files can be optimized one by one using command line tools, and you might decide to start showing the data only past a certain scale with a SLD scale dependency.