\section{Introduction}

The quantum state of a system is prone to decoherence via inhomogeneous dephasing due to variations among the system's constituents.
These variations include spatial inhomogeneities, primarily nonuniform external fields \cite{Lepoutre2018b,Munowitz1986,Kotler2011}, environmental imperfections \cite{Pingault2014,Cogan2018}, finite temperature effects such as a thermal velocity distribution \cite{finkelstein2019,finkelstein2018,Kaplan2002,Levine2018,Saffman2016}, and fabrication infidelities in engineered systems \cite{Senellart2017,Krantz2019,Devoret2013}. Slow temporal fluctuations and shot-to-shot variations, also relevant to single-constituent systems, may as well manifest as inhomogeneous dephasing  \cite{Levine2018,Manovitz2019,Aharon2013,Bermudez2012}. In quantum information processing, these variations limit the qubit coherence time.

In some cases, such as in dual-color magic-wavelength optical traps, it is possible to minimize inhomogeneous dephasing by introducing an additional field that induces the exact same inhomogeneity and thus directly balances the differential phase shifts \cite{Harber2002,Ye2008,Radnaev2010,Lahad2019,Hilton2019}.
In more general cases, inhomogeneous dephasing can be mitigated by pulse-based protocols, namely echo sequences and various dynamical decoupling methods~\cite{Viol1999,Khodjasteh2005,Souza2011,Rui2015}. In these protocols, the system regains coherence only at discrete, pre-determined times. Several protocols of continuous dynamical decoupling have been studied~ \cite{Fanchini2007,Tan2013,Bermudez2012,Golter2014}, typically relying on strong resonant driving of the qubit transition. Such driving is not compatible with all systems and requires additional measures for eliminating the sensitivity to drive noise.

Here we propose and demonstrate an alternative approach for mitigating inhomogeneous dephasing. Our scheme employs off-resonance fields for continuously dressing the target state with an additional \emph{sensor} state, which has an opposite and enhanced sensitivity to the same source of inhomogeneity. This admixture of a portion of the sensor state forms a protected, dephasing-free state.
Notably, the dressing field is a continuous wave, providing continuous protection and not limiting the extraction of quantum information to predefined times.

Experimentally, we focus on optical excitations in thermal atomic ensembles. Storing and processing quantum information as collective excitations in atomic ensembles
offers long coherence times \cite{Zhao2009,Dudin2013,Yang2016,Korber2018,Katz2018}, noiseless memory protocols \cite{finkelstein2018,kaczmarek2018}, remote entanglement \cite{Julsgaard2001,Krauter2013,Lettner2011}, and the generation of non-classical light \cite{Duan2001,Boyer2008}. Nevertheless, whether they are warm or ultracold, the atoms are not stationary; their inhomogeneous velocity distribution leads to dephasing of collective states formed by excitations with non-zero momentum transfer \cite{Whiting2017,Saffman2016,Levine2018}. Such motional dephasing is a prevailing decoherence mechanism, yet it is rarely mitigated by means beyond actual cooling and trapping of the atoms. While some pulsed spin-echo techniques have been proposed \cite{Moiseev2001} and demonstrated \cite{Rui2015}, a robust and continuous protocol for protection from motional dephasing has not been realized to date.

In the following, we first present a general analysis of the protection scheme with either one or two dressing fields. We identify the minimal requirements and conditions for optimal protection. We show that introducing two dressing fields reduces not only the sensitivity to frequency variations, but also to fluctuations in the dressing fields, known as drive noise \cite{Tan2013,Stark2018,Spielman2018}. We then turn to the case study of motional dephasing in atomic ensembles. Our experimental realization is based on light storage in a fast ladder memory (FLAME) \cite{finkelstein2018}, using the retrieval efficiency as a measure of atomic coherence~\cite{Phillips2001}. We demonstrate complete cancellation of motional dephasing, prolonging the memory lifetime, and verify the scaling of optimal conditions for protection. We further confirm that the double-dressing scheme, in contrast to the single-dressing scheme, inflicts no qubit phase rotation and is thus less sensitive to drive noise, providing robust continuous protection.

\begin{figure}[ht!]
    \centering
    \includegraphics[width=\columnwidth]{sketch_V7.pdf}
    \caption{\textbf{Continuous protection of a quantum state.} (a) Without protection, the transition frequency of the qubit $\mid\downarrow\rangle-\mid\uparrow\rangle$ is shifted by $\delta$ due to some source of inhomogeneity, while an auxiliary sensor state $|S\rangle$ experiences an opposite and possibly larger shift $-s\delta$. A single field with Rabi frequency $\Omega$ and detuning $\Delta$ dresses the qubit with the sensor state. (b) A protected qubit is formed under conditions of optimal dressing. The transition frequency of the protected qubit is slightly altered by the mean light shift. (c) Applying two dressing fields of equal intensities and opposite detunings further eliminates the mean shift and reduces the sensitivity of the protected qubit to dressing noise. (d) Optimal dressing ratio $\Omega/\Delta$ as a function of the sensitivity parameter $s$ for the single dressing (orange) and double dressing (purple) schemes. }
    \label{fig:Intro}
\end{figure}
\subsection*{Continuous protection using a sensor state}
Consider a qubit comprising the states $\mid\downarrow\rangle$ and $\mid\uparrow\rangle$, whose transition frequency experiences an inhomogeneous shift $\delta$ [Fig.~\ref{fig:Intro}(a)]. We require that one of these states, here chosen to be $\mid\uparrow\rangle$ without loss of generality, can be coupled by an external dressing field to a third state $|S\rangle$, which acts as a sensor. We further require that the overall transition frequency between $\mid\downarrow\rangle$ and $|S\rangle$ be sensitive to the same source of inhomogeneity, such that it experiences an opposite and possibly larger shift, $-s\delta$, where $s$ is the sensitivity factor. For example, in the case of motional dephasing, $s$ is the ratio between the Doppler widths of the transitions $\mid\downarrow\rangle-|S\rangle$ and $\mid\downarrow\rangle-\mid\uparrow\rangle$, and, in the case of magnetically sensitive transitions, $s$ is the ratio between their magnetic moment differences.
If the above requirements are satisfied, the dressing field can be tuned to protect the qubit from inhomogeneous dephasing.

We begin by considering a single dressing field with Rabi frequency $\Omega$ and detuning $\Delta$ [Fig.~\ref{fig:Intro}(b)]. The states $\mid\uparrow\rangle$ and $|S\rangle$ are mixed, forming dressed states. With respect to $\mid\downarrow\rangle$, the states $|S\rangle$ and $\mid\uparrow\rangle$ are inhomogeneously shifted by  \mbox{$[\Delta+\delta(s-1)]/2\pm\sqrt{[\Delta+(s+1)\delta]^2/4 +\Omega^2}$}. Optimal protection is achieved when the transition frequency of the qubit ($\mid\downarrow\rangle-\mid\uparrow\rangle$) is insensitive to variations in $\delta$ to first order, which occurs when
\begin{equation}
\label{single-dress cond}
\Omega^2/\Delta^2=s/(s-1)^2.
\end{equation}
This condition is shown in Fig.~\ref{fig:Intro}(d) by the orange line. Since the ratio $\Omega^2/\Delta^2$ determines the magnitude of mixing between the bare states, the larger the sensitivity factor $s$ is, the smaller the portion of $|S\rangle$ admixed into $\mid\uparrow\rangle$ at optimal protection.

For $s\gg1$, the dressed state frequency-shift takes the form of a light shift:  $\Omega^2/(\Delta+s\delta)\approx (\Omega^2/\Delta)[1-s\delta/\Delta+O(s\delta/\Delta)^2]$. It is the first-order term $(\Omega/\Delta)^2s\delta$ which counteracts the inhomogeneous shift $\delta$ when condition \eqref{single-dress cond} is fulfilled. Given a distribution of inhomogeneous shifts $\delta$ of width $\sigma$, a dressing field detuning $\Delta\gg s\sigma$ allows efficient protection for practically all $\delta$. This regime has been recently used for narrowing and enhancement of spectral lines \cite{Lahad2019,finkelstein2019}.

A caveat of the single-dressing scheme is the mean, zeroth-order shift $\Omega^2/\Delta$ added to the transition frequency of the protected qubit. This shift rotates the qubit phase in time, exposing it to noise in the drive field $\Omega$. To avoid this, two dressing fields at opposite detunings $\pm\Delta$ should be used [Fig.~\ref{fig:Intro}(c)]. The mean shifts they induce are opposite in sign and thus cancel. In fact, as all even-order terms cancel, a protection up to third order in $\delta/\Delta$ is achieved.

The above arguments carry on to any $s$. For a general solution with more than one dressing field, we use the Magnus expansion of the time-dependent Hamiltonian and find the condition under which the transition frequency is insensitive to variations in $\delta$ to first and second order:
\begin{equation}
\label{double-dress cond}
    J_0(2\sqrt{2}\Omega/\Delta)=(s-1)/(s+1),
\end{equation}
where $J_0$ is the zeroth-order Bessel function (see Supplementary Information for the full derivation). Equation \eqref{double-dress cond} converges to Eq.~\eqref{single-dress cond} for $s\gg 1$, but $\Omega/\Delta$ scales more favorably with $s$ for small $s$, such that less Rabi frequency and less mixing are required with double dressing [see Fig.~\ref{fig:Intro}(d)].

\begin{figure}[t!]
\centering
\includegraphics[width=\columnwidth]{New_Fig2_motionalV2.pdf}
\caption{\textbf{Motional dephasing and protection of a spin wave.} (a) Position-velocity representation of a collective excitation (spin wave), generated, {\it{e.g.}}, by light storage in a two-photon bi-chromatic transition. The collective state has a spatially dependent phase with a wavelength $\Lambda$ (top), clearly apparent when the real part of the spin-wave is averaged over all velocities (purple, bottom). (b) Thermal atomic motion (mean velocity $v_\mathrm{T}$) at time $t=2\Lambda/(\piv_\mathrm{T})$ after the initial excitation. Averaging over velocities results in decoherence of the collective state. (c) In the presence of a single dressing field, atoms with different velocities experience different light shifts, and the spatial coherence is maintained up to second order in velocity. The zeroth-order term introduces a global phase shift, evident when comparing the spin wave (purple) to its original form (gray). (d) When two dressing fields are applied, the protection is up to third order in velocity (thus covering a wider velocity range), and no global phase shift is introduced. Calculations are performed with $s=110$, $\Delta=4v_T/\Lambda$, and optimal $\Omega$ from Eqs.~(1) and (2).
}
\label{fig:dephasing}
\end{figure}
\subsection*{Motional dephasing}
Initializing an ensemble of atoms in a collective state is typically done by optical excitation. When the wavevector $k=2\pi/\Lambda$ of the collective excitation is nonzero, a spatial phase $e^{i(2\pi /\Lambda )z}$ is imprinted on the atoms, as illustrated in Fig.~\ref{fig:dephasing}(a). In multi-photon transitions, the excitation wavevector is the vectorial sum of the participating fields' wavevectors. For a single excitation, the collective state can be written as a spin wave $|W\rangle=\frac{1}{\sqrt{N}}\sum_j{e^{i(2\pi /\Lambda )z_j} \mid\uparrow\rangle_j\langle\downarrow\mid_j |G\rangle}$, where $z_j$ is the position of atom $j$ at storage time $t=0$, and $|G\rangle=\prod_j\mid\downarrow\rangle_j$ is the initial collective ground state.
When the atoms are displaced by thermal motion, atoms with different thermal velocities $v_j$ carry the original phase to different positions $z_j+v_jt$ along the spin wave [Fig.~\ref{fig:dephasing}(b)]. The resulting state is
\begin{equation}\label{Dephasing_eqn}
    |W'(t)\rangle=\sum_j{e^{i(2\pi /\Lambda )(z_j-v_j t)} \mid\uparrow\rangle_j \langle\downarrow\mid_j |G\rangle},
\end{equation}
whose overlap with $|W\rangle$, quantified by the squared coherence $\mathcal C(t)\equiv |\langle W|W'(t)\rangle|^2$, determines the retrieval efficiency of light into the phase-matched direction $\eta(t)=\eta(0)\mathcal{C}(t)$ \cite{Whiting2017}. Here $\eta(0)$ is the combined efficiency of the read and write processes. For pure motional dephasing, $\mathcal C(t)=e^{-(t/T_{\mathrm{inhom}})^2}$, where  $T_{\mathrm{inhom}}=\sqrt{2}/\sigma=\Lambda/(\sqrt{2}\piv_\mathrm{T})$ is the inhomogeneous dephasing time, and $v_\mathrm{T}$ is the atomic thermal velocity.

To apply the continuous protection scheme in the case of motional dephasing, we dress the spin-wave with a sensor state which has large and opposite sensitivity to velocity. This is achieved by an optical transition whose wavevector is larger than $2\pi/\Lambda$. The frequency shift of the transition $\mid\downarrow\rangle-\mid\uparrow\rangle$ due to the dressing field depends on the atom velocity via the Doppler effect and, for optimal dressing, exactly cancels the bare motional dephasing, rendering a velocity-insensitive state. Remarkably, as shown in Figs.~\ref{fig:dephasing}(c,d), although the atoms are constantly in motion during the evolution of the state, the spin-wave correlations between position and phase are largely maintained. This unique time-evolution is animated in Supplementary Video 1.

