Protoplasmic streaming in plant cells is directly visible in the cases of *Chara corallina* and *Nitella flexilis*, and this streaming is understood to play a role in the transport of biological materials. For this reason, related studies have focused on molecular transportation from a fluid mechanics viewpoint. However, the experimentally observed distribution of the velocity along the flow direction *x*, which exhibits two peaks at *V*_{x} = 0 and at a finite *V*_{x}(≠0), remains to be studied. In this paper, we numerically study whether this behavior of the flow field can be simulated by a 2D stochastic Navier–Stokes (NS) equation for Couette flow in which a random Brownian force is assumed. We present the first numerical evidence that these peaks are reproduced by the stochastic NS equation, which implies that the Brownian motion of the fluid particles plays an essential role in the emergence of these peaks in the velocity distribution. We also find that the position of the peak at *V*_{x}(≠0) moves with the variation in the strength *D* of the random Brownian force, which also changes depending on physical parameters such as the kinematic viscosity, boundary velocity, and diameter of the plant cells.

## I. INTRODUCTION

A circular flow called protoplasmic streaming is directly observed in the cells of specific plants, such as *Chara corallina* and *Nitella flexilis*, in which the cell size is relatively large, ranging from a few hundred micrometers to 1 mm.^{1–5} Such streaming inside cells is naturally considered to play a role in the transportation of biological materials.

The driving force of the flow is known to be molecules moving along actin filaments; hence, these molecules are called molecular motors.^{6–11} These molecular motors transport chlorophyll, which is very large, and drive the flow in plant cells. Interestingly, the speed of the flow in cells is closely related to the size of the plant.^{12} Moreover, the mechanism for the transportation of biological materials is understood to be the same as that in animal cells.^{13,14}

Therefore, protoplasmic streaming has attracted substantial attention both in scientific fields and in the context of agricultural technology.^{12} Kamiya and Kuroda observed the position dependence of the flow speed in a section vertical to the longitudinal direction of *Nitella* cells via optical microscopy in 1956^{5} [Fig. 1(a)]. This position dependence of the flow speed was later precisely measured via particle tracking velocimetry by Kikuchi-Mochizuki^{15} who reported results compatible with the simulation data obtained by Goldstein *et al.* using coupled Navier–Stokes (NS) and advection–diffusion equations.^{16,17} Goldstein *et al.* assumed a spiral flow as a boundary condition on the wall in their simulations and provided insight into the role of this spiral flow in molecular transportation. This velocity field was later shown to be compatible with experimental data obtained through magnetic resonance velocimetry.^{18,19} Niwayama *et al.* also simulated streaming in the case of *Chara corallina* using the 3D NS equation in a method called the moving particle semi-implicit method, in which the spiral flow is neglected, and reported results almost identical to those of Goldstein *et al.*^{20} Their original motivation was to use particle image velocimetry to measure the streaming velocity in the case of *Caenorhabditis elegans* embryos in which the mechanism of transportation is slightly different from that in the cases of *Chara corallina* and *Nitella flexilis*.^{20}

In 1974, Mustacich and Ware observed the distribution of the velocity *V*_{x} along the flow direction *x* by means of a laser light scattering technique called laser Doppler velocimetry and found two different peaks in the velocity distribution: one at *V*_{x} = 0 and the other at a finite *V*_{x}(≠0)^{21–23} [Fig. 1(b)]. Shortly thereafter, the velocity distribution was again measured using the same technique by Sattelle and Buchan^{24} who similarly detected two different peaks. Figures 2(a) and 2(b) show the experimental data extracted from Refs. 21–23, where the data are represented by solid lines approximating the data points. The horizontal axis represents the frequency of the laser light, which is proportional to the fluid velocity, and the shape and position of the second peak depend on the scattering angle; note that the corresponding velocities in these figures are identical to each other for data obtained at the same point inside the cell. The velocity of the second peak was reported to be 60 *μ*m/s^{21} and 72 *μ*m/s^{22,23} for the data shown in Figs. 2(a) and 2(b), respectively. The biological implications of the existence of the second peak are currently unclear; however, it is possible that the peak in a relatively high-velocity region may be closely related to an enhancement of some biological function, such as transportation or mixing.

However, the experimentally obtained velocity distribution has not yet been numerically verified. Although the peak at zero velocity is expected to be caused by the Brownian motion of the fluid molecules, as noted in Ref. 24, the peak at a finite velocity has yet to be explained. Clearly, a microscopic perspective is effective for studying this problem; therefore, to this end, we adopt Langevin simulation, which is a technique for simulating the Brownian motion of small particles.^{25–32}

The peak in the velocity distribution at zero velocity can be naturally understood from the fact that the fluid at the central part of a cell is expected to have a slow speed compared to the fluid at the wall;^{16,17} thus, the fluid in the central region is expected to be influenced by random Brownian forces. In contrast, the fluid close to the cell wall is strongly influenced by the activation forces of molecular motors; in other words, thermal fluctuations are suppressed by contact with the motors and the cell wall. On the other hand, fluid that is separated from the cell wall is not influenced by such boundary conditions, and the speed of the fluid is expected to continuously decrease toward the central region of the cell. Therefore, no intuitive explanation is available for the existence of the second peak at a relatively high velocity.

In this paper, we numerically solve the NS equation with a random Brownian force for flow fields in a square region by regarding twisting flows as straight flows along the longitudinal direction. This 2D NS equation is considered as a Langevin equation or a stochastic differential equation because it includes a random force. In this paper, we combine two different techniques:^{33} NS simulation for continuum fluids^{34,35} and Langevin simulation for particles.^{25–32} This simulation approach in combination with the NS equation is new and, thus, is not comparable to standard techniques for the NS equation without random forces; therefore, we carefully check the dependence of the results on parameters including spatial and temporal discretizations.

It will be shown that all qualitatively different simulation results can be obtained by merely varying the strength *D* of the random Brownian force and that two different peaks appear in the velocity distribution at intermediate values of *D*.

## II. STOCHASTIC NAVIER–STOKES EQUATION

The symbols for the variables and constants used in this paper are listed along with their units and descriptions in Table I.

Symbol . | Unit . | Description (assumed typical value) . |
---|---|---|

ψ | m^{2}/s | Stream function |

ω | 1/s | Vorticity |

$V\u2192\u2009=\u2009(Vx,Vy)$ | m/s | Velocity vector |

V_{B} | m/s | Velocity at the boundary (50 × 10^{−6}) |

$\eta \u2192\u2009=\u2009(\eta x,\eta y)$ | m/s^{2} | Gaussian random force |

$g\u2192\u2009=\u2009(gx,gy)$ | 1 | Gaussian random number of mean 0 and deviation 1 |

D | m^{2}/s^{3} | Strength of random force (2 × 10^{−7} ≤ D ≤ 4 × 10^{−6}) |

D_{dif} | m^{2}/s | Diffusion constant |

τ_{e} | s | Macroscopic relaxation time |

μ | Ns/m^{2} | Viscosity |

ν | m^{2}/s | Kinematic viscosity (1 × 10^{−4}) |

n_{X} | 1 | Total number of lattice points on one edge (100 ≤ n_{X} ≤ 300) |

n_{T} | 1 | Total number of iterations per time step |

α | 1 | Positive number for a unit change in length |

β | 1 | Positive number for a unit change in time |

γ | 1 | Positive number for a change in n_{X} |

δ | 1 | Positive number for a change in n_{T} |

d | m | Cell diameter (500 × 10^{−6}) |

Δx | m | Lattice spacing |

Δt | s | Discrete time step (8 × 10^{−9}) |

E | Parameter set (ν, V, D) | |

S | Parameter set (ν, V, D, Δx, Δt) | |

a | m | Diameter of a lump of fluid particles |

Exp(E) | Experimental data corresponding to E = (ν, V, D) | |

R_{e} | 1 | Reynolds number |

P_{e} | 1 | Péclet number |

S_{c} | 1 | Schmidt number |

Symbol . | Unit . | Description (assumed typical value) . |
---|---|---|

ψ | m^{2}/s | Stream function |

ω | 1/s | Vorticity |

$V\u2192\u2009=\u2009(Vx,Vy)$ | m/s | Velocity vector |

V_{B} | m/s | Velocity at the boundary (50 × 10^{−6}) |

$\eta \u2192\u2009=\u2009(\eta x,\eta y)$ | m/s^{2} | Gaussian random force |

$g\u2192\u2009=\u2009(gx,gy)$ | 1 | Gaussian random number of mean 0 and deviation 1 |

D | m^{2}/s^{3} | Strength of random force (2 × 10^{−7} ≤ D ≤ 4 × 10^{−6}) |

D_{dif} | m^{2}/s | Diffusion constant |

τ_{e} | s | Macroscopic relaxation time |

μ | Ns/m^{2} | Viscosity |

ν | m^{2}/s | Kinematic viscosity (1 × 10^{−4}) |

n_{X} | 1 | Total number of lattice points on one edge (100 ≤ n_{X} ≤ 300) |

n_{T} | 1 | Total number of iterations per time step |

α | 1 | Positive number for a unit change in length |

β | 1 | Positive number for a unit change in time |

γ | 1 | Positive number for a change in n_{X} |

δ | 1 | Positive number for a change in n_{T} |

d | m | Cell diameter (500 × 10^{−6}) |

Δx | m | Lattice spacing |

Δt | s | Discrete time step (8 × 10^{−9}) |

E | Parameter set (ν, V, D) | |

S | Parameter set (ν, V, D, Δx, Δt) | |

a | m | Diameter of a lump of fluid particles |

Exp(E) | Experimental data corresponding to E = (ν, V, D) | |

R_{e} | 1 | Reynolds number |

P_{e} | 1 | Péclet number |

S_{c} | 1 | Schmidt number |

### A. Discretization of the stochastic Navier–Stokes equation for the stream function

