text:
# 1 Introduction

## 1 Introduction

Radiative effects accompany any processes with charged light particle scattering. These effects are exactly calculable in QED and they should be incorporated as an essential step in the data analysis. For elastic processes, which are the subject of the present report, our starting point for the radiative effect calculation is the technique proposed by Bardin and Shumeiko  and later developed for polarized particles as well as for specific elastic and exclusive processes in papers . For other approaches, please see Refs. . In some cases the radiative effects can be calculated analytically. However, often the analytical integration over the photonic phase space is complicated because of specifics of detector geometry or is not desirable because of experimental setup (for example, detector resolution and kinematic cuts). In these situations an approach based on Monte Carlo simulation of the radiative effects is the most adequate. Monte Carlo generator RADGEN  for deep inelastic processes was developed on the basis of corresponding analytical code POLRAD  and successfully applied by several experimental collaborations analyzing polarized and unpolarized observables in deep inelastic electron–nucleon scattering.

The reason for careful implementation of the radiative correction procedure in elastic $`ep`$–scattering is twofold. First, elastic $`ep`$–processes are commonly used for detector calibration, and JLab CLAS detector is one of the examples. Second, a new generation of experiments on elastic $`ep`$-scattering provide data with very high accuracy (see, for instance, Refs. ), and set new precision requirements to data analysis and, in particular, to the radiative correction procedure. In this report we present a new Monte Carlo generator for elastic unpolarized $`ep`$-scattering events taking into account QED radiative effects. The main requirement for construction of the generator is appearance of simulated events (Born and radiative) as well as all event distributions with respect to the kinematic variables in accordance with their contribution to the total cross section.

The contribution to the total cross section from the Born process

$$e(k_1)+p(p_1)e^{}(k_2)+p^{}(p_2)$$
(1)
is described by the transfer momentum squared $`Q^2=(k_1k_2)^2`$ and the azimuthal angle $`\varphi `$. The cross section does not depend on the angle $`\varphi `$ up to acceptance effects. To describe the phase space of the radiative process

$$e(k_1)+p(p_1)e^{}(k_2)+p^{}(p_2)+\gamma (k),$$
(2)
three new kinematic variable have to be considered. We choose them as the proton transfer momentum squared $`t=(k_1k_2k)^2`$, the inelasticity $`v=(p_2+k)^2M^2`$ and the azimuthal angle $`\varphi _k`$ between the planes $`(𝐪,𝐤_\mathrm{𝟏})`$ and $`(𝐤_\mathrm{𝟏},𝐤_\mathrm{𝟐})`$. This set of the variables defines the four–vectors of all final particles in any frame. Usually Lab system is used with OZ axis along the beam direction and plane OZX parallel to the ground level.

Thus, the strategy for simulation of one event can be defined as follows.

* For the fixed initial energy and $`Q^2`$ the non–radiative and radiative parts of the total cross section are calculated. Note, that together with Born subprocess, the contributions of loop effects and soft photon radiation are also included into the non–radiative part.
* The channel of scattering is simulated for this event in accordance with partial contributions of these two positive parts into the total cross section.
* For any process $`\varphi `$ is simulated uniformly from $`0`$ to $`2\pi `$.
* For the radiative event the kinematic variables $`t`$, $`v`$ and $`\varphi _k`$ are simulated in accordance with their calculated distribution (see discussion below).
* The four–momenta of all final particles in the required system form are calculated.
* If the initial $`Q^2`$ has not a fixed value but is instead simulated according to some probability distribution (for example, the Born cross section) then the cross sections have to be stored for reweighting. The $`Q^2`$ distribution is simulated over the Born cross section, and realistic observed $`Q^2`$ distribution is calculated as sum of weights, they are ratios of the total and Born cross sections.

We note that separation of the bremsstrahlung process into the radiative and nonradiative parts necessarily requires to introduce an additional separating parameter. Usually this quantity is associated with photon energy resolution in the detector. In our program the minimal inelasticity $`v_{min}`$ is chosen as this parameter.

## 2 Non-radiated and radiated contributions

The radiative events are generated in accordance with their contribution to the observed total cross section, that can be separated into non–radiative and radiative parts:

$$\sigma _{obs}=\sigma ^{nonrad}(v_{min})+\sigma ^{rad}(v_{min}).$$
(3)
Here and later we define $`\sigma d\sigma /dQ^2d\varphi `$. The first part includes the Born cross section $`\sigma _0`$ (Fig.1 (a)), loop effects (Fig.1 (b,c)) contributions as well as soft photon radiation, which is restricted by the inelasticity value $`v<v_{min}`$:

$$\sigma ^{nonrad}(v_{min})=\sigma _0e^{\delta _{inf}}(1+\delta _{VR}+\delta _{vac}+\delta ^{add}(v_{min}))+\sigma _R^{add}(v_{min}).$$
(4)
The explicit expressions for $`\sigma _0`$, $`\delta _{inf}`$, $`\delta _{VR}`$ can be found in . The dependence of the non-radiative part of the cross section on $`v_{min}`$ appears in the following way:

$$\delta ^{add}(v_{min})=\frac{2\alpha }{\pi }\mathrm{log}\left(\frac{v_{max}}{v_{min}}\right)\left[1\mathrm{log}\left(\frac{Q^2}{m^2}\right)\right],$$
(5)
$$\sigma _R^{add}(v_{min})=\frac{\alpha ^3}{2S^2}\underset{t_1^{}}{\overset{t_2^{}}{}}𝑑t\underset{i}{}\left[\theta _i\frac{_i(t)}{t^2}4\theta _i^B_{IR}\frac{_i(Q^2)}{Q^4}\right].$$
(6)
Here $`S=2k_1p_1`$, $`v_{max}=SQ^2M^2Q^2/S`$, $`_i(t)`$ are the squared combinations of the electric and magnetic elastic form factors. The quantities $`\theta _i`$ as well as $`\theta _i^B`$ are the analytical functions  of kinematic invariants for the radiative and Born subprocesses, respectively.

The radiative part of the total cross section can be presented as an integral over the real photonic phase space in the following way:

$$\sigma ^{rad}(v_{min})=\frac{\alpha ^3}{2S^2}\underset{t_1}{\overset{t_2}{}}𝑑t\underset{i}{}\frac{_i(t)}{t^2}\mathrm{\Theta }_i(t),$$
(7)
where

$$\mathrm{\Theta }_i(t)=\underset{v_1}{\overset{v_{max}}{}}𝑑v\theta _i(t,v)=\underset{v_1}{\overset{v_{max}}{}}𝑑v\underset{0}{\overset{2\pi }{}}𝑑\varphi _k\theta _i(t,v,\varphi _k),$$
(8)
and

$`v_1`$ $`=`$ $`max\{v_{min},{\displaystyle \frac{(tQ^2)(\sqrt{t}\sqrt{4M^2+t})}{2\sqrt{t}}}\},`$ (9)
$`t_{1,2}`$ $`=`$ $`{\displaystyle \frac{2M^2Q^2+v_{max}(Q^2+v_{max}\sqrt{(Q^2+v_{max})^2+4M^2Q^2})}{2(M^2+v_{max})}},`$ (10)
$`t_{1,2}^{}`$ $`=`$ $`t_{1,2}(v_{max}v_{min}).`$ (11)
The radiative photonic phase space $`t`$\- and $`v`$-variables are presented in Figure 2.

Notice that in (6) as well as in (7) the structure functions $`_i`$ depend on only one integration variable, namely, on $`t`$. Therefore to apply our program to different fits or models of nucleon form factors in (6) and (7), the numerical integration over the variable $`t`$ is used. On the other hand, to speed up the process of event generation, integration of the expressions $`\theta _i`$ in (8) over the other photonic variables $`v`$ and $`\varphi _k`$ has been performed analytically because the quantity $`\theta _i(t,v,\varphi _k)`$ appears as a result of the calculation of the matrix element squares with the real photon emission and do not depend on any model for elastic electron–nucleon scattering.

## 3 Simulation of radiative events

To simulate $`v`$ and $`\varphi _k`$ in accordance with the analytically calculated distributions without any numerical integration over $`t`$, the first step is simulation of the variable $`t`$. Then taking into account the simulated value $`t`$, the distribution function for $`v`$ can be constructed. Finally, $`\varphi _k`$ is simulated according to the simulated values of $`t`$ and $`v`$. If the event generation is performed successfully, the event distributions over these three variables has to follow the corresponding probability densities:

$$\begin{array}{cc}\rho (t)=\frac{1}{N_t}\underset{i}{}_i(t)\mathrm{\Theta }_i(t)/t^2,\hfill & N_t=\underset{t_1}{\overset{t_2}{}}𝑑t\underset{i}{}_i(t)\mathrm{\Theta }_i(t)/t^2\hfill \\ \rho (v)=\frac{1}{N_v}\underset{i}{}_i(t_g)\theta _i(t_g,v),\hfill & N_v=\underset{i}{}_i(t_g)\mathrm{\Theta }_i(t_g)\hfill \\ \rho (\varphi _k)=\frac{1}{N_{\varphi _k}}\underset{i}{}_i(t_g)\theta _i(t_g,v_g,\varphi _k),\hfill & N_{\varphi _k}=\underset{i}{}_i(t_g)\theta _i(t_g,v_g)\hfill \end{array}$$
(12)
As an example, the simulated distributions of the kinematic variables $`t`$ and $`v`$ together with their probability densities for representative kinematics at CEBAF energies are compared in Figure 3.

In summary, we presented the first results of the FORTRAN code ELRADGEN, which is a Monte Carlo generator of elastic unpolarized $`ep`$–scattering that includes electromagnetic radiative events.


format: text

text:
# A backtracking survey propagation algorithm for 𝐾-satisfiability

## 1 Introduction

Recently a new heuristic algorithm for solving constrain-satisfaction problems has been proposed . The algorithm is based on survey propagation and it has been successfully applied to the K-sat problem  and to the coloring problem . The algorithm is well performing in the random instance case, where the factor graph  associated to an instance of the problem is locally cycle-less (i.e. it has the structure of a tree).

The properties of this survey propagation algorithm have been studied in details in the case of the random 3-sat problem  in the limit of very large number ($`N`$) of literals . Here the relevant parameter is

$$\alpha =\frac{M}{N},$$
(1)
where $`M`$ is the total number of clauses.

The $`01`$ law  should be valid: for $`\alpha <\alpha _c`$ all random systems (with probability one when $`N`$ goes to infinity) are satisfiable. On the other hand, for $`\alpha >\alpha _c`$ no large random system (i.e. with probability one when $`N\mathrm{}`$) is satisfiable. An heuristic analysis  suggests that $`\alpha _c=\alpha ^{}4.267`$ where $`\alpha ^{}`$ can be analytically computed using the survey-propagation equations defined later.

Numerical experiments  shows that for very large $`N`$, e.g. $`N=10^6`$, the survey propagation algorithm works up to $`\alpha _S4.252`$, that is definite less, but very near to $`\alpha _c`$ with a time that is approximatively linear in $`N`$. In other words in the region $`\alpha _S<\alpha <\alpha _c`$ the problem is satisfiable for large $`N`$, but we are unable to find a solution to the problem using the survey decimation algorithm.

The aim of this note is to present a backtracking version of the survey propagation algorithm that is somewhat slower, but it able to solve problems where the original survey propagation algorithm does not work. In section II we present some general considerations on the backtracking strategies in the context of survey propagation. In section III we recall the basic definitions and results of the survey propagation approach for random instances. Finally in the last section we present the details of a simple implementation of backtracking.

## 2 The basic approach

The basic idea behind survey propagation is rather simple. In a satisfiability problem there are $`L`$ boolean variable (or literals) $`\sigma (i)`$ ($`i=1,L`$) that may be true or false. We have to assign a truth value to these literals in such a way that a given formula (i.e. the sat-formula) is satisfied.

If we assign a value to one literal (e.g. the one indexed by $`k`$) we remain with a problem with $`L1`$ literals. If we make all the correct assignments, we start from a problem with $`N`$ literals and recursively we reduce it to a null problem. Of course making the correct assignments is a non-trivial task. An heuristic approach is the following.

Let us consider a given problem with $`L`$ variables and let us call $`𝒩`$ is the total number of different truth value assignments that make the sat formula satisfied. In a similar vein we could define $`𝒩(i)_T`$ and $`𝒩(i)_F`$ the total number of different truth value assignment of the remaining $`L1`$ literals that make the sat-formula satisfied after that if we assign respectively a true value or false value to the $`i^{th}`$ literal. We can thus define