Apart from protecting against dephasing of the spin wave (when $|W\rangle\rightarrow |W'\rangle$), as quantified by the light retrieval efficiency, we are interested in protecting a general qubit $\alpha |G\rangle + \beta |W\rangle$. This qubit is formed, for instance, by mapping a photonic qubit $\alpha |0\rangle + \beta |1\rangle$, where $|0\rangle$ and $|1\rangle$ are photon number states, onto atomic collective states. Decoherence of $|W\rangle$ would naturally lead to decoherence of the qubit, but so will global fluctuations of the relative phase of $\beta/\alpha$ \cite{Schmidt-Eberle2019}. The double dressing scheme, which does not introduce a global phase that vary with drive noise, is therefore advantageous for robust protection of a qubit.

\begin{figure}[t!] \includegraphics[width=\columnwidth]{New_Fig3_Storage_V3.pdf}
\caption{\textbf{Experimental realization of continuous protection from motional dephasing.} (a) Light storage scheme: Counter-propagating pulses of a strong control field and a weak signal field with wavevectors $k_\mathrm{c}$  and $k_{\mathrm{s}}$ couple the ground state $\mid\downarrow\rangle$ to the excited state $\mid\uparrow\rangle$ through an intermediate state $|i\rangle$. A dressing beam with wavevector $k_{\mathrm{d}}$ couples $\mid\uparrow\rangle$ to the sensor state $|S\rangle$ with Rabi frequency $\Omega$ and detuning $\Delta$. (b) The signal and control pulses generate a collective excitation with wavevector $2\pi/\Lambda = k_{\mathrm{s}}-k_\mathrm{c}$. The dressing field wavevector $k_{\mathrm{d}}$ has an opposite direction and a larger magnitude. (c) Typical data traces when the incoming and retrieved signals are either directly measured (top) or interfered with reference pulses to measure the phase of the signal (bottom).\label{fig:Exp_setup}}
\end{figure}

\subsection*{Experimental realization} To demonstrate and study the continuous protection against motional dephasing, we probe the coherence decay of a collective excitation in a hot atomic vapor, utilizing storage and retrieval of light. We employ the FLAME protocol \cite{finkelstein2018} to store 1.8-ns (FWHM) signal pulses (with $\sim 0.5$ photons per pulse) as electronic-orbital excitations  [Fig.~\ref{fig:Exp_setup}(a)] in a 5-mm-long cell of thermal $^{87}\mathrm{Rb}$ vapor at 98$^{\circ}$C. The $\mid\uparrow\rangle$ and $\mid\downarrow\rangle$ states reside in the levels $5S_{1/2}$ and $5D_{5/2}$ respectively. The two states are coupled via a two-photon transition, detuned by 1.4~GHz from the intermediate level $5P_{3/2}$ [$|i\rangle$ in Fig.~\ref{fig:Exp_setup}(a)]. After some storage time $t$, we send a second control pulse and retrieve a signal pulse from the collective excitation [Fig. 3(c) top]. The memory efficiency $\eta(t)/\eta(0)$ is determined from the retrieved power (relative to the input power), and the temporal decay of the squared coherence is calculated from $\mathcal{C}(t)=\eta(t)/\eta(0)$. To measure the phase difference between incoming and retrieved signals, we store slightly longer signal pulses (4-ns FWHM) and interfere them with reference pulses shifted by $2\pi\cdot 380$ MHz, which do not traverse the atomic vapor [Fig.~\ref{fig:Exp_setup}(c) bottom].

The wavevectors for the signal and control transitions are respectively $k_{\mathrm{s}}=2\pi/0.78~\mathrm{\mu m}^{-1}$ and $k_\mathrm{c}=2\pi/0.776~\mathrm{\mu m}^{-1}$, giving an excitation wavelength $\Lambda=2\pi/(k_{\mathrm{s}}-k_\mathrm{c})= 151~\mathrm{\mu m}$ [Fig.~\ref{fig:Exp_setup}(b)]. The Gaussian distribution of velocities with a thermal velocity $v_T=188$~m/s results in an inhomogeneous dephasing of the collective state over about 100 ns. Protection during storage is obtained by a field with wavevector $k_{\mathrm{d}}=2\pi/1.274~\mathrm{\mu m}^{-1}$, propagating along the direction of the signal field and weakly dressing the state $\mid\uparrow\rangle$ with the state $|S\rangle=28F_{7/2}$. We use an electro-optic modulator (EOM) to generate the double dressing configuration (see Methods). The sensor transition $\mid\downarrow\rangle-|S\rangle$ has an opposite and enhanced velocity-sensitivity $s=(k_{\mathrm{s}}-k_\mathrm{c}+k_{\mathrm{d}})/(k_{\mathrm{s}}-k_\mathrm{c}) \approx 110$ compared to that of the two-photon transition $\mid\downarrow\rangle-\mid\uparrow\rangle$, such that weak dressing is sufficient to form a velocity-insensitive collective state.
The bare homogeneous decoherence rate of the collective state $\gamma=2\pi\cdot 1.38$ MHz is limited by the $5D_{5/2}$ radiative lifetime and by the transit time of atoms through the signal beam waist. In the presence of a dressing beam, we measure an increase of $\le 25$\% in $\gamma$, which we attribute to resonant scattering due to imperfect extinction of the carrier frequency by the EOM and to residual excitation of atoms in the tail of the velocity distribution for small dressing field detunings (see Methods).

\begin{figure}[tb!]
    \centering
    \includegraphics[width=\columnwidth]{FigTime_resub.pdf}
    \caption{\textbf{Protection of a collective excitation from inhomogeneous dephasing.} The squared coherence is calculated for different storage times by measuring the memory efficiency. Solid lines are fit to a model accounting for both homogeneous and inhomogeneous dephasing. With double dressing, the decay is predominantly homogeneous; the dashed line marks the homogeneous component of the fitted decay model, with surrounding shaded area representing fit uncertainty. Inset: The inhomogeneous component of the decoherence $\mathcal{C}(t)/e^{-\gamma t}$, where the homogeneous component $e^{-\gamma t}$ is normalized out. The dressing field parameters are $|\Delta|=2\pi\cdot350$ MHz, $\Omega=2\pi\cdot34$ MHz (double dressing), and $\Omega=2\pi\cdot28$ MHz (single dressing). \label{fig:Time}}

\end{figure}

To quantify the effect of protection, we compare the decay of coherence with no dressing field, with a single dressing field, and with a double dressing field; these are shown in Fig.~\ref{fig:Time}. In the absence of dressing fields, the decay of efficiency with time has both exponential and Gaussian components, corresponding respectively to homogeneous and inhomogeneous dephasing. With a single dressing field, the inhomogeneous component is partially removed, while with the double dressing field, it vanishes almost completely. As a result, we observe an increase of over an order of magnitude in retrieved signal at long storage times.

The residual inhomogeneous dephasing is revealed by normalizing the squared coherence $\mathcal{C}(t)$ by the homogeneous decoherence $e^{-\gamma t}$, as plotted in the inset of Fig.~\ref{fig:Time}. With a double-dressing protection field, the inhomogeneous dephasing time is substantially elongated. Importantly, this implies that a dramatic improvement in the coherence times can be obtained for the broad range of systems that enjoy low homogeneous decoherence rate. The remaining, non-linear dependence of the qubit transition frequency on the inhomogeneous shift $\delta$ results in non-Gaussian dephasing \cite{Sung2019}. On long time scales, the coherence decays as a power-law $\mathcal C(t) \propto t^{-2/n}$, where $n$ is the first non-vanishing order of $\delta$ in the qubit transition frequency. However on a short time scale, the dephasing is well approximated by a Gaussian decay (see Supplementary Information).

 \begin{figure}[t!]
\centering
\includegraphics[width=\columnwidth]{Fig5_resub_sept20.pdf}
\caption{\textbf{Dependency of coherence time and accumulated phase on the dressing parameters.} (a) Inhomogeneous dephasing time $T_{\mathrm{inhom}}$, with double-dressing protection as a function of dressing Rabi frequency $\Omega$, for different detunings $\Delta$ (in $2\pi\cdot$MHz). Shaded gray line is a theoretical model, where the dephasing rates due to linear and cubic dependencies on $\delta$ are summed in quadrature. (b) The ratio $\Omega_{\mathrm{opt}}/\Delta$, where $\Omega_{\mathrm{opt}}$ is the Rabi frequency maximizing the gain for each $\Delta$, compared to the model prediction $2^{-3/2}J_0^{-1}[(s-1)/(s+1)]\approx s^{-1}\approx 0.095$ (solid line). (c) Phase accumulated due to the dressing fields during a storage time of 40 ns. Here $\Delta=2\pi\cdot480$ MHz. Solid lines are linear fits. For double dressing, the added phase is zero (up to the experimental uncertainty) and independent of the dressing intensity.}
\label{fig:DressPower}
\end{figure}

To further study the optimal conditions for continuous protection, we vary the dressing parameters $\Omega$ and $\Delta$ and measure the retrieval efficiency at $t=250$ ns, from which we infer the dephasing time $T_{\mathrm{inhom}}$ (see Methods). Figure \ref{fig:DressPower}(a) shows the improvement in $T_{\mathrm{inhom}}$ for the double-dressing protection scheme. The data agrees with the condition $|\Delta|>s\sigma$ for efficient protection (here $s\sigma=2\pi\cdot140$ MHz). Substantial improvement in dephasing time is obtained over a large range of dressing powers and is further enhanced with larger detuning. For the largest detuning measured, we compare the data to the theoretical model (shaded gray line). Indeed, the improvement in $T_{\mathrm{inhom}}$ is quadratic in $\Omega$ and reaches an optimal value at the protection condition \eqref{double-dress cond}, where it is limited by the third-order term in the inhomogeneous shift. The latter can be made negligible by increasing the detuning. We further validate the protection condition \eqref{double-dress cond} in Fig.~\ref{fig:DressPower}(b) by comparing the predicted ratio $\Omega/\Delta$ to the measured ratio $\Omega_{\mathrm{opt}}(\Delta)/\Delta$. Excellent agreement is found for large detunings.

An important property of any protection or dynamical decoupling technique is the qubit phase introduced by the driving fields and its dependence on their intensities. Measurements of this dependence for both single and double dressing are shown in Fig.~\ref{fig:DressPower}(c). For single dressing, we find indeed that the phase shift scales linearly with the mean light shift $\Omega^2/\Delta$ \cite{Parniak2019}, yielding a sensitivity to drive noise of $\sim$ 0.4 mrad/$(2\pi\cdot$MHz$)^2$. This is in contrast to the double dressing protection, where we measure essentially no added phase and no dependence on the dressing intensity or on the storage time.

\subsection*{Potential application to other systems}
Motional dephasing occurs across various platforms, where it limits the performance of quantum information applications. One prominent example is that of neutral atoms trapped in tweezer arrays and interacting via Rydberg excitations. The coherence time and bell-state preparation fidelity in several recent realizations of this system are limited by motional dephasing \cite{Levine2018,Graham2019}. Our scheme can be adapted to mitigate this effect by using a uniform field for globally dressing either the ground or the Rydberg state with an auxiliary meta-stable state or by implementing velocity-insensitive bichromatic entangling gates. An entirely different application is a single photon source based on four-wave mixing in atomic vapor. Whether they rely on Rydberg interactions or not \cite{Ripka2018,Lee2016}, the brightness and purity of these  sources are limited as well by motional dephasing, which can be suppressed by implementing our scheme. Here, when it comes to stochastic events such as the heralded spontaneous generation of photons, a protection scheme that is continuous is essential. For collective Rydberg excitations, a three-photon excitation to the Rydberg state can be employed with an intermediate (two-photon) transition that is velocity sensitive (\emph{i.e.}, with larger sensitivity than in our experiment); this seemingly peculiar choice would yield a state which is protected from motional dephasing and, at the same time, has a large Rydberg admixing, such that it inherits the desired Rydberg-level qualities, namely strong dipolar interactions.

Beyond the context of motional dephasing, our scheme can be applied to magnetically sensitive qubits, such as trapped ions and NV centers. Here again, one should identify a sensor transition with opposite and possibly larger magnetic moment difference relative to the qubit transition. For example, we suggest protecting an optical qubit in a $^{88}\mathrm{Sr}^+ $ ion encoded in the $\ket{S_{1/2};m_j{=}1/2}-\ket{D_{5/2};m_j{=}1/2}$ basis by dressing the ground state with the long-lived $\ket{D_{5/2};m_j{=}-3/2}$ state. This choice provides a sensitivity factor $s=6$ owing to the different land\'{e} $g$-factors of the states involved and thus entails low admixing of the sensor level into the qubit state.

\subsection*{Discussion} The presented protection scheme is akin to continuous dynamical decoupling but distinct in several aspects. Our protecting fields do not couple the two states comprising the qubit, but rather couple one of them to an external state. They are preferably off-resonant, and the overall protection improves with further detuning. In standard continuous dynamical decoupling, on the other hand, the protecting fields drive the qubit transition, often resonantly, exposing it directly to drive noise $\delta_\Omega$ when the field amplitude fluctuates $\Omega \rightarrow \Omega+\delta_\Omega$. Moreover, if $\delta$ fluctuates in time, the fluctuation bandwidth over which the protection is effective is limited by the available Rabi frequency \cite{Stark2017}. In our scheme, the bandwidth is fundamentally limited by the detuning (or  the effective Rabi frequency $\sqrt{4\Omega^2+\Delta^2}$), which can be made much larger given the same limited power, especially with large sensitivity $s$ (for which $\Delta \gg \Omega$)

In Table \ref{tab:summary}, we present the effectiveness of our protection schemes in terms of the residual (higher-order) sensitivities to inhomogeneous shift O$(\delta^2)$,O$(\delta^3)$; to Rabi frequency fluctuations O$(\delta_\Omega)$; and to cross-terms of the two O$(\delta\delta_\Omega)$. For a single dressing field, an enhanced sensitivity $s\gg1$ earned from utilizing a sensor state enables the use of a protection field that is far-detuned, thus reducing sensitivity to drive noise to $\delta_\Omega/\sqrt{s}$. Further suppression of drive noise is achieved with double dressing, which we find to be insensitive to drive noise to leading order O($\delta_\Omega$) [also supported by the measurements in Fig. \ref{fig:DressPower}(c)]. Our analysis additionally shows that the double dressing removes sensitivity to inhomogeneous shifts up to O$(\delta^3)$. The remaining sensitivity scales inversely with Rabi frequency or with detuning, but we note that it scales unfavourably with the sensitivity factor $s$. The dephasing time $T_{\mathrm{inhom}}$ is proportional to the inverse of the leading residual sensitivity term when substituting $\delta=\sigma$.

For a sensor state with sensitivity equal to that of the qubit state, $s=1$, the single-dressing protection condition [Eq. \eqref{single-dress cond}] requires a resonant field $\Delta=0$, which results in two, equally mixed, protected states \cite{Stark2018,Spielman2018}. The double dressing scheme also generates two protected states in this case (see Supplementary Information), without the need for resonant driving. The $s=1$ case resembles the case of a two-level system, where the condition \eqref{double-dress cond} provides for intrinsic dynamical decoupling in bichromatic entangling gates with trapped ions \cite{Sutherland2019}.

\begin{table}[tb!]
\caption{\textbf{Leading residual sensitivities of the qubit transition frequency:} first-order in the drive noise ($\delta_\Omega$), second- and third- order in the inhomogeneous shift ($\delta^2$, $\delta^3$), and the cross term (${\delta_\Omega\cdot\delta}$). For the double dressing scheme, exact numerical factors depend on the experimental realization, including the phase between the two fields. These factors and their exact derivation are given in the Supplementary Information.  }\label{tab:summary}
\begin{tabular}{ |c||c|c||c|c|c|c|  }
\colrule
 & \multicolumn{2}{c||}{{single dressing}} & \multicolumn{4}{c|}{{double dressing}} \\ [0.5ex]

& O($\delta_\Omega$) &  O($\delta^2$) & O($\delta_\Omega$) &  O($\delta^2$) & O($\delta \delta_\Omega$) &  O($\delta^3$) \\ [0.5ex]
\colrule
& & & & & &\\ [-2ex]
 $s=1$	& $\delta_\Omega$ & $\frac{1}{2}\frac{\delta^2}{\Omega}$
 & 0 & 0 & $\sim1\frac{\delta\delta_\Omega}{\Omega}$ & $\sim\frac{1}{2}\frac{\delta^3}{\Omega^2}$\\[1ex]

 $s\gg 1$ &	$\frac{2}{\sqrt{s}}\delta_\Omega$ & $\sqrt{s}\frac{\delta^2}{\Omega}$
 & 0	& 0 & $2\frac{\delta\delta_\Omega}{\Omega}$ & $\sim s\frac{\delta^3}{\Omega^2}$\\[1ex]
\colrule
 \end{tabular}
 \end{table}

In conclusion, we have introduced and demonstrated a new scheme for protection from inhomogeneous dephasing, which is continuous, efficient, and robust to drive noise. The scheme is particularly suitable for protection from motional dephasing of a collective state stored in an atomic ensemble, where limited solutions were suggested to date. The minimal requirements outlined in this paper may be found across numerous systems where a multi-level structure exists. We have discussed a few examples for such systems, where this scheme can be applied to enhance the performance of quantum sensors, sources and gates.

Our experimental demonstration confirms the validity of the scheme by eliminating the inhomogeneous dephasing of a collective excitation in a gas of thermal atoms. Remarkably, this is achieved not through time-reversal of the process of dephasing or of the direction of atomic motion, but through continuously maintaining the original position-dependent phases of atoms which move randomly to different positions.

The compatibility of this scheme to various protocols of quantum information processing, including gates and metrology, requires further research. It could provide an important tool for these applications operating with either single-constituent qubits or with collective quantum states of ensembles.

\subsection*{Methods} \emph{Experimental design. --}
The setup comprises a 780 nm distributed Bragg reflector (DBR) diode laser, serving as the signal beam, and a 776 nm external cavity diode laser (ECDL) amplified by a tapered amplifier (TA), serving as the control beam. The signal laser is offset-locked to a master ultra-stable fiber laser using a fast beat-note detector. The control laser is locked to a two-photon absorption feature in a reference cell when overlapped with the master laser, where the latter is frequency shifted by a fiber electro-optic phase modulator (EOPM). The signal field is amplitude modulated in time by two fiber electro-optic amplitude modulators (EOAM) to carve a Gaussian pulse of 1.8 ns FWHM, with a combined extinction ratio of 1:3000. The control field is amplitude modulated by two Pockels cells (PCs), generating pulses of 2.5 ns FWHM with an extinction ratio of 1:1000. The repetition rate of the experiment is set by that of the PCs to 100~kHz. After the modulators, the control beam is passed through a tilted filter  (Semrock LL01-780-12.5) functioning as a 776 nm bandpass, filtering out other frequencies that might be produced in the TA. Both beams are passed through single mode fibers (SMF), aligned with each other in a counter-propagating geometry, and overlapped at the center the vapor cell. The signal beam is focused down to $w_0=85~\mu$m, while the control beam is focused down to $w_0=200~\mu$m. Both beams are $\sigma^+$ polarized, and we optically pump all atoms to the stretched state, a combination which guarantees purely orbital transitions and greatly simplifies the multi-level structure of the atomic vapor. The optical pumping is realized by a `pump' and a `repump' at 795 nm, both $\sigma^+$ polarized, and separated from  each other by 6.8 GHz, such that the pump (repump) is resonant with the $F=2 \rightarrow F'=2$ ($F=1 \rightarrow F'=2$) transition of the rubidium D1 line.  The pump (repump) beam power at the vapor cell is 300~mW (200~mW). The pump beams are 1.2 mm wide and are directed at a small angle with respect to the control beam.
The 5-mm-long $^{87}$Rb vapor cell is anti-reflection coated for 780-1064 nm. It is heated to 72$^\circ$C at its coldest spot and 98$^\circ$C at its hottest spot using two electrical current heaters, to set a Rb density of $6.5\times10^{11}$~cm$^{-3}$ and an optical depth OD$\approx5$. We obtain OD$\approx8.5$ with continuous optical pumping.
After the cell, the signal beam is passed through a polarizing beam-splitter and two 780 nm bandpass filters to filter out any residual 776 nm and 795 nm components. It is then coupled to a SMF acting as a spatial filter, removing most of the spatially incoherent fluorescence emitted from the cell at 780 nm. The SMF is coupled either to a fast linear avalanche photo detector (APD) with bandwidth of 1~GHz, or to a single photon counting module (SPCM) connected to a time tagger with time bins of 100~ps.

\emph{Dressing fields generation.--} The dressing beam is produced by a 1274 nm ECDL and is then passed through both an EOAM and an EOPM. For the double-dressing protection scheme, the EOPM is modulated by a square-wave electronic signal (70 ps rise/fall time), generating dominant first-order side bands comprising 78\% of the total output power and a vanishingly small  carrier ($0.2\%$). The dressing beam is then amplified by an O-band booster optical amplifier (Thorlabs BOA1130S) and further amplified by a TA. The dressing laser frequency is stabilized to a wavelength meter with 1 MHz resolution (High-Finesse WS8-IR1). The dressing beam is combined with the signal beam on a dichroic mirror and focused to $w_0=200~\mu$m in the vapor cell. It is also $\sigma^+$ polarized, thus coupling only to the stretched state.

\emph{Homogeneous decoherence in the experiment.--}
With a single dressing field and, more significantly, with a double dressing field, we measure a small increase in the homogeneous decoherence rate (of up to 25\%). The first contribution to this increase is the scattering rate of atoms at the tale of the inhomogeneous distribution. For the intermediate detunings we use in the experiment $|\Delta| \lesssim 3\sigma$, there is a small but non-vanishing fraction of atoms whose transition frequency is near-resonant with the dressing field. The overall scattering rate due to this contribution is on the order of $1.2{-}6.3\cdot10^6 ~\mathrm{s}^{-1}$ in the range of parameters used in the experiment and scales inversely with $\Delta^2$.

For the double dressing scheme, the experimental realization (as detailed above) includes a single laser tuned to resonance and a fast square-wave phase modulation, whose imperfections result in a carrier intensity of up to $\Omega^2_\mathrm{carrier}/\Omega^2=0.3(2)\%$. This corresponds to a calculated resonant scattering rate of up to $0.6\cdot 10^5~ \mathrm{s}^{-1}$ (at the optimal Rabi frequency of $\Omega=2\pi\cdot 35$ MHz).

\emph{Evaluation of $T_{\mathrm{inhom}}$ .--}
To obtain the inhomogeneous dephasing time, we measure the memory efficiency for a long storage time of 250 ns. The results of this measurement are presented in the Supplementary Information (Fig. S2). We extract $T_{\mathrm{inhom}}$ from the efficiency by fitting to a Gaussian decay model $e^{-(t/T_{\mathrm{inhom}})^2}$ after normalizing out the homogeneous decoherence $e^{-\gamma t}$. We include the slight increase in $\gamma$ in the presence of the dressing field by measuring $\gamma$ for several dressing powers at $\Delta=2\pi\cdot350$ MHz (as in Fig. 4 in the main text) and calibrate our model for resonant scattering $\gamma=\gamma_0+\Gamma\Omega^2/\Delta^2$ accordingly. We verify the validity of the model across a wide range of dressing field parameters as presented in Fig. S2 (Supplementary Information). Larger available laser power and faster modulation of the dressing field would enable operation at larger detunings which would greatly suppress such scattering.

\emph{Phase measurements.--} In order to measure the phase between the incoming and retrieved signals, we generate two consecutive signal pulses and pass them through an acousto-optical modulator (AOM). The diffracted pulses, shifted by $2\pi\cdot380$ MHz relative to the signal, will act as a reference. They are thus coupled into an optical fiber acting as a delay line. The time difference between the consecutive pulses, the memory storage time, and the time delay in the delay line are all set to be identical. This allows part of the incoming signal, which was not stored due to limited efficiency, to interfere with the first reference pulse, and the retrieved signal to interfere with the second reference pulse. The phase difference between these two pulses, acquired in a single shot within tens of ns, is thus insensitive to interferometer drifts, which occur over a much longer time scale. We then average this value over hundreds of repetitions.

\emph{Author contributions--} R.F., O.L., O.D., S.K., and E.P. contributed to the experimental design, construction, data collection and analysis of this experiment. I.C., R.F., O.L., and E.P. developed the theoretical framework supporting the experiment. E.P. and O.F. supervised the entire project. All authors discussed the results and contributed to writing
the manuscript.

\emph{Competing interests}-- The authors declare no competing interests.

\emph{Acknowledgements--} We acknowledge financial support by the Israel Science Foundation and ICORE, the
European Research Council starting investigator grant QPHOTONICS 678674, the Pazy Foundation, the Minerva
Foundation with funding from the Federal German Ministry
for Education and Research, and the Laboratory in
Memory of Leon and Blacky Broder. I.C. acknowledges support from Marie Skodowska-Curie grant agreement no. 785902.

\newpage
\bibliographystyle{latexmkrc}


\section{Introduction} \label{S1}

Within a sublayer of wall turbulence of an incompressible fluid, one-point statistics such as the mean velocity and the variances of some velocity fluctuations vary logarithmically with the distance from the wall. This logarithmic scaling is unusual, contrasting to power laws and exponential laws found in many other systems. Hence, in wall turbulence, we are to study scaling laws of two-point statistics.

The configuration is as follows. We take the $x$--$y$ plane at the wall. The $x$ direction is that of the mean stream. While $U(z)$ denotes the mean velocity at a distance $z$ from the wall, $u(z)$ and $w(z)$ denote velocity fluctuations in the streamwise and the wall-normal directions. Turbulence is homogeneous in the streamwise direction. Its thickness $\delta$ is a constant. The two points considered here are those separated by a streamwise distance $r$.

Asymptotically in the limit of high Reynolds number, there is a sublayer at $z/\delta \rightarrow 0$ such that the momentum flux $\rho \langle -uw(z) \rangle$ is constant at a value of $\rho u_{\tau}^2$ \cite{my71}. Here $\rho$ is the mass density, $u_{\tau}$ is the friction velocity, and $\langle \cdot \rangle$ denotes an average. Even in an actual case over a smooth or rough wall at a high but yet finite Reynolds number, this constant-flux sublayer is still a good approximation for a range of distances $z$.

Throughout the constant-flux sublayer, the friction velocity $u_{\tau}$ serves as a characteristic velocity. Since there is no constant in units of length, the mean velocity $U$ obeys a relation $\partial U/\partial z \varpropto u_{\tau}/z$ \cite{ll59, s48, my71}. Then,
\begin{subequations} \label{eq1}
\begin{equation} \label{eq1a}
\frac{U(z)}{u_{\tau}} = c_U + d_U \ln \! \left( \frac{\delta}{z} \right)  \!
\ \ \mbox{with} \ \
d_U = - \frac{1}{\kappa} .
\end{equation}
Here $c_U$ is an integration constant. The von K\'arm\'an constant $\kappa$ appears to be universal. Its estimate of $0.39 \pm 0.02$ is common among various configurations of wall turbulence, e.g., pipe flows, channel flows, and boundary layers \cite{mmhs13}.

The same scaling exists for the variance of streamwise velocity fluctuations $\langle u^2(z) \rangle$. According to the attached-eddy hypothesis of Townsend \cite{t76}, i.e., a model of a random superposition of energy-containing eddies that are attached to the wall,
\begin{equation} \label{eq1b}
\frac{\langle u^2(z) \rangle}{u_{\tau}^2} = c_{u^2} + d_{u^2} \! \ln \! \left( \frac{\delta}{z} \right) \!  .
\end{equation}
This law has been confirmed recently by means of laboratory experiments and field observations \cite{mmhs13, hvbs12}. As for its constants, while $c_{u^2} \simeq 1.4$--$1.8$ in pipe flows is distinct from $c_{u^2} \simeq 2.0$--$2.5$ in channel flows and boundary layers, $d_{u^2} \simeq 1.2$--$1.3$ is common among them \cite{mmhs13, hvbs13, vhs15, ofsbta17, smhffhs18, hvbs12}.

Other scaling laws are also known. An example is the local rate per unit mass of the energy dissipation $\varepsilon (z)$. By equating its average $\langle \varepsilon \rangle$ to the mean rate of the energy production $\langle -uw \rangle \partial U/\partial z = u_{\tau}^3/\kappa z$ at each distance $z$ in the constant-flux sublayer \cite{ll59},
\begin{equation} \label{eq1c}
\frac{\langle \varepsilon(z) \rangle}{u_{\tau}^3/z} = c_{\varepsilon}
\ \ \mbox{with} \ \
c_{\varepsilon} = \frac{1}{\kappa} .
\end{equation}
\end{subequations}
This is in accordance with the attached-eddy hypothesis \cite{m17}, albeit possibly not exact with a discrepancy of $\pm 10$\% in the actual flow \cite{lm15}.

Being analogous to such one-point statistics, two-point statistics would exhibit some scaling laws. These are expected to offer much more information about the wall turbulence.

We are to use cumulants $\langle \alpha^n \beta^m \rangle_c$ of random variables $\alpha$ and $\beta$ at $n$ and $m = 0$, $1$, $2$, ... \cite{my71, ks77}. They are related to usual moments $\langle \alpha^l \beta^k \rangle$ at $l$ and $k = 0$, $1$, $2$, ....  For example \cite{c51},
\begin{subequations} \label{eq2}
\begin{align}
\label{eq2a} \langle \alpha          \rangle_c &=   \langle \alpha           \rangle                                                                  , \\
\label{eq2b} \langle \alpha  \beta   \rangle_c &=   \langle \alpha   \beta   \rangle - \langle \alpha         \rangle   \langle \beta \rangle         , \\
\label{eq2c} \langle \alpha^2\beta^2 \rangle_c &=   \langle \alpha^2 \beta^2 \rangle
                                                 -  \langle \alpha^2         \rangle   \langle        \beta^2 \rangle
                                                 -2 \langle \alpha   \beta   \rangle^2                                                        \nonumber \\
                                               & -2 \langle \alpha^2 \beta   \rangle   \langle \beta          \rangle
                                                 -2 \langle \alpha           \rangle   \langle \alpha \beta^2 \rangle
                                                 +8 \langle \alpha           \rangle   \langle \alpha \beta   \rangle   \langle \beta \rangle \nonumber \\
                                               & -6 \langle \alpha           \rangle^2 \langle        \beta   \rangle^2
                                                 +2 \langle \alpha^2         \rangle   \langle        \beta   \rangle^2
                                                 +2 \langle \alpha           \rangle^2 \langle        \beta^2 \rangle.
\end{align}
\end{subequations}
At $n+m \ge 2$, each moment $\langle \alpha^n \beta^m \rangle$ is contaminated nonlinearly with lower order moments. If all of such contamination is removed, the result is the cumulant $\langle \alpha^n \beta^m \rangle_c$. It is also identical to $\langle (\alpha - \langle \alpha \rangle )^n (\beta - \langle \beta \rangle)^m \rangle_c$. For a sum of independent random variables $\alpha_1$ and $\alpha_2$, each cumulant is identical to the sum of cumulants of the variables, i.e., $\langle (\alpha_1+\alpha_2)^n \beta^m \rangle_c = \langle \alpha_1^n \beta^m \rangle_c + \langle \alpha_2^n \beta^m \rangle_c$.

This linear character of cumulants would lead to simple scaling laws in the constant-flux sublayer.  Actually from the attached-eddy hypothesis \cite{m17},
\begin{subequations} \label{eq3}
\begin{equation} \label{eq3a}
\frac{\langle u^n(x,z)u^m(x+r,z) \rangle_c}{u_{\tau}^{n+m}} = C_{u^nu^m} \! \! \left( \frac{r}{z} \right) + D_{u^nu^m} \! \! \left( \frac{r}{z} \right) \! \ln \! \left( \frac{\delta}{z} \right)
\end{equation}
and
\begin{equation} \label{eq3b}
\frac{\langle \varepsilon^n(x,z) \varepsilon^m(x+r,z) \rangle_c}{(u_{\tau}^3/z)^{n+m}} = C_{\varepsilon^n \varepsilon^m} \! \! \left( \frac{r}{z} \right) \! .
\end{equation}
\end{subequations}
Here $x$ has become a dummy parameter because the turbulence is homogeneous in the streamwise direction. The functions $C_{u^nu^m}$, $D_{u^nu^m}$, and $C_{\varepsilon^n\varepsilon^m}$ are not yet determined but are related to constants of the one-point statistics. For example, $c_{u^2}$ in Eq.~(\ref{eq1b}) is identical to $C_{uu}(0)$ as well as to $C_{u^2}(0) = C_{u^2u^0}(0)$, while $c_{\varepsilon}$ in Eq.~(\ref{eq1c}) is identical to $C_{\varepsilon}(0) = C_{\varepsilon \varepsilon^0}(0)$.

We note that Eq.~(\ref{eq3a}) has been derived without adding any assumption to the original hypothesis of Townsend \cite{t76}. Although previous studies added assumptions \cite{pc82, dnk06, mizuno18, mm19}, they are not consistent with that hypothesis if their results do not satisfy Eq.~(\ref{eq3a}), aside from whether they appear reasonable or not \cite{m17,m19}.

The scaling laws of such two-point statistics are studied here. With some mathematics (Sec.~\ref{S2}), Eq.~(\ref{eq3}) is reproduced as an extension of Eq.~(\ref{eq1a}) for the mean velocity $U$. Then, by using data obtained from our experiments of boundary layers (Sec.~\ref{S3}), we confirm Eq.~(\ref{eq3}) in Sec.~\ref{S4}. Since its functions $C_{u^nu^m}$, $D_{u^nu^m}$, and $C_{\varepsilon^n\varepsilon^m}$ are not dependent on the above mathematics, they are discussed in terms of the attached-eddy hypothesis (Sec.~\ref{S5}). Finally, we conclude with remarks in Sec. \ref{S6}.

\section{Theory} \label{S2}

Our theory is to extend the scaling law of Eq.~(\ref{eq1a}) for the mean velocity $U$. Within the constant-flux sublayer at $z/\delta \rightarrow 0$, the local gradient $\partial U/ \partial z$ depends only on the local parameters $z$ and $u_{\tau}$ \cite{ll59, s48, my71}. Thus, we have used $\partial U/ \partial z \varpropto u_{\tau}/z$ to obtain Eq.~(\ref{eq1a}). Although  $\gamma \ne 1$ for $\partial U^{\gamma}/ \partial z \varpropto u_{\tau}^{\gamma}/z$ might appear equally plausible, this is not invariant under a Galilean transformation to add a constant to all of $U$ \cite{o01}.

The scaling law for any other statistics is required to be consistent with that for the mean velocity $U$. We rely on this requirement to constrain the former scaling via cumulants of the total streamwise velocity $U+u$. Another basis of our theory is that any nondimensional function is required to be described by nondimensional parameters alone.

With use of the friction velocity $u_{\tau}$, the total streamwise velocity $U+u$ is nondimensionalized as
\begin{equation} \label{eq4}
\alpha = \frac{U(z)+u(x,z)}{u_{\tau}}
\ \ \mbox{and} \ \
\beta = \frac{U(z)+u(x+r,z)}{u_{\tau}}.
\end{equation}
The distribution of $\alpha$ is described completely by its characteristic function $\langle e^{is\alpha} \rangle$ with a nondimensional parameter $s$ ranging from $-\infty$ to $+\infty$. By definition \cite{my71, ks77}, the cumulants $\langle \alpha^n \rangle_c$ are obtained from
\begin{subequations} \label{eq5}
\begin{equation} \label{eq5a}
\ln \langle e^{is\alpha} \rangle = \sum_{n=1}^{\infty} \langle \alpha^n \rangle_c \frac{(is)^n}{n!}
\end{equation}
or from
\begin{equation} \label{eq5b}
 \langle \alpha^n \rangle_c = \left. \frac{\partial^n}{\partial (is)^n} \ln \langle e^{is\alpha} \rangle \right\vert_{s=0} .
\end{equation}
For consistency with $\partial U/ \partial z \varpropto u_{\tau}/z$, we impose
\begin{equation} \label{eq5c}
\frac{\partial}{\partial z} \ln \langle e^{is\alpha} \rangle = - \frac{\phi (s)}{z} .
\end{equation}
This relation is still invariant under the aforementioned Galilean transformation, which affects only a linear term of $\langle \alpha \rangle_c = \langle \alpha \rangle = U/u_{\tau}$ at $n=1$ in Eq.~(\ref{eq5a}). Since $\phi$ is a nondimensional function, it does not depend on $z$ that is not nondimensional. From Eqs.~(\ref{eq5b}) and (\ref{eq5c}),
\begin{equation} \label{eq5d}
\frac{\partial}{\partial z} \langle \alpha^n \rangle_c = - \left. \frac{1}{z} \frac{\partial^n \phi (s)}{\partial (is)^n} \right\vert_{s=0} = - \frac{\phi^{(n)}(0)}{i^n z}.
\end{equation}
\end{subequations}
We replace $\phi^{(n)}(0)/i^n$ with a constant $d_{(U+u)^n}$. Because of $\langle \alpha^n \rangle_c = \langle (U+u)^n \rangle_c/u_{\tau}^n$,
\begin{equation} \label{eq6}
\frac{\partial}{\partial z} \frac{\langle [U(z)+u(z)]^n \rangle_c}{u_{\tau}^n} = - \frac{d_{(U+u)^n}}{z}.
\end{equation}
The integration of Eq.~(\ref{eq6}) leads to the logarithmic laws of Eq.~(\ref{eq1a}) for $\langle U+u \rangle_c = U$ via $d_{U+u} = d_U = -1/\kappa$ and of Eq.~(\ref{eq1b}) for $\langle (U+u)^2 \rangle_c = \langle u^2 \rangle_c  = \langle u^2 \rangle$ via $d_{(U+u)^2} = d_{u^2}$. Such a law is obtained also for $\langle (U+u)^4 \rangle_c =\langle u^4 \rangle_c = \langle u^4 \rangle - 3 \langle u^2 \rangle^2$ and so on. However, the corresponding law for a moment $\langle u^l \rangle$ is usually not simple because $\langle u^l \rangle$ is contaminated nonlinearly with cumulants $\langle u^n \rangle_c$ of orders $n < l$ as inferred from Eq.~(\ref{eq2}).

To extend our theory into two-point cumulants, we use $\alpha$ and $\beta$ from Eq.~(\ref{eq4}). For their joint distribution \cite{my71, ks77}, the characteristic function is $\langle e^{is \alpha + it \beta} \rangle$. The cumulants $\langle \alpha^n \beta^m \rangle_c$ are obtained from
\begin{subequations} \label{eq7}
\begin{equation} \label{eq7a}
\ln \langle e^{is \alpha + it \beta} \rangle = \! \! \! \! \sum_{n+m=1}^{\infty} \! \! \! \langle \alpha^n \beta^m \rangle_c \frac{(is)^n}{n!} \frac{(it)^m}{m!}
\end{equation}
or from
\begin{equation} \label{eq7b}
 \langle \alpha^n \beta^m \rangle_c = \left. \frac{\partial^{n+m}}{\partial (is)^n \partial (it)^m} \ln \langle e^{is \alpha + it \beta} \rangle \right\vert_{s=t=0} .
\end{equation}
For consistency with Eq.~(\ref{eq5c}), we incorporate the separation of the two points $r$ as
\begin{equation} \label{eq7c}
\frac{\partial}{\partial z} \ln \langle e^{is \alpha + it \beta} \rangle = - \frac{\varphi (s,t,r/z)}{z} .
\end{equation}
Here $\varphi$ is a nondimensional function of the nondimensional parameters $s$, $t$, and $r/z$. They are independent of one another and also of the other parameter $z$. From Eqs.~(\ref{eq7b}) and (\ref{eq7c}),
\begin{align} \label{eq7d}
\frac{\partial}{\partial z} \langle \alpha^n \beta^m \rangle_c & = \left. - \frac{1}{z} \frac{\partial^{n+m}\varphi (s,t,r/z)}{\partial (is)^n\partial (it)^m } \right\vert_{s=t = 0} \nonumber \\
                                                                                & =         - \frac{\varphi^{(n,m,0)} (0,0,r/z)}{i^{n+m} z}.
\end{align}
\end{subequations}
We consider the cases of $n+m \ge 2$ and use the function $D_{u^nu^m}(r/z)$ in place of $\varphi^{(n,m,0)} (0,0,r/z)/i^{n+m}$. Because of $\langle \alpha^n \beta^m \rangle_c = \langle u^n(x,z)u^m(x+r,z) \rangle_c/u_{\tau}^{n+m}$,
\begin{equation} \label{eq8}
\frac{\partial}{\partial z} \frac{\langle u^n(x,z)u^m(x+r,z) \rangle_c}{u_{\tau}^{n+m}} = - \frac{D_{u^nu^m}(r/z)}{z} .
\end{equation}
The integration of Eq.~(\ref{eq8}) leads to the logarithmic law of Eq.~(\ref{eq3a}).

For an extension into any other quantity, a joint distribution between $\alpha$ for this quantity and $\beta$ for $(U+u)/u_{\tau}$ is imposed to satisfy Eq.~(\ref{eq7}). Since Eq.~(\ref{eq7}) is reduced to Eq.~(\ref{eq5}) if $m=0$, the quantity itself satisfies Eq.~(\ref{eq5}) and then Eq.~(\ref{eq7}). We thereby obtain its own scaling laws.

The laws for fluctuations of the spanwise velocity are logarithmic as in the case of the streamwise velocity $u$. For the wall-normal velocity $w$, since it is equal to $0$ at the wall \cite{t76}, we adopt $\phi \equiv \varphi \equiv 0$ in Eqs.~(\ref{eq5}) and (\ref{eq7}). Hence, the one-point cumulants are constants. The two-point cumulants are functions of $r/z$ alone. Both of them are independent of $\ln (\delta /z)$. The same laws have been derived from the attached-eddy hypothesis \cite{t76,m17}.

We also consider the local rate of the energy dissipation $\varepsilon$. Although the dissipation is due to the fluid viscosity $\nu$, its value does not affect statistics of $\varepsilon$ such as the average $\langle \varepsilon \rangle$ if the Reynolds number is high enough \cite{ll59, t35, igk09}. The rate  $\varepsilon$ is nondimensionalized as
\begin{equation} \label{eq9}
\alpha = \frac{\varepsilon (x,z)}{u_{\tau}^3/z} \ \ \mbox{and} \ \ \beta = \frac{\varepsilon (x+r,z)}{u_{\tau}^3/z}.
\end{equation}
The mean rate of the energy production $u_{\tau}^3/\kappa z$ is enhanced in the limit $z/\delta \rightarrow 0$ of the constant-flux sublayer. With respect to this, the dissipation rate $\varepsilon$ needs to remain finite. We accordingly adopt $\phi \equiv \varphi \equiv 0$ in Eqs.~(\ref{eq5}) and (\ref{eq7}). Via integration, Eq.~(\ref{eq5d}) at $n=1$ yields Eq. (\ref{eq1c}) if the integration constant $c_{\varepsilon}$ is assigned to be $1/\kappa$, while Eq.~(\ref{eq7d}) yields Eq.~(\ref{eq3b}).

Thus, without invoking the attached eddies, we have reproduced the functional forms of Eqs.~(\ref{eq1}) and (\ref{eq3}). Such a scaling law exists simply because the constant-flux sublayer has no parameter except for the distance $z$ and the friction velocity $u_{\tau}$. The laws are limited to cumulants. Only their laws are extended systematically from that for the mean velocity $U$. If allowed by the condition at the wall surface, the law becomes logarithmic. For the existence of these scaling laws, an attached eddy is unnecessary, albeit useful to discussing their functions $C_{u^nu^m}$, $D_{u^nu^m}$, and $C_{\varepsilon^n \varepsilon^m}$ (see Sec.~\ref{S5}).

\begin{figure*}
\begin{center}
\rotatebox{90}{
\resizebox{5.50cm}{!}{\includegraphics*[9.9cm,1.9cm][18.5cm,29.2cm]{fig1.ps}}
}
\caption{\label{f1} One-point statistics (a) $U(z)/u_{\tau}$, (b) $\langle -uw(z) \rangle_c/u_{\tau}^2 = \langle -uw(z) \rangle/u_{\tau}^2$, (c) $\langle u^2(z) \rangle_c/u_{\tau}^2 = \langle u^2(z) \rangle/u_{\tau}^2$, (d) $\langle u^4(z) \rangle_c/u_{\tau}^4 = [ \langle u^4(z) \rangle - 3 \langle u^2(z) \rangle^2]/u_{\tau}^4$, and (e) $\langle \varepsilon(z) \rangle_c/(u_{\tau}^3/z) = \langle \varepsilon(z) \rangle/(u_{\tau}^3/z)$ against $z/\delta_{99}$ for $U_{\infty} = 6$\,m\,s$^{-1}$ (circles) and $12$\,m\,s$^{-1}$ (squares shifted horizontally by two units). The filled symbols lie in the constant-flux sublayer. To these, solid lines are regression fits of Eq.~(\ref{eq1a}), (\ref{eq1b}), or (\ref{eq3a}). We provide $\pm 2 \sigma$ errors, albeit not including those for $u_{\tau}$ and $z-z_{\rm wt}$.}
\end{center}
\end{figure*}

\section{Experiments} \label{S3}

Experiments of turbulent boundary layers were done in a wind tunnel of the Meteorological Research Institute. We use coordinates $x_{\rm wt}$, $y_{\rm wt}$, and $z_{\rm wt}$ in the streamwise, spanwise, and floor-normal directions. Their origin $x_{\rm wt} = y_{\rm wt} = z_{\rm wt} = 0$ is on the center of the floor at the upstream end of the test section of the tunnel. Its size is ${\mit\Delta}x_{\rm wt} = 18$\,m, ${\mit\Delta}y_{\rm wt} = 3$\,m, and ${\mit\Delta}z_{\rm wt} = 2$\,m. The cross section ${\mit\Delta}y_{\rm wt} \times {\mit\Delta}z_{\rm wt}$ is the same upstream to $x_{\rm wt} = -4$\,m.

Upon the entire floor from $x_{\rm wt} = -4$\,m to $+18$\,m with an interval of ${\mit\Delta}x_{\rm wt} = 0.1$\,m, spanwise rods of diameter $3.0$\,mm were set as roughness. It displaces the zero plane of the wall turbulence $z=0$ from the floor surface $z_{\rm wt}=0$ \cite{my71}, for which we assume $z - z_{\rm wt} = -1.5 \pm 1.5$\,mm \cite{losda03}.

The incoming flow velocity $U_{\infty}$ was set at $6$ or $12$\,m\,s$^{-1}$. Over some range of distances $z_{\rm wt}$ at $x_{\rm wt} = +14$ m and at $y_{\rm wt} = 0$\,m, where the turbulence had been well developed and had become almost independent of the position $x_{\rm wt}$, we measured the streamwise velocity $U+u$.

We used a hot-wire anemometer made up of a constant-temperature system (Dantec, 90C10) and of a single-wire probe (Dantec, 55P04). The wire was of platinum-plated tungsten, $5$\,$\mu$m in diameter, $1.25$\,mm in sensing length, and oriented to the spanwise direction. Its resistance overheat ratio was set at $0.80$. We calibrated the anemometer before and after each series of the measurements.

The anemometer signal was low-pass filtered and then digitally sampled. For each pair of $U_{\infty}$ and $z_{\rm wt}$, the sampling frequency $f_s$ was set as high as possible, provided that noise was still negligible at around $f_s$ where the energy spectrum of the signal had decayed substantially \cite{po97, mhk07}. The filter cutoff was at $f_s/2$.

The total length of the data at each of the distances $z_{\rm wt}$ from $45$ to $105$\,mm was as large as $3.2 \times 10^8$ for $U_{\infty} = 6$ m\,s$^{-1}$ under $f_s = 16$\,kHz and $4.0 \times 10^8$ for $U_{\infty} = 12$\,m\,s$^{-1}$ under $f_s = 44$\,kHz.  At the other distances $z_{\rm wt}$, we individually obtained $4.0 \times 10^7$ data.

During these measurements, we monitored the flow conditions such as the temperature. They are used to estimate the fluid viscosity $\nu$.

From our data along time $t_{\rm wt}$, spatial information is obtained via Taylor's hypothesis of $x = -Ut_{\rm wt}$. Despite some known problems \cite{dn08, dj09}, we rely on this hypothesis up to a large separation of $r/z \gtrsim 10^2$. The reason is the value of $\langle u^2 \rangle /U^2$. It was small enough, i.e., $\lesssim 0.03$, at all the distances $z_{\rm wt}$.

The local rate $\varepsilon$ of the energy dissipation is obtained as $15\nu (\partial_x u )^2$ by assuming local isotropy of the turbulence. For the calculation of the derivative $\partial_x u$, we use the four-point finite difference. This is a surrogate of the true rate, but its result for $\langle \varepsilon (x, z) \varepsilon (x+r, z) \rangle_c$ is reliable except at smallest separations $r$ around and below the Kolmogorov length $\eta = (\nu^3 / \langle \varepsilon \rangle)^{1/4}$ \cite{cgs03}.

These calculations are made for individual segments of length $10^7$ of our data. Among segments, statistics exhibit scatters. They originate in variations of experimental conditions, calibration uncertainties, and incomplete convergence due to a limited sampling time. After removing segments that are too noisy for some uncertain reason, we use these scatters to estimate the final errors in a standard manner \cite{br03}.

Supplementary short measurements were also done. To estimate the boundary layer thickness $\delta_{99}$, i.e., a distance $z$ at which $U$ is $99$\% of its maximum, we measured $U$ with an interval of ${\mit\Delta}z_{\rm wt} = 10\,\mbox{mm} \simeq 0.03\delta_{99}$. Furthermore, to obtain the momentum flux $\rho \langle -uw \rangle$, we measured $u$ and $w$ by utilizing a crossed-wire probe of the anemometer (Dantec, 55P53). Its wires were $1$\,mm in separation and oriented at $\pm 45^{\circ}$ to the streamwise direction. The other settings were the same as for the single-wire probe. We estimate the errors as described above about the long measurements.

\section{Results} \label{S4}

Figure \ref{f1} shows one-point statistics semilogarithmically as a function of $z/\delta = z/\delta_{99}$. Their parameters are summarized in Table \ref{t1}. For these and other following results, $\pm 2 \sigma$ errors are given as typical uncertainties \cite{br03}.

The constant-flux sublayer is observed in Fig.~\ref{f1} at least from $z/\delta_{99} \simeq 0.14$ to $0.28$ (filled symbols). Throughout this range, $U$ in Fig.~\ref{f1}(a) is logarithmic, $\langle -uw \rangle_c =  \langle -uw \rangle$ in Fig.~\ref{f1}(b) is constant, and $\langle u^2 \rangle_c = \langle u^2 \rangle$ in Fig.~\ref{f1}(c) is logarithmic.

\begin{table}[tbp]
\begingroup
\squeezetable
\caption{\label{t1} Parameters for $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$: boundary layer thickness $\delta_{99}$, friction velocity $u_{\tau}$, viscosity $\nu$, Reynolds number $\delta_{99} u_{\tau}/\nu$, aerodynamic roughness $z_0$ and von K\'arm\'an constant $\kappa$ of $U(z)/u_{\tau} = \ln (z/z_0)/\kappa$ as well as $C_{u^n}(0)$ and $D_{u^n}(0)$ of Eq.~(\ref{eq3a}). The uncertainties are $\pm 2 \sigma$ errors. We also provide ranges of the Kolmogorov length $\eta$ and of the total sampling time $T$ among distances $z$ in the constant-flux sublayer.}
\begin{ruledtabular}
\begin{tabular}{lccc}
                             & Unit             & $U_{\infty} = 6$\,m\,s$^{-1}$ & $U_{\infty} = 12$\,m\,s$^{-1}$ \\  \hline
$\delta_{99}$                & mm               & $392   \pm  5$                & $403   \pm 5$                  \\
$u_{\tau} $                  & mm\,s$^{-1}$     & $259   \pm  1$                & $512   \pm 0$                  \\
$\nu$                        & mm$^2$\,s$^{-1}$ & $15.0  \pm  0.0$              & $14.7  \pm 0.0$                \\
$\delta_{99} u_{\tau}/\nu$   & $10^3$           & $6.74  \pm  0.08$             & $14.0  \pm 0.2$                \\
$\eta$                       & mm               & $0.27$ to  $0.29$             & $0.16$ to $0.17$               \\
$z_0$                        & mm               & $0.22  \pm  0.07$             & $0.19  \pm 0.07$               \\
$\kappa$                     &                  & $0.35  \pm  0.02$             & $0.36  \pm 0.02$               \\
$C_{u^2}(0)=C_{uu}(0)$       &                  & $2.52  \pm  0.15$             & $2.59  \pm 0.16$               \\
$D_{u^2}(0)=D_{uu}(0)$       &                  & $1.25  \pm  0.07$             & $1.19  \pm 0.08$               \\
$C_{u^4}(0)=C_{u^2u^2}(0)$   &                  & $2.08  \pm  1.18$             & $1.48  \pm 1.14$               \\
$D_{u^4}(0)=D_{u^2u^2}(0)$   &                  & $-5.01 \pm  0.48$             & $-4.47 \pm 0.46$               \\
$T U_{\infty} / \delta_{99}$ & $10^6$           & $0.28$ to  $0.31$             & $0.20$ to $0.24$               \\
\end{tabular}
\end{ruledtabular}
\endgroup
\end{table}
\begin{figure}[htbp]
\resizebox{7.84cm}{!}{\includegraphics*[3.0cm,1.7cm][17.5cm,28.4cm]{fig2.ps}}
\caption{\label{f2} Two-point velocity cumulant in the constant-flux sublayer $\langle u(x,z) u(x+r,z) \rangle_c/u_{\tau}^2 = \langle u(x,z) u(x+r,z) \rangle/u_{\tau}^2$ against $r/z$ for (a) $U_{\infty} = 6$\,m\,s$^{-1}$ and (b) $12$\,m\,s$^{-1}$. The solid arrow indicates an increase in $\delta_{99}/z$. The gray areas are examples of $\pm 2 \sigma$ errors, albeit not including those for $u_{\tau}$ and $z-z_{\rm wt}$. The panel (c) shows $C_{uu}(r/z)$ and $D_{uu}(r/z)$ of Eq. (\ref{eq3a}). For $U_{\infty} = 12$\,m\,s$^{-1}$, we provide $\pm 2 \sigma$ errors including those for $u_{\tau}$ but not for $z-z_{\rm wt}$. }
\end{figure}
\begin{figure}[bp]
\rotatebox{90}{
\resizebox{5.50cm}{!}{\includegraphics*[3.0cm,7.2cm][16.7cm,26.5cm]{fig3.ps}}
}
\caption{\label{f3} Two-point velocity cumulant $\langle u(x,z) u(x+r,z) \rangle_c/u_{\tau}^2$ $= \langle u(x,z) u(x+r,z) \rangle/u_{\tau}^2$ against $\ln (\delta_{99}/z )$ at $r/z = 0.1$, $0.3$, $1.0$, $3.0$, and $10.0$ for (a) $U_{\infty} = 6$\,m\,s$^{-1}$ and (b) $12$\,m\,s$^{-1}$. The filled symbols lie in the constant-flux sublayer. To these, solid lines are regression fits of Eq.~(\ref{eq3a}). Although we provide $\pm 2 \sigma$ errors on all the data in the same manner as in Figs.~\ref{f2}(a) and \ref{f2}(b), none of them are discernible. }
\end{figure}
\begin{figure}[htbp]
\resizebox{7.84cm}{!}{\includegraphics*[3.0cm,1.7cm][17.5cm,28.4cm]{fig4.ps}}
\caption{\label{f4} Same as in Fig.~\ref{f2} but for $\langle u^2(x,z) u^2(x+r,z) \rangle_c/u_{\tau}^4 = [\langle u^2(x,z) u^2(x+r,z) \rangle - \langle u(z)^2 \rangle^2 - 2 \langle u(x,z) u(x+r,z) \rangle^2]/u_{\tau}^4$. }
\end{figure}
\begin{figure}[tbp]
\rotatebox{90}{
\resizebox{5.50cm}{!}{\includegraphics*[3.0cm,7.2cm][16.7cm,26.5cm]{fig5.ps}}
}
\caption{\label{f5} Same as in Fig.~\ref{f3} but for $\langle u^2(x,z) u^2(x+r,z) \rangle_c/u_{\tau}^4 = [\langle u^2(x,z) u^2(x+r,z) \rangle - \langle u(z)^2 \rangle^2 - 2 \langle u(x,z) u(x+r,z) \rangle^2]/u_{\tau}^4$. }
\end{figure}
\begin{figure}[tbp]
\resizebox{7.84cm}{!}{\includegraphics*[3.0cm,1.7cm][17.5cm,28.4cm]{fig6.ps}}
\caption{\label{f6} Two-point cumulant of the rate of the energy dissipation in the constant-flux sublayer $\langle \varepsilon (x,z) \varepsilon (x+r,z) \rangle_c/(u_{\tau}^3/z)^2$ $= [\langle  \varepsilon (x,z) \varepsilon (x+r,z) \rangle - \langle \varepsilon (z) \rangle^2 ]/(u_{\tau}^3/z)^2$ against $r/z$ for (a) $U_{\infty} = 6$\,m\,s$^{-1}$ and (b) $12$\,m\,s$^{-1}$. The solid arrow indicates an increase in $\delta_{99}/z$. The panel (c) shows $C_{\varepsilon\varepsilon}(r/z)$ of Eq.~(\ref{eq3b}). We provide examples of $\pm 2 \sigma$ errors as in Figs.~\ref{f2} and \ref{f4}. The circles and squares indicate the Kolmogorov length $\eta$. The insets are for large values of $r/z$.}
\end{figure}

The inner bound of the sublayer depends on roughness of the wall \cite{t76, hvbs13}. Specifically for our roughness, we have observed its direct effect at $z \lesssim 40$\,mm, i.e., $\lesssim 0.10 \delta_{99}$, by shifting measurement positions slightly in the $x$ direction. The outer bound depends on the turbulence itself. Our estimate of $z/\delta_{99} \simeq 0.28$ is larger than those of $0.15$ in most studies \cite{mmhs13, hvbs13, vhs15, mm19}. Since the sublayer at any finite Reynolds number is an approximation (Sec.~\ref{S1}), no unique definition exists about its bound. An estimate similar to ours is actually found in the literature \cite{ofsbta17}.

We have used $\langle -uw \rangle^{1/2}$ as the friction velocity $u_{\tau}$. The results for the von K\'arm\'an constant $\kappa$ in Table~\ref{t1} are small with respect to the standard value of $\kappa = 0.39 \pm 0.02$ \cite{mmhs13}. Since it yet lies within $\pm 2\sigma$ errors of our results, we have adopted $\pm 2\sigma$ as a typical level of the uncertainties. Then, $c_{u^2} = C_{u^2}(0)$ and $d_{u^2} = D_{u^2}(0)$ in Table~\ref{t1} are consistent with those in Sec.~\ref{S1} summarized from the literature \cite{mmhs13, hvbs13, vhs15, ofsbta17, smhffhs18, hvbs12}.

The logarithmic scaling is also observed in Fig.~\ref{f1}(d) for $\langle u^4 \rangle_c = \langle u^4 \rangle - 3 \langle u^2 \rangle^2$. We expect this from Eq.~(\ref{eq3a}) at $r = 0$ and $n + m = 4$. The results for $C_{u^4}(0)$ and $D_{u^4}(0) $ in Table \ref{t1} are consistent within $\pm 2\sigma$ errors between the cases of $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$ (see also Appendix \ref{apA}).

We have $\langle u^4 \rangle_c = \langle u^4 \rangle - 3 \langle u^2 \rangle^2 < 0$ in the constant-flux sublayer. Its streamwise fluctuations $u$ are known to be sub-Gaussian, i.e., $\langle u^4 \rangle < 3 \langle u^2 \rangle^2$. The spanwise and the wall-normal fluctuations are super-Gaussian \cite{ff96}.

The present and some other data exhibit $\langle u^4 \rangle_c / \langle u^2 \rangle_c^2 \simeq -0.3$ \cite{hvbs13, vhs15, smhffhs18, ff96, mm13}. Since $\langle u^4 \rangle_c$ is not so significant, $\langle u^4 \rangle = \langle u^4 \rangle_c + 3 \langle u^2 \rangle_c^2$ could be approximated by $3\langle u^2 \rangle_c^2 = 3\langle u^2 \rangle^2$, to which Eq.~(\ref{eq1b}) for $\langle u^2 \rangle$ is applicable \cite{mm13}. Nevertheless, the exact law is Eq.~(\ref{eq3a}) for the cumulant $\langle u^4 \rangle_c$.

Finally in Fig.~\ref{f1}(e), the ratio of $\langle \varepsilon \rangle_c = \langle \varepsilon \rangle$ to $u_{\tau}^3/z$ lies at around $1/\kappa = 2.6$ (dotted lines), which corresponds to the standard value of $\kappa = 0.39$ \cite{mmhs13}. Although that ratio appears to vary in contrast to the law of Eq.~(\ref{eq1c}) or (\ref{eq3b}) \cite{lm15}, a further study is desired because $\langle \varepsilon (x, z) \varepsilon (x+r, z) \rangle_c$ satisfies Eq.~(\ref{eq3b}) in a range of separations $r$ (see below and also Appendix \ref{apB}).

To confirm that the spatial and temporal resolutions of our experiments were high enough for an estimation of $\varepsilon$, we consider $\langle (\partial_x u)^4 \rangle / \langle (\partial_x u)^2 \rangle^2 \simeq \langle \varepsilon^2 \rangle / \langle \varepsilon \rangle^2$. As for $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$, its values in the constant-flux sublayer are $7.7$--$7.8$ and $8.5$--$8.6$. Given the microscale Reynolds numbers $\langle u^2 \rangle / \nu  \langle (\partial_x u)^2 \rangle^{1/2}$ of $400$--$410$ and $570$--$580$, they are consistent with results of the previous studies \cite{sa97}.

We are now to study scaling laws of two-point statistics for Eq.~(\ref{eq3a}) at $n = m = 1$ and $2$ as well as for Eq.~(\ref{eq3b}) at $n = m = 1$. These statistics are shown in Figs.~\ref{f2}--\ref{f6} as a function of $r/z$ or $\delta/z = \delta_{99}/z$.

Figure \ref{f2} shows $\langle u(x,z) u(x+r,z) \rangle_c/u_{\tau}^2$ against $r/z$ at several wall-normal distances $z$ lying in the constant-flux sublayer. Here $\langle u(x,z) u(x+r,z) \rangle_c = \langle u(x,z) u(x+r,z) \rangle$ is a velocity correlation that would offer basic information about the wall turbulence. While it is usual to study this correlation as a function of $r/(\nu/u_{\tau})$ \cite{ghhlm05} or $r/\delta$ \cite{mswc07, dn11a}, we have adopted $r/z$ from Eq.~(\ref{eq3a}).

If $r/z$ is less than a few times $10^{-2}$, $\langle u(x,z) u(x+r,z) \rangle_c$ at each of $\delta_{99}/z$ is almost constant at about the value of $\langle u^2(z) \rangle_c$. With an increase in $r/z$, it decays toward $0$. It is no longer persistent if $r/z$ exceeds a few times $10^1$.

If $r/z$ is fixed, $\langle u(x,z) u(x+r,z) \rangle_c$ becomes large with an increase in $\delta_{99}/z$ (an arrow). This is due to the logarithmic law of Eq.~(\ref{eq3a}). Actually in Fig.~\ref{f3}, data points of $\langle u(x,z) u(x+r,z) \rangle_c$ make up linear functions of $\ln (\delta_{99}/z)$. The same has been observed in a laser Doppler anemometry of a similar flow \cite{mizuno18}. By fitting Eq.~(\ref{eq3a}) to our data \cite{br03}, we calculate the functions $C_{uu}(r/z)$ and $D_{uu}(r/z)$. Those in Fig.~\ref{f2}(c) for $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$ collapse individually to single curves,  the shapes of which are to be discussed in Sec.~\ref{S5}.

Figure \ref{f4} shows $\langle u^2(x,z) u^2(x+r,z) \rangle_c/u_{\tau}^4$ against $r/z$. Since $\langle u^2(x,z) u^2(x+r,z) \rangle_c = \langle u^2(x,z) u^2(x+r,z) \rangle - \langle u(z)^2 \rangle^2 - 2 \langle u(x,z) u(x+r,z) \rangle^2$ at $r=0$ is reduced to $\langle u^4(z) \rangle_c = \langle u^4(z) \rangle - 3 \langle u^2(z) \rangle^2$ in Fig.~\ref{f1}(d), this is a correlation of some non-Gaussian component of the velocity fluctuations $u$.

The two-point cumulant $\langle u^2(x,z) u^2(x+r,z) \rangle_c$ at each of $\delta_{99}/z$ is constant up to $r/z \simeq 10^{-1}$ and is persistent up to $r/z \simeq 10^1$. If $r/z$ is fixed, it is dependent on $\delta_{99}/z$. The reason is the logarithmic law of Eq.~(\ref{eq3a}) as confirmed in Fig.~\ref{f5}. Its functions $C_{u^2u^2}(r/z)$ and $D_{u^2u^2}(r/z)$ in Fig.~\ref{f4}(c) are consistent between the cases of $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$. These are analogous to results for $\langle u(x,z) u(x+r,z) \rangle_c$ in Figs.~\ref{f2} and \ref{f3}, although the shapes of the curves are entirely different.

Figure \ref{f6} shows $\langle \varepsilon (x,z) \varepsilon (x+r,z) \rangle_c/(u_{\tau}^3/z)^2$ against $r/z$. Here $\langle \varepsilon (x,z) \varepsilon (x+r,z) \rangle_c = \langle  \varepsilon (x,z) \varepsilon (x+r,z) \rangle - \langle \varepsilon (z) \rangle^2$ is a correlation of the dissipation rate $\varepsilon$. While it is usual to study this correlation as a function of $r/\eta$ \cite{cgs03, po97, mthk06}, we have adopted $r/z$ from Eq.~(\ref{eq3b}).

There is an enhancement of $\langle \varepsilon (x,z) \varepsilon (x+r,z) \rangle_c$ at the smallest separations $r$ \cite{cgs03, po97, mthk06}. They lie in the dissipative range, which extends from the Kolmogorov length $\eta$ (circles or squares) by a factor of $20$--$30$ (dotted arrows). Since the fluid viscosity $\nu$ is not negligible, Eq.~(\ref{eq3b}) does not hold there. The enhancement of $\langle \varepsilon (x,z) \varepsilon (x+r,z) \rangle_c$ is to be described rather by the small-scale intermittency \cite{igk09, sa97}.

Above that dissipative range, there exist the inertial and the energy-containing ranges. Throughout these two, $\langle \varepsilon (x,z) \varepsilon (x+r,z) \rangle_c/(u_{\tau}^3/z)^2$ is independent of $\delta_{99}/z$ in accordance with Eq.~(\ref{eq3b}). It persists up to a large separation \cite{cgs03, po97, mthk06}, i.e., $r/z \simeq 10^1$, as predicted originally by Landau \cite{ll59, mthk06}. We calculate the average at each of $r/z$ to obtain the function $C_{\varepsilon \varepsilon}(r/z)$. Those in Fig.~\ref{f6}(c) for $U_{\infty} = 6$ and $12$ m\,s$^{-1}$ collapse to a single curve at least at separations $r$ in the inertial and the energy-containing ranges.

\section{Discussion} \label{S5}

Having confirmed the scaling laws of Eq.~(\ref{eq3}), we discuss their functions $C_{uu}$, $D_{uu}$, $C_{u^2u^2}$, $D_{u^2u^2}$, and $C_{\varepsilon \varepsilon}$ in terms of the attached-eddy hypothesis \cite{t76,m17}.

\subsection{Scaling laws from attached eddies} \label{S5a}

Figure \ref{f7}(a) is a schematic of the attached eddies. They have various finite sizes, a common shape that is extending from the wall, and a common characteristic velocity $u_{\tau}$. If $\mbox{\boldmath{$x$}}_{\rm e} = (x_{\rm e}, y_{\rm e}, h_{\rm e})$ is the most extended position of an eddy, $h_e$ is used as its size. It induces the streamwise velocity $u_e$ and the energy dissipation $\varepsilon_e$ at any position $\mbox{\boldmath{$x$}} = (x, y, z)$ as
\begin{equation}
\label{eq10}
\frac{u_{\rm e}(\mbox{\boldmath{$x$}})}{u_{\tau}}                 = f_u             \! \left( \!  \frac{\mbox{\boldmath{$x$}} - \mbox{\boldmath{$x$}}_{\rm e}}{h_{\rm e}} \! \right) \!
\ \ \mbox{and} \ \
\frac{\varepsilon_{\rm e}(\mbox{\boldmath{$x$}})}{u_{\tau}^3/h_e} = f_{\varepsilon} \! \left(  \! \frac{\mbox{\boldmath{$x$}} - \mbox{\boldmath{$x$}}_{\rm e}}{h_{\rm e}} \! \right) \! .
\end{equation}
Here $f_u$ and $f_{\varepsilon}$ are nondimensional functions that take finite nonzero values in a finite volume. As for the velocity function $f_u$, we impose a free-slip wall condition, i.e., $f_u \ne 0$ at $z = 0$.

The number of eddies of size $h_e$ per unit area of the wall is $N_e h_e^{-3} dh_e$  \cite{t76}. Here $N_e$ is a constant. Apart from $N_e$ and $h_e$, no quantity affects such a number density. Upon the wall, the distribution of the eddies is random and independent. They could overlap one another because we use them not as realistic organized structures but only as bases to describe the flow \cite{t76, mm19}.

The entire flow is a superposition of the eddies. Since they are random and independent, a sum of their cumulants is identical to the cumulant of this flow. The law for $z/\delta \rightarrow 0$ is used as a law for the constant-flux sublayer, in accordance with its asymptotic character (Sec.~\ref{S1}, see also Appendix \ref{apB}).

From this attached-eddy hypothesis, Eq.~(\ref{eq3}) is derived by simplifying the original derivation \cite{m17}. We again use nondimensional functions of nondimensional parameters such as $r/z$.

The velocity cumulant $\langle u^n(x,z) u^m(x+r,z) \rangle_c$ at any distance $z \le \delta$ is given by an integration from $h_e = z$ to $\delta$,
\begin{subequations} \label{eq11}
\begin{equation} \label{eq11a}
\frac{\langle u^n(x,z)u^m(x+r,z) \rangle_c}{u_{\tau}^{n+m}} = \! N_e \! \! \int^{\delta}_z \! \! \frac{dh_e}{h_e} I_{u^nu^m} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right) \! .
\end{equation}
Here $I_{u^n u^m}$ is a cumulant of eddies of particular size $h_e$. As in the case of Eq.~(\ref{eq2}), it is related to moments $J_{u^l u^k}$ and $J_{u^l} = J_{u^l u^0}$. For example,
\begin{equation} \label{eq11b}
I_{uu} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right) \! = \! J_{uu} \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right) \! - J_{u}^2 \! \left( \! 0, \frac{z}{h_e} \! \right) \!.
\end{equation}
We use $\mbox{\boldmath{$r$}} = (r,0,0)$ to define $J_{u^lu^k}$ via an integration at that distance $z$,
\begin{equation} \label{eq11c}
J_{u^lu^k} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right)
\! = \! \!
\iint \! \frac{dx_e}{h_e} \frac{dy_e}{h_e} f_u^l \! \! \left( \! \frac{\mbox{\boldmath{$x$}}                         - \mbox{\boldmath{$x$}}_e}{h_e} \! \right) \!
                                           f_u^k \! \! \left( \! \frac{\mbox{\boldmath{$x$}} + \mbox{\boldmath{$r$}} - \mbox{\boldmath{$x$}}_e}{h_e} \! \right) \! \! .
