# Two-time Lagrangian velocity correlation function for particle pairs in two-
dimensional inverse energy-cascade turbulence

Tatsuro Kishi<EMAIL_ADDRESS>Takeshi Matsumoto
Sadayoshi Toh  Division of Physics and Astronomy, Graduate School of Science,
Kyoto University, Kitashirakawa Oiwaketyo Sakyoku, Kyoto 606-8502, Japan

###### Abstract

We numerically investigate a two-time Lagrangian velocity correlation function
(TTLVCF) for particle pairs in two-dimensional energy inverse-cascade
turbulence. We consider self similarity of the correlation function by means
of incomplete similarity. In this framework, we propose a self-similar form of
the correlation function, whose scaling exponents cannot be determined by only
using the dimensional analysis based on the Kolmogorov’s phenomenology. As a
result, the scaling laws of the correlation function can depend on the initial
separation. This initial-separation dependency is frequently observed in
laboratory experiments and direct numerical simulations of the relative
dispersion, which is directly related to the correlation function, at moderate
Reynolds numbers. We numerically verify the self-similar form by direct
numerical simulations of two-dimensional energy inverse-cascade turbulence.
The involved scaling exponents and the dependencies on finite Reynolds-number
effects are determined empirically. Then, we consider implication of the
scaling laws of the correlation function on the relative dispersion, i.e. the
Richardson–Obukhov $t^{3}$ law. Our results suggest a possibility not to
recover the Richardson–Obukhov $t^{3}$ law at infinite Reynolds number.

††preprint: APS/123-QED

##  I Introduction

Velocity correlation function is fundamental to characterize turbulence. We
can understand dynamical couplings between two points and two times in
turbulence through the correlation function [1]. Eulerian and Lagrangian
velocity correlation functions have different characteristic time scales from
each other, which are essential to develop two-point closure approximations to
the Navier-Stokes equations without ad-hoc parameters [2, 3, 4]. These direct-
interaction approximations or Lagrangian renormalized approximations provide
consistent results with Kolmogorov’s 1941 (K41) phenomenology [5] and also
with the 2D analog [6, 7, 8]. In particular, the success of these closures
lies in applying the approximation not to the Eulerian velocity correlation
but to the Lagrangian velocity correlation.

The most general form of the second-order Lagrangian velocity correlation
function can be written in terms of Kraichnan’s generalized velocity notation
[3] as

$Q^{L}_{ij}(\bm{a},s_{1}|t_{1};\bm{b},s_{2}|t_{2})\equiv\langle
v_{i}(\bm{a},s_{1}|t_{1})v_{j}(\bm{b},s_{2}|t_{2})\rangle.$ (1)

Here $v_{i}(\bm{a},s|t)$ is the $i$-th component of the Lagrangian velocity
measured at time $t$ of a Lagrangian particle passed through a point $\bm{a}$
at time $s$ and $\langle\cdot\rangle$ denotes an ensemble average. Since the
general form, Eq.(1), is too hard to tackle, a majority of theoretical,
numerical and experimental investigations on the Lagrangian velocity
correlation are devoted to the abridged form of Eq.(1) by setting
$s_{1}=s_{2}$ and $t_{2}=s_{2}$ [9, 10, 11, 12, 13, 14], namely,

$Q^{L}_{ij}(\bm{a},s|t;\bm{b},s|s)=\langle
v_{i}(\bm{a},s|t)v_{j}(\bm{b},s|s)\rangle.$ (2)

It should be noticed that $v_{j}(\bm{b},s|s)$ coincides with the Eulerian
velocity at a point $\bm{b}$ at time $s$. Thus, the correlation Eq.(2) is
between the Lagrangian velocity and the Eulerian velocity.

An abridged form of Eq.(1) but involving only the Lagrangian velocity can be

$Q^{L}_{ij}(\bm{a},s|t_{1};\bm{b},s|t_{2})=\langle
v_{i}(\bm{a},s|t_{1})v_{j}(\bm{b},s|t_{2})\rangle,$ (3)

where the measuring times $t_{1}$ and $t_{2}$ are different from the labeling
time $s$. There are few studies on the Lagrangian correlation function Eq.(3)
that is an essential ingredient to solve unsteady problems of turbulence such
as turbulent diffusion and mixing [15]. In a notable study of the correlation
Eq.(3) [16], the authors performed a direct numerical simulation (DNS) to
obtain the Lagrangian velocity correlation function Eq.(3). However, to
analyze the computed correlation function, they had to resort to the theory
developed for the Lagrangian-Eulerian correlation Eq.(2). This may be not only
because a theory for the Lagrangian correlation Eq.(3) is not developed, but
also because its simple characterization remains to be done. By the simple
characterization, we mean answers to the following cascading questions: does
the Lagrangian correlation function Eq.(3) have a self-similar form?; if this
is so, what is the self-similar form?; if the self-similar form is a power-law
function, what are scaling exponents? In this paper, we address these
questions with phenomenological theory beyond the dimensional analysis and
direct numerical simulations. Our simulation here is limited to two-
dimensional (2D) energy inverse cascade turbulence, but the theory is
applicable both to two and three dimensions.

One of the difficulties in these questions is that the Lagrangian velocity
correlation function Eq.(3) is intrinsically dependent on both two times,
$t_{1}$ and $t_{2}$. To illustrate our approach, let us show a color map of
the Lagrangian correlation function as a function of $t_{1}$ and $t_{2}$ in
Fig.1. For reasons described shortly below, we do not consider the Lagrangian
correlation Eq.(3). Instead, we study the Lagrangian velocity increment or,
equivalently, the relative velocity between two Lagrangian particles whose
Lagrangian labels are $\bm{a}$ and $\bm{a}+\bm{r}_{0}$ at time $s$,

$\delta v_{i}(\bm{a},\bm{r}_{0},s|t)\equiv
v_{i}(\bm{a}+\bm{r_{0}},s|t)-v_{i}(\bm{a},s|t)$ (4)

and its correlation

$C^{L}_{ij}(r_{0},t_{1},t_{2})=\langle\delta
v_{i}(\bm{a},\bm{r}_{0},s|t)\delta v_{j}(\bm{a},\bm{r}_{0},s|t)\rangle,$ (5)

where $r_{0}=|\bm{r}_{0}|$. We call the Lagrangian correlation Eq.(5) two-time
Lagrangian velocity correlation function (TTLVCF) in this paper (on the left
hand side (lhs) of Eq.(5), we omit the dependence on $\bm{a},s$ for
simplicity). The TTLVCF shown in Fig.1 is numerically computed in 2D energy
inverse-cascade turbulence. The details will be explained in Sec.III in this
paper.

To characterize the TTLVCF shown in Fig.1, one way is to look at it along the
diagonal line through the origin, $t_{1}=t_{2}$. The line is parallel to the
“$T$” axis written in Fig.1. The other way is obviously to study it along the
lines perpendicular to the diagonal line, i.e, $t_{1}+t_{2}=c_{p}$, where
$c_{p}$ is a positive constant. These lines are parallel to the “$\tau$” axis
written in Fig.1. Accordingly, a correlation time can be defined for each
line. From Fig.1, it can be observed that a correlation time along a
perpendicular line grows as we increase the constant $c_{p}$. This sort of
non-stationary behavior is not present in the two-point Eulerian velocity
correlation, whose correlation time is constant due to the statistical
stationarity. Unlike the Eulerian one, the TTLVCF has more than one degree of
freedom. It implies that the scaling law of the TTLVCF cannot be obtained by
dimensional analysis. For this kind of problems, the incomplete similarity
[17] provides a framework to specify possible self-similar forms. In this
study, by using both the incomplete similarity and DNS, we propose a self-
similar form of the TTLVCF shown in Fig.1

Figure 1:  Color map of the TTLVCF, Eq.(5), with $T_{B}=3.4T_{\eta}$ for
$\mathrm{Re}_{\alpha}=160$, where $T_{B}$ and $T_{\eta}$ are the Batchelor
time and the Kolmogorov dissipation time scale, respectively. The infrared
Reynolds number $\mathrm{Re}_{\alpha}$ is defined in the main text. See Sec.
III for details.

Now let us explain why we consider the correlation of the Lagrangian velocity
increment, Eq.(4). The TTLVCF is directly related to the relative dispersion
as follows. The relative separation of a particle pair, whose Lagrangian
labels are $\bm{a}$ and $\bm{r}_{0}$ at time $0$, is written as

$\bm{r}(t)=\bm{r}_{0}+\int_{0}^{t}\delta\bm{v}(t^{\prime})~{}dt^{\prime}.$ (6)

The squared separation can then be written in terms of the TTLVCF as

$\langle
r^{2}(t)\rangle=r_{0}^{2}+2\int_{0}^{t}\bm{r_{0}}\cdot\langle\delta\bm{v}(t_{1})\rangle~{}dt_{1}+\int_{0}^{t}\int_{0}^{t}\langle\delta\bm{v}(t_{1})\cdot\delta\bm{v}(t_{2})\rangle~{}dt_{1}dt_{2},$
(7)

where $r=|\bm{r}|$. The turbulent relative dispersion has been widely
investigated since the pioneering work by Richardson [18], who first predicted
that the lhs of Eq.(7) grows as $\langle r^{2}(t)\rangle\propto t^{3}$ over an
appropriate time interval, see Ref.[19] as a review. The $t^{3}$ prediction
can be also derived by naively using the K41 dimensional analysis in the
inertial subrange [20, 21] and therefore, the $t^{3}$ law is referred to as
the Richardson–Obukhov law.

In the K41 framework, statistics on the relative separation, whose best
studied quantity is the second moment $\langle r^{2}(t)\rangle$, is considered
to be independent of the initial separation, $r_{0}$, as long as $r(t)$ is in
the sufficiently wide inertial subrange. This implies that, when we plot
$\langle r^{2}(t)\rangle$ starting from various $r_{0}$’s, the curves become
independent of $r_{0}$ and eventually collapses to the only one curve
independent of $r_{0}$, which is proportional to $t^{3}$. A tendency toward
such an asymptotic state is indicated by DNS in three dimensions at high
Reynolds numbers [22, 23, 24]. By contrast, one has never clearly observed the
$r_{0}$-independence in laboratory experiments in three dimensions [25, 26]
and in two dimensions [27, 28, 29, 30]. The situation in numerical simulations
of the 2D inverse energy-cascade turbulence [31, 32, 33] is similar to that of
the 3D turbulence [34, 35, 22, 36, 23]. In this sense, the Richardson–Obukhov
law is not verified satisfactorily by observations. Of course, with much
higher Reynolds number, a cleaner $t^{3}$ law irrespective of $r_{0}$ may be
observed. In this paper, we parametrize finite Reynolds-number effect on the
Richardson–Obukhov law by analyzing the TTLVCF Eq.(5). Specifically, we will
study Reynolds-number dependence of the correlation function. Then, through
Eq.(7), we argue that $\langle r^{2}(t)\rangle$ depends on $r_{0}$ at finite
Reynolds numbers and infer its asymptotic form at infinite Reynolds number.

In particular, at moderate Reynolds numbers it is known that the $t^{3}$
scaling behavior of $\langle r^{2}(t)\rangle$ is observed for a certain
selected initial separation, see [31, 37, 33] for further discussion. This
special initial separation is around the Kolmogorov dissipation length for
both 2D and 3D and thus dependent on the Reynolds number. Our argument on the
$r_{0}$-dependence reveals the nature of this special initial separation.

We organize the paper as follows. In Sec. II, we make ansatz for the scaling
laws of the TTLVCF shown in Fig.1 by adapting the incomplete similarity [17].
The method enables us to explore qualitatively scaling laws that deviate from
the K41 dimensional analysis. In particular, we take the finite Reynolds
number effect and the dependence on $r_{0}$ into consideration. It should be
noticed that this method does not depend on specific dimensions. We can apply
the method to both 2D and 3D turbulence.

Next, in Sec.III, we verify the ansatz and determine quantitatively the
involved scaling exponents by comparing with DNS data of the 2D inverse
energy-cascade turbulence. We use a modified Navier-Stokes equations which
have a superviscosity and a hypodrag. We have already estimated the artificial
effects from these numerical terms and verified that these effects can be
ignored for the statistics of particle pairs in turbulence in Ref. [33]. We
estimate the scaling exponents as a function of Reynolds number by using DNSs
with four different Reynolds numbers. Subsequently, we infer asymptotic values
of the scaling exponents at infinite Reynolds number by extrapolating them
from those at finite Reynolds numbers. There are several reasons for selecting
the 2D system; detailed numerical studies are more feasible; the Eulerian
velocity is intermittency free [38, 39] and therefore we factor out the
intermittency effects on the Lagrangian statistics. Of course, careful
discussion and further investigation are required when one applies the method
used here to the 3D system. We discuss a possibility to justify its validity
to 3D turbulence in Sec. V.

In Sec.IV, we discuss implications of the self-similar ansatz of the TTLVCF on
the Richardson–Obukhov law. Furthermore, we explain why the $t^{3}$ scaling
law of $\langle r^{2}(t)\rangle$ is observed for a selected initial separation
at moderate Reynolds numbers. Concluding remarks are made in Sec.V.

##  II Incomplete self-similarity of the Lagrangian correlation and scaling
exponents

In this section, we present a scaling ansatz of the TTLVCF for particle pairs.
It should be noted that the argument below is independent of specific spatial
dimensions. Therefore, we expect that the ansatz is meaningful for both 2D and
3D turbulence.

As we discussed in Sec.I, we consider the Lagrangian correlation function,
$C^{L}_{ii}(r_{0},t_{1},t_{2})$, of the velocity difference Eq.(4), where we
use the Einstein summation convention of the repeated indices. We write
$C^{L}_{ii}(r_{0},t_{1},t_{2})$ as $C^{L}(r_{0},t_{1},t_{2})$ below and term
it the TTLVCF as well.

We consider a statistically steady, homogeneous, and isotropic turbulent
Eulerian velocity field, and therefore we deal with an external forcing which
leads to such a statistical state. In the following scaling argument, we
ignore effects of the external forcing on the TTLVCF. Later, in Sec. V, we
discuss the effects when we analyze the DNS data of the TTLVCF.

