DOC: clarify average EEG reference behavior and update tutorial note

mne/_fiff/reference.py

@@ -404,6 +404,7 @@ def set_eeg_reference(
         Array of reference data subtracted from EEG channels. This will be
         ``None`` if ``projection=True``, or if ``ref_channels`` is ``"REST"`` or a
         :class:`dict`.
+
     %(set_eeg_reference_see_also_notes)s
     """
     from ..forward import Forward

mne/utils/docs.py

@@ -4103,6 +4103,14 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75):
     EEG signal by setting ``ref_channels='average'``. Bad EEG channels are
     automatically excluded if they are properly set in ``info['bads']``.
 
+    For sensor-space analyses, if the original reference electrode is not
+    present as a channel, the resulting average reference is biased because
+    the subtracted signal is divided by n channels instead of n + 1
+    (including the reference). For a correct average reference, add a
+    zero-filled reference channel with :func:`~mne.add_reference_channels`
+    before calling :func:`~mne.set_eeg_reference`. For further discussion,
+    see Kim et al. (2023, doi:10.3389/frsip.2023.1064138).
+
 - A single electrode:
     Set ``ref_channels`` to a list containing the name of the channel that
     will act as the new reference, for example ``ref_channels=['Cz']``.

tutorials/preprocessing/55_setting_eeg_reference.py

@@ -141,6 +141,14 @@
 raw_avg_ref = raw.copy().set_eeg_reference(ref_channels="average")
 raw_avg_ref.plot()
 
+# .. note::
+#     When computing an average reference, if the dataset does not include
+#     the original reference electrode, the average will be biased because the
+#     missing channel is ignored. To avoid this, first add a zero-valued
+#     reference channel using :func:`mne.add_reference_channels`. This ensures
+#     that all intended channels, including the original reference, contribute
+#     correctly to the average.
+
 # %%
 # .. _section-avg-ref-proj:
 #