\end{equation}
\end{subequations}
The first-order moment $J_u$ is not necessarily equal to $0$ \cite{pc82}. Although this was neglected in the original studies \cite{t76,m17}, their final results hold without any correction.

By using $\zeta = z/h_e$ and hence $d \zeta /\zeta = -dh_e /h_e$, we rewrite Eq.~(\ref{eq11a}) as a function of $z/\delta$ and $r/z$,
\begin{subequations} \label{eq12}
\begin{equation} \label{eq12a}
\! \frac{\langle u^n(x,z)u^m(x+r,z) \rangle_c}{u_{\tau}^{n+m}} = N_e \! \! \int^{1}_{z/\delta} \! \frac{d \zeta}{\zeta} I_{u^nu^m} \! \! \left( \frac{r}{z} \zeta, \zeta \right) \! . \! \!
\end{equation}
The free-slip wall condition implies that $I_{u^nu^m}(r \zeta /z, \zeta )$ at $ \zeta \rightarrow 0$ is not necessarily equal to $0$.  There is a function $\hat{I}^{(0)}_{u^nu^m}(r/z)$ such that $I_{u^nu^m} \rightarrow \hat{I}^{(0)}_{u^nu^m}$. By rewriting $I_{u^nu^m}$ as $I_{u^nu^m} - \hat{I}^{(0)}_{u^nu^m} + \hat{I}^{(0)}_{u^nu^m}$ in Eq.~(\ref{eq12a}) and by taking the limit $z/ \delta \rightarrow 0$, we derive Eq.~(\ref{eq3a}). Its functions $C_{u^nu^m}$ and $D_{u^nu^m}$ are
\begin{align}
\label{eq12b} &C_{u^nu^m} \! \! \left( \frac{r}{z} \right) \!
 = N_e \! \! \int^1_0 \! \frac{d \zeta}{\zeta} \! \left[ I_{u^n u^m} \! \! \left( \frac{r}{z} \zeta, \zeta \right) \! - \hat{I}^{(0)}_{u^nu^m} \! \! \left( \frac{r}{z} \right) \right] \! , \\