Let us first change time variables from $t_{1}$ and $t_{2}$ to the average
time, $T$, and the relative time, $\tau$ as

$T\equiv\frac{t_{1}+t_{2}}{2},\quad\tau\equiv t_{1}-t_{2}.$ (8)

The new variables are useful because the TTLVCF is symmetric with respect to
the $T$-axis as shown in Fig.1.

Now we present the scaling ansatz for the TTLVCF. Our arguments are given
thereafter. The TTLVCF can be this form:

$C^{L}(r_{0},T,\tau)=\varepsilon T\Phi(r_{0},T,\tau),$ (9)

which we assume to be valid at an appropriate time interval. Here
$\varepsilon$ is the energy dissipation rate. The non-dimensional function
$\Phi(r_{0},\tau)$ includes the deviation from the dimensional analysis. By
using the idea of incomplete similarity [17], we argue that
$\Phi(r_{0},T,\tau)$ can be written as

$\Phi(r_{0},T,\tau)=G\left(\frac{T_{B}}{T}\right)^{\gamma}g^{L}\left(\frac{\tau}{T_{B}^{\beta}T^{1-\beta}}\right),$
(10)

where $G$ is a non-dimensional, non-zero constant and $g^{L}(X)$ is a non-
dimensional self-similar function with $g^{L}(0)=1$. Here the two scaling
exponents, $\beta$ and $\gamma$, appear. They are not determined by
dimensional analysis. With the ansatz Eq.(10), the width of the ridge along
the diagonal line $t_{1}=t_{2}$ shown in Fig.1 is given by
$T(T_{B}/T)^{\beta}$. Later in Sec.III, by comparing with DNS data, we will
show that the exponents take the following functional forms

$\displaystyle\beta$
$\displaystyle=\beta_{0}+\tilde{\beta}\left(\frac{T_{\eta}}{T_{B}},\frac{T_{L}}{T_{B}}\right),$
(11) $\displaystyle\gamma$
$\displaystyle=\gamma_{0}+\tilde{\gamma}\left(\frac{T_{\eta}}{T_{B}},\frac{T_{L}}{T_{B}}\right).$
(12)

Here, we introduce three time scales, $T_{B}$, $T_{\eta}$, and $T_{L}$, which
are given by

$T_{B}\equiv\left(\frac{r_{0}^{2}}{\varepsilon}\right)^{1/3},T_{\eta}\equiv\left(\frac{\eta^{2}}{\varepsilon}\right)^{1/3},T_{L}\equiv\left(\frac{L^{2}}{\varepsilon}\right)^{1/3},$
(13)

where $T_{B}$ is the Batchelor time associated with the initial separation
$r_{0}$ [21], $T_{\eta}$ is the smallest time scale of turbulence associated
with the smallest length scale, $\eta$, such as Kolmogorov dissipation length,
and $T_{L}$ is the largest time scale of turbulence associated with the
largest length scale, $L$, such as the integral scale. Moreover, $\beta_{0}$
and $\gamma_{0}$ are the asymptotic exponents at infinite Reynolds number, and
therefore their values are independent of $T_{\eta}$ and $T_{L}$.

Now let us explain how we reach the scaling ansatz, (9)–(12), of the TTLVCF.
Our argument here follows the self-similar analysis of Ref.[17]. First, we
specify the system of all the governing parameters of the correlation
function, $C^{L}(r_{0},T,\tau)$. It depends on the average time, $T$, the
relative time, $\tau$, the initial separation of particle pairs, $r_{0}$, the
energy dissipation rate or the average energy flux in the inertial subrange,
$\varepsilon$, the smallest length scale of turbulence such as the Kolmogorov
length, $\eta$, and the largest length scale of turbulence such as the
integral scale, $L$. Taking them into account, we rewrite the arguments of the
TTLVCF as

$C^{L}(r_{0},T,\tau,\varepsilon,\eta,L)=\langle\delta
v_{i}(\bm{a},\bm{r}_{0},s=0|t_{1})\delta
v_{i}(\bm{a},\bm{r}_{0},s=0|t_{2})\rangle,$ (14)

Here, we take average over the Lagrangian label $\bm{a}$. Hence we omit the
dependence on $\bm{a}$. We set the labeling time to zero, i.e., $s=0$ and,
$t_{1}$ and $t_{2}$ are measured from this time origin. We also omit
dependence on $s$ on the lhs of Eq.(14). For the 2D inverse energy-cascade
turbulence, we can use the characteristic length of the drag as $L$ instead of
the integral scale and use the energy flux cascading inversely in the inertial
subrange as $\varepsilon$. In this case, we can explicitly write down $L$
dimensionally by using the drag coefficient and the energy flux. This may be
an advantage of the 2D inverse energy-cascade turbulence.

Second, we apply the Buckingham $\Pi$-theorem [17] to Eq.(14) by assuming that
the independent dimensions are $\varepsilon$ and $T$. This means that all the
other governing parameters are non-dimensionalized by $\varepsilon$ and $T$.
This leads to an expression with the dimensionless function, $C_{*}^{L}$, as

$C^{L}(r_{0},T,\tau,\varepsilon,\eta,L)=\varepsilon
T~{}C^{L}_{*}\left(\frac{\tau}{T},\frac{r_{0}}{\varepsilon^{1/2}T^{3/2}},\frac{\eta}{\varepsilon^{1/2}T^{3/2}},\frac{L}{\varepsilon^{1/2}T^{3/2}}\right).$
(15)

This is rewritten by using the time scales (13) as

$C^{L}(r_{0},T,\tau,\varepsilon,\eta,L)=\varepsilon
TC^{L}_{*}\left(\frac{\tau}{T},\left(\frac{T_{B}}{T}\right)^{3/2},\left(\frac{T_{\eta}}{T}\right)^{3/2},\left(\frac{T_{L}}{T}\right)^{3/2}\right).$
(16)

Third, we consider intermediate asympotics of the time scales, $T_{\eta}$,
$T$, and $T_{L}$, and reduce the number of the arguments on the rhs of
Eq.(16). As is clear from our choice of the independent dimensions, we assume
that $T_{\eta}$ and $T_{L}$ are sufficiently separated and that

$T_{\eta}\ll T\ll T_{L}.$ (17)

This intermediate time interval for $T$ is a Lagrangian counterpart of the
inertial subrange of the Eulerian velocity statistics. We call this time
interval the inertial subrange in this paper. In this inertial subrange, we
assume that the correlation function becomes independent of $T_{\eta}$ and
$T_{L}$. This implies that the complete similarity holds for $T_{\eta}$ and
$T_{L}$. Hence, with a dimensionless function $C^{L}_{**}(\zeta,\xi)$, Eq.(16)
is simplified as

$C^{L}(r_{0},T,\tau,\varepsilon)=\varepsilon
TC^{L}_{**}\left(\frac{\tau}{T},\frac{T_{B}}{T}\right).$ (18)

Fourth, let us also assume that the initial separation, $r_{0}$, is
sufficiently small. Namely, we consider that the average time is much larger
than $T_{B}~{}(\propto r_{0}^{2/3})$:

$(T_{\eta}\ll)\quad T_{B}\ll T\quad(\ll T_{L})$ (19)

To discuss behavior of $C^{L}$ in this time range, for simplicity, we rewrite
the dimensionless times as

$\zeta\equiv\frac{\tau}{T},\quad\xi\equiv\frac{T_{B}}{T}.$ (20)

The additional asymptotics (19) implies $\xi\to 0$. Now there are two
possibilities for the asymptotic behavior of $C^{L}_{**}(\zeta,\xi)$ as
$\xi\to 0$ [17]:

  1. (i)

the limit of $C^{L}_{**}(\zeta,0)$ exists and is finite and non-zero,

  2. (ii)

no finite limit of $C^{L}_{**}(\zeta,0)$ exists, or the limit is zero if it
exists.

We do not know a priori which case holds unless the full functional dependence
of $C^{L}(r_{0},T,\tau,\varepsilon)$ was obtained theoretically from the
Navier-Stokes equations. It is necessary to study data for small $\xi$, which
is obtained from DNS or laboratory experiment in order to conclude which case
is valid [17].

Now let us discuss implications of each case. In the case (i), a scaling law
for $C^{L}(\zeta,0)$ is consistent with the K41 phenomenology. In other words,
we can determine the scaling relations for $C^{L}(\zeta,0)$ by dimensional
analysis: complete similarity. In this case, we can estimate $\xi$-dependence
by the Taylor series

$C^{L}_{**}(\zeta,\xi)=C^{L}_{**}(\zeta,0)+A_{1}(\zeta)\xi+O(\xi^{2}),$ (21)

where $A_{1}(\zeta)=(\partial C^{L}_{**}/\partial\xi)_{\xi=0}$. Therefore, in
the case (i), the scaling law for $C^{L}(r_{0},T,\tau)$ is as follows:

$C^{L}(r_{0},T,\tau)=\varepsilon TA_{0}\left(\frac{\tau}{T}\right)+\varepsilon
T_{B}A_{1}\left(\frac{\tau}{T}\right)+O(\xi^{2}),$ (22)

where $A_{0}(\tau/T)=C^{L}_{**}(\tau/T,0)$.

Furthermore, when the intermediate asymptotics (17) is insufficient, we assume
that complete similarity for $C^{L}(\zeta,0)$ still holds but the constants
$A_{i}$ depend on $T_{\eta}$ and $T_{L}$ [40]. It should be noted that at this
situation, $\xi$ has a lower bound and therefore does not approach zero. Under
this assumption, in the inertial subrange, Eq.(22) may be modified as

$C^{L}(r_{0},T,\tau)=\varepsilon
TA_{0}\left(\frac{\tau}{T},\frac{T_{\eta}}{T_{L}}\right)+\varepsilon
T_{B}A_{1}\left(\frac{\tau}{T},\frac{T_{\eta}}{T_{L}}\right)+O(\xi^{2}).$ (23)

It should be noted that the width at $T$ of the ridge along the diagonal line
$t_{1}=t_{2}$ shown in Fig.1 is given by $T$ in the inertial subrange for any
Reynolds number.

On the other hand, in the case (ii), a scaling law for $C^{L}(\zeta,\xi)$ has
non-trivial scaling exponents which cannot be determined by dimensional
analysis: incomplete similarity. When $\xi$ is sufficiently small, as a
natural self-similar form suggested in [17], we propose

$C^{L}_{**}(\zeta,\xi)=G\xi^{\gamma}g^{L}\left(\frac{\zeta}{\xi^{\beta}}\right),$
(24)

where $g^{L}(X)$ is a dimensionless function and $g^{L}(0)=1$. Here, $G$ is a
non-zero constant factor and independent of $T_{B}$. It should be noticed that
Eq.(24) is consistent with the case (ii) since the function $g^{L}(X)$ is
bounded. The scaling exponents, $\beta$ and $\gamma$, are determined either by
the Navier-Stokes equations (or, more precisely, closure equations for the
TTLVCF) or by comparison with experimental data. As well as the case (i), when
the intermediate asymptotics (17) is insufficient and then $\xi$ is finite, we
expect that the incomplete self-similar form (24) still holds but the scaling
exponents $\beta$ and $\gamma$ may depend on $T_{\eta}$, $T_{L}$ and $T_{B}$.
The functional dependency of $\beta$ and $\gamma$ on $T_{B}$ are not obvious
like Eq.(23) because the Taylor expansion for $\xi$ does not exist.

As we will show in Sec.III, the case (ii) yields a better agreement with DNS
data of the 2D inverse energy-cascade turbulence than for the case (i).
Therefore, we conclude that the case (ii) holds for the 2D inverse energy-
cascade turbulence. Finally, we arrive at Eqs (9) and (10), which is supposed
to hold under the conditions (17) and (19).

Having obtained the ansatz for the temporal inertial subrange, we now consider
finite-Reynolds number effects. The argument below is heuristic and should be
justified experimentally. At finite Reynolds numbers, we assume that the above
self-similar form (9) is applicable. However, we assume that the scaling
exponents, $\beta$ and $\gamma$, are dependent on $T_{\eta}$ and $T_{L}$ as in
the Eqs (11) and (12). In the next section we show that these hypothetical
formulae of the exponents are useful to fit the DNS data obtained at finite
Reynolds numbers and to infer the asymptotic behavior of $C^{L}$.

It should be noted that undetermined scaling exponents such as $\beta$ and
$\gamma$ do not appear in scaling relations for the two-time Eulerian
correlation function $C^{E}(r,t_{1},t_{2})$ of the velocity increments, which
is defined by

$C^{E}(r,t_{1},t_{2})=\langle\delta u_{i}(\bm{x},\bm{r},t_{1})\delta
u_{i}(\bm{x},\bm{r},t_{2})\rangle,$ (25)

where $r=|\bm{r}|$ and the Eulerian velocity increment is given by $\delta
u_{i}(\bm{x},\bm{r},t_{1})=u_{i}(\bm{x}+\bm{r},t_{1})-u_{i}(\bm{x},t_{1})$. We
used spatial homogeneity and isotropy.

Because of the statistically steady state, $C^{E}(r,t_{1},t_{2})$ does not
depend on $T$ and can be written with all the governing parameters by

$C^{E}(r,t_{1},t_{2})=C^{E}(r,\tau,\varepsilon,\eta,L).$ (26)

According to the $\Pi$-theorem, there exists a dimensionless function,
$C^{E}_{*}$, by regarding $\varepsilon$ and $r$ as the independent parameters,
such that $C^{E}$ has the form,

$C^{E}(r,\tau,\varepsilon,\eta,L)=\varepsilon^{2/3}r^{2/3}C^{E}_{*}\left(\frac{\tau}{\varepsilon^{-1/3}r^{2/3}},\frac{\eta}{r},\frac{L}{r}\right).$
(27)

Furthermore, when we consider that $r$ is in the inertial subrange,

$\eta\ll r\ll L,$ (28)

we assume that, as $\eta\to 0$ and $L\to\infty$, the dependence on $\eta$ and
$L$ can be ignored. Then $C^{E}$ has a reduced form,

$C^{E}(r,\tau,\varepsilon)=C_{2}\varepsilon^{2/3}r^{2/3}g^{E}\left(\frac{\tau}{\varepsilon^{-1/3}r^{2/3}}\right),$
(29)

