You are right, the proper time is ill defined for massless particles or null vectors if you prefer. However you can still consistently define a parameter and take the derivative with respect to it to obtain the tangent vector for the case ds^2=0.

It should be discussed in any GR book, or take a look in this discussion:

https://is.gd/InKLaY

The proper time isn't ill-defined-it's zero, for a massless particle. It's a Lorentz invariant quantity, so there's no issue of it being ill-defined.

The 3-velocity of a relativistic particle, however, is ill-defined, since it isn' t invariant under Lorentz transformations; that's why it doesn't make sense discussing its value.

So the answer to the question in the text is No: While the norm of the photon's 4-velocity is equal to c in any inertial frame (that's why a massless particle doesn't have a rest frame), its direction depends on the Lorentz transformation that relates the frames. So the 3-velocity depends on the frame, as might be expected.

The answer to the main question goes as follows:

Massless particles ARE described by 4-vectors, just like massive particles are.

The equations of motion for both are d2xμ/dλ2=0, where λ parametrizes the worldline. The massive particle is distinguished from the massive particle by the constraints.

The constraint for a massive particle is p2=m2, where p is the energy-momentum 4-vector; while it's p2=0 for a massless particle. This constraint, that is invariant under Lorentz transformations, in turn, can be imposed in a fully covariant way by realizing that it is equivalent to local reparametrization invariance of the trajectory-something that's, unfortunately, not as well-known as it should be.

The point is that there's a gauge invariance involved, that's the same for the massive particle, too; just the constraints are different. GR doesn't have anything to do with it (though, for some reason, it's often discussed in such texts). The relativistic particle is a constrained system-that's the point.

This is discussed here: http://fma.if.usp.br/~amsilva/Livros/Zwiebach/chapter5.pdf for instance.

though it stops too soon.

The idea is to write the Lagrangian of the relativistic particle, -m[(dx/dλ)^2](1/2)

as (-1/2e(λ))(dx/dλ)2-(m2e(λ)/2), by introducing the auxiliary field e(λ). The two Lagrangians are equivalent, since the field e(λ) isn't dynamical, since the Lagrangian doesn't contain its derivatives with respect to the parameter λ- it imposes the constraint about the mass. The spurious issue of the first Lagrangian not having a useful limit, when m->0, is thus resolved by using the second Lagrangian, that does. It's the same thing as resolving the spurious singularities of coordinate systems.

In this way, the massless particle is described by the Lagrangian (1/2e(λ))(dx/dλ)2 in a redundant way, in terms of its coordinates in spacetime and the auxiliary field-called the einbein-and there's a gauge choice that must be made, in order for the initial value problem to admit a unique solution and one such gauge choice is setting e(λ)=1. This is discussed, for instance, here: https://inspirehep.net/literature/4176

So, Yes, the conjugate momenta are given by pμ=dxμ/dλ-the point is that they're subject to a constraint and λ isn't a physical quantity, since it's not invariant under reparametrization of the trajectory, that's a symmetry of the theory.

The whole GR and SR is not presenting a universe with billions of galaxies.

Article How the Quantum Universe Works