\label{eq12c} &D_{u^nu^m} \! \! \left( \frac{r}{z} \right) \!
 = N_e                                                                                \hat{I}^{(0)}_{u^nu^m} \! \! \left( \frac{r}{z} \right)         \! .
\end{align}
\end{subequations}
Because of $l \ge 1$ for $I_{u^nu^m} (r \zeta /z, \zeta ) - \hat{I}^{(0)}_{u^nu^m}(r/z) \varpropto \zeta^l$ at each of $r/z$ in the limit $\zeta \rightarrow 0$, the integration for $C_{u^nu^m}$ in Eq.~(\ref{eq12b}) is convergent. The divergent component has been removed as $D_{u^nu^m}$ in Eq.~(\ref{eq12c}). Likewise, at $n=1$ and $m=0$, Eq.~(\ref{eq12}) could yield the law of Eq.~(\ref{eq1a}) for the mean velocity $U(z)$ \cite{pc82},

The functional form of Eq.~(\ref{eq3a}) implies that it is always affected by the turbulence thickness $\delta$. Such an effect is through largest eddies of $h_e = \delta$, which contribute to all of nonzero cumulants \cite{m19}.

For the local rate of the energy dissipation $\varepsilon$, its cumulant is given by
\begin{subequations} \label{eq13}
\begin{equation} \label{eq13a}
\! \frac{\langle \varepsilon^n(x,z) \varepsilon^m(x+r,z) \rangle_c}{(u_{\tau}^3/z)^{n+m}} \!
=
N_e \! \! \int^{\delta}_z \! \! \frac{dh_e}{h_e} I_{\varepsilon^n \varepsilon^m} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right) \! . \! \!
\end{equation}
Here $I_{\varepsilon^n \varepsilon^m}$ is a cumulant of eddies of size $h_e$ that has been multiplied by a factor of $(z/h_e)^{n+m}$. The cumulant itself is related to moments $J_{\varepsilon^l \varepsilon^k}$ and  $J_{\varepsilon^l} = J_{\varepsilon^l \varepsilon^0}$. For example,
\begin{equation} \label{eq13b}
I_{\varepsilon \varepsilon} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right) \! = \frac{z^{2}}{h_e^{2}} \left[ J_{\varepsilon \varepsilon} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right) \!
                                                                                - J_{\varepsilon}^2           \! \! \left( \!             0, \frac{z}{h_e} \! \right) \! \right] .