$$P_T(i)=\frac{𝒩(i)_T}{𝒩},P_F(i)=\frac{𝒩(i)_F}{𝒩}.$$
(2)
It may be useful to define also

$$P(i)=\mathrm{max}(P_T(i),P_F(i)).$$
(3)
Let us suppose to we know the quantities $`P_T(i)`$ and $`P_F(i)`$. If we recursively decrease the size of the problem from $`N`$ to 0 by assigning at each step the literal having the maximal value of $`P(i)`$, we are going to find a correct assignment, if it exist. While it is out of question to compute exactly these two functions, it is possible to get some reasonable estimates of them. The survey equations (that will be described in the next section) provides an heuristic and efficient way to estimate these two functions.

Generally speaking, giving an heuristic for estimating the functions $`P_T(i)`$ and $`P_F(i)`$ it is possible to use this information to recursively decrease the size of the systems. This is the essence of the survey propagation method. Indeed in many cases the survey propagation equations tell us an estimate of the probabilities $`P_T(i)`$ and $`P_F(i)`$ that is the best that can be obtained by a local analysis (i.e. an analysis such that for given $`i`$ $`P_T(i)`$ can be computed with a fixed error in a time that is of order 1).

Unfortunately a local analysis is not exact and the survey estimations obtained in this way may be not fully correct: it is possible that choices done at the initial stage turn out to be wrong. This may be checked by computing the increase in the number of configurations that satisfy the sat formula if at the step $`L`$ we leave free the $`k^{th}`$ literal that was blocked at the step $`M`$. This quantity will be denoted by $`1/I(k)`$. If the literal $`k`$ has been frozen at the previous step (i.e. $`M=L1`$), we have that

$$P(k)=I(k).$$
(4)
If after freezing the other literals, when the difference $`LM`$ becomes large, $`I(k)`$ decreases, this is a signal that a wrong choice has been done at the step$`M`$ and we would like to backtrack this choice. An algorithm, like survey propagation, that is able to compute the quantities $`P(k)`$ is usually able to compute also the quantities $`I(k)`$.

These observations suggest a more general algorithm than pure decimation. In our problem we may have $`NL`$ literals frozen (i.e their value has been assigned) and $`L`$ literals free. Our aim is to start from $`L=N`$ and to arrive to $`L=0`$, with the highest possible value of $`𝒩`$ at the final destination. We can consider two kind of moves: the usual decimation moves ($`LL1`$) and the backtracking moves $`LL+1`$. It seems natural in this approach to select for a decimation move the variable with largest $`P(i)`$, where $`i`$ labels a variable that is not assigned, and to select for a backtracking move the variable with smallest $`I(k)`$. When

$$P_M\underset{i}{\mathrm{max}}P(i)>I_m\underset{k}{\mathrm{min}}I(k),$$
(5)
a combined move, where we assign a value to the variable at $`i_M`$ (i.e. $`P(i_M)=P_M`$) and we deassign a value to the variable at $`k_m`$ ($`I(k_m)=I_m`$), leads to a new configuration with the same value of $`L`$ but with an higher value of $`𝒩`$ (at least in the generic case where the points $`i_M`$ and $`k_m`$ are far away on the factor graph and the corresponding literals are not correlated).

## 3 A survey of the survey propagation equations

The survey propagation equations are described in details in the literature . Here we present only the basic ideas.

For a given problem it is convenient to consider the set of all the solutions of the problem, i.e. of all the assignments that satisfy the sat formula. We are we are interested not only to know the number of all the solutions, but also to know other properties of this set.

At this end it is convenient to say that two solutions are adjacent if their Hamming distance is less than $`ϵN`$, where $`ϵ`$ is an appropriate small number.

It has been suggested that in the limit of large $`N`$:

1. In the interval $`\alpha <\alpha _d`$ ($`\alpha _d3.86`$ for the 3-sat ) the set of all solutions is connected, i.e. there is a path of mutually adjacent solutions that joins any two solutions of the set.
2. In the interval $`\alpha _d<\alpha <\alpha _c`$ ($`\alpha _c3.92`$ for the 3-sat ) the set of all the solutions breaks in an large number of different disconnected regions that are called with many different names in the literature  (states, valleys, clusters, lumps…). Roughly speaking the set of all the solutions can be naturally decomposed into clusters of proximate solutions, while solutions belonging to different clusters are not close.
3. Only in the interval $`\alpha _b<\alpha <\alpha _c`$ ($`\alpha _b3.92`$ for the 3-sat; in other models $`\alpha _b`$ and $`\alpha _d`$ may coincide) there are literals $`\sigma `$ that take the same value in all the legal configurations of a region. We say that these literals form the backbone of a region. It is important to realize that a given literal may simultaneously belong to the backbone of one region and not belong to the backbone of an other region.
4. Only the interval $`\alpha _m<\alpha <\alpha _c`$ ($`\alpha _m4.153`$ for the 3-sat model) the simple survey approach is correct . For $`\alpha _b<\alpha <\alpha _m`$ the simple survey approach is likely to be an approximation and a more refined analysis is needed in order to get exact results.

The last region is the hardest from the computational point of view. Our analysis is done in this region and we will concentrate our efforts on the properties of the backbone. Only in the third region the survey propagation equations have a non-trivial solution, i.e. with surveys that do not have $`s_I(i)`$ equal to one for all $`i`$. The appearance of a non-trivial solution of the survey equation marks the entrance in the third region.

The main quantity we will study is $``$, i.e. the number of regions with different backbones. As far as for large $`N`$ there is an exponential large number of configuration in each region, $``$ is much smaller than $`𝒩`$. One usually introduce the entropy  $`S`$ and the complexity  $`\mathrm{\Sigma }`$ defined as

$$𝒩=\mathrm{exp}(S),=\mathrm{exp}(\mathrm{\Sigma }).$$
(6)
In the rest of this note we will consider $``$ and not $`𝒩`$ for two reason:

* It is much simpler to compute $``$ that $`𝒩`$.
* It can be argued that the quantity $``$ is more interesting than $`𝒩`$ because in the limit $`\alpha \alpha _c^{}`$, $``$ is a quantity of order 1, while $`𝒩`$ is exponentially large (nolo sumere, nondum acerbum est). In other words the entropy density ($`S/N`$) does not vanish in the limit, while the complexity density ($`\mathrm{\Sigma }/N`$) vanishes.

In a similar way the functions $`P_T(i)`$ and $`P_F(i)`$ will be defined as the ratio of the appropriate quantities $``$ (we must substitute $`𝒩`$ with $``$ in eq. (2)). One hopes that this does not make a great difference.

In the survey approach our aim is to compute $``$ and the surveys $`\stackrel{}{s}(i)`$ where $`\stackrel{}{s}`$ is a three component vector $`\stackrel{}{s}=\{s_T,s_I,s_F\}`$ (the sum of the components is one). The quantity $`s_T(i)`$ is the fraction of backbones that contain the literal at $`i`$ with value true, $`s_F(i)`$ is the fraction of backbones that contains the literal at $`i`$ with value false $`s_I(i)`$ is the fraction of backbones that do not contain the literal at $`i`$. If for any $`i`$ $`s_I(i)=1`$ backbones are not present and the problem is easy.

From the point of view or our paper the surveys are extremely interesting because with very good approximation we have that

$$P_T(i)=1s_F(i),P_F(i)=1s_T(i)$$
(7)
More precisely it can be shown that $`P_T(i)=1s_F(i)+O((1s_F(i))^2)`$ that is good enough for us because in most of the cases we are interested to those values of $`i`$ such that $`P_T`$ (or $`P_F`$) is very near to 1.

Unfortunately it is not possible to write local equations directly for the surveys, but we have to introduce the cavity surveys , i.e. the value that the survey at $`i`$ would have if we remove from the problem the clause $`c`$ (under the conditions that the literal at $`i`$ enters in the clause $`c`$). On large random problem the cavity surveys $`\stackrel{}{s}(i,c)`$ satisfy simple equations that are called the survey propagation equations .

The solution of these equations may be used to compute both $``$ and the surveys $`\stackrel{}{s}(i)`$.

## 4 The survey algorithm

### 4.1 Survey decimation algorithm

The survey decimation algorithm has been proposed  for finding the solution of the random K-satisfiability problem .

We start by solving the survey propagation equation and to estimate the function $`P(i)`$. The main step in the decimation procedure consists in starting from a problem with $`L`$ variables and in arriving to a problem with $`L1`$ variables where $`\stackrel{}{s}(i)`$ is fixed to be true (or false), i.e. the literal at $`i`$ has an assigned value. We have already remarked that if $`P(i)`$ is near to 1, the second problem it is easier to solve: it has nearly the same number of solutions of the belief equations and one variable less. (We assume that the complexity can be reliable computed by solving the survey propagation equations).

The decimation algorithm proceeds as follows. We reduces by one the number of variables choosing the node $`i`$ in the appropriate way, e.g. by select the literal with maximal $`P(i)`$. In a slightly less efficient but faster version we select $`R`$ literals with maximal $`P`$, where $`R=fN`$ and $`f`$ is a small number. If we are not too near the critical point and $`f<.01`$ the results are only weakly dependent on $`f`$.

We recompute the solutions of the survey equations and we reduce again the number of variables; it is convenient to monitor the complexity. At the end of the day different things may happen:

* We arrive to a negative complexity (in this case the reduced problem should have no solutions and we are lost).
* If in spite of this warning we go on with the algorithm, we usually find that the reduced problem contains a contradiction and in this case we are definitively lost. Sometimes we can find a contradiction also for positive complexity, but this happens always in the region of very small complexity density.
* Other troubles arise when the survey equations do not have anymore an unique solution and the usual iterative procedure for finding the solution of the survey equations does not converge. It is possible that one could find some way to escape from the troubles, but this case has never been studied in details.
* After a sufficient large number of iterations the survey equation do not have any more a non-trivial solution and the standard iteration algorithm converges to the trivial solution ($`s_T(i)=s_F(i)=0`$, $`s_I(i)=1`$). If this happens the reduced problem is easy and it can be solved using other algorithms.
* An other possible happy end is that the survey equations converge to a zero entropy solution where (depending on the site) or $`s_T(i)=1`$ or $`s_F(i)=1`$. In this case we do not need other algorithms to solve the reduced problem.

A careful analysis  of the results for large, but finite $`N=O(10^6)`$, shows that the algorithm works in the limit of infinite $`N`$ up to $`\alpha _S(f)`$ where in the limit of $`f0`$ we have $`\alpha _S4.252`$, that is definite less, but very near to $`\alpha _c`$.

It is interesting to note that for $`\alpha <\alpha _S`$ the survey decimation algorithm takes a time that is polynomial in $`N`$ and using a smart implementation the time is nearly linear in $`N`$. The survey algorithm is an incomplete search procedure that may not be able to find a solution to a satisfiable instance. This actually happens with a non-negligible probability e.g. for sizes of the order of a few thousands of variables also when $`\alpha <\alpha _S`$, however in the limit $`N\mathrm{}`$, it should work as soon $`\alpha <\alpha _S`$. In the interval $`\alpha _S<\alpha <\alpha _c`$ the decimation procedure leads to a regime of negative complexity so that the algorithm does not work.

Unfortunately at the present moment the value of $`\alpha _S`$ can be obtained only analyzing how the argument works on a finite sample and we are unable to compute analytically $`\alpha _S`$: it is a rather difficult task to understand in details why for $`\alpha _S`$ it so near to $`\alpha _c`$.

### 4.2 Backtracking surveys

In the backtracking version of the survey propagation algorithm we follow the general strategy outlined in the introduction.

In order to solve a sat problem with $`N`$ literals we consider auxiliary problems where $`L`$ literals are free and $`NL`$ literals have a value assigned. Obviously for given $`L`$ we have $`2^{NL}`$ of such problems.

Our aim is to start from $`L=0`$ and to arrive to to $`L=N`$, with an assignment to the literals such that the sat formula is satisfied. An heuristic argument suggests that a reasonable strategy is to keep the highest possible value of $`𝒩`$ at each intermediate step.

We can consider two kind of moves: the usual decimation moves ($`LL+1`$ ) and the backtracking moves ($`LL1`$). In the simplest version of the the algorithm we select for applying the decimation rule the literal with largest $`P(i)`$, where $`i`$ labels a literal that is not assigned, and we select for a backtracking move the literal with smallest $`I(k)`$. The computation of $`P(i)`$ is standard and $`I(k)`$ can be just computed in the same way as function of the surveys of the nodes adjacent to the $`k^{th}`$ node.