where $C_{2}$ is a universal constant related to the Kolmogorov constant and
the function $g^{E}(X)$ satisfies $g^{E}(0)=1$. This is consistent with
Kolmogorov’s phenomenology. In this way, the scaling law of the Eulerian
velocity correlation function can be determined by the dimensional analysis
thanks to the statistical stationarity. This is different from the TTLVCF.
However, it should be noted that the sweeping effect by large-scale advection
of eddies [41, 42] may be more dominant than the Kolmogorov time scale
$\varepsilon^{-1/3}r^{2/3}$. If that is the case, the Eulerian correlation
function may be different from the the scaling law given in Eq.(29), see
Ref.[43, 1] for review.

Before we leave this section, we summarize the assumptions that we made to
arrive at the scaling law (10), (11), and (12) with their physical meaning and
how to validate them.

(Assumption 1) Governing parameters of the TTLVCF

They appear in the arguments of the TTLVCF on the lhs of Eq.(14). It should be
noticed that we do not include parameters related to the external force. We
assumed here that they are not relevant for the sake of argument. However,
this may not be valid for the small-scale forcing in 2D turbulence. We will
later discuss this point in Sec. V.

(Assumption 2) Scale separations of time, Eqs.(17) and (19)

We assume that the Batchelor time $T_{B}$ is in the inertial range to arrive
at main result (10).

(Assumption 3) Vanishing dependence on $\bm{T_{\eta}}$ and $\bm{T_{L}}$ in the
intermediate asymptotics

This is related to Assumption 2 above. It is a natural assumption in the
inertial range as made in the K41 theory. However, it is not obvious whether
the assumption can be applied to the Lagrangian statistics. To validate this
point, we will perform DNS with different Reynolds numbers.

(Assumption 4) Persistent dependence of $\bm{T_{B}}$ in the intermediate
asymptotics

This is also related to Assumption 2 above. In our argument from Eq.(14) to
Eq.(19), we assume that the effects of the initial separations remain for a
long time and therefore that the statistics of particle pairs keep being
dependent on their histories. In the following sections, we will investigate
this memory effects of the initial separations by means of DNS of 2D inverse
energy-cascade turbulence.

(Assumption 5) Functional form of the TTLVCF at short times, Eq.(24)

Equation (24) is a standard candidate of the functional forms in the
incomplete self-similarity framework. The validation of this assumption with
DNS is our principal goal of the next section III.

(Assumption 6) Finite Reynolds-number effects on the exponents $\bm{\beta}$
and $\bm{\gamma}$

In Eqs.(11) and (12), we proposed that specific functional forms of $\beta$
and $\gamma$ where $T_{\eta}$ and $T_{L}$ re-appear, in spite of Assumptions 2
and 3. This is an empirical form we find with DNS to parameterize finite
Reynolds number effects, as we will see in the next section.

##  III Numerical experiments

###  III.1 Numerical details

We perform DNS of the 2D inverse energy-cascade turbulence in order to
numerically verify the ansatz made in Sec.II. We suppose that the Eulerian
velocity field, $\bm{u}(\bm{x},t)$, follows the 2D Navier-Stokes equations
with forcing, hyperviscous, and hypodrag terms. We numerically solve the
equations in terms of the vorticity,

$\frac{\partial\omega}{\partial
t}+\left(\bm{u}\cdot\nabla\right)\omega=(-1)^{h+1}\nu\Delta^{h}\omega+\alpha\Delta^{-q}\omega+f,$
(30)

where $\omega$ is vorticity filed,
$\omega(\bm{x},t)=\partial_{x}u_{y}(\bm{x},t)-\partial_{y}u_{x}(\bm{x},t)$.
The hyperviscous, and hypodrag terms are the first and second ones on the
right hand side (rhs) of Eq.(30) and $f$ is an external forcing term. For the
2D inverse energy-cascade turbulence, the smallest and largest time scales can
be explicitly described by the viscous coefficient, $\nu$, and the drag
coefficient, $\alpha$, respectively as below:

$T_{\eta}\equiv\left(\frac{\nu}{\varepsilon^{h}}\right)^{\frac{1}{3h-1}},\quad
T_{L}\equiv\left(\frac{1}{\alpha\varepsilon^{q}}\right)^{\frac{1}{3q+1}}.$

This is an advantage of the 2D energy inverse-cascade turbulence, because the
integral time scale of the 3D turbulence cannot be explicitly described.

The forcing term, $f({\bm{x}},t)$, is given in terms of the Fourier
coefficients,
$\hat{f}(\bm{k},t)=k^{2}\varepsilon_{in}/[n_{f}\hat{\omega}^{*}(\bm{k},t)]$,
where $\hat{f}$ denotes the Fourier transform of the function $f(\bm{x},t)$.
The energy input rate is denoted by $\varepsilon_{in}$, and $n_{f}$ denotes
the number of the Fourier modes in the following forcing wavenumber range. We
select the coefficients, $\hat{f}({\bm{k}},t)$, as non-zero only in high wave
numbers, $\bm{k}$, satisfying $k_{f}-1<|\bm{k}|<k_{f}+1$. Thus, the energy
input rate is maintained as a constant in time. Numerical integration of Eq.
(30) is performed via the pseudospectral method with the 2/3 dealiasing rule
in space and the 4-th order Runge–Kutta method in time. The setting and our
numerical method are identical to those used in [44, 45]. The typical
wavenumber of the hypodrag is dimensionally estimated as
$(\alpha^{3}/\varepsilon)^{1/(6q+2)}$, which is termed as the frictional wave
number, $k_{\alpha}$. Here we use the infrared Reynolds number,
$\mathrm{Re}_{\alpha}\equiv k_{f}/k_{\alpha}$, as proposed by Vallgren [46] in
order to quantify the span of the inertial subrange. In Table 1 we list the
parameters of simulations used in the study.

$N^{2}$ | $\delta x$ | $\delta t$ | $\nu$ | $h$ | $\alpha$ | $q$ | $k_{f}$ | $\varepsilon_{in}$ | $N_{p}^{2}$ | $\varepsilon$ | $\sigma_{\varepsilon}$ | $L$ | $u_{\mathrm{rms}}$ | $\mathrm{Re}_{\alpha}$ | $T_{\eta}$ | $T_{L}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$1024^{2}$ | $0.006$ | $0.002$ | $1.82\times 10^{-38}$ | 8 | $35$ | $1$ | $249$ | $0.1$ | $2048^{2}$ | $0.019$ | $2.9\times 10^{-4}$ | $0.38$ | $0.5$ | $40$ | $0.091$ | $1.1$
$2048^{2}$ | $0.003$ | $0.001$ | $4.664\times 10^{-43}$ | 8 | $35$ | $1$ | $496$ | $0.1$ | $2048^{2}$ | $0.019$ | $2.9\times 10^{-4}$ | $0.37$ | $0.5$ | $80$ | $0.057$ | $1.1$
$4096^{2}$ | $0.0015$ | $0.001$ | $1.05\times 10^{-47}$ | 8 | $35$ | $1$ | $997$ | $0.1$ | $4096^{2}$ | $0.019$ | $2.6\times 10^{-4}$ | $0.36$ | $0.5$ | $160$ | $0.036$ | $1.1$
Table 1: Parameters of numerical simulations: $N^{2}$, $\delta x=2\pi/N$,
$\delta t$, $\nu$, $h$, $\alpha$, $q$, $k_{f}$, $\varepsilon_{in}$, and
$N_{p}^{2}$ denote the number of grid points, grid spacing, size of the time
step, hyperviscosity coefficient, order of the Laplacian of the
hyperviscosity, hypodrag coefficient, order of the inverse Laplacian of the
hypodrag, forcing wavenumber, energy input rate of the forcing, and, number of
the Lagrangian particles, respectively. Turbulent characteristics:
$\varepsilon$, $\sigma_{\varepsilon}$, $L$, $u_{\mathrm{rms}}$, $Re_{\alpha}$,
$T_{\eta}$, $T_{L}$, and $N_{p}^{2}$ denote mean of the resultant energy flux
in the inertial subrange, standard deviation of the resultant energy flux,
integral scale, root-mean-square velocity, infrared Reynolds number, viscous
time scale, and, drag time scale, respectively.

To obtain the Lagrangian statistics, we employ a standard particle tracking
method. The flow is seeded with a large number of tracer particles. The number
of particles, $N_{p}^{2}$, for each simulation is described in Table 1. The
particles are tracked in time via integrating the advection equation,

$\frac{\mathrm{d}}{\mathrm{d}t}\bm{x}_{p}(t)=\bm{u}(\bm{x}_{p}(t),t),$ (31)

where $\bm{x}_{p}(t)$ denotes the particle position vector. The numerical
integration of Eq. (31) is performed using the Euler method. The velocity
value at an off-grid particle position is estimated by the fourth-order
Lagrangian interpolation of the velocity calculated on the grid points.

In Eq.(30), we use the hyperviscosity, $h=8$, rather than the normal
viscosity, $h=1$, for DNSs. This is because the hyperviscosity extends the
inertial subrange for a given spatial resolution. We confirmed that the
hyperviscosity does not affect the particle-pair statistics in Ref. [33].

First of all, let us consider to what extent the assumptions on the time
separations, $T_{\eta}\ll T\ll T_{L}$ (Eq.(17)) and $T_{B}\ll T$ (Eq.(19)),
made in Sec.II hold in our DNS. In the DNS, $T_{L}/T_{\eta}\lesssim 10^{2}$.
Certainly, this poses limitations on studying whether the asymptotic behavior
of the TTLVCF, Eq.(9), is valid. In theory, if $T_{\eta}\ll T_{B}\ll T_{L}$,
then the particle pairs may be hardly influenced by neither the viscosity nor
the large scale drag from the beginning of the relative diffusion. However, in
practice, due to the limited scale separation, $\xi=T_{B}/T$ may not become
sufficiently small in our DNS, as $T=(t_{1}+t_{2})/2$ increases while
satisfying $T_{\eta}<T_{B}<T<T_{L}$. Therefore, it is inevitable to consider
that the numerically obtained TTLVCF, $C^{L}(r_{0},T,\tau,\varepsilon)$,
depends on $T_{\eta}$ and $T_{L}$ even if the large $T_{B}$ condition,
$T_{\eta}<T_{B}<T_{L}$, is satisfied.

Given these practical limitations, it is useful to relax the large $T_{B}$
condition and to consider the case $T_{B}<T_{\eta}$, which we call the small
$T_{B}$ condition. Obviously under the small $T_{B}$ condition, we cannot
ignore viscous effects on particle-pair statistics. However, the value of
$\xi=T_{B}/T$ can become smaller as the average time $T$ increases in
$T_{\eta}\ll T\ll T_{L}$ than under the large $T_{B}$ condition. Some previous
studies investigate a particle-pair statistics under the small $T_{B}$
condition [27, 32, 47, 48]. Of course, it is not obvious that the two
different conditions give the same asymptotic behavior of
$C^{L}(r_{0},T,\tau,\varepsilon)$ as $\xi=T_{B}/T\to 0$. Therefore, we
investigate dependencies on $T_{\eta},T_{L}$ and $T_{B}$ for both conditions
in the following subsections.

More specifically, we investigate the two scaling exponents, $\beta$ and
$\gamma$, appeared in our ansatz (9). For this purpose, we decompose the
TTLVCF $C^{L}(r_{0},T,\tau,\varepsilon)$ into two parts:

$C^{L}(r_{0},T,\tau,\varepsilon)=C_{d}^{L}(T,T_{B})C_{p}^{L}(T,\tau,T_{B}),$
(32)

where $C_{d}^{L}(T,T_{B})$ corresponds to the TTLVCF along the diagonal line
$t_{1}=t_{2}$, that is, $C_{d}^{L}(T,T_{B})\equiv
C^{L}(r_{0},T,\tau=0,\varepsilon)$. The other part $C_{p}^{L}(T,\tau,T_{B})$
corresponds to the TTLVCF along a line $t_{1}+t_{2}=2T$, which is
perpendicular to the diagonal line. Its value at $\tau=0$ is normalized:
$C_{p}^{L}(T,\tau=0,T_{B})=1$. If the ansatz is correct,
$C_{p}^{L}(T,\tau,T_{B})=g^{L}(\tau/[T^{1-\beta}T_{B}^{\beta}])$.

In what follows, the values of the exponents, $\gamma$ and $\beta$, are
estimated from numerically calculated $C_{d}^{L}(T,T_{B})$ and
$C_{p}^{L}(T,\tau,T_{B})$, respectively, as we vary $T_{B}$ and $T_{\eta}$. We
consider first the large $T_{B}$ condition ($T_{B}>T_{\eta}$) and then the
small $T_{B}$ condition.

###  III.2 Large $\bm{T_{B}}$ condition: $\bm{T_{\eta}\ll T_{B}\ll T\ll
T_{L}}$

Figure 2:  Color maps of time evolution of the TTLVCF
$C^{L}_{ii}(r_{0},t_{1},t_{2})$ defined in Eq.(5) with
$T_{B}=(r_{0}^{2}/\varepsilon)^{1/3}=3.5T_{\eta}$. (a)
$\mathrm{Re}_{\alpha}=40$, (b) $\mathrm{Re}_{\alpha}=80$, (c)
$\mathrm{Re}_{\alpha}=160$. The white lines indicate $t_{i}=T_{B}$ ($i=1,2$).
Here the time axes, $t_{i}/T_{\eta}$ ($i=1,2$), span from $0$ to
$T_{L}/T_{\eta}$.

Let us consider the scaling laws of $C^{L}(r_{0},T,\tau,\varepsilon)$ under
the large $T_{B}$ condition. In DNS, although this condition,
$T_{\eta}<T_{B}<T_{L}$ is satisfied, Assumption 2 is not sufficient even in
our largest simulation with $\mathrm{Re}_{\alpha}=160$. Figure 2 shows color
maps of $C^{L}(r_{0},T,\tau,\varepsilon)$ in terms of the original time
variables $t_{1}$ and $t_{2}$ with $T_{B}=3.5T_{\eta}$ for the three values of
$\mathrm{Re}_{\alpha}$. We observe that the width at $T$ of the ridge along
the diagonal line (the region where $C^{L}(r_{0},T,\tau,\varepsilon)$ remains
large) becomes wider as the average time $T=(t_{1}+t_{2})/2$ increases. We
also observe that qualitatively this tendency appears to be independent of
$\mathrm{Re}_{\alpha}$.

