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 Vx = 0 and at a finite Vx(≠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 Vx(≠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.

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 19565 [Fig. 1(a)]. This position dependence of the flow speed was later precisely measured via particle tracking velocimetry by Kikuchi-Mochizuki15 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 

FIG. 1.

(a) The flow velocity V inside a cell. (b) The normalized velocity distribution h(Vx) along the x-direction.

FIG. 1.

(a) The flow velocity V inside a cell. (b) The normalized velocity distribution h(Vx) along the x-direction.

Close modal

In 1974, Mustacich and Ware observed the distribution of the velocity Vx 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 Vx = 0 and the other at a finite Vx(≠0)21–23 [Fig. 1(b)]. Shortly thereafter, the velocity distribution was again measured using the same technique by Sattelle and Buchan24 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/s21 and 72 μm/s22,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.

FIG. 2.

Plots of experimental velocity distributions obtained via laser Doppler velocimetry in (a) Ref. 21 and (b) Refs. 22 and 23. The horizontal axis represents the frequency of the laser light, which is proportional to the fluid velocity.

FIG. 2.

Plots of experimental velocity distributions obtained via laser Doppler velocimetry in (a) Ref. 21 and (b) Refs. 22 and 23. The horizontal axis represents the frequency of the laser light, which is proportional to the fluid velocity.

Close modal

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 fluids34,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.

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

TABLE I.

List of symbols with units and descriptions. The numbers inside the parentheses are the assumed typical values used as inputs to the simulations.

SymbolUnitDescription (assumed typical value)
ψ m2/s Stream function 
ω 1/s Vorticity 
V=(Vx,Vy) m/s Velocity vector 
VB m/s Velocity at the boundary (50 × 10−6
η=(ηx,ηy) m/s2 Gaussian random force 
g=(gx,gy) Gaussian random number of mean 0 and deviation 1 
D m2/s3 Strength of random force (2 × 10−7D ≤ 4 × 10−6
Ddif m2/s Diffusion constant 
τe Macroscopic relaxation time 
μ Ns/m2 Viscosity 
ν m2/s Kinematic viscosity (1 × 10−4
nX Total number of lattice points on one edge (100 ≤ nX ≤ 300) 
nT Total number of iterations per time step 
α Positive number for a unit change in length 
β Positive number for a unit change in time 
γ Positive number for a change in nX 
δ Positive number for a change in nT 
d Cell diameter (500 × 10−6
Δx Lattice spacing 
Δt Discrete time step (8 × 10−9
E  Parameter set (ν, V, D
S  Parameter set (ν, V, D, Δx, Δt
a Diameter of a lump of fluid particles 
Exp(E Experimental data corresponding to E = (ν, V, D
Re Reynolds number 
Pe Péclet number 
Sc Schmidt number 
SymbolUnitDescription (assumed typical value)
ψ m2/s Stream function 
ω 1/s Vorticity 
V=(Vx,Vy) m/s Velocity vector 
VB m/s Velocity at the boundary (50 × 10−6
η=(ηx,ηy) m/s2 Gaussian random force 
g=(gx,gy) Gaussian random number of mean 0 and deviation 1 
D m2/s3 Strength of random force (2 × 10−7D ≤ 4 × 10−6
Ddif m2/s Diffusion constant 
τe Macroscopic relaxation time 
μ Ns/m2 Viscosity 
ν m2/s Kinematic viscosity (1 × 10−4
nX Total number of lattice points on one edge (100 ≤ nX ≤ 300) 
nT Total number of iterations per time step 
α Positive number for a unit change in length 
β Positive number for a unit change in time 
γ Positive number for a change in nX 
δ Positive number for a change in nT 
d Cell diameter (500 × 10−6
Δx Lattice spacing 
Δt Discrete time step (8 × 10−9
E  Parameter set (ν, V, D
S  Parameter set (ν, V, D, Δx, Δt
a Diameter of a lump of fluid particles 
Exp(E Experimental data corresponding to E = (ν, V, D
Re Reynolds number 
Pe Péclet number 
Sc Schmidt number 

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 AC′, are assumed to be periodic. This square region is extracted and shown in Fig. 3(c) to clarify the 2D nature of the domain.

FIG. 3.

(a) Flow directions in the so-called indifferent zone on the surface and section AA′–CC′ of the cylinder at the center. (b) Section AA′–CC′ and the flow directions on the surface of the cylinder. (c) 2D simulation domain corresponding to section AA′–CC′, where VB denotes the fixed velocity considered as a boundary condition.

FIG. 3.

(a) Flow directions in the so-called indifferent zone on the surface and section AA′–CC′ of the cylinder at the center. (b) Section AA′–CC′ and the flow directions on the surface of the cylinder. (c) 2D simulation domain corresponding to section AA′–CC′, where VB denotes the fixed velocity considered as a boundary condition.

Close modal

The boundary condition given by the velocity VB 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 equation34,35 with a random Brownian force is given by

ωt=Vω+νΔω+×η(t)z,ω=Δψ,
(1)

where V=(Vx,Vy,0) is the fluid velocity obtained from the stream function ψ and ω is the third component of the vorticity ω=×V such that

Vx=ψy,Vy=ψx,Vz=0,ω=×V=(0,0,ω).
(2)

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 m2/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

ηiμ(t)ηjν(t)=2Dδijδμνδ(tt),
(3)

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:

ρVt+VV=p+μΔV+ρη(t),
(4)

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

Vt=VVρ1p+νΔV+η(t),V=0.
(5)

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 V=0 is exactly satisfied in Eq. (1); therefore, Eq. (1) is easier to solve numerically than the original NS equation for V 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

Hi,j(t)=tt+Δtηi,j(t)dt,
(6)

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 Hi,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

Hi,j(t)=0,Hi,j2(t)=2DΔt,
(7)

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 Hi,j(t) in Eq. (6) as

Hi,j(t)=ηi,j(t)Δt
(8)

and substitute this Hi,j(t) into the second expression in Eq. (7), we obtain a finite value

ηi,jμ(t)=2D/Δt
(9)

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

ωi,j(t+Δt)ωi,j(t)ΔtVωi,j(t)+νΔtΔωi,j(t)       +2DΔt(×gi,j(t))z,ωi,j=Δψi,j,
(10)

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

ηi,jΔt=2DΔtgi,j(t).
(11)
FIG. 4.

(a) 2D regular square lattice of size N = nX × nX, where a lattice site is represented by (i, j) and the lattice spacing is Δx in both the i and j directions. (b) The stochastic variable H(t) 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.

FIG. 4.

(a) 2D regular square lattice of size N = nX × nX, where a lattice site is represented by (i, j) and the lattice spacing is Δx in both the i and j directions. (b) The stochastic variable H(t) 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.

Close modal

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

ΔtVω=ΔtψyωxψxωyΔt4(Δx)2ψi,j+1ψi,j1ωi+1,jωi1,jψi+1,jψi1,jωi,j+1ωi,j1.
(12)

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

Δωi,j(1/Δx)2ωi+1,j+ωi1,j+ωi,j+1+ωi,j14ωi,j,
(13)

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

(×gi,j(t))z(gi+1,jygi1,jygi,j+1x+gi,j1x)/(2Δx).
(14)

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 Ddif related to the temperature T by means of the Einstein–Stokes–Sutherland formula Ddif = kBT/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 VB [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

2Dτe2(=Ddif)=kBT6πμa,
(15)

where μ(=ρν) is the viscosity, a is the size of a fluid particle or a group of particles in the fluid, and kB 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 Ddif 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

ψi,j(+1)(t)ψi,j()(t)+Aψi+1,j()(t)+ψi1,j()(t)+ψi,j+1()(t)+ψi,j1()(t)+(Δx)2ωi,j(t)4ψi,j()(t),
(16)

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.

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

ωi,nX=2(Δx)2ψi,nX1+VBΔx,ψi,nX=0(onΓ1),ωi,1=2(Δx)2ψi,2+VBΔx,ψi,1=0(onΓ3),
(17)

where the velocity V at the boundary is given by

V=(VB,0)onΓ1,V=(VB,0)onΓ3.
(18)

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

FIG. 5.

(a) The boundary condition ψ = 0 at the boundaries Γ1 and Γ3. (b) The velocity V=(Vx,Vy) is fixed to V=(VB,0) at Γ1 and to V=(VB,0) at Γ3. The lattice sites (i, nX − 1) and (i, 2) in (a), close to Γ1 and Γ3, respectively, are used to enforce the boundary conditions for the variable ω in Eq. (17).

FIG. 5.

(a) The boundary condition ψ = 0 at the boundaries Γ1 and Γ3. (b) The velocity V=(Vx,Vy) is fixed to V=(VB,0) at Γ1 and to V=(VB,0) at Γ3. The lattice sites (i, nX − 1) and (i, 2) in (a), close to Γ1 and Γ3, respectively, are used to enforce the boundary conditions for the variable ω in Eq. (17).

Close modal

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 ψ + f0 are exactly equivalent for any constant f0 in the sense that both ψ and ψ + f0 correspond to the same velocity configuration. Therefore, if ψ1,nX is nonzero such that ψ1,nX=c0 at (1, nX) ∈ Γ1, then ψ can be replaced by ψ + f0 with f0 = −c0, and hence, we have ψ1,nX=0 [Fig. 5(a)]. It is also easy to check that ψi,nX=0(i>1) because of the boundary condition Vy = 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,nX on the boundary Γ1 is obtained by means of Taylor expansion and the second expression in Eq. (1) as follows:

ψ(x,yΔy)=ψ(x,y)ψy(x,y)Δy+122ψy2(x,y)Δy2+=0+VBΔx12ω(x,y)Δx2+for(x,y)Γ1,
(19)

where Δy = Δx is assumed and Vx = −∂ψ/∂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

ωnX+1,j=ω1,j,ψnX+1,j=ψ1,j,ω1,j=ωnX,j,ψ1,j=ψnX,j.
(20)

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

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:

1m=α1αm,1s=β1βs.
(21)

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

νm2/s,VBm/s,dm,τs,Dm2/s3,Δxm,Δts.
(22)

The kinematic viscosity coefficient ν is explicitly included in Eq. (10). The velocity VB 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 m2/s3 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 τ:

Δx=dnX,Δt=τnT,
(23)

where nX is the total number of lattice points on one edge of the lattice and nT 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 (nT) is excessively large (small), then the simulation will not converge; hence, Δt0 should be fixed to Δt0Δtcr. This Δtcr 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

E=(ν,V,d),S=(ν,V,D,Δx,Δt).
(24)

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.

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, nX, and nT introduced in Sec. II C. These scale transformations are defined in terms of positive numbers α, β, γ, and δ(>0) such that

m,s,nX,nTαm,βs,γnX,δnT.
(25)

The first two transformations, (m, s) → (αm, βs), represent changes to the units of length and time. The third transformation, nXγnX, represents a change in lattice size, which corresponds to a change in the lattice spacing Δxγ−1Δx in Eq. (23). The final transformation, nTδnT, represents a change in the time step Δtδ−1Δt by means of Eq. (23). γ and δ for nX and nT 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 nX by α can be used because of the relation Δx (m) = d/nX (m) = α−1d/nX (αm) = d/(αnX) (αm). Indeed, this relation implies that a unit transformation by α can be understood as a change in the lattice size nX. However, in this case, it is impossible to observe the dependence of the results on the lattice size nX 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

ωi,jωi,j+Δt(Δx)2ψω+νΔt(Δx)2ω+2DΔtΔxg(1/s),
(26)

where (Δx)2ψω on the RHS represents the spatial discretization in Eq. (12) and (Δx)2ω and (Δx)1gz represent the discretizations of the Laplacian Δω and the rotation ×gz, 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/nX) (m) and Δt(= τ/nT) (s) are replaced by α−1γ−1Δx (αm) and β−1δ−1Δt (βs), and we also have ν (m2/s) = α−2βν [(αm)2/βs], and D (m2/s3) = α−2β3D [(αm)2/(βs)3]. The units of ψ and ω are m2/s and 1/s, respectively; therefore, we have ψ (m2/s) = α−2βψ [(αm)2/βs] and ω (1/s) = βω (1/βs). From these expressions and Eq. (26), we obtain

ωi,jωi,j+γ2δ1Δt(Δx)2ψω+γ2δ1νΔt(Δx)2ω+2γ2δ1DΔtΔxg(1/βs),
(27)

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δ−1D under the scale transformation (nX, nT) → (γnX, δnT), 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

ωi,jωi,j+γ2δ1Δt(Δx)2ψω+γ2δ1νΔt(Δx)2ω+γ2δ12DΔtΔxg(1/βs);
(28)

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) = α−1d (α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:

(ν,V,D,Δx,Δt)α2βν,α1βV,α2β3γ2δ1D,α1γ1Δx,β1δ1Δt.
(29)

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 Vx distribution and

  • the dependence of the normalized Vx on y

for (ω1, ψ1) are identical to those for (ω2, ψ2).

Thus, we have proven the following statement:

  • The solution (ω, ψ) of Eq. (10) remains unchanged if and only if the parameters S = (ν, V, D, Δx, Δt) scale in accordance with Eq. (29) under the scale transformations in Eq. (25).

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(E1) and Exp(E2) 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 S1 = (ν1, V1, D1, Δx1, Δt1) and S2 = (ν2, V2, D2, Δx2, Δt2) are called equivalent if there exists a set of positive numbers α, β, γ, and δ (all >0) such that

(ν1,V1,D1,Δx1,Δt1)=α2βν2,α1βV2,α2β3γ2δ1D2,        α1γ1Δx2,β1δ1Δt2.
(30)

This equivalence is denoted by S1S2. It is easy to confirm that S2S1 if S1S2 because α and β can be inverted to α−1 and β−1. S1 = S2 if and only if α = β = γ = δ = 1. A set of parameters equivalent to S = (ν, V, D, Δx, Δt) is written as S¯=(ν¯,V¯,D¯,Δx¯,Δt¯).

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

Finally, in this subsection, we introduce the notation

Exp(E)S0

for

E=(ν,V,d),S0=(ν0,V0,D0,Δx0,Δt0),
(31)

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) ≃ S0 is that “the experimental data Exp(E) are successfully simulated with the parameters S0.”

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(Ee) data, while the remaining parameters can be fixed to the parameters S0 used to simulate certain existing experimental data Exp(Ee,0), which are not always identical to Exp(Ee). This process is described in more detail in the following statement:

  • (B)

    Let Ee,0 = (νe,0, Ve,0, de,0) be a set of parameters that characterize Exp(Ee,0), and let S0 be a set of parameters given by S0 = (ν0, V0, D0, Δx0, Δt0). In this situation, if Exp(Ee,0) ≃ S0, then for any experimental data Exp(Ee), with Ee = (νe, Ve, de) and De, and for any given set of (nX, nT), there exists a unique Dsim such that Exp(Ee) ≃ (ν0, V0, Dsim, Δx0, Δt0).

Statement (B) indicates that only one parameter, Dsim, must be varied to simulate arbitrary experimental data Exp(Ee). Proof of statement (B) and further details are provided in the  Appendix.

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–32.

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

ε=1×109,
(32)

for a time step of the NS equation such that

1/Nij1ωij(t+Δt)ωij(t)<ε
(33)

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

1/Nij1ψij(t+Δt)ψij(t)<ε
(34)

from which the convergent configuration {ψ} is obtained. N(=∑ij1 = nX × nX) 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

1/Nij1ψij(+1)(t)ψij()(t)<ε
(35)

at each time step Δt.

Note that ε is related to τ in Eq. (22). Indeed, ε determines the total number of iterations nT, which satisfies nT = τ/Δ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 nT 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

Q=(1/n)k=1nQk,
(36)

where Qk is the kth sample corresponding to the kth 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 Qk in Eq. (36) from the convergent configuration {ω,ψ}k. The total number n of samples lies in the range of 1 × 104n ≤ 2 × 104 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, 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 can then be calculated using this mean value configuration ψ.

  • {V}k may first be calculated from the kth configuration {ω,ψ}k, and V can then be calculated as the mean value of the {V}k(k=1,n).

This is simply because of the commutativity of the calculation of V 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 nT, 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.

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 Dsim = 0, and Fig. 6(b) shows the mean ψ values of 1000 convergent configurations for Dsim = 400 and the V results calculated using this mean ψ. Figures 6(c) and 6(d) show the snapshots of two different convergent configurations for the same Dsim = 400. The parameters other than Dsim 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 (⇔Dsim = 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 Dsim, 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 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.

FIG. 6.

(a) Snapshot of the normalized stream function ψ and the velocity V obtained for Dsim = 0. (b) The mean values of 1000 convergent configurations for Dsim = 400. [(c) and (d)] Snapshots of convergent configurations for Dsim = 400. Small cones represent V, and the colors represent ψ. In the case of Dsim = 0, no random Brownian force is assumed, and a uniquely determined configuration appears. In contrast, in the case of Dsim = 400, all of the convergent configurations are slightly different, two examples of which are illustrated by the snapshots in (c) and (d). If Dsim = 400 is decreased to Dsim = 0, all of the possible configurations continuously deform to that depicted in (a), while if Dsim = 400 is continuously increased to Dsim = 4000, for example, then the vectors in each snapshot become more clearly deformed except in the boundary region; however, no deviation can be observed in the mean values of ψ or the velocity with respect to those depicted in (b) for Dsim = 400.

FIG. 6.

(a) Snapshot of the normalized stream function ψ and the velocity V obtained for Dsim = 0. (b) The mean values of 1000 convergent configurations for Dsim = 400. [(c) and (d)] Snapshots of convergent configurations for Dsim = 400. Small cones represent V, and the colors represent ψ. In the case of Dsim = 0, no random Brownian force is assumed, and a uniquely determined configuration appears. In contrast, in the case of Dsim = 400, all of the convergent configurations are slightly different, two examples of which are illustrated by the snapshots in (c) and (d). If Dsim = 400 is decreased to Dsim = 0, all of the possible configurations continuously deform to that depicted in (a), while if Dsim = 400 is continuously increased to Dsim = 4000, for example, then the vectors in each snapshot become more clearly deformed except in the boundary region; however, no deviation can be observed in the mean values of ψ or the velocity with respect to those depicted in (b) for Dsim = 400.

Close modal

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

Vx=12Δxψi,j+1ψi,j1,Vy=12Δxψi+1,jψi1,j,
(37)

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

FIG. 7.

(a) The distribution h(Vx) of |Vx| along the x-direction. (b) The distribution h(V) of the length of V. The lattice size is fixed to nX = 100, and the strength Dsim of the random Brownian force is varied from Dsim = 0 to Dsim = 4000 in the simulation units. The dashed lines are normalized distributions exp(Vx2/(2c2)) in (a) and V exp(−V2/(2c2)) in (b) expected from Maxwell–Boltzmann distribution corresponding to the ideal gas. In these expressions, c = 0.22 is assumed, and the normalization factor in V exp(−V2/(2c2)) is dropped. Due to the violation of the equipartition of energy, the distribution h(Vx) in (a) deviates from the ideal gas behavior even at the largest Dsim = 4000.

FIG. 7.

(a) The distribution h(Vx) of |Vx| along the x-direction. (b) The distribution h(V) of the length of V. The lattice size is fixed to nX = 100, and the strength Dsim of the random Brownian force is varied from Dsim = 0 to Dsim = 4000 in the simulation units. The dashed lines are normalized distributions exp(Vx2/(2c2)) in (a) and V exp(−V2/(2c2)) in (b) expected from Maxwell–Boltzmann distribution corresponding to the ideal gas. In these expressions, c = 0.22 is assumed, and the normalization factor in V exp(−V2/(2c2)) is dropped. Due to the violation of the equipartition of energy, the distribution h(Vx) in (a) deviates from the ideal gas behavior even at the largest Dsim = 4000.

Close modal

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

The distribution h(Vx) is obtained by constructing a histogram in which the velocity range 0|Vx|vxmax is divided into 100 subranges, and the velocity Vx 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=2VB, where VB denotes the velocity at the boundaries Γ1 and Γ3. The reason for this choice of vxmax is that vx,(i)max in the ith convergent configuration is expected to vary with i instead of remaining constant. The factor of 2 in 2VB 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(Vx) is also independent of procedures (i) and (ii) for the calculation of Vx. The height of h(Vx) is normalized such that the maximum height is equal to 1 for each Dsim, and the horizontal axis |Vx| is similarly normalized using the maximum |Vx| satisfying h(Vx) ≠ 0 for each Dsim. 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(Vx) is flat for Dsim = 0; this flat h(Vx) is compatible with the expectation from Eq. (A13). For a sufficiently large Dsim, the nonflat h(Vx) has a peak only at Vx = 0. However, for intermediate Dsim values, another peak appears at Vx ≠ 0. If all the ensemble configurations for Dsim ≠ 0 were identical to that in Fig. 6(a) for Dsim = 0, then h(Vx) 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(Vx) 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(Vx) 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(=(Vx,0)) but also V(=(Vx,Vy)), which has a small nonzero Vy component, can be detected. This drop of h(V) at V → 0 in Fig. 7(b), for sufficiently large Dsim, is also consistent with the drop observed in the Maxwell–Boltzmann distribution V exp(−V2/(2c2)) expected in the ideal gas. In this expression, the constant c(= 0.22) is obtained from the peak position V of h(V) for Dsim = 4000 in Fig. 7(b). The simulation result h(V) for Dsim = 4000 is almost consistent with the dashed line as expected, even though Dsim is finite. However, the corresponding h(Vx) considerably deviates from the dashed line exp(Vx2/(2c2)) in Fig. 7(a). The reason for this deviation is that the equipartition of energy, Vx2=Vy2, is obviously violated due to the effect of boundary velocity, at least for finite Dsim. Nevertheless, the ideal gas behavior observed in h(V) for sufficiently large Dsim is also reasonable because Dsim 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 Dsim(≠0) are located between exp(Vx2/(2c2)) (dashed line) expected from the ideal gas and h(Vx) = 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/cm2). 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 m2/s to νe,0 = 1.5 × 10−4 m2/s. Thus, we assume that νe,0 = 1 × 10−4 m2/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 Ve,0(= 50 μm/s) at the boundary and a diameter of de,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 Ee,0. The corresponding α0 and β0 values are obtained via Eq. (A1) using the parameters S0 shown in Table III.

TABLE II.

Physical parameters Ee,0 = (νe,0, Ve,0, de,0) assumed in the simulations, expressed in physical units, and the parameters α0 and β0 for the conversion to the simulation units.

νe,0 (m2/s)Ve,0 (μm/s)de,0 (μm)α0β0
1 × 10−4 50 500 1 × 10−6 1 × 10−1 
νe,0 (m2/s)Ve,0 (μm/s)de,0 (μm)α0β0
1 × 10−4 50 500 1 × 10−6 1 × 10−1 
TABLE III.

The parameters S0 used in the simulations; these values are given in the simulation units. γ0nX is written as nX0 for simplicity.

ν0[(α0m)2β0s]V0[α0mβ0s]Δx0 (α0m)Δt0 (β0s)nX0
1 × 107 8 × 10−8 100 
ν0[(α0m)2β0s]V0[α0mβ0s]Δx0 (α0m)Δt0 (β0s)nX0
1 × 107 8 × 10−8 100 

The lattice spacing Δx0 (α0m) is calculated using Eq. (A11). For Δx0 = 5 α0m, we have Δxe,0=Δx0α01=5×106m.

There is no need to mention that the parameters S0 in Table III are also obtained from the parameters Ee,0 in Table II and α0 and β0. Indeed, it is easy to check that the kinematic viscosity ν0=1×107α02m2/β0s is obtained from the relation νe,0(m2/s)=1×104β0/α02(α02m2/β0s)=ν0(α02m2/β0s). The velocity V0 (α0m/β0s) and d0 (α0m) are also obtained through the relations Veα01β0(α0m/β0s)=V0 and deα01=d0(α0m), which are not explicitly used in the simulations.

Each value of Dsim shown in Table IV is fixed as an input to the simulations, and the corresponding physical strength De,0 can be obtained as De,0=Dsimα02β03. The parameter τ0 in Table IV is estimated using the relation τ0 = nTΔt0, where τ0 is expected to satisfy τ0nTΔtcr, as discussed in Sec. II C. Using these De,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 kBT = 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 = nTΔt0 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.

TABLE IV.

The assumed parameters Dsim and De,0, the results τ0 estimated as τ0 = nTΔt0, and the diameter a of a fluid particle as estimated using Eq. (15).

Dsim[(α0m)2(β0s)3]De,0 (m2/s3)nTΔt0 (β0s)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[(α0m)2(β0s)3]De,0 (m2/s3)nTΔt0 (β0s)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 Ddif in Eq. (15) can be estimated by assuming that nTΔt0 is τe as follows: Ddif=2De,0(nTΔt0β0)2. For De,0 = 4 × 107 m2/s3 in Table IV, we have Ddif ≃ 4 × 10−11 m2/s = 4 × 10−7 cm2/s, which is almost 10 times smaller than the estimated value of Ddif = kBT/6πμa ≃ 5.3 × 10−6 cm2/s obtained with μ = 0.1 m2/s and a = 4 × 10−10 m. Note, however, that this value of 5.3 × 10−6 cm2/s is comparable to the value of 10−5 cm2/s reported in Ref. 16. The deviation between the estimates Ddif=2Dτe2 and Ddif = kBT/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 nTΔt0β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 Sc, which is calculated as the ratio of the Reynolds number Re and the Péclet number Pe such that Sc = Pe/Re. The Reynolds number is evaluated to be Re = 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 m2/s, respectively. The Péclet number is similarly evaluated to be Pe = Vd/Ddif ≃ 10 for Ddif = 5.3 × 10−10 m2/s. If we assume that Ddif=4×1011m2/s(=2De,0(nTΔt0β0)2), we obtain Pe ≃ 100, which is closer to the estimate of 102–103 given in Ref. 16. Therefore, we obtain Sc ≃ 2 × 105 for Ddif = 5.3 × 10−10 m2/s and Sc ≃ 2 × 106 for Ddif = 4 × 10−11 m2/s. Note also that Sc can be obtained directly from its definition, Sc = ν/Ddif. Thus, the relatively large Sc 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 term16,17 but also to the viscosity term. Nevertheless, the peaks in the velocity distribution essentially originate from this small thermal fluctuation effect.

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

FIG. 8.

(a) Distributions h(Vx) obtained on lattices with sizes ranging from nX = 100 to nX = 300. (b) The dependence of Vx on y. The parameters Dsim and Δx0 are varied in accordance with Eq. (29), and Dsim, Δx0, and the other parameters d0, V0, ν0, and Δt0 shown in (a) and (b) are all expressed in the simulation units. Dashed lines in (b) denoted by Exp. are normalized experimental data reported in Refs. 15 and 18.

FIG. 8.

(a) Distributions h(Vx) obtained on lattices with sizes ranging from nX = 100 to nX = 300. (b) The dependence of Vx on y. The parameters Dsim and Δx0 are varied in accordance with Eq. (29), and Dsim, Δx0, and the other parameters d0, V0, ν0, and Δt0 shown in (a) and (b) are all expressed in the simulation units. Dashed lines in (b) denoted by Exp. are normalized experimental data reported in Refs. 15 and 18.

Close modal

The height of the histogram h(Vx) is normalized such that the maximum height is equal to 1 for each Dsim, as in Fig. 7(a), while the horizontal axes |Vx| for all Dsim are normalized using a constant value equal to the maximum |Vx| for Dsim = 200 and a lattice size of nX = 100 satisfying h(Vx) ≠ 0.

We find that h(Vx) with respect to |Vx| is almost independent of either nX or Δx. This finding supports not only the correctness of statement (B) but also the Δx independence of the results. In fact, if h(Vx) with respect to |Vx| 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 Vx 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 Vx is found simply by averaging Vx from the convergent configurations using Eq. (37). Thus, the nontrivial behavior of two distinct peaks in h(Vx) is not always reflected in the dependence of Vx on y. The parameters used in the simulations are shown in Fig. 8(b).

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 nT, 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(Vx) with respect to |Vx| is defined in the same way as in Fig. 8(a).

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

FIG. 9.

(a) Distributions h(Vx) with respect to |Vx| for different combinations of Δt0 and Dsim. (b) nTΔt0 vs Δt0. In (a), the parameters Δt0 and Dsim are scaled to δ−1Δt0 and δ−1Dsim, where Δt0 is varied in the range 2 × 10−8Δt0 ≤ 32 × 10−8 with the simulation units (1 β0s = 0.1 s).

FIG. 9.

(a) Distributions h(Vx) with respect to |Vx| for different combinations of Δt0 and Dsim. (b) nTΔt0 vs Δt0. In (a), the parameters Δt0 and Dsim are scaled to δ−1Δt0 and δ−1Dsim, where Δt0 is varied in the range 2 × 10−8Δt0 ≤ 32 × 10−8 with the simulation units (1 β0s = 0.1 s).

Close modal

Finally, in this subsection, we check whether nTΔt0 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 nTΔt0 on ε. The results shown in Fig. 9(b) indicate that nTΔt0 is roughly independent of ε for sufficiently small Δt0.

In this subsection, we demonstrate how to use statement (B) to discuss the dependence of the normalized velocity distribution on the physical parameters νe, Ve, and de. From the perspective of statement (B), the purpose of the simulations in Sec. III B is to find S0 = (ν0, V0, D0, Δx0, Δt0) with a suitable D0 and the set of parameters listed in Table III. Indeed, from the simulation results, we find that D0 = 400 in the second line of Table IV is suitable because the shape of h(Vx) in Fig. 7(a) is relatively close to the experimental data reported in Refs. 21–23. Since |Vx| 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 Dsim = 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 S0, we can perform simulations of Exp(Ee) characterized by Ee, which is different from Ee,0.

The parameters are shown in Table V, where νe, Ve, and de are the elements of the experimental data E(i) (i = 1, 2, 3). E(1), E(2), and E(3) are different from Ee,0 in Table II only in terms of νe, Ve, and de, respectively (indicated by underlines). Note that νe/νe,0 = 2, Ve/Ve,0 = 2, and de/de,0 = 2 imply that νe = 2νe,0 = 2 × 10−4 m2/s, Ve = 2Ve,0 = 100 μm/s, and de = 2de,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, Ve,0, and de,0 in Table II depending on their size.2 

TABLE V.

Experimental data E(i) (i = 1, 2, 3), τe, De, and the corresponding parameters α, β, γ, δ, and Dsim assumed for the simulation parameters Se(i) (i = 1, 2, 3). These parameters are determined by the ratio Dsim/D0, where D0 = 400 from the second line of Table IV.

E(i)νeνe,0VeVe,0dede,0τeτe,0DeDe,0αα0ββ0γγ0δδ0DsimD0
(1) 2̲ 12 12 14 
(2) 2̲ 12 14 116 
(3) 2̲ 14 14 
E(i)νeνe,0VeVe,0dede,0τeτe,0DeDe,0αα0ββ0γγ0δδ0DsimD0
(1) 2̲ 12 12 14 
(2) 2̲ 12 14 116 
(3) 2̲ 14 14 

We assume that the macroscopic relaxation time τe in Eq. (15) is proportional to the inverse kinematic viscosity νe139 and the area Ae as follows:

τeAe/νe.
(38)

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

DekBTμaτe2de2νe,
(39)

where Ae 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 Dsim/D0 via statement (B). Since Exp(Ee,0) is simulated with S0 = (ν0, V0, D0, Δx0, Δt0), from statement (B), we have

(ν0,V0,D0,Δx0,Δt0)=α02β0νe,0,α01β0Ve,0,α02β03γ02δ01De,0,  α01γ01Δxe,0,β01δ01Δte,0,
(40)

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

(ν0,V0,Dsim,Δx0,Δt0)=α2βνe,α1βVe,α2β3γ2δ1De  α1γ1Δxe,β1δ1Δte,
(41)

using the parameters α, β, γ, and δ for the scale transformation between S0 and Se. Therefore, from Table V, we have

α02β0=2α2β,α01β0=α1βforE(1),α02β0=α2β,α01β0=2α1βforE(2),α02β0=α2β,α01β0=α1βforE(3),
(42)

which imply

α/α0=2,β/β0=2forE(1),α/α0=1/2,β/β0=1/4forE(2),α/α0=1,β/β0=1forE(3).
(43)

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

DsimD0=αα02ββ03γγ02δδ01DeDe,0.
(44)

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

Dsim=1600(E(1)),25(E(2)),100(E(3))
(45)

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 Ee,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 Re(=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 Re.

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 S0 correspond to Exp(Ee,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 Dsim = 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, Ve, De), such as those in Table IV, are identical to or located between the peaks in Fig. 7(a).

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

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

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.

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

  • (B)

    Let Ee,0 = (νe,0, Ve,0, de,0) be a set of parameters that characterize Exp(Ee,0), and let S0 be a set of parameters given by S0 = (ν0, V0, D0, Δx0, Δt0). In this situation, if Exp(Ee,0) ≃ S0, then for any experimental data Exp(Ee), with Ee = (νe, Ve, de) and De, and for any given set of (nX, nT), there exists a unique Dsim such that Exp(Ee) ≃ (ν0, V0, Dsim, Δx0, Δt0).

To prove statement (B), we first fix the parameters α and β using the parameter sets S0 and Ee such that

α=νeν0V0Ve,β=νeν0V0Ve  2.
(A1)

Indeed, from the expressions νe (m2/s) = νeα−2β [(αm)2/βs] and Ve (m/s) = Veα−1β (αm/βs), we have ν0 = νeα−2β and V0 = Veα−1β, which lead to Eq. (A1). Note that α and β in Eq. (A1) correspond to the unit transformation between Ee and S0 and are not always identical to α0 and β0 for the unit transformation between Ee,0 and S0.

The assumption Exp(Ee,0) ≃ S0 implies that there exist parameters α0, β0, γ0, and δ0 for the scale transformations m → α0m, s → β0s, nXγ0nX, and nTδ0nT such that

(ν0,V0,D0,Δx0,Δt0)=α02β0νe,0,α01β0Ve,0,α02β03γ02δ01De,0,  α01γ01Δxe,0,β01δ01Δte,0.
(A2)

The parameter De is assumed to be given in addition to Ee for Exp(Ee), and it is not always identical to De,0. The parameters Δxe and Δte are defined as

Δxe=γΔx0α,Δte=δΔt0β,
(A3)

where Δx0 and Δt0 are given by

Δx0=deγnXα1=deγnXν0νeVeV0(αm),
(A4)
Δt0=τeδnTβ1=τeδnTν0νeVeV02(βs).
(A5)

The parameter γ/γ0 is obtained from the constraint that Δx0 in Eq. (A4), expressed in units of αm for the simulation of Exp(Ee) and expressed in units of α0m for the simulation of Exp(Ee,0), is identical to Δx0; thus, we have

γγ0=dede,0νe,0νeVeVe,0.
(A6)

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,

δδ0=τeτe,0νe,0νeVeVe,02,
(A7)

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

Using De, Δxe, and Δte, we define a set of parameters Se such that

Se=(νe,Ve,De,Δxe,Δte).
(A8)

The strength Dsim of the random force is fixed to

Dsim=α2β3γ2δ1De,
(A9)

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

(ν0,V0,Dsim,Δx0,Δt0)=α2βνe,α1βVe,α2β3γ2δ1De,  α1γ1Δxe,β1δ1Δte.
(A10)

The relation in Eq. (A10) implies that Se ≡ (ν0, V0, Dsim, Δx0, Δt0), which means that Exp(Ee) can be simulated by (ν0, V0, Dsim, Δx0, Δt0) 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 Ee. 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 Δx0 in S0 used in Eq. (A3) is defined in terms of the experimental diameter de,0 (m) and the lattice size nX0(=γ0nX), similar to Eq. (A4), such that

Δx0=de,0nX0α01=de,0nX,0ν0νe,0Ve,0V0(α0m)
(A11)

in the simulation units for Exp(Ee,0). Since the diameter de,0 appears in Δx0, it is not explicitly included in S0 and Se. The parameter Δxe on the RHS of Eq. (A10) is not experimental and is simply defined by Eq. (A3). The final parameter in Se, i.e., Δte, is also defined by Eq. (A3).

The implications of statement (B) should be emphasized. The meaning of the equivalence between Se and (ν0, V0, Dsim, Δx0, Δt0), as expressed by Se ≡ (ν0, V0, Dsim, Δx0, Δt0), is that any experimental data Exp(Ee) characterized by the parameter Ee can be simulated with a single set of parameters S0 = (ν0, V0, D0, Δx0, Δt0) if D0 is replaced with Dsim. To simulate other experimental data Exp(Ee), it is sufficient to replace Dsim with Dsim such that Se ≡ (ν0, V0, Dsim, Δx0, Δt0).

Se ≡ (ν0, V0, Dsim, Δx0, Δt0) simply implies that the simulation results with parameter set Se are identical to those with parameter set (ν0, V0, Dsim, Δx0, Δt0); hence, this equivalence does not always imply that the real experimental data characterized by Ee are exactly the same as the simulation results obtained with (ν0, V0, Dsim, Δx0, Δt0). 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(Ee,0) ≃ S0. 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(Ee,0) ≃ S0 is that the set of parameters in S0 is already given. Using the parameters in S0 and Ee, we obtain α and β via Eq. (A1) for Exp(Ee).

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,

Vt=VV+νΔV+η(t).
(A12)

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

V=2VBdy,0,
(A13)

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

1.
N.
Kamiya
, “
Cytoplasmic streaming in giant algal cells: A historical survey of experimental approaches
,”
Bot. Mag., Tokyo
99
,
441
496
(
1986
).
2.
N.
Kamiya
and
K.
Kuroda
, “
Dynamics of cytoplasmic streaming in a plant cell
,”
Biorheology
10
,
179
187
(
1973
).
3.
M.
Tazawa
, “
Motive force of the cytoplasmic streaming in Nitella
,”
Protoplasma
65
,
207
222
(
1968
).
4.
N.
Kamiya
and
K.
Kuroda
, “
Measurement of the motive force of the protoplasmic rotation in Nitella
,”
Protoplasma
50
,
144
147
(
1958
).
5.
N.
Kamiya
and
K.
Kuroda
, “
Velocity distribution of the protoplasmic streaming in Nitella cells
,”
Bot. Mag., Tokyo
69
,
544
554
(
1956
).
6.
D.
Houtman
,
I.
Pagonabarraga
,
C. P.
Lowe
,
A.
Esseling-Ozdoba
,
A. M. C.
Emons
, and
E.
Eiser
, “
Hydrodynamic flow caused by active transport along cytoskeletal elements
,”
Europhys. Lett.
78
,
18001
(
2007
).
7.
S.
Klumpp
,
T. M.
Nieuwenhuizen
, and
R.
Lipowsky
, “
Movements of molecular motors: Ratchets, random walks and traffic phenomena
,”
Physica E
29
,
380
389
(
2005
).
8.
R.
Lipowsky
,
Y.
Chai
,
S.
Klumpp
,
S.
Liepelt
, and
M. J. I.
Müller
, “
Molecular motor traffic: From biological nanomachines to macroscopic transport
,”
Physica A
372
,
34
51
(
2006
).
9.
T.
Kawakubo
,
T.
Kobayashi
, and
S.
Sakamoto
, “
Drift motion of granules in chara cells induced by random impulses due to the myosin-actin interaction
,”
Physica A
248
,
21
27
(
1998
).
10.
F.
Jülicher
,
A.
Ajdari
, and
J.
Prost
, “
Modeling molecular motors
,”
Rev. Mod. Phys.
69
(
4
),
1269
1281
(
1997
).
11.
R. D.
Astumian
, “
Thermodynamics and kinetics of a Brownian motor
,”
Science
276
(
5314
),
9217
9922
(
1997
).
12.
M.
Tominaga
,
A.
Kimura
,
E.
Yokota
,
T.
Haraguchi
,
T.
Shimmen
,
K.
Yamamoto
,
A.
Nakano
, and
K.
Ito
, “
Cytoplasmic streaming velocity as a plant size determinant
,”
Dev. Cell
27
,
345
352
(
2013
).
13.
B. B.
McIntosh
and
E. M.
Ostap
, “
Myosin-I molecular motors at a glance
,”
J. Cell Sci.
129
,
2689
2695
(
2016
).
14.
M.
Tominaga
and
K.
Ito
, “
The molecular mechanism and physiological role of cytoplasmic streaming
,”
Curr. Opin. Plant Biol.
27
,
104
110
(
2015
).
15.
K.
Kikuchi
and
O.
Mochizuki
, “
Diffusive promotion by velocity gradient of cytoplasmic streaming (CPS) in Nitella internodal cells
,”
PLoS One
10
,
e0144938
(
2015
).
16.
J.-W.
Meent
,
I.
Tuval
, and
R. E.
Goldstein
, “
Nature’s microfluidic transporter: Rotational cytoplasmic streaming at high Péclet numbers
,”
Phys. Rev. Lett.
101
,
178102
(
2008
).
17.
R. E.
Goldstein
,
I.
Tuvalk
, and
J.-W.
van de Meent
, “
Microfluidics of cytoplasmic streaming and its implications for intracellular transport
,”
Proc. Natl. Acad. Sci. U. S. A.
105
,
3663
3667
(
2008
).
18.
J.-W.
van De Meent
,
A. J.
Sederman
,
L. F.
Gladden
, and
R. E.
Goldstein
, “
Measurement of cytoplasmic streaming in single plant cells by magnetic resonance velocimetry
,”
J. Fluid Mech.
642
,
5
14
(
2010
).
19.
R. E.
Goldstein
and
J.-W.
van de Meent
, “
A physical perspective on cytoplasmic streaming
,”
Interface Focus.
5
,
20150030
(
2015
).
20.
R.
Niwayama
,
K.
Shinohara
, and
A.
Kimura
, “
Hydrodynamic property of the cytoplasm is sufficient to mediate cytoplasmic streaming in the Caenorhabiditis elegans embryo
,”
Proc. Natl. Acad. Sci. U. S. A.
108
,
11900
11905
(
2011
).
21.
R. V.
Mustacich
and
B. R.
Ware
, “
Observation of protoplasmic streaming by laser-light scattering
,”
Phys. Rev. Lett.
33
,
617
620
(
1974
).
22.
R. V.
Mustacich
and
B. R.
Ware
, “
A study of protoplasmic streaming in Nitella by laser Doppler spectroscopy
,”
Boiophys. J.
16
,
373
388
(
1976
).
23.
R. V.
Mustacich
and
B. R.
Ware
, “
Velocity distributions of the streaming protoplasm in Nitella flexilis
,”
Boiophys. J.
17
,
229
241
(
1977
).
24.
D. B.
Sattelle
and
P. B.
Buchan
, “
Cytoplasmic streaming in Chara corallina studied by laser light scattering
,”
J. Cell Sci.
22
,
633
643
(
1976
), available at https://jcs.biologists.org/content/22/3/633.
25.
D. S.
Lemons
and
A.
Gythiel
, “
Paul Langevin’s 1908 paper “On the Theory of Brownian Motion” [“Sur la théorie du mouvement Brownien,” C. R. Acad. Sci. (Paris) 146, 530–533 (1908)]
,”
Am. J. Phys.
65
(
11
),
1079
1081
(
1997
).
26.
W.
Brenig
,
Statistical Theory of Heat
(
Springer-Verlag Berlin Heidelberg
,
1989
).
27.
R.
Metzler
and
J.
Klafter
, “
The random walk’s guide to anomalous diffusion: A fractional dynamics approach
,”
Phys. Rep.
339
,
1
77
(
2000
).
28.
G. G.
Batrouni
,
G. R.
Katz
,
A. S.
Kronfeld
,
G. P.
Lepage
,
B.
Svetitsky
, and
K. G.
Wilson
, “
Langevin simulations of lattice field theories
,”
Phys. Rev. D
32
,
2736
2747
(
1985
).
29.
A.
Ukawa
and
M.
Fukugita
, “
Langevin simulation including dynamical quark loops
,”
Phys. Rev. Lett.
55
,
1854
1857
(
1985
).
30.
K.
Höfler
and
S.
Schwarzer
, “
Navier-Stokes simulation with constraint forces: Finite-difference method for particle-laden flows and complex geometries
,”
Phys. Rev. E
61
,
7146
7160
(
2000
).
31.
H.
Koibuchi
, “
Langevin simulation of the inter-quark potential in SU(2) lattice gauge theory
,”
J. Phys. G: Nucl. Phys.
13
,
1463
1468
(
1987
).
32.
S.
Nagahiro
,
S.
Kawano
, and
K.
Kotera
, “
Separation of long DNA chains using a nonuniform electric field: A numerical study
,”
Phys. Rev. E
75
(
1-5
),
011902
(
2007
).
33.
B.
Uma
,
T. N.
Swaminathan
,
R.
Radhakrishnan
,
D. M.
Eckmann
, and
P. S.
Ayyaswamy
, “
Nanoparticle Brownian motion and hydrodynamic interactions in the presence of flow fields
,”
Phys. Fluids
23
,
073602
(
2011
).
34.
M. E.
Taylor
,
Partial Differential Equations III: Nonlinear Equations
, 2nd ed. (
Springer
,
New York
,
2010
), Chap. 17, pp.
511
614
.
35.
G.
Lukaszewicz
and
P.
Kalita
,
Navier-Stokes Equations: An Introduction with Applications
(
Springer
,
2015
).
36.
W. T.
Coffey
and
Y. P.
Kalmykov
, “
On the calculation of the macroscopic relaxation time from the Langevin equation for a dipole in a cavity in a dielectric medium
,”
Chem. Phys.
169
,
165
172
(
1993
).
37.
Y.
Feldman
,
A.
Puenko
, and
Y.
Ryabov
, “
Dielectric relaxation phenomena in complex materials
,” in
Fractals, Diffusion, and Relaxation in Disordered Complex Systems
, Advanced Chemical Physics Vol.133, edited by
W. T.
Coffey
and
Y. P.
Kalmykov
(
Wiley-Interscience
,
NJ
,
2006
).
38.
V. I.
Arkhipov
, “
Hierarchy of dielectric relaxation times in water
,”
J. Non-Cryst. Solids
305
,
127
135
(
2002
).
39.
L. I.
Zaichik
,
V. A.
Pershukov
,
M. V.
Kozelev
, and
A. A.
Vinberg
, “
Modeling of dynamics, heat transfer, and combustion in two-phase turbulent flows: 1. Isothermal flows
,”
Exp. Therm. Fluid Sci.
15
,
291
310
(
1997
).
40.
S.
Succi
,
The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond
, Numerical Mathematics and Scientific Computation (
Clarendon Press
,
Oxford
,
2001
).
41.
N.
Metropolis
,
A. W.
Rosenbluth
,
M. N.
Rosenbluth
,
A. H.
Teller
, and
E.
Teller
, “
Equation of state calculations by fast computing machines
,”
J. Chem. Phys.
21
,
1087
1092
(
1953
).
42.
D. P.
Landau
, “
Finite-size behavior of the Ising square lattice
,”
Phys. Rev. B
13
,
2997
3011
(
1976
).
43.
Intel Fortran Compilers for Windows and Linux are used on Windows and Linux PCs. The random numbers are generated by using a Fortran source code Mersenne Twister from home page http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/emt.html.