Solving this system of ODEs does not seem feasible, either symbolically or numerically. However, it is feasible to provide approximate solutions for t < 100.

For convenience, simultaneously solve the first three ODEs numerically.

ode1 = x''[t] + 3  x'[t] Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t]  
  Sqrt[2/3]])^2/4] + Sqrt[3/2] Exp[-x[t] Sqrt[2/3]]  (1 - Exp[-x[t]  
  Sqrt[2/3]]) == 0;
ode1IC = {x'[0] == -0.008226306418212731, x[0] == 5.630991866033891};

ode2 = a'[t]/a[t] == Sqrt[(x'[t])^2/6 + (1 - Exp[-x[t] Sqrt[2/3]])^2/4];
ode2IC = a[0] == 1;

ode3 = \[Tau]'[t] == 1/a[t];
ode3IC = \[Tau][149.4517772937791] == 0;

sol123 = NDSolve[{ode1, ode2, ode3, ode1IC, ode2IC, ode3IC}, 
  {x[t], \[Tau][t]}, {t, 0, 500}] // Flatten

Plot[Evaluate[{x[t], \[Tau][t]} /. sol123], {t, 0, 500}, PlotRange -> All, 
  AxesLabel -> {t, "x,\[Tau]"}, LabelStyle -> {10, Bold, Black}]

The small oscillations in x gradually die away with increasing t. a is not included, because it is not needed in subsequent calculations. Next, compute and plot f, which appears in ode4.

f = ((7.41193*10^6)^2 - ((3/2 + 1/2*((Sqrt[3/2]  Exp[-x[t] Sqrt[2/3]] 
  (1 - Exp[-x[t] Sqrt[2/3]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))^2 + 
  1/2 ((Exp[-Sqrt[8/3] x[t]] (1 -  Exp[Sqrt[2/3] x[t]] 
  (1 - Exp[-Sqrt[2/3] x[t]]))/(3/4 (1 - Exp[-x[t] Sqrt[2/3]])^2))))^2 - 1/4)
  /(\[Tau][t])^2);

Plot[f /. sol123, {t, 0, 145}, AxesLabel -> {t, "x,\[Tau]"}, 
  LabelStyle -> {10, Bold, Black}]

For larger values of t, f is strictly negative and oscillates between -10^-17 and -10^-30. The analysis below focuses on t < 100. As seen in the plot just above, f is approximately constant for t < 20, and ode4 produces solutions oscillating at a frequency Sqrt[f /. sol123 /. t -> 0], i.e., 7.41193*10^6. It is for this reason that NDSolve takes seemingly forever to solve ode4. On the other hand, ode4 can be solved approximately using the WKB approximation, if desired, for t < 100, except near f = 0. Here, an approximate symbolic solution exists:

FindRoot[f /. sol123, {t, 36}][[1, 2]]
(* 36.7625 *)
Series[f /. sol123, {t, %, 1}] // Normal
(* -0.015625 - 7.44662*10^12 (-36.7625 + t) *)

The constant term is neglibile and can be ignored, allowing ode4 to be solved in terms of Airy functions.

as = DSolveValue[u''[t] + %[[2]] u[t] == 0, u[t], t] /. {C[1] -> 1, C[2] -> 0}
(* AiryAi[2.62237*10^-9 (-2.73757*10^14 + 7.44662*10^12 t)] *)

(The second solution, deleted here, grows exponentially for large t.)

Plot[as, {t, 36.7625 - .002, 36.7625 + .002}, AxesLabel -> {t, u}, 
  LabelStyle -> {10, Bold, Black}]

The WKB approximation, which I shall not take the time to apply here in detail, is described in https://en.wikipedia.org/wiki/WKB_approximation. From it, the amplitude of the oscillations varies as f^(-1/4) for t < 36.5. Evaluating f at t = 0 and t = 36.5 indicates that the amplitude of oscillations in u grows by about a factor of

((f /. sol123 /. t -> 36.5)/(f /. sol123 /. t -> 0))^(-1/4)
(* 2.28105 *)