Figure 3:  (a) Time evolution of $C_{d}^{L}(T,T_{B})$ normalized by
$(r_{0}/T_{B})^{2}$ for various $T_{B}$’s at $\mathrm{Re}_{\alpha}=40$ (dashed
dotted), $80$ (dashed), and $160$ (solid). (b) Logarithmic local slope (LLS)
of $C_{d}^{L}(T,T_{B})$ at $\mathrm{Re}_{\alpha}=40$ (dashed) and $160$
(solid). The filled circle on each line indicates the position of the maximum.
(c) Value of the exponent $\gamma$ suggested by the maximum value of the LLS
plotted as a function of $(T_{B}-T_{\eta})/T_{L}$ at $\mathrm{Re}_{\alpha}=40$
(red) and $160$ (green). The blue solid line corresponds to
$1.34[(T_{B}-T_{\eta})/T_{L}]^{0.5}$, which is determined by the least square
fit in the range $0.039\leq(T_{B}-T_{\eta})/T_{L}\leq 0.16$.

Now we focus on the behavior of the correlation function along the diagonal
line. Figure 3(a) shows time evolution of
$C_{d}^{L}(T,T_{B})=C^{L}(r_{0},T,\tau=0,\varepsilon)$ for various $T_{B}$’s.
Obviously, it indicates that a scaling exponent, if it exists, depends on
$T_{B}$. Figure 3(b) shows logarithmic local slopes (LLSs) of
$C_{d}^{L}(T,T_{B})$. If the ansatz (9) is valid, the LLS becomes $1-\gamma$
(constant). We see that a narrow plateau region for each LLS. As $T_{B}$
approaches $T_{\eta}$ from above, we observe that it becomes wider and that
the value of the plateau region becomes closer to $1$ which corresponds to
$\gamma=0$. For further quantification, we infer the value of $\gamma$ for
each curve from the maximum value of the LLS. The data shown here indicate the
dependence on $T_{B},T_{\eta}$ and $T_{L}$. To circumvent this, we now use the
empirical form for $\gamma$ given in Eq.(12).

Figure 3(c) shows the maximum values of the LLSs, which we regard as $\gamma$
in Eq.(12); $\gamma$ is dependent on $T_{B},T_{\eta}$ and $T_{L}$. The
horizontal axis of Fig.3(c) is set to $(T_{B}-T_{\eta})/T_{L}$. We find
empirically this combination of the independent variables,
$(T_{B}-T_{\eta})/T_{L}$, to make the data points collapse onto a single
curve. The first observation concerns the behavior as
$(T_{B}-T_{\eta})/T_{L}\to 0$ (when $T_{B}$ approaches $T_{\eta}$ from above):
the exponent $\gamma$ seems to approach $0$. However, this limit $T_{B}\to
T_{\eta}$ violates the large $T_{B}$ condition, $T_{\eta}\ll T_{B}\ll
T_{\eta}$. The second observation is about the behavior under the large (but
smaller than 1) $(T_{B}-T_{\eta})/T_{L}$ range, which is consistent with the
large $T_{B}$ condition. In this range, we observe that the master curve
becomes independent of $\mathrm{Re}_{\alpha}$ as shown in Fig.3(c). Our best
fit function to the curve for the exponent $\gamma$ is

$\gamma\left(\frac{T_{\eta}}{T_{B}},\frac{T_{L}}{T_{B}}\right)=\overline{\gamma}_{0}\left(\frac{T_{B}-T_{\eta}}{T_{L}}\right)^{1/2},$
(33)

which is plotted as a solid line in Fig.3(c). Here $\overline{\gamma}_{0}$ is
a constant estimated about $1.34\pm 0.01$ by the fitting.

Now we come back to the ansatz (9) leading to
$C^{L}(r_{0},T,\tau=0,\varepsilon)\propto T^{1-\gamma}$. The functional form
of the exponent (33) indicates that $\gamma\to 0$ at infinite Reynolds number.
This implies that the K41 scaling,
$C^{L}(r_{0},T,\tau=0,\varepsilon)=C_{d}^{L}(T,T_{B})\propto T$, is recovered
under the sufficient scale separation. However, for this recovery, the
exponent $1/2$ in Eq.(33) suggests that we need an enormously large
$\mathrm{Re}_{\alpha}$. For example, in order to get the value of $\gamma$
valid for one effective figure, $\gamma\sim 0.01$, we may need
$T_{L}/T_{\eta}\sim 10^{5}$ (in our DNS here $T_{L}/T_{\eta}\sim 30$ at most),
which may correspond to $\mathrm{Re}_{\alpha}\sim 10^{6}$.

On the other hand, at small values of $(T_{B}-T_{\eta})/T_{L}$, $\gamma$
deviates from the relation (33) as shown in Fig. 3(c). Let us suppose that the
deviation persists at larger Reynolds numbers. Then $\gamma$ may have a
negative limit value as $T_{B}\to T_{\eta}$ (approaching $T_{\eta}$ from
above). If we extrapolate the deviation to $(T_{B}-T_{\eta})/T_{L}=0$ with
linear decrease, the limit value of $\gamma$ is about $-0.25$. We cannot
conclude whether the deviation remains at sufficient large Reynolds numbers
from our DNS.

In summary of the result for the TTLVCF along the diagonal line, our
simulation data suggest that the scaling law of $C_{d}^{L}(T,T_{B})$ at
sufficiently large $\mathrm{Re}_{\alpha}$ is,

$C^{L}(r_{0},T,\tau=0,\varepsilon)=C_{d}^{L}(T,T_{B})=G\varepsilon
T_{B}^{\overline{\gamma}_{0}\sqrt{\frac{T_{B}-T_{\eta}}{T_{L}}}}~{}T^{1-\overline{\gamma}_{0}\sqrt{\frac{T_{B}-T_{\eta}}{T_{L}}}}.$
(34)

where $G$, which is the constant appeared in Eq.(9), is estimated as $G\sim
80$ from the compensated plot of Fig.3(a) by $T_{B}^{\gamma}T^{1-\gamma}$ (the
compensated plot is not shown).

It is noted that $T_{L}$ is kept constant and $T_{\eta}$ is changed when
$\mathrm{Re}_{\alpha}$ is increased in our DNS. Thus, the limit
$(T_{B}-T_{\eta})/T_{L}\to 0$ is consistent with $T_{B}\to T_{\eta}$ in this
study. On the other hand, if we can change both values of $T_{\eta}$ and
$T_{L}$, the limit $(T_{B}-T_{\eta})/T_{L}\to 0$ indicates two states: One is
$T_{B}\to T_{\eta}$ and the other is $T_{L}\to\infty$ while keeping
$T_{B}-T_{\eta}$ constant, where the large $T_{B}$ condition can hold. Hence,
if we change the value of $T_{L}$ and fix the value of $T_{B}-T_{\eta}$, the
data points of $\gamma$ may collapse onto another curve different from the
former one at small $(T_{B}-T_{\eta})/T_{L}$, which indicates the latter
state. Nevertheless, we assume that the asymptotic value for the limit
$(T_{B}-T_{\eta})/T_{L}\to 0$ is the same for the two states.

Next, we consider the behavior of the TTLVCF along the lines perpendicular to
the diagonal line, which is given by $C_{p}^{L}(T,\tau,T_{B})$. Figure 4 shows
sectional views of the color map shown in Fig.2 for various sections given by
the lines $t_{1}+t_{2}=2T$. It should be noted that the curves shown in Fig.4
are normalized by $C_{d}^{L}(T,T_{B})$. Hence they are the graphs of
$C_{p}^{L}(T,\tau,T_{B})$ as a function of the relative time
$\tau=t_{2}-t_{1}$ at a fixed average time $T$.

We first notice that the typical width of the peak of
$C_{p}^{L}(T,\tau,T_{B})$ centered at zero relative time $\tau=0$ is given by
the dissipation time scale $T_{\eta}$ initially, i.e., for small average time
$T$. Then, the width becomes larger and larger as the average time increases.
At large average times, $T\lesssim T_{L}$, the function
$C_{p}^{L}(T,\tau,T_{B})$ decreases exponentially as shown in the insets of
Fig. 4. Moreover, the data indicate that $C_{p}^{L}(T,\tau,T_{B})$ decreases
faster than exponential at $\tau\sim T_{L}$.

Figure 4:  Normalized correlation function, $C_{p}^{L}(T,\tau,T_{B})$, defined
in Eq.(32) as a function of the relative time $\tau$ with $T_{B}=3.5T_{\eta}$
for various $T$’s. The average time varies in $T_{\eta}<T<0.54T_{L}$ and the
corresponding curves are colored from black to yellow. The three panels
correspond to (a) $\mathrm{Re}_{\alpha}=40$, (b) $\mathrm{Re}_{\alpha}=80$,
and (c) $\mathrm{Re}_{\alpha}=160$. Two vertical dashed lines in each panel
show $\tau/T_{B}=\pm T_{\eta}/T_{B}$, respectively. The horizontal axis spans
in $-T_{L}/T_{B}\leq\tau/T_{B}\leq T_{L}/T_{B}$. The insets show the same
plots as the outsets but in the lin-log coordinates.

To quantify the decay of $C_{p}^{L}(T,\tau,T_{B})$, we use an $n$-th decay
time scale, $\tau_{1/n}(T)$, defined as

$C_{p}^{L}(T,~{}\tau=\tau_{1/n}(T),~{}T_{B})=\frac{1}{n}.$ (35)

If $\tau_{1/n}(T)$ is power-law such as $\tau_{1/n}(T)\propto T^{1-\beta}$ and
the scaling exponent, $\beta$, is independent of the value of $n$, then
$C_{p}^{L}(T,\tau,T_{B})$ has the self-similar form of the ansatz (9), that
is, $C_{p}^{L}(T,\tau,T_{B})=g^{L}(\tau/[T_{B}^{\beta}T^{1-\beta}])$.

If $\beta$ depends on $n$, we can still expect that $C_{p}^{L}(T,\tau,T_{B})$
has a self-similar form in a certain interval of $\tau$. Hereafter, this
$n$-dependent $\beta$ is denoted by $\beta_{1/n}$. The interval is determined
by the value of $\tau_{1/n}(T)$.

Figure 5 shows $\tau_{1/n}(T)$ with $n=2,8$, and $32$ for various $T_{B}$’s.
With a small $n$ such as $n=2$, we probe the behavior in the vicinity of the
peak of $C_{p}^{L}(T,\tau,T_{B})$, and, with a large $n$ such as $n=32$, we
characterize the behavior in the tail region of $C_{p}^{L}(T,\tau,T_{B})$. For
$n=2$ as shown in Fig. 5(a), $\tau_{1/2}(T)$ strongly depends on $T_{B}$. This
is because $\tau_{1/2}(T)$ is smaller than $T_{B}$ for almost all $T$s.
Nevertheless, there may be a power-law behavior in a certain range of $T$. On
the other hand, for larger $n$ such as $n=8$ and $32$, as shown in Fig.5(b)
and (c), the power law behavior of $\tau_{1/n}(T)$ becomes clearer and
$\tau_{1/n}(T)\propto T^{1-\beta_{1/n}}$ holds at a certain time interval of
$T$. The scaling exponents, $\beta_{1/n}$ appear to become independent of
$T_{B}$ and the scaling region becomes larger as increasing
$\mathrm{Re}_{\alpha}$. These observations lead us to conclude that the ansatz
(9) is a reasonable description of the function $C_{p}^{L}(T,\tau,T_{B})$.

However, as shown in the insets of Fig. 5, the LLSs are too noisy to determine
the value of $\beta_{1/n}$ accurately. The noise may be suppressed as we
increase massively the number of particle-pair samples. Instead, here we use
compensated plots of Fig. 5 to estimate the value of the scaling exponent
$\beta_{1/n}$. The compensation is based on the self-similar variable
$\zeta/\xi^{\beta}$ in Eq.(24), which is the argument of the function $g^{L}$.
If the self-similarity is valid at $\tau=\tau_{1/n}(T)$ with the exponent
$\beta_{1/n}$, the self-similar variable

$\left.\frac{\zeta}{\xi^{\beta_{1/n}}}\right|_{\tau=\tau_{1/n}(T)}=\frac{\tau_{1/n}(T)}{T_{B}^{\beta_{1/n}}T^{1-\beta_{1/n}}}\equiv
D_{n,\beta_{1/n}}(T)$ (36)

becomes constant which neither depends on $T$ nor $T_{B}$. For each $n$, we
plot $D_{n,q}(T)$ by varying $q$ and find $q_{*}$ that gives the widest flat
region as a function of $T$. We regard this $q_{*}$ as $\beta_{1/n}$. We show
$D_{n,\beta_{1/n}}(T)$ in Fig. 6 for $n=2,8$, and $32$. These compensated
plots are less noisy than the LLSs, but they still have tiny oscillations. As
increasing $\mathrm{Re}_{\alpha}$, we observe that $D_{n,\beta_{1/n}}(T)$
becomes independent of $T_{B}$ except for $T_{B}=1.4T_{\eta}$, in particular,
for $n=32$ as shown in Fig. 6 (c). This indicates that the ansatz (9) is
reasonable for $C_{p}^{L}(T,~{}\tau,T_{B})$, albeit that the numerical data is
noisy.

Figure 5:  $n$-th decay time scale, $\tau_{1/n}(T)$ as a function of the
average time $T$ for (a) $n=2$, (b) $n=8$, (c) $n=32$ at
$\mathrm{Re}_{\alpha}=40$ (dashed dotted), $\mathrm{Re}_{\alpha}=80$ (dashed),
and $\mathrm{Re}_{\alpha}=160$ (solid). The insets show the LLS of
$\tau_{1/n}(T)$ shown in the ousets.