Once we know how to select the literals to decimate and the literals to backtrack, the most important decision is if we do a decimation or a backtracking move.

The simplest possibility is do decimations moves and backtracking moves at a fixed ratio. This algorithm depends on the parameter $`r`$, i.e. the ratio between the number of backtracking moves and the total number of moves (obviously we must have $`r<0.5`$). This choice has the advantage that the total number of moves is still bounded by a number proportional to $`N`$ (more precisely by $`N/(0.5r)`$). If it works, it is very efficient.

We have done a few numerical experiments and we have found that the introduction of the backtracking moves with the simple strategy outline before definitely increases the range of applicability of the survey equation algorithm.

In fig. 1. we show the results for the $`\mathrm{\Sigma }`$ as function of $`L`$ for a 3-sat problem with $`N=10^5`$ and $`\alpha =4.25`$ for the algorithm with $`f=10^3`$ with three values of $`r`$ (0, .25 and .4). In fig.2 we plot the quantity

$$F(L)=\frac{\underset{i\text{non-frozen}}{}|s_T(i)s_T(i)|}{L}$$
(8)
where the sum is done only on the $`L`$ non-frozen nodes. The vanishing of $`F`$ implies that the solution of the survey equation is trivial.

For $`r=0`$, i.e. the standard survey algorithm, the complexity arrives to zero with $`F0`$ and the survey algorithm fails to solve the problem. The further evolution (not shown) gives a negative complexity and after some more iterations we find a contradiction and the algorithms must stops (as expected). The fact that the survey algorithm fails for a value of $`\alpha `$ very near, but smaller than $`\alpha _S`$ is due to fluctuations in the value of $`\alpha _S`$ that are present for the finite value of $`N`$ and vanish when $`N\mathrm{}`$.

On the contrary both for $`r=.25`$ and $`r=.4`$ we arrive to a situation where the complexity and $`F`$ simultaneously go to zero, i.e. to the trivial solution of the survey equations. The backtracking algorithm performs better that the non-backtracking one.

The backtracking algorithms is able to find an easy reduced problem also in the case where the non-backtracking algorithm fails. It is interesting to note that increasing $`r`$, the complexity at fixed $`L`$ increases signaling a more efficient algorithms. Indeed there are other cases (e.g. for slightly larger $`\alpha `$) where the $`r=.25`$ algorithm fails and the $`r=.4`$ is successful.

The same phenomenon happens in a consistent way at different values of $`f`$ and $`N`$: also this simple version of the backtracking algorithm is more performing that the non-backtracking one for difficult problems.

More carefully and systematic investigations are needed in order find out the dependence of $`\alpha _A`$ on the backtracking ratio $`r`$ and to study how near we can go to the critical point ($`\alpha _c`$) using this method.

It is rather likely that the backtracking method in its present form does not work too near to $`\alpha _c`$. Indeed the first decrease of the complexity (for $`L`$ near $`N`$) is weakly sensitive to the backtracking; it is likely using a survey propagation algorithm we must enter in a region of negative complexity if we are sufficiently near to $`\alpha _c`$. In principle the algorithm is able to recover and to increase the complexity when we decrease $`L`$: this exactly what happens in a few cases for $`N=O(10^3)`$.

Unfortunately, when we are too near to $`\alpha _c`$ and we are deep in the region of negative complexity, the iterative procedure for finding the solution of the survey propagation equations does not converge anymore and the algorithm stop to work. The non-convergence of the iterative procedure for finding the solution of the survey propagation equations is a sign of the existence of multiple solutions of the survey propagation equations. It is likely that the next step to do is to implement a version of the algorithm that bypass this difficulty. It would also be very interesting to do a more careful analysis of different strategies in doing backtracking steps. The relative efficiency of backtracking in non-random instances is also a problem that one should investigate carefully.

## Acknowledgements

It is a pleasure for me to thank David Forney for a stimulating discussion.


text:
# New exact treatment of the perturbed Coulomb interactions

## Abstract

A novel method for the exact solvability of quantum systems is discussed and used to obtain closed analytical expressions in arbitrary dimensions for the exact solutions of the hydrogenic atom in the external potential $`\mathrm{\Delta }V(r)=br+cr^2`$, which is based on the recently introduced supersymmetric perturbation theory.

One of the challenging problems in nonrelativistic quantum mechanics is to find exact solutions to the Schrödinger equation for potentials that are used in different fields of physics. In particular, the perturbed Coulomb potentials represent simplified models of many situations found in atomic, molecular, condensed matter and particle physics. There has been much interest in obtaining analytical solutions of such potentials in arbitrary dimensions. These problems have been studied for years and a general solution has not yet been found.

Such class of potentials,

$$V(r)=a/r+br+cr^2,$$
(1)
are possible candidates for the quarkonium potential as has been indicated by the quarkonium spectroscopy. In the special case of $`c=0`$ and $`b>0`$ such potentials reduce to the well known charmonium potential. Apart from its relevance in heavy quarkonium spectroscopy, this class of potentials with $`c=0`$ has important applications in atomic physics. The Stark effect in a hydrogen atom in one dimension is given exactly by the charmonium like potential ($`b`$ being the electric field parameter). The more general class of these potentials with $`c>0`$ is also relevant in atomic physics. This could be interpreted as the potential seen by an electron of an atom exposed to a suitable admixture of electric and magnetic fields. In addition, nuclei in the presence of an electron background form a system which is important for condensed matter physics and for laboratory and stellar plasmas. The potential between two nuclei embedded in such a plasma is approximately Coulomb plus harmonic oscillator, which corresponds to $`b=0`$ in (1).

As the exact form of such interactions are being unknown to a great extent, it is thus desirable to study the general analytical properties of a large class of potentials in (1). In connection with this, the analyticity of the energy levels for these kind of potentials was investigated rigorously by many authors using different theories - in relation to their potential applications in spectroscopic problems.

In this letter, we introduce an alternative, simple formalism for an algebraic solution of the Schrödinger equation with the perturbed Coulomb potential and find exact solutions in $`N`$-dimesional space. The new formalism is based on the supersymmetric quantum mechanics and uses the spirit of perturbation theory, in which there is no room for the drawbacks encountered in the calculations with the perturbation theories available in the literature. The work introduced in this letter also clarifies, within the powerful frame of the present formalism, that the potential in (1) is in fact an exactly solvable shape invariant potential unlike the claim in a recent work  where the authors stated that the supersymmetric quantum mechanics yields exact solutions for a single state only for the perturbed Coulomb potential in (1) which was treated as a quasi-exactly solvable potential.

Let us now introduce the formalism. The reduced radial wave Schrödinger equation for a spherically symmetric potential in $`N`$-dimensional space reads

$$\frac{\mathrm{}^2}{2m}\frac{\mathrm{\Psi }_n^{\prime \prime }(r)}{\mathrm{\Psi }_n(r)}=V(r)E_n,V(r)=V_{ES}(r)+\mathrm{\Delta }V(r),V_{ES}(r)=V_0(r)+\frac{\mathrm{\Lambda }\left(\mathrm{\Lambda }+1\right)\mathrm{}^2}{2mr^2},$$
(2)
where $`V_0`$ is one of the exactly solvable potential, such as the Coulomb and harmonic oscillator potentials used through the present work, $`n=0,1,2,\mathrm{}`$ being the radial quantum number, and $`\mathrm{\Lambda }=(M3)/2`$ with $`M=N+2\mathrm{}`$. We see that the radial Schrödinger equation in $`N`$ dimensions has the same form as the three-dimensional one. Consequently, given that the potential has the same form in any dimension, the solution in three dimensions can be used to obtain the solution in arbitrary dimensions simply by the use of the substitution $`\mathrm{}\mathrm{\Lambda }`$. In the above equation, $`\mathrm{\Delta }V`$ stands for the small perturbation $`\left(=br+cr^2\right)`$ assuming the potential parameter $`a`$ is large when compared to other coupling parameters $`b`$ and $`c`$ in (1) with the consideration of perturbed Coulomb interactions while for the case of large $`c`$ and relatively small $`a`$ and $`b`$ the perturbing potential involves $`\mathrm{\Delta }V=a/r+br`$ leading to a perturbed harmonic oscillator potential. This significant feature will be used later in this work to test our results.

Using the spirit of the standard perturbation theory, we put forward here a more general formalism to treat quantum mechanical perturbation problems efficiently in a simple framework, which has been already discussed in detail through the recent works  involving different applications of the model. Start with

$$\mathrm{\Psi }_n(r)=\chi _n(r)\varphi _n(r),$$
(3)
in which $`\chi _n`$ is the known normalized eigenfunction for the exactly solvable potential in (1) including the barrier term and $`\varphi _n`$ is a moderating function due to the perturbing potential. Substituting (3) into (2) yields

$$\frac{\mathrm{}^2}{2m}\left(\frac{\chi _n^{\prime \prime }}{\chi _n}+\frac{\varphi _n^{\prime \prime }}{\varphi _n}+2\frac{\chi _n^{}}{\chi _n}\frac{\varphi _n^{}}{\varphi _n}\right)=VE_n.$$
(4)
With the new definitions,

$$W_n=\frac{\mathrm{}}{\sqrt{2m}}\frac{\chi _n^{}}{\chi _n},\mathrm{\Delta }W_n=\frac{\mathrm{}}{\sqrt{2m}}\frac{\varphi _n^{}}{\varphi _n}$$
(5)
one arrives at

$$\frac{\mathrm{}^2}{2m}\frac{\chi _n^{\prime \prime }}{\chi _n}=W_n^2\frac{\mathrm{}}{\sqrt{2m}}W_n^{}=V_{ES}ϵ_n,$$
(6)
where $`ϵ_n`$ is the eigenvalue for the exactly solvable potential of interest, and

$$\frac{\mathrm{}^2}{2m}\left(\frac{\varphi _n^{\prime \prime }}{\varphi _n}+2\frac{\chi _n^{}}{\chi _n}\frac{\varphi _n^{}}{\varphi _n}\right)=\mathrm{\Delta }W_n^2\frac{\mathrm{}}{\sqrt{2m}}\mathrm{\Delta }W_n^{}+2W_n\mathrm{\Delta }W_n=\mathrm{\Delta }V\mathrm{\Delta }ϵ_n,$$
(7)
in which $`\mathrm{\Delta }ϵ_n`$ is the correction term to the energy due to $`\mathrm{\Delta }V`$, and $`E_n=ϵ_n+\mathrm{\Delta }ϵ_n`$. Subsequently, Eq. (4) reduces to

$$\left(W_n+\mathrm{\Delta }W_n\right)^2\frac{\mathrm{}}{\sqrt{2m}}\left(W_n+\mathrm{\Delta }W_n\right)^{}=VE_n.$$
(8)
As the complete spectrum and wave functions are known in the analytical form for the exactly solvable potentials appeared in (6), one needs here to solve (7) to obtain in a closed form the corrections to both, energy and wave function.

Proceeding with the perturbed Coulomb potential in arbitrary dimensions,

$$V(r)=\left[\frac{a}{r}+\frac{\mathrm{\Lambda }\left(\mathrm{\Lambda }+1\right)}{2mr^2}\right]+br+cr^2,$$
(9)
and set the superpotential, for the ground solutions,

$$W_{n=0}(r)=\sqrt{\frac{m}{2}}\frac{a}{\left(\mathrm{\Lambda }+1\right)\mathrm{}}\frac{\left(\mathrm{\Lambda }+1\right)\mathrm{}}{\sqrt{2m}r},$$
(10)
leading to the first part of the potential in the bracket, and from the literature, corresponding normalized wave function and energy are

$$\chi _{n=0}(r)=N_0r^{\mathrm{\Lambda }+1}\mathrm{exp}\left[\frac{ma}{\left(\mathrm{\Lambda }+1\right)\mathrm{}^2}r\right],ϵ_{n=0}=\frac{ma^2}{2\mathrm{}^2\left(\mathrm{\Lambda }+1\right)^2}.$$
(11)
For the perturbing potential, $`\mathrm{\Delta }V=br+cr^2`$, the physically acceptable unique choice is

$$\mathrm{\Delta }W_{n=0}(r)=\sqrt{c}r,$$
(12)
which satisfies Eq. (7) from where one readily sees that

$$\mathrm{\Delta }ϵ_{n=0}=\frac{M\left(M1\right)b\mathrm{}^2}{4ma},b=\frac{2a\sqrt{2mc}}{\left(M1\right)\mathrm{}}.$$
(13)
It is stressed that such solutions in general have constraints on the potential parameters as appeared in (13). These constraints differ for each eigenvalue, and hence various solutions do not correspond to the same potential and are not orthogonal.

