diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 5016803b1fb1e4c6482cde91fae555d4d3e28c9d..7f453179641e2ba16b30e3172ddd7853245a1d2f 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -62,9 +62,9 @@ INLINE static void convert_gpart_vel(const struct engine* e,
   ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 /**
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index d6b6f36ba23cee02690091611fb485ad54b47b4d..d47c96fbf32e1ee00346888aaf2e8afabc22abc3 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -101,9 +101,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 8f503230e1b2014c506e4c476f8597b5b77f85a1..3f2af41dc7f0cc8f60992a15a0f09f3c90f764fe 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -113,9 +113,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/hydro/GizmoMFM/hydro_io.h b/src/hydro/GizmoMFM/hydro_io.h
index 6ce5a23eeed66a1f6a7205d4683ecaeb91e05259..59d579f70cd4aedc728dbf42038eff78d4c507d5 100644
--- a/src/hydro/GizmoMFM/hydro_io.h
+++ b/src/hydro/GizmoMFM/hydro_io.h
@@ -158,9 +158,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/hydro/GizmoMFV/hydro_io.h b/src/hydro/GizmoMFV/hydro_io.h
index 7b1c651698bafad2740e179f77e4f5c3c7d1f677..92e4378f071cb71678929716be86588a3405f40e 100644
--- a/src/hydro/GizmoMFV/hydro_io.h
+++ b/src/hydro/GizmoMFV/hydro_io.h
@@ -158,9 +158,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 833103e7ab00738070cfe3ef8fc4ce4a2b8d6348..879255640fc1a1d6a06a666c80d3860c9c31ab64 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -127,9 +127,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/hydro/MinimalMultiMat/hydro_io.h b/src/hydro/MinimalMultiMat/hydro_io.h
index 5542b93ae43d88c93ed2e2ebddc1a6baef82a557..61301157e12b625f3ca1fb6cf5d5238d82403c9f 100644
--- a/src/hydro/MinimalMultiMat/hydro_io.h
+++ b/src/hydro/MinimalMultiMat/hydro_io.h
@@ -130,9 +130,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
index 4f168ed81961438b996f7973c9b04a8525ecf039..78967faec218f0efffbb624c4e8d25af214aad94 100644
--- a/src/hydro/PressureEnergy/hydro_io.h
+++ b/src/hydro/PressureEnergy/hydro_io.h
@@ -131,9 +131,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 /**
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index f2247c234aa46ad881266417a6147db87f807afe..8c11bf6e334e18b10217e90f6573a42e40880955 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -125,9 +125,9 @@ INLINE static void convert_part_vel(const struct engine* e,
   hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
 
   /* Conversion from internal units to peculiar velocities */
-  ret[0] *= cosmo->a2_inv;
-  ret[1] *= cosmo->a2_inv;
-  ret[2] *= cosmo->a2_inv;
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
 }
 
 INLINE static void convert_part_potential(const struct engine* e,
diff --git a/src/space.c b/src/space.c
index aaea061e8ecc865785d2453f6ea7c9955b2e7185..3efb1080c37bb8057ec6712a67370b8ac3a2b39e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2354,7 +2354,7 @@ void space_first_init_parts_mapper(void *restrict map_data, int count,
   const struct cosmology *cosmo = s->e->cosmology;
   const struct phys_const *phys_const = s->e->physical_constants;
   const struct unit_system *us = s->e->internal_units;
-  const float a_factor_vel = cosmo->a * cosmo->a;
+  const float a_factor_vel = cosmo->a;
 
   const struct hydro_props *hydro_props = s->e->hydro_properties;
   const float u_init = hydro_props->initial_internal_energy;
@@ -2429,7 +2429,7 @@ void space_first_init_gparts_mapper(void *restrict map_data, int count,
   const struct space *restrict s = (struct space *)extra_data;
 
   const struct cosmology *cosmo = s->e->cosmology;
-  const float a_factor_vel = cosmo->a * cosmo->a;
+  const float a_factor_vel = cosmo->a;
   const struct gravity_props *grav_props = s->e->gravity_properties;
 
   for (int k = 0; k < count; k++) {
@@ -2486,7 +2486,7 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count,
 #endif
 
   const struct cosmology *cosmo = s->e->cosmology;
-  const float a_factor_vel = cosmo->a * cosmo->a;
+  const float a_factor_vel = cosmo->a;
 
   for (int k = 0; k < count; k++) {
     /* Convert velocities to internal units */
diff --git a/theory/Cosmology/coordinates.tex b/theory/Cosmology/coordinates.tex
index e3a22eaae025e6911e8aca92d6d29bf5fa82bf21..bc593606026217f345e2f77d90cee3f6632b3cef 100644
--- a/theory/Cosmology/coordinates.tex
+++ b/theory/Cosmology/coordinates.tex
@@ -40,19 +40,29 @@ $\Psi \equiv \frac{1}{2}a\dot{a}\mathbf{r}_i^2$ and obtain
   -\frac{\phi'}{a},\\
   \phi' &= a\phi + \frac{1}{2}a^2\ddot{a}\mathbf{r}_i'^2.\nonumber
 \end{align}
-Finally, we introduce the velocities $\mathbf{v}' \equiv
-a^2\dot{\mathbf{r}'}$ that are used internally by the code. Note that these
-velocities \emph{do not} have a physical interpretation. We caution that they
-are not the peculiar velocities, nor the Hubble flow, nor the total
-velocities\footnote{One additional inconvenience of our choice of
+Finally, we introduce the velocities
+$\mathbf{v}' \equiv a^2\dot{\mathbf{r}'}$ that are used internally by
+the code. Note that these velocities \emph{do not} have a physical
+interpretation. We caution that they are not the peculiar velocities
+($\mathbf{v}_{\rm p} \equiv a\dot{\mathbf{r}'} =
+\frac{1}{a}\mathbf{v}'$), nor the Hubble flow
+($\mathbf{v}_{\rm H} \equiv \dot{a}\mathbf{r}'$), nor the total
+velocities
+($\mathbf{v}_{\rm tot} \equiv \mathbf{v}_{\rm p} + \mathbf{v}_{\rm H}
+= \dot{a}\mathbf{r}' + \frac{1}{a}\mathbf{v}'$) and also differ from
+the convention used in \gadget snapshots
+($\sqrt{a} \dot{\mathbf{r}'}$) and other related simulation
+codes\footnote{One additional inconvenience of our choice of
   generalised coordinates is that our velocities $\mathbf{v}'$ and
   sound-speed $c'$ do not have the same dependencies on the
   scale-factor. The signal velocity entering the time-step calculation
-  will hence read $v_{\rm sig} = a\dot{\mathbf{r}'} + c = \frac{1}{a} \left(
+  will hence read
+  $v_{\rm sig} = a\dot{\mathbf{r}'} + c = \frac{1}{a} \left(
     |\mathbf{v}'| + a^{(5 - 3\gamma)/2}c'\right)$.}.
-This choice implies that $\dot{v}' = a \ddot{r}$. Using the SPH
-definition of density, $\rho_i' =
-\sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
+% This choice implies that $\dot{v}' = a \ddot{r}$.
+
+Using the SPH definition of density,
+$\rho_i' = \sum_jm_jW(\mathbf{r}_{j}'-\mathbf{r}_{i}',h_i') =
 \sum_jm_jW_{ij}'(h_i')$, we can follow \cite{Price2012} and apply the
 Euler-Lagrange equations to write
 \begin{alignat}{3}