First, we present the computational domain extracted from the cylindrical body shown in Fig. 3(a). The arrows in Fig. 3(a) illustrate the flow directions in the indifferent zone indicated by the dashed line on the surface. We consider the section *AA*′–*CC*′, which is also visualized in Fig. 3(b), with arrows on the boundaries *AA*′ and *CC*′, where the flow direction is modified to be parallel/antiparallel along the longitudinal direction. The arrows inside the square represent the velocity on the surface of the opposite side, which is not included in or differs from the square domain. The computational domain is the flat square region with boundaries *AA*′ and *CC*′, where the other two boundaries, *AC* and *A*′*C*′, are assumed to be periodic. This square region is extracted and shown in Fig. 3(c) to clarify the 2D nature of the domain.

The boundary condition given by the velocity *V*_{B} is simply the same as that for Couette flow. The real 3D flow is modified to this 2D flow for simplicity, and the flow direction in the 2D domain is obtained by modifying the flow direction on the surface of the cylinder, as stated above. Because of this modification of the velocity direction, we can determine whether the origin of the peaks in the velocity distribution lies in the spiral flow. It is also possible to investigate whether the peaks are related to the 3D nature of the flow.

Here, we should comment on the reason why the shape of the boundary *AA*′–*CC*′ is assumed to remain unchanged. Cell surfaces composed of soft biological materials may exhibit shape deformations that can be directly measured, for example, in the case of *Caenorhabditis elegans* embryos in which cytoplasmic streaming is also expected.^{20} However, in the case of plants such as *Chara corallina* and *Nitella flexilis*, the situation is different; the cell surface is relatively hard, and fluctuations can be neglected.

The continuous form of the NS equation^{34,35} with a random Brownian force is given by

where $V\u2192\u2009=\u2009(Vx,Vy,0)$ is the fluid velocity obtained from the stream function *ψ* and *ω* is the third component of the vorticity $\omega \u2192\u2009=\u2009\u2207\xd7V\u2192$ such that

The physical meaning of each term in Eq. (1) will be given below in reference to the equation for the velocity $V\u2192$. The parameter *ν* in Eq. (1) is the kinematic viscosity coefficient, where *ν* ≃ 1 × 10^{−6} m^{2}/s in the case of water at room temperature. The symbol $\eta \u2192(t)=(\eta x,\eta y,0)$ represents Gaussian white noise or a Gaussian random force corresponding to the Brownian motion of the fluid particles or a lump of fluid particles. The components of $\eta \u2192(t)$ are assumed to satisfy

where ⟨⋯⟩ denotes the expectation value, *D* is called the strength of the random force, and the subscript *i* denotes the fluid position. In Eq. (3), we introduce $\eta \u2192(t)$ in a discrete form because the NS equation is discretized on a square lattice in this numerical study. No confusion is expected between the symbol for the kinematic viscosity coefficient *ν* and the superscript of the Gaussian random force *η*^{ν}(*t*).

Here, we comment on the reason why the NS equation for the stream function is used instead of the NS equation for the velocity field. Indeed, it is easy to check that Eq. (1) is obtained from the following NS equation:

where *ρ* and *p* are the density and pressure of the fluid, respectively, and *μ* = *ρν* is the viscosity. The fluid is assumed to be Newtonian. The NS equation in Eq. (4) has the form of the standard equation of motion per unit volume for a fluid of density *ρ*. The second term on the left-hand side (LHS), which is called the advection term, arises from the fact that the fluid particles are moving with velocity $V\u2192$. This term is very small compared with the other terms in the case of protoplasmic streaming; however, we include it for completeness. The first term on the right-hand side (RHS) represents the force from the pressure *p*; the negative sign appears by definition. The second term on the RHS, defined by the Laplace operator, represents the force from the viscosity of the fluid. The final term on the RHS is the random Brownian force *η*(*t*), determined by Gaussian random numbers, on which detailed information will be given below.

Equation (4) can be conveniently modified by multiplying both sides of the equation by *ρ*^{−1}, and by additionally incorporating the condition $\u2207\u22c5V\u2192\u2009=\u20090$, we obtain

By multiplying this standard NS equation from the left by the rotation ∇×, we obtain the NS equation in Eq. (1). The NS equation in Eq. (1) is used instead of Eq. (5) because the condition $\u2207\u22c5V\u2192\u2009=\u20090$ is exactly satisfied in Eq. (1); therefore, Eq. (1) is easier to solve numerically than the original NS equation for $V\u2192$ given in Eq. (5) for the case of protoplasmic streaming.^{30} This is the reason why the NS equation for the stream function, i.e., Eq. (1), is used instead of the NS equation for the velocity field, i.e., Eq. (5).

To obtain the discrete form of the NS equation in Eq. (1) on a 2D regular square lattice [Fig. 4(a)], we introduce the quantity

where *η*_{i,j}(*t*) denotes the random force on the fluid particle at lattice site (*i*, *j*) at time *t* [Fig. 4(b)]. Note that since the representation of the lattice site is changed from *i* to (*i*, *j*), $\eta \u2192i\u2192\eta \u2192i,j$ accordingly. This $H\u2192i,j(t)$ is still considered a stochastic variable. From the expression in Eq. (6) and the relation in Eq. (3), it is easy to obtain

which are typical characteristics of stochastic variables. The first relation comes from the fact that $\eta \u2192ij$ corresponds to Gaussian white noise with a mean value of zero, and the time integral and the expectation value operation ⟨⋯⟩ are assumed to be commutative. If we rewrite *H*_{i,j}(*t*) in Eq. (6) as

and substitute this *H*_{i,j}(*t*) into the second expression in Eq. (7), we obtain a finite value

for the random Brownian force. This expression is used in the discrete version of Eq. (1), which is given by

where $g\u2192i,j(t)\u2009=\u2009(gi,jx(t),gi,jy(t),0)$ and the components of $g\u2192i,j$ are given by Gaussian random numbers with mean 0 and variance 1. This $g\u2192i,j(t)$ is related to $\eta \u2192i,j$ as follows:

On the RHS of the first expression in Eq. (10), the spatial discretization of the second term with respect to the lattice spacing *Δx* is given by

This term makes almost no contribution to the flow because the velocity is low (no higher than ∼100 *μ*m/s) in the case of protoplasmic streaming. Thus, the results are expected to be independent of this term, although we include it in the equation for our simulations. The discrete form of the Laplace operator *Δ* acting on *ω*_{ij} is given by

and *Δψ*_{i,j} in the second expression in Eq. (10) has almost the same discrete form. The discrete form of the final term is given by

Notably, $2D\Delta t$ in Eq. (10) effectively corresponds to the deviation of the random Brownian force. As determined through dimensional analysis, $(2D\Delta t)2\Delta t\u2009=\u20092D(\Delta t)2$ is the diffusion constant *D*_{dif} related to the temperature *T* by means of the Einstein–Stokes–Sutherland formula *D*_{dif} = *k*_{B}*T*/6*πμa*, which is identified with 2*D*(*Δt*)^{2}. Here, we introduce the notion of the macroscopic relaxation time *τ*_{e}, which is the time required for the fluid to equilibrate from the resting state to a stationary state compatible with the boundary condition given by the velocity *V*_{B} [Fig. 3(c)], and this *τ*_{e} is independent of whether the initial state is the resting state or a random state.^{36–39} According to this definition, *τ*_{e} is proportional to the area (or volume, more generally), in sharp contrast to the standard relaxation time, which is the mean time required for a molecule to return to its original position from a disturbed position. We replace *Δt* with *τ*_{e} because *Δt* is a numerically introduced quantity; thus, we have

where *μ*(=*ρν*) is the viscosity, *a* is the size of a fluid particle or a group of particles in the fluid, and *k*_{B} is the Boltzmann constant. Note that *a* is larger than the size of a molecule such as water because it is obtained by assuming *τ*_{e}, which is not a microscopic quantity. The actual value of *D* assumed in the simulation and its relation to the *D*_{dif} value reported in Ref. 16 will be discussed in Sec. III.

The second equation of Eq. (10), which is Poisson’s equation, is numerically solved by the convergent configuration of the iterations such that

where the superscript *ℓ* is an integer denoting fictitious time or the number of iterations and the constant *A* is the acceleration coefficient fixed to *A* = 1. This technique is generally called the successive-over-relaxation (SOR) technique and is equivalent to the Gauss–Seidel method when *A* = 1. The convergence criteria will be discussed in Sec. III.

### B. Boundary conditions

The boundary conditions for the variables *ω* and *ψ* at the boundaries Γ_{1} and Γ_{3} [Fig. 5(a)] are given by

where the velocity $V\u2192$ at the boundary is given by

Note that the third component of $V\u2192$ is henceforth assumed to be zero and is thus neglected in all expressions for simplicity.

The reason why the stream function *ψ* can be fixed to *ψ* = 0 on Γ_{1} and Γ_{3} in Eq. (17) is as follows: The function *ψ* is not uniquely fixed in the domain because the velocity is given by the first-order differentials in Eq. (2). For this reason, *ψ* and *ψ* + *f*_{0} are exactly equivalent for any constant *f*_{0} in the sense that both *ψ* and *ψ* + *f*_{0} correspond to the same velocity configuration. Therefore, if $\psi 1,nX$ is nonzero such that $\psi 1,nX\u2009=\u2009c0$ at (1, *n*_{X}) ∈ Γ_{1}, then *ψ* can be replaced by *ψ* + *f*_{0} with *f*_{0} = −*c*_{0}, and hence, we have $\psi 1,nX\u2009=\u20090$ [Fig. 5(a)]. It is also easy to check that $\psi i,nX\u2009=\u20090\u2009(i>1)$ because of the boundary condition *V*_{y} = 0 on Γ_{1} in Eq. (18). On Γ_{3}, we also have *ψ*_{i,1} = 0 because of the symmetry argument under a rotation by *π* around the *z*-axis perpendicular to the domain. Evidently, no gravitational force is considered, and no asymmetry is expected in the random Brownian force; therefore, this rotational symmetry of *ψ* is naturally expected.