\end{equation}
The moment $J_{\varepsilon^l \varepsilon^k}$ is defined as
\begin{equation} \label{eq13c}
J_{\varepsilon^l \varepsilon^k} \! \! \left( \! \frac{r}{h_e}, \frac{z}{h_e} \! \right)
\! = \! \!
\iint \! \frac{dx_e}{h_e} \frac{dy_e}{h_e} f_{\varepsilon}^l \! \! \left( \! \frac{\mbox{\boldmath{$x$}}                         - \mbox{\boldmath{$x$}}_e}{h_e} \! \right) \!
                                           f_{\varepsilon}^k \! \! \left( \! \frac{\mbox{\boldmath{$x$}} + \mbox{\boldmath{$r$}} - \mbox{\boldmath{$x$}}_e}{h_e} \! \right) \! \! .
\end{equation}
\end{subequations}
We use $\zeta = z/h_e$ to rewrite Eq.~(\ref{eq13a}),
\begin{subequations} \label{eq14}
\begin{equation} \label{eq14a}
\frac{\langle \varepsilon^n(x,z)\varepsilon^m(x+r,z) \rangle_c}{(u_{\tau}^3/z)^{n+m}}
 =
 N_e \! \! \int^{1}_{z/\delta} \! \frac{d \zeta}{\zeta} I_{\varepsilon^n \varepsilon^m} \! \! \left( \frac{r}{z} \zeta, \zeta \right) \! .