From (12), the moderating function

$$\varphi _{n=0}(r)=\mathrm{exp}\left[\frac{\sqrt{2m}}{\mathrm{}}\stackrel{r}{}\mathrm{\Delta }W_{n=0}(r)\right]=\mathrm{exp}\left[\frac{b\left(M1\right)}{4a}r^2\right].$$
(14)
Therefore, the full wave function $`\mathrm{\Psi }_{n=0}`$ for the potential in (9) takes the form

$$\mathrm{\Psi }_{n=0}(r)=\chi _n(r)\varphi _{n=0}(r)=N_0r^{\left(\frac{M1}{2}\right)}\mathrm{exp}\left[\frac{2ma}{\left(M1\right)\mathrm{}^2}r\frac{b\left(M1\right)}{4a}r^2\right].$$
(15)
Finally, the exact ground state energy is

$$E_{n=0}=ϵ_{n=0}+\mathrm{\Delta }ϵ_{n=0}=\frac{2ma^2}{\mathrm{}^2\left(M1\right)^2}+\frac{M\left(M1\right)b\mathrm{}^2}{4ma}.$$
(16)
The results in Eqs. (15) and (16) are exact and agree with those in .

As noted earlier, the potential in (9) behaves also like a harmonic oscillator in case for large $`c`$ and relatively small $`a`$ and $`b`$. This enables one to check out explicitly the results obtained above and the reliability of the formalism introcuded. For this consideration, we set the superpotentials

$$W_{n=0}(r)=\sqrt{c}r\frac{\left(\mathrm{\Lambda }+1\right)\mathrm{}}{\sqrt{2m}r},\mathrm{\Delta }W_{n=0}(r)=\frac{b/2}{\sqrt{c}},$$
(17)
which yield, through the use of either Eqs. (6) and (7) together, or (8) alone,

$$ϵ_{n=0}=\frac{\mathrm{}\sqrt{c}}{\sqrt{2m}}\left(2\mathrm{\Lambda }+3\right),\mathrm{\Delta }ϵ_{n=0}=\frac{b^2}{4c},$$
(18)
which are exactly equivalent to those found for the perturbed Coulomb case in (16). Similarly, using (17) one readily arrives at (15).

For the generalization, we further make clear that the present technique is also applicable to excited states, for which one needs to use the shape invariance property and the relation between supersymmetric partners

$$V^+(r,\alpha _0)=V^{}(r,\alpha _1)+R(\alpha _1),$$
(19)
where $`\alpha _0=\mathrm{\Lambda }`$ and $`\alpha _1=f\left(\alpha _0\right)=\mathrm{\Lambda }+1`$ are the position independent parameters while $`V^\pm `$ are the supersymmetric partners

$$\left(W_{n=0}+\mathrm{\Delta }W_{n=0}\right)^2\frac{\mathrm{}}{\sqrt{2m}}\left(W_{n=0}+\mathrm{\Delta }W_{n=0}\right)^{}=V^{}E_{n=0}^{},$$
(20)
which depend upon superpotentials in (10) and (12) for the case $`a>b,c`$. Bearing in mind that the replacement of $`\mathrm{\Lambda }`$ with $`\mathrm{\Lambda }+1`$ in $`V^{}(r,\mathrm{\Lambda })`$ leads to $`V^{}(r,\mathrm{\Lambda }+1)`$ and keeping the potential parameters $`b`$ and $`c`$ fixed as in  while allowing the change in $`a`$, see Eq. 13, we obtain by the straightforward calculations

$`E_n`$ $`=`$ $`E_n^{}{\displaystyle \frac{b^2}{4c}}+{\displaystyle \frac{\mathrm{}\sqrt{c}}{\sqrt{2m}}}\left(2\mathrm{\Lambda }+3\right),`$
$`E_n^{}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{n}{}}}R(\alpha _k)={\displaystyle \frac{\mathrm{}\sqrt{c}}{\sqrt{2m}}}\left[\left(2\mathrm{\Lambda }+3\right)2\left(n+\mathrm{\Lambda }\right)+3\right],`$
$`E_n`$ $`=`$ $`{\displaystyle \frac{b^2}{4c}}+{\displaystyle \frac{\mathrm{}\sqrt{c}}{\sqrt{2m}}}\left[2\left(n+\mathrm{\Lambda }\right)+3\right],n=0,1,2,\mathrm{},`$ (21)
which agrees with  in $`N`$-dimensional space and with  in three dimensions. Furthermore, one can easily construct the bound $`n^{th}`$ state wave functions from the ground state wave function using the supersymmetric definition

$$\mathrm{\Psi }_{n+1}^{}(r,\alpha _0)A^+(r,\alpha _0)\mathrm{\Psi }_n^{}(r,\alpha _1),A^+(r,\alpha _0)=\frac{\mathrm{}}{\sqrt{2m}}\left[\frac{d^2}{dr^2}+\frac{\mathrm{\Psi }_{n=0}^{}(r,\alpha _0)}{\mathrm{\Psi }_{n=0}(r,\alpha _0)}\right].$$
(22)
We conclude with two remarks. First, the present formalism can be generalized to all the polynomial forces

$$V_m(r)=Ar^{2m}+Br^{2m1}+\mathrm{}+F/r+G/r^2$$
(23)
as an alternative but simple treatment to the other works , and the references therein. With the consideration of the perturbed Coulomb problem we have clarifed that the explicit solution of the related Schrödinger equation with the potential family in (23) remains feasible in an almost complete parallel with their $`m=1`$ predecessors. Along this line, the works are in progress.

The second remark we wish to make is that apart from the inherent interest one has in the existence of exact solutions, the results reported here are likely to be useful in perturbation calculations for the excited state energies and wave functions, particularly if the method of other perturbation theories appeared in the literature for evaluating second- and higher-order corrections can be extended to cover such cases. In this context, we believe that the simple form of our new approach to such problems opens a new direction of development towards many practical applications yet to be constructed and appreciated.


text:
# 1 Some Concepts of Nonstandard Analysis

## 1 Some Concepts of Nonstandard Analysis

The nonstandard analysis theory of turbulence will be presented in this paper. The nonstandard analysis means the mathematic fundament of the nonstandard picture of turbulence. Therefore, some concepts of nonstandard analysis should be introduced in the first place.

On the basis of mathematic logic, A.Robinson  proved in the sixties of the last century that a real number system $`R`$ can be expanded into the hyperreal number system $`R^{}`$. Besides all real numbers (i.e., standard numbers), $`R^{}`$ contains also hyperreal numbers (i.e., nonstandard numbers). The infinitesimal $`\epsilon `$ and infinite $`L`$ are elementary nonstandard numbers. There are other nonstandard numbers, for example, $`\xi \pm \epsilon ,\xi \epsilon ,\xi L,\xi \pm L,\epsilon \epsilon ,LL,`$ etc. ($`\xi `$ stands for any real number).

It is known that real numbers one-to-one correspond to the points on a real line, but on the real line there are no such points which correspond to nonstandard numbers. In other words nonstandard numbers can not be indicated by the points on a real line. The concept of monad needs to be presented to indicate the nonstandard numbers. For every real number $`\xi `$, there exists a monad that is composed of the numbers infinitely close to $`\xi `$ and $`\xi `$ belongs also to the monad. In fact, there exist infinite monads containing the real number $`\xi `$. The size of all these monads is infinitesimal $`\epsilon `$. The numbers in this monad, except $`\xi `$, are all nonstandard numbers of which $`\xi `$ is the standard part.

So any point on a real line is expanded into a monad. The real number, represented by the point on the real line, is the standard part of the nonstandard numbers contained in the monad. And a real line, formed by the points which indicate only real numbers,is expanded into a hyperreal line composed of monads. The set of the standard parts of all these monads composes the real line mentioned above. Every monad on a hyperreal line is also called a point (i.e.,standard point). Clearly every such standard point is not the absolute point which has the length of absolute zero, but the monad whose length is infinitesimal $`\epsilon `$.

Every monad must contain only one real number. If a monad contains the real number $`\xi `$, we call the monad as $`\xi `$-monad.

The physical world is in a hierarchical structure. This hierarchical structure can be described by the concepts of nonstandard mathematics. For example, there are three levels: level $`\alpha `$ (with infinitesimal $`\epsilon `$, infinite $`\epsilon ^1`$), level $`\beta `$ ( with infinitesimal $`\epsilon ^3`$, infinite $`\epsilon `$) and level $`\gamma `$ ( with infinitesimal $`\epsilon ^1`$, infinite $`\epsilon ^3`$). Any monad on a hyperreal line in level $`\alpha `$ is also a whole hyperreal line in level $`\beta `$. However, the whole hyperreal line in level $`\alpha `$ is only a monad of the hyperreal axis in higher level $`\gamma `$. From the viewpoint of hierarchical structure, the meanings of infinite and infinitesimal are not fixed. It is understandable in physics. Physicists, in fact, can even determine how large infinite and infinitesimal are in a specific physical problem, and often take a certain number as infinite or infinitesimal. In physics any number (not absolute zero) is taken as infinite, infinitesimal or finite in various specific condition.

The coordinate systems in two levels will be set up for the convenience of description. Let space coordinates be $`(x_1,x_2,x_3)`$, time $`t`$ in level $`\alpha `$ and space $`(x_1^{},x_2^{},x_3^{})`$, time $`t^{}`$ in level $`\beta `$, the range of values of $`(x_1^{},x_2^{},x_3^{},t^{})`$ be $`[0,L_1),[0,L_2),[0,L_3)`$ and $`[0,T)`$ respectively. Here $`L_1,L_2,L_3`$ and $`T`$ are the infinites of space and time, respectively, in level $`\alpha `$. For any time-space point $`(\stackrel{}{x},t)`$ in level $`\alpha `$, the monad of this point is expressed by $`(\stackrel{}{x},t)`$ and the dimensions of the monad are $`(\epsilon _1,\epsilon _2,\epsilon _3,\epsilon _t)`$. Here $`\epsilon _1,\epsilon _2,\epsilon _3,\epsilon _t`$ are infinitesimals of space and time, respectively, in level $`\alpha `$. Moreover, assume that:

$$\epsilon _1L_1=1,\epsilon _2L_2=1,\epsilon _3L_3=1,\epsilon _tT=1$$
(1)
The space-distance between any two points

$$\mathrm{\Delta }l=\sqrt{\left(\mathrm{\Delta }x_1\right)^2+\left(\mathrm{\Delta }x_2\right)^2+\left(\mathrm{\Delta }x_3\right)^3}$$

in level $`\alpha `$, and

$$\mathrm{\Delta }l^{}=\sqrt{\left(\mathrm{\Delta }x_1^{}\right)^2+\left(\mathrm{\Delta }x_2^{}\right)^2+\left(\mathrm{\Delta }x_3^{}\right)^2}$$

in level $`\beta `$. The actual length (i.e., the length observed from the angle of level $`\alpha `$) of $`\mathrm{\Delta }l^{}`$ is:

$$\sqrt{\left(\mathrm{\Delta }x_1^{}\epsilon _1^2\right)^2+\left(\mathrm{\Delta }x_2^{}\epsilon _2^2\right)^2+\left(\mathrm{\Delta }x_3^{}\epsilon _3^2\right)^2}$$
(2)
If $`\epsilon _1=\epsilon _2=\epsilon _3=\epsilon `$ the actual length of $`\mathrm{\Delta }l^{}`$ is:

$$\mathrm{\Delta }l^{}\epsilon ^2$$
(3)
Similarly, the actual value (i.e., the value observed from level $`t`$) of $`\mathrm{\Delta }t^{}`$ is:

$$\mathrm{\Delta }t^{}\epsilon _t^2$$
(4)
And the actual lengths of the whole space and time axes in level $`\beta `$ are:

$$L_1\epsilon _1^2=\epsilon _1,L_2\epsilon _2^2=\epsilon _2,L_3\epsilon _3^2=\epsilon _3,T\epsilon _t^2=\epsilon _t$$
(5)
They are just the dimensions of a monad in level $`\alpha `$.