Note that the expressions for *ω* in Eq. (17) are well known and that the expression for $\omega i,nX$ on the boundary Γ_{1} is obtained by means of Taylor expansion and the second expression in Eq. (1) as follows:

where *Δy* = *Δx* is assumed and *V*_{x} = −*∂ψ*/*∂y* in Eq. (2) is used. The expression for *ω*_{i,1} on Γ_{3} in Eq. (17) can be obtained in the same manner.

On the boundaries Γ_{2} and Γ_{4}, periodic boundary conditions along the horizontal or *i* direction are assumed such that

These conditions imply that the lattice sites (1, *j*) on Γ_{2} and (n_{X}, *j*) on Γ_{4} are adjacent to each other [Fig. 5(a)].

### C. Physical and simulation units

In the actual process of protoplasmic streaming, length and time are measured in units of m and s, respectively, while the corresponding values in the simulations are *α* m and *β* s, where *α* and *β* are positive numbers. We use *α*m and *β*s to denote the corresponding so-called simulation units. The transformation rules from m and s to *α*m and *β*s are given as follows:

The numerical results should be independent of the values of *α* and *β*, which will be discussed in greater detail later. We consider the physical parameters

The kinematic viscosity coefficient *ν* is explicitly included in Eq. (10). The velocity *V*_{B} is necessary in the boundary conditions illustrated in Figs. 3(c) and 5(b). The third parameter, *d*, is the diameter illustrated in Fig. 5(b). The fourth parameter, *τ*, is set equal to the macroscopic relaxation time *τ*_{e} introduced in Sec. II A and numerically corresponds to the equilibration time of the fluid [see Fig. 4(b)]. This length of time should be no less than *τ*_{e}; otherwise, no equilibrium will be achieved in the numerical simulations. *D* in Eq. (15) must be expressed in units of m^{2}/s^{3} because each term in Eq. (10) must have units of 1/s, which is clear from the fact that *ω* is expressed in units of 1/s.

The final two parameters in Eq. (22), *Δx* and *Δt*, are necessary only in the simulations; however, these parameters do have indirect counterparts in the actual experimental phenomena. Specifically, these parameters have the following relations to the physical parameters *d* and *τ*:

where *n*_{X} is the total number of lattice points on one edge of the lattice and *n*_{T} is the total number of iterations. The parameter *τ* is simply related to the convergence of the time evolution corresponding to a set of random Brownian forces {*η*_{ij}(*t*)}, and it physically corresponds to the macroscopic relaxation time, which represents the typical time scale of the phenomenon, as mentioned above. This convergence is controlled by the small parameter *ε*, which will be discussed later in Sec. III. Therefore, exact information on *τ* or the macroscopic relaxation time *τ*_{e} is unnecessary, at least in simulations. Indeed, as long as *ε* is sufficiently small, then a configuration that is randomized by Brownian forces can correctly converge or reach equilibrium. Moreover, if *Δt* (*n*_{T}) is excessively large (small), then the simulation will not converge; hence, *Δt*_{0} should be fixed to *Δt*_{0} ≤ *Δt*_{cr}. This *Δt*_{cr} is considered the maximal time step satisfying the convergence criterion. Thus, the temporal discretization or time evolution is subtle compared to the spatial discretization. Both of the discretizations will be discussed in greater detail in Sec. III.

Here, we introduce the symbol *E* to represent the experimental parameters and the symbol *S* to represent the simulation parameters such that

In principle, the parameter *D* should be included in *E*; however, it is unknown for the flows under study. Therefore, *D* is instead included in *S* as an input for the simulations.

### D. Invariance under unit transformations

To discuss the invariance of the simulation results obtained using the discrete NS equation in Eq. (10) under unit transformations, we use the notion of scale transformations for the parameters m, s, *n*_{X}, and *n*_{T} introduced in Sec. II C. These scale transformations are defined in terms of positive numbers *α*, *β*, *γ*, and *δ*(>0) such that

The first two transformations, (m, s) → (*α*m, *β*s), represent changes to the units of length and time. The third transformation, *n*_{X} → *γn*_{X}, represents a change in lattice size, which corresponds to a change in the lattice spacing *Δx* → *γ*^{−1}*Δx* in Eq. (23). The final transformation, *n*_{T} → *δn*_{T}, represents a change in the time step *Δt* → *δ*^{−1}*Δt* by means of Eq. (23). *γ* and *δ* for *n*_{X} and *n*_{T} are necessary in addition to *α* and *β* for m and s because the simulation results should be independent of the lattice spacing *Δx* and the time step *Δt*. For a change in *Δx*, for example, the scaling of *n*_{X} by *α* can be used because of the relation *Δx* (m) = *d*/*n*_{X} (m) = *α*^{−1}*d*/*n*_{X} (*α*m) = *d*/(*αn*_{X}) (*α*m). Indeed, this relation implies that a unit transformation by *α* can be understood as a change in the lattice size *n*_{X}. However, in this case, it is impossible to observe the dependence of the results on the lattice size *n*_{X} without affecting other parameters that also depend on the length unit, which is why *γ* is necessary in addition to *α*. Moreover, as will be shown below in further detail, the rescaling of *Δx* always affects the force strength *D*, which implies that we need to modify *D* as well as *Δx* to observe the dependence of the results on *Δx*. This behavior arises from the fact that the order of *Δx* in the Brownian force term is different from that in the other terms; this difference is simply due to the difference in the order of spatial differentials such as *Δ* and ∇. The dependence of the results on *Δt* suffers from the same problem, as will also be studied in detail later. In this case, the difference in the order of *Δt* originates from the fact that the Brownian force is represented by a stochastic variable, as described in Sec. II A.

First, we rewrite Eq. (10) by explicitly including the lattice spacing *Δx* such that

where $(\Delta x)\u22122\psi \cdots \u2009\omega \cdots \u2009$ on the RHS represents the spatial discretization in Eq. (12) and $(\Delta x)\u22122\omega \cdots \u2009$ and $(\Delta x)\u22121gz\cdots \u2009$ represent the discretizations of the Laplacian *Δω* and the rotation $\u2207\xd7g\u2192z$, respectively. The symbol (1/s) represents the overall unit of the terms, which, in Eq. (26), are written with the physical units (m) and (s).