\end{equation}
Since $I_{\varepsilon^n \varepsilon^m}$ has a factor of $(z/h_e)^{n+m} = \zeta^{n+m}$ with $n+m \ge 1$ as in the case of Eq.~(\ref{eq13b}), the integration of Eq. (\ref{eq14a}) is convergent even in the limit $z/\delta \rightarrow 0$. Thus, Eq. (\ref{eq14a}) yields Eq.~(\ref{eq3b}). Its function $C_{\varepsilon^n \varepsilon^m}$ is
\begin{equation} \label{eq14b}
C_{\varepsilon^n \varepsilon^m} \! \! \left( \frac{r}{z} \right) = N_e \! \! \int^1_0 \! \frac{d \zeta}{\zeta}I_{\varepsilon^n \varepsilon^m} \! \! \left( \frac{r}{z} \zeta, \zeta \right) \! .
\end{equation}
\end{subequations}
That at $n=1$ and $m=0$ corresponds to the law of Eq. (\ref{eq1c}) for the average $\langle \varepsilon (z) \rangle$.

\begin{figure}[tbp]
\rotatebox{90}{
\resizebox{5cm}{!}{\includegraphics*[5.4cm,3.8cm][18.6cm,26.7cm]{fig7.ps}}
}
\caption{\label{f7} Schematic of (a) attached eddies, (b) $I_{uu}(r/h_e, z/h_e)$,  (c) $I_{u^2u^2}(r/h_e, z/h_e)$, and  (d) $I_{\varepsilon\varepsilon}(r/h_e, z/h_e)$. The solid arrow indicates an increase in $r/h_e$. The dotted lines are examples of $\hat{I}^{(0)}_{u^nu^n} (r/z)\varpropto D_{u^nu^n} (r/z)$ in Eq.~(\ref{eq12c}). They are to reproduce $D_{uu} \ge 0$ and $D_{u^2u^2} \le 0$ in Figs.~\ref{f2}(c) and \ref{f4}(c). As for integrations of $I_{u^nu^n}-\hat{I}^{(0)}_{u^nu^n}$ in Eq.~(\ref{eq12b}) and of $I_{\varepsilon\varepsilon}$ in Eq.~(\ref{eq14b}), the gray areas are examples to reproduce $C_{uu} (r/z)$, $C_{u^2u^2} (r/z)$, and $C_{\varepsilon \varepsilon} (r/z)$ in Figs. \ref{f2}(c), \ref{f4}(c), and \ref{f6}(c).}
\end{figure}

