Вот этот свободный 4-вектор и нужно преобразовывать в другие СО по обычному правилу преобразования ковариантных векторов
Э, нет! Это незаконный трюк. Ваш "свободный вектор" определён, вообще говоря, не в точке, а на некой гиперповерхности

. А для всей гиперповерхности разом Вы компоненты матрицы Якоби преобразования не посчитаете, ибо они будут разными в разных точках. Ладно, допустим, что в данном конкретном случае Вам повезло: Мировая линия частицы пересекает гиперповерхность

только в одной точке. А если у Вас две частицы, то что будете делать?
А с чего Вы взяли, что интегральная энергия является константой? Нет. В общем случае, она является 0-компонентой 4-векторного переменного поля (4-импульса). Просто в случае сохраняющегося 4-импульса - это поле ковариантно постоянно (т.е. все его ковариантные производные равны нулю). Так что в случае с двумя (и более) частицами ничего страшного не случится.
Возникает вопрос: каким образом Вам удалось получить правильный ответ?
Секрет прост. Вот эта формула для энергии правильная:

,
а вот эта неправильная:

.
Не совсем так. Давайте повторим рассуждения Фока и посмотрим в чем может быть "зарыта собака".
Пусть имеется физическая система с тензором энергии-импульса

в пространстве-времени Минковского (в произвольной криволинейной системе координат). Возьмем некоторое несингулярное (это уже от меня лично - не хочу заморачиваться случаем изотропного вектора) ковариантное векторное поле

и образуем свертку с ТЭИ

. Найдем ковариантную дивергенцию этой свертки

Проинтегрируем эту ковариантную дивергенцию по пространственноподобной гиперповерхности (3-объему), причем предполагаем что координата с индексом 0 как всегда времениподобна, а остальные пространственноподобны

Рассмотрим более подробно левую часть этого выражения. С учетом определения ковариантной дивергенции векторного поля получим

Выделим в этом интеграле слагаемое с производной по времениподобной координате

где

- пробегает 1,2,3. Во втором интеграле используем теорему Гаусса

Из всех этих формул следует

Таким образом, если система сосредоточена в конечной области (вне которой

), ТЭИ сохраняется (

) и выполнены уравнения Киллинга для вектора

, то мы имеем сохранение интеграла слева со временем

Далее в игру вступает тот факт, на который указывали Вы epros: вследствие линейности уравнений Киллинга для их решений выполнен принцип суперпозиции - если мы имеем два независимых решений

то решением является также и их произвольная линейная комбинация

, где

- произвольные константы. Т.е. какую бы мы комбинацию не подставили в последний интеграл, он всегда будет равен нулю - фактически сохранение интеграла со временем
не зависит от выбора Киллингова вектора. Это (а также форма интеграла в декартовых координатах) позволяет определить интеграл как сохраняющийся 4-импульс.
Отсюда же следует некорректность Ваших примеров: Вы использовали физические системы с несохраняющимся ТЭИ. В этих случаях, как мне кажется, определение Фока не работает.
...Мне более интересно, как Вы ухитрились получить нужный Вам ответ для интеграла

. Похоже, что при преобразовании вектора

Вы просто забыли про лоренцево сокращение размеров частицы в

раз, так что произошла "взаимная компенсация ошибок".
Я просто не довел вычисления до конца. Там получается плотность массы со сложным аргументом и интеграл от нее не равен

.
...Надеюсь, что теперь Вы уже не скажете:
Я просто знаю правильный ответ, а потому или поверьте мне на слово или проверьте.

Сейчас скажу это еще более уверенно. Благодаря Вам я разобрался с некоторыми темными для меня моментами.