Figure 6:  Compensated plots of $\tau_{1/n}(T)$ by $T^{1-\beta_{1/n}}$ for (a)
$n=2$, (b) $n=8$, (c) $n=32$ at $\mathrm{Re}_{\alpha}=80$ (dashed), and
$\mathrm{Re}_{\alpha}=160$ (solid). The values of $\beta_{1/n}$ are determined
in such a way that each of the compensated graphs has the widest flat region.

Now let us specify the empirical form of the exponent $\beta$ given in
Eq.(11). Figure 7 shows the measured values of $\beta_{1/n}$ with $n=2,8$ and
$32$ as a function of $(T_{B}-T_{\eta})/T_{L}$ that is the same independent
variable used in Fig.3 (c).

For $n=2$, the outset of Fig. 7(a) indicates that $(T_{B}-T_{\eta})/T_{L}$ is
not appropriate since the data points are still scattered. This leads us to
search for a more suitable self-similar variable for $\beta_{1/2}$, which is
found to be $T_{B}/T_{L}$ as shown in the inset of Fig. 7(a). This implies
that $\beta_{1/2}$ depends only on $T_{B}$ and $T_{L}$, but not on $T_{\eta}$.
Empirically we now fit the collapsed curve obtained in the inset of Fig. 7(a)
with a function $\check{\beta}(T_{B}/T_{L}+a)^{b}$ with constants
$\check{\beta},a$, and $b$. Our result is
$\beta_{1/2}\propto[T_{B}/T_{L}+\omega_{1}]^{0.4}$, where the constant
$\omega_{1}$ takes zero or a non-zero small value possibly in a range,
$0\leq\omega_{1}\lesssim 0.01$. The fitted functions are shown in the inset of
Fig. 7(a).

On the other hand, for $n=8$ and $32$, the exponents $\beta_{1/8}$ and
$\beta_{1/{32}}$ are dependent on $T_{B}$, $T_{\eta}$, and $T_{L}$. Our best
fit is $\beta_{1/n}\propto[(T_{B}-T_{\eta})/T_{L}+\omega_{2}]^{0.4}$ for $n=8$
and $32$, as shown in Figs. 7(b) and (c). Here, $\omega_{2}$ is a constant in
a range $0\leq\omega_{2}\lesssim 0.01$. The accurate values of $\omega_{1}$
and $\omega_{2}$ cannot be determined from the data shown in Fig. 7. This is
because the data are noisy and also the Reynolds numbers are not sufficiently
large for studying the behavior in $T_{B}/T_{L}\to 0$. Nevertheless, it is
obvious that the behavior of $\beta_{1/2}$ is different from the others. On
the other hand, for larger $n$’s such as $n=8$ or $32$, the behaviors of
$\beta_{1/n}$ are similar to each other. Therefore, these results suggest that
the exponent $\beta$ in the ansatz has two different self-similar forms
depending on $\tau\lesssim T_{B}$ and $\tau\gg T_{B}$ at sufficiently large
Reynolds numbers. Specifically, we infer from the data

$\beta\left(\frac{T_{\eta}}{T_{B}},\frac{T_{L}}{T_{B}}\right)=\begin{cases}\displaystyle\left(\frac{T_{B}}{T_{L}}+\omega_{1}\right)^{0.4}\equiv\beta_{1}&\text{for}\quad\tau\lesssim
T_{B},\\\ &\\\
\displaystyle\left(\frac{T_{B}-T_{\eta}}{T_{L}}+\omega_{2}\right)^{0.4}\equiv\beta_{2}&\text{for}\quad\tau\gg
T_{B},\end{cases}$ (37)

where $\omega_{1}$ and $\omega_{2}$ are $\mathrm{Re}_{\alpha}$ independent
constants, which may be zero. Accordingly, the function $g^{L}$ in the ansatz
(9) can be given by

$g^{L}\left(\frac{\tau}{T_{B}^{\beta}T^{1-\beta}}\right)=\begin{cases}\displaystyle
g^{L}_{1}\left(\frac{\tau}{T_{B}^{\beta_{1}}T^{1-\beta_{1}}}\right)&\text{for}\quad\tau\lesssim
T_{B},\\\ &\\\ \displaystyle
g^{L}_{2}\left(\frac{\tau}{T_{B}^{\beta_{2}}T^{1-\beta_{2}}}\right)&\text{for}\quad\tau\gg
T_{B},\end{cases}$ (38)

where $g^{L}_{1}$ and $g^{L}_{2}$ are self-similar functions.

Now we discuss the limit of $\beta$ as $T_{\eta}\to 0$ and $T_{L}\to\infty$.
Here, $\omega_{1}^{0.4}$ in Eq.(37) is the limit of $\beta_{1}$ as
$T_{B}/T_{L}\to 0$. Similarly, $\omega_{2}^{0.4}$ is the limit of $\beta_{2}$
as $(T_{B}-T_{\eta})/T_{L}\to 0$. Let us suppose $\omega_{2}=0$. Then the K41
scaling law is recovered at $\tau\gg T_{B}$ at sufficiently large Reynolds
numbers. It is impossible to determine the accurate value of $\omega_{2}$ from
Fig. 7. It appears that $\omega_{2}=0.01$ is the best fitted value judging
from Fig. 7(c) though $\omega_{2}=0$ is not ruled out. Both values
$\omega_{1}=0.01$ and $\omega_{1}=0$ seem equally good as in the case for
$\omega_{2}$. In order to determine the accurate values of $\omega_{1}$ and
$\omega_{2}$, we need to perform DNSs at much larger Reynolds number and with
much larger number of the particle pairs.

Figure 7:  Scaling exponents, $\beta_{1/n}$, evaluated from the compensated
plots of Fig. 5 for (a) $n=2$, (b) $n=8$, and (c) $n=32$ at
$\mathrm{Re}_{\alpha}=80$ (red) and $\mathrm{Re}_{\alpha}=160$ (green). The
inset of (a) shows the same plots in the outset, but the horizontal axis is
changed to $T_{B}/T_{L}$ from $(T_{B}-T_{\eta})/T_{L}$. The orange solid line
shows $[T_{B}/T_{L}]^{0.4}$. The blue solid line shows
$[T_{B}/T_{L}+0.01]^{0.4}$. For (b) and (c), the orange solid line shows
$[(T_{B}-T_{\eta})/T_{L}]^{0.4}$. The blue solid line shows
$[(T_{B}-T_{\eta})/T_{L}+0.01]^{0.4}$.

In Eq.(38), the two different behaviors of $g^{L}$ are inferred from those of
$\beta$. We now demonstrate that the two forms are consistent with the DNS
data. Figure 8 shows $C_{p}^{L}(T,\tau,T_{B})$ as a function of
$\tau/[T^{1-\beta_{1/n}}T_{B}^{\beta_{1/n}}]$ for $n=2$ in Fig.8(a) and $n=32$
in Fig.8(b). Figure 8(a) for $n=2$ is plotted in lin-lin coordinates, which
means that we mainly observe the regions where $C_{p}^{L}(T,\tau,T_{B})$ is
large. On the other hand, Fig. 8(b) for $n=32$ is plotted in lin-log
coordinates, which means that we mainly observe the regions where
$C_{p}^{L}(T,\tau,T_{B})$ is small. The master curve in Fig.8(a) corresponds
to $g_{1}^{L}$ and the one in Fig.8(b) corresponds to $g_{2}^{L}$ in Eq.(38).
Here we assume $\beta_{1}=\beta_{1/2}$ and $\beta_{2}=\beta_{1/32}$. Compare
the collapsed curves in Fig.8 to those shown in Fig.4 without taking any
appropriate similarity variable. Furthermore, let us assume that the rescaled
functions are exponential, namely $g_{1}^{L}(X)=\exp(-k_{1}X)$ and
$g_{2}^{L}(X)=\exp(-k_{2}X)$. This assumption is consistent with Fig.8. We can
estimate the constants as $k_{1}\sim 2.3$ and $k_{2}\sim 3.0$ from Fig.8,
though these values are also slightly dependent on $T_{\eta},T_{L}$, and
$T_{B}$. The exponential forms will be used to estimate the Richardson
constant in Sec.IV.

Figure 8:  Normalized correlation function, $C_{p}^{L}(T,\tau,T_{B})$ for
$T_{B}=2.2T_{\eta}$ rescaled as a function of (a)
$\tau/[T_{B}^{\beta_{1/2}}~{}T^{1-\beta_{1/2}}]$ and (b)
$\tau/[T_{B}^{\beta_{1/32}}~{}T^{1-\beta_{1/{32}}}]$, where $\beta_{1/2}=0.64$
and $\beta_{1/{32}}=0.76$, for $T_{B}<T<0.54T_{L}$ at
$\mathrm{Re}_{\alpha}=160$. The colors of the curves change from black to
yellow as the average time $T$ increases, which is similar to Fig.4. Black
dashed lines show (a)
$\exp(-2.3|\tau|/[T_{B}^{\beta_{1/2}}~{}T^{1-\beta_{1/2}}])$ and (b)
$\exp(-3.0|\tau|/[T_{B}^{\beta_{1/32}}~{}T^{1-\beta_{1/32}}])$. (c) Same as
(a) but for $T_{B}=3.5T_{\eta}$, where $\beta_{1/2}=0.42$. (d) Same as (b) but
for $T_{B}=3.5T_{\eta}$, where $\beta_{1/32}=0.30$.

###  III.3 Small $T_{B}$ condition: $\bm{T_{B}<T_{\eta}}$

Now we consider the scaling laws of $C^{L}(r_{0},T,\tau,\varepsilon)$ under
the small initial separation condition, $T_{B}<T_{\eta}$, where particle pairs
may be strongly influenced by small-scale effects caused by the viscosity and
the forcing. On the other hand, the condition $T_{B}\ll T_{L}$ is met more
easily than in the previous large $T_{B}$ condition. Hence, we expect that
$C^{L}(r_{0},T,\tau,\varepsilon)$ is independent of large-scale effects such
as the drag. Moreover, the $t^{3}$ scaling law for $\langle r^{2}(t)\rangle$
has been observed under this condition in many previous studies for both 2D
and 3D as mentioned in Sec.I. We also investigate the reason why the $t^{3}$
scaling is observed even at moderate Reynolds numbers only for a tuned initial
separation $r_{0}$.

In this subsection, we repeat what we have done in the previous subsection.
Therefore, we describe only the diffences. In terms of the correlation along
the diagonal line, the exponent $\gamma$ is determend as we did in Fig. 3(b),
see Fig. 9(a) and (b). In what follows we write the exponents with check in
the small $T_{B}$ condition. In Fig. 9(c), we show the measured
$\check{\gamma}$ as a function of $T_{B}/T_{\eta}$. This choice of the
variable yields a curve independent of $\mathrm{Re}_{\alpha}$, which we fit
with

$\check{\gamma}(T_{B},T_{\eta})=\ln\left(\frac{T_{B}}{T_{\eta}}\right)-\upsilon.$
(39)

Here the constant $\upsilon$ is $0.25\pm 0.02$ which is determined by a least-
square method. Therefore, at $T_{B}<T_{\eta}$, the scaling law of
$C_{d}^{L}(T,T_{B})$, can be,

$C_{d}^{L}(T,T_{B})=G\varepsilon
T\left(\frac{T_{B}}{T}\right)^{\check{\gamma}}=GT^{1-\ln\left(\frac{T_{B}}{T_{\eta}}\right)+\upsilon}.$
(40)

Figure 9:  (a) Time evolution of $C_{d}^{L}(T,T_{B})$ for various $T_{B}$’s at
$\mathrm{Re}_{\alpha}=40$ (dashed dotted), $80$ (dashed), and $160$ (solid).
(b) LLSs of $C_{d}^{L}(T,T_{B})$ at $\mathrm{Re}_{\alpha}=40$ (dashed) and
$160$ (solid). Filled circles on the curves show positions of their maximum
value. (c) Values of the exponent $\check{\gamma}$ as a function of
$T_{B}/T_{\eta}$, which are measured by the maximum value of the LLSs at
$\mathrm{Re}_{\alpha}=40$ (red) and $160$ (green). The blue solid line shows
$\ln(T_{B}/T_{\eta})-0.25$.

Next, we focus on the other part of the correlation $C_{p}^{L}(T,\tau,T_{B})$
defined in Eq.(32) and the scaling exponent, $\check{\beta}$. Figure 10 shows
$n$-th decay time, $\tau_{1/n}(T)$ and Fig. 11 shows its compensated plots by
$T^{1-\check{\beta}_{1/n}}$. Here, the exponent $\check{\beta}_{1/n}$ is
selected in the same way as in the previous large $T_{B}$ case. It should be
noted that what we show in Fig.11 is not $D_{n,\check{\beta}_{1/n}}(T)$
defined in Eq.(36), but $\tau_{1/n}(T)/T^{1-\check{\beta}_{1/n}}$. The plotted
variable $\tau_{1/n}(T)/T^{1-\check{\beta}_{1/n}}$ is less dependent of
$T_{B}$ than $D_{n,\check{\beta}_{1/n}}(T)$, which is consistent with the
$T_{B}$ independent behavior of the $n$-th decay time shown in Fig.10.

Figure 10:  Time evolution of the $n$-th decay time scale, $\tau_{1/n}(T)$ for
(a) $n=2$, (b) $n=8$, (c) $n=32$ at $\mathrm{Re}_{\alpha}=160$. The insets
show the same plots as the outsets but at $\mathrm{Re}_{\alpha}=80$.

Figure 11:  Compensated graphs of $\tau_{1/n}(T)$ by $\check{\beta}_{1/n}$ for
(a) $n=2$, (b) $n=8$, and (c) $n=32$ at $\mathrm{Re}_{\alpha}=80$ (dashed) and
$\mathrm{Re}_{\alpha}=160$ (solid). The Insets show $\beta_{1/n}$ as a
function of $T_{B}/T_{\eta}$ at $\mathrm{Re}_{\alpha}=80$ (red) and
$\mathrm{Re}_{\alpha}=160$ (green).