\subsection{Implication for attached eddies} \label{S5b}

By utilizing $D_{uu}(0) = D_{u^2}(0)$ and $D_{u^2u^2}(0) = D_{u^4}(0)$, we constrain the value of $N_e$. From Eq.~(\ref{eq12c}),
\begin{subequations} \label{eq15}
\begin{equation} \label{eq15a}
\frac{D_{u^4}(0)}{N_e} + \frac{2D_{u^2}^2(0)}{N_e^2} = \hat{I}^{(0)}_{u^4}(0) + 2 \! \left[ \hat{I}^{(0)}_{u^2}(0) \right]^2 \! \ge 0.
\end{equation}
This is because $\hat{I}^{(0)}_{u^4}(0)$ and $\hat{I}^{(0)}_{u^2}(0)$ are cumulants of orders $4$ and $2$, i.e.,  $\langle \alpha^4 \rangle_c + 2 \langle \alpha^2 \rangle_c^2 = \langle (\alpha - \langle \alpha \rangle)^4 \rangle - \langle (\alpha - \langle \alpha \rangle)^2 \rangle^2 = \langle [ (\alpha - \langle \alpha \rangle)^2 - \langle (\alpha - \langle \alpha \rangle)^2 \rangle ]^2 \rangle \ge 0$. Then,
\begin{equation} \label{eq15b}
N_e \le - \frac{2 D_{u^2}^2(0)}{D_{u^4}(0)}
\ \ \mbox{if} \ \
D_{u^4}(0) < 0 .
\end{equation}
\end{subequations}
The results for $D_{u^2}(0)$ and $D_{u^4}(0)$ in Table \ref{t1} yield $N_e \lesssim 0.6$. Here $N_e$ serves as a number of attached eddies of size $h_e$ per each volume of $h_e^3$. It is thereby inferred that the eddies are rather sparse so that those of similar sizes do not  overlap one another, although each of them contains many smaller eddies.