Under the scale transformations in Eq. (25), *Δx*(= *d*/*n*_{X}) (m) and *Δt*(= *τ*/*n*_{T}) (s) are replaced by *α*^{−1}*γ*^{−1}*Δx* (*α*m) and *β*^{−1}*δ*^{−1}*Δt* (*β*s), and we also have *ν* (m^{2}/s) = *α*^{−2}*βν* [(*α*m)^{2}/*β*s], and *D* (m^{2}/s^{3}) = *α*^{−2}*β*^{3}*D* [(*α*m)^{2}/(*β*s)^{3}]. The units of *ψ* and *ω* are m^{2}/s and 1/s, respectively; therefore, we have *ψ* (m^{2}/s) = *α*^{−2}*βψ* [(*α*m)^{2}/*β*s] and *ω* (1/s) = *βω* (1/*β*s). From these expressions and Eq. (26), we obtain

where the common factor *β* is eliminated from both sides. In the case of *γ* = 1 and *δ* = 1, nothing is changed except that the units are changed from m and s to *α*m and *β*s. The problem is the case of *γ* ≠ 1 or *δ* ≠ 1, where the factor $\gamma 2\delta \u22121$ in the final term is different from the factor *γ*^{2}*δ*^{−1} in the second and third terms. However, if *D* transforms as *D* → *γ*^{2}*δ*^{−1}*D* under the scale transformation (*n*_{X}, *n*_{T}) → (*γn*_{X}, *δn*_{T}), then we obtain a common factor of *γ*^{2}*δ*^{−1} in all three of these terms on the RHS of Eq. (27). In this case, we have

therefore, the convergent numerical solution is expected to remain unchanged. Indeed, in such a stationary or equilibrium configuration, the term *ω*_{i,j}(*t* + *Δt*) on the LHS is expected to be identical to the first term *ω*_{i,j}(*t*) on the RHS; hence, the common factor *γ*^{2}*δ*^{−1} in the remaining terms can be dropped.

The velocity *V* for the boundary conditions and the diameter *d* are included in the parameters *E* or *S* in Eq. (24), and their scaling properties under unit transformation are given by *V* (m/s) = *α*^{−1}*βV* (*α*m/*β*s) and *d* (m) = *α*^{−1}*d* (*α*m). Thus, under the scale transformations in Eq. (25), the RHS of Eq. (10) remains unchanged in the equilibrium configuration if the parameters *S* = (*ν*, *V*, *D*, *Δx*, *Δt*) scale as follows:

Next, we introduce the notion of equivalence in the simulation data. The simulation data obtained by solving Eq. (10) are denoted by (*ω*, *ψ*), while the experimental velocity data are denoted by Exp(*E*) because the experimental data Exp are characterized by the parameters *E* in Eq. (24). We define the term *equivalent* below.Two sets of simulation data (*ω*_{1}, *ψ*_{1}) and (*ω*_{2}, *ψ*_{2}) are *equivalent* if the following conditions are satisfied:

The histogram of the normalized velocity

*V*_{x}distribution andthe dependence of the normalized

*V*_{x}on*y*

for (*ω*_{1}, *ψ*_{1}) are identical to those for (*ω*_{2}, *ψ*_{2}).

Thus, we have proven the following statement:

Here, we comment on the dependence of the simulation results on *Δx*, as mentioned above. The parameter *γ* for rescaling *Δx* was introduced in Eq. (25) to elucidate this dependence. However, we find from statement (A) that *Δx* alone cannot be changed without affecting the results. The dependence of the results on *Δt* shows the same behavior as the dependence on *Δx*. This nonstandard situation with regard to *Δt* arises from the discrete form of the random Brownian force in Eq. (10) and is typical of such a discrete Langevin equation.^{25–32}

The experimental data Exp(*E*) can also be grouped into equivalent classes in the same way. Specifically, two different sets of experimental data Exp(*E*_{1}) and Exp(*E*_{2}) are considered equivalent if conditions (i) and (ii) are satisfied.

A similar notion of equivalence can be introduced for the parameters *S* in Eq. (24). Two sets of parameters *S*_{1} = (*ν*_{1}, *V*_{1}, *D*_{1}, *Δx*_{1}, *Δt*_{1}) and *S*_{2} = (*ν*_{2}, *V*_{2}, *D*_{2}, *Δx*_{2}, *Δt*_{2}) are called *equivalent* if there exists a set of positive numbers *α*, *β*, *γ*, and *δ* (all >0) such that

This equivalence is denoted by *S*_{1} ≡ *S*_{2}. It is easy to confirm that *S*_{2} ≡ *S*_{1} if *S*_{1} ≡ *S*_{2} because *α* and *β* can be inverted to *α*^{−1} and *β*^{−1}. *S*_{1} = *S*_{2} if and only if *α* = *β* = *γ* = *δ* = 1. A set of parameters equivalent to *S* = (*ν*, *V*, *D*, *Δx*, *Δt*) is written as $S\xaf=(\nu \xaf,V\xaf,D\xaf,\Delta x\xaf,\Delta t\xaf)$.

Notably, two different parameter sets *S*_{i} (*i* = 1, 2) produce the same solution (*ω*, *ψ*) for the discrete NS equation in Eq. (10) if these *S*_{i} (*i* = 1, 2) are equivalent; in other words, two solutions (*ω*_{1}, *ψ*_{1}) and (*ω*_{2}, *ψ*_{2}) are equivalent if the corresponding parameter sets *S*_{i} (*i* = 1, 2) are equivalent. In this sense, the solution to Eq. (10) depends only on the equivalent class of parameters $S\xaf$. If one of the parameters is transformed as *ν* → *ν*_{1}(≠*ν*) in $S\xaf$, then the parameters $(\nu \xaf1,V\xaf,d\xaf,D\xaf,\Delta x\xaf,\Delta t\xaf)$ are not equivalent to the original $S\xaf$.

Finally, in this subsection, we introduce the notation

for

which means that the experimentally observed data of the velocity distribution are equivalent to the simulation data in the sense defined above. The meaning of the expression Exp(*E*) ≃ *S*_{0} is that “the experimental data Exp(*E*) are successfully simulated with the parameters *S*_{0}.”

### E. Unique solution to the Navier–Stokes equation

The problem that we would like to clarify is how many parameters are sufficient to simulate the real experimental data Exp(*E*). We must consider this problem because the solution to Eq. (10) depends on many parameters *S* = (*ν*, *V*, *D*, *Δx*, *Δt*), even though only their equivalent classes are meaningful. One possible answer is that only the parameter *D* must be varied to enable the simulation of arbitrary Exp(*E*_{e}) data, while the remaining parameters can be fixed to the parameters *S*_{0} used to simulate certain existing experimental data Exp(*E*_{e,0}), which are not always identical to Exp(*E*_{e}). This process is described in more detail in the following statement:

- (B)
Let

*E*_{e,0}= (*ν*_{e,0},*V*_{e,0},*d*_{e,0}) be a set of parameters that characterize Exp(*E*_{e,0}), and let*S*_{0}be a set of parameters given by*S*_{0}= (*ν*_{0},*V*_{0},*D*_{0},*Δx*_{0},*Δt*_{0}). In this situation, if Exp(*E*_{e,0}) ≃*S*_{0}, then for any experimental data Exp(*E*_{e}), with*E*_{e}= (*ν*_{e},*V*_{e},*d*_{e}) and*D*_{e}, and for any given set of (*n*_{X},*n*_{T}), there exists a unique*D*_{sim}such that Exp(*E*_{e}) ≃ (*ν*_{0},*V*_{0},*D*_{sim},*Δx*_{0},*Δt*_{0}).

Statement (B) indicates that only one parameter, *D*_{sim}, must be varied to simulate arbitrary experimental data Exp(*E*_{e}). Proof of statement (B) and further details are provided in the Appendix.

## III. SIMULATION RESULTS

### A. Computational procedure and Langevin simulation technique

The simulation program is written in Fortran, and the Gaussian random numbers are generated by a Box–Muller transformation of uniform random numbers.^{43} The flow field is exactly the same as that for Couette flow if the Brownian force $\eta \u2192$ is zero, and whether the results obtained under the condition $\eta \u2192\u2009=\u20090\u2192$ are compatible with the expected solution will be determined below. In the case of $\eta \u2192\u2009\u2260\u20090\u2192$, the validity of the technique depends on whether the discrete Langevin dynamics with the term in Eq. (11) are meaningful. This problem has been shown to be meaningful in particle physics in Refs. 25–32.

The convergent configuration of the variable {*ω*} for the first equation in Eq. (10) is obtained using the small number,

for a time step of the NS equation such that

for *ω*, and the same small number *ε* is also assumed for the variable *ψ* such that

from which the convergent configuration {*ψ*} is obtained. *N*(=∑_{ij}1 = *n*_{X} × *n*_{X}) in Eqs. (33) and (34) represents the total number of vertices or the size of the lattice, and the subscript *ij* denotes a particular lattice site. The convergence of the SOR technique in Eq. (16) for Poisson’s equation, which is the second equation in Eq. (10), is given by the same number *ε* such that

at each time step *Δt*.

Note that *ε* is related to *τ* in Eq. (22). Indeed, *ε* determines the total number of iterations *n*_{T}, which satisfies *n*_{T} = *τ*/*Δt* and, hence, depends only on *τ* for a fixed *Δt*. From this relation, it is clear that *ε* depends on *τ*. If *ε* is excessively large, the iterative solution processes for the NS equation and the Poisson equation will not converge. By contrast, if *ε* is excessively small, then *n*_{T} will be very large, resulting in time-consuming simulations. The number *ε* in Eq. (32) is considered sufficiently small because if it is increased by a factor of 10, i.e., replaced with 10*ε*, and the results will remain unchanged. The dependence on the time step *Δt* will be checked separately inSec. III D.

The mean value of a physical quantity *Q* is calculated as

where *Q*_{k} is the *k*th sample corresponding to the *k*th convergent configuration {*ω*,*ψ*}_{k} and *n* denotes the total number of samples. The configuration {*ω*,*ψ*}_{k}, corresponding to a given set of Gaussian random forces {*g*}_{k}, is obtained by iterating the time step *Δt* and solving the Poisson equation in Eq. (10) once for every time increment of *Δt*. Thus, by repeating these two steps, we obtain the sample *Q*_{k} in Eq. (36) from the convergent configuration {*ω*,*ψ*}_{k}. The total number *n* of samples lies in the range of 1 × 10^{4} ≤ *n* ≤ 2 × 10^{4} for all simulations. Note that the formula given in Eq. (36) for the mean value is exactly the same as that for the Monte Carlo simulation technique for statistical mechanical models,^{41,42} and these two techniques are known to be equivalent.^{25–31} This is the reason why we call Eq. (1) [or Eq. (10)] the stochastic NS equation.

We note that the mean value of the velocity $V\u2192$, for example, is independent of the order of the calculations, which is as follows:

The mean value

*ψ*may first be calculated from the configurations {*ω*,*ψ*}_{k}(*k*= 1,*n*), and $V\u2192$ can then be calculated using this mean value configuration*ψ*.${V\u2192}k$ may first be calculated from the

*k*th configuration {*ω*,*ψ*}_{k}, and $V\u2192$ can then be calculated as the mean value of the ${V\u2192}k\u2009(k\u2009=\u20091,n)$.

This is simply because of the commutativity of the calculation of $V\u2192$ using Eq. (2) and the mean value calculation using Eq. (36).

Regarding the relation to the lattice Boltzmann method (LBM), the stochastic NS equation is obtained simply by including the Brownian force term in the NS equation; hence, this definition is slightly different from that of the LBM, which is a technique for describing fluid flow from the viewpoint of particle mechanics.^{40} Nevertheless, the stochastic NS equation technique can be regarded as a special case of the LBM because all terms, including the random Brownian force term, are understood to represent forces acting on the fluid particles. However, we do not delve into the details of this problem because our interest is simply focused on reproducing the experimentally observed velocity distribution.

Here, we comment on the dependence of the results on the initial configuration of the variables {*ω*, *ψ*}. On the boundaries Γ_{1} and Γ_{3}, the variables are fixed to certain assumed values in accordance with the boundary conditions. Inside the flow region, two possible initial configurations can be assumed for {*ω*,*ψ*}_{k}: *ω* = *ψ* = 0 and the convergent configuration {*ω*,*ψ*}_{k−1}, where {*ω*,*ψ*}_{0} is *ω* = *ψ* = 0. The results, including *n*_{T}, which is the total number of iterations in the sense of the mean values, are independent of these initial configurations, implying that the macroscopic relaxation time *τ*_{e} is also independent of the initial conditions.

### B. Normalized velocity distribution

First, before presenting the velocity distribution, we show snapshots of the flow field in Figs. 6(a)–6(d). Figure 6(a) shows the snapshot obtained for *D*_{sim} = 0, and Fig. 6(b) shows the mean *ψ* values of 1000 convergent configurations for *D*_{sim} = 400 and the $V\u2192$ results calculated using this mean *ψ*. Figures 6(c) and 6(d) show the snapshots of two different convergent configurations for the same *D*_{sim} = 400. The parameters other than *D*_{sim} are the same as those used for the velocity distribution, which will be presented below. In the graphics, the flow velocity is represented by small cones, and the stream function is normalized to −1 ≤ *ψ* ≤ 0 and represented by a gradient between two different colors, with *ψ* reaching its maximum value of *ψ* = 0 at the boundaries Γ_{1} and Γ_{3}. If the random Brownian force is neglected (⇔*D*_{sim} = 0), the flow field is uniquely determined by Eq. (A13). The snapshot in Fig. 6(a) is compatible with this expected solution. In contrast, for a nonzero *D*_{sim}, the flow field exhibits fluctuations even if it is convergent. The graphic in Fig. 6(b) appears to be almost the same as that in (a); however, the velocity distributions are different, as will be shown below. In Figs. 6(c) and 6(d), $V\u2192$ and *ψ* are different from those in Fig. 6(b) and fluctuate from one convergent configuration to another. We emphasize that these fluctuating configurations are ensemble configurations, which are thermally fluctuating in the statistical mechanical sense, and are not always identical to experimentally observed configurations in general. However, they are understood to be some of the possible configurations, and it is interesting that such vortex configurations are included among the ensemble configurations. In some specific cases, a vortex configuration can dominate; however, we do not address this problem in detail. The experimentally observable quantity that we numerically study in this paper is the velocity distribution, which will be presented below.

Now, the main results are presented. In Figs. 7(a) and 7(b), we plot the distribution (or normalized histogram) *h*(*V*_{x}) of the absolute velocity |*V*_{x}| along the *x*-direction (the longitudinal direction) and the distribution *h*(*V*) of the magnitude *V* of the velocity vector $V\u2192\u2009=\u2009(Vx,Vy)$. The discrete expressions are

which are the standard discrete forms corresponding to the first-order differentials in Eq. (2).

The lattice size is fixed to *N* = 10 000, with *N* = *n*_{X} × *n*_{X} and *n*_{X} = 100. According to statement (B), the simulation results are expected to depend only on the strength of the random Brownian force, *D*_{sim}; hence, *D*_{sim} is varied from *D*_{sim} = 0 to *D*_{sim} = 4000 in the simulation units. The lattice size dependence will be presented in Sec. III C.

The distribution *h*(*V*_{x}) is obtained by constructing a histogram in which the velocity range $0\u2264|Vx|\u2264vxmax$ is divided into 100 subranges, and the velocity *V*_{x} in Eq. (37) at lattice point (*i*, *j*) is counted in the histogram for every convergent configuration. The maximum velocity $vxmax$ is fixed to $vxmax\u2009=\u20092VB$, where *V*_{B} denotes the velocity at the boundaries Γ_{1} and Γ_{3}. The reason for this choice of $vxmax$ is that $vx,(i)max$ in the *i*th convergent configuration is expected to vary with *i* instead of remaining constant. The factor of 2 in 2*V*_{B} for $vxmax$ is sufficiently large because the fluctuations in $vx,(i)max$ are relatively small. Such an assumption regarding the maximum velocity for velocity normalization is unnecessary if the mean value of *ψ* is calculated first, in accordance with procedure (i) described in Sec. III A, because the velocity distribution is obtained from a single configuration corresponding to the mean *ψ* in that case. The histogram *h*(*V*_{x}) is also independent of procedures (i) and (ii) for the calculation of *V*_{x}. The height of *h*(*V*_{x}) is normalized such that the maximum height is equal to 1 for each *D*_{sim}, and the horizontal axis |*V*_{x}| is similarly normalized using the maximum |*V*_{x}| satisfying *h*(*V*_{x}) ≠ 0 for each *D*_{sim}. The histograms *h*(*V*) with respect to *V* in Fig. 7(b) are normalized in the same way.

From Fig. 7(a), we find that *h*(*V*_{x}) is flat for *D*_{sim} = 0; this flat *h*(*V*_{x}) is compatible with the expectation from Eq. (A13). For a sufficiently large *D*_{sim}, the nonflat *h*(*V*_{x}) has a peak only at *V*_{x} = 0. However, for intermediate *D*_{sim} values, another peak appears at *V*_{x} ≠ 0. If all the ensemble configurations for *D*_{sim} ≠ 0 were identical to that in Fig. 6(a) for *D*_{sim} = 0, then *h*(*V*_{x}) would be expected to be flat. Therefore, configurations including vortices, such as those shown in Figs. 6(c) and 6(d), are understood to be the reason for the peaks observed in *h*(*V*_{x}) in numerical studies, even though no such vortex characteristics are apparent in the mean value configuration in Fig. 6(b). We should note that the flow field shown in Fig. 6(b) is understood to be slightly different from the observable configurations in experimental measurements of *Nitella* cells due to the imposed simplifications. The distribution *h*(*V*) in Fig. 7(b) drops to zero, *h*(*V*) → 0 in the limit of *V* → 0, in contrast to *h*(*V*_{x}) in Fig. 7(a). This drop is reasonable because the fluid is always moving and there are no fluid particles with zero velocity, $V\u2192(t)\u2009=\u20090\u2192$, for all *t*. The drop in *h*(*V*) at *V* → 0 corresponds to the same drop that is visible in the experimentally reported data in Ref. 21. Indeed, with the light scattering technique, not only $V\u2192(=\u2009(Vx,0))$ but also $V\u2192(=\u2009(Vx,Vy))$, which has a small nonzero *V*_{y} component, can be detected. This drop of *h*(*V*) at *V* → 0 in Fig. 7(b), for sufficiently large *D*_{sim}, is also consistent with the drop observed in the Maxwell–Boltzmann distribution *V* exp(−*V*^{2}/(2*c*^{2})) expected in the ideal gas. In this expression, the constant *c*(= 0.22) is obtained from the peak position *V* of *h*(*V*) for *D*_{sim} = 4000 in Fig. 7(b). The simulation result *h*(*V*) for *D*_{sim} = 4000 is almost consistent with the dashed line as expected, even though *D*_{sim} is finite. However, the corresponding *h*(*V*_{x}) considerably deviates from the dashed line $exp(\u2212Vx2/(2c2))$ in Fig. 7(a). The reason for this deviation is that the equipartition of energy, $\u27e8Vx2\u27e9\u2009=\u2009\u27e8Vy2\u27e9$, is obviously violated due to the effect of boundary velocity, at least for finite *D*_{sim}. Nevertheless, the ideal gas behavior observed in *h*(*V*) for sufficiently large *D*_{sim} is also reasonable because *D*_{sim} is proportional to the temperature *T* as a result of the Einstein–Stokes–Sutherland formula; hence, for sufficiently large *T*, the random Brownian force is expected to be very large compared to the other interactions within the fluid. The most important point to note is that the simulation results at finite *D*_{sim}(≠0) are located between $exp(\u2212Vx2/(2c2))$ (dashed line) expected from the ideal gas and *h*(*V*_{x}) = 1 (○) expected from the exact solution of the Couette flow.

We will now discuss the parameters used in the simulations in detail. The viscosity is considered to be almost 100 times greater than that of water. In Ref. 2, the viscosity is reported to be 0.5 ≤ *μ* ≤ 1.5 (dyn s/cm^{2}). Therefore, if the density *ρ* is assumed to be the same as that of water, then the kinematic viscosity *ν*_{e,0} ranges from *ν*_{e,0} = 0.5 × 10^{−4} m^{2}/s to *ν*_{e,0} = 1.5 × 10^{−4} m^{2}/s. Thus, we assume that *ν*_{e,0} = 1 × 10^{−4} m^{2}/s. In the experiment conducted by Kamiya,^{5} a cell with a diameter of 0.46 mm was used for the velocity measurements, and 50 *μ*m/s was observed at the boundary. For this reason, a velocity of *V*_{e,0}(= 50 *μ*m/s) at the boundary and a diameter of *d*_{e,0}(= 500 *μ*m) are assumed, as shown in Table II. The parameters *α*_{0} and *β*_{0} in Eq. (21) that are used in the simulations are also shown in Table II. These values, *α*_{0} = 1 × 10^{−6} and *β*_{0} = 1 × 10^{−1}, imply that the simulation units for length and time are 1 *μ*m and 0.1 s, respectively. The first three parameters in Table II are collectively denoted by *E*_{e,0}. The corresponding *α*_{0} and *β*_{0} values are obtained via Eq. (A1) using the parameters *S*_{0} shown in Table III.

ν_{e,0} (m^{2}/s)
. | V_{e,0} (μm/s)
. | d_{e,0} (μm)
. | α_{0}
. | β_{0}
. |
---|---|---|---|---|

1 × 10^{−4} | 50 | 500 | 1 × 10^{−6} | 1 × 10^{−1} |

ν_{e,0} (m^{2}/s)
. | V_{e,0} (μm/s)
. | d_{e,0} (μm)
. | α_{0}
. | β_{0}
. |
---|---|---|---|---|

1 × 10^{−4} | 50 | 500 | 1 × 10^{−6} | 1 × 10^{−1} |

$\nu 0\u2009[(\alpha 0m)2\beta 0s]$ . | $V0\u2009[\alpha 0m\beta 0s]$ . | Δx_{0} (α_{0}m)
. | Δt_{0} (β_{0}s)
. | n_{X0}
. |
---|---|---|---|---|

1 × 10^{7} | 5 | 5 | 8 × 10^{−8} | 100 |

$\nu 0\u2009[(\alpha 0m)2\beta 0s]$ . | $V0\u2009[\alpha 0m\beta 0s]$ . | Δx_{0} (α_{0}m)
. | Δt_{0} (β_{0}s)
. | n_{X0}
. |
---|---|---|---|---|

1 × 10^{7} | 5 | 5 | 8 × 10^{−8} | 100 |

The lattice spacing *Δx*_{0} (*α*_{0}m) is calculated using Eq. (A11). For *Δx*_{0} = 5 *α*_{0}m, we have $\Delta xe,0\u2009=\u2009\Delta x0\alpha 0\u22121\u2009=\u20095\xd710\u22126\u2009m$.

There is no need to mention that the parameters *S*_{0} in Table III are also obtained from the parameters *E*_{e,0} in Table II and *α*_{0} and *β*_{0}. Indeed, it is easy to check that the kinematic viscosity $\nu 0\u2009=\u20091\u2009\xd7\u2009107\u2009\alpha 02m2/\beta 0s$ is obtained from the relation $\nu e,0\u2009(m2/s)=1\u2009\xd7\u200910\u22124\beta 0/\alpha 02\u2009(\alpha 02m2/\beta 0s)=\nu 0\u2009(\alpha 02m2/\beta 0s)$. The velocity *V*_{0} (*α*_{0}m/*β*_{0}s) and *d*_{0} (*α*_{0}m) are also obtained through the relations $Ve\alpha 0\u22121\beta 0\u2009(\alpha 0m/\beta 0s)\u2009=\u2009V0$ and $de\alpha 0\u22121\u2009=\u2009d0\u2009(\alpha 0m)$, which are not explicitly used in the simulations.

Each value of *D*_{sim} shown in Table IV is fixed as an input to the simulations, and the corresponding physical strength *D*_{e,0} can be obtained as $De,0\u2009=\u2009Dsim\alpha 02\beta 0\u22123$. The parameter *τ*_{0} in Table IV is estimated using the relation *τ*_{0} = *n*_{T}*Δt*_{0}, where *τ*_{0} is expected to satisfy *τ*_{0} ≤ *n*_{T}*Δt*_{cr}, as discussed in Sec. II C. Using these *D*_{e,0} and *τ*_{0} values, we can estimate the diameter *a* (m) of a fluid particle using Eq. (15), where the temperature is assumed to be *T* = 300 K and *k*_{B}*T* = 4.1 × 10^{−21} Nm. Although the time scale *τ* in Eq. (15) is expected to be smaller than *τ*_{0}, we use the results of *τ*_{0} = *n*_{T}*Δt*_{0} to calculate *a*. Nevertheless, the value of *a* calculated in this way should be larger than the size of a water molecule, i.e., 4 × 10^{−10} m, and indeed, we find that almost all the data in Table IV satisfy this condition.

$Dsim\u2009[(\alpha 0m)2(\beta 0s)3]$ . | D_{e,0} (m^{2}/s^{3})
. | n_{T}Δt_{0} (β_{0}s)
. | a (m)
. |
---|---|---|---|

200 | 2 × 10^{−7} | 8.7 × 10^{−3} | 3.6 × 10^{−9} |

400 | 4 × 10^{−7} | 7.0 × 10^{−3} | 2.8 × 10^{−9} |

800 | 8 × 10^{−7} | 7.0 × 10^{−3} | 2.8 × 10^{−9} |

1200 | 1.2 × 10^{−6} | 1.1 × 10^{−2} | 7.5 × 10^{−10} |

4000 | 4 × 10^{−6} | 2.6 × 10^{−2} | 4 × 10^{−10} |

$Dsim\u2009[(\alpha 0m)2(\beta 0s)3]$ . | D_{e,0} (m^{2}/s^{3})
. | n_{T}Δt_{0} (β_{0}s)
. | a (m)
. |
---|---|---|---|

200 | 2 × 10^{−7} | 8.7 × 10^{−3} | 3.6 × 10^{−9} |

400 | 4 × 10^{−7} | 7.0 × 10^{−3} | 2.8 × 10^{−9} |

800 | 8 × 10^{−7} | 7.0 × 10^{−3} | 2.8 × 10^{−9} |

1200 | 1.2 × 10^{−6} | 1.1 × 10^{−2} | 7.5 × 10^{−10} |

4000 | 4 × 10^{−6} | 2.6 × 10^{−2} | 4 × 10^{−10} |

The diffusion constant *D*_{dif} in Eq. (15) can be estimated by assuming that *n*_{T}*Δt*_{0} is *τ*_{e} as follows: $Ddif\u2009=\u20092De,0(nT\Delta t0\beta 0)2$. For *D*_{e,0} = 4 × 10^{7} m^{2}/s^{3} in Table IV, we have *D*_{dif} ≃ 4 × 10^{−11} m^{2}/s = 4 × 10^{−7} cm^{2}/s, which is almost 10 times smaller than the estimated value of *D*_{dif} = *k*_{B}*T*/6*πμa* ≃ 5.3 × 10^{−6} cm^{2}/s obtained with *μ* = 0.1 m^{2}/s and *a* = 4 × 10^{−10} m. Note, however, that this value of 5.3 × 10^{−6} cm^{2}/s is comparable to the value of 10^{−5} cm^{2}/s reported in Ref. 16. The deviation between the estimates $Ddif\u2009=\u20092D\tau e2$ and *D*_{dif} = *k*_{B}*T*/6*πμa* is expected to shrink in one of the following two possible cases: either *a* may be considered equal to the radius of a group of water molecules and, hence, should be larger than the radius of a single water molecule, or *τ*_{e} may be slightly larger than *n*_{T}*Δt*_{0}*β*_{0}. If either of these conditions is satisfied, then the two estimates are almost compatible, and the notion of *τ*_{e} as adopted in Eq. (15) is reasonable. Clearly, the first condition, at least, is quite reasonable.

Now, let us comment on the Schmidt number S_{c}, which is calculated as the ratio of the Reynolds number R_{e} and the Péclet number P_{e} such that S_{c} = P_{e}/R_{e}. The Reynolds number is evaluated to be R_{e} = *Vd*/*ν* = 5 × 10^{−5} because the velocity, diameter, and kinematic viscosity are assumed to be *V* = 50 × 10^{−6} m/s, *d* = 500 × 10^{−6} m, and *ν* = 1 × 10^{−4} m^{2}/s, respectively. The Péclet number is similarly evaluated to be P_{e} = *Vd*/*D*_{dif} ≃ 10 for *D*_{dif} = 5.3 × 10^{−10} m^{2}/s. If we assume that $Ddif\u2009=\u20094\u2009\xd7\u200910\u221211\u2009m2/s(=\u20092De,0(nT\Delta t0\beta 0)2)$, we obtain P_{e} ≃ 100, which is closer to the estimate of 10^{2}–10^{3} given in Ref. 16. Therefore, we obtain S_{c} ≃ 2 × 10^{5} for *D*_{dif} = 5.3 × 10^{−10} m^{2}/s and S_{c} ≃ 2 × 10^{6} for *D*_{dif} = 4 × 10^{−11} m^{2}/s. Note also that S_{c} can be obtained directly from its definition, S_{c} = *ν*/*D*_{dif}. Thus, the relatively large S_{c} implies that the viscosity force is larger than the diffusion force, which is activated by thermal fluctuations. Therefore, the effect of thermal fluctuations is relatively small compared not only to the advection term^{16,17} but also to the viscosity term. Nevertheless, the peaks in the velocity distribution essentially originate from this small thermal fluctuation effect.

### C. Lattice size dependence

In this subsection, we show that the simulation results are independent of the lattice size *n*_{X}. As mentioned in Sec. II D, we need to change not only *n*_{X} but also *D*_{sim} to observe this dependence. Figure 8(a) shows the histograms *h*(*V*_{x}) with respect to |*V*_{x}| obtained on lattices with sizes ranging from *n*_{X} = 100 to *n*_{X} = 300. In these simulations, the parameters *D*_{sim} and *Δx*_{0} are scaled to *γ*^{2}*D*_{sim} and *γ*^{−1}*Δx*_{0}, as indicated in Eq. (29). The results remain unchanged when the lattice size is modified from *n*_{X}(= 100) to *γn*_{X}, where 1 ≤ *γ* ≤ 3 [see Eq. (25)]. We start with *γ* = 1 for *D*_{sim} = *γ*^{2}200 and *Δx*_{0} = *γ*^{−1}5 (*α*m) with *Δt*_{0} = 4 × 10^{−8} *β*s. The parameters *D*_{sim} and *Δt*_{0} can also be replaced with *D*_{sim} = *γ*^{2}400 and *Δt*_{0} = 8 × 10^{−8} *β*s, respectively, which are identical to those plotted in Figs. 7(a) and 7(b). The reason *D*_{sim} = *γ*^{2}400 and *Δt*_{0} = 8 × 10^{−8} *β*s are replaced with *D*_{sim} = *γ*^{2}200 and *Δt*_{0} = 4 × 10^{−8} *β*s, respectively, is that if we start with *Δt*_{0} = 8 × 10^{−8} *β*s, the simulation does not converge for the case of *γ* = 3.

The height of the histogram *h*(*V*_{x}) is normalized such that the maximum height is equal to 1 for each *D*_{sim}, as in Fig. 7(a), while the horizontal axes |*V*_{x}| for all *D*_{sim} are normalized using a constant value equal to the maximum |*V*_{x}| for *D*_{sim} = 200 and a lattice size of *n*_{X} = 100 satisfying *h*(*V*_{x}) ≠ 0.

We find that *h*(*V*_{x}) with respect to |*V*_{x}| is almost independent of either *n*_{X} or *Δx*. This finding supports not only the correctness of statement (B) but also the *Δx* independence of the results. In fact, if *h*(*V*_{x}) with respect to |*V*_{x}| did depend on *Δx*, we would need to consider either statement (B) to be incorrect or the results to depend on *Δx*.

Additionally, the dependence of *V*_{x} on *y*, plotted in Fig. 8(b), is almost linear, and this linear behavior is the same as that of the trivial solution in Eq. (A13). We should note that this linear behavior is different from the previously experimentally observed behaviors.^{15,18} The reason for this deviation is considered to be the simplification of the model, as stated in the Introduction and in Sec. II A. The *y* dependence of *V*_{x} is found simply by averaging *V*_{x} from the convergent configurations using Eq. (37). Thus, the nontrivial behavior of two distinct peaks in *h*(*V*_{x}) is not always reflected in the dependence of *V*_{x} on *y*. The parameters used in the simulations are shown in Fig. 8(b).

### D. Discrete time-step dependence

In this subsection, the dependence of the results on *Δt* is checked, as mentioned in Sec. II D. This subsection contributes to the numerical Proof of statement (B). There are two possible origins, *τ* and *n*_{T}, for the change in *Δt*, as described in Eq. (23). However, in sharp contrast to the case of the diameter *d* for *Δx*, the relaxation time *τ* is not clear. Nevertheless, to check the scaling properties in Eq. (29), it is sufficient to observe the dependence of the results on *Δt* or *δ*. The normalization of the histograms *h*(*V*_{x}) with respect to |*V*_{x}| is defined in the same way as in Fig. 8(a).

In Fig. 9(a), the histograms *h*(*V*_{x}) of velocity vs |*V*_{x}| are plotted. The parameters *D*_{sim} and *Δt*_{0} are scaled to *δ*^{−1}*D*_{sim} and *δ*^{−1}*Δt*_{0}, respectively. The results are independent of *Δt*_{0}, which is varied in the range 2 × 10^{−8} ≤ *Δt*_{0} ≤ 32 × 10^{−8} in the simulation units (1 *β*_{0}s = 0.1 s). This range of *Δt*_{0} corresponds to a *δ* range of 1 ≤ *δ* ≤ 16. Thus, the results plotted in Fig. 9(a) confirm that the corresponding scaling properties, such as *δ*^{−1}*D*_{sim} and *δ*^{−1}*Δt*_{0}, are correct. This completes the numerical verification of statement (B).

Finally, in this subsection, we check whether *n*_{T}*Δt*_{0} depends on the *ε* used for the convergence criteria in Eqs. (33) and (34), which has previously been fixed to *ε* = 1 × 10^{−9} for all simulations. Here, the value is additionally set to *ε* = 1 × 10^{−7}, *ε* = 1 × 10^{−8}, and *ε* = 1 × 10^{−10} to assess the dependence of *n*_{T}*Δt*_{0} on *ε*. The results shown in Fig. 9(b) indicate that *n*_{T}*Δt*_{0} is roughly independent of *ε* for sufficiently small *Δt*_{0}.

### E. Dependence on physical parameters

In this subsection, we demonstrate how to use statement (B) to discuss the dependence of the normalized velocity distribution on the physical parameters *ν*_{e}, *V*_{e}, and *d*_{e}. From the perspective of statement (B), the purpose of the simulations in Sec. III B is to find *S*_{0} = (*ν*_{0}, *V*_{0}, *D*_{0}, *Δx*_{0}, *Δt*_{0}) with a suitable *D*_{0} and the set of parameters listed in Table III. Indeed, from the simulation results, we find that *D*_{0} = 400 in the second line of Table IV is suitable because the shape of *h*(*V*_{x}) in Fig. 7(a) is relatively close to the experimental data reported in Refs. 21–23. Since |*V*_{x}| in Fig. 7(a) is normalized, the position of the second peak can only be compared to those in Figs. 2(a) and 2(b). The peak position of the result of *D*_{sim} = 400 is approximately 0.6, while that of Fig. 2(b) is approximately 0.5; the deviation from the peak in Fig. 2(a) is clearly larger. However, the positions of the second peaks in Figs. 2(a) and 2(b) move to the right if the high-frequency part, where the intensity is close to zero, is removed. Although the relative intensity of the second peak in the numerical data is higher than those in Figs. 2(a) and 2(b), we consider that the existence of the second peak is clearly reproduced by this simplified 2D model. Thus, using the parameters *S*_{0}, we can perform simulations of Exp(*E*_{e}) characterized by *E*_{e}, which is different from *E*_{e,0}.

The parameters are shown in Table V, where *ν*_{e}, *V*_{e}, and *d*_{e} are the elements of the experimental data *E*_{(i)} (*i* = 1, 2, 3). *E*_{(1)}, *E*_{(2)}, and *E*_{(3)} are different from *E*_{e,0} in Table II only in terms of *ν*_{e}, *V*_{e}, and *d*_{e}, respectively (indicated by underlines). Note that *ν*_{e}/*ν*_{e,0} = 2, *V*_{e}/*V*_{e,0} = 2, and *d*_{e}/*d*_{e,0} = 2 imply that *ν*_{e} = 2*ν*_{e,0} = 2 × 10^{−4} m^{2}/s, *V*_{e} = 2*V*_{e,0} = 100 *μ*m/s, and *d*_{e} = 2*d*_{e,0} = 1 mm. All of these values are meaningful in the engineering viewpoint because these parameters in plant cells such as *Chara corallina* and *Nitella flexilis* are not always uniquely determined but distributed around the values *ν*_{e,0}, *V*_{e,0}, and *d*_{e,0} in Table II depending on their size.^{2}

E_{(i)}
. | $\nu e\nu e,0$ . | $VeVe,0$ . | $dede,0$ . | $\tau e\tau e,0$ . | $DeDe,0$ . | $\alpha \alpha 0$ . | $\beta \beta 0$ . | $\gamma \gamma 0$ . | $\delta \delta 0$ . | $DsimD0$ . |
---|---|---|---|---|---|---|---|---|---|---|

(1) | $2\u0332$ | 1 | 1 | $12$ | 2 | 2 | 2 | $12$ | $14$ | 4 |

(2) | 1 | $2\u0332$ | 1 | 1 | 1 | $12$ | $14$ | 2 | 4 | $116$ |

(3) | 1 | 1 | $2\u0332$ | 4 | $14$ | 1 | 1 | 2 | 4 | $14$ |

E_{(i)}
. | $\nu e\nu e,0$ . | $VeVe,0$ . | $dede,0$ . | $\tau e\tau e,0$ . | $DeDe,0$ . | $\alpha \alpha 0$ . | $\beta \beta 0$ . | $\gamma \gamma 0$ . | $\delta \delta 0$ . | $DsimD0$ . |
---|---|---|---|---|---|---|---|---|---|---|

(1) | $2\u0332$ | 1 | 1 | $12$ | 2 | 2 | 2 | $12$ | $14$ | 4 |

(2) | 1 | $2\u0332$ | 1 | 1 | 1 | $12$ | $14$ | 2 | 4 | $116$ |

(3) | 1 | 1 | $2\u0332$ | 4 | $14$ | 1 | 1 | 2 | 4 | $14$ |

We assume that the macroscopic relaxation time *τ*_{e} in Eq. (15) is proportional to the inverse kinematic viscosity $\nu e\u22121$^{39} and the area *A*_{e} as follows:

Using this relation, we obtain the ratio *τ*_{e}/*τ*_{e,0} in Table V. Consequently, from Eq. (15), we have

where *A*_{e} is replaced with $de2$, the viscosity *μ* is proportional to *ν*_{e}, and the temperature is assumed to be constant.

We will now explain how to obtain the values of *D*_{sim}/*D*_{0} via statement (B). Since Exp(*E*_{e,0}) is simulated with *S*_{0} = (*ν*_{0}, *V*_{0}, *D*_{0}, *Δx*_{0}, *Δt*_{0}), from statement (B), we have

using the parameters *α*_{0}, *β*_{0}, *γ*_{0}, and *δ*_{0} for the scale transformation between *S*_{0} and *S*_{e,0}. This is the assumption component, or the trivial case of statement (B), and is exactly the same as Eq. (A2). We also have

using the parameters *α*, *β*, *γ*, and *δ* for the scale transformation between *S*_{0} and *S*_{e}. Therefore, from Table V, we have

which imply

From Eqs. (A6) and (A7), we have *γ*/*γ*_{0} and *δ*/*δ*_{0} as listed in Table V for *E*_{(1)} and *E*_{(2)}. Thus, we obtain *D*_{sim}/*D*_{0} using these values of *α*/*α*_{0}, *β*/*β*_{0}, *γ*/*γ*_{0}, and *δ*/*δ*_{0} and by Eqs. (40) and (41), as follows:

Therefore, the experimental data corresponding to Exp(*E*_{(i)}) (*i* = 1, 2, 3) can be simulated with

and with the parameters in Table III. Thus, we expect the peak position for *E*_{(1)} (*E*_{(2)} or *E*_{(3)}) to move to the left (right) of the peak position for *E*_{e,0}. Moreover, the simulation results for Exp(*E*_{(i)}) (*i* = 1, 2, 3) are located between, or are only slightly different from, the curves in Figs. 7(a) and 7(b). This is why we regard the results in Sec. III B as our main results, which we emphasize in this subsection. Note that the results are true only if the assumption regarding *τ*_{e} in Eq. (38) is true. It should also be noted that the results in Eq. (45) are qualitatively reasonable. Indeed, the Reynolds number R_{e}(=*Vd*/*ν*) is decreased for *E*_{(1)} and increased for *E*_{(2)} and *E*_{(3)}, and consequently, the second peak position is expected to move in opposite directions depending on the variation of R_{e}.

Finally, we must note the implications of statement (B) and its supporting analyses in this subsection. Our main result that all qualitatively different simulation results can be obtained merely by varying the strength *D* of the random Brownian force, as stated in the Introduction, is limited in the sense that this is true only if the simulation results obtained with the parameter set *S*_{0} correspond to Exp(*E*_{e,0}). In fact, this assumption is not always exactly satisfied because the experimental data in Figs. 2(a) and 2(b) are not exactly the same as the simulation data for *D*_{sim} = 400 in Fig. 7(a). However, for the second peak position, as mentioned above, these experimental and simulation data are almost the same; therefore, statement (B) indicates that the second peaks of all the experimental normalized velocity distributions corresponding to different (*ν*_{e}, *V*_{e}, *D*_{e}), such as those in Table IV, are identical to or located between the peaks in Fig. 7(a).

## IV. SUMMARY AND CONCLUSION

In this paper, we study the flow fields associated with protoplasmic streaming in plant cells such as *Chara corallina* and *Nitella flexilis* by means of the stochastic or Langevin Navier–Stokes (NS) equation, which is a 2D equation for describing incompressible viscous flows with random Brownian forces. The study focuses on the experimentally observed distribution of the velocity along the flow direction, which exhibits two distinct peaks. To clearly illustrate the role of the random Brownian force assumed in the NS equation, the computational model is simplified such that the twist of the flows is neglected and 2D Couette flow is assumed.

From a dimensional analysis of the Langevin NS equation, we find that the normalized velocity distribution depends only on the strength *D* of the random Brownian force. This finding is numerically verified in detail in Sec. III B. Furthermore, in Sec. III E, we extract the reasonable finding that the position of the second peak moves to the right or left in accordance with the variation in the physical parameters, i.e., the kinematic viscosity, diameter, and boundary velocity. If the kinematic viscosity is decreased, for example, then the peak position is expected to move to the right or a higher-velocity region. This phenomenon is consistent with the intuitive understanding of the Reynolds number in the sense that a decrease in the kinematic viscosity is equivalent to an increase in velocity.

The results in this paper imply that the spiral flow and the 3D nature of real protoplasmic streaming are not essential for the emergence of the two peaks in the velocity distribution, although the shapes of the peaks are expected to be influenced by these experimentally observed characteristics. Rather, the random Brownian forces, represented by Gaussian random numbers, are confirmed to be the origin of the peaks.

As noted in Sec. III B, the dependence of *V*_{x} on *y* is almost linear and is slightly different from the experimentally reported results. One reason for this deviation is that compared with real flows, the simulation model used in this paper is simplified in many respects. In particular, this simulation model is a 2D model, and the twist of the flows and the interactions between the fluid and biological materials are neglected, as mentioned above. These neglected components should be included in the model for further fluid mechanics studies on protoplasmic streaming.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

H.K. acknowledges Dr. Kazuhiko Mitsuhashi for reminding him of this interesting field. This work was supported, in part, by a Collaborative Research Project of the Institute of Fluid Science (IFS), Tohoku University, and a Collaborative Research Project of the National Institute of Technology (KOSEN), Sendai College. H.K. also acknowledges Professor Jean-Yves Cavaille of INSA Lyon for the support from the IFS project and encouragement. V.E. acknowledges president Dr. Hiroshi Fukumura of Sendai KOSEN for the warm hospitality provided during a four-month stay from 2019 to 2020, which was supported, in part, by JSPS KAKENHI Grant No. JP17K05149.

### APPENDIX: PROOF OF STATEMENT (B) IN SEC. II E

Statement (B) in Sec. II E is as follows:

- (B)
Let

*E*_{e,0}= (*ν*_{e,0},*V*_{e,0},*d*_{e,0}) be a set of parameters that characterize Exp(*E*_{e,0}), and let*S*_{0}be a set of parameters given by*S*_{0}= (*ν*_{0},*V*_{0},*D*_{0},*Δx*_{0},*Δt*_{0}). In this situation, if Exp(*E*_{e,0}) ≃*S*_{0}, then for any experimental data Exp(*E*_{e}), with*E*_{e}= (*ν*_{e},*V*_{e},*d*_{e}) and*D*_{e}, and for any given set of (*n*_{X},*n*_{T}), there exists a unique*D*_{sim}such that Exp(*E*_{e}) ≃ (*ν*_{0},*V*_{0},*D*_{sim},*Δx*_{0},*Δt*_{0}).

To prove statement (B), we first fix the parameters *α* and *β* using the parameter sets *S*_{0} and *E*_{e} such that

Indeed, from the expressions *ν*_{e} (m^{2}/s) = *ν*_{e}*α*^{−2}*β* [(*α*m)^{2}/*β*s] and *V*_{e} (m/s) = *V*_{e}*α*^{−1}*β* (*α*m/*β*s), we have *ν*_{0} = *ν*_{e}*α*^{−2}*β* and *V*_{0} = *V*_{e}*α*^{−1}*β*, which lead to Eq. (A1). Note that *α* and *β* in Eq. (A1) correspond to the unit transformation between *E*_{e} and *S*_{0} and are not always identical to *α*_{0} and *β*_{0} for the unit transformation between *E*_{e,0} and *S*_{0}.

The assumption Exp(*E*_{e,0}) ≃ *S*_{0} implies that there exist parameters *α*_{0}, *β*_{0}, *γ*_{0}, and *δ*_{0} for the scale transformations m → *α*_{0}m, s → *β*_{0}s, *n*_{X} → *γ*_{0}*n*_{X}, and *n*_{T} → *δ*_{0}*n*_{T} such that

The parameter *D*_{e} is assumed to be given in addition to *E*_{e} for Exp(*E*_{e}), and it is not always identical to *D*_{e,0}. The parameters *Δx*_{e} and *Δt*_{e} are defined as

where *Δx*_{0} and *Δt*_{0} are given by

The parameter *γ*/*γ*_{0} is obtained from the constraint that *Δx*_{0} in Eq. (A4), expressed in units of *α*m for the simulation of Exp(*E*_{e}) and expressed in units of *α*_{0}m for the simulation of Exp(*E*_{e,0}), is identical to *Δx*_{0}; thus, we have

This expression indicates that *γ* is uniquely determined because the parameters that appear in this expression are already uniquely given. The uniqueness of the parameter *δ* can be understood from a similar expression obtained through the same procedure,

although we must assume that the macroscopic relaxation time *τ*_{e} is a well-defined quantity in the target experiments corresponding to *E*_{e,0} and *E*_{e}.

Using *D*_{e}, *Δx*_{e}, and *Δt*_{e}, we define a set of parameters *S*_{e} such that

The strength *D*_{sim} of the random force is fixed to

which is unique because the quantities on the RHS are all uniquely given. Thus, we have proven that

The relation in Eq. (A10) implies that *S*_{e} ≡ (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}), which means that Exp(*E*_{e}) can be simulated by (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}) from statement (A), thus concluding the Proof of statement (B).