Although the compensated data shown in Fig.11 oscillate getting stronger for
large $n$, we observe a plateau region for each graph. It is noticeable that
the values of the plateaus depend on $Re_{\alpha}$. This tendency is not
present in the large $T_{B}$ case as shown in Fig.6 (ignoring data for the
smallest $T_{B}$). If the self-similarity given in Eq.(9) is valid for
$C_{p}^{L}(T,\tau,T_{B})$, the values of the plateaus should become
independent of $Re_{\alpha}$. Therefore, Fig.11 suggests two possibilities:
one is that $C_{p}^{L}(T,\tau,T_{B})$ is not self-similar; the other is that
$C_{p}^{L}(T,\tau,T_{B})$ is self-similar but with yet another time scale,
$T_{X}$.

The second possibility is more likely, although numerical evidence is
marginally convincing as we will see. With the hypothetical time scale
$T_{X}$, a similarity variable for $C_{p}^{L}(T,\tau,T_{B})$ can be made as
$\tau/[T_{X}^{\check{\beta}}T^{1-\check{\beta}}]$. Therefore, the similarity
function $g^{L}$ under the small $T_{B}$ conditions is likely the similar form
to Eq. (38). The difference is that we just replace $T_{B}$ in the similarity
variable by $T_{X}$. Moreover, as we discussed with Fig.11, the levels of the
plateaus of the vertical axis, which we denote $E_{n,\check{\beta}_{1/n}}(T)$
are independent of $T_{B}$. Regarding change in the numerical values of the
plateau levels as we vary $n$, we observe
$E_{2,\check{\beta}_{1/2}}\simeq(1/3)E_{8,\check{\beta}_{1/8}}\simeq(1/5)E_{32,\check{\beta}_{1/32}}$,
yielding $\tau_{1/n}(T)=\tau_{1/2}(T)\log_{2}n$. This implies that
$C_{p}^{L}(T,\tau,T_{B})$ for $T_{B}<T_{\eta}$ decays exponentially in all the
range of $\tau$. Therefore, $g^{L}$ under the small $T_{B}$ condition is
likely to have two self-similar forms such as

$g^{L}\left(\frac{\tau}{T_{X}^{\beta}T^{1-\beta}}\right)=\begin{cases}\displaystyle\exp\left[{-\check{k}_{1}\left(\frac{\tau}{T_{X}^{\check{\beta}_{1}}T^{1-\check{\beta}_{1}}}\right)}\right]&\text{for}\quad\tau\lesssim
T_{B},\\\ &\\\
\displaystyle\exp\left[{-\check{k}_{2}\left(\frac{\tau}{T_{X}^{\check{\beta}_{2}}T^{1-\check{\beta}_{2}}}\right)}\right]&\text{for}\quad\tau>T_{B},\end{cases}$
(41)

where $\check{\beta}_{1}$ and $\check{\beta}_{2}$ are exponents for the two
self-similar regimes of $\tau$ and $\check{k}_{1}$ and $\check{k}_{2}$ are
constants. We consider that the exponent $\check{\beta}_{1}$ for small $\tau$
is represented by $\beta_{1/2}$ and $\check{\beta}_{2}$ for large $\tau$ is by
$\beta_{1/8}\simeq\beta_{1/32}$. Equation (41) is analogous to Eq.(38) under
the large $T_{B}$ condition. We cannot find a simple form of
$\check{\beta}_{1/n}$ unlike the large $T_{B}$ case.

Figure 12:  Normalized correlation function, $C_{p}^{L}(T,\tau,\varepsilon)$,
with $T_{B}=0.47T_{\eta}$ as a function of rescaled variable (a)
$\tau/[T_{B}^{1-\check{\beta}_{1/2}}~{}T^{\check{\beta}_{1/2}}]$ and (b)
$\tau/[T_{B}^{1-\check{\beta}_{1/32}}~{}T^{\check{\beta}_{1/{32}}}]$ under the
small $T_{B}$ condition. Here the exponents are $\check{\beta}_{1/2}=0.82$ and
$\check{\beta}_{1/{32}}=0.96$. Different curves have different $T$s in
$T_{B}<T<0.54T_{L}$ at $\mathrm{Re}_{\alpha}=160$. The black dashed lines
correspond to (a)
$\exp(-2.7|\tau|/[T_{B}^{\check{\beta}_{1/2}}~{}T^{1-\beta_{1/2}}])$ and (b)
$\exp(-4.8|\tau|/[T_{B}^{\check{\beta}_{1/32}}~{}T^{1-\beta_{1/32}}])$.

We show now that the two different scaling behaviors given in Eq.(41) are
consistent to the DNS data. Figure 12 shows $C_{p}^{L}(T,\tau,\varepsilon)$ as
a functions of $\tau/T^{\check{\beta}_{1/n}}$ for $n=2$ and $n=32$. More
precisely, since we do not know $T_{X}$, we use $T_{B}$ to non-dimensionalize
$\tau/T^{\check{\beta}_{1/n}}$ in the horizontal axis of Fig.12. All the
curves in each panel collapse to one curve with a suitable rescaling of
$\tau$. The similarity function $g^{L}$ can be fitted with an exponential
curve as shown in Fig. 12.

Let us summarize the results in this section. We have considered numerically
scaling behavior of $C^{L}(r_{0},T,\tau,\varepsilon)$ in comparison with the
ansatz (9) under two conditions, $T_{B}>T_{\eta}$ and $T_{B}<T_{\eta}$. The
DNS data of $C^{L}(r_{0},T,\tau,\varepsilon)$ are consistent to the ansatz for
both conditions. The difference between the two conditions is in the
functional forms of the exponents,
$\gamma\left(\frac{T_{\eta}}{T_{B}},\frac{T_{L}}{T_{B}}\right)$ (Eqs.(33) and
(39)), and $\beta\left(\frac{T_{\eta}}{T_{B}},\frac{T_{L}}{T_{B}}\right)$
(Eq.(37)), although $\beta$ for the latter case was not identified.
Furthermore, $\gamma$ and $\beta$ are probably continuous at $T_{B}=T_{\eta}$.
In particular, under the large $T_{B}$ condition, our DNS data indicate that
$\gamma$ and $\beta$ do not approach zero as $T_{B}\to T_{\eta}$. If this is
not a finite Reynolds number effect, the non-zero limits of the exponents
imply that the TTLVCF has exponents that deviate from the K41 dimensional
analysis. Consequently, the squared relative separation $\langle
r^{2}(t)\rangle$ disagrees with the Richardson–Obukhov law $t^{3}$, even if we
take $r_{0}\to 0$ at infinite Reynolds number. We will discuss this point in
the next section.

##  IV Implications on the Richardson–Obukhov law

Finally, we consider implications of the above scaling behaviors of
$C^{L}(r_{0},T,\tau,\varepsilon)$ on the relative separations of particle
pairs. The second moment of the relative separation $r(t)$ can be reduced to

$\displaystyle\langle r^{2}(t)\rangle$ $\displaystyle=$ $\displaystyle
r_{0}^{2}+2\int_{0}^{t}\bm{r_{0}}\cdot\langle\delta\bm{v}(t_{1})\rangle~{}dt_{1}+\int_{0}^{t}\int_{0}^{t}\langle\delta\bm{v}(t_{1})\delta\bm{v}(t_{2})\rangle~{}dt_{1}dt_{2}$
(42) $\displaystyle\sim$ $\displaystyle
r_{0}^{2}+2\int_{0}^{t/2}~{}dT\int_{0}^{2T}~{}d\tau~{}C^{L}(r_{0},T,\tau,\varepsilon)+2\int_{t/2}^{t}~{}dT\int_{0}^{2(t-T)}~{}d\tau~{}C^{L}(r_{0},T,\tau,\varepsilon)$

where the average time is $T=(t_{1}+t_{2})/2$ and the relative time is
$\tau=t_{2}-t_{1}$. We also assume
$\bm{r}_{0}\cdot\langle\delta\bm{v}(t)\rangle=0$ by taking the direction of
the initial separation $\bm{r}_{0}$ being randomly and isotropically
distributed, and use the symmetry of $C^{L}(r_{0},T,\tau,\varepsilon)$ with
respect to the diagonal line $t_{1}=t_{2}$.

First, we consider the scaling law of $\langle r^{2}(t)\rangle$ under the
large $T_{B}$ condition, $T_{\eta}\ll T_{B}\ll T_{L}$, at sufficiently large
Reynolds numbers. Under this condition, $C^{L}(r_{0},T,\tau,\varepsilon)$ has
the self-similar forms (9) and (10). The exponent $\gamma$ is given by
Eq.(33). The self-similar function $g^{L}$ and the other exponent $\beta$ take
two different forms as given in Eqs.(37)–(38), depending on $\tau\lesssim
T_{B}$ or $\tau\gg T_{B}$. More precisely, taking into consideration the self-
similarity, $g^{L}$ may take the following forms,

$g^{L}\left(\frac{\tau}{T_{B}^{\beta}T^{1-\beta}}\right)=\begin{cases}\displaystyle
g^{L}_{1}\left(\frac{\tau}{T_{B}^{\beta_{1}}T^{1-\beta_{1}}}\right)&\text{for}\quad\tau\leq
c_{1}T_{B}\left(\frac{T}{T_{B}}\right)^{1-\beta_{1}},\\\ &\\\ \displaystyle
g^{L}_{2}\left(\frac{\tau}{T_{B}^{\beta_{2}}T^{1-\beta_{2}}}\right)&\text{for}\quad\tau\geq
c_{2}T_{B}\left(\frac{T}{T_{B}}\right)^{1-\beta_{2}},\end{cases}$ (43)

where $c_{1}$ and $c_{2}$ are constants which determine the time to switch
from $g^{L}_{1}$ to $g^{L}_{2}$. These constants are dependent on Reynolds
number. The transition time is around $T_{B}$, that is,
$c_{1}T_{B}(T/T_{B})^{1-\beta_{1}}\sim
T_{B},~{}c_{2}T_{B}(T/T_{B})^{1-\beta_{2}}\sim T_{B}$ for $T_{B}\ll T\ll
T_{L}$. From this, we can estimate the values of $c_{1}$ and $c_{2}$:

$c_{1}\sim\left(\frac{T_{B}}{T_{L}}\right)^{1-\beta_{1}},\quad
c_{2}\sim\left(\frac{T_{B}}{T_{L}}\right)^{1-\beta_{2}}.$ (44)

Therefore, $c_{1}$ and $c_{2}$ are very small constants at sufficiently large
Reynolds numbers. Now, we substitute the self-similar forms (43) and calculate
the integrals under the condition $t\gg T_{B}$, then we obtain,

$\displaystyle\langle r^{2}(t)\rangle\sim$
$\displaystyle~{}\left(\int_{0}^{c_{1}}g_{1}^{L}(x)dx\right)t^{3-\gamma-\beta_{1}}+\left(\int_{c_{2}}^{\infty}g_{2}^{L}(x)dx\right)t^{3-\gamma-\beta_{2}}$
(45)
$\displaystyle+\int_{T_{B}}^{t/2}dT~{}T^{2-\gamma-\beta_{2}}\int_{\infty}^{2(T/T_{B})^{\beta_{2}}}g^{L}_{2}(x)~{}dx+\int_{t/2}^{t}dT~{}T^{2-\gamma-\beta_{2}}\int_{\infty}^{\frac{2(t-T)}{T_{B}^{\beta_{2}}T^{1-\beta_{2}}}}g^{L}_{2}(x)~{}dx.$
(46)

We should consider which term becomes dominant at large Reynolds numbers. To
calculate it in more detail, we use the functional forms
$g^{L}_{1}(x)=e^{-k_{1}x}$ and $g^{L}_{2}(x)=e^{-k_{2}x}$ observed in Fig.8.
Thereby, we can calculate Eqs.(45)-(46) as follows:

$\displaystyle\langle r^{2}(t)\rangle\sim$
$\displaystyle~{}(1-e^{-c_{1}})t^{3-\gamma-\beta_{1}}+e^{-c_{2}}t^{3-\gamma-\beta_{2}}$
$\displaystyle-\Gamma\left(\frac{3-\gamma-\beta_{2}}{\beta_{2}},2k_{2}\left(t/2T_{B}\right)^{\beta_{2}}\right)-t^{3-\gamma-\beta_{2}}\int_{0}^{1/2}(1-z)^{3-\gamma-\beta_{2}}\exp\left[-\frac{2k_{2}t^{\beta_{2}}z}{T_{B}^{\beta_{2}}(1-z)^{1-\beta_{2}}}\right]dz.$
(47)

Here, $\Gamma(a,x)$ is the upper incomplete gamma function defined by
$\Gamma(a,x)=\int_{x}^{\infty}z^{a-1}e^{-z}dz$.

Now we consider conditions to recover the Richardson-Obukhov law, $\langle
r^{2}(t)\rangle\propto t^{3}$. It is known that the upper incomplete gamma
function, $\Gamma(a,x)$, has the following asymptotic series: $\Gamma(a,x)\sim
x^{a-1}e^{-x}[1+\frac{a-1}{x}+\frac{(a-1)(a-2)}{x^{2}}+\cdots]$ as
$x\to\infty$ [49]. With these asymptotic formulae, we have
$\Gamma((3-\gamma-\beta_{2})/\beta_{2},2k_{2}(t/2T_{B})^{\beta_{2}})\propto
t^{3-\gamma-\beta_{2}}\exp[-2k_{2}(t/2T_{B})^{\beta_{2}}]$, as $t\to\infty$.
Here we assume $\beta_{2}\not=0$ at infinitely large Reynolds number, i.e.,
$\omega_{2}\not=0$ in Eq.(37). Furthermore, the last term of Eq.(47) can be
estimated as,

$t^{3-\gamma-\beta_{2}}\int_{0}^{1/2}(1-z)^{3-\gamma-\beta_{2}}\exp\left[-\frac{2k_{2}t^{\beta_{2}}z}{T_{B}^{\beta_{2}}(1-z)^{1-\beta_{2}}}\right]\leq
Ct^{3-\gamma-2\beta_{2}},$ (48)

where $C$ is a constant. Therefore, considering $\beta_{1}>\beta_{2}$ as
obtained in DNS, the dominant power-law scaling at large $t$ and large
$Re_{\alpha}$ is given by

$\langle r^{2}(t)\rangle\sim t^{3-\gamma-\beta_{2}}.$ (49)