Now the limit of tending to zero ($`\mathrm{\Delta }x_i0`$) in level $`\alpha `$ is not tending to absolute zero but infinitesmal $`\epsilon `$. Here $`\mathrm{\Delta }x_i`$ can not go into the inner of the zero monad. In the inner of the zero monad in level $`\alpha `$ exists $`\mathrm{\Delta }x_i^{}`$ instead of $`\mathrm{\Delta }x_i`$. In other words, the process of tending to zero does not mean tending to absolute zero but the dimension of a monad in the same level. The standard part of any number $`\omega `$ is denoted by $`St\omega `$. To show clearly the meanings of $`\mathrm{\Delta }x_i0`$ and $`\mathrm{\Delta }t0`$ mentioned above, we do not use the notation $`\mathrm{\Delta }x_i0`$ and $`\mathrm{\Delta }t0`$. Instead, we use

$$St\mathrm{\Delta }x_i0,St\mathrm{\Delta }t0$$
(6)
i.e.,

$$\underset{St\mathrm{\Delta }x_10}{lim}\mathrm{\Delta }x_1=\epsilon _1,\underset{St\mathrm{\Delta }x_20}{lim}\mathrm{\Delta }x_2=\epsilon _2,\underset{St\mathrm{\Delta }x_30}{lim}\mathrm{\Delta }x_3=\epsilon _3,\underset{St\mathrm{\Delta }t0}{lim}\mathrm{\Delta }t=\epsilon _t$$
(7)
In other words, the definitions of ($`\mathrm{\Delta }x_i\epsilon _i`$) and ($`\mathrm{\Delta }t\epsilon _t`$) are, respectively, ($`St\mathrm{\Delta }x_i0`$) and ($`St\mathrm{\Delta }t0`$). Similarly, the limit of tending to infinite does not mean tending to absolute infinite but the nonstandard number $`L`$, i.e.,

$$\underset{St\mathrm{\Delta }x_1\mathrm{}}{lim}\mathrm{\Delta }x_1=L_1,\underset{St\mathrm{\Delta }x_2\mathrm{}}{lim}\mathrm{\Delta }x_2=L_2,\underset{St\mathrm{\Delta }x_3\mathrm{}}{lim}\mathrm{\Delta }x_3=L_3,\underset{St\mathrm{\Delta }t\mathrm{}}{lim}\mathrm{\Delta }t=T$$
(8)
Moreover, in standard case, a function is expressed by $`f(x_1,x_2,x_3,t)`$ which represents the function value at the standard point $`(x_1,x_2,x_3,t)`$. However, in the nonstandard analysis, a point $`(x_1,x_2,x_3,t)`$ is a monad within which there are still infinite nonstandard points $`(x_1^{},x_2^{},x_3^{},t^{})`$. Then a function should be expressed by $`f(x_1,x_2,x_3,t,x_1^{},x_2^{},x_3^{},t^{})`$ which is the function value at the nonstandard point $`(x_1^{},x_2^{},x_3^{},t^{})`$ contained in the monad $`(x_1,x_2,x_3,t)`$. Still we define the partial derivatives in form:

$$\frac{f}{x_i^{}}=\underset{St\mathrm{\Delta }x_i^{}0}{lim}\frac{f\left(x_i^{}+\mathrm{\Delta }x_i^{}\right)f\left(x_i^{}\right)}{\mathrm{\Delta }x_i^{}},\frac{f}{t^{}}=\underset{St\mathrm{\Delta }t^{}0}{lim}\frac{f\left(t^{}+\mathrm{\Delta }t^{}\right)f\left(t^{}\right)}{\mathrm{\Delta }t^{}}$$
(9)
$$\frac{f}{x_i}=\underset{St\mathrm{\Delta }x_i0}{lim}\frac{f\left(x_i+\mathrm{\Delta }x_i\right)f\left(x_i\right)}{\mathrm{\Delta }x_i},\frac{f}{t}=\underset{St\mathrm{\Delta }t0}{lim}\frac{f\left(t+\mathrm{\Delta }t\right)f\left(t\right)}{\mathrm{\Delta }t}$$
(10)
$$\left(i=1,2,3\right)$$

Definition (9)-(10) is only a definition without the ordinary meaning of derivative in the standard case. The difference of $`f\left(x_i+\mathrm{\Delta }x_i\right)f\left(x_i\right)`$ in (10) is the abbreviation of $`f(x_i+\mathrm{\Delta }x_i,x_j,x_k,t,x_i^{},x_j^{},x_k^{},t^{})f(x_i,x_j,x_k,t,x_i^{},x_j^{},x_k^{},t^{})`$, the increment of the function value between two monads of the points $`\left(x_i+\mathrm{\Delta }x_i\right)`$ and $`\left(x_i\right)`$. The distance between the points

$$\underset{St\mathrm{\Delta }x_i0}{lim}(x_i+\mathrm{\Delta }x_i,x_i^{})$$

and

$$(x_i,x_i^{})$$

is infinite from the angle of level $`\beta `$. Therefore, $`\frac{f}{x_i}`$ and $`\frac{f}{t}`$ in (10) have the meanings of an average. If one insists on the concept of ordinary derivative in standard case, $`\frac{f}{x_i}`$ and $`\frac{f}{t}`$ in (10) have no sense.

## 2 Nonstandard Model Of Turbulence

On the basis of the above-mentioned concepts of nonstandard analysis a nonstandard description of turbulent field can be presented as follows.

> Assumption 1: Global turbulent field is composed of standard points, and every standard point is yet a monad. Each monad possesses the internal structure, namely a monad is also composed of infinite nonstandard points (so called interior points of the monad).

It is well known that a space is formed from the points. In fluid mechanics, fluid mechanists always take the fluid element (fluid particle) as point. In fact, a fluid particle has the volume not to be absolute zero but to be a micro volume, which is taken as infinitesimal from macroscopic viewpoint and infinite from microscope. People generally think of the whole fluid particle as uniform. Namely the fluid particle is taken to be equal to a point in reality. Now, in the case of turbulence, the motion of fluid varies so fast that a fluid particle can not be uniform as a whole but becomes a monad. From now on in this paper, the meanings of a fluid particle, a monad on a whole and a standard point are identical with each other and so are the meanings of a fluid particle in the lower level, an interior point of a monad and a nonstandard point.

Since the point of the turbulent field becomes a monad with internal structure, the motion features of different interior points of the monad are different. So it is obvious that there is a flow or flow field in the inner of the monad. Hence, a fluid particle (monad) is made up of nonstandard points. The volume of these nonstandard points is not absolute zero, and these nonstandard points, in fact, are the fluid particles in the lower level. Now the particles in the lower level are thought to be uniform. They are points (nonstandard points) actually. Moreover, a lot of fluid molecules are included still in the fluid particle of the lower level.

Therefore, two kinds of flow fields exist in turbulence. One is the global turbulent field and the other is the flow field in a monad ( called monad field in this paper). They are flow fields in two different levels. Now the global turbulent field is not composed of the points in which there are no structures, but of monad fields. So there are three levels in turbulent field: molecular movement, monad field and global field.

> Assumption 2: The flows in monad fields are controlled by the Navier-Stokes equations..

Then the motion-equations of the flows in monad fields can be written as follows (for incompressible fluid, unsteadiness and three dimensions):

$$\frac{U_i}{\left(\epsilon ^2x^{}\right)_i}=0$$
(11)
$$\frac{U_i}{\epsilon _t^2t^{}}+\frac{U_iU_j}{\left(\epsilon ^2x^{}\right)_j}=\frac{1}{\rho }\frac{\sigma _{ij}}{\left(\epsilon ^2x^{}\right)_j}$$
(12)
Here $`U_i`$ is velocity component in $`i`$-direction, $`\rho `$ density of fluid, $`\sigma _{ij}`$ stress tensor and $`\sigma _{ij}=\sigma _{ji}`$. The independent variables of these functions in (11)-(12) are $`(x_1,x_2,x_3,t,x_1^{},x_2^{},x_3^{},t^{})`$. $`\epsilon _t,\epsilon _1,\epsilon _2,\epsilon _3`$ are infinitesimals (i.e., the linear dimensions of a monad) of time $`t`$, space $`x_1,x_2`$ and $`x_3`$, respectively. Let $`\epsilon _t=\epsilon _1=\epsilon _2=\epsilon _3=\epsilon `$ here and later, then the governing equations becomes:

$$\frac{U_i}{x_i^{}}=0$$
(13)
$$\frac{U_i}{t^{}}+\frac{U_iU_j}{x_j^{}}=\frac{1}{\rho }\frac{\sigma _{ij}}{x_j^{}}$$
(14)
The governing equations (11)-(12) show the internal structure of a monad.

For actual turbulence the possible linear dimensions of flow fields in two levels could be estimated roughly as follows:

Global turbulent field: infinite $`L=\epsilon ^110^{1.5}cm`$

infinitesimal $`\epsilon 10^{1.5}10^2cm`$

(linear dimension of a monad)

Monad field: infinite $`\epsilon 10^{1.5}10^2cm`$

infinitesimal $`\epsilon ^310^{4.5}10^5cm`$

(linear dimension of the particle of lower level)

According to this estimation, one space monad of the turbulent field contains approximately $`10^910^{10}`$ nonstandard points (the particles in the lower level). If the fluid is gas, there are approximately $`10^510^6`$ gas molecules in one nonstandard point of a monad field. This number of gas molecules is enough for stable mean values in the statistical average.

> Assumption 3: Turbulent field is continuous.

Here the meanings of continuity are: firstly the global turbulent field is composed continuously of standard points (i.e., monad fields); secondly the flows in the monad fields are continuous; and thirdly the physical quantities on the interface between two infinitely close monads are continuous as well. The monads of

$$\underset{St\mathrm{\Delta }x_i0}{lim}(x_i+\mathrm{\Delta }x_i,x_j,x_k,t)$$

and

$$(x_i,x_j,x_k,t)$$

are taken as two infinitely close monads in space. And so are in time. The mathematic expression of the third is

$$\underset{St\mathrm{\Delta }x0}{lim}U(x+\mathrm{\Delta }x,0)=U(x,L),\underset{St\mathrm{\Delta }t0}{lim}U(t+\mathrm{\Delta }t,0)=U(t,T)$$
(15)
There is an even more important meaning of the continuity. That will be stated in Assumption 6.

## 3 Uncertainty Of Measurement

To make a physical measurement, some time-space point $`(x_1,x_2,x_3,t)`$ has to be assigned first, at which the measurement is carried out. The results of the measurement are taken as the measuring data of this point $`(x_1,x_2,x_3,t)`$. According to the standpoint of nonstandard analysis, the point $`(x_1,x_2,x_3,t)`$ is not an absolute point but a monad, in which there are infinite nonstandard points. Generally speaking, the motion features of different nonstandard points are different. So the meaning of saying “measuring at point $`(x_1,x_2,x_3,t)`$” becomes ambiguous. Obviously, if there are fields in two different levels in the system being studied, there must exist the new meaning of a measurement at any point of the field in higher level.

> Assumption 4: When a measurement at any point (monad) $`(x_1,x_2,x_3,t)`$ in a physical field is taken, the operation of the measurement will act randomly on one interior point (nonstandard point) of the point $`(x_1,x_2,x_3,t)`$.

This assumption states clearly that one measurement operation at some standard point (monad) in a field will always take place at an interior point of such a monad. But the measurer can not determine which interior point is acted on by the measurement operation. What actually happens is random. The measurer knows only the standard point (monad) in which the interior point acted on by one measurement operation is contained. It is surprising that when a measurement is carried out there is no way to determine the exact position of the measured object. That is the uncertainty of measurement. At first sight it seems that the uncertainty of measurement stems from the rough measuring skill and that the circumstance conditions of the measurement are not exactly controllable etc.. After further thinking, it is comprehended that the uncertainty rises essentially from limited ability and means of the description of the physical world in hierarchical structures. Apparently, in their quest of the laws of the universe, people need not (and can not in fact) describe all motion states of every molecule and electron in the universe. Conversely, even though an exact and perfect description of the motion states of every molecule and electron in the universe is given, the laws of the motion of the whole universe are still not understood. It is clearly that because the physical world has a hierarchical structure, therefore any description and measurement of physical phenomena should be related with the hierarchical structure. When the effect of hierarchical structure can not be ignored, the uncertainty of measurement occurs. The measurement and observation of turbulence can be taken for examples of this case.

An operation of measurement will act randomly on various interior points in a monad. By what probability does the measurement operation act on each interior point of a monad? In turbulence it is thought that one operation of measurement will take place in equiprobability at each interior point of a monad.

> Assumption 5: When a measurement at any point (monad) of a turbulent field is made, the operation of the measurement will act in equiprobability on various interior points of the monad. This assumption is called the equiprobability assumption.

