Blender 4.0+
As already mentioned in the comments, the solution for current Blender versions is as follows:

The Corners of Edge node returns (in the case that the edge is connected to two faces) the index of the face corner of the adjacent faces selected with Sort Index:
(Image from the official documentation)
The Face of Corner node uses this index to select the respective face.
In this example, I save the index of both adjacent faces in X and Y of a vector for convenient further processing.
Blender 3.6 or below
As I see it, this is only a question of the correct interpolation between the attribute domains.
If you want to select the two faces connected with an edge, then you can solve this as follows:

- Select the desired edge and capture this selection with
Capture Attribute as a float value in the domain Edge.
- Then interpolate from the domain Face Corner to Faces (the previously stored selection in Edge is transferred to Face Corner, and then to the domain Face) and store this value also in the domain Face.
Finally, you can use this selection (only two faces are selected) to fetch the corresponding values from these faces.
You can do this, for example, by separating the faces with this selection and Separate Geometry, and then fetching them from there with Sample Index (index $0$ and $1$):

But...
The solutions shown above are funny when it comes to finding out only one edge and its adjacent faces.
But what if you want to get the adjacent faces for all edges of a mesh, i.e. do exactly what you would expect from a non-existing node Faces of Edge?
Then the task becomes a lot more fun and the solution might look something like this:

The trick here is that the node Edge Angle fortunately gives you back the signed angle of an edge. This value is somewhere between $-180°$ and $180°$.
An angle of $0°$ means that the two adjacent faces are planar, an angle greater than $0°$ means that the two faces are convex, and an angle less than $0°$ means that the faces are concave.
Now by mapping this angle into a range of $90° - 180°$ (without Clamp), I get the bisected angle between the two faces.
If I now rotate the normal of an edge (which already represents the interpolated direction of the normals of both faces) with the direction vector of the edge as the Axis with this angle, this resulting vector is always exactly perpendicular to the normals of one of the two faces.
And now I simply sample the closest face and use exactly this vector scaled with a very small value like $0.000001$ and added to the position of the respective edge as Sample Position. Since an edge has at least one face (hopefully), I repeat this with the inverted angle as well.
In this way I can relatively reliably determine per edge the indices of both adjacent faces (and consequently also the normals of these faces).
Finally, I save these two indices into a 2D Vector called "faces_of_edge" to make them easily reusable.

...of course this solution works only for meshes whose edges have a maximum of two faces.
Ah, and by the way: This solution can theoretically be adapted from version 3.1 (not tested), which is also cool.