diff --git a/doc/changes/dev/13943.bugfix.rst b/doc/changes/dev/13943.bugfix.rst
new file mode 100644
index 00000000000..9ac1b6c21d3
--- /dev/null
+++ b/doc/changes/dev/13943.bugfix.rst
@@ -0,0 +1 @@
+Fix channel ordering bug in :func:`mne.viz.plot_compare_evokeds` when ``axes="topo"``, by `Lifeng Qiu Lin`_.
\ No newline at end of file
diff --git a/mne/evoked.py b/mne/evoked.py
index 17048bf0193..c1d62a1df52 100644
--- a/mne/evoked.py
+++ b/mne/evoked.py
@@ -1609,7 +1609,7 @@ def _get_evoked_node(fname):
     return evoked_node
 
 
-def _check_evokeds_ch_names_times(all_evoked):
+def _check_evokeds_ch_names_times(all_evoked, inplace=False):
     evoked = all_evoked[0]
     ch_names = evoked.ch_names
     for ii, ev in enumerate(all_evoked[1:]):
@@ -1618,7 +1618,8 @@ def _check_evokeds_ch_names_times(all_evoked):
                 raise ValueError(f"{evoked} and {ev} do not contain the same channels.")
             else:
                 warn("Order of channels differs, reordering channels ...")
-                ev = ev.copy()
+                if not inplace:
+                    ev = ev.copy()
                 ev.reorder_channels(ch_names)
                 all_evoked[ii + 1] = ev
         if not np.max(np.abs(ev.times - evoked.times)) < 1e-7:
diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index a62d2379f03..63aca378b59 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -2858,9 +2858,7 @@ def plot_compare_evokeds(
         for evk in evokeds[cond]:
             _validate_type(evk, Evoked, "All evokeds entries ", "Evoked")
     # ensure same channels and times across all evokeds
-    all_evoked = sum(evokeds.values(), [])
-    _check_evokeds_ch_names_times(all_evoked)
-    del all_evoked
+    _check_evokeds_ch_names_times(sum(evokeds.values(), []), inplace=True)
 
     # get some representative info
     conditions = list(evokeds)
@@ -3019,6 +3017,7 @@ def plot_compare_evokeds(
     if not do_topo:
         # add vacuous "index" (needed for topo) so same code works for both
         axes = [(ax, 0) for ax in axes]
+        assert len(axes) == 1
         if np.array(picks).ndim < 2:
             picks = [picks]  # enables zipping w/ axes
     else:
@@ -3130,7 +3129,8 @@ def click_func(
     c_func = None if do_topo else combine_func
     all_data = list()
     all_cis = list()
-    for _picks, (ax, idx) in zip(picks, axes):
+    # We need to truncate axes because of a possible additional ax for the legend
+    for ax, idx in axes[: len(picks)]:
         data_dict = dict()
         ci_dict = dict()
         for cond in conditions:
@@ -3145,7 +3145,7 @@ def click_func(
                 combine,
                 c_func,
                 ch_type=ch_type,
-                picks=_picks,
+                picks=picks[idx],
                 scaling=scalings,
                 ci_fun=ci_fun,
             )
diff --git a/mne/viz/tests/test_evoked.py b/mne/viz/tests/test_evoked.py
index c1e88acc1ef..5335f163bad 100644
--- a/mne/viz/tests/test_evoked.py
+++ b/mne/viz/tests/test_evoked.py
@@ -442,6 +442,16 @@ def test_plot_compare_evokeds(evoked):
     # test defaults
     figs = plot_compare_evokeds(evoked)
     assert len(figs) == 3
+    # test arbitrary ordering of channels is handled correctly, in topo mode
+    evoked_subset = evoked.copy().pick(["MEG 0113", "MEG 0112"])
+    fig_ordered = plot_compare_evokeds(evoked_subset, axes="topo")
+    evoked_reordered = evoked_subset.copy()
+    evoked_reordered.reorder_channels(["MEG 0112", "MEG 0113"])
+    figs_reordered = plot_compare_evokeds(evoked_reordered, axes="topo")
+    assert_allclose(
+        fig_ordered[0].axes[0].lines[0].get_ydata(),
+        figs_reordered[0].axes[0].lines[0].get_ydata(),
+    )
     # test passing more than one evoked
     red, blue = evoked.copy(), evoked.copy()
     red.comment = red.comment + "*" * 100
@@ -463,6 +473,19 @@ def test_plot_compare_evokeds(evoked):
         yvals = line.get_ydata()
         assert (yvals < ylim[1]).all()
         assert (yvals > ylim[0]).all()
+    # test that the channels are aligned when many evoked
+    # are passed in different orders, in topo mode
+    evoked_subset = evoked.copy().pick(["MEG 0113", "MEG 0112"])
+    evoked_reordered = evoked_subset.copy()
+    evoked_reordered.reorder_channels(["MEG 0112", "MEG 0113"])
+    # catch warnings when testing misalignment on purpose
+    with pytest.warns(RuntimeWarning, match="Order of channels differs"):
+        figs = plot_compare_evokeds(
+            dict(orig=evoked_subset, reordered=evoked_reordered), axes="topo"
+        )
+    assert_allclose(
+        figs[0].axes[0].lines[0].get_ydata(), figs[0].axes[0].lines[1].get_ydata()
+    )
     # test plotting eyetracking data
     plt.close("all")  # close the previous figures as to avoid a too many figs warning
     info_tmp = mne.create_info(["pupil_left"], evoked.info["sfreq"], ["pupil"])