In terms of its rigor, this proof is insufficient because the macroscopic relaxation time *τ*_{e} is not always explicitly given in actual experimental data corresponding to *E*_{e}. In such a case, the expressions in Eqs. (A5) and (A7) are meaningless. Therefore, this component is studied numerically in Sec. III. The problem to be numerically clarified is whether the scaling properties expressed in Eq. (A10) are correct for the cases of *α* = *β* = 1, *γ* ≠ 1, and *δ* ≠ 1. It is sufficient to check only the case of *δ* ≠ 1; however, the case of *γ* ≠ 1 will also be checked. It is numerically shown in Sec. III that these scaling properties are correct. Thus, we assume that statement (B) is correct.

The lattice spacing *Δx*_{0} in *S*_{0} used in Eq. (A3) is defined in terms of the experimental diameter *d*_{e,0} (m) and the lattice size *n*_{X0}(=*γ*_{0}*n*_{X}), similar to Eq. (A4), such that

in the simulation units for Exp(*E*_{e,0}). Since the diameter *d*_{e,0} appears in *Δx*_{0}, it is not explicitly included in *S*_{0} and *S*_{e}. The parameter *Δx*_{e} on the RHS of Eq. (A10) is not experimental and is simply defined by Eq. (A3). The final parameter in *S*_{e}, i.e., *Δt*_{e}, is also defined by Eq. (A3).

