The main problem is in the image warping parameters, they are not giving you the right gridded image.
About the projection itself, you are using the Lambert Conical projection in the right way. You have just not information about the datum transformation. You can't assume that geographic coordinates mapped in 1944 to the Bessel ellipsoid don't need a datum transformation to be transformed to the WGS84 ellipsoid.
About warping the image from the geographic coordinates of the corners with a linear transformation seems to me a wrong process, beyond that the result looks better than the previous one.
I have loaded the original .jpg file and digitized your Ground Control Points:
They seems good. But when you warp the image, you are forcing computation to first degree polynomials (-order 1), which is something like performing a linear transformation in the image (the overabundance of points is used only in the compensation). This is the result with a rectangular grid superimposed:
Since eight GCPs are enough to a second order polynomial, this is the result with -order 2 parameter (or without it, because it will be the default for 8 GCPs):
Seems better, but the scanner has its deformations. This is the result forcing the Thin Plate Spline transformation (-tps):
That's the best that we can do with those GCPs. Let's see how it is georeferenced:
Those are the differences between the Porta Sempione mapped in 1944 and today.
The Monte Mario datum based on the Hayford (International 1909) ellipsoid was adopted in 1940, but this map explicitly says that it is projected from the Bessel ellispsoid.
From Geodetic datums of the Italian cadastral systems:
In the first decades of the twentieth century, the Italian Institute of Military Geography (Istituto Geografico Militare; I.G.M) developped four geodetic network, all on the Bessel 1841 ellipsoid; the Genova 1902, the old Monte Mario, the Castanea delle Furie 1910 and the Guardia Vecchia datums, for northern, central and southern Italy and Sardinia, respectively ( Mori, 1922 ).
The authors of the same paper derived a set of unified transformation parameters, from these datums to WGS84:
dX= +656.5 m; dY= +138.2 m; dZ= +506.5 m; rotX= –5.187 arc sec; rotY= +2.540 arc sec; rotZ= –5.256 arc sec; scale factor= –12.61 ppm (according to ‘coordinate frame rotation’; average horizontal error, except Sicilia and Southern Italy: 2.5 m; maximum horizontal error: 5.6 m.)
We need to invert the rotation parameters because the convention for the towgs84 parameter is "position vector" instead of "coordinate frame":
gdalwarp -r near -tps -co COMPRESS=PACKBITS -t_srs "+proj=lcc +lat_0=45.90 +lon_0=14 +lat_1=45.90 +x_0=800000 +y_0=601000 +k_0=0.998992911 +ellps=bessel +towgs84=656.5,138.2,506.5,5.187,-2.540,5.256,-12.61 +units=m +no_defs" ./output_gcp.tif ./output_warped.tif
You will check the accuracy for the whole image, it is really very good.