Here the power-law exponent of $\langle r^{2}(t)\rangle$ is related to those
of the TTLVCF $C^{L}(r_{0},T,\tau,\varepsilon)$. In particular, it involves
$\beta_{2}$, which implies that $\langle r^{2}(t)\rangle$ is affected by the
correlation $C^{L}(r_{0},T,\tau,\varepsilon)$ far from the diagonal line.

We have found the empirical form of $\gamma$ as a function of
$T_{B},T_{\eta}$, and $T_{L}$, which is given in Eq.(33). It suggests that
$\gamma\to 0$ as $Re_{\alpha}\to\infty$. The similar form of $\beta_{2}$ given
in Eq.(37) indicates that $\beta_{2}\to\omega_{2}^{0.4}$ as
$Re_{\alpha}\to\infty$. As we discussed in Sec.III, with our DNS data we are
not able to conclude whether $\omega_{2}$ vanishes or not. However, at the
practically accessible Reynolds numbers, $\beta_{2}$ is not zero as indicated
by Fig.7. Therefore, now including the constant factor, the Richardson–Obukhov
law is modified at finite Reynolds numbers to

$\langle
r^{2}(t)\rangle=\frac{2G\varepsilon}{k_{2}(3-\gamma-\beta_{2})}T_{B}^{\gamma+\beta_{2}}t^{3-\gamma-\beta_{2}}+\mathrm{(subleading~{}terms)},$
(50)

where we take $c_{2}\sim 0$. Although $G$ and $k_{2}$ are slightly dependent
on $T_{\eta},T_{L},$ and $T_{B}$, they are estimated as $G\varepsilon\sim 1.5$
and $k_{2}\sim 3.0$ from the DNS data.

Now let us consider the numerical value of the Richardson constant $g_{R}$
involved in the Richardson–Obukhov law $\langle
r^{2}(t)\rangle=g_{R}\varepsilon t^{3}$. To evaluate $g_{R}$, we substitute
the values of $G,\varepsilon$ and $k_{2}$ in Eq.(50) by assuming that they do
not change much at infinite Reynolds number. We also assume that the modified
exponent $3-\gamma-\beta_{2}$ approach $3$ at infinite Reynolds number, i.e.,
$\gamma\to 0$ and $\beta_{2}\to 0$ as $Re_{\alpha}\to\infty$. Then the
Richardson constant is estimated as $g_{R}=2G/k_{2}\sim 2\times 10^{1}$, which
is distinctly different from $g_{R}=0.5$ and $3.8$ obtained in previous
experimental and numerical studies, respectively[27, 32]. This discrepancy of
$g_{R}$ is not surprising since the measurements in the previous studies were
done under the small $T_{B}$ condition.

Figure 13:  Mean-squared relative separation of particle pairs, $\langle
r^{2}(t)\rangle$, for various initial separations at
$\mathrm{Re}_{\alpha}=160$ (colored lines) and the scaling law (50) (black
lines). Two vertical dashed lines show the dissipation time scale $T_{\eta}$
(left) and the integral scale $T_{L}$ (right).

Here, we validate the scaling law (50) by our DNS data. Figure 13 shows
$\langle r^{2}(t)\rangle$ as a function of $t$ for various initial separations
at $\mathrm{Re}_{\alpha}=160$ with the scaling laws (50). The scaling law (50)
roughly coincides with the DNS data, especially at larger $T_{B}$. Strictly
speaking, the values of scaling exponents $3-\gamma-\beta_{2}$ in Eq. (50) are
slightly larger than those of DNS data. This is because we assumed that the
self-similar function $g^{L}$ takes the forms (43) at finite Reynolds numbers.
Contrary to this assumption, self-similar function $g^{L}$ may take more
complicated forms at finite Reynolds numbers. In other words, the scaling
exponent $\beta$ may be not binary but continuous function of $\tau$. In this
view, the scaling exponent for $\langle r^{2}(t)\rangle$ is slightly smaller
than $3-\gamma-\beta_{2}$ because $\beta$ may be monotonically decreasing
function of $\tau$.

Although the scaling laws (50) is good agreement with the DNS data for larger
$T_{B}$, this is not the case for smaller $T_{B}$ such as $T_{B}=1.4T_{\eta}$.
This deviation does not mean failure of the scaling laws (50) because we
assume the scale separation, $T_{\eta}\ll T_{B}\ll T_{L}$. The scaling law
(50) is still useful to understand the qualitative tendency of $\langle
r^{2}(t)\rangle$ at smaller $T_{B}$ than $T_{\eta}$ too.

We then consider the scaling law of $\langle r^{2}(t)\rangle$ under the small
$T_{B}$ condition. By doing an analogous calculation to that of the large
$T_{B}$ condition, we have

$\langle
r^{2}(t)\rangle=\frac{2G\varepsilon}{\check{k}_{2}(3-\check{\gamma}-\check{\beta}_{2})}T_{X}^{\check{\gamma}+\check{\beta}_{2}}t^{3-\check{\gamma}-\check{\beta}_{2}}+\mathrm{(subleading~{}terms)},$
(51)

at large $t$ and large $Re_{\alpha}$. This is similar to Eq.(50) for the large
$T_{B}$ condition. A crucial difference between Eqs.(51) and (50) is that the
exponent $\check{\gamma}$ is negative as seen from Eq.(39). This enables one
to tune $T_{B}$ for given $T_{\eta}$ such that
$-\check{\gamma}(T_{B},T_{\eta})-\check{\beta}_{2}(T_{\eta}/T_{B},T_{L}/T_{B})=0$
under the small $T_{B}$ condition (our DNS data suggest that
$\check{\beta}_{2}$ is generally positive). Consequently, we observe $\langle
r^{2}(t)\propto t^{3}$, the same scaling exponent as the Richardson–Obukhov
law. In contrast, this sort of tuning leading to $t^{3}$ is not possible under
the large $T_{B}$ condition since $\gamma$ is always positive, see Eq.(33).

Indeed, it is known that $\langle r^{2}(t)\rangle\propto t^{3}$ can be
observed even at moderate Reynolds number by tuning the initial separation
$r_{0}$ which satisfies the small $T_{B}$ condition $T_{B}\leq T_{\eta}$. See,
for example, Refs. [27, 32, 37, 28, 30, 33] in the 2D energy inverse-cascade
turbulence. Specifically, with the tuned initial separation $r_{0}$, $T_{B}$
is close to $T_{\eta}$. In those circumstances, the equivalent of the
Richardson constant can be given by
$\check{g}_{R}=4GT_{X}^{\check{\gamma}+\check{\beta}_{2}}/(3\check{k}_{2})$
from Eq.(51). It is noted that the value of $G$ in this range is strongly
dependent on $T_{B}$ and $T_{\eta}$. The Richardson constant measured in the
previous experimental and numerical studies [27, 32] with the tuned initial
separation should be therefore compared to $\check{g}_{R}$. However, we do not
compare it quantitatively since we can not determine the accurate value of
$T_{X}$.

Let us now argue that the nature of the $t^{3}$ law with the tuned initial
separation is different from that of the Richardson–Obukhov law. In general
terms, by the Richardson–Obukhov $t^{3}$ law, it is understood that the
$t^{3}$ law holds irrespective of the value of the initial separation $r_{0}$,
provided that the inertial subrange is sufficiently wide. Strictly speaking,
one should add a condition that $r_{0}$ is inside the inertial subrange [21].
The large $T_{B}$ condition which we have considered conforms to the added
condition. From Eq.(50), the Richardson–Obukhov law corresponds to $\gamma=0$
and $\beta_{2}=0$, and the resultant $t^{3}$ law does not depend on $T_{B}$,
or equivalently $r_{0}$. For the sake of the argument, let us relax the added
condition. Now we consider the small $T_{B}$ condition. From Eq.(51), The
$t^{3}$ law with the tuned initial separation corresponds to
$-\check{\gamma}-\check{\beta}_{2}=0$ and the resultant $t^{3}$ law has the
prefactor $T_{X}^{\check{\gamma}+\check{\beta}_{2}}$, which potentially
depends on $T_{B}$. Therefore, the $t^{3}$ law observed at a given Reynolds
number (however large) by tuning the initial separation is different from the
Richardson–Obukhov $t^{3}$ law. The agreement of the power-law exponents is
coincidental. In this sense, the $t^{3}$ scaling of $\langle r^{2}(t)\rangle$
observed at moderate Reynolds numbers is a different state from the complete
similarity for $\langle r^{2}(t)\rangle$, which is consistent with the
dimensional analysis naively using the K41 phenomenology.

It is interesting that such a coincidence do not occur under the large $T_{B}$
condition. Then, in this condition, can we say anything about observability of
the bona fide Richardson–Obukhov law? As far as our DNS data suggest, the
exponents $\gamma$ and $\beta_{2}$ do not vanish under the large $T_{B}$
condition. Consequently, the Richardson–Obukhov law is not observable at the
current Reynolds numbers. It should be noted that this is caused not by the
intermittency effects, but by correlation of the Lagrangian velocity.
Extrapolation of the data suggests a possibility that $\beta_{2}$ does not
vanish at infinite Reynolds number as we discussed. This implies that the no
matter how large the Reynolds number is, the Richardson–Obukhov law is not
observable.

We have presented here a framework to study the Richardson–Obukhov law by way
of the self-similarity of the TTLVCF. It can be adapted to the 3D turbulence.
The $t^{3}$ law with the tuned initial separation is also known in the 3D
case, see, for example, [26, 35, 22, 36, 23]. Our analysis in the 3D case will
be reported elsewhere.

##  V Concluding remarks

We have investigated the two-time Lagrangian velocity correlation function
(TTLVCF) for particle pairs with the incomplete self-similarity and the DNS of
the 2D energy inverse-cascade turbulence. First, we have made the self-similar
ansatz (9) of the correlation function by using the idea of incomplete
similarity. The ansatz includes the Bachelor time, the Kolmogorov dissipation
length, and the integral length as similarity variables, meaning that finite
Reynolds number effects and the initial separation dependence are encoded. The
ansatz is characterized by the two scaling exponents, $\beta$ and $\gamma$,
and the one-variable function $g_{L}$. The exponent $\gamma$ concerns the
equal-time correlation along the diagonal line through the origin shown in
Fig.1. The other exponent $\beta$ concerns how the correlation decreases along
the direction perpendicular to the diagonal line. Since the ansatz is an
example of the incomplete self-similarity, the two exponents cannot be
determined by dimensional analysis.

In order to verify the ansatz, we have performed DNS of the 2D inverse energy-
cascade turbulence and calculated the TTLVCF by varying parameters such as the
Batchelor time $T_{B}$ and the dissipation time $T_{\eta}$. We split the DNS
study into two parts: the large and small $T_{B}$ conditions. For both
conditions, we showed that the ansatz describes the DNS results reasonably
well. Then we measured the values of the two exponents and the functional form
of $g_{L}$ from the DNS data which are in some cases too noisy to obtain
reliable measurements. The measurements indicated that the exponents depend on
$T_{B}$ and $T_{\eta}$. In theory, we assumed that they are independent. This
dependence of the exponents is empirically determined as Eqs.(33), (37), and
(39). The function $g^{L}$ is determined as an exponential function. Moreover,
our results indicate that at finite Reynolds numbers, the correlation in
general has correction described by non-zero $\gamma$, non-zero $\beta$ and
the function $g^{L}$ to the K41 dimensional analysis for both large and small
$T_{B}$ conditions.

We next considered the limit of these empirical relations at infinite Reynolds
number. The extrapolation of the relations obtained at moderate Reynolds
numbers was subject to uncertainty. Thus, we could not determine the accurate
values of the scaling exponents, $\gamma$ and $\beta$ at infinite Reynolds
number. Especially, we observed that $\beta$ may not approach zero, which
suggests a possibility that the TTLVCF is not consistent with the K41
dimensional analysis at infinite Reynolds number under the large $T_{B}$
condition.

Finally, we have considered the relationship between the scaling law of the
TTLVCF and the Richardson–Obukhov $t^{3}$ law for the second moment of
relative separation via the integral (42). With the asymptotic argument, we
found that the Richardson–Obukhov law is not recovered for finite $T_{B}$ at
finite Reynolds numbers under the large $T_{B}$ condition. Instead, the
scaling law of the squared separation of particle pairs is modified as
Eq.(50). The modified scaling exponent, $3-\gamma-\beta_{2}$, is determined by
the scaling behavior of the TTLVCF on the diagonal line $\tau=0$ shown in
Fig.1 and far from the diagonal line. Namely, the scaling exponent is
influenced by correlations between two particles not only at a same time but
also at quite different times. Moreover, using the scaling law of the TTLVCF,
we explained why we, nevertheless, observe the $t^{3}$ scaling at moderate
Reynolds numbers with a special initial separation under the small $T_{B}$
condition because $\gamma$ and $\beta_{2}$ take a negative and positive
values, respectively. Therefore, we concluded that the physics of this
$t^{3}$-scaling behavior is different from that of the Richardson–Obukhov law.

In this paper, we assumed that forcing effects are negligible. The external
forcing is limited to small scales for 2D. In fact, the Eulerian statistics in
the Fourier space such as the energy spectrum or the energy flux is influenced
by the forcing only in the vicinity of the forcing scales[44, 50, 45]. Hence
the influence is considered as local. This may be the reason why the
empirically found functional forms of the scaling exponents, $\beta$ and
$\gamma$ depend only $T_{\eta},T_{L}$ and $T_{B}$ given by Eqs. (37) and (33)
under the large $T_{B}$ condition. Strictly speaking, we can neglect the
forcing effects if correlation between the forcing and the Lagrangian
velocity, $\langle f_{i}\delta v_{j}\rangle$, rapidly decays in time. Here
$f_{i}$ is the forcing increment between two Lagrangian particles and $\delta
v_{j}$ is the relative velocity between them. We speculate that the cross
correlation rapidly decays because the characteristic time scales of the
forcing and the velocity in the forcing scale are small.

Here, we discuss whether or not the scaling laws (10) and the modified
Richardson–Obukhov law (50) can be applied to 3D turbulence. First, the
scaling laws (10) could be suggested by only using the dimensinal analysis,
which can be applied to both 2D and 3D turbulence. Namely, the special
properties of 2D turbulence such as the inverse cascade of energy flux are not
necessary. This makes us convincing that this scaling laws (10) can be
implimeted to 3D turbulence. On the other hand, more careful investigation is
needed for the modified Richardson–Obukhov law (50). This is because this is
derived by using the DNS results such as Eq. (43).