So there are two kinds of averages. One is that one measurement is just taking average over a large number of molecules contained in some interior point of a monad, while the other average over the motions of all interior points in a monad has to be taken in order to get the mean values of the physical quantities over the standard point(monad) of the turbulent field. According to the equiprobability assumption, the second average formula can be given. Let $`U(\stackrel{}{x},t,\stackrel{}{x^{}},t^{})`$ represent some physical quantity at an interior point $`(\stackrel{}{x^{}},t^{})`$ of the monad $`(\stackrel{}{x},t)`$ and $`\stackrel{~}{U}(\stackrel{}{x},t)`$ express the average of $`U(\stackrel{}{x},t,\stackrel{}{x^{}},t^{})`$ over all interior points of the monad. Then it follows that

$$\stackrel{~}{U}(\stackrel{}{x},t)=\frac{1}{T}_0^T𝑑t^{}\frac{1}{L^3}_0^L𝑑x_1^{}_0^L𝑑x_2^{}_0^L𝑑x_3^{}U(\stackrel{}{x},t,\stackrel{}{x^{}},t^{})$$
(16)
Where the monad $`(\stackrel{}{x},t)(x_1,x_2,x_3,t)`$, the interior point of the monad $`(\stackrel{}{x^{}},t^{})(x_1^{},x_2^{},x_3^{},t^{})`$, and T, L are infinites of time and space respectively. Let $`T=L_1=L_2=L_3=L`$ here and later. Obviously, $`\stackrel{~}{U}`$ is only a function of the standard point $`(\stackrel{}{x},t)`$.

If the monad fields of a global flow field are uniform, there is no difference between saying “one measurement acting on the interior point of a monad” and “one measurement acting on the point (monad)”. Then the global field is laminar. Exactly,the essential feature of turbulence is that its particle (standard point) has internal structure, i.e., internal flow field (monad field). That could be taken as the definition of turbulence.

Because of the uncertainty of measurement, one measurement at the point $`(\stackrel{}{x},t)`$ will act randomly on any interior point. Therefore, the results of many measurements carried out at the same point(monad) under constant circumstance conditions will indicate randomly the motion states of various interior points of the monad. Surely fluctuation of data of the measurements must appear, provided that there are considerable different states among various interior points of the monad. This is the mechanism of fluctuation appearing in a turbulent field.

Thus, a turbulent field is composed of monads with interior structure (monad field). The flows in the monad fields are controlled by the governing equations of (13)-(14). Still continuity of turbulent field exists. Therefore, turbulence is also regular flow. The appearance of the fluctuation of a turbulent flow stems from the uncertainty of measurement and observation in the turbulent field.

To keep the concepts mentioned above in mind, the order of magnitude of fluctuation in a turbulent field can be estimated roughly. Let $`\mathrm{\Delta }x^{}`$ be finite length in the coordinate system of a monad, $`U0\left(1\right)`$ any physical quantity and its fluctuation $`uU\left(x^{}+\mathrm{\Delta }x^{}\right)U\left(x^{}\right)\frac{U}{x^{}}\mathrm{\Delta }x^{}`$. The actual rate of change (i.e., observed from the angle of the global field level) of function U is $`\frac{U}{\epsilon ^2x^{}}`$. In one case $`\frac{U}{x^{}}0\left(\epsilon ^2\right)`$, $`u0\left(\epsilon ^2\right)`$, no obvious fluctuation occurs. In another case $`\frac{U}{x^{}}0\left(1\right)`$, and $`\frac{U}{\epsilon ^2x^{}}0\left(\frac{1}{\epsilon ^2}\right)`$ showing that there is excess rate of shear strain in the field. This will lead to instability in the flow field. So the case should be excluded. Then in the third case $`\frac{U}{x^{}}0\left(\epsilon \right)`$ only, and $`\frac{U}{\epsilon ^2x^{}}0\left(\frac{1}{\epsilon }\right)`$ with fluctuation $`u0\left(\epsilon \right)`$. So there is notable fluctuation in order of magnitude $`u0\left(\epsilon \right)`$ in turbulence provided that there exists any set of interior points in a monad, and on the set $`\frac{U}{x^{}}0\left(\epsilon \right)`$.

## 4 Conservation Equations

Generally speaking, in fluid mechanics a small volume $`\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3`$ is taken and the flux of conservative physical quantity through the boundary planes of the volume is computed. Then let $`\mathrm{\Delta }x_10,\mathrm{\Delta }x_20,\mathrm{\Delta }x_30`$, the small volume will tend to the point $`(x_1,x_2,x_3)`$. Thus, conservation equations at the point $`(x_1,x_2,x_3)`$ are obtained. The boundary planes of the small volume are absolute planes without thickness. However, in nonstandard case a small volume is still taken. By Use of the same method as in standard case but let $`St\mathrm{\Delta }x_10,St\mathrm{\Delta }x_20,St\mathrm{\Delta }x_30`$, the small volume $`\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3`$ will become the monad of $`(x_1,x_2,x_3)`$. Now the boundary planes of the small volume are not absolute planes but very thin layers with thickness of $`\epsilon `$ (the linear dimension of the monad). Also conservation equations at this monad are obtained. Surely the equations are those of mean quantities over the monad in which infinite interior points are contained.

Now these conservation equations are given in the following.

Let $`\mathrm{`}\mathrm{`}\mathrm{"}`$ express the average operation over $`(x_1^{},x_2^{},x_3^{},t^{})`$, $`\mathrm{`}\mathrm{`}^i\mathrm{"}`$ over $`x_i^{}`$, $`\mathrm{`}\mathrm{`}\mathrm{"}`$ over $`t^{}`$, etc., in a monad.

Through the pair of the planes perpendicular to $`x_1`$ axis the mass inflow per unit time is:

$$\stackrel{~}{\stackrel{~}{\overline{U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(x_1\right)\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

the mass outflow

$$\stackrel{~}{\stackrel{~}{\overline{U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(x_1+\mathrm{\Delta }x_1\right)\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

the net inflow

$$\frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

Through the other two pairs of the planes perpendicular to the two other axes the net mass inflow per unit time are, respectively,

$$\frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{U_2}}^2}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

and

$$\frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{U_3}}^3}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

Because of the law of mass conservation, it follows that

$$\frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}+\frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{U_2}}^2}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_3}+\frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{U_3}}^3}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2}=0$$

Here the average of $`\mathrm{`}\mathrm{`}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{"}`$ is taken over the plane of $`\mathrm{\Delta }x_1\mathrm{\Delta }x_2`$, which tends to $`\mathrm{`}\mathrm{`}12\mathrm{"}`$ when $`St\mathrm{\Delta }x_10,St\mathrm{\Delta }x_20`$. And so are the others. Let $`St\mathrm{\Delta }x_10,St\mathrm{\Delta }x_20,St\mathrm{\Delta }x_30`$, then it is obtained that

$$\frac{}{x_1}\stackrel{~}{\overline{U_1}}^{123}+\frac{}{x_2}\stackrel{~}{\overline{U_2}}^{123}+\frac{}{x_3}\stackrel{~}{\overline{U_3}}^{123}=0$$

i.e.,

$$\frac{\stackrel{~}{U_1}}{x_1}+\frac{\stackrel{~}{U_2}}{x_2}+\frac{\stackrel{~}{U_3}}{x_3}=0$$
(17)
Similarly, the equations of momentum conservation can be obtained too.

Through the pair of the planes perpendicular to $`x_1`$ axis, the net inflow of the momentum in $`x_1`$ direction per unit time is

$$\rho \stackrel{~}{\stackrel{~}{\overline{U_1U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(x_1\right)\mathrm{\Delta }x_2\mathrm{\Delta }x_3\rho \stackrel{~}{\stackrel{~}{\overline{U_1U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(x_1+\mathrm{\Delta }x_1\right)\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

$$=\rho \frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{U_1U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

Through the other two pairs of the planes perpendicular to the other two axes, the net inflows of the momentum in $`x_1`$ direction per unit time are

$$\rho \frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{U_2U_1}}^2}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

and

$$\rho \frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{U_3U_1}}^3}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

respectively.

The force in $`x_1`$ direction exerting on the pair of the planes perpendicular to $`x_1`$ axis is:

$$\stackrel{~}{\stackrel{~}{\overline{\sigma _{11}}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(x_1\right)\mathrm{\Delta }x_2\mathrm{\Delta }x_3+\stackrel{~}{\stackrel{~}{\overline{\sigma _{11}}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(x_1+\mathrm{\Delta }x_1\right)\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

$$=\frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{\sigma _{11}}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

The forces in $`x_1`$ direction exerting on the other two pairs of the planes perpendicular to the other two axes are

$$\frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{\sigma _{21}}}^2}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

and

$$\frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{\sigma _{31}}}^3}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

respectively.

Then the increments of the momentum in $`x_1`$ direction in the small volume per unit time are

$$\rho \left\{\stackrel{~}{\overline{U_1}}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(t+\mathrm{\Delta }t\right)\stackrel{~}{\overline{U_1}}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\left(t\right)\right\}\frac{\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3}{\mathrm{\Delta }t}$$

$$=\rho \frac{}{t}\stackrel{~}{\overline{U_1}}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3$$

Because of the law of momentum conservation, it follows that

$$\rho \frac{}{t}\stackrel{~}{\overline{U_1}}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3}+\rho \frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{U_1U_1}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}+\rho \frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{U_2U_1}}^2}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_3}+\rho \frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{U_3U_1}}^3}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2}$$

$$=\frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{\sigma _{11}}}^1}^{\mathrm{\Delta }x_2\mathrm{\Delta }x_3}+\frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{\sigma _{21}}}^2}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_3}+\frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{\sigma _{31}}}^3}^{\mathrm{\Delta }x_1\mathrm{\Delta }x_2}$$

Let $`St\mathrm{\Delta }x_10,St\mathrm{\Delta }x_20,St\mathrm{\Delta }x_30`$, so $`\mathrm{`}\mathrm{`}\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3\mathrm{"}`$ tend to $`\mathrm{`}\mathrm{`}123\mathrm{"}`$, etc. Then it follows that

$$\frac{}{t}\stackrel{~}{\overline{U_1}}^{123}+\frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{U_1U_1}}^1}^{23}+\frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{U_2U_1}}^2}^{13}+\frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{U_3U_1}}^3}^{12}$$

$$=\frac{1}{\rho }\frac{}{x_1}\stackrel{~}{\stackrel{~}{\overline{\sigma _{11}}}^1}^{23}+\frac{1}{\rho }\frac{}{x_2}\stackrel{~}{\stackrel{~}{\overline{\sigma _{21}}}^2}^{13}+\frac{1}{\rho }\frac{}{x_3}\stackrel{~}{\stackrel{~}{\overline{\sigma _{31}}}^3}^{12}$$

i.e.,

$$\frac{}{t}\stackrel{~}{U_1}+\frac{}{x_1}\stackrel{~}{U_1U_1}+\frac{}{x_2}\stackrel{~}{U_2U_1}+\frac{}{x_3}\stackrel{~}{U_3U_1}$$

$$=\frac{1}{\rho }\left(\frac{}{x_1}\stackrel{~}{\sigma _{11}}+\frac{}{x_2}\stackrel{~}{\sigma _{21}}+\frac{}{x_3}\stackrel{~}{\sigma _{31}}\right)$$

Similarly, the momentum equations in $`x_2`$ and $`x_3`$ directions can be written too. The general form of momentum equations is:

