diff --git a/sgkit/stats/pc_relate.py b/sgkit/stats/pc_relate.py
index eb06a2b99..f19bdeb34 100644
--- a/sgkit/stats/pc_relate.py
+++ b/sgkit/stats/pc_relate.py
@@ -107,6 +107,7 @@ def pc_relate(
 
     :data:`sgkit.variables.pc_relate_phi_spec`: (S,S) ArrayLike
     pairwise recent kinship coefficient matrix as float in [-0.5, 0.5].
+    The dimensions are named ``samples_0`` and ``samples_1``.
 
     References
     ----------
@@ -173,5 +174,7 @@ def pc_relate(
     phi = gramian(centered_af) / gramian(stddev)
     # NOTE: phi is of shape (S x S), S = num samples
     assert phi.shape == (call_g.shape[1],) * 2
-    new_ds = create_dataset({variables.pc_relate_phi: (("sample_x", "sample_y"), phi)})
+    new_ds = create_dataset(
+        {variables.pc_relate_phi: (("samples_0", "samples_1"), phi)}
+    )
     return conditional_merge_datasets(ds, new_ds, merge)
diff --git a/sgkit/tests/test_pc_relate.py b/sgkit/tests/test_pc_relate.py
index 3817325ae..1977f971c 100644
--- a/sgkit/tests/test_pc_relate.py
+++ b/sgkit/tests/test_pc_relate.py
@@ -122,7 +122,7 @@ def test_pc_relate__identical_sample_should_be_05() -> None:
     g.call_genotype.loc[dict(samples=8)] = g.call_genotype.isel(samples=0)
     phi = pc_relate(g)
     assert phi.pc_relate_phi.shape == (n_samples, n_samples)
-    assert np.allclose(phi.pc_relate_phi.isel(sample_x=8, sample_y=0), 0.5, atol=0.1)
+    assert np.allclose(phi.pc_relate_phi.isel(samples_0=8, samples_1=0), 0.5, atol=0.1)
 
 
 def test_pc_relate__parent_child_relationship() -> None:
@@ -164,7 +164,7 @@ def simulate_new_generation(ds: xr.Dataset) -> xr.Dataset:
         ds["pc_relate_phi"]
         .to_series()
         .reset_index()
-        .pipe(lambda df: df.loc[df.sample_x >= df.sample_y]["pc_relate_phi"])
+        .pipe(lambda df: df.loc[df.samples_0 >= df.samples_1]["pc_relate_phi"])
         .pipe(
             pd.cut,
             bins=[p for phi in [0, 0.25, 0.5] for p in [phi - 0.1, phi + 0.1]],