The implications of statement (B) should be emphasized. The meaning of the equivalence between *S*_{e} and (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}), as expressed by *S*_{e} ≡ (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}), is that any experimental data Exp(*E*_{e}) characterized by the parameter *E*_{e} can be simulated with a single set of parameters *S*_{0} = (*ν*_{0}, *V*_{0}, *D*_{0}, *Δx*_{0}, *Δt*_{0}) if *D*_{0} is replaced with *D*_{sim}. To simulate other experimental data Exp($Ee\u2032$), it is sufficient to replace *D*_{sim} with $Dsim\u2032$ such that $Se\u2032$ ≡ (*ν*_{0}, *V*_{0}, $Dsim\u2032$, *Δx*_{0}, *Δt*_{0}).

*S*_{e} ≡ (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}) simply implies that the simulation results with parameter set *S*_{e} are identical to those with parameter set (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}); hence, this equivalence does not always imply that the real experimental data characterized by *E*_{e} are exactly the same as the simulation results obtained with (*ν*_{0}, *V*_{0}, *D*_{sim}, *Δx*_{0}, *Δt*_{0}). The latter problem is related to the fundamental problem of whether the Langevin NS simulation can successfully simulate real physical flows. In this paper, we assume that it can; this is the implication of the assumption that Exp(*E*_{e,0}) ≃ *S*_{0}. However, we emphasize that this assumption is true only because the experimentally observed peaks in the velocity distribution can be reproduced, which is the main result in this paper, as has been shown. Another implication of the assumption Exp(*E*_{e,0}) ≃ *S*_{0} is that the set of parameters in *S*_{0} is already given. Using the parameters in *S*_{0} and *E*_{e}, we obtain *α* and *β* via Eq. (A1) for Exp(*E*_{e}).

Although the kinematic viscosity coefficient *ν* appears in the NS equation given in Eq. (10), statement (B) indicates that the simulation results depend only on *D*, which is understood from the original NS equation in Eq. (5) for the velocity field without the pressure term,

This equation contains two parameters, *ν* and *D*, in the second and final terms, respectively. The first term can be neglected for protoplasmic streaming; this term is irrelevant to the following discussion, although it is included in Eq. (A12). If the final term $\eta \u2192$ is not present, it is easy to confirm that the solution is

which satisfies the boundary conditions in Eq. (18) and is independent of *ν*. Thus, the question is whether this solution is also expected to satisfy Eq. (A12) and to be independent of *ν* in the presence of $\eta \u2192(t)$. Statement (B) indirectly answers this question and shows that the results depend only on *D* in the presence of $\eta \u2192$, although this statement does not always imply that the results are independent of *ν*.