To explain $C_{uu}$, $D_{uu}$, $C_{u^2u^2}$, $D_{u^2u^2}$ and $C_{\varepsilon\varepsilon}$ observed in Figs. \ref{f2}(c), \ref{f4}(c), and \ref{f6}(c), we discuss the likely shapes of $I_{uu}$, $I_{u^2u^2}$, and $I_{\varepsilon\varepsilon}$. They are illustrated schematically in Figs. \ref{f7}(b)--\ref{f7}(d).

The functions $C_{uu}$, $C_{u^2u^2}$ and $C_{\varepsilon\varepsilon}$ in Eqs.~(\ref{eq12b}) and (\ref{eq14b}) are due to eddies of wall-normal sizes $h_e$ that are comparable to the observing distance $ z$ \cite{t76}. Since those in Figs.~\ref{f2}(c), \ref{f4}(c), and \ref{f6}(c) are persistent up to almost the same separation of $r/z \simeq 10^1$, the streamwise to wall-normal size ratio of an eddy is about $10^1$.

The functions  $D_{uu}$ and $D_{u^2u^2}$ in Eq.~(\ref{eq12c}) are due to wall-adjacent portions of the eddies \cite{t76}. Since those in Figs.~\ref{f2}(c) and \ref{f4}(c) are again persistent up to $r/z \simeq 10^1$, the streamwise size of such a portion is proportional to the distance $z$.

From the shapes of $C_{u^2u^2}$ and  $D_{u^2u^2}$ observed in Fig. \ref{f4}(c), we expect $I_{u^2u^2} < 0$ for all pairs of $r/h_e$ and $z/h_e$ (see Fig.~\ref{f7}). The sub-Gaussianity of the fluctuations $u$, i.e., $\langle u^4 \rangle_c < 0$ in Fig.~\ref{f1}(d), extends over the entire ranges of streamwise separations $r$ and of wall-normal distances $z$ in each of the eddies.

The function $C_{\varepsilon\varepsilon}$ in Fig.~\ref{f6}(c) differs between the cases of $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$ at $r \lesssim 30 \eta$ (dotted arrows). This dissipative range is where the attached-eddy hypothesis holds no longer \cite{m17, mm19}. From Eq.~(\ref{eq10}), the Kolmogorov length $\eta$ is obtained as $(\nu^3/\varepsilon_e)^{1/4} \varpropto (\nu /u_{\tau})^{3/4} h_e^{1/4}$, which is not proportional to the eddy size $h_e$. To describe motions in the dissipative range, we require some other class of small eddies, e.g., vortex tubes \cite{igk09, sa97, mhk07}, albeit not essential to statistics like $\langle u^n(x,z) u^m(x+r,z) \rangle_c$ that are almost constant there (Figs.~\ref{f2} and \ref{f4}).

\section{Concluding Remarks} \label{S6}

For the constant-flux sublayer of wall turbulence, the logarithmic scaling of $\langle u^n(x,z) u^m(x+r,z) \rangle_c$ in Eq.~(\ref{eq3a}) and the nonlogarithmic scaling of $\langle \varepsilon^n(x,z) \varepsilon^m(x+r,z) \rangle_c$ in Eq.~(\ref{eq3b}) have been studied experimentally. We have done experiments of boundary layers and have obtained those two-point cumulants at several distances $z$ from the wall (Sec.~\ref{S3}). The results in Figs.~\ref{f2}--\ref{f6} are consistent with Eq.~(\ref{eq3a}) at $n=m=1$ and $2$ as well as with Eq.~(\ref{eq3b}) at $n=m=1$ in the inertial and the energy-containing ranges of the separations $r$ (Sec.~\ref{S4}).

The mathematical reason for such a scaling law is that the constant-flux sublayer has only two local parameters, i.e., the distance $z$ and the friction velocity $u_{\tau}$. Under Galilean invariance, only laws for cumulants are extended systematically from the law of Eq.~(\ref{eq1a}) for the mean velocity $U(z)$. The logarithmic or nonlogarithmic character of the scaling is determined by a condition at the wall surface (Sec.~\ref{S2}). Even if those two points lie in the $y$ or $z$ direction \cite{m17}, we could rely on the same reasoning.

Having confirmed the scaling laws of Eq.~(\ref{eq3}), we have related them to the attached eddies of Townsend \cite{t76,m17}. Since the eddy number $N_e$ is rather small, eddies of each size $h_e$ do not overlap one another. The streamwise size of an eddy is about $10$ times its wall-normal size $h_e$. Within the individual eddies, the streamwise fluctuations $u$ are sub-Gaussian (Sec.~\ref{S5}).

These characteristics of the attached eddies might have been optimized to maximize the momentum flux $\rho \langle -uw \rangle$. Actually from recent applications of variational calculus to convective systems \cite{hcd14, mks18b}, it has been inferred that their scaling laws and eddy structures are optimal for their heat transfer. Such a study is equally desired for momentum transfer in wall turbulence.

We also desire to study cases other than the boundary layers. As noted in Sec.~\ref{S1}, among various configurations of the wall turbulence, $c_{u^2} = C_{u^2}(0)$ of Eq.~(\ref{eq1b}) is not common \cite{mmhs13, hvbs13, vhs15, ofsbta17, smhffhs18, hvbs12}. Then, $C_{u^nu^m}(r/z)$ of Eq.~(\ref{eq3a}) and possibly $C_{\varepsilon^n \varepsilon^m}(r/z)$ of Eq.~(\ref{eq3b}) are not universal. Nevertheless, since $d_{u^2} = D_{u^2}(0)$ is common, $D_{u^nu^m}(r/z)$ is expected to be universal. The reason would be that only $D_{u^nu^m}$ is determined locally in the constant-flux sublayer (Sec.~\ref{S2}).

Thus far, we have studied two-point cumulants, which are useful if wall turbulence is modeled as a superposition of finite-size motions like those of the attached eddies. Although it is usual to study energy spectra \cite{mm19}, their wavelengths do not correspond to any particular size of the individual motions \cite{m19}.

Finally, we remark on the turbulence thickness $\delta$. Wall turbulence is known to include long organized structures with streamwise lengths that exceed more than $10$ times the thickness $\delta$ \cite{mswc07, dn11b}. The ratio $r/\delta$ has accordingly been used as a parameter of two-point statistics \cite{mswc07, dn11a}. However, in our results and in some others \cite{m19, mizuno18}, those structures are not discernible. They are meandering so that their total lengths do not appear in one-dimensional correlations and spectra. Furthermore, the structures lie essentially away from the wall and do not affect any law of the constant-flux sublayer at least at $z/\delta \rightarrow 0$ \cite{mm19}. We have considered such asymptotic laws alone (see also Appendix \ref{apB}). Given the limit $z/\delta \rightarrow 0$, since the distance $z$ is always not equal to $0$, $r/\delta$ does not deserve to be a parameter. The true nondimensional parameters are $r/z$ and $\delta/z$ as demonstrated here for cases of the scaling laws of Eq.~(\ref{eq3}).

 \vspace{-3mm}

\begin{acknowledgments}
This work was supported in part by KAKENHI Grants No. 17K00526 and No. 19K03967.
\end{acknowledgments}

\begin{figure*}
\begin{center}
\rotatebox{90}{
\resizebox{5.50cm}{!}{\includegraphics*[9.9cm,1.9cm][18.5cm,29.2cm]{fig8.ps}}
}
\caption{\label{f8} Same as in Fig.~\ref{f1} but for (a) $\langle u^3(z) \rangle_c = \langle u^3(z) \rangle$, (b) $\langle u^5(z) \rangle_c = \langle u^5(z) \rangle -10 \langle u^3(z) \rangle \langle u^2(z) \rangle$, (c) $\langle u^6(z) \rangle_c = \langle u^6(z) \rangle -15 \langle u^4(z) \rangle \langle u^2(z) \rangle -10 \langle u^3(z) \rangle^2 +30 \langle u^2(z) \rangle^3$, (d) $\langle u^7(z) \rangle_c = \langle u^7(z) \rangle -21 \langle u^5(z) \rangle \langle u^2(z) \rangle -35 \langle u^4(z) \rangle \langle u^3(z) \rangle +210 \langle u^3(z) \rangle \langle u^2(z) \rangle^2$, and (e) $\langle u^8(z) \rangle_c = \langle u^8(z) \rangle -28 \langle u^6(z) \rangle \langle u^2(z) \rangle -56 \langle u^5(z) \rangle \langle u^3(z) \rangle - 35 \langle u^4(z) \rangle^2 +420 \langle u^4(z) \rangle \langle u^2(z) \rangle^2 +560 \langle u^3(z) \rangle^2 \langle u^2(z) \rangle -630 \langle u^2(z) \rangle^4$ nondimensionalized with use of $u_{\tau}$. The solid lines are regression fits of Eq.~(\ref{eq3a}).
}
\end{center}
\end{figure*}
\begin{table}[tbp]
\begingroup
\squeezetable
\caption{\label{t2} Parameters $C_{u^n}(0)$ and $D_{u^n}(0)$ of Eq.~(\ref{eq3a}) for $U_{\infty} = 6$ and $12$\,m\,s$^{-1}$. The uncertainties are $\pm 2 \sigma$ errors.}
\begin{ruledtabular}
\begin{tabular}{lcc}
                             & $U_{\infty} = 6$\,m\,s$^{-1}$ & $U_{\infty} = 12$\,m\,s$^{-1}$ \\  \hline
$C_{u^3}(0)$                 & $-2.17 \pm 0.33$              & $-1.99 \pm 0.34$               \\
$D_{u^3}(0)$                 & $0.64  \pm 0.13$              & $0.60  \pm 0.13$               \\
$C_{u^5}(0)$                 & $7.35  \pm 4.13$              & $5.35  \pm 5.11$               \\
$D_{u^5}(0)$                 & $3.75  \pm 1.58$              & $4.04  \pm 1.97$               \\
$C_{u^6}(0) = C_{u^3u^3}(0)$ & $-74.4 \pm 24.6$              & $-63.4 \pm 23.5$               \\
$D_{u^6}(0) = D_{u^3u^3}(0)$ & $79.5  \pm 9.8$               & $71.7  \pm 9.2$                \\
$C_{u^7}(0)$                 & $82    \pm 117$               & $122   \pm 151$                \\
$D_{u^7}(0)$                 & $-255  \pm 46$                & $-255  \pm 59$                 \\
$C_{u^8}(0) = C_{u^4u^4}(0)$ & $2790  \pm 950$               & $2460  \pm 910$                \\
$D_{u^8}(0) = D_{u^4u^4}(0)$ & $-2260 \pm 380$               & $-2070 \pm 350$                \\
\end{tabular}
\end{ruledtabular}
\endgroup
\end{table}