$$\frac{\stackrel{~}{U_i}}{t}+\frac{\stackrel{~}{U_iU_j}}{x_j}=\frac{1}{\rho }\frac{\stackrel{~}{\sigma _{ij}}}{x_j}$$
(18)
The equations $`\left(17\right)\left(18\right)`$ are about the mean quantities over a monad and called the conservation equations. Here no special mean-methods are assigned. But the mean $`\mathrm{`}\mathrm{`}\mathrm{"}`$ will always be taken as the mean-method showed in (16) in this paper.

## 5 Fundamental Equations of Turbulent Flows

Still instantaneous quantities $`U_i`$ and $`P`$ are decomposed into two parts such that

$$U_i=\stackrel{~}{U_i}+u_i,P=\stackrel{~}{P}+p,\stackrel{~}{u_i}=0,\stackrel{~}{p}=0$$
(19)
Here the quantities with $`\mathrm{`}\mathrm{`}\mathrm{"}`$ are called mean quantities; in this regard, we will adopt the mean-method of (16) i.e.,

$$\stackrel{~}{U}=\frac{1}{T}_0^T𝑑t^{}\frac{1}{L^3}_0^L𝑑x_1^{}_0^L𝑑x_2^{}_0^L𝑑x_3^{}U,\overline{U}=\frac{1}{T}_0^T𝑑t^{}U,\stackrel{~}{U}^i=\frac{1}{L}_0^L𝑑x_i^{}U$$
(20)
$$\stackrel{~}{U}^{ij}=\frac{1}{L^2}_0^L𝑑x_i^{}_0^L𝑑x_j^{}U,\stackrel{~}{U}^{ijk}=\frac{1}{L^3}_0^L𝑑x_i^{}_0^L𝑑x_j^{}_0^L𝑑x_k^{}U$$

The relations of (19) are the definitions of $`u_i`$ and $`p`$. The $`u_i`$ and $`p`$ have the same the order of magnitude as the turbulent fluctuation mentioned in Section 3. We will call the $`u_i`$ and $`p`$ in (19) as fluctuant quantities. To note that the $`\stackrel{~}{U_i}`$ and $`\stackrel{~}{P}`$ are mean quantities over a monad, and $`u_i`$ and $`p`$ are meaningful to a monad (or in a monad). They have different meanings with those in the existent theory of turbulence.

Before derivation of the fundamental equations of turbulent flows, it is necessary to present the assumption of close property between two monads when they are infinitely close to each other, and analyze the mutual relation of the two monads in motion characteristics.

> Assumption 6: In both the value and structure of function, physical function, defined on the interior points of the monads of a turbulent field, is very close between two monads, when these two monads are infinitely close to each other.

Considering any physical quantity $`U_i0\left(1\right)`$, its fluctuation $`u_i`$, the mathematical expressions of this assumption are as follows:

$$\left(a\right)\delta u_i0\left(\epsilon ^2\right),\left(b\right)\delta \left(u_iu_j\stackrel{~}{u_iu_j}\right)0\left(\epsilon ^4\right)$$
(21)
$$\left(0x_i^{}<L,0t^{}<T\right)$$

Here the operator $`\delta `$ is defined as

$$\delta _iU=\underset{St\mathrm{\Delta }x_i0}{lim}\left[U(x_i+\mathrm{\Delta }x_i,x_i^{})U(x_i,x_i^{})\right],\delta _tU=\underset{St\mathrm{\Delta }t0}{lim}\left[U(t+\mathrm{\Delta }t,t^{})U(t,t^{})\right]$$
(22)
It can be got from (a) in (21) that

$$\delta \left[\frac{U}{\epsilon ^2x_i^{}}\frac{U}{x_i}\right]0\left(\epsilon ^2\right)$$
(23)
Where $`\frac{U}{\epsilon ^2x_i^{}}`$ is the actual rate of change or figure slope of the function $`U`$ at the nonstandard point $`x_i^{}`$, and $`\frac{U}{x_i}`$ is the mean value of $`\frac{U}{\epsilon ^2x_i^{}}`$ over $`x_i^{}`$ (see(25)). Clearly the condition (23) and $`\left(a\right)`$ in (21) show that the figure shape of the function $`U`$ has very little difference between two infinitely close monads. Similarly the condition $`\left(b\right)`$ in (21) shows that the figure shape of the function $`u_iu_j`$ has very little difference between two infinitely close monads.

Assumption 6 could be called “Continuous Assumption”. It is well known that with a few exceptions, physical quantities vary usually continuously with time and space. In the case of standard analysis, a function is defined on the standard points. It is characterized only by its value. The continuity of a physical quantity only requires that the function representing the physical quantity be continuous in mathematics. On the other hand, in the nonstandard case a point becomes a monad and a function is defined at the interior point of the monad. There exists a function-structure on the whole monad. So the function is characterized by not only its value but also its structure on the monad. The continuity of a physical function will mean that in both the values and figure-shape of the function , the variation of the physical function with monads is continuous. Thus, there is very little difference between two infinitely close monads in both value and structure of the physical function. In other words, the Assumption 6 is the nonstandard expansion of the physical function continuity in the case of standard analysis. Moreover, the Assumption 6 is reasonable and natural, because the circumstance conditions imposed on two infinitely close monads are very close to each other.

In a monad, we consider that stress tensor for Newtonian fluid is

$$\sigma _{ij}=P\delta _{ij}+\mu \left(\frac{U_i}{x_j^{}}+\frac{U_j}{x_i^{}}\right)\frac{1}{\epsilon ^2}$$
(24)
Using (20) and(17),(15) (Assumption 3) we can give the following computation:

$$\frac{1}{\epsilon ^2}\stackrel{~}{\frac{U_i}{x_j^{}}}^j=\left(\frac{U_i}{x_j}\right)_{x_j^{}=0}$$
(25)
Then

$$\stackrel{~}{\sigma _{ij}}=\stackrel{~}{P}\delta _{ij}+\mu \left(\frac{\stackrel{~}{\overline{U_i}}^{kl}}{x_j}\right)_{x_j^{}=0}+\mu \left(\frac{\stackrel{~}{\overline{U_j}}^{mn}}{x_i}\right)_{x_i^{}=0}$$

$`{\displaystyle \frac{\stackrel{~}{\sigma _{ij}}}{x_j}}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{P}}{x_i}}+\mu {\displaystyle \frac{}{x_j}}\left({\displaystyle \frac{\stackrel{~}{\overline{U_i}}^{kl}}{x_j}}\right)_{x_j^{}=0}+\mu {\displaystyle \frac{}{x_j}}\left({\displaystyle \frac{\stackrel{~}{\overline{U_j}}^{mn}}{x_i}}\right)_{x_i^{}=0}`$
$`=`$ $`{\displaystyle \frac{\stackrel{~}{P}}{x_i}}+\mu {\displaystyle \frac{}{x_j}}{\displaystyle \frac{\stackrel{~}{U_i}}{x_j}}+\mu {\displaystyle \frac{}{x_j}}\left({\displaystyle \frac{\stackrel{~}{\overline{u_i}}^{kl}}{x_j}}\right)_{x_j^{}=0}`$
By use of (a) in (21) it follows that

$$\left(\frac{^2\stackrel{~}{\overline{u_i}}^{kl}}{x_j^2}\right)_{x_j^{}=0}0\left(\epsilon \right)$$
(26)
To omit them, it is obtained that

$$\frac{\stackrel{~}{\sigma _{ij}}}{x_j}=\frac{\stackrel{~}{P}}{x_i}+\mu \frac{^2\stackrel{~}{U_i}}{x_j^2}$$
(27)
Then the equations (18) can be written as

$$\frac{\stackrel{~}{U_i}}{t}+\frac{\stackrel{~}{U_iU_j}}{x_j}=\frac{1}{\rho }\frac{\stackrel{~}{P}}{x_i}+\nu ^2\stackrel{~}{U_i},^2=\frac{}{x_j}\frac{}{x_j}$$
(28)
(17) is written as

$$\frac{\stackrel{~}{U_i}}{x_i}=0$$
(29)
Here $`\nu =\frac{\mu }{\rho }`$ the kinematic viscosity.

Note that now the mean method in the equations (28) and (29) is defined by (20). Then by use of the relation (19), the decomposition of (29) and (28) can be done as

$$\frac{U_i}{x_i}\frac{u_i}{x_i}=0$$
(30)
$$\frac{U_i}{t}+\frac{U_iU_j}{x_j}\frac{u_i}{t}\frac{\left(U_iU_j\right)_{fl}}{x_j}$$

$$=\frac{1}{\rho }\frac{P}{x_i}+\frac{1}{\rho }\frac{p}{x_i}+\nu ^2U_i\nu ^2u_i$$
(31)
Here $`u_i,p`$ and $`\left(U_iU_j\right)_{fl}`$ are the fluctuations of the velocity, pressure and $`\left(U_iU_j\right)`$ respectively. The independent variables of the functions in the equations (30)-(31) are $`(x_1,x_2,x_3,t,x_1^{},x_2^{},x_3^{},t^{})`$. Based on Assumption 6 the order of magnitude of the terms concerned with fluctuation quantities is $`0\left(\epsilon \right)`$, but the order of magnitude of the terms concerned with instantaneous quantities is $`0\left(1\right)`$. So when (30)-(31) are split into two parts having different order of magnitude, the following can be obtained:

Equations, in which terms are in the order of magnitude $`0\left(1\right)`$, of instantaneous quantities are:

$$\frac{U_i}{x_i}=0$$
(32)
$$\frac{U_i}{t}+\frac{U_iU_j}{x_j}=\frac{1}{\rho }\frac{P}{x_i}+\nu ^2U_i$$
(33)
Equations, in which terms are in the order of magnitude $`0\left(\epsilon \right)`$, of fluctuation quantities are:

$$\frac{u_i}{x_i}=0$$
(34)
$$\frac{u_i}{t}+\frac{\left(U_iU_j\right)_{fl}}{x_j}=\frac{1}{\rho }\frac{p}{x_i}+\nu ^2u_i$$
(35)
Where $`\left(U_iU_j\right)_{fl}`$ in (35) is:

$`\left(U_iU_j\right)_{fl}`$ $`=`$ $`U_iU_j\stackrel{~}{U_iU_j}`$
$`=`$ $`\stackrel{~}{U_i}u_j+\stackrel{~}{U_j}u_i+u_iu_j\stackrel{~}{u_iu_j}`$
Substitution of the relations into (35) produces

$$\frac{u_i}{t}+u_j\frac{\stackrel{~}{U_i}}{x_j}+\stackrel{~}{U_j}\frac{u_i}{x_j}+\frac{u_iu_j}{x_j}\frac{\stackrel{~}{u_iu_j}}{x_j}=\frac{1}{\rho }\frac{p}{x_i}+\nu ^2u_i$$
(36)
Equations (28)(29),(32)(33) and (34)(35) are the governing equations about mean, instantaneous and fluctuant quantities respectively. They can be called as the fundamental equations of turbulent flows.

It is found easily that these equations have the same form as those in the existent theory of turbulent flows. But there are essentially differences between the two. These differences are:

1.Different methods to get them. In the existent theory of turbulence the mean equations are obtained by the average of the instantaneous equations (i.e., the Navier-Stokes equations). Then the mean equations are subtracted from the instantaneous equations to get the fluctuation equations. But in this paper, the mean equations will FIRST be obtained in view of the physical conservation laws, and then by virtue of the close property assumption, the instantaneous and fluctuant equations are obtained respectively.

2.The terms in the equations of the existent theory are the limits, such as

$$\underset{\mathrm{}x0}{lim}\frac{f\left(x+\mathrm{}x\right)f\left(x\right)}{\mathrm{}x}\underset{\mathrm{}t0}{lim}\frac{f\left(t+\mathrm{}t\right)f\left(t\right)}{\mathrm{}t}$$

But the terms in the equations of the nonstandard theory are the partial derivatives defined in (10). The equations in two theories are conceptually different. The former is based on, so called, the frame of $`\delta \epsilon `$, the latter is out of the frame of $`\delta \epsilon `$. Therefore, the equation in nonstandard picture of turbulence is the new kind of equation.

3. The physical quantities in the equations of the existent theory, such as velocity, pressure etc., are all functions of the standard point $`(\stackrel{}{x},t)`$. However, those in the nonstandard case are all functions of the point $`(\stackrel{}{x},t,\stackrel{}{x^{}},t^{})`$, actually functions of the nonstandard point $`(\stackrel{}{x^{}},t^{})`$ with $`(\stackrel{}{x},t)`$ showing which monad contains it. It is reasonable that the physical quantities are taken as the functions of the standard point $`(\stackrel{}{x},t)`$ in laminar flow, but not in turbulent flows, because a standard point becomes a monad in which there is an interior structure (i.e., the monad flow field) in the case of turbulence. In fact, there is no definition of instantaneous and fluctuant quantities of a whole monad, but there is the definition of instantaneous and fluctuant quantities of the interior points of the monad. So for the global field,the Navier-Stokes equations in existent meanings do not exist in principle. Therefore, though (32)-(33) are the same as the Navier-Stokes equations in forms, they are different. The fundamental equations obtained in this paper are governing equations about the physical quantities defined at those interior points of various monads. These interior points have the same nonstandard coordinates $`(\stackrel{}{x^{}},t^{})`$. But the fundamental equations in the existent turbulent theory are governing equations about the physical quantities defined at the standard points of the global flow field. There is an essential difference between them.

4. The average in the existent theory usually means the ensemble average, but the distribution function can not be given. So the closure problem continues to exist. Otherwise, it is definite in the nonstandard model to take average over all interior points of a monad by the equiprobability assumption.

5. It should be pointed out that the function values at some standard points are given for the initial and boundary conditions of the fundamental equations in the existent theory. However, only the function values at some nonstandard points can be taken as the initial and boundary conditions for the instantaneous and fluctuation quantities in the fundamental equations in the nonstandard case. Because a turbulent field always starts and develops from walls or laminar flows and evolves from laminar flows or static fluids, the function values on walls or laminar flows can be taken as the boundary conditions, and the conditions of static fluids or laminar flows can be taken as the initial conditions. No slip condition is still kept. Thus, the initial and boundary conditions of the fundamental equations obtained from the nonstandard analysis can be given in principle.

## 6 On Closure Problem

Substitution of the relation $`U_i=\stackrel{~}{U_i}+u_i`$ into (28) produces

$$\frac{\stackrel{~}{U_i}}{t}+\frac{\stackrel{~}{U_i}\stackrel{~}{U_j}}{x_j}+\frac{\stackrel{~}{u_iu_j}}{x_j}=\frac{1}{\rho }\frac{\stackrel{~}{P}}{x_i}+\nu ^2\stackrel{~}{U_i}$$
(37)
This is the same in form as famous Reynolds equation. It is not closure, and neither is equation (36) in which there is term $`\stackrel{~}{u_iu_j}`$. The closure problem stemming from nonlinearity is very difficult in turbulence research. Now in the nonstandard theory, the turbulence is composed of monad fields, and the point-average is adopted. Then by the close property between two infinitely close monads, we can give the reasonable closure methods as follows.

By the relation (21)(Assumption 6) and the order of magnitude of fluctuation $`u_i0\left(\epsilon \right)`$, it follows that

$$\delta \left(u_iu_j\right)=u_i\delta u_j+u_j\delta u_i0\left(\epsilon ^3\right)$$

$$\delta \stackrel{~}{u_iu_j}0\left(\epsilon ^3\right)$$

$$\frac{\stackrel{~}{u_iu_j}}{x_j}0\left(\epsilon ^2\right)$$

and

$$\delta \left(u_iu_j\stackrel{~}{u_iu_j}\right)0\left(\epsilon ^4\right)$$

also,it follows that

$$\frac{\left(u_iu_j\stackrel{~}{u_iu_j}\right)}{x_j}0\left(\epsilon ^3\right)$$
(38)
There are three choices for the closure of equations:

(A) Choice one:

The term $`\frac{\stackrel{~}{u_iu_j}}{x_j}0\left(\epsilon ^2\right)`$ is neglected from equations (37). The closed equations of the mean quantities will be obtained as follows.

$$\frac{\stackrel{~}{U_i}}{x_i}=0$$
(39)
$$\frac{\stackrel{~}{U_i}}{t}+\frac{\stackrel{~}{U_i}\stackrel{~}{U_j}}{x_j}=\frac{1}{\rho }\frac{\stackrel{~}{P}}{x_i}+\nu ^2\stackrel{~}{U_i}$$
(40)
(B) Choice two:

To note that

$$\frac{\left(u_iu_j\stackrel{~}{u_iu_j}\right)}{x_j}0\left(\epsilon ^3\right)\left(38\right)$$

and the other terms in equation (36) have the order of magnitude of $`0\left(\epsilon \right)`$. Hence, the term (38) can be neglected from (36). Then the equation (36) becomes

$$\frac{u_i}{t}+u_j\frac{\stackrel{~}{U_i}}{x_j}+\stackrel{~}{U_j}\frac{u_i}{x_j}=\frac{1}{\rho }\frac{p}{x_i}+\nu ^2u_i$$
(41)
With $`\stackrel{~}{U_i}=U_iu_i`$, this equation changes to

$$\frac{u_i}{t}+U_j\frac{u_i}{x_j}+u_j\frac{U_i}{x_j}2u_j\frac{u_i}{x_j}=\frac{1}{\rho }\frac{p}{x_i}+\nu ^2u_i$$
(42)
Under certain initial and boundary conditions, the equations (32)-(33) are solved and $`U_i`$ are obtained. Then the equations (42) with (34) becomes also closed, the fluctuations $`u_i`$ and $`p`$ can be obtained after the equations (34)(42) are solved. Finally the mean quantities $`\stackrel{~}{U_i}=U_iu_i`$ and $`\stackrel{~}{P}=Pp`$.

(C) Choice three:

Hence, the term (38) can be neglected from (36) and (37). Then it follows that

$$\frac{\stackrel{~}{U_i}}{x_i}=0$$
(43)
$$\frac{\stackrel{~}{U_i}}{t}+\frac{\stackrel{~}{U_i}\stackrel{~}{U_j}}{x_j}+\frac{u_iu_j}{x_j}=\frac{1}{\rho }\frac{\stackrel{~}{P}}{x_i}+\nu ^2\stackrel{~}{U_i}$$
(44)
$$\frac{u_i}{x_i}=0$$
(45)
$$\frac{u_i}{t}+\stackrel{~}{U_j}\frac{u_i}{x_j}+u_j\frac{\stackrel{~}{U_i}}{x_j}=\frac{1}{\rho }\frac{p}{x_i}+\nu ^2u_i$$
(46)
Obviously, the equations (43)-(46) are closed. We can obtain $`\stackrel{~}{U_i},\stackrel{~}{P},u_i,p,`$ from solving these equations.

Surely if $`u_i=0,p=0`$, i.e., no fluctuations occur, the equations (28)(29) and (32)(33) are reduced to ordinary Navier-Stokes equations and the field becomes a laminar field.

## 7 Some Remarks on the Numerical Calculations

To take the $`\frac{U}{t}`$ for example, in the standard case we have

$$\left(\frac{U}{t}\right)_{st}=\underset{\mathrm{\Delta }t0}{lim}\frac{U(x_1,x_2,x_3,t+\mathrm{\Delta }t)U(x_1,x_2,x_3,t)}{\mathrm{\Delta }t}$$

but in the nonstandard case by the definition of $`\left(\frac{U}{t}\right)_{nst}`$ in this paper, it follows that

$$\left(\frac{U}{t}\right)_{nst}=\underset{St\mathrm{\Delta }t0}{lim}\frac{U(x_1,x_2,x_3,t+\mathrm{\Delta }t,x_1^{},x_2^{},x_3^{},t^{})U(x_1,x_2,x_3,t,x_1^{},x_2^{},x_3^{},t^{})}{\mathrm{\Delta }t}$$
(47)
In form, there is no difference between the two cases when the discretization of time-space is taken for the numerical computation. However within a grid, obtained from the discretization of $`\left(\frac{U}{t}\right)_{st}`$, there is not structure (i.e., the whole grid should be uniform). But the interior structure is permitted through the grid obtained from the discretization of $`\left(\frac{U}{t}\right)_{nst}`$. This diference between two cases is very important.

Surely in the nonstandard case, $`U(x_1,x_2,x_3,t+\mathrm{\Delta }t,x_1^{},x_2^{},x_3^{},t^{})`$ and $`U(x_1,x_2,x_3,t,x_1^{},x_2^{},x_3^{},t^{})`$ are the values of physical quantity $`U`$ at the interior point $`(x_1^{},x_2^{},x_3^{},t^{})`$ of two different monads in time, that is monad of $`\left(t+\mathrm{\Delta }t\right)`$ and monad of $`\left(t\right)`$. Though the fundamental equations work for any interior point $`(x_1^{},x_2^{},x_3^{},t^{})`$ of every monad, it is better that let $`x_1^{}=x_2^{}=x_3^{}=t^{}=0`$ because the point $`(x_1,x_2,x_3,t,0,0,0,0)`$ is located on the boundary of the time-space of a turbulent field, when $`(x_1,x_2,x_3,t)`$ is the boundary monad of the field.

Moreover, $`U(x_1,x_2,x_3,t+\mathrm{\Delta }t,0,0,0,0)`$ and $`U(x_1,x_2,x_3,t,0,0,0,0)`$ are hyperreal numbers. In the numerical calculations, their standard parts will be taken.

Therefore, after the discretization, $`\frac{U}{t}`$ becomes:

$$\frac{StU(x_1,x_2,x_3,t+\mathrm{\Delta }t,0,0,0,0)StU(x_1,x_2,x_3,t,0,0,0,0)}{\mathrm{\Delta }t}$$

For simplicity it can also be written as

$$\frac{U(x_1,x_2,x_3,t+\mathrm{\Delta }t,0,0,0,0)U(x_1,x_2,x_3,t,0,0,0,0)}{\mathrm{\Delta }t}$$
(48)
provided it is kept in mind that the $`U`$ in (48) represents the standard parts of them.

There are similar discussions for the space partial derivertives of physical quantities.

After the $`U(x_1,x_2,x_3,t,0,0,0,0)`$ and $`u(x_1,x_2,x_3,t,0,0,0,0)`$, the instantaneous and fluctuant quantities respectively, are obtained, then

$$\stackrel{~}{U}(x_1,x_2,x_3,t)=U(x_1,x_2,x_3,t,0,0,0,0)u(x_1,x_2,x_3,t,0,0,0,0)$$
(49)
That is the mean value over the small volume of $`\mathrm{\Delta }x_1\mathrm{\Delta }x_2\mathrm{\Delta }x_3\mathrm{\Delta }t`$, but it is taken as the mean value over the point $`(x_1,x_2,x_3,t)`$.

Moreover, the mean value $`\stackrel{~}{U}`$ obtained by the methods mentioned above, is the mean over the point(monad) and still random oscillatory. The random oscillation of the mean values $`\stackrel{~}{U}`$ ought to stem from the unavoidable random disturbances in real turbulent fields. If the average of these mean values over a finite time period(or a finite space range) is taken once again, the final average results could be compared with the measured average values of physical quantities over corresponding time period(or space range).

## 8 Conclusions

In a view of the concepts of the nonstandard analysis and the fact that the physical world has a hierarchical structure, one nonstandard picture on turbulent flows has been presented in this paper. The key points of this picture are as follows:

1. There exist two kinds of fields in different levels. One is the global turbulent field composed of the standard points (monads), and the other is the monad field. And as a whole a monad is not already uniform. Now the interior structure occurs in a monad; namely, different interior points of the monad possess different characteristics of motion.

2. The flows in a monad field are governed by the Navier-Stokes equations.

3. If there are two kinds of fields in different levels, one operation of measurement at any point (actually a monad) of the field in a higher level will act randomly on one interior point of the monad field. The measurer can determine the monad, rather than its interior point on which the measurement operation acts.

The position of measured object can not be determined exactly, but over a range. Evidently there exists some uncertainty of measurement. That uncertainty stems from the description and observation of the physical phenomena in a hierarchical structure.

4. The flows in a turbulent field are also continuous and regular, while they are seen intuitively in disorderly and irregular form. Fundamentally, the disorder, irregularity and fluctuation occur because of the uncertainty of measurement. Hence, we can probably say that the disorder of turbulence stems from two sources: the first is the fluctuation because of the uncertainty of measurement. This kind of fluctuation is thought to be “real turbulent fluctuation”. The second is the violent unsteadiness and unhomogeneity of the field. The values of $`\stackrel{~}{U_i}`$ and $`\stackrel{~}{P}`$ vary fast with time and space. It is due to the flow-instability as the Reynolds number is very high. In this case, the small random disturbances will be amplified.

5. There are two kinds of averages. The values of the measurement acting on one interior point of a monad, merely stand for the mean results over a lot of fluid molecules contained in this interior point. The other is the average over all interior points in a monad. The second average results indicate the mean values of the physical quantities over the monad (the standard point) in the global field. By the equiprobability assumption, the mean formula of the second can be given in (20).

6. As the expansion of the fact that the physical functions are continuous in the standard case, the assumption of close property between two infinitely close monads shows that the variation of interior structure as well as values of the physical function $`U(\stackrel{}{x},t,\stackrel{}{x^{}},t^{})`$ with $`(\stackrel{}{x},t)`$ (monad) is continuous too. In other words, there exists only very little difference in characteristics of the function $`U(\stackrel{}{x},t,\stackrel{}{x^{}},t^{})`$ between two infinitely close monads. According to this assumption and the conservation laws, the governing equations of mean, instantaneous and fluctuant quantities in turbulence (i.e., the fundamental equations) are obtained, and the reasonable closure methods are given too. It should be noted that the point-average is used, instead of the Reynolds average, therefore the closure problem is easily overcome in the nonstandard case.

7.The new kind of equations of turbulent motion is presented in this paper. These equations are out of the frame of $`\delta \epsilon `$ and valid for the points, which are permitted to be not uniform. These points, in fact, are monads.