Geographic lat/lon coordinates are already angular measures, so you can't compute a simple arctangent with the coordinates based on {delta-y,delta-x} as if they were Cartesian values.

Using a port of the US Geodetic Survey FORTRAN code to solve the inverse (aka reverse) geodetic problem, and using the coordinate values in your Java code, I got the following:

p1                      |p2                       |gcbearing
POINT ( 2.3582 48.86124)|POINT ( 2.35922 48.86215)|36.4845280
POINT ( 2.3582 48.86124)|POINT ( 2.35956 48.86054)|127.9567007

The difference between the bearings is 91.472 degrees.

You cannot rely on a Web Mercator canvas for whether the intersection "looks like" a right angle because of the significant distortion in Web Mercator, even at only 49N.

I don't have any experience in JTS, but I would expect it has some geodetic functions, in addition to Cartesian ones.

Geographic (lat/lon) coordinates are already angular measures, so you can't compute a simple arctangent with the coordinates based on atan2(delta_y,delta_x) as if they were Cartesian values.

Using a port of the US Geodetic Survey FORTRAN code to solve the inverse (aka reverse) geodetic problem, and using the coordinate values in your Java code (assuming WGS84 datum), I got the following:

p1                     |p2                      |gcbearing
POINT (2.3582 48.86124)|POINT (2.35922 48.86215)|36.4845
POINT (2.3582 48.86124)|POINT (2.35956 48.86054)|127.9567

The difference between the bearings is 91.4722 degrees.

Reprojecting the points to WGS84 UTM Zone 31N and using Cartesian bearing, I got:

p1                            |p2                            |bearing
POINT (452926.994 5412229.237)|POINT (453002.658 5412329.763)|36.9681
POINT (452926.994 5412229.237)|POINT (453026.088 5412150.583)|128.4401

The difference here is 91.4720 degrees (since Paris is on the left side of the UTM cylinder the bearings are slightly different, as well)

Finally, reprojecting into a custom WGS84 datum flavor of the Europe Albers Equal Area Conic projection (central meridian 10.0, standard parallels 43.0,62.0), I got:

p1                            |p2                             |bearing
POINT ( -553711.962 2098386.8)|POINT (-553627.667 2098480.823)|41.8774
POINT ( -553711.962 2098386.8)|POINT (-553621.972 2098298.257)|134.5356

with a difference of 92.6582 degrees.

You cannot rely on a Web Mercator canvas for whether the intersection "looks like" a right angle because of the significant Y-axis distortion in Web Mercator, even at only 49N.

I don't have any experience in JTS, but it apparently only has Cartesian (not geodetic) functions, so you'd need to reproject to compute bearing changes (based on the above I'd recommend using UTM for relatively small areas).