Stochastic Fluid Dynamics Simulations of the Velocity Distribution in Protoplasmic Streaming

Protoplasmic streaming in plant cells is directly visible in the cases of \textit{Chara corallina} and \textit{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(\not=\!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 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(\not=\!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.

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 one at a finite V x ( = 0) [21][22][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 experimental data extracted from Refs. [21][22][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 the 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. (a) Refs. [22,23] Intensity 0 100 200 Frequency [Hz] (b) Intensity Ref. [21] FIG. 2. Plots of experimental velocity distributions obtained via laser Doppler velocimetry in (a) Ref. [21] and (b) Refs. [22,23]. The horizontal axis represents the frequency of the laser light, which is proportional to the fluid velocity.
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][26][27][28][29][30][31][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 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 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][26][27][28][29][30][31][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.     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 as-sumed 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 random Brownian force is given by where V = (V x , V y , 0) is the fluid velocity obtained from the stream function ψ and ω is the The physical meaning of each term in Eq. (1) will be given below in reference to the equation for the velocity V . 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 η(t) = (η x , η 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 η(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 η(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 . 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 ∇ · V = 0, we obtain  [30]. This is the reason why the NS equation for the stream function, i.e., Eq.
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), η i → η i,j accordingly. This H i,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 ሺ݅ǡ ݆ሻ is understood to have an expectation value of zero; however, its square is finite (Eq. (7)). The force η(t) randomly fluctuates inside a narrow square of width ∆t in (b), and its time integral H(t) can be intuitively understood as an impulse. The macroscopic relaxation time τ e is considerably longer than the width ∆t.
which are typical characteristics of stochastic variables. The first relation comes from the fact that η ij 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 H i,j (t) = η i,j (t)∆t (8) 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 i,j (t) = (g x i,j (t), g y i,j (t), 0) and the components of g i,j are given by Gaussian random numbers with mean 0 and variance 1. This g i,j (t) is related to η i,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∆t in Eq. (10) effectively corresponds to the deviation of the random Brownian force. As determined through dimensional analysis, ( √ 2D∆t) 2 ∆t = 2D(∆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 2D(∆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][37][38][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 the results section.
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 the results section.

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 at the boundary is given by Note that the third component of V 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 ψ 1,n X is nonzero such that ψ 1,n X = c 0 at (1, n X ) ∈ Γ 1 , then ψ can be replaced by ψ + f 0 with f 0 = −c 0 , and hence, we have ψ 1,n X = 0 ( Fig. 5(a)).
It is also easy to check that ψ i,n X = 0 (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 ω i,n X 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 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.
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 the results section. 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 the results section.
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 Sect. II C. These scale transformations are defined in terms of positive numbers α, β, γ, and δ(> 0) such that m, s, n X , n T → αm, βs, γn X , δn T .
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 . 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 Section II A.
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 γ 2 δ −1 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.
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 as follows: Two sets of simulation data (ω 1 , ψ 1 ) and (ω 2 , ψ 2 ) are equivalent if the following conditions are satisfied: (i) The histogram of the normalized velocity V x distribution and (ii) the 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 Sect. II C 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 in 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][26][27][28][29][30][31][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 . Finally, in this subsection, we introduce the notation 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.

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 η is zero, and whether the results obtained under the condition η = 0 are compatible with the expected solution will be determined below. In the case of η = 0, 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][26][27][28][29][30][31][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) 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 ten, i.e., replaced with 10ε, the results will remain unchanged. The dependence on the time step ∆t will be checked separately in Sec. 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 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][26][27][28][29][30][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 , for example, is independent of the order of the calculations: (i) The mean value ψ may first be calculated from the configurations {ω, ψ} k (k = 1, n), and V can then be calculated using this mean value configuration ψ. 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.
which are the standard discrete forms corresponding to the first-order differentials in Eq. (2).
The lattice size is fixed to N = 10000, with N = n X ×n X and n X = 100. According to  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 (t) = 0, 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 (= (V x , 0)) but also V (= (V x , V y )), 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 /(2c 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(−V 2 x /(2c 2 )) in Fig. 7(a). The reason for this deviation is that the equipartition of energy, V 2 x = V 2 y , 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(−V 2 x /(2c 2 )) (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   TABLE II 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 the following Table III.
The lattice spacing ∆x 0 [α 0 m] is calculated using Eq. (A11). For ∆x 0 = 5 α 0 m, we have ∆x e,0 = ∆x 0 α −1 0 = 5 × 10 −6 m. 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 Table IV is fixed as an input to the simulations, and the corresponding physical strength D e,0 can be obtained as D e,0 = D sim α 2 0 β −3 0 . 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 Section 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.
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 D dif = 2Dτ 2 e 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 = V d/ν = 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 = V d/D dif ≃ 10 for D dif = 5.3×10 −10 m 2 /s. If we assume that D dif = 4×10 −11 m 2 /s(= 2D e,0 (n T ∆t 0 β 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. Exp. [18] Exp. [15] FIG . 8. (a) Distributions h(V x ) obtained on lattices with sizes ranging from n X = 100 to n X = 300.  [15,18].

C. Lattice size dependence
In this subsection, we show that the simulation results are independent of the lattice size n X . As mentioned in Section 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 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 Section 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). 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 . 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 Section 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][22][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) 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 Table II only in terms of ν e , V e , and d e , respectively (indicated by underlines). Note that ν e /ν e,0 = 2,  Table II depending on their size [2].   Table IV.
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 d 2 e , 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 Section 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 (= V d/ν) 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 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 Section III B. Furthermore, in Section 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 the final part of Section 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.

V. DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. (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

ACKNOWLEDGMENTS
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 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 .

(A2)
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 ∆x e = γ∆x 0 α, ∆t e = δ∆t 0 β, where ∆x 0 and ∆t 0 are given by 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 S e = (ν e , V e , D e , ∆x e , ∆t e ).
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 the results section. 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 shown that these scaling properties are correct. Thus, we assume that statement (B) is correct.
The , it is sufficient to replace D sim with D ′ sim such that S ′ e ≡ (ν 0 , V 0 , D ′ sim , ∆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.
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 η 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 η(t). Statement (B) indirectly answers this question and shows that the results depend only on D in the presence of η, although this statement does not always imply that the results are independent of ν.