Given the self-similar form Eqs.(9) and (10) of the TTLVCF, one would like to
“derive” it from the Navier–Stokes equations using only plausible assumptions.
More precisely, we propose to use it as an input to a set of integro-
differential equations (closure equations) for the Lagrangian correlation
function obtained by a closure approximation such as direct-interaction
approximation [2, 3, 4, 51]. One standard procedure in the last step of the
closures is to substitute certain self-similar forms for the correlation
function and the linear response function to those in the closure equations
and then to study consistency of the self-similar forms with the closure
equations. By input, we mean to input the ansatz studied here into closure
equations of, for example, a direct-interaction approximation. This may give
analytically functional forms of the scaling exponents,
$\beta(T_{\eta}/T_{B},T_{L}/T_{B})$ and $\gamma(T_{\eta}/T_{B},T_{L}/T_{B})$
and their limits at infinite Reynolds number.

Closure approximations have been applied to study the Richardson–Obukhov
$t^{3}$ law, see e.g., [9, 16]. However, these studies have used one-time
Lagrangian velocity correlation function given by Eq.(2) in Introduction,
which is different from the TTLVCF $C^{L}(r_{0},T,\tau,\varepsilon)$ studied
here. In fact, the TTLVCF is unexplored with the Lagrangian renormalization
approximation [4] and perhaps other closure approximations 111Y. Kaneda
(private communication). Therefore, the results in this study play an
important role to develop a new avenue of closure theories.

Another approach can be to develop a stochastic model of turbulent relative
dispersion using the ansatz we have obtained here. Recently, continuous time
random walk (CTRW) models [53, 54] developed for the relative dispersion.
These models are constructed to be consistent with the Richardson–Obukhov law.
It is possible to modify these models to have the self-similar properties of
the TTLVCF obtained in this paper. Building such a model corresponds to
incorporating effects of time correlations [47, 55] and finite propagation
speed of the relative diffusion [56, 57]. We will report a stochastic modeling
based on the ansatz elsewhere.

###### Acknowledgements.

We are grateful for stimulating discussion with Yukio Kaneda. The numerical
computations in this work were performed at the computer facility of the
Yukawa Institute for Theoreical Physics at Kyoto University. This study was
supported by the Research Institute for Mathematical Sciences at Kyoto
University and by Kakenhi grant (A) No. 19H00641 from JSPS.

## References

  * He  _et al._ [2017] G. He, G. Jin, and Y. Yang, Space-time correlations and dynamic coupling in turbulent flows, Annual Review of Fluid Mechanics 49, 51 (2017).
  * Kraichnan [1959] R. H. Kraichnan, The structure of isotropic turbulence at very high Reynolds numbers, Journal of Fluid Mechanics 5, 497 (1959).
  * Kraichnan [1965] R. H. Kraichnan, Lagrangian-History Closure Approximation for Turbulence, Physics of Fluids 8, 575 (1965).
  * Kaneda [1981] Y. Kaneda, Renormalized expansions in the theory of turbulence with the use of the Lagrangian position function, Journal of Fluid Mechanics 107, 131 (1981).
  * Kolmogorov [1941] A. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Dokl. Akad. Nauk SSSR 30, 301 (1941).
  * Kraichnan [1967] R. H. Kraichnan, Inertial Ranges in Two-Dimensional Turbulence, Physics of Fluids 10, 1417 (1967).
  * Leith [1968] C. E. Leith, Diffusion approximation for two-dimensional turbulence, Physics of Fluids 11, 671 (1968).
  * Batchelor [1969] G. K. Batchelor, Computation of the Energy Spectrum in Homogeneous Two-Dimensional Turbulence, Physics of Fluids 12, II (1969).
  * Kraichnan [1966] R. H. Kraichnan, Dispersion of Particle Pairs in Homogeneous Turbulence, Physics of Fluids 9, 1937 (1966).
  * Gotoh and Kaneda [1991] T. Gotoh and Y. Kaneda, Lagrangian velocity autocorrelation and eddy viscosity in two-dimensional anisotropic turbulence, Physics of Fluids A: Fluid Dynamics 3, 2426 (1991).
  * Gotoh  _et al._ [1993] T. Gotoh, R. S. Rogallo, J. R. Herring, and R. H. Kraichnan, Lagrangian velocity correlations in homogeneous isotropic turbulence, Physics of Fluids A: Fluid Dynamics 5, 2846 (1993).
  * Kaneda  _et al._ [1998] Y. Kaneda, K. Gotoh, and T. Ishihara, Taylor Expansions and Padé Approximations of Lagrangian and Eulerian Two-Time Velocity Correlations in Turbulence, Journal of the Physical Society of Japan 67, 1075 (1998).
  * Kaneda  _et al._ [1999] Y. Kaneda, T. Ishihara, and K. Gotoh, Taylor expansions in powers of time of Lagrangian and Eulerian two-point two-time velocity correlations in turbulence, Physics of Fluids 11, 2154 (1999).
  * He  _et al._ [2009] G.-W. He, G. Jin, and X. Zhao, Scale-similarity model for Lagrangian velocity correlations in isotropic and stationary turbulence, Physical Review E 80, 066313 (2009).
  * Monin and Yaglom [1975] A. S. Monin and A. M. Yaglom,  _Statistical Fluid Mechanics, Mechanics of Turbulence_ , Vol. II (MIT Press, 1975).
  * Ishihara and Kaneda [2002] T. Ishihara and Y. Kaneda, Relative diffusion of a pair of fluid particles in the inertial subrange of turbulence, Physics of Fluids 14, L69 (2002).
  * Barenblatt [2003] G. I. Barenblatt,  _Scaling_ (Cambridge University Press, 2003).
  * Richardson [1926] L. F. Richardson, Atmospheric Diffusion Shown on a Distance-Neighbour Graph, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 110, 709 (1926).
  * Salazar and Collins [2009] J. P. Salazar and L. R. Collins, Two-Particle Dispersion in Isotropic Turbulent Flows, Annual Review of Fluid Mechanics 41, 405 (2009).
  * Obukhov [1941] A. Obukhov, On the distribution of energy in the spectrum of turbulent flow, Izv. Akad. Nauk SSSR, Ser. Geogr. Geofi 5, 453 (1941).
  * Batchelor [1950] G. K. Batchelor, The application of the similarity theory of turbulence to atmospheric diffusion, Quarterly Journal of the Royal Meteorological Society 76, 133 (1950).
  * Sawford  _et al._ [2008] B. L. Sawford, P. K. Yeung, and J. F. Hackl, Reynolds number dependence of relative dispersion statistics in isotropic turbulence, Physics of Fluids 20, 065111 (2008).
  * Buaria  _et al._ [2015] D. Buaria, B. L. Sawford, and P. K. Yeung, Characteristics of backward and forward two-particle relative dispersion in turbulence at different Reynolds numbers, Physics of Fluids 27, 105101 (2015).
  * Buaria  _et al._ [2016] D. Buaria, P. K. Yeung, and B. L. Sawford, A Lagrangian study of turbulent mixing: forward and backward dispersion of molecular trajectories in isotropic turbulence, Journal of Fluid Mechanics 799, 352 (2016).
  * Ouellette  _et al._ [2006] N. T. Ouellette, H. Xu, M. Bourgoin, and E. Bodenschatz, An experimental study of turbulent relative dispersion models, New Journal of Physics 8, 109 (2006).
  * Ott and Mann [2000] S. Ott and J. Mann, An experimental investigation of the relative diffusion of particle pairs in three-dimensional turbulent flow, J. Fluid Mech 422, 207 (2000).
  * Jullien  _et al._ [1999] M.-C. Jullien, J. Paret, and P. Tabeling, Richardson Pair Dispersion in Two-Dimensional Turbulence, Physical Review Letters 82, 2872 (1999).
  * Rivera and Ecke [2005] M. K. Rivera and R. E. Ecke, Pair dispersion and doubling time statistics in two-dimensional turbulence, Phys. Rev. Lett. 95, 194503 (2005).
  * von Kameke  _et al._ [2011] A. von Kameke, F. Huhn, G. Fernández-García, A. P. Muñuzuri, and V. Pérez-Muñuzuri, Double cascade turbulence and richardson dispersion in a horizontal fluid flow induced by faraday waves, Phys. Rev. Lett. 107, 074502 (2011).
  * Rivera and Ecke [2016] M. K. Rivera and R. E. Ecke, Lagrangian statistics in weakly forced two-dimensional turbulence, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 013103 (2016).
  * Boffetta and Sokolov [2002a] G. Boffetta and I. M. Sokolov, Relative Dispersion in Fully Developed Turbulence: The Richardson’s Law and Intermittency Corrections, Physical Review Letters 88, 094501 (2002a).
  * Boffetta and Sokolov [2002b] G. Boffetta and I. M. Sokolov, Statistics of two-particle dispersion in two-dimensional turbulence, Physics of Fluids 14, 3224 (2002b).
  * Kishi  _et al._ [2020] T. Kishi, T. Matsumoto, and S. Toh, Non-Kolmogorov scaling for two-particle relative velocity in two-dimensional inverse energy-cascade turbulence, Physical Review Fluids 5, 054601 (2020).
  * YEUNG and BORGAS [2004] P. K. YEUNG and M. S. BORGAS, Relative dispersion in isotropic turbulence. Part 1. Direct numerical simulations and Reynolds-number dependence, Journal of Fluid Mechanics 503, 93 (2004).
  * Biferale  _et al._ [2005] L. Biferale, G. Boffetta, A. Celani, B. J. Devenish, A. Lanotte, and F. Toschi, Lagrangian statistics of particle pairs in homogeneous isotropic turbulence, Physics of Fluids 17, 1 (2005).
  * Bitane  _et al._ [2013] R. Bitane, H. Homann, and J. Bec, Geometry and violent events in turbulent pair dispersion, Journal of Turbulence 14, 23 (2013).
  * Kellay and Goldburg [2002] H. Kellay and W. I. Goldburg, Two-dimensional turbulence: a review of some recent experiments, Rep. Prog. Phys. 65, 845 (2002).
  * Paret and Tabeling [1997] J. Paret and P. Tabeling, Experimental Observation of the Two-Dimensional Inverse Energy Cascade, Physical Review Letters 79, 4162 (1997).
  * Boffetta  _et al._ [2000] G. Boffetta, A. Celani, and M. Vergassola, Inverse energy cascade in two-dimensional turbulence: Deviations from gaussian behavior, Phys. Rev. E 61, R29 (2000).
  * Barenblatt [2014] G. I. Barenblatt,  _Flow, Deformation and Fracture: Lectures on Fluid Mechanics and the Mechanics of Deformable Solids for Mathematicians and Physicists_, Cambridge Texts in Applied Mathematics (Cambridge University Press, 2014).
  * Kraichnan [1964] R. H. Kraichnan, Kolmogorov’s Hypotheses and Eulerian Turbulence Theory, Physics of Fluids 7, 1723 (1964).
  * Tennekes [1975] H. Tennekes, Eulerian and Lagrangian time microscales in isotropic turbulence, Journal of Fluid Mechanics 67, 561 (1975).
  * Wallace [2014] J. M. Wallace, Space-time correlations in turbulent flow: A review, Theoretical and Applied Mechanics Letters 4, 022003 (2014).
  * XIAO  _et al._ [2009] Z. XIAO, M. WAN, S. CHEN, and G. L. EYINK, Physical mechanism of the inverse energy cascade of two-dimensional turbulence: a numerical investigation, Journal of Fluid Mechanics 619, 1 (2009).
  * Mizuta  _et al._ [2013] A. Mizuta, T. Matsumoto, and S. Toh, Transition of the scaling law in inverse energy cascade range caused by a nonlocal excitation of coherent structures observed in two-dimensional turbulent fields, Physical Review E 88, 053009 (2013).
  * Vallgren [2011] A. Vallgren, Infrared Reynolds number dependency of the two-dimensional inverse energy cascade, Journal of Fluid Mechanics 667, 463 (2011).
  * Scatamacchia  _et al._ [2012] R. Scatamacchia, L. Biferale, and F. Toschi, Extreme Events in the Dispersions of Two Neighboring Particles Under the Influence of Fluid Turbulence, Physical Review Letters 109, 144501 (2012).
  * Biferale  _et al._ [2014] L. Biferale, A. S. Lanotte, R. Scatamacchia, and F. Toschi, Intermittency in the relative separations of tracers and of heavy particles in turbulent flows, Journal of Fluid Mechanics 757, 550 (2014).
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York, 1964).
  * Boffetta and Musacchio [2010] G. Boffetta and S. Musacchio, Evidence for the double cascade scenario in two-dimensional turbulence, Physical Review E 82, 016307 (2010).
  * Kida and Goto [1997] S. Kida and S. Goto, A Lagrangian direct-interaction approximation for homogeneous isotropic turbulence, Journal of Fluid Mechanics 345, 307 (1997).
  * Note [1] Y. Kaneda (private communication).
  * Thalabard  _et al._ [2014] S. Thalabard, G. Krstulovic, and J. Bec, Turbulent pair dispersion as a continuous-time random walk, Journal of Fluid Mechanics 755, R4 (2014).
  * Bourgoin [2015] M. Bourgoin, Turbulent pair dispersion as a ballistic cascade phenomenology, Journal of Fluid Mechanics 772, 678 (2015).
  * Eyink and Benveniste [2013] G. L. Eyink and D. Benveniste, Diffusion approximation in turbulent two-particle dispersion, Phys. Rev. E 88, 041001(R) (2013).
  * Ogasawara and Toh [2006] T. Ogasawara and S. Toh, Model of Turbulent Relative Dispersion: A Self-Similar Telegraph Equation, Journal of the Physical Society of Japan 75, 083401 (2006).
  * Kanatani  _et al._ [2009] K. Kanatani, T. Ogasawara, and S. Toh, Telegraph-Type versus Diffusion-Type Models of Turbulent Relative Dispersion, Journal of the Physical Society of Japan 78, 024401 (2009).